A Numerical Assessment on the Inﬂuences of Material Toughness on the Crashworthiness of a Composite Fuselage Barrel

: In the present work, a numerical study on the dynamic response of a composite fuselage barrel, in relation to crashworthiness, has been investigated. The aim of this work is to investigate the inﬂuence of the material fracture toughness on the capability of a composite fuselage barrel to tolerate an impact on a rigid surface. Three di ﬀ erent material conﬁgurations with di ﬀ erent intra-laminar fracture energy values were considered to take into account variations in material toughness. Indeed, the dynamic behaviour of the analysed fuselage barrel has been numerically simulated by means of the FE (Finite Element) code Abaqus / Explicit. The e ﬀ ects of intralaminar fracture energy variations on the impact deformation of the barrel has been evaluated comparing the numerical results in terms of displacements and damage evolution for the three analysed material conﬁgurations.


Introduction
In recent years, several steps forward in increasing civil aircraft safety levels have been taken in relation to crashworthiness during an emergency landing. In this context, crashworthiness has become of fundamental relevance for the design and certification of an aircraft. The main aim in designing crashworthy aircrafts is to ensure the ability to absorb as much energy, deformations and breakage as possible during an impact event, without compromising the occupants' safety. Hence, the structure, under deformations and/or breakage conditions, must preserve the living space of passengers and crew, allowing escape routes from the aircraft [1]. Indeed, the entire structure should be able to guarantee the transfer to passengers of acceleration which can be tolerated by the human body [2][3][4][5]. In the past, the assessment of the crashworthiness performance of aircrafts was verified exclusively by several costly full scale experimental testing campaigns, as in the case of the Airbus A320 [6,7], Boing 707 [8][9][10], Boing 737 [11,12] and NAMC YS [13,14]. According to the literature, the experimental drop tests have allowed to identify the structural components which are generally involved in the absorption of energy during an impact event: circumferential frames close to the vertical supports and fasteners where large structural deformations occurred [15,16].
Nowadays, the use of advanced numerical tools can contribute to reduce the number of expensive experimental testing campaigns. Numerical developments focused on crash simulations of civil aircraft fuselages. Due to the crash phenomenon and fuselage structure's inherent complexity, numerical models with increasing levels of accuracy were progressively introduced by several authors. Basic models, using concentrated masses positioned at different locations and beams were used to start

Theoretical Background: Material Damage Model
Among several criteria for the prediction of the initiation and evolution of intra-laminar damage for structures made of unidirectional fibre-reinforced composite materials, the Hashin criteria are probably the most adopted ones. These criteria, based on the separation of the failure modes by introducing separated equations, allow to clearly identify the damage of the matrix and of the fibre under tension or compression stress conditions. Indeed, according to Hashin failure criteria [44], in Equations (1)-(4), the limit values for the onset of the damage for the fibre traction (F ft ), fibre compression (F fc ), matrix traction (F mt ) and matrix compression (F mc ) are, respectively, introduced.
In Equations (1)-(4) σ 11 , σ 22 , σ 12 are the components of the effective stress tensor, respectively, along fibre direction, matrix direction and shear while X T , X C , Y T , Y C , S L and S T are, respectively, the fibre tensile, fibre compressive, matrix tensile, matrix compressive and shear strengths in longitudinal and transversal directions. The evolution of the damage for each failure mode is modelled according to the bilinear law schematically represented in Figure 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW  3 of 16 probably the most adopted ones. These criteria, based on the separation of the failure modes by introducing separated equations, allow to clearly identify the damage of the matrix and of the fibre under tension or compression stress conditions. Indeed, according to Hashin failure criteria [44], in Equations (1-4), the limit values for the onset of the damage for the fibre traction (Fft), fibre compression (Ffc), matrix traction (Fmt) and matrix compression (Fmc) are, respectively, introduced.
In Equations (1)(2)(3)(4) , , are the components of the effective stress tensor, respectively, along fibre direction, matrix direction and shear while XT, XC, YT, YC, SL and ST are, respectively, the fibre tensile, fibre compressive, matrix tensile, matrix compressive and shear strengths in longitudinal and transversal directions. The evolution of the damage for each failure mode is modelled according to the bilinear law schematically represented in Figure 1. According to Figure 1, it is possible to identify the onset and evolution damage phases for each damage mode. In particular, point A of Figure 1 identifies the condition of Hashin's criterion satisfaction.
At this condition, the element has absorbed an amount of energy represented by the yellow triangle area. Beyond this condition, the evolution of the damage and the corresponding loss of stiffness of the material is simulated by introducing a parameter allowing to reduce linearly the characteristics of the undamaged material. The phase of the damage evolution is represented by the segment in Figure 1 connecting point A to point C representing the completely damaged element condition. Actually, the element is considered completely damaged when it has absorbed an energy represented by the sum of the yellow and blue triangles areas. Indeed, this global area is representative of the fracture toughness of simulated material system.
In order to evaluate the degradation status of the element, a material degradation coefficient di is introduced which is evaluated independently for matrix and fibre under traction and compression conditions. The degradation coefficient di definition is introduced in Equation (5).  According to Figure 1, it is possible to identify the onset and evolution damage phases for each damage mode. In particular, point A of Figure 1 identifies the condition of Hashin's criterion satisfaction.
At this condition, the element has absorbed an amount of energy represented by the yellow triangle area. Beyond this condition, the evolution of the damage and the corresponding loss of stiffness of the material is simulated by introducing a parameter allowing to reduce linearly the characteristics of the undamaged material. The phase of the damage evolution is represented by the segment in Figure 1 connecting point A to point C representing the completely damaged element condition. Actually, the element is considered completely damaged when it has absorbed an energy represented by the sum of the yellow and blue triangles areas. Indeed, this global area is representative of the fracture toughness of simulated material system.
In order to evaluate the degradation status of the element, a material degradation coefficient di is introduced which is evaluated independently for matrix and fibre under traction and compression conditions. The degradation coefficient di definition is introduced in Equation (5).
where δ i,eq is equivalent displacement; δ 0 i,eq is equivalent displacement on set damage; δ t i,eq is equivalent displacement on full damage (Point C in Figure 1). Moreover, f c , f t , m c , m t are fibre failure compression, fibre failure tension, matrix failure compression and matrix failure tension.
In Equation (6), the definition of the maximum equivalent displacement, reached at point C of Figure 1, is introduced.
In Equation (6), σ 0 i,eq and δ 0 i,eq are, respectively, the equivalent stress and displacement at the Hashin's limit condition (point A), and G c is the material fracture toughness, i.e., the area of the global triangle in Figure 1. In Table 1, the relations adopted to evaluate the equivalent stress and displacement are introduced.

12
Fibre Compression In Table 1, L C and < > are, respectively, the element characteristic length and the Macauley bracket operator [45].

Fuselage Barrel Test Case: Geometrical Description and Finite Element Model
As already mentioned, in this paper, a composite fuselage barrel has been adopted as a numerical test case to assess the influence of the fracture toughness characteristics on the mechanical behaviour of a composite structure undergoing a crash event. In this paragraph, the numerical test case is introduced by providing details on the geometry, materials, boundary conditions and on the finite element discretisation.

Geometrical Description
The geometry of the fuselage section adopted as a test case for the impact numerical analyses is shown in Figure 2. The considered fuselage barrel does not take into account geometric variations related to the tail plane, nose and wing attachments. This assumption has allowed to neglect the stress concentration effects at a global level while stress concentration effects were considered at sub-components level.
The fuselage is composed of different structural sub-components (shown with different colours in Figure 2). The overall dimensions of the fuselage section are reported in Figure 2, together with an indication on the impact angle between the fuselage and the impacted rigid plane, which has been set to 3 degrees. Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 16  Figure 3 allows to highlight the material systems adopted for the different subcomponents of the fuselage. Three different material systems were adopted: namely a fibre-reinforced composite material system with long unidirectional fibres (UNIDIRECTIONAL UD), a woven fabric material system and AL2025 aluminium. The related material mechanical properties are introduced in the following Tables 2-3-4.    Figure 3 allows to highlight the material systems adopted for the different subcomponents of the fuselage. Three different material systems were adopted: namely a fibre-reinforced composite material system with long unidirectional fibres (UNIDIRECTIONAL UD), a woven fabric material system and AL2025 aluminium. The related material mechanical properties are introduced in the following Tables 2-4.  Figure 3 allows to highlight the material systems adopted for the different subcomponents of the fuselage. Three different material systems were adopted: namely a fibre-reinforced composite material system with long unidirectional fibres (UNIDIRECTIONAL UD), a woven fabric material system and AL2025 aluminium. The related material mechanical properties are introduced in the following Tables 2-3-4.    For the aluminium plates, a thickness of 8 mm was chosen. Table 5 shows the stacking sequences for the composites' made sub-components. A discrete coordinate system has been introduced in the numerical model for each composite barrel sub-component. Fibres are oriented according to the fuselage longitudinal direction, while the normal axis has been chosen according to the element normal direction.

Material System
As already mentioned, three different material systems configurations, characterised by different fracture toughness energies, were considered for the simulations in order to assess the influence of fracture characteristics on the evolution of the damage and hence on the dynamic response of the whole fuselage. Each configuration, identified with a number (I,II,III), is associated to different values of the fracture toughness energies for the matrix and fibre under tensile and compression conditions for both the UNIDIRECTIONAL UD and the woven composite material systems. Moreover, tensile and compression fracture energy values for the woven fabric were assumed identical for all the considered configurations.
As it can be observed from Table 6, Configuration II can be considered the most toughened configuration while configurations I is characterised by the lowest toughness. Finally, configuration III is characterised by an intermediate/low toughness. Moreover, for all the material configurations, the degradation of the matrix toughness in both tension and compression was not taken into account. In fact, the small matrix toughness variations during propagation have a significant effect only on inter-laminar damages which have not been considered in the frame of the presented application.

Finite Element Description
The finite element model of the fuselage section is introduced in Figure 4A,B. The entire model consists of 1,976,157 nodes and 995,858 elements. The elements used for the aluminium components are discretised with three-dimensional elements with eight nodes and reduced integration scheme available in the Abaqus library (SC8R); the elements have a dimension of about 10 mm × 10 mm. Planar shell elements with rigid body constraint were adopted to model the impact plate between representing the impacted rigid ground. With reference to the computational grid of the composite components of the fuselage section, eight nodes continuum shell elements with a reduced integration scheme were adopted. These elements are general-purpose shells allowing finite membrane deformation and large rotations, and thus, they are suitable for nonlinear geometric analysis. These elements include the effects of transverse shear deformations and the effects of thickness change. Connections between the sub-components were simulated with tie-constraints allowing to couple a master surface with a slave one by taking into account with all the degrees of freedom. The connections between the vertical stanchions and the passenger floor/frames were simulated by fastener elements allowing to introduce breakage criteria based on maximum force in each reference direction.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 16 Connections between the sub-components were simulated with tie-constraints allowing to couple a master surface with a slave one by taking into account with all the degrees of freedom. The connections between the vertical stanchions and the passenger floor/frames were simulated by fastener elements allowing to introduce breakage criteria based on maximum force in each reference direction.

Boundary Conditions
The fuselage section, considered in the frame of the numerical simulations, has a mass of approximately 940 kg with several mass points added to simulate external non-structural components: row of seats and occupants and balancing mass. The balancing mass has been opportunely chosen in order to create an unbalancing effect between the right and the left side of the fuselage barrel during the simulated drop test. This unbalancing effect has been introduced in order to assess the effects of the variations in material toughness on lateral damage distribution during the impact simulation. A general contact type "all-with-self" available in Abaqus/Expicit was adopted with a friction coefficient of 0.3.
In addition, the weight of the three seats and the two dummies was applied in the centre of gravity of the seats and rigidly connected in the points of attachment of the seats to the passenger floor. Figure 5 shows the points of connection of the various masses.
approximately 940 kg with several mass points added to simulate external non-structural components: row of seats and occupants and balancing mass. The balancing mass has been opportunely chosen in order to create an unbalancing effect between the right and the left side of the fuselage barrel during the simulated drop test. This unbalancing effect has been introduced in order to assess the effects of the variations in material toughness on lateral damage distribution during the impact simulation. A general contact type "all-with-self" available in Abaqus/Expicit was adopted with a friction coefficient of 0.3.
In addition, the weight of the three seats and the two dummies was applied in the centre of gravity of the seats and rigidly connected in the points of attachment of the seats to the passenger floor. Figure 5 shows the points of connection of the various masses. The impact simulation of the fuselage section with the ground has been carried out by considering a rigid plane bounded in space and by applying an initial velocity to the entire fuselage section equal to 9900 mm/s. The initial velocity has been evaluated by considering the fuselage section dropping from a height of 5000 mm with respect to the rigid ground surface.

Numerical Results
In this section, the numerical results obtained for the three analysed configurations according to Table 6 are introduced. The different configurations were compared in terms of maximum displacement along the drop direction and in terms of qualitative deformation of the whole fuselage section.
For the three different configurations, the history of the vertical displacement of the control points A and B, as identified in Figure 6, was evaluated. The control points are positioned on the crossbeams of the passenger floor in the middle of the beams where the seat rails are fixed, as shown in Figure 6. The impact simulation of the fuselage section with the ground has been carried out by considering a rigid plane bounded in space and by applying an initial velocity to the entire fuselage section equal to 9900 mm/s. The initial velocity has been evaluated by considering the fuselage section dropping from a height of 5000 mm with respect to the rigid ground surface.

Numerical Results
In this section, the numerical results obtained for the three analysed configurations according to Table 6 are introduced. The different configurations were compared in terms of maximum displacement along the drop direction and in terms of qualitative deformation of the whole fuselage section.
For the three different configurations, the history of the vertical displacement of the control points A and B, as identified in Figure 6, was evaluated. The control points are positioned on the crossbeams of the passenger floor in the middle of the beams where the seat rails are fixed, as shown in Figure 6.  Figure 7A shows the vertical displacements of the control points A and B obtained for the configuration I which is the configuration with the lowest toughness. As shown in Figure 7A, for this configuration, the maximum value of the vertical displacement at control point was reached (about 250 mm). Moreover, a discrepancy in the displacements between the two points can be noticed. This is generated by the aforementioned unbalancing effect of the balancing mass between the left and the right side of the fuselage. However, the configuration with the lowest toughness is able to ensure a  Figure 7A shows the vertical displacements of the control points A and B obtained for the configuration I which is the configuration with the lowest toughness. As shown in Figure 7A, for this configuration, the maximum value of the vertical displacement at control point was reached (about 250 mm). Moreover, a discrepancy in the displacements between the two points can be noticed. This is generated by the aforementioned unbalancing effect of the balancing mass between the left and the right side of the fuselage. However, the configuration with the lowest toughness is able to ensure a quasi-symmetric distribution of the damage, and deformations between the left and the right side of the fuselage are small despite of the unbalancing effect of the introduced mass. Figure 7B-D shows the deformed structure at the maximum point of deflection. In particular, Figure 7B shows the frontal view of the deformed fuselage section. Figure 7C shows the isometric view, while Figure 7D shows a zoomed view of the deformed and broken structure in the cargo area. Indeed, the cargo area can be considered the most stressed one during a crash event. This consideration leads to pay special attention to this area in a crashworthy design in order to maximise the dissipation of the kinetic energy.  Figure 7A shows the vertical displacements of the control points A and B obtained for the configuration I which is the configuration with the lowest toughness. As shown in Figure 7A, for this configuration, the maximum value of the vertical displacement at control point was reached (about 250 mm). Moreover, a discrepancy in the displacements between the two points can be noticed. This is generated by the aforementioned unbalancing effect of the balancing mass between the left and the right side of the fuselage. However, the configuration with the lowest toughness is able to ensure a quasi-symmetric distribution of the damage, and deformations between the left and the right side of the fuselage are small despite of the unbalancing effect of the introduced mass. Figure 7B-D shows the deformed structure at the maximum point of deflection. In particular, Figure 7B shows the frontal view of the deformed fuselage section. Figure 7C shows the isometric view, while Figure 7D shows a zoomed view of the deformed and broken structure in the cargo area. Indeed, the cargo area can be considered the most stressed one during a crash event. This consideration leads to pay special attention to this area in a crashworthy design in order to maximise the dissipation of the kinetic energy.   Figure 8A shows the vertical displacements of the control points A and B obtained for the configuration II which is the configuration with the highest toughness. As shown in Figure 8A, for this configuration, the minimum value of the vertical displacement at the control point is reached (about 220 mm). Moreover, the maximum discrepancy in the displacements between the two points (A and B) can be noticed. Indeed, the configuration with the highest toughness amplifies the asymmetry in damage distribution and deformations between the left and the right side of the fuselage due to the unbalancing effect of the introduced mass. Figure 8B-D shows the deformed structure at the maximum point of deflection. In particular, Figure 8B shows the frontal view of the deformed fuselage section. Figure 8C shows the isometric view, while Figure 8D shows a zoomed view of the deformed and broken structure in the cargo area. Indeed, due to the toughening effect, the damage seems to be less distributed in the whole fuselage and more concentrated in the cargo area if compared to configuration I ( Figure 7D). Figure 8.B-D shows the deformed structure at the maximum point of deflection. In particular, Figure 8.B shows the frontal view of the deformed fuselage section. Figure 8C shows the isometric view, while Figure 8.D shows a zoomed view of the deformed and broken structure in the cargo area. Indeed, due to the toughening effect, the damage seems to be less distributed in the whole fuselage and more concentrated in the cargo area if compared to configuration I (Figure 7.D).  Figure 9.A shows the vertical displacements of the control points A and B obtained for the configuration III, which is the configuration with the intermediate low toughness. As shown in Figure  9.A, for this configuration, as expected, the minimum value of the vertical displacement at control point is almost identical to the one observed for configuration I (about 250 mm). The same considerations can be repeated for the magnitude of the discrepancy in the displacements between the two points (A and B) due to the unbalancing effect of the introduced mass. Figure 9.B-D highlights the damaged structure with particular focus on the maximum deflection point. Moreover, in Figure 9.B, a frontal view of the deformed fuselage section is reported. In Figure  9.C, an isometric view of the considered fuselage section is exhibited, while a zoomed view of the damaged cargo area of the structure is shown in Figure 9.D. Indeed, for configuration III, the cargo area seems to have a reduced damaged concentration if compared with configuration II.  Figure 9A shows the vertical displacements of the control points A and B obtained for the configuration III, which is the configuration with the intermediate low toughness. As shown in Figure 9A, for this configuration, as expected, the minimum value of the vertical displacement at control point is almost identical to the one observed for configuration I (about 250 mm). The same considerations can be repeated for the magnitude of the discrepancy in the displacements between the two points (A and B) due to the unbalancing effect of the introduced mass.

Configuration Comparison
In order to appreciate the effect of the fracture toughness on the mechanical behaviour of the fuselage undergoing the crash event, the displacement patterns were superimposed in Figure 10. Indeed, in Figure 10.A and Figure 10.B, respectively, the vertical displacements of control point A and point B obtained for the three analysed configurations were superimposed. As shown in Figure  10A,B, for all the configurations, the slopes at the beginning of the crash event are identical straight lines. In fact, at the beginning of the crash event, although the structure deforms and the damage has   Figure 9B, a frontal view of the deformed fuselage section is reported. In Figure 9C, an isometric view of the considered fuselage section is exhibited, while a zoomed view of the damaged cargo area of the structure is shown in Figure 9D. Indeed, for configuration III, the cargo area seems to have a reduced damaged concentration if compared with configuration II.

Configuration Comparison
In order to appreciate the effect of the fracture toughness on the mechanical behaviour of the fuselage undergoing the crash event, the displacement patterns were superimposed in Figure 10. Indeed, in Figures 10A and 10B, respectively, the vertical displacements of control point A and point B obtained for the three analysed configurations were superimposed. As shown in Figure 10A,B, for all the configurations, the slopes at the beginning of the crash event are identical straight lines. In fact, at the beginning of the crash event, although the structure deforms and the damage has started, the effects of the damage evolution were found almost negligible. Indeed, the effects of the in-plane fracture toughness were found more significant later during the crash event almost at the maximum deflection.
Configurations I and III, characterised by a similar value of fracture energies for both materials, both show a maximum deflection of about 250 mm for control point A and about 230 mm for control point B. For both these configurations, a significant variation of stiffness has been found at 0.022 s for control point A ( Figure 10A) and 0.027 s for control point B ( Figure 10B). On the other hand, configuration II, characterised by the highest toughness values, shows a maximum deformation of about 220 and 180 mm, respectively, for control points A and B. These maximum values, as expected, are lower if compared to configurations I and III.
Finally, as already remarked, the configuration with the highest toughness, differently from configurations I and III, amplifies the asymmetrical distribution of the damage between the left side (control point A) and the right side (control point B) of the fuselage induced by the unbalancing effect of the introduced mass. Moreover, in Figure 11.A, a comparison of damage energy dissipation is reported for the three fracture toughness energies. Configuration II dissipates a lower amount of energy if compared to the other configurations. Configurations I and III dissipate a relevant amount of the total energy as damage energy; consequently, a higher value of the maximum displacement was found for this configurations, compared to configuration II. Moreover, in Figure 11, a comparison of damage energy dissipation is reported for the three fracture toughness energies. Configuration II dissipates a lower amount of energy if compared to the other configurations. Configurations I and III dissipate a relevant amount of the total energy as damage energy; consequently, a higher value of the maximum displacement was found for this configurations, compared to configuration II. Moreover, in Figure 11.A, a comparison of damage energy dissipation is reported for the three fracture toughness energies. Configuration II dissipates a lower amount of energy if compared to the other configurations. Configurations I and III dissipate a relevant amount of the total energy as damage energy; consequently, a higher value of the maximum displacement was found for this configurations, compared to configuration II.

Conclusions
In the presented work, a numerical study on the influence of the in-plane toughness on the dynamic behaviour of a complex composite fuselage barrel, undergoing a crash event, was attempted. In order to evaluate the effects of intralaminar fracture toughness, three different material systems characterised by different toughness (low toughness, high toughness and intermediate/low toughness) were considered, both for the unidirectional fibre-reinforced composite and for the woven fabric one.
As a result of the performed numerical study, a relevant influence of the in-plane toughness on the global dynamic response of the fuselage barrel was found. Actually, the configuration characterised by the highest toughness showed the lowest vertical deflection and the most significant damage in the cargo area. Moreover, this configuration seems to be more sensitive to the unbalancing of the mass in the lateral direction producing the most significant asymmetry in damage distributions between the left and the right side of the fuselage. On the other hand, the configurations characterised by the lowest in-plane toughness showed the maximum vertical deflection and a more distributed damage evolution in the whole structure, leading to less significant damage in the cargo area. Finally, a low sensitivity to the unbalancing of the mass in the later direction was found.

Conclusions
In the presented work, a numerical study on the influence of the in-plane toughness on the dynamic behaviour of a complex composite fuselage barrel, undergoing a crash event, was attempted. In order to evaluate the effects of intralaminar fracture toughness, three different material systems characterised by different toughness (low toughness, high toughness and intermediate/low toughness) were considered, both for the unidirectional fibre-reinforced composite and for the woven fabric one.
As a result of the performed numerical study, a relevant influence of the in-plane toughness on the global dynamic response of the fuselage barrel was found. Actually, the configuration characterised by the highest toughness showed the lowest vertical deflection and the most significant damage in the cargo area. Moreover, this configuration seems to be more sensitive to the unbalancing of the mass in the lateral direction producing the most significant asymmetry in damage distributions between the left and the right side of the fuselage. On the other hand, the configurations characterised by the lowest in-plane toughness showed the maximum vertical deflection and a more distributed damage evolution in the whole structure, leading to less significant damage in the cargo area. Finally, a low sensitivity to the unbalancing of the mass in the later direction was found.