RD-E: 4900 Bird Strike on Windshield

Using SPH elements to simulate a bird hitting a windshield.

This example will introduce how to simulate a bird hitting a windshield.

Options and Keywords Used

/SPHGLO
/SPHCEL
/MAT/LAW6 (HYDRO or HYD_VISC)
/MAT/LAW36 (PLAS_TAB)
/MAT/LAW27 (PLAS_BRIT)
/PROP/TYPE34 (SPH)
/PROP/TYPE11 (SH_SANDW)
/INIVEL/TRA
Boundary condition (/BCS)
Rigid body (/RBODY)
Interface (/INTER/TYPE7 and /INTER/TYPE12)

Input Files

The input files used in this example include:

Bird Strike on Windshield: <install_directory>/hwsolvers/demos/radioss/example/49_bird_strike_on_windshield/*

Model Description

A bird weighing about 1.93kg (4 pounds) with an initial speed of 80m/s impacts a windshield of a plane. The bird model is extremely simplified using a cylinder (diameter: 106mm; height: 195mm), it represents all major options used in a bird strike analysis. The edge of the plane model is fixed using a rigid body. Boundary conditions are applied on the main node.

Units: mm, ms, kg, KN, GPa

The bird model has following characteristics.

For the bird /MAT/LAW6 is used which describes a fluid material and is based on water properties /MAT/HYDRO. Pressure is computed using Equation of State provided by definition of /EOS/POLYNOMIAL option.

Material Properties
Initial density: 9.82242e-7 $[\frac{k g}{m m^{3}}]$
C1 Hydrodynamic pressure law coefficient: 2.24
/PROP/TYPE34 (SPH) (face centered cubic): mp: 1.725149e-4 $[kg]$
qa: 2e-30
qb: 1e-30
order: 0
h: 6.286 $[mm]$

Note:

In property, mp is the mass of each particle. The bird model consists of 11190 particles, so the bird weighs 11190*1.725149e-4 = 1.93kg (4 pounds).
For no damping qa was set 2e-30 and qb to 1e-30.
To simulate a bird strike, the compression phase of the bird will be also considered; therefore, the initial compression wave should properly propagate through the bird. Energy conservation is also important, so that the bird communicates the most realistic amount of energy to the target, thus ensuring that the impact is reasonably conservative.
In property, $h$ is smoothing length:(1) $h = {(\frac{m_{p} \sqrt{2}}{ρ})}^{\frac{1}{3}} = {(\frac{1.725149 \cdot 10^{- 4} \sqrt{2}}{9.82242 \cdot 10^{- 7}})}^{\frac{1}{3}} = 6.286 [m m]$
The smoothing length, $h$ is highly recommended to use the minimal distance between each particle and its closest neighbor.
SPH correction order is used to satisfy the consistency conditions and to ensure a better representation when the particles are not so well organized than into a hexagonal compact net. High order value is very expensive. In this example, order is set to 0.

For the windshield in this example, composite material with matrix glass and plastic fabric is used. For both fabrics /MAT/LAW27 (PLAS_BRIT) (composite with isotropic layers) is used.

Glass Matrix

Material Properties
Initial density: 1.74e-6 $[\frac{k g}{m m^{3}}]$
E-Modules: 61 [GPa]
Poisson ratio: 0.3
Yield stress, a: 0.045 [GPa]
Plasticity hardening, b: 0.05
$ε_{t 1}$ (Tensile failure strain in principle strain dir. 1): 0.01
$ε_{m 1}$ (Maximum tensile failure strain in principal strain dir. 1): 0.03
element deleted at $ε_{f 1}$ (Tensile strain for element deletion in principal strain dir. 1): 0.06
Figure 3.

Fabric Plastic

Material Properties
Initial density: 8e-7 $[\frac{k g}{m m^{3}}]$
E-Modules: 3 [GPa]
Poisson ratio: 0.3
Yield stress, a: 1e30 [GPa]

/PROP/TYPE11 (SH_SANDW):

Material Properties
Thickness of windshield: 9 [mm]
Number of layers, N: 9
I_thick: 1; Figure 4. Windshield consisting of composite material with matrix glass and plastic fabric

For metal parts (aluminum) in the airplane, /MAT/LAW36 (PLAS_TAB) is used.

Material Properties
Initial density: 2.7E-6 $[\frac{k g}{m m^{3}}]$
E-Modules: 73.1 [GPa]
Poisson ratio: 0.33

For metal parts connection (bolt) in the airplane, /MAT/LAW2 (PLAS_JOHNS) and /PROP/TYPE 3 (BEAM) are used.

/MAT/LAW2 (PLAS_JOHNS)
Initial density: 7.8E-6 [ $[\frac{k g}{m m^{3}}]$
E-Modules: 210 [GPa]
Poisson ratio: 0.28
Yield stress, a: 0.58 [GPa]
Hardening parameter, b: 0.6
Hardening exponent, n: 1
Strain rate coefficient: 0.017

/PROP/TYPE 3 (BEAM)
Area: 50.265 [mm²]
Iyy = Izz (YY and ZZ section inertia, bending): 201.061 [mm⁴]
Ixx (XX section inertia, torsion): 403.12 [mm⁴]

Model Method

Bird modeling
The bird model consists of 11190 particles, the distribution follows a hexagonal compact net. The cell distribution must be cubic face centered. The diameter of each particle is about 5.0 mm.

Figure 5. Perspective View of the Cubic Centered Faces Net
SPH simulation setup
In order to define the SPH bird, the following keywords are needed:
- /PROP/TYPE34 (SPH) to define the bird SPH cells property.
- /SPHCEL to define SPH cells for bird. Use bird part ID in this keyword.
- /SPHGLO to define the SPH parameters for the simulation.
  The default value of 0.25 for $α_{s o r t}$ is used. This coefficient controls the search for neighbors, so that for each particle more than the actual neighbors are found. This reduces the computational time.
Contact between SPH bird and structure and impact setup
- Interface TYPE7 is selected to simulate the contact between SPH bird and windshield structure. The windshield structure is selected as main surface and the SPH bird as secondary nodes.
- I_gap=1 (use variable gap + gap scale correction + considering mesh size). The gap is ½ particle diameter + ½ thickness of the target. This is the physical value for contact.
- The stiffness between bird and windshield structure are quite different. In order to get results closer to reality, set I_stf=0 and Stfac=0.1. This means interface stiffness equals to 0.1 times stiffness of main side (structure).
- /INIVEL/TRA is used to set bird part with initial velocity 80m/s strike on the structure.
Connection between the windshield and fuselage
- For the weld connection between windshield and fuselage /PROP/TYPE13 spring beam in combination with interface TYPE2 tied contact are used.
Other remarks
The animation out uses the H3D format via /H3D. If /ANIM output is used, the /ANIM/VERS/44 option is needed to get the best display of SPH particles.

Results

Figure 6 shows von Mises stress of the windshield at time 3ms ~ 8ms.

The Interface Force between the bird and the windshield is shown in Figure 7.

Conclusion

The mesh size and the number of particles of the bird model seem to be a good compromise between the accuracy of the problem and the CPU calculation time.
For interface parameters, default values could be used. In TYPE7 contact interface, the interface stiffness I_stf has to be set to zero.
In the model a simplified hydro viscous law is used, which is well suited to describe (represent) the most important properties of the bird model.
The smoothing length h should be set to a value larger than the inter-cell distance (to avoid excessive de-cohesion), but not too large to avoid excessive energy absorption by the particles of the bird.
The differences between the real tests and the simulations are due to the used material models and their parameters, boundary conditions, target material laws and many other sources.