A Multi-Scale Numerical Method for the Study of Size-Scale Effects in Ductile Fracture

The use of a stress-strain constitutive relation for the undamaged material and a traction-separation cohesive crack model with softening for cracking has been demonstrated to be an effective strategy to predict and explain the size-scale effects on the mechanical response of quasi-brittle materials. In metals, where ductile fracture takes place, the situation is more complex due to the interplay between plasticity and fracture. In the present study, we propose a multi-scale numerical method where the shape of a global constitutive relation used at the macro-scale, the so-called hardening cohesive zone model, can be deduced from meso-scale numerical simulations of polycrystalline metals in tension. The shape of this constitutive relation, characterized by an almost linear initial branch followed by a plastic plateau with hardening and finally by softening, is in fact the result of the interplay between two basic forms of nonlinearities: elasto-plasticity inside the grains and classic cohesive cracking for the grain boundaries.


Introduction
The cohesive zone model (CZM) is widely applied in the computational mechanics community to model cracking in quasi-brittle materials at the macro-scale [1][2][3][4].Following this approach, the heterogeneous material composition is replaced by a homogenized continuum and the resistance to crack propagation is simulated by the action of cohesive tractions acting in the process zone along the crack path.At the same time, the CZM can be used at the meso-scale to depict the phenomenon of intergranular separation [5,6] in polycrystalline materials.In this instance, the nonlinearity is concentrated in the constitutive relation of the interfaces along the grain boundaries, whereas the grains are supposed to behave linear elastically.
In the majority of metals, on the other hand, large plastic deformation with significant elongation of the grains is observed, as, e.g., during metal forming.Even in this case, characterized by a high ductility, however, the load-displacement curve presents softening after the plastic regime, if the test is carried out up to very large displacements (see Figure 1).This softening, impossible to be explained as the result of purely plastic deformations, is the physical result of the formation of small micro-cracks and micro-cavities that coalesce, grow and finally lead to specimen failure.Figure 1.Typical example of micro-cracks coalescence in metals leading to the softening branch and to the final failure of a specimen tested in uniaxial tension (taken from [7,8]).
Another interesting problem where there is a strong interplay between plasticity of the continuum and grain boundary decohesion is represented by the MgCa0.8alloy investigated in [9].This brittle Magnesium alloy used for medical implants contains Calcium as alloying element.Besides the grain refinement, which is of minor importance since the material fully recrystallizes during extrusion, Calcium increases the strength and reduces the tendency to corrosion.Though the alloying content of Calcium is less than the maximal solubility of Calcium in Magnesium, degenerate eutectics formed between the grains.These very brittle eutectics break into little fragments during subsequent extrusion.The presence of these interfaces along grain boundaries promotes brittle intergranular cracking in addition to ductile plastic deformation of the grains.
These experimental evidences suggest that the interplay between the two elementary forms of nonlinearity related to grain boundary decohesion and plasticity inside the grains is responsible for the failure behavior of the material at the meso-scale.In this article, it will be shown that the simultaneous use of the CZM and plasticity in a meso-scale model of a polycrystalline metal leads to size-scale dependent homogenized stress-relative displacement relations.These relations can be used at the macro-scale, in an upscaling procedure towards a multi-scale analysis of the material response.Hence, interface constitutive laws for the simulation of crack propagation in ductile materials at the scale of the structural element/specimen can be defined based on the results established at the meso-scale.At the macro-scale, a numerical algorithm and an analytical limit analysis approach are proposed and used to simulate bending tests.The role played by the initial notch length and the size-scale effects on the mechanical response in metallic materials are investigated.Numerical predictions are also compared with experimental results, showing a remarkable good agreement.

Constitutive Relations for the Grain Boundaries and for the Grains
Grain boundaries are sometimes considered as perfectly bonded interfaces.However, these material discontinuities are often the source of damage and cracking in polycrystalline materials.Hence, imperfect bonding is a more realistic assumption [10,11].Following the 2D implementation of interface elements with a nonlocal CZM discussed in [10], a generalization to 3D FE topologies has been carried out in [9], where the analysis of a complex 3D set of grains and interfaces has been proposed.To this purpose, new user-defined subroutines for linear and quadratic interface elements compatible with brick and tetrahedra continuum elements have been implemented in the Finite Element Analysis Program FEAP [12].
The basic operations of an interface element are described in details in [10,[13][14][15].Essentially, this is a finite element with zero-thickness that is used to impose tractions on the grain boundaries dependent on the relative displacements of the grains sharing the interface.In the kinematic part of the element, the relative opening and sliding displacements of the opposing nodes of a common interface are computed.The cohesive tractions are determined in terms of these relative displacements according to a specified CZM relation, see also [1,3,4].In the present study, the CZM by Tvergaard [16] is adopted, since it is frequently used to model crack propagation in the microstructure of metals.The normal and the tangential tractions, σ and τ, are computed as functions of the normal and tangential relative displacements g N and g T in the local frame defined by the normal unit vector to the average plane of the interface element, n, and by the tangent vectors t 1 and t 2 : T Tc τ γσ (λ) where The normal and the resultant tangential tractions are shown in Figure 2 vs. g N , and g T In this formulation, l Nc and l Tc denote the critical separations in the normal and tangential directions corresponding to complete decohesion.The parameter σ p defines the Mode I peak CZM traction, whereas the Mode II peak CZM traction is τ p = γσ p .The variable λ physically represents a measure of the total relative displacement experienced by the interface.
Due to the nonlinearity of the CZM, the numerical simulations are carried out in an implicit FE code.Hence, a consistent linearization of the weak form of the finite element is performed to compute its tangent stiffness matrix and the residual vector, to be used within a Newton-Raphson iterative procedure.Regarding the constitutive relation for the bulk, a standard von Mises elasto-plastic formulation is considered.

Interplay between Intergranular Fracture and Grain Plasticity and Derivation of Homogenized Relations for the Upscaling Procedure
The basic interplay between intergranular cohesive fracture and grain plasticity is herein investigated in reference to a system consisting of two blocks joined by a cohesive interface and subjected to a uniaxial tensile test.The lower side is restrained to the displacements in the axial direction, whereas the upper side is subjected to an imposed vertical displacement δ.The two blocks are discretized with quadratic brick elements and quadratic interface elements, see Figure 3b.The uniaxial test is carried out under displacement control.An elasto-plastic constitutive model with hardening (E = 152 GPa, ν = 0.3, σ y = 80 MPa, k H = 120 MPa) is considered for the bulk (see Figure 3a).Regarding the nonlinear CZM at the interface, the formulation in [16] is used with a fracture energy G F = 0.1 N/mm.Different values of σ p are explored, see Figure 3a.For σ p ≤ σ y , the behavior of a brittle polycrystalline material is simulated and fracture is induced by intergranular cracking.The nonlinear response is therefore ruled by the CZM, see Figure 4a.
For σ p > σ y , the response for σ < σ y is due to the CZM initial stiffness, which is not infinite for numerical reasons, and the elasticity of the system (see Figure 4b).Therefore, the CZM contribution to the overall deformation before the peak load should not be interpreted as a real grain boundary separation.For higher stresses, plasticity takes place and modifies the global response (see Figure 4b).The interesting aspect is the appearance of a softening branch after hardening, which is now related to a real grain boundary separation.Therefore, numerical results shown in Figure 4 pinpoint the strong interplay taking place between grain boundary decohesion and grain plasticity, depending on the material parameters.These results are confirmed by numerical simulations on polycrystalline materials, see Figure 5.It is interesting to point out that the numerical curves obtained according to the interplay between CZM and plasticity can be interpreted as a new CZM by plotting the homogenized stresses vs. an axial displacement obtained by removing the elastic component from the imposed total displacement δ.An example of this procedure is shown in Figure 6 for one of the cases analyzed in Figure 4.The initial part of the obtained traction-separation law depends on the shape of the CZM used for the interface.For σ > σ y , on the other hand, hardening takes place due to plasticity.Finally, an almost linear softening of the cohesive traction is observed for large separations, again due to the CZM formulation.By increasing the sample size and assuming that σ p is size-scale independent, the traction-separation law is not scale invariant, although the fracture energy, the peak cohesive traction and the critical Imposed axial displacement,  Restrained boundary nodes opening displacement are the same in all the simulations.Here, the size variation refers to the geometrical change of the size of the specimen, keeping constant the grain size and the material properties.This size-scale effect takes place due to the interplay between plastic deformation, which depends on the number of grains in the plastic regime, and cohesive cracking, which depends on the number of grain boundaries involved in the crack pattern at failure.In general, the maximum separation corresponding to the onset of linear softening increases with the sample size, as also observed in experimental tests [17,18].

Numerical Algorithm
A numerical algorithm for a macro-scale analysis is herein proposed with reference to the three-point-bending test (TPB) geometry shown in Figure 7.In this approach, both fracturing and plastic compressive phenomena are assumed to be fully localized along the mid-span cross-section of the beam, whereas the two symmetric beam portions behave elastically.The mechanical behavior of the symmetry cross-section is described by means of the previous meso-scale constitutive model discussed in Section 2. The meso-scale constitutive model was obtained in pure uniaxial tension (Mode I) on specimens where the grain microstructure was explicitly modeled.It has been shown that such meso-scale simulations lead to stress-inelastic displacement curves whose shapes are consistent with those used in [19] for the hardening cohesive/overlapping zone model.Therefore, the meso-scale results can be used as an equivalent cohesive zone model for characterizing the constitutive behavior of the single dominant crack at the macro-scale, where the material microstructure is not taken into account as a material heterogeneity.Since the problems at the macro-scale are characterized by pure Mode I deformation, the meso-scale constitutive model can be applied to all the finite element nodes used in the macro-scale computational model.More complex loading scenarios leading to Mixed Mode deformation would require additional testing at the meso-scale for different Mixed Mode loading, to be consistent in the upscaling procedure.The same constitutive relation is adopted for both tension and compression.From a kinematical point of view, in fact, when the elastic range is overcome in compression, a fictitious interpenetration is assumed to take place, according to the overlapping crack model proposed in [19,20], and similar to the fictitious crack opening considered in tension.More in details, a linear-elastic stress vs. strain relation is assumed for the undamaged material, whereas a hardening stress vs. displacement law, followed by a linear softening tail, is adopted for the damaged and/or yielded material (see Figure 8).The latter law is consistent with the numerical results obtained in Figure 6.For practical purposes, the hardening branch of the cohesive law can also be well described by a Ramberg-Osgood-like relationship: where μ > 0 is the hardening exponent (μ = 0 corresponding to perfect plasticity), σ y is the yielding stress; w is the opening or overlapping relative displacement; w r corresponds to the peak stress and determines the shape of the stress vs. separation law.Beyond the ultimate displacement, the stresses decrease linearly and vanish in correspondence with the critical value for the opening or overlapping relative displacement, w cr .For the application of the numerical scheme, the symmetry cross-section of the specimen is subdivided into finite elements by n nodes.Consequently, cohesive and overlapping stresses are replaced by equivalent nodal forces by integrating the corresponding distributed tractions over the element side.Such nodal forces depend on the nodal opening or overlapping displacements according to the cohesive or overlapping laws shown in Figure 8b.The nodal forces, F, acting along the mid-span cross-section are computed as follows: where {F} is the vector of nodal forces, [K w ] is the matrix of the coefficients of influence for the nodal displacements, {w} is the vector of nodal displacements, and {K P } is the vector of the coefficients of influence for the applied load P. All the coefficients of influence are computed a priori with a linear-elastic finite element analysis.In addition, [K w ] and {K P } are independent of the specimen size, mesh and slenderness ratio being the same.When fracturing and plastic yielding compression take place, the following equations can be considered (see Figure 9, where the central portion of the TPB specimen is shown): The same hardening law is used for tension, Equation ( 8), and compression, Equation (10).Equations ( 6)-( 10) constitute an algebraic system of (2n) equations with (2n + 1) unknowns, i.e., the elements of the vectors {w} and {F} and the applied load, P. In this case, there are two alternative possible additional equations needed for solving the system: we can set either the force in the fictitious crack tip, m, equal to the tensile yielding force, or the force in the fictitious overlapping zone tip, p, equal to the compressive yielding force.In the numerical scheme, we select the situation that is closer to the critical condition.The driving parameter of the process is the tip that in the considered step has reached the limit resistance.Only this tip is moved when passing to the next step.The two fictitious tips advance until they converge to the same node.Finally, at each step of the algorithm it is possible to calculate the deflection, δ, as follows: where {D w } is the vector of the coefficients of influence for the nodal displacements, and D P is the coefficient of influence for the applied load.In addition, {D w } and D P are also independent of the specimen size, mesh and beam slenderness being kept constant.It is worth noting that the same numerical algorithm can be profitably used to study different specimen geometries and loading conditions, as, for instance, the compact tension test (CT).The only difference regards the elastic coefficients entering Equations ( 6) and (11).The geometry and the discretized model used to compute the elastic coefficients are shown in Figure 10.

Limit Analysis for the Three-Point-Bending Test
The numerical simulations carried out with the algorithm proposed in the Section 3.1 terminates when the two fictitious tips converge to the same node.On the other hand, for an exhaustive analysis, the entire post-peak load vs. deflection curve should be determined.To this purpose, a limit analysis approach based on the proposed constitutive laws can provide asymptotic approximation to the post-peak response, which can be relevant in the analysis of ductile-to-brittle transitions.
Once the tensile yielding strength and the compressive yielding strength have been achieved at the bottom and the top beam edges, respectively, fracture and plastic compressive processes are assumed to initiate.The supposed limit situation is shown in Figure 11a.The limit stage of the deformation and fracture process may be considered as that of two rigid parts connected by the hinge A in the mid-span cross-section.The equilibrium of each part is ensured by the external load, the support reaction, the closing cohesive forces and the opening overlapping forces, as evidenced in Figure 11b.For the sake of simplicity, a linear-hardening with subsequent vertical drop relationship is assumed for the stress vs. displacement law (dashed curve in Figure 8b) instead of the power-law hardening expressed in Equation ( 5) followed by the linear softening branch (solid curve in Figure 8b).Several investigations focusing on the effect of the shape of the traction-separation law on the resulting fracture behavior (see [21,22] for metals and [23] for concrete-like materials), in fact, came to the conclusion that such an effect is relatively weak.On the contrary, the cohesive/overlapping dissipated energy plays a fundamental role in determining the overall response.Therefore, a simplified law can be assumed for the limit analysis computations, provided that the area beneath the stress vs. displacement curve is equal to that of the more complex law assumed for the corresponding numerical simulations.The same mechanical parameters (cohesive/overlapping energy, yielding stress, ultimate stress) are adopted for both tension and compression.The geometrical similitude of the triangles ABC and A'B'C' in Figure 11a provides: which gives: The rotational equilibrium around point A is possible for any beam configuration only if the moment of external force and support reaction is equal to the moment of cohesive forces and overlapping forces (see Figure 11b): which gives: According to Figure 11a, the total extension of the cohesive and overlapping zones must be less than or equal to the initial ligament of the beam, so that the limit condition is: where ξ 0 is the ratio between the initial crack length, a 0 , and the beam depth, W. From Equation ( 12), we determine the lower bound to limit analysis in terms of deflection: Similarly, a closed form solution based on the limit analysis approach can be obtained also for the compact tension test scheme.More details are given in [19].

Three-Point-Bending Test of HY80 Steel
In this section, the experimental load vs. displacement curves from Zhu and Joyce [24] are compared with our numerical simulations.The material used for tests was HY80 steel.The mechanical properties, derived from tensile tests on cylindrical specimens with a 25 mm initial gauge, are: σ y = 630 MPa, σ u = 735 MPa, hardening exponent 0.1, and Young's modulus E = 207 GPa.
The above mechanical parameters have been employed in the numerical simulations for the cohesive/overlapping laws.The cohesive/overlapping energy is assumed to be 750 N/mm, w cr = 1.81 mm and w r = 0.3w cr .It is worth noting that the ratio w r /w cr has only a slight influence on the pre-peak behavior whereas the mechanical response, especially the post-peak performance, is significantly determined by the value of the cohesive energy [19].
The comparison between experimental results and numerical simulations is shown in Figure 12, where a very good approximation is obtained for all the considered initial crack lengths.The numerical simulations are terminated when the two fictitious tips converge to the same node, so that the tails of the load vs. deflection curves are described just by means of limit analysis.In addition, the numerical load-deflection curve for a 0 /W = 0.00 is also plotted in Figure 12, which demonstrates a very good continuity between numerical algorithm and limit analysis.The points A, B, C, D, E and F indicate the lower bounds of limit analysis for the different initial crack ratios, ξ 0 = a 0 /W.It may be presumed that all the experimental load-deflection curves asymptotically tend to the limit analysis curve as shown in Figure 12.Comparison between numerical simulations and experimental results [24] for TPB tests with different initial crack lengths (the points A, B, C, D, E, and F represent the lower bound of limit analysis for different initial crack lengths).

Compact Tension Test of DIN 22NiMoCr37 Steel
The material for the compact tension tests (CT) herein considered was obtained from a large forged, quenched and tempered ring segment of ferritic steel DIN 22NiMoCr37.All the experimental data are taken from the Euro-Fracture Dataset [17].Three different specimens scaled in the ratio 1:2:4 have been considered, as described in Table 1.According to the CZM derived in Section 2, size-dependent parameters should be assumed for the constitutive laws.The tensile properties were determined from cylindrical specimens having a diameter of 6 mm and an initial gauge length of 25 mm.They are: σ y = 470 MPa, σ u = 616 MPa, and E = 200 GPa.In the numerical and analytical simulations such values are assigned to the 1T CT specimen, which has the dimensions closer to those of the tensile sample.For such a specimen, the cohesive/overlapping energy is assumed equal to 3200 N/mm, and the Ramberg-Osgood hardening exponent to 0.12.Then, scale-dependent cohesive parameters are assumed for the other specimens.In addition to the size effects on the fracture energy, justified by the meso-scale model presented in Section 2, scale-dependent yielding and ultimate strengths have also been considered, consistently with the experimental investigations on the size effects in mild steel reported in [18].In particular, a variation for the yielding strength with a slope of −0.15 in the bi-logarithmic diagram σ y vs. W, is introduced.Such a scaling law determines σ y = 425 and 385 MPa for 2T and 4T CT specimens, respectively.As regards the scale effect on the cohesive/overlapping energy, it is determined by repeated solutions until the numerical and limit analysis predictions best fit the experimental results (Figure 13).In addition, in these cases, the numerical simulations are terminated when the two fictitious tips converge to the same node, and the asymptotic tails of the load vs. deflection curves are obtained from the limit analysis approach.The best fitting of the experimental results is obtained by assuming the values G F = 4100 and 6400 N/mm for specimens 2T and 4T CT, respectively.The parameter w cr is determined from σ y and G F once the shape of the cohesive law and the ratio w r /w cr are fixed.The cohesive/overlapping constitutive laws adopted for the three considered specimens are shown in Figure 14.In the investigated compact test examples, the material is very ductile as confirmed by the experimental force-displacement curves shown in Figure 13.Consistently, numerical results predict the propagation of the fictitious crack tip only, while the real crack tip does not propagate in the considered load range.Since a steady-state crack growth condition was not achieved, we cannot argue about the behaviour of energy dissipation and possible size-scale effects in that stage.However, we explored another simulation where the fracture energy was reduced in order to have both fictitious and real crack growths.For this case study, the diagram in Figure 15 shows that the dissipated energy rate, defined as the dissipated energy per unit of fictitious crack surface, is an increasing function of the fictitious crack tip position until the real crack tip starts propagating.Later on, the dissipated energy rate remains constant because the steady-state condition of crack growth is suddenly achieved.This particular result cannot exclude other cases when the size of the process zone changes during real crack tip propagation, leading to a further increase in the dissipated energy rate.Therefore, our proposed model can be used to predict the evolution of the energy dissipation also before the onset of real crack propagation mainly due to plastic deformation.

Conclusions
The hardening cohesive crack model is able to predict size-scale effects in ductile fracture taking place in metals, as firstly shown in [19] and further demonstrated with additional examples in the previous section.In the present study it has been shown that the shape of the hardening cohesive crack model can simply be deduced by the fundamental interplay between the following forms of nonlinearity at the meso-scale: grain plasticity and intergranular cracking in the polycrystalline microstructure of metals.Although a detailed calibration of the meso-scale parameters has not yet been proposed, since all the numerical results at the meso-scale refer to standard elasto-plastic constitutive parameters cohesive zone model properties of steel, the obtained stress-inelastic displacement curves are found to be in line with those assumed in [19].Therefore, this opens the possibility to relate the shape of the hardening cohesive/overlapping crack model postulated in [19] to standard mechanical parameters governing the behaviour of the material at the meso-scale.A detailed parameters identification is also ongoing and we prospect a strategy based on performing tensile tests on metallic samples inside the chamber of a scanning electron microscope.The observation of the evolution of cracks and plasticity in the material will permit to identify meso-scale material parameters at the meso-scale, instead of considering them as free parameters to fit the macroscopic response.Therefore, a multi-scale computational method where the shape of the hardening cohesive crack model (equivalent constitutive law) is computed from meso-scale finite element simulations is an effective approach to be pursued.

Figure 2 .
Figure 2. (a) Local frame introduced at the interface element level; (b) and (c) Mixed-Mode CZM tractions used in meso-scale simulations.

Figure 3 .
Figure 3. (a) Constitutive relations for the bulk and the interface and (b) mesh of the 3D uniaxial tensile test.

Figure 6 .
Figure 6.Hardening cohesive laws obtained from the results in Figure 4, for different specimen size-scales (the inelastic displacement w is measured in mm).

Figure 7 .
Figure 7. Finite element nodes along the mid-span cross-section of the TPB specimen.

Figure 8 .
Figure 8. Constitutive laws according to the cohesive/overlapping crack model: (a) linear-elastic stress vs. strain relationship for the undamaged material, and (b) hardening stress vs. displacement law for the damaged and/or yielded material (followed by a linear softening tail).

Figure 9 .
Figure 9. Nonlinear cohesive and overlapping stress distributions along the mid-span cross-section of the TPB specimen.

Figure 11 .
Figure 11.Limit analysis for TPB test: (a) limit condition of fracture with cohesive and overlapping forces; (b) cohesive and overlapping force distributions along the symmetry cross-section.

Figure 12 .
Figure 12.Comparison between numerical simulations and experimental results [24] for TPB tests with different initial crack lengths (the points A, B, C, D, E, and F represent the lower bound of limit analysis for different initial crack lengths).

Figure 15 .
Figure 15.Dissipated energy per unit of fictitious crack surface versus the relative crack tip position.