A Numerical Assessment on the Influences of Material Toughness on the Crashworthiness of a Composite Fuselage Barrel

Abstract

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 influence of the material fracture toughness on the capability of a composite fuselage barrel to tolerate an impact on a rigid surface. Three different material configurations with different 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 effects 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 configurations.

1. 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 exploring the mechanical response of a fuselage during a crash phenomenon. Subsequently, models with non-linear springs and concentrated masses were introduced [17,18].

Today, the use of finite element simulations with explicit methodologies allows to mimic, with good accuracy, the evolution of the crash phenomenon of a fuselage during a fall, allowing to identify both the damages and their development [19,20,21,22,23].

The high complexity of the phenomenon under examination is amplified by the complexity of the modern aircrafts. In particular, the huge size of the fuselage structure, the number of different sub-components characterised by several sources of non-linearity (geometry, materials, contacts, etc.) make the validation of numerical models very complicated.

In addition to all these conditions, new aircraft generations are characterised by a further source of complexity represented by the use of fibre-reinforced composite materials. Composite fibre-reinforced material, by their nature, are subjected to intra-laminar and inter-laminar damages, which are not easily numerically analysed [24].

The use of fibre-reinforced composite materials, even with the intrinsic advantages in terms of weight-performing ratio, is source of further complexities in terms of components crash behaviour. Indeed, the absence of a material plastic phase [25] and the presence of different interacting mechanisms, such as fibre and/or matrix breakage, delamination and fibre matrix de-bonding, for absorbing impact energy [26,27,28] makes the crash behaviour even more complex with respect to conventional metallic materials. Furthermore, a crash events response is a highly geometrical and materially nonlinear event [29]. Hence the crash behaviour of a complex composite structure is strongly influenced by the failure characteristics and strength of its material system. Consequently, the use of effective theories for the prediction for the onset and evolution of the damage becomes mandatory to understand and describe the complex mechanism of breakages in composite materials during a crash event.

In the literature, several theories are proposed to accurately describe the failure mechanisms for laminated aeronautical structures [30,31,32,33,34,35,36]. The theory of continuum damage mechanics (CDM) is mainly adopted for the assessment of intralaminar damages. According to this theory, the stiffness reduction can be physically compared to the presence of distributed micro-fractures or defects. The presence of micro damage, phenomenologically, can be represented by a variation of internal state variables of the material. Alternative methodologies for the simulation of degradation of the material properties (MPDM) based on a thermo-dynamic model can be found in [37,38,39,40]. These methodologies demonstrated to be very effective in conjunction with finite element modelling.

In general, the intra-laminar damages evolution in composite materials and then in complex composite structures is strictly influenced by fracture toughness for the four fracture modes (matrix and tensile and compressive fibre). Moreover, the experimental determination of these fracture toughness-related parameters is rather complex [41,42,43].

The purpose of this paper is to numerically evaluate the influence of these fracture toughness-related parameters on the mechanical behaviour of a composite fuselage section during a crash event. The Abaqus/Explicit platform has been adopted to perform the numerical simulations of the mechanical behaviour of the fuselage barrel subjected to a crash event on a rigid ground. The structural subcomponents of the fuselage were detailed discretised by adopting finite elements with a three-dimensional formulation. Three different configurations with material systems characterised by different in-plane fracture toughness values were analysed and compared to each other to assess the effects on the energy absorption capabilities of the structures and on the consequent deformations and stresses distributions.

In Section 2, the theoretical background is introduced, while in Section 3, the description of the adopted finite element model is provided. In Section 4, the numerical results obtained from the different analysed configurations are introduced and discussed.

2. 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.

$$F_{ft}={\left(\frac{{\stackrel{⏜}{\sigma }}_{11}}{X_T}\right)}^2+{\left(\frac{{\stackrel{⏜}{\sigma }}_{12}}{S_L}\right)}^2=1$$

(1)

$$F_{fc}={\left(\frac{{\stackrel{⏜}{\sigma }}_{11}}{X_C}\right)}^2=1$$

(2)

$$F_{mt}={\left(\frac{{\stackrel{⏜}{\sigma }}_{22}}{Y_T}\right)}^2+{\left(\frac{{\stackrel{⏜}{\sigma }}_{12}}{S_L}\right)}^2=1$$

(3)

$$F_{mc}={\left(\frac{{\stackrel{⏜}{\sigma }}_{22}}{2S_T}\right)}^2+\left[{\left(\frac{Y_C}{2S_T}\right)}^2-1\right]\cdot \frac{{\stackrel{⏜}{\sigma }}_{22}}{Y_C}+{\left(\frac{{\stackrel{⏜}{\sigma }}_{12}}{S_L}\right)}^2=1$$

(4)

In Equations (1)–(4) ${\stackrel{⏜}{\sigma }}_{11}$, ${\stackrel{⏜}{\sigma }}_{22}$, ${\stackrel{⏜}{\sigma }}_{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.

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).

$$d_i=\frac{{\delta }_{i,eq}^t\left({\delta }_{i,eq}-{\delta }_{i,eq}^0\right)}{{\delta }_{i,eq}\left({\delta }_{i,eq}^t-{\delta }_{i,eq}^0\right)};{\delta }_{i,eq}^0\le {\delta }_{i,eq}\le {\delta }_{i,eq}^t;i\in \left(f_c,f_t,m_c,m_t\right)$$

(5)

where ${\delta }_{i,eq}$ is equivalent displacement; ${\delta }_{i,eq}^0$ is equivalent displacement on set damage; ${\delta }_{i,eq}^t$ 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.

$${\delta }_{i,eq}^t=\frac{2G_c}{{\sigma }_{i,eq}^0}$$

(6)

In Equation (6), ${\sigma }_{i,eq}^0\mathrm{and}{\delta }_{i,eq}^0$ 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. Equivalent stress and displacement definitions.

Load Condition	Equivalent Stress	Equivalent Displacement
Fibre Tension	$\frac{L_c(\langle {\sigma }_{11}\rangle \langle {\epsilon }_{11}\rangle +{\sigma }_{12}\cdot {\epsilon }_{12})}{{\delta }_{ft,eq}}$	$L_c\sqrt{\langle {\epsilon }_{11}\rangle {}^2+{\epsilon }_{12}^2}$
Fibre Compression	$\frac{L_c\langle -{\sigma }_{11}\rangle \langle -{\epsilon }_{11}\rangle }{{\delta }_{fc,eq}}$	$L_c\langle -{\epsilon }_{11}\rangle$
Matrix Tension	$\frac{L_c(\langle {\sigma }_{22}\rangle \langle {\epsilon }_{22}\rangle +{\sigma }_{12}\cdot {\epsilon }_{12})}{{\delta }_{mt,eq}}$	$L_c\sqrt{\langle {\epsilon }_{22}\rangle {}^2+{\epsilon }_{12}^2}$
Matrix Compression	$\frac{L_c(\langle -{\sigma }_{22}\rangle \langle -{\epsilon }_{22}\rangle +{\sigma }_{12}\cdot {\epsilon }_{12})}{{\delta }_{mc,eq}}$	$L_c\sqrt{\langle -{\epsilon }_{22}\rangle {}^2+{\epsilon }_{12}^2}$

, the relations adopted to evaluate the equivalent stress and displacement are introduced.

In Table 1, L_C and < > are, respectively, the element characteristic length and the Macauley bracket operator [45].

3. 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.

3.1. 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.

3.2. Material System

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

Table 2. UNIDIRECTIONAL (UD) material system properties.

Fibre Composite Material
Young’s Modulus, E11 [MPa]	137,500
Young’s Modulus, E22 [MPa]	8200
Shear Modulus, G12 [MPa]	3950
Shear Modulus, G13 [MPa]	3950
Shear Modulus, G23 [MPa]	3950
Poisson’s Ratio, ν12 = ν13 = ν23	0.35
Fibre Tensile Strength F1t [MPa]	1890
Fibre Compressive Strength F1c [MPa]	1008
Matrix Tensile Strength F2t [MPa]	86.5
Matrix Compressive Strength F2c [MPa]	112
In-Plane Shear Strength, S12 [MPa]	95
Out-Plane Shear Strength, S23 [MPa]	100
Density [ton/mm3]	1.9 e−9
Ply Thickness [mm]	0.129

,

Table 3. Woven fabric material system properties.

Woven Fabric Material
Young’s Modulus, E11 [MPa]	55,000
Young’s Modulus, E22 [MPa]	55,000
Shear Modulus, G12 [MPa]	3363
Shear Modulus, G13 [MPa]	3363
Shear Modulus, G23 [MPa]	3363
Poisson’s Ratio, ν12 = ν13 = ν23	0.30
Fibre Tensile Strength F1t [MPa]	650
Fibre Compressive Strength F1c [MPa]	650
Matrix Tensile Strength F2t [MPa]	650
Matrix Compressive Strength F2c [MPa]	650
In-Plane Shear Strength, S12 [MPa]	150
Out-Plane Shear Strength, S23 [MPa]	150
Density [ton/mm3]	1.97 e−9
Ply Thickness [mm]	0.25

and

Table 4. Aluminium Al 2024 mechanical properties

Aluminium Al2025
Young’s Modulus, E [MPa]	70,000
Poisson’s ratio, ν12 = ν13 = ν23	0.33
Yield stress [MPa]	369
Ultimate Tensile Stress [MPa]	469
Density [ton/mm3]	2.7 e−9

.

For the aluminium plates, a thickness of 8 mm was chosen.

Table 5. Sub-components stacking sequences.

Stacking Sequence
Skin	[90/45/0/45]s
Stringer	[45/45/0/0/90/0/0]s
Frames	[90/45/0/45/-45/90/45/45/0/-45/45/-45]s
Components in Woven Fabric Material	[0/45/0/-45]s

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.

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. Fracture toughness energies.

		${\mathit{G}}_{1\mathit{C}}^{\mathit{T}}[\mathbf{kJ}/{\mathbf{m}}^2]$	${\mathit{G}}_{1\mathit{C}}^{\mathit{C}}[\mathbf{kJ}/{\mathbf{m}}^2]$	${\mathit{G}}_{2\mathit{C}}^{\mathit{T}}[\mathbf{kJ}/{\mathbf{m}}^2]$	${\mathit{G}}_{2\mathit{C}}^{\mathit{C}}[\mathbf{kJ}/{\mathbf{m}}^2]$
I	UNIDIRECTIONAL UD	24	10	3	6
I	Woven	8	8	8	8
II	UNIDIRECTIONAL UD	44	16	3	6
II	Woven	12	12	12	12
III	UNIDIRECTIONAL UD	26	12	3	6
III	Woven	6	6	6	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.

3.3. 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.

3.4. 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.

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.

4. 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.

4.1. Configuration I

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.

4.2. Configuration II

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).

4.3. Configuration III

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.

Figure 9B–D highlights the damaged structure with particular focus on the maximum deflection point. Moreover, in 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.

4.4. 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 10A and Figure 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 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.

5. 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.