A Gurson model in association with a cohesive zone model on the investigation of cleavage-ductile transition in metals

The rupture in metals can occur by cleavage, where all process is controlled by stresses, or by ductile fracture, which takes into account the damage caused by nucleation and growth of voids.The process is then dependent on stresses and strains. The Linear Elastic Fracture Mechanics, widely used in engineering practice, is based on the assumption that the first process prevails, which occurs only under certain conditions. Consideration of the second fracture process is not so well disseminated. In this work, two methodologies are considered to take into account the cleavage-ductile transition. One is based on the Tvergaard-Hutchinson's cohesive model and the other is based on the Gurson-Tvergaard-Needleman's ductile damage model. The two methodologies are considerably different and, in this work, initially the relationships between the two models are established. Then, the conditions for a transition from cleavage to ductile fracture are determined and discussed. Most of the results are presented based on crack growth resistence curves obtained for different material parameters. A strip in mode I rupture is considered firstly. It is shown that, depending on the yield stress and other factors, the two fracture modes can coexist. Also, even when only cleavage is occurring, it is affected by interactions with voids. Lastly, the present simulations are compared withCompact Tension experimental results.Results considering the coupling between the two fracture models presented a better fitting with experiments than other simulations where the coupling is not considered.


INTRODUCTION
Cohesive zone models (CZMs) are extensively explored in the study of cleavage fracture.The concept of this model is based on the principles developed by Dugdale [1] and Barenblatt [2].They considered tractions with opposite modules regarding applied stresses at the crack tip front, known as cohesive zone.The fracture process in the CZM takes place at the crack line, relating tractions to displacements at the coesive zone.Among the existing CZMs, some models developed over the years stand out, such as the constitutive laws formulated by Xu and Needleman [3], Tvergaard and Hutchinson [4], Zhang and Paulino [5], Roth [6], among others.The Tvergaard-Hutchinson [4] trapezoid model will be the basis for this work regarding this matter.In the model, the forces applied on the body cause the opening of the cohesive zone.Tractions increase until reaching a maximum value that will remain constant for a certain opening, and then they will finally decrease when the opening further increases.When tractions approach zero at the final stage, all of the material toughness is then consumed, completing the rupture process.
Ductile failure in continuum damage mechanics is established by softening through degradation of mechanical properties of the material [7].Unlike cleavage rupture, the ductile fracture considers plastic strains and not only stresses generated in the process as governing variables of the process.Among the various formulations explored for material damage, this work focuses on the formulation developed by Gurson [8], which takes into account voids as a degradation variable.Voids nucleate by debonding or fracture of second phase particles, which would grow until later agglomeration, triggering the phenomenon of coalescence.Most of the voids are nucleated in the loading process, even though an initial fraction of porosity may be considered.Tvergaard [9] states that the most significant part of the voids are formed in the early stages of deformation.
Over the years, many modifications have been made to Gurson initial model.Tvergaard & Needleman [9,10] proposed the insertion of new parameters to modify the effect of hydrostatic pressure at all stages of deformations, modeling problems that comprise shear bands in a more precise way.This new formulation became known as the Gurson-Tvergaard-Needleman (GTN) law of continuous damage, and it is the model used in the present work.Other modifications such as the ones made by Zhang & Niemi [11] and Nahshon & Hutchinson [12], were not taken into account in this work.
It is essential to establish that considerably less energy is dissipated during cleavage crack propagation compared to ductile failure [13], so one fracture mode can fundamentally change the other when both are considered together.Over the years, a few studies investigated crack propagation by using a sort of combination between CZM and Gurson's modified model.Tvergaard and Hutchinson [19] studied the reduction of the peak stress in the CZM through a plastic strain effect from the acceleration of void growth and nucleation, with the help of the GTN model.Tvergaard [35] made predictions of ductile fracture by CZM making use of a special interface element.In this element, unlike the initial CZM proposed by Tvergaard and Hutchinson [19] where the initial width of the interface element is zero, the separation law is calculated based on a background element with a non-zero width indicating the porosity of the material, where the GTN is applied.However, this approach would only allow ductile fracture analysis, considering that the separation law would no longer be based on cleavage.Uthaisangusuk et al. [30] analyzed the evolution of voids and cleavage in multiphase steels by RVE calculations, using GTN and CZM simultaneously, finding a good comparison with experimental results.In a more recent study, Hutter et al. [13] proposed the combination of the two rupture models for a complete analysis of the crack propagation.An exponential CZM by Xu and Needleman [31] was implemented in combination with a non-local modified Gurson model [24].They found that this kind of approach captures qualitative effects of corresponding experiments such as the cleavage initiation in front of a stretch zone, the formation of secondary cracks and possible crack arrest.Nevertheless, no analysis has been done associating the interference of one model functionality with the other.
The cleavage-ductile transition will be investigated in this work through numerical analyses of elastic-plastic materials modelled by finite elements.In particular, it will be explored how mechanical and fracture properties as well as inertial properties affect transition between the two fracture modes.This is a subject not fully understood or explored and it is a fundamental issue in the fracture ofmetals.The rupture analysis will be done mostly by the use of the crack growth resistance curves (R-curves) [9][10] [17][19] [28][14].In section 2 the methodology is described.Applications, comparisons and discussion are presented in section 3. Finally, in section 4 conclusions are made.

FAILURE MODELS
In this section the two types of failure models considered in this work, cleavage and ductile fracture, are briefly described as well as computational implementation.Aspects related to definition of the crack size when the ductile model is used are also explored.

Cleavage fracture by Cohesive Zone Model
Cleavage fracture is considered by a cohesive zone model based on the Tvergaard-Hutchinson proposal [4].Traction vs separation relation is shownin Figure 1.λ is a parameter obtained as a function of the cohesive zone separation and λ 1 and λ 2 appear as form factors of the separation law.
In Equation 1, Δ n and Δ t are respectively the normal and tangential openings and δ n and δ t are their respective critical values.Thus, as observed in Figure 1, λ will increase with crack opening, characterizing the rupture when λ=1.Unlike the model disseminated by Needleman [15], this factor will cause tangential displacements also to be accounted for in determining the traction in the cohesive zone.Thus, according to Tvergaard [4] as long as λ increases monotonically, tractions are expressed by: where the potential φ of consumed fracture energy (or distribution of fracture energy) is: G c is the energy consumed per unit area in the cleavage fracture process and can be obtained by the area below the ( ) σ λ function, Figure 1, as:

GTN-Model for ductile fracture
The development of the ductile fracture is modelled using the formulation of Gurson [8], modified by Tvergaard [9] [10], and Tvergaard and Needleman [17], which included the parameters 1 α and 2 α to change the effect of the hydrostatic pressure at all stages of deformation (GTN-Model).
The yield condition in the GTN-Model is: (5) Where f is the current void volume fraction, p the hydrostatic stress, s ij the deviatoric part of the stress tensor, y σ the yield stress and ij δ the Kronecker delta.The void growth rate ( f  ) will depend on the value of a critical void volumetric fraction C f , which will dictate the beginning of the void coalescence: Where n f  is the rate of nucleated voids, g f  is the growth rate and c f  the rate of void coalescence.Their values are expressed by: (1 ) where    is the plastic strain;   is the equivalent plastic strain and ̇ is the rate of it.N f is the volumetric fraction of nucleated voids; N ε is an equivalent deformation value around which the voids are formed and N s is the corresponding standard deviation.
Needleman and Tvergaard [17] have been using the parameters of in their studies for different types of materials.It is observed that these values offer a good range for analysis since nucleation does not occur abruptly.
In Equation 10, C f and ε ∆ are the parameters that control the void coalescence and U f is the void volume fraction corresponding to the rupture in the absence of hydrostatic stress and it is approximated by: Latin American Journal of Solids and Structures, 2021, 18(3), e356 5/18 Coalescence occurs in a volumetric void fraction between 0.1 and 0.2 [18].Tvergaard and Needleman [9][10] [17][20] previously used an C f from 0.15 to 0.2 in numerical analysis, which will be taken into consideration in the first phase of the applications in this work.However, higher values of C f become very critical when considering the cleavage-ductile transition, given that the cleavage failure will dominate the process before the void volume reaches C f .

Numerical implementation
The use of cohesive models in finite element simulations is based on the insertion of interface elements between volumetric elements of the mesh.The opening is controlled by the separation law established in the CZM, as seen in Figure 2.
Where T e is the traction vector and Γ T,e the surface of the cohesive element.In the present work, four Gauss points are used to integrate Equation 13.The term can be introduced in the Principle of Virtual Work, which can be written for the body as follows: : U represent displacements in the volume body Ω; Δ are the crack openings; ρ is the material specific mass; B are the volume forces; F are the forces on surface Γ F ; T are the cohesive tractions and σ are the Cauchy stresses.δ defines a virtual variation.Equation 14 is solved by a Central Difference Method when inertial contributions are taken into consideration (second term, on the left).Otherwise the equation is solver by a Newton-Raphson procedure.
The ductile model is implemented through the volumetric constitutive equations.The increment in Cauchy stress is obtained by: Latin American Journal of Solids and Structures, 2021, 18(3), e356 6/18 where Δσ c is the increment of Cauchy stress tensor in a corotational system, C is he Hooke tensor, c D is the strain rate tensor and ln,c E  the logarithmic strain rate, both calculated also on the corotational system.Total stresses are calculated as: where Δσ can be obtained from Δσ c (Equation 15) by rotation through the rotation matrix of the polar decomposition.
Integration of Equation 15 is done by a predictor-corrector type of algorithm, as in small deformation algorithms.Details can be found for instance in Bittencourt and Creus [32] and Bittencourt [33].Plastic strains must be also obtained from equations described in section 2.2.Details can be found in Cunda [7].
Bi-linear 4-node isoparamentric finite elements are used.4 Gauss points are used to integrate deviatoric part of stresses and 1 Gauss point is used to integrate the volumetric part in Equation 15.Hence, the two fracture models are coupled with the CZM working through the interface elements and the GTN-model acting on the constitutive law of the volumetric elements.

Ductile propagation
In order to analyze the transition between the two processes, it is important to define a relationship for the propagation of damage at the crack tip in the form of propagation of the crack itself ( a ∆ ).Differently from cleavage failure, in the continuous damage model there is not a discrete macroscopic crack.Propagation can be inferred following the porosity parameter in the volumetric finite elements.
In previous studies dealing with cleavage-ductile transition, the use of cells [21] in finite element meshes were widely used, mainly in dynamic analysis [22] [23].However, this kind of methodology has been used less frequently in recent years since the stresses that trigger the cleavage cannot be precisely defined.Also, it is desirable to have a model where the damage propagates continuously and homogeneously in the material, and not through a discrete distribution of cells.
Needleman and Tvergaard [29] studied cases of ductile-cleavage transition, defining a boundary in front of the crack tip where . The length of this region, a ∆ , is added to the crack length.In more recent studies, such as Linse et al.
[24] and Hütter et al. [13], crack growth amount a ∆ was measured as: where W is the total solid length, 0 A the initial crack length, 0 f the initial void volume fraction, u f the ultimate void volume fraction and * f the current void fraction volume.This definition is used in the present applications, considering only regions where * C f f ≥ .

APPLICATION AND RESULTS
Applications, comparisons and discussions are presented in this section.It is divided in two parts.First a strip in pure mode I is analyzed and then a Compact Tension test is simulated and compared to experiments.

Stripcase
Analyzed problem is a strip in plane strain, shown in Fig. 3. Dimensions of the strip are 0.003 vs.0.001 mm.The elastic-plastic solid has an initial tensile yield stress σ y , with an uniaxial behavior defined as follows: Where E is the Young's Modulus, h the strain hardening modulus.An aluminum piece with E=70000 MPa, ν=0.33 (Poisson's ratio) is investigated.Different values of h andσ y were considered.
In the middle of the strip a crack 1.4 x 10 -4 mm long is placed in x direction, as indicate by the red horizontal line in Figure 3. Prescribed nodal displacements of 3 x 10 -5 mm are applied at all the upper nodes of the piece in y direction and prescribed displacement of -3 x 10 -5 mm are applied at all bottom nodes of the piece, in the same direction.Therefore the crack is under pure mode I fracture.In Figure 3 is possible to see the used mesh, as well as the refined part along the fracture line.The mesh has 5920 elements with 6118 nodes, where the size of the finite elements near the crack tip is 2 x 10 -6 mm.Cohesive interface elements are placed only in the plane of the crack.The CZM toughness c G was considered equal to 0.001 N/mm, with a max σ equal to 1750 MPa (except when mentioned otherwise).The values of 1 λ and 2 λ used were respectively 0.15 and 0.5, according to Tvergaard previously studies.Therefore, it is possible to obtain (Equation 4) a normal critical opening n δ equal to 8.4656x10 -7 mm.Below three distinct situations are considered: first pure cleavage; then it is investigatigated the effects of the damage model on cleavage behavior; finally the transition cleavage-ductile is studied.

Pure Cleavage
First, cases with pure cleavage rupture are investigated, or only the CZM is acting in the failure process.In Fig. 4 and Fig. 5, R-curves (normalized applied energy J/G c vs.normalized crack propagation Δa/δn) for distinct values of r (σ max /σ y ) are presented.Different normalized strain hardening modulus h/E are also considered.As noticed, the higher the r value, the greater the material toughness.In other words, for crack propagation to occur, the increase in r will require more energy to create new surfaces due to plastic dissipation.There was practically no gain in toughness for r below 2.5, as also observed by Tvergaard [4].The peak stress required in the fracture process is around 2.5 σ y and then lower stresses will cause the crack to propagate without full development of the plastic zone.
Another important factor is that the crack propagation becomes slower with the reduction of the strain hardening modulus h.On the other hand, the computational analyzes were more stable with the increase of this parameter.

The influence of the GTN properties on the CZM
In the next applications, the ductile model is introduced.The effects of void nucleation followed by void growth are analyzed, allowing to investigate how much this process affects cleavage rupture in the body when the two models are coupled.In Figure 6, a delay was observed in the crack propagation with the beginning of nucleation.This is due to the fact that, with the increase in the void volume, stresses in the solid decrease, with part of the applied fracture energy (J) being consumed in the process.
The nucleation of voids occurred more slowly with the increase of ε N .Such occurrence takes place because ε N causes the nucleation law to change or the higher its value, the slower the nucleation phenomenon.It was possible to obtain the variation of normalized stressesσ 22 /σ y and void volume fraction f at the crack tip during growth and after coalescence (Figure 7 and Figure 8).The figures are a function of the equivalent plastic deformation   .For a given stress, a greater equivalent plastic strain   was found with the reduction of the initial impurity/porosity level, implying that despite the cleavage fracture being controlled only by stresses, it is indirectly affected by these aspects.In previous studies, Tvergaard and Hutchinson [19] reported that a decrease in the maximum traction σ max of the cohesive law should be taken when the effective plastic strain exceeds a critical value.This would be necessary because the cohesive zone model does not account for the reduction in material resistance due to plastic deformation.This is an important issue, mainly because σ max and G c are the most significant parameters in the CZM.In the present work this correction of σ max value is not done.

Transition between cleavage-ductile modes
The transition between fracture modes will be investigated below.For the examples that follow, a symmetric version of the previous mesh is used.This means that only the part above the crack line is considered, being all nodes at the crack line being blocked in y direction, with exception of the nodes of the initial crack surface.The GTN-Model parameters used in the examples are: f U =0.6, f N =0.04, ε N =0.4,s N =0.17, h=E/100, r=5, f 0 =0.001,Δε = 0.6, f C =0.15 α 1 =1.5 and α 2 =1.0.
For low values of initial yield stress, r=5, the transition occurred more quickly, with the propagation of damage accumulating in front of the crack tip and generating a crack blunting, as seen in Figure 9.The increase in the initial yield stress will delay the propagation of the damage, making the transition to take place later.For a high enough yield stress, the region where the damage is formed cannot keep up with the crack propagation by cleavage, resulting in a low density of voids in front of the crack tip.It will cause the solid to fail predominantly by cleavage, even if the material damage slightly changes the speed of crack propagation.This event can be seen in Figure 10, for r=3.5.However, still for low r values, in the course of propagation the hydrostatic stress at the crack tip may reach high enough values to trigger ductile fracture, even if there is nota considerable density of voids in the region.Then the ductile fracture is no longer triggered by the void evolution and starts to be controlled by hydrostatic stress.This aspect is shown in Fig. 10, where even for low f values (a), ductile rupture still takes place due to high values of pressure (b).It was also found that during the transition, in some cases, the two rupture processes occurred simultaneously for a brief period.The tendency of this event and its duration was increased by the beginning of coalescence in the material.
The numerical solution for a few cases stopped due to convergence problems during the ductile regime.The numerical problems are related to elements close to the crack tip that are near the failure and are being excessively deformed by the "weakening" effect of the damage process in the GTN-Model.The same kind of problem was reported by Linse et al. [24] and by Samal et al. [25] in highly refined meshes.
The reduction of σ max in the CZM caused the crack to propagate more quickly, delaying the transition to the ductile process, as seen in Figures 12 and 13, for the same r.This is an interesting aspect, indicating that the ratio r alone does not uniquely define propagation.However, ratior tends to be the dominant factor in the transition.For instance for r values below 2.68, 2.62 and 2.52 for respectively σ max equal to 1750 MPa, 1400 MPa and 1000 MPa, propagation always occurs only by cleavage.Conversely, propagation always occurs only by ductile fracture for r>5.5, approximately.
Latin American Journal of Solids and Structures, 2021, 18(3), e356 13/18 The R-curves for cases where inertial terms are considered are shown in Figure 15.The impact velocity on the following dynamic examples was 3000 mm/s.Transitions occurred in different energy ranges, compared to quasi-static cases.Therein, r values less than 2.2 cause the crack to propagate exclusively by cleavage, and for r values above 4.3, the rupture was completely ductile.The brittle-ductile transition occurred for r values inside this interval.Figure 15 shows a comparin between dynamic and quasi-static cases for few r combinations.In all cases, ductile and cleavage propagation was delayed by the inertial effect.The effect is expected, see for instance Freund [34], as stresses increase more slowly due to inertial effects.It is important noticing also that the Latin American Journal of Solids and Structures, 2021, 18(3), e356 14/18 transitions occurred at the same values of Δa, for each r, regardless of inertial contributions.The curves do not differ initially, with a tendency for higher energy dissipation in the dynamic case.

C(T) case -comparison with experiments
Below, simulations of a C(T) (Compact Tension) specimen are performed and compared with experimental and numerical data obtained by the studies in the references [26][27] [28].
The dimensions of the solid used for the test are: W=50 mm, a 0 /W=0.59where W is the full length of the piece, and a 0 the initial crack length, see Figure 16.The material is a ferric structural steel with the German classification StE460.It has the following properties: E=210000 MPa, ν=0.3, σ y =470 MPa.The hardening properties used by the authors in the experiments were not specified, so simulations in the present work are made for cases with and without strain hardening in order to find a better fit to the experiments.
The CZM parameters followed the same data given by the references [26][27] [28]: 567, where D is a scale factor [27], and it is widely used as a definition of the characteristic size L of finite elements in numerical analyzes, with L min =D/2.Therefore, D=0.2 mm and L min =0.1 mm.The characteristic cohesive lengthδ n , Equation 4, is then 0.05 mm.
The GTN model properties determined by Brocks et al. [26] are: f U =0.19, f N =0.2, ε N =0.3, s N =0.1, f 0 =0.0025, f C =0.021.The finite element mesh built for this simulation is shown in Figure 16.A prescribed displacements u LL is applied at the loading point, Figure 16.In Figure 17 the corresponding force at the loading point versus u LL is shown, together with results from the literature [26,27].Considering the present work, two different simulations were performed, as seen in Figure 17.One considering an elastic-perfect plastic material, and a second using a linear strain hardening of h=E/100.The simulation without hardening presented a sudden drop in the force soon after the beginning of unstable crack propagation.This same behavior is found in studies by Hütter [13].The case with hardening presented the best fitting with experimental behavior.The fitting is considerably better than results presented in Sigmund and Brooks [27].It must be emphasized here that in the present simulations the cohesive model and the GTN model are both considered.In simulations done in Sigmund and Brocks [27], only the GTN model was considered.Cleavage propagation started for u LL =0.394 mm.As the cohesive zone is activated before propagation begins, the difference between Sigmund and Brooks [27] and the present simulations can be explained by the consideration of the cleavage process.This application demonstrates the practical importance of coupling both fracture processes.
Figure 18 shows the R-curves obtained from simulations with linear strain hardening and J-integral [16] calculated according to ASTM-E1820-11, Equation 19: where K is the stress intensity factor, A pl the area under the force (F) vs. displacement (u LL ) curve, B N the thickness, b 0 the uncracked ligament (W-a 0 ), νis the Poisson's ratio and η pl a parameter given by η pl =2+0.522b0 /W.

Figure 2
Figure 2 Illustration of the cohesive element in the initial state and after loading.Cohesive zone openings are stored in a vector as shown below:

Figure 3
Figure 3 Finite element mesh.Refined part along the crack line.

Figure 4 R
Figure 4 R-curves for pure cleavage failure with h=E/300.

Figure 5 R
Figure 5 R-curves for pure cleavage failure with h=E/100 and an elastic-perfect plastic solid.

Figure 6 R
Figure 6 R-curve for CZM tests considering nucleation and growth of voids.(a) Elastic-perfect plastic material.(b) h=E/100.
Growth and coalescence effects: In this second part, growth and coalescence are also considered.The parameters used are f N =0.05, ε N =0.3, s N =0.1, h=E/100 and r=5.Three different values of impurities/initial porosity f 0 =0.01, 0.005, 0.001, two values of Δε= 0.25, 0.6 and also two values of f C =0.1, 0.15 are considered.

Figure 9
Figure 9 Propagation of damage near the crack tip for r=5.Elements wheref >f C are ruptured.

Figure 10
Figure 10 Distribution of damage (a) and Pressure field (b) near crack tip for r = 3.5.Elements with p>1000 MPa are ruptured by ductile mode.

Figures 11 to 15
show R-curves depicting where the transition occurs.In Figures 12 and 13 different σ max are also considered.In Figures 14 and 15 inertial terms are added.In general, initially the curves follow the behavior observed in the pure cleavage examples, where the greater the r value, the slower the cleavage fracture propagation.Also for greater r values, ductile rupture tends to dominate the propagation process, as expected.

Figure 11 R
Figure 11 R-curves for quasi-static tests with different values of r.

Figure 12 R
Figure 12 R-curves for a quasi-static tests with different values of r andσ max =1400 MPa.

Figure 13 R
Figure 13 R-curves for a quasi-static tests with different values of r andσ max =1000 MPa.

Figure 14 R
Figure 14 R-curves for dynamic tests with different values of r .

Figure 15
Figure 15 Comparison of R-curves for dynamic and quasi-static tests.

Figure 16
Figure 16 Symmetric mesh used for simulation of the C(T) test -(a) Full mesh.(b) Refined mesh around the crack tip.Dimensions in mm.

Figure 17
Figure 17 Force versus load-line displacement

Figure 18 R
Figure 18 R-curves for the CT specimen.