Perturbation solutions for the finite radially symmetric Stefan problem

International Journal of Thermal Sciences 104 (2016) 386e395

Contents lists available at ScienceDirect

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

Perturbation solutions for the finite radially symmetric Stefan problem Florica Ioana Dragomirescu a, *, Kathrin Eisenschmidt b, Christian Rohde a, Bernhard Weigand b a b

University of Stuttgart, Institute of Applied Analysis and Numerical Simulation, Pfaffenwaldring 57, 70569 Stuttgart, Germany University of Stuttgart, Institute of Aerospace Thermodynamics, Pfaffenwaldring 31, 70569 Stuttgart, Germany

a r t i c l e i n f o

a b s t r a c t

Article history: Received 15 December 2014 Received in revised form 19 January 2016 Accepted 20 January 2016 Available online 19 April 2016

The beginning of the ice-growth from a nucleus in the freezing process of a single supercooled droplet is modelled through a two-phase Stefan problem which is characterized by an unknown boundary separating the two distinct phases. A semi-analytical solution for the iceewater interface is presented using a perturbation method followed by an iterative procedure. The method offers a quantitatively good approximation of the interface position as well as the temperature distribution in the ice and in the water phase in comparison with the considered reference numerical results. Furthermore, it is less expensive in terms of time and memory resources. The results obtained for the initial radially symmetric ice growth of the small ice particle can be used as subgrid models in a direct numerical simulation of the dendritic ice growth in a supercooled droplet. © 2016 Elsevier Masson SAS. All rights reserved.

Keywords: Stefan problem Perturbation solution Supercooled droplet

1. Introduction Understanding the process of ice formation in supercooled water is a problem of significant importance due to its common presence in our environment and its importance in technical areas. Water is supercooled if its temperature is lower than the melting point, i.e. Tm ¼ 273.15 K, under atmospheric pressure. The supercooled liquid is in a metastable state, which implies that a small disturbance can initiate a rapid and unstable solidification process. Predicting this process for different environmental conditions allows a better understanding of its microphysics and will help to improve numerical models of precipitation dynamics. According to [1] the freezing process of a supercooled droplet can be described in four stages: a supercooling stage, a recalescence stage, a freezing stage and a tempering stage. The whole freezing process may take tens of seconds [1], which, of course, depends on the ambient conditions. During the recalescence stage supercooling drives rapid crystal growth from the nuclei until the droplet has reached its equilibrium freezing temperature. The duration of the recalescence stage is extremely short compared to the other three stages. In models which describe this stage, ice starts from the inside of the droplet. The subsequent freezing stage is determined by heat conduction and it ends when the droplet is completely frozen.

* Corresponding author. E-mail address: [email protected] (F.I. Dragomirescu). http://dx.doi.org/10.1016/j.ijthermalsci.2016.01.019 1290-0729/© 2016 Elsevier Masson SAS. All rights reserved.

This stage is typically modelled with ice starting from the outside. In the present paper we focus on the beginning of the recalescence stage. A supercooled droplet may contain several nuclei. However, the present model is based on the assumption of radial symmetry and is, therefore, limited to one nucleus at the droplet center. In this work we determine the growth rate of the ice particle from the initial nucleus up to a stable radius [2] (see Eq. (2)) as well as the corresponding temperature distribution in the ice and in the water phase. The duration of the ice particle growth up to the same radius has been evaluated from numerical simulations and is shown in Table 1. We denote by tstable the dimensionless time that is needed for the interface to reach this stable radius. The evolution of the iceewater interface is governed by a Stefan problem which describes a time-dependent heat conduction problem involving a liquidesolid phase change. It is characterized by an unknown boundary separating the two distinct phases. Thus, the interface must be considered as a part of the solution. One of the parameters characterizing the freezing process is the Stefan number defined throughout the paper as the ratio of the sensible heat to the latent heat, i.e.

St ¼ cw ðTm  T∞ Þ=L;

(1)

where cw represents the specific heat capacity, L the specific latent heat, Tm and T∞ the melting and the ambient temperature, respectively.

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395 Table 1 dim Values of dimensional time tstable and non-dimensional tstable as needed for the interface to reach Rdim (see Eq. (2)) and Rstable, respectively. stable T∞/[K]

dim =½s tstable

Rdim =½m stable

tstable/[e]

Rstable/[e]

243.15 248.15 258.15 263.15

4.9834e-09 1.0989e-08 7.0686e-08 2.6814e-07

4.2659e-08 5.0386e-08 8.1749e-08 1.2125e-07

0.6603e-05 0.1456e-04 0.9366e-04 0.3680e-03

0.0043 0.0050 0.0082 0.0121

387

in the mathematical modelling is described in detail in Section 3. Perturbation solutions for the freezing behaviour for two temperatures (high and low supercooling) (see Table 2) are presented in comparison with numerical results in Section 4. In order to validate the method, two additional temperatures have been considered. Section 5 summarizes the results.

2. Mathematical model for radially symmetric ice growth Much work has been conducted with the purpose of obtaining analytical or semi-analytical approaches as well as numerical approximations for this type of problems (e.g. [3e9]). Exact solutions for classical two-phase Stefan problems can be obtained only in very special cases, e.g. for infinite domains. This is due to the inherent nonlinearity introduced by the condition on the moving interface. Classical similarity solutions of moving front problems, going back to the original work of Stefan, can be found for very simple geometry and boundary conditions in e.g. [10]. McCue et al. [11] also used a similarity solution for a small-time limit of a melting process, which, however, could not satisfy all boundary conditions. Chadam et al. [12] and Herrero et al. [13] gave solutions for one-phase problems. In this work we investigate the influence of heat conduction in the ice phase on its growth. Font and Wu [8,14,15], who investigated melting problems, presented solutions that are restricted to small Stefan numbers. Weinbaum et al. [16] and Savovic et al. [17] also presented solutions for small Stefan numbers only. Due to the high supercooling, which is a requirement for the homogeneous nucleation and freezing that we want to investigate, the Stefan numbers are much higher in our case. Following [14,16,18], the coupled partial differential equations for the water and the ice phase can be linearized by introducing a perturbation parameter ε ¼ Stn. A solution using ε ¼ St1/3 as perturbation parameter has been constructed in [3] with a modified Stefan condition. This solution closely matches the numerical results, also presented in [3], but deteriorates for higher Stefan numbers. The present paper extends this work solving the problem for the whole range of realistic Stefan numbers. Our purpose is to obtain the temperature distributions in both the water and the ice phase as well as to track the location of the interface between the ice and the water using inexpensive methods in terms of time and computer resources. A semi-analytical approach based on Perturbation Techniques followed by an Iterative Procedure (denoted by PT-IP in the following) is provided. Numerical results obtained using the algorithm described in [3] are considered here as Reference Numerical Solution (RNS). The reference numerical results have been validated for a one-dimensional phase change problem with heat conduction in the ice and in the water phase against an analytical solution [19]. The outline of this paper is as follows: In the next section the mathematical formulation of the problem for the radially symmetric case is given. The solution of the free boundary problem is defined not only by the temperatures in the ice and water phase, but also by the interface position. The perturbation approach used

In this section, we present the mathematical model describing the radially symmetric freezing of a single droplet in a supercooled environment. The assumption of radial symmetry holds for the initial growth of a small nucleus. This small nucleus of critical size Rc is of the order 1e-9m and it starts to grow spherically inside a supercooled droplet up to the radius that limits the stable growth (Fig. 1). This non-dimensional stable radius Rstable is defined by Rstable ¼ Rdim =Rd with Rdim estimated in [2] as stable stable

Rdim stable ¼



 2s 4k Tm  T∞ 7þ i : Lr kw Tm

Here s represents the surface tension, r the density and ki and kw the heat conductivity of ice and water, respectively. The water droplet is supposed to be spherical with a given radius Rd ¼ 1e5 m. After exceeding Rstable instabilities are not damped any more and dendrities might form. We consider a geometrical setting as shown in Fig. 1. The model we present in the following can serve as a subgrid model for direct numerical simulations of the dendritic ice particle growth inside a supercooled water droplet. The movement of the iceewater interface is described by the Stefan condition which relates the release of latent heat at the interface to the heat flux into the water and into the ice phase. This determines the interface velocity. The temperature fields in the ice, Ti, and in the water phase, Tw, are described by the heat conduction equations for both regions

! vTi 2 vTi v2 Ti ¼ ai 0 0 þ 02 ; 0  r 0  R0 ðt 0 Þ; r vr vt 0 vr

Table 2 Selected values of parameters for different ambient temperatures. T∞/[K]

St

a/[e]

Rc Rd =½

TGT/[K]

qffiffiffiffiffi ai aw =½

243.15 248.15 258.15 263.15

0.4402 0.3444 0.1960 0.1292

1.683e-4 2.032e-04 3.428e-04 5.172e-4

2.197e-4 2.7575e-04 4.6164e-04 6.050e-4

250.1677 254.7206 262.0103 264.6013

3.2553 3.1065 2.9440 2.8927

(2)

Fig. 1. Geometry of the droplet formation.

(3)

388

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395

! vTw 2 vTw v2 Tw ¼ aw 0 þ 02 ; R0 ðt 0 Þ  r 0  Rd ; r vr 0 vt 0 vr

(4)

where ai and aw are the thermal diffusivities of ice and water. The lowering of the melting temperature due to the curvature of the ice particle's interface, which is described by the GibbseThomson effect, has to be taken into account. Thus, the temperature at the interface in the radial symmetric case reads [2]


TGT ðR0 ðt 0 ÞÞ ¼ Tm 1 

 2s : 0 0 rLR ðt Þ

(5)

The boundary condition at the interface of the solid and the liquid phase is written as

Tw ðt 0 ; R0 ðt 0 ÞÞ ¼ Ti ðt 0 ; R0 ðt 0 ÞÞ ¼ TGT ðR0 ðt 0 ÞÞ:

(6)

A lot of publications investigating non-planar problems (like [6,17,18,20,21]) do not account for the GibbseThomson effect. Since the size of the nucleus that initiates the homogeneous freezing is of the order of nanometers, the reduction of the melting temperature plays an important role in our case. King et al. [20], Back et al. [9,22], account for the non-equilibrium effect of kinetic undercooling by a generalized GibbseThomson equation. This prevents blow-up of the solution, e.g. in the case of inward melting where a singularity appears in the center, but is not necessary in our case as the freezing starts from a nucleus with a finite size. Font et al. [14] derived a general GibbseThomson equation from thermodynamic principles. However, they concluded from their investigations that the process (melting in their case) is not dominated by the differences in density or heat capacity of the two phases. According to their estimation of the influence of the different terms, we see that in our case, the terms containing the differences in heat capacity and density are much smaller than the dominant term which contains surface tension, curvature, and the latent heat of fusion. The influence of the non-dominant terms even decreases quickly with time. Therefore, we use the standard GibbseThomson Eqs. (5) and (6). Corresponding to the standard GibbseThomson condition, we use the standard Stefan condition which neglects the changes in heat capacities and densities between both phases [6,16]. The Stefan condition determines the interface velocity and is given by

rL

dR0 vTw  vT  ¼ kw 0  0 0 þ ki 0i  0 0 : dt 0 vr R ðt Þ vr R ðt Þ

(7)

Herein r denotes the density, which is assumed to be the same for both phases. The temperature at the watereair interface is set to the ambient temperature

Tw ðt 0 ; Rd Þ ¼ T∞ :

(8)

For r 0 ¼ 0 the temperature of the ice must be finite: 0  Ti ðt 0 ; 0Þ < ∞. The nucleation rate is evaluated for a constant temperature, i.e. the ambient temperature. Therefore, the initial condition for the temperature in both phases reads

Tw ð0; r 0 Þ ¼ Ti ð0; r 0 Þ ¼ T∞ :

t ¼ t1 ; x ¼

r r  RðtÞ ; h¼ ; RðtÞ 1  RðtÞ

(10)

the resulting set of equations describing the evolution of the modified temperature fields in the ice and the water has the following form [3].

x

dR vv vv a v2 v R þ R2 ¼ i ; dt vx vt aw vx2

0  x  1;

dR vu vu v2 u ð1  RÞðh  1Þ þ ð1  RÞ2 ¼ 2 ; dt vh vt vh

(11)

0  h  1;

(12)

with the boundary and initial conditions

uðt; 0Þ ¼ vðt; 1Þ ¼ a; uðt; 1Þ ¼ 1; vðt; 0Þ ¼ 0;

(13)


Rc Rc Rc þ ; vð0; xÞ ¼ x ; uð0; hÞ ¼ h 1  Rd Rd Rd

(14)

where

a¼

2sTm : ðTm  T∞ ÞrLRd

(15)

The Stefan condition (7) transforms accordingly into


 dR 1 vu  uð0Þ ki 1 vv  vð1Þ ¼ St ;    dt Rð1  RÞ vh 0 kw R2 vx 1 R2 R2 Rc Rð0Þ ¼ ; Rd (16) with St defined by (1). Values for a and the GibbseThomson temperature are given in Table 2. The ambient temperatures T∞ ¼ 243.15 K and T∞ ¼ 263.15 K have been chosen to validate our method PT-IP (see Section 4 below) for high and low supercooling, respectively. Both ambient conditions are realistic conditions in the atmosphere. Two other temperature values are selected for the sake of completeness. For small times the motion of the ice front in Eq. (11) follows the relation Rft1/2. Assuming R(t)~t1/2, i.e. dR2/dt is a constant, Eq. (11) can be solved analytically. By the procedure from [23], we get the following form for v(t,x)


vðt; xÞ ¼ TðtÞ$XðxÞ ¼ C$ð1 þ b$tÞl $M

 1l b ; 2;  x2 ; 2 2A

where M(a,b,x) represents the confluent hypergeometric function defined as [24].

a aða þ 1Þ x2 aða þ 1Þða þ 2Þ x3 þ þ… Mða; b; xÞ ¼ 1 þ x þ b bðb þ 1Þ 2! bðb þ 1Þðb þ 2Þ 3! and l is the separation constant, A ¼ ai/aw, C and b are constants depending on the initial conditions. However, this solution does not fulfil all the boundary conditions.

(9)

The following dimensionless quantities transform the problem into a planar one [10]: Qw=i ¼ ðTw=i  Tm Þ=ðT∞  Tm Þ, r ¼ r 0 =Rd , RðtÞ ¼ R0 ðt 0 Þ=Rd , u ¼ Qw r, v ¼ Qi r, t1 ¼ aw t 0 =R2d . After the transformation to fixed boundaries using the linear mappings

3. Perturbation solution The dependence on a small parameter given by the Stefan number, St (see Table 2), indicates that perturbation methods can

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395

be used to linearize the partial differential equations and to obtain an approximate analytical solution of this problem [25]. The perturbation technique herein relies on an expansion representation of the solution in terms of an expansion parameter different from the Stefan number. Then the Perturbation Technique Solution (PTS), which is written in an analytical form and satisfies the initial conditions exactly and Eqs. (11) and (12) approximatively, is used as a starting point in an iterative procedure in order to obtain a better agreement with the RNS. The asymptotic behaviour of a singularly perturbed two-phase Stefan problem due to slow diffusion in one of the two-phases was investigated in [26]. Their results emphasized the necessity to introduce a corrected interface condition which accounts for the boundary layer at the moving interface. Mitchell and O'Brien [7] also presented a full asymptotic analysis for an 1D Stefan problem describing solvent diffusion in glassy polymers, providing small and large asymptotic solution for the problem. A small and largertime behaviour of the solution for a Stefan problem with a constant heat source at the moving boundary is discussed in [27]. The significant difference to our previous work [3] is in fact the Stefan number. If it is small, the supercooling is low. So we are not able to describe a process which starts by homogeneous nucleation. It is convenient to rename the (small) Stefan number by the perturbation parameter ε as ε ¼ St and write the Stefan condition (16) in the form Pε[R,u,v] ¼ 0 with

P ε ½R; u; v :

 dR vu   ε R 0  að1  RÞ dt vh
 ki vv   ð1  RÞ  a : vx 1 kw

¼ R2 ð1  RÞ

X ∞ 2 Rc ð1Þn sinðnpxÞeM 1 ðnpÞ t a Rd np n¼1

(17)

(18)

and




X ∞ 2 Rc 1 sinðnphÞeM 2 ðnpÞ t ; uðt; hÞ ¼ ð1  aÞh þ a þ 2 a np Rd n¼1 (19) with M i , i ¼ 1,2, defined by

.



. aw ðRc =Rd Þ2 ; M 2 ¼ 1 ð1  ðRc =Rd ÞÞ2 :

(20)

Inserting the perturbation ansatz



Rc Rðt; εÞ ¼ þ εR1 ðtÞ þ ε2 R2 ðtÞ þ O ε3 ; Rd

(22)



uðt; h; εÞ ¼ u0 ðt; hÞ þ εu1 ðt; hÞ þ ε2 u2 ðt; hÞ þ O ε3

(23)

into Eqs. (11)e(16) and equating terms of the same order, we obtain a divergent asymptotic expansion of the solution R(t), qualitatively different from the RNS [3]. The first order expansion behaviour is linear, which contradicts the validated reference numerical solution (Fig. 2a). For a second order approximation the numerical values become extremely large and do not offer a valid solution (Fig. 2b). Hence, we introduce a new expansion parameter εsSt. We start by considering the temperatures defined by Eqs. (18) and (19) as initial approximations for the temperature in the water and the ice phase. We can motivate this choice by the good agreement of the modified temperatures in the ice and the water phase with the RNS values at tstable. As it is visible in Figs. 3 and 4 the non-dimensional temperature distributions are almost identical. The approximation of the temperature in the water phase described by Eq. (19) is the same as the 0th-order solution for the short-time scale case obtained in [3]. Moreover, the expressions of v and u represent the 0th-order approximations, v0(t,x) and u0(t,h), in Eqs. (22) and (23), respectively. To construct an asymptotic expansion we assume 2 X



εj Rj ðtÞ þ O ε3 ;

(24)

j¼0




M 1 ¼ ai



vðt; x; εÞ ¼ v0 ðt; xÞ þ εv1 ðt; xÞ þ ε2 v2 ðt; xÞ þ O ε3 ;

Rðt; εÞ ¼

The problem Pε[R,u,v] ¼ 0 depends on the small, positive parameter ε. The aim of the perturbation method will be to determine the behaviour of the solution R(t) with respect to ε for t2(0,tstable). Let us start with the unperturbed problem P0[R,u,v] ¼ 0 for ε ¼ 0. In case of a regular perturbation problem for small, non-zero values of ε, one typically obtains a convergent expansion of the solution with respect to ε consisting of the unperturbed solution and higher order corrections. The Stefan condition reduces in this case to dR/ dt ¼ 0. Taking into account the initial condition in Eq. (16), we get R(t) ¼ Rc/Rd. The analytical solutions of Eqs. (11) and (12), with the corresponding boundary conditions, take the form

vðt; xÞ ¼ ax þ 2

389

(21)

which is valid for a new expansion parameter εsε ¼ St with R0 ðt ¼ 0Þs0. We consider for some d>0 the transformation ε ¼ St1ε given by

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ε ¼ ln 1 þ ðd  StÞ :

(25)

From Eq. (25) we obtain St ¼ d  ðe2ε  1Þ and we write





¼: St þ O ε3 : St ¼ d  2ε þ 2ε2 þ O ε3

(26)

The Stefan condition (16) transforms in terms of St into

dRðt; εÞ vu  ¼ StðRðt; εÞ 0  ð1  Rðt; εÞÞuð0Þ dt vh  k vv  i ðð1  Rðt; εÞÞ 1  ð1  Rðt; εÞÞvð1ÞÞÞ þ O ðε3 Þ: vx kw (27) Rðt; εÞ2 ð1  Rðt; εÞÞ

such that ε < ε and  The value of the parameter d is chosen St  St  < tol, with the tolerance tol ¼ 2$102. This gives

.

. d ¼ 2 104 Rc Rd :

(28)

In Table 3 values of ε for the considered ambient temperatures are provided. With the assumptions in Eqs. (18) and (19) and M i , i ¼ 1,2 defined in Eq. (20), the corresponding temperature gradients in the water and the ice, respectively, evaluated at h ¼ 0 and x ¼ 1, take the form


X ∞ 2 vu  Rc ¼ 1  a þ 2  a eM 2 ðnpÞ t ; 0 vh Rd n¼1

(29)

390

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395

Fig. 2. Behaviour of R(t) for the asymptotic expansion (21)e(23) (first and second order approximations).

Fig. 3. Modified temperature distributions in the ice phase (a) and water phase (b), respectively, for T∞ ¼ 243.15 K at tstable ¼ 0.66039100e05.

Fig. 4. Modified temperature distributions in the ice phase (a) and water phase (b), respectively, for T∞ ¼ 263.15 K at tstable ¼ 0.3680e03.

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395 Table 3 Numerical values of the parameter ε for various values of the ambient temperature T∞ and d defined by Eq. (28). T∞/[K]

  St  St =½

ε=½

243.15 248.15 258.15 263.15

0.0105 0.0061 0.0017 0.0011

0.192 0.1614 0.1065 0.0917

Z~t Zt Ri ðtÞ ¼ 



 Gi ~t e

Fð~sÞd~s

0

391

Zt d~t e 0

Fð~sÞd~s ; Ri ðt ¼ 0Þ ¼ 0; i  1;

0

(35) where


X ∞ 2 vv  Rc ¼aþ2 a eM 1 ðnpÞ t : 1 vx Rd n¼1

(30)

We want to simplify these expressions in order to be able to give an explicit analytical solution. In both Eqs. (29) and (30), series of expðM i ðnpÞ2 tÞ, i ¼ 1,2, occur. The numerical values for the  vv  will rapidly parameter M 1 are very large such that the series in vx 1 tend to zero. This quasi-steady state assumption is confirmed by numerical calculations of the temperature field in the ice phase. The quasi-steady behaviour of the modified temperatures u and v was also observed in [3] by transforming the time scale (t ¼ t,St). The parameter M 2 is not small, having values around 1. For tztstable, we need a larger number of terms in order to obtain an accurate approximation of the series. However, an estimation of the integral occurring in terms of the infinite series shows an error of only 1.5% up to the time tstable [3]. Thus, we ignore this term and the Stefan condition simplifies to



 dRðt; εÞ  ¼ d  2ε  2ε2  O ε3 Rðt; εÞ2 1  Rðt; εÞ Rðt; εÞ  a : dt (31) Inserting the perturbation ansatz (24) into Eq. (31) and equating terms of the same order, the 0th-order boundary value problem is found to be

ε0 :

2 R0 1

 R0

 Rc ¼ d R0  a ; R0 ðt ¼ 0Þ ¼ dt Rd

dR0

dR 2 2 0 þ k; gðtÞ ¼ 2R0  3R0 ; f ðtÞ ¼ 3R0  2R0 dt

(36)

 g1 ðtÞ ¼ 2 R0  a ; hðtÞ ¼ 1  3R0 ;

(37)

g2 ðtÞ ¼ R1

dR1 2 dR gðtÞ þ R1 0 hðtÞ þ g1 ðtÞ þ 2R1 : dt dt

The procedure of successively determining Ri ðtÞ, i  1, with Ri ðt ¼ 0Þ ¼ 0 can be applied to any order, always resulting in a linear first order ordinary differential equation. 4. An iterative procedure The perturbation analysis gives us the starting solution (PTS) that will be used in the iterative numerical scheme. The approximate position RðtÞ of the iceewater interface in the PTS is used as an initial guess to solve the heat conduction Eqs. (11) and (12) and to obtain the modified temperature distributions v and u. The interface position is then recalculated via the Stefan condition. The initial boundary value problem defined by the Stefan condition is solved numerically using the Runge-Kutta-Fehlberg method by making use of the symbolic processing capabilities of Maple [28] to quickly code the differential equations. Iterations stop once a predetermined level of tolerance is reached, i.e.

eMax  etol ¼ 2$104 (32)

(38)

(39)

with

    max Rkþ1 ðtÞ  Rk ðtÞ:

and, for i  1

eMax ¼

dRi dR 2  f ðtÞRi þ gi ðtÞ ¼ 0⇔ i  FðtÞRi þ Gi ðtÞ ¼ 0; εi : R0 1  R0 dt dt (33)

By Rk we denote the solution obtained at iteration k. The case k ¼ 0 corresponds to the PTS solution R. The iterative scheme is summarized in Algorithm 1 below.

t2ð0;tstable Þ

(40)

where f(t) and gi(t), i1 are analytical functions depending on Rj ðtÞ, ji1 and

FðtÞ ¼

f ðtÞ

2 R0 1

 R0

; Gi ðtÞ ¼

gi ðtÞ ; 2 R0 1  R0

Ri ðt ¼ 0Þ ¼ 0:

Therefore, the 0th-order growth of the interface is a solution of the following equation 3 2   R R  0 þ ð1  aÞ 0 þ að1  aÞR0 þ a2 ð1  aÞlnR0  a ¼ dt þ C; 3 2 (34)

with C being a constant that is determined by the initial condition. This has also been presented in [3] for d ¼ 1 with the assumption R0 > a. The functions Ri ðtÞ, i  1 can also be analytically expressed, i.e.

Remark 1. (i) Eqs. (11) and (12) can be written in the general form




vw vw v vw ¼ xm xm s2 x; t; w; s1 x; t; w; vx vt vx vx 
vw þ s3 x; t; w; vx

(41)

392

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395

Fig. 5. Absolute error of the solution with respect to RNS.

for 0  t  tstable and 0  x  1, x ¼ x,h, a corresponding initial condition w(x,t0) ¼ w0(x) and associated boundary conditions in Eq. (13) for w ¼ v,u, respectively. This requires the coefficient functions s1, s2, s3 to be given in closed form. In our case R determines the functions s1, s2, s3, but it is only given in a discrete way. A spline interpolation procedure is used to obtain such closed form.

k , is also evaluated between. The absolute error at iteration k, Eabs in these cases (Fig. 5) using

(ii) A numerical routine [28] based on a second order (in space and time) centered implicit finite difference scheme is used to solve Eqs. (11) and (12). (iii) The accuracy of the scheme is improved by considering fivepoint backward/forward differences for the evaluation of the partial derivatives at the interface for the modified temperatures v and u, respectively, in Eq. (16). This will offer a more accurate approximation of the spatial derivatives compared to a three-point difference.

where RRNS is given as the reference numerical solution (RNS) obtained using the numerical procedure described in [3]. Numerical values of the physical parameters for the considered ambient temperatures in the numerical calculations are listed in Table 2. All analytical results are compared with the RNS values. The comparison is also extended to the iterative procedure results showing a qualitatively good agreement. The error occurring in the numerical evaluations in the PTS can be associated on the one hand with the difference St  St  (Table 3). A numerical error is also induced by considering the vanishing series in Eqs. (29) and (30), however, the numerical data improve after we apply the iterative procedure.

We test the evolution of the solution after the iterative procedure for realistic high and low supercooling cases. However, the procedure is valid for all values of the ambient temperature in-

    k Eabs ðtÞ ¼ Rk ðtÞ  RRNS ðtÞ;

Fig. 6. Modified temperature distribution in ice and water at tn ¼ 0.693,105: RNS and PT-IP at T∞ ¼ 243.15 K.

(42)

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395

Fig. 7. Modified temperature distribution in ice and water at tn ¼ 0:360$103 : RNS and PT-IP at T∞ ¼ 263.15 K.

Fig. 8. Position of the moving interface: comparisons between the RNS, PT-IP and IG-IP for t2(0,tstable(RNS)).

393

394

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395

Table 4 Comparison of the CPU time and the time tstable to reach Rstable for all schemes. T∞/[K]

CPU time/[s]

243.15

PT-IP 78.92 PT-IP 35.58 PT-IP 125.80 PT-IP 246.21

248.15 258.15 263.15

IG-IP 132.28 IG-IP 68.74 IG-IP 194.17 IG-IP 405.02

tstable/[e] FD5 433.2 FD5 75.6 FD5 142.54 FD5 835.8

PT-IP 0.660e-5 PT-IP 0.1456e-4 PT-IP 0.9366e-4 PT-IP 0.368e-3

IG-IP 0.660e-5 IG-IP 0.1456e-4 IG-IP 0.9366e-4 IG-IP 0.368e-3

FD5 0.8e-5 FD5 0.17e-4 FD5 0.878e-4 FD5 0.345e-3

As part of the solution, the modified temperature distributions in the water and the ice phase are also of importance, in particular at the time when the interface position reaches Rstable. These fields represent the initial condition in a direct numerical simulation of the overall dendritic growth of the ice particle. A quasi-steady temperature distribution is observed in the ice phase (Figs. 6a and 7a). In the water phase, for very small values of the spatial variable, i.e. r close to Rstable, a different behaviour is observed (Figs. 6b and 7b). The iterative solution is based on an initial guess of the iceewater interface and the PTS provides in our case a good approximation of this function. Let us consider a different initial guess for the position of the ice front, RðtÞ ¼ Rc =Rd þ t 1=2 , which is chosen taking into account the iceewater interface behaviour and the initial condition R(t ¼ 0) ¼ Rc/Rd. The iterative procedure in this case, denoted IG-IP, also converges (Fig. 8), yet more CPU time is needed to achieve the target error level eMaxetol (Table 4). The numerical solution of Eqs. (11)e(16) may also be achieved by using a finite difference scheme. Following [8] we employ a semi-implicit scheme (FD5) for the system: we discretize implicitly for v and u and explicitly for R [8]

nþ1  vnþ1  vnl vv vnþ1 vv vl l1 z z lþ1 ; ; vt vx Dt 2Dx

v2 v vx2

z

vnþ1  2vnþ1 þ vnþ1 i lþ1 l1 Dx2

(43)

;

5. Conclusions A perturbation technique followed by an iterative procedure was considered in this paper in order to provide an approximate solution of the radially symmetric ice growth in a supercooled droplet. The result is a numerical scheme that can be easily implemented for evaluating the interface position and the modified temperature distributions in the ice and water phase. The analytical solution, expressed with two correction terms in the perturbation expansion, offers a good agreement with the existing reference numerical results. The iterative solution rapidly converges towards the numerical solutions such that the computational effort is significantly reduced. There are no restrictions imposed on the range of validity in terms of ambient temperatures/ Stefan numbers, the error is of the same order for all tested cases. The transformation in the expansion parameter in order to get an analytical expression of the phase growth function is not unique. Here, with this transformation, the parameter d can be any expression/value that allows a good approximation of the Stefan number using this transformation, but, the initial error in the analytical expression of the interface position will be influenced by this choice. However, the iterative algorithm will converge towards the benchmark numerical results for all d and all corresponding analytical evaluations of the starting approximation phase growth function satisfying the initial condition. Also the predicted modified temperature distributions in the ice and water phase are in good agreement with the numerical distribution of the corresponding temperatures. The comparison with the modified temperature distribution obtained by the numerical calculations in [3] confirms the quasi-steady behaviour of the temperature in the ice phase. Acknowledgements The authors gratefully acknowledge financial support from DFG within the Collaborative Research Center SFB-TRR 75.

where l ¼ 1,..,M and n ¼ 1..,N. An analogous discretization is used for u. Then Eqs. (11) and (12) can be solved at each time step n. The phase front growth R is obtained from Eq. (16) using the explicit discretization

dR Rnþ1  Rn z : dt Dt

The PT-IP method offers a quantitatively sufficient approximation of the interface as well as the modified temperature distributions in the ice and in the water phase in comparison with the considered numerical reference results and it is less expensive in terms of time and memory resources.

(44)

A five-point backward/forward difference in the evaluation of the partial derivatives at the interface of the modified temperatures v and u, respectively in Eq. (16) is also used in this case. A comparison of all computational times to reach Rstable and also the value of tstable for each case are compared in Table 4 for a time step Dt ¼ 109. We notice a difference for the FD5 scheme in the values of tstable, the time needed to reach Rstable in the numerical procedure. Mesh refinements are possible in order to correct this value, however, computational time costs become much larger. The iterative method provides a solution very close to RNS in a reasonable time, making the refinement process superfluous in our case.

References [1] J.P. Hindmarsh, A.B. Russell, X.D. Chen, Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet, Int. J. Heat. Mass Transf. 46 (2003) 1199e1213, http://dx.doi.org/10.1016/S00179310(02)00399-X. [2] W.W. Mullins, R.F. Sekerka, Morphological stability of a particle growing by diffusion or heat flow, J. Appl. Phys. 34 (2) (1963) 323329, http://dx.doi.org/ 10.1063/1.1702607. [3] Eisenschmidt, K., Rauschenberger, P., Rohde, C., Weigand, B., Modelling of freezing processes in supercooled droplets on sub-grid scale, ILASS-Europe 2013, 25th European Conference on Liquid Atomization and Spray Systems, Chania, Greece, 1e4 September (2013). [4] J. Caldwell, K. Yuen-Yick, A brief review of several numerical methods for onedimensional Stefan problems, Therm. Sci. 13 (2) (2009) 61e72, http:// dx.doi.org/10.2298/TSCI0902061C. [5] M.K. Bernauer, R. Herzog, Implementation of an X-FEM solver for the classical two-phase Stefan problem, J. Sci. Comput. 52 (2) (2011) 271e293, http:// dx.doi.org/10.1007/s10915-011-9543-x. [6] S.W. McCue, B. Wu, J.M. Hill, Classical two-phase stefan problem for spheres, Proc. R. Soc. A: Math. Phys. Eng. Sci. 464 (2012) 2055e2076, http://dx.doi.org/ 10.1098/rspa.2007.0315. [7] S.L. Mitchell, S.B.G. O'Brien, Asymptotic, numerical and approximate techniques for a free boundary problem arising in the diffusion of glassy polymers, Appl. Math. Comput. 219 (2012) 376e388, http://dx.doi.org/10.1016/ j.amc.2012.06.026.

F.I. Dragomirescu et al. / International Journal of Thermal Sciences 104 (2016) 386e395 [8] F. Font, T.G. Myers, S.L. Mitchell, A mathematical model for nanoparticle melting with density change, Microfluid Nanofluid 18 (2015) 233e243, http:// dx.doi.org/10.1007/s10404-014-1423-x. [9] J.M. Back, S.W. McCue, T.J. Moroney, Including nonequilibrium interface kinetics in a continuum model for melting nanoscaled particles, Sci. Rep. 4 (2014) 7066, http://dx.doi.org/10.1038/srep07066. [10] H. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Clarendon Press, Oxford, 1959. [11] S.W. McCue, B. Wu, J.M. Hill, Micro/nanoparticle melting with spherical symmetry and surface tension, IMA J. Appl. Math. 74 (2009) 439e457, http:// dx.doi.org/10.1093/imamat/hxn038. [12] J. Chadam, S.D. Howison, P. Ortoleva, Existence and stability for spherical crystals growing in a supersaturated solution, IMA J. Appl. Math. 39 (1987) 1e15, http://dx.doi.org/10.1093/imamat/39.1.1. [13] M.A. Herrero, J.J.L. Velasquez, On the melting of ice balls, SIAM J. Math. Anal. 28 (1997) 1e32, http://dx.doi.org/10.1137/S0036141095282152. [14] F. Font, T.G. Myers, Spherically symmetric nanoparticle melting with a variable phase change temperature, J. Nanopart. Res. 15 (2013) 2086, http:// dx.doi.org/10.1007/s11051-013-2086-3. [15] B. Wu, S.W. McCue, P. Tillman, J.M. Hill, Single phase limit for melting nanoparticles, Appl. Math. Model. 33 (2009) 2349e2367, http://dx.doi.org/ 10.1016/j.apm.2008.07.009. [16] S. Weinbaum, L.M. Jiji, Singular perturbation theory for melting or freezing in finite domains initially not at the fusion temperature, J. Appl. Mech. 44 (1) (1977) 25e30, http://dx.doi.org/10.1115/1.3424008. [17] S. Savovic, J. Caldwell, Finite difference solution of one-dimensional Stefan problem with periodic boundary conditions, Int. J. Heat. Mass Transf. 46 (2003) 2911e2916, http://dx.doi.org/10.1016/S0017-9310(03)00050-4. [18] L.M. Jiji, S. Weinbaum, Perturbation solutions for melting or freezing in annular regions initially not at the fusion temperature, Int. J. Heat Mass Transf. 21 (1978) 581e592, http://dx.doi.org/10.1016/0017-9310(78)90055-8. [19] V. Alexiades, A.D. Solomon, Mathematical Modelling of Melting and Freezing Processes, Hemisphere, Washington, 1993. [20] J.R. King, J.R. Evans, Regularization by kinetic undercooling of blow-up in the ill-posed Stefan problem, SIAM J. Appl. Math. 65 (2005) 1677e1707, http:// dx.doi.org/10.1137/04060528X. [21] G.B. Davis, J.M. Hill, A moving boundary problem for the sphere, IMA J. Appl. Math. 29 (1) (1982) 99e111, http://dx.doi.org/10.1093/imamat/29.1.99. [22] J.M. Back, S.W. McCue, M.H.N. Hsieh, T.G. Moroney, The effect of surface tension and kinetic undercooling on a radially-symmetric melting problem, Appl. Math. Comput. 229 (2014) 41e52, http://dx.doi.org/10.1016/ j.amc.2013.12.003. [23] G. Brenn, Concentration fields in evaporating droplets, Int. J. Heat. Mass Transf. 48 (2) (2005) 395e402, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2004.07.039. [24] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover Publications, New York, 1972. [25] J. Kevorkian, J.D. Cole, Perturbation methods in applied mathematics, Appl. Math. Sci. 34 (1981). Springer-Verlag. [26] J. Struckmeier, A. Unterreiter, A singular-perturbed two-phase Stefan problem due to slow diffusion, Appl. Math. Lett. 14 (2) (2001) 217e222, http:// dx.doi.org/10.1016/S0893-9659(00)00139-7. [27] J.R. King, D.S. Riley, Asymptotic solutions to the Stefan problem with a constant heat source at the moving boundary, Proc. R. Soc. A Math. Phys. Eng. Sci. 456 (2000) 1163e1174, http://dx.doi.org/10.1098/rspa.2000.0556. [28] Maple 18 (March, 2014), Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario.

Nomenclature

Latin aw: thermal diffusivity of the water, [m2/s] ai: thermal diffusivity of the ice, [m2/s] cw: specific heat capacity, [J/(kgK)] eMax: max. diff. between Rkþ1 and Rk, [e] etol: level of tolerance in the iterative procedure, [e] i,w: ice, water, indices kw: heat conductivity of the water, [W/(mK)] ki: heat conductivity of the ice, [W/(mK)] L: specific latent heat, [J/kg] Rdim : dimensional maximum spherical ice particle radius, [m] stable Rstable: maximum spherical ice particle radius, [e] r 0 : radial coordinate, [m] R0 : interface position, [m] R: nondim. interface position, [e] R: nondim. interface position in the modified approach, [e] r: radial coordinate, [e] Rc: initial nucleus radius, [m] Rd: water droplet radius, [m] Rk: interface position at iteration k, [e] St: Stefan number,St ¼ cw(TmT∞)/L, [e] St: Stefan number in the modified approach, [e] dim : duration for reaching Rdim , [s] tstable stable tstable: duration for reaching Rstable, [e] t 0 : time, [s] t: non-dimensional time, [e] T: temperature, [K] Tm: melting temperature, [K] T∞: ambient temperature, [K] TGT: GibbseThomson temperature, [K] tol: tolerance in approximating St with St, [e] u,v: modified temperatures in the water and in the ice phase, [e] Greek

a: surface tension coefficient; a ¼ 2sTm =ðTm  T∞ ÞrLRD , [e] d: artificial parameter, [e]

ε: expansion parameter in the standard approach, [e] ε: expansion parameter in the modified approach, [e] h,x: coordinates in the water and in the ice phase, [e] r: density, [kg/m3] s: surface tension, [N/m] Q: non-dimensional temperature, [e] Acronyms

RNS: Reference Numerical Solution PTS: Perturbation Technique Solution PT-IP: Perturbation Technique followed by the Iterative Procedure IG-IP: arbitrary Initial Guess followed by the Iterative Procedure

395