Numerical Simulation of the Anti-Icing Performance of Electric Heaters for Icing on the NACA 0012 Airfoil

Abstract

Ice accretion is a phenomenon whereby super-cooled water droplets impinge and accrete on wall surfaces. It is well known that the icing may cause severe accidents via the deformation of airfoil shape and the shedding of the growing adhered ice. To prevent ice accretion, electro-thermal heaters have recently been implemented as a de- and anti-icing device for aircraft wings. In this study, an icing simulation method for a two-dimensional airfoil with a heating surface was developed by modifying the extended Messinger model. The main modification is the computation of heat transfer from the airfoil wall and the run-back water temperature achieved by the heater. A numerical simulation is conducted based on an Euler–Lagrange method: a flow field around the airfoil is computed by an Eulerian method and droplet trajectories are computed by a Lagrangian method. The wall temperature distribution was validated by experiment. The results of the numerical and practical experiments were in reasonable agreement. The ice shape and aerodynamic performance of a NACA 0012 airfoil with a heater on the leading-edge surface were computed. The heating area changed from 1% to 10% of the chord length with a four-degree angle of attack. The simulation results reveal that the lift coefficient varies significantly with the heating area: when the heating area was 1.0% of the chord length, the lift coefficient was improved by up to 15%, owing to the flow separation instigated by the ice edge; increasing the heating area, the lift coefficient deteriorated, because the suction peak on the suction surface was attenuated by the ice formed. When the heating area exceeded 4.0% of the chord length, the lift coefficient recovered by up to 4%, because the large ice near the heater vanished. In contrast, the drag coefficient gradually decreased as the heating area increased. The present simulation method using the modified extended Messinger model is more suitable for de-icing simulations of both rime and glaze ice conditions, because it reproduces the thin ice layer formed behind the heater due to the runback phenomenon.

1. Introduction

Ice accretion is a phenomenon whereby an ice layer is formed on a solid surface due to the adhesion of super-cooled water droplets and ice particles. This icing phenomenon has been observed in various industrial apparatus and it causes severe accidents during machine operation. In aircraft, the occurrence of icing affects many components and severely decreases aerodynamic performance and flight safety. For instance, icing on the wing surface changes the wing shape and surface roughness, which negatively affects aerodynamic performance. Moreover, icing on measuring instruments also increases the risk of aircraft accidents.

Generally, icing phenomena strongly depend on many factors. Between them, three factors, namely the liquid water content (LWC), ambient temperature, and droplet size, are considered important for predicting icing. Here, LWC represents the amount of condensed water per volume unit, whereas the medium volumetric diameter (MVD) represents the typical droplet size. The ambient temperature affects the ice characteristics. Thus, glazed and rime ice form under relatively high and low ambient temperature conditions, respectively. Rime ice refers to super-cooled droplets that impinge on the wing and freeze instantaneously, glaze ice refers to the water film created by an impinged droplet running back along the wing and freezing, and the runback phenomenon is the liquid film flows downstream without icing. The former and latter icing phenomena occur at low (less than −15 °C) and relatively higher temperature (more than −10 °C) conditions, respectively. Additionally, relatively large droplets, including super large droplets, are well known to exhibit complex droplet behavior, such as droplet splash and bounce [1,2].

With regard to accident prevention, many experimental and numerical studies have been conducted. From the viewpoint of financial and temporal cost, the development of a prediction method for icing is urgently needed. Therefore, many models for predicting the ice shape have been proposed. The Glen Research Center [3] developed the LEWICE code, which estimates the icing amount, position, and shape. Subsequently, the LEWICE code was expanded into the three-dimensional LEWICE3D code [4], which has been widely used to predict both aircraft icing and icing on the rotating blades [5]. Similar studies on the development of icing code have been conducted worldwide [6,7,8]. De- and anti-icing methods have been investigated not only for airfoils [9], but also wind turbines [10].

Notably, accidents have occurred 583 times from 1982 to 2000 [11]. Although the total number of icing accidents is relatively small, it is not negligible. Hence, the Federal Aviation Administration (FAA) [12] has provided regulations for avoiding icing during cruise flight. However, the decrease in aerodynamic performance by the icing has been reported, even under unregulated conditions [13]. Therefore, the development of anti- and de-icing devices is very important to ensure that aircraft safety is not compromised by icing.

To develop de- and anti-icing devices for aircraft wings, several techniques have been proposed and implemented in aircraft. Bleed air [14,15,16] prevents icing by heating up the wing surface using warm air supply from the compressor. However, the system and design are complicated. De-icer boots [17] deform the rubber on the wing surface by supplying compressed and decompressed air, but this is only effective for a specific ice thickness. Anti-freezing liquid [18,19] has also been used, although it is not environmentally efficient.

Recently, electric heaters have been introduced to aircraft, owing to their easy setting, environmental efficiency, and promotion of aircraft electrification. Al-Khalil et al. [20] conducted experimental and numerical investigations in order to assess the anti-icing performance of an electric heater on the NACA 0012 airfoil and reported the effect of different heating temperatures. Reid et al. [21] conducted numerical simulation for an electric-thermal de-icing system and obtained results that are in excellent agreement with experimental data. Additionally, they reported the temporal variation of surface temperature in the de-icing process. Bu et al. [22] conducted an icing simulation for an airfoil considering the energy conservation law with the heat flux from the wing surface and demonstrated the validity of their model by comparing the obtained results to the experimental results. Zhou et al. [23] conducted an experimental study and applied a dielectric barrier discharge (DBD) plasma actuator as a de- and anti-icing device. They confirmed that the heating effect of the DBD plasma actuator is effective under glaze-icing conditions, and even prevents icing close to the leading edge under rime-icing conditions. However, icing still occurs behind the heater, owing to water film runback. Mu et al. [24] suggested a mathematical model comprising water film runback dynamics, energy balance theory, and conjugate heat transfer model for an electrothermal heater. Their results are in good agreement with experimental results. Asaumi et al. [25] conducted an experimental study using a simple model to evaluate the anti-icing performance of NACA 0013 airfoil with a heater. They reported that the heater temperature threshold for achieving the anti-icing effect is 2–5 degrees Celsius. More recently, simulations of anti-icing by heating have been conducted using a commercial solver, not only for airfoils [26], but also for wing turbines [27,28].

As previously mentioned, many experimental and numerical simulation studies have been conducted on de- and anti-icing devices. The effectiveness of thermal heaters has been clarified, and the heater input has been optimized. However, optimization is required with regard to the heater area. In this study, a new model is proposed in order to predict the occurrence of icing with a heater. Two-dimensional icing simulations were conducted using the NACA 0012 airfoil to investigate the effect of the heating area on the aerodynamic performance. The extended Messinger model (EMM) was modified (MEMM) [29] to consider the heat transfer from the airfoil surface in order to simulate the ice growth. Through this study, the influence of the heating area on residual ice formation and on the drag and lift of the airfoil are elucidated.

The numerical scheme is described and validated in Section 2. Section 3 presents the simulation of icing with a heater for NACA 0012.

2. Numerical Scheme and Validation

The numerical simulation is based on the Euler–Lagrange method. The computation comprises four steps: (1) grid generation; (2) computation of the flow field; (3) computation of droplet trajectories; and, (4) thermodynamics computations. Each step is computed using in-house code and it is described below.

2.1. Grid Generation

The computational targets were NACA 0012 airfoil for the icing simulation and NACA 0013 airfoil for the validation. Figure 1 shows the computational domain system and grid. In this study, C-type grids and an overset grid system were employed. Interpolation between the main grid and sub-grid was performed using the Lagrangian interpolation method. These grids were generated based on Hermite polynomials [30]. A main grid of $221\times 71$ points was used to simulate the entire flow field around the airfoil, and the sub-grid of $301\times 51$ points had sufficiently high resolution to correctly obtain the ice shape around the leading edge and boundary layer on the airfoil. The total number of grid points was approximately $26,000$. The convergence of the ice shape by the grid points had preliminarily been confirmed through our previous investigation [31].

2.2. Flow Field Computation

The flow around the airfoil is assumed to be two-dimensional, compressible, and fully turbulent. The governing equations are the continuity, Navier–Stokes, energy, and transport equations of the turbulent kinetic energy k and its dissipation rate $\epsilon$. The Kato–Launder k-$\epsilon$ turbulent model [32] with a wall function was employed to suppress the over-production of turbulent eddy viscosity around the leading-edge region. The governing equations are discretized using a second-order upwind TVD (total variation diminishing) scheme [33] for the inviscid terms, a second-order central difference scheme for the viscous terms, and an LU-ADI (lower upper-alternating direction implicit) scheme [34] for time integration. The first grid points from the airfoil surface are in the range of 50 to 300 in wall units. The $L_2$ norm is employed as a convergence criterion.

2.3. Droplet Trajectory

In the droplet trajectory computation, we assume that the size and concentration of droplets are sufficiently small. Additionally, the collision and splitting of the droplets are assumed to be negligible. Thus, the trajectory was calculated based on the Lagrangian approach using the one-way coupling method. In other words, the droplet motion is affected by the flow field, whereas the droplet does not affect the flow field. The motion equation of a droplet is a simplified Basset–Boussinesq–Oseen (BBO) equation and it is expressed, as follows:

$$\begin{array}{c}\hfill \frac{d\overrightarrow{U_p}}{dt}=\frac{3}{4}{C}_D\frac{{\rho }_g}{{\rho }_d}\frac{1}{d_d}\overrightarrow{U_r}\left|\overrightarrow{U_r}\right|.\end{array}$$

(1)

where $\overrightarrow{U_p}$ is the droplet velocity and $\overrightarrow{U_r}$ is the relative velocity between the droplet and the surrounding flow; ${\rho }_g$ and ${\rho }_d$ are the gas (air) and droplet density, respectively; $d_d$ is the droplet diameter, and $C_D$ is the drag coefficient. Because we assume that the droplet does not deform or rotate, the Schiller model [35] is used for the drag coefficient of a sphere; it is defined as follows:

$$\begin{array}{c}\hfill C_D=\frac{24}{Re}\left(1+0.15R{e}_d^{0.687}\right),\end{array}$$

(2)

where $R{e}_d$ is the droplet Reynolds number based on the droplet diameter $d_d$, relative velocity $U_r$, fluid density ${\rho }_f$, and fluid viscosity ${\mu }_f$.

2.4. Thermodynamics

In the thermodynamic computation, which is, the calculation of ice growth, the weak coupling method is used because the time scales of the flow field and icing are significantly different. The EMM [29] is used as the icing model and is expressed, as follows:

$$\begin{array}{ccc}\hfill \frac{\partial T_i}{\partial t}& =& \frac{k_i}{{\rho }_i{C}_{pi}}\frac{{\partial }^2{T}_i}{\partial y_i^2},\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \frac{\partial T_w}{\partial t}& =& \frac{k_w}{{\rho }_w{C}_{pw}}\frac{{\partial }^2{T}_w}{\partial y_w^2},\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill {\rho }_i\frac{\partial B_i}{\partial t}+{\rho }_w\frac{\partial B_w}{\partial t}& =& {\dot{m}}_{im}+{\dot{m}}_{in}-{\dot{m}}_{es},\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\rho }_i{L}_f\frac{\partial B_i}{\partial t}& =& k_i\frac{\partial T_i}{\partial y_i}-k_w\frac{\partial T_w}{\partial y_w}.\hfill \end{array}$$

(6)

Here, T, k, $L_f$, B, t, y, and $\rho$ are the temperature, thermal conductivity, latent heat of solidification, thickness, time, wall normal, and density, respectively; subscripts i and w denote the ice and water, respectively. Additionally, ${\dot{m}}_{im}$, ${\dot{m}}_{in}$, and ${\dot{m}}_{es}$ denote the mass flow rates of impingement, runback-in, and evaporation or sublimation, respectively. The EMM is based on the mass and energy conservation law of the ice and water and the equation of phase change at the interface between the ice and water [31]. In the EMM, the temperature profile is not calculated explicitly and the heat transfer is computed based on the assumption of a linear temperature profile for ice and water. The present MEMM is in accordance with this assumption.

Figure 2 shows the EMM is widely used in icing simulations and the framework of the thermodynamics model of the MEMM. Here, there are four heat-receiving terms (aerodynamic heating $Q_a$, kinetic energy of incoming droplets $Q_k$, heat brought in by runback $Q_{in}$, and latent heat release (only for the rime ice condition) $Q_l$) and four heat-releasing terms (convection $Q_c$, cooling by the incoming droplets $Q_d$, evaporation and sublimation $Q_{es}$, and radiation $Q_r$). In the EMM model, the ambient temperature $T_{air}$ and freezing temperature $T_f$ are used as the wall surface and water film temperatures, respectively.

In contrast, in the modified EMM model developed in the present study, the contact temperature $T_c$ is used as the wall surface and water film temperatures $T_w$ in order to consider in the heating effect from the wall. The contact temperature $T_c$ at the interface between wall and the water film is evaluated using the impinging droplet and runback water temperatures, as

$$\begin{array}{ccc}\hfill T_w& =& \frac{{\dot{m}}_{im}{T}_{air}+{\dot{m}}_{in}{T}_{in}}{{\dot{m}}_{in}+{\dot{m}}_{in}},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill T_c& =& \frac{\sqrt{{\rho }_w{c}_w{k}_w}{T}_w+\sqrt{{\rho }_{wall}{c}_{wall}{k}_{wall}}{T}_{heater}}{\sqrt{{\rho }_w{c}_w{k}_w}+\sqrt{{\rho }_{wall}{c}_{wall}{k}_{wall}}},\hfill \end{array}$$

(8)

respectively, where $\dot{m}$, $\rho$, c, k, and T are the mass flow rate, density, specific heat, heat conductivity, and temperature, respectively, and the subscripts $im$, $in$, $air$, w, $heater$, and $wall$ denote the impingement, runback in, air, water, heater, and wall, respectively. Accordingly, the wall surface temperature is $T_c$. If $T_c$ is lower than the freezing point, icing (phase change) occurs, and the EMM is used to calculate ice layer growth; if $T_c$ is higher than the freezing point, icing does not occur. The runback water temperature was considered as $T_c$ in order to model heat conduction from the heater, where it was estimated by substituting $T_{heater}$ into $T_{air}$ in the unheated region. On the other hand, the runback water temperature is assumed to be equal to the freezing point in the EMM.

For the ice thickness computation, essentially the EMM is employed, although the energy brought in by runback water $Q_{in}$ is estimated using $T_c$ of the upstream neighboring computational cell, as follows;

$$\begin{array}{c}\hfill Q_{in}={\dot{m}}_{in}{C}_w(T_c-T_f).\end{array}$$

(9)

We then reconstruct the grids and compute the flow without ice-shape smoothing because the smoothing artificially affects the icing location. Note that the one-shot computation estimated ice shape.

2.5. Validation of Flow and Heating Simulation

The numerical results are compared with experimental results to validate the simulation, as described below.

First, the aerodynamic performance of clean airfoil (without icing and heating) is compared with the experimental results obtained by Sheldahl and Klimas [36] and Harris [37]. The present study simulated the aerodynamic performance of a NACA 0012 airfoil at Reynolds number $Re=2,800,000$ and Mach number $Ma=0.17$, where the $Re$ is defined by the chord length and the uniform flow velocity and $Ma$ is defined by the sonic and ambient air speeds. In addition, we simulated the flow with a $1.5$ times finer grid system. Figure 3 shows the lift and drag coefficients ($C_L$ and $C_D$, respectively), which are normalized by the dynamic pressure and chord length, with a different angle of attack (AoA) in comparison to those for $Re=2,000,000$, $Re=3,000,000$, and $Re=5,000,000$. Sheldahl and Klimas [36] used a clean airfoil and Harris [37] used a clean airfoil with tripping wire. The present results for lift coefficients show reasonable agreement with experimental results from Harris [37] without grid dependency. For the drag coefficient, the Reynolds number dependency was weak for AoA $\le 4°$. Because the numerical results obtained in this study are in reasonable agreement with the experimental values at an AoA $\le 4°$, we focus on a case with a small AoA ($0°\le$ AoA $\le 4°$), as described below. Additionally, Figure 4 shows pressure coefficients $C_p$ defined as

$$\begin{array}{c}\hfill C_p=\frac{p-p_0}{0.5\rho U^2},\end{array}$$

(10)

where p is the static pressure on the airfoil, while $p_0$ and U are the static pressure and uniform flow velocity, respectively, at the inlet. This clearly indicates the validity of the present numerical simulation for the flow field. In addition, the results obtained for the $1.5$ times finer grid system in all directions were plotted. $C_p$ and $C_L$ show no grid dependency. On the other hand, $C_D$ shows grid dependency and it seems to converge toward the experimental results. The $y^{+}$ values for the first grid locations from the wall were distributed from 50 to 300 in the present simulation and from 30 to 70 for the finer-grid case. Therefore, the wall function model should be changed to converge more closely toward the experimental result.

Second, the heater temperatures for the heating of NACA 0013 are compared. The experimental data were obtained from Asaumi et al. [25], who conducted an anti-icing experiment using a heater under constant heat flux. We conducted an icing simulation to reproduce their experiment for validation; the numerical conditions are presented in

Table 1. Numerical conditions for heater validation [25].

Airfoil	NACA 0013
Chord length	[m]	0.12
Angle of attack	[°]	0
Freestream velocity	[m/s]	30
Freestream temperature	[°C]	−8
Reynolds number	[-]	285,000
Median volume diameter (MVD)	[$\mathsf{\mu }$m]	58
Liquid water content (LWC)	[g/m3]	2.8
Initial droplet temperature	[°C]	−8
Total droplet number	[-]	1,000,000
Exposure time	[s]	300
Ambient pressure	[kPa]	101.325
Heating condition	Constant heat flux
Heating region	0–48% chord
Wall material	Titanium
Wall thickness	[$\mathsf{\mu }$m]	6

. The initial velocity and temperature of the droplet were equal to those under the main stream condition; the heater was set from the leading edge to $48\%$ of the chord length c (i.e., $48\%c$). To reproduce the constant heat flux condition, we adjusted the heater temperature $T_{heater}$ so that the temperature gradient (represented by $T_c$ and $T_{heater}$) in the wall-normal direction inside the wall is kept constant. The total number of droplets is 1 million, where the droplets are randomly arrayed over 1 chord width at the 7-chord-upstream location from the leading edge; thus, the formed ice thickness was modified by the LWC and MVD at the inlet. Although the arrangements of the heater and thermometer were different from those in the experiment, we can still compare the numerical results with the experimental results because the heater was very thin in the experiment. As shown in Figure 5, the heater temperature $T_c$ is in reasonable agreement with the experimental results. Hence, these results support the validity of the proposed numerical scheme.

3. Anti-Icing Performance of Heater for NACA 0012

This section describes the anti-icing simulation with the NACA 0012 airfoil. The numerical conditions are presented in

Table 2. Numerical conditions for NACA 0012 [38].

Airfoil	NACA 0012
Chord length	[m]	0.53
Angle of attack	[°]	4
Freestream velocity	[m/s]	58.1
Freestream temperature	[°C]	−27.8
Reynolds number	[-]	2,800,000
Median volume diameter (MVD)	[$\mathsf{\mu }$m]	18.0
Liquid water content (LWC)	[g/m3]	1.3
Initial droplet temperature	[°C]	−27.8
Total droplet number	[-]	1,000,000
Exposure time	[s]	480
Ambient pressure	[kPa]	95.61
Heating condition	Constant temperature
Heating wall temperature $T_{heater}$	[°C]	10.0
Heating region	1–12% chord (Case 1–Case 12)
Wall material	Aluminum
Wall thickness	[mm]	10

, which refers to an experimental study by Olsen et al. [38]. The free stream and droplet temperatures were set to −27.8 °C, which is the rime ice condition. Here again, rime icing is that the impinged super-cooled droplets freeze instantaneously, whereas the glaze icing is that a water film runs back along the wing and freezes. The total number of droplets is 1 million, where the droplets are randomly arrayed over 1 chord width in the 7-chord-upstream location from the leading edge, as shown in Figure 6. As an anti-icing measure, we installed a heater around the leading edge, as shown in Figure 6. The heating region was varied, which resulted in twelve different cases. Each case was named in accordance with the heating region. For example, the heating region of $x=0$–$1\%c$ is Case 1, while that of $x=0$–$2\%c$ is Case 2.

Figure 7 shows $C_L$ and $C_D$ normalized by the same values in the case without icing as a function of heating area. The normalized values are defined, as follows:

$$\begin{array}{c}\hfill C_L^{*}=\frac{C_L-C_{Lclean}}{C_{Lclean}}\times 100[\%],\end{array}$$

(11)

$$\begin{array}{c}\hfill C_D^{*}=\frac{C_D-C_{Dclean}}{C_{Dclean}}\times 100[\%],\end{array}$$

(12)

where the subscript clean shows the value without icing. These coefficients exhibited different behavior for the heating region. As the heating area increased, $C_L$ decreased, exhibited a negative peak at $4\%c$ of the heating area (Case 4), and then recovered. Moreover, $C_D$ gradually decreased and converged to 0%. Below, four cases, namely the clean airfoil without icing, Case 1, Case 4, and Case 10, are considered to discuss the flow fields.

Figure 8 shows the ice thicknesses. On the suction surface, the ice area decreases as the heating area increases. On the pressure side, the round ice forms very close to the leading edge in Cases 1 and 2. In contrast, in Case 4, the ice area spreads downstream, owing to the runback. In Cases 6 to 10, the ice area decreases, owing to the increased heating area.

Figure 9 shows the streamline and ice shape. Figure 9a shows the case without icing, where the flow smoothly moves downstream along the airfoil shape. According to the previous icing simulation by Hayashi and Yamamoto [31], icing mainly occurs close to the leading edge and is termed “round ice”. Figure 9b shows the heating in Case 1 ($1\%c$ heating). Icing did not occur around the leading edge where the wall was heated; therefore, the stagnation region became larger than that of clean airfoil. Round ice appeared close to the leading edge on the pressure side. On the suction side, a thin ice layer formed owing to the runback phenomenon that is a liquid film flow into downstream without icing. In addition, the ice layer behind the heater was similar to that in the experimental results [39,40]. Subsequently, the ice was melted by the heater (runback water), and the water flowed downstream and froze again owing to the heat loss to the airfoil wall and surrounding air. Interestingly, a thin ice layer formed downstream from the heater because the heater temperature was 10 °C and the runback water did not instantly freeze as soon as it passed the end of the heater. In the $4\%c$ heating (Case 4, Figure 9c), compared with Case 1, an ice layer on the suction side appeared downstream, and large flow separation was not observed. On the pressure side, the ice was thinner when compared with Case 1 and covered the surface. Flow separation was also not observed. In Case 10 (Figure 9d), as the heating area increased, the ice layers on both sides were thinner and formed further downstream.

Figure 10 shows the static pressure and ice shape. For the clean airfoil, as shown in Figure 10a, positive and high pressure appeared around the stagnation point, while small and negative pressure appeared on the suction side. In Case 1, the positive pressure increased at the stagnation point, owing to the blockage effect by the round ice on the pressure side. On the suction side, owing to the large flow separation observed in Figure 9b, a low-pressure region was created from the edge of the ice layer. In Case 4, the positive pressure region expanded, because the ice layer formed downstream on the pressure side, compared with Case 1, whereas the negative pressure on the suction side decreased, and its area became small. In Case 10, the pressure distribution is identical to that of clean airfoil. Although pressure variation was observed close to the ice edge, this variation was small.

Figure 11 shows the pressure coefficients on the airfoil surface for Cases 1, 4, and 10, along with those of the clean airfoil. The upper and lower blanches are the pressure and suction sides, respectively. For the clean airfoil, the pressure coefficient peaked at $x/c\approx 0$: on the pressure side, then gradually decreased and converged to a slightly negative value. On the suction side, the pressure coefficient suddenly decreased and gradually increased. In Case 1, the pressure coefficient dropped lower than that of clean airfoil, whose regions correspond with the flow separations caused by the icing. The decrease in pressure caused an increase in both lift and drag, as shown in Figure 7, and the decrease in $C_p$ near the ice ridge was also reported in an experimental study [41]. In Case 4, the pressure coefficient peaked at the stagnation point close to the leading edge, but another peak (or discontinuity profile) appeared downstream of the leading edge, at the edge of the ice. Moreover, the suction peak on the suction surface was attenuated by the formed ice, and this resulted in the decrease in lift. For $x/c\gtrsim 0.2$, the pressure coefficient converged to that of clean airfoil. In Case 10, the profile agrees with that of clean airfoil, except for the second peak at $x/c\approx 0.15$ because the heating area is wide and the icing layer is thin. Accordingly, the lift and drag contribution are small. In general, the lift decreases owing to icing, especially in experiments (e.g., [39,40,42]), and the present results are inconsistent with this tendency. In actual scenarios, the ice surface works as a rough surface because the shape of ice is three-dimensional and complex, and a hump-like ice formation appears behind the heater [40]. Accordingly, the icing shape deteriorates the aerodynamic performance of airfoil and the flow separation differs significantly from that in the present two-dimensional simulation. The present results are obtained using the two-dimensional simulation without considering the surface roughness and the icing helps in improving the aerodynamic performance. The lift is increased because of the large flow separation on the suction side for the small heating area cases and the thickened wing thickness for the large heating area cases. Moreover, the hump-like ice formed behind the heater [40], owing to large amounts of runback water, is smaller than that formed in the real situation. Therefore, it implies that there is scope to improve the icing model for large amounts of runback water created by the heater.

4. Conclusions

In this study, we simulated the occurrence of icing with a heater on a NACA 0012 airfoil. A modified extended Messinger model was developed in order to consider surface heating and runback water temperature. A combination of Euler–Lagrange approaches was adopted in the icing simulation. The flow field around the airfoil was calculated using the Eulerian approach, while the droplet trajectory was obtained using the Lagrange approach. The heater behavior was validated by considering experimental results, and reasonable agreement was obtained. For NACA 0012, an icing simulation with a heater was conducted under the rime ice condition. The heater was installed around the leading edge, and the effects of the heating area on icing and aerodynamic performance were investigated.

The present simulation method using the modified extended Messinger model is more suitable for use in de-icing simulations for both rime and glaze ice conditions, because it reproduces the thin ice layer behind the heater due to the runback phenomenon. It is well known that the extended Messinger model is suitable for both rime and glaze icing conditions as compared with the original Messinger model. Therefore, our proposed modification of the extended Messinger model is important in predicting the effect of a heater on icing.

The results obtained in this study reveal that the lift and drag depend on the heater area. In the case of the heating area of 4% of the chord length, the lift coefficient deteriorated, becoming smaller than that of the clean airfoil. However, increasing the heating area improved the lift coefficient by up to 4%. The large ice vanished owing to the heater, although a thin ice layer formed behind the heater. On the other hand, the drag decreased and converged to that in the case of the clean airfoil. The present study shows that the de-icing effect created by the heater can be simulated using the modified extended Messinger model. Although the computational target was the NACA 0012 airfoil in this study, the results suggest the possibility of using the heater for de-icing in other icing conditions and situations. Therefore, a more comprehensive study using a numerical method that is more robust against, for example, larger AoA cases and other airfoil shapes, will be pursed in future work.