The study of crack propagation in industrial components plays a critical role in ensuring their safety and reliability, particularly in high-stakes sectors such as aerospace1–4. In this field, where mechanical loads and performance requirements are exceptionally demanding, accurate modeling of defect initiation and growth is essential to prevent catastrophic failures.
To address this challenge, various numerical and theoretical frameworks have been developed. The so-called extended finite element method (XFEM) allows for the representation of crack discontinuities without remeshing5, by introducing localized enrichment functions. This approach proves effective in capturing local singularities and the branching of complex cracks. In parallel, energy-based criteria relying on the J and G integrals, originating from the formalism of linear elastic fracture mechanics (LEFM), quantify respectively the energy release rate and the material’s toughness under quasi-static elastic conditions6. These methods are primarily grounded in LEFM theory, which is well suited for brittle materials or small-scale defects. However, when dealing with elasto-plastic materials, the limitations of these models become readily apparent.
In fact, plasticity induces significant energy dissipation and leads to the development of residual stresses, thereby complicating the analysis and prediction of crack propagation. One of the major challenges lies in accurately predicting the crack growth direction, especially in dissipative environments.
To address these issues, the G-theta method offers an effective solution. This approach combines the computation of an energy quantity characteristic of fracture in dissipative media with an accurate prediction of the crack propagation path.
In this context, the numerical tool Z-cracks, developed through a collaboration between ONERA7 and TRANSVALOR and implemented within the Z-set finite element simulation suite, provides a generic and efficient framework for the study of 3D cracks. This module enables the computation of stress intensity factors (SIFs) for static and mixed-mode crack configurations, as well as the simulation of crack propagation, while integrating robust 3D adaptive remeshing techniques and efficient multi-threaded post-processing for SIF extraction.
In this paper, we investigate the effectiveness of the G-theta method integrated within Z-cracks. By applying this method to various cracking scenarios, we analyze its convergence properties, particularly in an elastoplastic context.
To numerically solve a fracture problem, the Z-cracks procedure involves several steps. The most challenging step consists of generating a conforming mesh around the crack during its propagation, which is a complex task due to the singularity of the stress and strain fields near the crack tip. First, a crack with a predefined geometric shape is inserted into an existing intact mesh. The intersection between this mesh and the crack mesh is then generated. The mesh of this intermediate structure is subsequently adapted to achieve a discretization of satisfactory quality, which is essential for accurately computing the stress intensity factors (SIFs) at the crack front. During this remeshing step, Z-cracks relies on external automatic adaptive remeshing tools (e.g. MeshGems), providing significant flexibility for finite element computations.
The prediction of crack propagation relies on the accurate evaluation of fracture-related quantities, particularly the stress intensity factors and the energy release rate G. A thermo-mechanical finite element analysis is conducted to estimate the stress state within the structure under a given loading. Based on this stress field, the SIFs are computed at multiple points along the crack front in order to determine the propagation angle, the advance increment, and the new position of each node along the front.
The advancement of the crack front is then governed by propagation laws (e.g. Paris-Elber8,9), which allow estimating the geometric evolution of the crack within the material. This evolution requires a continuous update of the finite element model’s topology.
The Z-cracks module automates this process by extracting fracture quantities and the bifurcation angle using the G-theta method. This method, proposed by Destuynder et al.10, relies on domain integrals and on the Lagrangian derivative of the potential energy with respect to a virtual crack extension velocity field, namely the θ field. We thus consider a reference domain Ω containing a crack with a front Γ 0. Let M be a point on the crack front and S the curvilinear abscissa along the front. The transformation F η of all points M in the domain Ω into a domain Ωη, due to crack propagation, is given by:
where η is a constant (η = 0 corresponds to the reference configuration) and θ(M) is the virtual crack displacement field.
Figure 1: Representation of a Crack within a volumetric domain Ω.
The virtual crack extension field θ is assumed to be tangent to the crack surface at the crack front, as well as to the boundaries of the body. This means that the transformations F η modify the position of the crack front Γ 0 without altering the body boundaries.
The energy release rate G(θ) for a virtual crack extension θ is given by the Lagrangian derivative of the potential energy W:
where εth denotes the thermal strain, and ε is the total strain.
where εae denotes the inelastic strain.
By integrating along the crack front Γ 0, with the front discretized using np control points associated with shape functions Nj(s), the energy release rate can be determined by solving the following system of equations with np unknowns:
Validation of the G-theta formulation is crucial to ensure that the modeling and evaluation of cracks in an elastoplastic context are reliable and consistent with material physics. To this end, the values of the energy release rate G were compared with those obtained for a purely elastic material. This comparison allows quantifying the contribution of inelastic residual stresses (thermal and viscoplastic) to the energy release in a realistic system where irreversible deformations play a significant role.
To strengthen the credibility of the approach, a validation phase was conducted using analytical and semi-analytical reference solutions in both 2D and 3D. These validations confirmed the method’s ability to reproduce the expected physical mechanisms in controlled configurations. This phase was then complemented by a numerical convergence analysis, evaluating the effects of mesh refinement and resolution parameters on the stability and accuracy of the results.
A sensitivity study on the calculation of the G value (see equation (2)) was conducted on a two-dimensional notched specimen subjected to uniaxial tension (cf. Figure 2(a)), with respect to the inner integration radius rb and outer integration radius ra. These two radii define the zone over which G is computed in both elastic and elastoplastic cases. The analysis consists of fixing one of the radii and observing the convergence of G as a function of the other. By setting rb = 0 and varying ra, the results for the purely elastic case show that G converges rapidly when applying the usual rule of choosing ra equal to three times the minimum element size of the mesh (cf. Figure 2(b)). In contrast, for the elastoplastic case, convergence of G is only achieved for ra values greater than 0.2 mm. This suggests that the outer radius ra should be at least three times larger than the estimated size of the plastic zone to ensure the reliability of the calculation.
left: (a)
right: (b)
Figure 2: (a) Single edge notched specimen subjected to a tensile load of 50 MPa and the corresponding mesh. (b) Comparison of the evolution of the G value for the elastic and elastoplastic cases as a function of the external integration radius ra.
In parallel, when the external radius ra is fixed, the results for G as a function of the internal radius rb show that the latter has no significant influence on G, provided that the distance between the two radii is sufficiently large. However, when ra is smaller than the size of the plastic zone, variations in rb affect the calculation of G.
This analysis refines the optimal selection of these radii, particularly for the elastoplastic case where the plastic zone extends beyond three mesh elements from the crack front, a length that previously defined ra in the elastic case.
The influence of the gradient near the crack front on the values of G was also investigated. The viscoplastic strain gradient varies over a region composed of a few elements from the crack front, thereby affecting the values of G depending on the selected parameter hmin, which corresponds to the minimum distance from the crack front beyond which the inelastic strain gradient is considered. The parameter hmin allows the exclusion of the gradient’s influence in the immediate vicinity of the crack front, where significant fluctuations may occur that could disrupt the evaluation of the mechanical field. These sensitivity studies were also carried out on an unstructured 2D mesh and under large deformation conditions, in order to assess the convergence of G values in both elastic and elastoviscoplastic cases.
This initial validation was then extended to three-dimensional cases using a planar geometry with dimensions 1 mm × 1 mm × 0.5 mm, subjected to a stress σ₂₂ = 110 MPa. Within this geometry, which features a symmetry plane at z = 0, an edge crack was inserted as illustrated in Figure 3(a) below.
The accumulated plastic strain field (epcum) was then computed at the surface of the geometry (z = 0.5 mm) as well as at its mid-plane (z = 0). The results show a slightly different field at these two locations (cf. Figure 3(b)). This difference observed in the plastic zone between the mid-plane and the surface can be explained by variations in the stress and strain fields through the thickness.
We then analyzed the convergence of G, both at the surface and at the core, by varying the external radius ra while keeping the internal radius rb constant at zero. The results of this analysis, shown in Figure 3(c), indicate that the difference in G between the core and the surface is negligible, in both the elastic and plastic cases, and can therefore be disregarded.
(a) - (b)
(c)
Figure 3: (a) The considered geometry, with a through-thickness edge crack, and the corresponding mesh. (b) The accumulated plastic strain field (epcum) showing the plastic zone at the surface and at the core of the geometry. (c) Results of the convergence analysis of G as a function of the external radius ra, at the surface and at the core, for both elastic and plastic materials.
However, for a through-thickness crack, the convergence of G differs between the surface and the core of the plate in the plastic case, as illustrated in Figure 4 below. Specifically, convergence is faster at the core than at the surface due to variations in the plastic zone through the thickness. In contrast, in the elastic case, G converges at the same rate at both the surface and the core, although the limiting values differ, which is not observed in the elastoplastic behavior.
Figure 4: Results of the convergence analysis of G at the surface and at the core as a function of the external radius ra (with rb = 0), for a through-thickness crack, considering both plastic and elastic material behaviors.
In a final step, we analyzed the effect of the interaction between two crack fronts on the value of G. The study focuses on a plate measuring 2 mm × 2 mm × 0.2 mm, subjected to a uniform tensile stress σ₂₂ = 150 MPa applied on its top and bottom faces. A through-thickness crack with a length of 2a = 0.5 mm is introduced (see Figure 5).
Convergence calculations of G at the core of the geometry were performed as a function of the external radius ra for both elastic and plastic behaviors. To reduce computation times, symmetry conditions were applied in two configurations. In the first case, the model was reduced to half of the geometry by exploiting symmetry with respect to the mid-plane z = 0, where the displacement condition U3 = 0 was imposed. In the second case, the reduction was extended to one quarter of the domain by introducing a second symmetry plane at x = 0, where the displacement condition U1 = 0 was applied in addition to U3 = 0 at the plane z = 0.
Figure 5: Results of the convergence analysis of G at the core (z = 0) as a function of the external radius ra (rb = 0, hmin = 0), in the case of two interacting crack fronts, for (a) elastic and (b) plastic behaviors.
The analysis results show that the convergence of G for the first crack front is influenced by the presence of the second. Indeed, in the first case, the values of G decrease significantly as the external integration radius approaches the second crack front (ra = 0.5 mm), thus disrupting the convergence of the results (see the blue curves in Figure 5).
In the second case, the proximity of the geometric boundary (ra = 0.25 mm) affects the results (see the orange curves in Figure 5), causing a decrease in the values of G beyond this limit. Although the two configurations are mechanically equivalent, the reduced integration zone beyond ra = 0.25 mm explains this divergence.
This study confirmed the relevance of the G-theta method for predicting crack propagation in materials with elastoplastic behavior. Integrated into the Z-cracks module of the Z-set suite, this approach enables the calculation of the energy released during crack growth, even in complex industrial contexts.
Thanks to its 3D adaptive remeshing and optimized post-processing, Z-cracks provides a high-performance solution for three-dimensional crack analysis. Z-set thus establishes itself as a powerful platform, ready to meet current and future challenges in fracture mechanics.
References