**2. Proposed Finite-Element-Analysis-Based Inspection Interval Optimization Approach**

Figure 1 shows the proposed FEA-based inspection interval optimization approach. The developed methodology consisted of the following five phases:


Phase 3.The uncertainty of mechanical parameters was included, and the design of experiments was created based on Latin hypercube sampling (LHS) scenarios.

Phase 4.A Monte Carlo simulation was performed by building a response surface model.

Phase 5. Risk assessment and inspection interval optimization of the tail wing structure were undertaken.

#### *2.1. Development of the Finite Element Model*

Modern aircraft structures, such as the wings and fuselage, utilize FRCPs. The composite laminates are highly susceptible to delamination, which is considered one of the major damage mechanisms in composites among all other types of damage. The delamination propagation can occur under the fatigue loading, thus causing stiffness and gradual strength degradation and leading to catastrophic failure of composite structures. Therefore, the characterization of the fatigue delamination resistance in composite materials is necessary for their damage tolerance design and reliability assessment. In aircraft composite structures, stringers are joined to the skin via adhesive bonding. During the flight, cyclic fatigue causes separation between the skin and stringer, which may result in delamination failure. Therefore, it is important to predict fatigue delamination initiation and propagation so that these composite panels do not fail prematurely. To enhance the fatigue life and reliability of composite structures, several researchers have investigated the delamination growth in composites under fatigue loading. A variety of models are proposed in the literature for the prediction of fatigue delamination growth at the coupon level. However, the prediction of fatigue delamination in complex geometries, such as aircraft wing structures, is still an open issue. The most popular approach for the characterization of fatigue crack

growth is Paris' law [17–19]. Paris' law relates the fatigue crack growth rate to the stress intensity factor and energy release rate. Equation (1) shows the basic form of Paris' law [20]:

$$dA/dN = \mathbb{C}f(G)^m\tag{1}$$

where *C* and m represent the fitting parameters and *f*(*G*) is a function of the energy release rate (*G*). In metals, the stress intensity factor variation (ΔK) is usually adopted for fatigue crack growth. On the other hand, in composite laminated structures, the energy release rate (Δ*G*) variation seems to provide a better description of experimental results, as shown in Equation (2):

$$\frac{dA}{dN} = \mathbb{C}(\Delta G)^{\mathfrak{m}} = \mathbb{C}(\mathbb{G}\_{\max} - \mathbb{G}\_{\min})\tag{2}$$

**Figure 1.** Schematic of the proposed FE-based inspection interval optimization approach.

The FE model implemented in this study and included in the FE code ABAQUS [21] was built on Equation (2) to estimate the fatigue crack growth rate.

In a numerical model, delamination is usually modeled using a damage mechanics or fracture mechanics approach in combination with Paris' law. Different approaches are found in the literature in the context of damage mechanics. Among these, the cohesive zone modeling approach [22–24] received considerable attention for delamination modeling. This approach provides good results in combination with Paris' law for simple geometries and test specimens. However, this approach is not suitable for complex geometries. On the other hand, the virtual crack closure technique (VCCT) [9,25] derived from fracture mechanics is found built into several commercial FE codes. This approach is suitable for implementation in complex geometries and can simulate fatigue delamination growth in

combination with Paris' law. The fatigue crack growth in an ABAQUS simulation employs VCCT criteria to compute the energy release rate. The crack propagation starts when the criterion based on Equation (3) is satisfied:

$$f = \frac{N}{c\_1 \Delta G^{\circ 2}} \ge G\_{\text{max}} > G\_{th} \tag{3}$$

where *N* is the current cycle number and *c*<sup>1</sup> and *c*<sup>2</sup> are fitting parameters that are experimentally determined. Once the onset criterion is satisfied, the fatigue delamination growth rate is governed by Paris' law, which is shown in Equation (4):

$$\frac{dA}{dN} = c\_3 (\Delta G)^{c\_4} \tag{4}$$

where *c*<sup>3</sup> and *c*<sup>4</sup> are fitting parameters.

The fatigue delamination finite element was first implemented for the double cantilever beam (DCB) composite specimen. The results of the numerical model were compared with the benchmark experimental study [23,26]. The composite specimen had a length of 254 mm and a width of 25 mm. The thickness of each arm was 3 mm. The initial delamination was induced in the specimen with a length of 51 mm and was located in the center of the plies, as shown in Figure 2. The specimen consisted of a T300/1076 carbon fiber epoxy composite with 24 unidirectional plies.

**Figure 2.** DCB specimen geometry with initial delamination.

The material properties used in the numerical model were taken from [26] and are shown in Table 1.



Continuum shell elements (SC8R) were used for the modeling of the DCB specimen. A single element is used throughout the thickness of each arm. The refined mesh was used in the damage propagation zone and the remaining region consisted of a coarse mesh. First, the load–displacement curve was computed for quasi-static loading. The simulation curve was computed and compared with the experimental benchmark case study result, as shown in Figure 3. The numerical model showed a linear response and a good correlation was found, both in terms of the initial stiffness and peak load between the numerical model results and the experimental benchmark case study.

**Figure 3.** Load vs. displacement curve of the DCB specimen [26].

For the fatigue crack propagation analysis, the triangular load cycle as applied in the experimental benchmark case study [23,26] was applied as a fatigue loading cycle. As expected, the delamination length increased rapidly at the start of the loading cycles. The delamination stopped propagating at around 3.7 million cycles. The numerical results computed in ABAQUS showed good agreement with the experimental benchmark results, as shown in Figure 4.

Once the fatigue delamination crack propagation model was established, the model was extended from the composite specimen to the KT-100 aircraft tail wing structure. The tail wing is located behind the main lifting surface of the aircraft and provides stability and control. The design of the KT-100 horizontal stabilizer consists of upper and lower skins. The upper and lower skins of the sandwich structure act as a spur structure. The upper and lower skins are bonded and tightened to make a box-shaped structure. Three ribs exist inside the box structure. Skin debonding from the rib is the most critical damage that can occur in an aircraft tail wing structure. Figure 5 shows skin-debonding damage. Therefore, in this study, skin/rib-debonding damage was modeled using the fatigue progressive damage model.

The geometry of the tail wing was provided by the Korean Airforce and imported from the 3D model in ABAQUS. To simplify the model, bolts and rivets were not considered in the geometry. Initial delamination of 32 mm was introduced between the middle rib and upper skin, as shown in Figure 6a. To avoid the loss of computational power, shell elements (S4R) were used to model the tail wing structure with 49,209 total elements. The mesh was kept refined in the crack propagation region and coarse in the remaining regions, as shown in Figure 6b.

**Figure 4.** Comparison of the numerical and experimental fatigue delamination propagation curve [23,26].

**Figure 5.** Skin-debonding damage in aircraft tail wing structure.

**Figure 6.** (**a**) Geometry of a tail wing with initial delamination and (**b**) the mesh configuration of the model.

The loads and boundary conditions are used as in an actual flight scenario. One side of the wing was kept fixed, and the shear forces were applied at the other end with the flight loading spectrum as a fatigue cycle provided by the Korean Airforce. Figure 7 shows the flight operating loading spectrum with a normalized amplitude.

The bond stat plot (crack propagation plot) output in ABAQUS represents the fatigue crack propagation. The blue color represents the unbonded region, while the red color represents the bonded region. Figure 8 represents the BDSTAT plot for the tail wing structure. The BDSTAT plot was computed between the upper skin and the middle rib of the composite structure.

**Figure 8.** BDSTAT plot representing the cracked region between the upper skin and the middle.

Different case simulations were carried out for the different values of the applied loading scenarios by using the flight operating load spectrum. The simulation was run for 8000 flight hours (FHs), and the crack propagation plots are analyzed. The first case scenario represented the actual loading condition with the magnitude of shear forces. The other three case scenarios were simulated to observe the fatigue crack propagation. Table 2 shows the case simulation scenarios.


**Table 2.** Case simulation scenarios.

The first case simulation consisted of the applied load of 1.849 KN (actual flight loading scenario) for the applied operating load spectrum. The simulation was carried out for 8000 FHs, and it was observed that the crack did not propagate. As the case represents a real loading scenario, it is customary that the crack propagation was slow, and it took more than 8000 FHs to propagate the crack. As the applied load was increased, the crack propagation region increased, as shown in Figure 9. Therefore, it was concluded that to observe the crack propagation, the applied loading must be increased to more than the real flight loading scenario.

**Figure 9.** Crack propagation plot for the different applied loading conditions.

Figure 10a shows the fatigue cracks propagation plot measured as a function of the number of cycles for case simulation 5. The cracks started to propagate around 153 flight times, and increase rapidly until 1436 FHs. From the 14,336 cycles over 3692 FHs, the crack propagation occurred linearly with a low crack propagation rate, while after 4692 cycles, the crack propagation increased rapidly up to 4783 cycles until it reached the defined crack threshold length of 10 mm. The graph trend of the current model results was compared with the crack propagation curve in the literature [27] for a composite. It was found that compared with the composite crack propagation trend in the literature, the current numerical model showed a similar trend of crack propagation.

**Figure 10.** (**a**) Crack propagation plot of the current numerical model and (**b**) a comparison with the literature trend of the composite fatigue crack propagation curve [27].

#### *2.2. Sensitivity Analysis and Design of Experiments (DOE)*

FRCPs subjected to cyclic loading are most likely to fail with fatigue delamination crack growth. Many factors can affect the fatigue delamination growth, such as the constituent material parameters, geometry, and environmental conditions. However, the most common practice is to check the effect of mechanical parameters on the fatigue delamination behavior of composite laminates. For this purpose, sensitivity analysis [28] was undertaken to assess the sensitivity of crack growth to uncertain mechanical parameters. Sensitivity analysis can be used as a tool to understand how uncertain inputs can affect a model's performance. Sensitivity analysis is a systematic assessment method that is generally performed to assess the impact of individual input uncertain parameters on the model output response. It is an essential part of every risk assessment analysis that seeks to learn things such as how the model outputs change with the change in inputs and how it affects the model output decisions. The sensitivity analysis enhances the overall confidence in the risk assessment. Further, it improves the prediction of the model by studying the model response to the change in input variables and by analyzing the interaction between the variables. The knowledge of these sensitive parameters can help to better comprehend the fatigue delamination behavior and can pinpoint the direction of an optimal composite design. In this study, the effect of varying the mechanical properties was studied to examine the fatigue delamination behavior. Table 1 shows that a total of eight mechanical properties and two fatigue parameters were considered to examine the effect of fatigue delamination behavior. The parameters were varied by +2% of the mean value, and the model output fatigue delamination growth was computed. It was observed that E1 and G23 showed the most sensitive behavior compared with the other parameters. Figure 11 shows the relative error plot.

After pinpointing the most sensitive mechanical inputs, the design of experiments (DOE) was carried out by using Latin hypercube sampling (LHS) [29]. For the generation of DOE scenarios, input factors should be defined. In the case of fatigue crack propagation in the tail wing structure, two input factors (E1 and G23) were chosen. The error variation of +2% was kept for random sampling for the generation of thirty DOE samples, as shown in Table A1 of Appendix A. Figure 12 shows the fatigue crack propagation curves combined into a single plot. It was observed that for a crack delamination length threshold of 10 mm, the number of flight times range from approximately 2500 to 4600.

**Figure 11.** Sensitivity analysis plot representing the relative error of each parameter with respect to the exact simulation parameters fatigue behavior.

**Figure 12.** Fatigue delamination behavior of DOE data generation scenarios.
