Sandwich shield subjected to bird impact: Use of surrogate models for influencing parameter analysis and shield behaviour understanding

Abstract

In this work, the behaviour of a sandwich shield subjected to a 1.82 kg bird impact at 175 m/s is studied using a finite element model. The 6 most influential design parameters are varied and their effects on the shield behaviour and on the target protection are assessed. First, we try to establish an engineer's visualization by varying parameters 2 × 2 using three 5-levels full-factorial design of experiments. These three 2D design of experiments enable us to visualize precisely the different effects of each parameter. Then a full sensitivity analysis (6D) is performed using a Latin Hypercube sampling to assess the possible interactions between parameters. Surrogate models are constructed using the Gaussian Process framework to follow the variation of the outputs in the 6D design space. These surrogate models are finally studied using two statistical methods: the Sobol′ method and the Morris method. The methodology developed in this study enables to improve the understanding of the behaviour of a shield under a soft body impact, as a first step towards a shield design tool.

Keywords

Sandwich shield impact design bird strike soft body impact parametric study Gaussian process

Introduction

During its flight, one of the major threats an aircraft can encounter is the collision with a bird. Such collisions are known to occur frequently, and a 2008 study by the European Aviation Safety Agency [1] estimates the occurrence at 186 per million flying hours. The possible damages of such an impact can be very diverse, due to the wide range of possible impacting scenarios, from multiple small flock birds (weighting approximatively 50 g) to heavy migratory birds (up to 4 kg). Moreover, the birds can impact all the forward facing structures, i.e. the nose, the windshield, the wings leading edge, the empennage, the engines, etc.

To ensure the protection of passengers and crew on a commercial aircraft, aviation authorities [2,3] require that the plane should be able to continue its flight and land safely after a 1.82 kg (4 lb) impact at operational speed at sea level (typically around 175 m/s). In the case of an impact on the aircraft nose, the main danger is the failure of the pressurised bulkhead, electrical systems and equipment situated behind the radome instruments. Thus, to meet the certification requirements, a shield is placed in front of the bulkhead, behind the instruments (cf. Figure 1).

Due to their high specific stiffness and the good absorption properties of their core [4], sandwich structures are the ideal candidates for such shields and are used for this application by most aircraft manufacturers. In this work, we focus on the study of the behaviour of such a sandwich shield under the certification impact (1.82 kg at 175 m/s).

The first studies about bird strike, in the late 1970s, were focused on understanding the behaviour of a bird impacting a plate at such speed. Using extensive testing, Barber et al. [5 –7] showed that the bird behaves similarly to a fluid during impact: the force transmitted to the target starts with a peak and is followed by a long plateau corresponding to a steady flow. They also showed that, during experiments, the bird can be substituted by an impactor made of gelatine with 10% porosity, with good improvements on repeatability. Consequently, these substitutes are used in tests in most latter studies.

Due to the great number of possible impacting scenarios (size and speed of bird, impacted structure, etc.), the literature on bird strike is quite large. In a 2011 review, Heimbs [8] lists more than 190 numerical studies on bird strike, and identifies for each the impacted structure (windshield, leading edge, etc.), the impact case (mass of bird, speed, geometry) and the bird modelling strategy used. Three different numerical strategies are commonly used to model a bird: Lagrangian, Arbitrary Lagrangian Eulerian (ALE) or smoothed particle hydrodynamics (SPH). These methods and their differences are presented in detail in [8].

Concerning the impactor shape and material laws, a study by Airoldi and Cacchione [9] showed that the simulations were closer to test data when the material laws used represented water with 10% porosity and when the impactor was modelled as a cylinder with hemispherical caps and a length to diameter ratio of 1.6. More complex shapes and material laws have also been studied in [10].

Only a few of these studies use a flat sandwich structure as target and focus on its behaviour. We present here these five papers:

Hanssen et al. studied a two-layer sandwich panel [11]. Using simulation with the ALE strategy, they were able to represent the observed experimental behaviour, including failure at the clamping bolts. In a second part of their work, they used their model to minimize the core thickness of a simple sandwich with 0.8 mm thick aluminium skins and a core made of 150 kg/m³ aluminium foam. The lighter sandwich able to stop the bird without front skin tearing had a 150-mm thick core.

Hohe et al. reproduced this load case [12] and showed that, using a graded core of three layers with increasing density, the front skin strain was spread on a greater surface, with the maximum strain being lower. They concluded that the use of a graded core could improve the resistance of the sandwich shield.

Liu et al. studied a sandwich with aluminium honeycomb core and a two-layer sandwich with the same core material and the same total height [13]. They created a finite element model to simulate the load case and were able to represent correctly both impacts.

Hedayati et al., in a numerical study, assessed the effect of the position of the middle skin in a two-layer sandwich, keeping the total height constant [14]. They showed that the minimal backward deflection is obtained for the shield with balanced cores (i.e. middle plate in the centre). In a second part, they showed that for cores with different densities, the best result was obtained with the lighter core facing the impactor.

The same year, Liu et al. studied the influence of skins and core thickness at constant mass on a simple sandwich with an aluminium foam core [15]. They showed that both the backward deflection and the energy absorption of the panel increase when the core height increases and the skins’ thickness decreases. Then, they studied the effect of the middle plate position on a two-layer sandwich with a constant total height. They showed that the minimal backward deflection is obtained with a first layer height equal to zero, i.e. a simple sandwich with a double front skin.

This latter conclusion seems to be in contradiction with the conclusion of [14]. This difference can be explained by the fact that those two studies use slightly different boundary conditions (respectively, riveted or clamped) and by the fact that in [14], no bonding is modelled between the skins and the cores.

All these results suggest that the behaviour of the shield is influenced not only by the core and skins properties but also by the structural coupling between the two. Thus, when designing a new shield, it seems important not to study the skins and core properties separately but to study all the shield design properties together.

According to this conclusion, in this work, the skins and core properties are all considered as design parameters of the shield, along with the different heights and thicknesses. We then study the influence of these parameters on the behaviour of the shield but also on the target protection. This paper is organised as follows:

In the next section, the case study is defined and the finite element model used is presented. Then, the different outputs which will be followed throughout the study are presented and, according to a previous study [16], the most influential parameters are determined.

In the ‘2D parametric studies’ section, these parameters are studied 2 × 2 using full factorial design of experiments (DOEs). The effect of each parameter and its interactions with other parameters are studied and physical interpretations are proposed for the observed phenomenon.

In Section ‘Expanding the parametric study in 6D’, these 2D studies are expanded in full 6D by adding 100 new simulation points. Surrogate models are then created to follow the output variations in the 6D design space using one of the main reference frameworks in the machine learning domain. Statistical methods are then used to analyse the effects of each parameter and the interactions with others. The results obtained are confronted to the observations made on the 2D studies of the previous section.

Case study description

The radome shields used today in the industry have a geometry which is strongly dependent on all the other surrounding systems, and thus they are different from one aircraft to the other, but these shields are usually flat sandwiches, with an area of about 1 m², directly supported by the protected bulkhead. These sandwiches are typically made of two aluminium skins and a core made of metallic foam or honeycomb. The skin thickness ranges from 1.5 to 5 mm and the core height is usually about 100 mm.

In this work, we will study the shield under the impact of a 1.82 kg bird at 175 m/s, representing the certification case. This impact represents an initial kinetic energy of 27.9 kJ.

Finite element model

In order to limit the number of geometrical and boundary condition parameters, the limit conditions and the shield geometry have been simplified compared to a real engineering case. Moreover, to be able to simulate a great number of different shield designs at an acceptable cost, the finite element model is kept as simple as possible, and the material laws used are kept generic.

Geometry

The geometry chosen is an 800-mm square sandwich, supported by a rigid frame with a 400-mm square aperture in its centre. The sandwich has two skins 3 mm thick and a core with a height of 100 mm. Due to symmetry conditions, only a fourth of the shield is modelled using ABAQUS/EXPLICIT (cf. Figure 2).

The support is modelled with rigid shell elements R3D4. The sandwich shield is modelled as one part, with no possible debonding between the skins and the core. Five layers of reduced integration brick elements (C3D8R) are used for the back skin, 10 layers of fully integrated elements (C3D8) are used for the core and five layers of reduced integration elements (C3D8R) are used for the front skin. During the parametric study, the core and skins thicknesses will change but the number of elements used through thickness is kept constant. The in-plane element size is 10 mm × 10 mm, for a total number of elements for the sandwich of 32,000. The impacting bird is modelled as a hemispherical-ended cylinder with a radius of 55 mm and a length of 220 mm. A total of 3490 Lagrangian reduced integration elements (C3D8R) are used, with enhanced hourglass control. The contact between the bird and the front skin, and between the back skin and the support, is modelled with ABAQUS general contact without friction. The total simulated time is 6 ms, enough to simulate the whole impact and the shield rebound.

Material models

In order to model the bird material behaviour, a tabulated equation of state representing the behaviour of water with porosity is used. The elements are deleted when their true strain becomes greater than 500%, in order to limit the decrease of the simulation time step due to highly distorted elements while allowing to represent the flow of the bird. The bird density is 955 kg/m³, which gives a 1.82 kg bird.

The shield design parameters being varied during the parametric study, we give here the values for the central point of the design space. This reference case represents a shield with aluminium skins and an aluminium honeycomb core.

The skin material is represented using an isotropic material law. A Young modulus of 72 GPa and a Poisson coefficient of 0.33 are used, with a density of 2800 kg/m³. The plastic behaviour is modelled by ABAQUS Johnson-Cook isotropic hardening (equation (1))

σ_{p} = A + B . {(ε_{p})}^{n}

(1)

where A = 300 MPa, B = 400 MPa and n = 0.27. This model is used because it is able to represent the plastic behaviour of many different metals with only three parameters [17]. No skin failure is modelled, but the maximal skin in-plane strain will be studied to ensure that this hypothesis stays valid.

To simulate the core behaviour, a generic material law is implemented using a VUMAT user routine. An elasto-perfectly plastic law is used for all directions, followed by a densification for the out-of-plane compression. Such a strain–stress curve has been shown to be typical for the cellular materials used as sandwich cores [4]. Uncoupled shearing and compressive behaviour is assumed. The material law parameters are adapted from [18]:

– In-plane modulus: $E_{x}$ = $E_{y}$ = 1 MPa.

– In-plane plateaus: $σ_{p x}$ = $σ_{p y}$ = 0.9 MPa.

– In-plane shear modulus: $G_{x y}$ = 5 MPa.

– In-plane shear plateau: $τ_{p x y}$ = 1.2 MPa.

– Out-of-plane shear modulus: $G_{x y}$ = $G_{y z}$ = 200 MPa.

– Out-of-plane shear plateaus: $τ_{p x y}$ = $τ_{p y z}$ = 1.2 MPa.

– Out-of-plane modulus: $E_{z c}$ = 1000 MPa.

For out-of-plane compression, a plateau followed by densification behaviour has been modelled using equation (2).

σ_{z} = σ_{p c} + C . (e^{\frac{ε}{ε_{d}}} - 1)

(2)

with

σ_{p c}

= 0.9 MPa the core crushing plateau, and C = 0.01 MPa and

ε_{d}

= 0.37 defining the densification phase. The core density is set to 50 kg/m³.

Outputs studied

With this generic finite element model, it is possible to simulate the impact of a bird on different shields and to study the effect of the shield definition on the target protection. The simulation of one design case takes approximatively 2 h on four CPUs, using the CALMIP (Calculateur Midi-Pyrénées) supercomputer. This relatively small computing time allows simulating numerous design points during the parametric study. Figure 3 presents the evolution of the reference shield shape during impact.

In such a parametric study, it is obviously impossible to analyse all the differences between all the finite element simulations, due to the important number of design cases tested. It is then necessary to define simple outputs that can be followed during the parametric study and that define the shield behaviour.

Here, we decided to follow four protection criterions, which could be used as design criterion in a real engineering case, and representing the protection capacity of the shield:

The maximum total force applied to the support $F_{S u p p}$ . The total force is computed at each time step as the sum of all the contact forces between the back skin and the support. This force is then filtered using a moving mean with a window of 0.15 ms and $F_{S u p p}$ is defined as the maximum along time. For the reference shield, $F_{S u p p}$ is placed in Figure 4 on the curve of the total force applied by the back skin on the support. The total force applied by the bird on the front skin is also given for comparison.

The maximum pressure applied to the support $P_{s u p p}$ . This pressure is computed using the contact forces between the back skin of the shield and the support. The contact forces are filtered using a moving average with a window of 0.15 ms along time and a Gaussian filter on a 3 × 3 nodes window along space. $P_{s u p p}$ is then defined as the maximum along time and space of the pressure obtained (cf. Figure 5 for the evolution of pressure field in the reference case.).

The maximum backward deflection of the back skin $δ_{b s}$ . This criterion is often used in the industry since an important backward deflection can endanger critical systems situated behind the pressurised bulkhead.

The maximum in-plan strain of the front skin $ε_{f s}$ . Even if no skin rupture has been modelled in this work, this output allows to differentiate two shield designs by knowing how close they are of a potential front skin tearing.

Using these four criterions, it is possible to assess the capacity of a shield to protect the target, but these outputs give no information about the shield deformation. The goal of this parametric study being to understand the behaviour of the shield and its influence on the target protection, it is necessary to have outputs describing how the shield is deformed.

In a previous work [16], we described a behaviour analysis tool allowing to decompose the deformation of a shield into three modes: indentation, bending and crushing. The main principle of this tool is to project at each time step the vector representing the shield deformation on a basis defined a priori and describing the three deformation modes. It has been shown that this decomposition describes efficiently the deformation of any shield with only a small residue (less than 10% in norm) [16].

Using this tool, it is then possible to extract the maximum norm of each deformation mode along time. These three behaviour criterions will give information about how the shield behaves and is deformed during impact without the need to analyse in detail each simulation.

Choice of parameters to study

The finite element model described here exhibits an important number of parameters. Thus, it is necessary to choose which parameters to take into account before conducting the parametric study. In [16], we conducted on this finite element model a screening analysis to identify the most influential parameters on the shield behaviour, using the deformation-based behaviour analysis tool presented above. This screening has been extended to identify influences on the shield behaviour (three outputs) and the four protection criterions (four additional outputs).

Fifteen parameters or group of parameters have been studied independently. Along with its reference value, each parameter has been assigned a minimal and maximal value. Thirty-one cases have then been simulated: the reference point and the minimal and maximal case for each parameter (one factor at a time DOE). For each of these simulations, we measured the difference of shield behaviour with respect to the reference case, in order to rank the 15 parameters (mean rank on the seven outputs). The main conclusions of this screening study are as follows:

The most influential parameter is clearly the core crushing plateau $σ_{p c}$ . Any change in the value of this parameter induces important change in the shield behaviour.

The second most influential parameter is the core out-of-plane shearing plateau $τ_{p c}$ .

Then the parameters rank as follows: the core height $H_{c}$ , the front skin yield strength $A_{f s}$ , the support aperture size $L_{a p}$ , the front skin thickness $t_{f s}$ , the elastic modulus of the front skin, the out-of-plane elastic modulus of the core and then the other parameters.

It appears that the in-plane properties of the core have very little influence on the shield behaviour, and that the back skin design has little influence compared to the front skin design.

The effects of the parameters on the shield behaviour can be strongly non-linear, and sometimes even non-monotonic.

According to these conclusions, this parametric study will be focused on the six most influential parameters: the core out-of-plane crushing plateau $σ_{p c}$ , the core out-of-plane shearing plateau $τ_{p c}$ , the core height $H_{c}$ , the support aperture size $L_{a p}$ , the front skin thickness $t_{f s}$ and the front skin yield strength $A_{f s}$ . Moreover, it is necessary to choose a DOE allowing us to study the possible non-linearities and interactions, as they were not studied in the screening study.

2D parametric studies

In order to be able to analyse precisely the effects of each parameter and to be able to visualize the results more easily, the six identified parameters are first studied in pairs: the two core out-of-plane plastic plateaus $σ_{p c}$ and $τ_{p c}$ , the two front skin parameters $A_{f s}$ and $t_{f s}$ and the geometric parameters $H_{c}$ and $L_{a p}$ . To be able to assess the non-linearities and the possible interactions, five levels by parameters are chosen. When varying the inputs pairwise, remaining parameters are kept at the baseline values (presented in ‘Finite element model’ section). The minimal and maximal parameter values used are presented in Table 1 and are the same as in [16], representing materials typically seen in the industry for such applications.

Table 1.

Minimal and maximal values of the six parameters studied.

		Lower limit	Reference	Upper limit	Units
$σ_{p c}$	Core crushing plateau	0.2	0.9	4	MPa
$τ_{p c}$	Core shearing plateau	0.2	1.2	4	MPa
$A_{f s}$	Front skin yield stress	100	300	500	MPa
$t_{f s}$	Front skin thickness	2	3	4	mm
$H_{c}$	Core height	50	100	150	mm
$L_{a p}$	Support aperture size	300	400	500	mm

The DOEs used for these three 2D series are shown in Figure 6. The reference case is indicated with a full red dot.

For the ‘core’ DOE, the choice was made to simulate more designs toward the lower bound of both $σ_{p c}$ and $τ_{p c}$ because the screening study showed that a decrease of those parameters has more effect on the shield behaviour than an increase.

Effect of the core design

For each of the 25 simulations of this DOE, the seven criterions (three behaviour criterions and four protection criterions) are measured. These results are given in Figure 7 as a surface reconstructed from the simulation points (black) using linearly interpolated scheme. The red dot indicates the reference case and the value of the outputs at this point is shown by a red line on the colour scales.

We can see on this figure that the effects of the core design are strongly non-linear and that an interaction exists between the two design parameters. Two main effects can be observed.

First, there is a strong effect of $σ_{p c}$ , which appears only when $σ_{p c}$ is smaller than 0.8 MPa. This effect can be observed in Figure 7 for indentation, crushing, $F_{S u p p}$ and $ε_{f s}$ . For a decreasing value of $σ_{p c}$ , the shield shows increasing values of indentation and crushing modes, associated with reducing values of $F_{S u p p}$ and $ε_{f s}$ . Thus, this strong crushing behaviour is associated with a smaller force transmitted to the support and a smaller strain of the front skin. This effect can be explained by the fact that the force transmitted to the support is linked with the force transmitted between the two skins by the core. With a low crushing plateau of the core, this force is necessarily smaller. Regarding $ε_{f s}$ , its decrease may be explained by the fact that a strong crushing behaviour induces an indentation which is less located at the centre of the shield.

Therefore, for a low core out-of-plane crushing plateau, the shield has a deformation dominated by indentation and crushing, with comparatively small bending, whereas for a higher core out-of-plane crushing plateau, the shield shows mainly indentation.

This change of behaviour can be seen clearly in Figure 8: the shield in (a) (with a low core crushing plateau) shows a strong crushing behaviour, whereas the shield in (b) (the reference case) shows mainly indentation. On this figure, only the outline of the core in the XZ plane is represented (with the second half reconstructed by symmetry for easier visualization). The undeformed profiles of the shields are represented by the dashed line and the rigid support is shown in black.

The second effect which can be observed in Figure 7 is the increase of bending, $P_{S u p p}$ , $δ_{b s}$ and $ε_{f s}$ for a low core out-of-plane shearing plateau associated with a strong core crushing plateau. The effect of $τ_{p c}$ can be explained by the fact that the core out-of-plane shearing plateau has a direct influence on the sandwich bending resistance. Thus, a decrease of $τ_{p c}$ will lead to an increase of bending, and this effect is enhanced for a high core crushing plateau, which increases the sandwich crushing resistance. This increase of bending leads to an increase of the back skin backward deflection $δ_{b s}$ , which can be seen clearly in Figure 8 where the shield (c) (with low shearing plateau and a high crushing plateau) shows an important backward deflection. Similarly, the shield (c) shows an indentation which is more located at the centre of the shield, which can explain the increase of the front skin maximum strain $ε_{f s}$ .

Another consequence of the increase of bending is the increase of the maximum pressure transmitted to the support $P_{S u p p}$ . This consequence is explained by the fact that the bending of the shield localises the contact between the back skin and the support. Thus, $P_{S u p p}$ increases without any increase of $F_{S u p p}$ .

In summary, the core design strongly influences the behaviour of the shield. A low core crushing plateau induces a behaviour dominated by the core crushing. On the contrary, a low core shearing plateau associated with a strong core crushing plateau leads to a shield behaviour with strong bending. For the reference case, the shield is deformed mainly by indentation.

Effect of the front skin design

The results for the 25 simulations studying the front skin design are presented in Figure 9. First, we can see that nearly all the outputs vary less than that in Figure 8, which shows that, using these ranges of variation, the core design has more influence than the front skin design. This conclusion is coherent with the conclusions of [16], where the core plateau parameters were identified as the most influent parameters. Only the front skin maximum strain $ε_{f s}$ is more influenced by the front skin design.

It is clear from Figure 9 that the front skin thickness $t_{f s}$ and its yield strain $A_{f s}$ play very similar roles. In fact, a shield with a thick skin and a low yield strain seems to behave similarly to a shield with a thin skin and a high yield strain. The only outlier design is the shield with $t_{f s} = 2 m m$ and $A_{f s} = 500 M P a$ , especially for the crushing output. When observing this finite element simulation in detail, it appears that this is due to the buckling of the front skin, under the compression stresses due to the sandwich bending. This front skin buckling does not appear on any other case simulated in this study. It may come from a thin skin submitted to high stresses, due to a high yield strain.

All the other designs show that a weak front skin (thin and with a low yield strain) increases the indentation and bending of the shield, while a strong front skin reduces it, with nearly no effect on the crushing behaviour. In fact, the front skin strength has an influence on the importance of the shield deformation, as can be seen in Figure 10, where the three shields shown are the two extreme cases and the reference.

As shown in Figure 9, it appears that the increase of the shield deformation due to a weaker front skin is associated with an increase of $P_{S u p p}$ (which is explained as in ‘Effect of the core design’ section by an increase of bending localizing the contact between the back skin and the support), but also with a decrease of $F_{S u p p}$ . This effect may come from the fact that, due to the important deformation of the shield, the impacting bird is slowed down more progressively by a weak front skin than by a strong one. Indeed, the maximum force transmitted to the support occurs at $t = 1.86 m s$ for the shield (a) in Figure 10, at $t = 1.62 m s$ for the shield (b) and at $t = 1.56 m s$ for the shield (c).

The increase of the shield deformation is also associated with an increase of the back skin backward deflection (clearly seen in Figure 10, shield (a)), and an increase of the front skin maximum strain. This effect is due to the strong indentation of the sandwich, localizing the front skin strain at its centre. We can note that the strong indentation also limits the capacity of the bird to flow radially. This phenomenon makes any front skin rupture dramatic since the loss of the front skin rigidity will induce a strong indentation, which will increase the front skin loading by limiting the bird radial flow.

Effect of the geometrical parameters

In the third 2D DOE, the influences of the core height $H_{c}$ and the support aperture size $L_{a p}$ are studied. The results of this DOE are given in Figure 11.

Here again, it appears that the variation ranges of the outputs are smaller than that for the core design study. Nevertheless, a clear influence of the geometrical parameters on indentation and bending can be seen in Figure 11.

First, it appears that a decrease of $H_{c}$ leads to a decrease of indentation and an increase of bending. The increase of bending is explained simply by the bending rigidity of a sandwich, which is directly influenced by its core height. This increase of bending is here associated with an increase of the backward deflection and the maximum pressure transmitted to the support, which is coherent with the phenomenon already identified.

Second, the increase of $L_{a p}$ seems to have the same effect as a decrease of $H_{c}$ on bending and $δ_{b s}$ . This is explained by the fact that this parameter represents the size of the unsupported area of the sandwich. The maximum pressure transmitted $P_{S u p p}$ seems to oscillate with $L_{a p}$ . This phenomenon does not have any physical justification and seems to be caused by the relative position of the support and back skin elements in the simulation, which change when the parameter $L_{a p}$ change. Similarly, no physical explanation has been found to justify the evolution of $F_{S u p p}$ .

Thus, it seems that the core height and the support aperture size have similar influences on the shield behaviour, with a stronger effect for low values of $H_{c}$ and high values of $L_{a p}$ . To illustrate this effect, the deformed shapes of three shields are represented in Figure 12. The shield (a) has a high core and a small aperture, the shield (b) is the reference and the shield (c) has a thin core and a large aperture.

Conclusion on the 2D studies

With these three ‘two factors at a time’ DOE, it was possible to explore in detail the influence of each of the six parameters chosen. The use of five levels per parameter allowed us to observe the non-linearity of the different effects, and studying the parameters in pairs showed that important interactions can exist between the different parameters.

These 2D studies confirmed that the most globally influent parameters are the core crushing and shearing plateaus, but they also showed that depending on the output studied, other parameters can be more influent. Different shield behaviours have been identified for different shield designs, and physical explanations for these different behaviours have been proposed. The links between the behaviour criterions (indentation, bending and crushing), and the target protection criterions ( $P_{S u p p}$ , $F_{S u p p}$ , $δ_{b s}$ , $ε_{f s}$ ) have been identified and explained when possible.

But all these conclusions have been obtained on a strongly reduced design space. In the full 6D space, only three 2D planes have been studied. Moreover, the results showed that interactions can exist between the different parameters. To study these interactions, and to assert the results obtained with these 2D studies, it is necessary to conduct a full 6D study, in which all the parameters are varied together.

Expanding the parametric study in 6D

A 6D study is far more complex to conduct than a 2D study. First, it is impossible to use a simple full factorial DOE, since the number of simulations needed increase exponentially (five levels for six parameters lead to $5^{6} = 15,625$ simulations). Second, the visualisation of the results is not possible in full 6D and can only be done using 2D slices, which makes the analysis of results complicated. It is then necessary to use statistical methods to measure the influence of the different parameters.

To handle these issues, the usual approach is to construct and train a surrogate model. This model is a mathematical function and its goal is to approximate the output studied everywhere in the design space and not only on the tested points. This approach is widely used in the field of engineering design, and described in detail in [19] and [20]. The main steps of this approach are as follows:

First, the choice of the design points to simulate in order to train the surrogate model. The model will interpolate the output between these known points, and thus the choice of this design of experiment is crucial.

Second, the construction, training and validation of the surrogate model. It is necessary to choose the form of the mathematical functions which will be used to construct the surrogate model, as it will constrain the precision of the model. The validation consists to check if the model is able to approximate efficiently the output. If not, it is possible either to change the form of the surrogate model or to return to the first step to add more known points.

Third, the analysis of the model obtained, in order to measure the effects of each parameter and their interactions, and to conduct physical interpretations. Global methods of analysis will be used to validate (or not) the observations made on the 2D studies in ‘2D parametric studies’ section.

Choice of the 6D DOEs

Usually, the choice of the DOE depends strongly on the form of the surrogate model to train, which is chosen depending on the a priori available knowledge of the output. Here, the 2D studies showed that the output behaviours can be strongly non-linear and with interactions. It is then necessary to have a surrogate model able to adapt to a wide range of behaviour. Moreover, in this work, seven different outputs are studied and thus seven different surrogate models will be constructed. As each simulation allows measuring the seven outputs, all seven surrogate models will be trained on the same set of known points. For these reasons, the DOE used needs to be as general as possible, to adapt to the different situations.

A lot of different methods exist to create a design of experiments [19,21,22]. We choose here to use a Latin Hypercube Sampling (LHS) strategy, because it allows to choose a priori the number of design points. In an LHS DOEs with $N_{s}$ points, each parameter has $N_{s}$ levels and takes each level once and only once. This type of DOE is of particular interest for computer experiments since no point is tested twice, and this property holds true if any parameter appears irrelevant. Constructing an LHS is simple (by random permutation of the rows of an identity matrix), but this approach can produce designs with very different uniformity. In order to obtain a good LHS, the original one can be optimized to increase the uniformity. For example, in [20], a strategy using the maximin [23] criterion is presented.

In this work, 73 sampling points have already been simulated during the 2D studies presented in ‘2D parametric studies’ section. In order to choose 100 new design points to be simulated, we adapted the strategy presented in [20] to optimize the maximin criterion of the complete 173 points DOE. The 100 new points chosen are then simulated and the seven outputs are measured for each case. The 173 known points can then be used to create the surrogate models.

Construction and validation of the surrogate models

As seven models have to be constructed, it would be possible to choose different frameworks for each, but for practical reasons we choose to create the seven surrogate models using the same framework. A lot of different frameworks of surrogate models exist [19,20,22], and we choose here to create our surrogate models using the Gaussian Process (GP) framework, described in detail in the book [24]. All the computations for the creation and analysis of these surrogate models are made on MATLAB using the GPML toolbox available with the book.

GP are used because they can reproduce many different behaviours. Moreover, their statistical framework allows us to estimate the approximation error at each point of the design space, which can be very useful to estimate the precision of the surrogate models, and possibly improve them by adding new simulation points.

Thus, for each of the seven output, a surrogate model is constructed using a Matérn ν = 3/2 kernel with Automated Relevance Determination, and a constant mean [24]. The nine hyper-parameters (the mean, the signal variance, the noise variance and the six length scales) are optimised using the log marginal likelihood maximisation function provided in the GPML toolbox [24]. In order to ensure a good convergence of the hyper-parameters optimisation, the data are first normalized using MATLAB zscore function, to set their means to zero and their variances to one in each direction.

The obtained surrogate models are then tested using the Leave-One-Out (LOO) method. This method consists of training the surrogate model on all the data points minus one, and then of testing it on the last data point. By repeating this process on each data point, a mean prediction error can be computed. This mean error is presented for each output in Figure 13, in percentage of the variation range of the output.

We can see from Figure 13 that the mean LOO error is lower than 5% for all the outputs studied. This shows that the surrogate models constructed are able to represent globally the evolution of the seven outputs. Thus, these models will be used to analyse the influence of each parameter.

Analysis of the full 6D design space

Methods used

In order to analyse a six dimensional space, we choose to use two complementary statistical methods: the Morris method [25] and the Sobol′ method [26].

The Morris method consists in measuring the elementary effect of each parameter around an initial point, using a one factor at a time DOEs. By using multiple initial points, evenly distributed in space, it is then possible to compute for each parameter $k$ , a mean effect $\bar{d_{k}}$ and its variance $S_{k}$ . Then, the greater $\bar{d_{k}}$ is, the more the parameter $k$ is influent, and the greater $S_{k}$ is, the more this influence is non-linear or with interactions with other parameters.

The Sobol′ method allows to compute the parameters’ sensitivity indices and their total sensitivity indices. The sensitivity indice represents the fraction of the total variance of the output contributed by the parameter individually, while the total sensitivity indice represents the influence of the parameter and all its interaction with other parameters. To estimate these indices, we use the MATLAB toolbox GSAT [27].

Thus, the Morris and Sobol′ methods both estimate the influence of each parameter, but in slightly different ways: the Morris method allows measuring the linear effect of a parameter and the variance of this effect, without being able to distinguish between non-linearity and interaction with other parameters, while the Sobol′ method measure the total influence of a parameter isolated and with interactions, without measuring the non-linearity of this influence. In this work, we used both approaches in order to obtain a maximum of knowledge about the output behaviour.

For both these approaches, it is necessary to use the unit hypercube as design space. Thus, the parameter space is normalized appropriately for all the surrogate models before any computation.

Results and interpretations

The results given by the Morris method are shown in Figure 14. For each of the seven surrogate models, we used 100,000 initial points pseudo-randomly distributed to estimate the mean effects $\bar{d_{k}}$ and their variance $S_{k}$ . In order to compare the influence of each parameter, they are placed in the $(|\bar{d}|, S)$ space for each output in Figure 14, where the dotted diagonal represents the line $|\bar{d}| = S$ . To visualize the sign of the mean effects, the parameters are represented by a • when the mean effect is positive and by a * when it is negative. A positive/negative effect means that the output increases/decreases when the parameters increases.

The results of the Sobol′ analysis are shown in Figure 15. For each surrogate model and each parameter, the corresponding sensitivity index and total sensitivity index are presented as superposed bars of different widths. Thus, the difference of height between the two bars of one parameter represents the sum of all the interactions with other parameters.

First, we can see from Figure 14 that the parameters are mainly placed near or above the $|\bar{d}| = S$ line, which shows important non-linearities or interactions. On the other hand, Figure 15 shows that the interactions between parameters are not so important. Thus, we can conclude that the effects of the different parameters are strongly non-linear. This conclusion is coherent with the observations in ‘2D parametric studies’ section.

By observing these two figures, it is clear that the more influential parameters are the core properties and then the core height. This confirms the parameter ranking obtained in [16] and the output ranges of variations observed in ‘2D parametric studies’ section. This is particularly true for two behaviour outputs (indentation and crushing), where the core out-of-plane crushing plateau $σ_{p c}$ seems to be the only influent parameter. This can be seen more clearly in Figure 16, where the 173 simulations are represented in the (indentation, $σ_{p c}$ ) and the (crushing, $σ_{p c}$ ) planes. There is a clear correlation between the output value and the $σ_{p c}$ parameter value.

Figure 16 also shows the strong non-linearity of the effect of $σ_{p c}$ , already measured by the Morris method, as can be seen in Figure 14 where the $σ_{p c}$ point has a strong value of $S$ . We can also see in Figure 16 the strong threshold effect of $σ_{p c}$ on the crushing output, already identified in ‘Effect of the core design’ section.

The Morris method also enables us to confirm the signs of the effects of the different parameters. For example, the effect of the core height $H_{c}$ is negative on bending and backward deflection $δ_{b s}$ , while positive on front skin maximum strain $ε_{f s}$ , which is coherent with the local observations of ‘Effect of the geometrical parameters’ section (Figure 12). These signs also show that the two core material parameters $σ_{p c}$ and $τ_{p c}$ have always opposed effects, which confirms that the two ‘extreme’ designs of core material are either a core with high crushing plateau and low shearing plateau or a core with low crushing plateau and high shearing plateau.

By observing Figure 15 and Figure 14, we can also see that the two front skin parameters (its thickness $t_{f s}$ and yield stress $A_{f s}$ ) have effects which are quite similar, and with same sign, on all the outputs. Once again, it confirms the local observations of ‘Effect of the front skin design’ section.

Conclusions on the full 6D study

In this section, the full 6D design space was studied using a set of 100 simulations spread evenly using a Latin Hyper-square DOE. Using these simulations and the 73 simulations already done in ‘2D parametric studies’ section, one surrogate model was constructed for each output studied.

After validation, it was possible to use these models to analyse the design space using two different methods: the Morris method and the Sobol′ method. We saw that these two methods are complementary and that using the two in parallel allowed us to measure both the non-linearity of the parameters effects and the interactions between parameters.

Using the results obtained, we saw that, on the studied design space, the parameters have effects which are often strongly non-linear, but with only small interactions. Then, the results were compared to the local observations of ‘2D parametric studies’ section. Both the ranking of parameters with respect of their influence and the signs of the parameter effects were coherent with local observations.

Thus, this 6D study allowed us to confirm globally the different behaviours observed in ‘2D parametric studies’ section. But it is important to note that, due to the strong non-linearity of the parameters effects, these behaviours can be false locally.

General conclusion

To tackle the problem of bird impact on sandwich shield (numerous parameters, time-dependant response, etc.), DOE and machine learning methods were used in this study with finite element simulations. After an initial screening to reduce the number of parameters to take into account, GP were used to create surrogate models in a 6D domain, using smart sampling. This method enables to analyse the effects of the six most influential parameters and their interactions on the shield behaviour (deformed shape and protection criterions) within a reasonable computation time.

The main conclusions of this study are the following:

The two most influential design parameters are clearly the core out-of-plane crushing plateau and the core shearing plateau. This means that during the design of a shield for soft impact, great care should be taken in the choice of the core material. We also saw that these two parameters have opposite effects . Thus, a shield with higher crushing plateau and lower shearing plateau will transmit more stresses (and more localized) to the back skin and the support, while reducing the area of core crushed and increasing backward deflection. Yet, when changing the core material, designers usually stay in the same category of cellular material (e.g. aluminium honeycombs), which means that changing both plateau parameters in the same direction (either increasing or decreasing).

The two skin parameters (thickness and material yield stress) have very similar effects. A stronger skin (i.e. higher thickness or higher yield stress) will tend to reduce the shield indentation and flexion, thus spreading more evenly the core crushing and reducing backward deflexion and front skin maximum strain. On the other hand, these two parameters have very different effects on the shield mass. For a lightweight shield, it is then more advantageous to increase the front skin yield strain while reducing its thickness. This conclusion shows that a soft impact problem is very different from a hard impact problem, where a softer front skin is usually beneficial [28]. This difference may come from the fact that, for a soft impact case, the strain gradients are smaller. This implies that the deformations inside the shield are more global and it appears that all the parts of the shield are working together at the same time and not one after the other (as can be the case in hard impact). Thus, the sandwich shield has to be studied as a whole and not as a succession of different layers.

The core height has a strong effect on the shield bending and backward deflection, which is coherent with usual sandwich behaviour. Thus, a higher core is usually beneficial for target protection.

For this work, we deliberately chose a quite simple finite element model in order to simulate numerous possible designs, and the conclusions have to be taken with precautions. In particular, the core material model used here is decoupled, which means that the real interactions between the two core plateau parameters are probably more important than calculated. Another limitation of this study is the fact that we did not take into account the possible front skin rupture. In the 173 configurations simulated, none showed front skin strains greater than 15%, which shows a posteriori that this hypothesis was valid for the material used. Nonetheless, we have to keep in mind that the front skin rupture would drastically change the shield behaviour, as losing the front skin rigidity would mean losing the ‘sandwich’ behaviour.

At last, thanks to this new understanding, general rules have been proposed to orient the design of better shields with similar boundary conditions. Moreover, the surrogate models created here are a first step towards a design tool to help engineers, as they can be used in an optimisation loop in order to find an optimal shield with only a few more simulations.

Of course, such a numerical study requires validation, and an experimental campaign is currently being conducted to validate the main conclusions.

Footnotes

Acknowledgements

The authors gratefully acknowledge the members of this project: Stelia, Airbus Group, Cedrem, Esteve, Ateca, Nimitech, Hutchinson, I2M and IRT Saint Exupéry. This work has benefited from access to the HPC resources of CALMIP (Calcul en Midi-Pyrénées).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors want to thank BPIFrance and the Région Occitanie for their financial support through the FUI project SAMBA (Shock Absorber Material for Bird-shield Application).

References

Dennis

Lyle

Bird strike damage & windshield bird strike final report. Farnham, UK: European Aviation Safety Agency, 2008.

CS-25 – Certification specifications and acceptable means of compliance for large aeroplanes, Section CS-25.631 – Bird strike damage. European Aviation Safety Agency, 2003.

14 CFR Part 25 - Airworthiness standards: Transport category airplanes, Section 25.631 - Bird strike damage. Federal Aviation Administration, Dept. of transportation, 2003.

Gibson

Ashby

MF.

Cellular solids: structure and properties. Cambridge: Cambridge University Press, 1999.

Barber

Taylor

Wilbeck

JS.

Characterization of bird impacts on a rigid plate: part 1. Air Force Flight Dynamics Laboratory, Technical Report AFFDL-TR-75-5, Dayton, US, January 1975.

Barber

Taylor

Wilbeck

JS.

Bird impact forces and pressures on rigid and compliant targets. Air Force Flight Dynamics Laboratory, Technical Report AFFDL-TR-77-60, Dayton, US, May 1978.

Wilbeck

JS.

Impact behavior of low strength projectiles. Wright-Patterson Air Force Base, Technical Report AFML-TR-77-134, October 1978.

Heimbs

Computational methods for bird strike simulations: a review. Comput Struct 2011; 89: 2093–2112.

Airoldi

Cacchione

Modelling of impact forces and pressures in Lagrangian bird strike analyses. Int J Impact Eng 2006; 32: 1651–1677.

10.

Hedayati

Ziaei-Rad

A new bird model and the effect of bird geometry in impacts from various orientations. Aerosp

Sci Technol 2013; 28: 9–20.

11.

Hanssen

Girard

Olovsson

, et al. A numerical model for bird strike of aluminium foam-based sandwich panels. Int J Impact Eng 2006; 32: 1127–1144.

12.

Hohe

Hardenacke

Fascio

, et al. Numerical and experimental design of graded cellular sandwich cores for multi-functional aerospace applications. Mater Des 2012; 39: 20–32.

13.

Liu

Shi

, et al. Dynamic response of bird strike on aluminum honeycomb-based sandwich panels. J Aerosp Eng 2012; 27: 520–528.

14.

Hedayati

Sadighi

Effect of using an inner plate between two faces of a sandwich structure in resistance to bird-strike impact. J Aerosp Eng 2016; 29; 04015020.

15.

Liu

Gao

, et al. Dynamic response of bird strike on aluminium foam-based sandwich panels. Int J Crashworthiness 2015; 20: 325–336.

16.

Wilhelm

Rivallant

Ferrero

J-F.

Study of the deformation of a sandwich shield subjected to bird impact: a behaviour analysis tool using vector decomposition. J Sandw Struct Mater 2019; 21: 689–706.

17.

Johnson

Cook

WH.

A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In: Proceedings of the 7th international symposium on ballistics. The Netherlands, 19–21 April 1983, pp.541–547.

18.

Kolopp

Alvarado

Rivallant

, et al. Modeling impact on aluminium sandwich including velocity effects in honeycomb core. J Sandw Struct Mater 2013; 15: 733–757.

19.

Queipo

Haftka

Shyy

, et al. Surrogate-based analysis and optimization. Prog Aerosp Sci 2005; 41: 1–28.

20.

Forrester

Sobester

Keane

Engineering design via surrogate modelling: a practical guide. John Wiley & Sons, 2008.

21.

Goupy

Creighton

Introduction aux plans d’expériences. 3ème éd. Paris: Dunod, 2006.

22.

Simpson

Lin

Chen

Sampling strategies for computer experiments: design and analysis. Int J Reliab Appl 2001; 2: 209–240.

23.

Johnson

Moore

Ylvisaker

Minimax and maximin distance designs. J Stat Plan Inference 1990; 26: 131–148.

24.

Rasmussen

Williams

CKI.

Gaussian processes for machine learning. Cambridge, MA: MIT Press, 2006.

25.

Morris

MD.

Factorial sampling plans for preliminary computational experiments. Technometrics 1991; 33: 161–174.

26.

Sobol′

IM.

Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comput Simul Second IMACS Semin Monte Carlo Meth 2001; 55: 271–280.

27.

Cannavó

Sensitivity analysis for volcanic source modeling quality assessment and model selection. Comput Geosci 2012; 44: 52–59.

28.

Kolopp

Rivallant

Bouvet

Experimental study of sandwich structures as armour against medium-velocity impacts. Int J Impact Eng 2013b; 61: 24–35.