ICE ACCRETION SIMULATIONS ON AIRFOILS

Ice shape predictions for a NACA0012 airfoil and collection e©ciency predictions for the Twin Otter airfoil are obtained and presented. The results are validated with reference numerical and experimental data. Ice accretion modeling mainly consists of four steps: §ow ¦eld solution; droplet trajectory calculations; thermodynamic analyses; and ice accretion simulation with the Extended Messinger Model. The models are implemented in a FORTRAN code to perform icing analyses for twodimensional (2D) geometries. The results are in good agreement with experimental and numerical reference data. It is deduced that increasing computational layers in calculations improves the ice shape predictions. The results indicate that collection e©ciencies and impingement zone increase with increasing droplet diameter.


INTRODUCTION
Ice accretion on airframes during §ight may cause great danger due to aerodynamic performance degradation and engine power loss.Hence, it is very crucial to simulate ice accretion in order to predict ice mass that accumulates on the surface and the regions which are prone to icing.Such a simulation will be useful to design and develop a de/anti-icing system for aircrafts and for airworthiness certi¦cation purposes.
In this study, ice shape predictions for a NACA0012 airfoil and collection e©ciency distributions for the Twin Otter airfoil are obtained.In the current approach, panel method is used for the §ow ¦eld solution.Droplet trajectories, collection e©ciency, and convective heat transfer coe©cient calculations are performed with the same approach as in the study of Ozgen and CanŠbek [1].The study of Myers [2] is used as a base for the ice accretion model.
The pioneering work of Messinger [3] represents an important foundation and a milestone in numerical ice accretion simulation, which introduces the Original Messinger Model that is a one-dimensional (1D), equilibrium energy balance.
With the increase in performance of digital computers in the 1970s, considerable amount of e¨ort was directed towards studying realistic geometries such as airfoils, wings, helicopter blades, and engine nacelles.The major contribution to the knowledge and the database was provided by NASA Lewis Research Center (USA), Defence Evaluation and Research Agency (DERA-UK), O©ce National d£Etudes et des Recherches A‚ erospatiales (ONERA ¡ France), Anti-Icing Materials International Laboratory (AMIL ¡ Canada), Ecole Polytechnique de Montreal (Canada), and Italian Aerospace Research Center (CIRA ¡ Italy).
The results of a collaborative e¨ort on icing research are published by Wright, Gent, and Gu¨ond [4].The report contains data obtained for airfoil geometries by three organizations including heat transfer coe©cients, collection e©ciency distributions, as well as ice shapes.The e¨ect of multistep computations is also studied and reported.
In their review paper, Gent, Dart, and Cansdale present the background and the status of the analyses developed to address the problem of icing on aircraft [5].The methods for the calculation of droplet trajectories and heat transfer coef-¦cients, ice accretion prediction, and aerodynamic performance degradation are outlined and discussed together with representative results.
Myers [2] presents a 1D mathematical model, addressing the shortcomings of the Original Messinger Model.Instead of setting the value of temperature to its equilibrium value, temperature distributions both in the ice and water layers are taken into account, bringing out the e¨ect of conduction in the ice/water interface.Including conduction has a cooling e¨ect at the ice/water interface, hence increasing ice production.Another shortcoming addressed by the Extended Messinger Model is that the freezing fraction decreases monotonically with increasing exposure instead of the ¦xed values for rime and glaze ice conditions producing a discontinuity of this parameter during transition in the Original Messinger Model.
Potapczuk and Bidwell [6] present a method for three-dimensional (3D) ice accretion modeling.Three-dimensional §ow ¦eld methods and droplet trajectory calculations are combined with 2D ice accretion calculations.
The e¨ect of supercooled large droplets (SLD) on icing has been reported by Wright and Potapczuk [7], where the methods used in the ice accretion software have been reviewed to account for SLD e¨ects.The numerical results obtained have been compared with the experimental data.The e¨ects of SLD on aircraft icing have attracted a lot of interest from researchers in the recent years and just very recently, icing conditions related to SLD have been included in the Certi-¦cation Speci¦cations for Large Airplanes by Federal Aviation Administration (FAA) in FAR-25 (FAA Regulations Part 25) [8] and European Aviation Safety Agensy (EASA) in CS-25 [9].
In a similar context, Papadakis et al. [10] have reported a comprehensive account of droplet impingement for small and large droplets on four di¨erent airfoils and two simulated ice shapes.Wind tunnel experiments have been supported by computations using the LEWICE 2D computer code.In general, good agreement was observed between experimental and numerical results for small droplets.For large droplets, however, the numerical results exhibited higher collection e©ciency levels and wider impingement limits.This discrepancy is attributed to the mass loss due to splashing experienced by large droplets during impingement.
Another phenomenon that has attracted a lot of interest in the recent years is icing in glaciated and mixed-phase conditions.It is reported that icing, particularly on aircraft engine components, may result in rollbacks, mechanical failures, and §ameouts.Villedieu, Trontin, and Chauvin [11] present numerical methods and results in glaciated and mixed-phase icing conditions using the ONERA 2D icing suite.Of particular interest are the models related to trajectory computations, heat and mass transfer during trajectory calculations, and impingement and accretion models.Again just very recently, meteorological conditions related to glaciated and mixed-phase conditions have been included in the Certi-¦cation Speci¦cations for Large Airplanes both by FAA in FAR-25 and EASA in CS-25.This manuscript summarizes the methods which are used in the developed computational tool and the results obtained for Appendix C and SLD icing conditions.The current computational tool has been in continuous development since 2007, ¦rst developed for Appendix C droplet conditions, then extended for SLD and most recently for glaciated and mixed phase conditions.The tool has also been extended for 3D wing and intake geometries, where parallel computation has been implemented for droplet trajectory calculations [12,13].Although most of the icing softwares that are in existence utilize the Original Messinger Model, the current computational tool employs the Extended Messinger Model.In the current study, the Extended Messinger Model has been extended to two and three dimensions, together with an extension to handle runback water.The results in the current manuscript include ice shape predictions for varying droplet median volume diameters (MVD), liquid water content (LWC), temperature, velocity, and exposure time.The e¨ect of SLD on ice shape calculation is also investigated.The collection e©ciency analyses are performed for di¨erent MVD and angles of attack (α).The comparisons with numerical and experimental data available in the literature are also presented.

METHODOLOGY
In this section, four modules used in the current approach are explained in detail.

Flow Field Solution
The 2D HessSmith panel method is used to determine the §ow velocities required for droplet trajectory calculations.The external velocity distribution is also provided, which is used in the boundary layer calculations to obtain heat transfer coe©cients.In this model, the airfoil is divided into quadrilateral panels, each having a constant strength source singularity element plus a vortex singularity that is constant for all panels.Using the §ow tangency boundary condition at the collocation points of the surface panels, the singularity strengths are obtained.This provides the velocity potential, in other words, the §ow velocity components at any point in the §ow ¦eld.The HessSmith panel method has also been extended to handle axisymmetric air intakes, where the mass §ow rate through the inlet also needs to be maintained in addition to §ow tangency condition at the control points.This requires a two-step computational procedure [14,15].

Droplet Trajectories and Collection E©ciency Calculations
The following assumptions are used for droplet trajectories: droplets are assumed to be spherical except for SLD where nonsphericity is also accounted for; the §ow ¦eld is not a¨ected by the presence of the droplets; and gravity and aerodynamic drag are the only forces acting on the droplets.
The governing equations for droplet trajectories are: (1) m y p = −D sin γ + mg ; (2) In the above equations, V x and V y are the §ow velocity components at the droplet location, while ' x p , ' y p , x p , and y p are the components of the droplet velocity and acceleration, respectively.Atmospheric density is denoted by ρ, while droplet cross-sectional area and drag coe©cient are represented by A p and C D , respectively.The droplet drag coe©cient is computed as a function of the droplet Reynolds number, Re = ρV rel d p /µ based on the droplet diameter d p , relative droplet velocity V rel and the atmospheric viscosity µ.
Drag coe©cients for droplets are calculated using empirical drag laws based on the droplet Reynolds number [5] for Appendix C droplets and droplet Reynolds number and eccentricity for SLD [16].Di¨erent drag laws are used for Appendix C droplets and SLD, where nonsphericity is accounted for in the latter.For the computation of the droplet trajectories, the Lagrangian approach is used where the droplet trajectories are obtained by integrating Eqs. ( 1) and ( 2) in time until the droplets impact the geometry.The particle impact pattern on the section determines the amount of water that impinges on the surface and the region subject to icing.The local collection e©ciency (β) is de¦ned as the ratio of the area of impingement to the area through which water passes at some distance upstream of the section.

Thermodynamic Analysis
In order to calculate the convective heat transfer coe©cients, a 2D Integral Boundary Layer Method is employed for both laminar and turbulent §ows [1].

Extended Messinger Model
Ice accretion on the geometry is found with the Extended Messinger Method.The ice shape prediction is based on phase change or the Stefan problem.The governing equations for the phase change problem are: energy equations in the ice and water layers, mass conservation equation, and a phase change condition at the ice/water interface [1]: (3) (4) In Eqs.(3)(6), θ and T are the temperatures; k w and k i are the thermal conductivities; C pw and C pi are the speci¦c heats; and h and B are the thicknesses of water and ice layers, respectively.On the other hand, ρ i and L F denote the density of ice and the latent heat of solidi¦cation of water, respectively.Ice density is assumed to have di¨erent values for rime ice, ρ r , and glaze ice, ρ g .The coordinate y is normal to the surface and ρ a is the LWC.In Eq. ( 5), ' m in and ' m e,s are impinging, runback, and evaporating (or sublimating) water mass §ow rates for a control volume, respectively.The boundary and initial conditions accompanying Eqs.(3)(6) are [2]: ice is in perfect contact with the wing surface: The surface temperature is taken to be the recovery temperature [5]: In the above expression, M = V ∞ /a ∞ , while the speed of sound is given by a ∞ = √ γRT a .Additionally, r is the adiabatic recovery factor: r = Pr 1/2 for laminar and r = Pr 1/3 for turbulent §ows, with Pr being the Prandtl number of air; the temperature is continuous at the ice/water boundary and is equal to the freezing temperature, T f : at the air/water (glaze ice) or air/ice (rime ice) interface, heat §ux is determined by convection, radiation, latent heat release, cooling by incoming droplets, heat brought in by runback water, evaporation, or sublimation, aerodynamic heating, and kinetic energy of incoming droplets; and wing surface is initially clean: In the current approach, each panel constituting the geometry is also a control volume.The above equations are written for each panel and ice is assumed to grow perpendicularly to a panel.
Rime ice growth is expressed with an algebraic equation from the mass balance in Eq. ( 5), since water droplets freeze entirely and immediately on impact: On the other hand, glaze ice thickness is obtained by integrating the ordinary di¨erential equation got by combining mass and energy equations over time.The di¨erential equation is: In this expression, Q c is the heat §ux by convection; Q e is the evaporation; Q d is the heat from incoming droplets; Q r is the radiation; Q a is the aerodynamic heating; Q k is the kinetic energy of incoming droplets; and Q in is the energy entering the control volume due to runback water.It is assumed that all of the unfrozen water passes to the neighboring downstream cell as runback water at the upper surface, while all water sheds at the lower surface [8].To calculate the glaze ice thickness, Eq. ( 7) is integrated numerically using a RungeKutta Fehlberg method.

Ice Shape Results
Ice shape predictions are obtained for the NACA0012 airfoil geometry.The results are compared with the numerical data reported in the literature obtained with LEWICE 2.0 and LEWICE 3.0 (software developed by NASA) and experimental data which are presented by Wright and Potapczuk [7].Chord length of the airfoil is 0.53 m and the angle of attack is 0 • for all the test cases; MVD, LWC, velocity, total temperature, and exposure time are varying for the test cases which are shown in Table 1.
In the calculations, multilayer calculation approach and the e¨ect of SLD are investigated.In multilayer calculation approach, exposure time is divided into segments.At the beginning of each time interval, the iced surface is considered as the new geometry to be exposed to icing and all the calculations are repeated.The e¨ect of SLD can be understood better by explaining the behavior of large droplets.Droplet breakup and splash are the phenomena which are characteristics of large droplets.The droplets with MVD greater than about 100 µm are considered as SLD.Inclusion of SLD e¨ects (breakup and splash) for the droplets with MVD lower than this value does not give correct results in the calculations.In the results, the cases for which the SLD e¨ects are included are shown as ¤SLD: on.¥In the same manner, the calculations performed without SLD e¨ects are stated as ¤SLD: o¨.¥  Figure 1 shows the parametric study of multilayer calculations with and without the inclusion of SLD e¨ects, respectively.The ice shapes obtained are compared with the experimental data.Although increasing number of layers does not change ice shape much, the result with 12 layers is slightly closer to the experimental data.The extent of ice and the ice shape, in general, is well-predicted, while there is a slight underestimation of the maximum ice thickness close to the leading edge.The ice shapes for this case are largely rime ice shapes due to low ambient temperature with relatively smooth contours but the presence of hornlike shape suggests that there is also glaze ice present, due to high LWC and high speed, in spite of the low ambient temperature.In Fig. 2, the e¨ect of SLD on ice shape prediction for 12 layers of calculation is shown.The prediction without SLD e¨ect can be said to be better when compared with the experimental data in terms of the symmetrical horn shape.This is expected since MVD for Test 1-22 is 40 µm which is not considered as an SLD case and inclusion of breakup and splash does not improve the ice shape prediction.
The best ice shape prediction obtained is without SLD e¨ects and 12 layers of calculation as shown above.This result is compared with numerical and experimental literature data in Fig. 3.Although LEWICE 2.0 predicts the horn  shape well, it overestimates the ice thicknesses in the horn regions.LEWICE 3.0 predicts a smaller ice mass and smoother ice shape.The current study captures the horn shape better than others, although it slightly underpredicts the ice thickness compared to the experimental result.
In Fig. 4, ice shape predictions for Test 1-1 are presented.As seen in Figs.4a  and 4b, increasing the number of calculation layers does not improve the results signi¦cantly, especially for the case where SLD e¨ects are included.In Test 1-1, MVD = 70 µm, which is not really an SLD case.This is backed up with Fig. 5, which shows that a better ice shape is obtained by excluding droplet breakup and splash e¨ects.The ice shapes for this case suggest that the conditions are Figure 6 shows the comparison of the results of the current study (12 layers of calculation, without SLD e¨ects) with reference data.All the results can be said to be similar to each other when compared with the ice shape obtained in the experiment.The impingement zone is predicted slightly wider in both current and reference results.However, the current study is slightly better than the others in terms of ice thickness and the overall ice shape.
The droplet diameter is increased to 160 µm for Test 1-4.In Fig. 7, it is clearly seen that increasing number of layers in the calculations does not improve the ice shape prediction for both with and without SLD e¨ects.This is mainly due to the fact that the ice shapes are almost pure-rime, where SLD e¨ects play very little role as the droplets freeze immediately upon impact with the geometry.
Moreover, Fig. 8 shows that including SLD e¨ects gives slightly closer ice shape prediction to experiment.
For 6 layers of calculation and the SLD: on case, the current study result is shown with the numerical and experimental data in Fig. 9. LEWICE 2.0 obtains a thicker ice mass when compared with the others in the horn regions.The current study and LEWICE 3.0 results are quite similar to each other, which predict the thickness successfully but miss the horn shapes.

Collection E©ciency Results
In this part of the study, collection e©ciencies are calculated for the Twin Otter airfoil pro¦le for varying droplet diameters and angles of attack.The validation study is performed for the §ow over a pro¦le with a chord length of 1.448 m and velocity of 78.25 m/s.Median droplet diameters vary between 11 to 168 µm.The analyses are done for α = 0 • and 4 • .Collection e©ciency results are validated with numerical and experimental data presented by Papadakis et al. [10].Liquid water content sprayed out during the experiment varies with the droplet diameter which is shown in Table 2.The droplets in the last two rows of the table are representative of SLD.  Figure 11 shows the collection e©ciency distributions for droplets with MVD of 11, 21, 79, 137, and 168 µm at α = 0 • .The current study results are compared with the numerical results obtained with LEWICE developed by NASA and also experimental data.Since the angle of attack is 0 • , maximum collection e©ciency value is at the leading edge where surface distance from highlight is 0 mm.The case corresponding to MVD = 11 µm (Fig. 11a) probably corresponds to the lower limit for the applicability of the drag coe©cient formulations because the droplet Reynolds number is very small, which is probably a contributing factor to the discrepancy of the current results and the literature data.For the case MVD = 21 µm (Fig. 11b), almost the same results are obtained from numerical and experimental literature data and the current study.Although Figure 10 De¦nition of surface distance from highlight in millimeters (signs) [10].the impingement zone is slightly underpredicted, the maximum β value is estimated well.When the other cases are investigated, it is seen that higher values of β are observed in the current study than the reference data, especially in the regions away from the leading edge.The same inference is also valid for LEWICE results, to a higher extent especially for MVD = 137 and 168 µm cases (see Figs. 11d and 11e).From Fig. 11, it can be deduced that collection e©ciency distribution and maximum values are captured fairly successfully for α = 0 • .In Fig. 12, collection e©ciency distributions for varying droplet diameters is seen.As expected, an increase in droplet diameter results in an increase in β as well.Moreover, impinge- ment limits both on lower and upper surfaces are wider for larger droplets.The reason for this is that larger droplets have higher inertia and, therefore, follow ballistic trajectories.This results in more particles impinging the surface and the impingement angle being higher, resulting in higher collection e©ciencies and wider impingement zones.When the angle of attack is increased to 4 • , collection e©ciency distributions are obtained as shown in Fig. 13.Due to angle of attack, maximum values of collection e©ciency shift to right, in other words, to the positions on the lower surface.

Curve shows the airfoil geometry
For α = 4 • cases, it is noticed that the maximum collection e©ciency values are almost the same as for α = 0 • cases.Current study results where MVD = 21 µm (see Fig. 13b) can be said to be the closest ones to the experimental results.Similar results to 0 • cases are obtained for the other droplet diameter cases.The location of the maximum β value is predicted almost the same, although the maximum value itself is found di¨erent.As in 0 • cases, impingement limits on lower and upper surfaces are calculated di¨erent as well.The current results indicate a slight overestimation of the collection e©ciencies away from the leading edge especially for MVD = 137 and 168 µm cases (see Figs. 13d  and 13e).
Figure 14 shows the collection e©ciencies for varying MVD values for angle of attack of 4 • .Similar to 0 • case, collection e©ciency values are higher and impingement zone is wider for larger MVD of droplets.
In this validation study, the e¨ects of angle of attack and droplet diameter on collection e©ciency are investigated.Maximum β values occur at the leading angle of attack is increased, ice accumulates more on the lower surface and maximum value position shifts to here.For a given angle of attack, increasing droplet diameter not only increases the maximum β value, but also extends the impingement zone both on the lower and upper surfaces.
The results suggest that there is a room for improvement regarding the droplet splash model.The state-of-the-art model was developed using experimental data obtained using a test setup and test article not exactly representing the current application.A new set of experimental data is de¦nitely needed.
It is also known that droplet breakup constitutes a secondary in §uence, except for very large droplets and high speeds [17].

CONCLUDING REMARKS
In the current study, ice shape predictions for a NACA0012 airfoil geometry and collection e©ciency calculations for the Twin Otter airfoil geometry are performed.
In the ¦rst validation study, it is seen that increasing number of layers of calculation improves ice shape prediction and a closer ice shape to experimental result is obtained.Moreover, for the droplets with median droplet diameter smaller than about 100 µm, including SLD e¨ects like breakup and splash, does not enhance the results unlike for the larger droplets.
In the collection e©ciency analysis, the e¨ects of angle of attack and droplet diameter are studied.It is deduced that the maximum collection e©ciency value occurs at the leading edge for zero angle of attack.The increase in angle of attack results in the shift of the position, where the maximum β value is observed, to the lower surface for positive angle of attack, although the value of this parameter itself seems to be fairly independent of angle of attack.Furthermore, for a given angle of attack, collection e©ciency becomes higher and the impingement zone widens when droplet diameter is increased.
The results for both ice shape prediction and collection e©ciency show that the current study results are in good agreement with the reference numerical and experimental data but there is a room for improvement in droplet splash modeling.It can be concluded that the current tool has a potential to be used for certi¦cation purposes as well as for the design of de/anti-icing equipment on aircraft.In order to assess the e¨ect of ice on aerodynamic performance, clean and iced geometries can be provided as inputs to a §ow solver employing Reynolds-averaged NavierStokes equations, so that the impact of ice on lift, drag, and moment characteristics can be brought out.

Figure 8 Figure 9
Figure8The SLD e¨ect on ice shape prediction for Test 1-4: 1 ¡ SLD: o¨; 2 ¡ SLD: on; 3 ¡ experimental data; and 4 ¡ clean airfoil Fig. 10 [10].It is clearly seen that highlight position is the leading edge of the airfoil.Negative values show the upper surface while positive values represent the lower surface.

Table 1
Test cases for ice shape predictions

Table 2
[10]LWC variation with MVD[10]Collection e©ciency distributions for di¨erent droplet diameters are shown in Figs.1014.The horizontal axis shows the surface distance from highlight (mm) which is de¦ned as in