A new methodology to predict damage tolerance based on compliance via global-local analysis

A BSTRACT . Over the years several design philosophies to address fatigue have been developed trying to combine structural safety and economy with aircraft manufacturing and operating processes. The safe-life approach consists of designing and manufacturing an aeronautical structure to be safe throughout its useful life. This approach results in factors that oversize structural elements to prevent possible failure and evidently incurs high design costs. Alternatively, a damage tolerance concept-based approach assumes that the structure, even when damaged, is able to withstand the actions for which it was designed until the detection of a crack due to fatigue or other defects during its operation. This research proposes a new methodology to address the damage tolerance problem in which two-dimensional global-local analysis at different levels of external requests will be made by means of compliance, aimed at finding a relationship between fatigue life and the Paris’ constant. Moreover, the study uses the BemCracker2D program for simulating two-dimensional crack growth. This methodology has proved to be an efficient alternative applicable to damage tolerance analysis.


INTRODUCTION
he National Transportation Safety Board [1] and Wanhill [2], have highlighted a number of documented studies to interpret the cause of airplane accidents due to fatigue processes such as those with the Comet and the Boeing 737-200. Sanford [3] states that all such structures are subject to fatigue, with cracks starting at the edges and growing until reaching critical size for brittle failure to occur, even when subjected to permanent loads only. Unlike T structures that are plastically overloaded, where large deformations occur before failure, those caused by fatigue occur suddenly, without warning (brittle failure). This proves to be a big problem as it makes it impossible to take preventive measures before the complete rupture, so information related to the variable load over time and, particularly, its effect on cracks is of fundamental importance for predicting the behavior of the structure as a whole. Over the years, several fatigue design philosophies have evolved trying to combine structural safety and economy in aircraft manufacturing and operating processes. The first approach, called safe-life, consists of designing and manufacturing an aircraft structure that is safe for its entire lifetime. To that end, the most extreme situations of foreseeable fatigue requirements arising during operation must be considered in the prototype tests. This methodology results in factors that oversize the structural elements in order to prevent the possibility of failure. It is an approach that evidently leads to high design costs and is not capable of guaranteeing security as to whether an unforeseen failure may occur during project life. Rationally, a new philosophy has been developed based on the concept of damage tolerance. In this methodology, it is assumed that the structure, even when damaged, is able to withstand the actions for which it was designed until the detection of a fatigue crack or other defects during its operation. The unit is then checked, repaired and put back into service until the end of its useful life. Palmberg [4] has been a pioneer in this new concept, performing a statistical analysis to control fatigue crack propagation and considering inspection intervals to keep the probability of complete failure low. Later, Wanhill [5,6] examined damage tolerance in the use of aluminum alloys for aircraft structural applications. Newman Jr. [7,8] has suggested that fatigue damage can be characterized by crack size. Schijve [9] proposed some aspects of the design, predictions and experiments associated with the same concept in aircraft structures. Barter and Molent [10] showed that the load cycles have a direct linear relationship with the logarithm of the crack size and that the largest cracks formed grow approximately exponentially (the so-called "main crack" methodology) [11] from small discontinuities inherent to the material, as soon as an aircraft enters into service [12]. Currently, the concept of damage tolerance is applied to aircraft with composite structures [13][14][15][16] in multiple crack analysis [17] and in shape optimization projects [18,19]. Probabilistic studies on damage tolerance are based on fabrication components [20] and fatigue life dispersion from an initial defect distribution [21]. Other works relate damage tolerance to computational methods, using Finite Element Method [22] and Extended Finite Element Method (XFEM) [23], Boundary Element Method (BEM) [24] and Dual Boundary Element Method (DBEM) [25]. In an aircraft fuselage, the evaluation of Fracture Mechanics (FM) parameters such as Stress Intensity Factor (SIF), number of load cycles and stress and displacement fields becomes difficult due to the complex nature of the panel details: brackets, shear clips, rivets, etc. On the other hand, knowledge of those parameters is of paramount importance for understanding the nature of the damage process, especially under the action of dynamic loads. Thus, designers are always looking for fast and reliable simulation methods that can produce accurate average data from these parameters to avoid damage processes and, consequently, the occurrence of accidents. Automation is seen as a key point, allowing for the evaluation of various analyzes as a means for conducting parametric studies and resulting in design optimization [26]. Among the several numerical methods for fracture modeling and analysis, DBEM has been consolidated and has the following advantages: simplified modeling of the crack area, direct SIF calculation, reduced execution times and accurate simulation of crack growth [27][28][29]. The behavior of a solid, discretizing only its contours, enables the analysis of the thousands of simulations necessary for a probabilistic study and, using DBEM, it is possible to study the defects, predicting the fatigue behavior, especially the damage process, multiple sites damage and reliability analysis, among others [30][31][32][33][34]. This work aims to find a solution for the damage tolerance problem detected in [35]. For that, two-dimensional globallocal analyzes will be carried out under different levels of external demands through compliance. The method aims to find a relationship between fatigue life and Paris' constants. Furthermore, a homemade software called BemCracker2D and its GUI BEMLAB2D version will be used for modeling and analysis of two-dimensional elastostatic problems involving cracks, which are based on the BEM and DBEM.

LITERATURE REVIEW
opez [36] presents an extensive review of the uncertainties involved in the monitoring of structural damage in aircraft addressing the existing methods developed for the problem of uncertainty in the areas of damage diagnosis, prognosis and control. Newman Jr. [8] predicts the fatigue life of various metallic materials under different loading conditions. This study made it possible to express crack growth as a function of the effective Stress Intensity Factor interval. The results obtained were compared with experiments on notched and unnotched specimens of aluminum and steel alloys. Fatigue is commonly represented by the study of [37]. Those authors made a fundamental L contribution to the modeling of cracks in solids subjected to cyclic loads from the Paris-Erdogan Law, also called da/dN, which relates the number of fatigue cycles with the crack size. Palmberg [4] is one of the pioneers when it comes to considering the concept of damage tolerance. That author performs a statistical analysis to control the propagation of fatigue cracks and considers inspection intervals in order to ensure that the probability of complete failure of the structure is always kept low. Pyo [38] deals with an alternative method for analyzing fatigue damage in aircraft structures. The author developed a methodology called Elastic Finite Element Alternating Method (EFEAM) to predict the maximum load capacity in cracked panels highlighting the Multi Site Damage (MSD) effect. Also in the same study, an analytical solution for a crack line in an infinite metal sheet is developed and as a result the author demonstrates that the classical LEFM approximation overestimates the load capacity. Jeong [39] presents a method for predicting the MSD threshold and also widespread damage to fuselages saying that the problem of widespread damage is the reduction in the residual strength of the structure below the tolerance while the MSD threshold refers to the point in useful life when there is fatigue coalescence (linkup) of two adjacent cracks still at the allowable stress level. The presented methodology determines a threshold value from the analysis of the combination of residual strength and fatigue crack growth results evaluated through laboratory tests. For that model, the results for the fatigue damage threshold were between 32,000 and 40,000 cycles and for the MSD threshold, about 70,000 cycles. Platz [40] points out that fatigue cracks in lightweight shells or panel structures can lead to major failures when used for sealing or cargo transport. In his research, that author investigates the application of piezoelectric systems applied to the surface of a thin cracked aluminum panel to reduce fatigue crack propagation. With the reduction of propagation, uncertainties in the strength of the structure, which remain even when the structure is used under damage tolerance conditions such as in an aircraft fuselage, can be reduced. Piezoelectrics act by inducing compressive forces at the crack tip to reduce cyclic SIF. As a result, it has been statistically highlighted using experimental samples that the crack propagation rate significantly reduces. Khan [41] has examined low cycle fatigue in aluminum Al 2024-T351 plate. Experimental analysis was performed for both monotonic and cyclic loading, using imaging technology to detect the crack start site. Breitbarth [42], based on biaxial tests with samples arranged in a cross, studied cracks in fuselage sections between the wing and tail of the aircraft, obtaining maximum SIF's values for metallic parts. The experimental results from digital images were compared with the analytical ones, obtaining studies of SIF's, J-Integral, plastic zone, and crack closure effects. Finally, the publications of [9] and [43] revealed that a fuselage subjected to compression testing to the Ultimate Limit State (ULS) and then fatigue tested had improved results due to residual compressive stresses. Those tensions were able to delay the development of cracks.

BASIC FORMULATION
his section presents the basic formulation of the main topics related to the methodology, namely, compliance and fatigue life.

Compliance
Compliance is the magnitude that represents the inverse of rigidity. Mathematically it is determined by: where is the displacement from an applied load P. In Fig. 1,

Fatigue Life
Fatigue life can be represented by the Paris' Law. This law relates the crack propagation rate (da/dN) to the variation of the Stress Intensity Factor (ΔK): where C and m are material constants, determined experimentally, a is the crack size, N is the number of loading cycles and K is the variation of the Stress Intensity Factor. Developing Eqn. (2), the classic condition of the number of cycles, which varies from an initial crack 0 a that can be detected by conventional methods, to a critical length c a , as determined by: where N is the number of cycles needed to increase the initial size crack and tot N is the expected number of cycles over the life. This variation between the number of cycles of the initial crack 0 N and the critical c N , when a safety factor (SF) is applied, is defined as the inspection interval to be adopted, mathematically represented by:

MATERIAL AND METHODS
his section presents the computational technique developed for optimizing the fatigue life of aircraft fuselage parts through compliance. The technique aims to find the relationship between the physical parameters of the Paris' Law material (C, m), such that the fatigue life equals that defined in the project. Moreover, a brief description will be made of the computational tools used in the global-local modeling and analysis process (BEMLAB2D and BemCracker2D softwares) in order to corroborate the presented methodology. This technique is based on the four steps of the following algorithm: 1) From the model in BEMLAB2D GUI by the user, the number of cycles required, or project, is defined (n*); 2) The algorithm computes the stress field in the macro analysis and locates the stress peak before reaching plastification, thus enabling the elastic analysis; 3) The algorithm positions the micro element at the stress peak found in step 2, and calculates the number of cycles for which compliance reaches 3C (N 3C ); 4) The algorithm, considering the initial fuselage defects (hole size and crack lengths), obtains a series of physical parameters of the material to be used (C and m) taking the number of cycles of the minimum instability compliance (N3C) to be that defined by the user (n*). For example, if the designer wants to obtain the instability at n*=6x10 12 , the optimization will show which series of physical parameters (C and m) the material needs to have so that it takes N 3C to the number of cycles defined in the design (n*), as illustrated in Fig. 2.

BemLab2D GUI
The BemLab2D software [44] is a graphic interface for manipulating two-dimensional models of boundary elements, allowing geometric information, boundary conditions and physical attributes to be managed in an efficient and userfriendly way. BemLab2D works both as a pre-processor when defining the geometric model of the problem, by associating physical attributes to the geometry and by generating the boundary elements mesh, and as a post-processor T when used to visualize the data provided by BemCracker2D, as deformed graphs, SIF's, number of cycles and crack propagation, among others. Fig. 3 illustrates these features. This graphical interface, for pre-and post-processing, is able to represent the macro and micro models by boundary elements, as detailed in the following sections.

BemCracker2D Program
BemCracker2D [45][46][47][48] is the processor program for elastostatic analysis of 2D problems. It is written in C++ language and entirely structured on the concepts of Object Oriented Programming (OOP) in order to perform analyzes using the Boundary Elements Method. Here, the module III (DBEM with propagation) of BemCracker2D is used, which consists in: Stress Analysis with standard BEM -SIF's Evaluation (J-Integral) -Direction/Correction of Crack Growth -Fatigue Life Evaluation (Paris' Law) and Coalescence of multiple cracks (linkup). This software was developed from the standard BEM modeling and the incremental analysis strategy for problems with cracks [49,50]. For analysis of Fracture Mechanics problems, BemCracker2D calculates elastic stresses using the conventional BEM and performs incremental analysis of the crack extension using the DBEM. SIFs are computed for each increment through the J-integral, the propagation direction by the three classical criteria in the literature, although here the maximum circumferential stress is used, and the crack growth rate by one modified Paris' equation, which uses an equivalent SIF considering fracture modes I and II.

COMPUTATIONAL TECHNIQUE
he technique developed is divided into two analyses: Macro, which shows the place where the stress peak occurs; and Micro, which shows the number of critical fatigue life cycles. In this case, there is a main program that performs the interaction between the global and local analyses, as well as the link between the BemLab2D interface and the processing by BemCracker2D.

Macro Analysis -pre-processing
Initially, the modeling of the fuselage is performed in BemLab2D, obtaining the data referring to the physical parameters and the boundary element mesh, such as the nodal coordinates, mesh topology, traction and displacement boundary conditions, etc. As an example, the model in Fig. 4 was adopted, which consists of a plate with normal P and shear Q external loads and displacement restriction in the other nodes. From the initial model in Fig. 4, the algorithm creates the internal points for the calculation of the stress field, by BemCracker2D. The generated set of internal points is only intended to obtain a uniform grid of points to enable the plots of stress field and, therefore, it is not a mandatory requirement of the BEM.

Macro Analysis -post-processing
After obtaining the internal stress fields of the macro model, by BemCracker2D, the stress fields σ x , σ y and τ xy are obtained, as illustrated, respectively, in Figs. 5, 6 and 7. With these stress fields, one can find the critical location, or a peak that considers the three stresses together. For this, the von Mises criterion is used -see Fig. 8. It is noteworthy that the von Mises criterion was adopted only to identify the peak, still in the elastic regime; there is no plasticity analysis, since the algorithm is recommended for the elastic regime. Thus, if other criteria such as Tresca are adopted, for example, there will be no significant change in the position of the peak. From the stress peak obtained in the macro model highlighted in Fig. 8, the micro model analysis is shown in the following item.

Micro Analysis
The micro model represented has a square shape with an edge size of 1 cm with initial flaws characterized by a central hole of radius r and two cracks with sizes L1 and L2 representing the upper and lower cracks, respectively. The micro-model size and the initial defects are among the main problems in a multi-scale analysis, especially in techniques in which physical quantities are transferred from the micro-scale to the global scale, as performed in homogenization techniques [51]. Numerous researchers have addressed this problem over the years [52][53][54][55][56]. In this work, the values used reflect the effective size of supposed initial faults, however, a probabilistic methodology with BEM -using BemCracker2D, allows the initial defects to be modeled under different sizes. Thus, the micro analysis is performed at the peak location of the stress field, since, as this is the most unfavorable situation, all other positions would naturally be met. Continuing with the example adopted in [35], the peak is shown in Fig. 8 and the positioning of the micro model is shown in Fig. 9 (a). In this case, the stresses are taken from the stress field at this point and applied on the right and upper edge with a fatigue load ratio (R) equal to 0.5 and displacement constraint on the left and lower edges, Fig. 9   From the micro model, BemCracker2d processes the incremental analysis and, at each crack advance, it brings as results the number of cycles and the mesh displacement, necessary for the compliance calculation from the average of the displacements of each edge and the respective tension on the considered edge (right or upper). Fig. 10 (a) shows the propagation path and the respective number of cycles and, and Fig. 10 (b) shows the mesh deformation. In Fig. 11, each point, from left to right, represents an increment of the crack from Fig. 10 (a). The initial compliance, before the crack propagation, is on the abscissa axis, and, with each propagation increment, the plate loses rigidity and exponentially increases the compliance versus cycles ratio. From these points, a curve fit corresponding to a spline is created, resulting in the curves in Fig. 12, from which it is obtained the number of cycles corresponding to 2x initial compliance and 3x initial compliance. The number of cycles referring to 3x initial compliance is considered unstable, as from this point onwards compliance tends to infinity. The results in Fig. 12 refer to Tip 1 of Crack 1 (Crack 1 Tip 1), see Fig. 9 (b). Fig. 13 presents the cycle number results for the four crack tips. For the purpose of damage tolerance, the smallest number of cycles that reach instability in 3x initial compliance (N 3C ) should be considered.

APPLICATION
his section presents two applications based on the model presented in [57], in order to demonstrate the accuracy of the proposed methodology, as well as to validate the robustness of the BemCracker2D program as part of it. For simplification purposes, the results of the macro and micro models regarding the mesh, stress field and deformed shape will not be shown, but only those related to the stress peak location, number of cycles, compliance and the respective curves. In both applications, P and Q represent the normal and shear loads, respectively; r is the radius of the central hole; L1 and L2 are the size of the top and bottom cracks, respectively.

A1 Application
Tab. 1 presents the macro and micro values for the first application, whose external request produces the location of the voltage peak as illustrated in Fig. 14. Next, the algorithm positions the micro element at the peak composed of the stresses arranged in Tab. 2 and executes the fatigue test to evaluate the respective number of cycles at each propagation crack until the last increment before reaching an edge. Finally, at each increment, there are the points for the construction of the fatigue life curve (N) versus the average compliance of the edges, as illustrated in Fig. 15. Thus, as this element was executed at the point of stress peak, these curves correspond to the worst situation, the N 3C minimums. T For this application, the values of C=5e-11 and m=2.5 of the Paris's constants were considered to calculate the number of cycles. However, by varying these values, the number of cycles also varies, thus making it possible to create the curve that relates the variables C and m with the number of cycles, as will be shown below.

Curve N(C, m) for A1 Application
Adopting C and m from Paris's Law in the form of a grid in the domain C = [5e- 11, 9.5e-11] and m = [2.7, 3.2], the points that relate them to the number of cycles are obtained, as illustrated in Fig. 16   Next, through the interpolation of the points obtained in Fig. 16, the surface that represents the function that relates C in with the number of cycles is obtained, Fig. 17. With that, to find the values of the constants that carry the number of minimum cycles to the design value, the intersection between the relationship curves with the surface of the required number of cycles must be found. For example, considering the number of cycles defined in the project as 1e+04, the intersection is shown in the red line, as shown in Fig. 18 (a). In Fig. 18 (b), in the C x m plane, the curve that shows the parameters C and m that the material must have for the number of cycles requested by the user is obtained, thus originating the graph in Fig. 19 showing that it supports 10,000 cycles.

A2 Application
Tab. 3 presents the macro and micro values for the second application, while the location of the stress peak is seen in Fig.  21. The stress peak is displayed in Tab. 4 and, finally, at each increment, there are the points for the construction of the fatigue life curve (N) versus average compliance of the edges, as illustrated in Fig. 22.    Curve N(C,m) for A2 Application Fig. 23 shows the domain C = [5e- 11, 9.5e-11] and m = [2.7, 3.2] relating them to the number of cycles. Fig. 24 (a) illustrates the surface that represents the function relating C and m with the number of cycles. Thus, considering the number of cycles defined in the project as 1e+04, for example, the intersection is shown in the red line. In Fig. 24 (b), in the C x m plane, the curve that shows the parameters C and m that the material must have for the number of cycles requested by the user is obtained, thus originating the graph in Fig. 25, showing that it supports 10,000 cycles. Furthermore, by modifying the value of the number of design cycles n* (black color surface of Fig. 24), it is possible to obtain the variation of C and m for each value of n* requested in design, as shown in Fig. 26.

Comparisons between A1 and A2 applications
Here, the curves C and m between applications A1 and A2 are compared. Tab. 5 shows the values of the variables for each application and, in Fig. 27, the respective curves are illustrated. In this table, one can observe the difference in the choice of Paris' values when we change the input variables P, Q, r, L1, L2.

FINAL REMARKS
he proposed methodology is an alternative to the classical damage tolerance analysis in which the evaluation of damage is done through the critical crack size. In this research, the critical size is disregarded and compliance is evaluated as the defining variable of instability. The use of the Boundary Elements Method (BEM), as a numerical method, has been essential because of its flexibility, precision and mesh simplicity. With the BEM, the proposed numerical fatigue life technique was programmed without complex re-meshing processing as in the FEM. The suggested formulation was useful to evaluate both the location of the stress peak and the compliance of the edges of the local analysis elements. This research presents the relationship between the Paris' constants (C and m) and damage tolerance, through curves relating the Paris' constants to the desired number of cycles. A function is presented with the parametric data of C and m related to the desired number of cycles in a conservative way to avoid instability. The automation of the technique presented in this research with the use of the computer programs BEMLAB2D and BemCracker2D enable the analysis of fatigue tolerance in two-dimensional geometries and under plain stress or strain. The programs allow damage analyses throughout the parametric data of C and m of Paris' law to ensure damage tolerance and avoid the undesired Limit State. This research is an initial contribution and more work is needed and encouraged. T T