Analysis of Electro-Thermal De-Icing on a NACA0012 Airfoil Under Harsh SLD Conditions and Different Angles of Attack

1. Introduction

Aircraft icing is a significant issue in aviation where an aircraft flies through a cloud of supercooled water droplets that freeze upon impact with its surfaces [1,2,3,4]. The most vulnerable areas prone to icing are wing leading edges, engine inlets, and windshield windows since these parts are directly exposed to incoming airflow and droplet impacts. Ice accumulation changes aerodynamic profiles, reducing performance, adding drag, and affecting stability and control [5,6,7]. Moreover, crushed ice fragments may damage internal components or enter the engine, reducing thrust or blocking the inlets [8,9]. Understanding when and how icing forms determines the variables that control its severity.

Several environmental and operational variables influence the extent of ice accretion. Airspeed, temperature, droplet size, liquid water content (LWC), and exposure time are some of the variables that affect the quantity of ice accretion under icing circumstances [10].The primary ice types are rime, glaze, and mixed: rime forms when small supercooled droplets freeze on impact at low temperature and low LWC, producing a rough, opaque, porous deposit; glaze forms when larger droplets and higher LWC near 0° allow water to spread and refreeze, yielding dense, transparent layers with runback horns; and mixed ice exhibits characteristics of both and can severely degrade aerodynamic performance [10,11,12]. Within this spectrum, one regime–Supercooled Large Droplet (SLD) icing–demands particular attention due to its operational risk.

SLD icing gained heightened attention in the aerospace industry after the fatal 1994 ATR-72 crash, which was attributed to ice buildup on the aircraft’s trailing-edge flaps leading to loss of control. This accident underscored the importance of understanding in-flight SLD icing and its effects on the aerodynamic performance of high-lift configurations. In response, research efforts have intensified to better understand SLD icing phenomena and associated hazards [13,14]. These risks, in turn, motivate robust and energy-efficient ice protection strategies on modern aircraft.

Commercial aircraft employ various Ice Protection System (IPS) technologies to prevent or remove ice. The most commonly used systems include engine bleed-air heating, electro-thermal heaters, pneumatic de-icing boots, and weeping wing (fluid) systems. Among these, hot-air and electro-thermal ice protection systems (ETIPS) are especially prevalent, with ETIPS gaining popularity due to their higher energy efficiency [15,16]. ETIPS are favored in modern aircraft designs because of their relative simplicity and energy efficiency [16], offering power savings of up to 33% compared to traditional thermo-pneumatic systems [17]. However, selecting how and where to apply heat introduces fundamental trade-offs between continuous prevention and periodic removal. IPS design strategies typically combine anti-icing and de-icing methods depending on aircraft requirements: some critical surfaces are continuously heated (anti-icing) while others allow ice to accumulate briefly before periodic removal (de-icing) [18]. A key trade-off exists between these approaches: anti-icing systems prevent ice from forming but consume more energy, whereas de-icing systems allow a thin ice layer to form and then shed it, using less energy overall [19,20]. Consequently, ETIPS are often designed to maximize efficiency by operating in a de-icing mode, which has been shown to be more energy-efficient than continuous anti-icing [21,22]. Furthermore, hybrid approaches that combine passive and active techniques (e.g., superhydrophobic coatings with electro-thermal heating) have proven effective in lowering energy consumption [23,24]. Because these design choices carry safety and energy implications, rigorous assessment methods are essential.

Developing and certifying an IPS can be costly, time-consuming, and even dangerous—especially if flight tests in natural icing conditions are required. Consequently, there is growing interest in Certification by Analysis (CbA), wherein numerical simulations supplement or replace some physical icing tests [25]. Accurate computational models are essential for designing efficient aircraft IPS. These simulations must capture complex phenomena including aerodynamic flow, droplet impingement, phase transitions (freezing and melting), heat transfer, and water film runback behavior [7,26,27,28]. This need has driven a mature ecosystem of icing simulation tools.

Several icing codes have been developed to estimate ice accretion on complex geometries. Representative examples include ANSYS FENSAP-ICE (commercial), NASA’s LEWICE (2D/3D; government-distributed), ONERA’s 2D/3D icing modules, and research/industry tools such as ICECREMO, CANICE, and PoliICE [7,29,30]. Other institution-specific frameworks—including the NRC Morphogenetic Code, IGLOO3D, PoliMIce, SIMBA/MESS3D, and ONICE3D—are not generally available to the public [31]. Availability and licensing therefore vary (commercial, government-distributed, or restricted), and most packages are not openly released or open-source. These tools typically employ Reynolds–Averaged Navier–Stokes (RANS) solvers, with the one-equation Spalart–Allmaras model commonly used for boundary-layer flows where icing originates [32]. Furthermore, time-resolved ice-shape evolution can be captured via multi-shot approaches combined with automatic remeshing, improving predictions of droplet behavior and aerodynamic performance [31].

Conjugate heat transfer (CHT)—the simultaneous solution of solid heat conduction in the multilayer leading edge (with embedded heaters) and fluid-side heat and mass transfer (convection, radiation, and freezing/melting at the interface)—is essential for realistic simulation of electro-thermal anti-icing and de-icing [7,33]. Recent work has shifted from steady approximations to unsteady, phase-change–aware formulations that resolve transient water-film/runback dynamics [7]. The Shallow-Water Icing Model (SWIM) provides an efficient surface-film description that captures runback and wall-temperature feedback during icing events [4]. Despite these developments, rigorous evaluations of electro-thermal de-icing performance across angles of attack and under severe SLD conditions remain limited. Understanding how droplet size and angle of attack influence de-icing effectiveness is critical for improving system reliability in extreme icing.

To address this gap, we investigate a cyclic electro-thermal ice protection system (ETIPS) originally designed for standard flight conditions, now evaluated under more realistic SLD environments. Using a validated RANS-based icing simulation framework that couples Eulerian droplet transport, shallow-water film/runback, phase change, and conjugate heat transfer, we simulate ice formation and de-icing for supercooled droplets larger than 40 $\mu$m (Langmuir D distribution [34]) across a range of angles of attack, while holding airfoil geometry and heater power fixed. The insights gained are intended to guide the design of more effective de-icing systems, improving both flight safety and operational efficiency.

2. Methodology

Icing and de-icing simulations are crucial for mitigating hazards from ice accretion on aerodynamic surfaces. This section presents a methodology for a coupled multiphase and conjugate-heat-transfer problem involving air (gas), supercooled liquid water droplets, and ice, thermally coupled to the solid surface. The influencing factors include atmospheric conditions (temperature, pressure, droplet size, liquid water content, surface roughness, and velocity) as well as the surface temperature. As shown in Figure 1, the simulation of the icing and de-icing process involves four main steps. In the first step, the airflow solution is computed, which simulates steady-state, compressible, viscous flow using the Finite Element Method (FEM) and the one-equation Spalart–Allmaras (S-A) turbulence model [31,35,36]. The second step is the solution of the droplet impingement problem, where an Eulerian formulation is adopted to treat droplet transport as a continuous-phase flow. The third step computes ice accretion, water runback, and ice surface deformation. The fourth step simultaneously solves heat conduction in the solid domain and the energy balance at the fluid–solid interface. At this stage, conduction from embedded heating pads through the material layers governs the heat flux at the solid surface. This workflow, detailed below, is implemented under the ANSYS FENSAP-ICE framework [35].

2.1. Governing Equations for the Icing/De-icing Model

For the dry-air flow, the conservation equations for mass, momentum, and energy are solved under the assumption of a one-way interaction, where droplets carried by the air do not affect the air velocity field. An Eulerian two-fluid approach is employed to calculate droplet impingement on 3D objects using results from the airflow solution. This approach solves for droplet velocities and water volume fractions, treating particle flow as a continuum. Based on inputs such as Liquid Water Content (LWC) and Median Volumetric Diameter (MVD), it generates a droplet distribution over the airfoil. Within an Eulerian two-fluid formulation, the dispersed water phase is governed by the volume-fraction transport (1) and momentum (2) equations, solved alongside the carrier-gas Euler/Navier–Stokes equations to model droplet transport:

$${\displaystyle \frac{\partial \alpha }{\partial t}}+\overrightarrow{\nabla }·\left(\alpha {\overrightarrow{V}}_d\right)=0$$

(1)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \left(\alpha {\overrightarrow{V}}_d\right)}{\partial t}}& +\overrightarrow{\nabla }·[\alpha {\overrightarrow{V}}_d\otimes {\overrightarrow{V}}_d]={\displaystyle \frac{C_{D}R{e}_d}{24K}}\alpha ({\overrightarrow{V}}_a-{\overrightarrow{V}}_d)\hfill \\ \hfill & +\alpha \left(1-{\displaystyle \frac{{\rho }_a}{{\rho }_d}}\right){\displaystyle \frac{1}{F{r}^2}}\overrightarrow{g}\hfill \end{array}$$

(2)

where ${\rho }_d$ is the density of liquid water, $\alpha$ is the water volume fraction and ${\overrightarrow{V}}_d$ the droplet velocity; for a given liquid water content (LWC), the water volume fraction $\alpha$ can be obtained as $\alpha =\mathrm{LWC}/{\rho }_d$. The first term on the right side of the equation represents the aerodynamic drag on droplets of mean diameter d, which is proportional to the relative velocity of the droplet in the air, its drag coefficient $C_D$, its Reynolds number $R{e}_d$, and an inertial parameter K, which are defined in equation (3) and equation (4) [35].

$$R{e}_d={\displaystyle \frac{{\rho }_{a}d{V}_{a,\infty }\parallel \overrightarrow{V_a}-\overrightarrow{V_d}\parallel }{{\mu }_a}}$$

(3)

$$K={\displaystyle \frac{{\rho }_d{d}^2{V}_{a,\infty }}{18L_{\infty }{\mu }_a}}$$

(4)

The second term represents the gravity and buoyancy forces, which depend on the local Froude number which is defined in equation (5):

$$Fr={\displaystyle \frac{\parallel V_{a,\infty }\parallel }{\sqrt{L_{\infty }g}}}$$

(5)

For simulating the splashing, bouncing, and breakup of water droplets, just for droplet distributions with a mean diameter greater than 40 $\mu m$ (SLD) [35], we consider appropriate physical models. For modeling the droplet breakup, we utilize the Pilch and Erdman correlation, which is based on a critical Weber number that must be exceeded for breakup to occur:

$$W{e}_{critical}=12\left(1+1.077O{h}^{1.6}\right)$$

(6)

Where the Ohnesorge number (Oh) is a dimensionless parameter that relates viscous forces to inertial and surface tension forces, defined as:

$$Oh={\displaystyle \frac{{\mu }_d}{\sqrt{d{\rho }_d{\sigma }_d}}}$$

(7)

Where ${\mu }_d$ is the droplet’s viscosity, d the droplet’s diameter, and ${\sigma }_d$ the droplet’s surface tension, and for Weber numbers below 13, the breakup mechanism is neglected. To predict whether a droplet will bounce or splash upon impact, we apply the Mundo model. This model modifies the droplet collection efficiency on a surface by determining whether all the impinging water bounces or if a portion of it splashes. In both cases, no splashed or bounced droplets are reintroduced into the computational domain. The condition of splashing or bouncing is determined using the Mundo parameter $k_m$, defined as:

$$k_m=k_c^{5/8}$$

(8)

Where $k_c$ is the Cossali parameter, which depends on the Weber and Ohnesorge numbers. Splashing will occur when:

$$k_c=O{h}^{-2/5}We\ge 540R^{-0.35}$$

(9)

Where R is the surface roughness height and is set to 500 $\mu m$ for all droplet sizes. The secondary-to-primary droplet diameter ratio ($d_s/d_i$) and the number of secondary particles ($n_s$) are calculated as:

$${\displaystyle \frac{d_s}{d_i}}=8.72e^{-0.028k_m}$$

(10)

$$n_s=1.676\times {10}^{-5}{\left(k_m\right)}^{2.539}$$

(11)

A system of two partial differential equations based on thin-film approximation is solved on solid surfaces. The first equation expresses mass conservation (equation (12)), and the second expresses conservation of energy (equation (14)).

$${\rho }_f\left[{\displaystyle \frac{\partial h_f}{\partial t}}+\overrightarrow{\nabla }·\left({\overline{V}}_f{h}_f\right)\right]=V_{\infty }\left(LWC\right)\beta -{\dot{m}}_{\mathrm{evap}}-{\dot{m}}_{\mathrm{ice}}$$

(12)

Where the subscript f identifies water film properties, meaning that $h_f$ is the film thickness and ${\overline{V}}_f$ is the velocity of the water film averaged across the thickness and assuming no-slip condition at the wall. The three terms on the right-hand side represent mass transfer by droplet impingement (source), liquid mass loss by evaporation (sink), and ice accretion (sink). Parameter $\beta$ (droplet collection efficiency) is defined in equation (13) [35]:

$$\beta =-{\displaystyle \frac{\alpha \overrightarrow{V_d}·\overrightarrow{n}}{\left(LWC\right)V_{\infty }}}$$

(13)

$$\begin{array}{c}{\rho }_f\left[{\displaystyle \frac{\partial h_f{c}_f{T}_f}{\partial t}}+\overrightarrow{\nabla }·\left({\overrightarrow{V}}_f{h}_f{c}_f{T}_f\right)\right]=\beta \left(LWC\right)V_{\infty }\left[e_d+{\displaystyle \frac{{\overrightarrow{V}}_d^2}{2}}\right]\hfill \\ -{\dot{m}}_{\mathrm{evap}}\left(L_{\mathrm{evap}}\right)+{\dot{m}}_{\mathrm{ice}}\left(L_{\mathrm{fus}}-c_{ice}{T}_{\mathrm{ice}}\right)\\ \hfill +\sigma \epsilon \left(T_{\infty }^4-T_f^4\right)-c_h\left(T_f-T_{\mathrm{rec}}\right)+Q_{\mathrm{anti}-\mathrm{icing}}\end{array}$$

(14)

Where $\overrightarrow{n}$ is the unit normal vector to the surface and the last three terms of the equation (14) are the radiative, convective, and conductive heat fluxes that could be supplied by a de-icing or anti-icing heater device. The parameters $\rho ,c_f,L_{\mathrm{evap}},c_{\mathrm{ice}},L_{\mathrm{fus}},\sigma ,\epsilon$ are physical properties of water and ice. The local wall shear stress $\tau$ and the convective heat flux $Q_h$ should be supplied by the airflow solver. The droplet impingement model provides local values of the collection efficiency $\beta$, droplet sensible energy $e_d$, and droplet impact velocity ${\overrightarrow{V}}_d$. $Q_{\mathrm{anti}-\mathrm{icing}}$ is the de-icing or anti-icing heat flux obtained from the solid heat conduction solver. The solid heat conduction solver handles domains composed of multiple materials with different thermal properties. Heat conduction is inherently unsteady, as heater pads are switched on and off sequentially and heat is transmitted through the various solid layers. To account for both sensible and latent heat, the governing equation can be written with an enthalpy-based term to represent phase change (equation 15) [31,35]:

$${\rho }_s{C}_p{\displaystyle \frac{\partial T}{\partial t}}=\nabla ·\left(K\left(T\right)\nabla T\right)+{\rho }_s{\displaystyle \frac{\partial H_{pc}}{\partial t}},$$

(15)

where ${\rho }_s$ is the solid density, $C_p$ the specific heat, $K\left(T\right)$ the thermal conductivity, and $H_{pc}$ the enthalpy term representing latent heat release during freezing and absorption during melting.

2.2. Numerical model setup and validation

To validate the unsteady de-icing solver, electro-thermal de-icing experiments conducted in the NASA Lewis Icing Research Tunnel (IRT) were selected [37]. The NACA 0012 airfoil, with a chord length of 0.9144 m, was used for the experiments. The leading edge of the airfoil consists of a multi-layered composite structure, which includes an erosion shield, two elastomer layers, a fiberglass/epoxy layer, and a silicone foam layer. Seven heater pads are installed on the leading edge of the airfoil, although they are not symmetrically placed due to manufacturing constraints. These heaters are offset by approximately 4.7625 mm towards the upper surface of the airfoil, as shown in Figure 2.

This study employs a Quasi-3D simulation with a computational grid that contains a single cell in the spanwise (+z) direction under periodic boundary conditions in order to effectively simulate a larger domain. High resolution of the airfoil boundary layer is required to model ice growth and accurately predict shear stress and convective heat transfer. Both the fluid and solid domain computational grids are well-structured and meet the Spalart-Allmaras (S-A) turbulence model demands, with $y^{+}\le 1$ [10]. The S-A model is commonly used in airfoil icing simulations since it has been found to perform well in predicting ice accretion and aerodynamics of iced airfoils [38]. It has been shown to agree well with experimental data for lift, stall angles, and pressure distributions for various ice shapes [39]. The benefit of the S-A model over the k-omega SST model is that, as a one-equation model, it is computationally more efficient. This is the reason for the choice of the S-A model in this study.

As shown in Table 1, a first cell height of $2\times {10}^{-6}$ m and a growth rate of 1.14 were applied uniformly across all computational grids for the fluid domain. To ensure mesh independence, four distinct grid resolutions—coarse, medium, fine, and very fine—were employed. The boundary conditions for the mesh independence study were as follows: freestream temperature $T_{\infty }=266.483\mathrm{K}$, freestream velocity $U_{\infty }=44.704\mathrm{m}/\mathrm{s}$, surface roughness of $0.0005\mathrm{m}$, and an angle of attack of $0^{\circ }$. In this study, only the airflow around the airfoil was simulated; the effects of droplet impact and ice accretion were not included. After refining the mesh to the fine grid, the results showed negligible change, with only a 0.69% variation. Consequently, the fine grid was selected for all subsequent simulations in the study.

For the simulations involving airflow, droplet impact, and ice accretion on the external airfoil surface, the fine mesh was employed across all relevant solvers. However, for the solid domain and to simulate heat transfer within the aircraft wing and between the heater layers, a different computational grid was used (Figure 3). This grid had a first element height of $6.7\times {10}^{-5}$ m and consisted of a total of 32,000 cells.

The accuracy of capturing temperature changes on heaters during the electrothermal de-icing cycles is minimally affected by the time step, and to reduce computational cost in this study, the time step is set to $\Delta t=0.5$ seconds [40].

The heater pads, positioned between the elastomer layers, distribute heat to the airfoil. The fiberglass/epoxy and foam layers minimize heat leakage to the internal structure of the wing. The properties of the different layers of the electro-thermal de-icing system, as well as the positions of all seven heating elements, are presented in Table 2 and Table 3, respectively.

For the heating cycle, all heaters operate in a specific on-off pattern. Each cycle lasts 120 seconds, and the total heating period consists of 5 consecutive cycles, totaling 600 seconds. In the experiment, the start of the heating cycle is precisely synchronized with the initiation of water spray onto the airfoil, meaning the onset of ice accumulation aligns exactly with the beginning of the heating process. To ensure that no ice forms on the leading edge of the airfoil, Heater A remains on continuously throughout the 120-second heating cycle. The other heaters operate cyclically, turning on and off during the heating cycle. The power of each heater and the duration of its operation within each 120-second heating cycle are specified in Table 4.

The different icing cases examined in this study are presented in Table 6. An important point to note is that for cases 1 to 5, where the mean diameter of supercooled droplets is considered to be 20 $\mu m$, a monodisperse distribution has been used, meaning all droplets are assumed to have the same diameter. However, for cases 6 to 8, according to Appendix C and O of the Federal Aviation Regulations (FAR) part 25 [41,42], supercooled large droplets (SLD) are defined as having diameters greater than 40 $\mu m$. Appendix O specifically addresses severe icing conditions, including freezing drizzle and freezing rain, which involve larger supercooled droplets that can create hazardous ice accretion beyond protected surfaces. These conditions are classified as harsh flight environments. For the harsh boundary conditions, supercooled droplets with an MVD of 50 $\mu m$, 100 $\mu m$, and 200 $\mu m$ have been used. Unlike the normal conditions, where all droplets were monodisperse, in harsh conditions a Langmuir "D" distribution has been used. The droplet diameters ($D_i$) corresponding to each class (i) and the associated dimensionless ratios and liquid water content weights are summarized in Table 5 [43].

Table 5. Langmuir “D” distribution Ratios and weight fractions.

Class i	Ratio ${\mathit{r}}_{\mathit{i}}={\mathit{D}}_{\mathit{i}}$/MVD	Weight Fraction [$\mathit{\Delta }{\mathit{w}}_{\mathit{i}}/{\mathit{w}}_{\mathit{t}}$] (%)
1	0.31	5
2	0.52	10
3	0.71	20
4	1.00	30
5	1.37	20
6	1.74	10
7	2.22	5

Table 5. Langmuir “D” distribution Ratios and weight fractions.

Class i	Ratio ${\mathit{r}}_{\mathit{i}}={\mathit{D}}_{\mathit{i}}$/MVD	Weight Fraction [$\mathit{\Delta }{\mathit{w}}_{\mathit{i}}/{\mathit{w}}_{\mathit{t}}$] (%)
1	0.31	5
2	0.52	10
3	0.71	20
4	1.00	30
5	1.37	20
6	1.74	10
7	2.22	5

Table 6. Different cases of icing conditions tested in the present work.

Case	${\mathit{T}}_{\mathit{\infty }}$	${\mathit{U}}_{\mathit{\infty }}$	LWC	MVD	AOA	Roughness Height	Icing Time
	[K]	[m/s]	[g/m3]	[$\mu$m]	[deg]	[$\mu$m]	[min]
Normal Flight Condition
Case 1	266.483	44.704	0.78	20	0	500	10
Case 2	266.483	44.704	0.78	20	-4	500	10
Case 3	266.483	44.704	0.78	20	4	500	10
Case 4	266.483	44.704	0.78	20	8	500	10
Case 5	266.483	44.704	0.78	20	12	500	10
Harsh Flight Condition
Case 6	266.483	44.704	0.78	50	0	500	10
Case 7	266.483	44.704	0.78	100	0	500	10
Case 8	266.483	44.704	0.78	200	0	500	10

Table 6. Different cases of icing conditions tested in the present work.

Case	${\mathit{T}}_{\mathit{\infty }}$	${\mathit{U}}_{\mathit{\infty }}$	LWC	MVD	AOA	Roughness Height	Icing Time
	[K]	[m/s]	[g/m3]	[$\mu$m]	[deg]	[$\mu$m]	[min]
Normal Flight Condition
Case 1	266.483	44.704	0.78	20	0	500	10
Case 2	266.483	44.704	0.78	20	-4	500	10
Case 3	266.483	44.704	0.78	20	4	500	10
Case 4	266.483	44.704	0.78	20	8	500	10
Case 5	266.483	44.704	0.78	20	12	500	10
Harsh Flight Condition
Case 6	266.483	44.704	0.78	50	0	500	10
Case 7	266.483	44.704	0.78	100	0	500	10
Case 8	266.483	44.704	0.78	200	0	500	10

As indicated by Figure 2 and Table 3, heater B is within the impact zone of supercooled droplets and is directly subjected to droplet impingement and ice growth, and thus validation of it with experimental data is extremely critical. On the other hand, Heater D, which is outside the impact zone, may be affected only by runback water flow from the de-icing system, and its simulated results must agree reasonably well with experimental data.

As illustrated in Fig. Figure 4, the performance evaluation of heaters B and D highlights the accuracy and efficiency of the electro-thermal de-icing model. For heater B, the experimental data indicate that the temperature reaches 270 K within 100 s. The results from this study show strong alignment with the experimental trend as well as with other models such as LEWICE [37]. However, the present experiment predicts slightly higher temperatures, with 273 K recorded at 100 s and 292 K at 110 s. The discrepancy could be attributed to differences in the estimation of the water film, which significantly contributes to the thermal behavior in de-icing. Heater B’s temperature then decreases in the cooling regime to around 273 K, closely matching the experimental observation of 271 K. For heater D, which is located outside the impingement zone, the calculated temperature fluctuations also agree well with the experimental data. The present study produced similar results, remaining in close agreement with the measured values, whereas LEWICE overestimated the minimum temperature, leading to a significant discrepancy with the experimental observations.

3. Results and Discussion

3.1. Simulation of Ice Accretion Without De-Icing

In this work, a multi-shot ice accretion simulation technique is utilized to simulate ice accretion on an airfoil. The technique is used by subdividing the overall icing time (e.g., 10 minutes) into smaller intervals (e.g., five 2-minute stages). First, ice accretion is simulated on a clean airfoil for the first interval (2 minutes). Based on the ice profile created on the airfoil’s leading edge, a new geometry for the airfoil is created and the mesh is recreated. The simulation is reiterated for the subsequent 2-minute time interval, considering the modified airfoil geometry resulting from the ice accreted previously. This is an iterative procedure for all the intervals until the entire 10-minute icing time is simulated.

The main advantage of the multi-shot method over a single-shot, 10-minute simulation is that it is more precise. With every stage, the mesh is redefined to reflect the present airfoil shape accurately, particularly in areas of primary interest close to the leading edge. The continuous mesh refinement guarantees that the most critical aerodynamic and thermal phenomena associated with ice accretion are simulated with higher accuracy. Besides, the geometry of ice formed at each stage has a significant influence on the following ice accretion, especially on the forefront, where surface roughness and asymmetry play the most significant roles.

Figure 5 (a) presents the local water droplet collection efficiency on the NACA0012 airfoil in electro-thermal de-icing simulation. For normal icing condition (if angle of attack is $0^{\circ }$), the collection efficiency reaches a maximum value of 0.54 at the leading edge and decreases to zero when the distance from the leading edge is increased ($\mathrm{S}/\mathrm{C}\ge -0.0313\mathrm{and}\mathrm{S}/\mathrm{C}\le 0.0313$), since the angle of attack (AOA) is $0^{\circ }$, the curves are symmetrical relative to the upper and lower surfaces of the airfoil, delineating the limits of droplet impingement for a chord position of $\pm 0.0313$. As shown in Figure 5 (b), with a transition from normal to harsh SLD conditions, there is a noticeable increase in the collection efficiency at the leading edge, particularly due to the larger supercooled droplet diameter. This leads to a significant expansion of the impingement region. Under normal flight conditions, the impingement limit typically encompasses heaters A, B, and C. However, under harsh conditions, if the MVD is equal to 50 or 100 $\mu m$, heaters A, B, C, D, and E fall within the impingement zone, and as the MVD further increases to 200 $\mu m$, all the heaters are affected.

Figure 6 illustrates the ice accretion profiles on the NACA 0012 airfoil (chord length of 0.9144 m) under the flight and icing conditions specified in Table 6. As shown in Figure 6 (a), with the increase in the angle of attack from -4° to 12°, the amount of ice formed on the leading edge of the airfoil moves more towards the lower surface of the airfoil. As the angle of attack increases, the accumulated ice on the airfoil increases on the lower surface. This expansion of the impingement zone directly influences ice accretion behavior. As depicted in Figure 6 (b), with the increase in the average diameter of the supercooled droplets, the ice thickness on the leading edge remains relatively unchanged, but the ice progression along the chord of the airfoil becomes more pronounced on both the upper and lower surfaces. This phenomenon is closely tied to droplet impact zone and collection efficiency, as the larger droplet diameters enhance droplet break-up, splashing, and bouncing, thereby intensifying the spread of ice formation.

3.2. Aerodynamics Degradation

Ice buildup on the leading edge of an airfoil significantly affects its aerodynamic performance by disrupting the smooth airflow, increasing drag, and reducing lift. he altered leading-edge geometry promotes premature boundary-layer separation; the flow becomes turbulent and forms a separated vortex downstream of the ice, particularly on the suction side, as show in Figure 7 for case 1. This separation slows down the airflow, raises pressure, and causes a sharp drop in lift. The accompanying rise in drag further reduces efficiency and increases fuel consumption.

Studies show that the lift coefficient can decrease by up to $48\%$ [6], while drag can increase by as much as 36% to 200% [44]. These changes, influenced by factors such as temperature and the duration of ice accretion, not only degrade aerodynamic performance but also create significant safety risks, including the potential for loss of control.

As shown in Table 7, under normal flight conditions, an increase in the angle of attack causes a decrease in the lift coefficient and an increase in the drag coefficient. The worst condition is achieved at an angle of attack of ${12}^{\circ }$, with the lift coefficient decreasing by 11.1% and the drag coefficient increasing by 48.2%. In harsh flight conditions, with the angle of attack fixed at $0^{\circ }$, the lift coefficient remains essentially unchanged as the median volume diameter (MVD) of supercooled droplets increases, whereas the drag coefficient rises markedly by $36.5\%$ at an MVD of 200 $\mu$m.

3.3. Analysis of electro-thermal de-icing system in normal icing conditions

When supercooled droplets in the air collide with the wings and body of an aircraft under normal conditions, they form a layer of ice. However, with the use of electro-thermal de-icing systems, heaters embedded beneath the wing surface warm the surface, causing the accumulated ice to melt and form a thin film of water. This thin water layer is pushed backward over the wing by airflow, and if these new areas are not thermally protected by the de-icing system, the water refreezes and forms ice. The concept of cyclic electro-thermal de-icing involves multiple heaters with varying power densities and lengths installed at different locations inside the wing, directly beneath its surface. These heaters operate cyclically, turning on and off in sequence. The de-icing process thus includes a combination of ice accumulation, ice melting, water film dynamics, and refreezing, which can lead to ice accumulation at various locations on the airfoil surface.

During the cyclic de-icing process, ice accretion can occur in various locations due to the presence of water from melting ice and the runback water film. In the first 100 seconds, when only heater A is on, the ice accretion process begins on the airfoil surface, which has a low initial temperature. As shown in Figure 8; immediately after heater A is turned on, the surface temperature remains below the freezing point ($T_m$=273.15 K). It takes about 5 seconds for the surface temperature to reach $T_m$, causing the supercooled water droplets (initially) to freeze instantly, without any runback water. As time passes, the heat power from heater A gradually raises the surface temperature, melting the ice layer and preventing further freezing in the heated area. Eventually, the surface temperature over heater A rises above the freezing point, preventing any ice formation in this area and creating a thin water film with a thickness of just a few microns. As a result, the volume of runback water generated increases, and it flows toward heaters B and C, leading to continued ice accumulation over heaters B and C (as shown in Figure 10 & Figure 12). Between 10 and 70s, the width of the water film continues to expand (Figure 9). By 70s, the surface temperature differences stabilize, and the ice thickness growth becomes steady (Figure 10), indicating that the thermal equilibrium of the icing process is approaching a steady state. During this period, although the width of the liquid film remains unchanged, the ice thickness over heaters B and C continues to increase significantly (as shown in Figure 10 & Figure 12).

Additionally, since the initial temperature of the entire airfoil surface is set to the static ambient temperature at the start of the simulation, the boundary layer temperature around the wing gradually rises due to electro-thermal heating. As shown in Figure 8, the surface temperature in areas outside the influence of water droplets and heaters (surfaces with no heater underneath) increases gradually over time due to convective heat transfer from the air, as also noted in the study by [40].

After 100 s, when heaters B and C are turned on (see timing diagram Figure 11), their heat raises the outer surface temperatures, initiating the melting of ice accumulated during the first 100 seconds over these heaters. As shown in Figure 9 & Figure 12, the width of the water film expands along the airfoil chord over areas D and E, indicating that the ice on heaters B and C is melting. At 110 seconds, heaters D, E, F, and G are turned on while heaters B and C are switched off. Within a short time, the surface temperatures in the areas of heaters B and C drop to the freezing point (Figure 4). Meanwhile, the water film flows downstream over heaters D, E, F, and G between 110 and 120 seconds, covering a wide area on the airfoil surface, as illustrated in Figure 9 and Figure 12.

Figure 9. Water film thickness distributions at different times in the de-icing process

Figure 9. Water film thickness distributions at different times in the de-icing process

Since the ice on these heaters melts, the liquid film flows beyond the areas of heaters F and G and begins to freeze again in unprotected regions, forming ice ridges, as shown in Figure 10 & Figure 12. Regarding the formation and movement of the water film, it is assumed that the water layer fully wets the surface. However, the continuous flow of the water layer outside the impingement limit may break into small rivulets or even separate droplets due to hydrodynamic instabilities in the water film [45].

3.4. Analysis of Electro-Thermal De-Icing System at Various Angles of Attack

As shown in Figure 13, for the airfoil surface related to heater A, for small angles of attack, ice accretion on the leading edge of the airfoil occurs only during the first 5 seconds when the surface temperature has not yet reached the ice melting temperature ($T_m=273.15K$). After this initial period, due to the continuous operation of heater A, the surface temperature rises above $T_m$, causing the average ice thickness to nearly reach zero. However, for larger angles of attack, such as $8^{\circ }$ or ${12}^{\circ }$, some ice accretion still occurs on the leading edge of the airfoil despite the heater being continuously active. Specifically, for $8^{\circ }$, ice accretion persists up to 100 s, while for ${12}^{\circ }$, it persists throughout the entire period of 600s. At higher angles of attack, the thermal load on the upper leading edge increases beyond the heater capacity. This indicates that even when the heater operates continuously, at ${12}^{\circ }$, the system is not able to fully melt the whole ice in this region. This can be interpreted as follows: as angle of attack ($\alpha$) increases, higher local velocities raise the convective heat-transfer coefficient $c_h$, enhancing cooling. The local energy balance, $q_{\mathrm{heater}}^{\prime \prime }<c_h(T_s-T_{\infty })+{\dot{m}}^{\prime \prime }\left[c_p(T_{\infty }-T_s)+L_f\right]$, may therefore not be satisfied at $\alpha ={12}^{\circ }$, so $T_s$ does not remain above $T_m$ and leading-edge ice persists throughout the 600 s window. At $\alpha =0^{\circ }$, the same heater overcomes these losses, and the surface becomes essentially ice-free after the first few seconds. Although the average ice thickness at $\alpha ={12}^{\circ }$ is small, even minor ice accretion can degrade aerodynamic performance and increase surface roughness.

Figure 10. Ice thickness distributions at different times in the de-icing process.

Figure 10. Ice thickness distributions at different times in the de-icing process.

Figure 11. Heater activation timing.

Figure 11. Heater activation timing.

For the surface above heater B, performance is best at $\alpha =0^{\circ }$, where the average ice thickness decays toward zero. Moving away from $0^{\circ }$ degrades performance—most notably at $4^{\circ }$—because the imposed thermal load during the accumulation phase of each 120 s cycle increases. In all nonzero angles, the first cycle shows the largest build-up, but subsequent cycles exhibit progressively lower peaks, indicating that the periodic heating is largely effective over time. Over heater C the system is effective at $\alpha =-4^{\circ },0^{\circ },4^{\circ }$: a sizable build-up appears in the first cycle, but the average thickness decreases as cycles proceed. At $\alpha =8^{\circ }$ and ${12}^{\circ }$ the trend reverses. Early cycles show modest accumulation, yet the average thickness grows in later cycles. This is consistent with runback water produced upstream (primarily above heater A) re-freezing near the heater-C footprint, where the available power density is insufficient to fully melt it during the cycle. A thin residual layer then remains at the end of each cycle and compounds over time. Improving high-$\alpha$ performance will require re-balancing power or modifying the heater layout in this region (e.g., extended length or higher local heat flux).

Figure 12. Film of water, ice profile, and temperature contours at different times. To improve visibility, the height of the water film and ice have been increased by a factor of 200 and 2, respectively.

Figure 12. Film of water, ice profile, and temperature contours at different times. To improve visibility, the height of the water film and ice have been increased by a factor of 200 and 2, respectively.

Above heater E, performance is weak at $\alpha =0^{\circ }$ and $4^{\circ }$; for higher angles the instantaneous accretion rate is lower than at small $\alpha$, but the average thickness still increases monotonically over time. Only at $\alpha =-4^{\circ }$ is a net reduction observed. Because this zone lies farther downstream on the suction surface, its aerodynamic penalty is smaller than that of the leading-edge regions, yet the persistent growth indicates that additional coverage or power would be beneficial. Finally, above heater G, higher angles ($4^{\circ },8^{\circ },{12}^{\circ }$) yield negligible accretion compared with $\alpha =-4^{\circ }$ and $0^{\circ }$. Some ice does accumulate over time, but the amounts remain small; at $\alpha =-4^{\circ }$, however, the end-of-cycle thickness is appreciable. For the thermally unprotected airfoil upper surface, which is located downstream of heater G, no ice accumulates during the de-icing cycle at an angle of attack of $0^{\circ }$. However, at other angles of attack, especially at $-4^{\circ }$, the amount of accumulated ice increases over time as expected, due to the increase in the runback water film. Nevertheless, the average ice thickness formed in this region remain small and aerodynamically insignificant over the simulated duration.

3.5. Analysis of Electro-Thermal De-Icing System in Harsh SLD Condition

Figure 14 compares normal icing (monodisperse) with SLD conditions represented by Langmuir D distribution (MVD=50 $\mu$m, 100 $\mu$m, 200 $\mu$m). Above heater A (continuous operation), only a brief < 5 s accretion occurs before the surface exceeds $T_m$. By 600 s the leading-edge average thickness is near zero in both normal and SLD conditions, indicating robust protection at this critical location. Because of the asymmetric internal layout, heaters B and C do not behave identically. Over heater B the lowest performance occurs for MVD=50$\mu$m: the average ice thickness remains roughly constant through time, which is undesirable given its proximity to the leading edge. For MVD $=200\mu$m, a large initial build-up appears in the first ∼ 100 s, but subsequent cycles reduce the average thickness to near zero by the end of the simulation. For the airfoil surface associated with heaters C the behavior for MVD $=200\mu$m is similar—high early accumulation followed by gradual removal to near zero. Unlike heater B, heater C handles MVD $=50\mu$m more effectively, reducing the average thickness throughout the cycles. Given that heater C covers part of the suction surface, this improved response is particularly important for aerodynamic performance. However, on heater E the response to MVD $=20\mu$m, $50\mu$m, and $100\mu$m is unsatisfactory: the average thickness remains nearly constant or increases with time. Although this region is farther from the leading edge (so its aerodynamic impact is smaller than for A/B/C), the monotonic growth indicates that additional coverage or power reallocation would help. On heater G the system performs well for MVD $=20\mu$m (normal icing conditions) and $50\mu$m, but poorly for $100\mu$m and $200\mu$m, where the average thickness is nearly steady or increasing.

In the thermally unprotected zone downstream of heater G, the response is acceptable for $20\mu$m, $50\mu$m, and $100\mu$m: after ∼ 100 s, runback film entering this region causes a small increase in average thickness, but the values remain low. In stark contrast, $200\mu$m droplets produce a sustained increase throughout the cycle, posing a clear hazard. With no installed heaters, ice can grow to significant thickness and mass, raising the risk of ice shedding and potential damage to downstream surfaces. The results argue for chordwise coverage (more/longer heaters) and/or higher local power density for electro-thermal de-icing in harsh SLD conditions. Consequently, alternative or hybrid strategies, such as combining passive and active methods, should be considered for better performance.

Figure 13. Comparison of average ice thickness at different angles of attack, Case 1 to Case 5.

Figure 13. Comparison of average ice thickness at different angles of attack, Case 1 to Case 5.

Figure 14. Comparison of average ice thickness at different supercooled droplet sizes: Case 1 and Cases 6 to 8.

Figure 14. Comparison of average ice thickness at different supercooled droplet sizes: Case 1 and Cases 6 to 8.

4. Conclusion

This study performed a comprehensive simulation of ice accretion and electro-thermal de-icing on an airfoil under various icing conditions. Accretion was first modeled for 10 minutes using a multi-shot approach that captured geometry changes due to ice build-up. The resulting aerodynamic penalties were significant: at an angle of attack 12°, the lift coefficient decreased by over 10% and the drag coefficient increased by almost 50%. Under SLD conditions with a median volume diameter (MVD) of 200 $\mu$m, the drag coefficient increased by more than 30%. These results confirm the detrimental impact of icing on performance and underscore the need for effective de-icing strategies.

A cyclic electro-thermal de-icing system (five 2-minute cycles) was then simulated. The model, which shows good agreement with available experimental and numerical data, reproduced heat-transfer processes during freezing, melting, and evaporation. When the leading-edge heater was activated during supercooled droplet impingement for less than 5 s, a thin ice layer formed due to subfreezing surface temperatures but subsequently melted. As the cycle progressed, a water film developed and remained stable for approximately 70 s until downstream heaters are activated, after which runback water refroze in unprotected regions, forming ice ridges. While the baseline assumed monodisperse droplets of (20 $\mu$m), larger SLD (>40 $\mu$m) would alter collection efficiency and impingement behavior, highlighting the sensitivity of de-icing outcomes to droplet size.

Performance varied with angle of attack and SLD conditions as a function of heater placement. At small angles of attack, leading-edge accretion was minimal under continuous operation of Heater A, with ice thickness approaching zero after an initial transient. At higher angles, however, ice accumulation persisted despite continuous heating, indicating limitations in fully eliminating accretion in critical aerodynamic regions. Under SLD conditions, Heaters A and C were relatively effective, whereas Heaters B, E, and G were less effective, particularly for larger droplets of over 50 $\mu$m. For 200 $/mu$m droplets, ice thickness increased substantially in thermally unprotected areas—an operational hazard due to both ice accumulation and potential shedding. These observations point to the value of additional or longer heaters to extend thermal protection under challenging conditions. Finally, as a prospective alternative and complement to electro-thermal configurations (revised heater sizing and placement), passive icephobic coatings should be explored to improve runback management and reduce refreezing.