Ballistic gelatin Lagrange Mooney-Rivlin material model as a substitute of bird in finite element bird strike case studies

* Corresponding author https://doi.org/10.1590/1679-78256215 Abstract Numerical simulation plays a crucial role in today’s aviation industry. Modern computers with the latest commercial FE codes have fueled the triumph of establishing simulation as a prerequisite of aerostructure part design from micro modeling of materials to large-scale structural analysis. In this current research paper, a Mooney-Rivlin material model of ballistic gelatin with Lagrange code is analyzed as a potential candidate for the computational bird model to simulate bird strike case studies. To investigate the practicability, the model is compared with other established Lagrange and SPH EOS models for both rigid and deformable body impacts along with experimental data found in the literature adopting explicit solver Ansys Autodyn. Despite some discrepancies found during rigid body impacts, deformable plate impacts confirm the robustness of the model with significantly faster computation time. Besides, the biggest criticism of the Lagrange model, mesh distortion problem as fluid during bird strike case studies is efficiently tackled by adopting the node erosion algorithm as an effective technique to solve Lagrange bird models without affecting the


INTRODUCTION
In impact engineering, a projectile is defined as a soft-body when the induced internal stresses during a collision phenomenon have greatly exceeded the projectile material strength but are significantly lower than the target material strength (Martin, 1990). For instance, when a bird colloids with an aircraft, the bird is considered as a soft body since the generated stress in the body is well above the bird material strength.
Based on the velocity of the projectile, a body impact event can be categorized into five different stages, namely: elastic, plastic, hydrodynamic, sonic, and explosive (Wilbeck, 1978). In case of the low velocity, when the developed stress in the projectile during an impact is much lower than the material strength, the projectile is considered as elastic. With the increased velocity, the projectile strength reaches the plastic regime and a further inclination of the velocity would cause the projectile to behave like fluid in the hydrodynamic regime, where the material density instead of the material strength controls the response of the projectile. Bird flows as fluid in this regime and does not bounce (Peterson and Barber, 1976).

Wilbeck's theory
The fundamental hydrodynamic theory of a soft body impacting a rigid target is laid by (Barber, Taylor & Wilbeck, 1974;Barber, Taylor & Wilbeck, 1978;Wilbeck and Barber, 1978). This consists of four main phases: i. initial shock (also known as Hugoniot pressure), ii. impact shock decay, iii. steady-state flow, and, iv. pressure decay, Fig. 1 (a) (Heimbs, 2011). A pressure profile is also associated of these events for better illustration of a soft-body impact, Fig. 1 (b) (Heimbs, 2011). According to the Wilbeck's theory, the flat-ended cylindrical projectile travels at initial velocity 0 u and impacts the rigid wall which develops a shock wave traveling back to the projectile body Fig. 1 (a) -i, with a velocity p u . This creates a sudden lose of pressure at the contact point and generates a release wave propagating perpendicular to the surface, Fig. 1 (a) -ii. After several reflections of these release waves, the material starts to flow steadily leading to constant pressure, Fig. 1 (a) -iii. Finally, the kinetic energy of the projectile is completely absorbed and the impact event terminates, Fig. 1 (a) -iv.  (Heimbs, 2011) Mathematical formulation of the Wilbeck's hydrodynamics theory applies the conservation of mass and momentum, which is defined as follow.
Furthermore, the total duration of the impact D t can be determined from:

Computational bird modeling
For bird-strike certification documented in FAA, regulations requires the aircraft manufacturers to present experimental data for safe flight operations (FAA, 2012). However, high costs involved in experiments and complex setups with lengthy preparation time often challenge the development of aircraft structural parts. Moreover, impacted structures during experiments become severely damaged and further reuse is not practicable (Plassard et al., 2015;Heimbs et al., 2018;Orlando et al, 2017).Therefore, commercial computer simulation codes (Abaqus, Ansys, LS-Dyna, MSC Dytran, etc.) are regularly used as alternatives of bird-strike experiments for early-stage research and development of various structural parts, materials, joints, etc. before final production.
In computational solid mechanics, fluid-structure interaction problems are solved adopting fluid as one of the four different models, namely, i. Lagrange, ii. Euler, iii. Arbitrary Lagrange Euler (ALE), and, iv. Smoothed Particle Hydrodynamics (SPH). For bird-strike related problems, almost 49% of researchers preferred the Lagrange model where 37% of researchers employed SPH to predict results (Heimbs, 2011). However, most of the recent works investigated the bird impact issues are carried out using the SPH method (Such as: Eren et al., 2017;Kim and Kim, 2017;Zhou et al., 2019b;Zhou et al., 2019c;Jang and Ahn, 2018;Zhang et al., 2018;Qiu et al., 2020;Yu et al., 2019). The main reason of adopting SPH is that it can predict fairly accurate results in comparison to experiments and the recent advancement of parallel computing has brought down the solution time significantly, which was previously the main obstacle of using this method.
SPH is mainly a meshless Lagrange method that was first introduced to solve astrophysics problems (Gingold and Monaghan, 1977). It has overcome some of the problems that the Lagrange method lacks behind. For instance, the high mesh distortion problem of the Lagrange method during impact often leads to energy error during computation and the solution terminates prematurely. Nonetheless, in Lagrange solver, highly distorted nodes can be eroded from the solution, which often results in well-predicted outcomes (Hedayati & Sadighi, 2015). Another advantage of using the Lagrange method is that it is computationally inexpensive, modeling is relatively simple and possesses well-defined contact algorithms (Hughes et al., 2013) where the SPH method requires special modeling consideration during the simulation and lacks sharp boundaries (Anghileri et al., 2005). Despite these contraries, both method allows for tracking of the material deformation and history (Grimaldi et al., 2013). For additional information, readers are referred to (Hedayati & Sadighi, 2015).

Constitutive models of bird material
Bird volume is assumed to be constituted of 85% -90% of water and the rest of the air mixture (also known as porosity), which represents the internal cavity of real birds (Niering, 1990). Therefore, in most numerical studies, the density of water ranges from 900 kg/m 3 to 950 kg/m 3 . Besides, in experimental tests, ballistic gelatin is regularly used as Latin American Journal of Solids and Structures, 2020, 17(6), e298 4/15 a replica of real birds which has similar density ranges (McCarthy et al., 2004;Dar et. al, 2013) and provides good agreements with 15% porosity compared with real chickens (Wilbeck and Rand, 1981). To simulate the bird material, for both gelatin and water, three different equations of state (EOS) are successfully adopted by researchers, namely, linear polynomial EOS (Such as: Hanssen et al., 2006;Ivančević & Smojver, 2011;Vijayakumar et al., 2015), Murnaghan EOS (Such as: Liu et al., 2014 and and shock EOS (Such as: Allaeys et al., 2017;Ziaei-Rad, 2011 and. However, for hemispherical bird models, Mie-Grüneisen form of shock EOS predicts better Hugoniot pressure with experiments than other two equations (Arachchige et al., 2020). Moreover, Murnaghan EOS cannot be used directly in some FE codes like Autodyn (Material Models in Autodyn, 2009).
To define the Mie-Grüneisen form of shock EOS, it is considered that particle velocity, p U and shock velocity, s U has an empirical linear relation which can be expressed as follow.
where 0 C is the sound velocity in the projectile and S is the material property. The final expression of pressure as a function of density can be derived as follow (Smojver & Ivančević, 2010).
where η is the nominal compressive volumetric strain which is equal to and m E is the internal energy per unit mass. Finally, a summary can be drawn from the above-mentioned literature that in most bird strike simulation cases, SPH bird models are always preferred along with EOS constitutive laws of water and gelatin. Nonetheless, a recent study suggests that the Mooney-Rivlin hyperelastic material model of ballistic gelatin can be a good alternative of EOS models to obtain realistic bird impact scenarios despite the difficulties of plotting the material constants (Zhou, 2017). Even though a series of successful validations with experiments are achieved, the main limitation of the study is that the Lagrange bird model along with EOS and Mooney-Rivlin constants cannot predict any outcome due to the mesh distortion problem.
To examine the potentiality of the gelatin Mooney-Rivlin material model as bird substitute material, further investigations are deemed necessary. Besides, with the latest development of commercial FE codes, a node erosion algorithm setup of the Lagrange model can be utilized to solve the mesh distortion problem by removing the highly distorted nodes while retaining their inertia during the solution process. Therefore, the main objective of this current research is to adopt the Lagrange bird model with gelatin Mooney-Rivlin hyperelastic constants and node erosion algorithm setup and compare the model for both rigid plate impacts and deformable plate impacts with Lagrange shock EOS and SPH shock EOS of gelatin and SPH shock EOS of water at different impact velocities. A research diagram is illustrated to comprehend the objective of the present investigations, Fig. 2. To carry out the computational research, an explicit solver Ansys Autodyn is employed. Final remarks and conclusions are made based on the numerical results of rigid plate impacts and deformable plate impacts including an experimental comparison with available literature data.

Bird geometry and material
In general, three different geometric bird models, namely straight-ended cylindrical, hemispherical-ended cylindrical, and ellipsoid are well accepted models to simulate the bird strike case events (Meguid et al., 2008). For the present numerical investigations, hemispherical-ended cylindrical bird model with a radius of 57 mm and a total length of 114 mm, Fig. 3 (a) is adopted. Both the mesh counts of the Lagrange model (total elements = 8960), Fig. 3 (b), and SPH model (Total SPH elements = 64297), Fig. 3 (c), are found to be optimal and further increment of elements would only result in more computation time. For experimental comparison, straight-ended cylindrical bird models are designed according to the geometric parameters found in the literature (Zhou et al., 2019a). For bird material models, Mie-Grüneisen EOS parameters of gelatin (Zhou et al., 2019a) and water (Hedayati and Jahanbakhshi, 2015) are considered, Table 1, for comparison with Mooney-Rivlin gelatin constants (Zhou, 2017), Table 2. It is important to note that both the Mooney-Rivlin and EOS gelatin models are characterized from the same type of ballistic gelatin (Zhou, 2017 andZhou et al., 2019a).

Target Plate material
For rigid body impacts, the target wall is modeled with linear elastic properties of structural steel assuming the density, ρ = 7850 kgm -3 , Young's Modulus, E = 2 × 10 5 MPa, and Poisson's ratio, ν = 0.3. For deformable plate impacts, Aluminum alloy 2024-T3/T351 is chosen with the Johnson-Cook material model and damage constants (Buyuk et al., 2009), Table 3, which can be further utilized for experimental comparisons.

Boundary Conditions
For rigid body impacts, a thick wall with a dimension of 500 mm × 92.5 mm × 500 mm is modeled and constrained as rigid. To mesh the wall, cubic solid elements are chosen and four different impact velocities : 95 ms -1 , 117 ms -1 , 145 ms -1 , and 175 ms -1 are given to the bird models to compare their pressure profiles generated at the central node which impacts with the wall at the beginning of the impact, Fig. 4 (a). For Lagrange models, both node erosion option and solution without eroding the nodes are set to run the simulations.
For deformable body impacts, Aluminum Alloy plate with a dimension of 800 mm × 1.4 mm × 800 mm is modeled and placed between two steel frames exposing an area of 600 mm × 600 mm to study the impact of different bird models, Fig. 4 (b). Both the target plate and the steel frames are meshed with cubic solid elements resulting in a total number of elements equal to 23136. The target plate is bonded with the steel frames while the frames are constrained as rigid.
For deformable plate impacts, two different impact velocities are chosen: 117 ms -1 and 175 ms -1 . For the impact velocity of 117 ms -1 , bird models producing deformation at selected nodes, Fig. 4 (c) are compared and generated equivalent von Mises stress and internal energy due to the impact at the target plates are examined. Since the damage predictions by the models are crucial to measure, a higher impact velocity of 175 ms -1 is adopted. Finally, a numerical setup for experimental comparisons is adopted from the literature (Zhou et al., 2019a).

Rigid plate impacts
At first, we will look at the pressure profiles of the node that first impacts with the rigid plate. From the illustration, Fig. 5, it is evident that, apart from the Lagrange Mooney-Rivlin model, all other models predict similar peak pressure (Hugoniot pressure) during the impacts. The predicted peak pressure by the Mooney-Rivlin model is significantly lower than other models, around 115%, which does not improve with increased impact velocity. Secondly, for higher impact velocities (145 ms -1 and 175 ms -1 ) the node is eroded after the impacts and provides a zero pressure line, which implies the elimination of the node from the solution.
To compare the effect of eroded nodes, it is noticeable that both eroded and non-eroded algorithm options provide the same peak pressure values for Lagrange EOS models. In general, the node erosion algorithm setup does not affect the solution since the node has already produced pressure on the object and continues to a steady-state regime before it is completely removed from the solver. Besides, for the impact velocity of 175 ms -1 , the solution of the Lagrange gelatin EOS model without nodal erosion algorithm cannot be converged fully due to the energy error problem, Fig. 5 (iv). However, the generated convergence problem can be tackled by introducing the node erosion technique.
For all impact cases, after reaching the peak pressure, the steady-state regime is found to be quite discontinuous for both the SPH EOS models and the Lagrange Mooney-Rivlin model; however, the steady-state pressure is much more disciplined for Lagrange EOS model. Overall, the pressure profiles of the SPH models are quite similar, which is something completely unlike between the Lagrange EOS and Mooney-Rivlin model.

Deformable plate impacts -I (Impact velocity = 117 m/s)
In comparison with the rigid plate impacts, deformable plate impacts are more practicable since the aerostructures are mostly made of deformable shells and plates. For the impact velocity of 117 ms -1 , nodal deformation results are plotted for the nodes placed at eight different locations on the target plate, Fig. 6. From the illustrations, it is observed that, for all bird models, the deformation predictions are very close to each other. Despite the similarities, SPH gelatin EOS predicts slightly higher values, which are only 2.5% to 4.5% higher than the Mooney-Rivlin model at peak stages of the deformation, occurs mainly at around 2 milliseconds by the models. Considering the area of the plate, these differences are very minor and as the time progresses, they are even minimized to a much lower percentage, around 1% to 2%. In most cases, the SPH water EOS and Lagrange Mooney-Rivlin model produce the same deformation outcomes.
Among the considered nodes, node I, which is located at the center of the impact, experiences the highest deformation during the impact. The predicted maximum deformation peaks are 82.2 mm, 81.8 mm, 81 mm, and 80.4 mm for SPH gelatin EOS, Lagrange gelatin EOS, SPH water EOS, and Lagrange Mooney-Rivlin model respectively. As the nodal distance from the impact center increases, the deformation of the nodes decreases consequently. For instance, node IV, which is located 84 mm away from the center horizontally and the farthest node examined for this study, exhibits the lowest deformation results, precisely, 56 mm, 55 mm, 54 mm and 53.2 mm for SPH gelatin EOS, Lagrange gelatin EOS, SPH water EOS, and Lagrange Mooney-Rivlin model respectively. It is also important to note that, in comparison with the SPH models, both the Lagrange models with distorted node erosion algorithm setup can successfully predict the deformation results without losing any accuracy. Besides, in case of retaining the distorted nodes in the solution, both the models fail to converge due to energy error as specified previously during the rigid plate's impact situation.
To compare the stress results generated in the plate due to different bird models, equivalent von Mises stress contours are plotted with an interval of 1 ms during the whole impact phenomenon, Fig. 7. At 1 ms, when almost half of the bird geometric model is dissipated, the predicted maximum equivalent stress results are found to be quite similar for all the models where the Lagrange Mooney-Rivlin model predicts slightly higher stress, around 5% more than others. As time progresses, at 2 ms, all the target plates experience the highest equivalent stress within the entire duration of the impact. The estimated maximum equivalent stress results are identical for both Latin American Journal of Solids and Structures, 2020, 17(6), e298 9/15 Lagrange gelatin EOS and Lagrange Mooney-Rivlin model, to be specific, 624.7 MPa and 627.4 MPa respectively. These predicted outcomes do not differ significantly with SPH models, only around 4.5%. At a later stage, for 3 ms, equivalent stress results are found to be quite different from each other. Since the dissipation of the bird geometric models does not occur with the same rate, therefore, the generated stress waves travel with different timing in the plate and develop different stress contour bands at this stage. Finally, at 4 ms, when the impact phenomenon is completely terminated, stress waves become weaker in the plate and the predicted values are appeared to be close with each other.
Next, we compared the absorbed energy plot of aluminum target plates due to the impact of different bird models, Fig. 8. From the illustration, it is evident that, despite the differences of bird's finite element models and material models, the absorbed energy is quite similar and the maximum internal energy generated by the plate is 5.56 kJ, 5.57 KJ, 5.31 KJ and 5.29 KJ for SPH gelatin EOS, Lagrange gelatin EOS, SPH water EOS, and Lagrange Mooney-Rivlin model respectively.

Deformable plate impacts -II ( Impact velocity = 175 ms -1 )
For the impact velocity of 175 ms -1 , some interesting observations can be made, Fig 9. Due to the high impact velocity of the bird models, the target aluminum plate is found to be completely damaged at the center and a hole is created with a round shape. However, the length of the perforation varies with the bird models. Both the Lagrange gelatin EOS and SPH gelatin EOS initiates similar perforation lengths of 243.2 mm and 245 mm respectively, which are almost 13% higher than the predicted perforated length of Mooney-Rivlin model. Nonetheless, SPH water EOS and Lagrange Mooney-Rivlin model generate almost similar perforation, to be exact, 219 mm and 215.5 mm respectively. Since complete damage is experienced on the plates, no further result plots are discussed for this section.

Computational time
To solve the bird impact problems, we have employed a powerful computer with AMD ryzen R7-3700x octa-core 5 GHz processor and 16 GB DDR4 RAM. However, to use parallel computing, Ansys Autodyn requires a high-performance computing (HPC) license with additional cost. Therefore, serial computing is used in this study. From the illustration, Fig. 10, it is evident that for both the impact velocity problems, adopting the Lagrange Mooney-Rivlin model provides a significantly faster solution while comparing with other bird models.
To investigate the impact results at 117 ms -1 , the Lagrange Mooney-Rivlin model requires only 43 minutes, which is around 125% and 159% faster than Lagrange gelatin EOS and SPH models respectively. A similar conclusion can be drawn for the impact velocity of 175 ms -1 , when Lagrange Mooney-Rivlin model takes about 27 minutes to compute the damage, and SPH models and Lagrange gelatin EOS model complete the computation taking 198 minutes (3 hours 18 minutes) and 107 minutes (1 hour 47 minutes) respectively. Nonetheless, we did not consider the comparison of computation time for rigid plate impacts since all the models provide almost similar solution time, varying maximum 8 to 10 minutes, despite the variation in modeling and generally, the computation time is much lower than deformable plate impacts.

Deformable plate impacts (Experimental Comparison)
Finally, we have compared the Lagrange Mooney-Rivlin bird model with the available experimental data found in the literature (Zhou et al., 2019a). To study the out of plane displacement results at different ballistic gelatin impact velocities, the authors have used 3D digital image correlation method (DIC) which is a non-contact optical method to capture the deformation of the target place accurately without the necessity of using strain gauges and accelerometers. From the illustration, Fig. 11 (i), a linear function of displacement and impact velocity can be observed for both experimental investigations and numerical studies using Lagrange Mooney-Rivlin model. Overall, an excellent correlation between the experimental data and numerical model is achieved. Furthermore, at an impact velocity of 122 ms -1 , Fig. 11 (ii), a comparison of central displacement profiles at the back of the aluminum target is illustrated. The measured discrepancy of peak displacement between the experimental data and Lagrange Mooney-Rivlin model is only 3%, which establishes the Lagrange Mooney-Rivlin model as an outstanding substitute of bird model and material. Besides, the model predicts better deformation profile than SPH EOS (Abaqus) outcome. Along with the deformation results, images of the gelatin deformation captured by the high-speed camera are also compared with the images of the Lagrange Mooney-Rivlin model, Fig. 12. From the illustration, it is evident that, the radial spread of the projectile is accurately generated by the model. Besides, the correlation between the numerical models are also found to be excellent. It is important to note that, as previous numerical case studies, node erosion algorithm setup is kept to predict the outcomes of the impact cases and found that without the erosion of the distorted nodes from the solutions, the solver provides an energy error with a premature state.

CONCLUSION
Numerical analysis is an essential tool to evaluate the bird strike case studies of aerostructures in the aviation industry. Therefore, it is crucial to model the bird with appropriate material and computational code. A large number of simulations are conducted in this study to establish the computational bird model with ballistic gelatin Mooney-Rivlin material properties along with Lagrange coding. To validate the efficiency of the model, other finite element methods and material modeling are utilized, namely, SPH gelatin EOS, Lagrange gelatin EOS, and SPH water EOS. Both rigid body impacts and deformable plate impacts are considered under different impact velocities. Besides, a node erosion algorithm setup is applied to remove the highly distorted node from the simulation during Lagrange modeling. Finally, numerical results are compared with experimental data from the literature. After the successful investigations, the following remarks can be drawn.
1. Despite the discrepancy found during rigid body impacts, the Lagrange Mooney-Rivlin model can accurately predict the outcomes of deformable plate impacts. The generated nodal deformation, von Mises stress contours and the absorbed energy of the target plate are almost identical to other models.
2. The Lagrange Mooney-Rivlin model can be precisely utilized for damage evaluation of aerostructures. Despite some differences observed during the damage predictions, it has generated a similar impact with SPH water EOS model.
3. By adopting the Lagrange Mooney-Rivlin model, bird strike solutions can be accomplished much faster compared with other models with similar accuracy. In particular cases, around 125% faster than Lagrange EOS and 159% faster than SPH EOS makes the model a computationally efficient code for further case studies.

4.
A close agreement between the experimental data and numerical predictions of the Lagrange Mooney-Rivlin model confirms the convenience of the computationl model as a substitute of bird in future finite element investigations.
5. Node erosion algorithm makes the Lagrange solver distortion-free, by removing the highly distorted nodes from the solution. After comparing all the numerical predictions and experimental results, it is evident that node erosion algorithm has no adverse effect during the solution. Rather, it helps the solver to predict accurate outcomes by providing an energy error-free solution while retaining the inertia of the nodes.