Gurson-based modelling of ductile damage and failure during cyclic loading processes at large deformation

Engineering Fracture Mechanics 160 (2016) 95–123

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Gurson-based modelling of ductile damage and failure during cyclic loading processes at large deformation D. Klingbeil a,⇑, B. Svendsen b, F. Reusch c a

Bundesanstalt für Materialforschung und -prüfung (BAM), Berlin, Germany RWTH Aachen University, Aachen, Germany c Beck Engineering Pty Ltd, Berlin, Germany b

a r t i c l e

i n f o

Article history: Received 21 September 2015 Received in revised form 19 February 2016 Accepted 14 March 2016 Available online 9 April 2016 Keywords: Ductile damage Combined non-linear hardening with cyclic loading Finite-element simulation Material parameter identification Residual stresses

a b s t r a c t Purpose is the formulation, numerical implementation, identification and application of a material model for ductile damage and failure during cyclic and non-proportional loading. The authors combined a hyperelasticity-based elasto-plastic model for non-linear isotropic as well as kinematic hardening with a modified Gurson model. Evolution strategy helped identify the model parameters for the high-strength steel 10MnMoNi5-5. The simulation of ductile failure in fracture mechanics specimens verified the model with respect to cyclic loading at two temperatures. The simulation of additional fracture mechanics applications validated the model as to the development of residual stresses at the crack tip under cyclic loads. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction The modelling of inelastic material behaviour at geometric imperfections such as cracks or notches requires in general the consideration of loading paths which may be locally non-proportional or cyclic even when the global loading is monotonic. Indeed, deformation localisation in tensile specimens, as well as crack tip blunting in fracture specimens, can cause a loading path change or reversal at a material point. As such, modelling of ductile damage and failure of fracture mechanics specimens and engineering structures should incorporate kinematic hardening effects. In the current work, the authors carried this out in the context of a modification of the Gurson–Tvergaard–Needleman (GTN) model [1–4]. Extensions of the original Gurson model [5] to the case of combined hardening have been proposed by [6], and to the case of non-spherical voids and combined hardening, by [7]. Lievers et al. [8] presented a formulation of kinematic hardening following [9] involving a modified GTN yield function. Likewise, [10] formulated an extension of the Gurson model to cyclic loading conditions based on the Chaboche–Lemaitre thermodynamical approach. Recently, a review about techniques of modelling ductile damage [11] stated the lack of a combination between ductile damage and kinematic hardening and its application to fracture mechanics specimens. The current work introduces a model that represents such an extension to the GTN model for void nucleation, growth and coalescence to the case of non-monotonic and non-proportional loading at large deformation [12,13]. It follows the approach of [7] and for simplicity, only the matrix material undergoes kinematic hardening. For monotonic and proportional loading, the model reduces to the original GTN formulation.

⇑ Corresponding author. E-mail address: [email protected] (D. Klingbeil). http://dx.doi.org/10.1016/j.engfracmech.2016.03.023 0013-7944/Ó 2016 Elsevier Ltd. All rights reserved.

96

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Nomenclature c b f f0 fn
f fc ff
fu h q1 , q2 r k z A X I DP E FE FP M RE T UE ln U E

en sn

eM c / k; l

rM rvM rH rA v

hardening parameter, kinematic hardening saturation parameter, kinematic hardening void volume fraction initial void volume fraction volume fraction of void-nucleating inclusions effective void volume fraction, GTN yield function critical value of f for void coalescence critical value of f for material failure
critical value of f for material failure hardening modulus, linear isotropic hardening, unit cell modelling Tvergaard parameters, GTN yield function hardening parameter, isotropic hardening saturation parameter, isotropic hardening accumulated magnitude of the inelastic rate of deformation strain-controlled void nucleation strength back stress tensor identity tensor local inelastic rate of deformation tensor small strain tensor local elastic deformation tensor local inelastic deformation tensor Mandel stress tensor elastic polar rotation tensor small stress tensor, Cauchy stress tensor elastic right stretch tensor elastic right logarithmic stretch tensor mean accumulated equivalent plastic strain for void nucleation standard deviation of distribution about en accumulated equivalent plastic strain in the matrix plastic multiplier Gurson yield function Lamé constants matrix yield stress von Mises equivalent stress hydrostatic stress yield stress for hydrostatic part of current yield condition triaxiality

As in any model identification, it is reasonable to (i) minimise the number of parameters, and (ii) identify as many parameters as possible independently of each other. The procedure described in the current work firstly identifies the hardening parameters (i.e., before damage begins), followed by the identification of the damage parameters at fixed hardening parameters. The two selected temperatures for testing were room temperature 25 °C (RT) and 300 °C. The use of a number of different experimental data sets minimises the non-unicity of the parameter identification. Therefore, the current work used the evolution strategy [14]. In the rear part of the present work, the model simulates the development of residual stresses at the crack tip during cyclic thermomechanical loading. In the case of the warm-pre-stress (WPS) effect for example, the thermomechanical loading history of a component influences the associated fracture behaviour significantly [15–17]. Such components often do not exhibit brittle fracture after cooling because the fracture toughness K IC (e.g. [18]) is larger than the crack resistance determined for this temperature by following the specified standard procedures [19]. This is due to the development of residual stresses at the crack tip. In turn, they come from significant plastic deformation at the crack tip and cause blunting of the initially sharp crack. The current analysis yields information about the plastic deformation and the ductile damage that occurs during the loading process and – after unloading – about the residual stress. The results from simulations are validated by comparison with measurements of residual and crack-opening stresses around the crack tip of a C(T)25 fracture mechanics specimen [15].

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

97

2. Review of the Gurson model and its extensions as to monotonic loading As is well-known from experiments on metal (e.g., [20]), nucleation, growth and coalescence of microscopic voids cause failure mainly in the ductile regime. Initially, the material contains voids determining an initial void growth volume fraction and inclusions determining the initial void nucleation volume fraction. For most materials, void nucleation is due to decohesion or fracture of (generally non-metallic) inclusions in the matrix material (e.g., [21]). On account of stress, the ratio of the hydrostatic to the deviatoric stress, i.e., the stress triaxiality v mainly influences void growth. As loading proceeds, voids continue to grow and begin to coalesce through shear banding between neighbouring voids. Hereby, a second population of microparticles greatly amplifies such shear banding and localisation as, e.g., [22] showed experimentally. Coalescence leads to accelerated ductile damage of the specimen as a whole, and to a macroscopic ductile crack and finally to failure. To formulate a material model that describes the effective behaviour of such metallic materials, Gurson [5] solved the mechanical equilibrium problem for a spherical representative volume element (RVE). In this model an ideal isotropic von Mises elasto-plastic matrix material with constant yield stress rM0 surrounds a single void of fractional volume f subjected to small-deformation loading conditions. In this fashion, he obtained the form




/ðT; f Þ ¼



rvM ðTÞ2 3 rH ðTÞ 2  ð1  f Þ ¼ 0 þ 2f cosh 2 r2M0 r2M0

ð1Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi of the yield function. Here, rvM ðTÞ :¼ 3 devðTÞ  devðTÞ=2 is the usual von Mises equivalent stress measure as based on the (in the RVE average or mean) stress T; devðTÞ :¼ T  rH ðTÞ I its deviatoric part, and rH ðTÞ :¼ I  T=3 the hydrostatic stress. For an undamaged material, i.e., f ¼ 0, (1) describes von Mises flow behaviour. For increasing values of f the von Mises yield surface cylinder in the principal stress space turns towards an ellipsoidal yield surface for the Gurson model. Assuming the matrix material being plastically incompressible, and ignoring any elastic volume changes, Gurson [5] obtained the evolution relation

f_ ¼ ð1  f Þ I  E_ P

ð2Þ

for f in terms of the small-deformation inelastic strain rate E_ P . The original Gurson model was completed by an associated flow rule, i.e.,

E_ P ¼ c /;T ;

ð3Þ

with c as the plastic multiplier. Using these micromechanical results, Gurson [5] then formulated his RVE on the assumption that plastic dissipation takes place only in the matrix material. Experimental observations show that in most metallic materials the final void volume fraction at ductile failure is much lower than the theoretical value of f ¼ 1 predicted by the of original Gurson model (e.g., [21]). This has to do with the fact that the original Gurson formulation [5], and in particular the evolution relation (2), is restricted to void growth and neglects additional processes such as void nucleation and coalescence. Related to the extended deformation-based form of (2) [1] gave an extension of the original Gurson model to the case of void nucleation,

Here,

f_ ¼ f_ growth þ f_ nucleation ¼ ð1  f Þ trðE_ P Þ þ AðeM Þ e_ M :

ð4Þ

( ) fn 1 ðeM  en Þ2 ffi exp  AðeM Þ ¼ pffiffiffiffiffiffiffiffiffiffi 2 s2n 2ps2n

ð5Þ

stands for a material property reflecting the strength of the strain-controlled void nucleation. This form relies on the idea that void nucleation due to decohesion and fracture of hard particles is statistical in nature and correlates with the Gaussian distribution of accumulated equivalent inelastic strain eM in the matrix about a critical value en for nucleation. In addition, f n represents the volume fraction of void-nucleating inclusions, en the mean equivalent plastic strain at which decohesion of the inclusions from the matrix begins and sn its standard deviation. Consistent with the formulation of the original model, the evolution relation for eM holds here,

e_ M ¼

T  E_ P 1 ¼ ð1  f Þ c /; rM ð1  f Þ rM

ð6Þ 1

following from (1) and (3). Note that the dependency on the term ð1  f Þ causes this relation to deviate from being associative. The development of eM also leads of course to isotropic hardening in the matrix material and a non-constant matrix yield stress rM in (1), depending now on eM . To extend the original model to the case of void nucleation, [2,3] utilised results from [23,24] to propose an extension of (1) which included void coalescence

98

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123




/ðT; rM ; f Þ ¼



r2vM ðTÞ 3 rH ðTÞ 2

 ðq1 f ðf ÞÞ ¼ 0  1 þ 2 q1 f ðf Þ cosh q 2 2 rM r2M

ð7Þ

that does not affect its functional form. This is referred to as the GTN (Gurson–Tvergaard–Needlemann) yield function in what follows. Tvergaard [23] introduced two additional parameters q1 and q2 that helped to obtain better agreement with numerical results. Perrin and Leblond [25] confirmed optimal values of q1 ¼ 1:5 and q2 ¼ 1:0 by analytical results: q1 ¼ 4=e  1:47 and q2 ¼ 1:0. The current work assumes these values as well. Other possibilities for q1 are q1 ¼ 1:25 [4] and q1 ¼ 1:35 [7]. The effective damage measure

(




f ðf Þ ¼

f 6 fc

f fc þ

f
u f c f f f c

ð8Þ

ðf  f c Þ f > f c

introduced in [24] accounts for the effect of accelerated damage development due to void coalescence upon the yield beha
viour. The relation ðf u  f c Þ=ðf f  f c Þ, accounts for the effect of void coalescence on the yield behaviour. Here, f f is the value at

which the failure in a material point occurs. The form (8) implies that f reaches the value f u when the material is totally

damaged, such that f u ¼ 1=q1 . Normally, one chooses a much smaller value for f u , to avoid numerical problems. If f ¼ f f holds, the material point fails. The restriction to f  1 and the recommendation in [3] to use f f  0:2 to 0.25 is consistent with the simulations of [26]. 3. Extension to non-monotonic loading Any constitutive model that involves damage for structures must be able to describe the effects of non-proportional loading on the material behaviour and damage. Progressive localisation of the deformation field in, e.g., homogeneous tensile specimens, or cracked specimens can cause loading path changes. The model of [5], as based on a yield surface which is symmetric with respect to the deviatoric stress planes, predicts negative void volume growth (i.e., healing) under hydrostatic pressure. This is at odds with the observation that, for constant triaxiality, the void volume fraction does not change during closed loading cycles [27]. Thus, the usage of the standard model for cyclic or non-monotonic loading requires an extension. Consider for example the case of localisation during monotonic loading. Early work (e.g., [28]) examined the influence of yield surface curvature using the deformation theory of plasticity and J 2 corner theory, or investigated shear banding using a model for von Mises plasticity with kinematic hardening (e.g., [29]). In the damage context, [6] gave a generalisation of the GTN model




/ðT; X; rF ; f Þ ¼



r2vM ðT  XÞ 3 rH ðT  XÞ

2  ðq1 f Þ ¼ 0 q  1 þ 2 q1 f cosh 2 2 rF r2F

ð9Þ

which is a direct generalisation of the GTN yield function (7) that counts for cyclic loading and kinematic hardening. In particular, the identification

jT  Xj

rF



jTj

ð10Þ

rM

gives the basis to determine the effective yield stress rF . This presumes that the magnitude of the plastic deformation rate be determined solely by that for isotropic hardening. As confirmed by [30], (9) leads to accelerated necking and failure and is not valid for general non-proportional and cyclic loading. A generalisation of the approach in [6] was considered by [31], who based their method on an evolution equation for the plastic spin as presented in [32]. In this regard [7] made a different approach based on the two measures R1 and R2 for isotropic hardening related to deviatoric and hydrostatic stress states, respectively. In the absence of kinematic hardening, these appear in the modified form

/ðT; R1 ; R2 ; f Þ ¼

r2vM ðTÞ R

2 1

 1 þ 2 q1 f cosh


 3 rH ðTÞ 2 f ¼0 2 R2

ð11Þ

of the Gurson yield function. Attempts to use different unit cells and different materials to determine R1;2 numerically proved to be very difficult and time-consuming. Consequently, they resorted to the upper Hashin–Strikman-limit of the yield stress of the porous material in the non-linear case [33–35] and determined R1;2 in approximate analytical form for small values of f . In any case, for realistic numerical simulations, one is still left here with the task of determining R1;2 as a function of the loading history. Using the ansatz of [6,7] generalised this procedure to the case of kinematic hardening. Again resorting to analytical solutions for the deviatoric and hydrostatic parts of X, they formulated the form

/ðT; X; R1 ; R2 ; f Þ ¼

r2vM ðT  rXÞ ðr rM0 þ ð1  rÞR1 Þ




 1 þ 2 f cosh 2

 3 rH ðT  rXÞ 2 f ¼0 2 ðr rM0 þ ð1  rÞR2 Þ

ð12Þ

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

99

for the Gurson yield function, with r as the fraction of combined hardening determined by kinematic mechanisms. Simulations on unit cells helped verify this approach. The application of this model to strongly isotropic hardening materials [27] showed a significant improvement in the prediction of the evolution of the void volume fraction compared to the original Gurson model. As mentioned in [36,27], elastic (or more generally stored-energy) effects may locally affect void growth under cyclic loading, something not considered in this model. In 2003, [8] presented an alternative extension of the GTN model following [9]. They worked with the modified GTN yield function of the form (9) in which the effective yield stress rF is given by rF ¼ ð1  bÞ rY þ b rM rather than by (10). Here, rY and rM are the initial and evolving yield stress of the matrix material, respectively. Further, b is the fraction of the combined hardening due to kinematic processes. This model served for parametric studies in order to determine the effect of kinematic hardening on the onset of macroscopic shear banding. 4. Cyclic loading of unit cells To gain insight into the role of isotropic and kinematic hardening in the process of void growth under cyclic loading conditions, we consider the cyclic loading of unit cells. Assuming a regular, periodic arrangement of voids in the matrix material and axisymmetric unit cells, the unit cell arrangement takes the honeycomb-like form shown in Fig. 1. The corresponding finite-element model of each unit cell is shown in Fig. 2a and b. The cell is subjected to cyclic deformation-controlled loading with a maximum logarithmic strain amplitude of 0.1 as shown in Fig. 2c [13,38]. In order to simulate conditions in the vicinity of a sharp crack tip in a fracture specimen [37], the triaxiality

v¼

rH 2ðs þ 1Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rvM 3 ðs  1Þ2

ð13Þ

in the cell is constrained to be constant during loading and equal to 1, with

s :¼

rradial raxial

ð14Þ

the ratio of the axial to radial stress in the cell. Relation (13) can be inverted to yield

rradial 3v  1 ¼ : raxial 3v þ 2

ð15Þ

To maintain geometric compatibility of the structure, the walls are constrained to remain straight during loading. The unit cells have a radius R of 1, a height h0 of 2, and contain a single void of initial radius 0.114. The matrix material is assumed to be a typical steel characterised by isotropic elastic (k = 121,154 MPa, l = 80,762 MPa), and von Mises inelastic behaviour with linear combined hardening. In particular, we have rM ðeM Þ ¼ 500 þ h eM (MPa) and X_ ¼ c DP (a special case of (24) below), and no damage. The results are to compare with the predictions of the GTN model for the same isotropic hardening behaviour and initial void volume fraction, neglecting void nucleation and coalescence. The distribution of accumulated plastic deformation in the vicinity of the void is shown in Fig. 3. In the isotropic as well as kinematic cases, the maximum value of eM away from the void is determined by the external load and reaches a value of 0.6. Close to the void, however, the reduced effective yield level upon load reversal caused by kinematic hardening results in an additional accumulation of plastic deformation in the matrix compared to solely isotropic hardening. For the purely kinematic case, the maximum of 2.6 lies far upon that for isotropic hardening. Fig. 4 shows the evolution of eM at point P close to the void boundary. In the isotropic hardening case, the slope of eM decreases, in contrast to the linear increase observed for purely kinematic hardening. This reduced amount of accumulated inelastic deformation for purely isotropic hardening results in reduced void growth. Fig. 5 compares results for the change in void volume fraction as a function of load cycle for purely isotropic hardening, purely kinematic hardening and the Gurson GTN model. In particular, the GTN model overestimates the amount of void growth after the first half cycle in comparison to the von Mises results. After load reversal, however, the purely kinematic hardening material near the void yields again before complete release, facilitating a reversal of the void growth during tension and in fact a slight reduction in the void volume below its initial value (Fig. 5). In contrast, no yielding of the purely isotropic material takes place during unloading. Consequently, the inelastic part of the initial void growth is frozen in, resulting in a net increase in f after one loading cycle. This is consistent with the results of [38]. As shown in [6], the interaction of neighbouring voids and shear banding is more pronounced for purely kinematic hardening, accelerates failure and leads to accentuated failure in comparison to purely isotropic hardening. 5. Current model The purpose of this section is to develop an extension of the Gurson model for ductile damage and failure applicable for arbitrary non-monotonic loading histories and arbitrary deformation. This is particularly important towards the application of such a model to predict failure at a sharp crack tip subjected to cyclic loading conditions (e.g., [12,13]). To this end, we

100

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 1. Regular arrangement of unit cells.

Fig. 2. Axis-symmetric finite-element model of a unit cell: (a) dimensions of a single unit cell and stress state, (b) finite-element mesh, and (c) cyclic deformation-controlled loading.

Fig. 3. Distribution of accumulated plastic deformation cycles.

eM in the matrix near the void for purely isotropic or purely kinematic hardening after 3 loading

accumulated plastic strain

M

3.0

isotropic hardening kinematic hardening 2.0

1.0

0.0 0.0

1.0

2.0

3.0

cycle Fig. 4. Evolution of

eM at point P shown in Fig. 3 close to the void for the cases of purely isotropic or purely kinematic hardening in the matrix.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

101

void volume fraction f

0.0018

isotropic hardening kinematic hardening Gurson model (GTN)

0.0016 0.0014 0.0012 0.0010 0.0008 0.0

1.0

2.0

3.0

cycle Fig. 5. Void volume fraction f vs. loading cycle in the unit cell labelled as ‘‘isotropic hardening” and ‘‘kinematic hardening”; compared to the GTN-modelbased RVE of same initial void volume fraction and isotropic hardening neglecting void nucleation and coalescence.

invoke form invariance in order to generalise the model formulation to large deformation by interpreting T as the Cauchy stress. In addition, the flow rule (3) generalises to

DP ¼ c /; M

ð16Þ

in this context as well. Here, DP ¼ symðF_ P F 1 P Þ represents the local inelastic rate of deformation, and F P the local inelastic deformation. Further, M represents the Mandel stress, which is related to the Kirchhoff stress K ¼ detðFÞ T via T T K ¼ F T E MF E  RE MRE ;

ð17Þ

the second form holding in the current context of metal plasticity (i.e., small elastic strain j ln U E j  1) and the polar decomposition F E ¼ V E RE ¼ RE U E of the local elastic deformation F E :¼ FF 1 P . In addition, X is now the back stress in the intermediate configuration. We work with the isotropic Hooke form

M ¼ k trðln U E ÞI þ 2 l ln U E

ð18Þ

of M in what follows. As usual, this is all in the context of the multiplicative elastoplastic decomposition of the deformation gradient [39], or more generally, in the context of the modelling of F P as a change of local reference placement [40]. In this context, the current model represents an extension of the Gurson model in the GTN form for the case of cyclic loading. To this end, we assume that the effect of kinematic hardening on the yield behaviour in the matrix degrades with
increasing f . Further, we follow [7] in introducing different yield stresses rM and rA for deviatoric and hydrostatic yield, respectively. In this case, (7) generalises to




/ðM; X; rM ; rA ; f Þ ¼



r2vM ðM  bðf ÞXÞ 3 rH ðMÞ

2  ðq1 f Þ ¼ 0; q  1 þ 2 q1 f cosh 2 2 rA r2M

ð19Þ

with


bðf Þ :¼ 1  q1 f ðf Þ

ð20Þ




(recall that q1 f ! 1 as f ! f f ). This form is based on the assumption that kinematic hardening processes primarily affect the yield behaviour of the matrix material, as well as on the idea that, in any case, rH ðMÞ  rH ðXÞ. Consequently, we can in fact neglect the dependence of the yield function on b rH ðXÞ without loss of physical generality. Consider next the evolution relation for rA . The form of this relation proposed here is based on the following observations: For strong kinematic hardening and isotropic softening in the matrix material, a dependence of the pressure-based term in / on rM would lead to an unrealistically large accelerated void growth. For materials that show a strong Bauschinger effect, the elastic and plastic ranges are difficult to separate experimentally.1 Further, the transition from the elastic range to the plastic range upon reloading is accompanied by significant hardening. Consequently, we assume that, after each unloading and loading path reversal, the saturation value for rA determines the hardening behaviour [12]. Cyclic loading leads to damage accumulation, as shown experimentally and numerically (e.g., [27]). We return to this issue below. For monotonic and proportional loading, the current model should reduce to the GTN model. 1

This had led some workers to eliminate the transition altogether and work simply with the stress as a function of the strain-rate.

102

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

I

X

m

III

A II

Fig. 6. The equivalent stress

A possible evolution relation for

r_ A ¼



rA determined by (21) as the limiting value of rV from (22).

rA meeting these requirements will be

r_ V if rV ¼ rA and r_ V > 0 0 if rV < rA or r_ V 6 0

ð21Þ

and in terms of the equivalent stress measure

rV ðrM ; XÞ :¼ rM þ rvM ðXÞ:

ð22Þ

Indeed, note that (21) represents a monotonically increasing function of the history-dependent maximum value of rV . Fig. 6 shows rA in principal stress space as the limiting stress for rV . Next we turn to the modelling of hardening in the matrix material. In the case of the back stress X, we assume that its evolution is given by the Armstrong–Frederik ansatz taking growth and saturation processes into account. Further, we generalise the approach of [41] to large deformation and consider a number of mk kinematic hardening mechanisms, so that

X¼

mk X Xa

ð23Þ

a¼1

holds. Here,

X_ a ¼ ca devðDP Þ  ba z_ X a

ð24Þ

with z_ ¼ jDP j the plastic arc length. As usual, ca governs the growth of kinematic hardening at the onset of plastic deformation, and ba determines the saturation behaviour. Finally, the yield stress rM in the matrix material may represent the following sum

rM ðeM Þ ¼ rM0 þ

Xmi

ka eM

r ð1  e a¼1 a

Þ

ð25Þ

of mi Voce-like mechanisms (e.g., [42]). Here, ka represents the saturation constant, and r a the magnitude of each mechanism. In particular, if ra is positive (negative), we have a hardening (softening) mechanism. With the current model now in hand, we apply it to the unit cell case from the previous section. Fig. 7 displays the thereby predicted evolution of the void volume fraction for 10 loading cylces. In contrast to the von Mises case in Fig. 5, the results here imply that purely kinematic hardening leads to much stronger void growth than purely isotropic hardening for the reasons discussed above. Indeed, kinematic hardening results in an enhanced triaxiality. This result is in agreement with that of [10], who show for unit cell calculations under both constant stress triaxiality and strain triaxiality that kinematic hardening leads to a faster void growth and failure. The systematic investigation of different triaxialities (v = 1.0, 2.0, 3.0) shows that under cyclic loading damage accumulation is a monotonically increasing function of v, as confirmed by [38]. Fig. 8 shows the void volume fraction f in a unit cell for three different triaxialities at the end of each loading cycle, corresponding to the curve min(f) in Fig. 7. The current extension of the Gurson model is qualitatively and quantitatively able to predict void growth as for nonmonotonic loading paths.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

103

Fig. 7. Development of void volume fraction in the unit cell from the previous section during 10 cycles predicted by the current extended Gurson model for the cases of purely isotropic and purely kinematic hardening.

6. Algorithmic formulation of the model Except for the inelastic anisotropy due to inclusion of the back stress, the model developed in the previous section is isotropic. As discussed in detail in [43–45], exponential backward-Euler integration (e.g., [46,47]) of the evolution relation for F P over a time interval ½t n ; t nþ1 , gives that the elastoplastic multiplicative decomposition of F, and the assumption of small elastic strains j ln U E j  1, result in the algorithmic split of the flow rule into

ln U E tr ¼ ln U E nþ1 þ tnþ1;n DP nþ1

ð26Þ

RE nþ1 ¼ Rnþ1 RTP nþ1 ¼ RE tr expðt nþ1;n W P nþ1 Þ;

ð27Þ

and

both to lowest order in j ln U E j, with t nþ1;n :¼ t nþ1  t n and

1 lnðC E tr Þ; 2 T ¼ F E tr F E tr ;

ln U E tr :¼ C E tr

ð28Þ

F E tr :¼ F nþ1 F 1 P n; represents the so-called trial value of ln U E nþ1 (with RE tr as the polar rotation of F E tr ). Although theoretically either RP nor the plastic spin W P could conceivably influence the evolution of X, we assume for simplicity here that RP ¼ I, in which case RE is modelled by R via (27)1 and we have the Green-Naghdi back-rotated configuration for X. In addition,

F E nþ1 ¼ RE nþ1 U E nþ1 ¼ Rnþ1 expðln U E nþ1 Þ

ð29Þ

then holds for the algorithmic local elastic deformation. Backward-Euler integration of the evolution relations (4), (6) and (24) over the time interval ½tn ; tnþ1 together with (26) and the yield condition / ¼ 0 as based on (19) yields the system of equations

r ln U E nþ1;n ¼ ln U E nþ1  ln U E tr þ cnþ1;n N nþ1;n ¼ 0;   2 devðNÞ  b1 jNj Y 1 ¼ 0; r Y 1 nþ1;n ¼ Y 1 nþ1  Y 1 n  cnþ1;n 3 nþ1;n .. .



 2 devðNÞ  bmk jNj Y mk ¼ 0; 3 nþ1;n n o 1 r eM nþ1;n ¼ eM nþ1  eM n þ cnþ1;n ð1  f Þ /; rM ¼ 0; nþ1;n n o 1 ¼ 0; r f nþ1;n ¼ f nþ1  f n  cnþ1;n ð1  f Þ trðNÞ þ ð1  f Þ A /; rM r Y mk nþ1;n ¼ Y mk nþ1  Y mk n  cnþ1;n

ð30Þ

nþ1;n

r c nþ1;n ¼ /nþ1;n ¼ 0; for the unknowns ðln U E nþ1 ; Y 1 nþ1 ; . . . ; Y mk nþ1 ; eM nþ1 ; f nþ1 ; cnþ1;n Þ, with cnþ1;n :¼ cnþ1 ðt nþ1  t n Þ the incremental plastic multiplier, znþ1;n ¼ cnþ1;n jN nþ1 j,

X a :¼ ca Y a in terms of the back ‘‘strain” Y a , and

ð31Þ

104

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 8. Void volume fraction in the unit cell during cycling predicted by the current extended Gurson model for purely kinematic hardening at three values of the triaxiality v.

N :¼ /; MX

ð32Þ

the normal to the yield surface. Because the backward-Euler form



rA nþ1 ¼ rA n þ

rV nþ1  rV n if rV nþ1 > rV n 0 if rV nþ1 6 rV n

ð33Þ

for the incremental update of rA from (21) via (22) is not continuous, rA is held constant at rA n during the solution of (30)
and then updated to rA nþ1 via rV nþ1 and (33) at the end. In addition, the approximation of the function (8) for f ðf Þ by a parabola avoids numerical problems and improves numerical robustness. On the basis of the relations (18), (22), (23), (25) and (31) for the stresses, the system (30) can be written in the compact form

r nþ1;n ¼ r n ðnþ1 ; F nþ1 Þ ¼ 0

ð34Þ

in terms of the set

nþ1 :¼ ðln U E nþ1 ; Y 1 nþ1 ; . . . ; Y mk nþ1 ; eM nþ1 ; f nþ1 ; cnþ1;n Þ

ð35Þ

of independent internal variables. This system can be solved via the usual predictor–corrector procedure (e.g., [48]) as based on the Jacobian

Jnþ1;n ¼ Jn ðnþ1 ; F nþ1 Þ :¼ r n; nþ1 ðnþ1 ; F nþ1 Þ:

ð36Þ

Adaptive time step methods helped to deal with convergence issues surrounding the treatment of the cosh-term in the yield function (19) and the nucleation term in the evolution relation (21). They were based on the minimum and maximum number of Newton iterations. Along with the solution (35), one obtains in this fashion the stress update

K nþ1 ¼ RE nþ1 M nþ1 RTE nþ1 ¼ Rnþ1 M nþ1 RTnþ1

ð37Þ

for the Kirchhoff stress via (17) and (18). Algorithmic linearisation of these with respect to F nþ1 yields

dK nþ1 ¼ ð@ aF nþ1 K nþ1 Þ½dF nþ1

¼ ðI  K nþ1 þ K nþ1 M IÞ½ð@ aF nþ1 Rnþ1 Þ½dF nþ1 RTnþ1 þ ðRnþ1  RTnþ1 Þðk I  I þ 2 l I  IÞð@ aF nþ1 ln U E nþ1 Þ½dF nþ1

ð38Þ

via (18) in terms of the tensor products ðA  CÞ½B :¼ ABC and ðA M CÞ½B :¼ ABT C on second-order Euclidean tensors. Here, the algorithmic derivative

ð@ aF nþ1 Rnþ1 Þ½dF nþ1 RTnþ1 ¼ fðdFÞRT  2 RðU ; C Þ½symðF T dFÞ F 1 gnþ1

ð39Þ

of Rnþ1 with respect to F nþ1 follows from the polar decomposition of F, and U ; C can be obtained, e.g., from the functional calculus (e.g., [49]) or via Padé approximation. Analogously, T 1 ð@ aF nþ1 ln U E nþ1 Þ½dF nþ1 ¼ ½J1 nþ1;n ln U ðln U E tr; F nþ1 Þ½dF nþ1 ¼ ½Jnþ1;n ln U ðln; C E tr Þ½symðF E tr dF E tr Þ

ð40Þ

1 1 follows from (28)1,2 and (30)1, with ½J1 nþ1;n ln U the upper left submatrix of Jnþ1;n , and dF E tr ¼ ðdF nþ1 ÞF P n . As with U ; C ; ln; C E tr can be obtained using the functional calculus.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

105

To implement the above algorithm in commercial programs such as ABAQUS via a user material interface such as UMAT, we must adapt the corresponding consistent tangent formulation as based on (38) to their input structure. To this end, consider the additive split

dF nþ1 ¼ Ddnþ1 F nþ1 þ W dnþ1 F nþ1

ð41Þ

of dF nþ1 , where Ddnþ1 is symmetric, and W dnþ1 skew-symmetric. On the basis of this split, (38) yields the form

@ dK nþ1 @Ddnþ1

¼ ðI  K nþ1 þ K nþ1 MIÞfI  V  2 ðR  F 1 ÞðU ; C ÞðF T  FÞgnþ1 þ ðRnþ1  RTnþ1 Þðk I  I T þ 2 l I  IÞ½J1 nþ1;n ln U ðln; C E tr ÞðF E tr  F E tr Þ:

ð42Þ

for the algorithmic derivative of K nþ1 with respect to Ddnþ1 . This last result divided by detðF nþ1 Þ represents the input variable DDSDDE in ABAQUS UMAT for the current model. As it turns out, the approximation

@ dK nþ1 @Ddnþ1

T  ðRnþ1  RTnþ1 Þðk I  I þ 2 l I  IÞ½J1 nþ1;n ln U ðln; C E tr ÞðF E tr  F E tr Þ

ð43Þ

to this tangent which neglects the contribution from the algorithmic variation of Rnþ1 with respect to F nþ1 as being small is also quite good [12]. 7. Material properties and specimen preparation The material considered in this investigation is the high-strength industrial steel, 10MnMoNi5-5. Table 1 lists its chemical composition. Table 2 gives an overview of the main properties, i.e., ductility and strength for 25 °C and 300 °C. The yield limit ReL in 10MnMoNi5-5 is well-defined up to about 100 °C. Above this temperature, one usually works with the 0.2%-strainbased yield stress Rp0:2 . The Lüders strain in 10MnMoNi5-5 at room temperature is about 1.5%. The used samples were part of a specimen cut out from the welded segments of a pressure vessel (wall thickness 340 mm, outside diameter 1840 mm, length 6 m). Klingbeil et al. [16] carried out a detailed series of material investigations for one of the two ring segments which were separated from this specimen. The manufacturing process applied the built-up material continuously in a helical form with the width of the welding bead overlapping half of the welding bead of 28 mm. After completion, the vessel was annealed for stress-relief for 20 h at 520 °C. The heat treatment during the application of subsequent welding layers (Fig. 9) substantially influenced the development of the ferrite-perlite microstructure. Thereby a formation of fringe crystals occurs with the directional hardening above the melting point curve. Directly under the melting point curve a coarse-grained structure develops, whilst a fine-grained structure develops in the greater distance. Instead of the anisotropy expected from the manufacturing process, the material shows a very isotropic behaviour regarding all material properties [16]. 8. Monotonic and cyclic test results A set of tension tests enabled the derivation of the yield curves shown in Fig. 10 that characterise the hardening behaviour of 10MnMoNi5-5. The Bridgman correction factor allowed the estimation of the true (Cauchy) stress vs. logarithmic strain after necking [50]. Details of the experimental set-up can be found in [16]. The available experimental data made it possible to separate the parameter identification into small and large deformation regimes. At small deformation, it is sufficient for the description of hardening to use two isotropic hardening mechanisms, and one kinematic hardening mechanism, to model the data. Including data at large deformation as well, however, requires an additional isotropic and kinematic hardening mechanism, such that mi = 3 and mk = 2. Due to the experimental restrictions, this large-deformation data is from monotonic tests. Besides tension tests, compression tests on short cylindrical specimens also provided data to determine yield curves at RT (Fig. 11). When properly corrected, data from compression tests at large deformation are more reliable. Hence, the experimental yield curves utilised in what follows base upon experimental data from both test types. The yield curve from compression tests demands a correction (shift) to the level of the curve from tension tests until uniform elongation is reached. Fig. 12 illustrates the resulting composite yield curves utilised for parameter identification. Application of data from tensile tests on notched round bars with notch radii of R = 4 mm (R4) and R = 10 mm (R10) verified the composite yield behaviour. Fig. 13 displays the experimental and simulation results for the reaction force F as a function of the change in diameter DD at the notch. The damage evolution, which initiates locally at the centre of the specimen in the notch, has no significant influence to the overall response of the specimen before reaching the load maximum. The accumulated plastic strain eM attains magnitudes of 50% (R10) and 90% (R4) at the notch. As such, it gets at magnitudes which are representative of the whole range of the hardening behaviour considered here. Note also that the triaxiality v reaches magnitudes of 0.8 (R10) and 1.2 (R4), which are much higher than in smooth round bars.

106

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Table 1 Chemical composition of 10MnMoNi5-5 taken from [16]. Element Mass (%)

C 0.099

Si 0.142

Mn 1.241

P 0.007

S 0.005

Cr 0.050

Ni 1.055

Element Mass (%)

Mo 0.587

V 0.002

Al 0.009

Cu 0.038

W 0.011

Co 0.029

Nb 0.002

Element Mass (%)

Sn 0.0014

Ti 0.005

As 0.005

Pb 0.003

Sb 0.0031

Table 2 Material properties of 10MnMoNi5-5 at two different temperatures: Young’s modulus E, Poisson’s ratio m, yield limit ReL (well-defined below 100 °C), yield stress Rp0:2 at 0.2%-strain (used above 100 °C) [16]. h (°C)

E (MPa)

m (–)

ReL (MPa)

Rp0:2 (MPa)

25 300

210,000 195,000

0.3 0.3

620 –

– 548

Fig. 9. Micrograph of the melting line in the cross section of a welding (V = 20:1) [16].

Fig. 10. Uniaxial tension test results for 10MnMoNi5-5 at 4 temperatures from (i) monotonic tensile tests (B8 40 according to DIN 50125) up to 100% strain, and (ii) strain-controlled cyclic tests with strain amplitudes between 1% and 3%.

For the determination of the yield curve at 300 °C, no compression test data was available. Since the simulation of notched specimens overestimates the force when based on tension test data, this data was corrected using the same factor as at RT (Fig. 14). Again, experimental data from notched specimens helped to verify the identified model parameters. Fig. 15 shows the averaged experimental data (R4: 3 experiments, R10: 3 experiments) and the simulation results obtained from the corrected yield curve and the yield curve from the tensile test. Next, we turn to the cyclic test results. There were data from strain controlled cyclic experiments on round bars at 25 °C, 100 °C, 200 °C and 300 °C [16]. Since higher strain levels lead to buckling, the strain amplitudes at higher temperatures had to be confined to 2%. Fig. 16 shows half cycles as a function of the accumulated strain at different temperatures.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 11. Possible measurement errors during tension (a) and compression (b) tests involve radius of curvature R, damage f and friction coefficient

107

l.

Fig. 12. Experimental yield curves from tension and compressions tests compared with the estimated combined curve including the Bridgman correction at RT. Ag represents the strain level in compression tests up to which no significant bulging of the specimen occurs.

Fig. 13. Comparison of simulation results as based on the estimated yield curve in Fig. 12 with experimental results for the reaction force F as a function of change in diameter DD on notched round bars with notch radii of 4 mm (R4) and 10 mm (R10).

108

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 14. As in Fig. 12 but now at 300 °C.

The reduction in yield stress and hardening level with increasing temperature is clearly visible. In addition, kinematic hardening effects are apparent at all temperatures. The cyclic test results at RT in Fig. 17 exhibit strong saturation towards a stable equilibrium hysteresis after only a small number of load cycles. At 300 °C this is the case not until 10 cycles, as shown in Fig. 18. 9. Hardening parameter identification To estimate starting values for the hardening parameter identification, consider the half-cycles of the cyclic tests discussed in the previous section. As these are all uniaxial and proportional, assume that this loading takes place in the 1-direction. Further, assume that T 11 and X 11 are the only components of interest of the Cauchy stress and back stress. DP is diagonal, and X_ 11 ¼ 2 c DP11  b jDP j X 11 holds for the evolution relation for the back stress from (24). This equation implies 3

that the initial increase of x :¼ 32 X 11 at yielding is given by c. Furthermore, the yield curve saturates to a stress value of rM0 þ c=b. These facts can be used to estimate c1 and b1 from the experimental data as shown in Fig. 19. At RT, for example, saturation occurs in 10MnMoNi5-5 after only two cyles. Note also that the slope c1 at the onset of plastic deformation is approximately equal the slope c01 after the first load reversal. The hardening behaviour at higher strain levels as illustrated in Fig. 17 is not asymptotically approaching a constant stress level, but approximates a straight line of a certain slope. An approximation of the kinematic hardening parameter c2 describes this characteristics and therefore the slope of the line. The assumption of constant isotropic hardening for large deformation allows the approximation of the ratio c2 =b2 with the aid of the yield curve, using b2 as a measure for the deviation from a straight line (Fig. 20): Table 3 summarises the parameter values estimated in this way for kinematic hardening. The same strategy is applied for 300 C. Due to the similarity of the forms of the curves in Fig. 16, the same estimates for the kinematic hardening parameters apply at the higher temperature.

Fig. 15. As in Fig. 13 but now at 300 °C.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

109

Fig. 16. Comparison of cyclic load test results for a strain amplitude of 1% at different temperatures. The test at 200 °C was carried out at a higher strain amplitude; for comparison with the other cases, only the loading phase is shown here.

Fig. 19 shows an equivalent estimate in case of isotropic hardening. In particular, microscopic shear-banding at the onset of yielding results in an initial isotropic softening as reflected in the value of r1 . This is followed by a second isotropic hardening mechanism r2 and k2 . Together, these lead at RT to a relatively quick dynamic recovery and saturation of the yield stress towards a constant value. In contrast, the increasing stress amplitude over a couple of load cycles at 300 °C in Fig. 18 indicates a continuing isotropic hardening for moderate deformation at this temperature in 10MnMoNi5-5. Consequently, a third isotropic hardening mechanism r3 and k3 is required to fit the data. This increases the number of isotropic hardening parameters at 300 °C to a total of six, in comparison to a total of four at RT (Table 4). Lastly, Table 5 lists the estimates for E and the initial yield stress used for the parameter identification. In particular, the value of rM0 at 300 °C differs from the value of Rp0:2 at this temperature given in Table 2. Indeed, the former value is estimated from the experimental results in Fig. 18. The parameter identification of the combined hardening model needs the solution of an optimisation problem and finding of a suitable procedure amongst possible strategies. The advantage of using non-gradient methods such as the evolution strategy [14] is easy implementation and the robustness of the solution. On the other hand, there is a higher computational cost for differentiable problems in comparison to gradient-type approaches. In the current context, this was not a major disadvantage since the homogeneous deformation state of the specimens allowed the parameter identification to be carried out for the most part without extensive finite-element simulation. Table 6 summarises the results for the parameter identification at each temperature as based on this approach. Comparison of the experimental results with those from the identified model in Figs. 21 and 22 demonstrates the robustness of the identified values for all kind of experiments. In particular, these results show that the back-stress X 1 is dominant at smaller deformation, whilst in the monotonic tension tests at larger deformation the contribution of X 2 becomes significant. This correlates with large differences in the magnitude of the corresponding parameter values as shown in Table 6. In

Fig. 17. Cyclic loading results at RT for strain amplitudes of 1%, 2% and 3%.

110

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 18. Cyclic load tests at 300 °C for strain amplitudes of 1%.

addition, the fit results account for the initial softening behaviour in r (Fig. 21). The subsequent moderate increase is followed by saturation at RT, whereas at 300 °C, hardening is present throughout the entire loading range. 10. Damage parameter determination Finite-element simulations of experiments on notched round bars helped to determine the damage parameters of the current model [12,13]. Comparison with experimental results of other fracture mechanics specimens from [16] verified the identification of the corresponding model. In particular, specific values of the triaxiality v characterise different fracture specimens as schematically illustrated in Fig. 23. In contrast to the hardening case, identification of the isotropic damage parameters requires finite-element simulation of all specimen geometries at the underlying constraint conditions. Due to the difficulty of investigating ductile failure in smooth specimens in a predictable, controlled manner, notched specimens better account for the identification of the damage parameters. Tensile tests on notched round bars with different notch radii exhibit a trixiality of about 1, with the sharpness of the notch determining the level of triaxiality. Higher triaxiality values of 2–3.5 develop at the crack tip of compact tensile specimens and centre cracked tensile specimens. Results from metallographic investigations provide estimates of the initial void volume fraction f 0 and the volume fraction f n of void-nucleating particles. Micrographs of different stages of the fracture process together with the accumulated plastic strain from finite-element simulations allow to assess the parameters en and sn . Since this has already been done for 10MnMoNi5-5 [16], we simply take advantage of these results. Electron-microscopy indicated a very small initial void volume fraction for this material, but similar to values identified for other high-strength heat-treated steels [16]. This complicates a reliable parameter adaptation due to the variation of the force maximum of tested notched round bar specimens. In this case it is assumed for the initial void volume fraction to be an order of magnitude below the volume fraction of void volume forming particles. Therefore, f 0 ¼ f n =10 holds. The smaller area fraction of non-metallic inclusions and their smaller diameter (oxides < 10 mm) [16] indicate smaller values for en

Fig. 19. Estimation of starting values for the hardening parameters from cyclic test data.

111

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 20. Estimation of kinematic hardening parameter c2 and b2 at RT.

and sn compared to other steels, since for smaller inclusions plastic deformation can lead to decohesion at an earlier time. Electron-microscopy restricted the order of magnitude of f n to f n < 0:005. As mentioned above, experimental as well as numerical results on notched bars were the basis to identify the critical void volume-fraction f c for coalescence as well as the void volume-fraction f f at ductile failure. We neglected any possible temperature-dependence of the damage properties. With a sufficiently high initial void volume fraction f 0 , the onset of ductile damage causes a reduction of the load maximum, which allows an additional estimate of f 0 . Initially, damage develops due to decohesion and void growth. The parameters f n , en and sn of the formerly introduced void nucleation model describe this behaviour. Moderately notched round bars subjected to displacement-controlled tensile loading show crack initiation at first in the centre of the specimen followed by ductile crack growth that causes a drop of the external load. This point of decreasing load as schematically displayed in Fig. 24, helps identify f c and f f . To determine the parameters, the data on notched round bars at RT and 300° C presented earlier is used. The results of this identification are presented in Figs. 25 and 26 and Table 7. Simulation results from the hardening parameter identification show the insensitivity of the load maximum. The point of ductile failure in the specimens is predicted very well with the identified parameter sets. At RT, the experimental results show almost no scatter. Therefore, the onset of the decrease in load is clearly identified (Fig. 25). At 300 °C, there was more scatter (Fig. 26). Consequently, the point of rapid load-drop in the identification was assumed to be given by the average of the experimental results.

11. Fracture specimens Experiments on ductile crack growth in side-grooved C(T) sg specimens (Fig. 27) [16] as well as numerical simulations were predestined to validate the currently identified model. Since side-grooved C(T) specimen geometry results in the suppression of transverse contraction at the crack tip [16], plane-strain conditions are justified. The specimens of W = 50 mm, a0 = 25 mm had a thickness of B = 25 mm, side grooves with a depth of 2.5 mm each, and a notch radius of 0.25 mm. The simulation results at both temperatures show a small deviation from the experimentally observed F—V LL curves taken from [16] (Figs. 28 and 29). The ramp-like form of the simulated curves arises because of element elimination during the simulation of ductile crack growth. Cracking starts when the first element in the ligament is fully damaged. Since the values of force and displacement at crack initiation show a strong scatter in the experiments, Table 8 summarises values that demonstrate the good agreement between the simulation and the experiments. The model correctly predicted crack-resistance curves from experiments up to a crack growth Da of 1 mm (Fig. 30). This good agreement with the experiments is an additional validation of the identified model. At 300 °C, crack growth is observed a small distance ahead of the crack tip, which agrees with the experiments. The simulated crack–resistance curve (Fig. 31) includes the accumulated length of all damaged elements. Again, the agreement between simulation and experiment is quite good. This completes model identification and validation. Assuming that its identification is robust, the model can be applied to completely different specimen geometries and structures consisting of the same material.

Table 3 Estimated values of the kinematic hardening parameters at two temperatures. h (°C)

c1 (MPa)

b1 (–)

c2 (MPa)

b2 (–)

25 300

100,000 100,000

300 300

1000 1000

3 3

112

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Table 4 Estimated values of isotropic hardening parameters at two temperatures. h (°C)

r 1 (MPa)

k1 (–)

r 2 (MPa)

k2 (–)

r 3 (MPa)

k3 (–)

25 300

400 200

200 200

200 100

50 50

– 100

– 1

Table 5 Youngs modulus E and initial yield stress

rM0 at two temperatures.

h (°C)

E (MPa)

rM0 (MPa)

25 300

210,000 195,000

620 500

Table 6 Identified hardening parameter values for 10MnMoNi5-5 at two temperatures. h (°C)

r 1 (MPa)

k1 (–)

r 2 (MPa)

k2 (–)

r 3 (MPa)

k3 (–)

25 300

406.0 140.8

213 3.45

167.0 96.7

62.1 46.7

– 112.1

– 1.11

h (°C)

c1 (–)

b1 (MPa)

c2 (–)

b2 (MPa)

25 300

101,000 63,500

287.0 244.8

1060 1170

2.6 3.12

Fig. 21. Comparison of cyclic experimental and parameter identification results at RT.

Fig. 22. Comparison of monotonic experimental and parameter identification results at RT.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 23. Schematic draw of the influence of triaxiality

113

v on crack initiation.

Fig. 24. Influence of the damage parameters on the force–displacement behaviour of round notched bars.

12. Experimental model validation Blumenauer et al. [15] carried out X-ray measurements of crystal lattice deformation around the crack tip in fracture mechanics specimens of 10MnMoNi5-5 subjected to cyclic loading. The measurement procedure used the W-differential method. Together with the elasticity constants for the crystal lattice [51], these measurements enabled the determination of residual compressive stress fields at the crack tip for different preloading-levels, and a detailed spatial resolution of stress gradients at the crack tip. The measurement took place at the free boundary in a thin layer up to 5 lm: [15] next to the surface. Initial simulations of the stress field at the crack tip in a C(T) specimen after repeated preloading and plane strain conditions (e.g., [12,13]) showed only qualitative agreement with these experimental results. The reason for the discrepancy was the poor approximation of the stress state at the free surface of the specimen using the plane strain assumption. Consequently, three-dimensional simulations had to resolve the problem. The finite-element model of the C(T) specimen under consideration (Fig. 32a) included elements of decreasing thickness towards the free surface (i.e., down to 0.05 mm). In the control area (Fig. 32a) elements with a smallest element edge-length of 0.1 mm discretised the ligament region located to the right of the (grey coloured) initial crack. Taking advantage of the specimen’s symmetry in the y- and z-directions, only one fourth of the specimen was modelled. Fig. 32b shows the history of experimental loading and measurements. The experimenters removed one or more layers of surface material twice from the specimen by electrolysis (wear out), thus, obtained stress results in the width (i.e., z) direction. As shown, the first load was 55 kN at room temperature. After unloading and removal of a layer of 0.5 mm thickness in the control area (20 mm 20 mm), the first measurement followed. The results showed only low residual stress characteristics. Then a second loading stage up to 78 kN into the range of crack initiation for this specimen took place. After unloading a second measurement took place without a removal of material. In order to check any effects of the removal procedure on the results, measurement 3 (Fig. 32b) occurred after the wear out of a second layer of material on the surface of the specimen. The removal procedure was also simulated by finite-element analysis through removing the corresponding elements and solving for the new equilibrium state. The resulting asymmetry of the specimen geometry was neglected. Fig. 33 compares measured residual compressive stresses in the crack opening direction T 22  T yy after the first loading cycle up to 55 kN with the finite-element solution. The measurements confirm the symmetry of the stress field [15]. But, the predicted stress is larger than the experimental, whereas the spatial extent of the residual stress field is smaller (Fig. 33). Fig. 34 shows the analogous comparison between the experimental and simulation results after the second loading stage. After removal of an additional layer, Fig. 35 gives the stress distributions of the residual stresses at measurement 3 (Fig. 32b). After the second loading stage, the measured stress maxima appear approximately 1 mm ahead of the crack tip, whereas the

114

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 25. Comparison of simulation by identified model with experimental results for notched round bars at RT.

Fig. 26. Comparison of simulation by identified model with experimental results for notched round bars at 300 °C.

predicted maximum is at the crack tip. It would appear that the much higher level of plastic deformation in 10MnMoNi5-5 at the crack tip results in strong blunting and a shift of the stress maximum away from the crack tip, something which texture and flow-line formation confirm. Apparently, the current finite-element simulation does not resolve this effect. Otherwise, however, the results of measurement and simulation are in good agreement for all loading stages. In each case, the highest residual stresses occur in the ligament adjacent to the crack tip and their spatial distribution takes the typical butterfly form relative to the crack tip. Comparing the results as a function of loading stage (e.g., Figs. 33 and 34), one sees that additional loading leads to ‘‘expansion” of the residual stress field further into the ligament. Comparing these results with those in Fig. 35, one sees that the layer removal process does not significantly affect the residual stress state. As Fig. 36 shows, the maximum values of T 22 are significantly higher inside the specimen than on the surface. This is due to the high hydrostatic stress in the interior, resulting in residual stress values up to approximately 1200 MPa. This stress reduces substantially towards the free surface, resulting in strong stress gradients in the material within a distance of approximately 0.1 mm. 13. Crack opening stresses of C(T) specimens The kind of hardening of the material in the ligament adjacent to the crack tip substantially influences the stress field around the crack tip during non-proportional loading. Indeed, since the equivalent plastic strain at the crack tip can achieve up to 20% during each load cycle and easily accumulate up to 100% after 5 cycles, substantial hardening can occur locally. To look into this in detail, consider the C(T)25 specimen without notches of Fig. 37. For such specimens, the ratio of the initial

115

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123 Table 7 Estimated and identified damage parameter values for 10MnMoNi5-5 at two temperatures. h (°C)

f 0 (–)

sn (–)

en (–)

f c (–)

f f (–)

25 300

0.0005 0.0005

0.1 0.1

0.05 0.05

0.025 0.025

0.2 0.2

Fig. 27. Compact tensile specimen C(T) sg with side grooves.

Fig. 28. Comparison of experimental [16] and simulated load–displacement curves for C(T) sg25 specimens at RT.

crack length to the specimen width, a0 =W, and the shape of the initial crack front, have significant effect on the force F measured during the experiment. Therefore, a finite-element simulation should reproduce the initial crack geometry measured by [16]. To take account of the asymmetry of the crack front, as illustrated in Fig. 37, symmetry only along the x–z-plane could be assumed. Initially, [16] carried out three-dimensional simulations of the propagation of such a crack. They mapped the crack front to a regular grid of nodes in the ligament by generating a set of two-dimensional meshes parallel to the x–y-plane along the z-axis. This technique affects the predictions for the development of stress and damage due to the constraint along the irregular mesh front. However, the present study modelled the curved crack front exactly resulting in a smoother stress distribution at the crack tip up to a simulated ductile crack propagation of Da = 2 mm. The following simulations subjected the specimen to a series of loading–unloading cycles under displacement controlled loading up to a maximum value of V LL = 0.5 mm. This was followed by force-controlled unloading down to 0 kN. All together, a total of 8 loading–unloading cycles were applied. The maximum displacement of V LL = 0.5 mm resulted in the maximum load level of 60 kN as observed in the experiment. On the other hand, this load produced a slightly lower load-line displacement in the simulations (V LL  0.48 mm) compared to the experiments. Fig. 38 displays experimental and predicted results of the force versus load-line displacement curves F—V LL . The experimental results show a clear hysteresis due to plastic deformation around the crack tip. The residual displacement after unloading increases with each load cycle, resulting in increasing V LL values in the unloaded state. Basically, the agreement of the experimental and simulation results is good. In the current cyclic loading context, the increase of V LL in the unloaded state (Fig. 39) is related to effects of (kinematic) hardening on the plastic flow behaviour near the crack tip, which result in

116

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 29. Comparison of experimental [16] and simulated load–displacement curves for C(T) sg25 specimens at 300 °C.

Table 8 Comparison of experimental and numerical results for crack initiation of a C(T) specimen. 25 °C

300 °C

F (kN)

V LL (mm)

F (kN)

V LL (mm)

Experiments 48.04 47.96 48.88 50.24 50.91

1.00 1.00 1.11 1.41 1.45

46.10 48.25 47.62 49.37 51.31

0.70 0.71 0.75 0.80 0.90

Average 49.20

1.19

48.53

0.77

FE simulation 51.14

1.06

50.42

0.83

Fig. 30. Comparison of experimental [16] and simulation results for crack resistance J as a function of crack growth Da for C(T) sg25 at RT. ‘averaged’: bestfit curve of test data.

increased accumulated plastic deformation in comparison to purely isotropic hardening. The overall discrepancy between predicted and measured values of V LL after unloading could be due to the finite-element mesh being too coarse, resulting in a somewhat stiffer response. Consider next the simulated spatial distribution of the crack-opening stress as for the loading cycles of Fig. 38. The 3D diagrams of each of Figs. 40–42 represent the stress for the 1st, 2nd and 8th loading cycle at maximum load, respectively [12,13,16]. The ligament has an area of 4 mm 12.5 mm, z = 0 is the centre of the specimen and the x-axis indicates the direction of crack propagation. The step-like curved crack front in the x–z-plane may be recognised clearly. For all three load cases, the stress maximum lies directly adjacent to the crack front. As discussed in the foregoing section, there is a stress

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

117

Fig. 31. Comparison of experimental [16] and simulated results for crack resistance J as a function of crack growth Da for C(T) sg25 at 300 °C.

Fig. 32. (a) C(T)25 specimen (a0 =W = 0.5) with control area in which the experimental measurements were carried out. (b) Procedure used to obtain experimental crack-opening and residual stress results at the crack tip [52].

Fig. 33. Experimental (left) and predicted (right) crack-opening stress T 22  T yy distribution at the crack tip after the initial loading cycle to 55 kN and the first surface layer removal via electrolysis (measurement 1: Fig. 32b).

gradient towards the free surface of the specimen along the z-axis which can be adequately resolved only by using a sufficiently fine mesh. The results clearly exhibit the distribution of the peak stress along the z-axis with the maximum located in the centre of the specimen and the minimum at the specimen boundary. During the 1st loading stage (Fig. 40), the ligament stress T 22 decreases monotonically along x after a certain distance away from the crack tip. Repeated loading leads to increased blunting and a reduction in the stress maximum (Figs. 41 and 42). The distribution of the crack-opening stresses during the 2nd and 8th loading has several similarities, whereas the stress distribution caused by the 1st loading is more discrete. The comparison of the crack-opening stresses during the 1st and during the 2nd loading stage clearly exhibits the effect of

118

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 34. Same as in Fig. 33 after the second load cycle to 78 kN (measurement 2: Fig. 32b).

Fig. 35. Same as in Fig. 34 after removal of an additional surface material layer via electrolysis (measurement 3: Fig. 32b).

Fig. 36. Simulated crack-opening stress T 22 in the ligament at the boundary of a C(T) specimen (measurement 3: Fig. 32b).

preloading. Increasing deviation from this behaviour occurs upon repeated loading of the material, which is then no longer free of residual stresses. Fig. 43 outlines the increase of equivalent plastic strain during cyclic loading. The graphic shows the growth of the 1%-strain-contour in the ligament for the 1st, 2nd and 8th loading cycle, respectively. Being unconstraint, the largest

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

119

Fig. 37. Initial crack geometry of the C(T)25 specimen for the finite-element simulation.

Fig. 38. Experimental and simulation results for the force F as a function of the corresponding load-line displacement V LL in a C(T)25 specimen during cyclic loading.

Fig. 39. Experimental and simulated values of V LL each after unloading.

accumulated plastic deformation occurs near the free surface of the specimen. On the other hand, the constraint reduces deformation in the interior of the specimen.

120

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 40. Simulated crack-opening stress T 22 in the ligament after the 1st loading stage.

Fig. 41. Simulated crack-opening stress T 22 in the ligament during the 2nd loading stage.

Fig. 42. Simulated crack-opening stress T 22 in the ligament during the 8th loading stage.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

121

Fig. 43. Development of plastic deformation at the crack tip with increasing load cycles.

Fig. 44. Evolution of the void volume fraction f at selected points along the crack front.

Fig. 45. Simulated residual stress T 22 in the ligament after the 1st unloading stage.

As shown in Fig. 44, beyond the development of plastic deformation, the asymmetry of the crack influences the state of triaxiality in the specimen and therefore the evolution of damage significantly. Here, the figure shows damage by void nucleation, growth and coalescence in terms of the evolution of the void volume fraction f over the number of loading cycles at various points along the crack front (z-axis). Note that the highest values occur at z ¼ 6:5 mm.

122

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

Fig. 46. Simulated residual stress T 22 in the ligament after the 8th unloading stage.

Consider next the residual compressive stress in the ligament at the crack tip in Figs. 45 and 46 at unloading. This stress is regarded as the main cause for the WPS effect [16,17]. The current simulation results imply that, after only one loading cycle (Fig. 45), this stress is already quite large, and of comparable magnitude as of the crack-opening stress (up to 1000 MPa). In the WPS case, this is related to specimen cooling in the LUCF path [16]. Further loading cycles lead to two effects. Firstly, the increase of plastic deformation (Fig. 43) and ductile damage (Fig. 44) reduces the number of possible sites for cleavage fracture. Secondly, Fig. 44 shows that the final void volume-fraction f f in an element at the crack tip is reached during the 7th loading cycle. When ductile crack propagation occurs, the crack grows in a place where the effect of preloading is smaller. In the process, the WPS-effect is also reduced which is in agreement with experimental investigations [15]. 14. Summary The presented procedure extends the Gurson–Tvergaard–Needleman ductile damage model by large deformation and hyperelasticity to deal with the effects of non-linear combined hardening during non-proportional and cyclic loading on damage development. Simulations of the behaviour of unit cells with voids motivated the enhancement of the Gurson yield function because they showed that kinematic hardening influences localisation on the mesoscopic level and therefore damage development. Implementation of the the material model into the finite-element program ABAQUS via the user interface UMAT enabled its identification using experimental data on hardening and damage processes for the high strength steel 10MnMoNi5-5 at room temperature and at 300 °C. The application of an evolution strategy solved the corresponding optimisation problem. Numerical data on smooth and notched round bars in comparison with experimental results verified the model identification by means of force–displacement and crack-opening resistance curves for different specimen geometries. Finally, the identified model showed its capabilities as to predict blunting effects at the crack tip during cyclic loading and, moreover, the development of the residual compressive stress state that decreases the likelihood of cleavage fracture. Formally, this analogously applies for the warm-pre-stress (WPS) effect in which thermomechanical cyclic loading leads to an apparent increase in fracture toughness. The introduced approach might serve as the basis for modelling the WPS effect in the future. Acknowledgements The authors thank the Federal Ministry of Education and Research of the Federal Republic of Germany for funding under Grant No. 1500 986 and Mr. A. Eberle of BAM, Berlin, for preparing the final version of this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

Chu C, Needleman A. Void nucleation effects in biaxially stretched sheets. J Engng Mater Technol 1980;102:249–56. Needleman A, Tvergaard V. An analysis of ductile rupture in notched bars. J Mech Phys Solids 1984;32:461–90. Tvergaard V, Needleman A. Analysis of the cup-cone fracture in a round tensile bar. Acta Metall 1984;32(1):43–68. Koplik J, Needleman A. Void growth and coalescence in porous plastic solids. Int J Solids Struct 1988;24:835–50. Gurson A. Continuum theory of ductile rupture by void nucleation and growth: part i yield criteria and flow rules for porous ductile media. J Engng Mater Technol 1977;99:2–15. Mear M, Hutchinson J. Influence of yield surface curvature on flow localization in dilatant plasticity. Mech Mater 1985;4(3-4):395–407. Leblond J-B, Perrin G, Devaux J. An improved Gurson-type model for hardenable ductile metals. Euro J Mech A/Solids 1995;4:499–527. Lievers WB, Pilkey AK, Worswick MJ. The co-operative role of voids and shear bands in strain localization during bending. Mech Mater 2003;35 (7):661–74. Tvergaard V. Effect of yield surface curvature and void nucleation on plastic flow localization. J Mech Phys Solids 1987;35(1):43–60. Besson J, Guillemer-Neel C. An extension of the Green and Gurson models to kinematic hardening. Mech Mater 2003;35:1–18. Besson J. Continuum models of ductile fracture: a review. Int J Dam Mech 2010;19:3–52.

D. Klingbeil et al. / Engineering Fracture Mechanics 160 (2016) 95–123

123

[12] Arndt S, Svendsen B, Klingbeil D. Modeling of the residual stress development at a crack tip with a damage model (in German). Tech Mech 1997;4:323–32. [13] Arndt S. Modeling of warm-prestressing in the steel 10MnMoNi5-5 using a damage model with combined hardening (in German) Ph.D. thesis. University of Magdeburg; 2000. [14] Rechenberg I. Optimization of thechnical systems according to principles of biological evolution (in German) Ph.D. thesis. Technische Universität Berlin; 1970. [15] Blumenauer H, Krempe M, Eichler B, Ude J. Micromechanical modelling of the warm prestress (WPS effect). In: 2nd European conference on the mechanics of materials (EUROMECH-MECAMAT). University of Magdeburg; 1998. p. 231–7. [16] Klingbeil D, Arndt S, Häcker R, Hünecke J, Reusch F. Evaluation of material mechanics during fast cooling operations in pressurized building components with postulated cracks (in German), Tech. Rep. 234. Federal. [17] Wallin K. Master curve implementation of the warm pre-stress effect. Engng Fract Mech 2003;70:2587–602. [18] Anderson T. Fracture mechanics: fundamental and applications. 2nd ed. London: CRC Press LLC; 1995. [19] ASTME399-90. Standard test method for plane-strain fracture toughness of metallic materials. ASTM-Handbook 3.01. [20] Bluhm J, Morrissey R. Analysis of the cup-cone fracture in a round tensile bar. First international conference on fracture, vol. 3. p. 1739. [21] Goods S, Brown L. The nucleation of cavities by plastic deformation. Acta Metall 1979;27:1–15. [22] Marini B, Mudry F, Pineau A. Experimental study of cavity growth in ductile rupture. Engng Fract Mech 1985;22:989. [23] Tvergaard V. Influence of voids on shear band instabilities under plane strain conditions. Int J Fract 1981;17:389–407. [24] Tvergaard V. Material failure by void coalescence in localized shear bands. Int J Solids Struct 1982;18:659–72. [25] Perrin G, Leblond J. Analytical study of a hollow sphere made of plastic porous material and subjected to hydrostatic tension application to some problems in ductile fracture of metals. Int J Plast 1990;6:677–99. [26] Anderson H. Analysis of a model for void growth and coalescence ahead of a moving crack tip. J Mech Phys Solids 1977;25:217–33. [27] Devaux J, Gologanu M, Leblond J, Perrin G. On continued void growth in ductile metals subjected to cyclic loadings. In: IUTAM conference on nonlinear analysis of fracture. Cambrige University Press; 1996. [28] Christofferson J, Hutchsinson J. A class of phenomenological corner theories of plasticity. J Mech Phys Solids 1979;27:465–87. [29] Hutchinson J, Tvergaard V. Shear-band formation in plane strain. Int J Solids Struct 1981;17:451–70. [30] Becker R, Needleman A. Effect of yield surface curvature on necking and failure in porous plastic solids. J Appl Mech 1986;53:491–9. [31] Lammering R, Pecherski B, Stein E. Theoretical and computational aspects of large plastic deformations involving strain-induced anisotropy and developing voids. Arch Mech 1990;42:347–75. [32] Paulun J, Pecherski R. On the application of the plastic spin concept for the description of anisotropic hardening in finite deformation plasticity. Int J Plast 1987;3:303–14. [33] Willis J. On methods for bounding the overall properties of nonlinear composites. J Mech Phys Solids 1991;39:73–86. [34] Ponte-Castaeda P. The effective mechanical properties of nonlinear isotropic materials. J Mech Phys Solids 1991;39:45–71. [35] Michel J, Suquet P. The constitutive law of non-linear viscous and porous materials. J Mech Phys Solids 1992;40:783–812. [36] Tvergaard V, Huang Y, Hutchinson J. Cavitation instabilities in a power hardening elasto-plastic solid. Euro J Mech A/Solids 1992;11:215–31. [37] Kuna M, Sun D. Three-dimensional cell model analysis of void growth in ductile materials. Int J Fract 1996;81:235–58. [38] Rabold F, Kuna F. Cell model simulation of void growth in nodular cast iron under cyclic loading. Comput Mater Sci 2005;32(3–4):489–97. [39] Lee E. Elastic–plastic deformation at finite strains. J Appl Mech 1969;36:1–16. [40] Svendsen B. On the modeling of anisotropic elastic and inelastic material behaviour at large deformation. Int J Solids Struct 2001;38:9579–99. [41] Chaboche J-L. Time-independent constitutive theories for cyclic plasticity. Int J Plast 1986;2:149–88. [42] Lemaitre J, Chaboche J-L. Mechanics of solids. Paris: Dumond; 1990. [43] Svendsen B. A thermodynamic formulation of finite-deformation elastoplasticity with hardening based on the concept of material isomorphism. Int J Plast 1998;14:473–88. [44] Svendsen B, Arndt S, Klingbeil D, Sievert R. Hyperelastic models for elastoplasticity with non-linear isotropic and kinematic hardening at large deformation. Int J Solids Struct 1998;35(25):3363–89. [45] Svendsen B. A logarithmic-exponential backward-euler-based split of the flow rule for anisotropic inelastic behaviour at small elastic strain. Int J Numer Meth Engng 2007;70(4):496–504. [46] Eterovic A, Bathe K. A hyperelastic-based large strain elastoplastic constitutive formulation with combined isotropic-kinematic hardening using logarithmic stresses and strain measures. Int J Numer Meth Engng 1990;30:1099–115. [47] Weber G, Anand L. Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids. Comp Meth Appl Mech Engng 1990;79:173–202. [48] Simo JC, Hughes TJR. Computational inelasticity. Springer-Verlag; 1998. [49] Silhavy M. The mechanics and thermodynamics of continuous media. Springer-Verlag; 1997. [50] Bridgman P. Studies in large plastic flow and fracture. McGraw Hill; 1952. [51] Eigenmann B, Macherauch E. X-ray investigation of the stress state in materials (in German). Mat-wiss u Werkstofftech 1995;26:148–60. [52] Krempe M. X-ray measurements of residual stresses in pre-deformed C(T)25 specimens (in German), Tech. Rep. 014/99. Institute for Material Technology and Testing; 1999.