Simulation of crack propagation based on eigenerosion in brittle and ductile materials subject to finite strains

In this paper, a framework for the simulation of crack propagation in brittle and ductile materials is proposed. The framework is derived by extending the eigenerosion approach of Pandolfi and Ortiz (Int J Numer Methods Eng 92(8):694–714, 2012. https://doi.org/10.1002/nme.4352) to finite strains and by connecting it with a generalized energy-based, Griffith-type failure criterion for ductile fracture. To model the elasto-plastic response, a classical finite strain formulation is extended by viscous regularization to account for the shear band localization prior to fracture. The compression–tension asymmetry, which becomes particularly important during crack propagation under cyclic loading, is incorporated by splitting the strain energy density into a tensile and compression part. In a comparative study based on benchmark problems, it is shown that the unified approach is indeed able to represent brittle and ductile fracture at finite strains and to ensure converging, mesh-independent solutions. Furthermore, the proposed approach is analyzed for cyclic loading, and it is shown that classical Wöhler curves can be represented.


Introduction
Engineering structures often consist of various different materials exhibiting a brittle or ductile behavior. As one example in the context of mechanized tunneling, consider the tools attached at the cutting wheel, namely the cutting disks or the chisels. These tools are usually composed of ductile steel with wear-resistant, brittle armoring, made of hard metal or metal matrix composites protecting the steel from wear. These armoring materials make in turn also use of microstructures consisting of brittle and ductile components, namely a ductile metallic matrix with embedded brittle ceramic particles. Whereas fracture at the macroscale has to be considered to evaluate potential failure of the tool, sub-critical crack propagation under cyclic loading is particularly important at the microscale as it is mostly responsible for wear in terms of surface spalling. This is just a representative example for a variety of engineering applications where brittle and ductile fracture play a role. However, most algorithms are rather specialized to either brittle or ductile fracture and thus, the connection of both is not necessarily straightforward, especially with view to ensuring mesh-independent simulations.
Brittle crack propagation has already been investigated by Griffith [7]. There, it was discovered that a crack propagates by an area increment in an elastic material if the elastic energy surpasses a certain energy level. This criterion is known as the Griffith criterion in which the constant energy level per crack area increment is defined as the Griffith energy release rate. This constant is equal to the surface energy that is needed to produce a certain surface. In this approach, it is assumed that this process works irreversibly and that the crack surfaces are stress-free. Later, Irwin [9] investigated metals with ductile material behavior and discovered that plasticity plays an important role for the crack propagation. For example, necking of tensile samples occurs before cracking [26] and the crack path changes. Furthermore, Irwin extended the original Griffith criterion for brittle materials to the Griffith-type criterion for ductile materials by additionally considering the energy dissipated through plasticity at the crack tip before the crack propagates.
A special form of crack propagation is associated with fatigue due to cyclic loading, which was recently investigated by Zenner and Hinkelmann [43]. Motivated by the investigation of train axles which broke in an incident in the nineteenth century [41], it was discovered that structures fail at loads under cyclic loading lower than the maximum load in a single load cycle. This results from macroscopic effects like macroscopic plastic deformations and sub-critical crack propagation as well as microscopic effects like plasticity at the microscale, void nucleation and microcracks as seen in Pineau et al. [30]. The basis for the analysis of this effect is the well-known Wöhler experiment which consists of series of single experiments to identify the number of cycles until failure for varied load amplitudes yielding the so-called Wöhler curve. With these, life time predictions and improvements of engineering structures can be made.
Although crack propagation governs phenomena important to engineering, like failure upon cyclic loading or surface spalling, its simulation is still a challenging task in engineering science. The reason for that is the requirement of continuity in the mechanical fields for the application of continuum mechanics which conflicts with the discontinuity imposed by the cracks. To circumvent numerical problems associated with the continuity requirements of the classical finite element method (FEM), which is mostly applied to engineering problems, multiple approaches have been made. For example, mesh-free methods, like e.g., the Peridynamics [17,35], do not necessarily require continuum mechanical descriptions and thus, a priori avoid these conflicts. However, the absence of a direct link to continuum mechanical descriptions for the bulk materials may pose different challenges. Other approaches are extensions of the finite Element framework, e.g., the incorporation of strong discontinuities within the elements [36] or through enriched shape functions in the sense of the extended FEM (XFEM) [5]. But these approaches are conceptionally not straightforward and technically demanding. A further alternative is the phase-field method as in Miehe et al. [21], where sharp cracks are smoothed by a continuous damage field. This approach has the disadvantage that this separate damage field requires usually an additional degree of freedom rendering it computationally demanding. Quite similarly, the introduction of gradient-enhanced damage formulations describes the crack in a smeared sense. Here, the need for the additional degree of freedom could be eliminated by the neighbored element method proposed in [11], see also [12] for finite strains. A further alternative in the context of smeared damage fields is based on the convexification of incremental variational formulations to allow for mesh-independent simulations, cf. e.g., [3]. However, the link of the smeared damage fields to failure criteria is still an open problem. Furthermore, cohesive-zone models as introduced in e.g., Barenblatt [4] and Dugdale [6] have been considered for the crack simulation using an FE framework (Hillerborg et al. [8]). There, interface elements describing cohesion during cracking are included between the regular elements where the crack propagates. Thus, the crack paths are enforced to follow the element edges, which may lead to erroneous crack patterns as shown in e.g., Schellekens and De Borst [33]. Another class of approaches is based on element erosion where single elements are removed from the calculation as soon as potential crack growth is detected within the element. This concept does only allow for an approximate representation of the crack geometry and usually requires very fine meshes. On the other hand, its implementation and handling is relatively straightforward and thus, it can be directly applied in a broad range of engineering problems. As a loss of ellipticity of the differential equations in structural boundary value problems is observed upon failure, mesh-dependent solutions, i.e., non-converging numerical schemes, can result therefrom. In order to avoid this, all the above-mentioned approaches for crack propagation make use of nonlocal enhancements of the originally local formulations. Whereas this is directly included in mesh-free, phase-field, or gradient-damage formulations, strong discontinuity approaches make use of nonlocal evaluations of stress-concentration factors. In the context of erosion algorithms, the first consistent nonlocal approach enabling mesh-independent calculations, the so-called "eigenerosion," was proposed in [27], based on the variational eigenfracture formulation proposed in [34]. This eigenerosion approach is based on the Griffith criterion for crack propagation and achieves its mesh independence from the regularization of the crack area by evaluating a specifically defined spatial -neighborhood. In Schmidt et al. [34], it was proven that -convergence of the regularized energy dissipation functional to the unregularized one is obtained as the length scale parameter ε goes to zero. This leads to an independence of the crack path, crack propagation and therefore of the overall structural response from the spatial discretization. The eroded elements are able to undergo eigendeformations for which no additional work has to be imposed in those.
Variations of this method within a mesh-free framework have been developed subsequently [16,24,28], for example, the simulation of high-impact loading and fragmentation. The eigenerosion approach has been originally developed for brittle crack propagation at small strains. One extension to ductile crack propagation at small strains with Drucker-Prager plasticity has been presented in Qinami et al. [31].
In this paper, we extend the eigenerosion approach to the geometrically nonlinear setting of finite deformations within the finite element framework. Furthermore, we propose an energy-based failure criterion of Griffith's type which can be directly used in the eigenerosion algorithm to extend its applicability to ductile crack propagation. In order to account for sub-critical cracks in metals, the finite strain J 2 elasto-plastic model of Simo and Miehe [38] is considered. Additionally, the elastic energy density is split into a tension and a compression part in order to ensure that the crack surfaces are stress-free under tension but still transfer forces under compression. Therefore, the asymmetric material behavior as seen in Miehe et al. [20] can be accounted for. This asymmetry is necessary for simulations in which the load direction changes, for example for the simulation of a specimen under cyclic loading, as it is essential for the calculation of the Wöhler curve.

Eigenerosion algorithm
Since the extensions proposed in this section are based on the original eigenerosion approach [27,34], we start by recapitulating the main equations governing the conceptual framework. The theory of brittle fracture is based on the work of Griffith [7]. There, the existence of a Griffith-type energy release rate is assumed. Herein, the potential energy U stored by imposing mechanical work depends on the displacement vector u and |C| denotes the area of the crack set C. Irreversibility of crack propagation and no healing of the material are taken into account. Consider a discrete time increment t, then the monotonicity constraint holds for the crack set C(t) at time t and the crack set C(t + t) at a later time step t + t. If the energy release rate G reaches the critical value of Griffith energy release rate G c , the crack propagates with the crack front velocity v. In that case, G − G c = 0 holds. Otherwise, if G − G c < 0, the crack rests and its tip velocity v becomes 0. Taking these requirements into account, the relations are obtained. Based on this, the energy-dissipation functional is constructed which, according to Mielke and Ortiz [23], has to be minimized with respect to the displacement field u and the crack set C under the monotonicity constraint in Eq. (2). Herein, the potential energy U competes with the crack resistance G c |C|. In order to achieve independency from the spatial discretization in finite element simulations and to obtain a converging numerical scheme, the energy dissipation functional is regularized by now considering the -neighborhood C of the crack set C within the influence radius > 0. Note that Schmidt et al. [34] has proven the -convergence of the regularized energy-dissipation functional F to the unregularized energy dissipation functional F as → 0, so that mesh independence of the algorithm is ensured. For the application in finite element simulations, an algorithmic treatment is required. The main component therein is the erosion of an element K as soon as the potential energy in this element exceeds the maximally tolerable energy which would be released through the formation of additional crack surface. This can be expressed by considering the differences of the energetic states in element K after the element erosion and before, which are expressed by the -operator. Accordingly, we define U K = −U K to be the difference of potential energy in element K after eroding the element, which is zero, and the potential energy before eroding the element, which is equal to U K . Therefore, − U K corresponds to the potential energy in element K if the element was not eroded. In the context of Griffith-type fracture criteria, this energy has to be compared with the maximally tolerable energy released through eroding the element, which is described by the term G c A K . Herein, A K denotes the effective crack area, which describes the additional, regularized crack surface resulting from eroding the element. Thus, the net energy gain − F K can be introduced as the amount of potential energy which does not dissipate from element erosion and which needs to be stored in the remaining elements through an appropriate algorithm. This net energy gain is thus defined as the difference of the total potential energy U K = − U K , which would have been stored in the element if no crack erosion was detected, and the maximally tolerable energy G c A K , i.e., Thus, the net energy gain per effective crack area corresponds to − F K / A K = G − G c . The definition of the effective crack area is actually crucial ingredient of the eigenerosion approach and differs from previous nonlocal erosion algorithms. It is defined as and represents only that part of the neighborhood area which would additionally be considered if the element domain K was eroded as shown in Fig. 1. In a finite element implementation, the effective crack area is computed using the same Gaußquadrature which is also applied for the integration of the element residual and tangent stiffness matrices. For this purpose, the set of integration points is determined by first selecting the ones whose distance to the center point of the elements in C and K is lower than the influence radius . This list of Gaußpoints is reduced by the ones in C ( -neighborhood of the eroded elements C), which are in turn defined as the ones being within a distance of from the respective element center points. Then, the area A K is derived by adding up the volumes of the obtained set of Gaußpoints and normalizing by the factor 2 . Note that the choice of the influence radius influences the crack propagation; has to be large enough to ensure mesh independence. In our implementation, the algorithmic strategy consists of solving mechanical equilibrium for the displacements u and the regularized crack field |C| in an iterative, staggered scheme. After the mechanical equilibrium equations are solved in the current time step t n , the net energy gain − F K is evaluated for each element K . If the net energy gain − F K of an element is larger than zero, this element becomes a possible candidate for erosion in the current time step. The element with the largest net energy gain − F K > 0 is eroded, which means that its resulting forces and stiffness due to the material law is set to zero, and the mechanical fields are solved again without proceeding in time. This procedure is repeated until no further element to be eroded is found and then, the algorithm moves to the next time step. That way, in each time step, mechanical equilibrium is achieved for the updated set of eroded elements representing the crack. In the original approach, the -neighborhood is defined as a sphere around the center of mass of each element. For higher-order shape functions, e.g., using a quadratic approximation of the displacements, the mid-nodes may be located far from the straight connections of the corner nodes if the element formulation is not reduced to straight-edge elements. This may move the center of mass in such a way that rather small neighborhoods will not necessarily reflect non-local effects from all sides of the element, and thus, imply unnecessarily large neighborhoods. In order to a priori avoid this, we consider the neighborhood as union of spheres around the integration points rather than the center of mass. Furthermore, we incorporate inertia effects using the Newmark scheme [25] in order to allow for the simulation of complete fracture and detaching fragments. As a result therefrom, the erosion of elements will imply a loss in mass. To keep its effect on the dynamics small, we distribute the mass in an uncoupled manner to the corner nodes of the eroded elements. Additionally, if an element is eroded, the list of Gaußpoints lying in the additional crack area (C ∪ K ) \C for each intact element is updated because it might be reduced due to the increased crack area. The full algorithm is shown in Fig. 2.

Formulation for brittle fracture
The crack propagation in the case of brittle materials has been investigated for the case of linearized kinematics in [27,28]. Here, we extend this approach to the geometrically nonlinear framework of finite deformations. Note that although brittle fracture is defined to occur after almost purely elastic deformations and thus, usually at small strains, a geometrically nonlinear setting may still be reasonable. Specific engineering problems may result in boundary value problems where large deformations occur although only small strains are observed. In these cases, a fully geometrically nonlinear framework a priori takes into account the additional internal forces resulting from considering equilibrium in the deformed configuration. This is in contrast to the geometrically linearized setting, where equilibrium in the undeformed configuration is considered. Note that in the following descriptions, we will solely focus on the differences compared to the linearized setting itself, which are relevant to the components of the eigenerosion algorithm, rather than explaining the whole geometrically nonlinear setting. For details with this regard, the reader may consider standard literature, e.g., Wriggers [42]. Despite the fact that constitutive equations may still be formulated under the assumption of small strains to be used within the finite deformation setting, we consider the general concept of finite strains to obtain a unified formulation more suitably connecting to the extension of ductile fracture in the following subsection. The constitutive equations are formulated in the left Cauchy-Green tensor, i.e., its isochoric part respectively. Herein, F with J = det(F) denotes the deformation gradient, which maps infinitesimal vector elements of the undeformed reference configuration to their counterparts in the deformed, actual configuration. For our numerical analysis presented later, the material law for Neo-Hookean solids according to Rivlin [32] is taken into account, where the volume-specific strain energy density is considered. Following the second law of thermodynamics and the standard argument of rational continuum mechanics, the Kirchhoff stress tensor is obtained by deriving the energy density with respect to the deformation tensor, i.e., Herein, the deviator operator is defined as dev(b) :=b − 1/3tr(b)I with I being the second-order identity tensor. In order to evaluate the net energy gain − F K for each element K , Pandolfi and Ortiz [27] proposed to derive the imposed energy − U K = 0.5 d T K k K d K directly from the element stiffness matrix k K and discrete element displacement vector d K , which is only directly possible for linear elasticity. For an extension to more general cases, we propose to compute the volume integral of the internal energy ψ over the element region K in order to capture the whole work that has been imposed in element K . For small strains, this approach coincides with the original approach. In our case, the energy associated with inertia is not involved in − U K because it is not dissipated due to crack propagation.

Formulation for ductile fracture
For the extension of the eigenerosion algorithm to ductile crack propagation, the rate-independent finite strain J 2 elasto-plasticity formulation with isotropic hardening from [37][38][39] is applied. For the numerical implementation, see Miehe et al. [19] and Klinkel [13]. In this material model, the deformation gradient is multiplicatively decomposed as into an elastic and a plastic part, F e and F p , respectively, cf. also Kröner [14,15]. Following this, the elastic left Cauchy-Green tensor and its spectral decomposition result in with the square root of the principal elastic stretches λ e i obtained as eigenvalues of b e , and their principal directions expressed by the eigenvectors n i . With the application of the principal logarithmic elastic strains e i = log(λ e i ), the elastic strain tensor is obtained in logarithmic representation so that the additive split of the strain ε = ε e + ε p into an elastic and plastic part ε e and ε p can be considered. The strain energy density is divided into an elastic part ψ e (ε e ) only depending on the elastic strains ε e and the plastic part ψ p (α) depending on the equivalent plastic strain α serving as internal variable. For the elastic part, the quadratic elastic strain energy density with the compression modulus κ and the shear modulus μ is taken into account. The Kirchhoff stress tensor can then be identified as The convex plastic dissipation potential is chosen as the superposition of linear-and exponential-type hardening, cf. Voce [40]. Herein, the hardening parameters have the following meaning: h exp specifies the degree of exponential hardening, h lin defines the slope of the superimposed linear hardening, y 0 denotes the initial yield strength, and y ∞ describes the plastic yield strength at the transition to the almost purely linear hardening. Hence, we obtain the hardening function We consider plastic incompressibility, which is typical for metal plasticity, so that the flow condition becomes the von Mises type Together with the plastic variable λ p ≥ 0, the Kuhn-Tucker condition is fulfilled, where λ p is connected to the evolution of the norm of plastic strains ||ε p || = λ p (22) and to the evolution of the equivalent plastic strainṡ It is numerically computed from φ = 0, see Eq. (20), using the Newton scheme. In the case of ductile crack propagation, the imposed energy is proposed to be which represents an extension of the elastic strain energy density by the plastic potential ψ p . This approach is conceptually in line with the idea in Irwin [9] who extended the Griffith energy release rate for brittle materials to a Griffith-type energy release rate for ductile materials by adding the energy that dissipates into plasticity near the crack tip before the crack propagates. Therefore, the Griffith-type energy release rate G c has to be large enough, in order to trigger plastic deformations before the crack propagation. Note that this parameter G c for ductile crack propagation should not be confused with the classical Griffith constant, which is connected to brittle fracture.

Regularization of plasticity
In order to avoid mesh dependence, which may already occur during the formation of plastic shear bands before fracture occurs, regularization is applied. Here, the regularization is achieved by extending the elasto-plastic material model to elasto-viscoplasticity of the Perzyna-type, cf. e.g., Perzyna [29], Junker et al. [10]. Since here the rate dependency of the material is small and not to be modeled, the viscous part is only introduced for numerical reasons. Instead of computing λ p from φ = 0, we assume the plastic parameter to be computed from with the Macaulay bracket (•) ± = ((•) ± |(•)|)/2 and the viscosity η. Here, the same evolution Eqs. (22) and (23) as for the previous elasto-plastic model are considered. Thereby, the plastic deformation becomes more diffuse because in contrast to the unregularized model, the plastic deformation is delayed in the affected elements so that the neighboring elements are subjected to higher loading and have to deform plastically as well. Note that Eq. (25) may lead to a violation of the yield criterion, i.e., φ > 0, if the viscosity η is chosen large enough. Hence, the parameter η controls the speed of plastic deformation. In order to adjust the thickness of the shear bands and thus, the diffusive character of the regularized model, the hardening function is modified to with the localization parameter reducing the yield stress with increasing equivalent plastic strain α. This is necessary to imply additional plastic deformation after plasticity once has taken place in the particular material point, because the viscosity η strongly smears the plastic deformations over the body and therefore weakens the plasticity on the corresponding elements. Hence, without the parameter , no clear shear bands can develop within the viscoplastic material model if the viscosity η is chosen large enough to regularize.
With increasing values of the localization parameter , the thickness of the bands of plastic deformations is reduced independently from the spatial discretization except if the thickness of this band is smaller than the element size. Hence, a low localization parameter leads to an increase in the plastic deformations in the elements within the shear bands. Furthermore, the value of this rate-independent parameter has to be chosen in line with the selected viscosity η in each specific application. For this regularized material model, the resulting plastic dissipation potential is modified accordingly to Note that δ has to be small enough to keep the convexity of the plastic dissipation potential ψ p for the occurring range of the equivalent plastic strains α. Due to the viscosity, the criterion for ductile crack propagation is extended through the modified imposed energy which additionally takes into account the viscous dissipation D vis = ηα 2 in order to capture the full dissipation of the viscoplastic part of this material model. In order to achieve results comparable to experiments, all material parameters have to be fitted again considering the two additional material parameters η and δ. Note that the special type of viscous regularization considered here can be exchanged by any kind of regularization, including alternative viscous formulations, gradient enhanced formulations or energy-relaxed approaches. The regularization is only required in order to avoid mesh dependency even before the crack evolves. However, it is of course not sufficient to guarantee mesh independence during crack evolution and thus, a suitable crack propagation approach-here the eigenerosion algorithm-has to be considered.

Compression-tension asymmetry
For the simulation of cyclic fatigue, a previous crack may close again. Then, it is assumed that the two sides of a crack shall only transfer compression forces due to their contact and no tensile forces. To account for this, the formulation has to be modified in order to only partly erode elements, i.e., the tensile part of the energies, and keep the compressive part in the cracked elements. Therefore, the split of the elastic energy density into a tensile and a compressive part similar to Miehe et al. [20] and Ambati et al. [2] is considered as Here, the variable D is binary and indicates if an element is either intact (D = 0) or eroded (D = 1). This variable D should not be confused with continuous damage variables commonly used in continuum damage mechanics approaches. Considering this split, the tensile and compressive parts of the energy density become with the Lamé parameter λ = κ − 2/3 μ. Herein, the strain tensor is assumed to be split into as in Miehe et al. [20]. Accordingly, the Kirchhoff stress tensor results in For intact elements, we assume the elastic energy density ψ e = ψ e + + ψ e − , the stress τ = τ + + τ − , and the material tangent C = C + + C − needed for the Newton-Raphson scheme, respectively, while for eroded elements, only the compressive parts remain, i.e., ψ e = ψ e − , τ = τ − , and C = C − holds. Hence, the total potential energy to be considered in the fracture criterion becomes The resulting criterion for crack propagation ensures that elements would rather erode under tension loads and than under compression and consistently represents the energy that is dissipated into plasticity and crack propagation. Herein, the term D vis is set to zero for elasto-plastic materials. Note that in eroded elements, no further plastic deformation is considered in order to avoid an unphysical response under further tension loads. The algorithmic implementation of the material law itself including the material tangent C follows Miehe and Lambrecht [18], Miehe et al. [20].

Numerical examples
In order to show the performance of the proposed approach for the simulation of crack propagation at finite strains, several numerical calculations are analyzed which are based on rather classical benchmark problems.
To enable a more directly connected investigation of brittle and ductile fracture, we focus our analysis on the crack propagation in a dogbone-shaped specimen under tension, cf. Miehe et al. [22] or Aldakheel [1], where the qualitative development of the crack path is well-known from experiments. There, an already qualitatively differing crack path is observed for brittle and ductile fracture. The analysis is concentrated on showing a meshindependent response in the crack path as well as in integral quantities such as the external load-displacement curves.
(a) (b) (c) Fig. 3 a Specimen for tension test with initial crack at time t 0 , b specimen for tension test without crack, and c considered geometric parameters; T denotes the thickness of the specimens

Brittle fracture
First, the performance of the algorithm is analyzed for brittle fracture. Two different scenarios are investigated: (i) a specimen with an initial crack and (ii) a specimen without an initial crack, where a small perturbation close to the center is considered. Both specimens are subject to an external displacement and then, the response under crack propagation is analyzed. Figure 3 shows a schematic illustration of the numerical tests and the chosen geometry parameters.

Plate with initial crack
First, the specimen as shown in Fig. 3a with a prescribed initial crack with the length a at the time t = 0 s is loaded with the constantly increasing displacementū = t · 0.05 mm/s. The specimen is discretized by 27-node hexahedral elements with quadratic ansatz functions with the different characteristic element sizes h = 0.22, 0.11, 0.065, 0.0325 mm. Note that the element size h represents the minimum in-plane edge length appearing in the specimen. In thickness direction, only one layer of elements is considered. The material parameters are chosen as κ = 164,883.0 MPa, μ = 76,100.0 MPa, G c = 2 N/mm and the density ρ = 7810 kg/m 3 . The resulting final crack path after running the simulation is shown in Fig. 4. As can be seen, the crack follows the expected, horizontal path and it does not depend on the chosen discretization. Admittedly, the path also follows the element edge lines and thus, the positive results may have been influenced by already arranging the mesh in line with the expected crack path. However, as shown later, reasonable cracks will also evolve in directions differing from the element alignment. In order to quantitatively assess the quality of the simulation results, Fig. 5a shows the resultant force at the Dirichlet boundary versus the applied boundary displacements for different initial crack lengths. As can be seen, for all scenarios a mesh-independent response is obtained. This confirms at least numerically mesh independency also for the proposed eigenerosion algorithm for finite deformations.

Plate without initial crack
In contrast to the previous example, where the crack propagation direction could have been influenced by the alignment of the initial crack, now a specimen without initial crack is considered, see Fig. 3b. Because in reality, failure initialization is governed by processes on lower scales, it cannot directly be covered in the To avoid this inconsistency, the specimen is slightly weakened in an area with the horizontal length of (W − 2d)/10 close to the center of the plate. On purpose, this perturbation is not perfectly centered to not artificially symmetrize the problem. Thereby, in principle, the crack may evolve in arbitrary directions from here. The perturbation is included by setting the critical Griffith-type energy release rate there to the small value G c = 0.0001 N/mm. Of course, as usual for brittle fracture, the crack is still expected to evolve in horizontal direction starting from the perturbed area. The simulation results show that the crack propagates reasonably in horizontal direction, see Fig. 6, no matter the considered discretization. Also the response of the load-displacement curve shown in Fig. 5b shows mesh independence.

Ductile crack propagation
In a second set of numerical experiments, ductile crack propagation is analyzed. For this purpose, the same virtual specimen as above without initial crack, but with weakened area in the center, is considered, cf. Fig. 3b. The material parameters shown in the table of Fig. 7b, which were adjusted to the stress-strain curve of a tension experiment of an iron alloy as shown in Fig. 7a, are applied. The critical Griffith-type energy release rate is chosen to be G c = 11 N/mm.

Elasto-plastic specimen
First, the unregularized elasto-plastic material formulation is considered. As a result of the simulation, for a prescribed, external displacement ofū = 0.2 mm, the crack has evolved according to the illustrations in Fig. 8. There, also the distribution of equivalent plastic strains is shown at this loading state for different element sizes. As it can be seen, localized shear bands appear diagonally in a 45 • angle as expected from the analytic work of Onat and Prager [26]. The crack occurs within these bands because here, the elements' imposed energy − U K is increased by the plastic potential ψ p in contrast to the elements without plastic deformation. In those, the imposed energy only consists of the elastic part ψ e and thus, is lower if the plastic deformation is large enough. With decreasing element size h, the thickness of the shear bands decreases because only one row of elements deforms plastically. When these elements start deforming plastically, their stiffness is reduced so that any further plastic deformation only occurs in those because their resistance against deformation is lower than the ones of the surrounding elements. This leads to different times t c at which crack propagation occurs. Hence, the localization of the plastic deformation influences the evolution of the crack and thus, mesh dependency is obtained, as can be seen in the mesh-dependent force-displacement curves in Fig. 9. However, despite this mesh-dependent quantitative response, the crack path follows the expected 45 • angle and appears to be rather mesh-independent, see Fig. 8. Thereby, it is shown that a reasonable crack evolution can also be obtained in discretizations where the element edges do not follow the crack path orientation. In contrast to the brittle material, the crack does not completely propagate within one time step until the two sides of the specimen are split, and a rather delayed crack evolution occurs. Because of this, finite jumps in the reaction force are observable in Fig. 9. Furthermore, the field of the equivalent plastic strains α and thereby the crack path loose their horizontal symmetry as shown in Fig. 8 due to the sensitivity of the plastic localization to the location of the weakened area, which is considered slightly shifted from the center of the specimen (Fig. 10).

Ductile crack propagation with regularization
In order to show qualitative as well as quantitative mesh independence of our erosion algorithm, the pathological localization of the plastic deformation needs to be cured such that the crack propagation will not be influenced by erroneous distributions of plasticity. Therefore, now the regularized elasto-plastic material formulation is considered. For the regularization, the viscosity η = 300 MPa s and the localization parameter δ = 0.2 are chosen. In contrast to the virtual tension tests of the unregularized elasto-plastic specimen, the resulting equivalent plastic strain α before crack evolution shown in Fig. 11 develops diagonally in both possible directions with a higher intensity in the diagonal from upper left to lower right. This effect results from the regularization which delays the plastic deformation and thus enables additional plastic zones. Furthermore, due to the regularization, the shear bands converge with decreasing element size h. Hence, the crack propagation behaves qualitatively mesh-independent as shown in the crack paths in Fig. 12 and quantitatively as observed  Fig. 13. From these results, we conclude that the viscous regularization if already sufficient to cure the mesh dependence of the localizing plastic fields, which enables the erosion algorithm to keep its mesh-independence properties.

Influence of regularization intensity
The regularization of the plastic deformations appears to be crucial to obtain a decent starting basis for the crack propagation. Therefore, here, the influence of regularization intensity is analyzed in more detail. For this purpose, we firstly consider the response at the material point level. We vary the localization parameter δ in order to separate the effects of both regularization parameters. Figure 14a demonstrates the effect on the hardening function. As can be seen, the slope of the hardening function is decreased with increasing δ. Due to this, the plastic dissipation potential ψ p , c.f. Figure 14b, is decreased as well. Second, the structural response is analyzed. For this purpose, the viscosity η and the localization parameter δ are varied for the weakened specimen and a homogeneous specimen where not even the perturbation is applied. Figures 16 and 17 show the simulation results for an external displacement ofū = 0.15 mm. As it can be observed, for a low viscosity, i.e., η < 1000 MPa s, diagonal shear bands appear, while for a high viscosity, i.e., η ≥ 1000 MPa s, rather smeared fields are obtained. Furthermore, the results of the simulations with very small viscosity η = 30 MPa s become even asymmetric and only show single bands due to the high impact of plastic localization. In contrast to this, an increase in the localization parameter δ leads to an increased localization and thus, to a decrease in the thickness of the shear band. However, as shown in Fig. 14a, the localization parameter δ = 0.05 only changes the results slightly. Therefore, simulations using values smaller than this can be interpreted as analysis where the localization term is not considered. Hence, the form of the localized plastic strains can be controlled by adjusting these two parameters. The force-displacement curves for the different values of the viscosity η and the localization parameter δ are shown in Fig. 15. The results show that the reaction force F increases with decreasing values of η and increasing values of δ. Especially for δ = 5 and δ = 20, a rapid decrease in the force F becomes visible because the plastic deformation develops quite fast compared to the one of the simulations with a lower localization parameter δ (Figs. 16, 17).

Simulation with cyclic loading
Ductile fracture comes often along with a sub-critical crack propagation in the sense that the crack does not necessarily go through the specimen in almost an instance. This effect is important in the context of fatigue, where sub-critical crack growth happens as a result of increasing the number of load cycles. In order to show that  26) and (27) considering the elasto-plastic material model the erosion framework proposed here enables the simulation of fatigue based on purely macroscopic effects, simulations with cyclic loadings F(t) = F max sin(2 π t) and varying force amplitudes F max are carried out on the specimen depicted in Fig. 18a. In these simulations, the geometry parameters from Fig. 3c are used again and the coarsest mesh from the previous simulations including the weakened area is considered. Note that out-of-plane displacements on the backside of the specimen are prohibited in order to avoid bending when the specimen is loaded under compression. Here, the Griffith-type energy release rate is varied and chosen as G c = 3, 6, 11, 22 N/mm; the material parameters from Fig. 7b are applied. Furthermore, we define structural failure as the situation where the specimen is completely separated into an upper and lower part. Note that the tension-compression split of the material formulations as described in Sect. 2.3.2 was used here to allow for a realistic response under cyclic loading. In Fig. 19, the force-displacement curve of one simulation is presented which shows the typical response due to plasticity. With every cycle, the plastic deformation grows until the imposed energy − U K in one element K surpasses the crack resistance so that the crack propagates. Later on, the crack propagates until the structure fails after a certain amount of cycles N c . If the load amplitude F max surpasses the maximum static force, the plate fails within one cycle. In order to be able to trigger cyclic failure after a certain amount of cycles, the maximum force amplitude F max has to be large enough to cause plastic deformation in at least one material point of the structure. Otherwise, no crack will develop. With increasing load amplitude F max , the number of cycles until structural failure N c decreases linearly in the double logarithmic space as shown in Fig. 18b. These curves show typical Wöhler curves containing the force amplitudes F max over the number of cycles until structural failure N c . The reasonable response observed in Fig. 18b is only achieved as a result of the relatively fine deviations from cycle to cycle shown in the zoomed illustration in Fig. 19b, which are hardly visible in Fig. 19a. Obviously, such results can only be obtained if a predictive, mesh-independent crack propagation algorithm as the one proposed is applied. Note that although here only low-cycle fatigue examples were analyzed, the algorithm would in principle also allow for high-cycle fatigue. However, then purely macroscopic simulations would not be sufficient since the small changes of crack propagation with increasing numbers of cycles would only happen at the microscale.

Conclusion
The aim of the work was to develop a mesh-independent framework which enables the simulation of brittle and ductile crack propagation at finite deformations. The framework was obtained by extending the eigenerosion approach from [27], which was originally proposed for brittle fracture at small deformations.
In numerical examples of brittle as well as ductile fracture, it was shown that reasonable crack propagation results were obtained which also converged for finer discretizations implying a mesh-independent approach. This was qualitatively shown in terms of comparing the obtained crack paths, as well as quantitatively, where the resulting force-displacement curves were compared. Specifically for the ductile material, well-known localization of plastic deformations, which evolve even before crack evolution starts, was observed. This localization was found to have a significant effect on subsequent crack growth as these localizations already initiate mesh dependence, even before cracks evolve. Whereas a rather moderate impact on mesh dependency regarding the shape of the crack path was observed, the impact on the quantitative mesh dependency in terms of force-displacement curves was significant. In order to enable an objective analysis of mesh independency of the proposed crack propagation algorithm, the localized plastic deformations were regularized using an algorithmically viscous formulation. Provided that regularized plastic deformations are taken into account, the proposed fracture algorithm was found mesh-independent, qualitatively as well as quantitatively. Furthermore, a specimen undergoing cyclic loading was simulated and typical Wöhler curves could be obtained as a result. In future works, we  Fig. 19 a Force-displacement curve for simulation with cyclic loading and amplitude F max = 790 N and critical Griffith-type energy release rate G c = 6 N/mm, and b zoomed-in illustration of the peak at the maximum of the implied force F plan to investigate the extended eigenerosion approach in hard metal or metal matrix microstructures. There, crack propagation through brittle as well as ductile phases has to be simulated where natural benefits of our unified approach can be exploited. Based on such simulations, the material design may be improved to provide new engineering tools with higher wear-resistance, which are characterized by an increased toughness against sub-critical crack propagation at the microscale and thus, an improved resistance against surface spalling.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.