Ice accretion simulation on multi-element airfoils using extended Messinger model

18 Pages·2008·0.85 MB·English

HeatMassTransfer(2009)45:305–322
DOI10.1007/s00231-008-0430-4
ORIGINAL
Ice accretion simulation on multi-element airfoils using extended
Messinger model
S. O¨zgen Æ M. Canıbek
Received:25October2007/Accepted:22July2008/Publishedonline:8August2008
(cid:1)Springer-Verlag2008
Abstract In the current article, the problem of in-flight to larger elements. The results also indicate that the multi-
ice accumulation on multi-element airfoils is studied layer approach yields more accurate results compared to
numerically. The analysis starts with flow field computa- theone-layerapproach,especiallyforglazeiceconditions.
tionusingtheHess-Smithpanelmethod.Thesecondstepis
thecalculationofdroplettrajectoriesanddropletcollection
efficiencies. In the next step, convective heat transfer
coefficient distributions around the airfoil elements are 1 Introduction
calculatedusingtheIntegralBoundary-LayerMethod.The
formulation accounts for the surface roughness due to ice Ice accumulation on parts of the airframe is one of the
accretion. The fourth step consists of establishing the fundamental problems of aviation. Ice growth on wings,
thermodynamic balance and computing ice accretion rates tail surfaces, fuselage and other items like the engine
usingtheExtendedMessingerModel.Atlowtemperatures intakes and pitot tubes result in severe performance deg-
and low liquid water contents, rime ice occurs for which radation, thus threatening flight safety. For example,
the ice shape is determined by a simple mass balance. At modification of the wing shape due to accumulated ice
warmer temperaturesand high liquid water contents,glaze results in reduced lift together with increased drag and
ice forms for which the energy and mass conservation weight.Iceformationoncontrolsurfacesresultsinserious
equationsarecombinedtoyieldasinglefirstorderordinary and often unpredictable degradations in the controllability
differential equation, solved numerically. Predicted ice of aircraft. If an airplane is to fly in icing conditions, it
shapes are compared with experimental shapes reported in must demonstrate that it can operate safely in conditions
theliteratureandgoodagreementisobservedbothforrime prescribed by Certification Authorities, like those defined
andglazeice.Iceshapesandmassesarealsocomputedfor in Federal Aviation Regulations, Part 25, Sect. 25.1419.
realistic flight scenarios. The results indicate that the Certification process may involve flight and/or laboratory
smaller elements in multielement configurations accumu- testing and numerical simulation. Numerical ice accretion
late comparable and often greater amount of ice compared simulationreduces (but never totallyreplaces)the demand
for flight and laboratory testing.
Efforts towards understanding the effects of ice on
S.O¨zgen(&)(cid:1)M.Canıbek
performance and flight mechanics started in the 1940s.
TurkishAerospaceIndustries,FlightSciencesDepartment,
These were mainly based on experiments and in-flight
MiddleEastTechnicalUniversityTechnopolis,
06531Ankara,Turkey testing. Among the pioneering works, the published work
e-mail:[email protected];[email protected] of Messinger [6] represents an important foundation and a
M.Canıbek milestone in numerical ice accretion simulation. With the
e-mail:[email protected] advent of digital computers in the 1970s, theoretical
research was directed towards representative geometries
S.O¨zgen
such as airfoils, wings and helicopter rotor blades. The
DepartmentofAerospaceEngineering,
MiddleEastTechnicalUniversity,06531Ankara,Turkey major contributors to the aircraft icing simulations are
123
306 HeatMassTransfer(2009)45:305–322
NASA Lewis Research Center (USA), Defence Research collection efficiency and ice accretion calculations. In this
Agency (DRA-UK), Office National d’Etudes et des context, a computer code is developed in FORTRAN
Recherches Ae´rospatiales (ONERA-France), Anti-Icing programming language. Inputs to the problem are the
Materials International Laboratory (AMIL-Canada) and ambient temperature T , freestream velocity V , liquid
a ?
Italian Aerospace Research Center (CIRA-Italy), each water content (LWC) of air q , droplet median volume
a
having developed an ice accretion simulation code. diameter d (MVD), total icing time t , angle of attack a
p exp
Cebeci et al. [1] describe a numerical method for com- and the airfoil geometry. The liquid water content (LWC)
puting ice shapes on airfoils and their effects on lift and istheweightofliquidwaterpresentinaunitvolumeofair,
drag coefficients. The Interactive Boundary Layer Method often expressed in g/m3. The droplet median volume
developed by Cebeci has been incorporated into the diameter(MVD),whichisoftengiveninlm(microns)isa
LEWICE code of NASA to improve the accuracy of ice term usedtodescribethedropletsize. Itis thedroplet size
shape predictions and to compute the performance char- at which one-half of the given volume consists of larger
acteristics of airfoils. droplets and one-half consists of smaller droplets.
A NASA Report [13] summarizes the results of a ten- The solution starts with the calculation of the pressure
year collaborative research on ice accretion simulation distribution around the given airfoil shapes using a panel
between NASA, DRA and ONERA. The report includes method.Thesamecalculationservestodeterminingairand
the descriptions of the codes developed by these institu- droplet velocities anywhere in the flow field. Droplets are
tions and the results obtained. The report also presents ‘‘fired’’ at a plane far upstream (10 chords) and their tra-
comparisonswithiceshapesobtainedexperimentallyinthe jectories are calculated by integrating the equations of
NASA Lewis Icing Research Tunnel with a 21’’ chord motion in differential form in two dimensions. Distance
NACA 0012 airfoil. between two adjacent particles is taken as 10-4 m in the
Mingione and Brandi [7] present results on ice shape present calculations. Impact locations and droplet collec-
simulation over multi-element airfoils. They describe and tion efficiencies on the airfoil surface are thus found.
compare different ways to solve the transient ice accretion Predictionoftheheattransfercoefficientdistributionplays
problem, i.e., single-step, multi-step and predictor-correc- animportantroleinicingpredictions.Characteristicsofthe
tor methods. viscous flow such as the skin friction coefficient distribu-
Inareviewpaper,Gentetal.[3]presentthebackground tion over the airfoil surface are determined by using
andthestatusofanalysesaddressingaircrafticingproblem. empirical relations either based on experimental results or
Methods for water droplet trajectory calculation, ice solutions of the Integral Boundary–Layer Equation. From
accretion prediction and aerodynamic performance degra- theobtaineddata,itispossibletocalculateconvectiveheat
dation are discussed and recommendations for further transfer coefficients using Reynolds’ Analogy. Depending
research are made. on parameters like the freestream velocity, ambient tem-
Myers [8] presents a one-dimensional mathematical perature, liquid water content and collection efficiency,
model, extending the original Messinger Model describing rime ice, glaze ice or mixed ice accumulates on the airfoil
ice growth. It is demonstrated that the model can also be surface. Rime ice prediction involves a simple mass bal-
extended to two and three-dimensions. A modified version ance since droplets freeze immediately on impact. Under
of the two-dimensional extension proposed by Myers is milder conditions glaze ice develops, involving a layer of
employed in the current study. water lying on top of a layer of ice. There is evidence that
Myersetal.[9]discussamathematicalmodelforwater glazeiceisalwaysprecededbyathinlayerofrimeiceand
flow in glaze ice conditions. Water flow can significantly the transition from rime to glaze ice is smooth, i.e.,
complicate the problem and can have a major impact on freezing fraction reduces smoothly from unity to the
finaliceshapes.Ithasbeenpointedoutthatpreviouscodes equilibrium value, which is less than unity. In the original
cannot deal adequately with this issue. The model is Messingermethod,thistransitionissudden,resultinginan
applied to ice accretion problem and results are presented underprediction of the ice thickness. Therefore, the
for ice growth and water flow driven by gravity, surface Extended Messinger Model used in this study reflects the
tension and constant air shear. physics of the problem better. Under gravitational and/or
Fortin et al. [2] propose an improved roughness model, aerodynamic forces, water layer may flow downstream
inwhichthewaterstateonthesurfaceisrepresentedinthe (called runback water) or may be shed.
form of beads, film or rivulets. The model is tested for The calculations may bedoneeither inone-layer mode,
severe icing conditions at six different temperatures cor- wheretheiceshapesarepredictedinonestepfortheentire
responding to dry, mixed and wet ice accretion. duration of t or in multi-layer mode where t is divi-
exp, exp
Present study is an effort to predict ice shapes combin- ded into segments (or layers). In the multi-layer mode,
ingestablishedapproachesforflowfield,droplettrajectory, flowfield, droplet trajectory and ice calculations are
123
HeatMassTransfer(2009)45:305–322 307
repeated for each layer. This approach allows the effect of
ice shapes on flowfield and droplet trajectories to be taken
into account, thus reflecting the physics of the problem
more realistically. It also allows cases with varying ambi-
ent and icing conditions to be treated, like climbing flight
whichisanoveltyofthecurrentstudy.Thisfeatureallows
icing computations to be performed for the entire flight
profile of an airplane. Another important aspect of the
present study is that, although multi-element airfoils have
been treated by other researchers in the past, an extensive
parametric study like the one presented in this study does
not exist to the authors’ best knowledge. Therefore, the
present study fills an important gap in the literature.
Section 2 describes the solution method, briefly
explaining the flowfield, droplet trajectory, droplet collec-
tion efficiency and convective heat transfer coefficient
calculations,andtheExtendedMessingerModel.Section 3
is devoted to code validation, where the ice shapes
obtainedinthecurrentstudyarecomparedtoexperimental
and numerical ice shapes reported in the literature. In
Section 4, several realistic icing scenarios are studied and
the resulting ice shapes for single and two-element con-
figurations are presented. Finally, Section 5 summarizes
the study and points out important conclusions.
2 Problem formulation and solution method
In this section, the method developed for ice accretion
calculations is summarized. A brief flowchart of the cal-
culationprocedureandthedevelopedprogramispresented
in Fig. 1.
2.1 Flow field solution: panel method
In order to determine the pressure distribution around the
airfoil and provide the air velocities required for droplet
trajectory calculations, a panel method [5] is employed in
this study. In this method, the geometry is discretized by
quadrilateral panels each being associated with a singu-
larity element of unknown but constant strength. The
developed code uses N panels to solve for N singularity
strengthsusingtheflowtangencyboundaryconditionatthe
surfaces. A velocity potential can then be constructed for
any point in the flow field using the calculated singularity
strengths. The velocity components at the given point are
the x-, y- and z-derivatives of this velocity potential. The
velocity distribution around the airfoil is also used for
boundary-layer calculations in order to determine the
convective heat transfer coefficients. Results of the devel-
oped code are compared with experimental data [13] in
Fig. 2forsingleandtwo-elementcases.Althoughthereisa Fig.1 Flowchartofthepresentcalculationprocedure
123
308 HeatMassTransfer(2009)45:305–322
a 3 • Droplet sizes are small, hence remain spherical,
PresentStudy • The flowfield is not affected by the presence of the
Experimental droplets,
• Gravity and aerodynamic drag are the only forces
2 involved.
These assumptions are valid for d B 500 lm. These
p
assumptions are safe, as droplet sizes of 25 lm or greater
Cp 1 are found in only 4% of encounters [4]. The governing
equations for droplet motion are:
mx€ ¼(cid:2)Dcosc; ð1Þ
p
0 my€ ¼(cid:2)Dsincþmg; ð2Þ
p
y_ (cid:2)V
c¼tan(cid:2)1 p y; ð3Þ
x_ (cid:2)V
p x
-10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 D¼1=2qV2CDAp; ð4Þ
x/c
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
(cid:2) (cid:3)2 (cid:2) (cid:3)2
b 3 V ¼ x_p(cid:2)Vx þ y_p(cid:2)Vy : ð5Þ
Presentstudy
Experimental Intheaboveexpressions,CDisthedropletdragcoefficient,
V and V are the components of the flow field velocity at
x y
the droplet location and x_ ; y_ ; x€ ; y€ are the components
2 p p p p
of the droplet velocity and acceleration. Atmospheric
density and droplet cross-sectional area are denoted by q
and A .
p
Cp 1 Thedragcoefficientsofthedropletsarecalculatedusing
the following drag law [3]:
C ¼1þ0:197Re0:63þ2:6(cid:3)10(cid:2)4Re1:38; Re(cid:4)3500;
D
0
C ¼ð1:699(cid:3)10(cid:2)5ÞRe1:92; Re[3500; ð6Þ
D
where, Re = qVd /l is the Reynolds number based on
p
droplet diameter d and relative velocity V, while l is the
p
-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 atmospheric viscosity. This parameter is calculated using
x/c Sutherland’s viscosity law as a function of ambient
temperature [12]. The pattern of droplet impact on the
Fig.2 Comparison of pressure distributions obtained from panel
method and experiment (NACA 23012 section for airfoil and flap). airfoil determines the amount of water that impacts the
aSingleelementairfoil,a=-0.27(cid:2);bairfoilwithanexternalflap, surface and the region subject to ice growth. The local
a=-1.05(cid:2)
collection efficiency is defined as the ratio of the area of
impingement to the area through which water passes at
some distance upstream of the airfoil. The local collection
slight disagreement of the pressure coefficients especially efficiency can be defined as:
ontheflap,theresultsareaccurateenoughforthepurposes
dy Dy
of this study. The mentioned disagreement can be attrib- b¼ o (cid:5) o; ð7Þ
ds Ds
uted tothe fact thattheflowmay beseparatedclose tothe
wheredy isthedistancebetweentwowaterdropletsatthe
trailing edge of the flap, a phenomenon that cannot be o
release plane and ds is the distance between the impact
modeled by the panel method used in this study.
locationsofthesametwodropletsontheairfoil,seeFig. 3.
ExamplesofparticletrajectoriesarepresentedinFig. 4for
2.2 Droplet trajectories and the collection efficiencies single and two-element configurations. Note that, for the
two-element configuration, the icing problem is influenced
Thefollowingassumptionsaremadefortheformulationof bymuchlargernumberofdropletscomparedtothesingle-
the equations of motion for the water droplets: element case.
123
HeatMassTransfer(2009)45:305–322 309
0.5 impact it compared to a case of a larger body. Greater
droplet size and airspeed increases the collection effi-
ciency,asbothoftheseparametersincreasedropletinertia.
0.25
Themoreinertia adroplethas, themoredifficultitwillbe
to deviate it from the body. On the other hand, angle of
m) 0
y ( attack determines the region of ice accumulation. The
highest collection efficiency occurs at the lowest angle of
∆s
-0.25 attack(inmagnitude)depicted,althoughthecorresponding
∆yo
region subject to icing is the lowest.
In the developed tool, droplet trajectory calculations
-0.5
-1 -0.5 0 0.5 1 consumemorethan90%oftheCPUtime.Inamulti-layer
x (m)
calculation,asthedropletcalculationsarerepeatedateach
Fig.3 Definitionofcollectionefficiency layer, total CPU time is proportional to the number of
layers. In this study, it is found that the best compromise
between computational time and accuracy is a four-layer
a1.2 approachregardlessoftheexposuretime,ashighernumber
of layers does not improve the accuracy in a significant
0.8
manner.
0.4
(m) 0 2.3 Calculation ofconvective heat transfer coefficients
y
-0.4
The current study employs an Integral Boundary Layer
-0.8 Method for the calculation of the convective heat transfer
coefficients.Thismethod enablescalculationofthe details
-1.2-2 -1 0 1 2 of the laminar and turbulent boundary layers fairly accu-
x(m) rately. Transition prediction is based on the roughness
b1.2 Reynolds number, Rek = qUkks/l, where ks is the rough-
ness height and U is the local airflow velocity at the
k
0.8 roughness height calculated from the following expression
[11]:
0.4
y(m) 0 Uk ¼2ks(cid:2)2(cid:4)ks(cid:5)3þ(cid:4)ks(cid:5)4þ1d2dUeks(cid:4)1(cid:2)ks(cid:5)3: ð8Þ
U d d d 6m ds d d
e a
-0.4
Intheaboveexpression,U istheflowvelocityoutsidethe
e
-0.8
boundary-layer at the roughness location and s is the
-1.2 streamwisedistancealongtheairfoilsurfacestartingatthe
-2 -1 0 1 2
stagnationpoint.Roughnessheightiscalculatedfromk ¼
x (m) s
ð4r l =q FsÞ1=3 [13], where r , q and l are the
w w w w w w
Fig.4 ParticletrajectoriesforaNACA4412airfoil,V?=92.6m/s, surface tension, density and viscosity of water,
a=4o.aone-element;btwo-element
respectively. Fraction of the airfoil surface that is wetted
by water droplets is denoted by F, while s denotes local
surface shear stress. The boundarylayerthicknessis given
Effects of key parameters like chord, droplet size, air- by [12]:
speedandangleofattackonthecollectionefficiencybare
illustrated in Fig. 5. Results in Fig. 5 suggest that smaller 315
d¼ h: ð9Þ
l
airframes, lower speed and larger drop sizes increase the 37
collection efficiency, hence possibility of ice formation.
Laminar momentum thickness is computed using
Smaller size increases collection efficiency as the body
Thwaites’ formulation [12]:
createsasmallerobstaclefortheincomingdropletsandthe
s
deviation of the droplets away from the body is not suffi- h2 0:45Z
cientforthemtoavoidit.Inotherwords,whenadropletis ml ¼ U6 Ue5ds: ð10Þ
e
heading for a small body, it is more likely that it will 0
123
310 HeatMassTransfer(2009)45:305–322
a 0.8 b0.8
c=0.5m dp=10µm
0.7 cc==12..00mm 0.7 ddpp==2400µµmm
0.6 0.6
0.5 0.5
β 0.4 β 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
-0.16 -0.12 -0.08 -0.04 0 0.04 -0.2 -0.16 -0.12 -0.08 -0.04 0 0.04
s/c s/c
c 0.8 d0.8
V =50m/s α= 0o
0.7 VViinnff==120000mm//ss 0.7 αα==−46oo
inf
0.6 0.6
0.5 0.5
β 0.4 β 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
-0.16 -0.12 -0.08 -0.04 0 0.04 -0.12 -0.08 -0.04 0 0.04 0.08 0.12
s/c s/c
Fig.5 Effectofvariousparametersonthedropletcollectionefficiency(NACA0012airfoil).aEffectofchord;beffectofdropletsize;ceffect
offreestreamvelocity;deffectofangleofattack
For laminar flow Rek B 600, the equation of Smith and St¼ Cf=2 ; ð13Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Spalding is adopted to calculate the convective heat Pr þ ðC =2Þ=St
t f k
transfer coefficient [3]:
where Pr = 0.9 is the turbulent Prandtl number. The
t
0:296kU1:435
h ¼ e ; ð11Þ roughness Stanton number is calculated from:
c qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
mR0sUe1:87ds Stk ¼1:92Rek(cid:2)0:45Pr(cid:2)0:8; ð14Þ
where k is the conductivity of air. This parameter is cal- where Pr = lC /k = 0.72 is the laminar Prandtl number.
p
culated by using viscosity computed from Sutherland’s TheskinfrictioniscalculatedfromtheMakkonenrelation:
viscosity law, assuming constant Prandtl number and spe- C 0:1681
f ¼ : ð15Þ
cific heat. Note that expression (11) is not dependent on 2 ½lnð864h=k þ2:568Þ(cid:6)2
t s
roughness.
For turbulent flow Re [600, the method of Kays and The turbulent momentum thickness is computed from:
k
Crawford is employed [3]. The turbulent convective heat 0 s 10:8
0:036m0:2 Z
transfer coefficient is evaluated from: ht ¼ U3:29 @ Ue3:86dsA þhtr; ð16Þ
h ¼StqU C ; ð12Þ e
c e p str
where C is the specific heat of air. The Stanton number where h is the laminar momentum thickness at transition
p tr
can be calculated from: location.
123
HeatMassTransfer(2009)45:305–322 311
Theboundary-layercalculationsstartattheleadingedge Table1 Parametervaluesusedinthecalculations
andproceeddownstream usingthe marchingtechniquefor
Symbol Definition Value
the upper and lower surfaces of the airfoil. Transition is
fixed at the streamwise location where Rek = 600, Cp Specificheatofair 1,006J/kgK
according to Von Doenhoff criterion. C Specificheatofice 2,050J/kgK
pi
C Specificheatofwater 4,218J/kgK
pw
e Saturationvaporpressureconstant 27.03
2.4 Extended Messinger model o
g Gravitationalacceleration 9.81m/s2
Theiceshapepredictionisbasedonthestandardmethodof ki Thermalconductivityofice 2.18W/mK
phase change or the Stefan problem. The phase change kw Thermalconductivityofwater 0.571W/mK
problemisgovernedbyfourequations:energyequationsin Le Lewisnumber 1/Pr
theiceandwaterlayers, massconservationequationanda LF Latentheatofsolidification 3.3449105J/kg
phase change condition at the ice/water interface [8]: LE Latentheatofevaporation 2.509106J/kg
L Latentheatofsublimation 2.83449106J/kg
oT k o2T S
¼ i ; ð17Þ Pr LaminarPrandtlnumberofair 0.72
ot qC oy2
i pi Pr TurbulentPrandtlnumberofair 0.9
t
oh k o2h e Radiativesurfaceemissivityofice 0.5–0.8
ot ¼qwCwpwoy2; ð18Þ lw Viscosityofwater 1.795910-3Pas
q Densityofrimeice 880kg/m3
r
oB oh
qi ot þqw ot ¼qabV1þm_in(cid:2)m_e;s; ð19Þ qg Densityofglazeice 917kg/m3
q Densityofwater 999kg/m3
w
oB oT oh r Stefan–Boltzmannconstant 5.6704910-8
qL ¼k (cid:2)k ; ð20Þ r
i F ot i oy woy r Surfacetensionofwater 0.072N/m
w
wherehandTarethetemperatures,k andk arethethermal
w i
conductivities,C andC arethespecificheatsandhand oh
pw pi Glazeice:(cid:2)k ¼ðQ þQ þQ þQ Þ
B are the thicknesses of water and ice layers, respectively. woy c e d r
In Eq. (19), qabV?, m_in and m_e;s are impinging, runback (cid:2)ðQaþQkþQinÞ aty¼Bþh:
andevaporating(orsublimating)watermassflowratesfora
ð23Þ
controlvolume(panel),respectively.InEq. (20),q andL
i F
oT
denotethedensityoficeandthelatentheatofsolidification
Rimeice:(cid:2)k ¼ðQ þQ þQ þQ Þ
of water, respectively. Ice density is assumed to have two ioy c s d r
different values for rime ice (q) and glaze ice (q ), see (cid:2)ðQ þQ þQ þQÞ aty¼B:
r g a k in l
Table 1.Thecoordinateyisnormaltothesurface.Inorder
ð24Þ
todeterminetheiceandwaterthicknessestogetherwiththe
temperaturedistributionateachlayer,boundaryandinitial 4. Airfoil surface is initially clean:
conditions must be specified. These are based on the fol- B¼h¼0; t¼0: ð25Þ
lowingassumptions [8]:
1. Iceisinperfectcontactwiththeairfoilsurface,whichis
In the current approach, each panel constituting the
takentobeequaltotheairtemperature,T inthisstudy:
a
airfoil is also a control volume. The above equations are
Tð0;tÞ¼T ¼T : ð21Þ
s a written for each panel and ice is assumed to accumulate
2. The temperature is continuous at the ice/water perpendicularlytoapanel.Thisisanextensionoftheone-
boundary and is equal to the freezing temperature: dimensional model described by Myers [8] to two-dimen-
sional, which is accomplished by taking mass and energy
TðB;tÞ¼hðB;tÞ¼T : ð22Þ
f terms due to runback water flow in the conservation
3. At the air/water (glaze ice) or air/ice (rime ice) equations into account, see Eq. (19). The differences
interface, heat flux is determined by convection (Q ), betweenthepresentmodelandtheonedescribedbyMyers
c
radiation (Q), latent heat release (Q), cooling by aresubtle.Themaindifferenceisrelatedtothemulti-layer
r l
incoming droplets (Q ), heat brought in by runback approach, where the exposure time is subdivided into
d
water (Q ), evaporation (Q ) or sublimation (Q), segments or layers. In this case, the boundary condition
in e s
aerodynamic heating (Q ) and kinetic energy of given in Eq. (25) is modified so that each layer except the
a
incoming droplets (Q ): first one starts with a non-zero ice thickness. Another
k
123
312 HeatMassTransfer(2009)45:305–322
differencewhichislessimportantisthatthepresentmodel ðQ þQ þQ Þ(cid:2)ðQ þQ þQ þQ Þ
hðyÞ¼T þ a k in c d e r ðy
takes the heat loss due to radiation into account. f k
w
(cid:2)BÞ:
2.4.1 Rime ice growth and temperature profile
ð31Þ
Rime ice thickness can be obtained directly from the mass Mass conservation Eq. (19) is integrated once to obtain
conservationEq. (19)aswaterdropletsfreezeimmediately the expression for water height, h:
on impact [8]: q bV þm_ (cid:2)m_ q
h¼ a 1 in eðt(cid:2)t Þ(cid:2) gðB(cid:2)B Þ; ð32Þ
q bV þm_ (cid:2)m_ q g q g
BðtÞ¼ a 1 in st: ð26Þ w w
q
r whereB andt aretheicethicknessandthecorresponding
g g
It has been shown that, for ice thicknesses less than time at which glaze ice first appears, respectively. When
2.4 cm (which is mostly the case), the temperature Eq. (32) is substituted into the phase change condition in
distribution is governed by [8]: Eq. (20), a first order ordinary differential equation for the
ice thickness is obtained:
o2T
¼0: ð27Þ (cid:2) (cid:3)
oy2 oB ki Tf (cid:2)Ts
q L ¼
g F ot B
Integrating the above equation twice and applying ðQ þQ þQ þQ Þ(cid:2)ðQ þQ þQ Þ
þk c e d r a k in :
conditions given in Eqs. (21) and (24) yields the w k
w
temperature distribution in the rime ice layer as:
ð33Þ
TðyÞ¼T
s During transition from rime ice to glaze ice, ice growth
ðQ þQ þQ þQÞ(cid:2)ðQ þQ þQ þQ Þ
þ a k in l c d s r y: rate must be continuous:
k
i (cid:5) (cid:5)
oB oB
ð28Þ ¼ at B¼B or t¼t ð34Þ
ot ot g g
rime glaze
2.4.2 Glaze ice growth
Using Eqs. (26) and (33) yields:
It has been shown that, if ice and water layer thicknesses
arelessthan2.4 cmand3 mm(whichisthecaseformost
(cid:2) (cid:3)
k T (cid:2)T
i f s
B ¼ ; ð35Þ
g ðq bV þm_ (cid:2)m_ ÞL þðQ þQ þQ Þ(cid:2)ðQ þQ þQ þQ Þ
a 1 in sub F a k in c d e r
applications), respectively, the temperature distributionsin q
t ¼ r B : ð36Þ
the ice and water layers are governed by [8]: g q bV þm_ (cid:2)m_ g
a 1 in sub
o2T o2h Inordertocalculatetheglazeicethicknessasafunction
¼0; ¼0: ð29Þ
oy2 oy2 oftime,Eq. (33)isintegratednumerically,usingaRunge–
Kutta–Fehlberg method.
After integrating above equation twice and employing
conditions(21)and(22),thetemperaturedistributioninthe
2.4.3 Energy terms
ice becomes:
T (cid:2)T
TðyÞ¼ f syþT : ð30Þ Theenergytermsappearingintheaboveequationsneedto
B s
be expressed in terms of the field variables. Although
The temperature distribution in the water layer is convective heat transfer (Q ) and latent heat (Q) are the
c l
obtained by integrating Eq. (29) twice and employing most prominent terms, all relevant energy terms are con-
conditions (22) and (23): sideredhere,andusedinthecomputerprogramdeveloped.
123
HeatMassTransfer(2009)45:305–322 313
Inthesubsequentformulation,T isthetemperatureatthe 2.4.4 Rime ice temperature distribution
sur
ice surface (for rime ice) or the water surface (glaze ice).
Equation (28) becomes:
• Convective heat transfer at the water surface (Q ):
c
Q þQ T
Q ¼h ðT (cid:2)T Þ: ð37Þ TðyÞ¼ 0r 1r syþT ; ð48Þ
c c sur a k (cid:2)Q B s
i 1r
• Cooling by incoming droplets (Q ):
d where
Q ¼q bV C ðT (cid:2)T Þ: ð38Þ
d a 1 pw sur a oB V2 V2
Q ¼q L þq bV 1þrh 1 þq bV C T
• Evaporative heat loss (Q ): 0r r F ot a 1 2 c2C a 1 pw a
e pa
Q ¼v e ðT (cid:2)T Þ; ð39Þ þh T þ4er T4þv e T þm_ C T ;
e e o sur a c a r a s 0 a in pw f
where v is the evaporation coefficient and e = 27.03. ð49Þ
e o
Evaporation coefficient is expressed as [8]: Q1r ¼qabV1Cpwþhcþ4errT13 þvse0þm_inCpw: ð50Þ
v ¼0:622hcLE; ð40Þ 2.4.5 Glaze ice temperature distribution and ice growth
e CpPtLe2=3 rate
where P is the total pressure of the airflow.
t Equation (31) can be written as:
• Sublimation heat loss (Qs): hðyÞ¼Q0þQ1Tf hþT ; ð51Þ
Qs ¼vseoðTsur(cid:2)TaÞ; ð41Þ kw(cid:2)Q1h f
Sublimation coefficient v is expressed as [8]: where
s
vs ¼0C:p6P2t2Lhec2L=3S: ð42Þ Q0 ¼qabV1V212 þrhc2VC12paþqabV1CpwTaþhcTa
þ4er T4þv e T þm_ C T ; ð52Þ
• Heat loss due to radiation (Q): r a e 0 a in pw f
r
Qr ¼4errTa3ðTsur(cid:2)TaÞ; ð43Þ Q1 ¼qabV1Cpwþhcþ4errT13 þvee0þm_inCpw: ð53Þ
where e is the surface emissivity and r is the Stefan- Equation (33) can be written as:
r
Boltzmann constant. oB T (cid:2)T Q þQ T
q L ¼k f s(cid:2)k 0 1 f : ð54Þ
• Aerodynamic heating term ðQaÞ: g F ot i B w kw(cid:2)Q1h
rh V2
Q ¼ c 1; ð44Þ Equation (35) can be written as:
a
2Cp k(cid:2)T (cid:2)T (cid:3)
i f s
where r is the adiabatic recovery factor (r ¼Pr1=2 for Bg ¼q L (cid:6)qabV1þm_in(cid:2)m_sub(cid:7)þ(cid:2)Q þQ T (cid:3): ð55Þ
laminar flow, r ¼Pr1=3 for turbulent flow). g F qr 0 1 f
• Kinetic energy of incoming droplets (Q ): 2.4.6 Freezing fractions and runback water
k
V2
Q ¼q bV 1; ð45Þ Freezing fraction foragivencontrolvolume (orapanelin
k a 1 2
this case) is the ratio of the amount of water that solidifies
• Energy brought in by runback water (Q ): totheamountofwaterthatimpingesonthecontrolvolume
in
plus the water that enters the panel as runback water.
Q ¼m_ C ðT (cid:2)T Þ; ð46Þ
in in pw f sur
q B
wherem_inisthemassflowrateoftheincomingrunback Rime ice: FF¼ r : ð56Þ
ðq bV þm_ Þt
water. a 1 in
• Latent heat of solidification (Ql): Glaze ice: FF¼qrBgþqg(cid:2)B(cid:2)Bg(cid:3): ð57Þ
oB ðq bV þm_ Þt
Q ¼q L : ð47Þ a 1 in
l r F ot
Runback water mass flow rate:
With these definitions, it is possible to express
m_ ¼ð1(cid:2)FFÞðq bV þm_ Þ(cid:2)m_ : ð58Þ
Eqs. (28), (31), (33) and (35) in terms of the airfoil out a 1 in e
surface temperature (T) and ambient temperature (T ) This becomes m_ for the neighboring downstream
s a in
only. panel. It is assumed that, all unfrozen water passes to the
123
314 HeatMassTransfer(2009)45:305–322
nextdownstreampanelfortheuppersurface.Forthelower a 0.1
surface,itisassumedthatalltheunfrozenwaterisshed[2]. Clean
Experimental
0.075 PresentStudy(4layers)
PresentStudy(1layer)
2.4.7 Evaporating or sublimating mass
0.05
Evaporating or sublimating mass is given as [11]:
m)0.025
0:7 (cid:4)p (cid:2)p (cid:5) y (
m_ ¼ h v;sur v;1 : ð59Þ
e;s C c P 0
pa 1
pv,sur and pv,? are the vapor pressures at the ice or water -0.025
surface and the ambient air, respectively. These are
computed from [11]: -0.05
-0.05 -0.025 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175
p ¼3386(cid:2)0:0039þ6:8096(cid:3)10(cid:2)6T(cid:1)2þ3:5579(cid:3)10(cid:2)7T(cid:1)3(cid:3); x (m)
v
ð60Þ b
0.1
T(cid:1)¼72þ1:8ðT (cid:2)273:15Þ: ð61Þ Clean
Experimental
0.075 PresentStudy(4layers)
DRA
NASA
ONERA
3 Code validation 0.05
m) 0.025
Inordertovalidatethedevelopedtool,theobtainedresults y (
are compared with experimental ice shapes [10] over a
0
NACA 0012 airfoil. Four test cases with significantly
varyingambienttemperaturesareselected(T = -27.8,-
a -0.025
19.8, -13.9 and -6.7(cid:2)C) for comparison. Numerical data
obtained by different research groups corresponding to -0.05
-0.05 -0.025 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175
these conditions are available in the literature [14], which
x (m)
is a reason why they were selected. The experimental and
numerical data have been reproduced from Wright et al. Fig.6 Comparison of ice shape predictions for NASA 27 case
(T =-27.8(cid:2)C)
[14].Geometricandflowconditionscorrespondingtothese a
cases are presented in Table 2.
InFig. 6,iceshapescorrespondingtoTa = -27.8(cid:2)Care predictions are similar and agree well with the experi-
presented. Figure 6a compares the ice shapes obtained in mental shape.
the current study with those reported by Olsen et al. [10]. InFig. 6b,obtainediceshapesarecomparedwiththose
The conditions are rime ice conditions as made evident by obtained numerically by DRA, NASA and ONERA,
theobtainediceshapes.Ascanbeobserved,thedeveloped respectively. All the codes (including the current one)
code predicts the experimental shape and the iced region predict similar shapes, all agreeing fairly well with the
very well. For the current study, the results are presented experimental one. As rime ice shapes are obtained from a
for one-layer and four-layer calculations. One-layer cal- simple algebraic equation, the good agreement observed
culations predict a slightly higher ice volume, but both here is expected.
InthecasepresentedinFig. 7,thetemperatureishigher,
but all other parameters remain the same as in Fig. 6.
Table2 Geometric and flow conditions for code validation Again, the results of the current study are presented using
calculations
one-layer and four-layer calculations. There is a marked
Variables Value difference between the results of two approaches as illus-
trated in Fig. 7a. One-layer calculations predict a typical
a,angleofattack 4(cid:2)
glaze ice shape, while four-layer calculations predict a
c,airfoilchord 0.53m
typical rime ice shape. However, it is the ice shape pre-
V ,freestreamvelocity 58.1m/s
? dicted bythe four-layerapproach thatagrees well with the
p ,ambientpressure 95610Pa
? experimental shape.
q,liquidwatercontent 1.3g/m3
a InFig. 7b,obtainediceshapesarecomparedwiththose
t ,exposuretime 480s
exp obtained numerically by DRA, NASA and ONERA,
d dropletdiameter 20lm
p respectively. The results of the current code and DRA
123