Abstract
Bird impact poses serious threats to military and civilian aircrafts as they lead to fatal structural damage to critical aircraft components. The exposed aircraft components such as windshields, radomes, leading edges, engine structure, and blades are vulnerable to bird strikes. Windshield is the frontal part of cockpit and more susceptible to bird impact. In the present study, finite element (FE) simulations were performed to assess the dynamic response of windshield against high velocity bird impact. Numerical simulations were performed by developing nonlinear FE model in commercially available explicit FE solver AUTODYN. An elastic-plastic material model coupled with maximum principal strain failure criterion was implemented to model the impact response of windshield. Numerical model was validated with published experimental results and further employed to investigate the influence of various parameters on dynamic behavior of windshield. The parameters include the mass, shape, and velocity of bird, angle of impact, and impact location. On the basis of numerical results, the critical bird velocity and failure locations on windshield were also determined. The results show that these parameters have strong influence on impact response of windshield, and bird velocity and impact angle were amongst the most critical factors to be considered in windshield design.
1. Introduction
High velocity bird impact is one of the most significant hazards to both civilian and military aircrafts [1]. When fighter military aircrafts are operated at lower altitude and higher speed, the probability of bird impact increases and proves lethal to the safety of pilot and critical aircraft components. Windshield is the exposed part of aircraft and prone to bird impact. In order to ensure the safety of aircraft, the windshield must be capable of withstanding high velocity impact threats. Design of an optimum impact resistant windshield is a challenge and requires extensive experimental testing. The advanced numerical techniques are being adopted as an effective tool to simulate the bird strike event and provide a substitute to excessive costly experimentations. Moreover, it allows analyzing the most stringent impact conditions that cannot be considered in the experiments and provides detailed insight to the impact process which is difficult to observe during experimental testing.
Several researchers such as Zang et al. [2], Samuelson and Sornas [3], and Boroughs [4] carried out finite element analysis to investigate the impact response of windshield against bird strike. McCarty et al. [5, 6] used Materially and Geometrically Nonlinear Analysis (MAGNA), a nonlinear finite element analysis program for designing of windshield and canopy of military aircrafts. Wang et al. [7–9] simulated the failure of aircraft windshield against bird impact by using a modified nonlinear viscoelastic constitutive model together with Zhu-Wang-Tang (ZWT) damage model for the PMMA-based windshield. The constitutive model and its failure criterion effectively predicted the failure of windshield under a range of impact velocities. The authors also examined different critical factors which affect the impact response of aircraft windshield against bird strike. The results proposed that material model for windshield, boundary conditions, mesh density, surrounding structure of windshield, and bird velocity are the critical factors that must be taken into account for FE analysis of aircraft anti-bird design. Guida et al. [10] numerically examined the influence of geometric parameters of windshield and impact parameters under high speed bird impact by using explicit FE code. The results showed that windshield panel dimensions, thickness, and curvature as well as bird velocity, size, and impact angle considerably affect the impact response of windshield. Liu et al. [11] applied smooth-particle-hydrodynamics (SPH) based approach to model the bird against impact on windshield structure by using explicit FE solver. The maximum displacements of various points on camber line of windshield were compared with experimental results and were found to be in good agreement. The results also indicated that SPH solver gives better comparison than Lagrangian solver. Zhu et al. [12] numerically studied the bird impact response on windshield by employing user-defined material subroutine in explicit FE solver. Numerical results of failure modes of the windshield, deformation, displacement, and strain curves of the measured points on the windshield agreed well with the experimental results. Yang et al. [13] carried out FE and experimental investigation to study the structure strength of windshield subjected to bird impact. Based on numerical scheme, the critical impact velocities on certain locations on windshield are determined and verified through experimental results.
The present work aims to numerically model the impact response of windshield and examine the effect of various factors that contribute to optimize the design of certain windshield. Several critical factors such as mass and shape of bird, impact velocity, angle of impact, and impact location were studied and their influence on the dynamic response of windshield was assessed. The numerical model was developed and implemented in explicit finite element hydro code ANSYS AUTODYN.
2. Numerical Modeling
2.1. Windshield Material Model
The windshield considered in this study is monolithic, uniform cross-section, PMMA-based aviation organic glass. The dynamic properties of this material have been extensively studied in the literature and several material models are available for modeling the impact response of this material [7, 14–17]. In this study, the elastoplastic material model along with maximum principal strain failure criterion was defined to predict damage and failure of windshield. The material model was implemented by incorporating isotropic hardening using Von Mises yield criterion together with rate-dependent Cowper Symonds plasticity law. Cowper Symonds strength model defines the yield strength of isotropic strain hardening of strain rate dependent material [10]. The yield surface can be defined as where is yield stress of the material, is yield stress at zero plastic strain, is strain hardening coefficient, is strain hardening exponent, and and are strain rate hardening coefficients. Von Mises is simple and convenient criterion to apply as it defines a smooth and continuous yield surface with good approximation at high stresses. At given principal stresses ,, and , the yield criterion is defined as
The maximum principal strain criterion implies that if the maximum tensile principal strain exceeds the prescribed limits, then material will instantaneously fail. Failure is predicted when either of the principal strains or , resulting from the principal stresses or , equals or exceeds the maximum strain corresponding to the yield strength of the material in uniaxial tension or compression. For yielding in tension the minimum principle strain would equal the yield strain in uniaxial tension. If the strains are expressed in terms of stress, then
The material properties of windshield used in the analysis are given in Table 1.
2.2. Bird Material Model
The actual bird is combination of flesh, blood, and bones, and it is difficult to implement the actual bird constitutive model in numerical program. Numerous researchers used different approaches to closely approximate the material response of bird. Some authors modeled the bird with elastoplastic material law along with certain failure criteria [12, 18–20], while others used equation of state approach for constitutive modeling with pressure volume behavior of water [21–24]. Bird is mostly composed of water. It behaves like a soft body and acts as a fluid on the structure during the impact. Water like hydrodynamic response by using Mie-Gruneisen equation of state (EOS) with negligible strength effects was implemented to model bird in this study. Mie-Gruneisen EOS correlates the material volumetric strength and pressure to density ratio, and it is easy to establish Gruneisen equation based on the shock Hugoniot [25] where where andare the initial and instantaneous densities of material, is the intercept, and is linear Hugoniot slop of shock velocity and particle velocity relationship. Gruneisen equation describes a linear relation between shock and particle velocity, where
Gruneisen form of equation of state based on shock Hugoniot is
It is assumed that constant and
For , this formulation gives a limiting value of compression because the pressure approaches to infinity as the term in (4) becomes zero and the pressure therefore becomes infinite and gives maximum density of . Also long before this regime is approached, the assumption of is not valid. Moreover, the assumption of a linear relationship between the shock velocity and the particle velocity does not hold for too large compression. And at high shock strengths, some nonlinearity in this relationship is apparent. To incorporate the nonlinearity in numerical model, two linear fits to the shock velocity and particle velocity relationship were established: one holding at low shock compressions defined by and other at high shock compressions where . The region between and is covered by a smooth interpolation between the two linear relationships [25]. The equation of state has been further improved to include a quadratic shock and particle velocity relationship where and are coefficients of slop, is Gruneisen parameter, and is internal energy. In the present numerical scheme, the water Gruneisen EOS parameters ,, and were used to model bird as soft projectile and other parameters were set to zero [25]. The density of the bird was taken as 900 kg/m3 as it has been taken by most of the researchers before.
2.3. Validation of Numerical Model
2.3.1. Model for Windshield
In order to validate the numerical model for windshield material, FE simulations of simple uniaxial tensile test were performed. A standard dumbbell-shaped tensile test specimen with dimensions of 200 × 20 × 2.67 mm and gauge length of 50.8 mm was built. The specimen was meshed with Lagrangian solid elements and the material properties were taken as described in Section 2.1. One end of the specimen was fixed while constant velocity load was applied at the other end. The rate of loading was adjusted according to rate of change of strain. The results from FE simulations were compared with available results [17] in terms of stress strain curves at different strain rates as shown in Figure 1. The results show that the model predicts the stress-strain behavior of PMMA at different strain rates with fair accuracy.
2.3.2. EOS for Bird
To validate the Gruneisen EOS implemented in FE solver, the Hugoniot and stagnation pressures were compared with theoretically obtained pressure values and Wilbeck [26] experimental results. Wilbeck experimental results carry significant importance in bird strike analysis, and many researchers used his experimental data for comparison of numerical results. For theoretical calculations, one-dimensional Hugoniot analysis was carried out in which bird impact is characterized in two stages. The first stage is initial shock called Hugoniot pressure which gives the maximum possible value of pressure during impact, and the other stage is steady state flow called stagnation pressure in which pressure stabilizes with time. The maximum pressure, that is, Hugoniot pressure, can be estimated by implying hydrodynamic impact theory and equations of conservation mass and momentum where ,,, and are the initial and instantaneous density and pressure values. is shock propagating velocity and is particle velocity behind shock. The particle velocity is assumed to be equal to the initial impact velocity of bird. The Hugoniot analytical pressure can be defined as
The shock velocity is function of initial impact velocity and can be obtained by solving the nonlinear equation [24]
where is velocity of sound in the material, is porosity of material (for 10% porosity ), and value of is determined experimentally. The variation of shock velocity with initial impact velocity is shown in Figure 2.
The steady flow pressure can be estimated by using Bernoulli relationship
For FE validation, a square rigid plate of 0.5 m side was modeled in Lagrangian grid fully constrained at sides. The bird is modeled as right cylinder of 0.18 m length and 0.06 m radius with Lagrangian solid elements and material properties described in Section 2.2. In the analysis, bird with initial impact velocity of 116 m/s was impacted against rigid plate, and values for Hugoniot and stagnation pressures at central point of impact were obtained as shown in Figure 3. In order to compare the results with experimental data, the velocity of bird was taken as 116 m/s used by Wilbeck in his experiments. The results were normalized by dividing pressure with stagnation pressure and time with total impact duration and compared with experimental and analytical results.
From analysis, the Hugoniot pressure was predicted to have a maximum value of 102.5 MPa and the stagnation pressure is 6.5 MPa, giving normalized values of 15.7 and 1, respectively. The analytical values of Hugoniot pressure and stagnation pressure were calculated as 86.3 MPa and 6 MPa, which after normalizing gives 14.3 and 1, respectively. The comparison of FE, analytical, and experimental values of normalized pressure is shown in Figure 4. The trend of the plot is consistent with experimental results where a sudden peak pressure value was observed at initial shock and the pressure then stabilized with time. The duration of pressure decay was also in accordance with experimental results. However, the FE simulation predicted peak pressure is higher than the experimental results. The overprediction in pressure peak can be credited to bird orientation during impact, because in real bird impact experiments the initial pressure peak is highly dependent on orientation of the bird as described by Hedayati and Ziaei-Rad [24] in their work. Moreover, the Hugoniot peak pressure occurs in a very short time about 4-5 μs and the capturing of the accurate pressure peak is a challenging experimental task. The modern high frequency state of the art pressure transducers must be employed to record the accurate pressure peak during impact event, and latest experimental data in this regard is required for better comparison of FE and experimental results.
However, for the present work, the bird model was further verified by impacting the bird on a flexible metallic plate and the resultant plastic deformation after impact was determined. Welsh and Centonze [27] experimental results were considered for this purpose in which the bird with initial velocity of 146 m/s was impacted on 6.35 mm thick T6061-T6 aluminum plate and corresponding plastic deformation of the plate was measured. FE simulations were performed in accordance with experimental parameters and the residual plastic deformation of the plate after impact was determined as shown in Figure 5. The FE predicted deformation of 43.68 mm was in close agreement with experimentally determined value of 41.275 mm. The validated material models for bird and windshield were then used to build a full scale model for bird-windshield impact problem.
2.3.3. 3D FE Model of Impact Problem
The windshield and bird were modeled with solid elements by using Lagrangian grid as shown in Figure 6. The windshield consists of 14,400 elements with 60 elements along the curvature, 40 elements along sides, and 4 through thickness elements. More refined elements distribution is adopted around the area of impact as most of the deformation takes place at this particular impact region. The 1.8 Kg bird was modeled as right cylinder of 60 mm radius and 180 mm length. The bird was modeled as soft body with 10 elements across radius and 20 elements along length of cylinder.
The elements in the mesh undergo severe distortion due to high rate of deformation. These distorted elements must be removed from the mesh in order to run program smoothly and to avoid reduced time step problem. The deletion of distorted elements is carried out by using geometric strain erosion model provided in AUTODYN. During the impact process, the bird highly deforms, perishes, and loses most of its mass due to elements removal process which causes inaccurate results of numerical calculations. This problem was overcome by employing the option of retaining inertia of eroded nodes in which the solver ascribes the removed cell mass to their nodal points. The contact between windshield and bird was defined by using Lagrange/Lagrange interaction option of AUTODYN. In this option, an automatic detection zone is defined around each interacting Lagrangian subgrid (independent or with itself). When a node enters into this zone, it is automatically repelled. The edges of the windshield are fully constrained to provide fixed boundary conditions. Four gage points 1, 2, 3, and 4 (Figure 6) were marked on the windshield central line to record the values of displacement, stress, and strain and to compare them with experimental results.
2.3.4. Experimental Confirmation
In the experimental results [11, 19], 1.8 Kg bird with velocity of 64.4 m/sec was impacted on windshield at four different locations (1, 2, 3, and 4). The displacement profiles for all the four locations were recorded during experiments. Numerical simulations were performed, and results were compared with experimental data and found in good agreement which prove the validity of numerical model as shown in Figure 7. This verified that numerical model was then further employed to investigate the effect of various parameters on windshield behavior.
(a)
(b)
(c)
(d)
3. Simulations
Numerical simulations were carried out to predict the dynamic response of windshield for range of impact velocities, impact angle and location of impact, and bird mass and shape. Figure 8 depicts various modes of deformation of windshield at different time intervals when impacted at location 2 with velocity of 64.4 m/sec. The windshield remains in the elastic state and keeps its original shape after impact. No sign of damage or failure was observed at this impact velocity. Also, at this velocity, the bird fails partially and slides along the surface of windshield.
When the velocity of bird is increased, the stresses escalate due to higher impact force and windshield tends to deform plastically. With further increase in impact velocity windshield shield reaches its elastic limit and suffers from permanent deformation leading to its complete failure. The equivalent stresses along cross-section of windshield for different velocities are shown in Figure 9. At higher velocity, the bird deforms severely and gets fragmented during sliding along the surface of windshield. Two different orientations of bird, that is, 15° and −15° from its axis, were studied to see the effect of impact angle on response of windshield and shown in Figure 10. The change of impact angle showed significant effect on normal displacement, stress, strain, and impact force. For 15° impact angle, partial bird failure occurs and whole bird body slides along the windshield surface in its way of impact.
At −15° angle, the bird transfers most of its energy to windshield giving higher values displacement and stress at point of impact. Most of the body of bird fails during this event.
Impact response for two additional bird geometries (hemispherical and ellipsoidal of similar length to diameter ratio) was also simulated and shown in Figure 11. For all bird shapes, the maximum length to diameter ratio and mass of the bird were taken constant. The simulations results show that impact The simulations results show that impact due to cylindrical-shaped bird produces slightly higher deformation on the windshield due to its increased contact area.
4. Results and Discussion
The impact response of windshield for various impact velocities was considered. It was noted that normal displacement at all gage locations increases with the increase of velocity. The time to reach maximum displacement at locations 3 and 4 increases as they lie farther from point of impact (location 2). The displacement time plots for different impact velocities at locations 2 and 4 are shown in Figures 12 and 13. Increase in velocity caused more deformation at the upper half of windshield because more of the bird mass slide and transferred more energy to the upper end. An instantaneous failure at the point of impact occurs when velocity was increased to 200 m/s.
Figures 14 and 15 show the plots for equivalent stress and corresponding strain at impact point. The amplitude of stress increases with the increase of velocity; at 64.4 m/s the value of maximum stress is 30 MPa which rises up to 50 MPa at 75 m/s velocity. At 100 m/s impact velocity, the stresses in the windshield become higher than the yield stress. On further increasing the velocity to 125 m/s, the ultimate stress limit is reached and upper end of windshield fails. At 200 m/s impact velocity, the windshield failed instantly at the point of impact in 1.9 ms.
Therefore, at velocity of 100 m/s and higher, the plastic deformation starts to prevail in windshield and it remains in the state of high stresses which also defines the critical impact velocity of the particular windshield. The modes of deformation and initiation of plasticity are shown in Figure 16. It can be seen that when the velocity is 100 m/s, sign of plastic deformation appears around the point of impact and upper end of windshield. The plastic deformation increases and windshield starts to fail at its upper end when the velocity is increased to 150 m/s. The further increase in velocity leads to major windshield failure at upper end and point of impact.
The influence of three different bird angles at the point of impact was examined. Figure 17 shows the plots for normal displacement at 64.4 m/s impact velocity from which it can be observed that at 15° angle of impact the maximum normal displacement is 3.5 mm which increases to 13 mm for 0° and 23 mm for −15° angle. The impact force also increases sharply with the change of impact angle as shown in Figure 18. At 15° angle, the maximum impact force of 9.9 kN is recorded at 5.6 ms which rises to 48.6 kN at 2.7 ms for −15° impact angle. Almost 5 times increase in peak impact force was observed due to change in impact angle from 15° to −15°.
The effect of impact angle on stress and strain history is shown in Figure 19. The maximum stress amplitude of 11.8 MPa and corresponding strain value of 0.00352 occur at 5.1 ms for 15° angle of impact. When the impact angle changes to −15°, the peak stress and strain value rise to 67.8 MPa and 0.0201 at 3.75 ms representing most severe impact conditions at 64.4 m/s velocity and windshield approaching its yield stress limit. Hence, impact angle is a critical parameter in determining the critical impact velocity for windshield.
(a)
(b)
The effect of bird mass on impact response was very obvious and shown in Figure 20. With the increase in bird mass, the corresponding value of displacement and stress increases because more kinetic energy of bird is transferred to windshield. This also shows that the critical impact velocity will be different for different bird masses to produce same deformation in the windshield.
(a)
(b)
The response of windshield impacted by three different shape birds was studied in this work. Similar displacement trends were observed for hemispherical- and ellipsoidal-shaped birds while cylindrical-shaped bird produces higher displacement as shown in Figure 21(a). The peak values of equivalent stress for hemispherical and ellipsoidal shapes are same and higher than cylindrical shape (Figure 21(b)).
(a)
(b)
5. Conclusions
The behavior of windshield against high speed bird impact was successfully simulated and the effect of various parameters on its dynamic response was studied. Bird impact velocity was among the most critical parameter having a strong influence on dynamic response. With the increase in bird velocity, windshield tends to deform plastically beyond its yield strength which finally leads to its major failure. For a range of velocities simulated in this study, there was a limiting impact velocity at which windshield suffers permanent plastic deformation and vulnerable to fail at certain crucial locations. The rearward fixed part of windshield was considered weakest at the critical velocity. For different impact angles, the response of windshield differs greatly. With less steep angle, most of the bird slides along the surface of windshield causing less damage. Steeper angle on the other hand produces high deformation and a plastic dimple is observed at the point of impact. The bird with higher mass proved more fatal to the windshield as they impact more kinetic energy to the structure. Although the shape of the bird did not show significant effect in this study, however, the bird with smaller length to diameter ratio and higher instantaneous contact area can affect the shock pressure and peak stress level in the structure. These critical factors can be parameterized together to predict the combined effect on impact response of windshield and can provide certain guiding principles for windshield design and optimization.
Acknowledgments
This work is supported by 973 Program (2012CB025904), NPU Foundation for Fundamental Research (NPU-FFR-JC201236), and Shaanxi Provincial Natural Science Foundation (2012JQ1003).