Finite element bird-strike modeling

Finite element bird-strike modeling 6.1

6

Introduction

Traditionally, experimental techniques were the only tool for the certification of aircraft components (Hedayati & Ziaei-Rad, 2013b; Hedayati, Ziaei-Rad, Eyvazian, & Hamouda, 2014). As explained previously in Chapters 4 and 5, experimental techniques are very costly (a single bird shooting test can cost several thousands of dollars) and also do not provide enough information about the structural response (Guida, Marulo, Meo, & Russo, 2013). The numerical techniques provide the designer with a wide range of useful data (e.g. stress distribution, displacements, 3D visual observations of structure deformation). Moreover, in numerical simulations, changing the material or geometry is possible without any additional costs. With the use of numerical simulations, the number of required tests is decreased significantly (Hedayati & Ziaei-Rad, 2013b; Hedayati & Ziaei-Rad, 2012). As specified in FAR 25.631, the experimental tests can be replaced by numerical analyses provided that the numerical models are validated with experimental tests on similar structures. For instance, this validation procedure was done for the upper cap zone of Airbus A380 (Faure, 2011; Heimbs, 2011). Due to the great advances in the efficiency of numerical methods, which is the result of improvement in computational hardware as well as solution algorithms, more and more aircraft manufacturers are using numerical methods to design novel bird-proof structures (Guida, 2013). Among the several available numerical methods, the finite element method (FEM) has obtained almost all the attention in bird-strike studies because of its abilities in analyzing structures with very complex geometries, material behaviors, and loading conditions (Hedayati & Ziaei-Rad, 2013b). The ongoing improvements in computer technology make it possible to carry out complex analyses in higher detail and within shorter timescales. Consequently, using numerical methods, rather than experimental tests, has become more and more justified every forthcoming year. The impact of a bird into an aircraft is an impact dynamics problem with a very short duration. Compared to quasi-static simulations, bird-strike modeling includes several complexities such as transient intense loads, fluid–solid interactions, strain rate effects, large deformations, severe element distortion, etc. (Wang & Yue, 2010). Since the duration of the event is very short (and, therefore, the materials are highly accelerated), the effect of material inertia cannot be ignored. Stress waves are so strong that their transmission and mutual interactions have to be taken into account. Due to the different natures of solids and fluids, proper modeling of their force interaction and impulse transfer poses a considerable problem. When the deformations are large, the impacted structure goes under plasticity (material behavior becomes nonlinear), and the effects of shear deformation and rotational inertia have to be considered in the plate or beam theories. Another complexity associated with bird-strike modeling is that high velocity

Bird Strike. http://dx.doi.org/10.1016/B978-0-08-100093-9.00006-6 Copyright © 2016 Elsevier Ltd. All rights reserved.

114

Bird Strike

impacts increase the strain rate which in turn increases the yield stress and the ultimate strength of the material (known as material nonlinearity) (Wang & Yue, 2010). Finite element modeling of bird strikes dates back to the late 1970s, when a few research studies were carried out to examine the response of the canopies and windshields of airplanes impacted by birds. Linear FE code IMPACT that was developed in 1977 was not able to yield acceptable results for bird strike (Heimbs, 2011). In 1978, the computer program MAGNA (Materially And Geometrically Nonlinear Analysis) was developed and 2 years later was applied to the impact of birds on the transparent components of aircraft (McCarty, 1980a,b). In the mid-1980s, in addition to the transparent components of aircraft, the bird-proof capability of other aircraft parts including laminated composites (Sun, 1972), horizontal stabilizers (McCarty, 1982,1983), and turbine engine fan blades (Hirschbein, 1982; Storace, Nimmer, & Ravenhall, 1984) were also investigated by means of finite element (FE). The computational power of very expensive computers at that time was very low when compared to a cheap personal computer today (Blair, 2008). Regarding this issue, in one of the first works undertaken on bird strike (Belytschko, Privitzer, Mindle, & Wicks, 1979), the authors stated that: “Simulations were first performed of bird-canopy impacts to determine the mesh refinement necessary to reproduce the magnitudes of experimentally observed deformations and the wave pattern of the observed displacement. It was found that only a very fine mesh proved satisfactory. This mesh could not be used in the computer simulations because of computer core-storage and cost limitations at our computer facility.” Therefore, the amount of work undertaken on bird strike in the 1980s was relatively low and remained so until the turn of millennium when the computing power of computers became sufficient for the accurate simulation of complex problems with relative ease (Blair, 2008). In order to properly model the structural deformation and material failure of both the bird projectile and the impacted structure, the implemented FE package must be able to handle the following areas (Guida et al., 2013): l

l

l

l

l

l

l

l

l

l

strain rate effects; excessive element distortion; hourglass elements; accurate prediction of load distribution in the interface between the bird and the target; modeling the fasteners (usually accomplished by defining the single node-to-node tie); reliable material behavior at very severe conditions (which usually happens at least at a few points of the FE model); geometrical and material nonlinearity; fast and accurate solutions; energy dissipation; and sufficient outputs.

The FEM is based on the idea of discretizing a complex geometry into a large number of small pieces with simple basic geometries (usually triangle and quadrilateral for 2D analyses, and hexahedral and tetrahedral for 3D analyses). The sufficiently small size for the elements of a FE model is the size that the normal or shear stress distribution on the surfaces of each element can be considered constant with acceptable accuracy. The responses of simple structural geometries under the action of the constant values of

Finite element bird-strike modeling

115

normal or shear stresses on their sides (in 2D analyses) or faces (in 3D analyses) is known. Therefore, considering the responses of each of the small elements to a given stress field, the response of the entire system can be obtained. The vertices of the constructing pieces (or elements) are known as the nodes of the FE model. In Fig. 6.1, a plate with two holes is shown pre- and post-discretization. Analyzing a physical problem in the FE packages is done in three main steps: l

l

l

Pre-processing: the FE model which is going to be solved is prepared. Solution: the prepared FE model is solved. Additional controlling parameters are usually used in this step to avoid solution instabilities. Post-processing: the quantities of interest are obtained from the solved problem, such as stress, strain, displacements, velocities, accelerations, resultant forces, visualization, etc.

The solution step is usually completely handled by the FE package, and therefore, the user does not have to deal with it. In some FE packages (e.g. ANSYS and ABAQUS), all the three steps are contained in a single graphical user interface (GUI). Whilst in others, the user has to use different programs for each step (e.g. in LS-DYNA which is the most prevalent program in bird-strike analysis, the solution program is separate from the pre- and the post-processing consoles). In some cases, the user may prefer using other more professional programs for the better handling of a specific step. For instance, some users implement the program HyperMesh for the more efficient discretization of very complex geometries. Pre-processing is the process in which all the required data by the solver program is prepared. In a pre-processing procedure, the following matters must be addressed: l

Creating the 2D or 3D CAD model of the physical object. In very complex geometries, the geometrical model can be prepared in a professional CAD program, and then exported into the FE pre-processing program. Figure 6.1 FEM discretization of a simple plate with holes.

(a)

A plate with two holes

(b)

The plate after discretization

116 l

l

l

l

l

Bird Strike

Specifying the material models that determine the behavior of each material type under different loading conditions. Specifying the material constants of each material. Selection of the suitable element formulation for different parts of the problem. Discretization of the model with well-shaped elements. Applying the loads and boundary conditions to the corresponding nodes.

In the solution step, first, the local stiffness matrix and the force vector of each element is constructed. By combining all the local stiffness and force matrices, a global stiffness matrix and a global force vector are constructed. The program then finds the displacements at the nodes by solving the set of equations: fF g ¼ ½ K  fd g

(6.1)

The stiffness matrix and force vector are usually very large. In fact, their size is equal to the number of nodes of the FE model. The solver programs usually first attempt to shrink the stiffness matrix and the force vector using matrix algebra procedures. After that, the FE solvers usually use advanced numerical methods to solve the huge set of equations simultaneously. The accuracy and speed of a particular FE solver depends on the efficiency of the shrinking and solving algorithms. Generally, the FE solvers solve a field problem by using either explicit or implicit techniques. If the load is high enough that the material under investigation is yielded or the displacements are larger than static displacements, the FE problem cannot be solved using Eq. 6.1 in a single step. In these conditions, the total duration of time through which the material is going to be deformed must be divided into very small time increments (or time steps). In an explicit analysis, depending on the amount of load or strain rates, the response of the material in each time increment can be different. At the end of each time increment, the geometry and material behavior of the structure is updated and the stiffness matrix elements are modified for the next increment. In the next time increment, the updated stiffness matrix is used for obtaining the new displacements. Similar to explicit analysis, the implicit FE solver uses the updated stiffness matrix in each time increment. The difference between the two analysis methods is that the implicit solver does Newton–Raphson iterations to ensure that the internal forces of the structure are in equilibrium with the externally applied loads. In contrast to the explicit method (in which the state of the system at a time step is calculated based on the state of the system in the previous time step), the implicit method finds the current state of the system by considering the state of the system at both the previous and the next time steps. This is why the implicit solutions usually give better accuracy. However, they are usually very time-consuming and therefore not suitable for problems with very large deformations. The implicit analyses are mainly useful for the transient problems which take a long time to occur. On the other hand, the explicit analyses are useful for high-velocity problems with very short durations in which the effect of shock waves is important (Guida, 2008). In some mechanical applications such as metal forming – when the process duration is neither short nor long – the

Finite element bird-strike modeling

117

advantages of using either the explicit or the implicit analyses are similar and consequently both of the analyses methods are in common use. Smaller time steps lead to greater accuracy, whilst larger time steps yield faster solutions. To obtain sufficiently accurate results, the time step must be as small as the time it takes for a shock wave to travel through the smallest element in the model. Therefore, the higher the velocities of the objects in a simulation are, the smaller the required time step will be. Time steps which are small enough ensure that all the energies related to the propagating shock waves are considered in the analysis (Guida, 2008). In a bird-strike problem with bird initial velocity of about 100 ms and with ordinary element sizes (e.g. the bird is constructed using about 5000–15000 elements), the automatically calculated time step is of the order 0.1–1 μs. For the problems in which the impacted structures are modeled with all the micro-geometrical details, the time step can be decreased to less than 1 ns. For a bird-strike problem which usually lasts for 1–10 ms, the small size of the resulting time steps means that the solver needs to solve the problem incrementally for more than thousands of times. Moreover, since the implicit analyses try to enforce the equilibrium of the system by Newton– Raphson iterations, the presence of a high amount of non-linearity (caused by material and/or geometrical nonlinearities, as well as large displacements) in the system makes the convergence in the implicit analysis more difficult. This is why using explicit analysis in a bird-strike problems is well-established. In the post-processing step, the user is able to obtain the required data from the results of the solver program after the completion of the solution. The most common results that are important for a bird-strike problem are pressure distribution in the interface surface, the permanent plastic deformations in the impacted component, the eroded elements of the structure, and the time-history visualizations. In the following paragraphs, the main sources of nonlinearity in a bird-strike simulation are described.

6.2

Structural nonlinearity

6.2.1

Geometric nonlinearity

Many engineering problems are solved by assuming that the strain and displacement i are related to each other through the equation εi ¼ @u @xi , for i ¼ 1, 2, 3. However, when the deformations of the problem becomes large, the relationship between the strain and displacement becomes non-linear and takes the form:
 1 @ui @uTi @ui @uTi + + εi ¼ 2 @xi @xi @xi @xi

(6.2)

Geometric nonlinearity usually occurs in very slender structures in flexure, e.g. shells and very slender beams. Tensile forces can also lead to geometrical nonlinearity in cables and inflatable membranes (Grimaldi, 2011).

118

6.2.2

Bird Strike

Material nonlinearity

In low stresses, the materials usually show a linear relationship between the applied load and the resulting displacement. In other words, the strain is related to the stress value proportionally. However, almost all the engineering materials show a nonlinear behavior if the applied stress exceeds their yield stress. In high-velocity events, the strain rate can also affect the stress-strain curve. Behavior of the materials at different temperatures is also very variable. There are other sources of nonlinearity in the behavior of the material, e.g. creep, viscoelasticity, and pre-stress. In all the aforementioned conditions, the linear material models cannot be used for FE modeling.

6.2.3

Force and displacement nonlinearity

In problems with small displacements, the applied loads or boundary conditions usually keep their orientations with respect to the edges of the studied structure. However, in some problems, such as those which include fluid–solid interactions or the problems in which the external load is applied through a large contact surface, the orientation and amount of the applied load is also dependent on the displacement of the interacting structure. In some problems, the boundary condition state can also be variable and dependent on the response of the system itself. Different examples of problems with nonlinear loading conditions are aerodynamic forces applied to an aircraft, a boat moving in a canal, the impact of two cars, and a bird-strike event. In the 1980s, the loading condition between the bird and the target was considered as a uniform normal pressure over a “bird impact footprint” (McCarty, 1980b; Heimbs, 2011). The pressure amount was set to the pressure the bird imposed on a rigid plate in its steady regime, and the Hugoniot phase was neglected. Later, Baughn and Graham (1988) considered the gradual spread of the bird over the target by increasing the “footprint” area with time (Heimbs, 2011). Although it was known for decades prior that the amount and the distribution of the impact loading heavily depends on the response of the impacted structure (McCarty, 1983; Engblom, 1980), it was only after the development of modern high-tech computers that modeling the actual fluid–solid interactions became possible.

6.3

Numerical approaches for bird strike

As stated above, in the early bird-strike simulations, the imposed load applied by the bird to the target was decoupled from the response of the target. In those simulations, the bird FE model was not created, and the effect of the birds was applied to the structure through temporal and spatial variation of load with several simplifying hypotheses (Wang & Yue, 2010). In later bird-strike studies, the bird model was also created, and its response was coupled to the response of the impacted structure by keeping the compatibility conditions in the contact area. The compatibility is maintained in the contact surface, if the applied load and displacement of the bird and the target structure are identical at each contact point. In the majority of the works published after 2000, the bird and target responses are coupled.

Finite element bird-strike modeling

119

In a bird strike, the bird behaves like a fluid, and in very severe situations is broken up into very small debris particles. Capturing this behavior of the bird with acceptable accuracy and simultaneously avoiding numerical instabilities has proven to be a big challenge (Hedayati & Ziaei-Rad, 2013b). Different numerical approaches have been presented for modeling a bird in a bird-strike analysis, which will be explained below. The first, and yet the most common, approach to model the bird is the Lagrangian method. The Lagrangian method is the default approach for discretizing solid parts in FE packages. This method uses the material coordinates as the reference, i.e. the nodes of the mesh are attached to the particles of the material. This model is able to follow the distortions of the bird material, and to some extent, the break-up of bird material into debris (Stoll & Brockman, 1997). However, there are several drawbacks in Lagrangian bird modeling. Since the nodes of the FE model are attached to the material, if the material is deformed greatly, some of the elements become highly distorted. As a result, the smallest dimension of the elements becomes very small, which in turn decreases the solution time step as well. The huge drop in the time step duration significantly increases the number of required time steps, and as a result the solution time increases (Hedayati & Ziaei-Rad, 2013b). A strategy to overcome the solution instabilities caused by small time steps is simply eliminating the highly distorted elements. However, if the number of distorted elements are relatively large compared to the total number of bird elements, the mass of the bird model drops significantly, which can lead to unrealistic results. Other techniques have also been proposed for dealing with the highly distorted element in the Lagrangian method, but have not been very efficient (Hedayati & Ziaei-Rad, 2013b). Alternative numerical modeling approaches have been proposed for overcoming the difficulties related to Lagrangian bird modeling (Jenq, Hsiao, Lin, Zimcik, & Ensan, 2007). Since at high pressures, the bird body behaves like a fluid, the Eulerian approach which is prevalent in modeling fluid dynamic problems would be helpful. In this method, a fixed void mesh is created in the space, and some of the cells are filled by bird material at the points where the bird must be present. As the bird material travels into the space, some cells become hollow and some others become filled with the bird material. The Eulerian approach requires a relatively fine mesh to yield sufficiently accurate results. This significantly increases the computational time required. A similar and more practical approach to modeling a bird strike is the Arbitrary Lagrangian Eulerian (ALE) method (which is basically similar to the Eulerian method). In the ALE method, the material is free to move with respect to the mesh. The only difference is that in the ALE method, the mesh is not necessarily fixed in the space and can move in the direction the projectile’s center of gravity moves. Therefore, the large number of void space elements in the Eulerian method can be reduced in the ALE method, and the void elements can be limited to only the regions around the bird model (Jenq et al., 2007). This method benefits from the stability of the Eulerian method and the low computational cost of the Lagrangian method. The interaction of the bird with the target surface is modeled by means of an ALE coupling algorithm. However, this method is relatively difficult to utilize, and the user must be familiar with a large set of controlling parameters to achieve acceptable results.

120

Bird Strike

The other proposed approach for bird modeling is smoothed particles hydrodynamics (SPH) which is a grid-less Lagrangian technique for modeling transient fluid motion using a pseudo-particle interpolation method (Hedayati & Ziaei-Rad, 2013b). This method was first developed by Monaghan in the 1970s for solving the movement of particles in the space under the interaction of each other, and with applications in hypervelocity impacts (Monaghan, 1992). Several years later, the SPH method was found to be also effective and accurate in predicting the deformation of fluids. In addition to the fluid computational problem, the SPH method can also be applied to solid mechanics problems in which deformations are very large, e.g. crash simulations (Grimaldi, Sollo, Guida, & Marulo, 2013). In the SPH method, the fluid is represented by a cloud of moving small particles, each one being an interpolation point, where all the fluid characteristics are known. An interpolation function (known as kernel function) is used to find the desired quantities for all the particles. The kernel function is active only over a limited zone for each interpolation point (Hedayati & Ziaei-Rad, 2013b). Each interpolation node is given a mass, and the values of state variables for the node are determined based on the mass of the node itself and the masses and distances of the adjacent nodes in its zone. Since the SPH method is mesh-less, the difficulties associated with severe element distortions do not exist. Moreover, no additional elements are required to be defined to represent the void space (whereas the Eulerian or ALE methods do) providing a decrease in computational time.

6.3.1

Lagrange

In the Lagrangian formulation, the mesh nodes are attached to the particles of the material. Therefore, each node represents one particle of the material under examination. The reference coordinate system in Lagrangian formulation is the material coordinate X. The motion of each particle in this formulation is described by (Alberto, 2006): x ¼ φðX, tÞ

(6.3)

where φ(X, t) is the function mapping the initial position of the material to its current position. The displacement of a point is defined as the distance between the initial and current positions of the point (Alberto, 2006): uðX, tÞ ¼ φðX, tÞ  X ¼ x  X

(6.4)

In this method, since the mesh nodes follow the material deformation, and therefore the boundary nodes remain on the external surfaces of the material, imposing the boundary conditions is simple. The other advantage associated with the Lagrangian method is the simplistic traceability of each material point. However, as described previously, using Lagrangian formulation for materials under very excessive deformations can lead to large distortions of the elements, which in turn increases the computational time (Guida, Marulo, Meo, Grimaldi, & Olivares, 2011). For having accurate results, the time

Finite element bird-strike modeling

121

step Δt must be at least smaller than the time required for a shock wave to travel through the smallest dimension of the element, lmin: Δt ¼

lmin c

(6.5)

where c is the speed of the sonic waves in the material. As a result of the excessive distortion in an element, while one of the dimensions of the element increases greatly, its other dimension (lmin) decreases to an unacceptable low value leading to very small time steps (Δt). The other problem that can arise from excessive element distortion is that the volume of some elements becomes non-physically negative, because they fold in on themselves, see Fig. 6.2. Several techniques have been presented to deal with the excessive mesh distortions in the Lagrangian method which include element removing, local mass-scaling, smallstrain simplifications, nodal masses (NM), and adaptive remeshing/rezoning. In the element removing technique, a critical failure strain is defined for the elements. The elements in which the effective strain exceeds the critical strain are removed. This technique has been shown to be very successful in avoiding numerical instabilities. However, the removal of a large number of elements also removes their mass and strain energy which causes violation in the principles of conservation of mass and energy. Moreover, since in a bird strike the majority of distorted elements are in contact with the target surface, removing them creates unrealistic severe fluctuations in the pressure profile read from the target surface, especially for coarse meshes (Nizampatnam, 2007). Mao, Meguid, & Ng (2008) suggested using “highly refined meshes” to remove the artificial noise in the contact force diagram in the Lagrangian simulations that include element deletion criteria. In addition to artificial fluctuations in the contact force, removing the distorted elements also decreases the average force applied by the bird to the target. In this regard, Castelletti and Anghileri (2004) stated that introducing a failure criterion for the bird model deteriorates the experimental–numerical correlation. Figure 6.2 Non-physical negative volume for an element.

(a)

Undeformed element Negative volume

(b)

Distorted element

122

Bird Strike

In the adaptive remeshing/rezoning technique, the regions that include several distorted elements are remeshed. This procedure increases the solution time and is a very complex task. In fact, if it is desired to do the process automatically, a very advanced remapping algorithm is required. In practice, the available remapping algorithms are not very successful in accurately mapping the original distorted mesh to the new wellshaped mesh and cause numerical errors, especially for complex geometries (Nizampatnam, 2007). Consecutive rezoning can lead to uneven external surfaces of a 3D object which make them unsuitable for applying boundary condition or contact. Some authors (e.g. Langrand, Bayart, Chauveau, & Deletombe, 2002) suggested using the “small-strain” simplification to overcome the severe element distortion or tangling problem. In this technique, the Jacobian matrix is not updated at each time step, and therefore the large structural distortions are neglected. Another method for avoiding the problems caused by distorted elements is using the mass of elements instead of their length to calculate the time step. In this case, the time step is called nodal and not elementary. To do this, the mass of distorted elements are increased artificially to keep the time step constant (local mass-scaling). However, in large strain problems such as bird strike, this technique may not be suitable because the final mass of the bird may become significantly large (Langrand et al., 2002). Advantages, disadvantages, and the relevant enhancements associated with the Lagrange formulation are listed below in Table 6.1.

6.3.2

Eulerian

For a better understanding of the Arbitrary Lagrange Eulerian (ALE) method, it is necessary to describe the Eulerian method. This method is mostly useful to fluid dynamics problems. In the Eulerian method, instead of material, the space is discretized. A mesh Table 6.1 Advantages and disadvantages of the Lagrangian formulation in bird-strike application Advantages

Disadvantages

Solution

Easy tracking of the time-history properties of each particle of material

Severe element distortion can increase the number of required time steps

Element erosion, remeshing, massscaling, small-strain

Simpler imposing of boundary conditions on the material interfaces

Element distortion can cause element tangling

Element erosion, remeshing, massscaling, small-strain

Low computational cost

Bird material cannot be split-up into debris

Simpler modeling and low number of analysis parameters that have to be dealt with by user

Finite element bird-strike modeling

123

consists of several stacked cells fixed in the space. Some of the cells of the mesh are initially filled by fluid and some are left hollow. As the fluid flows in the space, it leaves some cells and fills some initially hollow cells. Since the mesh does not follow fluid deformation, severe mesh distortions do not exist. Unlike the Lagrangian method, tracking the time history of material points as well as material interfaces is very difficult to be performed in the Eulerian method (Salehi, Ziaei-Rad, & Vaziri-Zanjani, 2010). The description of the Eulerian motion is the opposite of that in the Lagrangian formulation. This method uses the spatial coordinate system to express the material coordinates. Eulerian mapping can be defined as the inverse of Lagrangian motion mapping (Alberto, 2006): X ¼ φ1 ðx, tÞ

(6.6)

Since the Eulerian mesh is fixed in the space, it must include not only the locations where the fluid initially exists, but also all the locations where the fluid might be present at a later time. This fact does not pose a big difficulty in prevalent fluid dynamics problems in which the fluid usually flows around fixed objects such as pipe walls, skyscrapers, or airfoils. However, if the Eulerian mesh is going to be used for following the path of high velocity objects, such as a bird, the required meshed volume is much larger than the volume of the bird itself. This is why using the Eulerian mesh for a bird-strike problem is computationally expensive and therefore very rare (Nizampatnam, 2007). The other problem which accompanies the Eulerian approach for bird strike is that tracking the history of material particles in the domain is not easy at all. However, in order to do this, the stress and strain tensors can be moved from one cell to another. This measure does not usually lead to accurate results (Nizampatnam, 2007). The decoupling of the material and the mesh grid also makes it hard to follow the interface of the fluid, which is necessary if it is intended to engage solid–fluid interfaces. Whilst in reality, the bird external surface is usually curved, the Eulerian grids can only provide simple sinuate cubic boundaries which in fact are the diffused boundaries of a real bird. The Eulerian formulation also has its own sort of numerical problems. As explained by Nizampatnam (2007): “most Eulerian solvers are based on low-order difference schemes which have built-in sources of error, such as numerical dispersion and numerical dissipation, due to truncation of the Taylor series approximations. These numerical errors can be minimized by resorting to higher-order differencing schemes, but this comes at the expense of increased computational effort”. For overcoming all the inaccuracies and difficulties explained above, the mesh grid must be refined greatly making the problem solution very time-consuming (Mao et al., 2008; Birnbaum, Francis, & Gerber, 1999). Advantages, disadvantages, and relevant enhancements related to the Eulerian formulation are listed in Table 6.2.

124

Bird Strike

Advantages and disadvantages of the Eulerian formulation in bird-strike application

Table 6.2

Advantages

Disadvantages

Enhancements

No element distortion

Diffusion in the material boundaries

Using sophisticated interface tracking implementations

Large time steps

Tracking of the material history is difficult

Transferring the stress and strain tensors from one cell to another

Numerical dispersion and numerical dissipation

Using higher-order differencing schemes

A wide-meshed domain is necessary which increases computational time

6.3.3

Arbitrary Lagrangian Method (ALE)

The ALE formulation is a combination of the Lagrange and Eulerian formulations in which the reference is set arbitrarily by the user, in order to capture the advantages of the two methods whilst minimizing the disadvantages. In the ALE method, the mesh is arranged to be independent of the fluid motion; therefore, problems such as element distortion and tangling would not easily happen. However, the grid can move or contract/ expand to follow boundary motions. Since the mesh can move with a speed usually equal to the average velocity of all the particles of the fluid, the size of the meshed domain can be greatly decreased compared to the Eulerian method (Heimbs, 2011). The LS-DYNA software is capable of dealing with ALE domains filled by more than one type of fluid (e.g. bird and air in a bird problem) (Lavoie et al., 2007). In the ALE method, the location of the bird material is evaluated at each time step by comparing its position with respect to the nodes of the Eulerian mesh. The ALE simulation consists of three main stages at each time step: Lagrangian phase, Eulerian phase, and a smoothing phase in between (Guida, 2008). Mesh distortion can also become a problem in the ALE method, and elements may get negative volumes. For example, as reported by Anghileri, Castelletti, and Mazza (2005), when modeling the bird impact against engine fan blades using the ALE method, the background mesh became severely distorted and negative volumes appeared. However, in general, the ALE method tolerates greater material distortion than the Lagrangian formulation does because of the smoothing phase. Possible interaction of the bird material with a Lagrangian structure is evaluated by tracking the relative distance between the bird and the target (Lavoie et al., 2009). The bird material existent in the Eulerian mesh applies load to a Lagrangian target through the ALE coupling algorithm. Langrand et al. (2002) reported that the pressure profile obtained from the ALE method is smoother than that obtained from the Lagrangian bird model. They believed that the smoothing stage between the Lagrangian and Eulerian phases in the ALE formulation filters the dynamic load applied to the impacted structure, and therefore, makes the bird behave in a less aggressive manner.

Finite element bird-strike modeling

125

In the ALE method, the reference coordinates are denoted χ . If the particle position
tÞ and the mesh motion is defined as
x ¼ φðχ , tÞ, the mesh is defined as χ ¼ φðX, displacement can be given by (Alberto, 2006):

ðχ , tÞ ¼
x X u

(6.7)

The relationship between the ALE coordinates and the material coordinates is given by:
tÞ, tÞ ¼ Ψ ðX,
tÞ χ ¼ φ1 ðφðX,

(6.8)


tÞ ¼ φ1 oφ. where Ψ ðX, In Fig. 6.3, the different ways the aforementioned numerical approaches discretize a fluid (and its background void mesh, in the cases of the Eulerian and ALE methods) and the way the numerical models of the fluid deform in the space is depicted. Obviously, compared to the Eulerian method, the ALE method can better capture the material boundaries due to the deformations of the background mesh. Since in the ALE method, the background mesh is allowed to deform in accordance to the deformation of the fluid, the transport of fluid particles between the different cells is not as necessary as in the Eulerian method. However, transport of the fluid particles

Before deformation

After deformation

SPH

Lagrangian

ALE

Eulerian

Figure 6.3 Undeformed and deformed elements (and the background void mesh, if applicable) in the following methods: (a) SPH, (b) Lagrangian, (c) ALE, and (d) Eulerian.

126

Bird Strike

between different cells of the ALE does occur. (In Chapter 8, the flow of bird particles in a deforming ALE background mesh is effectively depicted in Fig. 8.15.) It must be kept in mind that since the number of ALE background elements remains constant throughout the bird deformation, when the bird spreads on the target, the total ALE structure has to be extended laterally to still include all the bird particles. Similar to Eulerian formulation, the accuracy of the ALE method significantly depends on mesh size. Therefore, after spreading of the bird model on the target, the accuracy of the results is greatly affected by the significant elongation of the elements in the lateral direction (Heimbs, 2011). Advantages, disadvantages, and relevant enhancements relating to the ALE formulation are summarized in Table 6.3. Coupled Eulerian Lagrangian (CEL) method, which is an approach similar to the ALE method, which is an approach similar to the ALE method, has been presented and implemented by Ivancˇevic´ and Smojver (2011). The Eulerian model in CEL analyses is usually represented by a stationary cube containing Eulerian finite elements. ABAQUS provides multi-material EC3D8R volume elements to model Eulerian problems, which may be completely or partially occupied by the Eulerian material (Ivancˇevic´ & Smojver, 2011). Besides the FE method (Lagrangian), finite volume method (Eulerian), and their combination (ALE), some numerical methods model the bird torso using discrete nodal points, which are known as discrete element methods (DEMs). These methods are very advantageous in situations in which the deformations are very large. Since the elements are not attached to each other, large distortion of the fluid does not pose a problem. The discrete element motions can be described using: (a) single rigid primitives (e.g. spheres); (b) single deformable bodies (which are also discretized using the FE method); and (c) single deformable bodies the deformations of which are calculated in an additional subroutine (Petrinic & Duffin, 2000). The two important DEMs are the SPH and the NM methods, which are going to be described in the following paragraphs. Table 6.3 Advantages and disadvantages of the ALE formulation in bird-strike application Advantages

Disadvantages

Enhancements

Smaller grid structure size compared to the Eulerian method

Element distortion can occur, which causes negative volumes

Finer mesh

Less element distortion compared to the Lagrangian method

Inaccurate results in high deformations

Finer mesh

More accurate contact modeling compared to the Eulerian method

Filtration of some portions of dynamic loading in EulerianLagrangian interaction

Finer mesh

Relatively good time-history data

Finite element bird-strike modeling

6.3.4

127

Smoothed particles hydrodynamics (SPH)

A fluid can be represented by several separate particles (but in interaction with each other) provided that the size of particles are large enough to include sufficient number of molecules so that the fluid properties can be considered uniform in each particle, and small enough to be able to show the gradual change of macroscopic fluid properties. Each SPH element is given a mass, the amount of which is determined by dividing the total mass of the fluid by the total number of the SPH elements. In addition to mass, each SPH element carries hydrodynamic and thermodynamic information of the fluid at that point (Guida, 2008). In the SPH method, the particles are described by: ðxi ðtÞ,mi ðtÞÞi2P

(6.9)

where P is the set of moving particles and xi(t) and mi(t) are, respectively, the position and mass of the moving particle i. The mass change of each particle is related to its movement by (Lacome, 2000): dmi ¼ r:V ðxi , tÞmi dt

(6.10)

where V is the velocity vector of the particle i. The quadrature formula can be written as: ð X   f ðxÞdx ¼ mj ðtÞf xj ðtÞ (6.11) Ω

j2P

Smoothing kernel is a valuable concept in the SPH technique. Before presenting the smoothing kernel, the auxiliary cubic B-spline function must be introduced: 8 3 3 > > 1  y2 + y3 y  1 > > < 2 4 θðyÞ ¼ α1  1 3 1y2 ð2  yÞ > > > > :4 0 y2

(6.12)

where α1 is a constant depending on the dimension and the slope of kernel function. In 10 2D problems: α1 ¼ . Lacome (2000) suggested the smoothing kernel as: 7π     1 xi  xj
W xi  xj , h ¼
θ (6.13) h h
where h
is the smoothing length (see Fig. 6.4). In initial SPH studies, the smoothing length was considered identical for all the particles and was unvaried throughout the fluid deformation. However, later studies revealed that it is better to consider different smoothing lengths for different particles depending on the number of particles close to

128

Bird Strike

Figure 6.4 Active domain around a particle in the SPH method (Grimaldi, 2011).

2hk particle, k

them (Alberto, 2006). The idea of variable smoothing length was first proposed and used by Benz (1988). The main concept of this idea is that it is necessary to keep enough but not too many numbers of particles in interaction with a particle. Therefore in very sparse locations, it is better to increase the smoothing length so that the change of variables in the fluid can be considered continuous. On the other hand, in highly populated regions, the interaction of a very large number of particles does not provide much additional accuracy compared to a reasonable lower number of interacting particles. In summary, the smoothing length is variable in time and space to avoid the numerical fracture caused by material expansion and the lengthy run times caused by material compression. Therefore, in the above equations, h
can be replaced by h(xi) (Goyal, 2013). Generally, a property A(xi) of the particle i is represented by its smoothed value Ah(xi) and is obtained by approximating the integral given in Eq. 6.11 (Alberto, 2006):    A xj  mj   W xi  xj , h
A ðx i Þ ¼ ρ xj j¼1 h

N X

(6.14)

The gradient of the property is found by applying the derivation operator on the smoothing length of Eq. 6.14. rAh ðxi Þ ¼

    A xj mj   rW xi  xj , h
ρ xj j¼1

N X

(6.15)

By replacing the general property A(xj) in Eq. 6.14 by mass density, the mass density of the fluid at point i is given by (Monaghan, 1989): ρðxi Þ ¼

N X   mj W xi  xj , h
j¼1

(6.16)

Finite element bird-strike modeling

129

The general equation of conservation of mass is (Monaghan, 1989): dρ ðxi Þ ¼ ρrV dt

(6.17)

On the other hand, by obtaining the gradient of mass density using Eq. 6.15 and inserting it in Eq. 6.17, the SPH approximation of conservation of mass can be obtained as (Monaghan, 1989): N X     dρ ðx i Þ ¼ mj v xj  vðxi Þ rWij dt j¼1

(6.18)

N X dρ ρ    ðx i Þ ¼ mj v xj rWij dt ρ j j¼1

(6.19)

or

The conservation of momentum for the SPH formulation may be written as: ( )   N X P xj dv Pð x i Þ ðxi Þ ¼ mj rWij  2 rWji dt ρ2i ρj j¼1

(6.20)

The neighbor search is a primary step in SPH simulations and requires high attention. It is important to determine what particles are in interaction with a particular particle at each time step (Hedayati, Sadighi, & Mohammadi-Aghdam, 2014). All the particles are given an influence domain which is a sphere with radius 2h. In the neighbor search step, all the neighbor particles that lie inside the influence domain of a particle at the time step are listed. In the search for finding the influencing particles for a set of N particles, the distance between the particle and N  1 other particles must be checked. For the total number of N particles, a total number of N ðN  1Þ distance calculations are needed, which can take a very long time for large models. To overcome this, the space covered by the particles is split into several stacked boxes (known as buckets). The neighbor particle search for each particle only takes place inside the bucket containing the particle and also its neighbor buckets. This algorithm is called bucket sort. After the hypothetical neighbors of a particle from the main bucket and its neighbor buckets are completely listed, the distances between the corresponding particle and the other particles are checked to see if the distances are lower than twice the smoothing length. Unlike the primary neighbor search in which the computational cost was proportional to N2, the total cost of the neighbor search using “bucket sort” is linear with the number of particles N (Lacome, 2001): The time rate of change for the smoothing length is: dh 1
¼ hr:V dt 3

(6.21)

130

Bird Strike

This equation is valid only if the mass of the particles in the neighborhood are kept constant. A sphere with radius of 2h including n particles has the mass of: 4 M ¼ n:m ¼ n:ρV ¼ n:ρ: π8h3 3

(6.22)

and the time rate of mass change is: dM dρ 32 32 dh3 ¼ n: : πh3 + n ρ π dt dt 3 3 dt

(6.23)

Due to the grid-less nature of the SPH method, the problems arising from severe mesh distortion or element tangling associated with the Lagrangian (and less severely with the ALE) methods do not exist. The time step is therefore kept constant throughout the fluid deformation (Heimbs, 2011). Moreover, due to the Larangian nature of both the bird model and the target, dealing with the fluid–solid coupling problem is much easier than the ALE method even with a much lower number of elements. Modeling fragmentation of the fluid particles in this method is easier than in both the Lagrangian and the ALE methods. Similar to the Lagrangian method, this method also provides tracking of the time history of material properties at different material points (Grimaldi, 2011). Nevertheless, this method also has some disadvantages. First, since in the SPH simulations, the interaction of a large set of particles must be checked at each time step (Audic, Berthillier, Bonini, Bung, & Combescure, 2000), the simulations are computationally demanding and require both high memory and CPU resources. The other problem associated with the SPH method is a lack of sharp boundaries (Castelletti & Anghileri, 2004) which makes it hard to apply boundary conditions especially after the deformations have started and particles are spread. The diffused boundaries can also affect the definition of fluid–solid interaction. In LS-DYNA SPH modeling, usually the contact type node-to-surface is used for defining the interaction between the bird and the target, where a set of nodes corresponding to the bird’s particles are defined as the slave and the target surface is defined as the master counterpart. The other problem of the SPH method is called “tension instability,” in which the particles under negative pressures (tension) form unphysical clusters which cause numerical collapse (Dyka & Ingel, 1995; Heimbs, 2011). Advantages, disadvantages, and relevant enhancements related to the SPH formulation are summarized in Table 6.4.

6.3.5

Nodal masses (MN) method

The other important type of DEM approach is NM in which the bird model is represented by a set of separate nodes each given a mass. Unlike the SPH method, the nodes of the NM method are not in interaction with each other. In a study by Anghileri et al. (2005), the NM method predicted the behavior of the bird in an unrealistic way. Using a dissipation mechanism (inclusion of damping in the contact algorithm at the

Finite element bird-strike modeling

131

Table 6.4 Advantages and disadvantages of the SPH formulation in bird-strike application Advantages

Disadvantages

Enhancements

No problem associated with element distortion or dangling

High computational cost

Using parallel computing

Good ability in dealing with FSI, compared to ALE, and Eulerian methods

Lack of sharp boundaries

Good ability in modeling fragmentation

Tension instability

Tracking the properties of the material points is simple

interface between the bird and the impacted surface) led to better results. The other drawback of the NM method is that the energy absorption due to the deformation of the bird is neglected (Anghileri & Bisagni, 2000).

6.3.6

Comparison of the numerical approaches

A comprehensive set of simulations was carried out using bird models created by the three methods: Lagrangian, ALE, and SPH. For all the simulations, the geometry of the bird was a hemispherical-ended cylinder which was then impacted to rigid normal and oblique targets. (The modeling process is explained in detail in Chapter 8.) To evaluate the accuracy of the results of each approach, the results were compared to the experimental results provided by Wilbeck (1978). The time sequences of bird model deformation on the rigid target can be compared between the three approaches in Fig. 6.5 for normal impact and in Fig. 6.6 for oblique impact. Moreover, the effect of mesh size in the pressure profiles obtained from the three approaches is depicted in Fig. 6.7a–c. In most cases, refining the mesh has led to smaller Hugoniot pressures (Fig. 6.7a–c). For fine mesh (12 elements through bird radius), all the three approaches predicted results close to the experimental data (Hedayati & Ziaei-Rad, 2011b). In another numerical/experimental study, Lavoie et al. (2009) stated that the Lagrangian method is inappropriate for bird-strike modeling due to pressure loss, mass loss, and inaccurate radial pressure distribution. The above-mentioned problems were resolved by using refined mesh (this however increased solution time). Lavoie’s ALE model gave acceptable results for all the evaluation criteria provided that the suitable values for solution parameters were used. Lavoie’s SPH model also provided acceptable results which were comparable to the ALE results, but with a lowered number of solution parameters to be dealt with. Finally, they suggested the SPH method as an accurate yet simple approach to simulate bird-strike problems.

Lagrange

SPH

ALE

0.45 ms

0.9 ms

1.2 ms

ALE

SPH

Lagrange

Figure 6.5 Evolution of bird deformation on normal rigid targets for different numerical approaches (Hedayati & Ziaei-Rad, 2011b). (Reprinted by permission of the publisher (American Institute of Aeronautics and Astronautics, Inc.).)

0.6 ms

0.9 ms

1.2 ms

Figure 6.6 Evolution of bird deformation on oblique rigid targets for different numerical approaches.

Finite element bird-strike modeling

3.5

133

´ 107 Coarse mesh Medium mesh Fine mesh Experiment

3

Pressure (Pa)

2.5 2 1.5 1 0.5 0 0.2

0.4

0.6

0.8

(a)

1

1.2

1.4

1.6

1.8

Time (s) 3.5

2 ´ 10−3

´ 107 Coarse mesh Medium mesh Fine mesh Experiment

3

Pressure (Pa)

2.5 2 1.5 1 0.5 0 0.2

0.4

0.6

0.8

(b)

1

1.2

Time (s)

1.4

1.6

1.8

2 ´ 10−3

Figure 6.7 Effect of mesh size on the pressure profiles obtained from the center of impact for: (a) Lagrangian; (b) SPH; and (Continued)

6.4

Bird material modeling

As stated in Chapter 4, the behavior of different internal parts of a bird in low-velocity impacts is neither uniform nor homogenous. Modeling bird behavior in such a situation is very complex, and depends on many parameters, e.g. bird age, sex, species, point of impact, etc. However, by increasing the shooting velocity, the bird tissues do not have enough shear strength against the generated high pressures and behave

134

Bird Strike

3

x 107 Coarse mesh Medium mesh Experiment

Pressure (Pa)

2.5

2

1.5

1

0.5

0 0.2

0.4

(c)

0.6

0.8

1 1.2 Time (s)

1.4

1.6

1.8

2 x 10−3

Figure 6.7, cont’d (c) ALE formulations in a bird strike (Hedayati & Ziaei-Rad, 2011b). (Reprinted by permission of the publisher (American Institute of Aeronautics and Astronautics, Inc.).)

as a homogenous jet of fluid. In this case, almost all of the bird torso portions have identical mechanical behavior, and a single material model can effectively predict their behavior in different situations (Wilbeck, 1978; Meguid, Mao, & Ng, 2008). Since Wilbeck (1978) showed that the behavior of a real bird impact in rigid target tests is similar to the behavior of a water jet with 10% porosity, in many numerical studies, a polynomial equation of state (EOS) with material constants associated with 10% porous water is used (Hedayati & Ziaei-Rad, 2011b). For modeling the material behavior of the bird, usually a “material model” along with an “equation of state” is used. A number of material models, e.g. null, isotropic elastic plastic hydrodynamic (IEPH), etc., are prevalent for water at high velocities. The null material model relates the stress and strain of the bird by: σ ij ¼ Pδij + 2ργ e_ij

(6.24)

where δij and e_ ij are the identity and the rate-of-deformation tensors, respectively, P is the fluid pressure, and γ is the dynamic viscosity coefficient. Although the shear strength of the bird material is usually neglected in high velocities, some authors (Brockman & Held, 1991; Cassenti, 1979) believe that superimposing a deviatoric material response (i.e. 2ργ e_ij) to the hydrodynamic constitutive law (i.e. Pδij ) will improve the numerical–experimental correlation (Airoldi & Cacchione, 2006). The IEPH material model is well suited for bird-strike modeling because it behaves as an elastic-plastic material at low pressures, and after the occurrence of a high velocity impact, it is governed by the pressure–volume relationship of the EOS. Therefore, a

Finite element bird-strike modeling

135

low shear strength value can be given to the bird allowing it to retain its shape until the impact. In an IEPH material model, the yield stress is defined as: σ y ¼ σ 0 + Eh ε p + ða1 + pa2 Þmax ½p, 0

(6.25)

where a1 and a2 are linear and quadratic pressure hardening coefficients, respectively. εp is the effective plastic strain which is given by: ðt  ε ¼ p

2 p p D D 3 ij ij

1

=2

dt

(6.26)

0

in which Dpij are the plastic components of the rate-of-deformation tensor. In Eq. 6.25, Eh is the plastic hardening modulus that is defined as: Eh ¼

Et E E  Et

(6.27)

where E and Et are the elastic and tangent moduli of the bird material, respectively.

6.5

Equations of state (EOS)

EOS is a relationship that relates state variables of a material in a physical condition. EOSs are useful in describing the behavior of both the fluids (gasses or liquids) and solids. The state functions in an EOS are usually temperature, pressure, density (volume), and internal energy. The EOS of ideal gases is described by the well-known analytical relationship: PV ¼ nRT J . mol:K However, no analytical relationship has yet been derived for the solids and liquids, and the proposed EOSs for them are semi-empirical formulas based on experimental data (Nizampatnam, 2007). This is why several competing EOS formulas have been proposed by different researchers which prompted Zukas (2004) to say: “Grown men sometimes descend into childish arguments over the merits of these EOS.” For bird-strike modeling, several EOSs have been proposed including linear, polynomial, tabulated, and Mie-Gru¨neisen formulas; these are now to be described in the following paragraphs. where R is the ideal gas constant and is equal to 8:3144621

6.5.1

Linear EOS

The simplest EOS used for bird modeling is the one in which pressure is related to the density in a linear way (Guida et al., 2013):

136

Bird Strike



ρ P¼K 1 ρ0

 (6.28)

where P is the current pressure in the fluid, K is the bulk modulus, ρ is the material current density, and ρ0 is the reference density at which the material has no pressure. The bulk modulus, K, measures the resistance to uniform compression. It is defined as the ratio of pressure increase needed to cause a given decrease in relative volume. Its SI unit is Pascal and its value for water is 2200 MPa. Another type of linear EOS used for modeling the hydrodynamic response of the bird is the Murnaghan EOS which has been used in several bird-strike studies including Vignjevic, Orłowski, De Vuyst, & Campbell, (2013), Liu, Li, and Gao (2014), Guo, Jia, and Hong (2012), and has the following form:  γ  ρ 1 P ¼ P0 + B ρ0

(6.29)

where P0 ¼ 0 is the reference pressure. B and γ are material constants and have to be determined experimentally. The values of B ¼ 128 MPa and γ ¼ 7:98 have been provided by McCarthy et al. (2005) for water.

6.5.2

Polynomial EOS

The most commonly used EOS for the water-bird is a polynomial of degree 3. This polynomial EOS for the bird model corresponds to a hydrodynamic, isotropic, and non-viscous constitutive law, and is given as follows:   (6.30) P ¼ c0 + c1 μ + c2 μ2 + c3 μ3 + c4 + c5 μ + c6 μ2 Ei where c0–c6 are the coefficients of the polynomial equation, Ei is the internal energy, and μ is the change in density during the impact: μ¼

ρ 1 ρ0

(6.31)

The coefficients are given by expressions based on the initial density, the speed of sound in the medium, and an experimental constant k. With a known and negligible initial equilibrium pressure, the values of the coefficients in Eq. 6.30 are given as follows: c1 ¼ ρ0 c20 c2 ¼ ð2k  1Þc1 c3 ¼ ðk  1Þð3k  1Þc1 c0 ¼ c4 ¼ c5 ¼ c6 ¼ 0

(6.32)

where ρ0 is the density of the medium (for the water ρ0 ¼ 1000 kg=m3 Þ, c0 is the speed of the sound in the medium (for water: c0, w ¼ 1483 kg=m3 ), and k is an experimental constant (for water k ¼ 2) (Ugrcˇicˇ, 2012).

Finite element bird-strike modeling

137

The coefficients of the polynomial equation (i.e. Eq. 6.30) have been obtained for porous and non-porous water in Brockman and Held (1991) and their validity in a birdstrike event has been evaluated for the porous models presented in Wilbeck (1978) and Cassenti (1979). Figure 6.8 compares the Hugoniot pressures of bird models with different apparent densities (or specific volumes) and with polynomial EOS with three sets of constants corresponding to the porosities: α ¼ 0; α ¼ 0:1; and α ¼ 0:15, taken from Wilbeck (1978) and Cassenti (1979). The corresponding values of ci are listed in Table 6.5 for different porosities. Comparison of the Hugoniot pressures of the numerical models having different degrees of porosities to their corresponding porous water counterparts calculated according to the methods described in Meyers (1994) 6e+08 Water Polynomial law for bird material with a = 0.00 Water, a = 0.10 Polynomial law for bird material with a = 0.10 Water, a = 0.15 Polynomial law for bird material with a = 0.15

5e+08

Pressure (Pa)

4e+08

3e+08

2e+08

1e+08

0

0.75

0.8

0.85 0.9 0.95 Relative specific volume

1

1.05

1.1

Figure 6.8 Hugoniot compressive curves of water-like materials and homogenized bird materials (Airoldi & Cacchione, 2006). (Reprinted by permission of the publisher (Elsevier).)

Table 6.5

values of ci in linear EOS for different porosities of water

α¼0

C1 ¼ 2250 Mpa C0 ¼ C2 ¼ C3 ¼ 0

α ¼ 0:1

C0 ¼ 0 C1 ¼ 511:7 Mpa C2 ¼ 8224 Mpa C3 ¼ 55150 Mpa

α ¼ 0:15

C0 ¼ 0 C1 ¼ 748:4 Mpa C2 ¼ 9622:1 Mpa C3 ¼ 36120 Mpa

138

Bird Strike

indicates that the interpolations can well predict the shock state of water-like materials with different porosities (Airoldi & Cacchione, 2006).

6.5.3

Tabulated EOS

The EOS for the bird material in the tabulated form has the form: P ¼ f ð εV Þ

(6.33)

where εV is the volumetric strain given by the natural logarithm of the relative volume V. The values of the tabulated EOS are listed in Table 6.6.

6.5.4

Mie-Gru¨neisen EOS

In 1903, Gustav Mie developed an intermolecular potential for deriving hightemperature EOSs for high-temperature solids (Mie, 1903). In 1912, Eduard Gru¨neisen extended Mie’s model to temperatures below the Debye temperature at which quantum effects become important (Gru¨neisen, 1912). This equation (also called Us  Up EOSs) describes a linear relationship between the shock and particle velocities. Gru¨neisen’s EOS with cubic shock velocity defines pressure for materials in compression as:

P¼"

 γ a ρ0 C2 μ1 + 1  0 μ  μ2 2 2 μ2 μ3  S3 1  ðS1  1Þμ  S2 μ+1 ðμ + 1Þ2

#2 + ðγ 0 + aμÞ E

(6.34)

and in tension: P ¼ ρ0 C2 μ + ðγ 0 + aμÞE

(6.35)

Table 6.6 Constants of tabulated equation of state for water with 10% porosity (Marulo & Guida, 2014)

1 2 3 4 5 6 7 8 9 10

εV

f

0.000 –0.105 –0.118 –0.128 –0.137 –0.154 –0.169 –0.183 –0.195 –0.217

0.000 2.37e8 4.25e8 5.86e8 7.27e8 9.72e8 1.18e9 1.37e9 1.54e9 1.84e9

Finite element bird-strike modeling

139

where C is the intercept of the vs  vp curve; S1, S2, and S3 are the coefficients of the slope of the vs  vp curve; γ 0 is the Gru¨neisen gamma; a is the first order correction to γ 0, and μ ¼ ρ=ρi  1 where ρ is the material density (LS-DYNA Keyword User’s Manual, v.9.71, 2006). The Mie-Gru¨neisen EOS’s parameters for water are C ¼ 1480; S1 ¼ 1:92; S2 ¼ 0; S3 ¼ 0; and γ 0¼0.1 (Chizari, Barrett, & Al-Hassani, 2009). Using the appropriate EOS is important. Zukas (2004) emphasizes that Mie-Gru¨neisen EOS is only applicable to solid materials that remain in the solid state throughout the impact event; that is, no phase change is allowed (Nizampatnam, 2007). However, many bird strike studies used this EOS for their numerical models.

6.6

Fluid–structure interactions

One of the main complications in modeling the impact problems in FE codes is proper modeling of the interaction between the projectile and the impacted structure. If the projectile is stiff and strong enough compared to the impacted structure (e.g. in the impact of a steel ball to a thin composite plate), its mechanical behavior in the FE code can be defined as rigid, due to the fact it does not suffer considerable deformation during impact. In that case, the impacted structure is penetrated locally and penetration depth can be calculated using the Hertizan contact law assuming the quasi-static behavior for the materials (Hou & Ruiz, 2007). However, the bird material is usually more deformable and weaker compared to the impacted structure. Moreover, due to the rapid deceleration at the point of impact, the material response of the bird models can be treated as a fluid. As a result, the interaction between the bird and the target surface can be regarded as a fluid–structure interaction. When a bird impacts a deformable target, the loads applied by target to the bird model depends on the stiffness and deformation of the target structure, and also the response of the bird to the generated loads. Therefore, proper simulation of interaction between the bird and the target is a crucial factor in prediction of the deformation and permanent damage of the target structure (Hedayati & Ziaei-Rad, 2011a). Different contact algorithms are implemented for the Lagrangian, ALE, and SPH approaches to bird strike. The algorithms which govern the interaction between the bird and the target are usually defined in different ways. The algorithms governing the bird/target interaction in Lagrangian and SPH approaches are called the “contact” algorithms. These algorithms are usually based on the penalty method which considers the differences of the mechanical properties of the target and the impactor. For the ALE approach, the interaction algorithm is usually called the coupling algorithm which is also based on the penalty method (Castelletti & Anghileri, 2008).

6.6.1

Contact algorithms for Lagrange

The interaction of the Lagrangian bird and the target is usually represented using a surface-to-surface or node-to-surface contact algorithm which are based on a penalty method whereby the extent of penetration of the slave set into the master set is determined by the penalty stiffness (Chuan, 2006). The penalty stiffness corresponds to an

140

Bird Strike

hypothetical spring acting between the slave and master surfaces to prevent their penetration into each other. The penalty stiffness of the slave and master segments can be obtained using the relationship: K slave : PSslave ¼ K master : PSmaster

(6.36)

where Kslave and Kmaster are the bulk moduli of the slave and master parts, respectively, and PSslave and PSmaster represent the penalty stiffness of the slave and master segments. Minor penalty stiffness usually leads to the excessive penetration of a slave segment into the master segment and therefore elimination of the contact. In theory, very high contact stiffness is favorable. However, very high penalty stiffness leads to ill-conditioning of the stiffness matrices which can cause numerical errors. The contact surfaces are usually defined using sets, which can include segments (element faces), nodes, beam elements, shell elements, polyhedral elements, and parts (Oasys LS-DYNA environment: 8.1, 2001). To model the bird contact with the target surface, the node-to-surface contact algorithm is the best choice. In this contact mode, the bird nodes are defined as the slave set, and the target surface (or body) is specified as the master set. During the solution, the nodes of the slave set are checked for not penetrating the master segment (Fig. 6.9). The other prevalent contact type for the Lagrangian bird model is the surface-tosurface contact mode (Fig. 6.9). In this contact type, it is not important which of the bird or the target are defined as slave or master sets. The only important point is that of the two surfaces in contact: one must be defined as the master that has a coarser mesh. In the surface-to-surface contact mode, at each time step, the nearest node of the slave surface to the master surface is found by checking the normal distance of each of the slave nodes to the target surface segments. If the nearest node can

Figure 6.9 Schematic view of node-to-surface contact type.

Finite element bird-strike modeling

141

be deemed “close” to the master segment according to the corresponding criteria, the node is translated to the master surface (Grimaldi, 2011). After that, the slave nodes are not allowed to travel toward the other side of the master segments. Shmotin et al. (2009) investigated the influence of the friction coefficient in the contact calculation between the bird and the metallic target structure with values from 0 to 1. Best results, compared to experimental results, were obtained with zero friction.

6.6.2

Contact algorithms for SPH

For the SPH bird models, there are two main contact algorithms for modeling the interactions in the model: particle-to-surface and particle-to-particle (Vignjevic et al., 2013). In the particle-to-particle contact algorithm, the distances of neighbor particles are checked to disallow the neighbor particles to penetrate each other. In this penalty-based method, the distance between each pair of particles is checked, and according to the distance, a contact force is applied to both of the particles along the hypothetical line connecting the centers of the two particles (Fig. 6.11a) (Vignjevic et al., 2013). In the particle-to-surface contact algorithm, at each time step, the normal distance between each of the SPH particles and the target surface segment is checked, and if it is larger than 0.5h, the contact has occurred (Fig. 6.11b). After the initiation of contact, a restoring penalty force is applied to the SPH particles as well as the target surface nodes. The direction of force at all the positions is perpendicular to the target surface at that point (Vignjevic et al., 2013).

Figure 6.10 Schematic view of surface-to-surface contact type.

142

Bird Strike

(a)

(b) Figure 6.11 SPH contact types: (a) particle-to-particle, (b) particle-to-surface (Vignjevic et al., 2013). (Reprinted by permission of the publisher (Elsevier).)

6.7

Hourglass control

Hourglass modes are nonphysical modes of deformation that produce zero strain and no stress. Hourglass modes occur only in under-integrated (single integration point) solid, shell, and thick shell elements (HOURGLASS, 2010). These modes are called hourglass modes because of the shape of the shell elements in 2D problems which resemble an hourglass. These modes also are given other names, e.g. keystone modes or spurious zero-energy modes (Nizampatnam, 2007). Consider a rectangular element as shown in Fig. 6.12. While the element has deformed to a great extent, the length of the perpendicular bisector to the sides of the elements is still unvaried. Therefore, the FE code considers no compressive/ tensional strain for the element in any direction. While some energy has been spent to cause such deformation, the internal energy of the element has not changed (due to no change in strains). This leads to “energy leak.”

X mode

Y mode

Figure 6.12 Hourglass modes of under-integrated solid elements.

Finite element bird-strike modeling

143

One way to avoid the formation of hourglass elements is by utilizing viscous damping or small elastic stiffness (e.g. by activating bulk viscosity control in LS-DYNA). However, this action might have negative effects on the stable global mode (Grimaldi, 2011) and may increase the solution time multiple times. Since the modes of hourglass deformation are orthogonal to the strain, the work done to avoid formation of hourglass elements might be neglected in the energy equation. However, the hourglass must always be controlled if under-integrated solid elements are used in the model (Grimaldi, 2011). In the absence of contact friction, the hourglass energy should not usually exceed 10% of internal energy to achieve reliable results in the simulation (Mao, Meguid, & Ng, 2009).

6.8 6.8.1

Bird geometry modeling Traditional bird models

The real birds have body shapes completely different from simple geometries. The inner part of the bird bodies is also very complex consisting of organs with different shapes and materials (Hedayati & Ziaei-Rad, 2013b). However, as the bird behaves like a fluid in a high velocity impact, and on the other hand, the density of all the bird organs are almost close (in the range of 950 kg/m3), in most bird-strike studies the bird material properties are considered uniform. Moreover, to simplify the bird modeling procedure and having consistent projectiles, several simple geometries have been proposed for the FE modeling of the birds (e.g. hemispherical-ended cylinder (Frischbier, 1997; Langrand et al., 2002; Hedayati & Ziaei-Rad, 2011b; McCarthy et al., 2004; Hedayati & Jahanbakhshi, 2015; McCarthy et al., 2005; Airoldi & Cacchione, 2006), ellipsoid (Richard, 2000; Hedayati & Ziaei-Rad, 2013a), and straight-ended cyclinder (Brockman & Held, 1991; Salehi et al., 2010). The hemispherical-ended cylinder geometry is the most common geometry, whilst the straight-ended cylinder is the less common one. Meguid et al. (2008) focused on the three configurations mentioned above at various length-to-diameter aspect ratios. The results showed that the aspect ratio is not very effective on the pressure profile or impulse diagrams. Moreover, they concluded that the initial contact area between the bird and the target has a very significant effect on the Hugoniot pressure value. The highest Hugoniot pressures belonged to the straight-ended cylinder, the hemispherical-ended cylinder, and the ellipsoid, respectively. Hedayati and Ziaei-Rad (2013a) modeled four SPH bird models (Fig. 6.13) namely sphere, straight-ended cylinder, hemispherical-ended cylinder, and ellipsoid, and impacted them to rigid targets from their axial and lateral sides. The results from the axial side impact were compared between the four geometries (and the experimental tests by Wilbeck (1978)) to find the most appropriate bird substitute geometry. Comparison of the four geometries was also made in the lateral impacts. The length-to-diameter ratio was set to two for all the geometries except the spherical impactor. Some shell elements (which played the role of pressure sensors) were attached to the center of impact using tied-node-to-surface contact type, and with

144

Bird Strike

Figure 6.13 Various geometries for the SPH bird model (Hedayati & Ziaei-Rad, 2013a). (Copyright © 2013 Elsevier Masson SAS. All rights reserved.)

specific distances to measure the pressure profile at different radii of the impacted surface. The pressure profiles were obtained by dividing the contact force profile for each shell element (or sensor) to its area. Deformation of different bird configurations at different time intervals is shown in Fig. 6.14. The deformations of the different models are very different at the initial moments of contact, but after a while they become relatively similar. In Fig. 6.15, pressure at the center of impact is compared between the Wilbeck experimental result (Wilbeck, 1978) and the four substitute bird configurations. The straightended cylinder (190 MPa), the sphere (168 MPa), the ellipsoid (120 MPa), and the hemispherical-ended cylinder (80 MPa) models had the higher Hugoniot pressures, respectively. Interestingly, the Hugoniot pressure obtained for the straight-ended FE model was very close to the theoretical Hugoniot pressure (198 MPa) which is also based on the straight-ended cylinder geometry (see Eq. 6.5, and Chapter 4, Fig. 4.5). The hemispherical-ended cylinder had the closest result to the experimental test by Wilbeck (1978), although it was still four times of that. As Wilbeck (1978) states the pressure transducers might have been unable to capture the very short time Hugoniot pressures of the impact due to their bandwidth limitations. An additional reason may be the difference between the initial contact area of the hemispherical-ended cylinder and the real bird. A bird model with a more realistic geometry would be helpful, which will be discussed in the following paragraphs. The peak pressure distribution on the target plate is compared between the four geometries in their axial impact in Fig. 6.16. The trend is similar for all the cases, i.e. the peak pressure is maximum at the center, and decreases by increase in the distance

Finite element bird-strike modeling

145

Figure 6.14 Comparison of the deformations of different bird geometries in the axial impact. (Reproduced from: Hedayati & Ziaei-Rad, 2013a. Copyright © 2013 Elsevier Masson SAS. All rights reserved.)

from the center of impact. As predictable, in straight-ended cylinder, a larger area of the target plate is affected by the very high pressure peaks. Therefore, it can be concluded that the straight-ended cylinder is the most damaging bird substitute in an impact.

6.8.2

Advanced geometry

The anatomic structure of birds includes several internal cavities, e.g. pneumatic bones, lungs, and especial air sacs that contribute to the complex bird’s makeup. Due to the reasons stated beforehand, researchers typically use various primitive

146

Bird Strike

2

x 108 Wilbeck (1978) experiment Ellipsoid Hemispherical−ended cylinder Sphere Straight−ended cylinder Reallistic bird (duck)

1.8 1.6

Pressure (Pa)

1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

1

2

3

4 Time (s)

5

6

7 x 10−4

Figure 6.15 Pressure profile at the center of impact for the four bird geometries, and the Wilbeck (1978) experiment. (Reproduced from: Hedayati & Ziaei-Rad, 2013a. Copyright © 2013 Elsevier Masson SAS. All rights reserved.)

geometries to represent the torso of the bird in numerical simulations. However, there are a few works in which the researchers have implemented more complex bird geometries which further resemble real bird shapes. McCallum and Constantinou (2005) built a complex multi-material bird model representing a Canadian goose weighing 3.6 kg. They compared the results of this multi-material bird model to the results of a hemispherical-ended cylinder. They observed that for the multi-material bird model with its neck extended, the target was pre-stressed due to the impact of the bird’s head and neck prior to the torso. Based on their results, McCallum and Constantinou (2005) asserted that modeling bird organs, e.g. head, neck, wing, etc. in addition to the bird torso is necessary when modeling large birds such as a Canadian goose. They stated that consideration of other body parts in bird models can have significant effects on damage initiation of aircraft components which can determine the final failure situation. In a similar effort to use more advanced bird models, Nizampatnam (2007) considered three types of bird models for bird impact against rigid targets: first, a hemispherical-ended cylinder made up of two distinct materials randomly distributed throughout the bird model. Second, a hemispherical-ended cylinder with uniform texture for most of the bird torso but with high density portions representing bird bones, as well as low density portions to represent the lungs and air sacs. Third, a bird model with a geometry more similar to real birds which included features such as neck, head, wing, lung, bone structure, and torso. Among the three afore-mentioned bird models,

Finite element bird-strike modeling

147

the more realistic bird model led to better results. Moreover, the multi-material bird model with realistic geometry was capable of predicting the damage initiation imposed by the bird head into the impacted structure. In a more recent study, Hedayati and Ziaei-Rad (2013b) attempted to better simulate a real bird body shape using CT scanning (Fig. 6.17). They imaged an anesthetized

2

×108 Ellipsoid Hemispherical-ended cylinder Sphere Straight-ended cylinder Wilbeck experiment

1.8

Pressure (Pa)

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0

0.5

1

1.5

Distance from the center of impact (cm)

Figure 6.16 Distribution of peak pressure on the rigid target plate for the tail side impact of birds with different geometries. (Reproduced from: Hedayati & Ziaei-Rad, 2013a. Copyright © 2013 Elsevier Masson SAS. All rights reserved.)

Figure 6.17 Setup of the anesthetized mallard bird before being imaged in the CT scan device.

148

Bird Strike

mallard bird using a medical CT scan device and recreated the 3D geometry of the bird using SPH elements. To create the SPH bird model, first the internal cavities of the bird were ignored and a uniform bird model without any internal cavity was created (Fig. 6.18c). All the SPH elements were given the material properties of water. Then the SPH elements located in the internal cavities of the bird model were removed (Fig. 6.18d), and then replaced by new SPH elements possessing the material properties of air (Fig. 6.18e). The resultant bird model consisted of 41,685 SPH elements with water material properties each weighing 0.0191 g and 7499 SPH elements with material properties of air each weighing 23 μg. The mass of each SPH particle was determined by dividing the total mass of the mallard bird by the number of SPH elements. The 3D view of the SPH bird model when impacting a rigid target from its tail side is shown in Fig. 6.19. The sensors (shell elements) installed on the target surface are also visible in Fig. 6.19. The sensors were installed at positions of: –2.5 cm, 0 cm, 2.5 cm, and 5 cm, on the target to capture the pressure plots at those locations. The air element was given a null material model with a linear polynomial EOS (which is given in Eq. 6.30). Since the air can be considered an ideal gas, the ideal gas law can be applied to it: P¼ρ

R T M

(6.38)

Figure 6.18 The procedure of constructing a multi-material SPH bird with realistic geometry: (a) a CT image slice; (b) the CT image is checkered; (c) the SPH points are placed on the vertices of the checkered image inside the periphery of the mallard bird; (d) the SPH elements located on the cavities are removed; (e) the cavities are filled by the SPH elements with the material properties of air. (Reproduced from: Hedayati & Ziaei-Rad, 2013b. Reprinted by permission of Taylor & Francis.)

Finite element bird-strike modeling

149

Figure 6.19 The 3D view of the SPH bird model when impacting a rigid target from its tail side.

R kJ ¼ 0:287 for air, the M kgK relationship P ¼ 84:1 ρ will be attained. Equalizing P ¼ C0 + C1 μ + C2 μ2 + C3 μ3 and P ¼ 84:1 ρ results in: C0 ¼ C1 ¼ 97:5 and C2 ¼ C3 ¼ 0: Distribution of the Hugoniot and the steady pressures on the rigid target plate are the most useful criteria for demonstrating the validity of a bird model. Comparison of the distribution of steady pressures on the target showed the excellent agreement of the numerical results obtained from the mallard bird model and the experimental results (Wilbeck, 1978). The agreement between the hemispherical-ended cylinder and the experimental result was less satisfactory. For the Hugoniot pressure distribution (Fig. 6.20), the mallard bird model shows much closer results to the experimental data than the hemispherical-ended cylinder does. The closer result of the mallard bird model can be attributed to the similarity of its tail geometry to that of a real bird. Although the experimental tests are always done in such a way that the birds are impacted to the target plate by their tail side (because of more stability of the bird in the launch tube) in real world bird strikes however, birds impact aircraft with their different body parts (see Fig. 6.21). The impact of birds with different body parts (e.g. head, tail, wings, and bottom side) can have different damaging effects on the aircraft component. To investigate the effect of impact orientation, four impact scenarios were considered for the mallard bird model (Hedayati & Ziaei-Rad, 2013b). The pressure peak distribution for the four impact scenarios (e.g. head, tail, wing, and bottom side) of the mallard bird model, the traditional (hemispherical-ended) bird model, and Wilbeck experimental tests are compared in Fig. 6.22. All the cases, except the bottom-side impact of a mallard model, show a maximum value at the center of impact, and the pressure peak decreases as the distance from the center of impact Considering the temperature 20°C and having:

150

Bird Strike

10

× 107

9

FE mallard tail side FE traditional bird model Wilbeck experiment

Maximum pressure (Pa)

8 7 6 5 4 3 2 1 0

−0.02

−0.01

0

0.01

0.02

0.03

0.04

0.05

Distance from center of Impact (m)

Figure 6.20 Hugoniot pressure distribution of the mallard bird model, hemispherical-ended cylinder, and experimental tests by Wilbeck (1978) in the impact from the bird tail-side. (Reproduced from: Hedayati & Ziaei-Rad, 2013b. Reprinted by permission of Taylor & Francis.)

Figure 6.21 Different bird orientations in a typical bird-strike event scenario.

increases. For the bottom-side impact, the high pressurized area is distributed over a larger area that shows its higher danger in bird-strike events. The head-side impact and the tail-side impact are the second and third most dangerous scenarios for the mallard bird model, and the wing-side impact shows the lowest pressures. Interestingly, in the wing-side impact, the existence of wings has led to central pressure peaks even lower than other points of the target. Again, the tail-side impact of the mallard bird model

Finite element bird-strike modeling

12

151

×107

Pressure (pa)

10

8 Bottom Head Tail Wing Traditional Wilbeck

6

4

2

0

−0.02

−0.01

0

0.01

0.02

0.03

0.04

0.05

Time (s)

Figure 6.22 The pressure peak distribution on the target for impacts from different orientations. (Reproduced from: Hedayati & Ziaei-Rad, 2013b. Reprinted by permission of Taylor & Francis.)

shows much closer results to the experimental tests than the hemispherical-ended cylinder does.

6.9

Differences in pressure readings

As previously observed in Chapters 4 and 5, the theoretical Hugoniot pressures are much higher than the experimental counterparts given in Wilbeck (1978). However, the theoretical and experimental steady pressures are almost in good accordance. Wilbeck (1978) contributed this huge difference between theoretical and experimental Hugoniot pressures to the inability of pressure gauges to capture very short-term peak pressures. Another reason may also be the effect of simplifications made in the derivation of the theoretical values. After the advent of numerical results, it became possible to evaluate which of the two theoretical and experimental-read Hugoniot pressures are representative of the actual Hugoniot pressure. However, the problem was not solved because the normalized Hugoniot pressures obtained in different numerical studies were shown to be in a wide range (between the experimental and theoretical values). The numerical studies on bird impact on rigid targets are listed in Table 6.7 and their Hugoniot and steady pressure values are compared to the experimental and theoretical values (shown in Fig. 6.23). The high scatter in numerical values can be seen in the figure. For example, it can be seen in Fig. 6.23a that at the impact velocity of about 100 m/s, the numerical study by Lavoie et al. (2009) has predicted a normalized Hugoniot pressure almost triple that of the Hugoniot pressure predicted by Tho and Smith (2011). The numerical result scatter is also visible in

152

Bird Strike

Table 6.7 List of FE simulations on bird impact against a rigid target as simulated by different academic research studies

1

2

3 4 5 6 7 8 9 10

11

12

13

Johnson and Holzapfel (116 m/s) (Johnson & Holzapfel, 2003) Airoldi and Cacchione (200 m/s) (Airoldi & Cacchione, 2006) Jenq (116 m/s) (Jenq et al., 2007) Meguid (225 m/s) (Meguid et al., 2008) Lavoie (116 m/s) (Lavoie et al., 2009) Mao (225 m/s) (Mao et al., 2008) Tho and Smith (116 m/s) (Tho & Smith, 2011) Mao (225 m/s) (Mao et al., 2009) Lavoie (95 m/s) (Lavoie et al., 2009) Smojver and Ivancˇevic´ (116 m/s) (Smojver & Ivancˇevic´, 2010) Ivancˇevic´ and Smojver (116 m/s) (Ivancˇevic´ & Smojver, 2011) Nishikawa (100 m/s) (Nishikawa, Hemmi, & Takeda, 2011) Smojver and Ivancˇevic´, 2012

Normalized Hugoniot pressure

Normalized steady pressure

Type of pressure recording

7.2

1.4

Not stated

6

1.4

Not stated

6.5

1

Not stated

2.6

0.3

Averaging

14

0.9

2.6

0.3

Center of impact Averaging

5.5

1.1

Not stated

2.6

0.3

Averaging

14

1.2

12.6

1.58

Center of impact Not stated

9

1.58

Not stated

16

2.13

12

1.58

Averaging over the initial area Not stated

(Reproduced from: Hedayati, Sadighi, et al., 2014. Copyright © 2014 Elsevier Masson SAS. All rights reserved.)

the steady pressures, but with less severity (Fig. 6.23b). Hedayati, Sadighi, and Mohammadi-Aghdam (2014) found out the main causes of these huge differences between the numerical results, and also which of the two values obtained from pressure transducers and theory, are in fact closer to the actual values. The normalized Hugoniot pressures obtained from the mallard and the hemispherical-ended cylinder bird models are compared to the theoretical and experimental values in different velocities in Fig. 6.24. Before discussing Fig. 6.24, it must

Finite element bird-strike modeling

153

20 18

Normalized Hugoniot pressure

16 14

* 12 10 8

Theory 1 2 3 4,6,8 5 7 9 10 11 12 13 Wilbeck test

6

*

4 2 50

100

150

(a)

200

250

300

Initial velocity (m/s)

2.5 Theory 1 2 3 4,6,8 5 7 9 10,11,13 12 Wilbeck tests

Normalized steady pressure

2

1.5

* 1

0.5

0 0

(b)

50

100

150

200

250

300

350

400

Initial velocity (m/s)

Figure 6.23 Comparison of: (a) normalized Hugoniot pressures; and (b) normalized steady pressures between the theoretical and experimental values, and the numerical results of the research listed in Table 6.7. (Reproduced from: Hedayati, Sadighi, et al., 2014. Copyright © 2014 Elsevier Masson SAS. All rights reserved.)

154

Bird Strike

22 Theory FE cylinder model (shell sensor) FE mallard model (shell sensor) FE cylinder model (averaging) Wilbeck experiment

20

Normalized Hugoniot pressure

18 16 14 12 10 8 6 4 2 50

100

150 200 Initial velocity (m/s)

250

300

Figure 6.24 Comparison of normalized Hugoniot pressures for the hemispherical-ended cylinder and the mallard FE models with the experimental and theoretical values. (Reproduced from: Hedayati, Sadighi, et al., 2014. Copyright © 2014 Elsevier Masson SAS. All rights reserved.)

be stated that different methods of measuring the pressure at the center of impact have been used in different numerical studies. In the first method, the contact force of the bird impact with the target surface is obtained and divided into the initial cross-section area of the bird model. The second method is to divide the contact force plot into the contact area between the bird and the target at any time. These two methods are called “averaging methods” in this book, because the obtained pressure is in fact the average of real pressure distribution imposed by the bird to the target. The averaging methods do not consider the severe pressure gradient at the point of impact. In the third method (implemented by Hedayati, Sadighi, et al., 2014), the contact force diagram between the bird and the “sensor surface” is obtained and divided by the area of the sensor itself. As can be seen in Fig. 6.24, the final method of pressure measurement has resulted in a good agreement of the Hugoniot pressures obtained from a hemispherical-ended cylinder to the theoretical values (if the pressure is read from sensors). The Hugoniot pressure of the mallard model is close to the experimental result, but the Hugoniot pressure of the hemispherical-ended cylinder is close to the experimental result, only if the pressure is obtained by averaging methods. However, using averaging methods for obtaining the pressure profile is not recommended because the actual Hugoniot pressure imposed by the hemispherical-ended cylinder is much larger.

Finite element bird-strike modeling

155

Theory FE cylindrical model (shell sensor) FE mallard model (shell sensor) Wilbeck experiments FE cylindrical model (averaging)

Normalized steady pressure

2

1.5

1

0.5

0 100

120

140

160

180 200 220 240 Initial velocity (m/s)

260

280

300

Figure 6.25 Comparison of normalized steady pressures for the hemispherical-ended cylinder and mallard FE models with experimental and theoretical values. (Reproduced from: Hedayati, Sadighi, et al., 2014. Copyright © 2014 Elsevier Masson SAS. All rights reserved.)

For the steady pressure, the situation is different. The obtained steady pressure obtained from the hemispherical-ended cylinder and the mallard bird model are both close to the experimental results (if the pressure is read from the shell sensors: see Fig. 6.25). Moreover, the difference between the steady pressures obtained by the two methods (sensors and averaging) for the hemispherical-ended cylinder is relatively lower (Fig. 6.25) than the difference observed between the Hugoniot pressures of the two methods (Fig. 6.24). This is due to the fact that in the steady regime, the pressure distribution is more uniform on the target, and as a result the averaging method does not lead to a very different steady pressure value.

6.10

Similarity law for bird strike

Employing small-scale experiments to replace full-scale bird-strike experimental tests has great advantages for reducing experimental costs and testing time. Yulong, Yongkang, & Pu (2008) investigated the model similarity law of aircraft structures under bird impact using dimensional analysis and similarity theory. Finite element models constructed with different scaling factors were numerically analyzed using ANSYS/LS-DYNA software. They concluded that for strain-rate-insensitive materials, the numerical results of the small-scale model and the full-scale model correlate well, and the similarity law can be applied to the bird-impact process. However, in high velocities, in which strain rate effects are determinative, using the similarity law proposed by Yulong et al. (2008) is not useful.

156

Bird Strike

References Airoldi, A., & Cacchione, B. (2006). Modelling of impact forces and pressures in Lagrangian bird strike analyses. International Journal of Impact Engineering, 32(10), 1651–1677. Alberto, C. (2006). Robust bird-strike modeling using LS-DYNA. Machine Design at University of Puerto RicoMayaguez Campus, Technical Comunicator at Ferreyros SAA, Peru. Anghileri, M., & Bisagni, C. (2000). Specific problems related to simulation of a bird impact against a turbofan inlet. In: Proceedings of the international crashworthiness conference (pp. 652–662). Anghileri, M., Castelletti, L. M., & Mazza, V. (2005). Birdstrike: Approaches to the analysis of impacts with penetration. In Impact loading of lightweight structures (pp. 63–74). Audic, S., Berthillier, M., Bonini, J., Bung, H., & Combescure, A. (2000). Prediction of bird impact in hollow fan blades. AIAA2000-3201, 1–7. Baughn, T. V., & Graham, L. W. (1988). Simulation of a birdstrike impact on aircraft canopy material. Journal of Aircraft, 25(7), 659–664. Belytschko, T., Privitzer, E., Mindle, W., & Wicks, T. (1979). Computer simulation of canopy-pilot response to bird-strike. Evanston, Illinois: Northwestern University, Civil Engineering. Benz, W. (1988). Applications of smooth particle hydrodynamics (SPH) to astrophysical problems. Computer Physics Communications, 48(1), 97–105. Birnbaum, N. K., Francis, N. J., & Gerber, B. I. (1999). Coupled techniques for the simulation of fluid-structure and impact problems. Computer Assisted Mechanics and Engineering Sciences, 6(3/4), 295–312. Blair, A. (2008). Aeroengine fan blade design accounting for bird strike. Dissertation, The University of Toronto. Brockman, R. A., & Held, T. W. (1991). Explicit finite element method for transparency impact analysis (No. UDR-TR-90-114). University of Dayton. Cassenti, B. (1979). Hugoniot pressure loads in soft body impact. In Proceedings of the 20th structural dynamics, and materials conference (pp. 241–248). Castelletti, L. M., & Anghileri, M. (2004). Multiple birdstrike analysis a survey of feasible techniques. In 30th European rotorcraft forum (pp. 495–505). Castelletti, L. M. L., & Anghileri, M. (2008). Birdstrike onto structures in rotational motion. In 26th international congress of the aeronautical sciences (pp. 14–19). ICAS. Chizari, M., Barrett, L. M., & Al-Hassani, S. T. S. (2009). An explicit numerical modeling of the water jet tube forming. Computational Materials Science, 45(2), 378–384. Chuan, K. C. (2006). Finite element analysis of bird strikes on composite and glass panels. Doctoral dissertation, National University of Singapore. Dyka, C. T., & Ingel, R. P. (1995). An approach for tension instability in smoothed particle hydrodynamics (SPH). Computers & Structures, 57(4), 573–580. Engblom, J. J. (1980). Coupled fluid/structure response predictions for soft body impact of airfoil configurations. In Emerging technologies in aerospace structures design structural dynamics and materials. ASME century 2 (pp. 209–223). Faure, J. M. (2011). Airbus dynamic analysis experiences using Abaqus/Explicit. In: 2011 SIMULIA customer conference, Barcelona, Spain. Frischbier, J. (1997). Bird Strike Capability of a Transonic Fan Blisk. In: Proceedings of the ASME Turboexpo, Orlando, FL, USA. Goyal, V. H. (2013). Smooth particle hydrodynamics for bird-strike analysis using LS-DYNA. American Transactions on Engineering & Applied Sciences, 2(2), 57–81.

Finite element bird-strike modeling

157

Grimaldi, A. (2011). SPH high velocity impact analysis: A birdstrike windshield application. (Doctoral dissertation, Università degli Studi di Napoli Federico II). Grimaldi, A., Sollo, A., Guida, M., & Marulo, F. (2013). Parametric study of a SPH high velocity impact analysis – A birdstrike windshield application. Composite Structures, 96, 616–630. Gru¨neisen, E. (1912). Theorie des festen Zustandes einatomiger Elemente. Annalen der Physik, 344(12), 257–306. Guida, M. (2008). Study, Design and Testing of Structural Configurations for the bird-strike compliance of aeronautical components. In University of Naples “Federico II”. Guida, M., Marulo, F., Meo, M., Grimaldi, A., & Olivares, G. (2011). SPH-Lagrangian study of bird impact on leading edge wing. Composite Structures, 93(3), 1060–1071. Guida, M., Marulo, F., Meo, M., & Russo, S. (2013). Certification by birdstrike analysis on C27J fullscale ribless composite leading edge. International Journal of Impact Engineering, 54, 105–113. Guo, Y., Jia, P., & Hong, G. (2012). Research on bird strike simulation of composite leading edge. AASRI Procedia, 3, 674–679. Hedayati, R., & Jahanbakhshi, M. (2015). Finite element analysis of an aluminum airplane stabilizer against birdstrike. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 1–10. Hedayati, R., Sadighi, M., & Mohammadi-Aghdam, M. (2014). On the difference of pressure readings from the numerical, experimental and theoretical results in different bird strike studies. Aerospace Science and Technology, 32(1), 260–266. Hedayati, R., & Ziaei-Rad, S. (2011a). Effect of impact orientation on bird strike analysis. International Journal of Vehicle Structures & Systems, 3(3). Hedayati, R., & Ziaei-Rad, S. (2011b). Foam-core effect on the integrity of tailplane leading edge during bird-strike event. Journal of Aircraft, 48(6), 2080–2089. Hedayati, R., & Ziaei-Rad, S. (2012). New bird model for simulation of bird strike on various layups used in transparent components of rotorcrafts. Journal of Aerospace Engineering, 27(1), 76–85. Hedayati, R., & Ziaei-Rad, S. (2013a). A new bird model and the effect of bird geometry in impacts from various orientations. Aerospace Science and Technology, 28, 9–20. Hedayati, R., & Ziaei-Rad, S. (2013b). Effect of bird geometry and orientation on bird-target impact analysis using SPH method. International Journal of Crashworthiness, 17(4), 445–459. Hedayati, R., Ziaei-Rad, S., Eyvazian, A., & Hamouda, A. M. (2014). Bird strike analysis on a typical helicopter windshield with different lay-ups. Journal of Mechanical Science and Technology, 28(4), 1381–1392. Heimbs, S. (2011). Computational methods for bird strike simulations: A review. Computers & Structures, 89(23), 2093–2112. Hirschbein, M. S. (1982). Bird impact analysis package for turbine engine fan blades. AIAA Paper. Hou, J. P., & Ruiz, C. (2007). Soft body impact on laminated composite materials. Composites Part A: Applied Science and Manufacturing, 38(2), 505–515. HOURGLASS (2010). LS-DYNA user’s manual. Ivancˇevic´, D., & Smojver, I. (2011). Hybrid approach in bird strike damage prediction on aeronautical composite structures. Composite Structures, 94(1), 15–23. Jenq, S. T., Hsiao, F. B., Lin, I. C., Zimcik, D. G., & Ensan, M. N. (2007). Simulation of a rigid plate hit by a cylindrical hemi-spherical tip-ended soft impactor. Computational Materials Science, 39(3), 518–526.

158

Bird Strike

Johnson, A. F., & Holzapfel, M. (2003). Modelling soft body impact on composite structures. Composite Structures, 61(1), 103–113. Lacome, J. (2000). Smooth particle hydrodynamics (SPH): A new feature in LSDYNA. In: Proceedings of the 6th international LS-DYNA users conference, session 7–3. Lacome, J. L. (2001). Smoothed particle hydrodynamics – Part II. In FEA information international news for the world-wide engineering community (pp. 6–11). Langrand, B., Bayart, A. S., Chauveau, Y., & Deletombe, E. (2002). Assessment of multiphysics FE methods for bird strike modelling – Application to a metallic riveted airframe. International Journal of Crashworthiness, 7(4), 415–428. Lavoie, M. A., Gakwaya, A., Ensan, M. N., & Zimcik, D. G. (2007). Validation of available approaches for numerical bird strike modeling tools. International Review of Mechanical Engineering, 1(4), 380–389. Lavoie, M. A., Gakwaya, A., Ensan, M. N., Zimcik, D. G., & Nandlall, D. (2009). Bird’s substitute tests results and evaluation of available numerical methods. International Journal of Impact Engineering, 36(10), 1276–1287. Liu, J., Li, Y., & Gao, X. (2014). Bird strike on a flat plate: Experiments and numerical simulations. International Journal of Impact Engineering, 70, 21–37. LS-DYNA Keyword User’s Manual. v.9.71. (2006). Livermore software. Mao, R. H., Meguid, S. A., & Ng, T. Y. (2008). Transient three dimensional finite element analysis of a bird striking a fan blade. International Journal of Mechanics and Materials in Design, 4(1), 79–96. Mao, R. H., Meguid, S. A., & Ng, T. Y. (2009). Effects of incidence angle in bird strike on integrity of aero-engine fan blade. International Journal of Crashworthiness, 14(4), 295–308. Marulo, F., & Guida, M. (2014). Design criteria for birdstrike damage on windshield. Advances in Aircraft and Spacecraft Science, 1(2), 233–251. McCallum, S. C., & Constantinou, C. (2005). The influence of bird-shape in bird-strike analysis. In: 5th European LS-DYNA users conference, Birmingham, UK. McCarthy, M. A., Xiao, J. R., McCarthy, C. T., Kamoulakos, A., Ramos, J., Gallard, J. P., & Melito, V. (2004). Modeling of bird strike on an aircraft wing leading edge made from fiber metal laminates – Part 2: Modeling of impact with SPH bird model. Applied Composite Materials, 11, 317–340. McCarthy, M. A., Xiao, J. R., McCarthy, C. T., Kamoulakos, A., Ramos, J., Gallard, J. P., & Melito, V. (2005). Modelling bird impacts on an aircraft wing – Part 2: Modelling the impact with an SPH bird model. International Journal of Crashworthiness, 10(1), 51–59. McCarty, J. E. (1982). 737 graphite-epoxy horizontal stabilizer certification. AIAA Paper. McCarty, R. E. (1980a). Computer analysis of bird-resistant aircraft transparencies. In: Survival and Flight Equipment Association, 17th annual symposium (pp. 93–97) Las Vegas, Nev., USA. McCarty, R. E. (1980b). Finite element analysis of F-16 aircraft canopy dynamic response to bird impact loading. In Structures, structural dynamics, and materials conference, Seattle, Wash., USA. McCarty, R. E. (1983). MAGNA computer simulation of bird impact on the TF-15 aircraft canopy. In Proceedings of the 14th conference on aerospace transparent materials and enclosures (pp. 974–1008). Scottsdale, AZ. Meguid, S. A., Mao, R. H., & Ng, T. Y. (2008). FE analysis of geometry effects of an artificial bird striking an aero engine fan blade. International Journal of Impact Engineering, 35(6), 487–498. Meyers, M. A. (1994). Dynamic behavior of materials. New York, USA: Wiley.

Finite element bird-strike modeling

159

Mie, G. (1903). Zur kinetischen Theorie der einatomigen K€ orper. Annalen der Physik, 316(8), 657–697. Monaghan, J. (1989). On the problem of penetration in particle methods. Journal of Computational Physics, 82(1), 1–15. Monaghan, J. J. (1992). Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics, 30, 543–574. Nishikawa, M., Hemmi, K., & Takeda, N. (2011). Finite-element simulation for modeling composite plates subjected to soft-body, high-velocity impact for application to bird-strike problem of composite fan blades. Composite Structures, 93(5), 1416–1423. Nizampatnam, L. S. (2007). Models and methods for bird strike load predictions. PhD Thesis. Oasys LS-DYNA environment 8.1: User Guide (2001). Blythe Valley Park, West Midlands: Oasys Ltd Petrinic, N., & Duffin, R. (2000). Discrete element modeling of soft body impact against rigid targets. 3rd B2000 users workshop, Enschede, The Netherlands. Richard, B. (2000). The development of a substitute artificial bird by the international bird strike research group for use in aircraft component testing. In ISBC25/ WP-IE3. Amsterdam: International Bird Strike Committee. Salehi, H., Ziaei-Rad, S., & Vaziri-Zanjani, M. A. (2010). Bird impact effects on different types of aircraft bubble windows using numerical and experimental methods. International Journal of Crashworthiness, 15(1), 93–106. Shmotin, Y. N., Chupin, P. V., Gabov, D. V., Ryabov, A. A., Romanov, V. I., Kukanov, S. S., & Saturn, N. (2009). Bird strike analysis of aircraft engine fan. In: Proceedings of the 7th European LS-DYNA users conference, Salzburg, Austria. Smojver, I., & Ivancˇevic´, D. (2010). Numerical simulation of bird strike damage prediction in airplane flap structure. Composite Structures, 92(9), 2016–2026. Smojver, I., & Ivancˇevic´, D. (2012). Advanced modelling of bird strike on high lift devices using hybrid Eulerian–Lagrangian formulation. Aerospace Science and Technology, 23(1), 224–232. Stoll, F., & Brockman, R. A. (1997). Finite element simulation of high speed soft-body impacts. In Proceedings of the 1997 38th AIAA/ASME/ASCE/AHS/ASC structure, structural dynamics, and materials conference (pp. 334–344). Kissimmee, FL, USA. Storace, A. F., Nimmer, R. P., & Ravenhall, R. (1984). Analytical and experimental investigation of bird impact on fan and compressor blading. Journal of Aircraft, 21(7), 520–527. Sun, C. T. (1972). An analytical method for evaluation of impact damage energy of laminated composites. In Composite Materials: Testing and Design (Fourth Conference), ASTM STP (Vol. 617, pp. 427–440). Tho, C. H., & Smith, M. R. (2011). Accurate bird strike simulation methodology for BA609 tiltrotor. Journal of the American Helicopter Society, 56(1), 12007-1–12007-10. Ugrcˇic´, M. (2012). Application of the hydrodynamic theory and the finite element method in the analysis of bird strike in a flat barrier. Scientific Technical Review, 62(3–4), 28–37. Vignjevic, R., Orłowski, M., De Vuyst, T., & Campbell, J. C. (2013). A parametric study of bird strike on engine blades. International Journal of Impact Engineering, 60, 44–57. Wang, F. S., & Yue, Z. F. (2010). Numerical simulation of damage and failure in aircraft windshield structure against bird strike. Materials and Design, 31(2), 687–695. Wilbeck, J. (1978). Impact behavior of low strength projectiles. Air Force Wright Aeronautical Labs, Air Force Materials Lab. OH: Wright-Patterson Air Force Base. Yulong, L., Yongkang, Z., & Pu, X. (2008). Study of similarity law for bird impact on structure. Chinese Journal of Aeronautics, 21(6), 512–517. Zukas, J. (2004). Introduction to hydrocodes. Elsevier.