Cross-Influence between Intra-Laminar Damages and Fibre Bridging at the Skin–Stringer Interface in Stiffened Composite Panels under Compression

In this paper, the skin–stringer separation phenomenon that occurs in stiffened composite panels under compression is numerically studied. Since the mode I fracture toughness and, consequently, the skin–stringer separation can be influenced by the fibre bridging phenomenon at the skin–stringer interface, in this study, comparisons among three different material systems with different fibre bridging sensitivities have been carried out. Indeed, a reference material system has been compared, in terms of toughness performance, against two materials with different degrees of sensitivity to fibre bridging. A robust numerical procedure for the delamination assessment has been used to mimic the skin–stringer separation. When analysing the global compressive behaviour of the stiffened panel, intra-laminar damages have been considered in conjunction with skin–stringer debonding to evaluate the effect of the fibre and matrix breakage on the separation between the skin and the stringer for the three analysed material systems. The latter are characterised by different toughness characteristics and fibre bridging sensitivities, resulting in a different material toughness.


Introduction
Layered composite materials are usually adopted to produce lightweight structures, able to fulfil specific stiffness and strength requirements, by opportunely stacking a sequence of layers with specific fibre orientations. Therefore, the material performance can be optimised by choosing appropriate fibre orientations according to loading and boundary conditions. However, the usage of composite structures is currently limited by a lack of knowledge of their damage behaviour, especially under compressive loading conditions [1][2][3]. In the case of Carbon-Fibre-Reinforced Polymer (CFRP) laminates, damages are generally triggered by inter-laminar fractures, such as delamination or debonding, that interact with intra-laminar failure mechanisms, such as matrix and fibre breakage. The delaminations in composite material structures have been extensively studied in the literature, both experimentally and numerically. In particular, experimental studies on the composites' inter-laminar fracture behaviour [4][5][6][7] are mandatory for the development of robust numerical procedures for damage evolution simulation. In [8], three different numerical techniques for the simulation of inter-laminar damage growth in Double-Cantilever Beam (DCB) specimens are proposed. The first two models are based on the Cohesive Zone Model (CZM), while the other model is based on the Virtual Crack Closure Technique (VCCT). Finally, the results have been compared against data from experimental tests. A numerical method, based on the VCCT, and an experimental study are presented also in [9] to examine the fracture behaviour of composite double-cantilever beam specimens. Finite element models for the investigation

VCCT-Based Numerical Procedure (SMXB)
In the following paragraphs, the theory behind the numerical models adopted in this work is described. The SMXB procedure, introduced in [13], for the simulation of the inter-laminar damage's evolution is briefly explained. It underlies the capability to evaluate crack propagation irrespective of the mesh and time-step size and to take into account the fibre bridging phenomenon. Finally, the adopted Hashin's criteria are introduced.

SMXB Methodology
Using the Ansys Parametric Design Languages (APDL), a new numerical methodology has been implemented in the Ansys ® Finite Element Model (FEM) simulation software. SMXB is described in detail in [12,13]. It is the result of an improved combination of the Virtual Crack Closure Technique (VCCT) and the Fail release (FR) method. In particular, the proposed methodology is able to overcome the main limitations of the standard VCCT-FR approaches [12] in terms of dependence on mesh and load-step size. The potential area of interest with the delamination growth is simulated by means of contact elements and fail-release conditions, which can be deactivated where the propagation criterion is fulfilled. The criterion used in SMXB is the Linear Power Law, which is described by Equation (1).
The term G j (with j = I, II, III) is the energy release rate, related to the failure mode j, while G jc is the corresponding critical value.
The delamination growth simulation strongly depends on the size of the elements. The delaminated area can generally be under-or overestimated. Through the combination of three separate and interacting moduli, the SMXB technique is able to change the local coordinate system related to each node, belonging to the delamination front, to properly determine the virtually released area. This allows it to equate the numerically evaluated delaminated area with the area that should be released to reach E d = 1. Furthermore, the peaks in the strain energy release rate that occur at corner nodes, due to the change in the delamination front's shape, are avoided by estimating the energy release rate on the delamination front's segments, rather than nodes. In Figure 1, a flowchart representing the SMXB workflow is shown. The detailed description of SMXB features is reported in [12,13].
Materials 2019, 12, x FOR PEER REVIEW 3 of 17 detail in [12,13]. It is the result of an improved combination of the Virtual Crack Closure Technique (VCCT) and the Fail release (FR) method. In particular, the proposed methodology is able to overcome the main limitations of the standard VCCT-FR approaches [12] in terms of dependence on mesh and load-step size. The potential area of interest with the delamination growth is simulated by means of contact elements and fail-release conditions, which can be deactivated where the propagation criterion is fulfilled. The criterion used in SMXB is the Linear Power Law, which is described by Equation (1).
The term Gj (with j = I, II, III) is the energy release rate, related to the failure mode j, while Gjc is the corresponding critical value.
The delamination growth simulation strongly depends on the size of the elements. The delaminated area can generally be under-or overestimated. Through the combination of three separate and interacting moduli, the SMXB technique is able to change the local coordinate system related to each node, belonging to the delamination front, to properly determine the virtually released area. This allows it to equate the numerically evaluated delaminated area with the area that should be released to reach Ed = 1. Furthermore, the peaks in the strain energy release rate that occur at corner nodes, due to the change in the delamination front's shape, are avoided by estimating the energy release rate on the delamination front's segments, rather than nodes. In Figure 1, a flowchart representing the SMXB workflow is shown. The detailed description of SMXB features is reported in [12,13].

Fibre Bridging Module
During delamination evolution, fibre-reinforced composite materials can show crack bridging, such as fibres crossing the open crack. Such a phenomenon, namely fibre bridging, increases the fracture toughness of the material as the crack length grows (see Figure 2) and, consequently, it is able to delay the delamination evolution. In [13], the SMXB procedure has been improved by implementing a modulus that is able to mimic the fibre bridging phenomenon by considering the GIc variations induced by fibre bridging (R-Curve, Figure 2c). Essentially, the GIc value changes as the constraints between nodes, representing the delaminated area, are released, as shown in Figure 2. Therefore, the fracture toughness grows according to the R-curve until a steady state value, which is representative of the maximum extension of the fibre bridging, is reached.

Fibre Bridging Module
During delamination evolution, fibre-reinforced composite materials can show crack bridging, such as fibres crossing the open crack. Such a phenomenon, namely fibre bridging, increases the fracture toughness of the material as the crack length grows (see Figure 2) and, consequently, it is able to delay the delamination evolution. In [13], the SMXB procedure has been improved by implementing a modulus that is able to mimic the fibre bridging phenomenon by considering the G Ic variations induced by fibre bridging (R-Curve, Figure 2c). Essentially, the G Ic value changes as the constraints between nodes, representing the delaminated area, are released, as shown in Figure 2. Therefore, the fracture toughness grows according to the R-curve until a steady state value, which is representative of the maximum extension of the fibre bridging, is reached. According to the theory underlying the fibre bridging modulus, the formula in Equation (1) should be written (where a is the crack extension) as:

Hashin Criteria and Material Property Degradation Rules
Hashin and Rotem in [23,24], on the basis of experimental evidence, propose for unidirectional composite materials the presence of two different intra-laminar mechanisms of failure that are dominated by the fibre phase and the matrix phase. A criterion for each of these two distinct failure modes has been introduced. In particular, under compressive load in the fibres' direction, a simple criterion of maximum stress is proposed, as shown in Equation (3).
In the case of tensile fibre failure, the expression in Equation (4) can be used.
In the case of transverse stress, the expression for predicting matrix failure under tensile stress is described by Equation (5), and, if the stress is compressive, Equation (6) should be considered.
According to Equations (3)-(6), σ1 and σ2 are the stresses in the fibre direction and normal to the fibre direction, respectively; and Xc, Yc, Xt, Yt and S are the fibre compressive, matrix compressive, fibre tensile, matrix tensile, and shear strengths, respectively.
The Damage Evolution Law, which defines the way the material degrades, is shown in Figure 3.  According to the theory underlying the fibre bridging modulus, the formula in Equation (1) should be written (where a is the crack extension) as:

Hashin Criteria and Material Property Degradation Rules
Hashin and Rotem in [23,24], on the basis of experimental evidence, propose for unidirectional composite materials the presence of two different intra-laminar mechanisms of failure that are dominated by the fibre phase and the matrix phase. A criterion for each of these two distinct failure modes has been introduced. In particular, under compressive load in the fibres' direction, a simple criterion of maximum stress is proposed, as shown in Equation (3).
In the case of tensile fibre failure, the expression in Equation (4) can be used.
In the case of transverse stress, the expression for predicting matrix failure under tensile stress is described by Equation (5), and, if the stress is compressive, Equation (6) should be considered.
According to Equations (3)-(6), σ 1 and σ 2 are the stresses in the fibre direction and normal to the fibre direction, respectively; and X c , Y c , X t , Y t and S are the fibre compressive, matrix compressive, fibre tensile, matrix tensile, and shear strengths, respectively.
The Damage Evolution Law, which defines the way the material degrades, is shown in Figure 3. Location A in Figure 3 represents Hashin's stress limit. When this location is reached, the material stiffness is reduced by means of degradation factor d. Following the onset of damage, a reduction in material stiffness occurs immediately. The constitutive relationship for a damaged material is given by Equation (7), where σ is the nominal Cauchy stress, ε is the total elastic strain and [D] d is the damaged elasticity matrix.  Location A in Figure 3 represents Hashin's stress limit. When this location is reached, the material stiffness is reduced by means of degradation factor d. Following the onset of damage, a reduction in material stiffness occurs immediately. The constitutive relationship for a damaged material is given by Equation (7), where σ is the nominal Cauchy stress, ε is the total elastic strain and [D]d is the damaged elasticity matrix.
For a transversely isotropic material in a plane stress state, the damaged elasticity matrix can be expressed as in Equation (8), where ( )( The terms ds, dm and df represent the shear, matrix and fibre damage variables, respectively, that can be calculated as in Equation (9), where λf + , λf -, λm + and λmare the fibre tension, fibre compression, matrix tension and matrix compression failure, respectively, evaluated from the effective stress σ. Two damage evolution approaches are available: instantaneous degradation and gradual degradation. The first approach, which has been considered in this work, instantaneously reduces the material stiffness based on damage variables. According to the second approach, the damage variables progressively grow.

Numerical Model
An I-stiffened reinforced composite panel, with induced delamination between the panel and the stringer, has been considered in this work. The geometrical description is shown in Figure 4. Twenty solid layered node elements have been considered to discretize the reinforcement and the skin underneath, while linear layered structural shell elements have been used to mesh the remaining skin panel, as shown in Figure 5, in order to reduce the computational cost [25].  For a transversely isotropic material in a plane stress state, the damaged elasticity matrix can be expressed as in Equation (8), The terms d s , d m and d f represent the shear, matrix and fibre damage variables, respectively, that can be calculated as in Equation (9), where λ f + , λ f − , λ m + and λ m − are the fibre tension, fibre compression, matrix tension and matrix compression failure, respectively, evaluated from the effective stress σ. Two damage evolution approaches are available: instantaneous degradation and gradual degradation. The first approach, which has been considered in this work, instantaneously reduces the material stiffness based on damage variables. According to the second approach, the damage variables progressively grow.

Numerical Model
An I-stiffened reinforced composite panel, with induced delamination between the panel and the stringer, has been considered in this work. The geometrical description is shown in Figure 4. Twenty solid layered node elements have been considered to discretize the reinforcement and the skin underneath, while linear layered structural shell elements have been used to mesh the remaining skin panel, as shown in Figure 5, in order to reduce the computational cost [25]. The model consists of 2730 solid elements and 840 shell elements. The overall number of nodes is 22,837.
In the framework of the performed numerical analyses, the panels' borders have been bounded, taking into account the effect of the potting regions in order to assure an even distribution of the load. The panel and the reinforcement have been linked, in the area where delamination is expected to grow, by using contact elements with "birth and death" features. Furthermore, node-to-surface contact elements have been positioned in the delaminated area to prevent penetration. Compressive displacements have been imposed on one side of the panel, and the other side has been clamped. Non-linear analyses have been performed, with variable time-steps, according to the SMXB procedure. Indeed, the numerical tool iteratively changes the time-step's size until E d = 1 is reached. Furthermore, no warping phenomena are expected to occur due to the symmetrical and balanced stacking sequence and the boundary conditions.   In the framework of the performed numerical analyses, the panels' borders have been bounded, taking into account the effect of the potting regions in order to assure an even distribution of the load. The panel and the reinforcement have been linked, in the area where delamination is expected to grow, by using contact elements with "birth and death" features. Furthermore, node-to-surface contact elements have been positioned in the delaminated area to prevent penetration. Compressive displacements have been imposed on one side of the panel, and the other side has been clamped. Non-linear analyses have been performed, with variable time-steps, according to the SMXB procedure. Indeed, the numerical tool iteratively changes the time-step's size until Ed = 1 is reached. Furthermore, no warping phenomena are expected to occur due to the symmetrical and balanced stacking sequence and the boundary conditions.
The mechanical properties of the used fibre-reinforced composite materials are listed in Table 1. Three different material models, with different GIc values, have been considered: • a reference material with GIc = 0.243 kJ/m 2 , which was obtained by means of the standard Solvay curing process; • a toughened material with GIc = 0.456 kJ/m 2 , which was obtained by means of a modified curing process aimed at improving the fracture toughness; and • a high fibre bridging material, with an increasing value of GIc, which was obtained by means of a modified curing process and characterised by high sensitivity to the fibre bridging phenomenon.    In the framework of the performed numerical analyses, the panels' borders have been bounded, taking into account the effect of the potting regions in order to assure an even distribution of the load. The panel and the reinforcement have been linked, in the area where delamination is expected to grow, by using contact elements with "birth and death" features. Furthermore, node-to-surface contact elements have been positioned in the delaminated area to prevent penetration. Compressive displacements have been imposed on one side of the panel, and the other side has been clamped. Non-linear analyses have been performed, with variable time-steps, according to the SMXB procedure. Indeed, the numerical tool iteratively changes the time-step's size until Ed = 1 is reached. Furthermore, no warping phenomena are expected to occur due to the symmetrical and balanced stacking sequence and the boundary conditions.
The mechanical properties of the used fibre-reinforced composite materials are listed in Table 1. Three different material models, with different GIc values, have been considered: • a reference material with GIc = 0.243 kJ/m 2 , which was obtained by means of the standard Solvay curing process; • a toughened material with GIc = 0.456 kJ/m 2 , which was obtained by means of a modified curing process aimed at improving the fracture toughness; and • a high fibre bridging material, with an increasing value of GIc, which was obtained by means of a modified curing process and characterised by high sensitivity to the fibre bridging phenomenon. The mechanical properties of the used fibre-reinforced composite materials are listed in Table 1. Three different material models, with different G Ic values, have been considered: • a reference material with G Ic = 0.243 kJ/m 2 , which was obtained by means of the standard Solvay curing process; • a toughened material with G Ic = 0.456 kJ/m 2 , which was obtained by means of a modified curing process aimed at improving the fracture toughness; and • a high fibre bridging material, with an increasing value of G Ic , which was obtained by means of a modified curing process and characterised by high sensitivity to the fibre bridging phenomenon.
The resistance curve, characterising the material as highly sensitive to the fibre bridging, is shown in the chart in Figure 6, and is compared with the constant values of G Ic . The resistance curve, characterising the material as highly sensitive to the fibre bridging, is shown in the chart in Figure 6, and is compared with the constant values of GIc.

Results
In this section, the results of the proposed numerical analyses are presented and discussed for all the considered material systems. The adopted numerical approach, introduced in Section 3, has been extensively investigated and validated by means of comparisons with experimental test results. These studies can be found in [13,21]; however, they are not reported in this work for the sake of brevity.
The compressive behaviour of the I-stringer panels, considering only the inter-laminar damages, is initially presented. Afterward, the Hashin's criteria have been added to the numerical simulations,

Results
In this section, the results of the proposed numerical analyses are presented and discussed for all the considered material systems. The adopted numerical approach, introduced in Section 3, has been extensively investigated and validated by means of comparisons with experimental test results. These studies can be found in [13,21]; however, they are not reported in this work for the sake of brevity.
The compressive behaviour of the I-stringer panels, considering only the inter-laminar damages, is initially presented. Afterward, the Hashin's criteria have been added to the numerical simulations, and the effect of the fibre/matrix breakage on the debonding evolution has been studied.

Debonding Simulations: Material Sensitivity Analysis
The global compressive behaviour of the analysed configurations, considering only the inter-laminar damage evolution, is reported in Figure 7, where the load versus applied displacement curves are compared for the three different analysed materials. It has been found that the delamination onset is influenced by the material's characteristics. Indeed, in the reference material model, the maximum attained load is 76.6 kN as compared to the 85.4 kN of the material that is highly sensitive to the fibre bridging phenomenon. On the other hand, the toughened material model reaches a maximum load of 90 kN.
In Figure 8, the results in terms of delaminated area versus applied displacements are introduced for the three analysed materials. Even if the fibre bridging phenomenon has been shown to delay the debonding evolution, the trend of the delaminated area does not change when considering the reference material and the material sensitive to the fibre bridging. Indeed, both are characterised by a stable separation evolution until a sudden collapse occurs, which corresponds to a complete loss of stiffness (maximum load in Figure 7). Debonding arises simultaneously in the reference material model and the configuration sensitive to the fibre bridging. In the toughened material configuration, skin-stringer separation onset is delayed and the debonding growth is more stable (as can be seen from the reduced curve slope) and seems to reach a steady state. This is because the starting value of G Ic is the same for the configurations that are characterised by a reference material system and a material system with high sensitivity to fibre bridging (0.243 kJ/m 2 according to the graph in the Figure 6), while the value is higher for the configuration characterised by the toughened material systems (the blue curve in the Figure 6), where there is a delay in delamination onset. Overall, it can be observed that the Mode I fracture toughness has a significant influence on the panel-reinforcement separation's progression. In particular, increasing the Mode I fracture toughness value stabilises the inter-laminar damage progression.  In Figure 8, the results in terms of delaminated area versus applied displacements are introduced for the three analysed materials. Even if the fibre bridging phenomenon has been shown to delay the debonding evolution, the trend of the delaminated area does not change when considering the reference material and the material sensitive to the fibre bridging. Indeed, both are characterised by a stable separation evolution until a sudden collapse occurs, which corresponds to a complete loss of stiffness (maximum load in Figure 7). Debonding arises simultaneously in the reference material model and the configuration sensitive to the fibre bridging. In the toughened material configuration, skin-stringer separation onset is delayed and the debonding growth is more stable (as can be seen from the reduced curve slope) and seems to reach a steady state. This is because the starting value of GIc is the same for the configurations that are characterised by a reference material system and a material system with high sensitivity to fibre bridging (0.243 kJ/m 2 according to the graph in the Figure 6), while the value is higher for the configuration characterised by the toughened material systems (the blue curve in the Figure 6), where there is a delay in delamination onset. Overall, it can be observed that the Mode I fracture toughness has a significant influence on the panel-reinforcement separation's progression. In particular, increasing the Mode I fracture toughness value stabilises the inter-laminar damage progression.   In Figure 8, the results in terms of delaminated area versus applied displacements are introduced for the three analysed materials. Even if the fibre bridging phenomenon has been shown to delay the debonding evolution, the trend of the delaminated area does not change when considering the reference material and the material sensitive to the fibre bridging. Indeed, both are characterised by a stable separation evolution until a sudden collapse occurs, which corresponds to a complete loss of stiffness (maximum load in Figure 7). Debonding arises simultaneously in the reference material model and the configuration sensitive to the fibre bridging. In the toughened material configuration, skin-stringer separation onset is delayed and the debonding growth is more stable (as can be seen from the reduced curve slope) and seems to reach a steady state. This is because the starting value of GIc is the same for the configurations that are characterised by a reference material system and a material system with high sensitivity to fibre bridging (0.243 kJ/m 2 according to the graph in the Figure 6), while the value is higher for the configuration characterised by the toughened material systems (the blue curve in the Figure 6), where there is a delay in delamination onset. Overall, it can be observed that the Mode I fracture toughness has a significant influence on the panel-reinforcement separation's progression. In particular, increasing the Mode I fracture toughness value stabilises the inter-laminar damage progression. The debonding evolution trends, at different applied compressive displacement values, are shown in Figure 9 for each analysed configuration. According to the figure, the debonding fronts grow with an inclination of about 45° with respect to the loading direction. The debonding evolution trends, at different applied compressive displacement values, are shown in Figure 9 for each analysed configuration. According to the figure, the debonding fronts grow with an inclination of about 45 • with respect to the loading direction.
In Figure 10, the load versus the out-of-plane displacement curves, representative of the buckling behaviour of the panels, are shown for each analysed configuration. Local buckling of the delaminated region can be appreciated at 30 kN. The delamination growth initiation load values, together with the buckling loads, are listed for each analysed configuration in Table 2. The delamination growth initiation load for the reference material is the same as that for the material with high sensitivity to fibre bridging. Indeed, for these two materials, the value of the critical energy release rate at the beginning of delamination growth is the same. In Figure 10, the load versus the out-of-plane displacement curves, representative of the buckling behaviour of the panels, are shown for each analysed configuration. Local buckling of the delaminated region can be appreciated at 30 kN. The delamination growth initiation load values, together with the buckling loads, are listed for each analysed configuration in Table 2. The delamination growth initiation load for the reference material is the same as that for the material with high sensitivity to fibre bridging. Indeed, for these two materials, the value of the critical energy release rate at the beginning of delamination growth is the same.

Cross-Influence between Intra-Laminar Damages and Fibre Bridging
The three different panel configurations (characterised by different material systems) have been further analysed by considering both the inter-laminar and the intra-laminar damage evolution. The influence of the fibre/matrix breakage on the delamination evolution can be observed in Figure 11, where the delaminated areas found with the SMXB procedure are compared with the delaminated areas found by using SMXB coupled with the Hashin's criteria, for all the considered material models. According to the charts in Figure 11, the debonding onset and evolution are the same as in the previously presented analyses up to the fibre/matrix damage onset. Indeed, the fibre's breakage arises along the delamination front, influencing the delamination evolution, as shown in Figure 12. The debonded areas that were found with the SMXB procedure are compared with the debonded areas that were found by using SMXB coupled with the Hashin criteria, at the same applied displacement value, which corresponds to the final state of each analysed material system found by using the SMXB approach coupled with the Hashin's criteria. In Figure 12, the fibre's failure is superimposed in red on the images reporting the delamination evolution status. According to Figure  12, it can be observed that the intra-laminar damages arising on the delamination front influence the delamination evolution, resulting in a reduction of the delaminated area. However, despite the reduction in inter-laminar damage, the intra-laminar damages critically reduce the structure's load carrying capability, as shown in Table 3 and Figures 13 and 14. In particular, in Table 3, the bucking load and applied displacement values, together with the intra-laminar and inter-laminar onset load and applied displacement values, are listed, for each analysed configuration. Figure 13 shows the load versus the out-of-plane displacement comparisons. According to Table 3 and Figure 13, the onset of local buckling in the delaminated region, which takes place before the onset of intra-laminar damage, is not affected by the intra-laminar damage. On the contrary, the appearance of the intralaminar damages modifies the global compressive behaviour of the panels, leading prematurely to failure before the complete development of the buckling phenomenon seen in the previously shown According to the charts in Figure 11, the debonding onset and evolution are the same as in the previously presented analyses up to the fibre/matrix damage onset. Indeed, the fibre's breakage arises along the delamination front, influencing the delamination evolution, as shown in Figure 12. The debonded areas that were found with the SMXB procedure are compared with the debonded areas that were found by using SMXB coupled with the Hashin criteria, at the same applied displacement value, which corresponds to the final state of each analysed material system found by using the SMXB approach coupled with the Hashin's criteria. In Figure 12, the fibre's failure is superimposed in red on the images reporting the delamination evolution status. According to Figure 12, it can be observed that the intra-laminar damages arising on the delamination front influence the delamination evolution, resulting in a reduction of the delaminated area. However, despite the reduction in inter-laminar damage, the intra-laminar damages critically reduce the structure's load carrying capability, as shown in Table 3 and Figures 13 and 14. In particular, in Table 3, the bucking load and applied displacement values, together with the intra-laminar and inter-laminar onset load and applied displacement values, are listed, for each analysed configuration. Figure 13 shows the load versus the out-of-plane displacement comparisons. According to Table 3 and Figure 13, the onset of local buckling in the delaminated region, which takes place before the onset of intra-laminar damage, is not affected by the intra-laminar damage. On the contrary, the appearance of the intra-laminar damages modifies the global compressive behaviour of the panels, leading prematurely to failure before the complete development of the buckling phenomenon seen in the previously shown SMXB analyses.    Moreover, the global compressive behaviour, considering the coupled intra-laminar and interlaminar damages, seems to be different both in the evolution of delaminated areas and in load trends, which are shown in Figure 14. As expected, the maximum load decreases when introducing the intralaminar damage, which anticipates a structural collapse with respect to the ideal case analysed with   Figure 13. Load vs. Out-of-plane displacement comparison.
Moreover, the global compressive behaviour, considering the coupled intra-laminar and interlaminar damages, seems to be different both in the evolution of delaminated areas and in load trends, which are shown in Figure 14. As expected, the maximum load decreases when introducing the intralaminar damage, which anticipates a structural collapse with respect to the ideal case analysed with the SMXB model. In Figure 15, the load versus displacement curves, related to the three analysed material systems with combined intra-laminar and inter-laminar damage evolution, are shown. From this figure, it is possible to appreciate that the Mode I fracture toughness values do not have a significant influence Moreover, the global compressive behaviour, considering the coupled intra-laminar and inter-laminar damages, seems to be different both in the evolution of delaminated areas and in load trends, which are shown in Figure 14. As expected, the maximum load decreases when introducing the intra-laminar damage, which anticipates a structural collapse with respect to the ideal case analysed with the SMXB model.
In Figure 15, the load versus displacement curves, related to the three analysed material systems with combined intra-laminar and inter-laminar damage evolution, are shown. From this figure, it is possible to appreciate that the Mode I fracture toughness values do not have a significant influence on the global behaviour of the analysed configurations under compression. However, the maximum attained load in the case of the reference material system remains lower than that of the other configurations, which have a higher Mode I fracture toughness.  The intra-laminar damage status, for the three analysed configurations, at the final states, are shown in Figures 16-18. Contour values over unity indicate a damaged region, while zero indicates an undamaged region (0, undamaged; 1, damaged; 2, completely damaged). The intra-laminar damage status, for the three analysed configurations, at the final states, are shown in Figures 16-18. Contour values over unity indicate a damaged region, while zero indicates an undamaged region (0, undamaged; 1, damaged; 2, completely damaged).   From Figures 16-18, as expected, no relevant differences between the three analysed material systems, in terms of intra-laminar damage distribution, can be noted. The intra-laminar damage caused a loss in the carrying capability of the panel, regardless of the inter-laminar damage onset load. Figure 19 shows the delaminated area versus the applied displacements for the three considered configurations. It is clear from the chart that an increase in Mode I toughness reduces the delaminated area even in the presence of intra-laminar damage. Indeed, as an example, considering the fixed applied displacement value of 0.55 mm, the delaminated area can be reduced by up to 16.5% by increasing the Mode I fracture toughness.  Figure 19 shows the delaminated area versus the applied displacements for the three considered configurations. It is clear from the chart that an increase in Mode I toughness reduces the delaminated area even in the presence of intra-laminar damage. Indeed, as an example, considering the fixed applied displacement value of 0.55 mm, the delaminated area can be reduced by up to 16.5% by increasing the Mode I fracture toughness.

Conclusions
In this paper, the influence of intra-laminar damages on the evolution of skin-stringer debonding in stiffened composite panels under compression has been investigated. Three different material systems have been analysed, taking into account the fibre bridging phenomenon. A validated robust numerical procedure, (SMXB), has been used for the numerical analyses. The proposed methodology is able to simulate the delamination evolution, in the presence of fibre bridging, with no mesh and time-step size dependency. First, a set of numerical simulations was carried out in which only inter-laminar damages were considered. It was found that, by increasing the Mode I fracture toughness, a significant decrease in delamination evolution can be obtained. In

Conclusions
In this paper, the influence of intra-laminar damages on the evolution of skin-stringer debonding in stiffened composite panels under compression has been investigated. Three different material systems have been analysed, taking into account the fibre bridging phenomenon. A validated robust numerical procedure, (SMXB), has been used for the numerical analyses. The proposed methodology is able to simulate the delamination evolution, in the presence of fibre bridging, with no mesh and time-step size dependency. First, a set of numerical simulations was carried out in which only inter-laminar damages were considered. It was found that, by increasing the Mode I fracture toughness, a significant decrease in delamination evolution can be obtained. In particular, the material systems with increased fracture toughness experience stable debonding growth and seem to reach a steady state condition. A second set of analyses was carried out by taking into account the evolution of intra-laminar damage in order to assess the influence of fibre/matrix breakage on the delamination evolution. The results show that the maximum achieved load is reduced and that the intra-laminar damages take place at the delamination front, influencing the delamination evolution and resulting in a reduction of the delaminated area. In conclusion, it was also observed that the value of the fracture toughness of the interface does not influence the fibre/matrix failure growth.