Investigation on the Mechanism of Heat Load Reduction for the Thermal Anti-Icing System

The aircraft ice protection system that can guarantee flight safety consumes a part of the energy of the aircraft, which is necessary to be optimized. A study for the mechanism of the heat load reduction in the thermal anti-icing system under the evaporative mode was presented. Based on the relationship between the anti-icing heat load and the heating power distribution, an optimization method involved in the genetic algorithm was adopted to optimize the anti-icing heat load and obtain the optimal heating power distribution. An experiment carried out in an icing wind tunnel was conducted to validate the optimized results. The mechanism of the anti-icing heat load reduction was revealed by analyzing the influences of the key factors, such as the heating range, the surface temperature and the convective heat transfer coefficient. The results show that the reduction in the anti-icing heat load is actually the decrease in the convective heat load. In the evaporative mode, decreasing the heating range outside the water droplet impinging limit can reduce the convective heat load. Evaporating the runback water in the high-temperature region can lead to the less convective heat load. For the airfoil, the heating power distribution that has an opposite trend with the convective heat transfer coefficient can reduce the convective heat load. Thus, the optimal heating power distribution has such a trend that is low at the leading edge, high at the water droplet impinging limit and zero at the end of the protected area.


Introduction
Ice accretion always occurs on aircraft wings and stabilizers when the aircraft is flying through super-cooled water droplet clouds. This phenomenon may have adverse impacts on flight safety where aerodynamic and maneuvering performances are threatened if no adequate protection is supplied [1]. To prevent the aircraft wings and stabilizers from ice accretion, all commercial aircrafts are equipped with ice protection systems [2]. Generally, the thermal anti-icing system is the most widely used for large commercial aircraft and it also accounts for about 20% energy consumption of the aircraft engines [3][4][5]. Therefore, it is critical to optimize the thermal anti-icing system to achieve energy conservation, which can lighten the take-off weight and extend the flying range [6,7].
To prevent ice accretion, the aircraft is equipped with many kinds of anti-icing systems, which can be classified as passive and active anti-icing systems. The passive systems utilize hydrophobic or icephobic materials [8,9] to coat the surfaces of the protected components. The hydrophobic or icephobic properties can keep the surfaces as clean as possible so that ice cannot easily accumulate [10][11][12]. However, in many cases, the passive systems alone are difficult to keep the surfaces free from ice revealed. Moreover, there is a lack of experimental investigations on the mechanism of the anti-icing heat load reduction. The influence of key parameters on the optimal heat flux distribution has not been fully studied. Therefore, in this paper, the mechanism of anti-icing heat load reduction was uncovered through both numerical simulation and experimental validations. The influences of the heating range, the surface temperature and the convective heat transfer coefficient on the anti-icing heat load reduction were also investigated. The outline of this paper is as follows. First of all, Section 2 presents the modeling and optimization methods of the thermal anti-icing system. Then, to demonstrate the applicability of the numerical models and the optimization methods, an experimental test is conducted and introduced in Section 3. Section 4 gives the numerical and experimental results and discusses the underlying mechanism of the anti-icing heat load reduction. The conclusion is presented in Section 5.

Modeling of the Anti-Icing Process
According to the Messinger model [32], the mass and energy conservation equations of a control volume for the runback water on the anti-icing surface can be formulated from Figure 1: Q w,in + (q w + q a + q n )dx = q evap + q conv dx + Q w,out (2) where . m in is the mass coming in from the previous control volume, . m out is the mass flowing into the next control volume, . m imp is the mass of water droplets impingement per unit length, . m evap is the mass of evaporation per unit length, q w is the heat flux converted from the kinetic energy of the water droplets impingement, q a is the aerodynamic heating flux, q n is the heat flux supplied by the thermal anti-icing system, q evap is the evaporative heat flux, q conv is the convective heat flux, Q w,in is the energy coming in from the previous control volume, and Q w,out is the energy flowing into the next control volume. The above mass terms are detailed as follows:  (5) where β is the local water droplet collection efficiency, u ∞ is the freestream air velocity, LWC is the liquid water content, α is the local convective heat transfer coefficient, c p,air is the specific heat of air, Pr is the Prandtl number, Sc is the Schmidt number, M water is the relative molecular mass of water, M air is the relative molecular mass of air and p e is the static pressure at the edge of the boundary layer. p v,s and p v,e are the saturation vapor pressure of water on the airfoil surface and at the edge of the boundary layer, respectively. The above energy terms are detailed as follows: m in · (T s,in − T r ) (6) Q w,out = c p,w · . m out · (T s − T r ) (7) q w = . m imp · u 2 ∞ 2 + c p,w · (T ∞ − T r ) (8) q a = α · γ · u 2 e 2c p,air (9) q evap = . m evap · r + . m evap · c p,w · (T s − T e ) (10) Energies 2020, 13, 5911 4 of 19 q conv = α · (T s − T e ) (11) T e = T ∞ + u 2 ∞ − u 2 e 2c p,air (12) where c p,w is the specific heat of water, T s,in is the temperature of the previous control volume, T s is the temperature of the current control volume, T r is the reference temperature (273.15K), T e is the temperature at the edge of the boundary layer, T ∞ is the freestream air temperature, γ is the recovery factor, u e is the velocity at the edge of the boundary layer and r is the latent heat of water evaporation.  (12) where c p,w is the specific heat of water, T s,in is the temperature of the previous control volume, T s is the temperature of the current control volume, T r is the reference temperature (273.15K), T e is the temperature at the edge of the boundary layer, T ∞ is the freestream air temperature, γ is the recovery factor, u e is the velocity at the edge of the boundary layer and r is the latent heat of water evaporation. Some parameters involved in the calculation of the anti-icing heat load, such as α, β and p e , are determined by the solution of flow field, which can be divided into two steps. The first step is to solve the air flow field by using ANSYS FLUENT CFD code. The Spalart-Allmaras model was selected as the turbulence model. The second step is to calculate the water droplet impingement characteristics by means of the Eulerian approach [33]. A UDF (User Defined Function) [34] was built to solely calculate the water droplets' movement based on the simulation results of the air flow field. To calculate the anti-icing heat load precisely, the heat conduction calculation was also taken into consideration [35,36].

Optimization Method of the Anti-Icing Heat Load
Based on the numerical simulation of the anti-icing process, an optimization method was proposed to optimize the heat flux distribution. For the convenience of practical application, the entire protected area was divided into multiple heating zones, as illustrated in Figure 2 whose heating power can be adjusted independently. The heating power of each zone can be treated as an individual variable and all these independent variables ultimately constituted the heating power distributions. By organizing a set of equations described in Section 2.1, the function for calculating the anti-icing heat load was established. The input of this function was the heating power distribution, and the output was the anti-icing heat load. The optimization goal is to find an optimal heating power Some parameters involved in the calculation of the anti-icing heat load, such as α, β and p e , are determined by the solution of flow field, which can be divided into two steps. The first step is to solve the air flow field by using ANSYS FLUENT CFD code. The Spalart-Allmaras model was selected as the turbulence model. The second step is to calculate the water droplet impingement characteristics by means of the Eulerian approach [33]. A UDF (User Defined Function) [34] was built to solely calculate the water droplets' movement based on the simulation results of the air flow field. To calculate the anti-icing heat load precisely, the heat conduction calculation was also taken into consideration [35,36].

Optimization Method of the Anti-Icing Heat Load
Based on the numerical simulation of the anti-icing process, an optimization method was proposed to optimize the heat flux distribution. For the convenience of practical application, the entire protected area was divided into multiple heating zones, as illustrated in Figure 2 whose heating power can be adjusted independently. The heating power of each zone can be treated as an individual variable and all these independent variables ultimately constituted the heating power distributions. By organizing a set of equations described in Section 2.1, the function for calculating the anti-icing heat load was established. The input of this function was the heating power distribution, and the output was the anti-icing heat load. The optimization goal is to find an optimal heating power distribution that can reduce the anti-icing heat load in the evaporative mode and the constraint condition can be expressed as follows: subject to . m out,end = 0 275K ≤ T i ≤ T max (15) where the subscript end denotes the end of the protected area. Equation (14) means that the runback water does not flow out of the protected area, ensuring the system operates in the evaporation mode. Equation (15) makes the protected area free of any ice. Moreover, the mechanical properties of the aircraft structure are greatly affected by temperature. The ultimate strength, the yield strength and the fatigue life of many composite materials decrease as temperature increases [37]. Furthermore, Energies 2020, 13, 5911 5 of 19 large temperature differences may cause non-uniform expansions in different parts of the structure, resulting in thermal stresses, added to the other imposed stresses [38]. As a consequence, it is recommended that the maximum temperature should be limited to the temperature at which the material is allowed to stand.
i m a x 275K T T   (15) where the subscript end denotes the end of the protected area. Equation (14) means that the runback water does not flow out of the protected area, ensuring the system operates in the evaporation mode. Equation (15) makes the protected area free of any ice. Moreover, the mechanical properties of the aircraft structure are greatly affected by temperature. The ultimate strength, the yield strength and the fatigue life of many composite materials decrease as temperature increases [37]. Furthermore, large temperature differences may cause non-uniform expansions in different parts of the structure, resulting in thermal stresses, added to the other imposed stresses [38]. As a consequence, it is recommended that the maximum temperature should be limited to the temperature at which the material is allowed to stand.  For the optimization problem, it is quite suitable to select the Genetic Algorithm (GA) [39]. The Genetic algorithm can deal with multiple individuals simultaneously in a population, which means For the optimization problem, it is quite suitable to select the Genetic Algorithm (GA) [39]. The Genetic algorithm can deal with multiple individuals simultaneously in a population, which means that it evaluates multiple solutions in the searching space, reducing the risk of falling into the local optimal solution. Moreover, the fitness function is not constrained by continuous differentiability and its domain can be set arbitrarily. In the present study, a heating power distribution under the evaporative mode was encoded into the initial population. By following the selection, mutation and intersection, the initial population was evolved into the high fitness population through generations. Subsequently, the optimal heating power distribution that can reduce the anti-icing heat load was obtained. The flow chart of the optimization method is illustrated in Figure 3.
distribution that can reduce the anti-icing heat load in the evaporative mode and the constraint condition can be expressed as follows: (13) subject to out,end 0 m =  (14) i m a x 275K T T ≤ ≤ (15) where the subscript end denotes the end of the protected area. Equation (14) means that the runback water does not flow out of the protected area, ensuring the system operates in the evaporation mode. Equation (15) makes the protected area free of any ice. Moreover, the mechanical properties of the aircraft structure are greatly affected by temperature. The ultimate strength, the yield strength and the fatigue life of many composite materials decrease as temperature increases [37]. Furthermore, large temperature differences may cause non-uniform expansions in different parts of the structure, resulting in thermal stresses, added to the other imposed stresses [38]. As a consequence, it is recommended that the maximum temperature should be limited to the temperature at which the material is allowed to stand.  For the optimization problem, it is quite suitable to select the Genetic Algorithm (GA) [39]. The Genetic algorithm can deal with multiple individuals simultaneously in a population, which means

Experimental Setup and Test Model
In the present study, a comprehensive experimental campaign was performed to validate the numerical optimized results. The experimental campaign was conducted in an ejector-driven and straight-flow icing wind tunnel, as shown schematically in Figure 4. The icing wind tunnel consisted of three parts, namely the ejector-driven wind tunnel system, the temperature adjustment system and the spraying system. The ejector-driven wind tunnel system was controlled by a pressure-regulating Energies 2020, 13, 5911 6 of 19 valve in front of the ejector, and the wind speed in the test section can be regulated from 0 to 90 m/s. The temperature adjustment system was manipulated by another pressure-regulating valve ahead of the cooling turbine, and the ambient temperature inside the insulated chamber can be adjusted from −25 to 0 • C. The spraying system can simulate the icing weather condition of which the liquid water content (LWC) varies from 0.5 to 2 g/m 3 and the median volumetric diameter (MVD) varies from 15 to 40 µm. Additional details of the icing wind tunnel are provided in [40].

Experimental Setup and Test Model
In the present study, a comprehensive experimental campaign was performed to validate the numerical optimized results. The experimental campaign was conducted in an ejector-driven and straight-flow icing wind tunnel, as shown schematically in Figure 4. The icing wind tunnel consisted of three parts, namely the ejector-driven wind tunnel system, the temperature adjustment system and the spraying system. The ejector-driven wind tunnel system was controlled by a pressure-regulating valve in front of the ejector, and the wind speed in the test section can be regulated from 0 to 90 m/s. The temperature adjustment system was manipulated by another pressure-regulating valve ahead of the cooling turbine, and the ambient temperature inside the insulated chamber can be adjusted from −25 to 0 °C. The spraying system can simulate the icing weather condition of which the liquid water content (LWC) varies from 0.5 to 2 g/m 3 and the median volumetric diameter (MVD) varies from 15 to 40 μm. Additional details of the icing wind tunnel are provided in [40]. An NACA0012 airfoil with a chord length of 0.3 m and a span length of 0.2 m was designed as a scaled test model. As shown in Figure 5, the test model comprised one main body and two spliced parts. The main body was 20 mm wide and made of polymethyl methacrylate (PMMA). The spliced parts on both sides were 40 mm wide and made of polyurethane foam. Near the front, the insides of the main body and the spliced parts were hollowed out for wiring and temperature sensor installation. Three parts were bolted together and mounted vertically in the test section of the icing wind tunnel. An NACA0012 airfoil with a chord length of 0.3 m and a span length of 0.2 m was designed as a scaled test model. As shown in Figure 5, the test model comprised one main body and two spliced parts. The main body was 20 mm wide and made of polymethyl methacrylate (PMMA). The spliced parts on both sides were 40 mm wide and made of polyurethane foam. Near the front, the insides of the main body and the spliced parts were hollowed out for wiring and temperature sensor installation. Three parts were bolted together and mounted vertically in the test section of the icing wind tunnel. The internal structure, as well as the material composition inside the main body, is illustrated in Figure 5. The outermost layer was the erosion shield, made of stainless steel. To observe and record the location and the distribution of the water film during the experiment, marks were curved on the erosion shields of each spliced part. Next to the erosion shield was the electric heating film layer. The heating film was made of nickel-chromium alloys with the polyimide tape covering its both sides. The size of the heating part was 20 mm × 5 mm. The resistance of each heating film was calibrated before the experiment, and 24 V DC power supplies with the precision of ±0.01V were applied to control the heating power of each heating film, of which the maximum value can be reached to 60 kW/m 2 . Since the NACA0012 airfoil was symmetrical, only one side needed to be investigated at the angle of attack (AOA) = 0°. Due to the size limitation of experimental conditions, the protected area on one side of the main body was divided into 8 heating zones to simulate various heating power The internal structure, as well as the material composition inside the main body, is illustrated in Figure 5. The outermost layer was the erosion shield, made of stainless steel. To observe and record the location and the distribution of the water film during the experiment, marks were curved on the erosion shields of each spliced part. Next to the erosion shield was the electric heating film layer. The heating film was made of nickel-chromium alloys with the polyimide tape covering its both sides. The size of the heating part was 20 mm × 5 mm. The resistance of each heating film was calibrated before the Energies 2020, 13, 5911 7 of 19 experiment, and 24 V DC power supplies with the precision of ±0.01V were applied to control the heating power of each heating film, of which the maximum value can be reached to 60 kW/m 2 . Since the NACA0012 airfoil was symmetrical, only one side needed to be investigated at the angle of attack (AOA) = 0 • . Due to the size limitation of experimental conditions, the protected area on one side of the main body was divided into 8 heating zones to simulate various heating power distributions. Each heating zone was covered by an electric heating film which is displayed in Figure 6. The heating zone lengths streamwise are presented in Table 1. On the other side of the main body, the protected area was covered by a larger electric heating film of which the size was 80 mm × 20 mm. Similarly, the leading edges of the spliced parts were also wrapped by the electric heating films, as shown in Figure 6. These heating films were used to keep as much of the test model clean of ice to prevent ice accretion from interference in the measurement on the main body. Material properties of each layer are presented in Table 2.
The internal structure, as well as the material composition inside the main body, is illustrated in Figure 5. The outermost layer was the erosion shield, made of stainless steel. To observe and record the location and the distribution of the water film during the experiment, marks were curved on the erosion shields of each spliced part. Next to the erosion shield was the electric heating film layer. The heating film was made of nickel-chromium alloys with the polyimide tape covering its both sides. The size of the heating part was 20 mm × 5 mm. The resistance of each heating film was calibrated before the experiment, and 24 V DC power supplies with the precision of ±0.01V were applied to control the heating power of each heating film, of which the maximum value can be reached to 60 kW/m 2 . Since the NACA0012 airfoil was symmetrical, only one side needed to be investigated at the angle of attack (AOA) = 0°. Due to the size limitation of experimental conditions, the protected area on one side of the main body was divided into 8 heating zones to simulate various heating power distributions. Each heating zone was covered by an electric heating film which is displayed in Figure  6. The heating zone lengths streamwise are presented in Table 1. On the other side of the main body, the protected area was covered by a larger electric heating film of which the size was 80 mm × 20 mm. Similarly, the leading edges of the spliced parts were also wrapped by the electric heating films, as shown in Figure 6. These heating films were used to keep as much of the test model clean of ice to prevent ice accretion from interference in the measurement on the main body. Material properties of each layer are presented in Table 2.  To obtain the surface temperature distribution without disturbing the heat transfer on the airfoil surface, several T-type thermocouples with the measurement range from −50 to 200 • C and the precision of ±0.5 • C, were installed beneath the heating film inside the main body. Each thermocouple was attached to a small piece of copper block which featured excellent thermal conductivity. Furthermore, the main body interior was filled with polyurethane elastic sealing foam which had a good adiabatic performance. As a consequence, the thermocouples were able to directly and precisely measure the surface temperature, which can be used to validate the numerical simulation and the optimization method.

Results and Discussions
According to the configuration of the icing wing tunnel and the test model, the computational domain of the numerical simulation was established and the boundary conditions were defined, which are displayed in Figure 7. After various tests, a structured grid with 3 million cells was used to calculate the flow field around the test model. The cells in the first boundary layer were refined to y + ≈ 1. The experimental conditions are listed in Table 3. To validate the numerical simulation and the optimization methods, the numerical results were compared with the experimental results.

Results and Discussions
According to the configuration of the icing wing tunnel and the test model, the computational domain of the numerical simulation was established and the boundary conditions were defined, which are displayed in Figure 7. After various tests, a structured grid with 3 million cells was used to calculate the flow field around the test model. The cells in the first boundary layer were refined to y + ≈ 1. The experimental conditions are listed in Table 3. To validate the numerical simulation and the optimization methods, the numerical results were compared with the experimental results.    The convective heat transfer coefficient plays an important role in the calculation of the anti-icing heat load. In the icing wind tunnel test, the experimental cases were conducted without the spraying system turned on. When the steady state conditions were achieved, the heating films were adjusted to a uniform temperature, typically in the range of 32 to 38 • C [41]. Roughly three minutes were required to record the temperature and heat flux distribution. The convective heat transfer coefficient can be derived by where q i is the power density of the ith heating film. The measurements of the convective heat transfer coefficient were carried out from Case 1 to Case 3. The comparisons of the results of the numerical simulation and the experiment are displayed in Figure 8. As can be seen, the distributions of the convective heat transfer coefficient in all cases show a decreasing tendency, which drops more at the leading edge and less backwards. This kind of distribution means that the leading edge has intensive heat transfer. The maximum error between the numerical and experimental results is less than 15%. In general, there is an acceptable agreement between the numerical and the experimental results.
Energies 2020, 13, 5911 9 of 19 coefficient were carried out from Case 1 to Case 3. The comparisons of the results of the numerical simulation and the experiment are displayed in Figure 8. As can be seen, the distributions of the convective heat transfer coefficient in all cases show a decreasing tendency, which drops more at the leading edge and less backwards. This kind of distribution means that the leading edge has intensive heat transfer. The maximum error between the numerical and experimental results is less than 15%. In general, there is an acceptable agreement between the numerical and the experimental results. Wrap distance(m)

Water Droplet Impinging Limit
The water droplet impinging limit in each case was measured by the rime ice method, which can form the thin rime ice on the surface of the test model. The pictures of the experimental results are shown in Figure 9, which are 17.5 mm in Case 1, 20 mm in Case 2 and 22.5 mm in Case 3. As can be seen, the water droplet impinging limits increase with the increase in the wind speed. There is more ice accretion at the leading edge and less at the impinging limit, which means the water droplet collection efficiency has a decreasing tendency along with the chord. The comparison between experimental results and numerical simulation of the water droplet collection efficiencies are shown in Figure 10. The results indicate that the numerical water droplet impinging limits are generally consistent with the experimental ones, which means the numerical simulation and the experiments are reliable.

Water Droplet Impinging Limit
The water droplet impinging limit in each case was measured by the rime ice method, which can form the thin rime ice on the surface of the test model. The pictures of the experimental results are shown in Figure 9, which are 17.5 mm in Case 1, 20 mm in Case 2 and 22.5 mm in Case 3. As can be seen, the water droplet impinging limits increase with the increase in the wind speed. There is more ice accretion at the leading edge and less at the impinging limit, which means the water droplet collection efficiency has a decreasing tendency along with the chord. The comparison between experimental results and numerical simulation of the water droplet collection efficiencies are shown in Figure 10.
The results indicate that the numerical water droplet impinging limits are generally consistent with the experimental ones, which means the numerical simulation and the experiments are reliable.

Surface Temperature
The heating power distributions used for the surface temperature validation are listed in Table  4. Figure 11 displays the numerical and experimental results of surface temperature distributions. As can be seen, the surface temperature distributions had humps in the middle of the protected area, which was caused by completely evaporating the runback water. There were small discrepancies between the experimental and numerical results in the low-temperature regions, while the differences became larger in the high-temperature regions. This was because the thermal conduction from the main body to the spliced parts of the test model was not considered in the numerical simulation. The maximum error between the numerical and experimental results was about 8.8%, and the average error was about 4.1%, which was acceptable.

Surface Temperature
The heating power distributions used for the surface temperature validation are listed in Table 4. Figure 11 displays the numerical and experimental results of surface temperature distributions. As can be seen, the surface temperature distributions had humps in the middle of the protected area, which was caused by completely evaporating the runback water. There were small discrepancies between the experimental and numerical results in the low-temperature regions, while the differences became larger in the high-temperature regions. This was because the thermal conduction from the main body to the spliced parts of the test model was not considered in the numerical simulation. The maximum error between the numerical and experimental results was about 8.8%, and the average error was about 4.1%, which was acceptable.

The Optimal Heating Power Distribution
To optimize the heating power distribution, the experiments started with the critical dry surface state which is the boundary between the dry surface mode and the evaporative mode. The dry surface mode was easily accessible in the experiment. Firstly, the heat power distribution was set up at a high

The Optimal Heating Power Distribution
To optimize the heating power distribution, the experiments started with the critical dry surface state which is the boundary between the dry surface mode and the evaporative mode. The dry surface mode was easily accessible in the experiment. Firstly, the heat power distribution was set up at a high level to ensure no runback water within the protected area. Then, from 8# to 1#, the power density of each heating film was reduced as low as possible until the critical dry surface state emerged. Then, the power density of 1# heating film was turned down and the power density of 2# heating film was turned up till the runback water covered 1# heating zone. After that, the power density of 1# and 2# heating films were turned down and the power density of 3# heating film was turned up till the runback water covered 2# heating zone, and so on till the runback water covered the entire protected area. Figure 12 shows the experimental heating power distributions with the runback water covering the different heating zones in each case. The corresponding anti-icing heat loads of the heating power distributions are also displayed in Figure 12. The results show that the anti-icing heat loads had a decreasing tendency as the runback water covered more protected areas. In Case 1, the anti-icing heat load had a minimum value when the runback water covered 4# heating zone and was reduced by 23.8% compared to the initial heating power distribution. In Cases 2 and 3, the anti-icing heat loads had minimum values when the runback water covered 5# heating zones and were reduced by 19.6% and 18.6%, respectively, compared to the initial heating power distributions. It can be seen from Table 2 and Figure 10 that the water droplet impinging limit of Case 1 was located in 4# heating zone while the water droplet impinging limits of Cases 2 and 3 were located in 5# heating zones. This indicates that the anti-icing heat load can be reduced to a minimum value as the runback water reaches the water droplet impinging limit. Considering that the PMMA has a poor temperature tolerance among the materials of the test model, the maximum temperature in Equation (15) is set to be 353 K [42]. The comparison of the optimized and the experimental results are displayed in Figure 13. The optimal heating power distributions of the optimization and the experiment were generally consistent in the trends, which were both low at the leading edge, high at the water droplet impinging limit and zero at the end of  Considering that the PMMA has a poor temperature tolerance among the materials of the test model, the maximum temperature in Equation (15) is set to be 353 K [42]. The comparison of the optimized and the experimental results are displayed in Figure 13. The optimal heating power distributions of the optimization and the experiment were generally consistent in the trends, which were both low at the leading edge, high at the water droplet impinging limit and zero at the end of the protected area. The anti-icing heat loads in the experiment were a little larger than those in the optimization because of the differences of the heating power distributions at the leading edge and the impinging limit. The difference between the optimizations and the experiments in the heating power distribution is that the experimental results are higher at the leading edge and lower at the water droplet impinging limit. As a matter of fact, during the experiment, when the heating power density at the leading edge was reduced excessively, the runback water film became unstable and easily turned into the rivulets, which failed the ice protection. However, the numerical simulation did not consider this process, which made the differences in the anti-icing heat load as well as the surface temperature.

Mechanism of the Anti-Icing Heat Load Reduction
According to the results of the optimization and the experiments, the optimal heating power distribution has such a characteristic that is low at the leading edge, high at the water droplet impinging limit and zero at the end of the protected area. This tendency can lead to a less anti-icing heat load. To find out the mechanism behind this consequence, it is necessary to figure out the composition of the anti-icing heat load. Integrating the control volume within the protected area, Equation (2) can be written as: where A refers to the protected area. The third term in Equation (16) should be zero because neither the inflow enthalpy nor the outflow enthalpy exists in the evaporative mode, except the enthalpy of water droplet impingement. According to References [24] and [31], q w and q a are usually much smaller than the others, which can also be neglected. Therefore, Equation (16) can be simplified as:

Mechanism of the Anti-Icing Heat Load Reduction
According to the results of the optimization and the experiments, the optimal heating power distribution has such a characteristic that is low at the leading edge, high at the water droplet impinging limit and zero at the end of the protected area. This tendency can lead to a less anti-icing heat load. To find out the mechanism behind this consequence, it is necessary to figure out the composition of the anti-icing heat load. Integrating the control volume within the protected area, Equation (2) can be written as: where A refers to the protected area. The third term in Equation (16) should be zero because neither the inflow enthalpy nor the outflow enthalpy exists in the evaporative mode, except the enthalpy of water droplet impingement. According to References [24] and [31], q w and q a are usually much smaller than the others, which can also be neglected. Therefore, Equation (16) can be simplified as: or where Q t is the total anti-icing heat load, Q evap is the evaporative heat load and Q conv is the convective heat load. It follows that the anti-icing heat load is mainly composed of the evaporative and convective heat loads. As can be seen from Equation (10), Q evap includes both the latent part and the sensible part. In the evaporative mode, supposing the thermal anti-icing system operates in a stable environment, the total amount of the collected water droplets remains constant. Consequently, Q evap does not change much even under different heating power distributions. As a result, the reduction in the anti-icing heat load can only be the decrease in Q conv . It is well known that Q conv is influenced by the heating range, the surface temperature as well as the convective heat transfer coefficient. Analyzing the influences of these factors is conducive to reveal the mechanism of the anti-icing heat load reduction.

Heating Range
The anti-icing heat loads reached the minimum value as the runback water exactly covered the water droplet impinging limit as shown in Figure 12. It indicates that the heating range is important to the anti-icing heat load reduction. To better investigate its influence, a uniform heating power distribution was experimentally tested. In this experiment, only the wet area covered by the runback water was heated, and the dry area was not heated. Hence, the heating range can be represented by the runback water length. Additionally, when the runback water length was less than the water droplet impinging limit, the heating range was fixed at the limit.
The experimental and the numerical results are displayed in Figure 14, which presents the relationship between the anti-icing heat load and the range of the runback water length. As can be seen, both the experimental and the numerical anti-icing heat loads had minimum values as the runback water lengths were at the water droplet impinging limits. The evaporative heat loads barely change with the heating power density. The trend of the convective heat load is similar to the anti-icing heat load, which has an inflection point corresponding to the impinging limit. This can be explained in two aspects. On one hand, when the heating power density was higher than the value at the impinging limit, the runback water length decrease. However, the heating range was fixed at the impinging limit, so the anti-icing heat load increased. On the other hand, when the heating power density was lower than the value at the impinging limit, both the runback water length and the heating range increased, as did the anti-icing heat load. Therefore, decreasing the heating range outside the water droplet impinging limit can effectively reduce the anti-icing heat load. As the heating range is equal to the impinging limit, the anti-icing heat load is the minimum. water droplet impinging limit as shown in Figure 12. It indicates that the heating range is important to the anti-icing heat load reduction. To better investigate its influence, a uniform heating power distribution was experimentally tested. In this experiment, only the wet area covered by the runback water was heated, and the dry area was not heated. Hence, the heating range can be represented by the runback water length. Additionally, when the runback water length was less than the water droplet impinging limit, the heating range was fixed at the limit. The experimental and the numerical results are displayed in Figure 14, which presents the relationship between the anti-icing heat load and the range of the runback water length. As can be seen, both the experimental and the numerical anti-icing heat loads had minimum values as the runback water lengths were at the water droplet impinging limits. The evaporative heat loads barely change with the heating power density. The trend of the convective heat load is similar to the antiicing heat load, which has an inflection point corresponding to the impinging limit. This can be explained in two aspects. On one hand, when the heating power density was higher than the value at the impinging limit, the runback water length decrease. However, the heating range was fixed at  Figure 15 displays the surface temperature distributions and the proportions of the convection and evaporation under the optimal heating power distributions. As can be seen, the evaporation has a higher proportion in the high-temperature regions, while the convection has a higher proportion in the low-temperature regions. This suggests that the temperature has different impacts on convection and evaporation. To better illustrate the influences of the surface temperature, the following parameters are selected as representative values: α = 200 W/(m 2 ·K), T e = 263.15 K, p e = 101.325 kPa. The changes in the evaporative and convective heat fluxes with the surface temperature are displayed in Figure 16. Noted that the evaporative heat flux can directly denote the evaporative mass transfer of the runback water according to Equation (10). According to Equations (4), (5), (10) and (11), the convective heat flux has a linear relationship while the evaporative heat flux has an exponential relationship with the surface temperature. As can be seen, when the same amount of water is evaporated simultaneously at both high and low temperatures, the temperature rise required in the high-temperature region is less than that in the low-temperature region. Moreover, the additional convective heat flux rise in the high-temperature region is less than that in the low-temperature region. For instance, as the evaporative heat flux increases from 0 to 30 kW/m 2 , the temperature will have an increment of about 49 K. However, as the evaporative heat flux increases from 60 to 90 kW/m 2 , the temperature rise is less than 10 K, which leads to a significant reduction in the convective heat load. Therefore, evaporating most runback water in the high-temperature region can surely reduce the convective heat load.

Surface Temperature
convective heat flux rise in the high-temperature region is less than that in the low-temperature region. For instance, as the evaporative heat flux increases from 0 to 30 kW/m 2 , the temperature will have an increment of about 49 K. However, as the evaporative heat flux increases from 60 to 90 kW/m 2 , the temperature rise is less than 10 K, which leads to a significant reduction in the convective heat load. Therefore, evaporating most runback water in the high-temperature region can surely reduce the convective heat load.

Convective Heat Transfer Coefficient
The surface temperature distribution is not only influenced by the heating power distribution, but also by the characteristics of the convective heat transfer. Supposed that a constant heating power density of 30 kW/m 2 is provided within the protected area, Te and pe are selected as representative values: Te = 263.15 K, pe = 101.325 kPa. According to Equations (10), (11) and (18), the changes of qevap and qconv as well as the surface temperature with the convective heat transfer coefficient are displayed in Figure 17a. As can be seen with the increase in the convective heat transfer coefficient, qconv tends to increase, while the surface temperature tends to decrease, so does qevap. The results show that qevap has a higher proportion in the region with the small convective heat transfer coefficient, while qconv has a higher proportion in the region with the large convective heat transfer coefficient.

Convective Heat Transfer Coefficient
The surface temperature distribution is not only influenced by the heating power distribution, but also by the characteristics of the convective heat transfer. Supposed that a constant heating power density of 30 kW/m 2 is provided within the protected area, T e and pe are selected as representative values: T e = 263.15 K, p e = 101.325 kPa. According to Equations (10), (11) and (18), the changes of q evap and q conv as well as the surface temperature with the convective heat transfer coefficient are displayed in Figure 17a. As can be seen with the increase in the convective heat transfer coefficient, q conv tends to increase, while the surface temperature tends to decrease, so does q evap . The results show that q evap has a higher proportion in the region with the small convective heat transfer coefficient, while q conv has a higher proportion in the region with the large convective heat transfer coefficient.
The surface temperature distribution is not only influenced by the heating power distribution, but also by the characteristics of the convective heat transfer. Supposed that a constant heating power density of 30 kW/m 2 is provided within the protected area, Te and pe are selected as representative values: Te = 263.15 K, pe = 101.325 kPa. According to Equations (10), (11) and (18), the changes of qevap and qconv as well as the surface temperature with the convective heat transfer coefficient are displayed in Figure 17a. As can be seen with the increase in the convective heat transfer coefficient, qconv tends to increase, while the surface temperature tends to decrease, so does qevap. The results show that qevap has a higher proportion in the region with the small convective heat transfer coefficient, while qconv has a higher proportion in the region with the large convective heat transfer coefficient.  Given the constant mass flow rate of the evaporation, of which m ʹ evap is equal to 0.01 kg/(m 2 •s), the changes of qevap and qconv, as well as the surface temperature with the convective heat transfer coefficient, are displayed in Figure 17b. As can be seen with the increase in the convective heat transfer coefficient, qconv tends to increase, while the surface temperature tends to decrease. However, due to constant m ʹ evap , the heating power density tends to increase. The results show that, when evaporating the same amount of water, the anti-icing heat load in the region with the small convective Given the constant mass flow rate of the evaporation, of which . m evap is equal to 0.01 kg/(m 2 ·s), the changes of q evap and q conv , as well as the surface temperature with the convective heat transfer coefficient, are displayed in Figure 17b. As can be seen with the increase in the convective heat transfer coefficient, q conv tends to increase, while the surface temperature tends to decrease. However, due to constant . m evap , the heating power density tends to increase. The results show that, when evaporating the same amount of water, the anti-icing heat load in the region with the small convective heat transfer coefficient is smaller than that in the region with the large convective heat transfer coefficient.
According to the influence of the surface temperature, evaporating the runback water in the high-temperature region can reduce the anti-icing heat load. The surface temperature can easily be raised in the region with the small convective heat transfer coefficient. For the flow field around the airfoil surface, the convective heat transfer coefficient is always large in the vicinity of the leading edge and small at the end of the protected area, due to the increase in the boundary layer thickness and the turbulence. The trend of the optimal heating power distribution is opposite to that of the convective heat transfer coefficient. Thus, the anti-icing heat load can be reduced.

Conclusions
In this paper, the optimal heating power distribution that leads to the less anti-icing heat load in the evaporative mode was investigated both in the numerical optimization and the experiment. By analyzing the influences of the key factors on the convective heat load, the mechanism of the anti-icing heat load reduction was revealed.
(1) The reduction in the anti-icing heat load is the decrease in the convective heat load. As the thermal anti-icing system operating in the evaporative mode, the evaporative heat load keeps nearly constant. The convective heat load can be varied with the surface temperature distribution which is influenced by the heating power distribution. Thus, there must be an optimal heating power distribution that can reduce the convective heat load as well as the anti-icing heat load. (2) The optimal heating power density obtained by the numerical optimization and the experiment has such characteristics that are low at the leading edge, high at the water droplet impinging limit and zero at the end of the protected area. These characteristics are mainly influenced by the heating range, the surface temperature and the convective heat transfer coefficient. In the evaporative mode, decreasing the heating range outside the protected area can reduce the anti-icing heat load effectively. As the heating range is decreased to the impinging limit, the anti-icing heat load has the minimum value. Due to the different impacts of the surface temperature on the evaporative and convective heat fluxes, it is better to evaporate the runback water in the high-temperature region, which can lead to a lesser additional convective heat load. The surface temperature distribution is affected by the convective heat transfer coefficient distribution of which the trend is high at the leading edge and decreases chordwise around the airfoil surface. As a consequence, the optimal heating power distribution has the opposite trend with the convective heat transfer coefficient distribution.
The present investigations on the mechanism of the anti-icing heat load reduction can provide valuable guidance for the design of the thermal anti-icing system, which can be used for wings, stabilizers and engine inlets of the aircraft as well as wind blades and cowlings of the wind turbines, or any component that has an aerodynamic profile. The configuration of the thermal anti-icing system can be improved according to the presented results, which can achieve the purpose of saving energy and reducing weight. However, the results obtained in this paper have certain limitations. Considering the anti-icing requirement in practical applications, the margins of both the anti-icing heat load and the protected area ought to be increased appropriately, and extra constraints such as maximum heating power density should also be taken into consideration.