Prediction of Fatigue Crack Growth in Metallic Specimens under Constant Amplitude Loading Using Virtual Crack Closure and Forman Model

: This paper considers the applicability of virtual crack closure technique (VCCT) for calculation of stress intensity factor range for crack propagation in standard metal specimen geometries with sharp through thickness cracks. To determine crack propagation rate and fatigue lifetime of a dynamically loaded metallic specimen, in addition to VCCT, standard Forman model was used. Values of stress intensity factor (SIF) ranges ∆ K for various crack lengths were calculated by VCCT and used in conjunction with material parameters available from several research papers. VCCT was chosen as a method of choice for the calculation of stress intensity factor of a crack as it is simple and relatively straightforward to implement. It is relatively easy for implementation on top of any ﬁnite element (FE) code and it does not require the use of any special ﬁnite elements. It is usually utilized for fracture analysis of brittle materials when plastic dissipation is negligible, i.e., plastic dissipation belongs to small-scale yielding due to low load on a structural element. Obtained results showed that the application of VCCT yields good results. Results for crack propagation rate and total lifetime for three test cases were compared to available experimental data and showed satisfactory correlation.


Introduction
One of the main tasks of advanced structural analysis is the consideration of the crack occurrence in structure [1], as well as structure lifetime estimation due to fatigue. Fatigue life in general consists of two phases: crack initiation period and crack growth period [2,3]. Related to that fact, estimation of the ability of the mentioned structures to resist growth of these cracks is of utmost importance for estimating the remaining life of the structure in its entirety. Some of the parameters that are commonly used to characterize and analyze the behavior of a structure containing a crack are energy release rate G, stress intensity factor K (SIF), Crack Tip Opening Displacement (COTD), and J -integral [4][5][6].
To predict the value of the mentioned crack parameters, analytical or numerical methods can be used.
One of the most common numerical method used to analyze complex structures is finite element analysis (FEA) [7,8]. In conjunction with additional calculations, it can be used to determine SIF for a crack and predict a lifetime of a structure in which that crack is contained. To predict the growth of an already existing crack in a structure, many different approaches can be used [2], including adaptive and moving mesh approaches [9][10][11][12]. The decision on which approach and fracture parameter will be used is dependent on a multitude of factors related to the material and geometry of a structure. One of the first constant amplitude loading crack propagation models that used SIF range to model crack growth was a model proposed by Paris, Gomez, and Anderson [13]. It is a simple model that has only two material parameters that are easily obtained by fitting the model to the available material data. The biggest shortcoming of that model was that it does not model the threshold region of crack growth or accelerated crack growth in the region where SIF range approaches values of critical stress intensity range K c . On top of that, it does not take into account the effects of stress ratio, so particular model parameters are required for a stress ratio of interest for a particular problem. These shortcomings of the Paris model are more or less accounted for by various proposed models that take into account some or all of these effects [2,4,14].
In this paper, crack propagation in metallic specimens containing sharp through thickness crack was analyzed using material parameters available in the literature, SIF range and Forman model in conjunction with the Barsom equation for crack propagation threshold [15]. SIF range for a dynamically loaded cracked element was calculated from results obtained from FEA in conjunction with the virtual crack closure technique (VCCT) [16,17]. VCCT was chosen since it can be easily implemented on top of any basic finite element (FE) code as it is a relatively simple FE post-processing method that does not require any special finite element types to determine the value of stress intensity factor at the crack tip. It is easy to implement in a wide range of basic 2D and 3D FE codes that use ordinary finite element types and linear-elastic analysis. Most frequently, it is used for fracture analysis of brittle materials that do not exhibit significant plasticity, i.e., for delamination of laminated composites. The groundwork for VCCT was first presented in the works of Hellen [18] and Parks [19], while the first application on finite element method (FEM) was done by Rybicki and Kanninen [17] and Rybicki et al. [20]. More recent books and research papers that cover VCCT theory and applications were published by many authors, including, but not limited to, Kuna [21], Krueger [16], Leski [22], Bonhomme et al. [23], Tavares et al. [24], etc.
The outline of this paper is as follows. Section 2 presents and overview of theoretical background used in the numerical algorithm presented in this paper, as well as some available analytical solutions for stress intensity factors, used for comparison and validation of VCCT results. Section 3 presents computer algorithm of the used numerical procedure in more detail. In Section 4, details on test cases used for validation are presented. Finally, in Section 5, results obtained from numerical algorithm are presented and compared to experimental data available from other published articles.

Crack Propagation Model
In real structures, fatigue crack growth in most cases occurs with variable amplitude time dependable loading. It is already well established that during these loadings various loading sequence effects can occur and significantly change the crack growth process. For example, crack tip plasticity during high load cycles in metals that exhibit significant plastic behavior can cause crack retardation effect, i.e., decreased crack propagation rate. This has an effect of the prolonged life of the structure made of such material [14,15]. Nevertheless, since it has been determined that sequence effects can often cancel each other out, and are often completely unpredictable, simple constant amplitude crack growth models can be used [2]. These simple models based on constant amplitude loading can be used to obtain results of sequences of constant amplitude loadings, representing real load sequence and neglecting mentioned sequence effects. The advantage is that this approach is relatively simple and gives a conservative estimate of life of a cracked structure.
The simplest, most widely used model is the Paris model. Since, as mentioned earlier, Paris model does not take into account many effects, a little more complex model can be used if enough material and load data are available. So, for the purpose of the numerical procedure presented in this paper, Forman model [2] will be used: When material constants for crack growth, C F and m F , and critical stress intensity K c are known, the Forman model can model the effect of load asymmetry and the effect of accelerated crack growth as the SIF range approaches the critical value of SIF range. In order to model the behavior of plates of various thicknesses and link them to more universal material property, critical stress intensity value was linked to plane-strain fracture toughness K Ic . Irwin presented a model to connect these two values [4]: and this model is used in a numerical procedure presented in this paper. This requires additional material property, yield stress σ y , to be known. Since it is always available for metals, this does not limit the applicability of the presented numerical procedure. To more accurately account for stopping of crack growth when SIF range approaches threshold value ∆K th , the threshold value of SIF range was modeled in such a way to account for load asymmetry with the use of stress ratio R, using Barsom Equation [15]: This requires threshold value of SIF for pulse loading ∆K th0 to be known. In case when ∆K th0 is not available, for stress ratio R ≥ 0.1 the threshold value can be computed using expression [25]:

Virtual Crack Closure Formulation for 4 and 8-Node Two-dimensional (2D) FEA
In numerical procedure presented in this paper, virtual crack closure technique (VCCT) was used to calculate SIF range based on FEA results for the models of the particular specimen containing sharp through thickness crack.
In the available literature, VCCT is also referred to as modified or virtual crack closure method. VCCT is derived from the crack closure method or two-step crack closure method that is based on Irwin's crack closure integral [16]. It assumes that crack extension from a + ∆a to a + 2∆a does not alter the state at the crack tip in a significant manner. Therefore, the displacements behind the crack tip at node B are approximately equal to the displacements behind the crack tip at node A when cracks tips are located at nodes C and B, respectively ( Figure 1a). In comparison to two-step crack closure method, which requires two finite element analyses, VCCT is faster to compute, i.e., more computationally efficient. The assumption from the VCCT enables us to calculate the work required to close the crack from one single FE analysis [16] and strain energy release rate G can be calculated by dividing the work with the created crack surface. When material constants for crack growth, CF and mF, and critical stress intensity Kc are known, the Forman model can model the effect of load asymmetry and the effect of accelerated crack growth as the SIF range approaches the critical value of SIF range. In order to model the behavior of plates of various thicknesses and link them to more universal material property, critical stress intensity value was linked to plane-strain fracture toughness KIc. Irwin presented a model to connect these two values [4]: and this model is used in a numerical procedure presented in this paper. This requires additional material property, yield stress σy, to be known. Since it is always available for metals, this does not limit the applicability of the presented numerical procedure. To more accurately account for stopping of crack growth when SIF range approaches threshold value ΔKth, the threshold value of SIF range was modeled in such a way to account for load asymmetry with the use of stress ratio R, using Barsom equation [15]: This requires threshold value of SIF for pulse loading ΔKth0 to be known. In case when ΔKth0 is not available, for stress ratio R ≥ 0.1 the threshold value can be computed using expression [25]:

Virtual Crack Closure Formulation for 4 and 8-Node two-Dimensional (2D) FEA
In numerical procedure presented in this paper, virtual crack closure technique (VCCT) was used to calculate SIF range based on FEA results for the models of the particular specimen containing sharp through thickness crack.
In the available literature, VCCT is also referred to as modified or virtual crack closure method. VCCT is derived from the crack closure method or two-step crack closure method that is based on Irwin's crack closure integral [16]. It assumes that crack extension from a + Δa to a + 2Δa does not alter the state at the crack tip in a significant manner. Therefore, the displacements behind the crack tip at node B are approximately equal to the displacements behind the crack tip at node A when cracks tips are located at nodes C and B, respectively ( Figure 1a). In comparison to two-step crack closure method, which requires two finite element analyses, VCCT is faster to compute, i.e., more computationally efficient. The assumption from the VCCT enables us to calculate the work required to close the crack from one single FE analysis [16] and strain energy release rate G can be calculated by dividing the work with the created crack surface. When 4-node two-dimensional (2D) elements are used to model a crack, strain energy release rate can be calculated by dividing the work required to close the crack ∆W with the newly created crack surface: where t is element thickness, ∆a is the length of the elements at the crack front, F xB and F yB are shear and opening forces at crack tip node B, and ∆u A and ∆v A are shear and opening displacements at a crack surface node A. In the expression (4), we can distinguish two parts of energy release rate. One is a result of mode I crack opening, and other is a result of mode II crack opening (in-plane shear mode), so we can write: where G I and G II are energy release rates for the corresponding crack opening modes.
In the case we model a crack with 8-node 2D elements (Figure 1b), the strain energy release rates for the corresponding crack opening modes can be calculated as: It has to be noted that, if VCCT is used for the simulation of crack propagation, a kinematically compatible node releasing scheme has to be used in order to obtain reliable results of strain energy release rate. If 8-node 2D elements are used, elementwise node release scheme has to be used, i.e., edge and midside nodes have to be released simultaneously [16,26].
For test cases in this paper, the calculation of SIF was necessary only for mode I crack opening, so only the calculation for the value of G I was used. In the case only small-scale yielding is present at the crack tip, SIF can be calculated from the value of known strain energy release rate G I [4] as: where K I is the value of SIF for mode I crack opening. Value denoted E* is different for the case of plane-stress or plane-strain condition. In the case of plane-stress condition, it is equal to elasticity modulus E, where in the case of plane-strain condition it can be calculated as where ν represents Poisson's ratio. If we calculate SIF only for mode I (opening mode), which is the main crack opening mode in examples used in this article, the following equations can be derived: Plane-strain : For the evaluation of the size of elements needed for the FE mesh on the crack path, known analytical solutions for SIF of the middle tension (M(T)) specimen, given by Irwin and Fedderson [14], are used: Fedderson : where normal stress σ is calculated as: where F represents specimen load, and t and w represent specimen thickness and specimen width, respectively.

Computer Algorithm
The numerical procedure presented in this paper was programmed in MATLAB and implemented in computer program VC-CGrow. As input data for this program, it is required to have available certain material data mentioned earlier, and results of crack surface displacements and crack tip forces acquired from FEA model. Based on this data, the program calculates SIF and crack propagation rate for a given load and given crack length. Since the program is based on VCCT and finite steps of crack propagation, SIF between crack length a and a + ∆a is linearly interpolated to give a more accurate value of SIF for the current crack length. Crack propagation rate and crack propagation threshold are calculated using Equation (1) and Equation (4), respectively. End of crack propagation and finite life (cycles at failure) are determined by achieving critical SIF value. It has to be noted that the implemented procedure does not include the effects of crack tip plasticity, crack closure effects or effects of overload. Those effects will be studied in future research and implemented in further development of numerical procedure presented in this article. A more detailed representation of the numerical procedure currently implemented in program VC-CGrow can be seen in Figures 2-4.

Test Cases and FEA Details
Considering the fact that the experimental data for crack propagation for standard test specimens was available in various published articles, standard test specimens were selected for the validation of the presented numerical procedure. These specimens also have known SIF solutions

Test Cases and FEA Details
Considering the fact that the experimental data for crack propagation for standard test specimens was available in various published articles, standard test specimens were selected for the validation of the presented numerical procedure. These specimens also have known SIF solutions and, as such, were suitable for the validation of SIF calculated by VCCT in this procedure. Therefore, the numerical procedure presented in this paper was tested on three different test cases of standard metal specimens containing sharp through-thickness cracks, subjected to constant amplitude loading. The first test case refers to the single edge cracked plate made of aluminum alloy 2024-T3 (test case 1) [27]. Second test case was middle tension specimen (M(T)), [28], made of martensitic structural steel 18G2A (Steel 1.0562, S355) [29] and third was compact tension specimen (C(T)), [28], made of welded pressure vessel steel A516 Grade 70 (Steel 1.0473, P355GH) with 45% overmatch weld [30]. All data regarding materials used in the preparation of test cases are presented in Table 1. Forman constants are either taken from appropriate articles or, if they were not available, obtained by fitting curves to available experimental data or Paris models. It should be noted that material fracture properties used in test case 3 were for weld material and not for the base metal. These values are denoted with ** in Table 1. In addition, since values of plane-strain fracture toughness were not available for steels used in this paper, values were estimated based on Charpy impact energy (CVN) for a given material using Roberts-Newton formula [4,31]: (18) and are denoted with * in Table 1. Specimen dimensions can be seen in Figures 5 and 6. Basic geometry data and loads of test specimens for particular test cases are shown in Table 2.
(a) (b) Figure 5. Sample geometry for: (a) test case 1, data from [27]; (b) test case 2, data from [29].        FEA was used to get the required data for various crack sizes, i.e., nodal forces and displacements needed for VCCT post-processing. Aluminum alloy specimens were modeled as a quarter model and steel specimens as a half model to alleviate some of the computational load. For the use of VCCT, crack front for all test cases was meshed with sufficiently small 8-noded 2D finite elements. Based on tests of VCCT calculated SIF accuracy in regard to FE size on the crack front, the size of 2 mm × 2 mm was chosen as a sufficiently small element size that gives reasonably accurate SIF results. A comparison of VCCT solutions for SIFs with various element sizes at the crack front and known analytical solutions can be seen in Figure 7. Elements smaller than 2 mm × 2 mm were unnecessary as they should not significantly contribute to the accuracy of calculated SIFs, only resulting in waste of available computer resources. Final FE meshes used in crack propagation tests are shown in Figure 8. In the end, results obtained from the numerical procedure integrated into the program VC-CGrow were compared with experimental data available in the literature [27,29,30].

Results
For test cases 1 and 2, calculated crack propagation rate was evaluated in comparison to experimental data available from the works of Mohanty, Verma, and Ray [27] and Skorupa and Skorupa [29]. For test case 1, it was found that there is a very good correlation with experimental data

Results
For test cases 1 and 2, calculated crack propagation rate was evaluated in comparison to experimental data available from the works of Mohanty, Verma, and Ray [27] and Skorupa and Skorupa [29]. For test case 1, it was found that there is a very good correlation with experimental data for load asymmetry ratio value of R = 0.1 (Figure 9). For test case 2, calculated crack propagation rate was compared to experimental data and was found to be very accurate for load asymmetry ratio of R = 0.15. From Figure 10, it can be seen that calculated values for crack propagation rate for load asymmetry ratio R = 0.5 were somewhat overestimated in comparison to experimental data available from Skorupa and Skorupa [29].
Metals 2020, 10, x FOR PEER REVIEW 11 of 16 was compared to experimental data and was found to be very accurate for load asymmetry ratio of R = 0.15. From Figure 10, it can be seen that calculated values for crack propagation rate for load asymmetry ratio R = 0.5 were somewhat overestimated in comparison to experimental data available from Skorupa and Skorupa [29].  Results for crack size versus the number of cycles, obtained with the presented numerical procedure, were compared to experimental data available in the before mentioned articles, for load asymmetry ratios given in Table 2. Results can be seen in Figures 11 and 12. Good correlation of crack growth curve with available experimental data was determined for both test cases. For test case 3, data for crack size vs. the number of cycles, calculated with numerical procedure implemented in program VC-CGrow, was compared to experimental data available from the article by Sarzosa et al. [30]. Good correlation with experimental data can also be seen in this case ( Figure 13).
Predicted lifetime of structural elements was compared to available data where possible. Experimental data for number of cycles to failure was available in the above-mentioned articles for test cases 1 and 3. Accuracy in comparison to data available from experiments was found to be 2.3% or less. Predicted number of cycles to failure for all test cases, as well as available experimental data, can be seen in Table 3.  was compared to experimental data and was found to be very accurate for load asymmetry ratio of R = 0.15. From Figure 10, it can be seen that calculated values for crack propagation rate for load asymmetry ratio R = 0.5 were somewhat overestimated in comparison to experimental data available from Skorupa and Skorupa [29].  Results for crack size versus the number of cycles, obtained with the presented numerical procedure, were compared to experimental data available in the before mentioned articles, for load asymmetry ratios given in Table 2. Results can be seen in Figures 11 and 12. Good correlation of crack growth curve with available experimental data was determined for both test cases. For test case 3, data for crack size vs. the number of cycles, calculated with numerical procedure implemented in program VC-CGrow, was compared to experimental data available from the article by Sarzosa et al. [30]. Good correlation with experimental data can also be seen in this case ( Figure 13).
Predicted lifetime of structural elements was compared to available data where possible. Experimental data for number of cycles to failure was available in the above-mentioned articles for test cases 1 and 3. Accuracy in comparison to data available from experiments was found to be 2.3% or less. Predicted number of cycles to failure for all test cases, as well as available experimental data, can be seen in Table 3.  Results for crack size versus the number of cycles, obtained with the presented numerical procedure, were compared to experimental data available in the before mentioned articles, for load asymmetry ratios given in Table 2. Results can be seen in Figures 11 and 12. Good correlation of crack growth curve with available experimental data was determined for both test cases. For test case 3, data for crack size vs. the number of cycles, calculated with numerical procedure implemented in program VC-CGrow, was compared to experimental data available from the article by Sarzosa et al. [30]. Good correlation with experimental data can also be seen in this case ( Figure 13).
Predicted lifetime of structural elements was compared to available data where possible. Experimental data for number of cycles to failure was available in the above-mentioned articles for test cases 1 and 3. Accuracy in comparison to data available from experiments was found to be 2.3% or less. Predicted number of cycles to failure for all test cases, as well as available experimental data, can be seen in Table 3.

Discussion
In this paper, the numerical procedure for calculating crack propagation was presented and results obtained with the use of the presented procedure was analyzed. The numerical procedure relies on FEA and VCCT method for the calculation of SIF range, and on available material data for the calculation of crack growth rate. Test cases included side cracked plate made from aluminum alloy, M(T) specimen made from 18G2A steel and a C(T) specimen made from A516 Grade70 welded steel, with crack progressing through weld material. These test cases were selected because of the availability of experimental data in published articles. Mentioned experimental data was used for the Figure 13. Crack propagation curve for test case 3 (compact tension (C(T)) specimen, A516 Grade70 welded steel).

Discussion
In this paper, the numerical procedure for calculating crack propagation was presented and results obtained with the use of the presented procedure was analyzed. The numerical procedure relies on FEA and VCCT method for the calculation of SIF range, and on available material data for the calculation of crack growth rate. Test cases included side cracked plate made from aluminum alloy, M(T) specimen made from 18G2A steel and a C(T) specimen made from A516 Grade70 welded steel, with crack progressing through weld material. These test cases were selected because of the availability of experimental data in published articles. Mentioned experimental data was used for the validation of the presented numerical procedure. Calculated crack propagation rate curves were compared to experimental data for two test cases, side cracked plate and M(T) specimen. In addition, crack propagation curves were validated against available experimental data for all test cases, and the calculated fatigue lifetime of cracked specimens was compared to experimental values for side cracked plate and C(T) specimen. Two-dimensional VCCT was shown to provide reliable results for SIF range needed in the calculation of crack growth rate. Calculated results for crack propagation and structural lifetime were compared to experimental data available in three different articles. Linear interpolation between SIF results for two VCCT steps did not have any noticeable impact on the calculated fatigue lifetime of the structure. In addition, negative impact of linear interpolation between two values of SIF calculated by VCCT cannot be seen or identified on any of the calculated crack propagation curves, given that the calculated curves closely follow the experimental data. Calculated fatigue lifetime was in line with experimental results. Since numerical procedure did not include various effects, such as plastic zone size, crack retardation effects, etc., further attempts are needed to include these effects in the existing numerical procedure. This should enable us to assess the influence of the mentioned effects on the prediction of fatigue lifetime for cracked specimens using the numerical procedure presented in this article.  Acknowledgments: This work has been financially supported by the University of Rijeka within the project: uniri-tehnic-18-42, "Investigation, analysis and modeling the behavior of structural elements stressed at room and high temperatures".

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. Nomenclature a 0 initial crack size a i crack size for cycle i a crack size a/w relative crack size C F , m F material constants for Forman model CVN Charpy impact energy da/dN crack propagation rate ∆a crack extension, length of the elements at the crack front ∆K stress intensity factor range for a load cycle ∆K th0 threshold stress intensity factor range for R = 0 ∆K thR threshold stress intensity factor range for specific load asymmetry ratio R ∆u shear displacement at crack surface node ∆v opening displacement at crack surface node ∆W work required to close the crack along one element side E elasticity modulus F specimen loading force F min minimum loading force in a cycle F max maximum loading force in a cycle F x shear force at the crack tip F y opening force at the crack tip G strain energy release rate G I strain energy release rate for crack opening mode I G II strain energy release rate for crack opening mode II K stress intensity factor K I stress intensity factor for crack opening mode I K II stress intensity factor for crack opening mode II K min minimum stress intensity factor in a load cycle K max maximum stress intensity factor in a load cycle K Ic critical stress intensity factor for plane-strain conditions, fracture toughness K c critical stress intensity factor N ef experimental results for cycles to failure N pf predicted cycles to failure R load asymmetry ratio σ normal stress σ UTS ultimate tensile strength σ min minimum normal stress in a cycle σ max maximum normal stress in a cycle σ y , σ 0.2 yield stress t specimen and element thickness u x-coordinate of crack surface node after load v y-coordinate of crack surface node after load w specimen width ν Poisson ratio