Abstract
During fatigue damage accumulation, cracks propagate through the material leading to catastrophic failure. As the cracks propagate, the natural frequency lowers, leading to a changing stress state. A new method has been developed where the damage accumulation rate is computed in the frequency domain using Linear Elastic Fracture Mechanics (LEFM), stress intensity, and the natural frequency. A finite element model was developed to predict the stress intensity and natural frequency during damage accumulation. Validation of the LEFM technique was done through comparison to experimental data. Reasonably good correlations between the FEM and the analytic model were achieved for the stress intensity and natural frequency.
1. Introduction
It is common knowledge that when a structure experiences vibration, a fatigue crack may eventually develop. When the crack becomes appreciably large, the natural frequency and mode shapes of the system will change [1]. Extensive health assessment work has been performed by utilizing the change in natural frequency as a damage indicator, prior to catastrophic failure. In order to make a meaningful life prediction of failure during random vibration, the change in natural frequency must be accounted for. It is possible to account for the change in the natural frequency using the virtual crack technique within the finite element method (FEM).
The growing role of virtual crack simulation during the design process of mechanical components induces engineers to improve the progress of life prediction using this technique. Industries that wish to improve the reliability of products by simulation of the dynamic response often use the FEM. This approach provides a mathematically stable environment and allows for modeling structures with complex geometries, both of which are essential for industrial application.
One may simulate the level of structural weakness using the virtual crack technique to predict the remaining useful life. For example, despite the presence of cracks in the wings of a typical aircraft, the aircraft may continue to be flight-worthy as long as the cracks do not exceed the damage tolerance. Nonetheless, calculating the remaining life of structures continues to be a challenge, especially when the vibration loading is random. Many methods have been developed for fatigue life evaluation based on a representation of the stress state, both in the time and the frequency domains [2–8]. Approaches using the power spectral density (PSD) and the root-mean-square acceleration, , were developed to obtain an equivalent maximum uniaxial stress or von Mises stress value to estimate the life of critical structures [4, 9]. Reference [3] contains an excellent survey of various equivalent fatigue damage models. Unfortunately, many of these models may not be fully suitable for a general use of the virtual approach to fatigue, especially when the frequency shifting or random dynamic loading is present.
The solution to the problem can be determined using the rate of frequency change (RFC) model [10]. This model characterizes multiaxial dynamic behavior of the mechanical system in the frequency domain. The stress intensity and natural frequency for each crack extension increment are required to fully predict time-to-failure. The RFC model and FEM are used to predict the life of structures experiencing random vibration. This method utilizes the frequency domain to model the crack growth. This approach incorporates the change in natural frequency and may offer reasonable predictions of the time-to-failure. Regardless of the structure complexity, the FEM feeds the RFC model with an updated natural frequency, , and a stress intensity factor, , for each virtual crack extension [10]. The model employs Linear Elastic Fracture Mechanics (LEFM) for fatigue crack propagation and accounts for the frequency shifting. The FEM development and its implementation into the RFC model are discussed in this paper. The full details of the RFC and experimental results are presented in [10].
The major advantage of the RFC approach is the ability to estimate time-to-failure in the frequency domain, where only the input power spectral density and damping factor are required. Monitoring the changes in the natural frequency and the stress intensity factor, regardless of the complexity of the geometry, may lead to a reasonable estimate of the remaining life of weakened structures. Experiments are conducted in parallel to validate the FEM models. Correlation between the analytical model and both FEA and experimental results is achieved.
2. Finite Element Method Development
In this study, single-edge -notched cantilever beam specimens were tested and modeled to validate the RFC-FEM coupled model, as shown in Figures 1 and 2. The beam specimens were fabricated from cold rolled 1018 steel. The cantilevered beam was 2.50 in long with a cross-sectional area of 1.250 × 0.1875 in2, as shown in Figure 2. The -notch detailed dimensions are depicted in Figure 3. The test sample was designed such that the first and second modes were closely spaced. The first bending mode in the vertical direction was 330 Hz. In the transverse direction, the first bending mode was also 330 Hz. Detailed explanation of this design is provided in [10]. The beam specimen was clamped by a four-bolt fixture, as shown in Figure 1. Each bolt was torqued to 30 ft-lb. Final failure was defined to occur when the beam tip made contact with a limit-bar located approximately 2.5 cm below the tip of the beam, as shown in Figure 1.
(a)
(b)
Notches in structures increase the localized stress concentration, which decreases the maximum load the structures can sustain. Hence, a criterion to evaluate the maximum load that a notched component can sustain is vitally important. The simplest and most frequent geometries are -notches. They are commonly observed in test samples and in notched structural components [11–13]. Studies show that the critical notch intensity factor as a function of notch angle can be used as a fracture criterion, provided that the -notches are sharp and plasticity is contained [11, 12].
In order to assess the validity of the RFC-FEM model, stress intensity factors and change in the frequency data were recorded for beams exposed to stationary, Gaussian, random vibration inputs. The input level was 0.0349 G2/Hz from 20 to 2000 Hz. The FEM was utilized to predict the beam modal response, track the natural frequency shift, and predict the stress intensity factor for mode , . Details of the particular vibration environments are discussed in [10].
The characterizes the local mode of the crack tip stress field in linear elastic material for mode loading condition where the principle load is applied normal to the crack plane. Mode of the crack surface displacement is considered to be the opening of the crack. Modes and are the sliding and tearing of the crack, respectively. Based on the FEM modal response analysis, the beam bending is the dominating mode. Therefore, only mode was considered in this study. The depends on the applied remote stress, the geometries of the structure, and the crack size as follows: where is the remote stress, is the crack size, and is the structure geometric factor. For simple structures, can be obtained from handbooks [14]. For complex structures, it is difficult to calculate analytically. Numerical methods such as FEM can be utilized to predict . Therefore, was determined from FEM using the -contour integral approach and updated as the virtual crack increased. This approach provided accurate results with surprisingly coarse mesh [13].
In designing the FEM mesh for fracture mechanics problems, the common focus is on the tip of the crack. The crack tip is a singularity point where the stress field becomes mathematically infinite. If the cracked region is modeled with conventional polynomial-based FEM, the mesh must be exceptionally dense around the crack tip. This approach may not be feasible and, depending on the problem complexity, may be costly and time consuming.
In FEM fracture mechanics, it is recommended to use 9-node biquadratic Lagrangian elements for two-dimensional problems and 27-node triquadratic Lagrangian elements in three-dimensional problems [13, 15]. The 8-node two-dimensional and 20-node three-dimensional elements are also common in crack problems. In this study, conventional 8-node fully elastic elements were used for the beam (away from the crack) and the fixture. For the area near the crack tip, 20-node elements were used and will be discussed in detail later. Multipoint constraints beam elements were utilized to model the torqued bolts used to clamp the cantilever beam, as shown in Figure 1. The preload value was applied to each bolt in the FEM model to produce accurate results.
The stresses near the crack tip are characterized by and expressed as follows [15]: For plane stress is zero and for plane strain is equal to . The displacement near the crack tip is calculated as follows [15]:where is the shear modulus. and are shown in Figure 4. It can be seen from the relationships above that the stress varies inversely with and the displacement varies proportionally with . The presence of a singularity of the stress at the tip was attained when approached zero.
During fatigue crack growth, both plane strain and plane stress will be present along the crack front [16]. However, the conditions ahead of a crack are neither plane stress nor plane strain but require treatment in three-dimensional case [13]. The material near the crack tip is at higher stresses than the surrounding material. Since there is no stress normal to the free surface, the material on the surface is in a state of plane stress. At the midplane of the crack, plane strain conditions exist and is .
In this investigation, the behavior of the stress and displacement near the crack tip were modeled using 20-node hexahedral fully elastic elements degenerating down to wedges. The wedge element was identical to the 20-node quadratic hexahedral element, except the nodes facing the crack site collapsed to form a wedge, as shown in Figure 5. Each of the red three-node edges was collapsed to one corner location, as shown in Figure 5, while the two middle nodes collapsed into one center location. The collapsed nodes were not merged into one node but tied together and became coincident. Subsequently, the collapsed nodes formed one line, which was the crack front. The middle nodes located at the faces orthogonal to the crack front were shifted by a quarter of the edge length closer to the crack tip. This modification enhanced the numerical accuracy without requiring significant mesh refinement to capture the crack tip stress field [13]. The final result was a “spider-web” meshing configuration at the crack tip region, as shown in Figure 6, which is also called the “spider-web” meshing technique.
As mentioned above, the stress state at the crack tip in a linear elastic material exhibits a mathematical singularity: where for plane stress, and for plane strain, where is the yield stress. The variable in this case is a first-order correction of the plastic zone size. In reality, the crack tip is surrounded by a zone where plastic deformation and material damage may occur. As can be seen from the FEM stress contour in Figure 7, the crack tip caused stress concentrations, where the stress gradient became larger as the crack tip was approached. The mesh was refined in the vicinity of the crack tip to extract accurate stresses. It is important to point out that LEFM is not accurate inside the plastic zone. However, LEFM may provide accurate results provided the plastic zone is small enough. It was necessary to take advantage of Paris law to predict the crack growth as a function of time. The crack growth rate can be calculated as follows [3]: where is the number of cycles, and are constant fatigue material properties, and and are the maximum and minimum stress intensity factors, respectively. To apply the Paris law, the stress intensity factor was calculated by means of the -contour integral first introduced by Rice [17]. The -contour integral is usually used in rate-independent quasistatic fracture analysis to characterize the energy release associated with crack growth. It can be related to the stress intensity factor if the material response is linear. The -contour integral is formulated as a path-independent line integral with a value equal to the decrease in potential energy per increment of crack extension in the material. The path independence implies that can be seen as a measure of the intensity of stresses and strains at the tip of the notch and crack [18]. Through the -contour integral,it is possible to calculate accurate values with coarse meshes without maintaining precise local stress values. It is important to point out that the -contour integral should be independent of the domain used; however, the -contour integral approximations from different rings in the web-mesh may vary due to the approximate nature of the finite element solution [19]. The strong variation in these estimates is commonly called domain or contour dependence. Therefore, the spider-web mesh was refined in the crack region with 18 contours, as shown in Figures 7 and 8. Although 18 contours might be excessive for this analysis, computing resources and the simplicity of the problem allowed this excess. Further study into the sensitivity of the results to the number of contours would be warranted. The crack front-line consisted of 21 nodes, which produced a 20-element-thick spider-web mesh, as shown in Figure 7.
The crack front was defined on the collapsed plane of the wedge elements, which is called the “edge plane.” The edge plane is orthogonal to the crack direction. Figure 7 is a magnified image of the -notch spider-web meshed region, where the center of the spider-web mesh is the crack front. Although not precisely true, the crack front is assumed to be a straight line. In the FEM model, the direction of the virtual crack extension was assumed to be in the downward vertical direction, which is denoted as the vector in Figure 8. The latter assumption is consistent with the experimental results obtained by Paulus et al. [10]. The virtual crack extension was updated in the FEM model manually based on Paulus et al. experimental results. Four virtual cracks extensions (, 0.0206, 0.0311, and 0.0516 in) were analyzed, as shown in Figure 9. The refined spider-web mesh at the crack tip was maintained for every crack extension.
(a)
(b)
(c)
(d)
3. Rate of Frequency Change Model
This section describes the implementation of the RFC model. Paulus et al. provided the derivation of the RFC model in a separate paper [10]. The model implementation is summarized in flowchart in Figure 10. The first step was to obtain the base excitation PSD . For an equivalent single degree of freedom system, the root-mean-square relative displacement can be calculated as follows [20]: where is the frequency and is the natural frequency for each crack extension. The is the equivalent viscous damping coefficient. Once the relative displacement is known for a given natural frequency, the root-mean-square far-field stress range can be calculated using equivalent strain energy as follows: where is the natural frequency for an uncracked uniform cantilever beam and is the Young elastic modulus. In this study, is the beam half height of the rectangular cross-section and is a simplified parameter that can be obtained from the first mode eigenvalue solution of a uniform uncracked cantilevered beam [10].
During failure, the natural frequency will change as the crack grows. Paulus et al. provided a general model for calculating the rate of natural frequency change, RFC, as follows [10]: where is total change in natural frequency and is an equivalent damage constant based on Rayleigh’s approximation. The constant represents a correction factor for the statistical distribution of the random vibration response, which is assumed to be narrowband, stationary, and Gaussian. Determination of can be accomplished statistically as discussed in detail by Paulus et al. [10]. Rayleigh’s approximation can be applied when the response is narrowband, as was the case for this analysis. The RFC model relates the crack growth rate to the change in the natural frequency of the structure. The major advantage of this model is that time-to-failure (TTF) estimation can be conducted in the frequency domain. However, the focus of this paper is on the rate of change in the natural frequency and not TTF prediction. The FEM allows the model to be extended to structures with complex geometry, where it may not be possible to compute the equivalent strain energy or the shift in the natural frequency. Figure 11 illustrates how the FEM can be combined with the RFC model to provide an approximate result instead of a closed-form solution.
Finally, the cracked structure TTF is obtained by integrating the inverse of the RFC with respect to the change in the natural frequency as follows: where and are the natural frequencies of the cracked structure at the initial life (when the crack is discovered) and end of life (when failure is reached), respectively.
4. Results
Computation of natural frequencies and the stress intensity factors as a function of the crack length were performed experimentally and numerically. In the first computation, the modal response analysis was performed for zero crack extension, . The experimental and numerical first bending mode values in the vertical direction were 330 and 351 Hz, respectively, as shown in Figure 12. In the transverse direction, the experimental and numerical first bending mode results were also 330 and 351 Hz, respectively. This analysis was repeated at varying crack depths from 0.010 to 0.114 in, as shown in Figure 13. Although a slight bias is apparent, the FEM modal analysis shows a linear relationship between the natural frequency and the crack depth, as shown in Figure 13. These results are in agreement with Paulus et al. assumption that the natural frequency varies linearly as a function of the crack depth [10]. The modal analysis also confirmed that the FEM model is correlated with Paulus et al. experimental results. Therefore, when analytical tools are unavailable or the structure geometry is complex, the rate of frequency shift from the FEM tool can be combined with the RFC model to calculate the remaining life, as shown in Figure 11.
The second computational analysis consisted of predicting the stress intensity factors for crack extensions: , 0.0206, 0.0310, and 0.0516 in, as shown in Figure 14. After , the beam started to experience excessive plasticity, which meant the use of the stress intensity factor was no longer possible, but the experiment was conducted to full failure. The stress intensity factor was calculated by applying a static load at the tip of the beam that was equivalent to the rms dynamic load obtained from the experiments performed by Paulus et al. [10]. For the zero crack case, both FEM and analytical methods were used to calculate the for each node along the crack front-line, which consisted of 21 nodes. The analytic stress intensity was determined using handbook calculations as described in [10]. The distribution of the stress across the cross-section was required to determine the stress intensity factor. This distribution of stress was determined using the stress concentration factor as determined by FEM. An average was calculated for each node on the crack front-line. Eighteen contours (rings) in the spider-web mesh were included in the calculations to ensure that the -contour integral was not domain dependent (or contour dependent). This is because the stress intensity factors have the same domain dependence features as the -contour integral. Numerical tests suggest that the estimate from the first ring of elements adjoining the crack front does not provide a high accuracy result [21].
(a)
(b)
(c)
(d)
The values obtained from FEM model were consistent with the analytical model, especially towards the center of the crack front-line for in, as shown in Figure 14. There was a slight deviation at the surface nodes. For the other crack extensions, the FEM and analytical values were in close proximity to each other and maintained the same trend. Ultimately, the RFC model required a root-mean-square stress intensity factor, , as shown in Figure 11. Thus, is calculated for both analytical and FEM analyses for each crack extension by averaging the nodal values and listed in Table 1. The values from the FEM and the analytical model are plotted as a function of the crack extension, as depicted in Figure 15. It is inferred from Figure 15 that the FEM values are close to the values obtained from analytical model. The increased as the shift in the natural frequency decreased (or as the crack depth increased).
5. Conclusion
This study provided a general virtual model that combined FEM with RFC. This model was used to predict change in the natural frequency, thus estimating fatigue life, using only frequency domain information. Execution of the model required only the input power spectral density, damping factor, and material properties. Integrating the FEM and the RFC model allows the model to be extrapolated to more complex geometries for which closed stress intensity values are not available.
The FEM further demonstrated the validity of the assumption of a linear relationship between crack depth and natural frequency over limited crack lengths. Although this assumption was used for the closed-form solution of the RFC model, use of the FEM allows this assumption to be relaxed in lieu of FEM results. The stress intensity quantities as a function of the crack growth are extracted from the FEM using the -contour integral. Reasonably good correlations between the FEM and the analytic model are achieved for the stress intensity and natural frequency. From the model, it can be deduced that the average stress intensity factor increased as the natural frequency decreased. Additional work is needed to conduct experimental and computational time-to-failure comparison.
The proposed model in this paper may be extended to accelerated life testing, virtual qualification, and reliability assessment. It can be used as a degradation model to analyze the relative severities of complex structures under harsh vibration environments. No explicit knowledge of the time history is needed. Thus, structural engineers could harness the flexibility of this model to reasonably predict the life-cycle when the only input is a PSD vibration profile. This approach may reduce the computation time and cost required to run a fully explicit FEM analysis.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.