A CFD–CSD coupling method for simulating the dynamic impact and expulsion of fragile foreign objects from the ‘inlet–bypass duct’ junction of a turboprop aircraft

Frequently appearing fragile hail, ice and other foreign objects pose great threats to advanced turboprop aircraft when they enter the intake system of the engine. However, it is difficult to simulate the effects of these objects as complicated aerodynamic and dynamic impact coupling phenomena are involved. Regarding a turboprop aircraft inlet with a bypass duct as the subject of the research reported in this article, a coupled Computational Fluid Dynamics–Computational Structural Dynamics (CFD–CSD) numerical simulation method is established based on an unstructured dynamic mesh and impact dynamics technology to solve the phenomena of ‘motion–collision–fragmentation–fragment motion’ for foreign objects such as fragile ice and hail entering the inlet. Based on airworthiness specifications, external hail and ice at the lip and inner lower wall of the inlet are modelled. Moreover, the dynamic motion of these foreign objects and their effects on the aerodynamic performance of the inlet system during the process are simulated and analysed in depth. The results show that high-speed hail breaks into small debris after collision, which may not cause a serious threat to the engine. In addition, the total pressure recovery coefficient and distortion rate are not heavily changed during the process. However, large pieces of ice at the lip of the inlet may lead to large-size fragments that impact the inner wall of the inlet, thus increasing the threat to the safety of the engine. Also, serious shielding and interference of large pieces debris with the inner flow in the dynamic process decrease the total pressure recovery coefficient while increasing the distortion rate, especially in the junction area between the inlet and the bypass duct. In particular, ice in the icing zone located in the inner lower part of the inlet should be of much more concern owing to the fact that ice in this region is more likely to fly directly into the main engine and cause a serious threat.


Introduction
Compared with other types of aircraft, turboprop aircraft (Marinus & Maison, 2016;Q.F. Wang et al., 2014) has been widely used in the field of military and civil aviation owing to their economy, safety, comfort and environmental friendliness. There are more than 5000 turboprop aircrafts in main and regional airlines. This indicates that the progress of aviation science and technology provides more opportunities and requirements for the development of new generations of turboprop aircraft.
Turboprop aircraft have good environmental adaptability (Ashenden & Marwitz, 1997). They can take off and land normally on airstrips such as clay runways, sand runways and runways covered by grass and snow, and they can also adapt to bad weather. That is to say, this CONTACT Mi Baigang mibaigang@163.com kind of aircraft ought to have enough tolerance to foreign objects that may appear on the runway or during flight, since these foreign objects may cause serious consequences when entering the inlet. In addition to affecting the efficiency of the inlet and reducing the performance of the engine, foreign objects may also damage the internal structure of the inlet or the engine and threaten flight safety. For this reason, in recent years, a kind of inlet with a bypass duct is gradually emerging and has been applied in a variety of advanced turboprop aircraft (Mi & Zhan, 2019) such as the DHC-8 turboprop aircraft and MA700 turboprop regional airliner in China. The bypass duct refers to a section of pipe connected to the outside of the engine at the end of the intake port of turboprop aircraft, which is mainly used to remove large foreign objects that cannot be swallowed by the engine. The application of this type of inlet-bypass duct is of great significance for improving the performance of advanced turboprop aircraft. Common foreign objects affecting turboprop aircraft include sand, components, birds, hail, ice, raindrops and dust (Lu et al., 2020;Marandi et al., 2014). Among these, fragile hail and ice have a high probability of occurrence. A complex 'motion-collision-fragmentation' phenomenon can easily be induced when they enter the inlet, which may seriously affect the performance of the turboprop aircraft's inlet. Therefore, foreign objects are of especial concern in current aircraft design processes. Rzadkowski et al. (2018) and Szczepanik et al. (2010) investigated, in a series of studies, the effects of foreign objects on the structure and aeroelastic behavior of rotor blades. These results well inform the design and evaluation of aircraft components faced with threats from various objects. L.M. Wang et al. (2020) and Mi (2020) have carried out several analyses of the aerodynamic effects of rigid foreign objects on the inlet-bypass ducts of turbo-aircraft. They developed a single CFD method to simulate the whole dynamic process, and found that the dynamic motion of rigid objects heavily reduced the performance of the inlet. However, the method cannot be applied to simulating more complex non-rigid foreign objects. Ignatyev et al. (2020) proposed a series of wind tunnel tests to evaluate the longitudinal steady and unsteady aerodynamic characteristics of a transport aircraft in icing conditions. The study showed that ice shapes lead to reduced stall Angle Of Attack (AOA), maximum lift and longitudinal damping; however, the ice studied was fixed to the surface of the aircraft and possible motion during the flight process was ignored. Costes and Moens (2019) presented a numerical method to investigate the flow around an NACA23012 airfoil with two ice shapes. They compared two turbulence models, the S-A model and the DRSM model, to evaluate their abilities to capture the detailed vortex shedding phenomenon of the airfoil with ice. This study did not consider the ice motion after it dropped from the airfoil, which limited the application of the method. Jin and Taghavi (2008) and Jin et al. (2011) used STAR-CCM+ code to evaluate the effect of flow angularity on an S-duct inlet with icing. Symmetrical and asymmetrical glaze ice shapes were computationally simulated on the inlet lip. The results showed that the flow angularity aggravated the low performance of inlets with icing. Although the study investigated the icing effects at different locations, it did not consider the motion of pieces of ice or their effects on inlet performance. Lee et al. (2020) developed a Large Eddy Simulation (LES) method to investigate the flows over two types of iced airfoils with three multi-elements. The results demonstrated that LES simulation can capture much more detailed flow phenomena on iced airfoils than general Reynolds Averaged Navier-Stokes (RANS) models. However, the high costs of calculation in the study could not be averted. Taghavi and Jin (2008) used a steady-state RANS calculation to simulate the effects of typical rime and glaze ice on the performance of the M2129 S-duct inlet. Different shapes of ice could cause different losses in the mass flow and total pressure recovery of the inlet. The conclusions of that paper were similar to other studies, and they only focused on the effect of ice shapes on the aerodynamic performance of the inlet. Most of the above studies only evaluated the change of aerodynamic characteristics of the inlet or aircraft after external freezing, ice growth on aircraft components and the influence of different ice types on the characteristics of the external flow field or inlet performance. However, few studies have focused on the dynamic motion of ice and hail in the inlet and their effects on the aerodynamic performance of the inlet, even though the separate analysis techniques of Computational Fluid Dynamics (CFD) and Computational Structural Dynamics (CSD) have been developed in depth and validated (Abadi et al., 2020;Ghalandari et al., 2019), which is mainly because of complicated coupling phenomena between internal flow and dynamic impact. The internal aerodynamic performance of the intake system is heavily affected by the foreign object's motion, and the motion trajectory and the possible impact on the intake wall are also determined by the inner flow. The nonlinear coupling process strongly challenges the simulation technique.
In view of the current research gap in the investigation of the expulsion characteristics of icy foreign objects in the inlet bypass system of turboprop aircraft, this paper innovatively studies a numerical simulation method for investigating the expulsion characteristics of fragile hail and ice in the inlet bypass system based on the nacelle model of turboprop aircraft. Foreign object models of hail and ice are established taking into consideration airworthiness specifications. The simulation calculation of the possible 'motion-collision-fragmentation' phenomenon after the iced foreign object enters the inlet system is realized by coupling the methods of CFD and CSD. Moreover, the influence of the foreign objects on the inlet performance and the threat to the engine are evaluated. The article provides a reference for the design and airworthiness research of new inlet systems for turboprop aircraft.

Modeling fragile foreign objects based on airworthiness specifications
For turboprop aircraft, 'fragile foreign objects' mainly refer to ice accretion and hail occurring during flight, although glassy foreign objects appearing during takeoff and landing can also be considered fragile foreign objects (Chen & Hutchinson, 2002). The probabilities of ice accretion and hail occurring are much higher. Therefore, these two eventualities are chosen for the proposed related simulation methods.

Typical inlet configuration of turboprop aircraft
The calculation model in this paper is a typical turboprop aircraft nacelle configuration, as shown in Figure 1. The bypass duct is located at the end of the inlet, and its outlet directly leads to the external atmosphere. The inlet area is 0.16 m2. Since the maximum swallowing capacity of the engine is 0.20 m × 0.11 × 0.01 m, the bypass duct of the inlet must eliminate foreign objects exceeding this size. The heart section between the front of the inlet and the main engine piping is called the 'characteristic section'. As the airflow in this section needs to enter the main engine through the flow-guide vane, it impacts the performance of the whole engine significantly. In the subsequent calculation and analysis, the total pressure recovery coefficient and total pressure distortion rate on this section are extracted to evaluate the performance of the inlet.

Properties of hail
Based on the national standard GB/T 27957-2011 (China Meteorological Administration, 2011), hail is described as hard solid precipitation that is spherical, conical or irregular. Based on the diameter, the hail can be divided into four levels: small hail (diameter d < 5 mm), medium hail (diameter 5 mm ≤ d < 20 mm), large hail (diameter 20 mm ≤ d < 50 mm) and outsize hail (D ≥ 50 mm). In this paper, hail is regarded as a regular spherical object, and its texture as transparent and clean without bubbles. Natural hail density ranges from 870 to 920 kg/m3. The average value, 895 kg/m3, is considered here.
The distribution of hail inhaled by a turboprop aircraft inlet can be tested with reference to airworthiness standard CCAR-33R2 (Li, 2011), which prescribes that, when an aircraft is flying in turbulence with a maximum flight altitude of 4500 m (15,000 feet), the engine shall not cause unacceptable mechanical damage or unacceptable power or thrust loss or require the engine to stop after inhaling heavy hail (of specific gravity 0.8-0.9) at the maximum real airspeed, under the condition of maximum continuous power. At this time, half of the hail should be randomly projected to the whole area in front of the inlet, and the other half should be projected to the key area in front of the inlet. The hail should be rapidly and continuously inhaled, and the quantity and size of hail should be determined as follows.
(1) For an engine with an inlet area of less than 0.064 m 2 , the diameter of a hail particle is 25 mm; (2) For an engine with an inlet area greater than 0.064 m 2 , one hail particle of diameter 25 mm and one of diameter 50 mm for every 0.0968 m 2 of inlet area or the remaining area.
The capturing area of the inlet is about 0.16 m 2 . According to the above standards, two uniform hail particles of diameters 25 and 50 mm can be arranged at the centerline of the inlet, as shown in Figure 2.
The hail is formed by the condensation of cumulonimbus, and then falls to the ground under the action of gravity. Therefore, the hail encountered in the flight has a certain falling velocity, and the influence of flight velocity is coupled when it is inhaled into the inlet. Considering the influence of the air drag, the final speed of the hail (Heymsfield et al., 2018;Zhou, 2011) V can be expressed as (1)  where ρ i , ρ a are the densities of the air and hail, respectively. D is the diameter of the hail, g is the acceleration due to gravity and C D is the drag coefficient of the hail. The flight altitude and Mach number of the aircraft in this paper are determined as 3 km and 0.338, respectively, then the falling velocity of the hail is 17.3 m/s (25 mm hail) and 24.5 m/s (50 mm hail), respectively. For the inlet in this paper, the velocity of the inner flow is basically constant in the range Mach 0.2-0.3 to ensure the stability of the performance. The relative horizontal velocity formed by the influence of suction on hail should also be between the internal velocity and the flight velocity. For simplicity, the average value of these two velocities is adopted. Therefore, the final absolute velocity of the hail should be the vector sum of the flying speed and the falling speed.

Properties of ice
The main forms of aircraft icing can be generally divided into three categories: dry icing, sublimation icing and droplet icing. One of the most common icing forms is droplet icing, which is the phenomenon that supercooled water droplets in the atmosphere freeze after hitting the fuselage surface. There are three main types of this icing shape: rime ice, glaze ice and mixed ice, as shown in Figure 3. Rime ice is mainly formed on the leading edge of the wing, and in the cloud layer with small supercooled water droplets. Its surface is fragile and it falls off easily. Glaze ice has a smooth surface, grows along the airfoil surface and distributes widely. It is formed owing to the high temperature and large size of supercooled water droplets. Mixed ice is one of the most common ice types, which occurs owing to the partial freezing of impinging supercooled droplets that deposit onto the growing interface as rime particles within a mushy layer. The mushy layer acts as a porous medium containing solid rime particles saturated by water. Underneath this, at the water-ice interface, the solidifying liquid then entrains the rime particles into a glaze matrix to form mixed ice (Janjua et al., 2015;Worster, 2000).
According to Article 77 of CCAR33-R2, the test regulation of engine ice absorption is 'The amount of ice should be the maximum amount of ice accumulated in a typical inlet fairing and in front of the engine due to a two-minute delay in opening the anti-icing system, or using a piece of ice with mass and thickness comparable to the size of the engine. The ice suction velocity should be able to simulate the velocity of the ice sucked into the engine inlet, and the engine should work at maximum cruise power and thrust state. The ice absorption test should be able to simulate the maximum continuous icing condition encountered at −4°C (25°F)'. The engine test specification (G.H.  also stipulates the requirements of engine ice swallowing: 'In [a] typical take-off (maximum), cruise and landing, there is one hail particle of 5 cm diameter and two of 2.5 cm diameter, with a density of 0.80 ∼ 0.90 g/cm3 for every 0.25 m 2 or a part of the intake area of the engine inlet section. During take-off and cruise, a large amount of ice formed in the inlet and leading edge may be swallowed. Its density is 0.8 ∼ 0.9 g/cm, and its size is 76 mm×229 mm×6.4 mm generally, and only one piece of ice is swallowed in each test.
When the aircraft flies in icing weather conditions, the leading edge of the inlet, the front fairing of the compressor and the guide vane of the first stage compressor may freeze. In addition, the leading edge and lip of the propeller of a turboprop aircraft may also freeze, as shown in Figure 4. With regards to the nacelle model studied in this paper, blade icing is ignored and only ice accretion at the inlet lip is considered because of ignorance of the influence of propeller power. Additionally, there is an ice accretion zone in the inlet of the engine ( Figure 5). Therefore, ice accretion and falling off at this part are also considered. The maximum expulsion capacity of the model studied in this paper is 203.2 mm×114.3 mm×9.525 mm ( Figure 6), and this size limitation is determined by the    Wang et al., 2020). The ice size required by the airworthiness standard is larger than this benchmark; therefore, the size of the selected ice model is the maximum expulsion size. The shape of the ice is modified according to the cross section of the three types of ice, the shape of the lip and the surface.
The ice falls off from the components; therefore, its initial velocity is the same as the fuselage. Under the scour of the air flow in the inlet, it falls off from the fuselage surface. At the initial time, the ice starts to tilt up under the action of air flow, and it starts to fall off when the force reaches a certain extent. Hence, the initial attitude of the ice is adjusted to an angle of 30°. At this time, the force is greater, and the movement trend is more obvious. Actually, this process has also been testified and validated in previous studies of the author (Mi, 2020;Mi & Zhan, 2019;L.M. Wang et al., 2020). In the beginning, the ice just rotates with the action of the aerodynamic forces. When the ice angle increases to about 30°, it starts to speed up to the inlet, and then the aerodynamic performance of the inlet begins to change significantly. Therefore, the initial stage of the ice motion is ignored to save computing costs. Figures 7 and 8 illustrate the modified ice calculation model for the ice accretion zone inside the lip and inlet.

Numerical simulation method
The fragile foreign objects move under the impact of the internal air flow after being sucked into the inlet. As the effect of their physical properties, these objects may break up after colliding with the inner wall of the inlet, as shown in Figure 9. From the perspective of numerical simulation, the original topological structure changes qualitatively, and this change of the topological structure relates to the movement of ice, which cannot be set in advance. Therefore, it is difficult to complete the simulation analysis of the process by using the CFD method   alone, and joint simulation needs to be completed by coupling CFD and CSD. The specific process is shown in Figure 10.
Specifically, the steps of simulation analysis are as follows.

Step 1: Foreign objects entering the inlet
Hail is inhaled from the outside of the inlet. It has horizontal velocity and falling velocity relative to the fuselage. Thus, dynamic mesh technology is adopted combined with Six-Degree-Of-Freedom (6-DOF) technology to simulate its entry process. The ice is attached to the inlet components. Relative to the fuselage, it is static. Under the constant scour of the air flow, ice tilts up and falls off from the components and moves into the inlet owing to the impact of aerodynamic force. Its motion attributes can also be described with 6-DOF equations; therefore, the simulation method is similar to hail.
(1) 6-DOF motion When the foreign object enters the inlet, it will move dynamically with coupled translation and rotation. In this paper, the unstructured dynamic mesh deformation where˙ v G is the movement velocity of the center of the control object, m is the mass of the motion control body and f G is the force vector related to gravity. The angular velocity vector of the body axis is easier to calculate, as shown in Equation (3): where L represents the inertia tensor, M B is the moment vector of the control object and ω B is the angular velocity vector.
The torque conversion mode in the inertial frame and body axis system is as follows: where R is the transformation matrix. Specifically, where C χ = cos(χ ),S χ = sin(χ ), and φ, θ, ψ represent Euler angles (around the Z-, Y-and X-axes), respectively.
Although the particle tracking technique has already been integrated in the Discrete Particle Method (DPM) model in ANSYS R Fluent software, the physical model is not the same as that in this study. The volume of the broken debris cannot be ignored, and secondary collisions among these pieces should also be considered; therefore, the 6-DOF method is applied to obtain the trajectory of the debris.
(2) Unstructured dynamic mesh technology When foreign objects move in the flow field, the unstructured mesh reconstruction technology in Fluent (Batina, 1991) is used to update the grid, which is more suitable for conditions of large deformation and large amplitude motion. The basic approach is to regenerate the grid in areas having large grid deformation within the set time steps. This approach is simple. However, it is necessary to ensure that the grid distortion rate is within the default range.

(3) Numerical simulation of unsteady aerodynamics
In the whole process of foreign objects expulsion, the flow field is constantly changing and highly unsteady. Hence, an unsteady aerodynamic numerical simulation method based on the Reynolds average is used to solve the problem. The turbulence model adopted is the SST k-ω two-equation model (Menter, 1994), and its basic equation is shown as Equation (5). The model has good adaptability for the simulation of small and medium separations, and is widely used in the calculation of aircraft flow fields.
Step 2: Collision judgement Solving the motion of the foreign object in the inlet by CFD, it is necessary to monitor the distance between the foreign object and the inner wall of the inlet in real time to judge or define the occurrence of collision. If there are multiple foreign objects, monitoring the distance between the foreign objects is also needed. When the minimum distance decreases continuously and is less than a set threshold δs, we consider the occurrence of collision. This threshold is greater than zero because the calculation will fail if the objects contact each other to form a singularity in the CFD calculation. (

1) Collision judgement of spherical hail
The collision judgement of spherical hail only needs to calculate the minimum distance between the hail sphere and the inner wall of the inlet and the threshold value d, as shown in Figure 11. When the distance decreases continuously and is less than the sum of the radius of the ball and the threshold, i.e. d ≤ r + δs, collision occurs.

(2) Collision judgement of ice
The shape of the ice is irregular; therefore, the collision judgement needs to calculate the distance of all the grid centers between the objects and the inner wall of the inlet, and take the minimum value, d min . When the minimum distance decreases continuously and is less than the threshold distance, i.e. d min ≤ δs, collision occurs, as shown in Figure 12.

Step 3: Foreign object status capture and flow field information storage before collision
Once it is judged that collision will occur in the next time step, it is necessary to suspend the CFD flow field calculation and save the current flow field information. In other words, the instantaneous attitude, position, translational and rotational velocity of the foreign object, as well as the pressure and velocity of other points in the flow field, are stored, i.e. transient flow field freezing. This information will be used to interpolate the initial flow field of the new foreign object motion simulation after the collision. In addition, in order to conduct the collision simulation, it is required to output the surface grid of the foreign object and the internal wall of the inlet, and reconstruct the geometry to form the collision model.

Step 4: Collision calculation
The finite mesh of the model outputted and reconstructed from the CFD result is generated, and the relevant property data of the foreign objects are reloaded. Setting the contact conditions, the collision process is simulated by the Ls-dyna program (Liu et al., 2014;Tian & Zhu, 2011). Information on the debris after collision is obtained, mainly including the volume, mass inertia, translational and rotational velocity, and surface geometry information.
The hail or ice impact process is a typical dynamic problem of large deformation on elastic-plastic materials. The basic kinetic equation is expressed as where F int represents the internal force vector, which is a nonlinear function of displacement u. F ext is the external force vector. M and C are the mass and damping matrices, respectively.
More detailed information on impact theory can be found in the users' manual for Ls-dyna software (Hallquist, 2003).

(1) Collision calculation method
There are three main models for Ls-dyna: the Finite Element (FE) method, the Smoothed Particle Hydrodynamics (SPH) method and the Arbitrary Lagrangian Eulerian (ALE) method. In general, the accuracy of the ALE model is low and the calculation time is long. The SPH method needs to use particle masses instead of a grid, which is not suitable for coupling with CFD calculations with grids in this paper. Therefore, the finite element method is used to solve the problem.
MAT24 material is used as the wall aluminum alloy model, and MAT_ISOTROPIC_ELASTIC_FAILURE (MAT13) is selected as the hail material in this paper. The contact in simulation is set as nodes-surface type, in which the slave surface is hail and the master surface is the engine wall. This setting mainly considers that, when hail and ice move near the wall, most contacts are point contacts. The collusion characteristics of hail and ice are calculated and analyzed in advance, and the threshold for collision is set as when the distance between the foreign object and the wall is 0.004 m.

Step5: Judgement of foreign object state and model reconstruction after collision
The equivalent size of debris after collision simulation is determined. If it is larger than the tolerance size of the inlet, a further calculation should be continued to track its trajectory. If it is smaller than the size limit, it can be considered that, even if it is inhaled into the main engine, it will have no significant impact and can be ignored. If fragments of the foreign object are larger than tolerance size, node information should be outputted and the model needs to be reconstructed to form a new foreign object model.

Step 6: Motion continuation of reconstructed foreign object fragments in the inlet
In order to transfer the parameters from a CSD calculation to CFD simulation, the new foreign object is meshed and loaded with the corresponding information on the foreign object, including mass inertia, position of the gravity center, initial velocity, etc. As a continuous motion and collision process, the CFD calculation at this time is based on the flow field before collision. Thus, it needs to use the flow field information in step 3 to interpolate to form the initial flow field in this step to calculate the motion of the new foreign object. Please note that, as the topologies of the foreign object before and after collision change strongly, the relative flow fields will also have been seriously affected (as shown in Figure 13). The complete configuration of the foreign object before the collision may be impacted into several small pieces, which will lead to a different case in which the required parameters of the flow field around the impacted debris are missing (dotted lines area in Figure 13). Therefore, before the next CFD simulation, the parameters of the flow field such as velocity, pressure and temperature should be initialized by interpolating the CFD data before the previous collision. In this study, a linear interpolation algorithm has been applied to calculate the parameters of P3 with those values of P1 and P2. This simple treatment can effectively provide initial values for the next CFD simulation.
The process is repeated until the foreign object completely breaks up into small pieces or enters the bypass duct or main engine.

Method validation for simulating the aerodynamic performance of the inlet without foreign objects
Firstly, the numerical simulation method is verified by the nacelle inlet model without foreign object interference. The unstructured mesh is generated by ANSYS ICEM TM CFD software and the grid is refined at the inlet system ( Figure 14). The size of the basic grid is six million nodes, and several other fine or coarse grids are used with total numbers of nodes of four, five, eight and ten million to validate the mesh independence. The main difference between these grids is the mesh size in the refined region of the inlet-bypass duct. The nacelle inlet surface grid is shown in Figure 15.   The Mach number selected for verification is 0.2. The periphery of the calculation domain is the far-field boundary condition. The nacelle inlet is set as the nonslip wall boundary. The outlet of the intake system is considered as the pressure outlet boundary. According to the test value, the static pressure here is set as 90527 Pa. Three parameters are used to evaluate the performance of the intake system: the outlet mass flow, the total pressure recovery coefficient and the total pressure distortion rate.
The calculation formula for the total pressure recovery coefficient is as follows: where P 1 is the average total pressure of the mass flow in the characteristic section. P 0 is the average total pressure of the mass flow of the free stream. The total pressure distortion rate is defined as follows: where P 1,max and P 1,min represent the maximum and minimum total pressure in the characteristic section, respectively. The calculated results within different grids are shown in Table 1. The results show that the basic grid with six million nodes can obtain small calculation errors and relative high efficiency, and the same grid generation strategy will be used to simulate all the other cases in the study.

Expulsion of spherical hail from the inlet
According to the airworthiness specification, four hail particles should be arranged at the entrance of the inlet in the model. In the process of elimination, there may be some complex phenomena, such as successive collision, collision with the inner wall, collision with each other and debris collision. This may lead to huge costs when using CFD/CSD co-simulation. This paper only establishes a general method. Therefore, only one large hail particle is set in the inlet to develop and test this method. The   Figure 16. The specific parameters of the hail are shown in Table 2.

(1) Flow field calculation before the first collision
The six-degree-of-freedom equation and the unstructured grid reconstruction technique are used to calculate the foreign object motion before collision. Figure 17 shows the calculated trajectory of the foreign object. When the calculation time is t = 0.004 s, the minimum distance between the hail particle and the lower wall is 0.0045 m, and the collision occurs. The location of the hail particle's gravity center is (11.92056,−4.4464, 0.1946364), the translational velocity is (69.97915, 0.002392371, −24.46934) and the rotational angular velocity is (−0.02305, −0.1474214,0). Since the initial rotation effect of the hail particle is not considered, the rolling effect caused by the air force is relatively small.

(2) First collision simulation
Based on the geometries of the foreign object and the inner wall of the inlet from the CFD data, a finite element calculation grid is generated. Inputting the mass inertia characteristic parameters of the foreign object, the collision process is calculated by Ls-dyna. The results are shown in Figure 18. Owing to the high-speed collision, the hail particle disperses into small pieces (other smaller pieces have been destroyed), and the separated larger pieces continue to move backward. The size of the small fragments is actually smaller than the swallowing size of the engine. Even if it is inhaled into the engine, it will not cause serious damage. However, in order to verify the calculation method in this paper, the simulation of the dynamic motion of this hail particle is continued. The maximum three-dimensional dimensions of the small fragments are (0.028, 0.024, 0.003 m), and the velocities in three directions are (55.257, 1.24, 16.49) ( Figure 19). Its mass is 0.000,69 kg, the center of gravity is (11.98, −4.446, 0.211) and the corresponding moments of inertia are I xx = 1.677e − 08kg · m 2 , I yy = 2.182e − 08kg · m 2 , Figure 17. Dynamic motion of a hail particle before the first collision: (a) initial position of the hail particle; (b) position of the first collision. Figure 18. Morphological changes of a hail particle before and after collision: (a) before collision; (b) after collision. I zz = 3.79e − 08kg · m 2 , I xy = 8.508e − 12kg · m 2 ,I xz = 1.053e − 12kg · m 2 , I yz = 1.864e − 11kg · m 2 respectively

(3) Continued flow field simulation after the first collision
After the first collision, the original spherical hail particle breaks into pieces, and the shape of the remaining larger fragments is shown in Figure 20(a). The model is remodeled and meshed to continue the flow field calculation. Figure 21 illustrates the initial and final foreign object positions in the CFD calculation after the first collision, after which the foreign object continues to move to the inlet by the action of air flow. The second collision occurs 0.016 s later at the entrance of the bypass duct. The fragments will continue to break into smaller ice fragments, and will be sucked into the bypass duct and then eliminated out of the inlet. At this time it is considered that the expulsion simulation of the hail can be ended. Figure 22 shows the motion trajectory of hail during the whole process. Figure 23 shows the effect of hail motion on inlet performance. As a whole, the existence of hail is a kind of disturbance to the inlet. When it moves into the inlet and impacts the inner wall, the influence on the flow field of the subsonic inlet is more obvious. The mass flow first decreases and then increases. The main reason for the increase is that the debris becomes smaller after the first collision, and the overall air flow shielding becomes smaller. However, compared with the case without foreign object interference, the mass flow is still reduced by about 11.9%. The variation law of the total pressure recovery coefficient is similar to that of the mass flow, and the maximum decrease is 0.93%. The distortion rate relates to the distance between the foreign object and the characteristic section. The smaller the distance, the faster the disturbance propagates and the greater the influence. When the fragment moves to the bifurcation region of the duct at the end of the inlet, the distortion rate noticeably increases, and is about 200% higher than that in the case without interference. This paper investigates only one hail particle in the inlet, and its influence is relatively small. If multiple hail particles are arranged according to airworthiness standards, the flow field shielding and  collision phenomena will be more complex, and their influence on the inlet performance will be more obvious.

Expulsion of irregular ice from the inlet
The calculation model and grid for the foreign object in the lower lip are shown in Figure 24. The ice without initial velocity starts to tilt and move inward under the action of aerodynamic forces. The parameters of the ice at the lower lip are given in Table 3.
(1) Flow field calculation before the first collision The coupled simulation method in the previous section is used to calculate the motion of ice without initial velocity. Firstly, the flow field before the first collision is calculated, and the results are shown in Figure 25. Since it is assumed that the ice is tilted 30°in the calculation, the ice is pushed to the inlet by aerodynamic forces. The ice may also rotate and roll owing to the uneven force caused by the shape. Before the first collision, the ice basically moves upwards and produces angular velocity around the    (11.685,45 m, −4.38, 0.208,9234 m) Inertia characteristics I xx = 0.000, 4792kg · m 2 , I yy = 0.000, 142kg · m 2 , I zz = 0.000, 6155kg · m 2 , I xy = 0, I xz = 0.000, 012kg · m 2 , I yz = 0 Initial velocity 0 center of gravity. The position of the center of gravity of the first collision is (11.795,08, −4.379,154, 0.396,411), the velocity and angular velocity of the first collision are (8.230,658, 0.111,8347, 13.086,44) and (−0.808,1403, 18.555,42, 0.087,0744), respectively.

(2) First collision calculation
The parameters of the foreign object model before collision are outputted and remeshed, and then the collision is simulated by Ls-dyna. Figure 26 shows the change of ice shape before and after the collision. The results show that the ice fragments impact the upper wall of the inlet with a low velocity and break into several fragments of different sizes. The largest fragment size is 136 mm × 24 mm × 5 mm, which is significantly smaller than the swallowing limit of the engine, but its subsequent motion is still calculated and analyzed. The mass of the largest fragment is 0.015 kg. The center of gravity is (11.842, −4.344, 0.416). The velocity is (V x = 8.6 m/s, V y = 0.146 m/s, V z = −1.78 m/s) and the moments of inertia are I xx = 2.275e − 05kg · m 2 , I yy = 7.356e − 07kg · m 2 , I zz = 2.327e − 05kg · m 2 , I xy = −2.22e − 16kg · m 2 , I xz = 2.168e − 7kg · m 2 , I yz = −3.469e − 18kg · m 2 respectively ( Figure 27).

(3) Flow field calculation after the first collision
The largest fragment is selected to reconstruct its geometric model. Owing to the lack of elements after collision, the shape of the largest fragment is heavily irregular, which makes the CFD modeling and meshing inconvenient. Therefore, the geometry of the fragment is scanned, small defects ignored and a simplified model established, as shown in Figure 28. The unstructured grid is shown in Figure 29.
The CFD method coupled with the 6-DOF equations are used to calculate the motion of ice fragments after the first collision. Figure 30 shows the instantaneous positions of the object at the initial time (t = 0.0282 s), the intermediate time (t = 0.0397 s, 0.0522 s) and the second collision time (t = 0.0687 s). The results show that the second collision occurs at the inner wall of the bypass duct entrance, i.e. the debris has entered the bypass duct and will be discharged from the inlet after the second collision. Hence, the calculation can be considered to be finished. The motion trajectory shows that the fragment    rotates to the inside under the action of air flow and collision speed, and the orientation of the foreign object changes dramatically in the whole movement process, especially in the fork area of the pipeline, which indicates that the flow mechanism in this area can be much more complicated.
The influence of ice movement and collision processes on the inlet performance after falling off from the inlet lip is shown in Figure 31. As the ice locates at the lower part of the inlet lip initially, and it is set at a certain angle, thus the intake effect is shielded at the initial time, and separated flow is formed on the leeward side. This leads to a small mass flow at the initial time. With the increase of the warping angle caused by movement of air forces, the shielding area and separation become serious. Therefore, the mass flow continues to decrease. The minimum mass flow is only 80% of that without interference, which leads to a smaller final total pressure recovery coefficient at the outlet of the intake system. But at this time, the foreign object is far away from the characteristic section, and the disturbance to the air flow has recovered well at the characteristic section. Hence, the distortion rate is not large. When the foreign object moves to the fork area, the disturbance is stronger when it is closer to the characteristic section, then the distortion rate notably increases, which is about 260% higher than that without interference. At this time, the influence on the inlet is reduced, and both the total pressure recovery coefficient and the mass flow are increased.

Expulsion of ice from the icing zone near the inner lower wall of the inlet
The ice on the inner lower wall of the inlet tilts up under the action of the air flow, and then moves inwards under the effect of the aerodynamic forces. The dynamic motion may be coupled translation and rotation. Figure 32 shows the calculation model and grid, and the calculation process is consistent with the previous section. The mass of the ice on the lower wall icing area is 0.124 kg. The center of gravity is (11.91782, −4.38, 0.20988), and the inertia is I xx = 0.000, 4435kg · m 2 , I yy = 0.000, 09766kg · m 2 , I zz = 0.000, 5088kg · m 2 I xy = 0, I yz = 0 I xz = 0.000, 03561kg · m 2 Figure 33 illustrates the calculated ice trajectory. It can be seen that the ice in the lower icing area is special. Instead of colliding with the upper and lower walls, it flys directly into the main engine and is not expelled through the bypass duct. This is the most dangerous situation, which needs to be of high concern. In the whole process, the ice curls up and floats towards the inlet. At the same time, owing to the uneven force, the ice rolls around the center of gravity and the rotation velocity around the local Y-axis is the highest. Figure 34 shows how the inlet performance changes during the whole process of ice movement. The calculation time is 0.032 s. When a foreign object locates in the icing area at the lower wall of the inlet, the mass flow noticeably decreases. With the different motion orientation of the foreign object, the perceived flow rate at the outlet of the intake system also    changes, and the difference of the flow rate in the calculation time reaches to 0.6 kg/s. In consequence, the total pressure recovery coefficient of the characteristic section also changes, and its maximum value decreases to 0.95, which is 95.8% of that without interference. The results show that the distortion rate increases with the presence of a foreign object, especially when the foreign object reaches the characteristic section, the disturbance to the air flow nearby increases and the distortion rate changes more dramatically. The maximum distortion rate reaches 9.5%, which is 389% higher than that without interference. The entry of this kind of foreign object directly into the engine not only seriously affects the safety of the engine, but also greatly reduces the performance of the inlet.

Conclusion
By coupling the unsteady CFD technique, the 6-DOF method and CSD impact theory, this paper establishes a numerical simulation method for the entry of two types fragile foreign object, ice and hail, into the nacelle inlet of an advanced turboprop aircraft with a new 'inlet-bypass duct' system in order to evaluate their expulsion characteristics and the performance effects on the inlet. The dynamic motions of the external hail particles, the ice near the lip of the nacelle, and the ice in the icing zone of the inlet, have been simulated in detail. Moreover, the total pressure recovery, mass flow rate and distortion rate of the 'inlet-bypass duct' system of the inlet during the expulsion process are also investigated. The results show the following.
(1) Hard hail particles have an initial falling velocity that is relatively fast when they enter the inlet system. Therefore, they are easily breaken into smaller pieces when they hit the inner wall and have a relatively small impact on the whole performance of the inlet. However, if a large number of hail particles enter into the inlet at the same time, they will greatly interfere with the performance of the inlet, even threatening flight safety.
(2) Ice at the lip is tilted and moves under the action of aerodynamic forces, which may produce large fragments beyond the endurance of the engine after collision. This leads to a serious block to the air flow in the inlet. Then the total pressure recovery coefficient and outlet mass flow decrease significantly compared with the non-interference situation, and the distortion rate increases significantly, especially when it moves to the fork position of the inlet-bypass duct/main engine inlet and the distortion rate increases by 260% compared with the non-interference situation, which seriously affects the inlet performance. (3) When the ice in the lower icing area is tilted, it may enter the main engine directly without collision, which is a great threat to the engine. The total pressure recovery coefficient decreases by 4.2% and the maximum distortion rate increases by 389% owing to the serious disturbance caused by the change of attitude. Therefore, the influence of this part of the ice should be monitored.
This paper only establishes a set of coupling simulation methods for calculating the expulsion characteristics of fragile foreign objects under typical working conditions. In fact, the dynamic expulsion process of a foreign object such as a fragile ice particle is much more complex. On the one hand, the properties of many kinds of ice vary owing to the physical mechanism of their condensation, which indicates that many more collision processes can appear, even on the same ice configuration, and further detailed and systematic analysis is proposed. On the other hand, the slipstream of a turboprop aircraft seriously interferes with the inlet, which may change the trajectory and expulsion characteristics of foreign objects. However, the corresponding numerical simulation methods for these phenomena perform consistently, i.e. the coupling method established in this  paper can be used for in-depth study. In future work, the coupled simulation method should be developed and validated by considering the influence of many more actual parameters. Moreover, combined with the actual probability of various foreign objects entering the intake system, the method can be further improved to support the aerodynamic optimization of the intake-bypass duct system.

Disclosure statement
The author declares that he has no financial or personal relationships with other people or organizations that could inappropriately influence his work, and that there is no professional or other personal interest of any nature or kind in any products, service and/or company that could be construed as influencing the position presented in, or the review of, the present manuscript.

Funding
This research was supported by the Fundamental Research Funds for the Central Universities, Northwestern Polytechnical University, Shaanxi province, China (grant number G2020KY05115) and the Natural Science Basic Research Program of Shaanxi, China (grant number 2021JQ-084).

Data availability statement
The data that support the findings of this study are available on request from the author.

Code availability
ANSYS Fluent and Ls-dyna codes were used in this research, and the coupling process was achieved by developing MATLAB code.