2.1 Numerical model

The proposed numerical model can determine the leading-edge ice thickness and successively find out the distorted aerofoil shape. Initially a control volume as Fig. 1 was considered onaerofoil surface to proceed with the model.

Figure 1. Mass and energy conservation within a control volume on aerofoil surface.

Within the control volume mass and energy are conserved as per the following equations [Reference Son, Oh and Yee15].

(1) \begin{align}{\dot m_{com}} + {\dot m_{in}} = {\dot m_{ice}} + {\dot m_{out}} + {\dot m_{eva}}\end{align}

(2) \begin{align}{\dot e_{com}} + {\dot e_{in}} = {\dot e_{ice}} + {\dot e_{out}} + {\dot e_{eva}} + {\dot e_{conv}}\end{align}

Where,

${\dot m_{com}}$ refers to the mass of impinging water that inflows onto the aerofoil surface.

${\dot m_{in}}$ is the mass of runback water towards the control volume that comes from the previous control volume.

${\dot m_{ice}}$ is the ice mass accreted on the surface.

${\dot m_{out}}$ is the runback water from the control volume under consideration.

${\dot m_{eva}}$ is the evaporated ice mass from the control volume.

Similarly, energy terms are applicable for the energy equations. All these terms in Equations (1) and (2) can be determined individually. As it is not possible to determine impinging water mass directly, it is calculated indirectly using collection efficiency [Reference da Silveira, Maliska and Estivamand Rafael Mendes26]. Where, collection efficiency is the ratio of perpendicular distance Δy between particles and distance ΔS on the surface where the impinging particles may stick Fig. 2.

Figure 2. Trajectory of the water droplet impinging on the surface (close up view).

Thus, the mass of liquid water impinging on the surface considering the particles are uniformly distributed in the air can be determined by the following equation.

(3) \begin{align}{\dot m_{com}} = \beta .LWC.{V_\infty }.\Delta S\end{align}

Where,

${V_\infty }$ is the free stream air velocity.

$\beta = \frac{{{\rm{\Delta }}y}}{{{\rm{\Delta }}S}}$ is termed as the collection efficiency or ice accumulation rate or catch rate on the surface.

$LWC$ is the liquid water content in the air, which depends on the temperature condition and cloud type.

As the calculation started the initial ${\dot m_{in}}$ is considered to be zero as there is no previous control volume to transfer mass. Evaporated water mass can be determined as per the following equation as it is the function of surface temperature [Reference Ozgen and Cambek27].

(4) \begin{align}{\dot m_{eva}} = \frac{{0.7}}{{{C_{Pa}}}}{h_c}\left( {\frac{{{P_{v,sur}} - {P_{v,\infty }}}}{{{P_\infty }}}} \right)\end{align}

Where,

${C_{Pa}}$ is the specific heat of air.

${h_c}$ is the convective heat transfer coefficient and it can be determined from the Nusselt number according to the following equation [Reference Anderson and Tsao9].

(5) \begin{align}{h_c} = \frac{{{k_a}{N_u}}}{d}\end{align}

Where, d measures the leading-edge diameter of the aerofoil. The stagnation line heat transfer coefficients of an NACA 0012 aerofoil were measured by Poinsatte [Reference Poinsatte28] at the IRT using the chord length as the characteristic length. In this study leading edge diameter of the aerofoil was considered as the characteristic length to cope up with the practice and Poinsatte expression can be re-written as per the following equation.

At 0° AoA,

(6) \begin{align}{N_u} = 1.10{R_e}^{0.472}\end{align}

Where,

${P_{v,sur}}$ and ${P_{v,\infty }}$ are the vapour pressure at the ice or water surface respectively for ambient air and can be computed from the following equations(Reference da Silveira, Maliska and Estivamand Rafael Mendes26).

(7) \begin{align}{P_{v,}} = 3386(0.0039 + 6.8096 \times {10^{ - 6}}{\bar T^2} + 3.5579 \times {10^{ - 7}}{\bar T^3})\end{align}

(8) \begin{align}\bar T = 72 + 1.8\left( {T - 273.15} \right)\end{align}

Ice mass and runback water mass can be determined using the freezing fraction which is the ratio of the freezing ice mass to the total ice mass as studied by Messinger [Reference Son, Oh and Yee15].

(9) \begin{align}f = \frac{{{{\dot m}_{ice}}}}{{{{\dot m}_{com}} + {{\dot m}_{in}}}}\end{align}

(10) \begin{align}{\dot m_{ice}} = f\left( {{{\dot m}_{com}} + {{\dot m}_{in}}} \right)\end{align}

(11) \begin{align}{\dot m_{out}} = \left( {1 - f} \right)\left( {{{\dot m}_{com{\rm{\;\;}}}} + {{\dot m}_{in}}} \right) - {\dot m_{eva}}\end{align}

From the energy conservation Equation (2), the energies coming from the air which is adjacent to the control volume and moves to the next control volume can be calculated using the following equations.

(12) \begin{align}{\dot e_{com}} = {\dot m_{com}}\left[\!\left( {{T_{sur}} - {T_i}} \right)\frac{{{V_\infty }^2}}{2}\right]\end{align}

(13) \begin{align}{\dot e_{in}} = {\dot m_{in}}\left[{c_{p,w}}\left( {{T_{sur\left( {I - 1} \right)}} - {T_i}} \right)\right]\end{align}

(14) \begin{align}{\dot e_{eva}} = {\dot m_{eva}}{[_{}}{c_{p,w}}\left( {{T_{sur}} - {T_i}} \right) + {L_v}]\end{align}

(15) \begin{align}{\dot e_{conv}} = {h_c}{[_{}}{T_{sur}} - \left( {{T_e} + \frac{{{r_h}{V_\infty }^2}}{{2{c_{p,a}}}}} \right)]\Delta S\end{align}

(16) \begin{align}{\dot e_{ice}} = {\dot m_{ice}}{[_{}}{c_{p,i}}\left( {{T_{sur}} - {T_e}} \right) - {L_f}]\end{align}

(17) \begin{align}{\dot e_{out}} = {\dot m_{out}}{[_{}}{c_{p,w}}\left( {{T_{sur}} - {T_i}} \right)\!]\end{align}

Where,

${c_{p,w}}$ and $\;{c_{p,i}}$ are the specific heat of water and ice, respectively [Reference Touloukian, Liley and Saxena29].

${L_{f}}$ and ${L_v}$ are the latent heat of fusion of ice and latent heat of evaporation of water, respectively.

While proceeding for the calculation all the terms of Equation (1) are substituted in Equation (2). The constant term values are obtained from the ambient conditions [30]. Thus, the final equation is left with only two unknown terms, surface temperature ${T_{sur}}$ and the freezing fraction, f. Finally, these two quantities are now calculated using an iterative process.

For initial iteration, the surface temperature is considered to be at 273.15K (freezing point) and freezing fraction f is calculated. If f< 0, it implies that the surface temperature is greater than freezing point, so f is set to be 0 and ${T_{sur}}$ is recalculated. On the other hand, if f >1, it means the ${T_{sur}}$ is less than 273.15K, so f is set to be 1 and ${T_{sur}}$ is calculated again. If two successive iterations give the ${T_{sur}}$ value close to each other, then only ${T_{sur}}$ and f are known and ${\dot m_{ice}}$ , are found using equations 1, 3, 4, 10 and 11 for each of the control volume.

Next step in the process is the calculation of the height of ice thickness for each surface panel of length $\Delta S$ . Height of ice growth on each surface panel is calculated based on the ice mass calculated in the previous step during a small-time interval Δt, using the following equation [Reference Son, Oh and Yee15].

(18) \begin{align}h = \frac{{{{\dot m}_{ice}}.\Delta t}}{{{\rho _{ice}}.\Delta S}}\end{align}

Where, ${\rho _{ice}}$ is the density of ice.

Assuming that ice grows only in the perpendicular direction to the aerofoil surface, these ice growth heights are added normal to the surface panels to produce a new profile. The ice growth becomes more as the time interval increases, and thus distorted aerofoil shape is generated due to ice accretion.

2.2 Computational fluid dynamics (CFD) model

An efficient design of aircraft anti-icing and de-icing devices and certificate of airworthiness requires the prediction of in-flight ice thickness and aerodynamic performance in varying ambient conditions. CFD offers the advantage of low cost and the ability to simulate realistic conditions. But the accuracy depends on the quality of the grid system [Reference Chi and Slater31] in resolving the geometry and hence a mesh validation was done comparing the experimental results for the proposed model. In this section advantages and disadvantages of various turbulence models are described and reasons for selecting k- $\varepsilon$ model for the present study are identified [Reference Shelil32, Reference Wilcox33].

1. (1) RANS Single equation model: SpalartAllamaras

If we want to observe the flow around an aerofoil or aircraft, Spalart Allamaras would be a great choice because it is well known and tested for this sort of application. but if we want to design further and determine the angle-of-attack at which the aircraft stall, it would be no longer a good choice.

1. (2) RANS two equation model:

   1. (a) Standard k- $\varepsilon$ model

      The standard k- $\varepsilon$ model is based on transport equations for the turbulence kinetic energy (k) and dissipation rate ( $\varepsilon$ ). For deriving the k- $\varepsilon$ the assumption is that the flow is fully turbulent and the effects of molecular viscosity are negligible. This model uses wall functions to analytically account for the fluid velocity in the viscous sub layer near the wall. This model offers good convergence and is not memory intensive and it is typically used for external flow with complex geometries.

1. (b) Realizable k- $\varepsilon$ model

   This model modifies the equation for epsilon and introduces the effect of mean flow distortion on turbulent dissipation. The model has improved performance for planar surfaces, round jets, rotation, recirculation and streamlines curvature. It can improve the boundary layer under adverse pressure gradients but can not do a magic as it’s still based on turbulent viscosity.

2. (c) RNG k- $\varepsilon$ model

   Results from the updated RNG turbulence model produces lower turbulence levels and can underestimate the value of k but it offers little or no advantage over the realisable k- $\varepsilon$ . So, the model is not convincible to use for complex geometries like ice accreted aerofoils.

3. (d) Standard k- $\omega$ model

   It is a popular two-equation model that pairs the k with specific rate of dissipation of kinetic energy ( $\omega$ ). It can model near wall interactions accurately than k- $\varepsilon$ models but has difficulty of convergence than k- $\varepsilon$ model. The model is very sensitive to inlet boundary conditions, which is disadvantage of the model compared to the k- $\varepsilon$ models.

4. (e) SST k- $\omega$ model

   This model is an enhancement of the original k- $\omega$ model and addresses the issue of sensitivity of freestream turbulence levels. It has superior performance in simulating the boundary layers with adverse pressure gradients but the performance is not very different from the realisable k- $\varepsilon$ model.

5. (f) Reynolds Stress Model (RSM)

   This is the most complete representation of physical turbulent flow. These models attempt to model the flow and terms directly in RANS equations. These models are based on six equations which represent turbulent stresses. These models can represent flow very well but are typically used for flows with extremely complex geometries or flows which have not been studied before because of their computational expense, sensitivity to initial conditions, amount of modeling and high-quality mesh requirement.

   Both steady and unsteady Spalart Allamaras, shear-stress transport, RNG k- $\varepsilon$ turbulence model was used by Y. Li et al. [Reference Chi and Li34] to predict the lift, drag and moment coefficients of a business jet aerofoil with rough and jagged rime ice and glaze ice in 2005 and it was reported that unsteady simulations did not provide better results whereas steady RANS, WIND and FLUENT simulations gave almost identical results for the same grid and order of accuracy. CFD solutions of NLF0414 and B575/767 aerofoil were generated with and without ice shape to demonstrate the developed model and grid by X. Chi et al. in 2002 for conditions like Mach number 0.29, T _s = −5 $^{\circ}$ C, R _e = 6.4 × 10⁶ with no-slip boundary conditions and no wall function [Reference Chi and Slater31]. Three different turbulence models as SpalartAllamaras, low Reynold number k- $\varepsilon$ model and SST k- $\omega$ model were explored by Kim et al. in 2013 [Reference Kim, Brown and Kreeger35] for aero-thermodynamic modeling of aircraft ice accretion phenomenon and they found acceptable correlation for all turbulence models at low angle-of-attack but there was noticeable differences in stall prediction and convergence to the steady-state among the three solvers. 2D numerical simulation study of DTU-LN221 aerofoil performance was studied by N. Shelil in 2020 [Reference Shelil32] using seven different CFD models such as RANS, Spalart Allmaras, standard k- $\varepsilon$ , RNG k- $\varepsilon$ , standard k- $\omega$ , SST k- $\omega$ and RSM and he reported that there was no general model that can be called the best over the range of study. Aerodynamic performance of NACA0012, NACA4412 and JOUKOWSKI (T = 12%) aerofoil, at various angle-of-attack (from 2 deg to 18 deg) and high Reynolds number were studied using realisable k-epsilon turbulence model and were compared with experimental results [Reference Jain and Mohammad36]. Results of the study identified that NACA 4412 gives better aerodynamic performance compared to NACA 0012 and JOUKOWSKI (T = 12%) aerofoil. Atmospheric ice accretion for different NACA series like NACA 0012, NACA 23012, NACA 643218 was studied by Talasila et al. [Reference Talasila and Ravi37] using a 2D model developed with solid works and the aerodynamic performance was investigated using standard k- $\varepsilon$ CFD model. They did the study to determine the NACA series, which is more efficient and resistant to icing. Results of their study found that six-series NACA aerofoil gives better efficiency at icing conditions. N.A. Khan et al. computationally studied the performance of a three-blade rotor of diameter 0.5m in ANSYS fluent and reported that the two equations realisable k- $\varepsilon$ turbulence model predicted the rates of round jets more accurately. Moreover, it has a great performance for flows involving separation, recirculation and boundary layer with strong adverse pressure gradient [Reference Khan and Islam38].

Studying the above literature and considering all advantages and disadvantages the CFD model used in this study is based on the RANS (Reynolds Average Navier Stoke) two-equation based standard k- $\varepsilon$ turbulence model for observing the flow filed around the aerofoil geometry. Much attention was given in producing a better grid system as the produced aerofoil shape was much complex than the base aerofoil. In the present study a 2D pressure-based (as it is used for small velocity up to 0.7M and in this study the maximum Mach number was 0.48M) turbulence model was considered.

1. (1) Governing equations

The motion of viscous fluid substances is described by the Navier-Stokes equation which is the governing equation of CFD. The continuity equation and momentum equation can be derived applying the mass and momentum conservation [Reference Zou39].

1. (a) Continuity equation

(19) \begin{align}\frac{{D\rho }}{{Dt}} + \rho \frac{{\partial {U_i}}}{{\partial {x_i}}} = 0\end{align}

1. (b) Momentum equation

(20) \begin{align}\rho \frac{{\partial {U_j}}}{{\partial t}} + \rho {U_i}\frac{{\partial {U_J}}}{{\partial {x_i}}} = - \frac{{\partial P}}{{\partial {x_j}}} - \frac{{\partial {\tau _{ij}}}}{{\partial {x_i}}} + \rho gj\end{align}

Where,

1^st term: Local change with time

2^nd term: Momentum convection

3^rd term: Surface force

4^th term: Molecular-dependent momentum exchange (diffusion)

5^th term: Mass force

And

(21) \begin{align}{\tau _{ij}} = -{{ {\rm{\mu }}\left( {\dfrac{{\partial {U_j}}}{{\partial {x_i}}} + \dfrac{{\partial {U_i}}}{{\partial j}}} \right) + 2/3{\delta _{ij}}\mu \dfrac{{\partial {U_k}}}{{\partial {x_k}}}}}\end{align}

For the purpose of simplification, the Navier-Stokes equations is

(22) \begin{align}\frac{{\partial \left( {\rho \varphi } \right)}}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}\left( {\rho {U_{i}}\varphi - {{\rm{\Gamma }}_\varphi }\frac{{\partial \varphi }}{{\partial {x_i}}}} \right) = {q_\varphi }\end{align}

1. (2) Assumptions

For the CFD simulations using the standard k- $\varepsilon$ turbulence model following assumptions were made:

1. (a) The flow is fully turbulent and the effects of molecular viscosity are negligible.

2. (b) Standard atmospheric pressure of 1.225atms was considered in order to make the wake effect zero.

3. (c) 6^th order of accuracy O(6) was considered for the computation.

   (23) \begin{align}{\varepsilon ^n}^{ + 1}-{\varepsilon ^n} \lt = {10^{ - 6}}\end{align}

4. (d) Standard wall function was considered.

5. (e) No slip boundary conditions were considered with roughness height = 0 and roughness constant = 0.5

6. (f) SIMPLE and 2^nd order upwind solution method was utilised for the solution.

1. (3) Geometry

   

   Figure 3. 2D aerofoil geometry with a selected C domain.

For the CFD computation the fluid domain has been selected by exhausting the periphery of the flow field. A 2D aerofoil geometry was considered with a C domain which is extended up to 12.5C at the radius and 25C (C is the chord of the aerofoil) from the aerofoil trailing edge to observe the flow at the vicinity and far away from the aerofoil as shown in Fig. 3.

1. (4) Mesh

   

   Figure 4. Generated unstructured mesh in the computational C domain.

   

   Figure 5. Inflation layers around the edges of the iced aerofoil (asymmetric NACA 2412 aerofoil).

   

   Figure 6. Close up view of the inflation layers around the edges (asymmetric NACA 2412 aerofoil).

The fluid domain was discretised into smaller elements based on a specific size. The unstructured mesh was applied in the fluid domain for better accuracy as the iced aerofoil has a complex geometry. A fine mesh was obtained with a total of 96,247 elements and 96,899 nodes with minimum and maximum cell sizes 0.00657 and 0.50m, respectively. The appropriate number of nodes can be determined by increasing the number of nodes so that further refinement does not change the results. So, a mesh validation was done and described in section IIIC of the present study to observe the effect of mesh quality. A body of influence with radius 3m, element size 5e-002m and local minimum size of 6.57e-003m was considered around the aerofoil. Smooth transition inflation was applied with maximum 10 layers, transition ratio of 0.272 and growth rate 1.2 for greater accuracy and better flow visualisation. The overall mesh and close-up mesh with inflation layer used for CFD study are presented in Figs. 4, 5 and 6, respectively.

1. (5) Boundary conditions

   

   Figure 7. Selection of domain and boundary conditions for the computational analysis.

   Table 1. Selected boundary conditions for ANSYS implementation

   

For the flow domain the applied boundary conditions are represented by Fig. 7 where the inlet velocity with free stream velocity is specified as inlet boundary and the x and y component of velocity is determined using the angle-of-attack. The upper and lower extents of the domain are termed as velocity inlet boundary conditions. A wall boundary condition is applied to the aerofoil and any ice accreted on aerofoil leading edge where no-slip boundary condition is applied and pressure outlet boundary condition is applied to the outlet presented by Table 1.

3.1 Ice accretion model validation

3.1.1 Validation for symmetric aerofoil

An experimental study was carried out by Hans and Palacios [Reference Hans and Palacios40] over an NACA 0012 aerofoil where two blade configurations of helicopter rotor were considered. Varying LWC, icing time, temperature and density ice accretion was considered over blades having a diameter of 2.7m and chord length of 0.533m presented by Fig. 8 and rotor rpm was varied between 200 to 1,200rpm. The temperature was maintained constantly between 0°C to −25°C inside a cold chamber at the test facility. The proposed ice accretion model had considered the following two experimental conditions for an NACA 0012 aerofoil as mentioned by Hans and Palacios [Reference Hans and Palacios40] listed in Table 2.

Figures 9 and 10 represent the comparison of the generated ice accreted aerofoil using the proposed numerical model with the experimental ice profiles for the above conditions. Though the numerically generated ice thickness almost appears to be the same with the experimental one, the leading-edge ice accretion for the proposed numerical model shows a slight deviation from the experimental ice shape referred from the open literature, which measures 6.26% and 3.89% for icing conditions 1 and 2 respectively at the stagnation point presented by Table 3.

Table 2. Icing conditions for the ice-shape models mentioned in Ref. [30].

Figure 8. Rotor-icing-test stand with the test blade of the experimental set up from reference literature [30].

Figure 9. Test Condition 1: Comparison between experimental ice shape and numerical ice shape from proposed method for symmetric aerofoil.

Figure 10. Test Condition 2: Comparison between experimental ice shape and numerical ice shape from proposed method for symmetric aerofoil.

3.1.2. Validation for asymmetric aerofoil

Ice accretion and aerodynamic performance analysis of a NACA 5-digit asymmetric aerofoil was validated with the experimental results studied by Lee et al. [Reference Lee and Addy41] in 2016 on a 2D NACA 23012 aerofoil (test conditions listed in Table 4) as a joint effort by NASA, ONERA and the university of Illionois. The objective of their research was to develop a data base for ice accretions on small-scale, large-scale 2D aerofoil and data base for the ice accretion effects on aerodynamics.

Figure 11 describes the comparison of the ice shape obtained from the present model with the experimental ice shape for a sub-scale model of NACA 23012 aerofoil referred by Lee [Reference Lee and Addy41]. Results of the study shows very close result for the stagnation point (4.51%) and at other points at the leading edge the difference was small. But significant deviation was observed at some other points where the numerical model overestimated the ice accumulation. However, as the model can estimate ice accretion close to the experimental data for symmetric aerofoil and ice thickness with minor difference for asymmetric (5-digit aerofoil), the model will be now applied to predict the ice accumulation over asymmetric (4-digit) aerofoil.

Table 3. Ice-shapes comparison with the reference literature

Table 4. IRT sub-scale model test conditions from Ref. [Reference Lee and Addy41].

Figure 11. Comparison between experimental ice shape and numerical ice shape from proposed method for asymmetric aerofoil.

3.2 Validation of the CFD model

3.2.1 For clean aerofoil (NACA 2412)

Shubham et al. in 2021 [Reference Shubham and Arnav42] had studied the performance of different asymmetric aerofoils NACA 2412, NACA 2414 and NACA 2415 for selecting the appropriate turbulence model for 2D simulation. They had investigated the aerodynamic performance of these aerofoils using standard k- $\varepsilon$ , standard k- $\omega$ , RST and transition k-kl-ω models (parameters of CFD model are listed in Table 5). Results of their study showed that the standard k-e model gives good efficiency for larger range of angle-of-attack (5–20 deg) while all other models had good efficiency from 2 to 5 deg only. So, the standard k- $\varepsilon$ CFD model utilised in this study was validated for NACA 2412 (which is our aerofoil of interest) and compared with this literature.

Table 5. Parameters for validation of the CFD model with asymmetric aerofoil

Figure 12. Comparison of lift force as a function of angle-of-attack for the present CFD model and from reference literature.

Figure 13. Comparison of drag force as a function of angle-of-attack for the present CFD model and from reference literature.

Figure 14. Wind-tunnel test section with aerofoil mounted for the experimental set up from reference literature [30].

Figure 15. Comparison of coefficient of Lift as a function of angle-of-attack for experimentally and numerically obtained ice shape with the base aerofoil.

Figure 16. Comparison of coefficient of drag as a function of angle-of-attack for experimentally and numerically obtained ice shape with the base aerofoil.

Figure 17. Dependency of mesh quality on the computational results.

Figures 12 and 13 describes the lift and drag force as a function of angle-of-attack for clean aerofoil (NACA 2412). From the graph of lift force for a range of angle-of-attacks from 0 deg to 22 deg shows that the present model shows only a small deviation at 8 deg angle-of-attack and at the stall angle-of-attack. For the present model the stall angle was 20 deg whereas for the reference literature stall angle was at 21 deg. The drag versus angle-of-attack graph shows almost similar nature except few places for both present model and reference literature. This validates that the present CFD model can be utilised for studying the aerodynamics of asymmetric iced aerofoil.

3.3. Mesh validation

With a view to observing the effect of mesh quality on the computational results, four different mesh qualities were considered. For a full scale of angle-of-attack, the lift coefficient was plotted and from Fig. 17, it was observed that with the increase of grid numbers, the lift curve almost converges to the experimental value. When an angle-of-attack of 16 deg was considered, it was observed that for mesh quality 2 and 4 the difference with the experimental result was minimum and was calculated as 3.09% and 3.035% respectively. Hence, for the proposed model, any one of the mesh qualities can be used. However, second grid numbers (96,000) were used throughout the study considering the computational time.

4.1 Solution

Ice accretion prediction and the aerodynamic effects of icing for symmetric aerofoils like NACA 0012 and NACA 0015 are widely available in the open literature both in form of experimental and numerical analysis. But much study on the asymmetric or cambered aerofoil had not been carried out to the best of author’s knowledge. But subsonic flight encompasses a pretty large range of airspeed/ airflow and NACA 2412 is a commonly used aerofoil for subsonic flight and the results are quite typical of aerofoil characteristics and it’s the simplest four-digit asymmetric aerofoil. Hence, in this study NACA 2412 was the chosen asymmetric aerofoil in order to observe the ice accretion phenomenon over cambered aerofoil. The main focus is to predict the ice accumulation over an asymmetric NACA 2412 (chord length of 1m) aerofoil and evaluate the aerodynamic performance degradation based on the aerofoil geometry, ambient conditions and flight conditions. Hence, it considers moderate to high altitude flights ranging from 10,000 to 25,000ft where air temperature varies from −10°C to −40°C and most of the icing phenomenon happens.

1. (1) Assumptions

To simplify the problem and ease of calculation the following assumptions were made:

1. (a) The freezing point of water and melting point of ice is the same.

2. (b) The aircraft surface or skin temperature is at freezing temperature for icing to be occurred.

3. (c) The analytical freezing fraction value ranges from 0.3 to 1.0 from Messinger’s heat balance analysis for effective ice accretion.

4. (d) Air flow around the aerofoil surface is considered to be turbulent and compressible.

5. (e) Ice accretion is considered for an aerofoil of a high-performance subsonic aircraft.

4.2 Analysis of icing parameters

For analysing the effect of each of the above-mentioned parameters, three reference cases were considered by changing each of the parameters at the same Reynolds number (Re = 6 $ \times {10^{6{\rm{\;\;}}}})$ . For different cloud conditions, the parameters were selected from the established standards and experimental data available from open literature.

Table 6. Ambient conditions for case I

1. (1) Case I

A high performance, subsonic aircraft moving at an airspeed of 122m/s within an altitude range of 12,700ft (3,898m approx.) through medium level Nimbostratus clouds [30]. At this altitude air contains a low amount of LWC, low catch rate but accumulation rate is higher as the air temperature is within −10°C. Other ambient conditions are listed in Table 6 according to the standard atmospheric chart (SI units).

As obtained from the standard atmospheric conditions, nimbostratus clouds have a horizontal coverage of several thousand feet (10–15NM) but small vertical coverage of about 3,000ft [43], flying within nimbostratus clouds will be affected by long duration of icing, which will result in rime ice accumulation. So, aircraft is supposed to be within icing clouds for about a duration of 227s (15NM/122m/s). Ice accumulation is calculated for the time duration and analytical freezing fraction values mentioned in Table 7.

1. (2) Case II

Table 7. Icing duration and intensity case

Table 8. Ambient conditions for case II

Table 9. Icing duration and intensity case

The very particular aircraft is now considered to be at an altitude of 17,000ft (5,395m approx) with a higher airspeed of about 142m/s. At this altitude medium level altocumulus clouds exist, which have with high concentration of SLD, higher catch rate and high LWC. This causes the worst icing case but for a short duration as cumulus clouds have limited horizontal coverage but extended vertical coverage of several thousand feet. Other ambient conditions at this altitude are listed in Table 8.

As a single cumulus cell extend for a horizontal range of 2NM to maximum 6NM [43], with an airspeed of 142m/s the aircraft is supposed to be within icing cloud for a duration of 80s (6NM/142m/s) on an average. So, icing duration will be longer if the aircraft is to fly through a number of cumulus cells contributing to the most critical icing conditions. So, ice prediction time and freezing fraction values are considered as per Table 9.

1. (3) Case III

If the same aircraft is now considered to fly at a higher altitude range of 22,000ft (about 6,800m) where the air temperature is within −30°C to −40°C and icing threat is present due to the presence of cirrocumulus clouds (20,000ft to 40,000ft) containing high LWC and high catch rate contributing to the glaze ice or clear ice. According to the standard altitude chart other ambient conditions are enlisted in Table 10.

Table 10. Ambient conditions for case III

As cumulus clouds usually have limited horizontal coverage (a single cumulus cell is extended up to 2–6NM), it covers a small duration of the flight time. So, icing time can be 56s (6NM/167m/s) to more. Hence, ice accretion prediction was made for the following freezing fractions and time durations of Table 11.

Table 11. Icing duration and intensity case

4.3 Analysis of relative importance of form drag versus skin friction drag

Form drag on an aerodynamic body is due to the shape and size of the body. On the other hand, skin friction drag is more complicated and is related to the process of transition (when the laminar boundary layer becomes turbulent). A realistic drag analysis should always consider a mixed boundary layer. For a streamlined body like an aerofoil; almost 20% drag is form drag and 80% drag is from skin friction drag [Reference Anderson44]. Again, orientation of the body also has some effects on the relative proportions of the skin friction and pressure difference between the front and back portion. For a streamlined body like aerofoil; the drag is dominated by the skin friction drag.

The drag coefficient C _d is sensitive to Reynolds number which is to be expected since both skin friction and flow separation are viscous effects. The skin friction drag coefficient is a combination of laminar skin friction drag coefficient over the forward part and turbulent skin friction drag coefficient over the remaining part. As a result, the boundary layer always starts from some distance from the leading edge and then transits to a turbulent boundary layer at some point downstream of the leading edge. A reasonable quantitative breakdown between the skin friction and form drag can be found in the results of Lombardi et al. [Reference Lombardi45]. Here using an accurate CFD technique at a Reynolds number 3×10⁶ they had calculated a total drag coefficient for NACA 0012 to be 0.00623 and a skin friction drag coefficient to be 0.00534, indicating that the form drag due to flow separation is 15% and this is reasonable as the skin friction is expected to be 80% and form drag is expected to be 20% of the total drag for a stream-lined body. However, at higher Reynolds number and hence higher velocity, increasing angle-of-attack and less streamlined body (adding some thickness, cambered or iced aerofoil) the transition point tends to move forward increasing the skin friction drag. On the contrary, as the body becomes less streamlined (more like a blunt body), the pressure drag becomes the dominant factor with the increase of angle-of-attack.

Assuming the aerofoil with a given surface (ice accretion-one of the factors that affect the location of transition) is in a flow at a free stream velocity V _∞, we can determine the location of transition point with the help of critical Reynolds number using the following equations [Reference Lee and Addy41].

(24) \begin{align}R{e_{{x_{cr}}}} = \frac{{{\rho _\infty }{V_\infty }{x_{cr}}}}{{{\mu _\infty }}}\end{align}

Where the assumption of critical Reynolds number = 500,000 is very conservative and more likely the actual value is closer to 1,000,000. If the velocity is increased, the transition point tends to move forward to the leading edge. The Reynolds number based on chord length is given by the following equation.

(25) \begin{align}R{e_c} = \frac{{{\rho _\infty }{V_\infty }c}}{{{\mu _\infty }}}\end{align}

Combining Equation (24) and (25) we can write that,

(26) \begin{align}\frac{{R{e_{{x_{cr}}}}}}{{R{e_c}}} = \frac{{{x_{cr}}}}{c}\end{align}

This location is the transition from the leading edge relative to the chord length. As the skin friction drag coefficient (c _f ) is always based on this distance, we can measure it using the following equation [Reference Lee and Addy41].

(27) \begin{align}{C_f} = \frac{{{x_{cr}}}}{c}{\left( {{C_{f,1}}} \right)_{laminar}} + {\left( {{C_{f,c}}} \right)_{turbulent}} - \frac{{{x_{cr}}}}{c}{\left( {{C_{f,1}}} \right)_{turbulent}}\end{align}

(28) \begin{align} \textrm{or},\; {C_f} = \frac{{{x_{cr}}}}{c}\frac{{1.328}}{{\surd R{e_{{x_{cr}}}}}} + \frac{{0.74}}{{R{e_c}^{1/5}}} - \frac{{{x_{cr}}}}{c}\frac{{0.74}}{{R{e_{{x_{cr}}}}^{1/5}}}\end{align}

Taking into account both the surfaces of the aerofoil, the net skin friction drag coefficient will be

(29) \begin{align}{C_f}\left( {total} \right) = 2{C_f}\end{align}

The measured aerofoil drag coefficient C _d includes both skin friction drag and pressure drag due to flow separation. So, the pressure drag coefficient can be obtained mathematically by subtracting the skin friction drag coefficient from the total drag coefficient.

5.1 Ice thickness at the stagnation point for different cases

1. Case I

2. Case II

3. Case III

The mass rate of ice accretion on a structure is affected by various factors like shape of aircraft surface (aerofoil shape, leading edge radius and aerofoil symmetry), collection efficiency, freezing fraction, time of ice accretion, aerofoil surface temperature, airspeed and meteorological conditions like liquid water content (LWC), droplet size, atmospheric temperature etc [Reference Heinrich, Ross, Zumwalt, Provorse and Padmanabhan4]. In this study the obtained ice thickness was much greater than the reference ice shape of section 3.1 obtained from previous literature [30]. The reason of such differences lies in the selection of various parameters mentioned in Table 12.

Table 12. Parametric differences between reference literature [30] and present study

5.2 Ice accretion for the worst cases

Using the proposed numerical model the worst cases of calculated leading edge ice thickness values are enlisted in Table 13.

The above results show the worst icing cases at the stagnation point in term of maximum thickness for freezing fraction value of 1 for each of the cases as the freezing fraction value was varied from 0.3 to 1.00. Though case II generated maximum (0.1282m) leading edge ice accretion at the stagnation point for maximum (700s) time of ice accretion, the free-stream velocity was maximum for case III (167m/s). Hence the time of ice accretion demonstrates a significant influence on aerofoil icing. Therefore, the ambient conditions like density of air, density of water, viscosity of air might have very small influence over ice accretion at the same Reynolds number.

Table 13. Stagnation point ice thickness

Figure 18. Ice shape for case I (Re = 6 $ \times {10^{6}})$ .

Figure 19. Ice shape for case II (Re = 6 $ \times {10^{6}})$ .

Figure 20. Ice shape for case III (Re = 6 $ \times {10^{6}})$ .

Figures 18, 19 and 20 reports that the ice shape accumulated on the NACA 2412 aerofoil resembles a circular shape which is different in shape than the reference literatures from Refs 30 and 43. The reasons lie in the fact that ice accumulation over aerofoil surfaces can be predicted through one layer mode or multi-layer mode. In one layer mode ice is predicted in one step through the entire duration of time which obviously overestimates the runback water and predicts the observed horns. The present study calculated ice thickness in multi-layer mode where time of ice accretion was divided into small segments. In this study, the height of ice growth for each surface panel of length $\Delta S$ . During a small time interval Δt, the height of ice growth on each surface panel is determined based on the ice mass calculated in the previous step. The flow- field, droplet trajectory calculations were repeated for each layer. This approach considered the effects of ice shapes on flow field trajectory and reflected the physics of the problem more realistically. The present model also treated the cases with varying ambient and flight conditions.

5.3 Aerodynamic performance analysis

The coefficient of lift and coefficient of drag of iced aerofoil was compared with the base aerofoil using the k- $\varepsilon$ turbulent model of CFD analysis and the following results were found.

1. Case I

From the above plot, the iced aerofoil shows a maximum 84.02% decrease in lift coefficient and 262% increase in drag coefficient at 6 deg angle-of-attack compared to the base aerofoil as shown by Figs. 21 and 22 respectively. And the maximum difference in ice thickness is 20% between upper and lower surface. The average increase in drag coefficient for the iced aerofoil is 158.28% compared to the base aedrofoil.

Figure 21. Coefficient of lift as a function of angle-of-attack for numerically obtained ice shape compared with the base aerofoil.

Figure 22. Coefficient of drag as a function of angle-of-attack for numerically obtained ice shape compared with the base aerofoil.

The lift curve shows a fall in coefficient of lift at 6 deg AoA, again a rise at 8 deg AoA and a fall at 10° AoA and then constant beyond 10° AoA. Such aerodynamic performance degradation of the iced aerofoil can be explained observing the velocity vector.

The velocity vector of Fig. 23 for base aerofoil shows still attached flow at 6 deg AoA generating more lift than drag. But the velocity vector of Fig. 24 for the iced aerofoil at 6 deg AoA shows a reverse flow at the ice accreted surface for both upper and lower surfaces causing the velocity to reach almost zero at the boundary layers. This reverse flow resulted in a very low lift coefficient and higher drag coefficient compared to the base aerofoil. This indicates that the aerofoil approaches to stall condition at only 6 deg AoA.

Figure 23. Velocity vector showing attached flow for base aerofoil at 6 deg AoA.

Figure 24. Velocity vector showing reverse flow for iced aerofoil at 6 deg AoA.

1. Form drag versus skin friction drag

An analysis of the skin friction drag and form drag for the icing condition of case I (velocity 122m/s) reveals the fact that for the base aerofoil of Figs. 25 and 26 and at lower angle-of-attack the drag was dominated by the skin friction drag which was about 65% of the total drag. And the skin friction drag tends to increase with the angle-of-attack but at a negligible amount whereas the form drag increases drastically as the body becomes less streamlined. At 16 deg angle-of-attack the drag was mostly form drag contributing 95% to the total drag.

Figure 25. Skin friction coefficient and form drag coefficient for with respect to the total drag coefficient for base aerofoil.

Figure 26. Relative percentage of the skin friction drag coefficient and form drag coefficient with he increase of angle-of-attack for base aerofoil.

Figure 27. Skin friction drag coefficient and form drag coefficient with respect to the total drag coefficient for iced aerofoil.

Figure 28. Relative percentage of the skin friction drag coefficient and form drag coefficient with the increase of angle-of-attack for iced aerofoil.

But the scenario was different for the iced aerofoil as can be seen from Figs. 27 and 28. As the surface roughness due to ice accretion gives rise to an early transition, the skin friction drag of the iced aerofoil (0.004287) was more compared to the base aerofoil (0.00410) but the total drag was dominated by the from drag (80.34% of the total drag) at zero deg angle-of-attack and it was due to the irregular size and shape of the iced aerofoil.

As the angle-of-attack increases the percentage of the skin friction drag to the total drag falls while the form drag increases rapidly.

2. Case II

The analysis of the results reports an average of 386.95% (Fig. 30) drag coefficient for the iced aerofoil compared to the clean one for the whole range of angle-of-attack. The analysis also reported that the lift coefficient falls to 89.713% (Fig. 29) of the clean aerofoil while the drag coefficient increased to 844% only at 2 deg angle-of-attack. However, the lift coefficient tends to increase again after 2 deg and showed the lowest possible degradation (reduction in lift coefficient is 77.42% and increase in drag coefficient was 199.77%.) at 10 deg angle-of-attack and finally the aircraft stalls at beyond 16 deg. Because due to flow separation at 2 deg angle-of-attack and reattachments beyond 2 deg.

Figure 29. Coefficient of lift as a function of angle-of-attack for numerically obtained ice shape compared with the base aerofoil.

Figure 30. Coefficient of drag as a function of angle-of-attack for numerically obtained ice shape compared with the base aerofoil.

Form drag versus skin friction drag

Figure 31. Skin friction coefficient and form drag coefficient for with respect to the total drag coefficient for base aerofoil.

Figure 32. Relative percentage of the skin friction drag coefficient and form drag coefficient with the increase of angle-of-attack for base aerofoil.

Figure 33. Skin friction coefficient and form drag coefficient for with respect to the total drag coefficient for iced aerofoil.

Figure 34. Relative percentage of the skin friction drag coefficient and form drag coefficient with the increase of angle-of-attack for iced aerofoil.

For case II (velocity 142m/s) the skin friction drag for both base aerofoil (at 4 deg AoA0.004508) and iced aerofoil (0.004514 at 4 deg AoA) represented by Figs. 31, 32, 33 and 34 respectively, was more than case I (0.0044 and 0.004506 respectively) because a higher angle-of-attack gives rise to an early transition increasing the skin friction drag.

Figure 35. Coefficient of lift as a function of angle-of-attack for numerically obtained ice shape compared with the base aerofoil.

Figure 36. Coefficient of drag as a function of angle-of-attack for numerically obtained ice shape compared with the base aerofoil.

Figure 37. Velocity vector showing flow separation for iced aerofoil at 4° AoA.

Figure 38. Skin friction coefficient and form drag coefficient for with respect to the total drag coefficient for base aerofoil.

Figure 39. Relative percentage of the skin friction drag coefficient and form drag coefficient with the increase of angle-of-attack for base aerofoil.

Figure 40. Skin friction coefficient and form drag coefficient for with respect to the total drag coefficient for iced aerofoil.

Figure 41. Relative percentage of the skin friction drag coefficient and form drag coefficient with the increase of angle-of-attack for iced aerofoil.

Figure 42. Comparison of skin friction drag as a function of angle-of-attack for three different icing intensities with different velocities.

Figure 43. Comparison of form drag as a function of angle-of-attack for three different icing intensities with different velocities.

For the iced aerofoil in this case with worst icing condition as described in section the drag was mostly form drag (93%–96% of the total drag) from the very beginning because of the large leading edge ice accumulation though the skin friction drag increases slightly with the angle-of-attack. The result of this study shows that the shape, size and orientation of the body has more effect on the total drag generated on an aerofoil.

3. Case III

From Fig. 35, the behaviour of iced aerofoil shows that the coefficient of lift falls only at 4° AoA as the region of high velocity separates from the aerofoil section, giving rise to the stall condition. However, for the iced aerofoil it increased again with angle-of-attack due to flow reattachment (Fig. 37) and again falls beyond 16° AoA. But the aerodynamic degradation of the iced aerofoil compared to the base aerofoil continued to be reducing with the increase in angle-of-attack as seen from Fig. 36.

Form drag versus skin friction drag

As can be seen from the Figs. 38 and 39, initially the drag was dominated by skin friction for the base aerofoil but as the angle-of-attack was increased pressure drag became prominent due to the orientation of the body. In terms of total percentage of drag skin friction drag yields only 6% of the total drag at maximum angle-of-attack and rest of the drag was generated due to the form drag.

For iced aerofoil of case III with a maximum velocity of 167m/s compared to case I and case II, the skin friction drag was the maximum at all angles of attack as seen from Fig. 40. At an angle-of-attack of 4 deg, skin friction drag was 0.004722, 0.004789 and 0.00493, respectively, for case I, II and III. As a total percentage of drag shown in Fig. 41 skin friction drag contributes to only 4.69% of the total drag at maximum angle-of-attack (18 deg) when the body becomes the less streamlined.

Figures 42 and 43 shows the comparison of skin friction drag and form drag for the three iced aerofoil. Results of the analysis show that with the increase of velocity the transition point moves forward giving rise to higher skin friction drag coefficient. Thus, ice case 3 with maximum velocity has the earliest transition and hence greater skin friction drag. Whereas case 2 is dominated by form drag mostly compared to case 1 and 3 because of highest ice accumulation and hence distorted aerofoil shape. Though the skin friction drag increases with angle-of-attack at higher angle-of-attack as the body becomes less streamlined drag is mostly dominated by the pressure drag for all three cases.

5.4 Applicable range and limitations of the proposed numerical model

1. a. By applying the present model a sensitivity analysis of roughness models can be made on the ice shape which is recommended for future extension of the present work.

2. b. As mentioned in various literatures, certain parameters affect the icing behaviour; this study considers free-stream velocity, freezing fraction and time of ice accretion as the prime parameters for varying ambient conditions. Effects of other parameters can be studied for further improvement of the model.

3. c. As it is not possible to determine impinging water mass directly, it is calculated indirectly using collection efficiency which is a major drawback of the current work.

4. d. Assuming that ice grows only in the perpendicular direction to the aerofoil surface, these ice growth heights are added normal to the surface panels to produce a new profile which is one of the limitations of the proposed method.

5. e. The present model calculates the ice thickness in multilayer mode that is unable to predict ice horn at different locations which can be obtained by one layer mode that overestimates the runback water and produce ice horns.

5.5 Limitations of the CFD model

1. a. A 3D simulations could have given better results than the 2D CFD study but at the cost of computational time and expenses.

2. b. RANS turbulence could have given better accuracy but it increases computational power and time.

3. c. The element size considered for the present CFD study was 8e⁻⁰⁰³m for meshing; a smaller element size would have given better accuracy.

4. d. Higher order of accuracy could be chosen for better CFD analysis.

5. e. As the k-epsilon model is based on average velocity concept, it can not produce turbulent wake.

6. f. The equation for in model, the k-epsilon model has the difficulty for solving epsilon, so it is postulated and that is why it is not the model with zero error.

7. g. Assessment and reduction of truncation error due to discretisation was neglected in the computational analysis.