Modeling stress state dependent damage evolution in a cast Al-Si-Mg aluminum alloy

Jump to: navigation, search

AbstractMethodologyMaterial ModelInput DataResultsAcknowledgmentsReferences


Internal state variable rate equations are cast in a continuum framework to model void nucleation, growth, and coalescence in a cast Al-Si-Mg aluminum alloy. The kinematics and constitutive relations for damage resulting from void nucleation, growth, and coalescence are discussed. Because damage evolution is intimately coupled with the stress state, internal state variable hardening rate equations are developed to distinguish between compression, tension, and torsion straining conditions. The scalar isotropic hardening equation and second rank tensorial kinematic hardening equation from the Bammann-Chiesa-Johnson (BCJ) Plasticity model are modified to account for hardening rate differences under tension, compression, and torsion. A method for determining the material constants for the plasticity and damage equations is presented. Parameter determination for the proposed phenomenological nucleation rate equation, motivated from fracture mechanics and microscale physical observations, involves counting nucleation sites as a function of strain from optical micrographs. Although different void growth models can be included, the McClintock void growth model is used in this study. A coalescence model is also introduced. The damage framework is then evaluated with respect to experimental tensile data of notched Al-Si-Mg cast aluminum alloy specimens. Finite element results employing the damage framework are shown to illustrate its usefulness. ®2000 Elsevier Science Ltd. All rights reserved.

Author(s): Mark F. Horstemeyer, J. Lathrop, A.M. Gokhale, M. Dighe

Corresponding Author: Mark F. Horstmeyer

Figure 1. Multiplicative decomposition of the deformation gradient. (click on the image to enlarge).

Figure 2. Schematic of fictitious material with increasing nucleation density and void growth in which the damage framework conceptually comprises. (click on the image to enlarge).

Figure 3. The damage framework encompasses the limiting cases shown by a single void growing in (a) and by damage accumulation from just nucleation in (b). (click on the image to enlarge).

Figure 4. Four representative damage parameters are illustrated as a function of stress triaxiality. Note the saturation that occurs at unity for \phi_2 and \phi_3. (click on the image to enlarge).

Figure 5. Comparison of how void growth and void nucleation relate to each other for the four representative damage parameters. For \phi_1-\phi_3 the void growth has units of volume and the nucleation particle density has units of number of particles per unit volume. For \phi_4, the units for void growth and nucleation are volume fractions and therefore unitless. (click on the image to enlarge).

Figure 6. Micrograph of fracture sites of silicon in tension. (click on the image to enlarge).

Figure 7. Micrograph of fracture sites of silicon in compression. (click on the image to enlarge).

Figure 8. Micrograph of fracture sites of silicon in torsion. (click on the image to enlarge).

Figure 9. The number of voids per unit area versus applied strain comparing the nucleation model and experimental data for tension, compression, and torsion loading conditions. (click on the image to enlarge).

Figure 10. Two different coalescence mechanisms observed in various materials. (click on the image to enlarge).

Figure 11. Damage evolution as a function of strain. The coefficient in the coalescence term is varied to illustrate its influence on total damage. (click on the image to enlarge).

Figure 12. Stress-strain comparison of constitutive model results and test data for compression, tension, and torsion for cast Al-Si-Mg aluminum alloy at ambient temperature and a strain rate of  1 X 10^{-4} /s. (click on the image to enlarge).


Kinematics in the continuum damage mechanics framework

The kinematics follows that in [4,5]. The material point motion is described by elastic straining, inelastic flow, and formation and growth of damage and is illustrated by the multiplicative decomposition of the deformation gradient shown in Fig. 1. The deformation gradient, \underline{F}, is decomposed into the isochoric inelastic, or plastic, (\underline{F}^p_d), dilational inelastic (\underline{F}^p_v), and elastic parts (\underline{F}^e) given by

 \underline{F} = \underline{F}^e\ \underline{F}^p_v\ \underline{F}^p_d.

Eq. (1) assumes that the motion of the body is described by a smooth displacement function. This precludes the initiation of discrete failure surfaces but still allows a continuum description of damage. The elastic deformation gradient,\underline{F^e}, represents lattice displacements from equilibrium. The inelastic deformation gradient,\underline{F}^p_d, represents a continuous distribution of dislocations whose volume preserving motion produces permanent shape changes. The volumetric inelastic deformation gradient, \underline{F}^p_v , represents a continuous distribution of voids causing the volume change of the material from that arises from inelastic deformation. It is assumed to have the form \underline{F}^p_v = \Phi \underline{1} , where \Phi is a function to be determined from kinematics (or conservation of mass).

The Jacobian of the deformation gradient is related to the change in volume or change in density for constant mass as

 J = det\ \underline{F}^p_v\ = \frac{V_2}{V_0} =\frac{\rho_2}{\rho_0}

and must be positive. The change in volume from the reference configuration (State 0) to the intermediate configuration (State 2) is  V_2 = V_0 + V_v assuming that the volume in State 0 equals that in State 1 because of inelastic incompressibility. The volume and density in the reference configuration are given by V_0 and \rho_0, respectively. In transforming the configuration from State 0 to State 2, an added volume from the voids, V_v), is introduced to the total volume, but the volume of the solid matter remains unchanged at its reference value, because the material is unstressed in this configuration. The intermediate configuration in State 2 then designates when elastic unloading has occurred. Damage, \phi), can be defined as the ratio of the change in volume of an element in the elastically unloaded state (State 2) from its volume in the initial reference state to its volume in the elastically unloaded state

 \phi = \frac{V_v}{V_2}.

From this definition, it follows that

 V_0 = (1-\phi)V_2,

where now the Jacobian is determined by the damage parameter, \phi, as

 J = det\ \underline{F}^p_v\ = \frac{1}{1-\phi}.  (5)

Consequently, the restriction that damage is assumed to produce isotropic dilatation gives the volumetric part of the deformation gradient as

 F^p_v = \frac{1}{(1-\phi)^{1/3}}\underline{I},

where  \Phi = (1- \phi)^{-1/3} . The velocity gradient associated with the deformation gradient, \underline{L} = \underline{F}\ \underline{F}^{-1} from Eq. (1) is given by

 \underline{L} = \underline{L}^e\ \underline{L}^p_v\ \underline{L}^p_d,

where  \underline{D} = (1/2)(\underline{L} + \underline{L}^T) and  \underline{W} = (1/2)(\underline{L} - \underline{L}^T) with analogous formulas holding for the elastic, volumetric plastic, and deviatoric plastic parts of the velocity gradients expressed as  \underline{L}^e = \underline{F}^e\ \underline{F}^{-1} ,  \underline{L}^p_v = \underline{F}^p_v\ \underline{F}^{p^{-1}}_v , and  \underline{L}^d_p = \underline{F}^e\ \underline{F}^p\ \underline{F}^{p^{-1}}\ \underline{F}^e\ \underline{F}^{e^{-1}} . The volumetric part of the velocity gradient is then given by

 \underline{F}^p_v\ \underline{F}^{p^{-1}}_v = \frac{\phi}{3(1-\phi)}\ \underline{I} ,

which defines the plastic volumetric rate of deformation as

 \underline{D}^p_v\  = \frac{\phi}{3(1-\phi)}\ \underline{I}.

Also note here that \underline{W}^p_v vanishes. The trace of the volumetric part, Eq. (9), is given as

 tr(\underline{D}^p_v) = \frac{\phi}{(1-\phi)}\ ,

so the damage parameter, \phi, directly relates to the volumetric rate of deformation. The elastic rate of deformation relates to the volumetric rate of deformation by the additive decomposition of the deformation rates similar to Eq. (7),

 \underline{D}^p_v = \underline{D} - \underline{D}^p_v - \underline{D}^p_d.

Similarly, the elastic velocity gradient can be decomposed into components like Eqs. (7) and (11), where the elastic spin equals the total spin when no plastic spin is prescribed. Recall that no volumetric component exists for the spin tensor, that is \underline{W}^p_v = \underline{0},

Now that the rate of deformation related to the damaged state is defined, we can describe damage in terms of void nucleation and void growth in the unstressed intermediate configuration. First, we let N equal the total number of voids in a representative continuum volume V_0 of material in the reference configuration (State 0) and let \eta^* be the number of voids per unit volume in the reference configuration; hence, \eta^* = N/V_0. The average void volume then is v_V = 1/N \sum_{i=1}^\N v_i, where v_i is the void volume from each particle that has nucleated a void. As such the volume of voids is given by

 V_v = \eta^*\ V_0\ v_V.

By combining this definition and inserting them into Eq. (3), the damage parameter, \phi, can be written as

 \phi = \frac{\eta^*\ V_0\ v_V}{V_0 + \eta^*\ V_0\ v_V} = \frac{\eta^*\ v_V}{1 + \eta^*\ v_V}.

This framework for damage was employed in [4]. If the number of voids per unit volume is defined in the intermediate configuration, then

 \phi = \frac{V_v}{V_2} = \frac{V_v}{N} \frac{N}{V_2} = v_V\ \eta ,


 \eta = \frac{N}{V_2} = \frac{N}{V_0} \frac{V_0}{V_2} = \eta^* \frac{V_0}{V_2}.

Recalling that the unstressed intermediate configuration has the volume V_2 = V_0 - V_Vand employing Eq. (3), there prevails

 \eta^* = \frac{\eta}{(1-\phi)}.

The density of voids is counted after the specimen is loaded to a certain strain level and then unloaded. From this point, the specimen is machined and the number counting of voids nucleated is performed representing the elastically unloaded intermediate configuration; hence, \eta is experimentally determined.

The damage framework is shown conceptually in Fig. 2. The number density of voids can change and growth of voids can occur independently or simultaneously. This framework is illustrated by the schematic in Fig. 3 for the limiting cases. One void can grow while others nucleate. A typical void growth model is assumed to have an initial void embryo of a size determined by optical micrographs or some other method. As such, the growth rule applies to both voids that are already present and those that are nucleating. These two types of voids would experience the same void growth in the damage framework. Because the void growth rule is initialized with a positive volume, the nucleated voids are assumed to be of the same size of the initial voids. Perhaps the most realistic embryo size for the newly nucleated site is the size of the second phase particle. The framework allows for this initial choice as well as others. For materials with second phases and pre-existing voids, one would anticipate that the average size of the second phase and average size of the preexisting voids be di€erent. Finally, nucleation is assumed to occur by decohesion of the particle-matrix interface or by particle fracture, and more than one void can be nucleated at a given particle.

Damage has also been defined in several other ways [5-12] as a combination of nucleation, \eta, and growth, v_V, terms. These different damage frameworks and the one used in this study, Eq. (17), will now be compared. For ease of reading, the damage parameters and their corresponding nucleation and growth terms will be identified with correlating numerical subscripts

 \phi_1 = \eta_1\ v_1,

 \phi_2= \frac{\eta_2\ v_2}{1+\eta_2\ v_2},

 \phi_3= 1 - exp(-\eta_3\ v_3),

 \phi_4= \eta_4 + v_4.

These four damage frameworks cover a large class of problems that have been addressed in the literature. The differences related to these damage frameworks have not been shown. Fig. 4 shows the damage evolution as a function of stress triaxiality illustrating the phenomenological differences. Fig. 5 shows the relation between void nucleation and growth in the four damage frameworks when the total damage is unity. \phi_4 is a linear relationship, and the other three are nonlinear to varying degrees. \phi_1 is most often referred to in the materials science literature and is representative of an intermediate configuration representation. As discussed in this section, \phi_2 is the appropriate form for a reference configuration analysis. \phi_3 is motivated from the grain nucleation model developed in [13-15]. The reason for the similarities in \phi_1, \phi_2, and \phi_3 damage frameworks is not obvious at first. When a binomial expansion is substituted in \phi_2 and an exponential expansion is substituted into \phi_3, the similarities are revealed. The two damage frameworks \phi_2 and \phi_3 are essentially \phi_1 with higher order terms.

The damage form, \phi_4 , is often used in the literature to separate the effects of void nucleation and growth [7,8,10-12,16]. In this case \eta_4 and v_4 are assumed to be normalized quantities and hence represent relative volume fractions. As such, \eta_4 is not a number density of voids as described in Eqs. (12)-(15), but is a volume fraction of cracked particles. Also, v_4 is a void volume fraction, not just a void volume. Although \phi_4 can be used to model damage in a phenomenological manner, this purely additive damage form lacks physical motivation in the opinion of the authors.

Material Model

Void nucleation, growth, and coalescence parameters

It follows that, the parameters for the void nucleation, growth, and coalescence terms are determined and explained, but first a description of the material and test specimens is warranted.

The cast A356 aluminum alloy considered in this study work hardenable aluminum matrix with the major second phase being silicon particles in the eutectic region. The aluminum alloy comprises 7% Si, 0.4% Mg, 0.01% Fe, 0.01% Cu, 0.01% Mn, 0.01% Sr, 0.01% Ti, and 0.01% Zn. The material was cast into plates A356 aluminum plates (20 cm x 14 cm x 5 cm) in iron chill molds on the top, bottom, and end of the casting cavity to simulate a permanent mold. A no-bake silica sand was used to create the sides of the plate, the riser, and the down sprue. A ceramic foam filter was used between the down sprue and the riser. A356.2 ingot was melted in an induction furnace. The melt was grain refined with titanium-boron, strontium modified, and degassed using a rotary degasser. The castings were poured between 950 and 977 K, and then cooled over a 16-h period. The plates were removed from the mold and then heat treated to a T6 anneal (solutionized at 810.8 K for 16 h, quenched in hot water at 344 K, and then aged for 4 h at 492.7 K). The microstructure contained aluminum-rich dendrite cells, equiaxed fine silicon particles distributed in the interdendritic regions, and sub-micron intermetallics.

Test specimens were machined from the chill end of the casting, where the amount of microporosity was measured to be low (<0.1% volume fraction), to determine the number of voids nucleated under tension, compression, and uniaxial cyclic boundary conditions. Test specimens were stopped at different strain levels, and then the specimens were sectioned in two pieces so that image analysis could be performed to quantify the void nucleation density. The resulting typical cast microstructure is shown in Fig. 6. Observe that the microstructure contains aluminum-rich dendrite cells, equiaxed fine silicon particles distributed in the interdendritic regions, sub-micron size Mg_2Si precipitates that are responsible for precipitation hardening (not observed by optical microscopy), and numerous Fe and Cu based minor intermetallic phases. After testing, the specimens were sectioned along the vertical plane containing the loading direction (torsional axis in case torsion test specimens), and metallographically polished using standard procedures.

Void nucleation

A series of tension, compression, and torsion experiments were performed to different strain levels to quantify void nucleation evolution of this cast A356 aluminum alloy. Figs. 6-8 show representative optical micrographs of specimens strained under tensile, compressive, and torsional loading conditions. In each micrograph, the loading direction is parallel to the height of the micrograph (see arrow). Observe that fractured Si particles are clearly observed under all three loading conditions. Further, the majority of the cracks-voids in the tensile test specimen (Fig. 6) are perpendicular to the loading direction, whereas the majority of the cracks-voids in the compression test specimen (Fig. 7) are parallel to the loading direction. In the torsion test specimen (Fig. 8), the cracks-voids are observed in all directions. In each specimen, the broken/debonded Si particles were counted, and their sizes, and orientations were measured by using interactive digital image analysis. These measurements were performed on more than 100 systematic random fields of view at 500 x in each specimen to obtain statistically reliable particle damage data. From these data, the fraction of damaged Si particles, their average size, and orientation distribution were computed.

The void nucleation rule of [17] is used to model the results from the cast A356 aluminum data under compression, tension, and torsion. The integrated form of the void nucleation rate equation is given by

 \eta(t) = C_{coeff} exp \left ( \frac{\epsilon(t)d^{1/2}}{K_{IC}f^{1/3}} \left \{ a \left [ \frac{4}{27} - /frac{J^2_3}{J^3_2} \right ] + b \frac {J_3}{J^{3/2}_2} + c \left \Vert \frac{I_1}{\sqrt{J_2}} \right \Vert \right \} \right ) ,

where \eta(t) is the void nucleation density, \epsilon(t) the strain at time t, and C_{coeff} is a material constant. The material parameters a, b, and c relate to the volume fraction of nucleation events arising from local microstresses in the material. These constants are determined experimentally from tension, compression, and torsion tests in which the number density of void sites is measured at different strain levels. The stress state dependence on damage evolution is captured in Eq. (21) by using the stress invariants denoted by I_1, J_2, and J_3, respectively. I_1 is the first invariant of stress (I_1 = \sigma_{kk}). J_2 is the second invariant of deviatoric stress (J_2 = \frac{1}{2} S_{ij}S_{ij}), where S_{ij} = \sigma_{ij} - \frac{1}{3} \delta_{ij}\sigma_{kk} . J_3 is the third invariant of deviatoric stress (J_3 = \frac{1}{3}tr(S_{ij})^3). The rationale and motivation for using these three invariants of stress is discussed in [17]. The volume fraction of the second phase material is f, the average silicon particle size is d, and the bulk fracture toughness is K_{IC}.

For the cast A356 aluminum alloy in our study, K_{IC} = 17.3 MPa m^{0.5}, d = 6 \mu m, and f = 0.07. The volume fraction and average size were determined from optical images of the sectioned test specimens. Fracture toughness tests were performed to determine K_{IC}. The stress state parameters were determined to be a = 615.46 GPa, b = 58.64 GPa, c = 30.011 GPa, and C_{coeff} = 90 GPa. In tension, compression, and torsion, specimens were strained to various levels and then stopped. The number of damaged sites were then counted. The details of determining the stress state parameters are given in Appendix A. Fig. 9 shows the nucleation model compared to the experimental results for compression, tension and torsion, respectively. Note that torsion incurred the highest number of voids nucleated followed by tension and then by compression. It is emphasized here that compression does indeed induce fracture sites in which damage accumulates.

Void growth

A crucial feature in determining the damage state, besides nucleation of voids, is void growth. Many void growth rules have been developed and studied [1,10] but none can comprehensively capture different levels of stress triaxialities, different hardening rates, different strain rates, and different temperature regimes. The damage framework allows for different void growth rules to included and evaluated. In the present study, we use the void growth rule [2]

 v_V = \frac{4}{3} \left ( R_0 exp \left [ \epsilon(t) \frac{\sqrt{3}}{2(1-n)} \times sinh \left ( \sqrt{3} (1-n) \frac{\sqrt{2}I_1}{3\sqrt{J_2}} \right ) \right ] \right )^3

In Eq. (22) the void volume grows as the strain and/or stress triaxiality increases. The material constant n is related to the strain hardening exponent and is determined from the tension tests. R_0 is taken to be the initial radius of the voids. As with most void growth models, the McClintock model allows voids to grow in tension, but not in compression or torsion. This complies with physical observations from measurements of this cast Al-Si-Mg aluminum alloy.


Another item related to damage is the phenomenon of void coalescence. Coalescence is the joining of voids either at the microscale or macroscale and has been observed to occur by two main mechanisms. The first mechanism [18] occurs when two neighboring voids grow together until they join as one, that is, as the ligament between them necks down to a point as illustrated in Fig. 10. Another mechanism occurs when a localized shear band occurs between neighboring voids [19,20], often referred to as the "void sheet" mechanism also shown in Fig. 10.

In the damage framework described in Eqs. (12)-(15), coalescence is restricted to the case where two voids grow together into one and no void sheet occurs. It arises naturally with the multiplicative relation between the nucleation and growth terms in Eq. (18). As Fig. 10 demonstrates, we start with two voids that are nucleated and independently grow until they join together. Then, one void emerges as they coalesce together. The coalescence event causes a discontinuous jump in the nucleation evolution and growth evolution but allows for continuous growth of total damage evolution, \phi. Although discontinuities occur in discrete regions for the nucleation and growth rules, the rate equations evolve as internal state variables at a higher length scale in the continuum where their effects are observed on macroscale effective quantities and thus are continuous functions.

In a phenomenological manner, we include a coalescence term in Eq. (17) as

 C = 1 + f(\eta, v_V),
 \phi = C \eta v.

In the limiting case when the function f(\eta,v_V) =0 in Eq. (23) equals zero, simple coalescence occurs and Eq. (18) results. When some function is used for f(\eta,v_V), microvoid linking is said to have occurred and the rate of damage is increased. The parametric trend of this effect is shown in Fig. 11. For implementation into the constitutive relations, we assume f(\eta,v_V) = C_{26}\eta v_V. It has been observed [1,21] that the microvoid sheet mechanism is related to particles initiating small voids in between two larger voids as the larger voids impose their influence on the surrounding region. As such, coalescence is a function of both nucleation and void growth. Obviously, other forms can be used depending on the material being analyzed.


The BCJ internal state variable plasticity model is modified to account for stress state dependent damage evolution. The pertinent observable and internal state variables are shown below in their rate form as

\stackrel{\circ}{\underline{\sigma}} & = \stackrel{\bullet}{\underline{\sigma}} - \underline{W}^e \underline{\sigma} - \underline{\sigma} \underline{W}^e \\
& = \lambda(1-\phi)tr(\underline{D}^e) \underline{I} + 2\mu(1-\phi)\underline{D}^e - \frac{\stackrel{\bullet}{\phi}}{1-\phi} \underline{\sigma}, \\

 \underline{D}^e = \underline{D} - \underline{D}^p_d - \underline{D}^p_v,

 \underline{D}^p_d = \sqrt{\frac{3}{2}} f(T) \times sinh \left [ \frac{ \sqrt{\frac{3}{2}} \left \Vert \underline{\sigma}^' - \sqrt{\frac{2}{3}}\underline{\alpha} \right \Vert - \{ R + Y(T) \} \{ 1 - \phi \} } {V(T) \{ 1 - \phi \} } \right ] \times \frac{\underline{\sigma}^' - \sqrt{\frac{2}{3}} \underline{\alpha} } { \left \Vert \underline{sigma}^' - \sqrt{\frac{2}{3}} \underline{\alpha} \right \Vert } ,

\stackrel{\circ}{\underline{\alpha}} & = \stackrel{\bullet}{\underline{\alpha}} - \underline{W}^e \underline{\alpha} - \underline{\alpha} \underline{W}^e \\
& = h(T) \underline{D}^p_d - \left [ \sqrt{\frac{2}{3}} r_d (T) \left \Vert \underline{D}^p_d \right \Vert + r_s(T) \right ] \sqrt{\frac{2}{3}} \left \Vert \underline{\alpha} \right \Vert \underline{\alpha}, \\

 R = H(T) \sqrt{\frac{2}{3}} \underline{D}^p_d - \left [ \sqrt{\frac{2}{3}} R_d(T) \left \Vert \underline{D}^p_d \right \Vert + R_s(T) \right ] R^2,

 \underline{W}^e = \underline{W} - \underline{W}^p

and are generally written as objective rates (\stackrel{\circ}{\underline{\sigma}},\stackrel{\circ}{\underline{\alpha}}) with indifference to the continuum frame of reference assuming a Jaumann rate in which the continuum spin equals the elastic spin (\underline{W} = \underline{W}^e). The elastic Lame` constants are denoted by \lambda and \mu. The elastic rate of deformation (\underline{D}^e) results when the total deformation (\underline{D}), which is defined by the boundary conditions, is subtracted from the deviatoric and volumetric components of the flow rule. The deviatoric plastic flow rule, \underline{D}^p_d, is a function of the temperature, the kinematic hardening internal state variable (\underline{\alpha}), the isotropic hardening internal state variable (R), the volume fraction of damaged material (\phi), and the functions f(T), V(T), and Y(T), which are related to yielding with Arrhenius-type temperature dependence. The function Y(T) is the rate-independent yield stress. The function f(T) determines when the rate dependence affects initial yielding. The function V(T) determines the magnitude of rate dependence on yielding. These functions are determined from simple isothermal compression tests with different strain rates and temperatures

 V(T) = C_1 exp \left ( \frac{-C_2}{T} \right ), Y(T) = C_3 exp \left ( \frac{C_4}{T} \right ), f(T) = C_5 \left ( \frac{-C_6}{T} \right ).

The kinematic hardening internal state variable, \underline{\alpha}, reflects the effect of anisotropic dislocation density, and the isotropic hardening internal state variable R, reflects the effect of the global dislocation density. As such, the hardening equations (24) and (25) are cast in a hardening-recovery format that includes dynamic and static recovery. The functions r_s(T) and R_s(T) are scalar in nature and describe the diffusion-controlled static or thermal recovery, while r_d(T) and R_d(T) are scalar functions describing dynamic recovery. The anisotropic hardening modulus is h(T), and the isotropic hardening modulus is H(T).

The hardening moduli and dynamic recovery functions account for deformation-induced anisotropy arising from texture and dislocation substructures by means of stress dependent variables. It was shown [22] that by using J_3 in the hardening equations the different hardening rates between axisymmetric compression and torsion (torsional softening) were accurately captured. Included in [23] is the feature of [24] that distinguishes between compression and torsion. Because damage evolves differently under compression, tension, and torsion for this cast Al-Si-Mg aluminum alloy, the stress arising from these different loading conditions is different. As such, the equations for the hardening and recovery variables are expressed as:

 r_d = \left \{ C_7 \left ( 1 - C_{19} \left [ \frac{4}{27} - \frac{J^2_3}{J^3_2} \right ] - C_{20} \frac{J_3}{J^{1.5}_2} \right ) \right \} \times exp \left ( \frac{-C_8}{T} \right ),

 H = \left \{ C_9 \left ( 1 + C_{19} \left [ \frac{4}{27} - \frac{J^2_3}{J^3_2} \right ] + C_{20} \frac{J_3}{J^{1.5}_2} \right ) \right \} \times exp \left ( \frac{-C_{10}}{T} \right ),

 r_s(T) = C_{11} exp \left ( \frac{-C_{12}}{T} \right ),

 r_d = \left \{ C_13 \left ( 1 - C_{19} \left [ \frac{4}{27} - \frac{J^2_3}{J^3_2} \right ] - C_{20} \frac{J_3}{J^{1.5}_2} \right ) \right \} \times exp \left ( \frac{-C_{14}}{T} \right ),

 H = \left \{ C_{15} \left ( 1 + C_{19} \left [ \frac{4}{27} - \frac{J^2_3}{J^3_2} \right ] + C_{20} \frac{J_3}{J^{1.5}_2} \right ) \right \} \times exp \left ( \frac{-C_{16}}{T} \right ),

 R_s(T) = C_{17} exp \left ( \frac{-C_{18}}{T} \right ).

Because attention is focused on stress state dependent damage in this study, we will reduce the BCJ equations to quasi-static strain rates for examining the phenomena at ambient temperature. As such, many of the BCJ constants will be zero. The constant values are given in Appendix A for the Al-Si-Mg cast material discussed in this paper.

Finally, a note here about continuum thermodynamics is in order since the BCJ model is motivated from internal state variable (ISV) theory [25]. Essentially, the ISVs reflect effects from lower length scale microstructural causes. In this model the ISV rate equations (28) and (29) and time rate of change of Eq. (17) are functions of the observable variables (temperature, stress state, and rate of deformation). In general, the rate equations of generalized displacements, or thermodynamics fluxes, describing the rate of change may be written as independent equations for each ISV or as derivatives of a suitably chosen potential function arising from the hypothesis of generalized normality [26]. An advantage of assuming generalized normality, although somewhat restrictive, is unconditional satisfaction of the Kelvin inequality of the second law of thermodynamics (nonnegative intrinsic dissipation), i.e.,

 \underline{\sigma} : \underline{D}^p_d - \underline{b} : \stackrel{\circ}{\underline{\alpha}} - \kappa \cdot \stackrel{\bullet}{R} - \Phi \cdot \stackrel{\bullet}{\phi} \geqq 0.

Here, the backstress, \underline{b}, is the thermodynamic force related to kinematic hardening; the isotropic stress, \kappa, is the thermodynamic force related to isotropic hardening; and the energy release rate, \Phi, is the thermodynamic force related to the kinematic damage variable, \phi. The damage ISV encompasses dissipation according to Eq. (38) from the void nucleation, growth, and coalescence terms.

Input Data

Model implementation into finite element code

A few comments are warranted in regard to the implementation of the plasticity-damage model. When damage \phi approaches unity, failure is assumed to occur. The goal is to implement the damage framework into a finite element code for solving complex boundary value problems, so failure occurs as \phi\rightarrow 1.0 within an element. Damage accumulation less than unity would be designated as failed material by engineers. In fact, it was stated [27] that a damage level of 50% is the limitation on the degraded elastic moduli. Practically, the total damage for final failure should perhaps be even less than 50%, but in applications [28] using the void growth rule [29,30], the damage goes rapidly to unity just after a few percent void volume fraction. More complicated functions for the damaged elastic moduli can be used such as that in [31], but these are computationally expensive.

To implement the model in a finite element code, we replace the deviatoric plastic rate of deformation with the total rate of deformation in the recovery terms of the hardening rate equations (28) and (29). This substitution makes the hardening rate equations directly solvable. Eqs. (28) and (29) become:

 \underline{\alpha} = h(T) \underline{D}^p_d - \left [ \sqrt{2/3}r_d(T) \left \Vert \underline{D} \right \Vert + r_s(T) \right ] \sqrt{2/3} \left \Vert \underline{\alpha} \right \Vert \underline{\alpha},
 R =  \sqrt{2/3} H(T) \left \Vert \underline{D}^p_d \right \Vert - \left [ \sqrt{2/3}R_d(T) \left \Vert \underline{D} \right \Vert + R_s(T) \right ] |R| R.

At the beginning of each time step, determine the values for J_2 and J_3 from the previous step to modify r_d, R_d, h, and H. Employed then is a radial return method to determine the plastic part of the strain by assuming the strain to be all elastic (i.e. \underline{D}^p_d =0 ). This gives the following trial values for the deviatoric stress and internal hardening variables:

\underline{\sigma}^*_{n+1}  & = \underline{\sigma}^'_{n} + \int\limits_{t_n}^{t_{n+1}} 2\mu(1-\phi) \underline{D}^' \, dt - \int\limits_{t_n}^{t_{n+1}} \frac{\phi \underline{\sigma}^'}{1+\phi} \, dt\\
& \approx \underline{\sigma}^'_n (1- \frac{\phi \Delta t}{1+\phi}) = 2\mu (1-\phi) \underline{D}^' \Delta t, \\

\underline{\alpha}^*_{n+1}  & = \underline{\alpha}_{n} - \int\limits_{t_n}^{t_{n+1}} (r_d \sqrt{2/3} \left \Vert \underline{D} \right \Vert + r_s )\sqrt{2/3} \left \Vert \underline{\alpha} \right \Vert \underline{\alpha} \, dt \\
& \approx \{1 - (r_d \sqrt{2/3} \left \Vert \underline{D} \right \Vert + r_s)\sqrt{2/3}\left \Vert \underline{\alpha} \right \Vert \Delta t \} \underline{\alpha}_n  , \\

R^*_{n+1}  & = R_{n} - \int\limits_{t_n}^{t_{n+1}} (R_d \sqrt{2/3} \left \Vert \underline{D} \right \Vert + R_s ) |R| r  \, dt \\
& \approx \{1 - (R_d \sqrt{2/3} \left \Vert \underline{D} \right \Vert + R_s)|R|\Delta t \}R_n. \\

Note that a differential equation of the form R = -kR has a solution that tends to zero, but if we take too large a time step, our approximate integration could oscillate about zero. To avoid this behavior, we limit the terms multiplying (al) and R_n in Eqs. (42) and (43) to be greater than or equal to zero.

Define the tensor \underline{\xi} = \underline{\sigma^*} - (2/3) \underline{\alpha}, such that the flow rule can be written as

 \underline{D}^p_d = \sqrt{3/2}f(T) \times sinh \left ( \frac{ \sqrt{3/2} \left \Vert \xi \right \Vert - (R+Y(T))(1-\phi)}{V(T)(1-\phi)} \right ) \times \frac{\xi}{\left \Vert \xi \right \Vert}.

By taking the norm of both sides, Eq. (44) can be inverted to give

 \Phi = \sqrt{\frac{3}{2}} \left \Vert \underline{\xi} \right \Vert - (1-\phi) \left [ R+Y+Vsinh^{-1} \left ( \frac{ \sqrt{2/3} \left \Vert \underline{D} \right \Vert }{f} \right ) \right ] = 0

If on evaluation of \Phi we find \phi \leqq 0, then the elastic assumption is valid. Use the trial values of \underline{\sigma}{\alpha}, and R as the actual values. Otherwise, look for a deviatoric plastic strain component \underline{D} such that

 \int\limits_{t_n}^{t_{n+1}} D^p \, dt = \frac{\gamma}{\left \Vert \xi \right \Vert} \xi, \int\limits_{t_n}^{t_{n+1}} \left \Vert D^p \right \Vert \, dt = \gamma.

This leads to corrections to the trial values of:

\sigma^'_{n+1}  & = \sigma^*_{n+1} - \int\limits_{t_n}^{t_{n+1}} 2\mu(1-\phi)D^p\, dt \\
& = \sigma^*_{n+1} - \frac{2\mu(1-\phi)\gamma}{\left \Vert \xi \right \Vert} \xi , \\

\alpha_{n+1}  & = \alpha^*_{n+1} + \int\limits_{t_n}^{t_{n+1}} h(T)D^p \, dt \\
& = \alpha^*_{n+1} + \frac{h\gamma}{\left \Vert \xi \right \Vert} \xi , \\

R_{n+1}  & = R^*_{n+1} + \sqrt{2/3} \int H(T) \left \Vert D^p \right \Vert \, dt \\
& = R^*_{n+1} + \sqrt{2/3} H \gamma . \\

Substituting these corrected values back into the inverted flow rule, it is easy to show that \Phi = 0 is satisfied by choosing \gamma as:

 \gamma = \frac{ \left \Vert \xi \right \Vert - \sqrt{2/3} (1-\phi)(Y+R+Vsinh^{-1}(\left \Vert D \right \Vert /f ))}{2\mu(1-\phi)+ \frac{2}{3} (h+H(1-\phi))}.

Eq. (49) is then used to correct the trial values. We then calculate the total effective strain as

 \epsilon_{N+1} = \epsilon_N = \sqrt{2/3} \gamma.

At this point we can calculate J_2 and J_3 from the corrected \sigma^'_{n+1} and calculate a new damage term \phi from updated nucleation, void growth, and coalescence from Eqs. (21)-(24). The updated damage equations become:

 \eta_{N+1} = C_{coeff} exp\left ( \frac{\epsilon_{N+1}d^{1/2}}{K_{IC}f^{1/3}} \left \{ a \left [ \frac{4}{27} - \frac{J^2_3}{J^3_2} \right ] + b \frac{J_3}{J^{3/2}_2} + c \left \Vert \frac{I_1}{\sqrt{J_2}} \right \Vert \right \} \right ),

 v_{N+1} = \frac{4}{3} \left ( R_0 exp \left [ \epsilon_{N+1} \frac{\sqrt{3}}{2(1-n)} \times sinh \left ( \sqrt{3}(1-n) \frac{\sqrt{2I_1}}{3\sqrt{J_2}} \right ) \right ] \right )^3,

 C_{N+1} = 1 + C_{26}\eta_{N+1}v_{N+1},

 \phi_{N+1} = \eta_{N+1}v_{N+1}C_{N+1}.

Finally, we add the pressure term to update the total stress as

 \sigma_{n+1} = \sigma^'_{n+1} + p_{n+1},


 p_{n+1} = \frac{1}{3} tr(\sigma_n)(1-\phi) + (1-\phi)\Delta t K tr(D).

Fig. 12 shows a comparison of the model to stress-strain data from compression, tension, and torsion tests at ambient temperature and quasistatic loading conditions (1 \times 10^{-4}/s). These curves reflect the inclusion of the void nucleation, growth, and coalescence terms. The hardening rate differences arising from these different global stress states is driven more by void nucleation than the growth and coalescence. It can be seen from Fig. 9 that the void nucleation rate is increasing as you go from compression to tension to torsion. The void nucleation relaxes the local dislocation density around the particles to relieve the local microstresses. As such, the global hardening rate is increasing as you go from torsion to tension to compression, the reverse order of the void nucleation rate.

A comment should be made regarding the plotting of the nucleation density of voids in a finite element code. First, physical measurements of the nucleation density is a number count per unit area, so in 2D finite element calculations, the void growth area must be used. Also, if a total number of voids nucleated is desired as opposed to a number density, then the area of the element must be included with the corresponding units conversion. The current measure is a number count per millimeter squared.


Model comparison to notch Bridgman tensile data

Finite element analyses were performed of notch tensile specimens employing the modified BCJ plasticity-damage model described in the previous sections. Quarter "plate" axisymmetric analyses using the finite element code ABAQUS were performed to study the effects of void nucleation, growth, and coalescence.

The typical progression of damage evolution for wrought materials in notch tensile tests arises from void growth at the center of the specimen because of the high stress triaxialities. Actually, stress triaxiality starts at the notch edge at the start of the deformation because of the stress concentration but moves fairly rapidly to the center as deformation proceeds. The stress triaxiality is nonuniform throughout the cross section of the specimen and reaches a level sufficient to start void growth when its peak reaches the center of the notch specimen. Fig. 13 shows contour plots of the stress triaxiality (hydrostatic stress/deviatoric stress) of a notch tensile test illustrating spatial movement of the peak stress triaxiality. Again, this is typical for wrought, ductile materials.

For ductile materials with brittle phases, such as the cast A356 aluminum alloy examined in this study, the final failure location does not occur at the center of the specimen. As the peak stress triaxiality increases in magnitude and moves toward the center of the specimen, voids have been nucleated by the fracture of the particles or by debonding of the particle-aluminum interface. As such, the damage increases by not only voids growing but by new voids initiating and then growing. Damage reaches a level of unity (unity is interpreted as defining a hole in the material the size of the finite element) before it reaches the center of the specimen like a wrought alloy. In a wrought alloy at first element failure, the highest level of plastic strain occurs at the notch edge, but the highest stress triaxiality occurs at the center of the specimen. As such, it can be concluded that the stress triaxiality drives the void growth as opposed to the plastic deformation, because damage evolves mainly at the center of the specimen. For this cast aluminum, the highest plastic strain occurs at the notch edge and the stress triaxiality occurs away from the edge similar to wrought materials, but Fig. 14 illustrates that the highest damage (which comprises void nucleation, growth, and coalescence) occurs at two different locations. The two different locations are aligned with the peak plastic strain and peak stress triaxiality as shown in Figs. 15 and 16, respectively. In Fig. 15, the dark contour shows that the plastic strain reached a level of 42% at the notch edge and around only 1% near notch center. In Fig. 16, the lightest contour shows that the highest negative tensile pressure (stress triaxiality equals the negative pressure over the deviatoric stress) occurs between the notch edge and specimen center. Fig. 17 shows contour levels of void nucleation in which the dark contour shows the peak level occurring near the highest at the notch edge and near the peak stress triaxiality. Fig. 18 shows contours of void volume throughout the specimen. Note in Fig. 18 that the peak void growth occurs at the peak stress triaxiality.

Coalescence was assumed to be unity for the above finite element analysis. Other analyses were performed in which the coalescence term was increased. The results essentially showed that as the coalescence term increased, the total damage level reached unity closer to the notch edge. The interpretation is that when coalescence increases, the void microlinking increases. This in turn makes the response look more brittle on a macroscale, since the more ductile the material, the closer to the center of the specimen would final failure occur.

Figure 13. Progression of stress triaxilities (hydrostatic stress/deviatoric stress) over the spatial domain of a notch specimen under tensile loading. (click on the image to enlarge).
Figure 14. Total damage contours just before final failure of a notch tensile specimen with a notch acuity to radius ratio of 0.117 for a cast A356 aluminum alloy. Note that the dark region indicates the peak damage level in two locations along the radial plane of symmetry. (click on the image to enlarge).
Figure 15. Plastic strain contours just before final failure of a notch tensile specimen with a notch acuity to radius ratio of 0.117 for a cast A356 aluminum alloy. Note that the dark region indicates the peak plastic strain level near the notch edge. (click on the image to enlarge).
Figure 16. Pressure contours just before final failure of a notch tensile specimen with a notch acuity to radius ratio of 0.117 for a cast A356 aluminum alloy. Note that the light region indicates the peak tensile hydrostatic stress between the notch center and notch edge. (click on the image to enlarge).
Figure 17. Void nucleation contours just before final failure of a notch tensile specimen with a notch acuity to radius ratio of 0.117 for a cast A356 aluminum alloy. Note that the dark region indicates the peak nucleation level at two locations. (click on the image to enlarge).
Figure 18. Void growth contours just before final failure of a notch tensile specimen with a notch acuity to radius ratio of 0.117 for a cast A356 aluminum alloy. Note that the dark region indicates the peak void growth level between the notch center and notch edge. (click on the image to enlarge).


The support and interest of USCAR/USAMP Lightweight Metals Group under the leadership of Richard Osborne is appreciated. Gerry Shulke provided the cast plates and Westmoreland Mechanical Testing and Research performed the tests. The work by M.F. Horstemeyer and J. Lathrop was performed at the US Department of Energy, Sandia National Laboratories under contract DE-AC04-94AL85000. The authors thank Dan Mosher, Doug Bammann, Ganesh Ramaswamy, Dave McDowell, and Jinghong Fan for their comments regarding this paper.


[1] W.M. Garrison, N.R. Moody, J. Phys. Chem. Solids 48 (1987) 1035-1074.

[2] F.A. McClintock, A criterion for ductile fracture by the growth of holes, ASME J. Appl. Mech. 35 (1968) 363-371.

[3] D.J. Bammann, M.L. Chiesa, G.C. Johnson, Modeling large deformation and failure in manufacturing processes, in: T. Tatsumi, E. Wannabe, T. Kambe (Eds.), Theoretical and Applied Mechanics, Elsevier, Amsterdam, 1996, pp. 359-376.

[4] L. Davison, A.L. Stevens, M.E. Kipp, Theory of spall damage accumulation in ductile metals, J. Mech. Phys. Solids 25 (1977) 11-28.

[5] D.J. Bammann, E.C. Aifantis, A damage model for ductile metals, Nucl. Eng. Des. 116 (1989) 355-362.

[6] T.W. Barbee, L. Seaman, R. Crewdson, D. Curran, Dynamic fracture criteria for ductile and brittle metals, J. Mater., JMLSA 7 (3) (1972) 393-401.

[7] A.L. Gurson, Continuum theory of ductile rupture by void nucleation and growth - 1. Yield criteria and flow rules for porous ductile media, J. Eng. Mater. Technol. 99 (1977) 2-15.

[8] A.L. Gurson, Porous rigid-plastic materials containing rigid inclusions-yield function, plastic potential, and void nucleation, in: D.M.R. Taplin (Ed.), Proc. Int. Conf. Fracture, vol. 2A, 1977, pp. 357-364. 46 M.F. Horstemeyer et al. / Theoretical and Applied Fracture Mechanics 33 (2000) 31-47

[9] E.P. Chen, M.E. Kipp, D.E. Grady, A strain-rate sensitive rock fragmentation model, in: K.P. Chong, J. Ward Smith (Eds.), Mechanics of Oil Shale, Elsevier, New York, 1984, pp. 423-456.

[10] V. Tvergaard, Material failure by void growth to coalescence, in: J.W. Hutchinson, T.Y. Wu (Eds.), Advances in Applied Mechanics, vol. 27, 1990, pp. 83-151.

[11] J.B. Leblond, G. Perrin, J. Devaux, An improved Gursontype model for hardenable ductile metals, Eur. J. Mech. A 14 (4) (1995) 499-527.

[12] E. Marin, A critical study of finite strain porous inelasticity, PhD Thesis, Georgia Institute of Technology, 1993.

[13] A.N. Kolmogorov, Bull. Acad. Sci. USSR, No. 3 (1937).

[14] M. Avrami, Kinetics of phase change, J. Chem. Phys. 7 (1939) 1103-1112.

[15] W.A. Johnson, R.F. Mehl, Processes of nucleation and growth, Trans. AIME (Iron & Steel Div.) 135 (1939) 416-442.

[16] D.L. McDowell, E. Marin, C. Bertoncelli, A combined kinematic-isotropic hardening theory for porous inelasticity of ductile metals, Int. J. Damage Mech. 2 (1993) 137-161.

[17] M.F. Horstemeyer, A.M. Gokhale, A void nucleation model for ductile materials, Int. J. Solids Struct., accepted.

[18] R. Garber, I.M. Bernstein, A.W. Thompson, Scripta Metall. 10 (1976) 341.

[19] T.B. Cox, J.R. Low Jr, An investigation of the plastic fracture of AISI 4340 and 18 Nickel-200 grade maraging steels, Metall. Trans. 5 (1974) 1457-1470.

[20] H.C. Rogers, The tensile fracture of ductile metals, Trans. AIME 218 (1960) 498-506.

[21] P.E. Magnusen, E.M. Dubensky, D.A. Koss, The effect of void arrays on void linking during ductile fracture, Acta Metall. 36 (6) (1988) 1503-1509.

[22] M.P. Miller, D.L. McDowell, The effect of stress state on the large strain inelastic deformation behavior of 304L stainless steel, J. Eng. Mater. Technol., Trans. ASME 118 (1) (1996) 28-36.

[23] M.F. Horstemeyer, D.L. McDowell, D.J. Bammann, Torsional softening and the forming limit diagram, in: K. Chen (Ed.), Analysis of Autobody Stamping Technology, SAE Tech. Ser. 960599, SP1134, 1995, pp. 81-89.

[24] D.J. Bammann, Modelling the temperature and strain rate dependent large deformation of metals, Appl. Mech. Rev. 43 (5) Part 2 (1990) 312-319.

[25] B.D. Coleman, M.E. Gurtin, Thermodynamics with internal state variables, J. Chem. Phys. 47 (1967) 597-613.

[26] J.R. Rice, Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity, J. Mech. Phys. Solids 9 (1971) 433-455.

[27] B. Budiansky, Thermal and thermoelastic properties of isotropic composites, J. Composite Mater. 4 (1970) 286-295.

[28] D.J. Bammann, M.L. Chiesa, M.F. Horstemeyer, L.I. Weingarten, Failure in ductile materials using finite element methods, in: N. Jones, T. Weirzbicki (Eds.), Structural Crashworthiness and Failure, Elsevier, Amsterdam, 1993, pp. 1-54 (Chapter 1).

[29] A.C.F. Cocks, M.F. Ashby, Intergranular fracture during power law creep under multixial stresses, J. Metal Sci. (1980) 395-402.

[30] A.C.F. Cocks, M.F. Ashby, On creep fracture by void growth, Prog. Mater. Sci. 27 (1982) 189-244.

[31] Y.H. Zhao, G.P. Tandon, G.J. Weng, Elastic moduli for a class of porous materials, Acta Mechanica 76 (1989) 105-130.

Personal tools

Material Models