An improved crack tip location algorithm using the principles of thermoelastic stress analysis

Abstract

A novel method of determining the crack tip location from a thermoelastic quadrature signal is presented. The method is utilised for crack tip locations within complex stress fields, namely within fastened aircraft lap joints. Coupled structural-thermal finite element modelling was undertaken to investigate the thermal response field around the crack tip location and develop the algorithmic principles. Experimental validation of the crack tip location was conducted using established crack mouth compliance techniques and optical measurements. The crack tip finding algorithm used the location of the maximum spatial gradient of the thermal field in the direction of crack growth. The method showed good accuracy when compared to traditional methods. Resultant crack growth rates were further verified using quantitative fractography.

Nomenclature

$a$

crack length ( ${\rm{mm}}$ )

${C_p}$

specific heat ( ${\rm{J}}/{\rm{kgK}})$

$d$

distance between max and min ( ${\rm{mm}}$ )

$E$

modulus of elasticity ( ${\rm{GPa}})$

$k$

thermal conductivity ( ${\rm{W}}/{\rm{mK}})$

$\Delta T$

change in temperature ( ${\rm{K}}$ )

${T_0}$

absolute temperature ( ${\rm{K}}$ )

$R$

load ratio

$\alpha $

coefficient of thermal expansion ( ${{\rm{K}}^{ - 1}}$ )

$\nu $

Poisson’s ratio

$\rho $

density ( ${\rm{kg}}/{{\rm{m}}^3}$ )

$\Delta {\sigma _{sum}}$

change in stress sum ( ${\rm{MPa}}$ )

${\sigma _{yy}}$

stress component in the loading direction ( ${\rm{MPa}}$ )

${\sigma _{xx}}$

stress component transverse to the loading direction ( ${\rm{MPa}}$ )

1.0 Introduction

Thermoelastic stress analysis (TSA) is showing great capability as an experimental tool for fatigue crack growth analysis [1–5]. The method is non-contact and allows for full-field stress measurement, meaning a large area of a structure can be assessed at once (compared to area-averaged point-data instruments such as strain gauges). In addition, TSA can be used for tracking of fatigue crack growth and evaluation of crack tip parameters such as the stress intensity factor (SIF) [1–4]. TSA can also be used to characterise and locate crack tips when there is no free-edge crack mouth (i.e. no provisions for a crack mouth opening displacement gauge), and is capable of observing multiple cracks at once for multiple site damage (MSD) assessments [6, 7]. The availability of compact and low-cost microbolometers further adds to the benefits of this technique. Use of a damage tolerance philosophy for certification of aircraft is commonplace. To ensure continued airworthiness, accurate fatigue crack predictions and the validity of the data supporting the predictions is paramount. Thus, continued development of TSA-based crack evaluation techniques is vital to realise the full benefits and potential cost savings for the aircraft industry.

Studies have shown the potential for using features in the TSA output signal(s) as references for the location of the crack tip. Thermal imaging during cyclic loading provides two signal measurements, which are a component in-phase with the remote loading and an out-of-phase component (or alternatively, a magnitude and phase angle). One method has been to take the point where the phase plot changes from positive to negative as an estimate for the crack tip [2]. The reasoning for this approach is usually based around arguments that the phase shift is due to plastic deformation of the material ahead of the crack tip or due to closure of the crack faces behind the crack tip. Other studies have suggested that the crack tip is located at the local minimum of the phase field [3, 8]. Several studies have taken the additional step to implement automated crack tip detection and tracking during live TSA measurement [4, 5, 8]. Thatcher et al. [4] utilise the point of maximum temperature response as the location of the crack tip. The work by Rajic and Brooks [5] uses the local minimum in the quadrature signal, rather than the phase, and continually re-positions the thermal camera with linear stages to centre and track the crack tip. Multiple past studies have noted that estimates for the crack tip location using phase/quadrature are consistently ahead of the actual crack tip position [3, 5]. The past research highlights the advantages of using TSA for crack tracking and evaluation, however, further work is required to understand and refine the process.

This paper discusses the ongoing development of a TSA-based crack growth evaluation method, including application to fastened lap joint structures. First presented is a detailed analysis of the thermoelastic stress field features surrounding a crack tip using coupled-field finite element (FE) modelling and TSA imaging. From this, the thermoelastic response features are quantified and the results used as an interpretation tool to make a better estimation of the crack tip location with the TSA field. The method is validated using standard test samples and then applied to a more complex situation, which is a fastened lap joint. The method presented here will be fundamental for ongoing work in evaluation of mixed-mode SIF and multi-site damage using TSA.

2.0 Experimental methods

2.1 Sample overview

All test samples were manufactured from 7075-T6 alclad aluminium sheet (AMS QQ-A-250/13). Static material properties and fatigue behaviour for this material are well established in the literature. For the validation case, a traditional compact tension geometry is used. The specific samples here are disk-shaped compact tension (DCT) as defined in ASTM E399. The samples were waterjet cut to size from 4.826mm (0.190”) thick material and a chevron notch was machined using a 0.5mm slitting saw. Specific dimensions of the validation samples are shown in Fig. 1. A total of four samples were tested with various cyclic load amplitudes to measure the crack tip location over a range of stress intensity factors. Samples were coated with a matte paint to minimise any effects from external infra-red reflections.

Figure 1. Schematic of DCT sample (units in mm).

2.2 Test machine, instrumentation and data acquisition

Samples were cyclically tested in a servo-hydraulic universal testing machine (Instron 1342) at 10Hz under constant amplitude load (R=0.1) until the sample catastrophically failed. The test machine has a maximum capacity of 250kN and was fitted with a 50kN load cell for the current testing. Load and actuator displacement data was periodically recorded though the machine data acquisition system every 5k cycles, for a duration of at least 50 cycles, and supplemented with custom instrumentation as follows.

An independent crack tip location measurement was included for the DCT samples in the form of a crack mouth opening displacement (CMOD) gauge. There is a well-established analytical relation between the crack length and the crack mouth opening displacement (compliance) for a DCT sample defined in ASTM E399. Strain gauges were included on all samples for validation of the FE modelling, comparing the measured strain to the FE strain field. The strain gauge data is used for the stress/strain calibration of the thermal imaging data, such that the FE thermal response can be experimentally validated. Calibrated fields are also beneficial for subsequent stress intensity factor evaluation, but that work is not discussed here. The DCT samples had FLA-2-23 strain gauges applied at several locations along the sample (constant 4mm above the crack path), allowing for a range of datapoints to be obtained for validation of the FE model over range of crack lengths. Strain gauge and CMOD gauge data was collected periodically every 5k cycles for a duration of at least 50 cycles.

For TSA data acquisition, a Xenics Gobi 640 microbolometer detector (temperature sensitivity <50mK) was used in conjunction with a laptop running the DSTG “MiTE” software [1]. The microbolometer was also mounted on a set of XY linear stages for precise field-of-view alignment of the sample being tested. The microbolometer lens was attached to the sensor body via a range of focal length extension tubes to further tune in the field-of-view to a desired area. For the DCT samples the resultant spatial resolution of the acquired TSA field was 200µm/pixel. The MiTE software also requires a load-correlated signal for TSA correlation. To achieve this, the direct load control signal from the Instron controller was run through a NI USB-6001 DAQ into the laptop running MiTE. TSA data was acquired periodically in intervals of approximately 5k cycles for an acquisition period of ∼500 cycles.

3.0 Finite element modelling

3.1 Modelling methodology

FE modelling of the DCT samples was utilised to understand the in-phase and out-of-phase (quadrature) thermal response field. This includes the associated features of the thermal response and the relationship with the physical location of the crack tip. Finite element modelling was undertaken using ANSYS 2020R2. The model uses a one-way structural-thermal coupling methodology. The structural analysis is solved first, then the resultant stress field is used to provide a heat generation input via the thermoelastic equation, which is passed to the thermal harmonic model. This allows for the heat-dissipation effects of different frequencies to be appropriately captured. The model surfaces are assumed to be perfectly insulated (no heat loss). The material properties for the 7075-T6 sheet are given in Table 1 and were taken from in house testing and published data (MMPDS-15). Both mechanical and thermal properties were required for the structural-thermal coupled analysis to simulate the thermoelastic effect. For the modelling presented in this paper, material plasticity is ignored as additional modelling (details omitted) showed limited effect of plasticity on the thermally evaluated stress field for the current loading regime. The DCT sample is modelled with 2D planar elements. The meshed model is shown in Fig. 2. A high refinement level is used at the crack tip (20µm elements) to capture the local stress gradients. Convergence was checked to verify the mesh refinement by assessment of the resultant strain field.

Table 1. Material properties for the FE model

Figure 2. Geometry and mesh, showing high mesh density at the crack tip.

3.2 Validation of remote strain field using strain gauges

Strain gauges were applied to the back of one of the DCT samples, and the cyclic minimum and maximum strain recorded at specific intervals throughout the duration of a test. The resulting strain values are plotted in Fig. 3 against the crack length. The strain gauge was located directly above the crack path at a crack length of 35mm. Hence the peak in measured strain is just prior to that point, followed by a swift decline for data points where the tip of the crack is beyond the location of the gauge. The vertical strain data was extracted from independent FE models with manually set crack lengths as an average over the gauge-area in the same location as the experimental strain gauge to create the FE strain gauge dataset. The correlation between the results suggests a physically realistic strain distribution is being achieved in the FE model. Thus, the crack tip strain field can be considered trustworthy for the purposes of this study.

Figure 3. FE model validation using strain gauges.

3.3 Validation of crack tip stress field using TSA

Cyclic loading of the DCT sample produces a thermal field that is measured by the TSA camera. For a linear-elastic, isotropic, and homogenous material, this vector temperature field (i.e. amplitude and phase) is converted to an equivalent stress field, using a simplified form of the heat transfer equation (assuming no heat transfer at surfaces):

(1) \begin{equation}\Delta {\sigma _{sum}}{\rm{\;}} = {\rm{\;}}\Delta \left( {{\sigma _{xx}}{\rm{\;}} + {\rm{\;}}{\sigma _{yy}}} \right){\rm{\;}} = {\rm{\;}} - \rho {C_{}}\Delta T/\alpha {T_0}\end{equation}

For the current work, the in-phase and quadrature components of the equivalent stress field are utilised. This means Equation (1) evaluated in-phase with the load signal, and 90 degrees out-of-phase from the load signal. Figures 4 and 5 show the equivalent stress sum field as evaluated from the thermal response in the near-tip region for a crack length of 29.6mm. The results are given as the in-phase and quadrature components of the equivalent stress sum for both the FE prediction (Fig. 4) and corresponding experimental measurements (Fig. 5). For reference, the crack tip is marked with a black cross. The in-phase plot shows a stress concentration delta (min load to max load) at the crack tip, while the quadrature plot contains features that are a product of the significant heat conduction caused by the corresponding high stress gradients. The features in the quadrature field reduce in size with increased loading frequency, and the in-phase response approaches the linear-elastic structural solution. As such, the accuracy of the TSA data is dependent on the local strain gradients and loading frequency. The FE results and experimental TSA data show consistent thermal response and feature distribution.

Figure 4. FE stress sum at 10Hz evaluated from the thermal response (a) in-phase and (b) quadrature.

Figure 5. Experimental TSA for the thermal response field at 10Hz: (a) in-phase and (b) quadrature.

Taking a path plot through the crack tip aligned with the direction of crack growth provides the data shown in Fig. 6 (in-phase and quadrature). Calibration of the TSA measurements was achieved by comparison of the experimental data with FE and the strain gauge values at several reference locations. The experimental calibration involves a scale factor that accounts for ambient conditions, the material properties, sample geometry and experimental setup. Details of the calibration are omitted here as only the location of the thermal features are required for finding the crack tip. In this case the FE simulations and experiments quantitatively match well. As such, the FE model of the DCT sample, including the near-tip stress field is considered valid for ongoing algorithm development purposes.

Figure 6. FE and experimental (calibrated TSA) stress sum line plot: (a) in phase field and (b) quadrature field.

4.0 Crack tip location algorithm development and validation

4.1 Quadrature feature distribution (QFD) method

The identification of crack tip locations from TSA is performed as a post-process operation, with the images acquired using MiTE being fed into a processing script. The first step in both methods is to scale the TSA image spatially. This is achieved by selecting reference locations at a known distance on the samples. The known distance can then be used to determine an overall spatial scale factor (i.e. mm/pixel). The initial algorithm employed used a constant distribution (proportion) between the minimum and maximum locations in the quadrature plot as the location of the crack tip. Figure 7 shows the quadrature response of the DCT samples, over a range of crack lengths, extracted as path plots in the direction of crack growth and centred at the crack tip. The locations of the minimum and maximum quadrature response values are constant with respect to the crack tip. When considering the total distance (d) between the two identifiable features, the FE modelling shows the crack tip consistently lies 0.1184d, behind the quadrature minimum, and 0.8816d ahead of the quadrature maximum.

Figure 7. FE quadrature path plots.

Numerically, this method is stable and reliable over a range of SIFs. However, a strong sensitivity to the exact location of the minimum and maximum quadrature values was observed. In the numerical case, the minimums and maximum are well defined, and can be programmatically found using simple min/max searches (see Fig. 4(b)). Difficulty arises when applying this method to experimental data, as the resolution of the sensor, general signal noise and the effects of surrounding structure/free edges often detract from the theoretical ideal stress field (see Fig. 5(b)). As such, when applying the QFD method to experimental data, an additional layer of curve fitting (surface fitting) was applied to minimise TSA noise. This allowed honing in on the best approximation of the ideal maximum and minimum locations, rather than just a single maximum or minimum pixel. Figure 8 illustrates the additional curve fitting applied to experimental data.

Figure 8. Feature fitting process (a) raw data, (b) cosine surface fit, (c) original and fitted data overlayed.

4.2 Maximum gradient method

For the revised maximum gradient method, more effort was dedicated to programmatically isolating the thermal features exclusively associated with the crack tip. There are occasionally hot-spots present in the thermal field that are unrelated to the crack tip and associated stress field. Common examples of these undesirable hot-spots are any free edges (heads of fasteners or edges of cracks), areas of high friction (under fastener heads or at fixture/grip locations), or areas of high reflection (chips/cracks in paint). The full acquired field is subjected to a progressive masking process prior to applying the final crack tip location algorithm (as shown in Fig. 9).

Figure 9. (a) Full quadrature response, (b, c) progressive masking, (d) spatial derivative and absolute maximum.

Once sufficient masking has been applied, the final crack tip location algorithm utilises the notion of a stress singularity at the crack tip. The stress singularity generates a concentration region for the in-phase thermal response. Due to associated heat conduction from the hotspot, this also leads to a concentration in the quadrature response. When the spatial derivative (in the direction of crack propagation) of the quadrature is calculated, the stress singularity will present as the point at which the magnitude of the derivative is maximised. This concept has been validated using the same FE data from the QFD method development, simply differentiating with respect to distance. This produces the data sets shown in Fig. 10. The plot clearly illustrates the divergence of the derivative of the quadrature field at exactly the physical location of the crack tip. It is expected that this method would be employable for both the in-phase and quadrature components of the TSA acquisition. However, previous work conducted in-house has shown that applying the method to the quadrature field yields the most accurate crack tip location. When applying the method to experimental data, the derivative of the entire masked field is taken in the direction of the crack growth. Then, simple min/max searches can be employed to locate the location of the maximum gradient (crack tip).

Figure 10. Differentiated quadrature field path plots, showing divergence at the crack tip.

4.3 Results

In Fig. 11, the results of applying the developed crack tip location methods to an individual DCT sample are presented along with the crack mouth compliance measurement (ASTM E 399 A6.4). It also shows results using previously published crack tip finding methods such as the phase minimum, phase zero-crossover, and quadrature minimum [2, 3 and 4]. The initial quadrature feature distribution method shows better correlation with the CMOD curve than the previously published methods, with a maximum relative error of 1.5% compared to the CMOD data. Similar observations were made for all other samples. The revised maximum gradient method is computationally more efficient, and results in a lower relative error when compared to the CMOD derived crack length data (max error of 0.75%). The QFD method is implementable for any crack length, provided the minimum and maximum locations of the quadrature plot are visible. The maximum gradient method can also be applied to any thermal field where the crack tip is located within the field-of-view. This includes situations with multiple cracks and near-fastener locations.

Figure 11. (a) Crack length measurements, (b) relative error.

5.0 Application to fastened joint samples

5.1 Specimen design

To demonstrate the applicability of the method for crack tip location in aircraft structures, two different fastened joint sample configurations have been manufactured and tested. Given the variety in aircraft structure, the two different samples allow for a broader understanding of the applicability of the method. These samples have been assessed using thermal response data and the QFD method, as well as an optical measurement using a high-resolution camera.

Initially, a double lap configuration is assessed. Secondary bending is minimised using this design which ensures that the fatigue cracks initiate as through-thickness cracks. It was also critical to achieve fatigue crack initiation on one of the visible outer-sheets of the double-lap joint. As such, the middle sheet of the joint was increased to three times the thickness of the outer sheets (outer sheet thickness = 0.063”, middle sheet thickness = 0.190”). Furthermore, one of the outer sheets was chevron-notched at a fastener hole to promote a consistent fatigue crack initiation location and orientation. The sample configuration is shown in Fig. 12(a), and the notch geometry highlighted in Fig. 12(c). The sample sections were waterjet cut and fastened together using eight HL18PB-5-5 Hi-Lok pins and HL70-5 collars. Fastener holes were reamed to produce a nominal 0.0015” interference fit. Three double-lap samples in total were tested and tracked through the course of the crack growth.

Figure 12. Schematic of the double-lap joint test sample.

Next, a riveted joint was assessed, in a 1.5 lap joint configuration. It consists of one load carrying member continuing through the joint, and one member terminating at the joint. This allows for significant fastener load transfer while minimising effects of out of plane bending. L-angle stiffeners are also included on the back face (tail side) of the joint to further prevent out of plane bending. This ensures fatigue cracks are initiated through the thickness of the critical sheet. Sample sections were waterjet cut from 0.063” thick material and a fatigue starter notch was machined into the edge of one of the critical fastener holes. The sections were then fastened together using 12x MS20470AD6-7 protruding head rivets. A schematic of the stiffened 1.5 lap joint sample is shown in Fig. 13. A total of five riveted joints were employed for this study.

Figure 13. Schematic of the 1.5 lap joint test sample.

5.2 Test methodology and instrumentation

Samples were cyclically tested in the same testing machine as the DCT samples at 10Hz under constant amplitude load (R=0.1). The test continued until the crack grew from the initiated fastener hole to the adjacent one. Double lap samples were tested at a maximum load of 30kN, and the 1.5 lap samples at a maximum load of 15kN. In each case, this represented a net section stress in the notched sheet that was approximately one-third of the of the yield stress of the material. Load and machine displacement data was periodically recorded though the machine DAQ every 5k cycles, for a duration of at least 50 cycles.

For TSA data acquisition, the same Xenics Gobi 640 microbolometer and “MiTE” software was used as the DCT samples. For the double lap joints the focal length was adjusted to achieve a spatial resolution of 80µm/pixel. In the case of the 1.5 lap joints, the spatial resolution of the TSA field was 135µm/pixel. The MiTE software was provided with the load signal from the Instron controller through a NI USB-6001 DAQ for correlation. TSA data was acquired periodically in intervals of approximately 5k cycles for an acquisition period of ∼500 cycles.

As an independent form of crack path validation, a digital single-lens reflex (DSLR) camera was utilised for the double lap joints. The DSLR line-of-sight is not perpendicular to the surface of the sample (as illustrated in Fig. 14) and so a projective correction process was applied. This aligned the co-ordinate locations obtained from the DSLR images with the TSA data capture. The correction was achieved through use of a standard checkerboard calibration plate and a MATLAB script to determine the transformation matrix. The crack tip locations are found (optically) from each captured image, and the resulting crack path geometry is then transformed via the matrix, to provide spatially calibrated crack path data.

Figure 14. Experimental setup for the double lap joint samples.

5.3 Results

Samples were tested and data recorded through MiTE software and through the high-resolution camera. Figure 15 shows an example of the quadrature thermal response of the double lap samples, as the fatigue crack grows from the initiated fastener hole to the adjacent one. The primary output from both the TSA maximum gradient crack tip location algorithm and the DSLR projection correction code is the X and Y co-ordinate locations for the crack tip(s). A curve-fitting process was used to smooth the crack position data independently for each of the coordinate sets. Noise was observed in the position measurements due to various reasons, including the limited pixel resolution of the thermal camera, motion of the test sample during testing, and quality of the paint coating. The raw results and the Savitsky-Golay (SG) fit (polynomial order 3, window width of 15 points) are shown for an example case (a stiffened 1.5 lap joint sample) in Fig. 16.

Figure 15. Experimental quadrature response of double lap joint.

Figure 16. Crack path fitting process (rivet data shown): (a) x-coordinate, (b) y-coordinate, and (c) final path data.

Once fitted, each of the double lap sample data sets and the 1.5 lap rivet sample data sets are plotted on spatial axes in Fig. 17(a) and (b). The results of the double lap samples show high level of curvature in the crack paths due to the interference fit of the Hi-Loks, yet all sample paths are tightly grouped. There is also strong agreement between the DSLR and TSA crack paths. In the case of the 1.5 rivet lap joints, the results show four sample crack paths tightly grouped and one outlier growing above the others. The crack paths in general stay closer to the mid plane between the two fasteners. This is likely due to the limited clamping force and hole-fill pressure of the rivets when compared to the Hi-Loks. The lower clamping force of the rivets means a reduced faying surface frictional load transfer, which in-turn leads to a straighter crack path. For the 1.5 lap rivet joints extensive DSLR crack tip location was not employed. At this stage, the maximum gradient method was considered dependable, having been verified with the double lap samples and the DCT samples. However, an image was taken of the final crack path to be compared with the TSA determined path as a form of limited validation. An example is shown in Fig. 18.

Figure 17. Results of the maximum gradient method for (a) double-lap samples and (b) 1.5 lap samples.

Figure 18. TSA crack path on image of final fracture.

For assessment of fatigue crack growth rate, the XY data is converted to total path length using a Pythagorean increment. The data in Fig. 19, shows the total crack length measured from the edge of the fastener hole. An initiation length of 4mm in the double lap and 5mm in the 1.5 lap was used. The corresponding number of cycles were zeroed-off to provide a comparison of the fatigue crack growth period only. Initiation periods for the samples shown in Fig. 19 ranged from 3,000 cycles to 300,000 cycles. In the case of the high-resolution camera, despite capturing a much higher resolution image, the optical method cannot distinguish the crack tip as well as the TSA method. The camera particularly struggles at small crack lengths, and this is evident by Fig. 19(a). However, both the TSA and optical data sets are in strong agreeance for crack lengths of approximately 7mm and longer.

Figure 19. Crack length vs cycles (a vs N) for the (a) Hi-Lok and (b) riveted samples.

As an independent form of crack growth rate validation, one of the double lap samples was disassembled and the fracture surface assessed under a scanning electron microscope (SEM). The striations were counted and averaged to provide a direct measurement of fatigue crack growth rate at certain points during the test. It should be noted that the fracture surface overall was non-uniform and locating regions of consistent striations in the direction of global crack growth (as shown in Fig. 20) was challenging. However, where clear images of striations were found, the measured crack growth rate aligns well with the data obtained by TSA and DSLR.

Figure 20. Quantitative fractography growth rates compared with TSA growth rates.

6.0 Conclusions

Finite element modelling has been used to develop a crack tip location algorithm from TSA experiments. The initial method used features in the quadrature component of the thermal response to locate the crack tip between peaks in the quadrature field. These peaks are due to heat conduction at the high stress gradient locations around the crack tip. The algorithm was tested using DCT samples under cyclic loading. The crack tip location estimates obtained from the TSA algorithm show good correlation with crack lengths determined using the crack mouth compliance method.

The method was further refined for computational efficiency and experimental robustness. The revised method utilised the notion of a stress singularity and divergence of the spatial derivative at the crack tip. Pre-notched double lap samples have been fatigue tested and the resultant crack paths were evaluated using TSA and optical methods. Crack paths were obtained and clearly show a high level of curvature in the double lap samples due to the high interference of the Hi-Lok fasteners. The riveted samples show more linear behaviour in the crack path. Optical measurements were obtained from high resolution photos of the double lap samples and the resultant crack lengths are in good agreeance with the TSA method of crack tip location. SEM images were also acquired at various intervals for one of the samples, and the directly measured crack growth rate aligns very well with the TSA and DLSR data.

The results demonstrated show the high value and validity of the TSA crack tip location method for fatigue crack growth and growth rate evaluation. The method works well even within complex stress fields, such as for fastened aircraft lap joints. The method developed here forms a fundamental basis for more advanced measurement using TSA such as evaluation of the stress intensity factor.