Dynamic Behavior and Heat Transfer Characteristics of Non-spherical Ice Crystals in High-Temperature Air Flow

，，，，，，，2*

1.School of Energy and Power Engineering，Xi’an Jiaotong University，Xi’an 710049，P.R.China；

2.Anti/De-Icing Key Laboratory，China Aerodynamics Research and Development Center（CARDC），Mianyang 621000，P.R.China

Abstract: This study investigates the dynamic behavior and heat transfer of ice crystals in high-temperature air flow after being ingested by an aero-engine.A numerical model under the Lagrange framework is built to simulate ice crystal motion，heat and mass transfer.Different ice crystal melting models and ice crystals with various shapes and sizes are also compared.The results show that when using the “naked” ice particle model，it takes longer to completely melt compared with the water cover model.It takes 0.04 s for a 20 μm particle to melt，but the melting rate of a “naked” ice particle is faster before 0.023 s.Under the same condition，the liquid water content of a spherical ice crystal is high，and that of a 20 μm spherical one is 49.05% at the outlet.The liquid water contents of ellipsoidal and hexagonal plate ice crystals are close，which are roughly 40% at the outlet.The smaller the size，the earlier the crystals start to melt，and the higher the liquid water content of the crystals.The liquid water content is 60.4% of a 20 μm ice crystal and only 15.5% of a 40 μm one at the outlet.

Key words：ice crystals；aero-engine；Lagrange method；heat and mass transfer；non-spherical

0 Introduction

Aircraft icing problems were recognized as early as the beginning of modern flight.The first reported aircraft icing can be traced back to 1920s［1］.Researchers first discovered and studied the phenomena of supercooled water droplets icing.Related prediction and protection have also been widely applied in practice［2］.However，since 1990s，there has emerged an increasing number of accidents caused by engine thrust loss at the altitude of 6.7 km even though supercooled water droplets rarely exist there.According to the accident investigation and analysis conducted by the relevant authorities，the main cause of accidents is the adhesion and freezing of ice crystals after being ingested by the engine［3］.

In recent years，intensive investigations have been conducted to resolve the problem of ice crystal icing inside engines via theoretical and numerical modelling，and experiment［4-5］.For the ice crystal melting and freezing process，Villedieu et al.［6］and Trontin et al.［7］proposed the ice crystal freezing model based on the supercooled water droplet freezing model，which divided the ice crystal melting phase transition into three processes，namely sublimation of supercooled ice crystal，ice crystal melting and evaporation of droplets after melting completed.Lou et al.［8］proposed three models of ice crystal melting，which are respectively the “naked”ice model，the water film cover model and the porous model.They compared the first two models and showed that the water film cover model had a faster melting rate.Norde［9］and Nilamdeen et al.［10］used the Euler method to simulate ice crystal freezing in airfoils and engines，considering the complex processes of ice crystal adhesion，fragmentation and erosion of ice crystals.Veres et al.［11-13］applied an engine icing computational tool to simulate the flow through a turbofan engine and assess the risk of ice accretion.They found that the ice accumulation was affected by the local wet bulb temperature，ice crystal speed and ice crystal liquid water content.Feulner et al.［14-15］used the Monte Carlo particle trajectory simulation method and studied the motion and crushing process of ice crystal particles ingested by the engine，and found that the ice crystal fragmentation mainly occurred on the low-pressure compressor blades.The internal flow field，ice crystal inertia and centrifugal effect of the compressor had great effects on the radial distribution of ice crystals.Nilamdeen et al.［16-17］extended the original supercooled water icing software and developed a mixed-phase icing calculation program for air，water droplets and ice crystals，and conducted the three-dimensional modeling on the ice crystal freezing in an engine to obtain the evolution process of ice shape on the surface of stator blades.

However，relevant research in China is still in its infancy stage.Especially，the motion，heat transfer and icing characteristics of ice crystal ingested in the turbofan engine are rarely studied.Tan［18］used the Euler method to simulate the freezing of airfoil ice crystals，used the Messinger model to obtain the ice shape，and found that freezing is easier to occur in a low-pressure environment.Bu et al.［19］and Guo et al.［20］modified the Messinger model to calculate the icing of NACA0012 airfoil under mixed-phase conditions.The results showed that liquid water content，airflow temperature and ice crystal adhesion effect have great influences on the amount and profile of icing.Jiang et al.［21］used the Lagrange method to study the motion，heat transfer and adhesion of ice crystals to the engine.They obtained ice profile of the surface area of stator.Guo et al.［22］simulated the movement and heat transfer of ice crystals in a typical large icy wind tunnel configuration based on Euler method，and clarified the thermal/force balance characteristics of medium ice crystals in large icy wind tunnels.Ma et al.［23］developed an airfoil icing model on the NNWICE platform，and analyzed the influence of ice crystal shape and flow parameters on the phase transition and adhesion of moving ice crystals.

Most of the reported studies simplify ice crystals as spherical particles.Ice crystals in clouds can undergo freezing and condensation，leading to complex shapes and size distributions.Therefore，in simulations of ice crystal freezing，it is crucial to take into account the impact of irregular particles which can affect the behavior of ice crystals.Naturally-formed ice crystals also have different melting phase transition processes under high-speed air flows，so it is necessary to compare and analyze different melting models.In this study，by using the Lagrange method，the dynamic behavior and heat and mass transfer of ice crystals after being ingested into an engine are investigated numerically.The heat transfer characteristics of ice crystals inside an aircraft engine are obtained.In the model，different melting models of crystals are compared，and the effects of shape irregularity of ice crystals and compressor blade thickness are examined.This study consolidates the foundation for further numerical simulation on ice crystal icing phenomena in aircraft engines.

1 Numerical Modeling

According to the two-phase flow theory，the icing condition of aircraft ice crystals is in the range of unidirectional coupled dilute two-phase flow，and the following assumptions are made.

（1）Only consider the effect of air flow on ice crystals，and neglect the effect of ice crystals to air flow.

（2）Ice crystals neither condense nor break in motion.

（3）Neglect internal temperature variation of an ice crystal，and the viscosity and density remain unchanged during the movement.

1.1 Motion and heat transfer model of ice crystals

In the engine flow channel，the ice crystal particles are affected by a variety of forces，including viscous forceFD，gravity and buoyancyFG，Magnus forceFM，Basset forceFB，etc.According to Newton’s second law

whereMiceis the mass of the ice crystal and duice/dtthe acceleration of the ice crystal.Miceis a function of the radius of ice crystals shown as

whereρiceis the ice density andRicethe radius of ice crystals.

Since ice density is much greater than that of air，and the size of an ice crystal is small，the viscous force is two to three orders of magnitude greater than gravity，only the viscous force is considered，and other forces such as gravity with smaller magnitude are negligible.Combined with the drag coefficient，the equation of motion can be expressed as

whereρa(bǔ)iris the air density，CDthe drag coefficient of ice crystals，uairthe airflow velocity，anduicethe ice crystal velocity.

Lawson et al.［24］measured and generalized the size and shape of ice crystal particles in cirrus clouds.The classification of the ice crystal shape is shown in Fig.1（a）.They are categorized into the following mathematical models：Spherical，ellipsoidal，columns，hexagonal plate ice crystals，as shown in Fig.1（b）.

Fig.1 Classification of the shape of ice crystal particles in clouds and simplifying mathematical models[24]

The sphericity is defined byΦ=Ap/A，whereApis the surface area of the sphere with the same volume as the ice crystal，andAthe actual surface area of the ice crystal.Φis used to measure the degree to which the non-spherical ice crystal approximates the spherical shape，and then further correct the drag coefficient and heat and mass transfer coefficient of the movement.By comparing theCDcorrection methods of various non-spherical particles proposed in the studies［9，25-26］，the empirical correlation of spherical，ellipsoidal and Columns is

whereA(Φ)=exp(2.328 8-6.458 1+2.448 6Φ2)；B(Φ)=0.096 4+0.556 5Φ；C(Φ)=exp(4.905-13.894 4Φ+18.422 2Φ2-10.259 9Φ3)；D(Φ)=exp(1.468 1+12.258 4Φ-20.732 2Φ2+15.885 5Φ3).

Calculation empirical correlation ofCDwith a hexagonal plate ice crystal

whereReis the relative Reynolds number.

The ice crystals that are sucked into the engine experience a temperature difference with the surrounding hot air.As a result，the outer surface of the ice crystal melts into a water film，which either falls off or covers the surface of the ice crystal and continues to evaporate.Due to this melting process，the diameter of the ice crystal particles decreases.This phase change melting process can be simplified into the following two idealized models.

The first model is the “naked” ice particle model，where the ice crystal particles are always in contact with air without any external water layer.In this model，once the outer surface melts，the resulting water film is blown away from the surface of the ice particles.The second model is the water cover model，where the water film resulting from melting covers the outer surface of the ice crystal particles，as shown in Fig.2.

Fig.2 Schematic diagram of two melting models

The initial stage of both models involves the sublimation process of supercooled ice crystals，with a heat transfer relationship as

wherehis the convective heat transfer coefficient on the surface of the ice crystal，Cpithe specific heat capacity of ice crystals，hmthe mass transfer coefficient on the surface of ice crystals，ρgathe density of water vapor in the air，ρgithe saturation vapor density on the surface of the ice crystal，Tairthe airflow temperature，Ticethe ice temperature，andLsgthe sublimation heat of ice.

For the “naked” ice model，after the ice crystals begin to melt，all the heat absorbed from the air is used for the melting process，and the energy balance equation is

whereLslis the melting heat of ice crystals.

For the water cover model，after the ice crystals begin to melt，the heat absorbed from the air is used for the melting of ice crystals and the evaporation of the water film.The energy conservation relationship therefore is

whereLlgis the heat of vaporization of water，Lslthe heat of melting of ice crystals，andMgthe mass of water film evaporated into vapor.

The mass of the ice-water mixture is the summation of those of ice and water，which is also equivalent to the mass of the ice-water mixture at the initial moment of melting minus the mass of evaporated water vapor

whereMmixis the mass of the ice-water mixture，Mwaterthe mass of water in the ice-water mixture，andMsmthe mass of the ice-water mixture at the initial moment of melting.

The equivalent densityρmixof ice-water mixtures can be determined by the mass of ice and water as

Assuming that the water film produced after the melting of ice crystals evenly covers the surface of the particles，the equivalent radiusRmixof the ice-water mixture can be calculated from the equivalent density

Similarly，the heat transfer of non-spherical particles is corrected，and the number of non-spherical particlesNuis calculated as［26］

wherePris the Prandtl number.

The mass transfer coefficient of ice crystal or ice-water mixture is obtained by heat and mass transfer analogy as

whereCpais the specific heat capacity of air，andLethe Lewis number.

The density of ice crystals is 919.1 kg/m3in this study，and the thermophysical parameters of water and air including density，viscosity，specific heat capacity and thermal conductivity are calculated according to temperature interpolation.

1.2 Engine flow channel model

This study uses the typical low-pressure stage compressor model of Ref.［27］for calculation.Fig.3 shows a five-stage compressor planar blade grid，with each stage composed of rotors and stators.

Fig.3 Schematic diagram of ice crystals passing through the interior of a low-pressure stage compressor[27]

The NACA airfoil family thickness distribution is adopted［28］.In order to facilitate the calculation，a polynomial expression is selected as

wherec=5% is the relative thickness，andxandyare dimensionless horizontal and ordinate coordinates of the blade profile line relative to the chord length，respectively.

Ref.［22］is used to perform triangulation on the inlet and outlet of the rotor.This provides the axial and circumferential velocities of airflow within the flow channel.Interpolation and calculation are then used to determine the temperature and pressure of the airflow at all stages.The calculation follows a discrete time and progressive space method，as depicted in Fig.4.

Fig.4 Calculation flow chart

2 Result Analysis

2.1 Comparison of two melting models

The “naked” ice particle model and the water-covered model are both simplified.It is necessary to compare the two models as they both exist in actual processes.

In Fig.5，a comparison is made between the melting processes of the “naked” ice particle model and the water-covered model under the same high-temperature air flow environment and with the same initial conditions of ice crystals.The specific running parameters are：constant airflow environment withUair=50 m/s，Tair=50 ℃，h=10%；Rice=20 μm，Tice=-20 ℃ for ice crystal initial conditions.The melting trends of the two models are similar，but the “naked” ice particle model takes longer to completely melt under the same working conditions.Specifically，the 20 μm particle takes 0.038 s to completely melt in the “naked” ice particle model.This is because，in the “naked” ice particle model，the radius of the ice crystal decreases over time，resulting in a reduction of the surrounding air convection heat transfer area.This，in turn，it leads to lower convective heat transfer during the later stages of melting compared to earlier stages.In contrast，when covered with a water film，the outside of the ice particle is wrapped with a water film，which keeps the heat exchange area and heat exchange rate essentially constant.Therefore，the melting time is shorter in the water-covered model.

Fig.5 Comparison of melting processes in the “naked” ice particle model and the water cover model

Furthermore，it can be observed that，during the initial melting stages before 0.023 s，the “naked” ice particle model actually melts slightly faster than the water-covered model.This can be attributed to the fact that，in the water-covered model，the presence of the water film means that some of the convective heat transfer must be used to evaporate the water film.This results in less heat being available for melting the ice crystals.

Therefore，in the early stages of melting，the presence of the water film slows down the melting rate due to the energy required for evaporation.However，in the middle and late stages of melting，the water film increases the convective heat transfer area such that it is much larger than that of the anhydrous film in the “naked” ice particle model.This leads to a more significant heat exchange，thereby accelerating the melting speed.Overall，the presence of the water film covering accelerates the melting of ice crystals.

2.2 Effect of initial conditions of ice crystals on motion heat transfer

Fig.6 shows the movement of ice crystals（equivalent volume diameter is 40 μm）of different shapes inside the compressor after being sucked in from the same position.The columns ice crystals collide in turn with the third and the fourth stage stator，with a liquid water content of 15.3% when colliding with the fourth stage stator.The hexagonal plate ice crystal collides with the third stage rotor and third stage stator with the liquid water content of 12.3%.

Fig.6 Trajectory of ice crystals with different shapes after inhalation

In Fig.7，the change in liquid water content along the axial position is shown after ice crystals of different shapes，but at the same initial position，are sucked into the compressor.The percentage content of liquid water mass of different shapes of ice crystals along the axial direction gradually increases.The liquid water content of spherical ice crystals is the largest at the same position，while the liquid water content of ellipsoidal and hexagonal ice plate crystals is similar.Specifically，at the outlet of the low-pressure compressor，the liquid water content of spherical，ellipsoidal，column-shaped，and hexagonal ice plate crystals is 49.05%，44.32%，29.74%，and 37.13%，respectively.

Fig.7 Changes in liquid water content after inhalation of ice crystals with different shapes

Fig.8 shows the speed variation of ice crystals along the axial positions after different shapes of ice crystals are sucked in from the same position.The gray dotted line is the axial velocity of air，and collisions between the first stage of rotor and stator are absent before 400 mm.The ice crystals continue to accelerate in the air flow，and their final speed approaches that of air.The subsequent velocity change follows the same trend as air.The occurrence of collision can also be found in the distribution of axial velocity，such as ellipsoidal ice crystals（red curve）violently colliding at the head of the third stage rotor.This results in greatly reduced axial velocity and an almost total loss of axial momentum due to a direct impact at the front end of the fifth stage rotor.

Fig.8 Changes in axial velocity of ice crystals with different shapes after suction

It is evident that there is a noticeable difference between spherical and non-spherical ice crystals in terms of their motion during the phase transition.Combining Figs.7 and 8，it can be found that when there is no collision effect，the speed of columnar particles increases the fastest，and that of spherical particles increases slower，because non-spherical particles have greater viscous force and greater acceleration in the air flow.Since spherical particles move more slowly in the air flow，the same axial displacement takes longer，which increases the heat and mass transferred between the hot air and the ice crystals，resulting in the maximum water content of spherical liquids and the lowest water content of columnar liquids.Due to the consistency of motion and melting laws，it can also be concluded that nonspherical characteristic indirectly affects the melting behavior through the particle trajectory and velocity in this study.

Fig.9 shows the variations in heat transfer coefficient for three ice crystal shapes in the engine runner.The heat transfer coefficient of the three ice crystal shapes are basically the same following the axial position trend with oscillatory fluctuations.The heat transfer coefficient of spherical ice crystals is large，and that of column ice crystals is small，attributed by the large difference between spherical ice crystals and air velocity.This difference not only affects the movement of ice crystals but also influences heat exchange.At an axial position of approximately 0.8 m，the ellipsoidal ice crystal had a violent collision，resulting in a significant difference between its axial velocity and air velocity，and a sudden increase in the heat transfer coefficient.

Fig.9 Changes in heat transfer coefficient after inhalation of ice crystals with different shapes

Fig.10 shows trajectories of ice crystals as they move within the compressor at various radii after being ingested.The larger the radius of the ice crystal，the stronger the inertia which results in upward drift of the trajectory.Ice crystals with a radius of 10 μm（yellow trajectory）collide with the third rotor and fourth stage stator where the liquid water content reaches 86.96%.As a result，the possibility of adhesion of 10 μm ice crystals on the fourth stage stator is greater.Ice crystals with a radius of 20 μm（purple trajectory）collide with the third stage stator，the fourth stage rotor and the fifth stage rotor.The liquid water content at the collision point is 45.53%，so the ice crystal with a radius of 20 μm is more likely to adhere to the third stator.Ice crystals with a radius of 100 μm（blue trace）collide with the first stage rotor，the first stator，and the third stage rotor，with the liquid water content still less than 10% at the outlet of the compressor.Therefore，ice crystals with a radius of 100 μm may not adhere nor freeze inside the low-pressure stage compressor.

Fig.10 Trajectory of ice crystals with different radii after inhalation

Fig.11 shows the axial temperature variation inside the compressor after inhalation of ice crystals with different radii and the same initial temperature.Ice crystals are warmed up from supercooled to 0 ℃along the axis before starting to melt.The smaller the ice crystals，the faster the temperature increases to reach 0 ℃.Ice crystals with a large size，such as those with a radius of 100 μm（purple curve），travel farther inside the compressor and enter the fifth stage of rotor and stator before they start melting.Due to their low liquid water content，the likelihood of adhesion and freezing is also low.

Fig.11 Temperature changes after inhalation of ice crystals with different radii

Figs.12 and 13 respectively show the melting rate and percentage variation of liquid water mass in the compressor after inhalation of ice crystals with different radii.The melting rates of ice crystals with different radii remain stable after melting starts.The liquid water content also increases slowly after melting starts.However，ice crystals with different radii begin to melt at different locations.The smaller the size，the earlier the crystals begin to melt.The melting rate for a larger size is faster，which is due to the larger surface area for heat transfer.

Fig.12 Melting rate of ice crystals with different radii changes after inhalation

The 80 μm particles exhibit a unique behavior in the front of the third stage stator（at approximately 0.6 m），which is due to their violent collision with the front of the blade.This collision results in a great decrease in the axial velocity of ice crystals，causing its speed to be greatly different from air flow speed，thereby enhancing heat exchange.Furthermore，due to the decreased velocity，the heat transfer time is longer under the same displacement，resulting in a sudden increase in the temperature，melting speed and liquid water content of the particle.Fig.13 shows that the liquid water content of 20，40，60，80，and 100 μm at the outlet of the low-pressure compressor is 60.4%，15.5%，6.3%，6.3%，and 1.9%，respectively.The liquid water content of the 60 μm ice crystal is less than 10% as it reaches the outlet.Therefore，ice crystals with size greater than 60 μm are difficult to adhere and freeze in the flow channel of the low-pressure compressor.

Fig.13 Percentage content of liquid water mass after inhalation of ice crystals with different radii changes

3 Conclusions

The study utilized a numerical model under the Lagrange framework to analyze the dynamic behavior and heat and mass transfer of ice crystals in the flow channel of an aeroengine after ingestion.The main conclusions are as follows：

（1）Two melting models of ice crystal particles were compared and analyzed.The results showed that the “naked” ice particle model took longer to fully melt compared to the water-covered model.The 20 μm particle took 0.04 s to melt，while the“naked” ice particle melted faster before 0.023 s.

（2）The liquid water content of spherical ice crystals was higher compared to other shapes under the same conditions.Specifically，at the outlet，the liquid water content of 20 μm spherical ice crystals was calculated to be 60.7%.On the other hand，the liquid water content of ellipsoidal and hexagonal plate ice crystals was close，with both shapes having a liquid water content around 40% at the outlet.

（3）The initial diameter affected motion and heat transfer of ice crystals in the compressor.The smaller the size，the earlier it started to melt，and the higher the liquid water content.The 20 μm crystal’s liquid water content at the outlet was as high as 60.7%，while the 40 μm one was only 17.8% under the same condition.The liquid water content at the outlet of the 60 μm ice crystal was close to 10%，enabling difficulty for the crystal to adhere to the ice.