Towards modelling of initial and final stages of supercooled water solidification

International Journal of Thermal Sciences 92 (2015) 150e161

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Towards modelling of initial and final stages of supercooled water solidification A. Criscione a, b, *, I.V. Roisman a, b, S. Jakirli
c a, b, *, C. Tropea a, b a b

€t Darmstadt, Germany Institute of Fluid Mechanics and Aerodynamics, Technische Universita €t Darmstadt, Germany Center of Smart Interfaces, Technische Universita

a r t i c l e i n f o

a b s t r a c t

Article history: Received 25 March 2014 Received in revised form 21 January 2015 Accepted 22 January 2015 Available online 3 March 2015

It is well known from experimental studies (e.g., Jung et al., 2012, [13]), that solidification of supercooled water occurs in two stages: the initial rapid, recalescent stage of crystallization and the final slower, quasi-isothermal stage. To date, it is not well understood how these stages can be modelled. In the present study it is demonstrated that both freezing stages of supercooled water can be mathematically modelled by utilizing exclusively the minimal (pertinent to the heat-diffusion description) model of solidification: the first unstable stage is physically well described by the dendritic-growth-approach, whereas the second stage of the freezing process is qualitatively and quantitatively defined by a stable planar solidification. The computational model, introduced initially in [7,31], has been appropriately upgraded with respect to the high-fidelity discretization of the thermal-energy equation as well as the normal-to-the-interface temperature gradient determination. At supercooling degrees DT higher than z5 K (corresponding to high supercooling range), the previously obtained numerical results related to the unstable freezing exhibit a distinct deviation from available experimental data. Possible reasons for this departure can be either the non-consideration of the thermal interaction between the growing needle-like dendrites or neglect of the kinetic effects influencing the growth at such high supercooling conditions. Two-dimensional computations concerning the growth of a needle-like dendrite surrounded by an array of needles indicate that they do not considerably influence the steady-state tip velocity of an isolated needle at higher supercooling. Therefore, the thermal energy equation has been extended through an appropriately derived term mimicking the kinetic undercooling phenomenon. The relevant computational validation results in a very good qualitative and quantitative agreement with experiments in the high-supercooling range as well, illustrating the model's predictive capability of describing both stages of crystallization: the first rapid, dendritic growth (pertinent to unstable freezing) and the second planar growth stage (corresponding to the stable freezing). © 2015 Elsevier Masson SAS. All rights reserved.

Keywords: Supercooled water Solidification Array of needle-like dendrites Kinetic effects Heat-diffusion-driven growth

1. Introduction There are many unanswered questions about the basic physics of supercooled water crystallization [2,8,22,30]. Water drops which exist in the liquid form at temperatures below 273.15 K are termed supercooled. They often occur in clouds located at altitudes which aircrafts have to pass during start and landing. The terminology used frequently in regard to the in-flight airframe icing classifies them as Supercooled Large Droplets (SLD). The SLD are defined in

* Corresponding authors. Institute for Fluid Mechanics and Aerodynamics, €t Darmstadt, Germany Technische Universita E-mail addresses: [email protected] (A. Criscione), s.jakirlic@sla. tu-darmstadt.de (S. Jakirli
c). http://dx.doi.org/10.1016/j.ijthermalsci.2015.01.021 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.

accordance with The World Meteorological Organization as those water drops with a diameter larger than 50 microns. Usually these drops impact the leading edge area of the airfoil. This may result in the formation of an ice layer leading to increased aerodynamic drag and weight, associated with a reduction in lift and thrust [29]. The in-flight airframe icing represents actually a significant aviation hazard, keeping in mind numerous serious difficulties in controlling this phenomenon [28]. Detection and identification of aircraft icing [4,23] is a topic of current interest, aiming at a better understanding of the physics underlying this phenomenon. It should be mentioned that icing is not exclusively an aviation problem. For example, wind turbines can also be significantly affected by icing. An ice layer on an operating turbine causes additional drive-train loads, exceeding often the design loads. For this reason many

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

turbines are shut down at icing onset. A restart is activated only after an inspection confirms that the ice no longer represents a danger to operation. This is regarded as a cumbersome practice, especially in remote locations and at night, limiting drastically the economy of wind power plants [11,16,17]. At the instant corresponding to the impact of the supercooled water drop onto a solid substrate, the crystallization of the water is initiated, causing possible growth of any free nuclei in the supercooled liquid, often creating dendritic-shaped crystals. When water crystallizes, latent heat of fusion is released. If an SLD is going to freeze only a small portion of the drop will freeze instantaneously, not necessarily more than enough to raise the temperature up to the crystallization/melt temperature of 273.15 K. In other words the drop freezes until the initial supercooling DT (representing basically the liquid temperature reduction below its freezing point without solidification) in the liquid is exhausted. Furthermore, progressive freezing takes place as the drop loses heat by evaporation and conduction. Accordingly, the solidification of supercooled water occurs in two stages [3,8,30]. Recently, Jung et al. [13] have experimentally investigated the freezing mechanism of an SLD residing on a superhydrophobic surface without external flow, Fig. 1. The drop has a volume of 5 ml and the initial supercooling corresponds to DT ¼ 15 K. Both the ambient and the surface have a constant initial temperature of 258.15 K. In the absence of any flow the pure diffusive evaporation rate of the supercooled water drop is negligibly small and a heterogeneous nucleation at the solideliquid interface, corresponding to lower critical energy barrier, initiates freezing. The solidification occurs in two stages: the snapshots aed in Fig. 1, recorded at 0, 8, 16 and 32 ms respectively, display the first stage of solidification, whereas the snapshots eeh taken at 0, 6, 13 and 20 s illustrate the progress of the second stage. The first crystallization phase is obviously much faster than the second one: the average freezing front velocity of the first stage corresponds approx. to 6.6 cm/s whereas the solidification front in the second stage has an average velocity value of 0.11 mm/s; this is a significant velocity difference. Important issues of discussion are the substantially different structural properties of the solidification process. The initially generated structure corresponds to a milky opaque appearance. In the second stage the solidification structure appears more clear-transparent; in the snapshot h (at 20 s) it can be recognized that the drop increases its initial volume as well. To date, it is not completely understood how these stages can be physically modelled.

151

Theoretically, the common approach to crystallization problems relies on the grade of the initial supercooling in the liquid. At low supercooling, the freezing process can be considered as a Stefan problem corresponding to thermal diffusion-driven growth [1,8]. The solideliquid interface is an active free boundary (actually a moving boundary) from which the latent heat of freezing is released during the phase transformation. This heat is conducted away from the interface through both the liquid and the solid phase. The heat equation is solved in each phase separately. Both temperature fields, in the liquid and in the solid phase, are coupled through two boundary conditions at the moving boundary. The first boundary condition is the constant melting/crystallization temperature at the interface. This temperature will be locally altered by an amount depending on the interfacial tension between the solid and liquid phases and the local curvature, in accordance with the GibbseThomson effect. The second boundary condition is represented implicitly by the velocity of the moving boundary. It is derived from a heat balance at the interface (in line with the Stefan condition); hence it depends on how fast the latent heat of solidification is removed from the interface. In the case of increasing supercooling, the literature reports that either the non-consideration of the thermal interaction between the growing needle-like dendrites [35] or neglecting the kinetic effects [34] can influence the growth at such high supercooling conditions. The transition point has been observed in experiments at a supercooling of approximately 5 K. The shape of the solideliquid interface depends on various parameters such as the initial degree of supercooling, the surface tension and its anisotropy. In case that the initial supercooling of the liquid is low and the heat flux into the solid is high, planar solidification is likely to occur. The latent heat of solidification released at the interface is mainly transferred into the solid phase. The more heat is conducted into the liquid the more unstable becomes the shape of the solideliquid interface. At high degrees of supercooling in the liquid phase, the planar solidification front transforms to a dendrite-shaped surface with side-branching growth. Even arrays of branchless dendrites (needle-like crystals) emerge at higher degrees of supercooling. Shibkov et al. [32e34] investigated the free growth of an ice crystal in pure supercooled water. The experiments were performed on a horizontal 200 mm thick water film stretched over area of 30 mm2. In order to seed an ice crystal in the supercooled water and to trigger the nucleation process they used a thin steel frosted rod. Various shapes of ice crystal patterns in the range between low (DT ¼ 0.1 K) and high (DT ¼ 14.5 K) supercooling have been

Fig. 1. SLD crystallizing on a superhydrophobic surface [13]: first (aed) and second (eeh) stage of solidification of a 5 ml pure water drop on a superhydrophobic surface without convective flow. The initial supercooling of the drop values 15 K. Both the ambient and the surface have also a constant initial temperature of 258.15 K. The origin and direction of the arrows, respectively, mark the instantaneous position and direction of motion of the solidification front.

152

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

observed: for a high supercooling an array of needle-like dendrites has been observed. It becomes evident that the shape of the solidification front depends strongly on the initial degree of supercooling in the water film. Furthermore, the shape of the dendritic front is determined by the balance between surface energy criterion and the efficiency of the interface in removing heat. Over the last decades, various computational approaches have been proposed and applied to simulate solidification especially in metals. However, the existing computational models describing the crystallization of pure water are not numerous. The direct solution of a time-dependent moving-boundary problem, which governs both the planar solidification and dendritic crystal growth, represents a significant challenge. To date, the phase-field method has been usually used in simulations describing crystallization phenomena to avoid the difficulties in tracking a sharp boundary [9,15,39]. In the present work, a computational model based on the level set technique [27] used in Refs. [6,7,31] has been appropriately modified, aiming at increasing the accuracy level of the discretization of the equation governing the thermal energy and the normal-to-interface derivative of the temperature field. The model describes the freezing mechanism under supercooled conditions, relying on the physical and mathematical description of the twophase moving-boundary approach (minimal model of solidification) considering kinetic effects. The relevant numerical algorithm is implemented into the open source software OpenFOAM® and is based on the algorithm presented in Ref. [7]. In the following section, a brief outline of the physics governing the freezing mechanism and its mathematical description is given. The computational algorithm used in this work is introduced in section 3. The computational results regarding the different freezing stages of supercooled water solidification are discussed in section 4. Finally, based on the present computational study, a physical interpretation of the two different stages is presented.

2. Mathematical modelling 2.1. Governing equations for a stable solidification Presently, the phenomenon of solidification is considered mathematically as a two-phase moving boundary case [1,8]. Both the liquid and the solid phases are active, i.e. the heat conservation is solved in both sub-domains. Hence, heat flows from the interface in both directions. We consider a domain D of pure water where the water is either in liquid (supercooled) or in solid phase, Fig. 2; here, T represents the temperature. The motion within the liquid region is not presently considered. The energy equation describing the

time-dependent heat conduction in both regions reduces correspondingly to:

vT ¼ V$ðks VTÞ; vt

x2Us

(1)

vT ¼ V$ðkl VTÞ; vt

x2Ul ;

(2)

rs cv;s

rl cv;l

where r defines the density, cv the volumetric heat capacity and k represents the thermal conductivity. The indices s,l denote the solid and the liquid phase, respectively. The physical properties in the solid and liquid phases are defined in Table 1.At the interface (denoted by X), two boundary conditions (BC) are needed to close the computational model. The first one (1st BC) describes the heat flux balance in a region close to the interface. This BC is also known as the Stefan condition. It expresses the dependence of the local normal velocity of the interface on the heat flux discontinuity within the interface region (the quantities with the subscript X relate to the solideliquid interface itself):


 rLvn ¼ kl VT 

X;l

   ks VT 

 X;s

$n;

x2X;

(3)

where L is the latent heat released by the phase change (Table 1), n is the unit normal vector directed towards the liquid region and vn denotes the magnitude of the solidification speed. Here, r ¼ rs ¼ rl is assumed.The second boundary condition (2nd BC) at the interface relates to the temperature distribution. In the case of a curved interface the temperature at the interface, TX, needs to be described by the GibbseThomson relation:

TX ¼ Tm  DTn ¼ Tm  Gk;

G¼

sTm ; rL

x2X;

(4)

with Tm describing the melting temperature of water, s representing the interfacial energy, k the curvature of the interface, G denoting the capillary constant and DTn is the capillary undercooling.The solideliquid interface is in dynamic equilibrium when the molecules attach and detach continually at equal rates. In this case of equilibrium the temperature of the interface, TX, is equal to Tm. In the case of a curved solideliquid interface the interface temperature is appropriately altered in order to account for the capillary undercooling, Eq. (4). In the published literature it is suggested that the molecules become stronger bounded to the interface if TX < Tm [1,8]. Thus, their number detaching per unit time decreases. The interface velocity, vn, increases, according to the difference between Tm and TX, until the point where the molecules in the liquid become sluggish and the rate of attachment decreases. Hence, for pure substances with low latent heats like metals [8], the velocity of the solideliquid interface can be approximated as a linear function of the kinetic undercooling, DTkin:

vn ¼ kkin ðTm  TX Þ ¼ kkin DTkin ;

(5)

where kkin represents the kinetic coefficient. Accordingly, the effect of the kinetic undercooling at the interface can be modelled as follows:

Table 1 Physical properties of both liquid and solidified water. Fig. 2. Schematic of the flow domain considered: the water is either in the liquid, supercooled phase or in the solid one. T(x,t) represents the temperature as a function of position and time. Us describes the region of solidified water, whereas Ul defines the region of liquid water. X is the interface between the two regions.

Water as

r [kg/m3]

L [kJ/kg]

k [W/(m K)]

cp [kJ/(kg K)]

Liquid Solid

1000 917

333 e

0.59 2.2

2.06 4.18

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

TX ¼ Tm  k1 kin vn :

(6)

Assuming that the latent heat of pure metals is comparable to that of pure water (L ¼ 398 kJ/kg for Aluminium, 368 kJ/kg for Magnesium and 333 kJ/kg for Water), the afore-mentioned approximation should be valid for water as well. Hence, in accordance with some references, the interface velocity can be approximated as a linear function of the initial supercooling. 2.2. From stable planar freezing interface to the unstable growth of needle-like dendrites The solideliquid interface, representing the phase-transition region where solid and liquid coexist, is characterized by a planar shape for most pure materials under ordinary freezing conditions, implying a constant Tm, Fig. 3 I. During the crystallization of supercooled water the interface between the solid and the liquid phase becomes unstable resulting in formation of small interface curvature perturbations. The supercooling acts towards an interface destabilization with the solidification rate depending on the degree of the initial-supercooling, which drives this process. On the contrary, the interfacial tension tends to stabilize the interface, damping small perturbations back in line with the GibbseThomson relation, Fig. 3 II. The balance between the stabilizing and destabilizing effect can be analyzed by the classical approach to morphological instability introduced by Mullins and Sekerka [24,25] in the context of directional solidification, Fig. 3 III. The analysis of the morphological instability yields the cutoff wavelength representing the largest possible wavelength for a stable interface:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2aTm scn pffiffiffiffiffiffiffiffiffi lc z z dd dc ; vn;0 L2v

(7)

where a is the thermal diffusivity, Lv ¼ Lr and vn,0 represents the characteristic velocity of the system. dd ¼ 2a/vn,0 and dc ¼ Tmscn/L2v represent the thermal diffusion and the capillary length scales, respectively. The wavelength of the fastest growing mode, lf, is defined as [10].

Fig. 3. Sequence of events, from planar freezing front to the growth of needle-like dendrites; see Sections 2.1 and 2.2 for meaning of the quantities and symbols.

lf ¼

pffiffiffi 3lc :

153

(8)

At this wavelength, the interfacial perturbations should experience the fastest growth into the supercooled liquid. Assuming that the thermal diffusivity, a, and the volumetric heat capacity, cv, are constant, the cutoff wavelength depends directly on vn,0. If the value of the interface velocity increases, the largest possible wavelength for a stable interface decreases. Thus, the cutoff wavelength gets smaller at higher supercooling for the liquid phase subjected to the condition which implies a constant temperature gradient in the solid. This theoretical statement is verified experimentally. Shibkov et al. [32e34] observed that when the surface of a very thin film of supercooled pure water is slightly touched by a metallic needle with a frosted tip (trigger), an ice crystal freely grows from the point of contact. The temperature in the crystal suddenly jumps to Tm. The temperature jump is caused by the release of the latent heat of crystallization. Regardless of the level of the initial supercooling in the water, the incipient structure of the crystal resembles a smooth circular disk at a constant temperature corresponding to Tm. Only when the radius of the initial nucleus is larger than the critical radius, Rc, it will grow towards R ¼ ∞, otherwise the nucleus will implode due to the system-dominating surface energy. The nucleus grows as a smooth spherical particle as far as the critical radius of instability, RI, is reached [1,8]. When the nucleus radius is larger than RI, the supercooling starts to destabilize the system: small bumps, i.e. instabilities at the interface start to grow into the supercooled water (Fig. 4). In this case, heat conduction in the solid would stabilize the particle. As the temperature in the nucleus is constant the conduction in the solid can be regarded as negligible. In a second step, the disk develops some bulges/bumps which grow into fingers representing the initial state of the unstable growth. These bumps are exposed to a steeper temperature gradient and evolve consequently into the deformation of fingers. During this transition process the primary fingers develop into crystals of different shapes: from isolated dendrites with sidebranching to an array of needles. This transition depends strongly on the initial supercooling degree of the pure water; hence, they are strictly correlated with the morphological cutoff wavelength. The spacing between bumps, Fig. 4, at the solideliquid interface resulting from the morphological instability theory of Mullins and Sekerka is in close agreement with the critical cutoff wavelength, Eq. (7). Thus, perturbations at the interface evolve into the formation of various

Fig. 4. A sketch of a cross section of a growing spherical nucleus. Between critical nucleation radius, Rc, and critical radius of instability, RI, the interfacial energy resists opposing growth of bumps (smooth circular disk). Corresponding temperature profile of the initial nucleus is shown. When the nucleus radius is larger than RI, the destabilizing effect (supercooling) starts to dominate the system: small bumps on the interface start to grow into the supercooled water.

154

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

polycrystalline structures, depending on the initial supercooling degree. Ivantsov [12] derived an analytical solution for the steady-state tip velocity and the tip radius of growing dendrites, under the following assumption: first, the kinetics is instantaneous, k1 ¼ 0, kin second, the surface tension is equal zero, s ¼ 0, and third, the crystal tip has a parabolic structure. This theory considers a single needle growing into an infinite half-space of supercooled liquid.
clet number, Pe, is used to establish a correlation between The Pe the tip velocity, vt, and the tip radius, rt:

Pe ¼

vt rt : 2a

(9)


clet number valid for the entire range of There is a unique Pe dimensionless supercooling, D ¼ DT/(L/cp), obtained by solving the Ivantsov function:

D ¼ Pe0 expPe0 E1 ðPe0 Þ:

(10)

The solution of this function represents a continuous family of parabolas/paraboloids for a given initial dimensionless supercooling D. One of the first important hypothetical principles regarding the tip-shape selection was the Marginal Stability Theory (MST) introduced by Langer and Müller-Krumbhaar [18e20] in their work on Universal Law of Crystal Growth. They analyzed the stable steady-state of the Ivantsov paraboloidal dendrite introducing the interfacial tension effect as a linearized perturbation function and found that the Ivantsov's continuum family of solutions may be divided into a stable and an unstable region. When interfacial energy, s, is present and kinetic supercooling is neglected, a new dimensionless parameter can be defined, representing the so-called controlling parameter for the operating point of the needle:

ε¼

Tm scn 2a dc dd lc ¼ 2 ¼ : 2prt rt L2n vn rt2

(11)

The selected dendrite corresponds to the point of marginal stability; hence the emerging tip radius, rt, corresponds to the cutoff wavelength, Eq. (7). Eq. (11) provides an additional relation between the tip growth velocity and the tip radius. The first relation
t number, Eq. (9). On the basis of is provided in terms of the Pecle the MST, a unique tip growth velocity as a single valued function of the supercooling can be calculated. The solid line in Fig. 5 represents the dimensionless steady-growth speed of freely growing ice crystals obtained theoretically as a function of the dimensionless supercooling [7]. The MST assumes a local thermal equilibrium at growing solideliquid interface which may not be the case at high supercooling degrees, DT  5 K, because of the slow water molecular transfer across the interface (sluggish interface kinetics). Following published literature, for higher supercooling degrees the speed of solidification should predominately be determined by the mechanism of molecular attachment kinetics at the solideliquid interface. A further possible reason for the deviation of the experiments from the MST at higher supercooling degrees than DT ¼ 5 K, Fig. 5, can be the thermal interaction between growing needles: in experiments [32] an array of needles growing into the liquid is observed; thus, the thermal field of the monitored crystal tip can be influenced by neighbouring needles, decreasing the steady tip velocity. In the limit of very small supercooling (DT / 0), Spencer and Huppert [35e38] derive analytical solutions which confirm the thermal influence of surrounding needles on the steady tip velocity of the observed crystal depending on the spacing between the needles. Both

Fig. 5. Ice crystals freely growing from supercooled pure water with tip velocity, vt, as a function of the initial supercooling, DT. The solid line represents the diffusional MST of Langer and Müller-Krumbhaar with the controlling parameter, ε* ¼ 0.025. Computational [7] and experimental results of Fujioka and Langer [21], Ohsaka and Trinh [26], and Shibkov et al. [32e34] show a significant deviation from the theoretical curve at high supercooling.

possible reasons for deviation, i.e. kinetic effects and thermal interaction of an array of needles, will be considered and discussed in section 4. 3. Numerical procedures The computational model used in Refs. [7,31] has been correspondingly modified towards fulfilment of an appropriate accuracy level with respect to the thermal energy equation discretization and the approximation of the temperature normal derivative close to the interface. The level set function, F, is defined as a signed distance function [27]: the zero level set function coincides with the solideliquid interface, whereas the signed value of the outer level set field represents the distance to the interface. The sign indicates the side of the interface one looks at, a negative sign of the outer level set field represents the solid phase while the positive sign indicates the liquid phase. Fig. 6 visualizes iso-surfaces of the level set field for the cross-section of a circular ice nucleus immersed in water. The heat transfer equations are solved in both, the liquid and solid phase independently from each other. At the solideliquid interface a constant value (Dirichlet BC) for the temperature field is imposed: the interface temperature is calculated according to the GibbseThomson relation, Eq. (4). For it, the curvature of the solideliquid interface, representing the zero level set

Fig. 6. Iso-surfaces of both the zero and the outer level set functions; the outer level set field constitutes the distance to the interface (zero level set function, F ¼ 0). The sign indicates the side of the interface, i.e. liquid phase, F > 0, and solid phase, F < 0.

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

function, is needed. In cells adjoining the interface, the curvature is calculated at the central node as the divergence of the unit normal vector, which is defined as n ¼ VF=jVFj. Subtracting the normal distance between the cell center and the interface (outer level set value) from the curvature radius, one obtains the interface curvature radius, whose inverse describes the curvature of the solideliquid interface [6,7,31]. An appropriate arrangement of the numerical stencil is applied, consisting of ghost-faces, gf, corresponding ghost-points, gp, and interface-points, ip (Fig. 7). It is necessary to fulfil the second boundary condition implying the energy balance at the interface. The interface temperature, TX, is imposed directly as a Dirichlet boundary condition for both the liquid and the solid temperature fields. Hence, the interface temperature is stored at ghost-points and interface-points. The two-dimensional temporal and spatial discretization of the energy equation in the liquid phase, using Euler and the central differencing schemes, leads to

rl cn;l

0 Ti;j  Ti;j

Dt

Vi;j

2 3   k Siþ1;j k Si;jþ1 k Si;j1 rl cv;l Vi;j l i;j l i;j l i;j d 2d i1;j 2 0 4 þkl Si;j þ iþ1;j þ i;jþ1 þ i;j1 5Ti;j Dt d1 d2 d1 d1 di;j di;j di;j 0 1 i;j1 i;jþ1 kl Si;j kl Si;j 1 2d0 d1 A iþ1;j Tiþ1;j þ i;j1 Ti;j1 þ i;jþ1 Ti;jþ1 þkl Si;j @ iþ1;j þ 2d2 d1 d2 d d d i;j

¼ rl cv;l

0 Ti;j

Dt

i1;j

Vi;j þkl Si;j



i;j

(12) 0 is the value stored at the grid node P at the previous where Ti;j i,j time step with V representing the cell volume. S is the face area, whereas d1 and d2 represent the characteristic lengths between cell center and ghost-point, Fig. 7. The indices gp,x and gp,y denote the ghost-points in x and y direction, respectively. In order to obtain the matrix format of the equation system discretized by employing a second-order accurate scheme the diffusive ghost-face term in Eq. (12) has to be decomposed as follows (it is demonstrated presently by means of the ghost-face oriented into x-direction; the y-direction-oriented ghost-face is regarded as a regular face within the liquid phase):

Fig. 7. Magnification of the ’i’ area in Fig. 6: Discretization of the thermal energy close to the interface; visualization of ghost faces, ghost and interface points. Direct imposition of Dirichlet boundary condition at the interface.



i;j

2d0 d2 2d0 d1 2d2 2d1 Tgp;x 2d2 d1 2d1 d2 (13)

In order to compute the temperature derivative normal to the interface, minimizing the grid anisotropy (which depends strictly on the mesh arrangement), the least-square-approach is used to approximate the temperature at a virtual node, ND, normal to the interface, i.e. perpendicular to the tangent line at the interfacepoint, ip:

vND ¼ v0ip  v0i;j þ maxðFic Þnip ;

Tiþ1;j  Ti;j iþ1;j Tgp;x  Ti;j i1;j ¼ kl S þ kl Si;j d2  d1 i;j d1 Ti;jþ1  Ti;j i;jþ1 Tgp;y  Ti;j i;j1 þ kl Si;j þ kl Si;j ; i;jþ1 i;j1 di;j di;j

155

(14)

v0ip and v0i;j represent the origin vector of the interface-point and the cell center, respectively with nip denoting the normal vector at the interface-point, Fig. 8; max(Fic) describes the maximal level set value among all cells bordering the interface, the so-called interface cells. In the two-dimensional case the following biquadratic approach,

TND ¼ Ti;j þ a0 xND þ a1 yND þ a2 ðxyÞND þ a3 ðxND Þ2 þ a4 ðyND Þ2 ; (15) is used to approximate the temperature at the point ND, Fig. 8. In order to compute unknown coefficients ai, the temperature is approximated at all center nodes of neighbouring cells (we recall that in a two-dimensional system at least five evaluation points have to be accounted for) leading to

Fig. 8. Approximation of the temperature in the virtual node (ND) using the least squares approach.

156

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

sðx; yÞ ¼

n X

aj fj ðx; yÞ ¼ A$a;

(16)

j¼1

where f(x,y) represents the basis function evaluated at each evaluation point. By minimizing the residuum of Eq. (16), the unknown coefficients are calculated as follows

a¼




AT A

1

AT sðx; yÞ:

(17)

Every entry in A, corresponding to an evaluation point, is weighted with the ratio between the distance of the closest evaluation point and the distance between the evaluation point under consideration and the interface cell center, Pi,j. By doing so, the evaluation points closer to the interface are more influential than the remote ones. Finally, the interface-normal derivative of the temperature field is approximated as follows:

  vT  T  Tip   $nz ND $n: vn dn;2 ip

(18)

The normal-to-interface velocity, vn, is estimated in accordance with the fulfilment of the heat flux balance close to the interface, Eq. (3). Within a narrow band around the interface, whose width is temporally adjusted to the maximum curvature of the interface, the velocity vn is appropriately expanded. This velocity is used to update the level set function within the narrow-band. Afterwards the outer level set function is re-initialized to enable the setting of a new band around the interface at the consequent time step. For detailed description of this update and re-initialization procedure we refer to [6,7,31]. The numerical algorithm has been validated along with, first, the analytical solution of the one-dimensional two-phase planar solidification and, second, the growth of an isolated needle-like dendrite. 4. Results and discussion

this influence at high supercooling degrees, a two-dimensional simulation setup has been appropriately designed (Fig. 9): a sinusoidal function as initial perturbation with amplitude A and wavelength ls is prescribed, with rt and vt representing the radius and velocity at the needle-tip, respectively. An adiabatic condition at the upper boundary of the solution domain, a symmetry-plane condition at the right and left sides (enabling the consideration of possible influence through surrounding needles) and a fixed temperature value, T ¼ 273.15 K, at the lower boundary close the simulation setup. DT describes the initial supercooling in the liquid phase denoting by F > 0. In the solid phase, defined by F < 0, the temperature is assumed to be constant at T ¼ 273.15 K. To assure a correct reproduction of the physics characterizing the system, an appropriately high spatial resolution (a Cartesian grid was used) corresponding to at least thirty control volumes was applied in order to resolve the curvature radius of the needle-tip correctly. The stroboscopic image, Fig. 10, exhibits a series of short samples illustrating the needle growth into the supercooled water; the initial supercooling in the liquid phase prescribed is DT ¼ 15 K. The initial amplitude of the sinusoidal function (ls ¼ 0.8 mm) corresponds to 0.2 mm. Hence, the initial tip radius is rt0 ¼ 8:5,108 m. The time interval between samples corresponds to Dt ¼ 0.25 ms. In Fig. 10 it may be observed that any side-branching effects will be suppressed by the thermal effect of the neighbouring needles (taken into account by applying the symmetry-plane boundary condition at the domain sides), having a significant impact on the surface structure of the growing needle. In Fig. 11 the steady-state tip velocity of the growing needle is plotted against the ratio of the wavelength of the perturbation function, ls, to the fastest-growing-mode wavelength, lf, for a constant value of initial supercooling (DT ¼ 15). The system characteristic time step, Dtc, is defined as the ratio of the steady tip radius of curvature to the tip velocity after a physical time step corresponding to 1,108 s. Hence, the steady-state of the system is reached when the transient deviation of the tip velocity is smaller then 0.1% for a physical time of 10  Dtc. The maximum value of the

The computational results pertinent to the growth of an isolated dendrite displayed in the authors' previous works [7,31] exhibit very good qualitative and quantitative agreement with analytical solutions as well as with the available experiments in the heatdiffusion-dominated region, Fig. 5. However, at supercooling degrees higher than 5 K, one may observe a specific departure of the experimental results from analytic solutions and computations. Possible reasons for this deviation could be sought either in the thermal interaction between growing needles or in the kinetic effects, being both neglected in analytical solutions and computations. An in-depth analysis of both influences is the subject of the following sections. 4.1. Growth of an array of needles In typical icing applications pertinent to high supercooling degrees the needle-crystals arrangement does not resemble a chaotic infinite collective but rather a kind of an array. Starting from an isolated dendritic shape at a low supercooling degree, the crystal shape evolves into the formation of needle arrays due to high destabilizing effects pertinent to higher supercooling degree. The steady-state tip velocity of an isolated needle cannot be applied if an array of needles is investigated, because the thermal boundary layer around the needle-tip should be influenced by the surrounding growing needles. At very low supercooling degrees, investigations reported in Refs. [35e38] show that the thermal interaction between needles has a significant impact on the steadystate tip radius and velocity. In order to investigate computationally

Fig. 9. Two-dimensional simulation setup referring to an array of growing needles; Sinusoidal function as initial perturbation with amplitude A and wavelength ls, rt and vt represent the radius and the velocity, respectively, at the needle-tip. Adiabatic upper boundary condition, symmetry-plane condition at the right and left side and fixed value condition (T ¼ 273.15 K) at lower boundary. DT describes the initial supercooling in the liquid phase (F > 0). In the solid phase (F < 0), temperature is assumed to be constant at T ¼ 273.15 K.

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

Fig. 10. Stroboscopic image of an ice needle growing into supercooled water; the needles displayed on the sides, representing the result of identical simulation, should illustrate a growing array. The initial supercooling in the liquid phase values 15 K, the initial amplitude of the sinusoidal function (with a wavelength of l ¼ 0.8 mm) amounts to 0.2 mm. Hence, the initial tip radius is rt0 ¼ 8:5,108 m. Time between samples values Dt ¼ 0.25 ms.

tip velocity (z17.5 cm/s) corresponding to the fastest-growingmode (lps/ffiffiffilf ¼ 1; we recall the onset of the unstable growth at ls ¼ lf = 3) agrees well with the analytical value being pertinent to the growth of an isolated needle (compare it with Fig. 5); however, it deviates substantially from the experimentally determined value (z6.6 cm/s) which has been determined for an array of needles according to Shibkov et al. (2005, [34]). This motivated finally the presently considered computational setup accounting also for a needle array. For higher ratios ls/lf (here the values 4, 9, 18, 27, 36 and 45 are adopted) a gradual departure from the theoretical value is obvious. Increasing systematically the wavelength of initial perturbation ls relative to that of the fastest growing mode lf causes a substantial decrease of the steady-state tip velocity of the dendrite. It means that the growth of a dendrite characterized by

Fig. 11. Steady-state tip velocity, vt, of a needle growing into supercooled water (DT ¼ 15 K) as a function of ratio between the wavelength of the initial perturbation function, ls, and the fastest-growing-mode wavelength, lf. The maximum tip velocity value at the fastest-growing-mode agrees well with the analytical value (growth of an isolated needle) but departs substantially from relevant experimental value (z6.6 cm/ s). At high supercooling degrees the thermal interaction between needles in an array is negligible. Hence, the freezing process in the first stage can be simplified as growth of an isolated needle-like dendrite.

157

such an enhanced initial perturbation is appropriately decelerated through the neighbouring, faster growing dendrites (whose wavelengths correspond to those of respective fastest-growing mode). It means furthermore that the tip velocity of a dendrite growing in accordance with the relevant fastest growing mode, corresponding to a theoretical velocity of an isolated needle, is not noticeably affected by a possible thermal influence through surrounding dendrites. Consequently, a conclusion can be drawn that at high supercooling degrees the thermal interaction between needles in an array is negligibly small compared to the thermal heat diffusion in the liquid phase due to supercooling. Accordingly, the kinetic effects resembling a specific rearrangement of molecules at the iceeliquid interface are recognized to be the most probable reason for the departure of the experimental results (of Shibkov et al., 2005, [34]) from those pertinent to the marginal stability theory (see the following section). It should also be emphasized that the initial amplitude, A, of the sinusoidal perturbation function has an insignificant influence on the steady-state tip velocity as can be observed in Fig. 12. In this figure the computational results pertinent to the case ls/lf ¼ 36 are displayed. The resulting steady state tip velocity follows closely the value corresponding to vt ¼ 1 cm/s (in accordance with Fig. 12). 4.2. Growth of a needle-like dendrite including kinetic effects Shibkov et al. (2005, [34]) pointed out that considering the effects pertinent to the kinetics-limited dendrite growth is of importance in the higher supercooling range (DT > 5). Unfortunately, there are no elaborated studies in the published literature about how such a phenomenon could be quantified in the framework of a mathematical model. Starting point of the present modelling activity is the work of Davis (2001, [8]) who suggested the velocity of the solideliquid interface at high supercooling to be approximated by a linear function of a fraction of the total undercooling originating from molecular kinetics, DTkin, Eq. (5). This approximation is valid for pure substances with low latent heat of crystallization such as pure metals and water. As mentioned above, the intensity of a coefficient accounting for the kinetic effects originating from the molecular rearrangement at the interface from liquid state to a solid state (denoted as kinetic coefficient kkin) is, to date, not well defined in the scientific literature. In the following, an empirical model for the kinetic coefficient is outlined. It represents a revised version of the model derived in Ref. [5]. In order to

Fig. 12. Steady-state tip velocity, vt, of a needle growing into supercooled water (DT ¼ 15 K) as a function of ratio between the initial amplitude, A, and wavelength, ls, of perturbation function (ls/lf ¼ 36). Initial amplitude of the sinusoidal perturbation function does not influence the steady-state of the needle-like dendrite.

158

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

quantify the model term mimicking the kinetics-limited growth the theoretical and experimental results displayed in Fig. 5 are comparatively analyzed. First of all, it is assumed that the steadystate dendrite tip velocity determined experimentally in [34], vn,exp, is directly proportional to

vn;exp fDTT :

(19)

with DTT representing the total undercooling (consisting of viscous and kinetic fractions, DTn and DTkin respectively) at the solideliquid interface, which is defined as

DTT ¼ DTn þ DTkin :

(20)

In the analytical MST solution, the surface tension and hence, the viscous undercooling is accounted for. In order to account for the kinetic undercooling, the ratio of the theoretical steady-state velocity (related to the viscous undercooling only) to the experimentally determined   one (influenced in addition also by the kinetic undercooling), Dvn  ¼ vMST =vn;exp , is computed assuming its proportionality to the ratio of the viscous (capillary) undercooling to the total one

DTT ¼

DTn : jDvn j

(21)

Inserting it into Eq. (20) and adopting a linear function for the kinetic undercooling, Eq. (5), following expression for the kinetic coefficient is obtained

kkin ¼

jDvn jvMST Lr : ð1  jDvn jÞksTm

(22)

The preliminary results obtained at high supercooling degrees by accounting for the kinetic undercooling in the present computational model exhibit good agreement with the experimental data. The only limitation when using Eq. (22) is that the kinetic coefficient is directly dependent on the theoretical steady-state velocity, hence, on the initial-supercooling degree. In order to find a coefficient value valid for all supercooling degrees considered, the relation between the kinetic coefficient and the theoretical steadystate velocity is introduced by applying the following relation: 2

kkin ¼ xðvMST Þ3 ;

(23)

Fig. 13. Schematic of the computational domain and mesh detail in the fine resolution region. ld represents the initial height of the parabolic needle.

the fine-resolution region. In order to assure a correct curvature reproduction, the curvature radius of the dendrite tip is resolved by approximately twenty grid cells. In order to fulfil the conditions at infinity, adiabatic boundary conditions for the temperature field have to be prescribed at the left, right and upper domain boundary, Fig. 13. After accounting for the kinetic effects by considering the present approach the computational results at high supercooling degrees follow closely the experimental results, Fig. 14. 4.3. Modelling the initial and final stages of the supercooled water freezing mechanism Schematics displayed in Fig. 15 along with the experimental images illustrate the mechanism of the supercooled water freezing pertinent to both characteristic freezing stages. The difference between the first and the second stage is basically caused by the direction of the heat flux. During the first stage, the latent heat of solidification is mainly added to the liquid due to high supercooling. Regardless of the level of the initial supercooling in the water, the incipient structure of the crystallization front is a smooth planar layer, Fig. 15a, at a constant temperature value corresponding to Tm, whose onset coincides with the substrate contact area (in

The value of the coefficient x representing a pre-exponential factor, determined, i.e. calibrated by relevance to the reference experiment by Shibkov et al. (see Fig. 5), amounts finally (p/11)2/3. Accordingly, the coefficient x represents a dimensional quantity whose unit corresponds to m1/3/(s1/3K). Under consideration of Eq. (5), the kinetic undercooling can be formulated as follows:

DTkin ¼

vn : kkin

(24)

Assuming that the solideliquid interface velocity, vn, corresponds to the theoretical value vMST obtained by neglecting the kinetic effects, the kinetic undercooling can be redefined resulting in following equation: 1

DTkin ¼

ðvMST Þ3 : x

(25)

Herewith, the appropriate quantification of the kinetic undercooling influence on the dendrite growth in the high supercooling range is provided. Corresponding computational validation is performed by adopting a solution domain whose schematic is depicted in Fig. 13. The parabolic needle-like dendrite is initialized within

Fig. 14. Steady state tip velocity, vt, dependence on initial supercooling DT for ls/lf ¼ 1: computational vs. MST-theoretical and experimental results; the results obtained by present empirical model accounting for kinetic undercooling at the solideliquid interface exhibit good agreement with reference experiments over the entire range of initial supercooling.

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

Fig. 15. Crystallization stages of supercooled pure water: in the first stage (left) the initial planar solidification front becomes unstable due to a supercooling degree of 15 K, small bumps at the interface grow to needles. A small proportion of the drop freezes instantaneously, not more than enough to raise the temperature to the crystallization temperature. In the second stage the cold substrate cools the drop and hence the water between the needles. It freezes smoothly (planar solidification).

accordance with heterogeneous nucleation, [13]). The planar front develops small bumps, Fig. 15b, due to the destabilizing effect of supercooling. These bumps at the interface start to grow into the supercooled water. The instability mechanism depends on the thermal boundary layer on the liquid side of the system. The more heat is conducted into the liquid the more unstable becomes the shape of the solideliquid interface. Or in other words, the smaller the cutoff wavelength, the largest possible wavelength for a stable interface. It corresponds to the fastest-growing-mode in line with Eq. (8). The heat conduction in the solid, which would expectedly stabilize the solidification front, is negligibly small at this stage of freezing. A small bump at the interface propagates further into the liquid, steepening the magnitude of the temperature gradient on the liquid side. Accordingly, the heat flux balance causes the bump to grow faster. During this transition process the small bumps at the solideliquid interface develop into crystals of different shapes. At

159

higher degrees of initial supercooling, DT  5 K, the solidification front progresses into arrays of branchless dendrites/needle-like crystals, Fig. 15c. The thermal interaction between needles does not influence the steady tip velocity, as demonstrated in the previous section. In contrast, kinetic effects have to be accounted for correspondingly; the appropriately formulated model yields results agreeing well with experimental findings. After the first stage the initial supercooling is exhausted. The drop temperature raises up to the crystallization temperature. The ratio of the sensible heat, cvDT, to the latent heat, L, reveals that only a portion (e.g. 18.8% for DT ¼ 15 K) of the drop volume is solidified after the first stage of crystallization. The experiments do not show any indication of dendrites during the first stage of crystallization, Fig. 1. For an initial supercooling of 15 K in the liquid, the critical wavelength for the dendritic growth corresponds to 2.6  108 m. The critical wavelength, lc, represents the largest possible wavelength for stable growth. It means, a perturbation with a smaller wavelength will be suppressed by the interfacial energy, while a perturbation with a larger wavelength will growth very fast into the supercooled liquid (dendritic growth). The wavelength of the fastest growth mode, defined by Eq. (8), represents the characteristic wavelength for the unstable growth. Dendrites exhibit a wavelength of the order of magnitude 108 m, while the drop in the experiments has a volume of 5 ml. By modelling the drop as a truncated sphere, the base diameter can be reasonably estimated and amounts approximately 2 mm. Indeed, the needlelike dendrites are too small to be visualized experimentally. At the beginning of the second stage, the temperature of the drop remains constant at the melting temperature, Tm. Due to an appropriately colder substrate in the experiments, Tsub ¼ 258.15 K, the heat is now removed via heat conduction into the substrate. Fig. 16 illustrates the mechanism of the freezing pertinent to the second stage. The heat is conducted from the solid phase (having higher heat conductivity) to the substrate. After the first stage the solidified water becomes colder and consequently cools down the water between the needles. Thus, the water begins to freeze

Fig. 16. (a) A sketch of the freezing mechanism in the second stage: heat from the solid phase, at Tm ¼ 273.15 K, is conducted to the cold substrate. Water between needle-like dendrites starts to freeze: red dotted lines visualize the possible time evolution of the solidification front. (b) A two-dimensional simulation setup, referring to freezing mechanism in the second stage, is depicted from the sketch; initial temperature in both the solid and liquid phase values Ts ¼ Tl ¼ 273.15 K. The lower boundary condition at the substrate, Tsub ¼ 258.15 K, represents the cooling substrate. Adiabatic upper boundary condition and symmetry-plane condition at the right and left side close the setup. The mesh resolution is defined by Dx, whereas h(t) describes one-dimensional growing solid/ice layer from the bottom/substrate to the top as a function of time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

160

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161

Fig. 17. Growing ice layer plotted against time; simulations converge to the analytical solution of the one-dimensional planar solidification [1]. After 20 s the ice layer reaches a height of 1.997 mm; in experiments [13] the water drop has a volume of 5 ml, assuming a spherical shape it leads to a diameter of 2.11 mm. The gain in volume during the phase change is not accounted in computations.

upwards from the solid cold boundary representing the substrate. The relevant solidification front velocity is considerably lower compared to the freezing velocity pertinent to the first stage. The second stage of freezing can be mathematically described as a onedimensional freezing, starting from a cold boundary, Fig. 15d. To verify this theory, a suitable two-dimensional simulation setup is configured, Fig. 16. The initial temperature in both the solid and liquid phases is Ts ¼ Tl ¼ 273.15 K. The lower domain boundary at Tsub ¼ 258.15 K represents the cooling substrate. An adiabatic condition at the upper domain boundary is prescribed and the symmetry-plane condition was adopted at the right and left domain sides. At this point, it should be mentioned that the freezing process is also influenced by the thermal diffusivity of the substrate affecting the freezing front velocity. Experiments by Jung et al. [14] confirm this behaviour. In the present work, the influence of the thermal diffusivity is neglected; a constant temperature of the substrate is assumed when simulating the second stage of freezing. In the experiments of Jung et al. [13], the volume of the supercooled water drop is 5 ml, i.e. if assuming a perfectly spherical shape of the drop the corresponding diameter would be equal to 2.11 mm. The duration of the second crystallization stage is approx. 20 s, see Fig. 1. Computations with different grid resolutions show, first, that the results converge unambiguously to the analytical solution (lower grid resolutions, implying Dx ¼ 5 and even 2.5 mm, have also been tested in some relevant configurations, see e.g. Refs. [7,31]) and, second, that the height of the growing ice layer after 20 s reaches the value of 1.997 mm, Fig. 17. The gain in volume and, consequently, in velocity of the solidification front during the phase change is negligibly small (only ca. 3% of the front velocity), being as such represented also in the present computations. Therefore, we can conclude that the governing equations describing mathematically the solidification mechanism resembling a heat-diffusionbased solidification model, section 2, can account for both freezing stages: first dendritic stage and second planar stage solidification. 5. Conclusions A computational model describing the solidification of supercooled water, introduced recently in Refs. [7,31], has been extended to account for both stages of crystallization e the first rapid,

recalescent stage and the second slower, quasi-isothermal stage e in conjunction with an appropriate numerical algorithm. The latter relates specifically to the discretization of the heat transfer equation and the temperature gradient determination at the solideliquid interface. Based on the present computational results, the different mechanisms underlying both freezing stages can be explained as follows: in the first stage the initial planar solidification front becomes morphologically unstable due to a high degree of supercooling. Small bumps/instabilities evolving at the interface propagate further into the liquid, causing a relevant steepening of the temperature gradient at the liquid side of the interface, contributing decisively to its rapid growth. During this transition process the small bumps at the solideliquid interface develop into crystals of different shapes. At higher degrees of supercooling, DT  5 K, the instabilities grow into a dense array of needles. A small fraction of the drop freezes instantaneously, not more than enough to raise up the temperature to a value corresponding to the crystallization, i.e. until the thermal energy rate originating from supercooling is exhausted. The crystallization front can be modelled as growing needle-like dendrites influenced by kinetic effects. In the second stage, the drop loses heat by evaporation and conduction. Within a sessile drop, the cold substrate cools the water between the needles and, accordingly, the crystallization front grows up planar from the cold boundary. The velocity of the solidification front is considerably lower than the freezing velocity in the first stage. The second stage can be mathematically described as a one-dimensional solidification process. Acknowledgements This research was supported by the German Scientific Foundation (DFG) in the framework of the SFB-TRR 75 collaborative  Tukovi
research center. Our special thanks go to Dr. Z. c (University of Zagreb) for his support in numerical implementation regarding the open-source software OpenFOAM® and Mr. D. Kintea for useful discussions about modelling kinetic undercooling. References [1] V. Alexiades, A.D. Solomon, Mathematical Modeling of Melting and Freezing Processes, Hemisphere Publishing Corporation, Washington, 1993. [2] M.B. Baker, M. Baker, A new look at homogeneous freezing water, Geophys. Res. Lett. 31 (2004) L19102. [3] S. Bauerecker, P. Ulbig, V. Buch, L. Vrbka, P. Jungwith, Monitoring ice nucleation in pure and salty water via high-speed imaging and computer simulations, J. Phys. Chem. C 112 (2008) 7631e7636. [4] F. Caliskan, C. Hajiyev, A review of in-flight detection and identification of aircraft icing and reconfigurable control, Prog. Aerosp. Sci. 60 (2013) 12e34. [5] A. Criscione, D. Kintea, S. Jakirli
c, I.V. Roisman, C. Tropea, A new approach for water crystallization in the kinetics-limited growth region, in: Proceedings of the 8th International Conference on Multiphase Flow, Jeju, South Korea, 2013. [6] A. Criscione, D. Kintea, Z. Tukovi
c, S. Jakirli
c, I.V. Roisman, C. Tropea, On computational investigation of the supercooled stefan problem, in: Proceedings of the 12th International Conference on Liquid Atomization and Spray Systems, Heidelberg, Germany, 2012. [7] A. Criscione, D. Kintea, Z. Tukovi
c, S. Jakirli
c, I.V. Roisman, C. Tropea, Crystallization of supercooled water: a level-set-based modeling of the dendrite tip velocity, Int. J. Heat. Mass Transf. 66 (2013b) 830e837. [8] S.H. Davis, Theory of Solidification, Cambridge University Press, Cambridge, 2001. [9] M. Do-Quang, G. Amberg, Simulation of free dendritic crystal growth in a gravity environment, J. Comput. Phys. 227 (2008) 1772e1789. [10] C. Godreche (Ed.), Solids Far from Equilibrium, Cambridge University Press, 1991. [11] M. Homola, M. Virk, W.T.,P. Niklasson, P. Sundsbo, Effect of atmospheric temperature and droplet size variation on ice accretion of wind turbine blades, J. Wind Eng. Ind. Aerod 98 (2010) 724e729. [12] G. Ivantsov, Temperature field around spherical, cylindrical, and needleshaped crystals which grow in supercooled melt, Dokl. Akad. Nauk. SSSR 558 (1947) 567e569.

A. Criscione et al. / International Journal of Thermal Sciences 92 (2015) 150e161 [13] S. Jung, M.K. Tiwari, N.V. Doan, D. Poulikakos, Mechanism of supercooled droplet freezing on surfaces, Nat. Commun. 3 (2012) 615, http://dx.doi.org/ 10.1038/ncomms1630. [14] S. Jung, M.K. Tiwari, D. Poulikakos, Frost halos from supercooled water droplets, Proc. Natl. Acad. Sci. U. S. A. 109 (2012) 16073e16078. [15] R. Kobayashi, Modeling and numerical simulations of dendritic crystal growth, Phys. D. Nonlinear Phenom. 63 (1993) 410e423. [16] A. Kraj, E.L. Bibeau, Measurement method and results of ice adhesion force on the curved surface of a wind turbine blade, Renew. Energy 35 (2010a) 741e746. [17] A. Kraj, E.L. Bibeau, Phases of icing on wind turbine blades characterized by ice accumulation, Renew. Energy 35 (2010b) 966e972. [18] J.S. Langer, R.F. Müller-Krumbhaar, Theory of dendritic growth-I. Elements of a stability analysis, Acta Metall. 26 (1978a) 1681e1687. [19] J.S. Langer, R.F. Müller-Krumbhaar, Theory of dendritic growth-II. Instabilities in the limit of vanishing surface tension, Acta Metall. 26 (1978b) 1689e1695. [20] J.S. Langer, R.F. Müller-Krumbhaar, Theory of dendritic growth-III. Effects of surface tension, Acta Metall. 26 (1978c) 1697e1708. [21] J.S. Langer, R.F. Sekerka, T. Fujioka, Evidence for a universal law of dendritic growth rates, J. Cryst. Growth 44 (1978) 414e418. [22] G. Lock, The Growth and Decay of Ice, Cambridge University Press, Cambridge, 2005. [23] J. Melody, T. Basar, W.R. Perkins, P. Voulgaris, Parameter identification for inflight detection and characterization of aircraft icing, Control Eng. Pract. 8 (2000) 985e1001. [24] W.W. Mullins, R.F. Sekerka, Morphological stability of a particle growing by diffusion or heat flow, J. Appl. Phys. 34 (1963) 323e329. [25] W.W. Mullins, R.F. Sekerka, The stability of a planar interface during solidification of a dilute binary alloy, J. Appl. Phys. 35 (1964) 444e451. [26] K. Ohsaka, E.H. Trinh, Apparatus for measuring the growth velocity of dendritic ice in undercooled water, J. Cryst. Growth 194 (1998) 138e142. [27] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988) 12e49.

161

[28] K. Petty, F. D.J, A Statistical Review of Aviation Airframe Icing Accidents in U.S., National Transportation Safety Board, Washington D.C., 2010. [29] M. Politovich, Aircraft icing, Enc. Atmos. Sci. (2003) 68e75. [30] H. Pruppacher, J. Klett, Microphysics of Clouds and Precipitations, second ed., Kluwer Academic Publishers, Netherlands, 1997. [31] P. Rauschenberger, A. Criscione, K. Eisenschmidt, D. Kintea, S. Jakirli
c, Z. Tukovi
c, I.V. Roisman, B. Weigand, C. Tropea, Comparative assessment of volume-of-fluid and level-set methods by relevance to dendritic ice growth in supercooled water, Comput. Fluids 79 (2013) 44. [32] A.A. Shibkov, Y.I. Golovin, M.A. Zheltov, A.A. Korolev, A.A. Leonov, Morphology diagram of nonequilibrium patterns of ice crystals growing in supercooled water, Phys. A 319 (2003) 65e79. [33] A.A. Shibkov, Y.I. Golovin, M.A. Zheltov, A.A. Korolev, A.A. Vlasov, Kinetics and morphology of nonequilibrium growth of ice in supercooled water, Crystallogr. Rep. 46 (2001) 496e502. [34] A.A. Shibkov, M.A. Zheltov, A.A. Korolev, A.A. Kazakov, A.A. Leonov, Crossover from diffusion-limited to kinetics-limited growth of ice crystals, J. Cryst. Growth 285 (2005) 215e227. [35] B. Spencer, H. Huppert, Steady-state solutions for an array of stronglyinteracting needle crystals in the limit of small undercooling, J. Cryst. Growth 148 (1995) 305e323. [36] B. Spencer, H. Huppert, On the solidification of dendritic arrays: an asymptotic theory for the directional solidification of slender needle crystals, Acta Mater. 45 (1997) 1535e1549. [37] B. Spencer, H. Huppert, On the solidification of dendritic arrays: selection of the tip characteristics of slender needle crystals by array interactions, Acta Mater. 46 (1998) 2645e2662. [38] B. Spencer, H. Huppert, The relation between dendrite tip characteristics and dendrite spacings in alloy directional solidification, J. Cryst. Growth 200 (1999) 287e296. [39] A.A. Wheeler, B.T. Murray, R. Schaefer, Computation of dendrites using a phase field model, Phys. D. 66 (1993) 243e262.