Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 17.1.1 - Eulerian Method for Ice Accretion on Multple-Element Airfoils
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00412
 Material Information
Title: 17.1.1 - Eulerian Method for Ice Accretion on Multple-Element Airfoils Droplet Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Hospers, J.M.
Hoeijmakers, H.W.M.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: Lagrangian simulations
Eulerian simulations
droplet dynamics
icing
 Notes
Abstract: A computational method is presented that computes ice accretion on multiple-element airfoils in specified icing conditions given the flow solution. The numerical simulation method (Droplerian) uses an Eulerian method to determine the droplet trajectories and distribution of the Liquid Water Content (LWC). This Eulerian method is based on a Lagrangian ice accretion method 2DFOIL-ICE described in Jacobs et al. (2008). Droplerian employs a Finite Volume Method for unstructured grids based on Koop (2008) and Kelleners (2007). Through the droplet velocities and Liquid Water Content at the surface of the airfoil configuration the droplet catching e ciency is calculated which is used as input for the icing model, based on Messinger’s model for ice accretion. The method includes multi-disperse droplet distributions with an arbitrary number of droplet bins and a droplet splashing model based on Honsek et al. (2008). Good agreement is found with experimental catching e ciencies and with ice shapes predicted by other computational methods. Agreement with the experimental ice shapes is fair. The catching e ciency predicted by both Droplerian and 2DFOIL-ICE is very similar.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00412
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1711-Hospers-ICMF2010.pdf

Full Text



7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Eulerian Method for Ice Accretion on Multiple-Element Airfoils


J.M. Hospers* and H.W.M. Hoeijmakers*

Faculty of Engineering Technology, group of Engineering Fluid Dynamics, University of Twente, the Netherlands
j.m.hospers @utwente.nl and h.w.m.hoeijmakers @ utwente.nl
Keywords: Lagrangian simulations, Eulerian simulations, droplet dynamics, icing




A computational method is presented that computes ice accretion on multiple-element airfoils in specified icing
conditions given the flow solution. The numerical simulation method (Droplerian) uses an Eulerian method to determine
the droplet trajectories and distribution of the Liquid Water Content (LWC). This Eulerian method is based on a
Lagrangian ice accretion method 2DFOIL-ICE described in Jacobs et al. (2008). Droplerian employs a Finite Volume
Method for unstructured grids based on Koop (2008) and Kelleners (2007). Through the droplet velocities and Liquid
Water Content at the surface of the airfoil configuration the droplet catching efficiency is calculated which is used as
input for the icing model, based on Messinger's model for ice accretion. The method includes multi-disperse droplet
distributions with an arbitrary number of droplet bins and a droplet splashing model based on Honsek et al. (2008).
Good agreement is found with experimental catching efficiencies and with ice shapes predicted by other computational
methods. Agreement with the experimental ice shapes is fair. The catching efficiency predicted by both Droplerian and
2DFOIL-ICE is very similar.


Nomenclature


Roman Symbols
D Drag force [N]
fD Drag force per unit mass [N/kg]
g Gravitational acceleration vector [m/s2]
1T Unit normal vector [-]
7 Local droplet velocity [m/s]
A Cross-sectional area [m2]
c Chord length [m]
CD Drag coefficient of a droplet [-]
d Droplet diameter [m]
f Liquid Water Content fraction for a single
droplet bin [-]
I Splashing term in momentum conservation equa-
tion [mk]
K Cossali's splashing parameter [-]
Ky Yarin and Weiss' splashing parameter [-]
Kc,dy Trujillo's splashing threshold for a dry surface [-]
Ky,cri, Yarin and Weiss' splashing threshold [-]
LWC Cloud liquid water content [kg/m3]
M Splashing term in mass conservation equa-
tion [k]
N Number of secondary droplets [-]
Nbin Total number of droplet bins [-]
Oh Droplet Ohnesorge number [-]
Rd Non-dimensional surface roughness [-]
Red Reynolds number based on the relative droplet
velocity [-]


s Airfoil coordinate, clockwise positive [m]
V Volume [m3]
We Droplet Weber number [-]
Greek Symbols
a Volume fraction of water contained in drop-
lets [-]
p Local catching efficiency [-]
p Local dynamic viscosity [Pa s]
0 Mass-loss-coefficient [-]
p Local droplet density [kg/m3]
r- Surface tension [N/m]
Subscripts
co Free stream quantity
a Air related quantity
d Droplet related quantity
s Secondary droplet related quantity
splash Splashing related quantity
w Water related quantity

Introduction

Aircraft icing has long been recognized as a serious flight
safety problem. According to Petty and Floyd (2004), 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 seri-
ous 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 charac-
teristics 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 dan-
gerous, 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 pro-
cess 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 accu-
rately 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, tempera-
ture, 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 man-
ufacturing 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 Univer-
sity of Twente also an ice accretion prediction method
has been developed. This method, designated 2DFOIL-
ICE (Jacobs et al. 2008; Dillingh and Hoeijmakers 2004,
2003; Snellen et al. 1997), 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 systems 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 formulated 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 Im) versus droplets in the normal-regime
(d < 50 Im) some additional effects have to be modeled
in order to accurately predict the trajectories of SLD
(Wright and Potapczuk 2004):


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Droplet deformation, resulting in a change in drag
force

Droplet break up, resulting in different droplet tra-
jectories

Droplet splashing and/or rebound, resulting in a dif-
ferent deposition rate

The present work describes the development of an Eu-
lerian 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 signif-
icant amounts of computer time. Furthermore, an Eu-
lerian 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 distribu-
tions, while also a model for splashing droplets has been
added. Numerical results from the Eulerian droplet track-
ing 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, espe-
cially for the multiple-element airfoil cases. Due to the
potential-flow 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 accepts 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 calcu-
late a droplet velocity and density field on a grid around
the airfoil. For the Eulerian method, discrete droplet tra-
jectories do not need to be calculated, reducing the com-
putational 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 associated with the impinging droplets. The local











rate of mass can be obtained from the local catching
efficiency P. 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 trajecto-
ries been straight lines. For a Lagrangian method, 3 is
calculated as:
dy
=d-7 (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, 3.









Figure 1: Definition of local catching efficiency P for
Lagrangian methods

For an Eulerian method, since we do not (necessar-
ily) 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 Pd and local droplet velocity Ud
can be derived:

Pd d.I
p= (2)
LWC Ud,o,

where the local droplet density Pd is actually the volume
fraction of water contained in the droplets a multiplied
with the local water density pw:

Pd = apw (3)

Equation (2) contains the impinging mass rate of liquid
water per square meter Pd d n (see Fig. 2), which can
be used directly as input for the icing model.



Ud





Figure 2: Droplet velocity Ud and surface normal ft

To calculate both Pd and Od the droplets are considered
as a second fluid phase. Solving the mass and momentum


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


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:

a- + PdUd = 0 (4)


SV (Pdd)Ud = PdfD + Pd(l


P) g (5)


where the only other source-term considered besides the
drag force fD is due to gravity. For the present case other
forces such as lift force and Basset history force are ne-
glected.
The drag force is expressed in terms of the drag coeffi-
cient CD:


f--
d
Pw -
Pw Vd


Pa Ua d (a Ud)AdCD
2pwVd


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

Paa U-d d
Red (7)
Pa
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 CD is derived from Langmuir and Blodgett
(1946):

-d= 1 + 0.0197Re.63 + 2.6* 10-4Re138 (8)
24 d d
which is valid for Red < 1000.
The resulting equations are solved numerically by em-
ploying a finite volume method based on similar nu-
merical flow simulation methods (Koop 2008; Kelleners
2007). This edge-based finite volume method uses an un-
structured grid and is suited for both two-dimensional and
three-dimensional domains, discretized using any combi-
nation 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 directed away from the
airfoil the following boundary conditions are applied:

IUd R>0: wall = 0 (9a)

where q is the variable of the considered governing equa-
tion; either Pd or Pd 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
1o0:

tdr < 0: < Vwall = ko + (Zwal ,o) *V40o (9b)

The resulting equations are solved for a multi-disperse
droplet distribution. For each droplet-bin in the distribu-
tion a separate set of equations is solved. The results from
each droplet-bin are combined in a single f, which is sub-
sequently 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:
Nbm,
J = 1.0 (10)
S
with Nbin the total number of droplet-bins, such that


Nb,.

fi-


Splashing Model


For traditional icing models (non SLD-regime), it is com-
mon to assume that every droplet that impinges on the
airfoil surface is fully deposited onto that surface, ei-
ther turning to ice, as runback-water or a combination
of ice and runback-water. However, several other deposi-
tion regimes are possible, such as splashing and rebound
of impinging droplets. One of the major differences be-
tween 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 air-
foil surface, while the remainder of the water will be
splashed away from the airfoil, possibly re-impinging
further downstream of the original point of impact. By
including a model for the mass-loss due to splashing,
a mass-loss-coefficient 0 < 0 < 1 can be determined,
which reduces the catching efficiency:

,Isplash = (1 0)f (12)


Models for droplet splashing often also predict the
number of secondary droplets and their velocity after


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


the splashing event (see Fig. 3), making it possible to
track the droplets further downstream, until possible re-
impingement.

N=3

d U d,, a d,






Figure 3: Variables related to a splash event

The basic model for mass-loss due to splashing from
Trujillo et al. (2000) causes a severe over-prediction of
the mass-loss coefficient. For improved performance, a
version of this basic model calibrated for mass-loss co-
efficients in SLD conditions by Honsek et al. (2008) is
used. The resulting mass-loss coefficient is expressed as:


(Ky) =3. (1 exp [-0.85 (Ky 17)1)


where Ky is the splashing parameter defined by Yarin and
Weiss (1995) as:

Ky =A-3/8(Oh-2/We)5/16

S3 LWC 5/16
S[Id(LWr /3-3 Oh-2/5We) 516 (14)

with the droplet Weber number We defined as:


Pa d2 d
We=-
o-d


and the Ohnesorge number Oh as:


Oh Pa
JPaOdd


Red
Red


In Eq. (13) the splashing threshold from Yarin and
Weiss appears; Ky,cri, = 17. Splashing only occurs for
Ky > 17 as illustrated in Fig. 4. Also note that, within 1%
accuracy, = 3 for Ky > 23.

For typical conditions (LWC z 0.5 10-3 kg/m3, Pa
1.225 kg/m3, Pa ` 1.789 10-5 Pa s, pw z 1000 kg/m3,
O-d 0.072 N/m, Ud ` 100 m/s, 0 m/s < a -Ud <
100 m/s) and diameters ranging from 10 Itm d <
1000 Im the following ranges for the relevant dimen-







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


From mass conservation and Eq. (17) the average size
of secondary droplets can be calculated:


- 1.0-
. 0.8-
-
0.6-
-
S0.4-
i 0.2-
o .o-


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


sionless numbers can be found:

We E [0, 1000]
Red e [0, 10000]
Oh e [0,0.1]
Ky e [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 average number of secondary droplets:


1 d I
N= 2 0.0437 K '' Kdry -44.92 (17)
22 Od t

where K is the splashing parameter defined by Cossali
et al. (1997), which is similar to Eq. (14):

K = Oh-- le (18)

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:


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


dPd,i IPd
d + V Pd,i Ud,i = -Mi
iBt


Trujillo et al. also determined relations for the average
secondary droplet velocities:


-ds t t d t-
Udr= 0.85 + 0.0025 arctan )
ffd-f [Udrt)


0.12 + 0.002 arctan Ud
Ud rT)


Ud,s n
IUd


(20a)


(20b)


dPd,i d,i Pa
at +V pd, di Ud,i = Pd,ifD,i+Pd,i 1 -- i-k
9t Pwl
(25)
where i = 1, 2..., Nbin indicates the considered droplet-
bin.


Results


In order to validate the Eulerian droplet tracking method
a suitable test case had to be selected. Because of the


3.8

- "


Ky,crit


For a Lagrangian method the information obtained
from Eqs. (13), (20) and (21) is sufficient to correct/p
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, corre-
sponding to the mass and momentum of the secondary
droplets, respectively.
From the local values ofpd and Od and from Eq. (13)
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 expe-
riences splashing:


M= J(pd df)> 0 : Splashing(22)
10: No splashing

The re-injected momentum then follows from the re-
injected mass rate (Eq. (22)) and Eq. (20), leading to
the following source-term, which again is only included
when splashing occurs in the considered finite-volume
element.

Sd .ft) Ud,,s MUd,s: Splashing
0 : No splashing

The calculation is performed using a multi-disperse
droplet distribution. The droplet-bin in which the flux-
terms from Eqs. (22) and (23) are included is chosen
such that the diameter is closest to the diameter of the
secondary droplets calculated in Eq. (21).
With the inclusion of terms M and I the following
equations are obtained:


0


ds =) d


Kc,dry = l h" '


I











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 pre-
diction code, the results from the Eulerian droplet track-
ing 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. Be-
sides 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. (2007) was se-
lected. They performed experiments with different droplet
median volume diameters (MVD) for a NACA-23012 air-
foil at 2.50 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.


Table 1: Conditions for selected cases
20 lm MVD 236 ~tm MVD
AoA 2.50
c 0.9144 m
LWC 0.19 g/m3 1.89 g/m3
Ua,1 78.23 m/s
T, 299 K
po 101330 Pa



A first assessment was made by performing a numer-
ical simulation for both the 20 tm 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 distribu-
tions 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.
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


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010



Table 2: 10-Bin droplet distributions for selected cases


Droplet
bin
1
2
3
4
5
6
7
8
9
10


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


Droplet size [im]
20 tm MVD 236 ~tm 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


is created for each vertex in the original triangular mesh.
The resulting catching efficiencies are shown in Figs. 7
and 8. For the 20 tm 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 diam-
eter 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 so-
lution is compared to the experimental result. However,
for the large MVD case of 236 Im, the improvements are
not so clear as for the small MVD case of 20 tm 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 Eule-
rian 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, es-
pecially around the nose of the airfoil (s z 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.
For the catching efficiencies calculated with the multi-
disperse Eulerian method shown in Figs. 7 and 8 cor-
responding 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 tm 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 re-
sulting ice shapes are shown in Fig. 9. The larger catching
efficiency of the 236 tm MVD case results in a larger ice













0.3-


0.2-


0.1- -


0.01
0.0 5 10 15 20 25 30 35 40 45

Droplet diameter, d [tm]
(a) 20 rtm MVD


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010



0.3-


0.2-


0.1-


0.0
0.0 200 400 600 800 1000
100 300 500 700 900 1100
Droplet diameter, d [ntm]
(b) 236 Ipm MVD


Figure 5: 10-Bin droplet distributions for selected cases


Dimensionless x-coordinate,

(a) Entire domain


-0.5




S0.0



-0
-0.5


0.0 0.5 1
Dimensionless x-coordinate,

(b) Close-up around the airfoil


Figure 6: Grid used for Eulerian computations


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


-0.2 -0.1 0.0
bottom nose top

Dimensionless airfoil-coordinate, s-s".
S
(a) Droplet distribution effects


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


-0.3 -0.2 -0.1 0.0
bottom nose top

Dimensionless airfoil-coordinate, s-s..
c
(b) Eulerian (Droplerian) and Lagrangian (2DFOIL) results


Figure 7: Calculated catching efficiencies, 20 tm MVD


44:







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


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


-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2
bottom nose -> top
Dimensionless airfoil-coordinate, s-s'..
c


1.0 Droplerian, 10 droplet-bins
S-. 2DFOIL-ICE, 10 droplet-bins
0.9 Experimental results


-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2
bottom nose -> top
Dimensionless airfoil-coordinate, s-ss.
c


(a) Droplet distribution effects (b) Eulerian (Droplerian) and Lagrangian (2DFOIL) results

Figure 8: Calculated catching efficiencies, 236 im MVD


accretion. The thickness of the ice accretions are similar
around the leading edge. The amount of energy imping-
ing 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 tm MVD case extends
further downstream of the leading edge, while the ice ac-
cretion for the 20 im MVD case is only present around
the leading edge.


- ice-shape, 20 vm MVD
- ice-shape, 236 vm MVD


0.04-

0.02-


-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, !


Figure 9: Calculated ice accretions

Splashing Effects The splashing model discussed in the


Table 3: Conditions for icing cases
AoA 2.50
c 0.9144 m
LWC 0.19 g/m3
t,,ac 78.23 m/s
To 268.15 K
po 101330 Pa




preceding 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 efficiencyf. 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 dra-
matic improvement is observed. Near the leading-edge
the under-prediction increased, but the match with exper-
imental results further downstream of the leading-edge
improved significantly.
For the icing shape calculated from the catching effi-
ciency 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 sim-
ilar to the icing shape determined without a splashing
model. The mass of ice that freezes near the leading-edge


0.08-











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 oc-
curs closer to the leading edge than for the case without
splashing.
Since the Eulerian model provides information of the
full droplet velocity and droplet density 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 [tm) is shown in Fig. 12. Droplets from the bins
with larger droplet diameters splash, creating small sec-
ondary droplets with a diameter close to 16 Itm. 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.
Figures 11-13 provide an image of the amount of water
in a single droplet bin only. The total (cumulative) amount
of water is shown in Fig. 14. Again a comparison is shown
between the cases with and without splashing model. A
small change can be seen in the region around the leading
edge caused by the re-injection of the secondary droplets
in the smaller bins. It can be noted that, since most of
these droplets are convected downstream of the airfoil,
this explains the decrease in catching efficiency noted in
Fig. 10(a).


Forces on Droplets

One of the most important parts of the ice accretion sim-
ulation method (for normal as well as for and SLD con-
ditions) is the calculation of droplet trajectories. These
trajectories determine the calculated amount of water on
the airfoil section and therefore the icing shape. Since
these trajectories are at the base of the method they are
focal point of this study. In the icing field it is common
for the trajectory calculations to only include the most
basic of the forces acting on the droplet: drag, gravity and
the buoyancy caused by the gravitational field. In order
to check the validity of the assumption that these are the
only significant forces for the droplet trajectory calcu-
lations a number of calculations have been performed
including other forces acting on the droplets. The analy-
sis of the equations for these forces is further explained
in van Eijkeren and Hoeijmakers (2010). The considered
forces are:


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Drag

Gravity

Buoyancy

Basset history force

Virtual mass

The remaining force terms are induced by rotation in the
flow which means that, for the potential-flow model used
in the current method, they will amount to zero.
Drag Force The drag force is accounted for using the
formulation described for the Eulerian model in Eq. 6
and Eq. 8.
Gravity Force In most common models the gravity is
combined with at least part of the buoyancy. For this study
the two have been separated. The gravity force consists
only of the force acting on the droplet by gravitational
acceleration:
f = pwVd (26)
Buoyancy Force The buoyancy force is caused by a
pressure gradient. In most common methods, only the
pressure gradient caused by the constant gravity field is
taken into account, resulting in:


f = -PaVdg


However, for this analyses the pressure gradient also takes
into account the local pressure gradient induced by the
flow field.
f = -VdVp (28)
Basset History Force Perhaps the most important ques-
tion here is whether or not the Basset history force has
to be accounted for. This term is based on the relative
speed at which the boundary layer adapts to changes in
the surrounding flow. Calculation of this term is particu-
larly complicated because it involves a time-integration
over the path of the droplet.


A t d( 0d)
fB = -3d,7r K(t r, r) dt dr
-co


where the kernel K has been chosen as the kernel for
non-creeping flow conditions by Mei et al. (1991):


4K (( r,) = 47 (t(t -r)2 3l -25
Td f Td
(30)
with rd = P- and fH = 0.75 + 0.2Red
Virtual Mass Force The virtual mass force is based
on the acceleration of the air surrounding the droplets.
Accelerating a droplet means that the air surrounding it








7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


- Droplerian, splashing
- Droplerian, no splashing
SExperimental results


0.08-- ice-shape, splashing
ice-shape, no splashing
0.06-


0.04-

0.02-


-0.02

-0.04-


-0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2
bottom +-- nose top

Dimensionless airfoil-coordinate, s-s'..
S
(a) Catching effiencies


-0.06--
-0.02


0.0 0.02 0.04 0.06 0.08 0.10 0.12
Dimensionless x-coordinate, 1
(b) Ice accretions
(b) Ice accretions


Figure 10: Results 236 tm MVD with splashing effects


(a) Without splashing model


(b) With splashing model


Figure 11: Calculated droplet distributions [kg/m3], 236 tm MVD, bin 10 (1046 pm)






7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


(a) Without splashing model








(b) With splashing model
Figure 12: Calculated droplet distributions [kg/m3], 236 wm MVD, bin 1 (16 pm)


(a) Without splashing model (b) With splashing model
Figure 13: Calculated droplet distributions [kg/m3], region near the leading edge, 236 tm MVD, bin 1 (16 Im)


W
0. OM 91,
0.000128
8.57e-05
4.28e-05
0.00


W
0. OM 91,
0.000128
8.57e-05
4.28e-05
0.00



























(a) Without splashing model


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


















(b) With splashing model


Figure 14: Calculated droplet distributions [kg/m3], region near the leading edge, 236 ~m MVD, Cumulative


has to accelerate as well. The virtual mass force therefore
depends on the relative acceleration:


2 ( dt dt )
m = ad(-T dg d (31)

Force Analysis A force analysis has been performed
using the potential flow-field and the droplet bins for
the 236 ~tm MVD case. The trajectories of the droplets
are calculated using all the previously mentioned forces.
The analyses reveals that when the droplets are within
one chord-length of the airfoil the dominant force is the
drag force. For the largest and smallest droplet bin this
is shown in Fig. 15, only the droplet trajectories corre-
sponding to the maximum of each force term are shown
(the extreme cases). When comparing the influence of
the individual forces the drag force, the buoyancy force
and the Basset history force appear to have the largest
influence. The virtual mass term appears to have virtually
no influence, its ratio is one thousandth of the drag force.
The most important part of the droplet trajectories is
the part close to the airfoil, this determines whether a
droplet hits or misses the airfoil surface. Examining the
magnitude of the forces the drag is the only relevant
force when looking at the smallest droplet size. How-
ever, for the largest droplet size the buoyancy force be-
comes increasingly important. The Basset force also be-
comes more important, but it is an order of magnitude
smaller than the drag and the buoyancy forces. Further-
more, when looking at the direction of the individual
forces, the Basset force works in the same direction as
the drag, meaning that the direction of the droplet is not
changed. This validates the generally accepted approach
in which the Basset history force is not taken into ac-
count. Furthermore, considering the actual trajectories
the difference between the trajectories calculated using


the traditional approach and the trajectories calculated
using all the forces is very small. The effect on the dis-
tribution of the amount of water is therefore minimal
and not justifying the extra effort needed to calculate the
Basset history force.


Conclusions

An Eulerian method to calculate ice-accretions on two
dimensional airfoils has been developed. The method
was based on and compared with a similar Lagrangian
method, which produced satisfactory results for small
droplet diameters. For larger (SLD) droplets, compared
with experimental data, both methods over-predict the
catching efficiency. To improve the matching with ex-
perimental catching efficiency data, a splashing model
has been added to the Eulerian method. Splashing is ac-
counted for by introducing a mass-loss-coefficient and by
re-injecting secondary droplets into the droplet bin corre-
sponding to the diameter closest to the secondary droplet
diameter obtained from the splashing model. The inclu-
sion of a splashing model accounts for a decrease in the
catching efficiency on the airfoil, bringing the catching ef-
ficiency prediction closer to experimental results. This is
a clear improvement compared to the case where splash-
ing was not accounted for. 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 catch-
ing 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 air-
foil, by using more advanced flow-models instead of the
currently applied potential-flow-model, might also im-







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


le-06


le-07


le-08


le-09


le-10


le-11


0.001 i i i i
Drag -
Virtual Mass -
Buoyancy -



le-05


le-06


le-07


le-08


le-09
0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
(distance from impact)/c
(a) Largest droplet bin (1046 prm)


0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05
(distance from impact)/c
(b) Smallest droplet bin (16 prm)


Figure 15: Calculated maximum forces on droplets


prove the results.
A detailed analysis of the droplet forces has also been
performed. The assumption that the drag and buoyancy
forces are the dominant forces appears to be valid, al-
though for larger droplets the Basset history force be-
comes more important. The effect of the inclusion of
additional force terms in terms of the droplet trajectories
is minimal while the impact on calculation/computational
time is very large.



Acknowledgements

This research was made possible by the European Union
sponsored project EXTICE. This project focuses on in-
flight SLD ice accretion on aircraft. International part-
ners include universities, aircraft manufacturers and re-
search institutes. The project includes: droplet impact
experiments and icing experiments, leading to improved
SLD-specific models and two and three dimensional ic-
ing results that can be used for validation. Numerical
simulations with these models and validation cases will
be performed.


References

G. E. Cossali, A. Coghe, and M. Marengo. The impact
of a single drop on a wetted solid surface. Experiments
in Fluids, 22(6):463-472, April 1997.


J. E. Dillingh and H. W. M. Hoeijmakers. Numeri-
cal simulation of airfoil ice accretion and thermal anti-
icing systems. In ICAS 2004, 24th international congress
of the aeronautical sciences, Yokohama, Japan, August-
September 2004.

J. E. Dillingh and H. W. M. Hoeijmakers. Simulation
of ice accretion on airfoils during flight. In FAA In-
Flight Icing' Ground De-Icing International Conference
& Exhibition, Chicago, IL, June 2003. 12 pages.

R. Honsek, W. G. Habashi, and M. S. Aub6. Eule-
rian modeling of in-flight icing due to supercooled large
droplets. Journal of Aircraft, 45(4):1290-1296, July-
August 2008.

S. J. Jacobs, J. M. Hospers, and H. W. M. Hoeijmak-
ers. Numerical simulation of ice accretion on multiple-
element airfoil sections. In ICAS 2008, 26th international
congress of the aeronautical sciences, Anchorage, AK,
September 2008.

P. H. Kelleners. An edge-based finite volume method
for inviscid compressible flow with condensation. PhD
thesis, University of Twente, Enschede, the Netherlands,
December 2007.

A. H. Koop. Numerical simulation of unsteady three-
dimensional sheet cavitation. PhD thesis, University of
Twente, Enschede, the Netherlands, September 2008.


S Drag -
Virtual Mass -
Buoyancy
Basset-





~-~-


~-----







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


I. Langmuir and K. B. Blodgett. A mathematical inves-
tigation of water droplet trajectories. Technical Report
5418, US Army Air Forces, 1946.

R. Mei, C. J. Lawrence, and R. J. Adrian. Unsteady
drag on a sphere at finite reynolds-number with small
fluctuations in the free-stream velocity. Journal of Fluid
Mechanics, 233:613-631, 1991.

Michael Papadakis, Arief Rachman, See-Cheuk Wong,
Hsiung-Wei Yeong, Kuohsing E. Hung, Giao T. Vu, and
Colin S. Bidwell. Water droplet impingement on sim-
ulated glaze, mixed and rime ice accretions. Technical
Report NASA/TM-2007-213961, NASA, October 2007.

K. R. Petty and C. D. J. Floyd. A statistical review
of aviation airframe icing accidents in the us. In llth
Conference on Aviation, Range, and Aerospace, Hyannis,
MA, October 2004.

M. Snellen, O. J. Boelens, and H. W. M. Hoeijmakers. A
computational method for numerically simulating ice ac-
cretion. In 15th AIAA Applied Aerodynamics Conference,
Atlanta, GA, June 1997. Technical Papers pt 2.

M. F. Trujillo, W. S. Mathews, C. F. Lee, and J. F. Peters.
Modelling and experiment of impingement and atomiza-
tion of a liquid spray on a wall. International Journal of
Engine Research, 1(1):87-105, 2000. ISSN 1468-0874.

D. F. van Eijkeren and H. W. M. Hoeijmakers. Influence
of the history term in a lagrangian method for oil-water
separation. In 7th International Conference on Multi-
phase Flow, Tampa, FL, 2010.

W.B. Wright and M.G. Potapczuk. Semi-empirical mod-
elling of sld physics. In 42nd AIAA Aerospace Sciences
Meeting and Exhibit, Reno, NV, January 2004.

A. L. Yarin and D. A. Weiss. Impact of drops on solid
surfaces: self-similar capilary waves, and splashing as a
new type of kinematic discontinuity. Journal of Fluid
Mechanics, 283:141-173, January 1995.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs