Abstract

The ice crystal supercooled droplet mixed phase icing problem is an important research direction in aircraft icing and has received more attention in recent years. The thermodynamic process of the water film after the ice crystals impact on the surface determines the final ice shape, which is an important part of the accurate prediction of aircraft icing. In this paper, a thermodynamic model of ice crystal supercooled droplet mixed phase icing is proposed based on the extended Messinger model, according to the results of flow field and particle trajectory calculations. In this model, the mass and energy conservation equations of ice crystals, supercooled droplets, and liquid water are considered. The equations take the process of ice crystal adhesion and erosion into account, and the solution method of the equations is given. Ice shapes are calculated under various ice crystal supercooled droplet mixed phase conditions and compared with experimental results to demonstrate the validity of the numerical method. The effects of ice crystal erosion rate, melting ratio, and adhesion coefficient on the calculation results are analyzed by a numerical method. The results show that the ice crystal erosion rate has little effect on the ice shape, while a larger melting ratio and adhesion coefficient lead to more ice accretion.

1. Introduction

Icing phenomenon is a noteworthy problem in modern civil aviation. It is generally accepted that aircraft icing is caused by supercooled droplets impacting on the surface and freezing [1, 2]. Aircraft icing will destroy the original aerodynamic shape of the aircraft, resulting in a decrease in lift and an increase in drag, affecting the flight performance of the aircraft. At present, a lot of research has been carried out on the icing problem caused by supercooled droplets, and a large number of results have been achieved [3, 4]. In the amendment to FAR25 Appendix O [5], supercooled large droplets and ice crystals are included in the aircraft airworthiness certification. There is a growing concern about icing problem caused by ice crystals [6, 7].

There is a great difference between ice crystal icing and supercooled droplet icing. Firstly, ice crystals exist in high-altitude clouds above 6700 m, with a maximum water content of 9 g/m3 and a maximum diameter of 200 μm [8]. The water content and diameter of the ice crystals are much larger than that of conventional supercooled droplets, and the meteorological conditions for icing between the two are quite different. Secondly, when the ice crystals impact on the dry surface of the aircraft, the ice crystals will not adhere to the surface and will not freeze due to the low surface temperature. However, when the ice crystals impact on the wet surface of the aircraft, they will transfer heat and mass with the water film on the surface, and thus, icing will occur. At the same time, rebound, adhesion, and erosion will occur during the ice crystal impact, and the adhesion process is related to the content of ice crystal and the thickness of the water film. Therefore, the mechanism of ice crystal icing is quite different from that of supercooled droplets [911]. Finally, the icing range of supercooled droplet is determined by surface water film flow processes, whereas ice crystal icing is more complex. Ice crystals will increase the mass flux of the surface water film during the icing process, leading to an increase in the range of water film flow. The increase of water film flow range will lead to more ice crystals impacting on the surface to participate in ice accretion. Ice crystal icing is therefore more extensive and hazardous than supercooled droplet icing. It is the fact that ice crystal icing is more severe and hazardous than supercooled droplets, which requires special attention.

At present, the icing problem caused by ice crystals has been studied by many researchers and numerous results have been obtained. Khalil et al. [12] carried out wind tunnel experiments on the performance of thermal anti-icing systems under mixed phase icing meteorological conditions, and significant erosion phenomenon was observed. Struk et al. [13] carried out rational research on ice crystal icing of many small samples using ice crystal experimental device, and icing and shedding processes were observed. Bartkus et al. [14] analyzed the mixed phase icing experiment of NACA 0012 airfoil and proposed two different types of ice accretion based on different surface energy balance. Baumert et al. [15] studied the mixed phase icing of NACA 0012 airfoil and cylinder through experiments. By adjusting the content of supercooled droplets and ice crystals, the icing conditions of different mixed phases were obtained and it was found that the lower the temperature, the lower the icing amount. In a numerical simulation, Wright et al. [16] added the mass and energy terms related to ice crystals to the thermodynamic model of supercooled droplet icing and obtained the thermodynamic model of mixed phase icing. Nilamdeen and Habashi [17, 18] developed a mixed phase icing model considering air, supercooled droplets, and ice crystals to simulate the icing growth. In the absence of experimental results, the effectiveness of the model was verified by using NACA0012 airfoil and cylinder. Trontin et al. [19] proposed a semiempirical model considering the effect of ice crystals on the mass and momentum conservation. By adjusting the parameters in the model, the predicted results of the model were consistent with the experimental data. Norde et al. [20] used the Eulerian method to calculate the trajectory of ice crystals and calculated the ice accretion of ice crystals after impact. Zhang et al. [21] studied the mixed phase icing of ice crystal based on a secondary development of Fluent. Rios and Cho [22] divided the ice accretion of ice crystals into two stages: the first stage is that the ice crystals form a water film on the blade surface, and the second stage is that the subsequent ice crystals are captured by the water film, and then, ice accretion was calculated using a modified Messinger thermodynamic model. Tsao et al. [23] developed a thermodynamic model to describe the icing of ice crystals on the unheated surface within the compression system of a turbofan engine. Through the analysis of the experimental data, it is believed that the model was able to capture some important qualitative trends of the ice crystal icing process in the target region of low-pressure compressor under different icing conditions.

In this paper, the numerical simulation method of the icing process under the mixed phase of ice crystal and supercooled droplet is studied. The physical phenomena in the process of ice crystal impacting on the surface is analyzed, and the physical model of ice crystal impacting is developed. Based on the improved Messinger model, the heat and mass transfer process between the ice crystal and water film is analyzed, and the phase transition model of icing after the ice crystal impact is built to simulate the ice shapes of mixed phase. The numerical simulation results are compared with the experimental results to verify the correctness of the method. Finally, the effect of different model parameters on the calculation results is investigated.

2. Numerical Simulation Methods

The numerical simulation of mixed phase icing of ice crystals and supercooled droplets mainly includes many modules such as mesh generation, flow field calculation, ice crystal/supercooled droplet trajectory calculation, and ice accretion model [24, 25]. In this paper, a two-dimensional structural grid is used to mesh around the airfoil. The flow field data is obtained by solving the N-S equation, and the standard SST model is selected as the turbulence model. Ice crystal/supercooled droplet trajectory calculation uses Lagrangian method to calculate single particle trajectory, and the surface collection coefficient is obtained by calculating a large number of particle trajectories in space. The collection coefficient and air flow field data are used as the input conditions of the ice accretion model. Mesh generation, flow field calculation, and ice crystal/supercooled droplet trajectory calculation methods have been relatively mature, and this paper describes the ice accretion model in detail.

2.1. Mixed Phase Icing Model

In order to simulate icing, it is necessary to establish control units on the surface, and subsequently, the mathematical model is constructed. The control units for icing calculations are distributed mainly on the icing surface. The upper surface of the control unit is in contact with the air, and the lower surface is a solid surface, which can be an impingement surface or an ice layer. The ice accretion model of mixed phase is extended based on the Messinger theory. Ignoring the heat conduction inside the ice layer and between the ice layer and the surface, the mass and energy conservation equations of the control unit are analyzed by using the control volume method.

2.1.1. Mass Conservation Equation

The mass conservation equation within the control unit can be established as follows: where is the mass flux of the liquid water flowing from the upstream control unit into the current control unit, is the mass flux of the supercooled droplets impacting on the surface, is the mass flux of the melted parts of the ice crystals impacting on the surface, is the mass flux of the unmelted parts of the ice crystals impacting on the surface, is the mass flux of water evaporated within the control unit, is the mass flux of liquid water flowing from the current control unit into the downstream control unit, and is the mass flux of liquid water frozen into ice.

In the mass equation, is the mass flux of the supercooled droplets impacting on the surface, that is, the mass flux of liquid water collected by the surface, which can be expressed by the following equation after the flow field and droplet trajectories have been calculated: where is the liquid water content, is the far-field inflow velocity, and is the droplet collection coefficient.

The mass flux of the melted parts of the ice crystals impacting on the surface is as follows:

The mass flux of the unmelted parts of the ice crystals impacting on the surface is as follows: where is the ice crystal content, is the ice crystal collection coefficient, and is the ice crystal melting ratio.

The mass flux of water evaporated within the control unit can be obtained by the heat and mass transfer analogy as follows: where is the convective heat transfer coefficient, which can be obtained by boundary layer integration method; is the specific heat capacity of air at constant pressure, which can be taken as ; is the air density; is the water vapor gas constant, which can be taken as ; is the Lewis number of the heat and mass transfer analogy criterion, whose value is related to the water content in the air, and can be taken as 1 under icing conditions; and and are the saturated water vapor pressure corresponding to the surface temperature of the control unit and incoming flow temperature , respectively, which can be expressed as follows:

when , When ,

is the phase change equilibrium temperature of water, which is 273.15 K.

As a result, the mass flux of flowing from the current control unit into the downstream control unit can also be obtained:

Under the assumptions of the Messinger model, the mass flux of water flowing into a particular control unit is equal to the mass flux of water flowing out of its upstream control unit.

2.1.2. Energy Conservation Equation

After defining the meaning of each item in the mass equation, the expression for the energy equation can be obtained. The energy conservation equation in the control unit is established as follows: where is the heat flux of the liquid water flowing from the upstream control unit into the current control unit, is the heat flux of the supercooled droplets impacting on the surface, is the heat flux of the melted parts of the ice crystals impacting on the surface, is the heat flux of the unmelted parts of the ice crystals impacting on the surface, is the heat flux conducted by the surface, is the heat flux of evaporation, is the heat flux of the liquid water flowing out of the current control unit, is the heat flux of liquid water frozen into ice, and is the heat flux of convective heat transfer.

As the temperature of water remains constant during the freezing phase change, the heat flux in the energy equation are expressed in terms of enthalpy in order to better simulate the icing process.

The heat flux of the supercooled droplets collected within the control unit consists of the internal and kinetic energy of the supercooled droplets.

Similarly, the heat flux of the ice crystals collected within the control unit includes both the internal and kinetic energy of the ice crystals. where is the specific heat capacity of ice at constant pressure, which can be taken as ; is the latent heat of ice melting, which can be taken as ; is the specific heat capacity of water at constant pressure, which can be taken as ; and is the velocity of supercooled droplets or ice crystals impacting on the surface, which can generally be approximated as the inflow velocity, namely, .

The heat flux taken away by evaporation can be described by the following expression: where is the latent heat of water evaporation, which can be taken as in the calculation.

The heat flux of the liquid water flowing from upstream control unit can be expressed as

The heat flux of the liquid water flowing from the control unit to the downstream control can be expressed as

In the icing calculation, the anti-icing and deicing system is generally not considered, and the surface is considered to be adiabatic, so the heat flux conducted by the surface is equal to 0:

The heat flux of convection heat transfer can be obtained from the convection heat transfer calculation equation: where is the convective heat transfer coefficient, which can be solved by the boundary layer integral method in icing calculation, and is the boundary layer recovery temperature, which can be expressed as follows: where is the recovery coefficient and is the airflow velocity at the outer boundary of the boundary layer.

The heat flux of liquid water frozen into ice is

2.2. Adhesion Effect of Ice Crystal

In the study of supercooled droplet icing, it is generally assumed that droplets fully adhere to the surface after impacting on the surface and participate in the process of icing phase change, but the physical process of ice crystal impingement is quite different from that of supercooled droplets. The adhesion of the ice crystal impingement varies depending on the surface condition. When there is no liquid water on the surface, ice crystals will rebound after hitting the surface, and no ice crystals will adhere to the surface. According to ice crystal experiments by Baumert et al. [15], only a small part of the rebound ice crystals will impact on the surface again, which has little effect on the whole icing process. Therefore, the rebound process can be considered that the ice crystal does not adhere to the surface and does not affect the movement of other ice crystals. When the ice crystals impact on the surface of water film, adhesion phenomenon will occur, especially under the mixed phase of ice crystals and supercooled droplet conditions. Trotin and Villedieu [26] analyzed the effect of liquid water content in the mixed phase on the adhesion of ice crystal impingement based on the experiment and given the corresponding empirical equations. According to this equation, the adhesion coefficient of the ice crystal under mixed phase conditions can be obtained. By introducing the adhesion coefficient that only takes into account the melting of ice crystals without supercooled droplets and the adhesion coefficient that takes into account the combined effect of supercooled droplets and ice crystals , the final adhesion coefficient used is the larger of the two: where where is the melting ratio of ice crystal, which is 0 when the ice crystal does not melt completely. According to the experiment under the ice crystal condition, it can be taken as . is the ratio of the mass of supercooled droplets impacting on the surface to the total mass of particles impacting on the surface. is the ratio of the mass of ice crystals impacting on the surface to the total mass of particles impacting on the surface. It can be found that the sum of and is 1, which can be expressed as follows:

and are constants. Considering that the lower the surface temperature is, the less ice crystals will adhere, especially on the dry surface with lower temperature, the ice crystals will not adhere. Thus, the constants in the adhesion model are related to the surface temperature.

From the above definition,the ratio of all the water in liquid form impacting on the surface to the total mass of the impingement particles can be found. Thus, if the ice crystal in the mixed phase does not melt, ; then, .

When ice crystals and supercooled droplets impact on the surface under mixed phase conditions, due to the low ambient temperature, it can be assumed that the ice crystals do not melt during the movement process, and the adhesion coefficient is simplified to . Considering the ice crystal adhesion effect, the mass flux of ice crystals that actually adhere to the icing surface and participate in the icing process is as follows:

2.3. Erosion Effect of Ice Crystal

Another major difference between the process of ice crystal and supercooled droplet impacting on the surface is that the erosion phenomenon will occur when the ice crystal impact on the surface of the ice layer. As ice crystals are brittle particles, part of the ice flies out under the impingement of ice crystals, resulting in the reduction of ice layer mass. The icing wind tunnel experiments have found that the erosion phenomenon has a strong influence on the icing process of ice crystals, which will lead to the reduction of ice amount. At present, the research on ice crystal erosion usually uses the erosion rate to characterize the effect of erosion on ice accretion. It can be found that the erosion rate is related to many parameters such as the velocity, angle, and mass of ice crystal impingement, and the relevant empirical relationship can be established according to the experimental data. Trontin and Villedieu [26] developed an erosion model based on the results of wind tunnel experiments with icing and obtained an empirical model of erosion rates in the following equation: where is the mass fraction of liquid water in the ice layer. , , , and .

Based on the existing wind tunnel experimental data, the empirical coefficient is obtained, and the mass fraction of liquid water in the ice layer is considered. The higher the liquid water content is, the more obvious the ice crystal erosion effect is. Empirical coefficients were obtained based on available wind tunnel experimental data, and the mass fraction of liquid water in the ice layer was taken into account. It was concluded that the higher the liquid water content is, the more obvious the ice crystal erosion effect is. The higher the content of solid ice crystals, the weaker the erosion effect. The erosion effect is infinitely greater when the mass fraction of liquid water reaches 0.6. At the same time, the model considers that the erosion effect is related to the tangential velocity of ice crystal impact, but this will lead to the abnormal calculation results near the stagnation point. Therefore, the model considers that the erosion effect increases with the increase of surface curvature, and the erosion rate is corrected by curvature.

According to the assumption of empirical equation, the mass flux of erosion is determined by the erosion rate and the impingement mass of ice crystals:

It can be seen that the erosion must be less than the icing on the surface.

2.4. Solution Method of Model

For a two-dimensional airfoil, the mass flux of liquid water flowing into the control unit at the stagnation point is considered to be 0. The liquid water within the control unit at the stagnation point is divided equally into two parts, which overflow to the upper and lower surfaces, respectively. There is a constraint relationship between the water flowing into and out of the control unit, that is, the mass flux of the water flowing out of the upstream control unit is equal to the mass flux of the water flowing into the downstream control body.

The specific method for solving the ice accretion model of mixed phase icing is as follows: the icing surface is divided into upper surface and lower surface from the stagnation point. First, the control body at the stationary point is calculated. Then, the downstream control units are solved sequentially for the upper and lower surfaces, respectively, until the surface water overflow limit is reached, and the solution of all control units is completed.

For each control unit, the ice accretion model is solved by an iterative process. It is assumed that the surface temperature is equal to the melting temperature . According to the energy conservation equation, the mass flux of frozen ice is obtained as follows:

The results within the control unit can be determined from the mass flux of frozen ice as follows: (1)If , it means that the mass flux of frozen ice is greater than the total mass flux of liquid water in the control unit, and all liquid water within the control unit is frozen into ice. Then, the assumption of is not tenable. At this time, , icing state is glazed ice. Set , can be obtained by the iterative calculation of the energy conservation equation(2)If , it means that the solid ice crystal impingement within the control unit are not completely melted, and the liquid water within the control unit is not completely frozen. Ice and water exist at the same time. The assumption of is valid, and the icing state is rime ice(3)If , it means that the solid ice crystal impingement within the control unit is melted, and there is no ice in the control unit; only liquid water is present. Then, the assumption of is not tenable. At this time, , the surface is in a wet state. Set , can be obtained by iterative calculation of the energy conservation equation

Through the above process, the ice accretion model is solved, and the surface temperature and the mass flux of frozen ice of each control unit on the surface are obtained. After solving the ice accretion model, the freezing rate can be defined as follows:

The mass flux of frozen ice is first solved from the equations of mass and energy conservation. The icing growth rate can be obtained by subtracting the mass flux of frozen ice of each control unit from the mass flux loss caused by erosion [26]. The ice shape can be obtained by converting icing growth rate into ice thickness and increasing it along the surface normal direction. The ice thickness calculation equation is as follows: where is the icing time and is the ice density.

3. Method Validation

In order to verify the correctness of the method in this paper, several calculation conditions are used [27]. As pure ice crystals do not lead to icing, the ice crystal supercooled droplet mixed phase icing conditions are used for calculation. The calculation model is NACA0012 with the chord length of 0.9144 m, so the same grid is used in the calculation. As shown in Figure 1(a), the two-dimensional structured grid is used to mesh the flow field grid around the airfoil. The far field distance is set at 10 m to ensure that the flow field is not affected by the boundary condition. The total number of grid nodes is 22000, of which the number of grid nodes on the airfoil surface is 200. The grid independence calculations are performed by selecting the grid heights of 0.05 mm, 0.01 mm, 0.005 mm, and 0.001 mm for the first layer, respectively. It is confirmed that when the grid height of the first layer is 0.005 mm, and the flow field calculation results changed slightly, and the grid meet the calculation requirements. In the flow field calculation, the results shown in Figure 1(b) are obtained by solving the N-S equations, and the standard SST model is selected as the turbulence model. Figure 1(c) shows the motion trajectory of 150 μm ice crystals calculated using the Lagrangian method [28], and the collection coefficients are shown in Figure 1(d). The calculated icing results will be compared with the experimental results and the calculation results of existing method in the reference [19].

3.1. Condition (1)

The inflow Mach is 0.163, the ambient pressure is 98 kPa, and the angle of attack is 0°. The droplet diameter is 20 μm,  g/m3, the ambient temperature is 266.18 K, and the icing time is 600 s. The calculated ice shape is compared with the experimental result and reference result as shown in Figure 2.

3.2. Condition (2)

The inflow Mach is 0.163, the ambient pressure is 98 kPa, and the angle of attack is 0°. The droplet diameter is 20 μm,  g/m3, the ice crystal diameter is 200 μm,  g/m3, the ambient temperature is 266.18 K, and the icing time is 600 s. The calculated ice shape is compared with the experimental result and reference result as shown in Figure 3.

3.3. Condition (3)

The inflow Mach is 0.165, the ambient pressure is 98 kPa, and the angle of attack is 0°. The droplet diameter is 20 μm,  g/m3, the ice crystal diameter is 150 μm,  g/m3, the ambient temperature is 260.65 K, and the icing time is 600 s. The calculated ice shape is compared with the experimental result and reference result as shown in Figure 4.

3.4. Condition (4)

The inflow Mach is 0.165, the ambient pressure is 98 kPa, and the angle of attack is 0°. The droplet diameter is 20 μm,  g/m3, the ice crystal diameter is 150 μm,  g/m3, the ambient temperature is 260.65 K, and the icing time is 600 s. The calculated ice shape is compared with the experimental result and reference result as shown in Figure 5.

Through the above comparison between the calculation results and the experimental results under the different conditions [27], it can be seen that the calculation results of this paper are in good agreement with experimental results, which show the feasibility of the calculation model. Compared with condition (1), condition (2) increased the ice crystal content. It can be seen that the experimental results of condition (2) decreased compared with condition (1), indicating that there is an erosion effect in addition to the adhesion effect of the ice crystals. Conditions (3) and (4) are glazed ice conditions and the calculations agree well with the experimental results. Compared with the ice accretion model in reference [19], it can be found that the model in this paper can effectively simulate ice shapes under different conditions, and the calculation results are closer to the experimental values under some conditions. From the perspective of engineering applications, the predicted ice shape results in this paper are more conservative and more meaningful for aircraft icing and anti-icing design.

It is worth noting that the calculation results for other conditions are very similar, but the ice shapes for condition (2) are larger than the experimental values. This is due to the higher ice crystal content and the higher icing temperature in condition (2), which makes the erosion effect more obvious. The erosion rate of 0.3 selected from the reference is not suitable for this condition, which indicates that it is necessary to study the effect of various model parameters on the icing characteristics.

4. The Effect of Different Model Parameters on the Icing Characteristics of Ice Crystal

In order to investigate the effect of different model parameters on the calculated results, the effects of erosion rate, melting ratio, and adhesion coefficient on the calculated ice shape are compared, respectively.

4.1. The Effect of Ice Crystal Erosion Rate on Icing Characteristics

The erosion effect is mainly reflected in the erosion rate in the calculation model, and the magnitude of the erosion rate is expressed as an effect on the ice crystal to icing rate. Therefore, this paper mainly studies the effect of erosion rate on ice shapes.

Validation condition (2) is chosen for the calculation model and conditions. According to the basic erosion rate in the reference being 0.3, the ice shapes for erosion rates of 0, 0.3, 0.5, 0.8, 1, and 1.2 are studied in this paper, respectively.

The ice shapes under different erosion rates are shown in Figure 6. As can be seen from the figure, with the increase of the erosion rate, the overall amount of ice gradually decreases. This is due to the enhancement of erosion, and the erosion effect on the icing is gradually reduced. Different erosion rates have a certain effect on the calculation results. The erosion effect will lead to the removal of a portion of the ice layer during the ice crystal impact process. The stronger the erosion effect, the smaller the amount of ice accretion. For condition (2), comparing the calculation results under different erosion rates, it can be seen that the calculation results obtained with larger erosion rates are closer to the experimental values.

The erosion mainly occurs on the upper and lower sides, and erosion near the stagnation point is relatively small. This is because according to the erosion calculation equation, erosion is related to the tangential velocity of the impingement, while the tangential velocity near the stagnation point is relatively small; hence, the erosion amount is relatively small.

4.2. The Effect of Ice Crystal Melting Ratio on Icing Characteristics

The ice crystal melting ratio is the percentage of ice crystals that have become liquid during the ice crystal movement. Different melting ratios of ice crystal will lead to changes in the physical process of ice crystal icing, which will have a great effect on the final ice shape. Therefore, the effect of melting ratios on the ice shape is mainly studied.

Validation condition (2) is chosen for the calculation model and conditions. Normally, the ice crystal melting ratio is 0 because the temperature of the ice crystals is the same as the ambient temperature and lower than 0 degrees centigrade. In certain conditions, especially when ice crystals enter the engine inlet, melting may occur. The ice shapes for melting ratios of 0, 0.1, 0.2, 0.5, and 0.8 are investigated, respectively, and the results are shown in Figure 7.

It can be seen from the figure that with the increase of the melting ratio, the overall amount of icing gradually increases and the icing range gradually expands. As the melting ratio increases, the proportion of liquid water in ice crystals gradually increases. More ice crystals will be involved in surface heat and mass transfer after impacting on the surface, which is equivalent to an increase in the total amount of liquid water collected on the surface. The liquid water on the surface flows to the upper and lower sides and gradually freezes under the influence of the ambient environment. Therefore, with the increase of melting ratio, the amount and range of ice accretion will increase.

When the melting ratio exceeds 0.5, the range of ice accretion will increase more and the amount of ice accretion will also be greater. This is because the impact range of ice crystals is larger than the impact range of droplets. When the melting ratio is small, only a small portion of the melted ice crystals collected outside the impact range of droplets. Although this portion of the melted ice crystals can form ice, under the impact of a large number of unmelted ice crystals, these small amounts of ice will fly away due to the erosion effect, making it difficult to generate an ice layer. When the melting ratio is relatively large, more melted ice crystals participate in the icing process, while the erosion effect of unmelted ice crystals decreases. More molten ice crystals enter the mixed phase icing phase transition, increasing the amount and flow range of liquid water, resulting in a significant increase in icing amount and icing range of ice crystals.

4.3. The Effect of Ice Crystal Adhesion Coefficient on Icing Characteristics

The ice crystal adhesion coefficient refers to the percentage of ice crystals that will adhere to the surface during the impingement of ice crystals. The difference of the adhesion coefficient will affect the amount of collected ice crystals, which in turn affect the icing process of ice crystals, and ultimately have a greater impact on the ice shape. Thus, the effect of the adhesion coefficient on the ice shape is mainly investigated.

Validation condition (2) is chosen for the calculation model and conditions. The adhesion coefficient can be calculated to be around 0.2 according to the equations. In this paper, the adhesion coefficients are set to 0.2, 0.5, 0.8, and 1, respectively, and the ice shapes are calculated as shown in Figure 8.

As can be seen from the figure, with the increase of the adhesion coefficient, the overall amount of icing gradually increases, but the icing range does not expand significantly. This is due to the fact that as the adhesion coefficient increases, it indicates that ice crystal impinging on the wet surface are more likely to adhere to the surface rather than splash out. The amount of ice crystals adhering to the surface gradually increases, and more ice crystals are involved in the heat and mass transfer on the surface. Under the action of the external environment, ice crystals and droplet icing, resulting in an increase of ice accretion. However, outside the droplet impingement area, due to the absence of liquid water, ice crystals do not adhere to the dry surface after impact. The adhesion coefficient is meaningless for dry surfaces, and no ice crystals participate in the icing process. Therefore, the increase in the adhesion coefficient will not lead to ice accretion outside the icing range of droplets.

5. Conclusion

The numerical simulation method of icing under the mixed phase of ice crystal supercooled droplet is studied, and a thermodynamic model of mixed phase icing is proposed. Based on the improved Messinger model, the mass and energy conservation equations of ice crystals and supercooled droplets impacting on the surface are established. The erosion and adhesion phenomena during ice crystal impingement are investigated in the model, and the solution method of the control equations is given. In order to verify the numerical calculation method in this paper, the ice shapes of several mixed phase conditions are calculated, and the results are compared with the experimental results, which prove that the model proposed in this paper is effective. The ice shapes of the mixed phase condition are investigated under different model parameters, and the effects of erosion rate, melting ratio, and adhesion coefficient on the ice shape are mainly studied. Based on the numerical simulation method in this paper, the following conclusions are obtained: (1)The ice accretion model of ice crystal supercooled droplet mix phase proposed in this paper can simulate the ice shape under different conditions, and the calculated results are in agreement with the experimental results(2)Under the mixed phase condition of ice crystal supercooled droplet, the erosion rate has a relatively small effect on the ice shape. The larger the erosion rate is, the smaller the ice shape is(3)The melting ratio in the process of ice crystal movement has a great impact on the ice shape. The higher the melting rate, the more liquid water film on the surface, and the larger the amount and range of ice accretion(4)The coefficient of adhesion of ice crystals has a certain influence on icing; the larger the coefficient of adhesion, the greater the amount of icing, but the range of icing is almost constant(5)The ice crystal adhesion coefficient has a certain effect on ice accretion. The larger the adhesion coefficient, the greater the amount of icing, but the icing range is almost unchanged

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 51806105), the Open Fund of the Key Laboratory of Icing and Anti/De-icing (Grant No. IADL20190309), the Postgraduate Research & Practice Innovation Program of Jiangsu Province, and the Scientific Research Foundation of Nanjing Institute of Technology (Grant No. YKJ202009).