A ductile failure model applied to the determination of the fracture toughness of welded joints. Numerical simulation and experimental validation

Engineering Fracture Mechanics 73 (2006) 2756–2773 www.elsevier.com/locate/engfracmech

A ductile failure model applied to the determination of the fracture toughness of welded joints. Numerical simulation and experimental validation I. Pen˜uelas *, C. Betego´n, C. Rodrı´guez Departamento de Construccio´n e Ingenierı´a de Fabricacio´n, Universidad de Oviedo, Edificio Departamental Oeste, Desp. 23, Bloque 5, 33203 Gijo´n, Spain Received 11 November 2005; received in revised form 28 April 2006; accepted 5 May 2006 Available online 19 June 2006

Abstract A ductile-failure model for analysing the fracture behaviour of welded joints has been implemented. Finite element analyses of mismatched welded joints have been performed using the computational cell methodology applied to SE(B) specimens. Different crack lengths, material mismatching, and widths of weld metal have been considered. Ductile parameters have been experimentally and numerically obtained. The influence of geometry and material mismatching on the fracture behaviour of cracked welded joints has been validated by means of a testing program. In addition, the experimental results have been explained through the crack tip constraint, which has been numerically determined.  2006 Elsevier Ltd. All rights reserved. Keywords: Ductile fracture; Constraint parameter; Welded joints

1. Introduction Ductile fracture of metallic materials involves micro-void nucleation and growth, and final coalescence of neighbouring voids to create new surfaces of a macro-crack. The ductile failure process for porous materials is often modelled by means of the Gurson model [1], which is one of the most widely known micro-mechanical models for ductile fracture, and describes the progressive degradation of material stress capacity. In this model, which is a modification of the von Mises one, an elastic–plastic matrix material is considered and a new internal variable, the void volume fraction, f, is introduced. Although the original Gurson model was later modified by many authors, particularly by Tvergaard and Needleman [2–4], the resultant model is not intrinsically able to predict coalescence, and is only capable of simulating micro-void nucleation and growth. This deficiency is solved by introducing an empirical void coalescence criterion: coalescence occurs when a critical void volume fraction, fc, is reached. The value of fc is *

Corresponding author. Tel.: +34 985181980; fax: +34 985181945. E-mail address: [email protected] (I. Pen˜uelas).

0013-7944/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2006.05.007

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2757

Nomenclature a crack length a0 initial crack length Da crack length increment A strain controlled nucleation rate B specimen thickness BL slope of the blunting line CMOD crack mouth opening displacement CTOD crack tip opening displacement e total elongation E Young’s modulus f current void volume fraction f* modified void volume fraction f0 initial void volume fraction fc critical void volume fraction fF void volume fraction at final failure fN void volume fraction of nucleating particles in the Gaussian distribution of the nucleation rate fu ultimate void volume fraction f_ growth void volume fraction growth rate f_ nucleation ; f_ n void volume fraction nucleation rate h weld semi-width H weld width J J-integral KI stress intensity factor in mode I of load l0 initial element size at the fracture process zone m mismatch ratio (level of mismatching) n strain hardening exponent p macroscopic hydrostatic stress P applied load q macroscopic von Mises effective stress q1, q2, q3 fitting parameters introduced by Tvergaard and Needleman Q constraint Q-parameter (triaxility parameter) r, h polar coordinates with origin at the crack tip rvoids void space ratio R R-curve (resistance curve) S Distance between the support cylinders Sij deviatoric components of the Cauchy stress tensor SN standard deviation in the Gaussian distribution of the nucleation rate T T-stress Tg T-stress due to geometry Tm T-stress due to material mismatching (ficticious) W specimen width Z reduction of area Greek symbols a proportionality parameter in the Ramberg–Osgood law bg geometrical constraint parameter (biaxiality parameter) bm constraint parameter due to material mismatching bT total constraint parameter

2758

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

components of the plastic strain rate tensor e_ pij ep equivalent plastic strain e_ p equivalent plastic strain rate e0 strain at yield stress e1, e2, e3 principal strains eN mean strain in the Gaussian distribution of the nucleation rate m Poisson’s ratio r0 yield stress r0b yield stress of the base metal r0w yield stress of the weld metal r1 current maximum principal stress rm macroscopic mean stress rR tensile strength  r flow stress of the matrix material U yield function of Gurson–Tvergaard–Needleman s normalised T-stress

selected beforehand or is numerically fitted from tension tests [5], and is introduced in the model as a constant value which depends on the material. This criterion was initially validated by several numerical studies [3,6]. However, later theoretical, numerical and experimental studies [7,8] suggested that the coalescence criterion should include some micro-structural information related to the void/ligament rate, the stress state and the geometry of the specimen, since different values of fc can be obtained from the same tension tests [9–14]. The plastic limit load proposed by Thomason [9] for void coalescence improves the prediction of ductile fracture. Accordingly to Thomason, on a surface of ductile fracture, micro-void coalescence is the result of the failure by plastic limit load (microscopic internal necking) of the matrix between voids. Thus, the localised deformation state of void coalescence is completely different to the homogeneous deformation state during void nucleation and growth. It is necessary to take into account both types of deformation in the analysis of ductile fracture, coalescence depending on the competition between these two types of deformation. As deformation begins, voids are small and the stress needed to follow a homogeneous deformation is lower than that required to follow a localised deformation. As the plastic deformation grows and the void volume fraction increases, the stress required for the localised deformation mode decreases. When this stress equals that needed for the homogeneous mode, a fork point is reached and coalescence takes place. Assuming that voids are initially spherical and that they evolve in a spherical way, a ductile failure model is formulated, based on the so-called complete Gurson model [13]. This complete model combines the Gurson–Tvergaard–Needleman model and the coalescence criterion proposed by Thomason, and it has been shown to give accurate predictions for any level of stress triaxiality, for both strain non-hardening and strain hardening materials. The implementation of such model requires the integration of the constitutive equations. This integration requires the calculation of the stress values and the different plastic variables in a finite strain increment, and different integration procedures have been developed, using explicit and implicit methods. In the former, the yield function and the plastic variables are evaluated at known stress states, and although it is not necessary to use an iterative procedure, it is common to restore the final stress and the plastic variables to the yield surface by means of an iterative correction, since this condition is not forced by integration [15]. In the implicit methods [16], stresses are unknown and so a system of non-linear equations should be solved, with iterative methods usually being employed for this purpose. Since the yield criterion is used for defining the system of equations to be solved, the obtained stresses directly satisfy the yield criterion. In this work, an implicit method for integrating the constitutive equations that describe the ductile failure process of a metallic material in a welded joint, is implemented. The model is based on the complete Gurson model and employs Runge–Kutta based methods, which have been shown to provide good results in elastoplasticity [17], to solve the resulting system of differential equations. The model has been implemented in FORTRAN and it has been introduced in the ABAQUS finite element commercial code by means of the

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2759

UMAT user subroutine. Variable update algorithm is based on the return mapping algorithms [18], and a robust iteration scheme is used. To simulate the extension of the crack, the computational cell methodology proposed by Xia and Shih [19– 22] is used. This methodology includes a realistic void growth mechanism and a micro-structural length scale physically coupled to the size of the fracture process zone. Thus, according to these authors, under Mode I ductile fracture occurs within a thin layer of material symmetrically located about the crack plane. This layer consists of cell elements of height l0, where l0 is associated with the mean spacing of the larger void-initiating inclusions. In this paper, progressive damage and subsequent macroscopic material softening in each cell is modelled using the complete Gurson model. Finally, when the void volume fraction reaches a specific value of final failure, the crack extends in a discrete manner by a distance of one cell. Within this context, finite element analyses of SE(B) specimens with mismatched welded joints have been performed. Different crack lengths, widths of the weld material, and mechanical material mismatching have been considered. Ductile parameters have been experimentally and numerically obtained. The influence of geometry and material mismatching on the fracture behaviour of cracked welded joints, which has been analysed in previous papers [23], has been here validated by means of a testing program, and the results have been explained through the crack tip constraint, which has been numerically determined. 2. Constitutive model: the complete Gurson model To describe the evolution of void growth and subsequent macroscopic material softening in computational cells, the yield function of Gurson [1], modified by Tvergaard [2,3] and Tvergaard and Needleman [4,24], is used in this work. This modified yield function is defined by an expression in the form:   q2 3  q2  p ; f Þ ¼ Uðq; p; r þ 2  q1  f   cosh  ð1Þ  ð1 þ q3  f 2 Þ ¼ 0   r 2r  is the flow stress of the matrix material (microscopic), where p = rm with rm the macroscopic mean stress, r  ¼ H ðep Þ, f is which relates to the equivalent plastic strain by means of an uniaxial equation in the form r * current void volume fraction, f is the modified void volume fraction (related to f), and q is the macroscopic von Mises effective stress given by rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3  ðS ij  S ij Þ q¼ ð2Þ 2 where Sij denotes the deviatoric components of the Cauchy stress. The constants q1, q2 and q3 are fitting parameters introduced by Tvergaard [2,3] to provide better agreement with the results of detailed unit cell calculations. For metallic materials the usual values of these constants are q1 = 1.5, q2 = 1, q3 ¼ q21 . The modified void volume fraction, f *, was introduced by Tvergaard and Needleman [4] to predict the rapid loss in strength that accompanies void coalescence, and is given by ( f if f 6 fc   f f ¼ ð3Þ c fc þ ffuF f  ðf  f Þ if f > fc c c where fc is the critical void volume fraction, fF is the void volume fraction at final failure, which is usually fF = 0.15, and fu ¼ 1=q1 is the ultimate void volume fraction.  and f. Thus the evolution law for the void volume The internal variables of the constitutive model are r fraction is given in the model by an expression in the form: f_ ¼ f_ growth þ f_ nucleation

ð4Þ

The void nucleation law implemented in the current model takes into account the nucleation of both small and large inclusions. The nucleation originated at large inclusions is stress controlled, and it is assumed that large inclusions are nucleated at the beginning of the plastic deformation, and are considered as the initial void volume fraction.

2760

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

The nucleation of smaller inclusions is strain controlled and, according to Chu and Needleman [25], the nucleation rate is assumed to follow a Gaussian distribution, which is f_ nucleationsmall particles ¼ A  e_ p ð5Þ where e_ p is the equivalent plastic strain rate, and  2 ! fN 1 ep  eN pffiffiffiffiffiffiffiffiffi  exp   A¼ 2 SN SN  2  p

ð6Þ

where eN is mean strain, SN is the standard deviation and fN is the void volume fraction of nucleating particles. The growth rate of the existing voids can be expressed as a function of the plastic strain rate in the form: f_ growth ¼ ð1  f Þ  e_ p ð7Þ kk

e_ pij

is the plastic strain rate tensor. where As was stated in the introduction, in the modified Gurson–Tvergaard–Needleman model the critical void volume fraction, fc, is a material constant. However, in the plastic limit load criterion proposed by Thomason [9] for void coalescence, fc is not a material constant, but the material response to coalescence. Thus, in the complete Gurson model [13,26,27] the coalescence criterion can be written in the form: 8    2 > bcoalescence r1 1 > p ffiffiffiffiffiffiffi > < a   1 þ  ð1  p  r2voids Þ ) no coalesc: occurs coalescence < r rvoids rvoids   ð8Þ  2 > > 2 > r1 ¼ acoalescence  1  1 þ bcoalescence pffiffiffiffiffiffiffi Þ ) coalesc: occurs  ð1  p  r : voids  rvoids rvoids r where the first term is related to the homogeneous deformation state, and the second one to the plastic limit load needed for coalescence. In the previous expression, r1 is the current maximum principal stress, acoalescence = 0.1, bcoalescence = 1.2, and rvoids is the void space ratio, which is given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3f

3  expðe þ e þ e Þ 1 2 3 4p pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rvoids ¼ ð9Þ expðe2 þe3 Þ 2

with f, the actual void volume fraction; and e1, e2 and e3, the principal strains. The complete Gurson model has been verified by Zhang et al. [13] for non-hardening materials (acoalescence = 0.1). In addition, Pardoen and Hutchinson [26] have shown that in the case of a Ramberg– Osgood material with a strain hardening exponent n, and considering that voids remain always spherical, better predictions are obtained with acoalescence ¼ 0:12 þ

1:68 ; n

bcoalescence ¼ 1:2

ð10Þ

These are the values used in this paper. 3. Ductile fracture in mismatched welded joints 3.1. Crack tip stress fields in welded joints When analysing welded joints, the heterogeneity of the joint is a key factor for understanding the fracture behaviour of any welded structure containing cracks. In an ‘‘ideal’’ weld, with a crack contained within the weld material and running along the material’s centre-line, parallel to the weld-base material interface, and where the effect of the heat affected zone is negligible [28–30], a two material idealisation of the weld structure can be taken into account. This consists of weld material with yield stress r0w and width 2h, which contains the crack, and base material with yield stress r0b. Fig. 1 shows this idealisation for a short-crack bend specimen. In this case, the level of material mismatching can be defined by the mismatch ratio m, in the form: r0w ð11Þ m¼ r0b

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2761

P

σ 0w W

σ 0b

a = W/2 a 2h

B=W

S = 4W Fig. 1. Two material idealisation of a three-point bend mismatched specimen.

with m > 1 referring to material strength overmatching, m < 1 to undermatching and m = 1 corresponding to a homogeneous specimen of weld metal. Using the finite element method and the slip-line field theory, Hao et al. [31], Joch et al. [32] and Burstow and Ainsworth [33] have demonstrated that material strength mismatching significantly affects both the stress fields and the crack resistance curves of tension and bend specimens [29,34,35]. That is, constraint is not only a function of geometry, but also of material mismatching. Thus, in addition to any geometrical constraint, constraint due to material mismatching should be taken into account [36]. Geometrical constraint can be quantified by means of both the elastic T-stress [37–40], which directly characterises the geometrical constraint effect, and the non-linear Q-parameter [41,42], which is a direct measure of the elastic–plastic stress fields that can be related with the HRR field [43,44], and it usually describes the deviation of the stress field from a reference stress state at a specified position ahead of the crack tip, thus Q¼

rhh  rRef hh r0

at

r  r0 ¼ 2; h ¼ 0 J

ð12Þ

where r0 is the yield stress of the material and rRef hh is the stress distribution for which T = 0. On the basis of a dimensional argument O’Dowd and Shih [41] demonstrated that, at least under small scale yielding (SSY), for each material Q and T are univocally related by expressions of the form: Q ¼ A  s þ B  s 2 þ C  s3

ð13Þ

where s = T/r0 is the normalised T-stress. The evolution of T with the load can be characterised by the socalled biaxiality parameter, or geometrical constraint parameter, bg; so that, for each geometry and at each instant, the value of the T-stress is defined by an expression in the form: bg  K I T ¼ pffiffiffiffiffiffiffiffiffi pa

ð14Þ

Geometries with positive bg values lead to positive T-stresses and raise the stress fields slightly. Conversely, geometries with negative bg values lead to negative T-stresses and lower the stress fields significantly. Constraint due to material mismatching has been analysed in depth by Burstow et al. [34], who used boundary layer formulations to investigate the crack tip stress fields in different mismatched cases. They also found that, even in the case of a null T-stress, T = 0, the development of the stress distribution around the crack tip depends not only on the applied load, but also on the level of mismatching, m, and the weld material width, 2h. In the case of overmatching (m > 1) the crack tip stress fields are lowered relative to those obtained in homogeneous weld material; conversely, in the case of undermatching (m < 1) the stress fields are raised relative to the homogeneous reference situation. In all cases, the higher the level of mismatching and the thinner the weld metal strip are, the more severe the former effects. In order to quantify constraint due to material mismatching, different methodologies, based on that utilised for the quantification of geometrical constraint, have been defined. Thaulow and co-workers [45–47] defined a procedure similar to the J–Q one. Thus, they established a J–Q–M formulation for the crack tip stress fields.

2762

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

In this three-parameter formulation, J is related with the load, Q with the geometrical constraint and M with the material constraint. Later, Burstow et al. [34] defined a normalised load parameter h Æ r0w/J that scales the plastic zone with the width of the weld region. Thus, for a given mismatch ratio, the crack tip stress fields depend only on this normalised load parameter. By means of this parameter, they also established a J–Q– M formulation in order to examine the crack tip stress fields, and related Q and m for the overmatching case. Thaulow and co-workers [47,48], also extended the model to interface cracks and three-material problems where the HAZ was taken into account. Detailed analyses of the crack tip stress fields have also been performed by Burstow et al. [34,35], who explored how both brittle and ductile fracture can be influenced by material mismatching, observing that fracture toughness for this kind of idealised overmatched welds is, in general, higher than that for undermatched specimens. On the other hand, Betego´n and Pen˜uelas [23] defined a procedure similar to the J–T one. Thus, they developed a J–Tg–Tm formulation where J is related with the load, Tg with the geometrical constraint and Tm with the material constraint. By means of finite element analyses of plane strain crack tip stress fields from homogeneous and heterogeneous modified boundary layer formulation, as well as homogeneous and mismatched full field solutions, they established a new constraint parameter, bm, for overmatched welded joints, that quantifies the material mismatching effect on the crack tip stress fields. In the case of complete specimens, both geometry and material mismatching affect the crack tip stress fields and they also defined a total constraint parameter bT, obtained by addition of bg and bm. In the case of undermatching, the constraint parameter bm cannot be defined. Although this three-parameter formulation of the stress fields in welded joints is defined for small scale yielding conditions, it can be extended to large scale yielding conditions for some configurations and enlarged the validity of other approaches. 3.2. Ductile fracture in welded joints Due to the nature of the welding process, the welded zone contains, in general, more flaws and defects than the surrounding base material, and it is a critical place for fracture. Although fracture can develop by mechanisms of cleavage and ductile tearing, in steels at room temperature under plain strain conditions, the crack usually grows by ductile tearing, and the fracture behaviour can be characterised by the ductile crack resistance curve. In addition, any typical steel welded joint contains both small and large inclusions. Thus, the fracture behaviour of welded joints can be analysed by means of the complete Gurson model described in the previous sections. Within the context of the computational cell methodology proposed by Xia and Shih [19–22], finite element analyses of three-point bend specimens with single edge cracks of different lengths (a/W = 0.1, 0.2, 0.5), and values of the weld semi-width h ranging from 1 mm to 10 mm, have been performed. The ABAQUS [49] finite element commercial code has been used. Simulations have been carried out under plane strain conditions, and they have accounted for large strains around the crack tip, and contacts. Four-node bilinear, hybrid elements (CPE4H) have been used. Table 1 gives details of the yield stress of the materials and the range of material mismatch conditions analysed. Fig. 2(a) shows two finite element meshes used to model the problem; also shown are the support cylinder and the cylinder where the load is applied, which are modelled as rigid bodies. The mesh represents one-half of the specimen, since symmetry has been applied. Fig. 2(b) shows details of the mesh at the crack tip region; layers of uniformly sized void-containing elements are also shown. Table 1 Material mismatch conditions analysed Configuration

m

r0w (MPa)

r0b (MPa)

Homogeneous 20% overmatched 40% overmatched 60% overmatched 20% undermatched 40% undermatched

1 1.2 1.4 1.6 0.8 0.6

625 625 625 625 625 625

625 520.8 446.4 390.6 781.3 1041.7

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2763

Fig. 2. Meshes for the three-point bend models: (a) general meshes and (b) the crack tip region and ductile cells (in grey).

Fig. 3. Resistance curves for short crack (a/W = 0.2) specimens with l0 = 0.1 mm, weld semi-width h = 1 mm and different levels of material mismatching, m.

The selected material parameters correspond to structural steels with moderate strength, hardening and toughness. All the materials are assumed to follow a Ramberg–Osgood law written in the form:  n e r r ¼ þa ð15Þ e 0 r0 r0 with e0 = r0/E, Young’s modulus E = 200 000 MPa, strain hardening exponent n = 10 and a = 1.56. The Poisson’s ratio is assumed m = 0.3.

2764

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

In the fracture process zone the micro-structural model parameters are: initial porosity f0 = 5.0 · 104, which includes the void volume fraction inherent to the welding process and the large inclusions volume fraction; Gaussian distribution of the nucleation rate of small inclusions with void volume fraction of nucleating particles fN = 0.002, mean strain eN = 0.3 and standard deviation SN = 0.1; fitting parameters of the yield

Fig. 4. Resistance curves for short crack (a/W = 0.2) specimens with l0 = 0.1 mm, material mismatching m = 1.6 and different weld semiwidths, h.

Fig. 5. Resistance curves for short crack (a/W = 0.2) specimens with l0 = 0.1 mm, material mismatching m = 0.6 and different weld semiwidths, h.

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2765

function of Gurson, Tvergaard and Needleman q1 = 1.5, q2 = 1.0 and q3 ¼ q21 ; void volume fraction at final failure fF = 0.15; and different initial element dimensions l0 = 0.05 mm, l0 = 0.1 mm and l0 = 0.2 mm. First of all, the results obtained from three-point bend specimen solutions for different mismatching and geometrical configurations are shown. Fig. 3 shows resistance curves for different levels of material

Fig. 6. Fracture toughness at two different crack growth values as a function of the constraint parameter.

Fig. 7. Resistance curves for long crack (a/W = 0.5) specimens with l0 = 0.1 mm, different levels of material mismatching, m, and different weld semi-widths, h.

2766

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

mismatching, m, a fixed weld semi-width, h = 1 mm, a short crack of a/W = 0.2, and l0 = 0.1 mm. Also shown are the values of the total constraint parameter, bT, for the different configurations. In the case of material strength overmatching, the reduction in constraint results in a significant increment in weld resistance to fracture. Both, fracture initiation and fracture propagation (R-curve gradient) are increased. In the case of material strength undermatching, the increment in constraint results in a significant reduction in weld resistance to fracture. Both, fracture initiation and fracture propagation (R-curve gradient) are reduced. In all cases, the higher the level of mismatching and the thinner the weld width are, the more severe the former effects. Figs. 4 and 5 show the resistance curves and the bT values for different weld semi-widths for the overmatched case m = 1.6 and the undermatched case m = 0.6, respectively, for l0 = 0.1 mm. For the overmatched configurations, the bT values are also shown. In addition, the smaller the ratio a/W is, and in consequence the lower the constraint is, the higher the fracture resistance. This crack size effect can be explained in terms of the geometrical constraint that affects the crack tip stress fields [50]. Fig. 6 shows, for the overmatched configurations, the J values, at crack extensions of 0.2 mm and 1.0 mm, respectively, as a function of the total constraint parameter. The effect of bT on the crack resistance curves is the same as on the crack tip stress fields. Thus, the more negative the constraint parameter is, the higher the J value at fracture initiation. At a crack extension of 1.0 mm this constraint effect is higher than at a crack extension of 0.2 mm, since bT also affects the gradient of the resistance curves. In geometries with long cracks, that is to say with positive or higher values of the constraint parameter, the influence of varying both material mismatching and the weld width is lighter than in the case of geometries with quite negative values of the constraint parameter. Fig. 7 shows the resistance curves for a three-point bend specimen with a long crack of a/W = 0.5, different mismatched situations and different weld semi-widths. Although the effect on the resistance curve of a 60% overmatched welded joint with h = 1 mm is quite noticeable, bigger weld semi-widths (h = 5 mm) lead to resistance curves close to the one of the homogeneous all weld case, and for h = 10 mm the resistance curve is close enough to a homogeneous curve to accept that the welded joint behaves as homogeneous all weld material. What is more, for a 40% undermatched welded joint, the effect is negligible for big semi-widths (h = 5 mm and h = 10 mm), but it is important for smaller semi-widths (h = 2 mm), where a profound reduction in fracture resistance is observed. As was stated in the introduction, the size of elements in the fracture process zone, l0, is associated with the mean spacing of the larger void-initiating inclusions. Thus, the fracture resistance curves obtained with different meshes, where the initial element size varies, have to be affected by this geometrical dimension l0. Fig. 8

Fig. 8. Resistance curves for short crack (a/W = 0.2) homogeneous specimens obtained with different cell element sizes.

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2767

shows the resistance curves for a homogeneous specimen obtained with cell elements of l0 = 0.05 mm, l0 = 0.1 mm and l0 = 0.2 mm, other micro-structural parameters remaining the same. The bigger the element size is, the higher the fracture resistance to initiation and propagation (higher gradient of the resistance curve), the differences between each curve and the closest one of being nearly 40%. This figure reveals the importance of a correct material characterisation in order to obtain realistic predictions of the ductile fracture behaviour of any material under analysis. Different normalisations, applied to modified boundary layer formulations, have been found in literature [13]. However, although they have been taken into account in this paper, no possible normalisation has been found for full field solutions. 4. Experimental procedure and validation In order to analyse the transferability of the numerical results and the validation of the ductile failure model implemented, it is necessary to study the real fracture behaviour of different pre-cracked welded joints by means of experimental techniques. Thus, specimens extracted from overmatched and undermatched welded joints, obtained with one weld material and two base materials (MB1 and MB2), with short and long cracks and different widths of the weld strip, have been tested within a complete experimental procedure. In previous sections, computational studies have been carried out of different welded joints, with geometrical and material configurations, as well as micro-structural fracture parameters, suitable for understanding the influence of both geometry and material mismatching in the fracture behaviour of welded joints. Thus, welds of very small semi-widths have been simulated. However, in practice it is not possible to obtain right welds when welding the thick plates necessary to assure plain strain conditions (B  20 mm) if a very small root gap (2 Æ h = 2 mm) is required. On this basis, joints with mean widths of h = 5 mm and h = 10 mm (that is, H = 2 Æ h = 10 mm and H = 20 mm) have been obtained, although the effect of material mismatching is, in these situations, less noticeable. The welding procedure has been developed in order to obtain welded joints with no mechanical or geometrical discontinuities, and a heat affected zone small enough to consider a two-material (base and weld material) idealisation of the welded joint. An appropriate joint geometry is also required for obtaining specimens where the crack grows within the weld material. Thus, in the case of welded joints with big widths, a multi-run open square butt weld, with full penetration welded from both sides, was utilised. However, in the case of small widths, due to accessibility problems, a multi-run single V-butt weld with backing strip and full penetration was used. The optimal operational and metallurgical conditions have been reached using a semiautomatic mechanically operated welding process with active shielding gas (MAG – 88% Ar + 12% CO2), with welding wire type E 70 S6, welding current 140–200 A, welding voltage 20–25 V, welding speed 20–30 cm/min, heat inputs 10 kJ/cm and inter-run temperature lower than 250 C. Coupons have been inspected with non-destructive tests including visual inspection of macrographs, X-rays and ultrasonics, in order to reject the coupons with high levels of defects or with no negligible heat affected zones. Fig. 9(a) shows a macrograph of an undermatched H = 20 mm coupon; Fig. 9(b) shows a macrograph of an overmatched H = 10 mm coupon. Destructive tests have also been conducted on samples of the coupons for qualification purposes. Thus, the verification of the mechanical properties of the welded joints by means of tensile, hardness and bend tests has been carried out. The inclusion content of the weld material has also been

Fig. 9. (a) Macrograph of an undermatched H = 20 mm coupon and (b) macrograph of an overmatched H = 10 mm coupon.

2768

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

determined using quantitative microscopic optical techniques. Two different populations of inclusions, small and large, have been observed. In addition, the initial void volume fraction inherent to the welding process is also determined. Thus, an initial porosity value f0 = 1.0 · 104, and a volume fraction of void-nucleating particles fN = 0.004, were considered. The tensile properties of the weld and base metals have been determined, according to the ASTM E8M-04 standard [51], using round specimens extracted from the weld and base regions parallel to the weld direction. Table 2 gives details of the tensile properties at room temperature of the materials used. The mechanical behaviour of these materials within the plastic strain hardening zone has been fitted by a Ramberg–Osgood law (see Eq. (15)). Single edge bend specimens (SE(B)) extracted normal to the weld direction, have been used for estimating both the crack growth resistance (J–R) curves and the fracture toughness of the welded joints. The specimens have been notched, fatigue pre-cracked, and sidegrooved. The notch was located in the centre line of the weld strip with the crack growing in the thickness direction. Fig. 1 shows the bend test specimen geometry. Table 3 gives details of the different specimen configuration tested, which consist of four overmatched, m = 1.4, and three undermatched, m = 0.8, welded joints with weld mean widths H = 10 mm and H = 20 mm. Two crack lengths corresponding to a short crack, a/W = 0.24, and a long crack, a/W = 0.5, have been tested. In all the configurations the specimen width and thickness were W = 18 mm and B = W, respectively, and the span was S = 4 Æ W. For overmatched configurations, also shown are the values of the total constraint parameter, bT, defined by Betego´n and Pen˜uelas [23]. Fracture tests were carried out at room temperature, in accordance with the ASTM E1280-05a standard [52], and the single specimen method for the J–R curve determination has been used. This method allows Table 2 Tensile properties at room temperature: r0 – yield stress, rR – tensile strength, e – total elongation, Z – reduction of area Material

r0 (MPa)

rR (MPa)

e (%)

Z (%)

Weld MB1 MB2

450 320 576

530 427 607

26 31 13

70 69 61

Table 3 Specimen configurations tested, values of the total constraint parameter (bT), and slope of the blunting line (J = BL Æ CTOD) Configuration

r0w (MPa)

r0b (MPa)

a/W

H (mm)

bT

BL

1_Weld-MB1 2_Weld-MB1 3_Weld-MB1 4_Weld-MB1 5_Weld-MB2 6_Weld-MB2 7_Weld-MB2

450 450 450 450 450 450 450

320 320 320 320 576 576 576

0.24 0.24 0.5 0.5 0.24 0.24 0.5

10 20 10 20 10 20 10

0.255 0.237 +0.035 +0.061 – – –

906 1130 1160 1240 1290 1410 1500

Final fracture Final crack tip

Fatigue

Spray painted Stable crack growth Initial crack tip Fatigue pre-cracked Side groove

Initial notch

Fig. 10. Fracture surface of a short crack and H = 10 mm overmatched specimen.

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2769

the crack length during test to be indirectly obtained from the specimen compliance variation. However, since this method leads to big inaccuracies when testing specimens of very ductile materials and short-cracked specimens, a colour marking technique of the crack front has been applied. This technique allows additional information to be obtained from a single test. Thus, two physical measurements of the crack growth are determined during each test: an intermediate painting sprayed crack length, and the final crack length. Fig. 10 shows a detailed fracture surface of a short crack and thin width (H = 10 mm) overmatched specimen. The amount of stable crack growth has been obtained by means of post-test examinations of the fracture surface, and crack length has been obtained, according to ASTM standards [52], by a weighted mean of nine through-thickness crack length values measured on the fracture surface. In a post-analysis of the testing results, the crack growth values, Da, indirectly obtained using the compliance method [41], have first been corrected to match the values physically measured on the fracture surface, and then interpolated between these values. This procedure leads to lower (J–R) curve dispersion. It is important to note that, whereas in the experimental curves, where the crack growth values are obtained using the compliance method, is included the effect of blunting at the crack tip, in physical measurements of the crack growth this effect is not included. Thus, and in order to compare both types of measurements, blunting at the crack tip has been added to the later ones. After correction, the resistance curves obtained for each configuration from the individual pairs (J, Da), have been numerically fitted by expressions in the form [52,53]: C2

J ¼ C 1  ðDaÞ

ð16Þ

Fig. 11 shows the crack growth resistance curves obtained for the different configurations detailed in Table 3. In these curves, two different zones are represented: an initial zone that corresponds to the blunting behaviour of each configuration, where the experimental results have been corrected in accordance with the standards; and a second zone, where the experimental results have been fitted with expression (16). In order to apply the proposed ductile model to simulate ductile tearing in SE(B) specimens, various model parameters must be determined. The first set of parameters corresponds to the stress–strain law of the matrix

Fig. 11. Fitted experimental resistance curves for the tested specimens.

2770

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

material. The law has been obtained from the tensile test already described, and has been fitted by a Ramberg– Osgood law with n ffi 7 and a = 1 for MB1 and weld metals, and n ffi 20 and a = 1 for MB2 metal. The second set of parameters consists on q1, q2 and q3, which are related to the hardening of the matrix material. In this work, q1 and q2 are fixed based on cell calculations by Kopling and Needleman [6], so that q1 = 1.5, q2 = 1 and q3 ¼ q21 . The third set of parameters, eN, SN, fN and f0, is related to the void nucleation and the initial porosity value. Both, the void volume fraction of nucleating particles, fN, and the initial void volume fraction, f0, have been determined using quantitative microscopic optical techniques. Thus, they have been obtained from the total volume of inclusions due to debonding particles. Two different populations of inclusions have been observed. The first one consists on particles of small mean diameter (about 0.6 lm) uniformly distributed. The second one consists on bigger globular inclusions (mean diameter 7 lm) randomly distributed in the sample material. In this paper, the uniformly distributed small particles are considered to be nucleating particles, so fN has been chosen to correspond to the volume fraction of these particles, fN = 0.004. The remaining nucleation parameters have been chosen eN = 0.3 and SN = 0.1, in accordance to Chu and Needleman [25]. The initial void volume fraction, f0, is obtained from the large particles population, since, as it was pointed out in Section 2, the debonding between the matrix and the inclusions occurs as soon as load begins, so that the large inclusions volume fraction equals the initial void volume fraction and f0 = 1.0 · 104. The porosity value that determines the final failure, fF, has been set to fF = 0.15. Finally, as fc is not a material parameter, but determined from Eq. (8), the only unknown parameter is the critical length, l0, which corresponds to the size of the cell elements near the crack tip, that is, to the crack tip mesh size. This critical length is related to the material micro-structure. However, in this paper, as all the other parameters have been fixed, the mesh size has been determined by adjusting the numerical predictions to match the experimental fracture resistance curve of certain configurations, in this case a/W = 0.5, H = 10 mm, overmatched. And it has been held fixed for the other configurations. An alternative procedure would consists on fixing this critical length as the mean inter-large-particle spacing and then determine the values of q1, q2 and q3 by fitting the numerical and experimental data. It is necessary to point out that, because of the parameter selection procedure, it is the complete set of values which describe the ductile behaviour of the material. Thus, if any parameter is changed, the numerically fitted parameters must be recalculated. Once the mechanical properties and the ductile parameters were known, the geometries and mismatched configurations tested have been numerically simulated and their resistance curves have been obtained. Finally, in order to analyse the transferability of the numerical results and the validation of the ductile failure model implemented, the numerical and experimental results have to be compared. To do this, the numerical resistance curves need to be corrected in order to include the blunting effect of each mismatch configuration and geometry. Thus, to the crack growth value determined numerically, one-half of the CTOD value is added (Dacorrected = Da + CTOD/2), following a similar procedure to that followed for the experimental physical values. The CTOD values are obtained from: (a) the CMOD calculated in simulations, ductile failure being taken into account, and (b) the CTOD obtained from numerical simulations where, in order to simulate blunting at the crack tip, finite element meshes with an inner radius at the crack tip have been used, ductile failure not being taken into account. For each configuration, these two numerical CTOD values have been found to be the same, and they agree with those obtained experimentally. Table 3 shows the slopes of the blunting lines for the different configurations analysed. From this table it can be observed that, in the case of overmatching, the lower the constraint is, the lower the slope of the blunting line. Besides, in the case of undermatching, the slope of the blunting line grows with the relative size of the crack and the width of the weld strip. Once corrected, numerical and experimental curves can finally be compared. Fig. 12 shows the resistance curves obtained numerically and experimentally. The former are represented by isolated points, and the latter by continuum and dot lines. It has to be noted that all configurations with a similar resistance curve in Fig. 11 (experimentally determined), have been represented in Fig. 12 as a unique curve. From Fig. 12, it can be observed that constraint, measured through the total constraint parameter bT, qualitatively determines the crack growth resistance curve, even for loads out of the range of application of bT. Thus, the more negative the constraint parameter is, see Table 3, the higher the fracture initiation, the steeper the fracture propagation, and the lower the slope of the blunting line are. That is to say, the lowest constraint configuration tested, which corresponds to a short-crack overmatched specimen with thin weld width, leads to

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2771

Fig. 12. Comparison of the experimental and numerical resistance curves (corrected with the blunting line) for the tested specimens.

the highest fracture resistance values, followed by the short-crack overmatched specimen with thick weld width. The remaining configurations, which correspond to positive or close to zero total constraint values, have very similar constraint values, and both fracture initiation and propagation are close enough to each other to be considered the same for all of them, as was pointed out earlier. Thus, in the case of short-crack undermatched specimens, whatever the weld width, as well as in the case of long-crack specimens, whatever the weld width and the mismatch nature, a single resistance curve is obtained. Such behaviours are shown by both numerical and experimental results. In addition, very good agreement has been found between both kinds of results for all the configurations analysed. These results confirm the high dependence on constraint of both fracture toughness and ductile propagation, and show the ability of numerical models for quantifying fracture behaviour. 5. Conclusions An implicit method for integrating the constitutive equations for materials which behave according to the complete Gurson model has been implemented. A Runge–Kutta based method for solving differential equations has been used within an error controlled sub-stepping procedure. The implemented model has been applied to mismatched welded joints, and the key role of the constraint effects on ductile fracture has been analysed. These effects have been quantified by means of a total constraint parameter, bT. In addition, a validating experimental programme has been carried out, and good agreement between numerically predicted and testing-obtained values is shown. Thus, the ability of numerical models for predicting the fracture behaviour of complete geometries with welded joints not only qualitatively but also quantitatively, has been shown. Acknowledgements The authors acknowledge the financial support of the Spanish Ministry of Science and Technology, project MAT2000-0602, and the access to ABAQUS under academic license from Hibbitt, Karlsson and Sorensen.

2772

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

References [1] Gurson AL. Continuum theory of ductile rupture by void nucleation and growth. Part I. Yield criteria and flow rules for porous ductile media. J Engng Mater Tech 1977;99:2–15. [2] Tvergaard V. Influence of voids on shear bands instabilities under plane strain conditions. Int J Fract 1981;17:389–407. [3] Tvergaard V. On localization in ductile materials containing spherical voids. Int J Fract 1982;18:157–69. [4] Tvergaard V, Needleman A. Analysis of cup-cone fracture in a round tensile bar. Acta Metall 1984;32:57–169. [5] Sun DZ, Kienler R, Voss B, Schmitt W. Application of micro-mechanical models to the prediction of ductile fracture. In: Alturi SN, Newman Jr JC, Raju IS, Epstein JS, editors. Fracture mechanics, 22nd symposium. ASTM STP 1131, vol. II. Philadelphia: ASTM; 1992. p. 368–78. [6] Koplik J, Needleman A. Void growth and coalescence in porous plastic solids. Int J Solids Struct 1998;24:835–53. [7] Steiglich D, Brocks W. Micromechanical modelling of damage and fracture of ductile materials. Fatigue Fract Engng Mater Struct 1998;21:1175–88. [8] Zhang ZL, Hauge M. On the Gurson parameters. In: Panontin TL, Sheppard SD, editors. Fatigue and fracture mechanics, vol. 29. American Society for Testing and Materials (ASTM STP 1321), 1998. [9] Thomason PF. Ductile fracture of metals. Oxford: Pergamon Press; 1990. [10] Thomason PF. Ductile fracture by the growth and coalescence of microvoids of non-uniform size and spacing. Acta Metall et Mater 1993;41:2121–34. [11] Thomason PF. A three-dimensional model for ductile fracture by the growth and coalescence of microvoids. Acta Metall 1985;33:1087–95. [12] Thomason PF. A view of ductile-fracture modelling. Fatigue Fract Engng Mater Struct 1998;21:1105–22. [13] Zhang ZL, Thaulow C, Ødegard J. A complete Gurson model approach for ductile fracture. Engng Fract Mech 2000;67:155–68. [14] Zhang ZL. A complete Gurson model. In: Aliadadi MH, editor. Nonlinear fracture and damage mechanics. Computational Mechanics Publications; 1998. [15] Sloan SW. Substepping schemes for the numerical integration of elastoplastic stress strain relations. Int J Numer Meth Engng 1987;24:893–911. [16] Ortiz M, Simo JC. An analysis of a new class of integration algorithms for elastoplastic constitutive relations. Int J Numer Meth Engng 1986;23:353–66. [17] Buttner J, Simeon B. Runge–Kutta methods in elastoplasticity. Appl Numer Math 2002;41:443–58. [18] Narasimhan R, Rosakis AJ, Moran B. A three-dimensional investigation of fracture initiation by ductile failure mechanisms in a 4340 steel. I. J Fract 1992;56:1–24. [19] Xia L, Shih CF. Ductile crack growth – I. A numerical study using computational cells with microstructurally based length scales. J Mech Phys Solids 1995;43:233–59. [20] Xia L, Shih CF. Ductile crack growth – II. Void nucleation and geometry effects on macroscopic fracture behaviour. J Mech Phys Solids 1995;43:1953–81. [21] Xia L, Shih CF. Ductile crack growth – III. Transition to cleavage incorporating statistics. J Mech Phys Solids 1996;44:603–39. [22] Xia L, Shih CF, Hutchinson JW. A computational approach to ductile crack growth under large scale yielding conditions. J Mech Phys Solids 1995;43:398–413. [23] Betego´n C, Pen˜uelas I. A constraint based parameter for quantifying the crack tip stress fields in welded joints. Engng Fract Mech, in press, doi:10.1016/j.engfracmech.2006.02.012. [24] Tvergaard V. Material failure by void growth to coalescence. Adv Appl Mech 1990;27:83–151. [25] Chu CC, Needleman A. Void nucleation effects in biaxially stretched sheets. J. Engng Mater Tech 1980;1028:249–56. [26] Pardoen T, Hutchinson JW. An extended model for void growth and coalescence. J Mech Phys Solids 2000;48:2467–512. [27] Pardoen T, Hutchinson JW. Micromechanics-based model for trends in toughness of ductile metals. Acta Mater 2003;51:133–48. [28] Gordon JR, Wang YY. The effect of weld mis-match on fracture toughness testing and analyses procedures. In: Schwalbe KH, Koc¸ak M, editors. Mis-matching of welds, ESIS 17. London: Mechanical Engineering Publications; 1994. p. 351–68. [29] Kirk MT, Dodds RH. The influence of weld strength mismatch on crack-tip constrain in single edge notch bend specimens. Int J Fract 1993;18:297–316. [30] Toyoda M, Minami F, Ruggieri C, Thaulow C, Hauge M. Fracture property of HAZ-notched weld joint with mechanical mismatching – Part I. Analysis of strength mismatching of welds on fracture initiation resistance of HAZ-notched joint. In: Schwalbe KH, Koc¸ak M, editors. Mis-matching of welds, ESIS 17. London: Mechanical Engineering Publications; 1994. p. 399–415. [31] Hao S, Schwalbe KH, Cornec A. The effect of yield strength mis-match on the fracture analysis of welded joints: slip-line field solutions for pure bending. I. J Solids Struct 2000;37:5385–411. [32] Joch J, Ainsworth RA, Hyde TH. Limit load and J-estimates for idealised problems of deeply cracked welded joints in plane-strain bending. Fatigue Fract Engng Mater Struct 1993;16:1061–79. [33] Burstow MC, Ainsworth RA. Comparation of analytical, numerical and experimental solutions to problems of deeply cracked welded joints in bending. Fatigue Fract Engng Mater Struct 1995;18:221–34. [34] Burstow MC, Howard IC, Ainsworth RA. The influence of constraint on crack tip stress fields in strength mismatched welded joints. J Mech Phys Solids 1998;46:845–72. [35] Burstow MC, Howard IC. Constraint effects on crack growth resistance curves of strength mismatched welded specimens. In: Schwalbe KH, Koc¸ak M, editors. Mis-matching of interfaces and welds. Geesthacht, FRG: GKSS Research Center Publications; 1997. p. 357–69.

I. Pen˜uelas et al. / Engineering Fracture Mechanics 73 (2006) 2756–2773

2773

[36] Kim YJ, Schwalbe KH. Numerical analyses of strength-mismatch effect on local stresses for ideally plastic materials. Engng Fract Mech 2004;71:1177–99. [37] Al-Ani AM, Hancock JW. J-dominance of short cracks in tension and bending. J Mech Phys Solids 1991;39:23–43. [38] Betego´n C, Hancock JW. Two-parameter characterization of elastic–plastic crack-tip fields. J Appl Mech 1991;58:104–10. [39] Bilby BA, Cardew GE, Goldthorpe MR, Howard IC. A finite element investigation of the effect of specimen geometry on the fields of stress and strain at the tips of stationary cracks. In: Size effects in fracture. I Mech E, London 1986. p. 37–46. [40] Du ZZ, Hancock JW. The effect of non-singular stresses on crack tip constraint. J Mech Phys Solids 1991;39:555–67. [41] O’Dowd NP, Shih CF. Family of crack-tip fields characterized by a triaxiality parameter – I. Structure of fields. J Mech Phys Solids 1991;39:989–1015. [42] O’Dowd NP, Shih CF. Family of crack-tip fields characterized by a triaxiality parameter – II. Fracture applications. J Mech Phys Solids 1992;40:939–63. [43] Hutchinson JW. Singular behaviour at the end of a tensile crack in a hardening material. J Mech Phys Solids 1968;16:13–31. [44] Rice JR, Rosengren GF. Plane strain deformation near a crack tip in a power-law hardening material. J Mech Phys Solids 1968;16:1–12. [45] Zhang ZL, Hauge M, Thaulow C. Two parameter characterisation of the near tip stress field for a bi-material elastic-plastic interface crack. Int J Fract 1996;79:65–83. [46] Zhang ZL, Hauge M, Thaulow C. The effect of T stress on the near tip stress field of an elastic–plastic interface crack. In: Advances in fracture research proceedings, ICF9, vol. 5, 1997. p. 2643–50. [47] Zhang ZL, Thaulow C, Hauge M. Effects of crack size and weld metal mismatch on the HAZ cleave fracture toughness of wide plates. Engng Fract Mech 1997;57:653–64. [48] Taulow C, Zhang ZL, Hauge M, Burget W, Memhard D. Constraint effects on crack tip stress fields for cracks located at the fusion line of weldments. Comput Mater Sci 1999;15:275–84. [49] ABAQUS 6.4., Hibbit, Karlsson and Sorensen, Inc., Pawtucket, 2003. [50] Betego´n C, Rodriguez C, Belzunce FJ. Analysis and modelisation of short crack growth by ductile fracture micromechanisms. Fatigue Fract Engng Mater Struct 1995;20:633–44. [51] ASTM E8M-04. Standard test method of tension testing of metallic materials [Metrics]. Book of Standards Volume 03-01. [52] ASTM E 1820-05a. Standard test method for measurement of fracture toughness. Book of Standards Volume 03-01. [53] Landes JD. The blunting line in elastic-plastic fracture. Fatigue Fract Engng Mater Struct 1995;18:1289–97.