SPH-FEM Design of Laminated Plies under Bird-Strike Impact

: Composite laminates can potentially reduce the weight of aircrafts; however, they are subjected to bird strike hazards in civil aviation. To handle their nonlinear dynamic behaviour, in this study, the impact damage of composite laminates were numerically evaluated and designed by means of smoothed particle hydrodynamics (SPH) and the ﬁnite element method (FEM) to simulate the interaction between bird projectiles and the laminates. Attention was mainly focused on the di ﬀ erent damage modes in various laminates’ plies induced by bird impact on a square laminated plate. A continuum damage mechanics approach was exploited to simulate damage initiation and evolution in composite laminates. Damage maps were computed with respect to di ﬀ erent ply angles, i.e., 0 ◦ , 45 ◦ and − 45 ◦ . The damage distributions were comparatively investigated, and then the ply design was considered for crashworthiness improvement. The results aim to serve as a design guideline for future prototype-scale bird strike studies of complex laminated structures.


Introduction
As an important soft-body impact, structural damage induced by bird strike has become an increasingly serious and catastrophic issue in the civil aviation industry [1,2], particularly threatening the safety of aircraft and passengers, e.g., loss of aircrafts and even the loss of lives. Therefore, there is a great need to improve the bird strike resistance of aircraft structures by developing high-performance materials and better crashworthiness designs [3]. Generally, bird strike is a short duration (generally in the millisecond range) and high-intensity loading event, with complex coupling interactions. Under the bird soft-body impacts, some critical parameters, e.g., geometry effects of birds [4][5][6][7], material's models of birds [8][9][10], initial velocities and impact angles of the bird projectile [10][11][12][13], etc., can influence the deformations and damage mechanisms of the front-facing structures.
Due to their fine mechanical properties, composite laminates are now used as key materials for aircraft's front-facing components, e.g., the fan blades of aeroengine structures. Composite laminates offer the potential for reducing the weight of civil aircrafts, and hence some civil aircraft models in China, e.g., C919, are seeking to use composite laminates in key components in the future. Meanwhile, bird strike certification compliance, e.g., FAR Part 25, should be considered [14,15]. For these applications, composite laminates undergo complex internal damage when they are subjected to foreign object damage (FOD) via impacts. In terms of composite materials under impact loadings, the laminated layup will be critically vital [16][17][18]. With much more complexity compared to metal materials, several failure modes in composite materials can be caused by impact loading, such as fiber fracture, matrix cracking, delamination, etc. For the damage and failure modes of composite materials upon FOD, influential factors mainly include structural geometries [15,19], impact velocity/locations [20][21][22], properties of the projectile [20], layer configurations [17,19,23,24], fiber and matrix properties [25], etc.
for the analytical relation equation as pv = RT. Let ρ be the current mass density; in order to model the hydrodynamic response of fluidic bird's material accurately, an EOS needs to be employed as p = f (ρ, E m ) (1) where E m is the internal energy per unit mass. Common forms of EOS that are extensively adopted include polynomial, Mie-Grüneisen and Murnaghan forms [2]. In this article, in order to predict the behaviour of birds well, a Mie-Grüneisen EOS is adopted as where p H and E H represent the Hugoniot pressure and specific energy (per unit mass), respectively. The Hugoniot pressure is referred to as the initial pressure peak in the contact point in a perpendicular impact, and is mathematically defined as p H = ρ 0 U 0 U s . It can be seen that the Hugoniot pressure is a function of the initial density of the impactor (ρ 0 ), the initial velocity of the impactor (U 0 ), and the shock wave velocity (U s ). Γ is the Grüneisen ratio defined as where Γ 0 represents a material constant, ρ 0 represents the reference density and ρ is the current density. The Hugoniot energy (E H ) can be related to the Hugoniot pressure (p H ) by where η = 1 − ρ 0 /ρ. From the above equations, we can yield The linear Mie-Grüneisen EOS (called a U s − U p EOS) was employed to describe a linear relationship between the shock and particle velocities. A typical U s − U p EOS can be given as The linear relationship between the linear shock velocity U s and the particle velocity U p can be defined by c 0 and s, as U s = c 0 + sU p . Under the above assumptions, the linear U s − U p Hugoniot form can be deduced as In the present simulation, SPH modelling was conducted in ABAQUS/Explicit by using PC3D element. The EOS parameters c 0 , s and Γ 0 are material constants with values of c 0 = 1483 m/s, s = 0 and Γ 0 = 0. Different substitute bird impactor geometries were proposed, e.g., the cylinder, the cylinder with hemispherical ends, the ellipsoid and the sphere [39], but there is no standardized artificial bird shape. The bird body here is modelled as a cylinder (aspect ratio = 2:1) projectile, with 200 mm length and 100 mm diameter. The bird strike velocity was fixed at 60 m/s. The projectile has an initial density of ρ 0 = 950 kg/m 3 , and the weight of the projectile model is ≈1.5 kg. A total of 473 particles were used to model the bird projectile, and the distances between the particles were approximately 10 mm.
For calibration of the material constants, an experiment was conducted on a 600 mm × 600 mm × 2 mm steel plate under a 84.2 m/s bird strike [15]. For convenience in studying the damage evolution, we decided to use a square plate. The impact damage of composite laminates can be efficiently modelled by means of the finite element method (FEM) [19,32,42]. For the fiber breakage and matrix cracking simulation, continuum damage mechanics (CDM) enables one to predict with fine accuracy the onset and growth of damages, by introducing degradation factors for the material properties d f , d m and d s for fiber, matrix and shear modes, respectively. The damaged elasticity matrix C d has the following form The degradation factors range from 0 to 1, corresponding to the intact and fully-damaged states, respectively. In the present simulations, the laminate edge length and projectile length are 500 mm and 200 mm, respectively. The initial impact radius is 50 mm and the initial distance between the projectile and the target is 100 mm, which defines the initial position of the projectile's centre of mass. The thickness of the composite laminated plate is 1.2 mm, with an initial laminate layup [0/45/0/−45] s . Namely, eight-ply laminate was considered, which is similar to [27]. The geometric details of the plate and the bird are as shown in Table 1. The plate was clamped, supported on the four edges. By increasing the mesh density until the mesh showed no significant change in results, the refined mesh density was determined by performing repeated analyses; in particular, it was found that the projectile meshes should be symmetrical. The final numerical model of the impacted plate consists of 2500 conventional shell finite elements (linear quadrilateral elements S4R in ABAQUS). Hashin's criteria formulation [43] was adopted for evaluating the failure modes, which can enable the prediction of fiber and matrix failure onsets in compression or tension for each mode.
The orthotropic elasticity, mass density and strength parameters of a unidirectional CFRP material are listed in Table 2 for the present simulations. In the table, E 11 represents the Young's modulus in the fiber direction; E 22 represents the Young's modulus in the direction perpendicular to the fibers; G 12 represents the shear modulus; ν 12 represents Poisson's ratios; X T represents the longitudinal tensile strength; X C represents the longitudinal compressive strength; Y T represents the transverse tensile strength; Y C represents the transverse compressive strength; and S 12 , S 13 and S 23 represent the shear strength.
For the FEM (CDM) modelling of the unidirectional composite lamina, the in-plane fracture energies [44]: longitudinal tensile fracture energy, longitudinal compressive fracture energy, transverse tensile fracture energy, and transverse compressive fracture energy, are needed, of which the values are listed in Table 3. Figure 1 shows an example of the shape of the cylinder artificial bird projectile made of gelatin jelly in bird strike experiments; and the final SPH-FEM numerical model is shown in Figure 2. The initial ply angles are given in Figure 3. Table 2. Mechanical properties of the CFRP lamina [45].

Property Type Notation Value
Mass density ρ (kg/m 3 ) 1600 Orthotropic elasticity  Orthotropic elasticity Table 3. In-plane fracture properties of the composite lamina.   Orthotropic elasticity Table 3. In-plane fracture properties of the composite lamina.

Results and Discussion
For such analysis of an extremely discontinuous processes, an explicit dynamic analysis was carried out using ABAQUS/Explicit, which is considered to be computationally efficient for the analysis of large models with relatively short dynamic response times. The projectile begins to initially contact the impacted plate at approximately 1.7 ms, and the total simulation duration was 10 ms.
Bird impacts with high kinetic energy, can produce severe laminar damage, e.g., fiber fracture and matrix cracking. In the failure modes of composite laminates, the fiber fracture can result in catastrophic failure of overall laminated structures; therefore, in this section, attention will be firstly focused on the fiber-tension failure mode if not specified. Then, some other critical failure modes will be also evaluated.

Results and Discussion
For such analysis of an extremely discontinuous processes, an explicit dynamic analysis was carried out using ABAQUS/Explicit, which is considered to be computationally efficient for the analysis of large models with relatively short dynamic response times. The projectile begins to initially contact the impacted plate at approximately 1.7 ms, and the total simulation duration was 10 ms.
Bird impacts with high kinetic energy, can produce severe laminar damage, e.g., fiber fracture and matrix cracking. In the failure modes of composite laminates, the fiber fracture can result in catastrophic failure of overall laminated structures; therefore, in this section, attention will be firstly focused on the fiber-tension failure mode if not specified. Then, some other critical failure modes will be also evaluated. Figure 4 illustrates the damage maps of the fiber tension of the first two plies: (a) the first ply at the impacted side; (b) the second ply at the impacted side. The center of the rectangular laminate is (0,0,0), and the center of the initial impact area is also (0,0,0) as the white dots (projectile particles) indicate in Figure 4. It can be observed that the initial impact location is the highly damaged area, while the damage propagations are along different directions in Figure 4a,b. Furthermore, the damage at the boundaries are very different. Damage also propagates along the boundary, and the severity is also rather high due to high tension stress in the clamped boundary. Results show that damage initiates and mainly propagates along the fiber direction of the ply. In 0 • -fiber ply, the damage is almost widely distributed from the impacted area to the boundary; however, for 45 • -fiber ply, intermittent regions appear in the damage from the impacted area to the boundary. This phenomenon indicates that bird impact damage may be locally distributed for such clamped 45 • ply. Figure 5 shows the ultimate damage maps (damage at 10 ms) of the first two plies in the rear side: (a) the first ply at the rear side; (b) the second ply at the rear side. The damaged regions have an approximately similar area, but along different directions because of the different fiber directions. The local contact regions are the highest damaged, while the damage along the boundaries are less distributed for both cases. This is because the boundary of the rear side experiences compressive deformations.

ms.
Bird impacts with high kinetic energy, can produce severe laminar damage, e.g., fiber fracture and matrix cracking. In the failure modes of composite laminates, the fiber fracture can result in catastrophic failure of overall laminated structures; therefore, in this section, attention will be firstly focused on the fiber-tension failure mode if not specified. Then, some other critical failure modes will be also evaluated.  indicate in Figure 4. It can be observed that the initial impact location is the highly damaged area, while the damage propagations are along different directions in Figure 4a,b. Furthermore, the damage at the boundaries are very different. Damage also propagates along the boundary, and the severity is also rather high due to high tension stress in the clamped boundary. Results show that damage initiates and mainly propagates along the fiber direction of the ply. In 0°-fiber ply, the damage is almost widely distributed from the impacted area to the boundary; however, for 45°-fiber ply, intermittent regions appear in the damage from the impacted area to the boundary. This phenomenon indicates that bird impact damage may be locally distributed for such clamped 45° ply.   Figure 6 compares the damage maps of the middle two −45 • plies (the ply 4 and 5). The impacted area experienced local damage along the −45 • direction, while the non-contacted zones had almost zero damage. The reasons can be divided into two portions. First, the middle layers have a smaller stress response than the front or rear layers do. Second, −45 • plies have a longer distance to the clamped boundary, i.e., larger area to bear the impact loading. According to this damage mechanism of the laminate, it can be expected that ±45 • plies will give more locally-distributed impact damage. It can be observed from Figures 4-6 that correlation exists between the ply angles and damage modes. Considering these results, a designed plies can be proposed, as in Figure 7, to possibly achieve more locally-distributed impacted damage modes, by using ±45 • plies. Figure 5 shows the ultimate damage maps (damage at 10 ms) of the first two plies in the rear side: (a) the first ply at the rear side; (b) the second ply at the rear side. The damaged regions have an approximately similar area, but along different directions because of the different fiber directions. The local contact regions are the highest damaged, while the damage along the boundaries are less distributed for both cases. This is because the boundary of the rear side experiences compressive deformations.  Figure 6 compares the damage maps of the middle two −45° plies (the ply 4 and 5). The impacted area experienced local damage along the −45° direction, while the non-contacted zones had almost zero damage. The reasons can be divided into two portions. First, the middle layers have a smaller stress response than the front or rear layers do. Second, −45° plies have a longer distance to the clamped boundary, i.e., larger area to bear the impact loading. According to this damage mechanism of the laminate, it can be expected that ±45° plies will give more locally-distributed impact damage. It can be observed from Figures 4-6 that correlation exists between the ply angles and damage modes. Considering these results, a designed plies can be proposed, as in Figure 7, to possibly achieve more locally-distributed impacted damage modes, by using ±45° plies.   The designed plies of the composite laminate consist of all ±45 • plies. Figure 8 gives the damage maps in the all plies (i.e., the envelope damage of each ply), for both the initial and designed ply angles. The white dotted boxes in Figure 8a,b have the same shape and area, which is for comparative purposes. The results indicate that bird-impact damage is globally-distributed for the initial laminate; but locally-distributed for the designed one. The highest-damaged regions also distribute along the clamped boundary in the initial plies. The results of these final damage maps indicate that the effect of ply angles can greatly influence the damage modes, which are mainly through the different damaged areas. Although energy that was slightly higher dissipated in the laminate with 0 • plies (as shown in Figure 9), attention had to mainly focus on the damaged areas for the purpose of future real applications. Furthermore, it can be seen that the maximum damage index shows local characteristics, remaining almost identical (≈0.56).  The designed plies of the composite laminate consist of all ±45° plies. Figure 8 gives the damage maps in the all plies (i.e., the envelope damage of each ply), for both the initial and designed ply angles. The white dotted boxes in Figure 8a,b have the same shape and area, which is for comparative purposes. The results indicate that bird-impact damage is globally-distributed for the initial laminate; but locally-distributed for the designed one. The highest-damaged regions also distribute along the clamped boundary in the initial plies. The results of these final damage maps indicate that the effect of ply angles can greatly influence the damage modes, which are mainly through the different damaged areas. Although energy that was slightly higher dissipated in the laminate with 0° plies (as shown in Figure 9), attention had to mainly focus on the damaged areas for the purpose of future real applications. Furthermore, it can be seen that the maximum damage index shows local characteristics, remaining almost identical (≈ 0.56).  Figure 10 describes the time histories of kinetic energies in the initial and designed models. In the impact process, the overall trends and curve shapes of the two kinetic energies show agreement. However, the kinetic energy curve of the designed case shows a right shift in the graph compared with the initial case. Such right shifting shows that the impact motion experienced a time lag in the designed case, which can be attributed to the changed laminated plies.  Figure 10 describes the time histories of kinetic energies in the initial and designed models. In the impact process, the overall trends and curve shapes of the two kinetic energies show agreement. However, the kinetic energy curve of the designed case shows a right shift in the graph compared with the initial case. Such right shifting shows that the impact motion experienced a time lag in the designed case, which can be attributed to the changed laminated plies.
After comparing the fiber-tension damage mode, it is also worthwhile examining the other damage modes. It was observed that very small, similar damage appeared in the fiber-compression and matrix-compression modes of both the initial and designed laminates. Figure 11 compares the envelope matrix-tension damage maps of the initial (a) and designed (b) laminates. It is observed that both the initial and ±45 • plies shared a similar damage map in this mode, with most of the regions damaged. Figure 10 describes the time histories of kinetic energies in the initial and designed models. In the impact process, the overall trends and curve shapes of the two kinetic energies show agreement. However, the kinetic energy curve of the designed case shows a right shift in the graph compared with the initial case. Such right shifting shows that the impact motion experienced a time lag in the designed case, which can be attributed to the changed laminated plies.  After comparing the fiber-tension damage mode, it is also worthwhile examining the other damage modes. It was observed that very small, similar damage appeared in the fiber-compression and matrix-compression modes of both the initial and designed laminates. Figure 11 compares the envelope matrix-tension damage maps of the initial (a) and designed (b) laminates. It is observed that both the initial and ±45° plies shared a similar damage map in this mode, with most of the regions damaged. For the plate surface facing the bird projectile, the fiber-compression damage is also critical, due to the possibility of fiber buckling failure under compressive deformation. Figure 12 displays the fiber-compression damage of the foremost plies (i.e., the ply directly contacted by the bird projectile) in the two laminates: (a) initial; (b) designed. Local damage distributions in the vicinity of the initially contacted regions can be observed for both models.  Figure 11. The envelope matrix-tension damage maps of the initial (a) and designed (b) laminates.
For the plate surface facing the bird projectile, the fiber-compression damage is also critical, due to the possibility of fiber buckling failure under compressive deformation. Figure 12 displays the fiber-compression damage of the foremost plies (i.e., the ply directly contacted by the bird projectile) in the two laminates: (a) initial; (b) designed. Local damage distributions in the vicinity of the initially contacted regions can be observed for both models.
The impact deflection can be considered as an indicator of the structural stiffness. Figure 13 depicts the deflections of the plate's centres in the two laminates. On the one hand, the peaks of the designed laminate were found to be a little larger than the initial one, which indicates that the ±45 • designed plies slightly decrease the structural local stiffness to a small degree. On the other hand, the curve of the designed case also shows a right shift, which agrees well with the results of the kinetic energies.
For the plate surface facing the bird projectile, the fiber-compression damage is also critical, due to the possibility of fiber buckling failure under compressive deformation. Figure 12 displays the fiber-compression damage of the foremost plies (i.e., the ply directly contacted by the bird projectile) in the two laminates: (a) initial; (b) designed. Local damage distributions in the vicinity of the initially contacted regions can be observed for both models. The impact deflection can be considered as an indicator of the structural stiffness. Figure 13 depicts the deflections of the plate's centres in the two laminates. On the one hand, the peaks of the designed laminate were found to be a little larger than the initial one, which indicates that the ±45° designed plies slightly decrease the structural local stiffness to a small degree. On the other hand, the   To further compare the inherent structural global properties for the purpose of the load-carry capability, the modal frequencies were calculated as shown in Figure 14 for the first 10 modes. The results indicate that the designed laminate exhibits higher modal frequencies overall, except the third mode. In addition, the fundamental frequencies are almost the same, which can validate the fine dynamic properties of the designed laminated plies. To further compare the inherent structural global properties for the purpose of the load-carry capability, the modal frequencies were calculated as shown in Figure 14 for the first 10 modes. The results indicate that the designed laminate exhibits higher modal frequencies overall, except the third mode. In addition, the fundamental frequencies are almost the same, which can validate the fine dynamic properties of the designed laminated plies.
Hereto, the SPH-FEM numerical simulations and analysis were conducted to examine the plies effect on the damage modes of a rectangular laminated composite, which is a flat structure. However, because the laminates were flat, slicing of the projectile did not occur. This aspect needs to be analysed for real aero-engine fan blades in the future.
To further compare the inherent structural global properties for the purpose of the load-carry capability, the modal frequencies were calculated as shown in Figure 14 for the first 10 modes. The results indicate that the designed laminate exhibits higher modal frequencies overall, except the third mode. In addition, the fundamental frequencies are almost the same, which can validate the fine dynamic properties of the designed laminated plies. Hereto, the SPH-FEM numerical simulations and analysis were conducted to examine the plies effect on the damage modes of a rectangular laminated composite, which is a flat structure. However, because the laminates were flat, slicing of the projectile did not occur. This aspect needs to be analysed for real aero-engine fan blades in the future.

Conclusions
In this article, the damage modes of impacted laminates with various ply angles were numerically investigated, by using SPH-FEM computation. It was found that the damage maps of the square laminates are sensitive to ply angles. However, the local maximum damage of the plate was less affected by the ply angle effect. As the laminated plies changed, the damage in the composite laminate changed from globally-distributed damage to a more locally-distributed mode by changing from 0 • ply to 45 • ply; while the boundary effect still exists with the deceased degree. Such a damaged region transition allows us to design the laminated plies, which is one of the critical factors for discussing the bird strike resistance of composite laminates. By analysing the damage in different plies, a ±45 • plies design was proposed to achieve more locally-distributed impact damage. The simulation results have validated such a design in terms of impact damage, indicating an improvement in crashworthiness. The results can serve as a design guideline for future prototype-scale bird strike studies of complex laminated structures.