7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Eulerian Method for Ice Accretion on MultipleElement 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 multipleelement 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 2DFOILICE 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 multidisperse 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
2DFOILICE 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 Crosssectional 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 Nondimensional 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 Masslosscoefficient []
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 19922000, 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 nonaerodynamic 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 leadingedge
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 twodimensional flow. Studies to
extend the twodimensional 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 multielement airfoil sections. In this method
the droplet trajectories are calculated using a Lagrangian
method. Ice accretion modeling is based on the Messinger
model, a quasisteady model that takes into account all
important mass and heat transfer processes that occur
when supercooled 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 SLDregime. Due to the size of droplets in the SLD
regime (d > 50 Im) versus droplets in the normalregime
(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 2DFOILICE. This development is
motivated by a Lagrangian computation, using several
droplet bins for multipleelement 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 monodisperse to multidisperse 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 2DFOILICE.
Eulerian Droplet Tracking
The 2DFOILICE method is based on Lagrangian droplet
tracking, in which the (potential) flow field is calculated
using a panelmethod. This has some limitations, espe
cially for the multipleelement airfoil cases. Due to the
potentialflow model in 2DFOILICE, the geometry of
multipleelement 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 2DFOILICE, but also Euler or full
NavierStokes). 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 threedimensional 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
=d7 (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 sourceterm 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 Ud 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* 104Re138 (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 edgebased finite volume method uses an un
structured grid and is suited for both twodimensional and
threedimensional domains, discretized using any combi
nation of element types. The flow field of the surrounding
air is used only as input, since oneway 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 multidisperse
droplet distribution. For each dropletbin in the distribu
tion a separate set of equations is solved. The results from
each dropletbin 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 dropletbins has to
equal one:
Nbm,
J = 1.0 (10)
S
with Nbin the total number of dropletbins, such that
Nb,.
fi
Splashing Model
For traditional icing models (non SLDregime), 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 runbackwater or a combination
of ice and runbackwater. 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 reimpinging
further downstream of the original point of impact. By
including a model for the massloss due to splashing,
a masslosscoefficient 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 massloss due to splashing from
Trujillo et al. (2000) causes a severe overprediction of
the massloss coefficient. For improved performance, a
version of this basic model calibrated for massloss co
efficients in SLD conditions by Honsek et al. (2008) is
used. The resulting massloss 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 =A3/8(Oh2/We)5/16
S3 LWC 5/16
S[Id(LWr /33 Oh2/5We) 516 (14)
with the droplet Weber number We defined as:
Pa d2 d
We=
od
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 103 kg/m3, Pa
1.225 kg/m3, Pa ` 1.789 105 Pa s, pw z 1000 kg/m3,
Od 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: Massloss 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 nondimensional
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 )
ffdf [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  ik
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 reinjected
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 reinjected mass rate can be computed. This
leads to the following term, which is only included when
the considered finitevolume is next to the wall and expe
riences splashing:
M= J(pd df)> 0 : Splashing(22)
10: No splashing
The reinjected momentum then follows from the re
injected mass rate (Eq. (22)) and Eq. (20), leading to
the following sourceterm, which again is only included
when splashing occurs in the considered finitevolume
element.
Sd .ft) Ud,,s MUd,s: Splashing
0 : No splashing
The calculation is performed using a multidisperse
droplet distribution. The dropletbin 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 2DFOILICE, which calculates a potential
flow field using a panelmethod in order to carry out
Lagrangian droplet tracking, a validated iceaccretion 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
potentialflow field as calculated by 2DFOILICE. 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 NACA23012 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 multidisperse 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 2DFOILICE and the Eulerian
Droplerian, with a monodisperse 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 dropletbins.
For the following simulations 10 dropletbins are used.
The 10bin 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 surfaceelements 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: 10Bin 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
multidisperse solution (Fig. 7(b) and Fig. 8(b)). This is
due to the inclusion of dropletbins with a smaller diam
eter than the MVD, leading to droplet trajectories that
better follow the streamlines. The inclusion of multiple
dropletbins 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
overprediction 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 timestep 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: 10Bin droplet distributions for selected cases
Dimensionless xcoordinate,
(a) Entire domain
0.5
S0.0
0
0.5
0.0 0.5 1
Dimensionless xcoordinate,
(b) Closeup around the airfoil
Figure 6: Grid used for Eulerian computations
 Droplerian, 10 dropletbins
 Droplerian, 1 dropletbin
 Experimental results
0.2 0.1 0.0
bottom nose top
Dimensionless airfoilcoordinate, ss".
S
(a) Droplet distribution effects
0.6 Droplerian, 10 dropletbins
2DFOILICE, 10 dropletbins
. Experimental results
0.3 0.2 0.1 0.0
bottom nose top
Dimensionless airfoilcoordinate, ss..
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 dropletbins
 Droplerian, 1 dropletbin
 Experimental results
0.5 0.4 0.3 0.2 0.1 0.0 0.1 0.2
bottom nose > top
Dimensionless airfoilcoordinate, ss'..
c
1.0 Droplerian, 10 dropletbins
S. 2DFOILICE, 10 dropletbins
0.9 Experimental results
0.5 0.4 0.3 0.2 0.1 0.0 0.1 0.2
bottom nose > top
Dimensionless airfoilcoordinate, sss.
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.
 iceshape, 20 vm MVD
 iceshape, 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 xcoordinate, !
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 2DFOILICE
and Droplerian. Implementing reinjection of droplets
is quite complicated for the Lagrangian method, due to
the calculation of the catching efficiencyf. Reinjection
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 leadingedge
the underprediction increased, but the match with exper
imental results further downstream of the leadingedge
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 leadingedge is very sim
ilar to the icing shape determined without a splashing
model. The mass of ice that freezes near the leadingedge
0.08
is similar, however, the amount of runback 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 1113 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 reinjection 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 potentialflow 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 timeintegration
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
noncreeping 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 iceshape, splashing
iceshape, 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 airfoilcoordinate, ss'..
S
(a) Catching effiencies
0.06
0.02
0.0 0.02 0.04 0.06 0.08 0.10 0.12
Dimensionless xcoordinate, 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.57e05
4.28e05
0.00
W
0. OM 91,
0.000128
8.57e05
4.28e05
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 flowfield 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 chordlength 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 iceaccretions 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 overpredict 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 masslosscoefficient and by
reinjecting 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 iceaccretion 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 flowmodels instead of the
currently applied potentialflowmodel, might also im
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
le06
le07
le08
le09
le10
le11
0.001 i i i i
Drag 
Virtual Mass 
Buoyancy 
le05
le06
le07
le08
le09
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
SLDspecific 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):463472, 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 DeIcing International Conference
& Exhibition, Chicago, IL, June 2003. 12 pages.
R. Honsek, W. G. Habashi, and M. S. Aub6. Eule
rian modeling of inflight icing due to supercooled large
droplets. Journal of Aircraft, 45(4):12901296, 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 edgebased 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 reynoldsnumber with small
fluctuations in the freestream velocity. Journal of Fluid
Mechanics, 233:613631, 1991.
Michael Papadakis, Arief Rachman, SeeCheuk Wong,
HsiungWei 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/TM2007213961, 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):87105, 2000. ISSN 14680874.
D. F. van Eijkeren and H. W. M. Hoeijmakers. Influence
of the history term in a lagrangian method for oilwater
separation. In 7th International Conference on Multi
phase Flow, Tampa, FL, 2010.
W.B. Wright and M.G. Potapczuk. Semiempirical 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: selfsimilar capilary waves, and splashing as a
new type of kinematic discontinuity. Journal of Fluid
Mechanics, 283:141173, January 1995.
