The Evaluation of Front Shapes of Through-the-Thickness Fatigue Cracks

: Fatigue failure of structural components due to cyclic loading is a major concern for engineers. Although metal fatigue is a relatively old subject, current methods for the evaluation of fatigue crack growth and fatigue lifetime have several limitations. In general, these methods largely disregard the actual shape of the crack front by introducing various simpliﬁcations, namely shape constraints. Therefore, more research is required to develop new approaches to correctly understand the underlying mechanisms associated with the fatigue crack growth. This paper presents new tools to evaluate the crack front shape of through-the-thickness cracks propagating in plates under quasi-steady-state conditions. A numerical approach incorporating simpliﬁed phenomenological models of plasticity-induced crack closure was developed and validated against experimental results. The predicted crack front shapes and crack closure values were, in general, in agreement with those found in the experimental observations.


Introduction
The evaluation of fatigue life and failure conditions of structural components is of permanent and primary interest for engineers. Over the past five decades, significant progress has been made toward the development of more appropriate fatigue crack growth models and life assessment procedures. Significant research effort has been directed to the study of the fatigue crack closure phenomenon, which was first introduced by Elber [1] to explain the experimentally observed features of fatigue crack growth in aluminium alloys. The number of publications grew rapidly since this pioneering study, and continues to grow. It is now commonly accepted that the contributions of various crack closure mechanisms, specifically plasticity-induced crack closure, roughness-induced crack closure, and oxideinduced closure, are significant, and these mechanisms are capable of explaining many fatigue crack growth phenomena, e.g., the influence of thickness on crack growth rates, retardation effects associated with overloads, or higher propagation rates of small cracks in comparison with long cracks [2].
It is well-established that for relatively long cracks propagating in a non-aggressive environment, the plasticity-induced crack closure dominates over the roughness-induced crack and oxide-induced closures. The plasticity-induced crack closure models rely on far fewer assumptions than the two other closure mechanisms. The first theoretical model was developed by Budianski and Hutchinson [3] based on the two-dimensional Dugdale stripyield model [4]. The theoretical results demonstrated that opening stress intensity factor is surprisingly high, and increases with an increase in the R ratio. All early crack closure models for plate components utilised both plane strain and plane stress simplifications, although real cracks are inherently three-dimensional (3D). To examine the thickness effect on crack propagation rates, empirical constraint factors were often used, demonstrating a stronger correlation with experimental results. With the advance of numerical methods and the increase in computational power, it became possible to eliminate these simplifications and study more realistic geometries, as well as various 3D effects [5,6].
In 3D problems, the order of the singularity at the intersection of the crack front with the free surface depends on Poisson's ratio and the intersection angle. From energy considerations, it follows that fatigue cracks have to preserve the 1/ √ r singularity. Therefore, the fatigue crack has to intersect the free surface at a critical angle, β cr , which is a function of Poisson's ratio. Several experimental studies have reported that, at least, the Mode I fatigue crack front is shaped to ensure the square root singular behaviour along the entire crack front. However, it seems that the effect of 3D corner singularity is not very significant in the presence of a sufficiently large crack front process zone [7]. This is because the 3D corner singularity effect is a point effect, and is very much localised. Therefore, it might be negated by the plasticity and damage formation near the surface. For example, in an experimental study of steel circular bars subjected to bending and torsion, the experimental intersection angles were found to be very different from the theoretically predicted critical angles [8].
Considerably less effort has been directed toward the study of the effects of the 3D corner singularity and elasto-plastic constraints on plasticity-induced crack closure. Generally, the direct 3D elasto-plastic simulations of fatigue crack growth demand much greater computational resources [9]. These simulations have many issues associated with the validation of the numerical solution and the accuracy of the obtained results. A number of factors affect the accuracy, which are difficult to control: the mesh refinement, the type of finite element, the crack advance scheme (which usually consists of releasing nodes ahead of the crack front), contact conditions, and the local criterion of crack front opening. Branco et al. [10] recently provided an exhaustive review concerning these aspects. The overall conclusion was that the direct numerical approaches are capable of describing the shape evaluation of fatigue cracks. However, the application of these approaches to particular problems can be quite cumbersome. Each problem needs a large effort to calibrate the solution and verify the results. These efforts are usually focused on the reduction in the number of finite elements, the number of simulations required in the analysis, or, eventually, the computation time, which cannot be considered to be of practical relevance [11].
In the present paper, a simplified procedure for the evaluation of the fatigue crack front shapes of through-the-thickness cracks propagating under the cyclic loading conditions is presented. The procedure is based on simplified methods for the evaluation of the plasticity-induced crack closure effect, namely the equivalent-thickness method introduced by Yu and Guo [12,13], as well as the analytical model developed by Kotousov et al. [14,15]. The outcomes of the simulation are compared with available experimental results obtained at the same propagation conditions for validation purposes. The paper is organised as follows: Section 2 addresses the method used to evaluate the crack front shape, as well as the models introduced to evaluate the crack closure along the crack front. Section 3 describes the finite element model developed to calculate the stress intensity factors along the crack front. Section 4 compares the predicted crack front shapes with those obtained experimentally for different materials and propagation conditions. The paper ends with some concluding remarks.

Crack Shape Simulation and Crack Closure Models
The main idea behind the evaluation of the steady-state shape of a fatigue crack front proposed in this paper is to select a curve from a parametric family that minimises the deviation of the fatigue driving force along the crack front. In other words, we first specified a possible parametric set of curves in the crack plane (e.g., parabolic, hyperbolic, or elliptical shapes) and then evaluated the local fatigue driving force using a finite element model, along with simplified plasticity-induced crack closure models. In this study, the local fatigue driving force was defined by the effective stress intensity factor range, ∆K eff , given by the formula: where U is the normalised load ratio parameter, or normalised effective stress intensity factor, which is often used to describe the effects of loading and geometry on crack closure, and ∆K is the traditional linear-elastic stress intensity factor range [16] defined by the maximum and minimum values of the stress intensity factor experienced for a given load cycle. The load ratio is therefore given by R = K min /K max . In the case of 3D problems, this normalised load ratio is not a constant, but rather a function of the position along the crack front, U = U(z). Thus, the local crack growth rate is a function of the effective stress intensity factor range, i.e., where K op (z) is the local opening load stress intensity factor, which corresponds to the minimum load at which the crack faces, at point z, which are fully separated. A number of sophisticated finite-element (FE) models were developed to evaluate U(z) for different geometries and loading conditions. However, as discussed above, these models have many limitations, and are quite difficult to apply in fatigue calculations. Below, we consider two simplified methods for the evaluation of the normalised load ratio, the equivalent-thickness model introduced by Yu and Guo [13], and the analytical model proposed by Kotousov et al. [14,15], which are addressed in Sections 2.1 and 2.2, respectively. These methods will be further incorporated into the 3D linear elastic finite element simulations to evaluate the shape of the through-the-thickness cracks. This evaluation will be performed via the corner singularity method [17], which is briefly presented in Section 2.3.

Equivalent-Thickness Model
For through-the-thickness cracked plates, She et al. [17] proposed defining the equivalent thickness based on a numerical analysis of the 3D distribution of the out-of-plane stresses and constraint factor, T z , which is defined as: where σ x , σ y , and σ z are the normal stresses. This method is illustrated in Figure 1. The equivalent thickness, 2h eq , for point P on the crack front is identified as the plate thickness, which leads to the same distribution of T z at the mid-plane.
deviation of the fatigue driving force along the crack front. In other words, we first specified a possible parametric set of curves in the crack plane (e.g., parabolic, hyperbolic, or elliptical shapes) and then evaluated the local fatigue driving force using a finite element model, along with simplified plasticity-induced crack closure models. In this study, the local fatigue driving force was defined by the effective stress intensity factor range, ΔKeff, given by the formula: where U is the normalised load ratio parameter, or normalised effective stress intensity factor, which is often used to describe the effects of loading and geometry on crack closure, and ΔK is the traditional linear-elastic stress intensity factor range [16] defined by the maximum and minimum values of the stress intensity factor experienced for a given load cycle. The load ratio is therefore given by R = Kmin/Kmax. In the case of 3D problems, this normalised load ratio is not a constant, but rather a function of the position along the crack front, U = U(z). Thus, the local crack growth rate is a function of the effective stress intensity factor range, i.e., where Kop(z) is the local opening load stress intensity factor, which corresponds to the minimum load at which the crack faces, at point z, which are fully separated. A number of sophisticated finite-element (FE) models were developed to evaluate U(z) for different geometries and loading conditions. However, as discussed above, these models have many limitations, and are quite difficult to apply in fatigue calculations. Below, we consider two simplified methods for the evaluation of the normalised load ratio, the equivalent-thickness model introduced by Yu and Guo [13], and the analytical model proposed by Kotousov et al. [14,15], which are addressed in Sections 2.1 and 2.2, respectively. These methods will be further incorporated into the 3D linear elastic finite element simulations to evaluate the shape of the through-the-thickness cracks. This evaluation will be performed via the corner singularity method [17], which is briefly presented in Section 2.3.

Equivalent-Thickness Model
For through-the-thickness cracked plates, She et al. [17] proposed defining the equivalent thickness based on a numerical analysis of the 3D distribution of the out-of-plane stresses and constraint factor, Tz, which is defined as: where σx, σy, and σz are the normal stresses. This method is illustrated in Figure 1. The equivalent thickness, 2heq, for point P on the crack front is identified as the plate thickness, which leads to the same distribution of Tz at the mid-plane. An empirical equation was suggested to evaluate the equivalent-thickness as follows: where z is the distance from the mid-plane and h is the half-thickness of the plate. The normalised load ratio parameter in this method can be calculated as follows: where κ is a function of the R ratio: (6) and α g is a global constraint factor, α g , defined by the formula: where ν is the Poisson's ratio and t is given by: with: where σ 0 is the flow stress. These empirical equations were extended to the corner, and surface cracks and were extensively validated using 3D finite element analyses.

Analytical Model for the Evaluation of Crack Closure
Another method for the evaluation of local plasticity-induced closure is based on a simplified 3D analytical model. In accordance with this model, the parameter U for Mode I loading under small-scale yielding conditions can be approximated from the following expression: where the fitting functions a, b and c can be written in the form: The above equations were obtained within the first-order plate theory based on the Budiansky-Hutchinson crack closure model [3,15]. The results, which correspond to the classical two-dimensional theories (or plane stress state, or plane strain state), can be obtained as limiting cases of very thin and very thick plates, i.e., when η → ∞ or η → 0, respectively. The details of the derivation of these equations can be found in the original paper by Codrington and Kotousov [14].

Corner Singularity Method
In this study, the evaluation of the steady-state front in the through-the-thickness cracks was carried out using the corner singularity method. First, we approximated the shape of the crack front by a two-parameter elliptical curve, which can be described as: where a and b are the major and minor semi-axes of an ellipse, respectively, as shown in Figure 2. 021, 11, x FOR PEER REVIEW 5 of 14 classical two-dimensional theories (or plane stress state, or plane strain state), can be obtained as limiting cases of very thin and very thick plates, i.e., when η → ∞ or η → 0, respectively. The details of the derivation of these equations can be found in the original paper by Codrington and Kotousov [14].

Corner Singularity Method
In this study, the evaluation of the steady-state front in the through-the-thickness cracks was carried out using the corner singularity method. First, we approximated the shape of the crack front by a two-parameter elliptical curve, which can be described as: where a and b are the major and minor semi-axes of an ellipse, respectively, as shown in Figure 2. The crack front tends to intersect the free plate surface at the critical angle, βc, when the plasticity effects are small. The critical angle is a function of the Poisson's ratio and the type of loading. We found that the critical intersection angle can be approximated by the following formula [18]: where ν is the Poisson's ratio. Typically, when the size of the plastic zone is greater than 1% of the plate thickness, the stress state near the vertex location is not controlled by the elastic singularity. In these cases, the plasticity effects become more important, and together with the vertex singularity effect, lead to greater critical angles for elastic-plastic materials. To find , we need to make sure that: where b is defined by: Substituting Equation (15) into Equation (12), we obtain: The crack front tends to intersect the free plate surface at the critical angle, β c , when the plasticity effects are small. The critical angle is a function of the Poisson's ratio and the type of loading. We found that the critical intersection angle can be approximated by the following formula [18]: where ν is the Poisson's ratio. Typically, when the size of the plastic zone is greater than 1% of the plate thickness, the stress state near the vertex location is not controlled by the elastic singularity. In these cases, the plasticity effects become more important, and together with the vertex singularity effect, lead to greater critical angles for elastic-plastic materials.
To find b, we need to make sure that: Substituting Equation (15) into Equation (12), we obtain: This equation meets the condition that the crack front intersects with the free surface at the critical angle given by Equation (13), and represents a parametric curve with one single parameter, a. Further, the steady-state condition of the crack propagation requires that the projection of the effective stress intensity factor along the crack propagation direction (x-direction, Figure 1) is constant for all points along the crack front. This condition cannot be satisfied exactly with any multi-parametric equation describing the possible crack front shapes. However, the shape that minimises the difference of the effective stress intensity factor along the crack front can be considered as the best approximation of the actual fatigue crack front shape.

Numerical Approach
This section describes the numerical model developed in this research to determine the stress intensity factor ranges along the crack front. The stress intensity factor ranges, along with the crack closure models described in the previous section, enabled the computation of the local fatigue driving force, which was used to obtain a steady-state crack front shape. The steady-state crack front shape was selected as the one producing the minimum deviation of ∆K eff along the crack front. This evaluation needs to be completed for each curve from the parametric set.
To reduce the computational overhead, we developed a simplified geometry by introducing adequate boundary conditions, capable of describing 3D effects near the crack front. Section 3.1 describes the details of the numerical modelling, and Section 3.2 addresses the boundary conditions considered in this paper. The last section, Section 3.3, is devoted to the validation of the stress intensity factor values obtained with the proposed approach.

Finite Element Model Description
The typical finite element geometry, developed here to study a through-the-thickness crack in an elastic plate, is shown in Figure 3. As can be seen, the rectangular cross-section geometries were modelled to evaluate the stress and displacement fields near the crack tip. The size of the finite element models is sufficient to avoid the effect of the finite boundaries on the stress state. By taking advantage of the symmetry conditions (i.e., XY symmetry, XZ symmetry, and YZ symmetry), only one-eighth of the crack problem was modelled. The height of the FE models taken was approximately ten times larger than the plate thicknesses. In accordance with the previous studies, this is sufficient to accurately describe the 3D effects near the crack front [19,20]. The FE models corresponding to different values of ( Figure 1) were meshed with linear 8-node hexahedral elements of type C3D8R. A reasonably uniform element grid with a structured mesh was considered. A denser mesh, with a spider-web pattern ( Figure  3c) was used near the crack front, where the stress gradients were expected to be maximum (Figure 3b), consisting of 5 concentric rings centred at the crack tip with a radial discretisation of 10° (Figure 3c). Thirty nodes along the plate half-thickness (Figure 3b) were used to define the crack front shape. The specimen was subjected to uniaxial loading The FE models corresponding to different values of a ( Figure 1) were meshed with linear 8-node hexahedral elements of type C3D8R. A reasonably uniform element grid with a structured mesh was considered. A denser mesh, with a spider-web pattern (Figure 3c) was used near the crack front, where the stress gradients were expected to be maximum (Figure 3b), consisting of 5 concentric rings centred at the crack tip with a radial discretisation of 10 • (Figure 3c). Thirty nodes along the plate half-thickness (Figure 3b) were used to define the crack front shape. The specimen was subjected to uniaxial loading applied at the bottom surface (i.e., at the XZ-plane with a Y-coordinate equal to H/2). The assembled mode is exhibited in Figure 3a. Further details about the modelling approach can be found in papers published by the present authors [19,21].
The numerical simulations were carried out using Abaqus/CAE 2020 (© Dassault Systèmes, 2019), assuming a homogeneous, isotropic, and linear-elastic behaviour. The mechanical properties inserted into Abaqus/CAE 2020 to perform the numerical simulations were the Young's modulus and the Poisson's ratio of the tested materials ( Table 1). The displacement field far from the crack tip was calculated in accordance with the William's solution using MATLAB R2020b, and the obtained results were applied for the boundary conditions. The 3D solutions of the J-integral were used to calculate the stress intensity factor near the crack front. One layer of elements surrounding the crack front was used to calculate the first contour integral. The additional layer of elements was used to compute the subsequent contours. The different contour solutions were approximately coincident after eight contours. The results from averaging contours five through eight was considered. A similar strategy, either in terms of mesh framework or simulation analysis, was carried out for all geometries and crack configurations studied in the present paper.

Boundary Conditions
The plane-stress displacements far from the crack tip were calculated in accordance with William's solution [22]: u y (r, θ) = r 2π Being: f I y (θ) = sin where r is the distance from the crack tip, θ is the angle measured from the symmetry line, K ∞ I is the remotely applied Mode I stress intensity factor, and k is Kolosov's constant for plane stress and plane strain conditions. The plane stress k value was considered in the boundary conditions, i.e., where ν is the Poisson's ratio. Bakker [23] showed that a cracked plate under plane stress undergoes a change to plane strain behaviour near the crack tip. He proved that the radial position, where the plane stress to plane strain transition takes place, strongly depends on the position in the thickness direction. The degree of plane strain is essentially zero at distances from the tip greater than five times the thickness, even in the middle plane of the plate [24].

Validation Study
The numerical results obtained for the maximum stress intensity factor are presented in Figure 4 as a function of the thickness for a Poisson's ratio of 0.3. The classical results for both the plane stress state and the plane strain state are also given in Figure 4. It is evident from Figure 4 that the stress intensity factor changes with the thickness of the plate until the thickness exceeds a critical value. In this particular problem, the results showed that the critical thickness is 25 mm. Once the thickness exceeds the critical dimension, the stress field in the vertex singularity region has a negligible impact on the behavior of the whole structure. The stress intensity factor becomes relatively constant in the sufficiently thick plate, and is equal to the value for plane strain conditions.

Crack Front Shape Evaluation and Comparison with Experimental Studies
The proposed method for the evaluation of the steady-state crack front shapes was compared against two independent experimental studies. The specimen geometries used in the experimental tests are exhibited in Figure 5, and were made of 6082-T6 aluminium alloy and polymethyl methacrylate (PMMA), separately. The main mechanical properties of both materials are listed in Table 1. The former (Figure 5a) consisted of a standard middle-crack tension specimen with a thickness of 3 mm [11,25]. The tests were conducted under constant-amplitude axial loading using a stress ratio equal to 0.25. Figure 6a shows an example of the typical fracture surfaces obtained in the tests. Fatigue cracks grew over a sufficiently large distance from the initial notch to ensure the quasi-steady-state conditions of propagation. The beach-marking technique was applied to mark the crack front at the fracture surface.

Crack Front Shape Evaluation and Comparison with Experimental Studies
The proposed method for the evaluation of the steady-state crack front shapes was compared against two independent experimental studies. The specimen geometries used in the experimental tests are exhibited in Figure 5, and were made of 6082-T6 aluminium alloy and polymethyl methacrylate (PMMA), separately. The main mechanical properties of both materials are listed in Table 1. The former (Figure 5a) consisted of a standard middle-crack tension specimen with a thickness of 3 mm [11,25]. The tests were conducted under constant-amplitude axial loading using a stress ratio equal to 0.25. Figure 6a shows an example of the typical fracture surfaces obtained in the tests. Fatigue cracks grew over a sufficiently large distance from the initial notch to ensure the quasi-steady-state conditions of propagation. The beach-marking technique was applied to mark the crack front at the fracture surface.
of both materials are listed in Table 1. The former (Figure 5a) consisted of a standard middle-crack tension specimen with a thickness of 3 mm [11,25]. The tests were conducted under constant-amplitude axial loading using a stress ratio equal to 0.25. Figure 6a shows an example of the typical fracture surfaces obtained in the tests. Fatigue cracks grew over a sufficiently large distance from the initial notch to ensure the quasi-steady-state conditions of propagation. The beach-marking technique was applied to mark the crack front at the fracture surface.   Regarding the latter (Figure 5b), the specimen geometry was made of polymethyl methacrylate. It had a rectangular cross-section (Figure 5b), with a thickness of 40 mm [26,27], and an initial straight notch at the middle of the specimen. The tests were conducted under four-point bending loading conditions using a stress ratio equal to 0. The crack front shape was evaluated in situ using a high-resolution digital camera. As in the previous case, fatigue cracks propagated over a sufficiently large distance from the initial notch to ensure the quasi-steady state conditions of propagation. An example of the crack front shapes observed in the experiments is exhibited in Figure 6b.  Figure 7a,b displays a comparison of the experimental crack front shapes and those obtained with the proposed methods for the 6082-T6 aluminium alloy and PMMA, respectively. Overall, the results showed that the equivalent-thickness method provides a satisfactory approximation for the fatigue crack propagation under small yielding conditions. Moreover, the experimental results confirmed that the angle at which the crack front intersects the free surface is greater than the proposed empirical equations in the sufficiently plastic materials. We think that the careful combination of the hyperbolic and elliptical functions might provide accurate crack front shape estimation in the presence of residual stresses or large crack closure effects. The good agreement demonstrated in the previous analysis confirmed the possibility of the accurate evaluation of stress intensity factors using the proposed approach in materials controlled by 3D corner singularity effects. Regarding the latter (Figure 5b), the specimen geometry was made of polymethyl methacrylate. It had a rectangular cross-section (Figure 5b), with a thickness of 40 mm [26,27], and an initial straight notch at the middle of the specimen. The tests were conducted under four-point bending loading conditions using a stress ratio equal to 0. The crack front shape was evaluated in situ using a high-resolution digital camera. As in the previous case, fatigue cracks propagated over a sufficiently large distance from the initial notch to ensure the quasi-steady state conditions of propagation. An example of the crack front shapes observed in the experiments is exhibited in Figure 6b. Figure 7a,b displays a comparison of the experimental crack front shapes and those obtained with the proposed methods for the 6082-T6 aluminium alloy and PMMA, respectively. Overall, the results showed that the equivalent-thickness method provides a satisfactory approximation for the fatigue crack propagation under small yielding condi-tions. Moreover, the experimental results confirmed that the angle at which the crack front intersects the free surface is greater than the proposed empirical equations in the sufficiently plastic materials. We think that the careful combination of the hyperbolic and elliptical functions might provide accurate crack front shape estimation in the presence of residual stresses or large crack closure effects. The good agreement demonstrated in the previous analysis confirmed the possibility of the accurate evaluation of stress intensity factors using the proposed approach in materials controlled by 3D corner singularity effects. This methodology can also be applied to conduct parametric studies associated with the main variables affecting the fatigue crack growth of through-the-thickness cracks. A subject that can be analysed with the developed approach is the effect of the stress ratio on crack closure values. Figure 8 plots the ratio of the opening stress intensity factor (Ko) to the maximum stress intensity factor (Kmax) along the crack front for both materials. As shown, the plane stress curve represents the upper limit, while the plane strain curve represents the lower limit. The values of Ko/Kmax are between two limiting cases, and decrease with an increase in the stress ratio. In addition, at lower stress ratios, the differences between the maximum and minimum values of Ko/Kmax are higher for PMMA and tend to be closer for the aluminium alloy. This methodology can also be applied to conduct parametric studies associated with the main variables affecting the fatigue crack growth of through-the-thickness cracks. A subject that can be analysed with the developed approach is the effect of the stress ratio on crack closure values. Figure 8 plots the ratio of the opening stress intensity factor (K o ) to the maximum stress intensity factor (K max ) along the crack front for both materials. As shown, the plane stress curve represents the upper limit, while the plane strain curve represents the lower limit. The values of K o /K max are between two limiting cases, and decrease with an increase in the stress ratio. In addition, at lower stress ratios, the differences between the maximum and minimum values of K o /K max are higher for PMMA and tend to be closer for the aluminium alloy. Figure 9 plots the variation in the K o /K max ratio at the crack surface obtained from the presented 3D FE simulations against previously published relationships based on experimental tests that incorporated plasticity-induced crack closure. Notably, the results of the presented procedure agree well with the outcomes of the experimental and theoretical studies reported in the literature [1,16,[27][28][29]. The variation between the presented method and published data decreases with an increase in the R ratio, as the size of the reverse plasticity zone (or monotonic plastic zone) becomes smaller in the fatigue crack growth rates. These results provide further support to and validation of the numerical technique outlined in this paper. on crack closure values. Figure 8 plots the ratio of the opening stress intensity factor (Ko) to the maximum stress intensity factor (Kmax) along the crack front for both materials. As shown, the plane stress curve represents the upper limit, while the plane strain curve represents the lower limit. The values of Ko/Kmax are between two limiting cases, and decrease with an increase in the stress ratio. In addition, at lower stress ratios, the differences between the maximum and minimum values of Ko/Kmax are higher for PMMA and tend to be closer for the aluminium alloy.   Figure 9 plots the variation in the Ko/Kmax ratio at the crack surface obtained from the presented 3D FE simulations against previously published relationships based on experimental tests that incorporated plasticity-induced crack closure. Notably, the results of the presented procedure agree well with the outcomes of the experimental and theoretical studies reported in the literature [1,16,[27][28][29]. The variation between the presented method and published data decreases with an increase in the R ratio, as the size of the reverse plasticity zone (or monotonic plastic zone) becomes smaller in the fatigue crack growth rates. These results provide further support to and validation of the numerical technique outlined in this paper.

Conclusions
In this paper, new numerical modelling tools capable of simulating the crack shape development of through-the-thickness fatigue cracks in finite plates were presented. The proposed approaches assume a pre-defined crack front shape, and include plasticity-induced crack closure. The methodology was successfully tested for cracked rectangular cross-section geometries when subjected to Mode I loading. The following conclusions can be drawn: 1. The maximum stress intensity factor becomes relatively constant in the sufficiently thick plates and is equal to the value obtained for plane strain state conditions. The plane strain fatigue models (2D) may lead to inaccurate predictions when applied to the analysis of fatigue crack growth of thin structural plates; 2. The proposed methodology leads to satisfactory crack front predictions, either for ductile materials or brittle materials. Moreover, it is sensitive to the plate thickness, enabling good results for both thin and thicker geometries. In addition, it is capable

Conclusions
In this paper, new numerical modelling tools capable of simulating the crack shape development of through-the-thickness fatigue cracks in finite plates were presented. The proposed approaches assume a pre-defined crack front shape, and include plasticityinduced crack closure. The methodology was successfully tested for cracked rectangular cross-section geometries when subjected to Mode I loading. The following conclusions can be drawn:

1.
The maximum stress intensity factor becomes relatively constant in the sufficiently thick plates and is equal to the value obtained for plane strain state conditions. The plane strain fatigue models (2D) may lead to inaccurate predictions when applied to the analysis of fatigue crack growth of thin structural plates; 2.
The proposed methodology leads to satisfactory crack front predictions, either for ductile materials or brittle materials. Moreover, it is sensitive to the plate thickness, enabling good results for both thin and thicker geometries. In addition, it is capable of dealing with different stress ratios; 3.
The opening stress intensity factor increases with increasing values of stress ratio, maximum stress intensity factor, and distance from the centre of the crack. Predicted values obtained by the proposed methodology are quite close to those found in the literature for the same propagation conditions.
The comparison with experimental results is encouraging, and demonstrates the validity of the underlying assumptions: (1) the crack front shape intersects the free plate surface at the critical angle; ad (2) the local stress intensity factor can be considered as the fatigue crack driving force, which leads to the formation of the crack front shape under high cycling loading. The above assumptions might not be correct in the case of large plastic effects near the crack tip. In this case, the plasticity-induced crack closure, which is significantly different along the crack front, will be the one of the most influential factors affecting the crack front shape.
Future work will be directed to the application of the proposed methodology to more complex problems in terms of geometry, loading scenario, and crack shape configuration. Lastly, the simplicity and speed of calculation of the proposed approach, compared to the current numerical solutions used for the same purpose, make it quite attractive for simulating the fatigue crack growth, in both practical applications and parametric studies.

Data Availability Statement:
The data presented in this study are available from the corresponding author, upon reasonable request. The data are not publicly available due to ethical restrictions.

Conflicts of Interest:
The authors declare no conflict of interest.