An Alternative Finite Element Formulation to Predict Ductile Fracture in Highly Deformable Materials

As presented by Cordeiro [1], engineering projects have become increasingly complex, and, thus, the accurate prediction of the behavior of mechanical and structural components is essential. In this context, an interesting topic of research is the analysis of ductile fracture in highly deformable materials, which requires sophisticated material models. The structural integrity prediction is crucial in many engineering applications such as nuclear reactors [2]. Oliveira and Andrade [3] address the extensive use of machines employing rotating shafts, as well as vehicles, emphasizing the effect of cracks and notches on the response of such components. Furthermore, Jesus [4] presents a study concerning the importance of investigating failure mechanisms, illustrating it with the aeronautical accident involving the EC225 helicopter in Norway, caused by the fracture of a metallic mechanical component. Since ductile damage involves a progressive deterioration of material strength and large plastic deformations, those engineering applications reinforce the need for a model capable of reproducing substantial levels of inelastic deformation, damage evolution, and, in extreme situations, component fracture.

According to Callister Jr and Rethwisch [5], two types of fractures can occur in metallic materials: ductile and brittle. In the case of ductile fracture, the material exhibits nonlinear behavior with considerable plastic deformation [6]. In addition, ductile fracture in a bar can be defined as follows [7]: small cavities form in the cross section; these cavities grow due to plastic behavior, approach each other, and coalesce until they form an elliptical crack. This crack has its major axis perpendicular to the direction of stress application. The crack continues to expand in the direction parallel to its axis until fracture finally occurs. Therefore, ductile damage also involves the evolution of microcavities (or internal voids) in the material [8,9]. Additional phenomena include void nucleation from second-phase particles and shear banding, which corresponds to the localization of plastic strain and damage along a narrow band [10], as well as void shearing [11,12]. As expected, considering additional mechanisms may improve the accuracy, but results in a more complex formulation. Brittle fracture, in turn, involves negligible plastic deformation and is not in the scope of the present work. An example of a finite-strain formulation developed for ductile-brittle transition can be found, for instance, in the work of Hütter et al. [13].

From the theoretical point of view, the Rice J-integral approach from fracture mechanics is employed to quantify the resistance of standard cracked specimens under elastoplastic regime. However, that global approach presents some major drawbacks: it strongly depends on the specimen geometry (it is not an intrinsic material property); it is restricted to preexistent cracks and simple geometries; and it cannot reproduce crack initial and propagation. To circumvent such limitations, several material constitutive equations and computational tools have been developed to simulate damage evolution, crack initiation, and propagation, as well as ductile fracture in the context of the finite element method (FEM). They can be classified into the following broad categories [14]: continuous approaches, discontinuous approaches, and alternative methodologies. In the first case, the material degradation process is assumed to be continuous and the objective is to quantify the effect of the void evolution on the mechanical response. As pointed out by Shakoor et al. [14], the advantage of continuous models is that the damaged zone can evolve, grow, and coalesce without any numerical difficulty. In the context of metallic alloys, the most used ductile damage formulation is probably the micromechanics-based Gurson–Tvergaard–Needleman (GTN) formulation and the phenomenological model of Tvergaard and Needleman [15]. The former is the result of the modifications proposed by Gurson [16] in the Gurson's porous plasticity model [17], including experimental adjustments, void nucleation, and coalescence, as well as material failure. Both formulations are compared, for instance, in the current state of art review presented by Fincato and Tsutsumi [18]. Basically, it can be concluded that, although the Lemaitre's model is thermodynamically consistent, only the GTN formulation fulfills the mass conservation principle and represents the plastic volume variation (which are important physical conditions), as well as involves model parameters with physical interpretations. In addition, the coupling between the level of damage and the elastic material stiffness, captured only by the Lemaitre formulation, is only relevant at the end of the fracture process. Another relevant aspect is that damage closure (or reduction in the void volume fraction) is directly accounted for by the GTN model. As an example of superior performance of the GTN model compared to the Lemaitre formulation in experimental tests, one can cite the work of Delgado-Morales et al. [19] with AISI 4340 steel alloy. Examples of GTN model applications for ductile fracture analysis can be found in Refs. [20–24], among many others.

Continuous damage models can be easily incorporated into a finite element computer code to simulate ductile damage evolution. In this case, the fracture can be reproduced by identifying the points and regions reaching the failure limit of the material (or the critical damage level). Once they are identified, two simple strategies can be employed: the fully damaged zone or the element erosion (or deletion) method. In the first case, a numerical failure limit is adopted (usually 90% of the failure porosity level in the GTN model), and the zone composed of the fully damaged points is assumed to behave as a crack, accounting for both initiation and propagation phenomena. The fully damaged zone technique is successfully used, for instance, in the works of Refs. [23,25], in the context of thermoviscoplastic materials (i.e., considering inelastic deformations, temperature changes, and strain rate dependence). In the second case, based on an appropriate fracture indicator, the damaged elements are removed from the mesh, which requires a modification in the numerical algorithm. Alternatively, an elastic behavior with a very small stiffness can be applied to the fractured elements.

Nevertheless, finite element formulations based on continuous models usually suffer from mesh dependence, as well as strain and damage localization. To avoid such limitations, the other two categories (discontinuous approaches and alternative methods) have been proposed. Regarding discontinuous formulations, the most used techniques are the enriched FEM and the cohesive zone models (CZMs). In the case of enriched finite elements, as the extended FEM (or X-FEM), the discontinuities (cracks or interfaces) can be modeled independently of the mesh, introducing a strong discontinuity in the displacement field. As a result, both crack initiation and propagation can be captured efficiently. On the other hand, enriched finite elements cannot reproduce the energy dissipation rate along the crack opening process, as the local load-bearing capacity is instantaneously lost once the displacement jump is introduced [14]. To model the energy dissipation rate in a more realistic way, CZMs are employed. The idea is that there is a finite size transition region between a fully damaged material and a bulk (or sound) material. Moreover, it is assumed that the energy dissipation is due to the breakage of atomic bounds and corresponds to the amount of energy required to create new surfaces, which is more physically based. Despite the advantages of discontinuous models, they are quite complicated to be implemented, requiring additional numerical techniques. Besides, reproducing cracks initiated at multiple sites remains a major limitation of discontinuous approaches, especially in 3D. Examples involving ductile fracture analysis with X-FEM and CZMs can be found in Refs. [26,27], among many others.

The third group corresponds to the alternative techniques, such as peridynamics, phase-field model, and meshless methods (MMs). According to the review of the peridynamics theory presented by Dorduncu et al. [28], an important drawback of usual finite elements is that the spatial derivatives required do not exist along the crack discontinuity. The partial differential equations are replaced by an integrand that does not contain any spatial derivatives of displacements. Nonlocal interactions are modeled through a network of bounds and damage is reproduced by eliminating those interactions. The main feature of peridynamics is modeling both damage initiation and evolution by bond breakdown. In the phase-field model, in turn, a variational approach to regularize the energy functional is employed and the crack surfaces are weakly described by a continuous variable, which represents the smooth transition between the bulk material and the discontinuous crack. With this method, complicated crack geometries can be reproduced in terms of initiation, propagation, coalescence, and branching, avoiding the need of tracking algorithms (as in X-FEM), cohesive elements, and remeshing techniques (which is costly and difficult in 3D problems). A review of phase-field models applied to ductile fracture analysis, including theories and computational aspects, can be found in Ref. [29]. Another way to deal with discontinuities is the MMs, in which the approximation is built from nodes only, eliminating some difficulties concerning the approximation dependence on the mesh. According to the review presented in Ref. [30], one of the advantages of MMs is that problems involving moving discontinuities (as crack propagation) can be easily modeled. On the other hand, as in the other approaches, despite the advantages of such alternative techniques, there are some important limitations: the peridynamic model is limited by the definition of nonlocal boundary conditions and surface effects [31]; the phase-field model results in a widening of the damaged zone as the crack approaches the external boundaries [31]; the mathematical rigor of the phase-field approach makes it difficult to include more sophisticated physical mechanisms [14]; and the definition of essential boundary conditions in meshless methods is not simple, as in traditional FEM [32]. In addition to the more complicated formulations, such alternative techniques in general results in a highly computational cost, which is a major disadvantage when compared to the other methods.

An alternative way to keep the advantages of continuous approaches and, at the same time, to alleviate the inherent mesh dependence and damage localization problems is the use of nonlocal models. According to Pascon and Waisman [25], the damage localization leads to ill-posed boundary-valued problems, requiring an additional physical length scale parameter to regularize the partial differential equations. Further details can be found, for instance, in the works of Refs. [33–35]. The most common procedure to avoid localization problems in continuous formulations is to use a nonlocal approach, in which the behavior of the material at some point depends not only on the local values of state variables but also on the corresponding values in a domain around such point. Examples of nonlocal formulations include the integral-type model [36], the gradient-enhancement formulation [37], the higher-order microstrain theory [24], and the phase-field regularization [38], already addressed. According to the review of continuous models presented by Besson [39], the integral-type and gradient-enhanced formulations are equivalent to each other and simple to implement. However, numerical convergence is not assured in the first case and the second strategy requires the evaluation of the gradient of the nonlocal variable at the Gauss point. Moreover, in the higher-order theories presented in Ref. [24], additional degrees-of-freedom are needed, which demands more computational effort and memory capacity. In the work of Ref. [25], it is demonstrated that using the porosity as a nodal degree-of-freedom together with a gradient-enhanced formulation reduces mesh dependence and provides more consistent porosity contours, since the singularity at the crack tip is attenuated in such models [31]. Additionally, it is shown that the number of iterations required along the simulation is considerably smaller for the nonlocal model, compared to the local formulation, in the slow strain rate situation. From the numerical point of view, the enrichment proposed in Ref. [25] requires only the modification of the porosity field equations, keeping the algorithm structure unchanged and avoiding additional numerical techniques.

Moreover, even with extensive literature available in the field, the development of numerical tools capable of predicting fracture in ductile materials, especially those highly deformable, remains a major challenge. Many numerical formulations within the context of large strain GTN models have been developed along the last decades in order to predict ductile fracture (see, among many others, the works of Refs. [39–45]). In such formulations, thermodynamic principles are employed, the strain decomposition is based on the Kröner–Lee split of the deformation gradient [46,47], and the stress measure adopted is usually the Kirchhoff or the Mandel stress tensor, which means that nonlinear continuum mechanics is required [48,49]. The resultant mathematical treatment is considerably more complex than the small strain regime counterpart.

Another relevant issue is that the standard numerical procedure includes the elastic prediction and plastic correction phase, which is a staggered scheme. In this case, the global variables (or degrees-of-freedom) are initially obtained and, if necessary, another system of equations is solved in order to find the internal variables, related to plastic strains and damage, for example. According to McAuliffe and Waisman [50], the accuracy of staggered schemes is questionable in the context of coupled constitutive models (like coupled damage–plasticity), as the accumulated numerical errors affect the results and may not be recovered in the following steps. Those authors also highlight that the global or monolithic solution (solving all the variables simultaneously) is more stable and accurate. Moreover, in usual elastoplasticity formulations, the consistency condition derived from the yield criterion is employed to set the material tangent stiffness regarding the stress–strain relationship along the elastoplastic regime. Another strategy, rarely exploited, is to embed the influence of the internal variable evolution on the material stiffness concerning the global variables. In the case of mixed formulations (in which not only nodal displacements are used as degrees-of-freedom), such an alternative strategy may be more convenient, since it directly incorporates the derivatives related to the internal variables into the derivatives of the global system of equations.

The purpose of this study is to develop an alternative finite element formulation to predict ductile damage and fracture in highly deformable materials. The proposed model is intended to be a simple and accurate way to reproduce ductile fracture and finite-strain levels by means of a nonlocal continuous model, avoiding any special and complicated numerical technique. To this end, the GTN constitutive model is employed together with a mixed formulation. A linear-order triangular finite element is adopted, and the degrees-of-freedom are current positions and nonlocal porosity at the nodes. Associative plastic flow rule, isotropic hardening, and porosity evolution based on growth, nucleation, and coalescence are also adopted. It should be emphasized that rate-independent response, isothermal, quasi-static, and plane-stress conditions are assumed. Another feature of the present work is the development of a finite-strain ductile damage model using a mixed finite element formulation based on a numerical technique of incorporating the nonlinear system of internal variables into the global system of equations regarding the degrees-of-freedom. In addition, the present formulation can be considered as a finite-strain extension of the model proposed in Ref. [25], in which the global system involves all the nodal and Gauss point variables (displacements, nonlocal porosity, and internal variables), which results in a larger system size and, thus, in a higher computational effort. Therefore, the numerical strategy of incorporating the internal variable derivatives into the main system is assessed in the current work within the context of finite-strain levels. The formulation proposed is also similar to the thermomechanical model of Håkansson et al. [42], in which an eight-node plane–strain isoparametric element is used together with reduced integration.

In this section, the adopted constitutive model is described, including measures of deformation and stress, yield criterion, ductile damage variable, plastic flow law, material hardening law, and elastic response.

The constitutive model employed in this work is restricted to homogeneous materials under isothermal conditions, a plane-stress state, and quasi-static loadings.

The numerical strategies adopted are provided in this section, comprising the finite element and the numerical solution of the governing equations.

The performance of the proposed numerical formulation is demonstrated by three mechanical problems involving highly deformable materials under ductile fracture: a double-notched specimen, the necking of a bar, and an inhomogeneous tension situation. For all the simulations, the material parameters adopted are provided in Table 1, corresponding to a structural steel with an initial yield stress of 312 MPa, as in the thermoviscoplastic GTN model of Pascon and Waisman [23]. Each problem has been analyzed with different discretizations in order to assess the effect of the mesh refinement on the mechanical response. The unstructured meshes have been generated in program Gmsh [66].

All the simulations have been carried out by means of displacement control, dividing the analysis into 10,000 steps due to the highly nonlinear character of the algebraic equations involved. Moreover, to avoid convergence problems, the damage level stops evolving once it reaches 90% of the porosity failure limit $(0.9fF=0.225)$⁠, even if the material remains in elastoplastic regime. In addition to the mesh refinement, the influence of the nonlocal parameter $aNL$ defined in Eq. (20) is studied in each numerical example, as in Ref. [41]. Nevertheless, in that work, the effect of the nonlocal parameter on the evolution of the porosity field is not addressed in detail. Only the reduced mesh dependence provided by the nonlocal formulation is shown.

One should highlight that, although isothermal and quasi-static conditions are considered in the present work, temperature and rate-dependent effects may be relevant along the ductile damage process (see, for instance, the work of Pascon and Waisman [23]). Moreover, in the postprocessing stage, the porosity field is obtained at the end of some specific steps, and the fully damaged zone defined by $f=0.225$ is considered to represent the crack initiation and propagation phenomena.

The first problem is the specimen with two asymmetrically rounded notches depicted in Fig. 1. The present numerical example has been proposed by Brokken [67] in order to evaluate the discrete cracking procedure developed for modeling the blanking process [68]. The left and top boundaries undergo a prescribed vertical displacement, whereas the right and bottom boundaries remain fixed. As the initial damaged zone is expected to take place in the round notches, the mesh refinement for all the discretizations adopted has been concentrated around such areas, see Fig. 2. The plane-stress conditions are ensured by using a small thickness of 1 mm.

According to Fig. 3, a more flexible response is obtained with mesh refinement and by reducing the nonlocal parameter $aNL$⁠. In all cases, one can observe the usual ductile material response: an initially linear elastic behavior, followed by strain-hardening in elastoplastic regime and a drop in the stress level due to the evolution of damage (or porosity). Mesh dependence seems to be less noticeable for the nonlocal formulation with $aNL=0.01$ (Fig. 3(b)). Regarding the most refined mesh, increasing the nonlocal parameter leads to an additional resistance to damage evolution and, thus, a less-pronounced porosity-induced softening behavior, providing a more ductile behavior (Fig. 3(d)). In the work of Mediavilla et al. [68], only two meshes composed of plane–strain bilinear quadrilateral elements are employed, showing a very different behavior especially in the crack propagation stage.

Considering the most refined mesh (366 nodes and 670 elements), the effect of the nonlocal parameter in the evolution of porosity field is analyzed. As shown in Fig. 4, the fully damaged zone $(f=0.225)$ starts at the notches and propagates toward the inclined line connecting them. A similar behavior has been reported in the work of Aldakheel et al. [44] considering a modified GTN model coupled with phase-field formulation. In the study by Mediavilla et al. [68], in turn, only the damaged zone in the upper notch evolves toward the other one. Moreover, it should be highlighted that the nonlocal parameter reduces the damage concentration in the elements outside the fully damaged zone, providing a smoother porosity field evolution, which is also captured in the nonlocal small strain GTN model of Pascon and Waisman [25]. As far as the authors know, the effect of eliminating the damage concentration along the element edges by means of a nonlocal and continuous finite-strain model has not been addressed in scientific literature yet.

The second example is the necking of a rectangular bar, depicted in Fig. 5. Due to the symmetry conditions, only one half of the specimen is discretized. In this case, three meshes are employed with a mesh refinement level concentrated around the center of the bar. In order to trigger the necking instability, the upper and bottom faces are restricted along the horizontal direction, instead of using a geometric imperfection. In addition, all the nodes along the central vertical line are restricted along the horizontal direction to reproduce the symmetry condition. The prescribed vertical displacement is applied to the upper face, whereas the bottom edge is fixed. As in the first example, the selected thickness in this case is 1 mm. The necking of a rectangular bar has been extensively studied in scientific literature considering finite-strain elastoplasticity coupled with ductile damage (see, for instance, [36,44,69–71]).

The convergence analysis in terms of the applied forces is provided in Fig. 6 together with the influence of the nonlocal parameter on the behavior of the most refine mesh. One can note that forces are practically the same for all the meshes considering the nonlocal formulations before the necking instability, indicating that the nonlocal model reduces mesh dependence. Once again, increasing the nonlocal parameter provides a higher ductility.

The differences comparing the meshes employed are more pronounced after the necking instability, which shows the difficulty in capturing the instable behavior. Such an effect is also observed in the nonlocal formulation presented by Nguyen et al. [63], in which a plane strain specimen undergoes a necking instability. However, the differences in that work are slightly larger for distinct values of the nonlocal parameter, regarding the unstable stage of the force–displacement graph. In the study by Leclerc et al. [36], only two values of the nonlocal parameter are employed and the effect on the force–displacement graph is not so clear.

As in the work by Aldakheel et al. [44], the initial failure occurs at the center of the specimen and propagates along an inclined direction (see Fig. 7), resembling the cup-cone failure mode (see, for instance, the experimental results of the study by Li et al. [72]). The nonlocal parameter also provides a smoother fully damaged zone, as well as delayed initiation and propagation phenomena. Comparing the two cases depicted in Fig. 7, the nonlocal formulation $aNL=0.1$ prevents the occurrence of a region with a small damage level (in red) surrounding the fully damaged zone. In addition, one can observe the highest transverse displacement levels around the center of the specimen, confirming that the strategy of restricting all the displacements of the upper and bottom edges is effective to capture the necking instability without using any geometric imperfection.

The third example is the specimen under a state of inhomogeneous tension depicted in Fig. 8. The geometric data have been extracted from the study by Miehe et al. [73]. According to those authors, the objective of the hyperbolic shape proposed is to induce the fracture initiation at the center of the specimen. The bottom face is restricted along the horizontal and vertical directions, while the upper face is restricted only in the horizontal direction with a prescribed vertical displacement increasing up to 0.775 mm, which corresponds to a nominal average strain of 10%. The selected thickness for the specimen is 0.1 mm.

For the nonlocal formulation, mesh convergence is achieved considering the stage before the fracture initiation, which is represented by the sudden drop in the force level (see Fig. 9). Similarly to the other examples, the increase in ductility provided the nonlocal parameter is more evident after the fracture initiation. In the formulation of Miehe et al. [73], two nonlocal parameters are employed: one for the phase-field regularization and another for the gradient-plasticity model. The adopted values in those works have a noticeable effect on the ductility and fracture initiation of the specimen under inhomogeneous tension. Nevertheless, it is difficult to calibrate such nonlocal parameters.

Regarding the most refined mesh employed (227 nodes and 424 elements) and the nonlocal formulation $aNL=0.1$⁠, the porosity evolution and some strain fields are provided in Figs. 10 and 11. As expected, the fully damaged zone starts at the center of the specimen and propagates along a slightly inclined direction toward both edges. In the work of Miehe et al. [65], the crack propagation occurs along a line with an inclination of 45 deg, following the shear banding. The difference is probably related to the distinct constitutive model adopted. Moreover, the strain localization at the center of the specimen can be observed (see Fig. 11). One should also highlight that the extremely high levels of maximum strains (up to 721.4%) reinforce the need for a finite-strain constitutive model. Even the transverse normal strain achieves considerably high values at the end of the simulation (260.2%). Another important aspect in Figs. 10 and 11 is the fact that the response is not symmetrical and, thus, the mesh should involve all the specimens and not only one half of it.

In this work, an alternative finite element formulation to predict ductile damage and fracture in highly deformable materials is presented. A large-strain GTN constitutive law for porous metallic materials is employed to describe plasticity-induced damage, together with a Lagrangian description. It is worth mentioning that plane-stress conditions are employed, rather than the usual assumption of plane strains. The numerical model adopted is based on a mixed triangular finite element of linear degree, whose degrees-of-freedom are current positions and nonlocal porosity levels at the nodes. A numerical strategy to capture the change in material stiffness along the evolution of the internal plastic variables is developed, instead of using the classical elastoplastic stress–strain relationship based on the consistency condition. The phenomena of crack initiation and propagation are reproduced via the strategy of the fully damaged zone based on the nonlocal porosity field. One should note that the present model is restricted to isotropic hardening, associative plastic flow rule, rate-independent material response, isothermal conditions, and quasi-static loading, as well as constant-strain triangular elements.

The performance of the proposed formulation is assessed by means of three numerical examples involving ductile fracture and finite plastic strain levels: a double-notched specimen under tension; the necking instability of a rectangular bar; and a specimen with variable width, resulting in a state of inhomogeneous tension. In the first case, the presence of the two round notches leads to a 45-deg propagation similar to a shear loading, even in a tension problem. In the second example, the present formulation is capable of capturing the necking instability and a final damaged zone resembling the cup-cone fracture, usually observable in ductile materials. In the third problem, a slightly inclined propagation surface is obtained, together with remarkably large strain levels.

All the problems have been simulated varying the number of elements and the value of the nonlocal parameter. Results show that convergence can be obtained with mesh refinement, and increasing the nonlocal parameter provides a more ductile response. One should highlight that, even for the most refined meshes, the local formulation leads to a strange final porosity field, which is not observed in experimental tests.

Some features of the present study can be cited. First, a constitutive model capable of reproducing ductile damage and large strains is essential for ductile fracture analysis of metallic materials. Providing the strain levels along the fracture process is not usual in scientific literature, although it is necessary in order to justify the finite-strain damage model. Second, the three numerical examples involving ductile fracture are solved with a nonlocal continuous model implemented in a mixed finite element formulation, avoiding the need for special numerical techniques. Finally, the effect of the nonlocal parameter on the mechanical response is addressed in detail. It is shown that the nonlocal parameter results in a smoother porosity field, regarding the elimination of the damage concentration along the element edges and, thus, reducing the damage gradient across the fully damaged zone. Such analysis within the context of finite-strain continuous model is rarely performed in the scientific literature.

Therefore, the present formulation is promising regarding analysis of ductile fracture under large strains. Future studies include extension to thermoviscoplasticity, phase-field coupling, cohesive zone models, more refined meshes, and improvements in the nonlocal strategy in order to reduce or eliminate mesh dependence.

The authors appreciate all the support provided by the Materials Engineering Department of Lorena School of Engineering of the University of São Paulo (DEMAR/EEL/USP), the Department of Civil and Environmental Engineering of the Bauru College of Engineering of the São Paulo State University (DEC/FEB/UNESP), and the financial aid granted by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), which allowed to acquire a workstation to carry out the numerical simulations.

There are no conflicts of interest.

The authors attest that all data for this study are included in the paper.

C =

right Cauchy–Green stretch tensor

E =

Green–Lagrange strain tensor

F =

deformation gradient

I =

second-order identity matrix

S =

second Piola–Kirchhoff stress tensor

Y =

vector of current nodal positions

$f$ =

effective porosity

$p$ =

hydrostatic pressure (mean stress)

$f¯$ =

nonlocal porosity

$E$ =

Young's modulus

$aNL$ =

nonlocal parameter

$f0$ =

initial porosity level

$fC$ =

coalescence limit of porosity

$fF$ =

failure limit of porosity

$rf$ =

scalar function for porosity evolution

$rκ$ =

scalar function defining the evolution of $κ$

$wf¯$ =

nonlocal porosity weighting function

$fext$ =

vector of external forces

$fint$ =

vector of internal forces

$rF$ =

force residuum vector

$rM$ =

residual vector regarding main variables

$rI$ =

residual vector regarding internal variables

$Dκ$ =

void nucleation function

$Meq$ =

equivalent Mandel stress

$Wp$ =

rate of plastic work (or plastic power)

$Dp$ =

plastic strain rate tensor

M_e =

elastic Mandel stress tensor

$RCp$ =

plastic flow tensor

$YM$ =

residual vector of main variables

$YI$ =

residual vector of internal variables

$fN,sN,εN$ =

nucleation parameters

$q1,q2$ =

GTN yield function parameters

$zΔλ,zκ,zf,zCp$ =

residual functions of internal variables

$K,ε0,n$ =

swift coefficients

$JMM,JII,JIM,JMI$ =

Jacobian matrices regarding main and internal variables

$κ$ =

equivalent plastic strain

$ν$ =

Poisson's ratio

$σY$ =

yield stress

$ϕ$ =

yield function

$φ$ =

pressure-dependent function