Abstract
A novel 3-D unsteady model of in-flight electrothermal deicing process is presented in this paper to simulate the conjugate mass and heat transfer phenomena of water film runback, phase change, and solid heat conduction. Mathematical models of water film runback and phase change are established and solved by means of a loosely coupled method. At the current time step, solid heat conduction, water film runback, and phase change are iteratively solved until the heat boundary condition reaches convergence, then the temperature distribution and ice shape at the moment are obtained, and the calculation of the next time step begins subsequently. A deicing process is numerically simulated using the present model following an icing tunnel experiment, and the results match well with those in the literatures, which validate the present model. Then, an in-flight deicing process is numerically studied to analyze the effect of heating sequence.
1. Introduction
The water droplet in clouds may remain in a liquid state even if the temperature is below freezing point due to the surface tension and a lack of condensation nucleus. Supercooled water droplet would solidify when impinging on windward surfaces of aircraft, which would cause the deterioration in the aerodynamic performance due to the change of aerodynamic configuration [1]. Aircraft icing would impose serious adverse effects on flight safety, which has long been recognized as an important issue to prevent.
In view of the serious threats that ice accretion would impose on flight safety, ice protection methods must be applied to prevent or control the ice accretion. Hot air anti-icing method is widely applied in commercial jets. The bleed air from engine compressor impinges on the structure to heat the surface, so that the water evaporates or stays in a liquid state rather than freezes. Electrothermal ice protection method uses heating pads, which are incorporated in the multilayer structure, to heat the surfaces. Electrothermal pads are available in both anti-icing mode with a constant power density and deicing mode with a periodic power density. Besides, some other deicing methods have been developed, such as the thermomechanical expulsion deicing method [2], while they are not widely applied. A large amount of bleed air is needed in hot air anti-icing system, which would significantly affect the performance of an engine, especially during the takeoff or landing period. Meanwhile, the jet flow would cause excessive temperature at the impingement point, which may damage the composite materials that are widely used in a modern aircraft [3] or even the aluminum skin. The electrothermal ice protection method has a higher heating efficiency and a lower power density, which results in the advantages in structure safety and energy efficiency. As the concepts of more-electric-aircraft and all-electric-aircraft are increasingly employed in a modern aircraft design, the electrothermal ice protection method is drawing more attention. The ice protection system of B-787 [4], which is based on the electrothermal method, serves as an example.
Since the icing tunnel experiments and the flight experiments are complex, expensive, and not able to cover all environment conditions, the investigation of a deicing process by an experimental method is quite limited and seldom reported to date. The numerical simulation method, which is time-efficient and cost-effective, becomes an important tool for the design of the ice protection system. The in-flight deicing process is a conjugate mass and heat transfer phenomenon, which is very complex and consists of a variety of coupled processes such as the air-droplet flow, water film runback and phase transition, and multilayer solid heat conduction. The phase state varies with both spatial locations and time. Due to the periodic power density, spatial distribution of heaters, and sustained droplet impingement, phase transition phenomena, such as evaporation, solidification, liquefaction, and sublimation, may coexist on the surface of the protection area, which makes it complex and typically unsteady.
The early study of electrothermal deicing concentrated on the solution of multilayer structure heat conduction. Stallabrass [5] numerically analyzed the heat conduction during the deicing process, focusing on the effects of multilayer materials on the distribution of skin temperature, especially the effects of inner and outer insulating layers. Then, the latent heat was introduced to simulate the phase change during the deicing process by Baliga [6]. As the research went in depth, Marano [7] and Roelke et al. [8] developed the “enthalpy method,” in which the unified energy equation was applied in the whole computational domain and the phase interface was determined by enthalpy distribution. Based on the enthalpy method, plenty of researches were conducted to simulate the phase transition process. Chao [9], Leffel [10], and Masiulaniec [11] applied the method to 2-D cases and analyzed the effects of surface curvature and heating pad spatial distribution. Chang [12, 13] numerically studied the 2-D deicing process by a finite volume method and analyzed the effects of the heating sequence and power density. Yaslik et al. [14] studied a 3-D deicing process, using a Douglas method to discretize the heat conduction equation. Xiao [15] applied a porous medium method to simulate the ice melting process. As is briefly reviewed above, researchers have done much work on the study of multilayer structure heat conduction and ice melting process based on enthalpy method, while the study of the in-flight deicing process, which involves sustained droplet impingement and water film runback, is quite rare.
Wright et al. [16] developed a 2-D deicing model, which was based on the classic Messinger mass and energy balance models [17]. The mass and energy balance equations took into account the input water caused by droplet impingement, but the water film runback mechanism was not taken into consideration. Habashi et al. [18–20] developed a conjugate model for the in-flight deicing process. The airflow and droplet impingement were calculated as initial solutions, and the solid conduction was solved, which was coupled with a water phase transition in a loosely coupled way. Finite element solvers were applied to the solution as the modules of the commercial software FENSAP-ICE, while the detailed solution method was not available in the open literature.
As electrothermal deicing method is increasingly employed in the new generation of aircraft, the study on the mechanism and simulation methods is urgently needed. While most current studies on the deicing process concentrated on the multilayer structure heat conduction and ice melting process, studies available on the in-flight deicing were quite rare.
This paper focuses on the conjugate heat transfer mechanism of the in-flight electrothermal deicing. A novel 3-D unsteady model is established based on the water film runback dynamic mechanism and phase transition thermodynamic model and solved by a loosely coupled method. Using the present model, the deicing process is numerically investigated, which contributes to a better understanding of electrothermal deicing so as to guide the design and optimization of an aircraft ice protection system.
2. Mathematical Model
The in-flight deicing process is a conjugate mass and heat transfer process, which involves air-droplet flow, solid heat conduction, water film runback, and the phase transition. The mass and energy balance, as is briefly shown in Figure 1, is determined by a variety of factors including convection, droplet impingement, evaporation, solidification, heat conduction, and film runback. Due to the periodic power density and spatial distribution of heaters, the phase transition process is typically unsteady.
2.1. Air-Droplet Flow
Since the volume fraction of droplet is very small, typically under 10−6 for icing conditions, the air-droplet flow can be solved by a one-way coupled method [21]. Reynolds-averaged Navier–Stokes equations (RANS) are applied to solve the air flow field. The mass and momentum equations are expressed as (1), and the energy equation is expressed as (2): where the Reynolds stress is defined by Boussinesq as where is the turbulent viscosity, is the average velocity, and is the turbulence kinetic energy.
The continuity and momentum equations of droplet flow are expressed as [22] where is the density of water, is the droplet volume fraction, is the velocity vector of droplet, is the velocity vector of air, is the external body forces exerted on droplet, such as the gravity or inertial force, and is the air-droplet momentum exchanger coefficient defined as where is the dynamic viscosity of air, is the diameter of droplet, and is the drag function. The Schiller and Naumann model is adopted here. is the relative Reynolds number and is given as
2.2. Water Film Runback and Phase Transition
The differential form of continuity equation for runback water film on the surface of ice protection area is expressed as where is the velocity of water film, is the water mass of droplet impingement, is the icing rate, and is the evaporation rate. The water film is incompressible, and the development of water film is quick due to a quite small thickness. Therefore, the unsteady water mass term is neglected, and the continuity equation is derived by integrating the differential form in a control volume. where is the total mass of water entering the current control volume from adjacent ones per unit time, is the total mass of water flowing out of the control volume, and the mass flux through a face can be expressed as where is the film thickness, is the edge length, and represents the unit vector normal to face .
The mass flux of droplet impingement is determined by local collection efficiency , as expressed as where is the velocity of droplet at far field, is the liquid water content, and is the area of control volume.
The evaporation rate is determined by the Chilton-Colburn analogy theory where the convective mass transfer coefficient is obtained from the convective heat transfer coefficient as shown below. where is the Schmidt number. The convective heat transfer coefficient is calculated according to the surface heat flux and the temperature difference. where is the recovery coefficient; is the air specific heat ratio, of which the value is 1.4; and is the Mach number of air flow. Then the evaporation rate is derived as where is the air molecular weight, is the air specific heat, is the vapor saturated pressure, is the total pressure, is the total temperature, and is the relative humidity, of which the value on water film surface is 1.
The momentum equation of water film runback is given by incompressible Navier–Stokes equation as where is the pressure, is the unit tensor, is the coordinate normal to surface, and is the gravitational acceleration.
The air-film boundary condition is defined as where is the shear stress, is the coefficient of tension, and is the water film surface curvature which is given by
The solid-film boundary is depicted under the nonslip boundary condition.
Previous studies suggest that the effects of pressure gradient should only be considered for water film with a large thickness [23]. During the deicing process, the water film is very thin; therefore, the terms of gradient and gravity are negligible, of which the effect is slight compared with other terms such as shear stress. The tangential velocity gradient is also very small; therefore, the momentum diffusion term is negligible [24]. The value of surface curvature is generally very small for film which entirely wets the surface, and the shear stress is the dominant factor at the air-film interface. Under such condition, the momentum equation and boundary condition of water film are simplified as
The velocity distribution normal to the surface is derived as
The average velocity is obtained by integrating the above equation along the film thickness direction, as expressed as
The energy balance of water film is described by the following differential equation:
Integrating the above equation, the energy balance in a control volume is expressed as where is the specific heat of water, and are the temperature at current time step and previous time step, is the deicing heat flux from solid, obtained from the calculation of solid heat conduction during the iterations, and is the energy carried by runback water entering the current control volume through the face , and it can be expressed as
is the energy of impinging water, which consists of two parts: the internal energy and the kinetic energy, as expressed below.
is the convective heat flux obtained by the convective heat transfer coefficient and the temperature difference between surface temperature and recover temperature .
The latent heat of solidification , the out-flow energy , and the latent heat of evaporation are related to the phase condition of the control volume, which are discussed as follows by considering the freezing fraction , namely the ratio of icing mass to the total entering water mass.
If and the surface is clean, which means currently that there is no ice accretion, the water in the control volume remains at a liquid state, icing rate is zero, and the energy of runback water and evaporation can be expressed as where is the latent heat of water evaporation.
If and there is ice accretion in the current volume, the ice melts at a liquid-solid mixed state and the control volume remains at the phase transition temperature. where is the melting rate and is the solidification latent heat.
If 0 < f < 1, part of water freezes, and the control volume is at the phase transition temperature.
If all the water freezes and the energy terms can be expressed as where is the specific heat of ice and is the sublimation latent heat.
2.3. Heat Conduction
The unsteady heat conduction through the multilayer materials is expressed as where is the enthalpy of the material, is the thermal conductivity, and is the heat source term. For electrothermal deicing process, the source term is determined by the spatial location of heaters and the power sequence.
The Dirichlet heat boundary condition of heat conduction is provided by the solution of water film runback and phase transition during the coupled solution. In return, the deicing heat flux is calculated and sent to the calculation of water film. The deicing heat flux is obtained at the boundary of a solid structure as shown in where is the normal vector of the boundary surface.
3. Solution Procedure
During the solution of the above unsteady conjugate heat transfer model, the water film runback and phase transition are coupled with the solid heat conduction, due to the fact that the solution of the solid heat conduction provides the interface heat flux which is needed in the water film energy balance equation; on the other hand, the boundary condition of the solid heat conduction is provided by the solution of the film runback and phase transition. A loosely coupled method is applied to solve the present model, in which both the water film runback and the heat conduction are iteratively calculated until the heat boundary condition reaches convergence, and during each iteration, the surface temperature and heat flux are exchanged, as briefly shown in Figure 2.
Considering that the ice thickness during deicing process is typically controlled at a very small value, the effect of the ice shape on air-droplet flow is slight. As a result, the steady airflow solution is computed as an initial condition and is assumed unchanged during the simulation. Besides, it is accurate as long as the ice thickness does not exceed a certain limit due to a protection failure. The RANS equations of airflow are discretized, using a finite volume method in a second order upwind scheme. To simulate the turbulent flow, the transition SST model, which is shown to obtain good results for wall-bounded flows, such as the airfoil, is utilized. A CFD solver FLUENT is used to solve the governing equations. The governing equations of droplet flow field are solved by the finite volume method using the User-Defined Scalar (UDS) transport equation. The droplet volume fraction and velocity components are set as the UDS, and the convective terms, diffusion terms, and source terms are defined by codes which are programmed using the User-Defined Functions (UDF). The solution of air-droplet flow is the initial condition of water film runback and solid heat conduction, and the data exchange is achieved by interpolation due to the difference between flow field mesh and solid mesh.
A loosely coupled method is applied to solve the deicing model. At the current time step, both the solid heat conduction and the water film runback are iteratively calculated until convergence. During each iteration, the surface temperature and heat flux are exchanged, and the boundary conditions are updated. The ice shape is obtained by the icing rate or melting rate, and the calculation of the next time step then begins. The flowchart is shown in Figure 3, and some brief introduction is provided below for the flowchart. (1)Solve the air-droplet flow and transfer the data by interpolation;(2)Loop all the surface control volumes and check them. If the input water mass is already known, assume an initial temperature, solve the mass and energy balance equations, update the value of temperature and phase state, and provide input conditions for adjacent volumes. Keep the calculation of water film runback and phase transition until the calculations of all control volumes are done; the temperature distribution at boundary surface is obtained;(3)Set Dirichlet boundary condition for solid heat conduction using the temperature distribution of step (2) and solve the heat conduction; calculate the deicing heat flux at boundary.(4)Update the deicing heat flux of water film energy equation using the data of step (3).(5)Repeat steps (2)–(4) until convergence, calculate ice accretion at the moment, and advance to the next time step.
4. Results and Discussion
Validations of the present model are conducted, and the results are compared with the experimental data. Then, a simulation of in-flight deicing process is conducted to optimize the heating sequence.
4.1. NASA Deicing Experiment in the Lewis Icing Research Tunnel
In order to perform the in-flight deicing simulation, a deicing experiment model is selected in the very rare records available in the open literature, which is conducted in the NASA Lewis Icing Research Tunnel (IRT) by Al-Khalil et al. [25].
The experiment model is a NACA 0012 airfoil with a chord of 0.914 m (36 in), and the environment conditions are listed in Table 1. The multilayer structure at the leading edge protection area is composed of four different layers, which are the erosion shield, the elastomer, the fiberglass, and the silicone foam, with a thickness of 0.2 mm, 0.56 mm, 0.89 mm, and 3.4 mm, respectively. The structure is shown in Figure 4, and the material properties are listed in Table 2. Seven heating pads are arranged in the elastomer layer, of which the heating sequence and power density are controlled independently. The length is 1.905 cm for heater 4; 2.54 cm for heaters 2, 3, 5, and 6; and 3.81 cm for heaters 1 and 7. The heating sequence is shown in Figure 5. Heater 4 acts as the heat blade, which keeps activated during the whole cycle to avoid ice accretion on the leading edge. The rear heating pads are activated alternately after 100 s.
The steady air-droplet flow was solved to obtain parameters such as the convective heat flux, the shear stress, and the local droplet collection efficiency. Then, the data were transferred to solid mesh by interpolation, and the unsteady deicing process was coupled solved. Structured grids were generated for the solution of air-droplet flow. To verify the mesh independence of the solution, three mesh files were applied, and the normal distance of the first layer is 0.01 mm (mesh a), 0.0075 mm (mesh b), and 0.005 mm (mesh c), respectively. The y+ at the ice protection area of all three mesh files was controlled around or lower than 1.
The simulated convective heat transfer coefficient curves are shown in Figure 6, which show good agreement with the experimental data and LEWICE solution. The results of three mesh files match well, indicating that the solutions are mesh independent. The convective heat transfer coefficient reaches the peak at the stagnation point and drops rapidly as it moves backwards due to the development of boundary layer. The experimental value is slightly larger than that of the simulated results, and the turbulence intensity of the icing tunnel test might be a possible explanation. The droplet collection efficiency was not recorded during the experiment, while it was simulated by NASA LEWICE code in the literature. The collection efficiency of the present model and LEWICE is shown in Figure 7, and the results match well. The value reaches its maximum at the stagnation point, and then the value decreases as it moves backwards. There is a droplet shadowed zone when the wrap distance exceeds 0.03 m, where the local collection efficiency is zero.
The simulated temperature of heater 3 is shown in Figure 8 with the experimental result and FENSAP simulated result provided in [18]. The temperature variation curve of each cycle is similar except for the beginning period when the solid structure starts to warm up. The result of the present model shows good agreement with the experimental data and FENSAP. The simulated ice shapes at 100 s, 110 s, and 120 s (the moment when heaters are turned on or off) are shown in Figure 9.
(a) t = 100 s
(b) t = 110 s
(c) t = 120 s
At the first stage (0–100 s), only the heat blade is activated. Due to the solid conduction and the runback liquid water, the temperature of heater 3 rises. Figure 9(a) shows that water remains liquid on the surface over the heat blade, and ice forms downstream where the heater is turned off. The runback ice covers part of the protection area. At the second stage (100–110 s), heaters 3 and 5 are activated and the temperature of heater 3 rapidly rises to 290 K. From Figure 9(b), it is observed that the ice over heater 3 melts and the runback ice area moves backward, due to the fact that the heating enlarges the water runback area. At the third stage (110–120 s), heater 3 is turned off and the temperature drops. Figure 9(c) shows that the runback ice, which forms during the second stage, melts and ice forms over heater 3 as it is turned off.
The simulated temperature of heater 3 is slightly higher than that of the experimental data. When heater 3 is off, the experimental temperature is around 270 K, and it is about 3 K lower than the simulated results. At this period, the surface is under a runback icing condition, which means a water-ice mixed state, and the temperature of such condition is set at 273.15 K (the freezing point) in the thermodynamic model during simulation. The possible reasons for the temperature difference between experiment and simulation are as follows: (1) The supercooled water film runback phenomenon, which is observed in the experiments [26–28], has not been considered in the present model. In simulation, the temperature of the runback water film would not be lower than the freezing point temperature (273.15 K). (2) The runback water might form beads or rivulet flow and does not completely wet the surface. Part of the surface is exposed to air convection, and the temperature is lower than that when the surface is completely wetted. (3) The temperature measuring point in the experiment is beneath the heating pad, and the precise location is not mentioned in the report; while in simulation, the average value of the heater area is calculated as the result. The average value might be larger when the heater is turned off, because the temperature at the margin is higher due to the heating of the adjacent heating pads.
4.2. Deicing Simulation under Different Heating Sequences
The above simulations validate the present model and the solution method. In this section, another deicing simulation is conducted under in-flight environment conditions to analyze the deicing performance of different heating sequences. The environment conditions are listed in Table 3. Compared with NASA icing tunnel experiment, the temperature is lower, and the air velocity and liquid water content are larger, which would lead to a severer ice protection condition.
The heating sequence is shown in Figure 10, with a cycle of 125 s. The heat blade (heating pad 4, at the leading edge) is kept activated during the whole cycle; heaters 3 and 5 are activated during 75–100 s, and other heaters are activated during 100–125 s. The ice protection structure model and the materials are the same with those in Section 4.2.
The convective heat transfer coefficient distribution is shown in Figure 11(a), and the contour of droplet collection efficiency is shown in Figure 11(b). The convective heat transfer intensity reaches its peak at the stagnation point and drops rapidly as it moves backwards due to the development of the boundary layer. Similarly, the droplet collection efficiency reaches its maximum at the stagnation point, and the value decreases as it moves backwards. There is no droplet impingement at the rear part.
(a) Convective heat transfer coefficient (W/m2K)
(b) Local collection efficiency
Figures 12 and 13 show the solid temperature distribution of the cross section and the ice shapes at t = 75 s, 100 s, and 125 s, which are the end time of the three stages. During the first stage (0–75 s), the structure temperature at the leading edge increases due to the activated heat blade, and the temperature of the heat blade reaches 289 K as is shown in Figure 12(a). The temperature of the adjacent area also increases due to heat conduction, while the value is lower than that of the heat blade and the heated area is limited, because the thermal conductivity of the structure materials is low. The water remains at a liquid state over the heat blade, as shown in Figure 13(a), and ice forms downstream. Compared with the result of Section 4.2, the ice accretion area and the ice thickness are larger in this case. All the ice protection area is covered by ice, because the air velocity and liquid water content are larger and the temperature is lower in this case. Later in the second stage (75–100 s), heaters 3 and 5 are turned on, ice starts to melt, and solid temperature increases. At 100 s, ice over the rear part of heaters 3 and 5 melts completely, and the surface temperature exceeds 273.15 K as shown in Figures 12(b) and 13(b). While there is ice left in the front of heaters 3 and 5, the accrete ice thickness is large there and the power is insufficient for all the accrete ice to melt. The ice continues to form on the nonheating surface area. In the third stage (100–125 s), the temperature at heaters 3 and 5 drops rapidly when the heaters are turned off, and runback water freezes over heaters 3 and 5, as shown in Figures 12(c) and 13(c). Ice over heaters 1, 2, 6, and 7 melts, and the temperature there increases to a high level since the convective heat dissipation there is relatively weak.
(a) t = 75 s
(b) t = 100 s
(c) t = 125 s
(a) t = 75 s
(b) t = 100 s
(c) t = 125 s
At the end of the second heating stage (t = 100 s), there is ice left unmelted in the front of the surface over heaters 3 and 5, which would lead to ice accretion on that area. According the discussion above, we know that the deicing power needed in the front is larger than that in the rear part. Therefore, the heating sequence is altered by extending the second stage to 30 seconds and shorten the third stage to 20 s, as shown in Figure 14. The total time of a cycle remains the same, while the power needed for a cycle is reduced, because the area of heaters 1, 2, 6, and 7 is larger. The simulated ice shape with the altered heating sequence is shown in Figure 15. Figure 15(a) shows that after the second stage the surface over heaters 3 and 5 is clean with all accrete ice melted. Figure 15(b) shows that after the whole deicing cycle, most part of the surface is clean, and only very slight ice accretion occurs at the location near the heat blade, which would be removed in the next cycle.
(a) t = 105 s
(b) t = 125 s
5. Conclusions
Based on the water film runback dynamics and energy balance theory, an unsteady conjugate heat transfer model for electrothermal deicing is established, and a loosely coupled solution method is developed. The model is applied in the simulation of deicing process, and the conclusions are as follows: (1)In-flight deicing process is very complex due to factors such as droplet impingement and water film runback. The present model is capable of simulating the in-flight deicing process. Simulation following an icing tunnel experiment has been conducted to validate the present model, and the results show good agreement.(2)The environment conditions would strongly affect the solid temperature distribution, the water film runback, and phase transition. A larger velocity or liquid water content would correspond to a larger runback icing range, and a higher power of longer heating duration is needed to perform the deicing process. Water remains at a liquid state over the heat blade, and ice forms on the surface where the heating power is insufficient.(3)Heating sequence is a key factor for the deicing performance. A proper heating sequence not only leads to a better deicing performance, but also saves energy. The optimization of heat sequence can be conducted by means of numerical simulation.(4)However, there are several factors not yet considered, including the ice shedding mechanism, the contact thermal resistance between multilayer materials, and the anisotropy of material properties.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work was supported by the National Natural Science Foundation of China (Grant no. 51206008).