ProCrackPlast: a finite element tool to simulate 3D fatigue crack growth under large plastic deformations

Many structural components and devices in combustion and automotive engineering undergo highly intensive cyclic thermal and mechanical loading during their operation, which leads to low cycle (LCF) or thermomechanical (TMF) fatigue crack growth. This behavior is often characterized by large scale plastic deformations and creep around the crack, so that concepts of linear-elastic fracture mechanics fail. The finite element software ProCrackPlast has been developed at TU Bergakademie Freiberg for the automated simulation of fatigue crack growth in arbitrarily loaded three-dimensional components with large scale plastic deformations, in particular under cyclic thermomechanical loading. ProCrackPlast consists of a bundle of Python routines, which manage finite element pre-processing, crack analysis, and post-processing in combination with the commercial software Abaqus . ProCrackPlast is based on a crack growth procedure which adaptively updates the crack size in finite increments. Crack growth is controlled by the cyclic crack tip opening displacement Δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta $$\end{document}CTOD, which is considered as the appropriate fracture-mechanical parameter in case of large scale yielding. The three-dimensional Δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta $$\end{document}CTOD concept and its effective numerical calculation by means of special crack-tip elements are introduced at first. Next, the program structure, the underlying numerical algorithms and calculation schemes of ProCrackPlast are outlined in detail, which capture the plastic deformation history along with the moving crack. In all simulations, a viscoplastic cyclic material law is used within a large strain setting. The numerical performance of this software is studied for a single edge notch tension (SENT) specimen under isothermal cyclic loading and compared to common finite element techniques for fatigue crack simulation. The capability of this software is featured in two application examples showing crack growth under mixed-mode LCF and TMF in a typical austenite cast steel, Ni-Resist. In combination with a crack growth law identified in terms of Δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta $$\end{document}CTOD for a specific material, the tool ProCrackPlast is able to predict the crack evolution in a 3D component for a given thermomechanical loading scenario.


Introduction
Nowadays, hot parts in combustion and automotive components like vessels, piping, turbo-chargers, exhaust pipes, etc., undergo highly intensive cyclic thermal and mechanical loading during their service phase, which leads to isothermal (LCF) or thermomechanical (TMF) fatigue crack growth. This failure behavior is accompanied by large scale plastic defor-mations and creep around the crack. Thus, established finite element (FEM) or boundary element simulation tools for 3D fatigue crack growth such as FRANC3D (Fracture 2016), ZenCrack (Zentech 2019), AdapCrack (Schöllmann et al. 2003), CrackTracer (Dhondt 2016), ProCrack (Rabold et al. 2013;Rabold and Kuna 2014) or BEASY (CM 2011), are not suitable, since they are based on linear elastic fracture mechanics. Therefore, the current procedure in practice with detected defects often consists of a prophylactic, cost-intensive replacement of the entire part. To resolve this issue, there is high demand for a pragmatic FEM-based crack growth simulation software to predict and assess lifetime and safety of three-dimensional structural components under complex thermomechanical loading conditions, if the intensity and extend of inelastically deformed zones cannot be neglected. This is the motivation for the presented research results.
Among the numerical simulation tools, the so-called Extended Finite Element Method (X-FEM) (Fries and Belytschko 2010;Moës et al. 2002) has been established in the last decades, which avoids the repeated remeshing of the structure by extending the shape functions of the cracked elements. The additional shape functions reflect the displacement jump across the crack faces and the singular crack tip fields. However, the method is limited to linear elasticity. Despite its attractions, applications of X-FEM to large scale yielding (LSY) fracture problems or LCF/TMF are difficult, since the required local crack-tip fields for enhancement are not known. Therefore, the X-FEM does not allow to calculate the fracture parameter CTOD with required accuracy.
In recent years, the phase field method has been developed to address fracture mechanics problems as well, see e. g. Ambati et al. (2015) and Keip et al. (2016). This approach considers the damage state of the material as an order parameter and is able to simulate crack propagation by minimizing an energy functional. However, this technique requires extremely fine meshes and needs further improvement to solve complex engineering problems as 3D LCF/TMF in irreversible materials.
For the assessment of high cycle fatigue crack growth, the cyclic stress intensity concept has been established for decades in the frame of linear elastic fracture mechanics (LEFM). In case of low cycle fatigue (LCF) with plasticity at room temperature and/or creep at high homologous temperatures, global fracture mechanics parameters J or C * are traditionally used or derived quantities. They can be determined from path-independent contour integrals, and most commercial FE programs offer this function. The disadvantage is that the validity of path independence, which is necessary for a reliable numerical evaluation, imposes very restrictive conditions on the material laws and on the loading scenario (e.g. monotonic loading, steady-state creep for the case of C * ), which are not fulfilled for large scale yielding (LSY) and under TMF conditions. The application of the cyclic J-integral ΔJ as advocated by Beesley et al. (2015), Ochensberger and Kolednik (2015), Dowling and Iyyer (1987), Tanaka (1983) and Wüthrich (1982) is also problematic since there is usually no global, locationindependent load reversal point for a thermally stressed component, which would be necessary to determine ΔJ . Its numerical calculation (Muhamad Azmi et al. 2017;Vormwald 2016) requires a relatively elaborate evaluation methodology in 2D and is hardly applicable to creep fatigue and TMF. For 2D cracks under mode I, Metzger et al. (2015) have shown that ΔJ and the cyclic crack tip opening ΔCTOD are equivalent quantities for describing crack propagation in isothermal LCF.
An alternative to the above mentioned global fracture parameters is the application of local loading parameters like the crack-tip opening displacement CTOD, which retain their physical validity under more general conditions. Numerous experimental studies (Laird and Smith 1962;Pelloux 1970;Chowdhury and Sehitoglu 2016) and numerical simulations of cyclic crack tip plasticity (Pippan et al. 2011;Tvergaard 2004) have shown that crack propagation in ductile metals under fatigue loading results directly from the irreversible deformation processes that occur at the crack tip. The ΔCTOD is a measure of these local plastic deformations. Many investigations at room temperature (Pelloux 1970;Ktari et al. 2014;Antunes et al. 2017;Vasco-Olmo et al. 2020) and at high temperatures (Kiyak et al. 2008;Schweizer 2013) have found a correlation between the measured crack propagation rate da/dN and the experimentally or numerically determined ΔCTOD in the following form: Although the ΔCTOD concept has a convincing physical meaning, its application to fatigue crack growth (CG) in engineering structures is rarely found in the literature. One reason for this lies in the difficulty to determine the CTOD by experimental measurements or by numerical computations with sufficient accuracy. Please note that the value of ΔCTOD has the same order of magnitude as the crack growth rate per cycle itself, i.e. values of a few micrometers, which even may dramatically vary during crack growth. In last decades, the digital image correlation technique (DIC) has been advanced so far that the CTOD can be measured with high resolution on a propagating crack, see exemplarily Vasco-Olmo et al. 2019). Unfortunately, the method is restricted to the surface of the body and primarily to laboratory samples.
The calculation of CTOD requires very fine FE discretizations and a flexible control of the element size at the crack tip in dependence on the solution variables. Therefore, most numerical studies in literature consider only 2D mode I problems and use a fine mesh of regular isoparametric elements along the prospective straight crack path. The crack propagation is numerically driven by the node release technique, which is controlled by a global quantity like ΔK eff and not a local physical variable. Another ambiguous issue is the proper definition of CTOD. Most researchers determine CTOD from the vertical displacement of the first node (Antunes et al. 2017(Antunes et al. , 2018aTinoco et al. 2019) or the second node (Pommier 2002) behind the crack tip, or use the 45 •secant intersection method with the crack profile, see e. g. Kiyak et al. (2008) and Kuna (2013). Many numerical issues have been addressed in numerous investigations of this kind, as there are: convergence studies regarding the influence of mesh refinement, required numbers of load cycles to stabilize solution, influence of crack face contact (called plasticity induced crack closure), etc. Here, we refer to a selection of representative papers by Solanki et al. (2004), Antunes et al. (2017Antunes et al. ( , 2018a, Cochran et al. (2011) andJiang et al. (2005).
Real 3D simulations of fatigue crack growth are relatively rare due to their high numerical effort. In a series of papers, Camas et al. (2019Camas et al. ( , 2020 studied the influence of the crack growth scheme on fatigue crack closure in CT aluminum specimens using the node release technique. However, node release technique is bound to the fixed size of elements placed along the entire crack path, which drives up the numerical expense. Moreover, any change or prediction of crack path direction is not possible with such simulations. More flexibility is achieved by a local re-meshing procedure during crack growth, whereby a fine mesh is only con-centrated around the crack tip and moved with it. Zuo et al. (2004) were the first researchers, who implemented such a FE-technique for crack growth analyses in an elastic-plastic material. Their code CRACK3D was used by Lan et al. (2007) and Sutton et al. (2005) to simulate ductile crack growth under monotonic loading in a thin aluminium specimen using a CTOD criterion. The group of Vormwald (Rossetti et al. 2012) applied the remeshing and mapping procedure in combination with Abaqus to simulate mode I fatigue crack growth using the ΔJ -integral concept.
Another very important aspect of the modeling of ΔCTOD is the realistic choice of the constitutive law. Nowadays it is undoubtedly accepted that advanced nonlinear kinematic hardening models are needed to capture properly the cyclic plastic behavior including strain ratchetting, which affect the evolution of cyclic plastic zone and crack opening.
Please note that nearly all FE simulations of fatigue crack growth rely on concepts of LEFM and are restricted to small scale yielding. An extension to large scale yielding and TMF has been published recently by the present authors (Gesell et al. 2023a, b). Extensive 2D analyses of single edge notch tension (SENT) specimens revealed that collapsed special crack tip elements are superior compared with commonly used regular quadrilateral 8-node elements. At the same level of accuracy of ΔCTOD, they require an about ten times coarser mesh and show less sensitivity w.r.t. element size for both stationary and propagating cracks.
The present work aims to develop a numerical procedure for fatigue crack growth in 3D structures, which is controlled by ΔCTOD under conditions of large scale yielding. Compared to most previous works, progress is attempted w.r.t. the following points: -The crack growth (CG) rate is formulated as function of ΔCTOD. An extended formulation for 3D mixed-mode situations is proposed. -CTOD and ΔCTOD are calculated by using special crack-tip collapsed elements, which lead to reduced number of elements, more efficiency and robustness. -The magnitude of crack extension Δa is directly and adaptively controlled by the ΔCTOD itself in combination with a CG law. Thus, referring to any other global loading parameters is avoided, as done in many crack growth schemes in the literature.
-After each crack extension step Δa, a remeshing and solution mapping are performed, i. e. there are no restrictions regarding size and amount of viscoplastic deformations. The size of the elements at the crack tip is adaptively controlled by CTOD .
The paper is organized as follows: The CTOD concept is recapitulated as basic fracture mechanical parameter in case of large elastic-plastic and/or creep deformations around the crack. In particular, the extension to three-dimensional and mixed-mode crack problems is introduced in Sect. 2. For the numerical determination of ΔCTOD along a spatial crack front, a new method based on collapsed isoparametric hexahedral elements (CTE) is established in Sect. 3. Its performance is compared with the up-to-date common standard, where regular hexahedral meshes are used to simulate fatigue crack propagation, using the example of a Single Edge Notch Tension (SENT) specimen. The main body of the paper in Sect. 4 gives insight into the structure of the developed software tool ProCrackPlast. It addresses the implemented fully automated incremental crack extension technique to carry on the plastic deformation history with the moving crack. Recommendations are suggested regarding the influence of numerical parameters like element size, crack extension etc. Finally, two application examples are presented in Sect. 5 to demonstrate the capability and performance of ProCrackPlast for 3D mixedmode crack configuration and a TMF crack problem.

State of the art
Nearly all experimental and numerical investigations on CTOD-concept presented in the literature for fatigue crack propagation are restricted to small scale yielding (SSY) and mode I loading only. Consequently, the identified crack propagation laws are exclusively based on mode I behavior. However, in the general case of a complex external loading and/or a curved crack surface in space, mixed loading arises along the crack front composed of modes I (tension), II (in-plane shear) and III (anti-plane shear = torsion). Therefore, for each point of the crack front, it is necessary to decide in which direction the crack propagates and how fast it propagates there. To answer this question, it is necessary to extend the existing CTOD-concept to the spatial mixedmode crack opening situation. An extensive literature search revealed that there are no satisfactory international studies and findings on this topic to date, see recent review by Wang et al. (2020). The idea of using the normal and shear components of the CTOD vector in mixed-mode loading goes back to Li (1989), who was able to use it to correlate both crack deflection and propagation rates with experiments in the context of small plastic zones (SSY) for 2D problems. Later, the group of Sutton et al. (2000aSutton et al. ( , 2000b) elaborated a CTOD-criterion for crack deflection and growth under monotonic mode I/II loading in plasticity and successfully validated it with experiments in 2D slices of aluminum. Similar 2D investigations and FE simulations on compact tension shear (CTS) specimens by Pirondi and Dalle Donne (2001) have confirmed that the CTOD-criterion can predict the direction and initiation value in ductile fracture. When mode I is dominant, normal stress rupture occurs in the direction of maximum circumferential stress, while a shear crack forms along the maximum slip plane when mode II loading prevails. Adaptive 3D elastic-plastic FE simulations were first presented in Lan et al. (2007), successfully recalculating the experiments of Sutton et al. (2000a) on ARCAN specimens and thus further substantiating the developed mixed-mode CTOD-concept. In the paper (Floros et al. 2019), different crack deflection criteria were analyzed with experiments from the literature using 2D cyclic elastic-plastic FE simulations, but exclusively for discrete crack lengths, without considering the growth process. The comparison showed very good results for the CTOD-criterion. In summary, there is no established generalized CTOD-concept for the assessment of crack growth under LCF/TMF loading for large plastic zones and (visco)plastic material behavior. Also, the extended CTOD-concept favored here has not yet been tested experimentally or numerically for 3D crack problems. A comprehensive solution of this problem is not the subject of the present computational approach. It would require elaborate 3D mixed-mode LCF/TMF experiments. Thus, in the absence of an experimental database, only heuristic, pragmatic concepts could be provided in this work, which should satisfy the following properties: i) Simple and robust numerical evaluation, ii) Transition to known criteria of LEFM in case the plastic zone becomes very small.  For a 3D crack the crack tip opening varies along the crack front. Figure 1 shows such a discrete point along the crack with the associated local coordinate system (x, y, z). The y-coordinate gives the normal of the crack surface and z the tangent to the crack front. The three modes of fracture are defined by the relative displacements of the upper (+) and lower (-) crack face w.r.t. the three coordinate axes. Mode I is the opening Δu y direction, mode II and mode III refer to the in-plane and out-of plane sliding Δu x and Δu z , respectively. Figure 2 shows the three components CTOD M of the displacement vector in the (x,y)-plane after deformation. Under cyclic loading they are a function of time t: The range ΔCTOD M is the difference between its maximum and minimum value during a stabilized load cycle. Therefore, it is usually taken from the ascending path of the last load cycle in Fig. 8, i. e. from point B* to point C: The advancement of each crack front point by an amount Δa is based on the crack growth law specific to the material. The direction of crack growth is controlled by a crack deflection criterion, giving the deflection angle φ.

Criteria for crack deflection
In LEFM, the criterion of maximum tangential stress (MTS) has been widely established to determine crack deflection under mixed-mode loading based on the well-known crack tip solution as a function of stress intensity factors K I and K II . In the considered generalized case of viscoplastic material behavior, comparable analytical solutions are lacking. Therefore, the angle φ of the maximal circumferential stress component σ θθ must be determined numerically from FE-analysis on a constant radius r around the crack tip in the local (r, θ)-plane: Looking at CTOD, an obvious assumption for crack deflection is that the crack will seek a direction where the mode I opening component attains a maximum, which follows also from the theory of local symmetry. Accordingly, the new crack orientation is perpendicular to the inclination angle φ of the CTOD-vector in the local (x, y)-plane as shown in Fig. 3. Thus, the deflection angle φ is calculated from the (x, y) components of CTOD at maximum loading as follows: Both the versions of crack deflection criteria are implemented in ProCrackPlast with MTS as the preferred one. As with most other 3D deflection definitions, an inclination towards the z-component would lead to geometrically discontinuous crack fronts by local twisting, and is excluded here.
2.4 Equivalent CTOD definition for 3D mixed-mode crack growth The transfer of crack propagation laws obtained from pure mode I tests to the 3D case of mixed-mode crack loading requires the formulation of an equivalent measure for CTOD and ΔCTOD. The question, then, is what value to put in a crack propagation law like Eq. (1) ? Which equivalent 3D value ΔCTOD eq causes the same crack propagation in mixed-mode case as in pure mode I ? In this work, three heuristic definitions of the equivalent ΔCTOD eq are proposed, which are composed of the three modes: Option 1: Option 2: Option 3: Option 1 assumes that only the mode I component determines crack growth and can also be set for pure mode I loading. Option 2 is based on the well-known fracture limit diagram of Richard et al. (2014) and defines ΔCTOD eq as the Euclidean norm with variable coefficients. An alternative energetic approach, proposed by Tanaka (1974) in LEFM, consists of the sumof-squares of K factors, which, converted to CTOD, corresponds to the sum of magnitude values. To allow the user of ProCrackPlast the greatest possible flexibility later on, all three variants were implemented with freely selectable coefficients α, β, γ to weight the impact of the individual modes I, II and III, respectively.

Numerical calculation of ΔCTOD along a 3D crack front
To calculate the ΔCTOD most accurately and efficiently, special FE-techniques are required. In Pro-CrackPlast the special technique, which has been successfully developed in Gesell et al. (2023a) for 2D problems, is extended to 3D crack configurations. This technique provides the crack tip blunting accurately with least mesh dependency and gives the three modes of ΔCTOD. It goes without saying that all simulations must be performed considering finite strains and rotations in order to properly capture the large distortions and the blunting of the crack tip. To avoid the interpenetration of crack faces, appropriate contact definitions must be introduced against each other in the general mixed-mode case of a complete model.

Special mesh design along crack front
Commonly, simple hexahedral finite elements HEX8 are used. They have the advantage of being easy to use and allow simplified crack propagation algorithms such as the node release technique. However, they require a very fine mesh and are therefore computationally intensive. In this work, 20-node elements HEX20 with quadratic shape functions and full integration scheme are considered more suitable to capture the high elasticplastic deformations at the crack tip, see Fig. 4a. A typical crack tip mesh with these elements along a crack front is shown in Fig. 5a. The new idea is to apply so-called collapsed elements, which have proved efficiently in ductile fracture analysis under monotonic loading, see Barsoum (1977) and Kuna (2013). Hereby, one face of a regular HEX20 element is degenerated into a single edge as shown in Fig. 4b, so that these nodes are on the same coordinates in the undeformed state. However, they can move independently when deformed. In this way, the shape functions of the element are distorted to produce a 1/r singularity in the strain fields along each radial direction and reflect a constant stress state, which corresponds to the crack tip singularity for an ideally plastic material. Of course, this is an approximation for a general elastic-viscoplastic material, which exhibits hardening and creep according to a power law leading presumably to asymptotic near-fields with a singularity between the elastic 1/ √ r and the 1/r type. The com-

Example problem single edge notch tension (SENT) specimen
At first, the stationary crack problem in a SENT specimen is analyzed under cyclic tensile loading. The specimen has a rectangular cross section of 4 × 18mm with a straight through initial crack of length a = 2.45 mm. This specimen was subject of intensive experimental LCF and TMF investigations in Gesell et al. (2023b). Here, only two load cases with different Rratio are considered at room temperature to compare the numerical performance of HEX20 and CTE elements, see Table 1. The schematic of the boundary conditions and applied load of the SENT specimen are shown in Fig. 7a. The specimen is cyclically loaded with a nominal stress σ max which acts on the clamping jaws of the testing machine. Therefore, in this area on the surface of the specimen, the displacements in the radial direction are set to zero. The vertical displacement component is fixed in the lower clamping and is specified at the upper clamping according to the applied stress. The crack experiences pure mode I under this loading and it grows along its initial plane. A typical FE-mesh of the SENT is shown in Fig. 7b. It is built as a small circular tube along the crack front with very fine discretization, whereas the remaining part of the body is automatically meshed with tetrahedral elements. The fine meshing along the crack-front region is constructed either with HEX20 or with CTE elements in the manner shown in detail in Fig. 5a and b, resp. The same size L e = 0.02 mm of elements around the crack front is chosen for both meshes.  Any cyclic elastic-(visco)plastic material law can be used in ProCrackPlast, either from the Abaqus family or as user-defined routine. Throughout this work, we use exemplarily the viscoplastic material model of Chaboche-type (Lemaitre and Chaboche 1990), which has been developed in previous research (Gesell et al. 2023b, a). It was formulated for large strains as an Abaqus UMAT to account for rate-dependent plastic deformations and creep occurring at the crack tip during high temperature loading. Realistic model parameters have been identified from test data of the austenitic cast iron EN-GJSA-XNiSiCr 35-5-2 (Ni-Resist-D5), which is commonly used for casing parts of gas turbines, exhaust manifolds and turbochargers (Kühn et al. 2017). Details of the material model are described in the appendix together with the specific parameters of Ni-Resist-D5 for the considered temperature of 20°and 700°.
In ProCrackPlast, any external thermomechanical loading of the structure with crack is allowed, which emulates a typical operation cycle during service (switch on/off, heating/cooling, etc.). Such a loading cycle is schematically depicted in Fig. 8, which is normalized by maximum load and may have any load ratio R = σ min /σ max . Before each crack growth step, this load cycle is multiply applied to ensure a stabilized stress-strain response at the crack tip, from where a reliable ΔCTOD can be measured. In agreement with experience in literature we found that two cycles are necessary and sufficient in this respect, see path A-B-C.  To illustrate the performance of the program for large scale yielding, the size of the plastic zone across the SENT specimen is shown in Fig. 9 for the load case RT0. The color plot represents the accumulated inelastic (at RT mainly plastic) strain. All values in gray lie above the initial yield strength of the material, i. e. the plastic region has reached the remote boundary (fully plastic yielding). The figure is a snap shot after four increments of crack growth. One can see the plastic wake behind the current crack tip.

Definition and determination of CTOD and ΔCTOD
These crack opening parameters are determined during the course of cyclic loading from the distance between a distinct pair of points, lying opposite on the upper and lower crack face in the undeformed configuration, cf. Fig. 2. The choice of these points has some arbitrariness. For hexahedral elements usually the first nodes on both crack faces behind the crack tip are chosen. In the case of HEX20 elements, these are the mid-edge nodes of the opposite first elements. When CTE elements are used, the opposite collapsed nodes of the first and last element (in the circumferential direction) of the rings around the crack tip are taken, see Fig. 6. For the symmetric SENT specimen, the mode I CTOD is simply 2 u y displacement of the selected points. Fig. 10 shows the obtained crack profiles in the center of the SENT specimen for both load cases at maximum and minimum load level. The points of CTOD definitions are marked. The cyclic ΔCTOD is obtained from the temporal course of the CTOD as the difference between the maximal and minimal value during the last ascending path B*-C of the load scheme in Fig. 8. The mathematical definition of ΔCTOD is also expressed by Eq. (3). Another visualization is given in the hysteresis loops of Fig. 11.
The crack profiles of Fig. 10 show a remarkable influence of the applied load ratio. For cyclic tensile loading R = 0, the crack does not close after remote unloading, which is hindered by the positive inelastic eigenstrains in the cyclic plastic zone. Contrary, the largest part of the crack is closed at maximal reverse load for the ratio R = −1. Only the first nodes behind the crack tip remain open. This is due to negative plastic strains at the very crack tip, which cause local tensile stresses and opening. This phenomenon has been observed as well by other researchers (Tvergaard 2004). The evolution of CTOD and ΔCTOD during cycling is better visualized in Fig. 11. The hysteresis loops show a pronounced ratchetting effect. However, the ΔCTOD values do not change much between first and second cycle, which indicates stabilization. The steep drop of CTOD for R = −1 is caused by the crack face contact.

Comparison of FE-techniques
A look on the crack profiles in Fig. 10 obtained with CTE and HEX20 meshing technique shows a rather rapid blunting / opening of the crack with CTE elements, whereas HEX20 elements need some distance to arrive at comparable CTOD level. If the first node behind the crack tip is used to determine CTOD for HEX20 elements, then a good agreement with CTE results is obtained. The second node would give much higher values. The absolute size of crack tip element L e = 0.02 mm must be seen in relation to the values of CTOD = 0.0036 mm (5× for RT0) and CTOD = 0.0015 (≈ 13× for RT-1). The hysteresis curves of both meshing types give comparable results and about the same ΔCTOD, provided the first node is used for HEX20 elements. To complete the picture, the distribution of CTOD and ΔCTOD is represented along the crack front of the SENT specimen as been calculated by HEX20 and CTE meshes for RT0. The agreement is quite good, see Fig. 12. The crack opening is maximal in the center and decreases towards the surface. Finally, the influence of the element size at crack tip was investigated. Instead of the previously used "fine" mesh, the same simulation RT0 was run with a five times coarser discretization having an element size of L e = 0.10 mm. The results are depicted in Fig. 13. The crack profile for the CTE elements keeps nearly unchanged, whereas the HEX20 profiles differ a lot, since the coarser elements cannot reproduce the blunting. Even the first node would yield too high CTOD values. The reason for the advantage of the CTE meshing is illustrated in Fig. 14, which shows the v.-Mises equivalent stress on both crack tip meshes deformed at maximum load for the RT0 case. The angular distribu-tion of stresses and plastic strains can be much better resolved by the fan-shaped CTE elements. In addition, the independently moving collapsed crack-tip nodes can better capture the blunting. More detailed comparative studies between regular and collapsed crack-tip elements have been reported for 2D in Gesell et al. (2023a), which have verified the higher robustness and efficiency of the CTE approach. Therefore, the CTE approach is also favored for 3D fatigue crack simulations, and implemented in ProCrackPlast.

PROCRACKPLAST software structure
Considering the different aspects of numerical crack growth simulation for large plastic deformation, presented in the previous sections, an incremental crack growth technique for 3D applications is developed. This is implemented in the crack growth simulation software, ProCrackPlast. ProCrackPlast has been created on the basis of ProCrack, a linear elastic crack growth simulation tool that was developed at TU Bergakademie Freiberg (Rabold et al. 2013;Kuna 2014, 2016;Ludwig et al. 2020). ProCrackPlast is operating in combination with the FE package Abaqus (Abaqus 2021) and consists of several program modules in Python for an automated pre-processing, FE analysis, and post-processing of the crack growth analysis. The first version of Pro-CrackPlast was developed for the Linux operating system and is compatible with Abaqus versions 2020 and 2021. For more details see the User Manual (Dude 2022).

Incremental crack growth algorithm
In ProCrackPlast, the growth of a crack in a 3D structure is simulated as a sequence of crack configurations called crack increments based on the associated crack growth law, which is formulated using the CTOD and ΔCTOD. The basic concept of Pro-CrackPlast is the sequential execution of the following steps for each crack increment: -Definition of the crack in the component at a geometric level for each increment. -Generation of a special FE mesh of the crack configuration. -FE analysis of the model for the applied loads and boundary conditions. The flow of process in ProCrackPlast is depicted in Fig. 15. The inputs required for the simulation include the Abaqus /CAE model of the crack-free component and a configuration file to define all the control parameters for the simulation. In the input CAE file, the component part, assembly, boundary conditions, load sequence and the required outputs should be defined. The most important control parameters in the configuration file include the definition of initial crack geometry, mesh parameters and crack growth law. All the parameters and their definitions are given in the User Manual (Dude 2022). The user must make efficient engineering judgments to set up these parameters in the configuration file prior to the simulation.
The software takes the input files and enters a loop as shown in the Fig. 15. Each cycle in this loop is an incremented crack configuration and the associated FE analysis. The main processes in the loop are shown in the orange boxes. The first process in the flowchart represents the generation of the FE model of the cracked component. In ProCrackPlast, the FE model of each crack configuration is generated as a Focused Global Model (FGM), the details of which are given in Sect. 4.2. This is done in the CAE module of the Abaqus package, which provides routines for generation of geometry and mesh. In the first crack increment loop, the initial crack is defined based on the associated parameters in the configuration file. For the subsequent crack growth increments, the current crack front is updated by extending it based on the used fracture mechanical concept.
The FE model is solved with Abaqus /Standard in the next step and the results are processed in steps three and four. In the post processing, the CTOD and the ΔCTOD for the three modes of crack tip displacements are calculated along the crack front. Furthermore, the crack increment and the crack deflection at each point of the crack front are determined. The details are given in Sect. 4.3.
The crack propagation is terminated when the equivalent CTOD value is below a threshold representing crack arrest or it is above a critical value which represents unstable crack growth. The parameters for the threshold and the critical value should be set in the configuration file. If the crack termination condition indicates further growth, the crack front coordinates for the incremented crack are calculated in the next step and the crack is updated in the next loop and the process continues. A distinct feature of this algorithm is the mapping of deformation history, the details of which are presented in Sect. 4.4.

Crack modeling in FGM
ProCrackPlast enables to model a 3D crack with a continuous crack front, which extends to the surface of the component. At the geometrical level, the crack is considered as a single 3D surface which is updated in each crack increment step. This crack surface is approx- Fig. 15 Crack growth algorithm in ProCrackPlast imated as a tessellation of geometric faces of quadrilaterals and triangles. The actual crack is generated at the FE level with this crack surface using the seam function in Abaqus /CAE. The initial crack in ProCrackPlast is a planar crack defined by the user in the configuration file. The crack front and the area of the initial crack are realized from the intersection of a user-defined ellipse and a line segment. For instance, a planar semi-elliptical surface crack in a cuboidal block (see Fig. 17) can be modeled as shown in Fig. 16. Defining this initial crack is illustrated here. First, a crack coordinate system (X,Y,Z) is considered such that the plane of the ellipse is parallel to the (X,Y)-plane of this coordinate system. Next, a line segment is defined, parallel to the X axis of the coordinate system. The location of this line segment can be adjusted with respect to the (X,Y,Z) coordinate system, to lie in the plane of the initial crack. Next, the center of the ellipse is defined and its semi-axes c and a, which lie parallel to the X and Y axes of the coordinate system, respectively. The ellipse intersects the block at two points, and this segment defines the crack front. The line segment and the elliptical segment span a planar surface. A separate geometry is generated from this, the crack surface part. The intersection of the crack surface part and the component defines the initial crack. The parameters to define the coordinate system, length of the segment and its offset, center and the dimensions of the ellipse are given in the configuration file.
The crack front is modeled as a polygon of straight line segments in ProCrackPlast. There are user input parameters in the configuration file to control the num- Fig. 16 Initial crack definition scheme in ProCrackPlast ber of crack front segments. The intersection of these segments form the crack front points where the fracture mechanical parameters are calculated. After each crack growth increment, a new set of crack front points is obtained relative to the current points as specified in Sect. 4.3. These points are combined with the old crack surface part to generate the new crack surface. The intersection of it and the component geometry gives the updated crack. The schematic of the crack in the cuboidal block under mode I loading after three increments is shown in Fig. 17.
Growing cracks often assume the shape of a varying curved crack front. Due to this feature, the extension of the crack front increases or decreases as the crack grows, depending on the convexity of the crack front. For example, in the surface crack growth given in Fig. 17, the crack front is elongated. The crack front segmentation is adapted in these conditions in Pro-CrackPlast. If the distance between two crack front points is too large, then the polygon segments do not approximate the crack front appropriately. Therefore, an additional crack front point is introduced in between, as soon as the distance exceeds a limit set by the user. Similarly, if the distance is shorter than a certain limit, then it is difficult to generate a valid HEX20 mesh in the focus region. Then, the point with the shortest distance to the neighbor is removed.

Mesh generation in FGM
A Focused Global Model is essentially a FE method to discretize the cracked component with a tubular focused mesh around the crack front embedded in a tetrahedral mesh of the component, see Fig. 18. Such a model is chosen due to the inability of the meshing tool to automatically generate structured elements with a focused mesh at the crack tip, especially for complex geometries.
To create the focused mesh region, a tubular geometric partition around the crack front is created in the component. To this end, the crack tip coordinate system in Fig. 1 is considered and the geometric profiles of the focused mesh are constructed along the crack front. The section profiles at each crack front point lie in the (x, y)-plane of the corresponding crack tip coordinate system and are oriented along the x axis conforming to the curvature of the crack front at that crack tip. A separate tubular shell part is created sweeping through these section profiles using the loft feature in Abaqus /CAE. The intersection of the tubular part with the component defines the geometric partition for the focused mesh region as shown in Fig. 17. In the tubular part, the end surfaces are extended to ensure an error-free intersection with the component.
The global region outside of the tubular partition in the FGM geometry is meshed with tetrahedral elements using the meshing function in Abaqus /CAE. The crack surface is defined as a seam before meshing and this feature generates the crack in the FE model. The seam crack has zero thickness, with the crack faces lying on top of each other. The region inside the tubular partition is meshed with hexahedral elements in the focused mesh pattern. The focused mesh is generated with Python functions. The nodal coordinates of the focused mesh are identified at each crack front point based on the crack tip coordinate system. The collapsed nodes lie at the origin of the crack tip coordinate system while the other nodes lie concentrically in the (x, y)plane of the coordinate system. This method ensures a very good conformance of the tetrahedral and the hexahedral element surfaces at the interface. However, the surface nodes of the focus part may not lie on the component surface unless the crack-tip plane (x, y) coincides with the component surface. In case they do not, the surface nodes are projected back onto the component surface. The HEX20 elements of the focused mesh are defined based on these nodes.
Four-noded (C3D4) or ten-noded (C3D10) tetrahedral elements can be used in the global region, while the focused mesh consists of twenty noded hexahedral (C3D20) elements alone. The displacement continuity on the interface between the dissimilar mesh types in Fig. 18 is enforced with suface-to-surface tie binding. Contact is applied between the FE crack surfaces, the definitions of which can be set up by the user in the configuration file. The user can also control the parameters for the FGM via the configuration file. The important ones include the radius of the tubular partition/focused region, the number of radial and circumferential elements in the focused mesh, the mesh density and corresponding bias at different regions of the FE model.

Propagation of the crack
The crack propagation scheme in ProCrackPlast basically integrates the given crack growth law and extends the crack considering either of the crack deflection criteria given in Sect. 2.3. It is designed to guarantee a continuous crack surface upon crack incrementation. As such, the crack extension is realized by extending the current crack front points only in the (x, y)-plane of the corresponding crack tip coordinate system as shown in Fig. 19.
In ProCrackPlast, the FE analysis is carried out for each crack increment imposing the applied load cycle. From the crack tip opening histories at the crack front points, the three modes of CTOD and ΔCTOD are identified for the corresponding point, as described in Sects. 2.2 and 3.3. Furthermore, the equivalent measures of CTOD eq and ΔCTOD eq are determined based on the user-defined option (see Sect. 2.4) from the three modes of CTOD and ΔCTOD.
The CTOD eq and ΔCTOD eq are used in the crack growth law to determine the crack extension, Δa i at each crack front point. Three specific forms of crack growth functions, g(ΔCTOD eq , CTOD eq ) based on CTOD /ΔCTOD are available in ProCrackPlast as given below. The user can further implement other types of crack growth laws based on the history of crack tip opening in ProCrackPlast, refer to Gesell et al. (2023b) for details. The crack growth per cycle, da/dN has the following generic form, where g(ΔCTOD eq ) is called crack growth function: da dN = g(ΔCTOD eq ).
(10) Fig. 18 Mesh details in the FGM of the block model Three variants of g(ΔCTOD eq , CTOD eq ) are available in ProCrackPlast: where C, m are material parameters.
where C, m, n are material parameters and R is the load ratio. C, n are temperature dependent.
GENERALIZED : g = C 1 (ΔCTOD eq ) m 1 +C 2 (T )(CTOD eq ) m 2 . (13) where, C 1 , C 2 , m 1 , m 2 are material parameters and C 2 is a temperature dependent coefficient. In the crack extension scheme of ProCrackPlast, the user needs to set a numerical parameter, Δa min , which dictates the minimal amount of crack advancement at each increment. To ensure an error-free generation of the Focused Global Model, Δa min is kept larger than the focus region partition radius R F . The user can input this parameter as an absolute value or a multiple of the CTOD of the previous increment. Δa min corresponds to the crack extension at the crack front point with the least equivalent CTOD along the crack front. It is further assumed that during this finite crack extension, the fracture parameters (CTOD, ΔCTOD) do not change along the crack front. With this assumption, the number of physical load cycles, ΔN , required for the Δa min in a crack increment step can be calculated  Figure 19 shows the crack extension, Δa i , and the crack deflection, φ i , at the crack front points. The numerical details of determining the crack deflection, φ i , with the MTS criterion Eq. (5) are omitted for brevity. This crack incrementation scheme sometimes results in a non-smooth crack front, due to numerical errors in the process. It has been observed that, if such perturbations occur, it worsens with further crack increments and sometimes results in aborting the simulation. To eliminate this effect, crack front smoothing schemes are available in ProCrackPlast, one based on linear regression and another based on spline interpolation. The smoothing option and the amount of smoothing can be set up through the configuration file. More details are given in the User Manual (Dude 2022).

Mapping of deformation and field variables to the new crack configuration
The inelastic deformation history effects are taken into account in ProCrackPlast by transferring the deformation and the internal state variables at each crack increment step. The mesh-to-mesh solution mapping feature in Abaqus is used to map the state variables of the viscoplastic material law, while a linear interpolation of the displacement field between the nodes of successive crack increment meshes is used for deformation mapping.

Mapping of state variables
In the mesh-to-mesh solution mapping algorithm, all the field variables at the nodes are directly interpolated between the meshes excluding the displacements. The Gauss point variables of the old mesh are first interpolated to the nodes, then averaged over the elements sharing a node, and finally interpolated with shape functions to the new mesh. This algorithm induces numerical diffusion to fields with large gradients, in particular when used between meshes with very different element size. Nonetheless, this method is used in ProCrack-Plast considering its simplicity and robust implementation in Abaqus . Hereby, it may happen that nonphysical, unacceptable values are assigned to the state variables. To overcome this problem, possibly necessary corrections to each state variable are performed before further analysis. These corrections were integrated into the UMAT subroutine (Gesell et al. 2023a).

Mapping of deformations
The boundary conditions and the load configuration of the new mesh, to which the state variables are to be mapped, have to match those of the old deformed mesh. To ensure this in ProCrackPlast, the mapping point is recommended to be chosen as 'B' in Fig. 8. With this choice, the mapping configuration 'B' consistently matches the initial configuration 'A' of the load cycle in the incremented crack growth step. Furthermore, the geometry of the new mesh for the incremented crack must match the deformed mesh in the current increment, except for the crack increment itself. This is especially relevant considering the residual opening along the crack faces and also for the better convergence of the mapped state variables. Considering the highly intricate deformation in the context of complex 3D cracks, the crack increment and remeshing are performed on the initial undeformed geometry (reference configuration) for every crack increment step, and the deformation history is sequentially added to the nodal coordinates of the corresponding mesh. The deformation mapping scheme in ProCrackPlast is depicted in Figs

Submodel analysis in ProCrackPlast
In industrial applications, the FE analysis of structures with a crack is most often restricted by the large size and complexity of the components including assemblies. Therefore, a submodel analysis is included in ProCrackPlast to extend its applicability to such large-scale problems. With this feature, a region around the crack and its potential growth zone can be given to ProCrackPlast as a submodel.  Kühn et al. (2017). This material is commonly used for the casing in gas turbines, exhaust manifolds and turbochargers. Since the crack growth law of this material is still under development, a first estimate is used in these simulations with a Paris type law Eq. (11). The parameters are chosen as m = 2 and C = 100/mm.

Crack growth in a three-point bending specimen
This example has been often used in literature to simulate fatigue crack growth under mixed-mode loading, but with elastic material behavior. Figure 22 shows the schematic of the three point bending specimen with boundary conditions and loading at room temperature. The specimen contains an initial straight surface crack of depth 0.6 mm, which is inclined at an angle of 60 • to the longitudinal edge of the specimen. A cyclic force of maximum 500 MPa (R=0) is applied at the center of the lower surface with two cycles as described in Fig. 8. On the upper surface, the bearings are modeled by displacement boundary conditions. Ten crack growth increments are simulated for this beam. The equivalent ΔCTOD eq for the crack growth law is chosen as option 2 in Sect. 2.4, i. e. Eq. (8) with α = β = γ = 1.0. The MTS criterion Eq. (5) is used for crack deflection. Figure 24 shows the evolution of the crack area in the CAE model (right) and the FGM mesh at the end of the simulation (left). It can be seen that the crack twists as it grows and becomes perpendicular to the longitudinal axis of the beam, i.e. perpendicular to the bending stresses. At the beginning, Fig. 22 Three point bending specimen with oblique crack. Dimensions are height = 6 mm, width = 5 mm, length = 28 mm the crack front in the middle of the beam is mainly stressed in mode I, while towards the side faces all three modes occur, see Fig. 23a. As the crack grows, the original mixed-mode loading evolves to a dominant mode I condition everywhere, see Fig. 23b. Figure 25 shows the equivalent inelastic strain contours in the beam after two crack growth increments. The crack tip at the mid-section experiences dominant mode I loading, while on the beam surface all the three modes occur, see the isolated focus mesh at the mid-section and the surface in Fig. 25.

TMF crack growth in a tube
TMF crack growth in a short tube under cyclic heating is presented in this example. This is a simplified structural model to mimic TMF crack growth occurring in turbo-charger hot parts. Figure 26 shows the schematic of the tube and the boundary conditions. The tube is 4 mm thick and contains an initial elliptical surface crack at the inner wall and aligned parallel to the axis. A section of the outer wall opposite to the crack side is fixed. The inner wall of the tube is cyclically heated between 200 • C and 700 • C by prescribing the temperature, see Fig. 27. Convective heat transfer conditions are applied on the outer wall with an ambient temperature of 200 • C. Otherwise, the component is mechanically load-free (no pressure). The crack experiences a pure mode I opening under this loading.
First, the heat transfer analysis is performed in the tube without crack, since the crack does not affect the heat flow. For this simple geometry a regular mesh is considered. The resulting temperature distribution in the tube at various instants of the heating cycle is shown in Fig. 28. The transient thermal fields from this simulation are used to drive the structural analysis of the tube (sequentially coupled thermo-mechanical analy-  . This is done by providing the *.odb-file for the thermal analysis as an additional input to ProCrack-Plast. The thermal gradients across the thickness of the tube induce a non-uniform stress state that causes the crack tip loading, see Fig. 29.
The uniform heating of the inner wall of the tube results in varied circumferential expansion across the thickness of the tube. As such, the region close to the inner wall expands more compared to the outer and this results in compression of the crack and the development of negative inelastic strains during the heating phase until point C. As the temperature is reduced down to point D, all thermal strains disappear. During this phase the elastic region pulls the embedded compressive inelastic zone back, resulting in tension at the crack tip and crack tip opening. The stress distribution in Fig. 29 along the radial ligament in front of the crack reflects these features. Further, the crack tip opening history resulting from this stress evolution, characterizes this problem as an out-of-phase TMF as evident from Fig. 30.
The crack growth through the thickness of the tube is documented in Fig. 31. The crack propagates more along the surface at the inner wall than into the depth of the tube, which is understandable due to the higher thermal tensile stresses at the surface. Further, the finite element mesh around the final crack configuration is also depicted in Fig. 31.

Conclusions
ProCrackPlast is a customized finite element software to simulate incremental crack growth in real 3D engineering components under large plastic (LCF) and thermomechanical (TMF) loading conditions. The special FE modelling technique used in ProCrackPlast, Focused Global Model, enables the automatic generation of structured meshes in complex crack geometries. The cyclic crack tip opening ΔCTOD is suggested as the appropriate fracture mechanical parameter to control fatigue crack propagation under large scale (visco)plastic deformations. For the accurate determination of ΔCTOD for a propagating fatigue crack, a special FE technique is developed. Here, for the first time, the effectiveness and robustness of collapsed crack-tip elements compared to regular hexahedral elements are demonstrated for 3D cracks. The size of these elements is automatically adapted to the magnitude of the CTOD and plastic zone size during the crack growth. The challenge of such simulations is to capture the history effects due to the irreversible inelastic deformations. This is accounted for in ProCrack- Plast by updating the deformation and by mapping the solution fields between the successive incremental FE meshes of the crack model. The use of any (visco)plastic constitutive law is possible and reasonable in combination with ProCrackPlast, provided it considers large strain setting and kinematic hardening. The two application examples show the potential of ProCrackPlast to predict the evolution of 3D cracks under mixed-mode and thermomechanical loading, when large inelastic deformations occur.
The main focus of the present work was on the development of suitable algorithms and numerical studies, whereby a simple crack growth law of Paris-type has been used. In cooperation with the Federal Institute of Materials Research and Testing (BAM), ΔCTODcontrolled crack propagation laws have been elaborated for the austenitic cast iron Ni-Resist (Gesell et al. 2023b). For this purpose, extensive isothermal LCF

Fig. 30
Temperature and CTOD history measured at the crack tip node in the mid-section of the initial crack and TMF experiments have been conducted for 20 • C -700 • C, and specific crack growth laws were identified by 2D numerical simulations. The results have shown that no unique relationship exists for all load ratios and temperatures. Therefore, three different groups of crack growth laws were specified, as presented here in Eqs. (7)-(9). These relationships can be used by others to match their experimental results. Further 3D simulations and experimental validation of ProCrack-Plast for this material are ongoing work.
The presented investigations have revealed a number of open issues regarding the application of the ΔCTOD concept to 3D mixed-mode crack configurations. A general extension of the approach and specific experimental investigations are needed and missing in the current literature. First, this concerns the formulation of an equivalent CTOD quantity that is representative for mixed-mode curved crack growth, and can be applied to generalize a mode I crack propagation law. Here, we suggest three different heuristic equivalent measures ΔCTOD eq , which however need to be validated by suitable experiments. Next, the criterion for the deflection angle of the crack under mixed-mode large scale plasticity requires further investigations. The geometric CTOD criterion suggested in Sect. 2.3 has proven to be dependent on the amount of plasticity and temperature. Therefore, the maximum tangential stress criterion is recommended for the situations at hand here.
Finally, it should be emphasized that software tools like ProCrackPlast form the prerequisite to analyze such crack problems and to calculate the crack tip parameters, i. e. to allow the interpretation of experimental results for improvement of the fracture mechanical concepts.