Yet Another Approach to Fatigue Crack Growth Simulation

: The analysis of a material that is subjected to variable loads is a complex subject and generally treated separately by fatigue and fracture mechanics. We present an attempt to extend the validity of conventional fatigue approach (here strain-life) in the scope fracture. This was achieved by introducing a zero thickness cohesive contact element coupled with a damage parameter that was developed from material observations of strain controlled fatigue experiments. The presented simulation framework results in a predictable crack growth direction on a compact tension specimen, although further experimentation is needed to validate the proposed approach.


Introduction
The lifespan of a cyclically loaded product can be divided up as follows. The product delivers the expected functionality for the majority of its useful life, without any noticeable changes on the macro scale. But if we look closer, on the scale of crystals and grain boundaries, microcracks emerge, join, coalesce and grow [1,2]. When microcracks reach a human-noticeable length scale, e.g., 1 mm, we refer to them as macrocracks. If left untreated, these cracks continue to grow, which leads to the product no longer functioning as intended and, finally, to complete rupture.
As this is a major safety issue, the fatigue analysis of products is very important. As we still lack general practical applications to describe the underlying micro-mechanical mechanism, we focus our attention on the macroscale. Here, a natural division is made between the safe-life (SL) [3] and the damage-tolerant (DT) [4] approaches. The SL approach is based on a product being safe to use up to the point of macrocrack nucleation. The DT approach, on the other hand, permits the existence of cracks and sets the limits of permissible damage. As the two approaches cover different aspects of fatigue life, their underlying approaches are inherently different. Typically, exponential decay functions describe the correlation between the cyclic load's amplitude and the number of cycles to failure in the SL approach. We refer to this as the conventional fatigue method, which is further discussed in Section 2.3. We could state that DT, on the other hand, focuses on the material's state in the vicinity of the notch tip, which is the topic of fracture mechanics. The topic of fatigue and fracture is broad and the underlying methods are not easily divisible into the scope of SL and DT approaches. Modeling in fatigue can be further divided into ultra-low, low, high and ultra-high cyclic fatigue. At very high cyclic loads the stress state exceeds the yielding point and the deformation is predominately plastic. The material is limited in the accumulation of load cycles and this is termed the ultra-low cyclic fatigue [5][6][7]. Here the distinction between the strain hardening and cyclic damage is not so clear and cyclic damage models can also be derived based on the dissipated plastic energy potential for monotonically loaded material [8]. Stress amplitude is lower in the low-cyclic regime. In the case of a proportional and non-varying cyclic load, some materials exibit a stabilised plastic dissipation curve also denoted as a hysteresis curve in the stress-strain space [9,10]. It is considered that the accumulation of plastic deformation is the predominant effect that leads to fatigue crack nucleation. The strain-life [11,12] and energy-life [13,14] approaches can be used model this regime. By further lowering the cyclic load amplitude we continue into the field of high-cyclic and further into ultra-high-cyclic fatigue [15]. Here the magnitude of elastic strain is dominant over the plastic one and environment, surface condition and other effects should be considered [16]. The goal of all fatigue approaches is to somehow state the number of load cycles before the crack appears. Analysis of fracture is thus the natural continuation. If we limit our scope to the macro level, thus jumping across the field micro-mechanics of cracks, a somewhat similar division is noticeable in fracture as is in the fatigue. The division is made on the basis of crack length e.g., short and long crack and finally a region of complete fracture, when the crack has grown to the critical length [17]. The fundamental mechanics of fracture is different from that of fatigue [18] where stress and strain are replaced by the gradient or speed of crack growth and stress intensity. Despite the fact that we have limited our approach to macro-fatigue observations, the complexity is obvious. In addition, the SL and DT approaches use incompatible fatigue models. This means that in the product-design phase, a fatigue analysis is limited to one of the two approaches, not only due to incompatibility, but also due to the fact that the process of building fatigue models is heavily dependent on time-consuming laboratory tests of materials.
This paper presents a modeling framework that tries to resolve the difficulties described above. The end goal is a discussion of an unifying approach to the total fatigue-life calculation. More explicitly the research is exploring the possibility of using conventional fatigue, here strain-life approach, to also model the process of fracture. This in term can result in the lowering of resources and shortening of time needed in the process of design validation.
The framework is built upon observations of cyclically loaded material, which are presented in the next section. The usability of fatigue calculations is very dependent on the speed of the achieved results [19]. Here, a computational technique known as direct cyclic analysis is used to calculate the state of the material. This subject is further developed in Section 2.2. A crack is de-facto damage and a means to describe it must be presented. Here, we applied the continuum-damage-mechanics approach in describing the effect of damage. This is now a well-evolved topic used to describe various damage processes. We present the basic idea in Section 2.4. The finite-element method (FEM) was used to simulate the crack's nucleation and propagation. The FEM is versatile, but built upon continuum mechanics. Therefore, direct modeling of the dis-continuum crack's nature is impossible and heavily mesh dependent. To circumvent this to some degree, we propose the use of cohesive contacts in the FEM mesh to model the presence of cracks. This topic is discussed in Section 2.5. Section 2.6, synthesizes all of the approaches in an algorithm, which is tested on two load cases, and the results and conclusion cover Sections 3 and 4. Throughout the paper we reference observations of a cyclically loaded, low-carbon S235 steel. Here, it is only used to provide a real database and was arbitrarily chosen.

Material Observations and Framework Idea
In fatigue, the cyclic material properties are extracted from a uniaxial, cyclically loaded specimen. A round, dog-bone specimen with the geometry shown in Figure 1 was manufactured from a hot-rolled S235 steel bar. The specimens were strain-controlled loaded with a ratio R = min max = −1. The testing was conducted on an MTS Landmark 100-kN servo-hydraulic tensile-testing machine, equipped with a MTS 634.11F-24 axial extensometer, both from MTS systems corporation, Eden Prairie, MN, U.S.A. The relevant results are collected in Table 1.   Figure 2 shows the evolution of the peak force for one of the specimens. The phenomenon is complex. Initially, the material either hardens or softens, which can be seen as a change in the peak force. The rate of change gradually reduces to a uniform response, which we describe as a quasi-stable cyclic response. Near the end of the experiment we observe another change in the peak force, as a decline up to the point of failure. To model this observation, certain simplifications have to be made. In the case of the S235 steel, the initial hardening/softening behavior appears to be strain-amplitude and strain-rate dependent. As this transition is relatively short compared to the total number of cycles to failure, we decided to dismiss it from the modeling. The quasi-stable response, on the other hand, accounts for the majority of the specimen's life. A combined isotropic-kinematic, hardening material model with three back-stress evolution formulations was used here [20]. The advantages and disadvantages of this model and the parameter estimation are described in [21]. Here, Table 2 lists only the parameters used for the model given in Equation (1), where C and γ are the material tangential hardening modulus and plastic strain p dependent memory coefficient.
Peak-force evolution for S235 steel at a = 0.5%. Our unified framework uses the FEM to calculate the material's response. The direct simulation of the cyclically loaded domain increases the demand for computational power significantly, and is therefore not applicable for real-world scenarios. Still, a means to evaluate this stable state is needed. Our framework is based on the direct cyclic algorithm to overcome this obstacle. An elaboration of this point is provided in Sections 2.2 and 2.3. Now we return to the last observed phase in Figure 2. The peak load's decline is associated with an accumulation of damage to the specimen. Here, microcracks that were evolving in the stabilized regime coalesce into macrocracks, which in turn reduces the load-bearing surface. With strain-controlled testing this is seen as a reduction in the nominal stress and the restive force. To model this phase, the continuum-damage-mechanics approach was used, see Section 2.4. Damage is seen as the crack's nucleation and growth. To incorporate the crack model into our framework, the cohesive-contact approach was used. It is derived from cohesive-elements modeling that is used in fracture simulations. Further discussion is given in Section 2.5.

Direct Cyclic Analysis
Despite the fact that fatigue loading is a dynamic process, some materials exhibit what is referred to as a stabilized cyclic response. In the FEM both implicit and explicit schemes can be used to simulate this cyclic response, including the initial hardening or softening. Due to the extensive computational effort, such approaches are only meaningful for a limited number of cycles. There is another numerical approach called direct cyclic analysis (DCA) [22]. This is a quasi-static analysis that optimizes the displacement function of the structure for each time increment of the load cycle. For this purpose, the system displacement function u(t) (Equation (2)) is described as a Fourier series .
For each iteration time, the non-linear equilibrium equation (Equation (3)) is solved.
where R, F and I describe the residual, cyclic and internal force vectors. If we denote c k as the correction coefficients of the displacement function and K as the system stiffness matrix, then the equilibrium equation can be written as Equation (4).
To be coherent with the displacement function, the residual force vector R is also written as a Fourier series on an element basis (Equation (5)) The iterative scheme progresses at the end of each loading cycle by solving for the displacement corrections and adding them to the next displacement coefficients (Equation (9)) To validate the proposed method, experimental observations were compared with the numerical direct cyclic results and are shown in Figure 3. For the DCA calculation, the combined isotropic-kinematic hardening material model was used. This is described in the previous section.

Conventional Fatigue Approach
The conventional fatigue approach refers to the laws derived from correlations between the cyclic load and the number of completed load cycles to failure. This branch of modeling dates back to the 19th century, when A. Wöhler [23] proposed an exponential decay law between the cyclic stress amplitude and the number of cycles to failure. Because this type of modeling is still used today, numerous derivations are available. Here, we mention the Coffin-Manson model [24,25] that correlates the cyclic strain amplitude with the number of cycles to the point when a noticeable crack appears (Equation (12)).
Here, a stands for the cyclic strain amplitude, E is the material's Young's modulus and S f an E f are the fatigue strength and strain coefficients. The model is a linear addition of the elastic and plastic strains that is decomposed from a . This is also the basis for the division into low-and high-cycle fatigue. In low-cycle fatigue, the plastic strain amplitude is dominant. The transition from low-to high-cycle fatigue occurs at approximately N f ≈ 10 5 cycles for steel, or when the elastic strain amplitude becomes dominant. To complete the model, two exponential material parameters c and b must be fitted to the measured material's response. The model parameters for S235 steel are given in Table 3 and the results are shown in Figure 4.  If a material is subjected to a variable-amplitude load, a model is needed to describe the effect of each load separately. Research shows [26,27] that the order of the cyclic load's amplitude has a major effect on the fatigue life, not just the number of occurrences. Because the accumulation laws of the N f calculation are not the main topic of this paper, the simplest version was used, i.e., the Palmgren-Miner linear-accumulation law (Equation (13)).
Here, n i is the number of cycles at a given load and N i is the number of cycles to failure N f , at the same load. Here, D f describes the damage accumulation in the interval [0, 1], where at D f = 1 a noticeable crack appears.

Continuum Damage Mechanics
Continuum damage mechanics (CDM) describes the mechanical properties of materials coupled with the damage. In essence, damage is derived on the basis of the thermodynamics of irreversible processes. One such process could be derived from the accumulating plastic strain under fatigue loading, which eventually leads to the onset of damage in the form of microcracks. With continued loading, the damage accumulates as the cracks grow to the meso and then to the macroscale, until finally complete rupture occurs. CDM incorporates the damage parameter in a fully coupled scheme with the material's stressstrain response. An example in Figure 5 shows the principle of the damage effect on the load-bearing surface [28,29]. If we limit the scope to isotropic damage, the direction-invariant damage is described by the scalar variable as Equation (14).
In tension, the increase in the damage and thus the decrease in the load-bearing surface can be understood as a change in the material's stiffness. In a uniaxial case, a simple constitutive model (Equation (17)) can be derived for the elastic strain part from the material's effective Young's modulus E.

Cohesive Contact
In the FEM both the cohesive-elements [30][31][32], and the cohesive-contact approaches define the interference layer of the domains with some sort of traction-separation law. The most significant difference is the thickness of the two approaches. The cohesive elements have a finite thickness, whereas the cohesive contact imposes none. Cohesiveelement and cohesive-zone modeling have been used extensively to model crack growth, both under monotonic loads [33][34][35] and fatigue loads [36][37][38]. In the regions with expected crack nucleation, the FEM mesh is separated and an interstitial cohesive element rebinds the mesh as a thin layer of added material. The constituted material model of this added layer is normally given as a traction-separation law, following the principles of linear-elastic fracture mechanics and virtual crack-opening displacement.
The development of cohesive-zone modeling for fatigue-crack growth appears to trail the formulations of the LEFM, as the complexity increases significantly in the transition from static brittle-crack to fatigue-crack modeling, resulting in numerous model approaches. Different FEM modeling schemes exist in fatigue-crack simulation, for example, element deletion or softening [39], node separation [40,41], XFEM [42] and intermittent cohesive elements. We must focus our attention on the principles that capture the nature of the phenomenon as closely as possible. Element deletion introduces an unreal volume and mass decrease, which is negligible in fatigue cracking. Node separation is abrupt as the simulated crack grows instantly on the element edges, which is de facto a serious meshdependency concern. The now widely accepted XFEM [43] approach minimizes the mesh dependency in the crack-growth simulation. Both fatigue-crack nucleation and propagation can be modeled inside this framework. Some results [44] indicate that the XFEM procedure leads to a numerical instability in arbitrary crack-nucleation-site prediction.
Lastly, the cohesive-element method of crack simulation introduces new fictional material between the edges of the failing FEM elements. As this is hard to model at the point of crack nucleation, normally, a beforehand global element separation and reconnection with cohesive elements is introduced [45][46][47]. This heterogeneous FEM mesh is a poor model of the actual structure as the cohesive elements introduce additional softening. Here, we propose to introduce the cohesive contact (CC) [22] instead of a cohesive element at the boundary of FEM elements only at the instance of the predicted failure (see Figure 6). The CC is modeled as a FEM contact when the interface is loaded in compression and as a linear elastic traction material when loaded in tension. The material model of the CC interface is an elastic constitutive matrix that relates the normal and shear interface loads with an appropriate displacement. In the local 3D coordinate system index 1 denotes the normal and indices 2, 3 the shear directions. The traction stress vector of the CC can be written as Equation (19).
Here, t i denotes the interface traction stress vector, K i,j is the stiffness matrix and δ i is the corresponding separation or displacement vector. If K ij is a full matrix it describes a fully coupled model. In this research, a simplified, uncoupled 2D model is used and the system reduces to Equation (20).
As the cohesive-contact models crack, a damage parameter must be introduced. The idea of isotropic damage D as noted in Equation (14) is simply introduced in the elastic stiffness matrix as Equation (21) To model the damage-evolution law we return to the phase-three material observations of a cyclically loaded specimen. The damage D accumulation is seen at the end of the specimen's life as a gradual reduction in the peak cycle force, see Figure 2. Testing was stopped on the basis of a large gradient in the peak-force drop. As such, the load does not fall to 0 N. The initial trend was approximated with a linear function to obtain the number of damage-accumulating cycles for each given nominal strain amplitude. One for each strain-amplitude load is shown in Figure 7. Here, the derived strain amplitude vs the number of cycles to failure is the dataset for the exponential-decay-damage model (see Equation (22) and Figure 8), which is an application of the same principle as the model described in Equation (12). As the observed law represents the effect of crack growth over a specimen area that is considerably larger that that of the CC, a scaling function like Equation (23) is used in order to calculate the number of cycles to failure for each CC.
where N D,CC represents the number of cycles to failure in the CC, A CC and A spec are the CC and specimen cross-section, respectively. The CC stiffness K cannot be arbitrarily set as it must capture the nature of softening between the elements as the crack progresses. The limits for the tensile CC stiffness are set based on [48] and are given as Equation (24).
Here, ν is the material's Poisson's ratio, E is the Young's modulus and L mesh is the characteristic mesh length.

Algorithm
The algorithm of the total fatigue-life calculation presented in the following section is based on the FEM solver Abaqus. The results are then iteratively post-processed using an external routine written in Python. Before diving into the algorithm, we must distinguish between two forms of damage. One is the relative damage up to the crack's nucleation, denoted D f , and the other is the crack damage D that influences the material's properties in CCs. Figure 9 presents the algorithm flow chart. The first operation calculates the stable cyclic response with a direct cyclic analysis. The selected variable is exported from the element-integration points. The field variable to export is damage-model specific, e.g., stress, strain. In the example presented later, the number of cycles to failure is calculated with the Coffin-Manson strain-amplitude approach. What will be presented below is a node-splitting routine for crack simulation. Therefore, the number of cycles to failure is extrapolated into nodes. Here, a distance-weighted average extrapolation was used to describe the result in nodes using its surrounding elements. The principle is shown in Figure 10 and the calculation in Equation (26). Here, N f ,j (e c ) stands for the cycles to the crack's initialization in the element, l j is the distance from element's centroid e c to the node n. The index j represents the node n touching element.   Supposing that the cyclic load function is deterministic and proportional to the algorithm jumps in the calculation by the minimum number of cycles calculated in nodes. The jump is thus made up to the first crack nucleation or D f = 1. The crack's direction and length are governed by the failing node and its neighbor node that has the smallest number of cycles or highest D f . Mesh dependency is thus warranted. The line that connects both nodes, the failed and its neighbor, is on the edge of two elements. The failed node is split and a cohesive contact (CC) is placed between the new surface boundary. The next iteration calculates the stable cyclic response now with an additional CC. First, a check is made in calculating the maximum permissible cycle number jump, which is governed by both nodes and the CC calculation. The strain amplitude in CCs is calculated as an average of the connected elements. The CC set damage jump ∆D and the CC calculated strain amplitude govern the permissible jump in ∆N D number of cycles by the relationship in Equation (22). The permissible node cycle number is calculated by a damage-accumulating function, for example, the simple Palmgreen-Miner linear damage accumulation. If the permissible number of cycles in CC is smaller than that in nodes, a CC damage routine executes, which also appends to the D f damage in nodes. If, on the other hand, the node cycle calculation is smaller, a new CC is introduced at the new fail node.

Application of the Proposed Framework on a Compact Tension Specimen
The geometry of the specimen, which resembles a compact tension specimen, is presented in Figure 11. Due to the current limitations of the proposed algorithm, the domain is reduced to a 2D domain with the implied assumption of a plain-strain state. The bottom pin hole is a rotational support and the load is applied on the top hole. The domain was meshed with linear 3-node plain strain CPE elements. Size of the mesh was 0.5 mm near the specimen notch with a gradual growth up to 5 mm elsewhere. Total number of elements is 4449. For the purposes of the presentation a damage step ∆D = 1 was used.  Figure 12 shows two results. In both cases an alternating load is applied, with magnitudes of 2 and 10 kN. In the initial phase of crack growth we observe the jaggedness as the crack progresses along the element edges. The 10-kN case also highlights the limitation of the proposed algorithm. As the crack grows, the load-bearing surface decreases. This results in a wide area of high stress in front of the crack tip, which leads to the concurrent failure of multiple nodes. This is seen as a deviation and zig-zagging of the crack's front.
The results in Figure 12 confirm that the proposed algorithm correctly predicts the crack's nucleation site and its propagation path. This is all that can be claimed in the current state of the algorithm, as no direct comparison is made with a laboratory experiment.

Limitations and Observations of the Proposed Framework
Additionally, we would like to address the limitations of the algorithm in its current state. Mesh dependency is one of the major concerns, as it influences the path of the crack's growth and also the calculation in the area of the crack's tip. The crack can grow only along the element edges. Therefore, triangular elements were used to increase the number of available paths (compared to a quadrilateral mesh). The node-splitting method results in a sharp geometry that represents the crack's front. This leads to a singularity point in the calculation of the FEM variables that cannot be remedied by a further reduction in the size of the mesh. By defining the crack's tip area with an elastic cohesive contact, this effect is minimized, but a further evaluation is necessary to evaluate what has been observed.
The idea of implementing a step-wise approach is based on the goal of reducing the computational time needed compared to the direct approach. As stated before, this is based on the assumption that when subjected to fatigue the material stabilizes and a proportional and in-phase load is applied. This results in the calculation of an average response up to the point of the predicted crack nucleation. As the crack nucleates with the introduction of the first cohesive contact, the step in the calculation must reduce to describe the characteristically transient phenomenon of crack growth. A direct approach would be more appropriate, but changes to the FEM algorithms mid analysis can be troublesome. Here, the FEM algorithm remains as DCA and a ∆D in the CC governs the step. The rate of the crack's growth is thus directly specified according to the current state in the CC and the size of ∆D, as it remains constant until the next step in the calculation. The ∆D parameter is thus balanced between the immediate growth of the crack on the length of the element with ∆D = 1 and the direct computation method with ∆D → 0.
As seen on Figure 12 the framework fails to predict the direction of fatigue crack growth if the load state in front of the crack tip is to severe. This causes rapid accumulation of damage in the surrounding elements that is extrapolated to the nodes. As the current CC fails and the framework proceeds in the selection of the next CC the choice of the crack tip neighbour node is arbitrary. With following iteration, this is seen as a zig-zaging of the crack, which is believed to be unrealistic and thus results at this point should be dismissed.

Conclusions
The proposed algorithm can be successfully used to simulate both fatigue and growth of a crack. The DCA method greatly reduces the time of the calculation up to the point of the crack's nucleation. The crack's nucleation site depends on the loading and the geometry conditions without the need for a direct specification. In the proposed scope, only the results of strain-controlled fatigue experiments were used, thus minimizing the need for specific crack-growth-rate experiment. The introduction of cohesive contacts is an elegant way to simulate damage-coupled crack propagation. This allows us to include the crack-closure effect and simulate possible crack branching. The proposed algorithm was numerically tested on a specimen shape similar to a compact test specimen. The results indicate that the crack's nucleation site and the propagation path are correct, but additional experiments are needed to validate the proposed approach. Highlights:

•
Conventional fatigue approach, here the strain-life method, can be used to simulate both fatigue and crack growth. • The direct cyclic analysis reduces calculation time in fatigue if the cyclic load is considered proportional. • The locally defined cohesive contact can simulate crack if coupled with a damaging parameter. • Results indicate that the crack grows in a predictable direction, but further experimental observations should be made to validate its speed. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: