Numerical Simulation of an Airfoil Electrothermal-Deicing-System in the Framework of a Coupled Moving-Boundary Method

A numerical method for the analysis of the electrothermal deicing system for an airfoil is developed taking into account mass and heat exchange at the moving boundary that separates the water film created due to droplet impingement and the ice accretion region. The method relies on a Eulerian approach (used to capture droplet dynamics) and an unsteady heat transfer model (specifically conceived for a multilayer electrothermal problem on the basis of the enthalpy theory and a phase-change correction approach). Through application of the continuous boundary condition for temperature and heat flux at the coupled movingboundary, several simulations of ice accretion, melting and shedding, runback water flow and refreezing phenomena during the electrothermal deicing process are conducted. Finally, the results are verified via comparison with experimental data. A rich set of data concerning the dynamic evolution of the distribution of surface temperature, water film height and ice shape is presented and critically discussed.


Introduction
Ice accretion may occur on the components of aircraft when the aircraft is flying through the clouds that are rich in supercooled droplets, which will threaten the flight safety of the aircraft [1,2]. Irregular ice accretion on the airfoil could destroy the aerodynamic configuration of the aircraft and cause earlier separation of airflow, which will lead to the decrease of lift-drag ratio and stalling angle of attack (AOA) [3]. Therefore, aircrafts are required to be equipped with anti-icing/deicing system to protect them from the risk of ice accretion [4,5]. The heat energy in the deicing system is usually generated by the engine bleed gas or electrothermal pads to remove the ice accretion on the surface of aircraft [6]. Because the electrothermal deicing system is convenient to control the heating law and improve the deicing performance, it has been widely applied in the helicopter blade, aircraft tail and some aircraft wing nowadays.
The early research about the electrothermal deicing system focus on the numerical simulation of heat transfer characteristic for multilayer electrothermal structure. At the beginning, the theory of heat transfer and phase change for airfoil structure and ice layer was preliminarily established for one-dimensional multilayer electrothermal pads [7][8][9]. Later, Chao [10], Wright [11] extended the simulation methods to two-dimensional rectangular deicing pads with the finite difference method. On this basis, Masiulaniec [12] performed electrothermal deicing simulation for curved configuration by applying coordinate transformation to the two-dimensional heat transfer equations and boundary conditions. Then Yaslik et al. [13] simulated the 3-D deicing process with the application of Douglas method for the discretization of heat transfer equations. Huang et al. [14] adopted an improved finite element method to simulate the heat transfer and phase change process of complex 3-D deicing cases.
There are abundant research results about the heat transfer and phase change process in the multilayer electrothermal structure and ice layer, which ignore the effect of droplet impingement and ice accretion on deicing characteristic. The related studies of the coupled ice accretion and heat transfer characteristics in the deicing process are quite rare. On the other hand, the numerical simulation method for the phase change [15,16] and heat transfer phenomena [17][18][19] have been developed greatly in the past few decades, which allows researchers to conduct further advancement on aircraft deicing issues.
Wright et al. [20] conducted experimental research for the electrothermal deicing process for NACA0012 airfoil with the heat flux term of ice accretion process calculated by Messinger's model [21], which became a part of LEWICE software. However, the ice shape transformation approach is needed in Wright's model to generate the ice shape around the airfoil, which will be time-consuming for complex configuration. Besides, the retained water is ignored in Messinger's model, which will cause computation deviation in the amount of runback water. Al-Khalil et al. [22] developed the simulation method of runback water, water film breakup and rivulets on the aircraft anti-icing surface and integrated it into ANTICE software. Reid et al. [23] simulated the unsteady conjugate heat transfer characteristics of electrothermal deicing system in FENSAP-ICE software. The results of the heat transfer module are in good agreement with analytical solution and the temperature response in the airfoil is in same tendency with experimental data, with the existence of some little discrepancy in the first cycle of periodic heating process as well as in the cooling stage of remaining cycle. Nevertheless, the water film surface and the upper boundary does not coincide and the surface projection approach that computes the heat transfer between these domains is not discussed in detail. Silva et al. [24] established the mathematical models for gaseous flow boundary, momentum or thermal boundary layer, water film flow and airfoil solid surface through the domain division approach to predict the surface temperature distribution and the mass flow rate of runback water. Due to only the anti-icing situation is considered in Silva's model, the mass rate of ice accretion and the heat transfer characteristic in the ice layer are ignored. Harireche et al. [25] solved the Myers model and unsteady heat transfer model with the explicit finite volume approach to simulate the runback water and ice accretion characteristic in anti-icing/deicing process, the results are in good agreement with the other codes and the corresponding code is implemented in ICECREMO2 software. However, only the water film flow and temperature distribution on the flat cases are considered and the simulation of anti-icing/deicing characteristic on the curved airfoil configuration is not conducted. The heat transfer in the ice layer is also simplified with linear assumption of temperature, which is not precise enough as the thickness of ice increasing in the deicing cases.
Nowadays, the electrothermal deicing technology is widely applied in current aircraft design, and the computational method and technology for fluid dynamics and thermodynamics have been applied in more and more research fields [26][27][28][29]. Thus, the ability for the comprehensively numerical simulation of aircraft electrothermal deicing characteristic is required. Through the review of current researches, it could be found that the simulation method of the coupled mass and heat transfer process for the aircraft deicing system still needs to be developed and improved, especially for the research of coupled method among external flow field, runback water film, ice layer and airfoil structure.
In current research, a novel coupled moving-boundary method is presented to simulate the coupled mass and heat transfer characteristic on the moving boundary caused by the ice growth or ice shedding in the deicing process, which is often simplified to a static coupled boundary located on the airfoil skin in most researches for the convenience of heat transfer simulation. The ice recognition method and corresponding boundary heat transfer approach are established to calculate the temperature distribution of the variable ice shape with static structure grid. The water film flow and ice accretion model on the upper surface of ice and airfoil are established to precisely simulate the physical characteristics of runback water and ice growth in deicing process. Finally, the complete simulation for electrothermal deicing process with the consideration of water film flow, ice accretion, ice melting and shedding as well as refreezing phenomena is achieved. The current models are validated by comparing the results with experimental data for the deicing cases, and the temperature response, water film flow and ice accretion characteristic are analyzed to study the deicing performance.

Mathematical Model
Before solving the coupled water film flow and heat transfer model, the dependent parameters (such as the local droplet collection efficiency and the surface convective heat transfer coefficient on the airfoil surface) are needed by solving the two-phase flow field of air and water droplets.

Air Flow Field Model
Because the liquid water content (LWC) of supercooled droplets in the air is relatively low, the influence of droplets on the air can be ignored, so that the air flow field and the droplets flow field can be decoupled. The Reynolds Averaged Navier-Stokes (RANS) equations are used to solve the air flow field. Based on the Boussinesq eddy viscosity assumption [30], the relationship of Reynolds stress with turbulent viscosity coefficient and average velocity gradient can be established, while the turbulent viscosity coefficient can be determined by the suitable turbulent model. Finally, the density, velocity, pressure and temperature distribution of air flow field can be derived by the solution of RANS equations as follows [31]: @v a @y @u a @y þ @v a @x @v a @x þ @u a @y 4 3 @v a @y À 2 3 @u a @x where q a ,V a ¼ u a ; v a ð Þ, e a , T a and p a correspond to the density, velocity vector, total energy, temperature and pressure of air. l L , l T are the laminar and turbulent viscosity coefficients, while k L , k T represent the laminar and turbulent thermal conductivity, respectively. s is the RANS stress matrix.

Droplets Flow Field Model
There are two approaches for the simulation of droplets motion. One is the discrete phase approach [32,33] which treats droplets as the group of single particle. Another is the continuous phase approach [34] that assumes the droplets have continuous physical properties, so the Eulerian method can be applied for the solution of droplets motion. In current research, the Eulerian two-phase flow method is used to calculate the droplets flow field. The governing equations of droplets flow field can be expressed as where q d ¼ a d q w corresponds to the apparent density of droplets, a d , q w are the volume density of droplets and the density of liquid water, respectively. V d ¼ ½u d ; v d T , C D , D d and F g ¼ ½0; g 1 À q a =q d ð Þ T are the velocity vector, drag coefficient, diameter and gravity vector of droplets. K d and Re d represent the droplet inertial parameter and Reynolds number based on the relative droplet velocity, respectively.
After the solution of the droplet flow field converges, the apparent density and velocity distribution of droplets on the airfoil surface can be obtained. Then the droplets impact characteristics on the airfoil surface are represented by the local droplet collection efficiency b, which is calculated as follows: where q d;1 , V d;1 are the apparent density and velocity of the freestream of droplets, respectively. And n is the normal vector on airfoil surface.

Coupled Water Film Flow and Heat Transfer Model
There are complex mass and heat transfer phenomena during the deicing process on the airfoil surface, such as ice accretion, melting and shedding, runback water and refreezing, as shown in the Fig. 1.The exist ice indicates the ice that has accreted on the airfoil surface before the current moment, whose heat transfer characteristic is solved by the solid heat transfer model. The accreted ice indicates the newly formed ice which is produced by the water film flow in the current time step. The solution data in the two models are exchanged through the surface temperature and heat flux on the coupled boundary. As ice grows and sheds, the boundary of exist ice moves and the external flow field is correspondingly changed. A dynamic ice recognition approach and coupled moving-boundary method are proposed to consider the effect of moving boundary of ice accretion in the deicing process. The following assumptions are applied for the simplification of the calculation: (1) The thickness and reduced Reynolds number of water film are sufficient small so the lubrication theory can be used [35]. (2) The water film flow is continuous and ignore the breakup of water film as well as the formation of beads and rivulet. (3) The thermal properties of airfoil structure material, water and ice are independent with temperature, the density of ice remains unchanged during melting process. (4) The liquid water formed by the ice melting is removed from the control volume and is not counted into the amount of runback water. (5) The mechanical characteristics in the ice layer is ignored, which means the ice shedding process is only subjected to the effect of temperature distribution.

Water Film Flow Model
When the surface temperature is relatively high, such as in the glaze ice accretion situation or in the deicing process, the supercooled droplets that impact on the surface will form water film. The governing equations about the mass and heat transfer of the water film flow on the substrate need to be established in order to determine the water film distribution and ice shape. The control volume (CV) of the water film flow model is shown in Fig. 2.
Generally, the thickness of the water film on the aircraft surface (typical thickness about 10 −6~1 0 −4 m) is much smaller than the chord length of the airfoil. Thus, the lubrication theory can be used to simplify the governing equations of thin water film. Through the magnitude analysis of the continuous equation and momentum equation and neglect the less important terms, the mass governing equation for the water film flow [36] is derived as where l w is the dynamic viscosity of water and g s , A s represent the gravity and air shear stress component along the tangential direction of surface curve, respectively. The _ m imp , _ m ev are calculated by where h cv , c p;a are the convective heat transfer coefficient and the air specific heat capacity. p s , p m represent the saturated vapor pressure and the average air pressure between freestream and airfoil surface. T aw is the temperature at water film-air interface.
The energy equation in the water film can be expressed as where T wf , p wf , s wf , V wf and k wf are the temperature, pressure, shear stress matrix, velocity vector and thermal conductivity of water film. Due to the thickness of water film is rather small, the transient term, convection term, stress work term and lateral conduction can be neglected, then the energy equation of thin water film is simplified as In the thin ice layer, the similar energy equation can be derived as where T acÀice represents the temperature in the accreted ice. For the convenient of description, the symbol T ðzÞ is used to represent T wf and T acÀice . According to the air temperature and substrate temperature, three kinds of water film flow and ice accretion types exist on the substrate, as shown in Fig. 3. For the type (a), the boundary condition at the accreted ice-substrate interface is where T s and _ Q s represent the temperature and the heat flux of substrate. At the water film-accreted ice interface, the boundary condition is the Stefan condition, which can be expressed as where T f , L f , k i , k w are the temperature of freezing point, the latent heat of ice, thermal conductivity of ice and water, respectively. By substituting the temperature conditions at z ¼ 0 and z ¼ b into the Eq. (17), the temperature distribution in the accreted ice can be derived as At the water film-air interface, the boundary condition is Q ev are the droplet impact kinetic energy, droplet cooling effect on the water film, aerodynamic heating effect, convective energy between water film and air as well as the evaporative energy. The energy terms are calculated by where R, V e , L ev are the temperature recovery factor, air velocity at the edge of boundary layer and the latent heat of evaporation. Similarly, the temperature distribution in the water film is computed as Substituting the temperature distribution in the water film and accreted ice into the heat flux boundary condition at the water film-accreted ice interface, the governing equation of the ice accretion height can be obtained as where q 0 and q 1 can be calculated by For the type (b), the impacting droplets freeze immediately without the forming of water film. The ice height and the temperature distribution can be computed as where _ m su , T ai represent the mass flux of sublimation and the temperature at accreted ice-air interface.
For the type (c), because the surface temperature of substrate exceeds the freezing point, the exist ice gradually melts and subsequent impacting droplets keep liquid on the surface. The water film height and the temperature distribution can be computed as

Unsteady Heat Transfer Model
The heat transfer domain of airfoil is the multilayer structure with the exist ice layer adhesive on the top of the airfoil skin. For the airfoil structure, the heat transfer process is governed by the unsteady solid heat transfer equation [37], while the simulation of the heat transfer process in exist ice layer is a little difficult because of the existence of phase change. Though the application of enthalpy theory [8], the unified heat transfer equation can be established for ice layer, and the phase change interface can be determined by the enthalpy distribution. The unsteady heat transfer equation based on the enthalpy theory can be written as [38]: where H is the total enthalpy. k ly and T ly is the thermal conductivity and temperature of each layer. _ Q h is the heat flux generated by the heater pads.
For airfoil structure, the total enthalpy is calculated by the following formula: where q ly , c p;ly is density and specific heat capacity of each layer in the airfoil structure.
For the ice layer, as temperature rises, the phase state of ice will change among ice, ice-water mixture and water. For the convenience of numerical simulation, the enthalpy theory assumes that the phase change occurs in a small temperature range DT m and the ratio of ice and water in the mixture changes linearly with temperature. The schematic diagrams of this assumption are shown in Figs. 4 and 5. In Fig. 4, the movement of phase change boundary in control volume (CV i ) which is in the ice-water mixture state is linear along the direction of phase change, i.e., the location of phase change boundary xðtÞ is a linear function. Thus, with the above assumptions the proportion of liquid water f in CV i can be expressed as: where T sm is the phase change temperature of ice and ice-water mixture state. The enthalpy value of CV i is calculated by: where c i is the specific heat capacity of ice.
According to the above assumption, the total enthalpy-temperature relationship in the ice layer is plotted in Fig. 5. Therefore, the total enthalpy in the ice layer can be calculated as follows: where T lm is the phase change temperature of ice-water mixture state and liquid water. In current research, the T sm is equal to 273.15 K and DT m is assumed to be 0.01 K.
For the phase change process, the material property and heat transfer characteristic are changing discontinuously. Thereby the temperature calculated by the enthalpy theory needs to be modified with phase change correction approach. Assumed that material is transformed from the phase state 1 to the  phase state 2 and the calculated enthalpy at the next time by the heat transfer equation isH nþ1 , the modified temperature T nþ1 is computed by where T 12 , H 12 are the phase change temperature and total enthalpy of phase state 1 and 2.

Coupled Moving-Boundary Method
In the process of deicing simulation, the coupled boundary between the water film domain and solid heat transfer domain moves with the variation of ice shape caused by the ice growth and shedding, which leads to three main difficulties in numerical simulation: The heat transfer calculation in the variable ice layer, the coupled method for water film flow and solid heat transfer and the simulation of changing external flow. The coupled moving-boundary method is proposed to deal with these problems.
The temperature distribution in the ice layer should be calculated with the unsteady heat transfer PDE (Partial Difference Equation) rather than the approximate function as the ice thickness increases, which requires the determination of solid heat transfer domain of variable ice layer. The ice recognition method is established for the variable layer domain based on the approach of piecewise linear approximation on the structure grid. Firstly, a background grid with small size is established on the top of airfoil structure. The points of intersection between current ice shape function f ice ðn; gÞ and background grid can be determined with Lagrangian interpolation approach. Then the ice shape between the adjacent points of intersection are approximated with linear configuration, as shown in Fig. 6. When the variation of ice height caused by ice growth and shedding exceeds the recognition standard (Db > b cell ), the solution domain of the ice layer is updated with the ice recognition again. Compared with the dynamic grid technology, this method has less computational cost and avoids producing grid distortion, which is rather common when the rough ice shape is accreted on the airfoil.
The temperature at the coupled boundary T s is interpolated with temperature of the cell on the top of ice layer (T p ), as shown in Fig. 6.
Then the water film flow governing equations are solved with T s to obtain the ice height, water film height as well as the temperature distribution. The heat flux q n ¼ ðq n ; q g Þ at the coupled boundary which is derived by the temperature distribution in the water film flow model is imposed as the boundary Figure 6: The diagram of ice recognition method and coupled variables calculation on the boundary condition of solid heat transfer domain. After an iterative step, the temperature in the water film and solid heat transfer domain are obtained. Because there exists a delay of data exchange between the water film domain and heat transfer domain, the iterative process is required before the convergence of temperature and heat flux on the coupled boundary. The maximum temperature of water film and solid heat transfer domain is chosen as the convergence criterion. That means if the difference between the maximum temperature at current iterative step and at previous iterative step (Max|T k+1 -T k |) is less than the sufficiently small value E the converged solutions at current time step are obtained. Assumed that ice grows along the normal direction of recognized coupled boundary, i.e. n s , the ice shape at next time can be obtained with the converged ice height: where r n is the ice shape position vector in Cartesian coordinate system at current time. Generally, the grids of water film flow model and unsteady solid heat transfer model are in different sizes, thus the variables at coupled boundary need to be interpolated with Lagrangian interpolated approach.
In actual flight, ice accretion will melt as the heating of the airfoil skin, and shed under the action of aerodynamic force or the inertial force. The mechanics characteristic of ice shedding is rather complicated, so ice shedding process is simplified by only considering the effect of ice layer melting in current simulation. That means if the phase-change line in the ice layer reaches at the top of the ice shape or fully exceeds the range of ice accretion, the ice shedding occurs. The external flow fields vary with the ice growth and shedding, affecting the water film flow and ice accretion characteristics. The dynamic process of external flow is approximately simulated with quasi-steady multistep approach to decrease the computational cost. When the variation of ice thickness Db e exceeds the updated standard of external flow fields Db crt , the grid of external flow fields is re-sized and the solutions of external flow fields are updated.
Above all, the flow chart of the coupled water film flow and heat transfer model based on the coupled moving-boundary is shown in the Fig. 7.

Numerical Discretization and Calculation Method
The air flow field is calculated with the CFD solver FLUENT [39]. The governing equations of air flow field are solved using Finite Volume Method (FVM), the convection term is discretized with the second-order upwind scheme, the pressure velocity coupling is dealt with the SIMPLE scheme [40], and the turbulent viscosity coefficient is calculated using the Spalart-Allmaras turbulence model [41]. After the solution of the air flow field converges, the density, velocity and pressure distribution on the grid are derived to solve the droplets flow field, while the convective heat transfer characteristic is obtained to solve the coupled water film flow and heat transfer model.
The droplets flow field is calculated with FORTRAN codes. The governing equations of droplets flow field are discretized using the FVM, the integral form of the governing equations on the control volume can be expressed as: where V cv and S are the area and side length of the control volume, n ¼ n x ; n y À Á is the normal vector on the side of the control volume. The discretization of above equations are accomplished with the ADI (Alternating Direction Implicit) method [31]. After the solution of the droplets flow field converges, the droplets impact characteristic on the airfoil surface is input into the water film flow model for the simulation of deicing process.
The water film flow and ice accretion model is computed with FORTRAN codes. The governing equation of water film flow is discretized with Finite Difference Method (FDM). The time term is discretized with first-order Euler implicit scheme, while the convection term is discretized with water film flux on the side of cell Q nþ1 kAE1=2 , as shown below:  where Dt and Ds k are the time step and space step. The water film flux on the side of cell is calculated with hybrid upwind and Lax-Wendroff scheme: The mass flow rate of ice accretion is computed with the explicit scheme, i.e., After the discretization, the algebraic linear equations can be solved using TDMA (Tridiagonal Matrix Algorithm) method [31].
The unsteady heat transfer model is calculated with FORTRAN codes. The FVM is used for the discretization of unsteady heat transfer equation. The integral form of the governing equations on the control volume can be expressed as: where Q ¼ k ly @T ly @x ; k ly @T ly @y . The normal heat flux at the side of the control volume is calculated with Jacobian transformation as follows: Q Á n ð Þ S ¼ k ly @T ly @n n x n x þ n y n y À Á þ @T ly @g g x n x þ g y n y À Á ! S where n x , n y , g x , g y are the Jacobian transformation coefficients which can be computed by the node coordinates of grid. The time term is discretized by Euler implicit format. After the discretization, the multi-diagonal sparse coefficient matrix linear equations are obtained and the LU-SGS (Lower-Upper Symmetric Gauss-Seidel) method which has the advantages of good stability, simple format, and high efficiency is used for solution of discretized linear equations [31].

Grid Independence Test
The deicing case of UH-1H helicopter blade in the literature [42] is used for the grid independence test. The experimental objective is UH-1H helicopter blade, whose airfoil is NACA0012 and has a chord of 1 m. The airfoil is a 7-layer composite structure, and the adjacent layers are bonded with adhesive, as shown in Fig. 8. In the zone close to the stagnation point of the leading edge, the noseblock is added to the substrate, which causes the structure is slightly different from the structure at the side of the blade. The thickness and material property of each layer are shown in the Tab. 1, where a ¼ k ly /ρ ly c p , ly is the thermal diffusivity. There are eight heater pads in the airfoil, each one has the length of 25.4 mm, which are also shown in Fig. 8. The experimental conditions are listed in Tab. 2.
The grids used in current research include the grid of external flow field as well as the grid of internal heat transfer domain, as shown in Figs. 9 and 10. The air/droplets flow field are simulated with external structure grid in Fig. 9. The quasi-steady assumption is applied for the calculation of the air/droplets flow field solutions of iced airfoil. That means if the variation of ice thickness Db e is less than the critical ice thickness Db crt , the air/droplets flow flied is simplified as steady flow. After Db e exceeds Db crt , the grid of external flow field is updated with the new configuration of iced airfoil. The data of unsteady heat transfer     Fig. 10, and the ice accretion calculated by the water film flow model is recognized on the background grid introduced in Section 2.3 (which is not shown in Fig. 10).
To analyze the mesh sensitivity, three kinds of different grid sizes are chosen for external flow field and internal heat transfer domain, respectively. The sum of local droplet collection efficiency b of the droplets flow field as well as the maximum temperature of the unsteady heat transfer domain at 20 s are output to analyze the deviation. The results are shown in Tabs. 3 and 4. It is noted that the deviation of the sum of b between Grid 3a and Grid 2a is obviously smaller than the deviation between Grid 1a and Grid 2a. Thus, the Grid 3a is chosen as the grid of external flow field. As shown in Tab. 4, the maximum temperature of the Grid 3b is in little difference with Grid 2b when compared to the accuracy improvement of results between Grid 1b and Grid 2b. Considering the computational time and costs, the size distribution of Grid 2b is used for the simulation of unsteady heat transfer process.

Validation of Phase Change Model
The 2D deicing test case about a rectangular multilayer domain from NASA/DRA/ONERA [43] is used to verify the current solid heat transfer and phase change method. The structure is shown in the Fig. 11 and the material properties are shown in Tab. 5.   The simulated results of temperature at ice-titanium interface by current model are compared with the data of NASA/DRA/ONERA codes, which are shown in Fig. 13. The temperature rises above the melting  Due to the effect of phase change latent heat, the temperature of liquid water will remain at freezing point for a period. This phenomenon is successfully simulated with current model, which is important for the evaluation of deicing and refreezing characteristics. It can be seen that the current simulation results are in good agreement with the calculation results of other codes, which verifies the current heat transfer model based on the enthalpy theory and phase change correction approach.

Deicing Case of a NACA0012 Airfoil
The coupled deicing case from experiment [42] is used for the validation of coupled water film flow and heat transfer model, which includes the phenomena of heat transfer in the airfoil and ice layer, ice melting and shedding, runback water flow and refreezing. The airfoil structure, material properties and the experimental conditions have been introduced in the Section 3.2. In the experiment, the blade was exposed in the ice wind tunnel for 5 minutes before heater pads are activated. Then the experimental ice shape is input into the codes to conduct the coupled water film flow and heat transfer simulation. Two different heating sequences are selected for the simulation of simultaneous heating process and asynchronous heating process. In heating sequence 1, all eight heater pads are turned on for 20 s and turned off for 60 s. In heating sequence 2, the heat pads 5~7 works for 20 s and then rest, while heater pads 3, 4, 8 works during 10 s~15 s and heater pads 1, 2 works during 15 s~20 s.
Firstly, the heating sequence 1 is used for current simulation. The temperature response in the position of heater 7 is compared with the experimental results and Leffel's simulation results to verify the models, which is shown in Fig. 14.
The temperature response of current models at the AB shield-Ice interface is consistent with the experimental result. When the temperature of ice exceeds the freezing point, the material property of ice layer transforms from ice to water due to the phase change phenomenon, leading to the variation of slope of the temperature response in Fig. 14. The temperature at the heater-insulator interface is in acceptable deviation with experimental result. For further analysis of the temperature difference, the standard deviation of temperature response among current simulation, experimental data and Leffel's simulation are computed, as shown in Tab. 6. At heater-insulator interface, due to the lack of consideration of ice shedding and runback water the Leffel's result overpredicted the temperature near the peak when compared with the experimental data, while current result underpredicted the maximum temperature. The reason may be attributed to the neglect of the material properties change related with temperature. The standard deviation of current simulation and experimental data is slight better than Leffel's simulation and experimental data, but the difference of both results is acceptable. It should be stated that the temperature response in the structure is affected by many factors, such as the variation of material properties caused by temperature change, boundary condition mutation caused by the ice shedding and runback water flow, as mentioned in literature [44], which cannot be totally precisely considered in the process of current simulation. Due to the deicing characteristic is mainly dependent on the surface temperature, the temperature response at abrasion shield-ice interface is a more important parameter for the evaluation of simulation result. The temperature response at abrasion shield-ice interface of current result and Leffel's result are both close to the experimental data while the standard derivation of current result is slightly smaller than Leffel's result. In generally, the results of current simulation are in acceptable deviation with the experimental data and can be used for further analysis of coupled deicing process.
The temperature distribution, the water film height and ice accretion height are calculated to analyze the deicing characteristics. The temperature contours of the airfoil in the heating stage of 0 s-20 s are shown in Fig. 15.  It can be seen that the ice layer still exists after heating for 5 s, and the top surface of the ice layer and the interface between the ice layer and abrasion shield are warmed up. At 12.5 s, the initial ice accretion has shed, then the runback water flows downstream and forms ice accretion in the unprotected zone. At 20 s, the range and height of ice accretion in the unprotected zone increase, and no ice grows on the surface of the heating zone. To further analyze the above process, the surface temperature, water film height and ice height at three moments are compared, as shown in the Figs. 16-18.  At 5 s, the temperature above the heater 5-7 is below the freezing point, which means that the ice layer doesn't shed yet. It can be seen that the surface temperature of the iced zone is slightly lower than the temperature of the surrounding clear zone. It is because the heat flux transferred by the ice layer is larger compared to the convective heat transfer process. At this time, the coupled boundary is located on the top of icing layer and there exists water film on the icing surface, which means the top surface of ice layer is in the state of phase change. Thus, the temperature of the ice surface is at the freezing point, which is higher than the ambient temperature, leading to the rise of temperature in the area near the outside surface of ice layer. The ice height at 5 s has also increased compared to the initial time, as shown in Fig. 18b. At 12.5 s, the initial ice layer has shed. It can be seen in Fig. 16 that the surface temperature of the heating zones is higher than the freezing point. This means that the supercooled droplets impacting on the surface of heating zone will form runback water and move downstream until reaching the unprotected zone. After leaving the heating zone, the surface temperature decreases below the freezing point. Thus, the temperature of runback water drops and starts to freeze. The main reason for the formation of two ice peaks at the outside of heating zone is that the runback water will freeze from the beginning of unprotected zone and stop freezing when the water film height decreases to zero. As this process proceeding, two ice peaks finally form at the outside of the heating zone in Fig. 18a. Observing the water film height and the ice height at 12.5 s, it can be found that after leaving the heating zone, the water film height gradually decreases because of the formation of ice accretion until all runback water freezes. At 20 s, the range and height of ice accretion in the unprotected zone continue to increase. The temperature contour at 20 s shows that the temperature at the top of the ice layer has increased due to the effect of heating. Thus, the freezing rate of runback water in these areas will reduce and the runback water will continue to move downstream and refreeze, which expands the range of ice accretion.
The heating pads are all switched off at 20 s and then the system would be cooled by the ambient temperature. The temperature contours of the airfoil in the cooling process is shown in Fig. 19. It can be found on the temperature contour of 30 s that after the heaters are turned off, the system gradually cools down, but ice accretion still doesn't occur on the protected zone. On the contrary, the thickness and range of ice accretion in the unprotected zone further increases. At 55 s, the ice shape in the unprotected zone remains basically unchanged, while the ice accretion occurs on the leading edge of the airfoil. When one The process of ice accretion has gone through three stages after the heaters are turned off. At the beginning of the cooling process, the surface temperature in the heating zone is still above the freezing point, thus the runback water will still move downstream and form ice in the unprotected zone, as shown by the results of 30 s in Figs. [20][21][22]. Subsequently, as the cooling process continues, the surface temperature in the heating zone gradually decreases below the freezing point. Therefore, the ice accretion emerges on the surface of heating zone and the water film gradually shrinks inward. Observing the ice height at 55 s in Fig. 22, it can be seen that there exists a small amount of ice accretion in the middle area between the stagnation point and the edge of heating zone, which is formed during the shrinkage  Fig. 20 shows the surface temperature in the heating zone basically reduces below the freezing point and Fig. 21 shows the water film only exists in a small area close to the stagnation point at 55 s. Therefore, the impacting supercooled droplets only freezes on the surface near the stagnation point and no longer moves downstream at the third stage. The reason why the surface temperature near the stagnation point is higher than surrounding area at the third stage is that the existence of water film on the top surface of ice layer keeps the temperature of coupled boundary higher than the dry zone.  From above results, it can be found that turning on all the heaters at the same time can achieve good deicing performance at the beginning. But if the surface heat flux generated by heater pads is not enough to completely evaporate the impacting supercooled droplets, the runback water will leave the heating zone and form ice ridge in the unprotected zone, which is difficult to be removed. As the deicing process continues, the height of ice ridge increases and will greatly affect the aerodynamic characteristics of the airfoil. A more reasonable way is to use different period for each heater, so that the runback water can periodically freeze and shed before leaving the heating zone, which will achieve a better deicing effect at the refreezing stage.
The heating sequence 2 is selected to verify the effect of asynchronous heating. Other conditions keep same with the case of heating sequence 1. The comparisons of distribution of surface temperature, ice height and water film height between sequences 1 and 2 are shown in Fig. 23.
In the heating process, the surface temperature above the heater pads 3, 4, 8 is close to the ambient temperature before 10 s because these heater pads are deactivated, as shown in Fig. 23a. Therefore, after the initial ice accretion sheds the runback water will move downstream and form ice accretion on the surface above heater pads 3, 4, 8, as shown in Figs. 23b and 23c. After the heater pads 3, 4, 8 are fired at 10 s, the surface temperature rises and gradually exceeds the freezing point. Then, the ice accretion on the lower surface totally sheds while the ice accretion on the upper surface does not shed, as shown in Fig. 23b. The reason is that in the upper surface the heater pads are not equipped on all iced zone, only the ice at the ice-abrasion shield above the heating zone melts and the remaining part still attaches on the unprotected zone. So the ice accretion on the upper surface will not fully shed. Thus, the ice accretion can only partially shed after the ice layer above the heating zone totally melts. After the heater pads 1, 2 are activated, the ice accretion on the lower surface sheds again. Then the runback water flows to the unprotected zone and forms ice accretion. Due to most ice accretion on the lower surface has been removed by ice shedding in asynchronous heating process, the height and range of ice ridge on the lower surface is obviously less than on the upper surface at the end of heating process, as shown in Fig. 23d.
In the cooling process, the surface temperature on the lower surface in heating sequence 2 is lower than heating sequence 1 at 30 s because that the heater pads 1~4 work less time in sequence 2. Thus, thicker ice accretion can be observed on the surface above heater pads 1~4 during the shrinkage stage of the water film. It is also noticed in Fig. 24c that the end points of water film are more forward to the leading edge on both side of airfoil surface at 30 s, which shows that the less activation time for the heater pads on the backward location would be benefit for the reduction of the amount of runback water on the unprotected zone. After one period ends, the amount of ice accretion on the lower surface in heating sequence 2 is less than in sequence 1 with less activation time of heater pads 1~4, which means that the deicing performance is improved by the optimization of the activation period of different heater pads.

Conclusion
In the current work, a novel coupled moving-boundary method is proposed for the simulation of dynamic electrothermal deicing process on the airfoil. Though the establishment of the ice recognition approach for variable ice shape and heat transfer condition at recognized coupled boundary, the coupled water film flow and heat transfer model is solved and the phenomena such as ice growth, ice melting and shedding, runback water and refreezing on the deicing surface are successfully simulated. The following conclusions can be obtained with the validation cases: 1. The interface temperature of deicing case on the rectangular multilayer domain from current model is in good agreement with the results of NASA/DRA/ONERA codes, which verifies the heat transfer The simulation results of the coupled water film flow and heat transfer model are compared with the experimental data of UH-1H helicopter blade, the temperature response at ice-abrasion shield interface is in good agreement with experimental data, while the result at heater-insulator interface is in acceptable deviation. 3. Based on the coupled moving-boundary method, the flow and icing characteristic of water film on the variable ice surface as well as the heat transfer and phase change characteristic of the variable ice layer in deicing process could be dynamically simulated. 4. The results of deicing process clearly show that the simultaneous activation of all heater pads will lead to the formation of ice ridge in the unprotected zone. By proper optimization of the activation period of different heater pads to periodically remove the ice accretion, better deicing performance can be achieved with less energy consuming. Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.