Eulerian Method for Ice Accretion on Multiple-Element Airfoil Sections Jacco M. Hospers1 and Harry W.M. Hoeijmakers2 University of Twente, PO Box 217, 7500 AE Enschede, the Netherlands

A computational method is presented that given the flow solution computes ice accretion on multiple-element airfoils in specified icing conditions. The numerical simulation method (Droplerian) uses an Eulerian method to determine the droplet trajectories and distribution of the Liquid Water Content (LWC). To solve the equations for the droplet trajectories and liquid water content distribution, Droplerian uses a Finite Volume Method for unstructured grids. Through the droplet velocities and Liquid Water Content at the surface of the airfoil configuration the droplet catching efficiency is calculated. The droplet catching efficiency and droplet velocities at the airfoil surface are input for the icing model, which is based on Messinger’s model for ice accretion. The method includes a multi-disperse droplet distribution with an arbitrary number of droplet bins and a droplet splashing model. For a single-element airfoil a good agreement is found with measured catching efficiencies and with the ice shapes predicted by other computational methods. For increasing droplet diameter the agreement with experimental results deteriorates. The application of the method to a three-element airfoil is described. The comparison of the catching efficiency predicted by both the Droplerian method and a Lagrangian method (2DFOIL-ICE) is good. The agreement of predicted ice accretions with available experimental data is reasonable.

Nomenclature D fD g n U A c CD d f I K Ky Kc,dry Ky,crit LWC M N Nbin Oh Rnd 1 2

Drag force [N] Drag force per unit mass [N/kg] Gravitational acceleration vector [m/s2 ] Unit normal vector [-] Local droplet velocity [m/s] Cross-sectional area [m2 ] Chord length [m] Drag coefficient of a droplet [-] Droplet diameter [m] Liquid Water Content fraction for a single droplet binh [-] i Splashing term in momentum conservation equation mkg 2 s2 Cossali’s splashing parameter [-] Yarin and Weiss’ splashing parameter [-] Trujillo’s splashing threshold for a dry surface [-] Yarin and Weiss’ splashing threshold [-] Cloud liquid water content [kg/m3 ] h i Splashing term in mass conservation equation mkg3 s Number of secondary droplets [-] Total number of droplet bins Droplet Ohnesorge number [-] Non-dimensional surface roughness [-]

PhD student, faculty of Engineering Technology, group of Engineering Fluid Dynamics Professor, faculty of Engineering Technology, group of Engineering Fluid Dynamics, senior member AIAA 1 American Institute of Aeronautics and Astronautics

Reynolds number based on the relative droplet velocity [-] Airfoil coordinate, clockwise positive [m] Volume [m3 ] Droplet Weber number [-]

Red s V We

Subscripts ∞ a d s splash w

Free stream quantity Air related quantity Droplet related quantity Secondary droplet related quantity Splashing related quantity Water related quantity

Symbols α Volume fraction of water contained in droplets [-] β Local catching efficiency [-] µ Local dynamic viscosity [Pa s] φ Mass-loss-coefficient [-] ρ Local droplet density [kg/m3 ] σ Surface tension [N/m]

Introduction

A

 icing has long been recognized as a serious flight safety problem. According to Ref. 1, in the period 1992–2000, airframe icing has been attributed to more than 50 accidents and incidents, claiming more than 800 lives, in the US alone. Icing occurs when super-cooled water droplets hit the aircraft, flying at a level where the temperature is at or below the freezing point. Ice accretion on the wing leading edge or on the tail plane can result in non-aerodynamic shapes and in serious degradation of the aerodynamic performance, such as a decrease in the stall angle, an increase in drag, a decrease in maximum lift, and altered moment characteristics of the aircraft. Also, ice accretion on parts of the engine nacelles or on propellers can cause dangerous situations when the accretions break free. Specifically Supercooled Large Droplets or SLD are relatively dangerous, the shapes of ice accretions caused by SLD often differ from that of conventional ice accretions in that they are positioned further downstream of the leading-edge of the airfoil and SLD ice accretions tend to grow very quickly. Numerical simulation of the ice accretion process provides an attractive method for determining the ice shapes on aircraft wings and evaluating a wide range of icing conditions. An ice accretion model that accurately predicts growth shapes on an arbitrary airfoil is valuable for analysis of the sensitivity of airfoils to ice accretion and for analysis of the influence of variables such as airspeed and angle of attack, pressure, temperature, humidity, droplet size, etc. on the accretion process. The predicted ice shapes can be used in wind tunnel and flight tests to assess aircraft performance and handling qualities degradation in icing conditions. Nowadays, it is common practice in the aircraft manufacturing industry to apply computational methods to predict ice accretion in two-dimensional flow. Studies to extend the two-dimensional ice growth method to three-dimensional flows are in progress at for example NASA GRC as well as at CIRA and ONERA. At the University of Twente also an ice accretion prediction method has been developed. This method, designated 2DFOIL-ICE,2–5 predicts the growth of ice on single and multi-element airfoil sections. In this method the droplet trajectories are calculated using a Lagrangian method. Ice accretion modeling is based on the Messinger model, a quasi-steady model that takes into account all important mass and heat transfer processes that occur when super-cooled water droplets strike an airfoil. The droplets either freeze immediately upon impact or freeze partly while the rest of the water runs back on the airfoil. SLD ice accretions can occur downstream of the anti-icing equipment fitted on many aircraft, presenting a little understood danger. SLD icing conditions are currently outside of the flight certification envelope, but a new certification process is being designed and there is a need for numerical tools to predict ice accretions in the SLD-regime. Due to the size of droplets in the SLD-regime (d > 50 µm) versus droplets in the normal-regime (d ≤ 50 µm) some additional effects have to be modeled in order to accurately predict the trajectories of SLD:6 • Droplet deformation, resulting in a change in drag force • Droplet break up, resulting in different droplet trajectories 2 American Institute of Aeronautics and Astronautics

• Droplet splashing and/or rebound, resulting in a different deposition rate The present work describes the development of an Eulerian droplet tracking method, used together with the icing model from 2DFOIL-ICE. This development is motivated by a Lagrangian computation, using several droplet bins for multiple-element airfoils, taking significant amounts of computer time. Furthermore, an Eulerian model can be more easily extended to 3D. Both the original Lagrangian droplet tracking model and the new Eulerian droplet tracking model have been extended from mono-disperse to multi-disperse droplet distributions, while also a model for splashing droplets has been added. Numerical results from the Eulerian droplet tracking model are compared to experimental data available from literature and results obtained with 2DFOIL-ICE.

Eulerian Droplet Tracking The 2DFOIL-ICE method is based on Lagrangian droplet tracking, in which the (potential) flow field is calculated using a panel-method. This has some limitations, especially for the multiple-element airfoil cases. Due to the potentialflow model in 2DFOIL-ICE, the geometry of multiple-element airfoils often has to be approximated. Only a single sharp edge is permitted on each element, requiring the addition of cove bounding streamlines. A second problem can be the process of finding droplet trajectories that impinge on one of the airfoil elements. This can become a time consuming task. To reduce both of these limitations an Eulerian droplet tracking method has been developed. The Eulerian method can obtain flow field data from any available flow model (e.g. potential-flow such as used in 2DFOIL-ICE, but also Euler or full Navier-Stokes). This flow field is used as input to calculate a droplet flow field on a grid around the airfoil. For the Eulerian method, discrete droplet trajectories do not need to be calculated, reducing the computational load, while allowing a detailed investigation of the droplet variables close to the airfoil. A further reason for developing this Eulerian method is the easy extension towards three-dimensional geometries. The output needed from the droplet tracking method should be suited as input for the icing model. For the icing model only two quantities related to the impinging droplets are needed: the rate of mass of water impinging locally on the airfoil surface, and the rate of kinetic energy that is contained within the impinging droplets. The local rate of mass can be obtained from the local catching efficiency β. From the velocity of the impinging droplets the rate of kinetic energy can be calculated. For Lagrangian methods the local catching efficiency is defined as the ratio of the rate of impinging mass divided by the rate of impinging mass had the droplet trajectories been straight lines. For a Lagrangian method, β is calculated as: β=

dy ds

(1)

This is illustrated in Fig. 1, the mass of the droplets between two droplet trajectories remains constant. Given a fixed value of dy, the smaller the contour element ds the larger the rate of impinging mass and therefore, β. For an Eulerian method, since we do not (necessarily) compute individual droplet trajectories, this approach cannot be applied. Using the liquid water content (LWC) of the cloud, the following relation depending only on the local droplet density ρd and local droplet velocity Ud can be derived: β=

ρd Ud · n LWC Ud,∞

(2)

where the local droplet density ρd is actually the volume fraction of water contained in the droplets α multiplied with the local water density ρw : ρd = αρw (3) Equation 2 contains the impinging mass rate of liquid water per square meter ρd Ud · n (see Fig. 2), which can be used directly as input for the icing model. To calculate both ρd and Ud the droplets are considered as a second fluid phase. Solving the mass and momentum balances for a discretized domain provides these local quantities. One of the source terms for the momentum balance will have to be the drag force, since this is the main driving force for the droplet phase. The equations for mass conservation and momentum conservation of the droplet phase are: ∂ρd + ∇ · ρd Ud = 0 ∂t 3 American Institute of Aeronautics and Astronautics

(4)

n

s

Ud ds dy

Figure 1. Definition of local catching efficiency β for Lagrangian methods

Figure 2. Droplet velocity U d and surface normal n

! ∂ρd Ud ρa + ∇ · (ρd Ud ) Ud = ρd f D − ρd 1 − g ∂t ρw

(5)

where the only other source-term considered besides the drag force f D is due to gravity. For the present case other forces such as lift force and Basset-history force are neglected. The drag force is modeled using the drag coefficient C D : fD =

D ρa |Ua − Ud | (Ua − Ud ) Ad C D = ρw Vd 2ρw Vd

(6)

where C D is usually a function of the Reynolds number based on the relative droplet velocity Red : Red ≡

ρa |Ua − Ud | d µa

(7)

The expression for the drag coefficient can range from an expression for small diameter droplets to special relations for deforming droplets (SLD diameter droplets). In the current model C D is derived from the Langmuir and Blodgett law:7 C D Red = 1 + 0.0197Re0.63 + 2.6 · 10−4 Re1.38 d d 24

(8)

which is valid for Red < 1000. The resulting equations are solved numerically by employing a finite volume method based on similar numerical flow simulation methods.8, 9 This edge-based finite volume method uses an unstructured grid and is suited for both two-dimensional and three-dimensional domains, discretized using any combination of element types. The flow field of the surrounding air is used only as input, since one-way coupling between the air and droplet phase is assumed. It can be obtained from any available flow solver. With respect to the boundary conditions, the boundary condition on the airfoil requires special attention. In order to calculate a catching efficiency the normal component of the droplet velocity at the airfoil is needed. However, when the normal component of the droplet velocity is negative (using the definition from Fig. 2), i.e. away from the airfoil, droplets should not be created from the inside of the airfoil. This leads to a boundary condition on the surface of the airfoil, which depends on the sign of the normal component of the droplet velocity. For a normal component of the droplet velocity away from the airfoil the following boundary conditions are applied: Ud · n ≥ 0 :

ψwall = 0

(9)

where ψ is the variable of the considered governing equation; either ρd or ρd Ud . In case the normal component of the droplet velocity is directed into the airfoil the values at the wall ψwall are calculated by extrapolation from the values in the first finite volume next to the airfoil surface ψ0 : Ud · n < 0 :

ψwall = ψ0 + (xwall − x0 ) · ∇ψ0

(10)

The resulting equations are solved for a multi-disperse droplet distribution. For each droplet-bin in the distribution a separate set of equations is solved. The results from each droplet-bin are combined in a single β, which is subsequently used in the icing model, using the fraction of the total LWC f . The sum of f for all droplet-bins has to equal one: Nbin X

fi = 1.0

i=1

4 American Institute of Aeronautics and Astronautics

(11)

with Nbin the total number of droplet-bins, such that β=

Nbin X

fi βi

(12)

i=1

Splashing Model For traditional icing models (non SLD-regime), it is common to assume that every droplet that impinges on the airfoil surface is fully deposited onto that surface. However, several other deposition regimes are possible, such as splashing and rebound of impinging droplets. One of the major differences between SLD and smaller droplets is thought to be caused by splashing. Due to splashing only part of the water contained in the droplet will be deposited on the airfoil surface, while the remainder of the water will be splashed away from the airfoil, possibly re-impinging further downstream of the original impact. By including a model for the mass-loss due to splashing, a mass-loss-coefficient 0 ≤ φ < 1 can be determined, which reduces the catching efficiency: β splash = (1 − φ) β

(13)

Models for droplet splashing often also predict the number of secondary droplets and their velocity after the splashing event (see Fig. 3), making it possible to track the droplets further downstream, until possible re-impingement. The basic model for mass-loss due to splashing from Trujillo, Matthews, Lee and Peters10 causes a severe overprediction of the mass-loss coefficient. For improved performance, a version of this basic model calibrated for mass-loss coefficients in SLD conditions by Honsek, Habashi and Aub´e11 is used. The resulting mass-loss coefficient is expressed as:   h  i 3.8  φ Ky = p 1 − exp −0.85 Ky − 17 (14) Ky where Ky is the splashing parameter defined by Yarin and Weiss12 as: Ky = Λ

−3/8



−2/5

Oh

We

5/16

 ! −3/8  5/16  3 LWC 1/3    Oh−2/5 We =  2 ρw

(15)

with the droplet Weber number We defined as: We = and the Ohnesorge number Oh as:

ρa |Ud |2 d σd

√ We µa Oh = p = Re d ρa σd d

(16)

(17)

In Eq. 14 the splashing threshold from Yarin and Weiss appears; Ky,crit = 17. Splashing only occurs for Ky > 17 as illustrated in Fig. 4. Also note that, within 1% accuracy, φ = √3.8 for Ky > 23. Ky

For typical conditions (LWC ≈ 0.5 · 10−3 kg/m3 , ρa ≈ 1.225 kg/m3 , µa ≈ 1.789 · 10−5 Pa s, ρw ≈ 1000 kg/m3 , σd ≈ 0.072 N/m, Ud ≈ 100 m/s, 0 m/s ≤ |Ua − Ud | ≤ 100 m/s) and diameters ranging from 10 µm ≤ d ≤ 1000 µm the following ranges for the relevant dimensionless numbers can be found: We ∈ [0, 1000] Red ∈ [0, 10000] Oh Ky

∈ [0, 0.1] ∈ [0, 100]

The velocity and the number of secondary droplets after the splashing event, which are only needed to follow secondary droplets, are taken from Trujillo et al.’s original expressions. They determined the following relation for the 5 American Institute of Aeronautics and Astronautics

√3.8 Mass-loss coefficient, φ

Ky

N=3 n

d

ds

Ud,s

Ud

Ky,crit = 17

1.0 0.8 0.6 0.4 0.2 0.0

Figure 3. Variables related to a splash event

0

10 20 30 40 50 60 70 80 90 100 Splashing parameter, Ky

Figure 4. Mass-loss coefficient as a function of splashing parameter

average number of secondary droplets:     !2   1  |Ud |  N= − Kc,dry  − 44.92 0.0437 K 22 Ud · n

(18)

where K is the splashing parameter defined by Cossali, Coghe and Marengo,13 which is similar to Eq. 15: K = Oh−2/5 We

(19)

The splashing threshold for a dry surface Kc,dry is defined by Trujillo et al. as a function of the non-dimensional surface roughness Rnd as: Kc,dry = 180R−0.035 (20) nd Trujillo et al. also determined relations for the average secondary droplet velocities: ! Ud,s · t Ud · t = 0.85 + 0.0025 arctan Ud · t Ud · n Ud,s · n Ud · t = 0.12 + 0.002 arctan Ud · n Ud · n

(21a)

! (21b)

From mass conservation and Eq. 18 the average size of secondary droplets can be calculated: ds =

 φ 1/3 N

d

(22)

For a Lagrangian method the information obtained from Eqs. 14, 21 and 22 is sufficient to correct β and follow the secondary droplet trajectories downstream. However, for an Eulerian method an extra step has to be taken if the secondary droplets are to be re-injected into the flow. A splashing boundary condition has to be defined, which at the airfoil surface adds a term to the flux in both the mass and momentum equations, corresponding to the mass and momentum of the secondary droplets, respectively. From the local values of ρd and Ud and from Eq. 14 the local re-injected mass rate can be computed. This leads to the following term, which is only included when the considered finite-volume is next to the wall and experiences splashing:     φ (ρd Ud · n) > 0 : Splashing M= (23)   0 : No splashing

6 American Institute of Aeronautics and Astronautics

The re-injected momentum then follows from the re-injected mass rate (Eq. 23) and Eq. 21, leading to the following source-term, which again is only included when splashing occurs in the considered finite-volume element.     φ (ρd Ud · n) Ud,s = MUd,s : Splashing I= (24)   0 : No splashing The calculation is performed using a multi-disperse droplet distribution. The droplet-bin in which the flux-terms from Eqs. 23 and 24 are included is chosen such that the diameter is closest to the diameter of the secondary droplets calculated in Eq. 22. With the inclusion of terms M and I the following equations are obtained: ∂ρd,i + ∇ · ρd,i Ud,i = −Mi ∂t

(25)

! ∂ρd,i Ud,i
ρa + ∇ · ρd,i Ud,i Ud,i = ρd,i f D,i − ρd,i 1 − g − Ii ∂t ρw

(26)

where i = 1, 2, . . . , Nbin indicates the considered droplet-bin.

Results

0.3 0.2 0.1 0.0 0.0

5

10 15 20 25 30 35 40 45 Droplet diameter, d [µm] (a) 20 µm MVD

Fraction of total LWC, f

Fraction of total LWC, f

In order to validate the Eulerian droplet tracking method a suitable test case had to be selected. Because of the availability of 2DFOIL-ICE, which calculates a potential-flow field using a panel-method in order to carry out Lagrangian droplet tracking, a validated ice-accretion prediction code, the results from the Eulerian droplet tracking method are compared to the results from 2DFOIL-ICE. This implies that the underlying flow field used for the Eulerian droplet tracking is also obtained from the potential-flow field as calculated by 2DFOIL-ICE. Besides this comparison between two computational results, a comparison with experimental impingement data was performed. For this comparison experimental data is available through an experiment of Papadakis et al.14 was selected. Papadakis performed experiments with different droplet median volume diameters (MVD) for a NACA-23012 airfoil at 2.5◦ angle of attack (AoA). The impingement data is presented in the form of a catching efficiency and the LWC distribution has been recorded so that it can be used in a multi-disperse simulation. Two cases were selected from the range of MVD’s, the smallest and the largest MVD, to investigate the ability to predict both normal and SLD impingement regimes. The selected cases and the corresponding conditions are shown in Table 1. A first assessment was made by performing a numerical simulation for both the 20 µm and 236 µm MVD, for both the Lagrangian 2DFOIL-ICE and the Eulerian Droplerian, with a mono-disperse droplet distribution. A single droplet class in which the droplets have a diameter equal to the MVD was used. To test the addition of multi-disperse droplet bins, a second calculation was performed. Papadakis et al. analyzed the measured droplet distributions and divided them into 10 and 27 binary droplet-bins. For the following simulations 10 droplet-bins are used. The 10-bin distributions for both MVD cases are shown in Table 2 and illustrated in Fig. 5. 0.3 0.2 0.1 0.0

0.0

200 100

400 300

800

600 500

700

1000 900 1100

Droplet diameter, d [µm] (b) 236 µm MVD

Figure 5. 10-Bin droplet distributions for selected cases 7 American Institute of Aeronautics and Astronautics

Table 2. 10-Bin droplet distributions for selected cases Droplet bin 1 2 3 4 5 6 7 8 9 10

Table 1. Conditions for selected cases 20 µm MVD AoA c LWC Ua,∞ T∞ p∞

236 µm MVD 2.5◦ 0.9144 m 0.19 g/m3 1.89 g/m3 78.23 m/s 299 K 101330 Pa

LWC [%] 5 10 20 30 20 10 3 1 0.5 0.5

Droplet size [µm] 20 µm MVD 236 µm MVD 3.850397 16.25037 9.390637 63.65823 13.80175 135.4827 19.60797 298.5197 25.4820 508.4572 30.73474 645.4684 35.19787 715.8689 38.32569 747.3936 40.66701 763.2455 44.36619 1046.767

3 y c

2 1 0 -1 -2 -3 -5

-4

-3 -2 -1 0 1 2 Dimensionless x-coordinate, cx (a) Entire domain

3

4

Dimensionless y-coordinate,

Dimensionless y-coordinate,

y c

The number of panels used to obtain the velocity field is equal to the number of surface-elements on the airfoil for the Eulerian method: 400. The domain used for the Eulerian calculation is a rectangular box, discretized with 8767 triangular elements, as shown in Fig. 6. A median-dual mesh was created, in which a control volume element is created for each vertex in the original triangular mesh.

0.5

0.0

-0.5 -0.5

0.0 1 0.5 Dimensionless x-coordinate,

1.5 x c

(b) Close-up around the airfoil

Figure 6. Grid used for Eulerian computations The resulting catching efficiencies are shown in Figs. 7 and 8. For the 20 µm MVD results it can be observed that the sharp impingement limits present in the mono-disperse solution are flattened and move outwards for the multi-disperse solution (Fig. 7(b) and Fig. 8(b)). This is due to the inclusion of droplet-bins with a smaller diameter than the MVD, leading to droplet trajectories that better follow the streamlines. The inclusion of multiple droplet-bins leads to a clear improvement when the solution is compared to the experimental result. However, for the large MVD case of 236 µm, the improvements are not so clear as for the small MVD case of 20 µm because the large droplets follow almost ballistic trajectories. The over-prediction for the large MVD case is thought to be due to splashing and rebound effects. When comparing the Lagrangian solution to the Eulerian solution it can be concluded that they are very similar. It should be noted that the catching efficiency that was computed with the Eulerian method is slightly lower than the catching efficiency from the Lagrangian method, especially around the nose of the airfoil (s ≈ 0). This is most likely due to numerical dissipation in the employed finite volume scheme. Overall, the agreement with the experimental results appears similar for the Eulerian and the Lagrangian results. 8 American Institute of Aeronautics and Astronautics

0.6

Droplerian, 10 droplet-bins Droplerian, 1 droplet-bin Experimental results

0.6 0.5 Catching efficiency, β

Catching efficiency, β

0.5

Droplerian, 10 droplet-bins 2DFOIL-ICE, 10 droplet-bins Experimental results

0.4 0.3 0.2 0.1

0.4 0.3 0.2 0.1

0.0 -0.3

-0.2

-0.1 bottom

0.0 nose

Dimensionless airfoil-coordinate, (a) Droplet distribution effects

0.1

0.0 -0.3

-0.2

top s−s stag c

-0.1 bottom

0.0 nose

Dimensionless airfoil-coordinate,

0.1 top s−s stag c

(b) Eulerian (Droplerian) and Lagrangian (2DFOIL) results

Figure 7. Calculated catching efficiencies, 20 µm MVD

For the catching efficiencies calculated with the multi-disperse Eulerian method shown in Figs. 7 and 8 corresponding ice accretion shapes were calculated. The impingement tests were performed at a temperature of 299 K, too high for any icing to occur. Therefore, the ice accretions have been calculated at a temperature of 268.15 K (-5 ◦ C) and at the LWC of the 20 µm MVD case: 0.19 g/m3 as shown in Table 3. The time-step used for the calculation of the ice accretions is 360 s. The resulting ice shapes are shown in Fig. 9. The larger catching efficiency of the 236 µm MVD case results in a larger ice accretion. The thickness of the ice accretions are similar around the leading edge. The amount of energy impinging on the airfoil surface only differs due to the kinetic energy of the impinging droplets. This amounts to a very small change for the area around the leading edge where the mass of ice is not limited by the mass of impinging water. The ice shape for the 236 µm MVD case extends further downstream of the leading edge, while the ice accretion for the 20 µm MVD case is only present around the leading edge. Splashing Effects The splashing model discussed in the previous section was implemented in both 2DFOIL-ICE and Droplerian. Implementing re-injection of droplets is quite complicated for the Lagrangian method, due to the calculation of the catching efficiency β. Re-injection is therefore only implemented in the Eulerian method. The results for the catching efficiencies are shown in Fig. 10(a). When comparing the catching efficiency with splashing from Droplerian with the catching efficiency without splashing and the experimental results, a dramatic improvement is observed. Near the leading-edge the under-prediction increased, but the match with experimental results further downstream of the leading-edge increased significantly. For the icing shape calculated from the catching efficiency with splashing, shown in Fig. 10(b), a decrease of the ice thickness can be noted downstream of the leading-edge. The icing shape near the leading-edge is very similar to the icing shape determined without a splashing model. The mass of ice that freezes near the leading-edge is similar, however, the amount of run-back water for the case with the splashing model is smaller, caused by the lower catching efficiency, so the edge of the horns occurs closer to the leading edge than for the case without splashing. Since the Eulerian model provides information of the full field a closer look at the droplet density distribution can be taken. The resulting droplet distributions for the largest droplet bin (bin 10, 1046 µm, see Table 2) are shown in Fig. 11. For this droplet bin there is no visual difference in the droplet distributions with and without splashing. Droplets that splash are removed from the flow (just as they would without splashing) and injected into a smaller bin. To illustrate this the result for the smallest droplet bin (bin 1, 16 µm) is shown in Fig. 12. Droplets from the bins with larger droplet 9 American Institute of Aeronautics and Astronautics

1.0 0.9

Droplerian, 10 droplet-bins Droplerian, 1 droplet-bin Experimental results

0.9 0.8 Catching efficiency, β

0.8 Catching efficiency, β

Droplerian, 10 droplet-bins 2DFOIL-ICE, 10 droplet-bins Experimental results

1.0

0.7 0.6 0.5 0.4 0.3

0.7 0.6 0.5 0.4 0.3

0.2

0.2

0.1

0.1

0.0 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 bottom

0.1

nose

Dimensionless airfoil-coordinate, (a) Droplet distribution effects

0.2

0.0 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0

top

bottom

s−s stag c

nose

Dimensionless airfoil-coordinate,

0.1

0.2

top s−s stag c

(b) Eulerian (Droplerian) and Lagrangian (2DFOIL) results

Figure 8. Calculated catching efficiencies, 236 µm MVD

diameters splash, creating small secondary droplets with a diameter close to 16 µm. These droplets are injected from the airfoil surface into the flow (Fig. 12(b) compared to Fig. 12(a)). Figure 13 provides a closer look around the leading edge showing a region of increased droplet concentration, with droplets being convected downstream.

Concluding Remarks An Eulerian method to calculate ice-accretions on two dimensional airfoils has been created. The method was based on a similar Lagrangian method, which produced satisfactory results for small droplet diameters. For larger (SLD) droplets, both methods over-predict the catching efficiency compared with experimental data. To improve the matching with experimental catching efficiency data, a splashing model has been added to the Eulerian method. Splashing is accounted for by introducing a mass-loss-coefficient and by re-injecting secondary droplets into the droplet bin corresponding to the diameter closest to the secondary droplet diameter obtained from the splashing model. The inclusion of a splashing model accounts for a decrease in the catching efficiency on the airfoil, bringing the catching efficiency results closer to experimental results. This is a clear improvement compared to the case without the splashing model. Due to the decreased catching efficiency a smaller ice-accretion is formed around the leading edge of the considered airfoil. The experimentally and numerically calculated catching efficiencies still differ, which is possibly caused by effects due to rebound and breakup. Models for these effects are considered for future improvements to the method. Improved prediction of the flow around the airfoil, by using more advanced flow-models instead of the currently applied potential-flow-model, might also improve the results.

References 1

Petty, K. R. and Floyd, C. D. J., “A statistical review of aviation airframe icing accidents in the US,” 11th Conference on Aviation, Range, and Aerospace, Hyannis, MA, October 2004. 2 Jacobs, S. J., Hospers, J. M., and Hoeijmakers, H. W. M., “Numerical Simulation of Ice Accretion on Multiple-Element Airfoil Sections,” ICAS 2008, 26th international congress of the aeronautical sciences, Anchorage, AK, September 2008. 3 Dillingh, J. E. and Hoeijmakers, H. W. M., “Numerical Simulation of Airfoil Ice Accretion and Thermal Anti-Icing Systems,” ICAS 2004, 24th international congress of the aeronautical sciences, Yokohama, Japan, August-September 2004. 4 Dillingh, J. E. and Hoeijmakers, H. W. M., “Simulation of ice accretion on airfoils during flight,” FAA In-Flight Icing/Ground

10 American Institute of Aeronautics and Astronautics

Dimensionless y-coordinate,

y c

0.08

ice-shape, 20 µm MVD ice-shape, 236 µm MVD

0.06 0.04

Table 3. Conditions for icing cases

0.02

AoA c LWC Ua,∞ T∞ p∞

0.0 -0.02 -0.04

2.5◦ 0.9144 m 0.19 g/m3 78.23 m/s 268.15 K 101330 Pa

-0.06 -0.02 0.0 0.02 0.04 0.06 0.08 0.10 0.12 Dimensionless x-coordinate, cx Figure 9. Calculated ice accretions

De-Icing International Conference & Exhibition, Chicago, IL, June 2003, 12 pages. 5 Snellen, M., Boelens, O. J., and Hoeijmakers, H. W. M., “A computational method for numerically simulating ice accretion,” 15th AIAA Applied Aerodynamics Conference, Atlanta, GA, June 1997, Technical Papers pt 2. 6 Wright, W. and Potapczuk, M., “Semi-Empirical Modelling of SLD Physics,” 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, January 2004. 7 Langmuir, I. and Blodgett, K. B., “A mathematical investigation of water droplet trajectories,” Tech. Rep. 5418, US Army Air Forces, 1946. 8 Koop, A. H., Numerical simulation of unsteady three-dimensional sheet cavitation, Ph.D. thesis, University of Twente, Enschede, the Netherlands, September 2008. 9 Kelleners, P. H., An edge-based finite volume method for inviscid compressible flow with condensation, Ph.D. thesis, University of Twente, Enschede, the Netherlands, December 2007. 10 Trujillo, M. F., Mathews, W. S., Lee, C. F., and Peters, J. F., “Modelling and Experiment of Impingement and Atomization of a Liquid Spray on a Wall,” International Journal of Engine Research, Vol. 1, No. 1, 2000, pp. 87–105. 11 Honsek, R., Habashi, W. G., and Aub´e, M. S., “Eulerian Modeling of In-Flight Icing Due to Supercooled Large Droplets,” Journal of Aircraft, Vol. 45, No. 4, July-August 2008, pp. 1290–1296. 12 Yarin, A. L. and Weiss, D. A., “Impact of drops on solid surfaces: self-similar capilary waves, and splashing as a new type of kinematic discontinuity,” Journal of Fluid Mechanics, Vol. 283, January 1995, pp. 141–173. 13 Cossali, G. E., Coghe, A., and Marengo, M., “The impact of a single drop on a wetted solid surface,” Experiments in Fluids, Vol. 22, No. 6, April 1997, pp. 463–472. 14 Papadakis, M., Rachman, A., Wong, S.-C., Yeong, H.-W., Hung, K. E., Vu, G. T., and Bidwell, C. S., “Water Droplet Impingement on Simulated Glaze, Mixed and Rime Ice Accretions,” Tech. Rep. NASA/TM-2007-213961, NASA, October 2007.

11 American Institute of Aeronautics and Astronautics

1.0 0.9

Droplerian, splashing Droplerian, no splashing Experimental results

0.08 y c

0.7

Dimensionless y-coordinate,

Catching efficiency, β

0.8

0.6 0.5 0.4 0.3 0.2 0.1 0.0 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 bottom

nose

Dimensionless airfoil-coordinate,

0.1 top s−s stag c

0.2

ice-shape, splashing ice-shape, no splashing

0.06 0.04 0.02 0.0 -0.02 -0.04 -0.06 -0.02 0.0 0.02 0.04 0.06 0.08 0.10 0.12 Dimensionless x-coordinate, cx

(a) Catching effiencies

(b) Ice accretions

Figure 10. Results 236 µm MVD with splashing effects

(a) Without splashing model

(b) With splashing model Figure 11. Calculated droplet distributions [kg/m3 ], 236 µm MVD, bin 10 (1046 µm) 12 American Institute of Aeronautics and Astronautics

(a) Without splashing model

(b) With splashing model Figure 12. Calculated droplet distributions [kg/m3 ], 236 µm MVD, bin 1 (16 µm)

(a) Without splashing model

(b) With splashing model

Figure 13. Calculated droplet distributions [kg/m3 ], region near the leading edge, 236 µm MVD, bin 1 (16 µm)

13 American Institute of Aeronautics and Astronautics