Numerical study of quasi-static and fatigue delamination growth in a post-buckled composite stiffened panel

In this work, an approach based on the Virtual Crack Closure Technique, included in the commercial finite element code ABAQUS, is adopted to study the propagation of delamination in composite structures under quasi-static and fatigue loads. The methodology, originally capable of simulating only delamination under quasi-static loads, has recently been extended introducing the possibility to analyze damage progression under fatigue load condition. The approach is assessed on simple specimens, Double Cantilever Beam and Mixed Mode Bending test, comparing the results with literature data. Afterwards, the behavior of a single-stringer specimen with an initial delamination is numerically investigated considering compressive loading conditions. At first, the single-stringer specimen is analyzed under quasi-static compressive load showing a clear correlation between local buckling phenomena and delamination growth. Then, a cyclic compressive load is applied such that the specimen switches between preand post-buckling conditions in a single load cycle. The outcomes of the numerical analyses are compared with the experimental data obtained from an experimental test campaign previously performed, showing the advantages of the adopted numerical technique but also the limitations that need to be addressed to properly analyze this phenomenon. Keyword: Fatigue; Crack Propagation; Finite Element Analysis; Post-buckling; Virtual Crack Closure Technique.


Introduction
Modern thin-walled aeronautical structures, such as fuselage and wing structures, make extensive use of composite stiffened panels where the stringers are typically joined to the skin through adhesive bonding or co-curing. In most cases the buckling load of these structures is much lower than their failure load [1]. However, the out-of-plane displacements in the post-buckling regime may become considerably large, and, due to their cyclic occurrence, may induce separation in the interface between the skin and the stringer. The lack of reliable numerical procedures able to accurately predict this phenomenon forces to consider the buckling load as limit load during the design of aerospace stiffened panels. Allowing the structure to buckle in some well-defined service conditions would result in an increase of the load carrying capability or in a reduction of the structural weight at the same limit load. However, in this case, particular attention must be paid to fatigue delamination initiation and propagation to ensure that the panel does not fail prematurely.
Although the problem has been extensively studied at the coupon level, the prediction of fatigue delamination propagation in post-buckled stiffened composite panels is still an open issue. Under fatigue loading conditions, delamination growth is usually characterized using the power law proposed by Paris et al. [2,3], which relates the crack growth rate to the stress intensity factor or equivalently to the energy release rate. This law is widely accepted in literature and overall provides an excellent approximation of the experimental results. Nowadays almost all available numerical methodologies used to simulate fatigue delamination growth rely on the Paris law or on one of its variations.
According to how the delamination is numerically represented and how the Paris law is implemented, the numerical approaches can be divided in two main categories, namely, damage mechanics and fracture mechanics [4].
In the context of damage mechanics, interface elements with an embedded cohesive law have been widely adopted for the simulation of delamination onset and propagation under quasi-static and impact load conditions [5][6][7][8]. Following the good results obtained for these problems, a few authors have extended cohesive formulation to simulate fatigue crack growth. However, inside the cohesive law, the energy release rate is not directly defined and this has led to the development of a variety of models that relate the cohesive damage variable to the crack growth rate defined by the Paris law [9][10][11][12]. The majority of these models provide good results in analyzing fatigue crack propagation in small coupons. However, it is difficult to apply them to more complex structural problems due to the large number of interface elements required to adequately represent the fatigue cohesive zone resulting in unfeasible computational times.
On the other hand, numerical approaches derived from the fracture mechanics are based on the direct application of the Paris law in conjunction with a methodology for the calculation of the energy release rate, such as the Virtual Crack Closure Technique (VCCT) [13][14][15]. This approach, developed for simulating delamination propagation under quasi-static load conditions [16], is now implemented into several commercial Finite Element (FE) codes. A number of issues with the original VCCT formulation, such as the problem of orthogonality between the crack front and the structural mesh or the bi-material interface simulation, have been identified and addressed by a variety of researchers across a large number of publications [17][18][19].
Although this methodology has still some limitations when dealing with fatigue delamination problems [20] and requires an initial damage, it currently represents an effective solution for the simulation of fatigue crack growth in structural components.
The aim of this work is to study the delamination growth behavior in a composite single-stringer specimen with an initial delamination, subjected to cyclic load in postbuckling regime, using a VCCT based approach recently introduced in the FE Code ABAQUS [21]. The Single-Stringer Compression (SSC) specimen analyzed in this work has been designed and tested in [22][23][24][25][26] to be representative of the response of large multi-stringer panels such as those of aircraft fuselage structures. In the building block approach, commonly used for aerospace composite structures, the specimen fits between the coupon specimen level and the large multi-stringer stiffened panel level.
The SSC specimen has been designed to have a small size and a relatively low manufacturing cost, while having a structural complexity on the level of multi-stringer panels. The intent is to improve the ease and affordability with which the response of stiffened panels can be experimentally determined, enabling the study of their failure behavior and damage tolerance in the post-buckling regime. This specimen provides the opportunity to verify and validate quasi-static and fatigue damage models due to its relatively complex geometry compared to coupon specimens, and its small size, which makes it computationally tractable.
In section 2 the theoretical background of the adopted numerical formulation are illustrated, while, in section 3, simple specimens such as the Double Cantilevered Beam (DCB) and the Mixed Mode Bending (MMB) test specimens are numerically analyzed to verify the effectiveness of the approach. Finally, in section 4, the model of the SSC specimen is presented together with the results of the quasi-static and fatigue numerical analyses and is compared with experimental data.

Theoretical Background
The Paris Law used to characterize the fatigue crack growth rate is given, in its most basic form, in equation (1): where C and m are fitting parameters and f(G) is a function of the energy release rate (G). The typical behavior of the crack growth rate is qualitatively shown in Figure 1. The curve in Figure 1 can be divided into three distinct regions. In region I, a sharp drop in the crack growth rate occurs when the energy release rate approaches the threshold value (Gth), which depends on the material and on the loading conditions. In region II the behavior is almost log-linear and is captured by the Paris law. In region III, when the energy release rate approaches the critical value (Gc), the crack growth rate drastically increases.
Many different functions of G can be found in literature [4,27] and there is no clear consensus on which function better describes the experimental data. Most of the formulations tend to adopt Gmax, or ΔG=Gmax-Gmin, respectively, the maximum value of the energy release rate and the variation of the energy release rate during the load cycle.
However, as pointed out by Pascoe et al. in [27], since the delamination growth is ultimately driven by the applied fatigue cycle, which cannot be described by a single parameter, any equation which use only one parameter is not able to fully capture the delamination growth behavior. Furthermore, the use of ΔG can result in an erroneous interpretation of experimental data since it violates the principle of similitude. Based on the analogy with the stress intensity factor variation (ΔK), adopted for fatigue crack growth in metal, according to [27], the parameter ((Gmax) 0.5 -(Gmin) 0.5 ) 2 seems to provide a better description of experimental results.
Despite these limitations, the large majority of experimental data available in literature are reported using the variation of the energy release rate during each load cycle, as shown in equation (2): The methodology adopted in this work and included in the FE code ABAQUS [21] is based on equation (2) for the calculation of the crack growth rate and on the VCCT equations to evaluate the energy release rate.
The VCCT is based on Irwin's assumption [28] that the strain energy released when a crack extends by a small amount Δa from a to a+Δa is equal to the work required to close the crack to its original length, a. In the framework of a finite element analysis, assuming that the extension Δa does not change significantly the state of stress along the crack front, it is possible to evaluate the energy release rate considering the forces acting on the node at the delamination front and the displacements of the nodes immediately behind it, as shown in Figure 2 for a 2D problem.
where X1 and Z1 are the shear and opening forces at crack tip (node 1), Δu2,3 and Δw2,3 are the relative shear and opening displacements at upper and lower crack face in the nodes behind the crack tip (node 2 and 3) and Δa is the increment in the crack length.
The propagation condition in a quasi-static analysis is reached when the sum of all the components of the energy release rate reaches the fracture toughness of the material.
Several criteria exist for the calculation of the fracture toughness in a general 3D case involving all the three opening modes. One of the most adopted methods is the Benzeggagh-Kenane criterion [29], shown in equation (4): III  c  Ic  IIc  Ic  I  II  III where GIc and GIIc are the fracture toughness for pure mode I and mode II, while η is a fitting parameter obtained from experimental data. The propagation criterion adopted in ABAQUS can be expressed as follow: where ftol is the release tolerance, whose default value is 0.2. It is possible to reduce this value to improve the accuracy of the analysis but this results in an increase in the computational times.
Once the propagation criterion is met, by default, the connection between the two coincident nodes on the delamination front is instantaneously released. It is also possible to define a ramp behavior in which the constrained force is linearly reduced to zero with the opening displacement, according to the value of the fracture toughness.
This non-default approach allows to improve the convergence ratio of the analysis and to better capture rounded delamination fronts. However, ramp debonding is not currently supported for fatigue analyses where only the instantaneous release is available.
The fatigue crack growth analysis capability in ABAQUS adopts the same VCCT equations to evaluate the energy release rate. The procedure requires to define the load history of a single load cycle. Because only constant amplitude fatigue loads can be simulated, the load cycle remains the same up to the total number of cycles specified by the user. The solution of each load cycle is obtained using the static nonlinear algorithm of the ABAQUS/Standard solver. The values of Gmin and Gmax, and therefore ∆G, are calculated using the VCCT equations when the structure is subjected respectively to the minimum and maximum load during the load cycle. The onset of fatigue delamination growth is determined with a criterion based on both Gmax and ∆G. The first part of the criterion consists in an exponential curve fit of the experimental fatigue crack onset data for the material being tested. In addition, the value of Gmax must exceed the threshold energy release rate of the material (Gth). The conditions that must be met for delamination to begin propagating are summarized in equation (6): where N is the current cycle number and c1 and c2 are fitting parameters experimentally determined. Once the onset criteria are satisfied, the fatigue delamination growth rate is governed by the Paris law shown in equation (7).
where c3 and c4 are fitting parameters.
The algorithm adopts a damage extrapolation technique to avoid the simulation of the entire load history. In particular, if at the end of the load cycle N, the onset criteria are satisfied in at least one point along the delamination front, the crack length is extended from the actual value aN to aN+∆N by releasing one node. The number of cycles required for the crack to grow is evaluated at each node along the crack tip reversing the Paris law, as shown in equation (8): where ∆Nj is the number of cycles required to release the node J, and ∆aNj is the crack extension that produces the release of the node J. The procedure releases at least one node at each load cycle by choosing the node with the minimum value of the cycles number evaluated with equation (8). This value also represents the number of cycles that are jumped in the following increment. The procedure is schematically illustrated in It is possible to accelerate this process defining a non-zero tolerance ∆DNtol which allows to release in the same cycle more than one node as long as the relation expressed in equation (9) is satisfied: The damage extrapolation tolerance determines how many nodes are released at each load cycle taking into account how close, in terms of number of cycles at failure (ΔNj), any other node on the delamination front is with respect to the weakest node (ΔNmin). A less stringent tolerance reduces the computational time, since more nodes can be released in a single load cycle, however the accuracy of the analysis decreases. The default value of ∆DNtol is 0.1.

Numerical Analysis of Coupon Tests
The numerical procedure described in the previous section is adopted at first to analyze two test cases, the Double Cantilevered Beam (DCB) and the Mixed Mode Bending

Analysis of DCB Test
The DCB specimen is a widely recognized standard to determine the mode I interlaminar fracture toughness [33] and mode I fatigue crack growth [34] for unidirectional fiber composite materials. The geometrical characteristics of the DCB specimen analyzed in this study are taken from literature [30] and shown in Figure 4.  The specimen is made of T300/1076 graphite/epoxy with 24 unidirectional plies and an initial delamination positioned in the center of the thickness. The properties of the material adopted in the numerical analyses are taken from literature [30] and shown in   The load-displacement curve of the quasi-static analysis performed on the DCB model is presented in Figure 6, along with the benchmark results from Krueger [30]. The model shows a linear response until the delamination begins to propagate at the peak reaction force. There is a good agreement between the numerical results and the reference data, both in terms of initial stiffness and peak load, which differs from the benchmark by less than 2%. In Figure 7  The delamination during the propagation shows a slightly rounded crack front, as expected from comparison with literature experimental results [30].
The same DCB specimen is then analyzed under fatigue loading conditions using the parameters shown in Table 2, taken from literature [31].
where Gpl is the Paris limit energy release rate when the crack growth rate reaches the unstable growth region. If the energy release rate is below Gth no delamination growth is assumed to occur, while, if it is above Gpl, quasi-static delamination growth is assumed to occur and the fatigue delamination growth analysis is terminated. The default value of 0.1 is adopted for the fatigue release tolerance ∆DNtol.
The analysis is divided into two steps: a pre-load step, where the opening displacement The delamination length measured at the center of the specimen as function of the number of cycles is shown in Figure 9 comparing it with benchmark results taken from literature [31]. As expected with a displacement controlled analysis, the delamination length rapidly increases in the first load cycles, then the growth rate decreases until the energy release rate becomes smaller than the growth threshold and the delamination stops propagating around 3.7 million cycles. From Figure 9, the agreement between the numerical results and the benchmark data can be appreciated, with the two curves deviating by less than 1%. The delamination front at different cycles is shown in Figure 10.

Analysis of MMB Test
The MMB specimen is one of the standard specimens to determine mixed-mode crack propagation for static and fatigue loading conditions [35]. The geometrical characteristics of the specimen and of the load fixture, considered in this study are taken from literature [32], and are shown in Figure 11 The specimen has 24 unidirectional plies of carbon epoxy IM7/8552 with a delamination in the middle of the thickness. The material properties are taken from literature [32] and are reported in Table 3.   The same mesh discretization adopted for the DCB is used to model the MMB  Figure 13 compared with benchmark data taken from literature [32]. A linear response is observed until the delamination propagation occurs and the reaction force decreases. When the crack front reaches the lever contact point at half span of the specimen, the change in the contact conditions leads to an increase in the reaction force.
The stiffness of the model as well as the behavior after the propagation agree well with the benchmark results. The deformed shape of the specimen and the delamination front at the end of the analysis are presented in Figure 14. In Figure 14 it can be noted, as already seen for the DCB, that the numerical analysis predicts the delamination front shape observed in experimental tests [32].
To analyze the specimen under fatigue loading conditions the parameters shown in Table 5, taken from literature [32], are adopted. The maximum displacements were determined to produce a maximum value of the energy release rate at the delamination front equal to 60% of Gc. A release tolerance ∆DNtol of 0.001 is adopted for the fatigue analysis.   Also in this case, there is a good match between the numerical results and the benchmark data, although a quite large deviation can be observed at the end of the analysis. In both cases, the differences begin to occur when the crack length reaches the mid-span length of the specimen, suggesting that a more accurate modeling of the load fixture including friction effects may be needed.

Single-Stringer Compression Specimen
In this section, the quasi-static and fatigue response of the SSC specimen, developed and tested in [22][23][24][25][26], is numerically investigated. The specimen consists of an omegashaped stringer co-cured with a skin panel. An initial delamination of 40 mm is created under one of the stringer flanges in the middle of the specimen, using a Teflon insert during the manufacturing process. The geometrical characteristics of the stiffened panel are shown in Figure 16.  for the stringer. Although the material system is the same as the MMB specimen analyzed in the previous section, the properties reported in [30] are slightly different from those shown in Table 4 due to two different versions of the same material distributed in Europe and in the United States. The material properties from [25], reported in Table 6, are adopted in this numerical simulation.    shown any significant improvement while resulting in an increase in the computational time. Figure 19 shows the obtained load-displacement curve and compares it to the experimental data measured during two quasi-static tests [24]. The initial trend of the obtained load-displacement curve reported in Figure 19 matches the experimental curves, although a difference in stiffness between the numerical results and the experimental data can be observed. This may be due to one or a combination of different types of imperfections. The test specimens contain geometrical imperfections such as variation in the potting length, alignment of the loaded surfaces of the potting, alignment of the specimen to the loaded surfaces and flatness of the specimen. In addition to these imperfections part of the difference may also be caused by the assumed boundary conditions to model the potting. Furthermore, from Figure 19 it is possible to notice that the numerical model quite underestimates the first load drop and therefore the delamination growth onset. After this point, at an applied displacement of around 0.4 mm, the numerical results and experimental data start to deviate from each other.
This is due to the interface properties adopted to simulate the crack propagation, which are measured performing experimental tests according to the ASTM standards [33][34][35] on coupons with delamination positioned between 0° plies. In the specimen under investigation, the initial delamination is located in the 45°/-45° interface resulting in higher values of the fracture toughness and higher delamination growth resistance. In Figure 22, the out-of-plane displacements obtained by the analysis are compared to the DIC data measured during the test just before the delamination starts to propagate [25]. As it can be noted, the numerical and experimental post-buckling shapes match very well. The skin as well as the skin portion under the stringer buckle with three halfwaves, however, the direction of the numerical displacements is opposite with respect to the experimental data. In Figure 23, the comparison in terms of out-of-plane displacements is presented after the opening of the delamination, just before the global collapse of the specimen. The post-buckling shape obtained from the analysis perfectly agrees with that one measured with the DIC. On the delamination side of the specimen, the skin presents a large single half-wave, which extends for almost the full length of the specimen, while on the opposite side exhibits a three half-waves mode shape.
After the quasi-static investigation, the SSC specimen is analyzed under fatigue loads.
The analysis is performed under load controlled condition to reproduce the experimental procedure. The two experimental tests were conducted at a maximum load of 23 kN, however, because this load is higher than the load at which the delamination starts to propagate according to the quasi-static analysis, it cannot be used for the fatigue propagation analysis because the delamination would propagate entirely after the first cycle. For this reason, a maximum load of 13.8 kN, just before the beginning of quasistatic crack propagation, is adopted. The fatigue analysis of the SSC specimen is performed using the fatigue crack growth approach adopted for the DCB and MMB specimens. A triangular load cycle is defined with maximum load equal to 13.8 kN and minimum load of 1.38 kN. Although, as shown in Figure 21, the mode-mixity changes along the delamination front and during the propagation, the fatigue parameters are not updated during the ABAQUS analysis and, therefore the parameters for 20% mixedmode, taken from literature [32] and reported in Table 7, are adopted for this problem.
This value of the mixed-mode is chosen because it represents a good approximation of the average mode-mixity through the whole analysis. The default value of 0.1 is adopted for the release tolerance ∆DNtol.  In Figure 24, the out-of-plane displacements of the specimen are reported at minimum and maximum load of the first cycle and at maximum load after 28836 cycles, when the fatigue analysis is terminated. From Figure 24, it is possible to appreciate how, during a single load cycle, the SSC specimen oscillates between pre-and post-buckling conditions. Furthermore, the increase of the out-of-plane displacements at the end of the analysis with respect to the first cycle is evident. This is due to the growth of the delamination which causes a reduction of the global stiffness of the specimen and an increase of the applied compressive displacement. At the beginning of the analysis, the flange of the stringer follows the buckling direction of the skin while, as the delamination length increases, right before the last load cycle, it starts to buckle in the opposite direction resulting in a rapidly increase of the energy release rate at the delamination front. The delamination front is shown in Figure 25 at different load cycles compared with the experimental image taken using an ultrasound system after 24,000 cycles. The comparison with experimental results shows a quite large underestimation of the delamination length. Several causes can be identified to explain the lower growth rate of the numerical analysis. The first and most important is that the numerical compressive load is much lower than the load applied during the test. The second issue regards the coefficients of the Paris law which are defined in advance and kept constant during the analysis and therefore do not take into account possible changes in the mode-mixity or stress ratio. Finally, more experimental data of the SSC specimens are needed to validate the numerical results. Test data from only two specimens were available, and they exhibited a substantial scatter in the observed fatigue delamination growth rates.

Conclusions
The numerical investigation of delamination growth under quasi-static and fatigue load using an approach based on Virtual Crack Closure Technique and available in the finite element code ABAQUS has been conducted. At first, Double Cantilever Beam and the Mixed Mode Bending specimens have been analyzed under quasi-static and fatigue loads to evaluate the capabilities of the numerical procedure. In both cases, the adopted methodology has produced good results in terms of crack length as function of the number of cycles compared to data taken from literature. Then, a single-stringer specimen has been analyzed under both quasi-static and fatigue loading conditions and the results have been compared with experimental data. The numerical approach has been proven to be capable of predicting delamination growth in a post-buckled singlestringer compression specimen, however, a few limitations have been encountered.
During the quasi-static analysis, the lack of experimental data regarding the value of the interface properties for delamination positioned between different oriented plies has led to an overestimation of the crack propagation. In the fatigue analysis, the use of a single set of Paris law parameters, related to specific values of load ratio and mode-mixity without taking into account their variation during the analysis resulted in an underestimation of the crack growth rate. Despite these limitations, the results obtained from the numerical analysis are qualitatively similar to the experimental data, showing the potential of the adopted numerical approach for the simulation of delamination growth under fatigue loading condition in complex structural problems.