Damage-Tolerant Design of Stiffener RunOuts : A Finite Element Approach

Modern aerostructures are predominantly of semi-monocoque construction characterized by a thin skin and stiffeners. The latest generation of large passenger aircraft also use mostly carbon-fibre composite material in their primary structure and there is a trend towards the utilization of bonding of subcomponents in preference of mechanical fastening. Current design philosophy requires that certain stiffeners are terminated, for example due to an intersecting structural feature or an inspection cut-out. In these circumstances, the loading in the stiffener must be diffused into the skin, leading to complex three-dimensional stressstates. The development and utilization of reliable virtual component testing, in the design of composite aerostructures, can potentially yield significant cost reductions. Such reliability requires a thorough understanding of the damage mechanisms and failure processes in realistic aerostructures, particularly in critical regions such as stiffener run-outs.


Introduction
Modern aerostructures are predominantly of semi-monocoque construction characterized by a thin skin and stiffeners.The latest generation of large passenger aircraft also use mostly carbon-fibre composite material in their primary structure and there is a trend towards the utilization of bonding of subcomponents in preference of mechanical fastening.Current design philosophy requires that certain stiffeners are terminated, for example due to an intersecting structural feature or an inspection cut-out.In these circumstances, the loading in the stiffener must be diffused into the skin, leading to complex three-dimensional stressstates.The development and utilization of reliable virtual component testing, in the design of composite aerostructures, can potentially yield significant cost reductions.Such reliability requires a thorough understanding of the damage mechanisms and failure processes in realistic aerostructures, particularly in critical regions such as stiffener run-outs.
When a stiffener is terminated, the loads which it carries must be transferred to the skin, making the design of the run-out region vital, hence improved design methodologies are required.Several studies of skin-stiffener failure [1][2][3][4][5][6][7][8] have been carried out.Falzon et al. [7,8] investigated the failure of realistic stiffener run-outs loaded in uniaxial compression.Different stacking sequences and skin thicknesses at the run-out region were tested and a wealth of complexity in the response and subsequent failure was reported.For all tests, failure initiated at the edge of the run-out and propagated across the skin-stiffener interface.It was found that the failure load of each specimen was greatly influenced by changes in the geometric features of these specimens.Falzon and Hitchings [8] used a Virtual Crack Closure Technique (VCCT) [9] to predict the crack growth characteristics of the modelled specimens and reported shortcomings in the quantitative correlation between the predicted and observed failure loads and modes.
Bisagni et al. [10] investigated, experimentally and numerically, the postbuckling response of hat-stiffened composite panels using cohesive zone models and predicted collapse loads which were in agreement with experiments.Krueger [6] used VCCT and a numerically effective shell/3D modelling approach to predict the delamination failure of skin-stiffener run-outs.Camanho et al. [11] implemented a cohesive element to numerically investigate the debond strength of skin-stiffener composite specimens using cohesive zone models and compared with experimental results.
In a recent study by the authors [12], a parametric numerical analysis was conducted to optimize the design of the run-out section to increase the crack growth stability under axial compression.Improved damage tolerance (stable crack propagation) was reported in the modified stiffener run-out design as compared to the baseline configuration.The modified design eventually failed catastrophically by interlaminar delamination, not bondline failure (debonding), which had not been considered in the numerical study.A more detailed analysis of different configurations, which accounts for delamination, was therefore undertaken in the work presented in this chapter.Building on the previous findings, the merits of a Compliant termination scheme are presented.
The research in this paper focuses on stiffener run-outs loaded in compression with a selection of stiffener termination schemes.These schemes are analysed numerically (Section 2) in order to compare the influence of the design on the energy release rate for debonding and delamination.The Baseline, Tapered and Compliant schemes were manufactured (Section 3) and tested to failure (Section 4).The experimental results are presented in Section 5, and compared to the predictions in Section 6.
The failure mode of the Tapered specimen configuration was found to be a combination of interlaminar and intralaminar failure [12].The Tapered stiffeners had an unexpected delamination in the flange between the bottom-most 0 o ply and above the 45 o ply.Furthermore, the delamination led to an intralaminar failure in the form of a matrix crack across the 0 o ply near the filler ends and continued delaminating between the filler and 0 o ply, Figure 19a.For this reason, a new finite element (FE) model was created in order to investigate these modes of damage, Figure 19b and c.

Skin-stiffener configurations
The main focus of this study was a new configuration to improve the crack growth stability of a stiffener run-out.The structural performance of a baseline skin-stiffener configuration under longitudinal compression, with geometry and dimensions shown in Figure 1a, was compared to that of a modified parametric configuration shown in b.The modified configuration has a widening flange towards the termination end of the stiffener but this added material is offset by the taper of the stiffener web.This results in a stiffener design with a similar overall weight to the baseline design.For the parametric configuration, various values of b, c and d were analysed.The materials used in this study were IM7/8552 carbon/epoxy pre-preg, with ply thickness 0.25 mm, for the skin and the stiffener, and FM300 adhesive film (0.15 mm thick) for the bondline.Material In previous test results with geometry similar to that in Figure 1a, the specimens failed by unstable debonding of the stiffener from the skin.Therefore, the different configurations in this study were assessed by comparing the energy release rates of the run-outs for a given displacement and for several debond lengths.

The FE model
All the FE simulations of the parametric study were carried out in ABAQUS [13] and the parameterized models were created in Python [14].The main model has five different parts, the skin, the adhesive between the skin and the stiffener, the parts of the stiffener and the filler.The material properties for IM7/8552 measured using standard tests and are shown in Table 1.A mesh sensitivity study was carried out to ensure that all results presented are mesh converged.

Element Idealisation
All parts were modelled with three dimensional hexahedral solid elements C3D8 to accurately capture stresses in the through-thickness direction.Also, solid elements are capable of modelling several layers of different materials for the analysis of laminated composites which is ideal for our numerical study.In consideration of computational cost, linear order elements were chosen and to ensure higher fidelity results whew required, the models executed using a mesh density based on a ability to define a detailed mesh convergence study.One of the major advantages of solid elements is the feature of defining several ply stack ups within a single element along the stacking direction with an improved element formulation, illustrated in the Figure 2 below.All the models were meshed with fully integrated solid elements (C3D8) rather using reduced integration elements (C3D8R), since the reduced integration elements were not able to capture damage variable calculations at all the nodes while delaminating, Figure 3. On further investigating this issue, it was observed that the problems arose because of hourglassing.Even though this particular problem could be partially mitigated by using the hourglass stiffness enhanced element formulation, the use of C3D8R elements in the numerical models was avoided for greater accuracy.Consequently, C3D8 elements were used in all the models.

Stacking Sequence
The skin consists of eight plies, while the stiffener consists of five plies.The stacking sequences used in [1] for the skin and the stiffener, are shown in Table 2. Figure 4 shows the ply orientations for the skin.The thickness of each ply is 0.25mm (double thickness) and the number of integration points is set to three within each ply/element through the thickness.

Mesh
The number and distribution of elements is a crucial step when using the FE method to solve structural problems.Figure 5 shows the element families that are most commonly used in stress analysis.The main difference among the element families is the geometry type that each of the family represents.As mentioned earlier, an 8-node linear brick element was used and a course mesh of the modified design is shown in Figure 6.

Constraints
Constraints were defined in the Assembly module [13] for the initial positions of instances.The type of constraint created for this model was a "tie constraint", which allows the user to fuse together two regions even though the meshes created for sectors may be dissimilar [8].For the current model, sixteen constraints were created between all the parts being in contact.In Figure 7 it can be seen how the master and the slave surfaces are displayed in the model.

Boundary conditions
In the current study, two boundary conditions were defined.The first condition, named BC-1, was a displacement/rotation type boundary condition.All displacements and rotations were set to zero.BC-1 was applied to one edge of the stiffener as shown in Figure 8.The second BC, named BC-2, was a displacement/rotation type BC as well, where all displacements and rotations were set to zero apart from U3, which was set to U3=1.U3 is the displacement in the axial direction and represents the testing machine's displacement.

Mesh sensitivity study
Finite element analysis provides an approximate solution and it can only guarantee that equilibrium is satisfied on average over an element.This does not mean that it will satisfy equilibrium over any smaller volume compared to a complete element.As a consequence, equilibrium is enhanced when the size of the element is decreased.Moreover, in regions of stress concentrations, it is necessary to increase the accuracy of the FE solution by either using elements with higher-order shape functions (p-refinement) or by using a finer mesh of elements (h-refinement).The goal that a designer needs to achieve is to select the best mesh density which is not prohibitively expensive to run and at the same will provide accurate and acceptable results [15].
In this study, three different meshes were used (a coarse mesh, an intermediate mesh and a fine mesh).13100 elements were used for the coarse mesh, 18700 elements for the intermediate mesh and 51000 elements for the fine mesh.Table 3 illustrates the three different meshes used, the strain energy and the running time of each model.According to the values of  Moreover, investigation was undertaken to assess the accuracy of using multiple ply orientations within a solid finite element.Abaqus® provides the option to make a partition through the ply thickness and define the material orientation.As a result, eight partitions were created for the skin and five for each stiffener using material orientations given in Table 2.The results obtained from both approaches were similar, as illustrated in Figure 9.

The Python script
In order to perform the parametric study, a large number of models were created and a script that generates these models was needed.This script was written in Python [14] because of the advantage of using the Abaqus scripting interface to generate automate repetitive models and execute.
The parametric script generates the stiffener run-out models with different crack lengths.
The script uses a basic stiffener run-out configuration and changes the design each time according to the parameter values.When a new design is generated, the script propagates the crack and the energy release rate is calculated for every step.

Energy Release rate along crack
After converging to an optimum mesh density, failure in the stiffeners was captured by consecutive models, propagating a delamination crack or a debonding crack, 1 mm at a time.The results of the parametric study are shown in Fig. 10, where the values of  Interlaminar and intralaminar failures were observed in the modified run-out stiffener.The stiffener had an unexpected delamination between the 0 0 ply and the 45 0 ply.After examination of the specimens it was assumed that the delamination led to an intralaminar failure in the form of a matrix crack across the 0 0 ply near the filler ends and continued delaminating between the filler and 0 0 ply.In Figure 11, the failed specimens clearly shows that delaminations and matrix cracks occurred in the 0 0 ply on both sides of the specimen.

Modeling the debonding of the stiffener
The finite element model of the stiffener run-out had the following features: 1. Skin -single part of (45 0 /-45 0 /0 0 / 90 0 )S laminate 2. Filler -made-up of 0 0 fibres 3. Left -0 0 lamina of the stringer 4. Left -single part of (0 0 / 90 0 /-45 0 /45 0 ) laminate of the stringer 5. Right -0 0 lamina of the stringer 6.Right -single part of (0 0 / 90 0 /-45 0 /45 0 ) laminate of the stringer The part descriptions are clearly illustrated in a front view of the skin-stiffener assembly, Figure 18.The interface between the skin and the stiffener was connected by using a surface-based cohesive layer, as seen in Figure 12.Cohesive zone modelling is generally used for the numerical simulation of interlaminar failure.Damage initiation is driven by a traction seperation law and the value of the maximum traction, t o , Figure 14(a).New crack surfaces are formed when the fracture toughness Gc is equal to the area surface under the tractionseparation curve.Considering the nature of predictions made by the parametric study, the model was analysed with cohesive interaction properties of adhesive FM 300 and listed in Table 4. Knn in Nmm -3  Kss in Nmm -3  Ktt in Nmm -3

Initial linear elastic behaviour
In order to accurately predict damage initiation, the interaction of traction components was taken into account and the quadratic stress criterion used.This criterion is included in Abaqus and was formulated based on Ye's criterion [13] including the interaction between traction components.According to this criterion, the damage initiates when a quadratic interaction function reaches a value of one.The quadratic interaction function is shown in the following equation.
where, t 0 n , t 0 s and t 0 t are the peak values of the stress when the separation is either purely normal to the interface or purely in the first or the second shear direction respectively.In the present numerical study, these values for IM7/8552 are 50 MPa in the normal direction and 100 MPa in the two shear directions [10], see Table 4.
When the initiation criterion is met, the cohesive stiffness degrades at a rate defined by the damage evolution model.The overall damage of the damage zone is represented by a scalar damage variable, D, and is implemented in the damage evolution model.After the damage initiation, the damage variable monotonically evolves from 0 to 1 on increasing the loading.
For the mix-mode definition of fracture energies, Benzeggagh-Kenane's criterion [16] (BK law) was used which defines the energy dissipated as in equation below, ( ) with, where Gn, Gs and Gt are the fracture toughness values in the normal and two shear directions respectively which were measured in-house for IM7/8552 and given in Table 4.In this study the value of the BK mode-mixity power parameter, η=1.6, was obtained from Maimi [17,18].

Response of the Numerical Model
The debonding failure initiation load predicted by the model was compared with the debonding loads predicted by using the energy release rate and the experimental results, Table 5.

Cohesive layer Prediction
Failure  Recalling the results from the energy release rate analysis, when the maximum energy release rate across the width of the stiffener was used, the predictions were closer to the experimental results.In addition, when the average energy release rate was used, the difference from the experimental results had an over-prediction of 8.5%.By using the cohesive zone model, the difference with the experimental results was 5.6%.
Also, the debonding growth along the width of the specimen correlates well with the predictions using strain energy release rate across the width of the stiffener as can be seen in Figure 15.Comparing the FE model pattern with the Tapered slop on the right, it is clearly observed that there is a correlation between the extent of damaged and the value of normalized strain energy release rate across the width of the specimen.

Modeling the debonding failure using VCCT
The implemented VCCT method in Abaqus standard, developed by Boeing and Simulia, suggested promising results, especially when the mismatched, tie-constrained, meshes had only 3% error comparing with pairing meshes.As can be seen in Figure 16(a) there was good correlation between the VCCT and the parametric study.On the other hand, the VCCT method could not capture the increase in normalised GT/GC arising from the geometric discontinuity at the edge of the flange.In order to capture the detail at the edge of the flange, the mesh resolution was increased and was biased towards these edges, Figure 17, and the results can be seen in Figure 16(b).Despite the good results, the size of the model and the time needed for execution made the use of VCCT impractical and was not used in further investigations.

2 nd iteration
The objectives defined in the first part of this study were to create two FE models for the baseline and the modified configuration.The baseline stiffeners failed due to debonding of the skin/stiffener interface, while the modified stiffeners failed by delamination between the 0 0 and 45 0 ply interface.Improved damage tolerance (stable crack propagation) was reported in the modified stiffener run-out design as compared to the baseline configuration.The modified design eventually failed catastrophically by interlaminar delamination, not bondline failure, which had not been considered in the numerical study.A more detailed analysis of different configurations, which accounted for delamination, was therefore undertaken.Building on the previous findings, the merits of a compliant termination scheme are presented.

Skin stiffener runout configurations
The structural performance of three different skin-stiffener configurations -Baseline (B), Tapered (T) and Compliant (C) -under longitudinal compression, with geometry and dimensions shown in Figure 18, was assessed.Compared to the Baseline stiffener (Figure 18a), the other two configurations have a widening flange towards the termination end of the stiffener but this added material is offset by the taper of the stiffener web (Tapered configuration, Figure 18b).The third configuration includes a curved tape (Compliant, Figure 18c).
The proposed Compliant design was developed by considering potentially beneficial local stiffness variations.This resulted in a stiffener design with a similar overall weight to the Baseline design.For the modified configurations, various values of b, c and d (see Figure 18b) were analysed [12].
The failure mode of the Tapered specimen configuration was found to be a combination of interlaminar and intralaminar failure [12].The Tapered stiffeners had an unexpected delamination in the flange between the bottom-most 0 o ply and above the 45 o ply.
2nd order polynomial  Furthermore, the delamination led to an intralaminar failure in the form of a matrix crack across the 0 o ply near the filler ends and continued delaminating between the filler and 0 o ply, Figure 19a.For this reason, a new FE model was created in order to investigate these modes of damage, Figure 19b and Figure 19c.

Energy release rate along crack
The three different configurations, Baseline, Tapered and Compliant (see Figure 18), were analysed for debonding and delamination growth stability.The results of this analysis are presented in Figure 20, where the values of GT = GΙ + GΙΙ + GΙΙΙ, the total strain energy release rate, have been normalized to the GT of the reference parametric stiffener for 0.5 mm of crack length.Figure 20 compares the normalised strain energy release rate for the Baseline, the Tapered and the Compliant stiffeners.The negative slope of the G(a) curve indicates crack growth stability, while a positive slope indicates instability (assuming constant fracture toughness).Recalling the failure modes obtained experimentally [12] the Baseline stiffener failed by debonding and the Tapered stiffener initially experienced debonding until it finally failed by delamination.This is in agreement with predictions, Figure 20.Consequently, both models were able to correctly describe these experimental results [12].In addition, the stability analysis for the Compliant stiffener predicts that this design will fail stably by debonding, Figure 20.

Energy release rate along the width of the crack tip
The strain energy release rate along the width of the crack tip was calculated for the Baseline, Tapered and Compliant configurations (Figure 18).Fracture initiation is expected when the GT exceeds the fracture toughness Gc for a given mixed-mode ratio GII / GT at each point along the crack tip.In other words, propagation at each point occurs when GT / Gc >1 [19,20].The interlaminar fracture toughness Gc can be calculated by using the equation 2 [20] The value of Gc is normalised to the width-average value for the Tapered specimen.It can be observed that the trend is similar for the Baseline and Tapered specimen types but is different in the centre of the Compliant stiffener.This is due to the difference in the web of the stiffeners.The curved taper has reduced the normalized strain energy release rate in the centre without affecting the trend in the flange.The maximum value of the energy release rate can be used to predict the load corresponding to the initiation of fracture using where P is the load at initiation of fracture, PFE is the load from the FE model, Gc the critical strain energy release rate (Equation 1.26), and GT is the strain energy release rate predicted by the FE model as defined previously.Two different predictions for P can be made: one using the maximum value of GT along the width, and another using the average, Table 6.

Experiments
The tests were carried out in an Instron testing machine, equipped with a 100 KN load cell, at a loading rate of 0.5 mm/min.Load and crosshead displacement were recorded continuously by a PC data logger connected to the load cell and the Instron machine at a sampling rate of 2 Hz.The specimens were aligned by careful measurement in the loading direction to avoid bending.AE sensors were used to identify and investigate failure, within the specimens, during testing.The AE equipment was manufactured by Physical Acoustic Corporation (PAC) and failure was monitored by AEwin software.Broadband (WD) sensors with an operating frequency range of 100 Hz to 1000 kHz were used and positioned in order to obtain the best results without affecting the specimens behaviour [21].
The Baseline stiffeners had an average failure load of 16.5 kN while the Tapered stiffeners had an average failure load of 17.7 kN and the Compliant an average failure load of 18 kN, Table 6.The fracture surfaces for selected specimens are shown in Figure 23, and the load versus displacement curves for selected specimens of the three stiffener designs are shown in Figure 22.The predicted loads (using Eq. 2) match well with the experimental values when the maximum G across the width is used, Table 6.
The acoustic emission signals (Figure 22) show that there was an increase in AE activity 0.1 mm before catastrophic failure for the Baseline specimen.For the Tapered specimen type, the increase in AE emission started about 0.05 mm before catastrophic failure and for the Compliant specimen, 0.02 mm. Figure 22 also shows the peak frequency during the tests for all specimen types.It can be observed that there is some very low-energy micro-cracking from the start of the test and this is possibly at the resin pots.The Tapered specimen configuration promoted a combination of failure modes including delamination and fibre bridging which preceded catastrophic failure.In addition, the Compliant stiffener, according to AE data and as visually observed (Figure 23c), suffered only from debonding.

Remarks
The strain energy release rate analysis yielded good results in the investigation of the runout design influence in debonding/delamination for stiffener terminations.The FE models accurately predicted the failure loads and failure modes for the specimens tested and the predictions were improved when the distributions of the strain energy release rate across the width was considered.The differences in the predictions using the average and the maximum energy release rates are shown in Table 6 and can be compared to the experimental failure loads.The load-displacement, as well as the peak frequencydisplacement plots (Figure 22), show that the Tapered design is slightly more damage tolerant than the Baseline one and this improved further with the Compliant design.The AE monitoring proved to be valuable in detecting and analysing the failure modes experienced by the specimens.

Conclusions
This study was based on the strain energy release rates for debonding and delamination and successfully predicted the failure loads for the three different specimen types.The predictions were more accurate when the maximum strain energy release rate across the width was used.It can be concluded that the variation of the energy release rate across the width should be considered when stiffener run-outs are designed.AE data recorded during skin-stiffener run-out compression tests proved useful to analyse the failure processes which take place in these specimens.The results show that in the design of skin-stiffener run-outs it is important to consider the possibility of failure modes other than debonding, and that compliant termination schemes offer the possibility of improved damage tolerance.

Figure 1 .
Figure 1.Designs and dimensions in mm of a) the baseline stiffener and b) the parametric stiffener.

Figure 4 .
Figure 4. Stacking sequence for the skin

Figure 9 .
Figure 9. Strain energy with composite lay-up (blue line) and with material orientations (red line)

Figure 10 .
Figure 10.Normalized energy release rates as a function of crack length (a) comparison between baseline stiffener design (Fig. 1a) and selected modified stiffener (Fig 1b with b = 3 mm, c = 10 mm and d =6.25 mm), (b) Influence of parameter b on G, (c) Influence of parameter c on G and (d) Influence of parameter d on G.

Figure 11 .
Figure 11.(a) Front view of failed specimen; (b) Exploded view showing the failed area; (c) Front view of bottom part showing 0 0 plies ; (d) Bottom view of the upper part showing delaminated 45 0 plies

Figure 12 .Figure 13 .
Figure 12.The parts of the FE model

Figure 14 .
Figure 14.Illustration of (a) Traction-separation law for cohesive zone models (b) Modified law to implement in FEM FM300 (adhesive)

Figure 15 .
Figure 15.Damage growth pattern compared with energy release rate predictions for tapered specimen

Figure 16 .
Figure 16.Comparing the results of the parametric study with the VCCT method (a) along the crack and (b) along the width of the stiffener.

Figure 17 .
Figure 17.The refined model that used in the VCCT method

Figure 19 .
Figure 19.a) Tapered stiffener after testing, b) FE model showing delamination path, and c) FE model of a specimen with boundary conditions.

Figure 20 .
Figure 20.Normalized strain energy release rates as a function of crack length; comparison between Baseline stiffener design, Tapered stiffener and Compliant stiffener with b=3 mm, c=10 mm and d=6.25 mm.

Figure 21 .
Figure 21.Normalized GT/Gc across the crack tip for crack a = 1 mm for the Baseline, Tapered and Compliant specimens.

Figure 22 .Figure 23 .
Figure 22.Loads and Peak frequencies versus displacement for a) the Baseline b) the Tapered and c) the Compliant stiffeners.A scale on the right hand side indicates the mode of failure typically associated with these peak frequencies [21].

Table 2 .
Stacking sequence for the skin and the stiffener

Table 3 ,
the solution is converged.Since the running time for the intermediate mesh was 35 minutes and the results appeared accurate and acceptable, this specific density of elements was selected for the rest of this study.

Table 3 .
Mesh Sensitivity Study results.

Table 4 .
Cohesive interaction properties

Table 6 .
The Imperial Data Acquisition (IDA) program was used to record load and displacement during the tests.Failure loads for the different specimen types, as well as the predicted failure loads using Eq. 2.