Automation of Aircraft Fatigue Life Estimation

Fatigue life is an extremely important factor during the exploitation of any machine. Its estimation requires an in-depth knowledge on the structure of the object, as well as conditions it will be exposed to. Although many industries, such as automotive (Reza Kashyzadeh, 2020), agricultural (Mattetti et al., 2012), or aerospace (Tavares & de Castro, 2017), use different standards for estimation and measurement of overall fatigue strength of a designed structure, the fundamental methods stay constant. The main idea is to predict the most probable time or number of cycles the structure can endure (Zhang et al., 2023). The prediction is mostly performed using data on statistical fatigue strength of a given material’s test sample. Depending on type of the analyzed material, expected scale and directions of loads, and other factors, multiple models of fatigue accumulation may be implemented.

In aviation, throughout the years, several methods of load prediction for the expected mission profiles were developed. In order to accelerate the initial design phase, engineers may use standardized load spectra. For instance, an analytically acquired block-spectra for sailplanes was published by Thielman, Franzmeyer in 1969 and was improved by work of Kossira, Reinke with experiment-based spectra in 1982 (Kossira & Reinke, 1989). However, as advised by works of Rodzewicz (Rodzewicz, 2007) and Kensche (Kensche, 2019), substantially better results may be achieved using dedicated spectra, based on flights of similar aircrafts with identical or comparable mission designs.

Given recent improvements in unmanned aerial vehicle (UAV) technology, as well as ongoing miniaturization of onboard inertial measuring units (IMU) and microcontroller units (MCU), it is possible to gather live data during the flights of aircrafts even lighter than 150 -g, such as Cinewhoops, which is a well-known device in the first-person view (FPV) drone-use community. Some of the popular operating systems for UAVs offer automated onboard sensor indications’ recording (Logs — Copter documentation; Blackbox logging internals | Betaflight), which reduces the burden of collecting in-flight load spectra. This brings a new possible approach in preparation of load estimations for aircraft-specific mission profiles. Devices operating on such software could be quickly deployed in both unmanned and manned aircrafts as a flight controller or an additional low-weight payload, respectively.

After gathering the data, a filtration must be applied to properly represent loads that could be reenacted on the analyzed structure by means of an experiment (Rodzewicz & Przekop, 2000) or a simulation (Maksimovic et al., 2021). In aircrafts, while some events and resulting accelerations may be treated as random (mainly gusts and landings), flight loads can be assumed as mostly deterministic (“Load Spectra”, 2009). This allows for use of time-based methods, instead of frequency-based ones, which are better suited for stochastic events. Historically in aviation the filtration was performed on the go, mostly due to memory limitations. The load history was stored as an array containing the number of discrete load level (LL) changes. Depending on the method used for detecting a cycle, as well as the format used for data storage, the array was called a halfcycle array or a full-cycle array. One of the most popular time-based filtration methods was the Rainflow Cycle Counting (RCC), as presented in work of Matsuichi, Endo (Matsuichi & Endo, 1968). It follows the locally extreme load values and counts the number of LL changes. The length of the half-cycle depends on previously detected cycles or the minimal LL encountered in a cycle. Its main advantage is robustness against irregularly increasing (or decreasing) LL, in which case classical load level difference counting would ignore low frequency high amplitude load changes.

Lastly, the influence of cycles endured is estimated using an application-specific material fatigue model. For metals and many other materials, estimation accuracy is satisfied by the Palmgren-Miner rule (Murakami, 2019), although other methods may be implemented for non-isotropic materials, such as textile composites (Lomov & Xu, 2015) or 3D prints. Combined with the half-cycle array, it enables a relatively low-cost computation-wise estimation of overall fatigue damage.

The main goal was to develop a software aiding fast and reliable preparation of a custom load spectra for a structure integrated into an aircraft. The processed data were to come from a real-time flight telemetry, collected by either ground station or saved in onboard memory. The data should be filtered and converted into a half-cycle array. The results of the software should consist of calculated load spectra as presented by Owczarek, Rodzewicz (Owczarek & Rodzewicz, 2009), as well as predicted fatigue lifetime of the structure.

The remainder of this paper is organized as follows: the forthcoming section contains assumptions made during the design phase of the software. Section 3 describes methods of problem analysis. Section 4 provides individual software module tests. The sample analysis results are discussed in Section 5. The paper is summarized with conclusions and further development possibilities in Section 6.

As stated in Section 1, the main goal of the project was to develop a software capable of aiding aircraft engineers in fatigue life calculations of designed structures. During the planning phase of the project, the Python 3 (Drake & Rossum, 2009) programing language was chosen, due to its multiplatform capabilities. The main assumptions of the project were as follows:

the available data come from a three-axis accelerometer set in the aircraft’s center of gravity, or is pre-processed to eliminate linear components of angular acceleration;

The forthcoming section describes the means of achieving the goals described earlier.

During the program execution, three main problems are resolved:

filtration of the available data,

The methods of solving the above are presented in the forthcoming subsections.

3.1.

Filtration of the available data

In terms of fatigue damage level calculation, it is useful to convert the direct log of accelerations into a Markov array (Kossira & Reinke, 1989). Historically in aviation it is a two-dimensional array with a size of 32 columns by 32 rows. In order to create it, the time-based data must be filtered by extreme values and transformed from accelerations to LLs (Rodzewicz, 2007). An example of the available data array is shown in Table 1.

Table 1.

Example of telemetry log used in calculations

Date	Time	Hdg(@)	Ptch	Roll	AccX(g)	RSSI(dB)	AccY(g)	AccZ(g)
2021-10-10	16:27:29.650	178.60	65.6	-17.0	-0.06	82	-0.02	0.05
2021-10-10	16:27:29.850	178.60	65.6	-17.0	-0.06	87	-0.02	0.05
2021-10-10	16:27:30.050	186.60	49.7	0.5	-0.09	70	-0.02	0.05
2021-10-10	16:27:30.250	186.60	49.7	0.5	-0.09	70	-0.09	0.14
2021-10-10	16:27:30.450	202.10	49.7	0.5	-0.09	80	-0.09	0.14
2021-10-10	16:27:30.650	202.10	28.9	8.5	-0.07	85	-0.11	0.21
2021-10-10	16:27:30.850	202.10	28.9	8.5	-0.07	85	-0.11	0.21
2021-10-10	16:27:31.050	233.70	23.1	28.5	-0.02	90	-0.08	0.19
2021-10-10	16:27:31.250	233.70	23.1	28.5	-0.02	90	-0.08	0.19
2021-10-10	16:27:31.450	273.10	40.2	28.5	-0.02	87	-0.08	0.19
2021-10-10	16:27:31.650	273.10	40.2	44.1	0.00	87	-0.03	0.17
2021-10-10	16:27:31.850	273.10	40.2	44.1	0.00	80	-0.03	0.17
2021-10-10	16:27:32.050	295.70	54.7	41.5	0.00	87	-0.01	0.09
2021-10-10	16:27:32.250	295.70	54.7	41.5	0.00	87	-0.01	0.09
2021-10-10	16:27:32.450	328.70	59.7	34.9	0.00	90	-0.01	0.09

The filtration begins with normalization of the logged accelerations. According to Assumption 9 described in Section 2, the available data should be appended with a calibration set of accelerations (aτxi)calib${(a_\tau ^{{x_i}})_{calib}}$, with x_i. meaning one of three main aircraft axis. The calibration should be performed by positioning the sensor in three perpendicular orientations, and so sensor readings for a = g acceleration in X, Y, and Z axes are acquired. Then, the normalization is performed on the data using the equation, 3.1 ntxi=atxi∑τ=0k(aτxi)calibk $$n_t^{{x_i}} = {{a_t^{{x_i}}} \over {{{\mathop \sum \limits_{\tau = 0}^k {{(a_\tau ^{{x_i}})}_{calib}}} \over k}}}$$ where t denotes the sample index and axi the acceleration in a chosen axis. It is worth noting that, despite the unusual notation of acceleration in Table 1, the load factors resulting from Eq. (3.1) are represented correctly in non-dimensional form. This means that, provided the calibration recording was performed the same way as the flight itself, the method will return correct values. A comparison between raw and normalized data is shown in Figures 3.1a and 3.1b.

Figure 3.1.

The next step of the filtration is computation of LLs based on load factors calculated earlier. According to Assumption 7, the critical load factors of the aircraft are known. With such information it is possible to calculate a (n) function. Kossira and Reinke (Kossira & Reinke, 1989) define LLs in a function of load factors as a discrete line between points (3, (n^x_i)_maxcrit) and (31, (n^x_i)_mincrit) in (LL, n) space. This means that the desired function is a line with factors of: 3.2 [aibi]=[3−31(nxi)maxcrit−(nxi)mincrit31−3−31(nxi)maxcrit−(nxi)mincrit⋅(nxi)mincrit] $$[\matrix{ {{a_i}} \cr {{b_i}} \cr } ] = [\matrix{ {{{3 - 31} \over {{{({n^{{x_i}}})}_{maxcrit}} - {{({n^{{x_i}}})}_{mincrit}}}}} \hfill \cr {31 - {{3 - 31} \over {{{({n^{{x_i}}})}_{maxcrit}} - {{({n^{{x_i}}})}_{mincrit}}}} \cdot {{({n^{{x_i}}})}_{mincrit}}} \hfill \cr } ]$$

Data from Figure 3.1b transformed by Eq. (3.2) for (n^x_i)_maxcrit = 4.796 and (n^x_i)_mincrit = −2.398 are shown in Figure 3.2.

Figure 3.2.

Knowing the LLs during the flight, it is possible to filter out non-extreme values from the log. The filtration may be performed by analysis of LL gradient sign. When the derivative of (n) function by data point index changes sign and is non-zero, the algorithm saves the argument of found extremum. It is highly important to not consider points with the derivative value of 0, as it introduces noise into the filtered data. This step of analysis enables matrix size reduction and, by extension, lowers computational costs of the subsequent steps. With such an approach it may be easier to implement the project into a super low-power device capable of long-term low-maintenance data collection onboard an aircraft.

The block flow diagram of extreme values filter (EVF) is shown in Figure 3.3 and an example of filtered data is shown in Figure 3.4.

Figure 3.3.

Figure 3.4.

With the filtered set of extreme values, it is possible to calculate the half-cycle array.

3.2.

Half-cycle array computation

Based on the advice provided in the work of Owczarek, Rodzewicz (Owczarek & Rodzewicz, 2009), the RCC (Tavares & de Castro, 2017) method was used during the preparation of a half-cycle array. The main idea behind RCC is to count “streams” flowing down the LL extremes time-wise. Each stream begins in a consecutive extreme LL (source) and is broken if any of the following conditions is met:

a source of a deeper extremum exists (LL₀ > LL_t for even streams [upper] and LL₀ > LL_n for odd streams [lower]),

or a stream of higher extreme LL already flows through the current stream position (the end of the stream is set to the LL of an older stream that has been met).

With the above set of rules it is possible to prepare a recursive algorithm for the RCC method. A block flow diagram of the algorithm and its results are shown in Figures 3.5 and 3.6, respectively.

Figure 3.5.

Figure 3.6.

The computed streams can then be saved as vectors of source and end LLs: 3.3 s→=[ LL(0),LL(n−1)] $$\vec s = [LL(0),LL(n - 1)]$$ where n is the complete number of steps flowed by a stream. In this form it is less computationally demanding to count the number of specific streams, as both values are used as column and row index of half-cycle array, respectively. An example of an active zone of the calculated half-cycle array is shown in Table 3.2.

Table 3.2.

Example of active zone in half-cycle array.

LL	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26
25	1	0	2	0	2	1	0	1	1	0	0	1	0	0	1	1	1	x	0
24	0	1	2	2	0	1	1	1	1	0	0	0	0	0	0	2	x	1	0
23	0	1	2	2	1	3	0	4	2	1	0	1	0	0	3	x	2	1	0
22	0	1	0	0	1	3	5	7	2	1	2	1	1	2	x	3	0	1	0
21	0	0	0	0	2	0	3	8	10	8	13	10	5	x	2	0	0	0	0
20	0	0	0	0	0	1	0	3	5	10	30	53	x	5	1	0	0	0	0
19	0	0	0	0	0	0	1	2	0	7	29	x	52	10	1	1	0	1	0
18	0	0	0	0	0	0	1	3	3	8	x	29	31	12	3	0	0	0	0
17	0	0	0	0	0	0	0	0	1	x	8	7	10	8	1	1	0	0	0
16	0	0	0	0	0	0	0	0	x	1	3	0	5	11	2	1	1	1	0
15	0	0	0	0	0	0	1	x	0	0	3	2	3	8	7	4	1	1	0
14	0	0	0	0	0	0	x	1	0	0	1	1	0	3	5	0	1	0	0
13	0	0	0	0	0	x	0	0	0	0	0	0	1	0	2	3	2	1	0
12	0	0	0	1	x	0	0	0	0	0	0	0	0	2	1	1	0	2	0
11	0	0	0	x	1	0	0	0	0	0	0	0	0	0	0	2	2	0	0
10	0	0	x	0	0	0	0	0	0	0	0	0	0	0	0	3	0	2	2
9	0	x	0	0	0	0	0	0	0	0	0	0	0	0	1	1	1	0	1
8	x	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1
7	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0

The resulting half-cycle array [T] can be used in fatigue life calculation, which is described in the forthcoming subsection.

3.3.

Fatigue lifetime calculation

In order to acquire an estimated time of the designed structure’s fatigue life, a specific deconstruction of the half-cycle array must be performed. First, the LLs must be converted back into load factors using the inverse linear function according to Eq. (3.2): 3.4 [ab]=[(nxi)maxcrit−(nxi)mincrit3−31(nxi)mincrit−(nxi)maxcrit−(nxi)mincrit3−31⋅31] \[[\begin{matrix} a\\ b\\ \end{matrix}]=[\begin{array}{*{35}{l}} \frac{{{({{n}^{{{x}_{i}}}})}_{maxcrit}}-{{({{n}^{{{x}_{i}}}})}_{mincrit}}}{3-31}\\ {{({{n}^{{{x}_{i}}}})}_{mincrit}}-\frac{{{({{n}^{{{x}_{i}}}})}_{maxcrit}}-{{({{n}^{{{x}_{i}}}})}_{mincrit}}}{3-31}\cdot 31\\ \end{array}]\]

Next, a load factor array is prepared by computing for each non-zero half-cycle array cell: 3.5 ni,ja=|nz(LLv)−nz(LLh)2| $$n_{i,j}^a = |{{{n_z}(L{L_v}) - {n_z}(L{L_h})} \over 2}|$$ 3.6 ni,jm=|nz(LLv)+nz(LLh)2| $$n_{i,j}^m = |{{{n_z}(L{L_v}) + {n_z}(L{L_h})} \over 2}|$$ where [n^a] represents the load factor amplitude matrix, [n^m] the mean load factor matrix, LL_h the column index of half-cycle array, and LL_v the row index of half-cycle array. However, it is important to emphasize that the exact values of the prepared arrays are independent from the half-cycle array computed in the preceding subsection. Both [n^a] and [n^m] matrices represent the corresponding load factors based only on column and row index. LL_h and LL_v are just sets of numbers between 1 and 32. The proper use of the half-cycle array values will be shown directly during fatigue damage level matrix computation.

According to Assumption 6 from Section 2, the mean stress value for the n = 1 situation is known. The computed matrices may then be scaled by a factor of known σ(n = 1), which results in the following arrays of mean stress matrix [σ^m] and stress amplitude matrix [σ^a]: 3.7 [σm]=[nm]⋅σ(n=1) $$[{\sigma ^m}] = [{n^m}] \cdot \sigma (n = 1)$$ 3.8 [σa]=[na]⋅σ(n=1) $$[{\sigma ^a}] = [{n^a}] \cdot \sigma (n = 1)$$

Having calculated stress matrices, it is now possible to compute the maximum number of cycles for each non-zero cell of the half-cycle array [T]. Due to σ(n = 1) being a result of a continuous function, it may not always fit the discrete material data available within the input data (Assumption 8).

The software compensates for lack of information using a linear approximation between the lines of material constant mean stress (as shown in Figure 3.7). The maximal number of cycles is calculated with: 3.9 Ni,j=10a(σi,jm)⋅σi,ja+b(σi,jm) $${N_{i,j}} = {10^{a(\sigma _{i,j}^m) \cdot \sigma _{i,j}^a + b(\sigma _{i,j}^m)}}$$ where both a and b are a function of σi,jm$\sigma _{i,j}^m$.

Figure 3.7.

Having calculated both the half-cycle array and the maximum number of cycles matrix, it is now possible to calculate the fatigue level matrix: 3.10 Di,j=Ti,jNi,j \[{{D}_{i,j}}=\frac{{{T}_{i,j}}}{{{N}_{i,j}}}\]

The overall fatigue level may then be calculated using the Palmgren-Miner rule as a sum of all elements of matrix D: 3.11 Dc=0.5∑i=132∑j=132Di,j=0.5∑i=132∑j=132Ti,jNi,j \[{{D}_{c}}=0.5\underset{i=1}{\overset{32}{\mathop \sum }}\,\underset{j=1}{\overset{32}{\mathop \sum }}\,{{D}_{i,j}}=0.5\underset{i=1}{\overset{32}{\mathop \sum }}\,\underset{j=1}{\overset{32}{\mathop \sum }}\,\frac{{{T}_{i,j}}}{{{N}_{i,j}}}\] and the hourly fatigue lifetime: 3.12 TFL=tDc \[{{T}_{FL}}=\frac{t}{{{D}_{c}}}\] where t indicates time period of the logged data for which the half-cycle array [T] was calculated.

The forthcoming section describes methods of testing the described algorithms.

It is useful to prove the validity of a new software’s submodules before performing complete testing. The developed tool consists of three components – the aforementioned EVF, stream generator (SG), and array generator (AG). Each of the submodules performs calculations according to the rules and equations specified in Section 3.

4.1.

Validity of Extreme Values Filter

The most important factor of data preparation is correctly filtering extreme values for use in half-cycle array. The EVF submodule should return only local extremes of LLs without considering regions of zero LL difference. Due to its implementation, the module is expected to find only the last extreme point before a change of derivative sign. Values returned are expected to be integers between 1 and 32.

The validation was performed on three sets of logged LLs by comparing the input with the results of the filtration. All three sections have shown correct behavior of computed polylines. Test 3 seems to be of significant importance as it contains multiple steps of little LL difference; however, the local extremes are found correctly. The tests are visualized in Figures 4.1–4.3.

Figure 4.1.

Figure 4.2.

Figure 4.3.

4.2.

Validity of Stream Generator

In order to validate the SG submodule, generated streams were imposed on points of extreme LLs from Subsection 4.1. The expected behavior of streamlines was specified in Subsection 3.2. The tests were visualized in a manner similar to that adopted in Subsection 4.1, as portrayed in Figures 4.4–4.6. They show the streams following the specified rules, with the exception of Test 3, where streams with sources at Points 180 and 194 are broken seemingly without a cause. This is, however, a limitation of the visualization method, as the global chart shows the streams are blocked by at least one older stream (as shown in Figure 4.7).

Figure 4.4.

Figure 4.5.

Figure 4.6.

Figure 4.7.

Knowing that each of the main submodules works correctly, it is possible to test the complete software by analyzing real-life data. Due to a lack of third party data for result comparison, the output was validated against the rules specified in previous sections.

During the tests, an aluminum sample with a static load of n_z = 1 characterized by a stress of σ(n_z = 1) = 37.298 MPa was analyzed. The model represents a pin of camera mounting in a typical sport quadcopter; however, the method should be applicable to any aluminum component of the structure. The considered critical load factors are (n^x_i)_mincrit = −2.398 and (n^x_i)_maxcrit = 4.796. The material characteristics are shown in Table 5.1. The analyzed mission profile was designed according to the Polish NSTS-06 rules of flight.

A log of computed LLs is visualized in Figure 5.1. The log does not comply with the assumption of continuous acceleration sampling; however, the available data are labelled with exact time stamps of measurements.

Figure 5.1.

The visualization of the generated RCC streams (Figures 5.2 and 5.3) was divided into odd and even streams for improved readability. The results show coherence with earlier tests, with acquired streams following the rules specified in Subsection 3.2. In order to further improve the readability of the following diagrams, index of sample was used instead of time stamp.

Figure 5.2.

Figure 5.2.

Based on the above results, the software calculated half-cycle array [T]. The active zone of the matrix is shown in Table 5.2. The results are on par with the expectations specified in Subsection 4.3. The highest values are reached near cells (18, 19) and (19, 18) corresponding to the load factor of n_z ≈ 1. The diagonal consists of only 0 values, which shows that no zero LL changes were found in the log. This proves further the correctness of the EVF submodule algorithm. The matrix has the same vertical and horizontal dimensions, which shows that streams were most probably generated correctly.

The active zone of the fatigue level matrix [D] is shown in Table 5.3. Evaluating it against the rules and expectations specified in Subsection 4.3, the array seems to be generated correctly. However, it is currently not validated against third party data, as such were not available during the software development. It is essential to perform further validation of the project in the near future.

The overall fatigue level of the analyzed object after 20 min of loads shown in Figure 5.1 is estimated to be D_c = 0.00154. This suggests the object fatigue lifetime to last about T_FL ≈ 220 h. The result may seem to be low; however, it is important to notice the tests’ high loads relative to a typical aircraft. This may generate high stresses in the structure, possibly introducing conditions for low cycle fatigue.

The main goal of the project was to develop a tool aiding engineers in fatigue lifetime calculations, based on real-life data collected during exploitation of an aircraft. The software was supposed to perform automatic transformation of telemetry logs to custom load spectra and estimated fatigue level, which was then to be recalculated to hourly fatigue lifetime. The input data were supposed to consist of flight log formatted according to the OpenTX telemetry standard and an array of discrete material fatigue data based on experiments.

The paper describes methods, algorithms, and validation of implemented software modules used in calculations, as well as the results of test data analysis. The test data consisted of accelerations imposed on an aluminum sample simulating a pin of a camera mounting on a typical sport UAV during flight according to the Polish NSTS-06 rules. The results show high coherence with all of expectations stated during the algorithm design phase, which suggests that the estimations have been computed correctly.

The results of the analysis show rather a low fatigue lifetime of the sample. However, it is important to emphasize that during the test flights, high-G acrobatics were abundantly performed, which significantly influences the overall load spectra of the aircraft. Moreover, the tests represent a relatively short sample size. In order to prepare full custom load spectra, the test flights should represent real-life situations as much as possible. To adequately estimate fatigue lifetime of the sample, further tests must be performed.

A strong advantage of this software is its multiplatform approach, with a focus on the requirements for the computational power and memory needed during calculations to remain low. This means the developed software could be used on mobile platforms for live data gathering and analysis. It would require integration of the pre-processing submodule, eliminating the effects of being placed outside of the center of gravity. This approach could provide a highly augmented view on fatigue levels in both small UAVs and manned aircrafts.

Further project development should also focus on extrapolation of half-cycle arrays. This approach could provide a view on a higher range of loads, sometimes very hard, problematic, or dangerous to invoke in real-life testing.