In-flight icing on unmanned aerial vehicle and its aerodynamic penalties

A numerical prediction of ice accretion on HQ309, and airfoils and its aerodynamic penalties is described. Ice accretion prediction on a three-dimensional (3D) swept wing is also presented. In ad-dition to air§ow and drop trajectory solvers, NRC£s (National Research Council) original, 3D, morphogenetic icing modeling approach has been used. The analysis was performed for a wide range of icing conditions identi¦ed in the FAA (Federal Aviation Administration) Appendix C icing envelope. They cover a range of drop sizes, air temperatures, and liquid water contents. For selected icing conditions, the resulting decrease in lift and increase in drag have been calculated.


INTRODUCTION
Flying through clouds containing supercooled liquid water can be hazardous, because of the risk of ice accretion on the leading edges of wings.The use of unmanned aerial vehicles (UAVs) to perform surveillance and reconnaissance operations has increased signi¦cantly in recent years.Consequently, the need for UAVs to operate in adverse weather conditions becomes increasingly crucial.Icing reduces lift and stall angle and increases aerodynamic drag, aircraft weight, and power requirements.While icing impacts the operation of all types of aircraft, icing on UAVs presents a unique set of challenges.The very features that make UAVs attractive to the operator ¡ high aerodynamic e©ciency, long endurance, small size, low power and remote control ¡ combine to make ice protection a challenge.In addition, UAVs are particularly vulnerable to icing, as a result of generally lower cruising altitudes and, consequently, higher moisture content.Their comparatively low air speed leads to longer exposure time and increases their vulnerability.Therefore, a better understanding of ice formation and its aerodynamic consequences on UAVs is needed.
The objective of this study is to use numerical simulation to investigate ice formation on a number of UAV airfoils and to evaluate its consequences for their aerodynamic performance.Ice accretion simulations were also performed on a 3D swept wing.It is assumed that the UAV airfoils §y through atmospheric icing conditions identi¦ed in FAA Federal Aviation Regulations (FAR) 25 Appendix C. Presently, there are very limited, publicly available research results dealing with the prediction and consequences of icing on UAVs [14].

NUMERICAL MODULE DESCRIPTION
There are three main physical processes that in §uence ice formation, each with a corresponding numerical module.Air §ow is determined by the NRC, in-house, unsteady computational §uid dynamics (CFD) §ow solver INS §ow, or by the German Aerospace Center (DLR) compressible FLOWer code.Drop trajectories are determined using a custom drop trajectory solver.Drop motion and freezing on the surface is calculated using NRC£s proprietary ¤morphogenetic¥ approach, which is based on a discrete formulation and simulation of ice formation physics.Below, the main characteristics of these three predictive tools are brie §y described.

Computational Fluid Dynamics Solvers
The in-house, block-structured, incompressible code INS §ow [5,6] was used for the HQ309, SD7032, and SD7037 airfoil §ow calculations.The code is applicable to small UAVs when the Mach number is less than 0.3.It was validated for various §ows in the low-Reynolds-number regime.The integral forms of the mass and momentum conservation laws are used in INS §ow.A fully implicit, secondorder, temporal di¨erencing scheme is used in the discretization, which makes the algorithm stable for large time steps.The discretization of the convective and di¨usive §uxes is carried out in a co-located variable arrangement, using a ¦nite-volume approach that is second-order accurate in space.The heat transfer is computed from the §ow temperature distribution, by solving the energy equation based on the §ow solution [7].The two-dimensional (2D) simulations were performed on C-type meshes with 801 × 97 grid points.The far-¦eld of the computational domain was set at about 25 chords from the airfoil.The computations were performed for angles of attack from 0 • to 9 • at 3 degree intervals and from 10 • to 18 • at 1 degree intervals.The LangtryMenter correlation-based γ-Re θt transition model [8] was used re §ecting the nature of low-Reynoldsnumber aerodynamics.The free-stream turbulence intensity was set as Tu = 0.001 and to re §ect ice surface coarseness for heat transfer computations, the roughness of the surface was assumed to be 2 mm.
Flow calculations for a 3D swept wing were performed using the DLR compressible FLOWer code [9].The Illinois Model (IM) swept wing has been considered.It is a version of the Common Research Model (CRM) and is fully described in [10].The FLOWer code was chosen because of potential extension of our capabilities to other larger wings or faster §ight speeds.Although the Reynolds number was still located in the transitional regime, it is believed that the transition process is very limited in this 3D swept wing case, in particular, when it is iced.As a result, the §ow was assumed to be fully turbulent and Menter£s SST (shear stress transport) turbulence model was applied, without transition model coupling.The NavierStokes solver FLOWer uses a block-structured computational domain around the aerodynamic con¦guration.The §ow equations were discretized based on the ¦nite-volume approach.A cell-centered formulation was used for the discretization.Convective §uxes were evaluated using a second-order accurate central di¨erencing scheme with scalar dissipation.The ¦rst-order accurate §ux-di¨erence upwind scheme was used for the convective §uxes of the turbulence equations at the ¦nest grid level.This scheme was used mainly because of issues related to the numerical stability and convergence.
The second-order scheme may have possibly improved the solution at some §ow conditions but, in addition, it was not available in used version of the CFD code.

Drop Trajectory Solver
The drop trajectory solver uses the §ow-¦eld velocity from the CFD §ow solver.The drop trajectories are computed using a Lagrangian approach.It is assumed that the water drops are monodisperse and spherical and that they do not disturb the air §ow.The only forces considered to act on the drops are aerodynamic drag and gravity.Initially, the drops are placed far upstream from the airfoil, on a line perpendicular to the free-stream direction.The initial drop velocity is set equal to the computed air §ow velocity at that point plus the drop terminal velocity.The §ow solver calculates the values of the air velocity on a grid, and these are passed to the drop trajectory solver.Because the drop trajectory code requires air §ow velocities at o¨-grid locations, the air velocity components at the drop locations are calculated as a weighted average of the values at the closest grid locations.The weightings are proportional to the inverse distance between the drop and the grid point.Computation of the drop trajectories determines the distribution of collision e©ciency and drop impingement direction along the airfoil, which is passed to the morphogenetic ice accretion model.

NRC Morphogenetic Ice Accretion Solver
The morphogenetic model is a discrete element, random walk model that emulates the motion and freezing of individual §uid elements arriving at the accretion surface.The model §uid elements may be imagined to consist of an ensemble of cloud drops, all of which undergo identical histories.A 3D rectangular lattice de¦nes the accretion domain.Each §uid element begins a stochastic motion downstream along the surface from its initial impact location.At each step in the process, a random number is generated and, according to its value, the element either freezes or moves along the surface.A number of simplifying assumptions have been made in the development of the present version of the morphogenetic model.One of these is that the §ow and drop trajectory solutions are obtained for a clean wing.This approach is justi¦able for small ice growth but for larger ice dimensions, it leads to inaccuracies.However, it is a goal of this paper to show that even a one-step approach to drop impingement gives at least qualitatively realistic results.The present authors are planning further model improvement by using a multistep approach to drop impingement.Because the morphogenetic approach to ice accretion growth is ¦ne-grained and multistep, it may partly compensate for the single-step approach to droplet impingement.
The present icing model is based on conservation of mass and energy for impinging, §owing, and freezing water on the airfoil surface.It is assumed that once a drop hits the airfoil surface without splashing, part of it may freeze in situ and part of it will §ow downstream along the surface under wind stress where further freezing may occur.The variation of the in¦nitesimal water mass §ux along the airfoil surface is, therefore, determined by the di¨erence between the impinging water mass §ux and the freezing mass §ux.It is assumed that the physical processes occurring on the ice surface are governed by a steady-state heat balance, which includes the convective, evaporative, aerodynamic heating and sensible heat §uxes [3].
The model is sequential, so that as soon as a particular §uid element freezes, the behavior of the next §uid element is considered.It is assumed that the mass of each §uid element is the mass of a frozen ice cube of dimensions 200 µm.The shape of the moving §uid element is not relevant in this model.However, its mass is the same as the mass of a spherical water droplet of diameter 241 µm.The total number of impinging §uid elements is determined by the free-stream velocity, liquid water content, §uid element size, duration of the icing event, and spatial distribution of the collision e©ciency.The results presented in the next sections were obtained in a morphogenetic model domain of 500 × 500 × 500 elements, equivalent to 10 × 10 × 10 cm in each direction.The §uid element size was chosen as a compromise between model precision and computational e¨ort.The full details of the morphogenetic ice accretion model were described in [11,12].

MODEL RESULTS AND DISCUSSION
Ice accretion and its aerodynamic penalties have been examined for a number of 2D and 3D geometries and a range of icing conditions.The authors started with an HQ309 airfoil and employed 2D CFD analysis to compute the §ow and aerodynamic penalties.However, a 3D version of the morphogenetic model was used to predict the ice formation.Next, ice accretion and its aerodynamic consequences were compared for three airfoils exposed to the same §ight and icing conditions.Finally, a 3D predictive analysis of the air §ow and ice accretion on a swept wing was performed.

Prediction of Ice Accretion and Its Aerodynamic Consequences on an HQ309 Airfoil
An HQ309 airfoil, which has been considered for UAVs, is used in the following analysis.The cruise speed was assumed to be 25.7 m/s with 3 degree angle of attack and an ambient air pressure corresponding to 1500 m above sea level.It was also assumed that the airfoil, with a chord of 0.47 m, travels a total distance of 32.2 km, which is given in FAR 25 Appendix C as a base value.This leads to §ight durations of 20.9 min.The choice of values for the icing parameters was also based on the Appendix C icing envelope.The aim was to cover the entire range of possible conditions.First, the cases with drop size (median volume diameter, MVD) 25 µm and static air temperature −2, −10, and −30 • C, associated with liquid water contents 0.46, 0.30 and 0.10 g/m 3 , respectively (as speci¦ed in the FAR 25 Appendix C) were considered.Figure 1 shows the predicted 2-centimeter wide ice accretion sections as seen from below the airfoil.Since the drop size was the same and trajectories were computed only for a clean airfoil and were not recalculated during ice growth, the distribution of the collision e©ciency is identical for all three cases.The di¨erences between the ice shapes results from the di¨erent air temperature and liquid water content.At an air temperature of −2 • C, there is insu©cient heat transfer to freeze all the impinging drops instantly.Consequently, they §ow along the airfoil surface before freezing.This leads to the formation of ice rivulets that are clearly visible in Fig. 1a.At an air temperature of −10 • C, the impinging drops freeze instantly, forming the rime accretion shown in Fig. 1b.The rough rime ice features arise because of the drop shadowing e¨ect of locally forming and protruding ice.These ice lobes grow preferentially by intercepting incoming drops, while preventing drops from impinging directly behind them.Since a further decrease of air temperature to −30 • C is associated with a lower liquid water content, the resulting mass of rime ice accretion for this case is decreased, but the main ice features remain much the same.
To perform a 2D computation of the aerodynamic penalties due to ice accretion, the predicted ice shape was averaged in the spanwise direction.Figure 2 depicts the velocity distribution around the leading edge of an HQ309 airfoil iced at −2 • C. In the grid generation for the CFD simulation, the ice surface was smoothed using a smoothing spline technique, allowing the predicted ice surface to increase or decrease up to 10% of the resolution (20 µm in this case) of the ice predicted values.Although the iced surface was smoothed, a pronounced horn structure and small-scale surface irregularities are nevertheless visible.In the concave regions, §ow separation and vortical §ow structures with a large vortical velocity component are observed as shown by the blue color in Fig. 2.This §ow macrostructure, combined with an irregular surface roughness, leads to a reduction in aerodynamic performance; lift decreases by 10% and drag increases by 143% for the −2 • C conditions.The changes in the aerodynamic performance of the iced HQ309 airfoil are depicted in Fig. 3.The lift and drag coe©cients are based on purely aerodynamic forces.Consequently, the ice weight is irrelevant in this analysis.However, the magnitude of the performance degradation is a function of ice shape, location, and extent.Since the accretions at −10 and −30 • C are similar, the main di¨erence being ice thickness, their aerodynamic characteristics also appear to be similar.When the angle of attack is less than 3 • , the lift coe©cient for the −10  and −30 • C cases is marginally greater than for the clean airfoil.For the −2 • C case, the lift coe©cient is smaller than for the clean airfoil as a result of the horn-like structure.For the same range of angle of attack, the drag coe©cient increases for all the iced airfoil cases.For larger angles of attack, the smallest aerodynamic penalty in the lift and drag coe©cients seems to be for the −2 • C ice accretion.This is because the drooped leading edge increases airfoil camber and thus delays §ow separation when compared with −10 and −30 • C cases.However, all §ows in the iced airfoil cases stall at lower angles of attack than for the clean airfoil.The calculated reduction in maximum lift coe©cient could lead to safety issues for an UAV.The increase in the drag coe©cient causes higher energy consumption, thereby limiting vehicle endurance and range.It should also be noted that the accumulated ice adds to the UAV£s weight; this could be important for some UAV applications.
Also, the in §uence of drop size on ice shape was examined (Fig. 4).This was done over a range of air temperature, with corresponding values of the liquid water content as speci¦ed in the FAR 25 Appendix C envelope.Larger drops have greater inertia, which is re §ected in their slower response to changes in streamline direction and a smaller curvature of the drop trajectories.This results in a greater maximum collection e©ciency and a greater impingement extent.It should be noted that the trajectories were computed for a clean airfoil; they were not recalculated during ice growth.Smaller drops are associated with a higher liquid water content.However, the impingement extent and intercepted §ux diminish for smaller drops.The displayed ice shapes are the result of these two opposing tendencies.In addition, warmer temperatures lead to §ow of water on the surface before freezing.This also in §uences the accretion shape.For the considered cases, only an air temperature of −2 • C produces glaze ice.The resulting ice shape is strongly in §uenced by the distribution of the convective heat transfer.At −10 • C, rime ice forms and the ice shape and extent are closely related to the drop impingement characteristics.As a result, the ice extent increases with drop size, but di¨erences in the maximum ice thickness are largely determined by the di¨erences in liquid water content.Since the ice accretion at −30 • C is rime, the main di¨erence from the −10 • C shapes, which are also rime, is the smaller ice thickness, resulting from the lower liquid water content under colder conditions.

Comparison Between Ice Accretions Forming on HQ309, SD7032, and SD7037 Airfoils
There was examined the aerodynamic performance degradation resulting from ice formation on selected airfoils exposed to identical §ow and icing conditions.
The following three airfoils that could be used for UAV applications were chosen: HQ309, SD7032, and SD7037.The parameter values assumed in the previous section were also used in this analysis, namely, airfoil chord 0.47 m, angle of attack 3 • , ambient air pressure corresponding to 1500 m above sea level, cruise speed 25.7 m/s, total §ight distance 32.2 km, and §ight duration 20.9 min.In addition, the following values of icing parameters were assumed: drop size 25 µm, air temperature −10 • C, and liquid water content 0.30 g/m 3 .Lift and drag coe©cients of clean and ice covered airfoils are also shown.The predicted ice shapes on the three airfoils, under the chosen §ow and icing conditions, are rather similar (Fig. 5).However, smaller airfoil thickness is associated with smaller ice mass.The ice masses for HQ309, SD7037, and SD7032 are 108, 118, and 120 g/m, respectively.For the speci¦ed conditions, ice forms as rime, there is no water §owing along airfoil surface and ice grows in the direction of the local drop impingement.Since the collection e©ciency at the stagnation point is almost the same for all three cases, the ice thickness there is very similar.However, the ice extent di¨erences between three considered cases are related to the fact that drop velocities have altered vertical component due to the §ow characteristics in airfoil vicinity.Consequently, the extent of the ice accretion on lower airfoil surface is the greatest for the HQ309 airfoil and the smallest for the SD7037 airfoil.
The changes in aerodynamic performance that result from ice accretion are also depicted in Fig. 5.The clean HQ309 and SD7037 airfoils exhibit similar aerodynamic characteristics for the simulated conditions.However, once ice is formed, this pairing seems to disappear.The obtained results suggest that the iced HQ309 airfoil has better aerodynamic performance than the SD7037 airfoil under cruise conditions (low angles of attack).However, the roles are reversed near the stall angle (high angle of attack).The SD7032 airfoil appears to have the best aerodynamic performance at both cruise and near-stall conditions.However, the authors suspect that these preliminary results could be sensitive to §ight and icing conditions.Additional analysis is needed in order to reach satisfactory, generalized conclusions.

Three-Dimensional Ice Accretions Forming on Illinois Model Swept Wing
The 3D, so-called IM was used in the numerical investigation of ice formation on a swept wing.This is a comparatively small swept wing consistent with what might be used in UAV con¦guration.Figure 6 shows some details of the IM swept wing geometry.The wing is fully described in [10].The following §ow parameter values were used: angle of attack 0 • , free-stream air pressure 10 5 Pa, and free-stream air velocity 49.1 m/s.This gives a Reynolds number of 6 • 10 5 based on the mean aerodynamic chord.Figure 6 also depicts three 2-centimeter wide domains, each perpendicular to the leading edge, where the ice accretion is analyzed.These three domains are centered at 80%, 50%, and 20% of the semispan measured from the root.Several drop trajectories for these domains are also depicted.It should be noted that the drop trajectory calculations initiated ¦ve chords upwind from the impingement area.Figure 6 shows only trajectories that are close to the wing.The discrepancies between ice accretion processes at di¨erent spanwise locations are driven by di¨erent local chord length and wing shape that in §uence Figure 6 Sketches of ice accretion forming on IM wing at three spanwise 2-centimeter wide regions centered at 80%, 50%, and 20% of semispan from the root and perpendicular to the leading edge.Samples of drop trajectories (in red) are also depicted air §ow, drop trajectories, and drop impingement area and water mass.All presented calculations were obtained for a drop diameter of 25 µm and an icing duration of 10.9 min based on a travel distance of 32.2 km.Ice accretion was calculated for air temperatures of −2, −10, and −30 • C in the three IM domains shown in Fig. 6.Following the Appendix C icing envelope, each value of liquid water content, 0.46, 0.30, and 0.10 g/m 3 , was associated with a particular static air temperature.They were −2, −10, and −30 • C, respectively.
Figure 7 shows these icing predictions for a static air temperature of −2 • C. At all locations, glaze conditions prevail, with impinging drops §owing along the ice surface before freezing.However, there are obvious di¨erences between the three predicted ice accretions.Close to the root, the greater wing thickness leads to an increase in the impinging mass §ux and, consequently, to an increase in total ice mass, since there is no shedding.The ice mass in the 80, 50, and 20 percent regions is 84, 123, and 164 g/m, respectively.Since the predicted convective heat transfer coe©cient decreases along the leading edge towards the root, the ice thickness also decreases.All unfrozen water moves downstream along the upper and lower surfaces of the wing and it freezes gradually.The ice accretion extent is greater closer to the root, because of the greater impinging mass and less e©cient water freezing.The images in the left column of Fig. 7 show the three-dimensionality of the ice accretion process, with ice rivulets aligned with  the external air §ow.In the right column of Fig. 7, three ice cross sections, perpendicular to the leading edge and 1 cm apart, are depicted using di¨erent colors.The images show the variability of the ice accretion shape, even on a small scale.
For a static air temperature of −10 • C, the model predicts rime ice (Fig. 8).A characteristic compact ice formation grows in the vicinity of the leading edge and discrete rough ice features are produced further back along the wing surface.These ice features are a consequence of the shadowing e¨ect that occurs when a growing ice lobe intercepts incoming drops and the region immediately behind it is shadowed.The 3D images and 2D cross sections show the randomness of these features.They also show the growth of ¤ice ¦ngers,¥ approximately perpendicular to the external air §ow, in agreement with experimental observations.The fact that the ice cross sections can di¨er signi¦cantly highlights di©culties of verifying ice accretion using 2D cross section.The ice mass increases towards the wing root and is 55, 80, and 107 g/m, respectively.
Icing simulations for a static air temperature of −30 • C are shown in Fig. 9. Rime ice forms as in the −10 • C case, but because of the lower liquid water content, the ice accretion features are less pronounced and the ice mass is smaller.In the three domains, the ice mass is 18, 27, and 36 g/m, respectively.

CONCLUDING REMARKS
A morphogenetic ice accretion model, coupled with a CFD solver, which also predicts surface convective heat transfer, and a drop trajectory solver, has been used to simulate ice accretion on a number of airfoils, HQ309, SD7032, and SD7037, and on a 3D IM swept wing.The numerical model was used to simulate ice accretion for a range of icing conditions speci¦ed in the FAR 25 Appendix C icing envelope.The results demonstrate the ability of the morphogenetic model to predict complex, 3D features of ice accretions under rime and glaze conditions, including ice accretions on a swept wing.
Aerodynamic analysis of the HQ309 airfoil suggests that for angles of attack less than 9 • , the di¨erences between the aerodynamic penalties for the glaze and rime cases are small.However, for angles of attack exceeding 9 • , the rough surface characteristics of rime ice lead to greater aerodynamic penalties than the smoother surface of glaze ice.This occurs despite the fact that glaze covers a larger area of the airfoil.Further analysis is needed to generalize these results to an even wider range of icing conditions and to identify an airfoil geometry that will minimize the aerodynamic consequences of ice accretion.
In the near future, the authors plan to perform validation experiments in the NRC Altitude Icing Wind Tunnel, with a focus on UAVs, including swept wing con¦gurations.The experimental ice accretion shapes will be used to validate the numerical codes.In order to increase the ¦delity of the numerical prediction, also, it is intended to take into account temporal variations of the air §ow, surface heat transfer, and drop trajectories as the ice accretion evolves with time.

Figure 2
Figure 2 Distribution of the streamwise velocity at the leading edge of the iced HQ309 airfoil at angle of attack 3 • , drop size 25 µm, and air temperature −2 • C. The surface of the original airfoil is the heavy black line

Figure 7
Figure 7 Ice accretion prediction on IM wing for an airspeed 49.1 m/s, angle of attack 0 • , drop size 25 µm, and air temperature −2 • C at three spanwise locations: (a) 80%, (b) 50%, and (c) 20% semispan from the root.Three-dimensional images (left column) show 2-centimeter wide accretion domains and 2D images (right column) show their three cross sections, 1 cm apart, depicted using di¨erent colors, cyan, orange and olive being furthest from the root

Figure 8 Figure 9
Figure 8The same as Fig.7but static air temperature is −10 • C