Dynamic Responses of the Aero-Engine Rotor System to Bird Strike on Fan Blades at Different Rotational Speeds

: To study the effect of a bird striking engine fan on the rotor system, a low-pressure rotor system dynamic model based on a real aero-engine structure was established. Dynamic equations were derived considering the case of the bird strike force which transferred to the rotor system. The bird strike force was obtained from the bird strike process simulation in LS-DYNA, where a smoothed particle hydrodynamics (SPH) mallard model was constructed using a computed tomography (CT) scanner, and ﬁnite element method (FEM) was used to simulate the bird strike on an actual fan model. The dynamic equations were solved using the Newmark- β method. The effect of rotational speeds on the rotor system dynamics after bird strike was investigated and discussed. Results show that the maximum bird impact force can reach 104 kN at 3772 r/min. Impact time is only 0.06 s, but the bird strike on fan blades lead to a transient shock on the rotor system. Under the action of transient shocks, the rotor system displacement in the horizontal and vertical directions increase sharply, and the closer the mass point is to the fan, the more it is affected; the vibration amplitude at the fan will increase 15 times within 0.1 s of the bird strike and will gradually decrease with the effect of damping. The dynamics of the rotor system changes from a stable single periodic motion to a complex irregular quasi-periodic motion after a bird strike, and the strike force excites the ﬁrst-order vibrational mode of the rotor system. This phenomenon occurs at all speeds when bird strikes occur. Bird strikes will cause resonance in the rotor system, which may cause damage to the engine. It was also seen that the bird strike force, and hence the effects on the rotor system, increases as the engine rotational speed increases; the peak force is larger and the number of peaks has increased. The impact force at 3772 r/min is 99.5 kN higher than at 836 r/min, and three additional peaks emerged. This effect is more reﬂected in the amplitude, and the overall vibration characteristics do not change. Combining the bird strike with the rotor dynamics calculation, the dynamic response of the aero-engine rotor system to bird strike is studied at different ﬂight stages, which is of guiding signiﬁcance for power evaluation of aero engines after bird strike.


Introduction
The bird strike is a common hazard of aircraft which has a significant effect on their operation safety. According to the statistics of the Federal Aviation Administration (FAA), wildlife strikes have damaged more than 271 aircraft and have led to more than 292 deaths worldwide between 1988 and 2019 [1]. As shown in Figure 1a, the wildlife strike frequency on aircraft has been constantly increasing since 1990. In 2019, 17,760 wildlife strikes took strikes took place, among which the proportion of bird strike was 94.4%. As shown in Figure 1b, the proportion of engine damage is the highest, reaching 26%, namely, the aero engine is a frequent occurrence area for bird striking [2]. When a bird impacts a jet engine, the fan is struck first, which may cause the failure of fan blades and induce secondary faults. Research on the effect of bird strike on the operation of an aero engine is of great importance to the working safety and reliability of the aero engine. In early bird strike studies, experimental methods were the most valid and popular means of gaining knowledge about bird strike. However, the experimental tests are expensive and difficult to perform [3]. Therefore, numerical simulation methods have been widely used and fully developed in the past 30 years [4]. Several simple geometries have been adopted to represent the shape of a bird. The geometry of the substitute bird models are of four main types: straight-ended cylinder, hemispherical-ended cylinder, ellipsoid, and sphere [5,6]. The hemispherical-ended cylinder has been the most commonly used geometry due to its more accurate results [7]. More recently, bird models with realistic geometries based on CT scan images have replaced the geometric model due to the more accurate results [8,9].
As for the target fan model, the initial approach was to use a plate with a thickness [19,20]. In more recent studies, fan blades with realistic geometries have become more common, whereafter transient response of rotating is obtained under bird-slicing impact [21][22][23]. Meanwhile, different cases of bird strikes with different bird mass and velocity, impact location and timing, and rotary speed have been investigated on the rotary engine [24,25]. In this paper, the SPH bird impactor model has been chosen and the finite element method on a fan model with realistic geometry is applied.
Bird strikes often have a dramatic effect on the aircraft's engines. The impact dynamics analysis under bird strike (which can be considered as foreign object impacts) can cause severe deformation in the fan blades, leading to impact-induced force generation [21]. Bird strike experiments on the engine primary compressor under different bird strike conditions were investigated by Liu et al., and a significant effect on the engine performance was found [24]. A bird strike creates a substantial impact force, which may apply significant imbalanced loads on the fan shaft, leading to problems occurring in the rotor performance [26]. In early bird strike studies, experimental methods were the most valid and popular means of gaining knowledge about bird strike. However, the experimental tests are expensive and difficult to perform [3]. Therefore, numerical simulation methods have been widely used and fully developed in the past 30 years [4]. Several simple geometries have been adopted to represent the shape of a bird. The geometry of the substitute bird models are of four main types: straight-ended cylinder, hemispherical-ended cylinder, ellipsoid, and sphere [5,6]. The hemispherical-ended cylinder has been the most commonly used geometry due to its more accurate results [7]. More recently, bird models with realistic geometries based on CT scan images have replaced the geometric model due to the more accurate results [8,9].
As for the target fan model, the initial approach was to use a plate with a thickness [19,20]. In more recent studies, fan blades with realistic geometries have become more common, whereafter transient response of rotating is obtained under bird-slicing impact [21][22][23]. Meanwhile, different cases of bird strikes with different bird mass and velocity, impact location and timing, and rotary speed have been investigated on the rotary engine [24,25]. In this paper, the SPH bird impactor model has been chosen and the finite element method on a fan model with realistic geometry is applied.
Bird strikes often have a dramatic effect on the aircraft's engines. The impact dynamics analysis under bird strike (which can be considered as foreign object impacts) can cause severe deformation in the fan blades, leading to impact-induced force generation [21]. Bird strike experiments on the engine primary compressor under different bird strike conditions were investigated by Liu et al., and a significant effect on the engine performance was found [24]. A bird strike creates a substantial impact force, which may apply significant imbalanced loads on the fan shaft, leading to problems occurring in the rotor performance [26].
Transient impacts on the rotor system can affect the stability of the whole system. Adding shock excitation to the bearing location of the rotor system can result in a significant impact on the rotor system [27]. The transient impact caused by airflow excitation has an obvious influence on the dynamic response of the rotor system [28]. Meanwhile, the engine rotates at high speed and its response is subjected to sudden unbalanced forces that have been studied [22]. Bird strike will not only apply transient impact load, but also induce rubbing between fan and casing [29]; then, a more common phenomenon has emerged. Transient shocks caused by rub-impact have been widely studied, and most studies have investigated the transient dynamic response and nonlinear characteristics of the system by introducing a rub-impact system model into the rotor system [30][31][32], the rub-impact model from the simple rotating static sub rubbing [33], development for the blade-casing rubbing [34][35][36], and the now widely used blade-coating rubbing [37]. Different models produce different transient shocks and have different magnitudes of impact on the rotor, but all studies could illustrate the high susceptibility of rotor systems to transient shocks, leading to changes in operating conditions. In most studies [33,34,[37][38][39], the impact load is often added to the dynamic equation as a function. The response of rotor dynamic is complicated when impact loads are added to the equation. However, as impact loads are often irregular, and hard to be described by the equations for some cases, this may lead to the difference between numerical results and actual responses of the rotor system. In this paper, to overcome the irregularity of bird impact force and the difficulty of fitting the equation, the impact force obtained by simulation is added to the equation by interpolation.
From the above-mentioned studies, it is evident that, even though in many previous studies transient shocks have been introduced into the rotor systems, only a limited number of models have directly incorporated bird impact forces into the dynamics calculations. This is despite the fact that the bird impact forces are also a distinct transient force that can have severe effects on rotor systems. Therefore, in this paper, based on the simulation of the bird strike engine fan process, the scope of the research is limited to the deformation of the fan blades due to bird strike, and the case of the blade out is excluded. The impact force generated by the bird strike on the system is obtained, and the dynamic response of the rotor system before and after the occurrence of the bird strike is analyzed, which provides some reference for the study of the rotor impact mechanics of the bird strike engine.

Bird Strike Modeling
According to the statistical report from FAA [1], during 1990-2019, mallards struck aircraft a total of 1146 times, making it one of the birds with the highest frequency of collisions with aero-engines. Therefore, the mallard was taken as the prototype bird model. The weight of the selected mallard was 820 g with 40 cm in body length and a 45 cm wingspan. The mallard bird was scanned by the CT scan device at a medical imaging center. A total of 1566 DICOM (Digital Imaging and Communications in Medicine) images were obtained. DICOM images and the SPH method were adopted to model the mallard. During the modeling process, the internal cavities of the mallard bird were considered, and no SPH elements were used in such areas. The final model composed of 41,685 SPH elements each weighing 0.0191 g ( Figure 2).
The bird was endowed with a Null material model with Gruneisen equation of state (EOS) [5,40]; the viscous stress was calculated by providing the intrinsic relationship of the model. The state equation was used to calculate the pressure, with the material initial density of 938 kg/m 3 [40]. Based on the EOS [41], the pressure of the compressed material is defined as follows:  The bird was endowed with a Null material model with Gruneisen equation of state (EOS) [5,40]; the viscous stress was calculated by providing the intrinsic relationship of the model. The state equation was used to calculate the pressure, with the material initial density of 938 kg/m 3 [40]. Based on the EOS [41], the pressure of the compressed material is defined as follows: For expansion of the material, the pressure is expressed as: Where C is the intercept of the velocity curve; S1, S2, and S3 are the slope coefficients of the velocity curve; is the Gruneisen constant; a is the first-order volume correction to ; E is the Elastic Modulus of material; and is the reference density. The parameter is defined by: where is density of the material. Based on Ref [42], the EOS's parameters of C = 1480, S1 = 1.92, S2 = 0, S3 = 0, and γ = 0.1 were used.
A type of aero turbofan engine fan system is demonstrated in Figure 3. The fan system consisted of 24 narrow chord blades. The height of each blade was 603.2 mm, and the initial torsion angle was 61.3°. The blade was connected to the hub with a dovetail and tongue-and-groove mating method. The entire fan system was discretized using 1,520,208 hexahedral solid elements (3 mm in size), see Figure 3. The material properties of Ti-6Al-4 V alloy with plastic dynamics model (MAT_PLASTIC_KINEMATIC) was employed for the fan blade (Table 1).  For expansion of the material, the pressure is expressed as: where C is the intercept of the velocity curve; S 1 , S 2 , and S 3 are the slope coefficients of the velocity curve; γ 0 is the Gruneisen constant; a is the first-order volume correction to γ 0 ; E is the Elastic Modulus of material; and ρ 0 is the reference density. The parameter µ is defined by: where ρ is density of the material. Based on Ref [42], the EOS's parameters of C = 1480, S 1 = 1.92, S 2 = 0, S 3 = 0, and γ 0 = 0.1 were used. A type of aero turbofan engine fan system is demonstrated in Figure 3. The fan system consisted of 24 narrow chord blades. The height of each blade was 603.2 mm, and the initial torsion angle was 61.3 • . The blade was connected to the hub with a dovetail and tongue-and-groove mating method. The entire fan system was discretized using 1,520,208 hexahedral solid elements (3 mm in size), see Figure 3. The material properties of Ti-6Al-4 V alloy with plastic dynamics model (MAT_PLASTIC_KINEMATIC) was employed for the fan blade (Table 1). The bird was endowed with a Null material model with Gruneisen equation of state (EOS) [5,40]; the viscous stress was calculated by providing the intrinsic relationship of the model. The state equation was used to calculate the pressure, with the material initial density of 938 kg/m 3 [40]. Based on the EOS [41], the pressure of the compressed material is defined as follows: For expansion of the material, the pressure is expressed as: Where C is the intercept of the velocity curve; S1, S2, and S3 are the slope coefficients of the velocity curve; is the Gruneisen constant; a is the first-order volume correction to ; E is the Elastic Modulus of material; and is the reference density. The parameter is defined by: where is density of the material. Based on Ref [42], the EOS's parameters of C = 1480, S1 = 1.92, S2 = 0, S3 = 0, and γ = 0.1 were used.
A type of aero turbofan engine fan system is demonstrated in Figure 3. The fan system consisted of 24 narrow chord blades. The height of each blade was 603.2 mm, and the initial torsion angle was 61.3°. The blade was connected to the hub with a dovetail and tongue-and-groove mating method. The entire fan system was discretized using 1,520,208 hexahedral solid elements (3 mm in size), see Figure 3. The material properties of Ti-6Al-4 V alloy with plastic dynamics model (MAT_PLASTIC_KINEMATIC) was employed for the fan blade (Table 1).   When a bird hits the fan blade, there is an uncertainty in the impact location, and different impact locations led to different damage effects [43]. To better reflect the effect of impact forces on the rotor system, a single impact location was selected for the simulations. The same applies to the selection of the impact orientations. Five impact positions: 1/6, 1/3, 1/2, 2/3, and 5/6 fan blade height, and 15 types of impact orientations were considered. Based on the previous studies [44], the bird impact force is pronounced at 1/2 blade height, which leads to significant impact on the rotor system. Therefore, 1/2 fan blade height as the impact position and the head side as the orientation were chosen, which is given in Figure 4.

Rotor System Modeling
The operation state of an aero engine may be affected by the impact force caused by a bird strike. Understanding the response characteristics of the rotor system under bird strike is useful for grasping the aero-engine operation state. The key to studying the dynamic response of an aero-engine rotor system is to develop a rotor dynamics model that can describe the effects of bird strikes on the fan blade. The schematic diagram of the aeroengine low-pressure rotor system is depicted in Figure 6. The low-pressure rotor system is supported by three rolling bearings (#1, #2, and #3), including a 1-stage fan and a 4-stage low-pressure turbine, and the bearing supporting mode is "0-2-1". In the simulation, the sucked-in velocity of the bird was 116 m/s. Some of the previous studies [5] on the bird impacting fan did not take the rotation of the fan into account, namely, the blades were stationary. As a result, two problems that made the simulation deviate the actual impact process arose: (1) the dynamic stiffening effect on the blade due to the centrifugal force was ignored; (2) the direction of the contact force was constant, as the relative velocity between the bird and the blade was not considered. For a more realistic simulation of the bird strike onto the fan, the rotor system was run at different rotational speeds. Bird strike occurs in several stages of flight [1]. The differences of fan speed at each stage leads to the different dynamic responses of the fan blade and thus the rotor system. According to the fan rotation speed in QAR (quick access recorder) data of an aero engine, four rotational speeds were 836 r/min, 3772 r/min, 3344 r/min, and 1984 r/min corresponding to five typical operating conditions of the ground idle stage, approaching stage, climbing stage, cruising stage, and descending phase. The noted speeds were selected as the initial fan rotational velocities. The bird model was verified by comparing the numerical results with Wilbeck's experimental results [41,45].
Comparing the experimental data of Wilbeck, the simulation data of the traditional model (sphere model), and the mallard model built in this paper, we monitored the pressure magnitude of the model after striking on the target. From Figure 5, when using the mallard model, the magnitude of the pressure values and the pattern of its variation with time are close to those of the Wilbeck's experiment compared to the conventional model. The pressure difference does not exceed 15%, other than the peak value, so the bird collision model can be considered to reflect the real bird collision process to some extent.

Rotor System Modeling
The operation state of an aero engine may be affected by the impact force caused by a bird strike. Understanding the response characteristics of the rotor system under bird strike is useful for grasping the aero-engine operation state. The key to studying the dynamic response of an aero-engine rotor system is to develop a rotor dynamics model that can describe the effects of bird strikes on the fan blade. The schematic diagram of the aeroengine low-pressure rotor system is depicted in Figure 6. The low-pressure rotor system is supported by three rolling bearings (#1, #2, and #3), including a 1-stage fan and a 4-stage low-pressure turbine, and the bearing supporting mode is "0-2-1".

Rotor System Modeling
The operation state of an aero engine may be affected by the impact force caused by a bird strike. Understanding the response characteristics of the rotor system under bird strike is useful for grasping the aero-engine operation state. The key to studying the dynamic response of an aero-engine rotor system is to develop a rotor dynamics model that can describe the effects of bird strikes on the fan blade. The schematic diagram of the aero-engine low-pressure rotor system is depicted in Figure 6. The low-pressure rotor system is supported by three rolling bearings (#1, #2, and #3), including a 1-stage fan and a 4-stage low-pressure turbine, and the bearing supporting mode is "0-2-1".
Appl. Sci. 2021, 11, x FOR PEER REVIEW 7 of 22 Figure 6. Schematic diagram of aero-engine low-pressure rotor structure and its support system.
To establish an accurate rotor dynamics model, three coordinate scanning and 3D modeling of the rotor shaft, disk, and other components were carried out: first, BEST-SCAN 751 (HOLON 3D, China) 3D scanner and HL-Scan software were used to obtain the point cloud data of the components, and the data was then imported into the software CATIA for constructing 3D models. Combined with the material parameters of each component, the mass distribution and center of mass position of the rotor system were obtained. Based on the obtained data points, the fan, the shaft, and the disk were simplified as several mass elements through the lumped-mass method, and a low-pressure rotor system model of the aero engine was established, as shown in Figure 7.
When using the lumped-mass method, the blade and the disk are treated as a mass point, so the bird impact force is considered to the disk, and its actual impact on the blade is similar; the mass of the disk is much larger than the mass of the mass point on the shaft, and the influence of its eccentricity is much larger than the eccentricity of the shaft, so the influence of concentrating the eccentricity on the disk on the calculation results is accepted. As there is no need to consider the connection method of each component in the calculation, considering it as a whole will have no effect on the calculation results; as only the x and y directional displacement of the rotor is considered in this paper, the torque is ignored in the simplification; and as the maximum speed is 3772 r/min and the speed is To establish an accurate rotor dynamics model, three coordinate scanning and 3D modeling of the rotor shaft, disk, and other components were carried out: first, BESTSCAN 751 (HOLON 3D, China) 3D scanner and HL-Scan software were used to obtain the point cloud data of the components, and the data was then imported into the software CATIA for constructing 3D models. Combined with the material parameters of each component, the mass distribution and center of mass position of the rotor system were obtained. Based on the obtained data points, the fan, the shaft, and the disk were simplified as several mass elements through the lumped-mass method, and a low-pressure rotor system model of the aero engine was established, as shown in Figure 7. Based on Newton's second law, the equations of motion for the established rotor system were derived as: where MR1 is the mass of Disk 1 of the rotor, MR2 is the mass of Disk 2, and MB1, MB2, MB3 are the masses of the concentrated points at the three rolling bearings, respectively. CR and CB are the damping coefficients of the disk and ball bearings, Kr is the stiffness of the shaft, supposing that the stiffness of each section of the rotor is identical. eR1 and eR2 are the eccentricities of Disk 1 and Disk 2, respectively. Bx and By are the decomposition of bird impact forces into the x and y directions, respectively. FBix and FBiy (i = 1, 2, 3) are the bearing support forces in the x and y directions. is the rotational speed of the rotor in rad/s. As this paper mainly studies the impact of bird strike on the fan rotor system, the bearings are regarded as an elastic system. The bearing stiffness in the rotor system was considered to be constant (Kb). The equations relating FBix and FBiy (i = 1, 2, 3), stiffness, and displacements could be given by: From the above equations, the bearing support force is proportional to the displacement, and the negative sign indicates that the reaction force points in the opposite direction of the displacement.
In our previous study [31], we used the same 3D scan method to build a '0-2-1' rotor test setup using the kinetic similarity principle (see Figure 8). The vibration response of the rotor system at 1500 rpm was measured using a test rig and is shown in Figure 9a. At a speed of 1500 r/min, the vibration waveform of the system shows obvious peaks at the peaks and valleys. The fundamental frequency amplitude is the main frequency of the When using the lumped-mass method, the blade and the disk are treated as a mass point, so the bird impact force is considered to the disk, and its actual impact on the blade is similar; the mass of the disk is much larger than the mass of the mass point on the shaft, and the influence of its eccentricity is much larger than the eccentricity of the shaft, so the influence of concentrating the eccentricity on the disk on the calculation results is accepted. As there is no need to consider the connection method of each component in the calculation, considering it as a whole will have no effect on the calculation results; as only the x and y directional displacement of the rotor is considered in this paper, the torque is ignored in the simplification; and as the maximum speed is 3772 r/min and the speed is low, the influence of the gyroscopic force is not obvious, therefore ignoring these two items has little influence on the calculation results.
As the focus of this research was the effect of bird strike on the rotor system in the x-y plane of radial vibration, other vibrational effects were not taken into account, and the following assumptions were considered: 1.
The force from the bird strike was taken into account on the disk.

2.
The unbalanced mass of the rotor system concentrated on the two rotor disks, and the mass of the rotor system was discretized into the five mass points on the shaft.

3.
The materials of the blade and the disk behave according to Hooke's law.

4.
The blade, disk, and shaft were considered to be rigidly connected to each other, ignoring the connections between them. 5.
The torsional vibration and gyroscopic moment of the system were neglected.
Based on Newton's second law, the equations of motion for the established rotor system were derived as: .. .. ..
where M R1 is the mass of Disk 1 of the rotor, M R2 is the mass of Disk 2, and M B1 , M B2 , M B3 are the masses of the concentrated points at the three rolling bearings, respectively. C R and C B are the damping coefficients of the disk and ball bearings, K r is the stiffness of the shaft, supposing that the stiffness of each section of the rotor is identical. e R1 and e R2 are the eccentricities of Disk 1 and Disk 2, respectively. B x and B y are the decomposition of bird impact forces into the x and y directions, respectively. F Bix and F Biy (i = 1, 2, 3) are the bearing support forces in the x and y directions. ω is the rotational speed of the rotor in rad/s.
As this paper mainly studies the impact of bird strike on the fan rotor system, the bearings are regarded as an elastic system. The bearing stiffness in the rotor system was considered to be constant (K b ). The equations relating F Bix and F Biy (i = 1, 2, 3), stiffness, and displacements could be given by: From the above equations, the bearing support force is proportional to the displacement, and the negative sign indicates that the reaction force points in the opposite direction of the displacement.
In our previous study [31], we used the same 3D scan method to build a '0-2-1' rotor test setup using the kinetic similarity principle (see Figure 8). The vibration response of the rotor system at 1500 rpm was measured using a test rig and is shown in Figure 9a. At a speed of 1500 r/min, the vibration waveform of the system shows obvious peaks at the peaks and valleys. The fundamental frequency amplitude is the main frequency of the rotor system, and its amplitude is much larger than other high frequencies. The numerical method is used to calculate the rotor response at 1500 r/min, as shown in Figure 9b, and it can be found that some of the results are the same, differing in units due to different model simplification treatments, which verifies the accuracy of the proposed rotor system model to a certain extent.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 22 rotor system, and its amplitude is much larger than other high frequencies. The numerical method is used to calculate the rotor response at 1500 r/min, as shown in Figure 9b, and it can be found that some of the results are the same, differing in units due to different model simplification treatments, which verifies the accuracy of the proposed rotor system model to a certain extent. rotor system, and its amplitude is much larger than other high frequencies. The numerical method is used to calculate the rotor response at 1500 r/min, as shown in Figure 9b, and it can be found that some of the results are the same, differing in units due to different model simplification treatments, which verifies the accuracy of the proposed rotor system model to a certain extent.

Bird Strike Force at Different Rotational Speeds
The bird strike simulation was performed in LS-DYNA to obtain the impact forces transmitted to the rotor system. The bird state at different times during the impact process is shown in Figure 10, and it shows that the bird body SPH model has been cut into particles after the strike.

Bird Strike Force at Different Rotational Speeds
The bird strike simulation was performed in LS-DYNA to obtain the impact forces transmitted to the rotor system. The bird state at different times during the impact process is shown in Figure 10, and it shows that the bird body SPH model has been cut into particles after the strike. The bird strike process was simulated at four rotational speeds, and the contact forces in the x and y directions (in coordinate system depicted in Figure 7) are extracted and shown in Figure 11.
(a) (b) Figure 11. The bird impact forces at four rotational speeds: (a) in the X direction, and (b) in the Y direction. Figure 11 shows that the peaks of the magnitude of the bird impact forces in both the x and y directions increase with the rotational speed of the rotor. Moreover, the force in the x direction is significantly larger than that in the y direction, which is caused by the The bird strike process was simulated at four rotational speeds, and the contact forces in the x and y directions (in coordinate system depicted in Figure 7) are extracted and shown in Figure 11.

Bird Strike Force at Different Rotational Speeds
The bird strike simulation was performed in LS-DYNA to obtain the impact forces transmitted to the rotor system. The bird state at different times during the impact process is shown in Figure 10, and it shows that the bird body SPH model has been cut into particles after the strike. The bird strike process was simulated at four rotational speeds, and the contact forces in the x and y directions (in coordinate system depicted in Figure 7) are extracted and shown in Figure 11.
(a) (b) Figure 11. The bird impact forces at four rotational speeds: (a) in the X direction, and (b) in the Y direction. Figure 11 shows that the peaks of the magnitude of the bird impact forces in both the x and y directions increase with the rotational speed of the rotor. Moreover, the force in the x direction is significantly larger than that in the y direction, which is caused by the  Figure 11 shows that the peaks of the magnitude of the bird impact forces in both the x and y directions increase with the rotational speed of the rotor. Moreover, the force in the x direction is significantly larger than that in the y direction, which is caused by the geometry of the fan blade.
The curve at 3772 r/min includes several peaks, which is due to the splitting of the bird body into multiple pieces by the fan blades during the impact process. The magnitude of each peak is highly dependent on the size of the colliding bird body. The number of peaks varies at different speeds, which is due to the fact that the number of blades cutting the bird body varies with the fan rotational speed. The faster the fan rotates, the more blocks the bird body is divided into and, as a result, the more peaks appear. Additionally, is it notable that the direction of the force changes in some curves, especially in the y direction at the speed of 1984 r/min. This is as the strike process between the bird body and the blade can be divided into two parts. As the fan is rotating, the bird hits the blade, which results in the blade slapping the bird body. As can be seen in Figure 10, after the bird was cut by the fan, the blades may hit the bird body, thus causing a change in the direction of the force in Figure 8.
After obtaining the curve of variation of Bx and By with time, the data were substituted into the dynamic equation by linear interpolation method to calculate the rotor dynamic responses.

Numerical Simulation of Rotor System Dynamic
The second-order differential equations of motion of the rotor system were solved with the Newmark-β method. The required data are listed in Table 2. Table 2. The geometrical and material parameters of rotor system. In Table 2, µ B is the clearance of the ball bearing for the non-dimensional calculations. To facilitate the presentation of the results, some data has been dimensionless. The dimensionless displacement parameter is defined as X = x µ B and Y = y µ B , the dimensionless speed parameter is given by  A total of 400 calculation cycles were performed. Considering the volatility of the calculation, in this research, the bird impact force was added to the calculations after the rotor system results were stabilized. The bird impact force was added at t = 3 s to the calculations when the rotor system had sufficient time to be stabilized. The displacement, velocity and acceleration data of each mass point of the rotor system were obtained. The orbit plot, the Poincaré section map, the phase diagram, the time domain map, and the spectrogram were used to illustrate the dynamic response of the rotor system after bird strike.

The Dynamic Responses of Rotor System Fan before and after Bird Strike
The dynamic responses of the rotor system fan before and after bird strike are demonstrated in Figure 12. Figure 12a depicts the variations of the rotor displacement in the x direction vs. time in the steady-state operational condition when the bird strike has not yet occurred (t = 2-2.5 s). The periodic fluctuations of displacement in the x direction come from the interaction of eccentric force, stiffness and damping of the system. The system runs smoothly during this duration and the maximum dimensionless displacements remains constant at 0.8. Similarly, in Figure 12c, the same periodic motion can be observed, but in the y direction. With an initial displacement in the y direction due to the gravitational force, the amplitude range is −1.68 to −2.98. The system operated smoothly before the bird strike occurs. Figure 12e shows the axial rotation of the system in the steady-state operation, with the axis of the fan moving in the range of approximately −0.8 to 0.8 in the x direction and −3 to −1.8 in the y direction. The axial rotation remains constant as a closed and relatively regular circle.
(m) (n) Figure 12. The dynamic response of rotor system fan before and after bird strike at 3772 r/min: (a) time history of displacement before bird strike in the X direction, (b) time history of displacement after bird strike in the X direction, (c) time history of displacement before bird strike in the Y direction, (d) time history of displacement after bird strike in the Y direction, (e) axis orbit before bird strike, (f) axis orbit after bird strike, (g) phase diagram before bird strike, (h) phase diagram after bird strike, (i) the Poincaré section before bird strike, (j) the Poincaré section after bird strike, (k) frequency spectrum before bird strike, (l) frequency spectrum after bird strike, (m) time-frequency graph before bird strike, and (n) time-frequency graph after bird strike. Figure 12g shows the phase diagram of the rotor operation, where the trajectory line is demonstrated as a closed limit ring, indicating that the system has a single periodic motion. There is also only one isolated point in the Poincaré section map in Figure 12i. Both diagrams in Figure 12g,i indicate that the rotor system undergoes a stable single periodic motion before bird strike. It can be seen from Figure 12k that there is no other frequency component, except the fundamental frequency in the spectrum graph, which confirms the correctness of the previous graph. It would be intuitive to investigate the effect of the bird strike on the noted graphs. Figure 12b,d depict the horizontal and vertical time history of the fan displacement after the occurrence of bird strike at = 3 s. In the x direction, due to the influence of the Figure 12. The dynamic response of rotor system fan before and after bird strike at 3772 r/min: (a) time history of displacement before bird strike in the X direction, (b) time history of displacement after bird strike in the X direction, (c) time history of displacement before bird strike in the Y direction, (d) time history of displacement after bird strike in the Y direction, (e) axis orbit before bird strike, (f) axis orbit after bird strike, (g) phase diagram before bird strike, (h) phase diagram after bird strike, (i) the Poincaré section before bird strike, (j) the Poincaré section after bird strike, (k) frequency spectrum before bird strike, (l) frequency spectrum after bird strike, (m) time-frequency graph before bird strike, and (n) time-frequency graph after bird strike. Figure 12g shows the phase diagram of the rotor operation, where the trajectory line is demonstrated as a closed limit ring, indicating that the system has a single periodic motion. There is also only one isolated point in the Poincaré section map in Figure 12i. Both diagrams in Figure 12g,i indicate that the rotor system undergoes a stable single periodic motion before bird strike. It can be seen from Figure 12k that there is no other frequency component, except the fundamental frequency in the spectrum graph, which confirms the correctness of the previous graph. It would be intuitive to investigate the effect of the bird strike on the noted graphs. Figure 12b,d depict the horizontal and vertical time history of the fan displacement after the occurrence of bird strike at t = 3 s. In the x direction, due to the influence of the applied bird impact excitation, the amplitude of the rotor displacement changes dramatically after t = 3 s, with a maximum dimensionless amplitude close to 15 in the x direction. In the y direction, the maximum displacement amplitude increases for around 3.7 times compared to the initial position, which is due to the larger bird impact force in the x direction. The difference between the x direction and y direction displacements before and after the bird strike demonstrates that the magnitude of the bird impact force has a direct effect on the displacement of the fan rotor. Looking at the plots from t = 3 s to t = 3.5 s, it can be seen that the amplitude of the rotor displacement is modulated after the onset of the shock, showing the multi periodic motion and the response amplitude gradually decreases under the effect of damping, and that the motion is eventually is smoothed out and goes back to the situation before the impact. As the research did not consider the deformation of fan blades and the fan blade-out events during bird strike, the operation of the rotor system goes back to normal after bird strike.
From the axial orbit diagram in Figure 12f, it can be seen that, after the bird strike, the fan motion has deviated significantly from circular axial orbit curves in Figure 12e. It shows multiple chaotic intersecting curves with intense and irregular variations. The amplitude of motion also increases significantly compared to the stable operation before bird strike occurrence. However, there exists a clear stable region in the middle of the axial diagram ( Figure 12h). Meanwhile, it can be seen from the phase diagram (Figure 12h) that there are multiple limit rings of motion during this period, that the situation is chaotic, and that several overlapping situations occur. At the end of bird strike phenomenon, a stable limit ring appears in the middle of the phase diagram ( Figure 12h) and Poincaré section map (Figure 12j). There are several points in the Poincaré section map due to the bird strike, which indicate multiple cycles of the motion which gradually converge to a fixed position.
From above-mentioned analysis, it can be concluded that a bird strike on an aeroengine operating at 3772 r/min causes a strong transient impact on the rotor system, and that the rotor motion will be greatly disturbed, affecting the normal operation of the rotor. However, due to the instantaneity of the impact, the rotor system gradually regains its stability due to damping effect after the shock. It can be seen in Figure 12l that after the rotor is impacted by the bird, the rotor rotation frequency still exists, even though it has a slight increase in amplitude compared to stable operation. A sharp peak is produced at frequency of 43.2 Hz. This frequency is independent of the fundamental frequency, which is due to the bird strike exciting the first mode of the rotor. The time-frequency graph of the dynamic response before and after bird strike in the x direction is demonstrated in Figure 12m,n. During t = 2-2.5 s, the frequency band in the graph is concentrated on the position of the rotational frequency, as the bird strike had not occurred at this time. After the bird strike occurring at 3 s, a clear abrupt change in frequency emerges. In addition to the fundamental frequency of 62.9 Hz, a frequency of 43.2 Hz has appeared, which is approximately equal to the first-order mode. The first-order response of the rotor system (43.2 Hz) was excited by the transient shock caused by the bird strike.
In summary, the rotor motion becomes extremely unstable after a bird strike; therefore, corresponding measures must be taken, otherwise it will cause great harm to the normal operation of the aircraft engine.

The Dynamic Response of the Rotor System
The previous section analyzed the dynamic response of the fan position where the bird strike force acts directly, while the bird strike affects the whole rotor system. Figure 13 shows the displacement of other mass points vs. the time. Looking back at the time history of displacement diagrams of the fan in Figure 12b, the displacement amplitude of the bearings and the turbine is significantly smaller. When the fan is not hit by a bird, the dimensionless displacement amplitude is 0.87. After bird strike, the displacement rises to 14.1, an increase of 12.23. In the other mass points of the rotor, the displacements are, respectively 0.23, 0.01, 0.21, and 0.07 from the position close to the fan hub to the end of the blade before bird strike. After bird strike, the maximum displacements increase significantly to 3.97, 1.41, 1.13, and 0.38, respectively, for the noted positions. As the bird strike force directly acts on the fan hub, the closer the bird impact point to the fan hub is, the larger the maximum displacements would be.
the blade before bird strike. After bird strike, the maximum displacements increase significantly to 3.97, 1.41, 1.13, and 0.38, respectively, for the noted positions. As the bird strike force directly acts on the fan hub, the closer the bird impact point to the fan hub is, the larger the maximum displacements would be. To study the response of the rotor system at different times, the position displacement of the rotor system at five typical instances is demonstrated in Figure 14. The effect of bird strikes on the horizontal displacement of the rotor system is obvious from the graph. At = 2 s, the bird strike has not yet occurred, and the displacement of the rotor is small. When the bird strike occurs at = 3 s, in a very short time, i.e., at = 3.01 s, the overall displacement level of the rotor increases significantly. Similar to the previous conclusion, the closer the impact position to the fan position is, the greater the magnitude of the displacements increase would be. At bearing 3, which is far from the fan, the variation of displacement is very small, and the maximum displacement does not exceed 0.38. In practical studies, special attention should be paid to the response at the bearing close to the fan. At = 3.5 s, the vibration amplitude starts to drop rapidly due to damping effect, reaching 75% at the fan. Finally, at = 4.5 s, the displacement curve is almost identical to that at the initial steady state. Therefore, bird strike is very likely to cause damage to the engine in an instant. To study the response of the rotor system at different times, the position displacement of the rotor system at five typical instances is demonstrated in Figure 14. The effect of bird strikes on the horizontal displacement of the rotor system is obvious from the graph. At t = 2 s, the bird strike has not yet occurred, and the displacement of the rotor is small. When the bird strike occurs at t = 3 s, in a very short time, i.e., at t = 3.01 s, the overall displacement level of the rotor increases significantly. Similar to the previous conclusion, the closer the impact position to the fan position is, the greater the magnitude of the displacements increase would be. At bearing 3, which is far from the fan, the variation of displacement is very small, and the maximum displacement does not exceed 0.38. In practical studies, special attention should be paid to the response at the bearing close to the fan. At t = 3.5 s, the vibration amplitude starts to drop rapidly due to damping effect, reaching 75% at the fan. Finally, at t = 4.5 s, the displacement curve is almost identical to that at the initial steady state. Therefore, bird strike is very likely to cause damage to the engine in an instant.

The Influence of the Rotational Speed
The dynamic responses of the fan of the rotor system under three other rotational speeds are demonstrated in Figure 15, then a comparison between the dynamic responses was performed in conjunction with the calculation at 3772 r/min in Figure 12. The time history variations of displacement (Figure 15a-c) have the same trend as that presented in Figure 12. Due to the different magnitude of the bird strike force, the maximum value of each displacement has obvious differences and increases with rotational speed. Similar patterns can be observed in the axis orbit diagram (Figure 15d-f). The range of the axis orbit motion increases significantly with the increase in rotational speed. It can be seen from the phase diagram shown in Figure 15g-ithat the higher the speed is, the higher the number of limit rings would be. It can be seen in the Poincaré section diagram (Figure 15jl) that, as the speed increases, more points appear in the diagram, so more periods of motion exist. From what was said above, it can be concluded that as the rotor speed increases, the impact of bird strike becomes more drastic and the time required to go back to the steady state becomes longer. From the frequency spectrum diagram shown in Figure  15m-o, the fundamental frequency of the rotor system is always present with or without bird strike, and that bird strike only causes small fluctuations in the amplitude. It is noteworthy to say that the vibrational peak occurred at around 43 Hz. As the natural frequency of the rotor is also around 43 Hz, the first-order mode of the system is excited by the bird strike, which can cause serious damage to the engine. In Figure 15p-r, the time-frequency graphs show the same results. At different speeds, first-order frequencies are excited except for the fundamental frequency. From what was said above, it can be concluded that the effect of bird strikes on the rotor system only differs in magnitude at different speeds, and there is no significant difference in vibration characteristics.

The Influence of the Rotational Speed
The dynamic responses of the fan of the rotor system under three other rotational speeds are demonstrated in Figure 15, then a comparison between the dynamic responses was performed in conjunction with the calculation at 3772 r/min in Figure 12. The time history variations of displacement (Figure 15a-c) have the same trend as that presented in Figure 12. Due to the different magnitude of the bird strike force, the maximum value of each displacement has obvious differences and increases with rotational speed. Similar patterns can be observed in the axis orbit diagram (Figure 15d-f). The range of the axis orbit motion increases significantly with the increase in rotational speed. It can be seen from the phase diagram shown in Figure 15g-ithat the higher the speed is, the higher the number of limit rings would be. It can be seen in the Poincaré section diagram (Figure 15j-l) that, as the speed increases, more points appear in the diagram, so more periods of motion exist. From what was said above, it can be concluded that as the rotor speed increases, the impact of bird strike becomes more drastic and the time required to go back to the steady state becomes longer. From the frequency spectrum diagram shown in Figure 15m-o, the fundamental frequency of the rotor system is always present with or without bird strike, and that bird strike only causes small fluctuations in the amplitude. It is noteworthy to say that the vibrational peak occurred at around 43 Hz. As the natural frequency of the rotor is also around 43 Hz, the first-order mode of the system is excited by the bird strike, which can cause serious damage to the engine. In Figure 15p-r, the time-frequency graphs show the same results. At different speeds, first-order frequencies are excited except for the fundamental frequency. From what was said above, it can be concluded that the effect of bird strikes on the rotor system only differs in magnitude at different speeds, and there is no significant difference in vibration characteristics.

Conclusions
The "0-2-1" rotor system model was modeled mathematically and combined with the bird strike force obtained from bird strike simulation in LS-DYNA; the responses of the rotor system under bird strike and at different rotational speeds were studied. The conclusions can be summarized as follows:

Conclusions
The "0-2-1" rotor system model was modeled mathematically and combined with the bird strike force obtained from bird strike simulation in LS-DYNA; the responses of the rotor system under bird strike and at different rotational speeds were studied. The conclusions can be summarized as follows: (1) In the process of the bird strike, due to the blade geometry, the impact force in the y direction is much smaller than that in the x direction at 3772 r/min, with a difference of about 99.5 kN. The same holds true for other operational speeds. As the bird body in the bird strike process is cut into multiple parts, there are several peaks in the bird impact force diagrams.
(2) At the climbing state, the dimensionless displacement in the x direction was 0.88 before bird strike, and it reached 14 after bird strike, which shows an increase of about 16 times. The motion type also changed from a single periodic type to irregular quasi-periodic type. The displacement of the rotor system also increased significantly, and the closer to the fan hub, the more obvious the increase in amplitude. Therefore, bird strike has a serious impact on the rotor operation.
(3) At different rotational speeds, the response of the rotor after bird strike is similar. The difference caused by speed is only reflected in the amplitude, as it leads to different impact forces. At 836 r/min, the displacement in the x direction increases from 0.05 before bird strike to 0.53 after bird strike, which is an increase of about 10 times. The increase is 16 times at 3772 r/min. The maximum value is much larger than that at a low rotational speed. Therefore, the higher the rotational speed, the greater the damage caused by a bird strike is.
(4) The frequency spectrums at all four speeds show a peak at 43 Hz. Combined with the natural frequency of the rotor, bird strike excites the first-order vibration of the rotor.
This paper does not simply analyze the process of bird impact on the fan blade, but combines the bird impact process with the study of rotor system dynamics to obtain the rotor system dynamics response when the rotor system is struck by a bird, which is of guiding significance for the dynamics analysis of aero engines under bird strike, and the conclusions obtained also provide ideas for the safety design of aero engines; for example, bearing force can be used to guide the design of bearing structure, which is very effective in suppressing the vibrations caused by bird strikes.

Data Availability Statement:
The data presented in this paper are available on request from the corresponding author.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.