- Email: [email protected]
Transient heating of an evaporating droplet with presumed time evolution of its radius
International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Transient heating of an evaporating droplet with presumed time evolution of its radius S.S. Sazhin a,⇑, P.A. Krutitskii b, I.G. Gusev a, M.R. Heikal a a b
Sir Harry Ricardo Laboratories, Centre for Automotive Engineering, School of Computing, Engineering and Mathematics, University of Brighton, Brighton BN2 4GJ, UK Keldysh Institute for Applied Mathematics, Department 4, Miusskaya Sq. 4, Moscow 125047, Russia
a r t i c l e
i n f o
Article history: Received 8 April 2010 Accepted 14 October 2010 Available online 1 December 2010 Keywords: Droplets Diesel fuel Heating Evaporation Moving boundary Analytical solution Stefan problem
a b s t r a c t New solutions to the heat conduction equation, describing transient heating of an evaporating droplet, are suggested, assuming that the time evolution of droplet radius Rd(t) is known. The initial droplet temperature is assumed to be constant or allowed to change with the distance from the droplet centre. Since Rd(t) depends on the time evolution of the droplet temperature, an iterative process is required. Firstly, the time evolution of Rd(t) is obtained using the conventional approach, when it remains constant during the timestep, but changes from one timestep to another. Then these values of Rd(t) are used in the new solutions to obtain updated values of the time evolution of the distribution of temperatures inside the droplet and on its surface. These new values of droplet temperature are used to update the function Rd(t). This process continues until convergence is achieved, which typically takes place after about 15 iterations. The results of these calculations are compared with the results obtained using the previously suggested approach when the droplet radius was assumed to be a linear function of time during individual timesteps for typical Diesel engine-like conditions. For sufficiently small timesteps the time evolutions of droplet temperatures and radii predicted by both approaches coincided. This suggests that both approaches are correct and valid. Similarly to the case when droplet radius is assumed to be a linear function of time during the timestep, the new solutions predict lower droplet temperature and slower evaporation when the effects of the reduction of Rd are taken into account. It is shown that in the case of constant droplet initial temperature, models both taking and not taking into account the changes in the initial droplet temperature with the distance from the droplet centre predict the same results. This indicates that both models are correct. Crown Copyright Ó 2010 Published by Elsevier Ltd. All rights reserved.
1. Introduction The conventional approach to the modelling of droplet heating and evaporation is based on the assumption that droplet radius remains constant during each timestep, but changes from one timestep to another due to the evaporation process [1,2]. In our previous paper [3], we took into account the time evolution of droplet radius Rd(t) assuming that Rd(t) is a linear function of time during the timestep. It was shown that this approach not only allows the reduction of the number of timesteps required for the calculation, but also leads to the prediction of lower droplet temperatures and longer evaporation times compared with the conventional approach. These new and in practice important results encouraged us to investigate this problem further and relax our assumption that Rd(t) is a linear function of time. The results of
⇑ Corresponding author. Tel.: +44 0 1273642677; fax: +44 0 1273642309. E-mail address: [email protected] (S.S. Sazhin).
our investigation are presented in this paper. Some preliminary results and their analysis were presented in [4]. Basic equations and approximations are presented and discussed in the next section. In Section 3, the solution based on the assumption that Rd(t) is a linear function of time, obtained in [3], is reproduced. The new solution, based on the assumption that Rd(t) is an arbitrary (but twice continually differentiable) function of time, but the initial droplet temperature Td0 does not depend on the distance from the droplet centre R, is derived in Section 4. In Section 5, the most general solution for arbitrary Rd(t) and Td0(R) is presented. The application of the new solutions to calculation of Rd(t) and T(R, t) of heated and evaporated droplets is described in Section 6. The results of calculation of droplet temperatures and radii for Diesel engine-like conditions using the conventional approach (Rd = const during the timestep), the previously suggested approach (Rd is the linear function of t during the timestep) and the new approach are presented, compared and discussed in Section 7. The main results of the paper are summarised in Section 8.
0017-9310/$ - see front matter Crown Copyright Ó 2010 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2010.10.018
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
1279
Nomenclature BM c fn G h h0 H(t) k L M(t) qn R Rd t T u U(t, R)
v vn W Y
Greek symbols evaporation rate defined by Eq. (9) Hn function defined by Eq. (11) j thermal diffusivity k, kn eigenvalues defined by Eq. (15) l(t) function defined by Eq. (5) l0(t) function defined by Eq. (12) m(t) function defined by Eq. (22) q density x(t, s) function defined by Eq. (A.13) X(t, s) function defined by Eq. (A.12)
Spalding mass transfer number specific heat capacity function introduced in Eq. (11) function introduced in Eq. (22) convection heat transfer coefficient R0 ðtÞR ðtÞ hðtÞ Rd ðtÞ 1 d 2jd kl function introduced in Eq. (5) thermal conductivity specific heat of evaporation function introduced in Eq. (3) function introduced in Eq. (16) distance from the droplet centre droplet radius time temperature TR function defined by Eq. (30) u R Td0 or function defined by Eq. (34) function defined by Eq. (13) function defined by Eq. (17) mass fraction
a
Subscripts d droplet e evaporation g gas l liquid n timesteps in the integration p constant pressure s surface v vapour 0 initial
for t 2 (0, te), R 2 (0, Rd(t)) with the boundary conditions
2. Basic equations and approximations
Let us assume that a stationary evaporating droplet is immersed into a homogeneous hot gas at temperature Tg. The droplet is heated by convection with the convection heat transfer coefficient h(t) depending on time t and droplet radius Rd(t), and cooled down due to evaporation. Rd(t) is a twice continuously differentiable function of time in the range 0 6 t 6 te, where te is the evaporation time. Both Rd(t) and h(t) are assumed to be known. Effects of thermal radiation are ignored, which is justified for typical Diesel fuel droplets and gas temperatures about or less than 1000 K, when the effect of radiation on droplet evaporation time is typically less than 1% [5]. The changes in the droplet temperatures are described by the heat conduction equation in the form [6,7]:
@T @ 2 T 2 @T ¼j þ @t @R2 R @R
! ð1Þ
for 0 < t < te, 0 < R < Rd(t), where j = kl/(cl ql) = const, kl is the liquid thermal conductivity, cl is the liquid specific heat capacity, ql is the liquid density, R is the distance from the centre of the droplet. Eq. (1) needs to be solved subject to the boundary condition:
@T kl ¼ hT g þ ql LR_ d ðtÞ; þ hT @R R¼Rd ðtÞ
ð2Þ
T is finite and continuous at R ? +0, Ts T(Rd(t), t) is the droplet’s surface temperature, L is the latent heat of evaporation. We took into account that R_ d ðtÞ dRd =dt 6 0. Eq. (2) is the energy balance condition at R = Rd(t). The initial condition is taken in the form T(t = 0) = Td0(R), where 0 6 R 6 Rd0 = Rd(t = 0). Let us rewrite boundary condition (2) in the form
@T h h q ¼ T g þ l LR_ d ðtÞ MðtÞ; þ T @R kl kl kl R¼Rd ðtÞ
ð3Þ
and introduce the new variable u = TR. Using this new variable, we can rewrite Eq. (1) as
@u @2u ¼j 2 @t @R
ð4Þ
@u þ HðtÞu ¼ lðtÞ; @R R¼Rd ðtÞ
ujR¼0 ¼ 0;
ð5Þ ð6Þ
where
HðtÞ ¼
hðtÞ 1 ; kl Rd ðtÞ
lðtÞ ¼ MðtÞRd ðtÞ
for t 2 (0, te). The initial condition is
uðRÞjt¼0 ¼ RT d0 ðRÞ
ð7Þ
for R 2 (0, Rd0). 3. Solution for the case when Rd(t) is a linear function Let us assume that Rd(t) is the linear function of t during the time step
Rd ðtÞ ¼ Rd0 ð1 þ atÞ;
ð8Þ
where Rd0 is the initial droplet radius. Remembering the physical background to the problem, the expression for a can be presented as [2]:
a¼
kg lnð1 þ BM Þ
ql cpg R2d0
;
ð9Þ
where BM = Yvs/(1 Yvs) is the Spalding mass transfer number, Yvs is the mass fraction of fuel vapour near the droplet surface, kg and cpg are gas thermal conductivity and specific heat capacity respectively. For Rd(t) defined by (8) the solution to Eq. (1), subject to the above boundary and initial conditions, was obtained in the form [3]:
" # 1 aRd0 R2 TðRÞ ¼ pffiffiffiffiffiffiffiffiffiffiffi exp 4jRd ðtÞ R Rd ðtÞ " # 1 X R l ðtÞ R ; Hn ðtÞ sin kn þ 0 Rd ðtÞ 1 þ h0 Rd ðtÞ n¼1
ð10Þ
1280
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
where
"
Hn ðtÞ ¼ Hn ð0Þ exp " exp
k2n R2d0
j a
k2n R2d0
j a
1 1 1 þ at
#
1 1 1 þ at 1 þ as
pffiffiffiffiffiffiffiffiffiffiffi
l0 ðtÞ l~ ðtÞ Rd ðtÞ exp
where
þ fn
# ds;
Z 0
t
Gðt; s; RÞ ¼ G0 ðt s; R Rd ðsÞÞ G0 ðt s; R þ Rd ðsÞÞ; pffiffiffiffi j x2 : G0 ðt; xÞ ¼ pffiffiffiffiffiffi exp 4jt 2 pt
dl0 ðsÞ ds ð11Þ
0 Rd ðtÞRd ðtÞ ; 4j
ð12Þ
l~ ðtÞ ¼ MðtÞR2d ðtÞ; MðtÞ ¼ khl T g þ qkll LR_ d ðtÞ; f n ¼ kvsinn kk2nk2 , n
v n ðnÞ ¼ sin kn n
ðn ¼ 1; 2; . . .Þ ! 1 sin 2kn 1 h0 2 ; ¼ 1 1þ 2 kv n k ¼ 2 2 2kn h0 þ k2n
ð13Þ ð14Þ
Expression G(t, s, R) can be presented in an alternative form: ( " # " #) pffiffiffiffi j ðR Rd ðsÞÞ2 ðR þ Rd ðsÞÞ2 exp : Gðt; s; RÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 4jðt sÞ 4jðt sÞ 2 pðt sÞ ð23Þ
Note that G(t, s, R = 0) = 0. m(t) is an unknown continuous function to be found later from one of the boundary conditions. Function v(R, t) is known as a single layer heat potential and it has the following properties for any continuous function m(t) [7,8]: (1) (2) (3) (4)
kn are positive solutions to the equation
k cos k þ h0 sin k ¼ 0;
ð15Þ
presented in ascending order,
h0 ¼
hðtÞ R0 ðtÞRd ðtÞ Rd ðtÞ 1 d kl 2j
Z t @ v ðR; tÞ mðtÞ @Gðt; s; RÞ þ ¼ m ð s Þ ds: @R R!Rd ðtÞ0 2 @R 0 R¼Rd ðtÞ
is assumed to be constant during the time step,
Hn ð0Þ ¼ qn þ l0 ð0Þfn ; qn ¼
Z
1 2
kv n k
W 0 ðnÞv n ðnÞdn:
3=2
W 0 ðnÞ ¼ Rd0 nT 0 ðnRd0 Þ exp
R0d ð0ÞRd0 2 n : 4j
ð17Þ
4. Solution for the case of arbitrary Rd(t) but Td0(R) = const The analysis of this section is based on the assumption that Td0(R) = Td0 = const. In this case we can introduce the new variable v = u R Td0 and rearrange Eq. (4) as:
@v @2v ¼j 2 @t @R
ð18Þ
@Gðt; s; RÞ 1 ¼ pffiffiffiffiffiffiffi @R 4 pjðt sÞ3=2 R¼Rd ðtÞ (
# ðRd ðtÞ Rd ðsÞÞ2 ðRd ðtÞ Rd ðsÞÞ exp 4jðt sÞ " #) ðRd ðtÞ þ Rd ðsÞÞ2 : ðRd ðtÞ þ Rd ðsÞÞ exp 4jðt sÞ
@v ¼ l0 ðtÞ; þ HðtÞv @R R¼Rd ðtÞ
ð19Þ
v jR¼0 ¼ 0
ð20Þ
for t 2 (0, te) and the initial condition
Since Rd(t) is a continuously differentiable function, we obtain in the limit s ? t 0:
@Gðt; s; RÞ 1 p ffiffiffiffiffiffiffiffiffiffi ffi : / O @R ts R¼Rd ðtÞ
mðtÞ
þ
Z
t
mðsÞ
0
Z t @Gðt; s; RÞ d s þ HðtÞ mðsÞGðt; s; Rd ðtÞÞds @R 0 R¼Rd ðtÞ
¼ l0 ðtÞ; or
mðtÞ 2
hðtÞ l0 ðtÞ ¼ T d0 HðtÞRd ðtÞT d0 þ lðtÞ ¼ Rd ðtÞT d0 þ lðtÞ: kl We look for the solution to the problem (18)–(21) in the form:
Z
0
t
mðsÞGðt; s; RÞds;
ð22Þ
ð26Þ
It follows from this equation, that there is an improper integral on the right hand side of Eq. (24). In view of Eqs. (24) and (22) we can rewrite the boundary condition (19) as:
þ
Z 0
t
mðsÞ
( ) @Gðt; s; RÞ þ HðtÞGðt; s ; R ðtÞÞ ds ¼ l0 ðtÞ; d @R R¼Rd ðtÞ
ð21Þ
for R 2 (0, Rd(t)), where
"
ð25Þ
2
for t 2 (0, te) and R 2 (0, Rd(t)) with the boundary conditions
v jt¼0 ¼ 0
This means that for any continuous function m(t) the potential satisfies Eq. (18) and boundary and initial conditions (20) and (21). Choice of the function m(t), satisfying integral Eq. (24), should be made in such a way that the remaining boundary condition (19) is satisfied as well. From Eq. (23) it follows that
v(R, t)
The first step in the derivation of this equation is based on the reduction of the problem (1) and (2) to the problem (4)–(7) with the corresponding initial conditions. Note that the condition h0 = const can be violated for t in the immediate vicinity of te for finite timesteps. This imposes a rather severe restriction on the choice of the timestep just before the droplet fully evaporates. Condition (8) allows us to use the above solution for individual timesteps but not for the whole period of droplet heating and evaporation. In what follows this assumption will be relaxed and we consider arbitrary Rd(t) provided that this function is twice continuously differentiable (natural condition remembering the physical background to the problem).
v ðR; tÞ ¼
ð24Þ
ð16Þ
1
0
It satisfies Eq. (18) for 0 < t < te and 0 < R < Rd(t); It satisfies Conditions (20) and (21); It is continuous at R ? Rd 0; For the derivative @ v (R, t)/@R the following limiting formula is valid:
ð27Þ where G(t, s, R) and its derivative with respect to R are defined by Eqs. (23) and (25). As follows from the definition of G(t, s, Rd(t)) (see Eq. (23)), in the limit s ? t 0 this function has the singularity:
1 Gðt; s; Rd ðtÞÞ / O pffiffiffiffiffiffiffiffiffiffiffi ts
1281
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
(cf. Eq. (26)). Therefore, the integral in Eq. (27) is defined as an improper integral. Eq. (27) is an integral equation of Volterra type. It has a unique continuous solution. A scheme for its numerical solution is described in Appendix A. This solution is then substituted into Eq. (22). The final distribution of temperature inside the droplet can be calculated from the following expression: pffiffiffiffi Z t j pffiffiffiffi 2R p 0 ( " # " #) 2 2 m ð sÞ ðR Rd ðsÞÞ ðR þ Rd ðsÞÞ pffiffiffiffiffiffiffiffiffiffiffi exp exp ds: 4jðt sÞ 4jðt sÞ ts
Tðt; RÞ ¼ T d0 þ
ð28Þ
Details of the numerical calculation of the integral on the right hand side of Eq. (28) are given in Appendix B.
Having substituted Eq. (34) into Eq. (4) and boundary and initial conditions (5)–(7), we obtain problem (18)–(21) for v(t, R) in which
l0 ðtÞ ¼ U 0R ðt; RÞ þ HðtÞUðt; RÞ R¼Rd ðtÞ þ lðtÞ:
ð35Þ
The solution of the latter problem is similar to the one discussed in Section 4. The expression for l0(t) contains
U 0R ðt; RÞR¼R
d ðtÞ
¼
Z
Reff 0
@G1 ðt; R; fÞ ðfT do ðfÞÞ df; @R R¼Rd ðtÞ
ð36Þ
where
@G1 ðt; R; fÞ 1 @ G ¼ ð ðt; R fÞ G ðt; R þ fÞ Þ 0 0 @R j @R R¼Rd ðtÞ R¼Rd ðtÞ 1 @ ¼ ðG0 ðt; R fÞ G0 ðt; R þ fÞÞ ; j @f R¼Rd ðtÞ ð37Þ
5. Solution for the case of arbitrary Rd(t) and Td0(R) Let us assume that an arbitrary twice continuously differentiable function Td0(R) is defined for 0 6 R 6 Rd0. This definition is extended for the whole range 0 6 R < 1:
8 > < T d0 ðRÞ when 0 6 R 6 Rd0 ; T d0 ðRÞ ¼ T out ðRÞ when Rd0 < R 6 Reff ; > : 0 when R > Reff ;
ð29Þ
G0 is the same as in Eq. (22). Remembering the latter equation we ðt;R;fÞ can rewrite the expression for @G1@R in a more explicit form:
( R¼Rd ðtÞ " # @G1 ðt; R; fÞ 1 ðR fÞ2 ¼ ðR fÞ exp p ffiffiffi ffi @R 4jt 4 pðjtÞ3=2 R¼Rd ðtÞ " #) ðR þ fÞ2 : ðR þ fÞ exp 4jt R¼Rd ðtÞ
where
Hence, we obtain an explicit expression for l0(t) in the form:
h io 1n Rd0 T d0 ðRd0 Þ þ ðR Rd0 Þ ðRT d0 ðRÞÞ0R R¼R T out ðRÞ ¼ ; d0 R
l0 ðtÞ ¼ pffiffiffiffi 3=2 4 pðjtÞ
1
Reff is the effective outer radius such that Reff > Rd0. Function Td0(R) defined by Eq. (29) is continuously differentiable in the range 0 6 R 6 Reff. Let us now introduce a new function U(t, R) defined as:
Uðt; RÞ ¼
Z
Reff
ðfT do ðfÞÞG1 ðt; R; fÞdf;
ð30Þ
0
where
G1 ðt; R; fÞ ¼
1
j
ð31Þ
(1) It satisfies Eq. (18) for 0 < t < te and 0 < R < 1; (2) It satisfies the boundary Condition (20) for 0 < t < te; (3) It satisfies the initial condition
RT d0 ðRÞ when 0 6 R 6 Reff ; 0
when R > Reff :
ð32Þ
The latter relation follows from the property of the delta-function:
lim
adelta !1
adelta pffiffiffiffi exp a2delta x2 ¼ dðxÞ: p
ð33Þ
We look for the solution to Eq. (4) in the form:
uðt; RÞ ¼ Uðt; RÞ þ v ðt; RÞ:
0
ðRd ðtÞ fÞ2 4jt
#
"
ð34Þ
i
l0 ð0Þ ¼ ðfT do ðfÞÞ0f f¼Rd0 þ Hð0ÞRd0 T d0 ðRd0 Þ þ lð0Þ:
ð38Þ
ð39Þ
Combining Eqs. (22) and (34) we can present the final solution to our problem in the form:
Tðt; RÞ ¼
Note that G1(t, R = 0) = 0. Function U(t, R) has the following properties [7,8]:
"
ðfT d0 ðfÞÞ ðRd ðtÞ fÞ exp
#) ðRd ðtÞ þ fÞ2 df ðRd ðtÞ þ fÞ exp 4jt ( " # Z Reff HðtÞ ðRd ðtÞ fÞ2 pffiffiffiffiffiffiffiffiffi ðfT d0 ðfÞÞ exp 4jt 2 pjt 0 " #) ðRd ðtÞ þ fÞ2 df þ MðtÞRd ðtÞ: exp 4jt
h
G0(t, x) is the same as in Eq. (22). Remembering the latter equation, the expression G1(t, R) can be presented in an alternative form:
Uðt; RÞjt¼þ0 ¼
(
Reff
In the limit t ? 0+ the expression for l0(t) is simplified to (see Appendix C):
½G0 ðt; R fÞ G0 ðt; R þ fÞ;
( " # " #) 1 ðR fÞ2 ðR þ fÞ2 exp : G1 ðt; RÞ ¼ pffiffiffiffiffiffiffiffiffi exp 4jt 4jt 2 pt j
Z
" ( " # pffiffiffiffi Z t 1 j mðsÞ ðR Rd ðsÞÞ2 pffiffiffiffiffiffiffiffiffiffiffi exp Uðt; RÞ þ pffiffiffiffi R 4jðt sÞ 2 p 0 ts " #) # ðR þ Rd ðsÞÞ2 ds ; exp 4jðt sÞ
ð40Þ
where m(s) is the solution to Eq. (27) with l0(t) given by Eq. (38), U(t, R) is given by Eq. (30). Note that taking into account the initial distribution of temperature along R is absolutely essential when the solution is applied to individual timesteps. In the solution described in the last two sections, however, the same formulae describe the time evolution of droplet temperatures during the whole period of their evaporation. It is anticipated that in this case the effect of the initial distribution of droplet temperatures is not important in most practical applications. Hence, the solution described in Section 4 is expected to be more practically important than the solution described in this section.
1282
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Note that although the analysis presented so far refers to stationary droplets, it can be generalised in a straightforward way to the case of moving droplets, based on the effective thermal conductivity model [1,2]. 6. Implementation of the new solutions into a numerical code In the solutions presented in the last two sections it was assumed that Rd(t) is known. From the point of view of the physical background to the problem, however, Rd(t) depends on the time evolution of the droplet temperature T(R, t), which is the solution to the problem. Hence, an iterative process is required. Firstly, the time evolution of droplet radius Rd(t) is obtained using the conventional approach, when it remains constant during the timestep, but changes from one timestep to another due to the evaporation process. Then these values of Rd(t) are used in the new solutions to obtain updated values of the time evolution of the distribution of temperatures inside the droplet and on its surface T(R, t). These new values of droplet temperature are used to update the function Rd(t). This process is continued until convergence is achieved, which typically takes place after about 15 iterations. 7. Results In Figs. 1–7 we compared the results of calculations of droplet surface temperatures and radii, taking into account the effects of evaporation, using the new integral solution for arbitrary Rd(t) but constant Td0 (Eq. (28)), the previously reported solution, based on the linear approximation of Rd(t) (Eq. (10)), and the conventional approach based on the assumption that droplet radius does not change during the timestep. Droplets are assumed to be those of n-dodecane (Mf = 170 kg/kmole), and ambient air is assumed to
be at the pressure of p = 30 atm = 3000 kPa (typical values for Diesel engine-like conditions). The results of calculations for Rd0 = 5 lm, and ambient air temperature 1000 K are shown in Fig. 1. One thousand timesteps were used in the conventional approach and the approach based on Solution (10). In the integral solution based on Eq. (28) the integral (A21) was approximated as the sum of 100 terms, and up to 15 iterations were used. At the first iteration the time evolution of the droplet radius was assumed to be the same as predicted by the conventional model. As follows from Fig. 1, the results predicted by the integral solution (28) and linear solution (10) practically coincide, which suggests that both approaches are correct and valid. Both these solutions predict lower droplet temperatures and longer evaporation times in agreement with the results reported in [3]. Note that deviations between the predictions of the integral and linear solutions were observed in the immediate vicinity of the time when the droplet completely evaporates. There were obvious numerical problems when we approached this time due to the fact that the time derivative of Rd becomes infinitely large. In practice the extrapolation, based on the assumption that the second derivative of Rd(t) is constant, was used for these times. This leads to small deviations between the predicted evaporation times. In the case shown in Fig. 1, the evaporation times predicted by the conventional model, linear solution, and integral solution were 0.595 ms, 0.622 ms and 0.628 ms respectively. That means that the difference between the evaporation times predicted by the linear and integral solutions was less than 1% and can be safely ignored in most practical applications (this error can be reduced further if required). The same comment applies to other cases considered below. The effect of the choice of the number of iterations on the prediction of the integral solution is illustrated in Fig. 2 for the same
Fig. 1. The plots of Ts versus time (a) and Rd versus time (b) using the conventional model (thick solid), integral model based on Eq. (28) (dashed) and linear model (thin solid) for a stationary n-dodecane (Mf = 170 kg/kmole) droplet with an initial radius 5 lm, evaporating in ambient air at a pressure of p = 30 atm = 3000 kPa and temperature 1000 K.
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Fig. 2. The same as Fig. 1 but for different numbers of iterations in the integral solution.
Fig. 3. The same as Fig. 1 but for a droplet with initial radius 50 lm.
1283
1284
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Fig. 4. The same as Figs. 1 and 3 but for a droplet with initial radius 100 lm.
Fig. 5. The same as Fig. 1 but for an ambient gas temperature equal to 800 K.
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
1285
Fig. 6. The same as Figs. 1 and 5 but for an ambient gas temperature equal to 1200 K.
Fig. 7. The plots of T versus n = R/Rd for a stationary n-dodecane (Mf = 170 kg/ kmole) droplet with initial radius 5 lm, evaporating in ambient air at a pressure of p = 30 atm = 3000 kPa and temperature 1000 K. The moments of time are indicated near the curves. The calculations were performed based on Eq. (28) (derived for constant initial distribution of temperatures inside droplets) and (40) (derived for the general case of the arbitrary distribution of the initial temperature inside the droplet).
case as shown in Fig. 1. This effect is shown only for the times when the deviation between the results predicted by the linear and integral solutions is maximal. For the first iteration, the time evolution of droplet radius is the same as predicted by the conventional model. The deviation of the corresponding droplet temperatures predicted by the integral and linear solutions appears to be quite noticeable. For the fifth iteration the droplet surface temperatures predicted by the integral and linear solutions
practically coincide up to t 0.45 ms. The corresponding plots of Rd(t), predicted by the integral solution, turned out to be closer to those predicted by the linear solution than those predicted by the conventional model. The closeness between the plots predicted by the linear and integral solutions improved as the number of iterations increased. However, even for the 15th iteration the deviation between the results remains visible, although not important for practical applications (cf. Fig. 1). For higher iterations the results are practically indistinguishable from those predicted by the 15th iteration. Interestingly, odd iterations predicted smaller Rd(t) and even iterations predicted larger Rd(t) compared with those predicted by the linear solution. At the qualitative level this could be related to the fact that a faster evaporation rate, assumed for the first iteration (conventional model), leads to a lower droplet surface temperature. At the second iteration, this lower droplet surface temperature leads to a slower evaporation rate etc. As to the computational efficiency of the new integral model, we note that for a PC Xeon 3000 Hz (the calculations were processed on one kernel only) with 2 GB RAM, the conventional approach requires 3586 s to calculate 1191 steps. Once these calculations have been completed, the integral model requires 453 s to calculate 15 iterations. This makes this model potentially suitable for incorporation into computational fluid dynamics (CFD) codes. The results, similar to those shown in Fig. 1, but for droplets with initial radii 50 and 100 lm are shown in Figs. 3 and 4 respectively. As can be seen from these figures, the plots of droplet surface temperatures and radii are largely unaffected by the initial droplet radii. This agrees with similar results reported in [3] (see Figs. 4–6 in that paper). The results, similar to those shown in Fig. 1, but for the gas temperatures 800 and 1200 K are shown in Figs. 5 and 6 respectively.
1286
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Comparing Figs. 1 and 5 one can see that a decrease in gas temperature by 20% leads to an increase in the evaporation time predicted by all models by more than 50%. The results predicted by the integral and linear solutions for gas temperature equal to 800 K, are much closer to those predicted by the conventional model than is the case for gas temperature equal to 1200 K. In the case of droplet surface temperature, they are hardly distinguishable, and we had to use an insert figure to zoom the difference between them. As in the case shown in Fig. 1, the results predicted by the integral and linear models for gas temperatures 800 and 1200 K practically coincide. Comparing Figs. 1 and 6 one can see that an increase in gas temperature by 20% leads to a decrease in the evaporation time predicted by all models of more than 25%. The deviation between the results predicted by the integral and linear models and those predicted by the conventional model is noticeably larger in the case shown in Fig. 6 than in the case shown in Figs. 1 and 5, in agreement with similar results reported in [3] (see Fig. 7 in that paper). Also, for the case shown in Fig. 6, the deviation between the results predicted by the integral and linear models is much more visible compared with the previous plots, especially at the final stage of droplet evaporation. The plots of T versus n = R/Rd for T(R, t = 0) Td0(R) = 300 K, Rd0 = 5 lm, Tg = 1000 K and various moments of time (indicated near the curves) are shown in Fig. 7. The calculations were performed based on Eq. (28) (the case of constant Td0(R)) and Eq. (40) (the general case of arbitrary Td0(R)). As can be seen from this figure, the predictions of Eqs. (28) and (40) coincide. This shows the correctness of both approaches to the problem. In agreement with the earlier models described in [2] and [3], Fig. 7 shows that the gradient of temperature inside droplets cannot be ignored at least at the initial stage of droplet heating and evaporation. The effect of non-constant initial distribution of droplet temperature on the time and space evolution of this distribution is illustrated in Fig. 8. Two cases are compared in this figure. In both cases, the initial droplet radii are assumed to be equal to 5 lm, and gas temperature is assumed to be constant and equal to Tg = 1000 K. In the first case the initial distribution of temperature was assumed to be independent on R (or n) and equal to 300 K, as in the case shown in Fig. 7. In this case the analysis was based on Eq. (28). In the second case, the initial distribution of droplet temperature was approximated as
T d0 ðRÞ ¼ 300 þ 10ðR=Rd0 Þ2 ¼ 300 þ 10ðnÞ2 ;
ð41Þ
and the analysis was based on Eq. (40).
Comparing the plots referring to both cases, shown in Fig. 8, one can see that these plots visibly converge with time. This can be related to the fact that increased droplet surface temperature in the general case leads to decreased convective heating of droplets. Hence the droplet surface temperature increases more slowly in the general case than in the case of constant initial temperature inside droplets. We appreciate that the errors associated with the conventional assumption that the droplet radii remain constant during the timestep can be comparable with or even larger than those associated with other effects, including uncertainties in gas temperature measurements, convection heat transfer coefficient approximations and effect of interactions between droplets in realistic sprays. The importance of the latter effect is discussed in our recent papers [9,10], but its analysis lies beyond the scope of this paper. 8. Conclusions Two new solutions to the heat conduction equation, describing transient heating of an evaporating droplet, are suggested, assuming that the time evolution of droplet radius Rd(t) is known. The initial droplet temperature is assumed to be constant or allowed to change with the distance from the droplet centre. The results turned out to be the simplest in the first case and the main focus of our analysis is upon these. Since Rd(t) depends on the time evolution of the droplet temperature, an iterative process is required. Firstly, the time evolution of Rd(t) is obtained using the conventional approach, when it remains constant during the timestep, but changes from one timestep to another. The droplet surface temperature in this case is obtained from the analytical solution to the heat conduction equation inside the droplet. It is assumed that this droplet is heated by convection from the ambient gas, and its radius remains constant during the timestep. Then these values of Rd(t) are used in the new solutions to obtain updated values of time evolution of the distribution of temperatures inside the droplet and on its surface. These new values of droplet temperature are used to update the function Rd(t). This process continues until convergence is achieved, which typically takes place after about 15 iterations. The results of the calculations of droplet surface temperature, using this approach, are compared with the results obtained using the previously suggested approach when the droplet radius was assumed to be a linear function of time during individual timesteps for typical Diesel engine-like conditions. For sufficiently small timesteps the time evolutions of droplet surface temperatures and radii, predicted by both approaches coincide. This suggests that both approaches are correct and valid. Similarly to the case when droplet radius is assumed to be a linear function of time during the timestep, the new solution predicts lower droplet temperatures and slower evaporation when the effects of the reduction of Rd are taken into account. It is shown that in the case of constant droplet initial temperature, models both taking and not taking into account the changes in initial droplet temperature with the distance from the droplet centre, predict the same results. This suggests that both models are correct. It is shown that the temperatures predicted by the models based on the assumption of constant initial droplet temperature, and the one taking into account the increase in this temperature with the distance from the droplet centre, tend to converge with time. Acknowledgements
Fig. 8. The same as Fig. 7 but for the general solution (Eq. (40) applied to the case when the initial distribution of temperature inside the droplet is given by Eq. (41)).
The authors are grateful to the European Regional Development Fund Franco-British INTERREG IVA (Project C5, Reference 4005) and the Royal Society (joint project with Russia, Reference JP 071790) for the financial support of this project.
1287
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Appendix A
where
A1. Numerical solution of the integral Eq. (27)
pffiffiffiffiffiffi xðt n ; sn ÞDt : g n ¼ Xðt n ; sn Þ Dt þ 2
Remembering Eqs. (23) and (25) we can rewrite Eq. (27) as:
mðtÞ þ
Z
t
0
1
mðsÞ pffiffiffiffiffiffiffiffiffiffiffi Xðt; sÞ þ xðt; sÞ ds ¼ 2l0 ðtÞ; ts
ðA11Þ
Rd ðtÞ Rd ðsÞ pffiffiffiffi þ jHðtÞ ts p 2 j " # ðRd ðtÞ ðRd ðsÞÞ2 ; exp 4jðt sÞ 1
Xðt; sÞ ¼ pffiffiffiffi pffiffiffiffi
Rd ðtÞ þ Rd ðsÞ pffiffiffiffi jHðtÞ ts " # ðRd ðtÞ þ ðRd ðsÞÞ2 : exp 4jðt sÞ 1
ðA12Þ
1
ðA13Þ
Functions X(t,psffiffiffiffiffiffiffiffiffiffi ) and ffi x(t, s) are continuous for s 2 [0, t]. Hence, the singularity 1= t s of the kernel in Eq. (A11) is presented in an explicit form. We look for the solution of Eq. (A11) for t 2 ½0; ^t, where ^t 6 t e is an arbitrary, but fixed positive constant. Let Dt ¼ ^t=N and tn = nDt, where N is the total number of timesteps, n = 0, 1, . . . , N is the number of the current timestep. Note that t0 = 0 and t N ¼ ^t. Discretisation of Eq. (A11) gives:
mðtn Þ þ
Xðt n ; sÞ ffi þ xðt n ; sÞ ds ¼ 2l0 ðtn Þ; mðsÞ pffiffiffiffiffiffiffiffiffiffiffiffi tn s t j1
j¼1
ðA17Þ
n1 1 X mðtj Þ þ mðtj1 Þ 1 þ g n j¼1 2 1 þ gn ( ) Xðt n ; sj Þ 2l ðtn Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðtn ; sj Þ Dt þ 0 : 1 þ gn t n sj
mðtn Þ ¼
xðt; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pðt sÞ 2 j
n Z X
¼ 2l0 ðt n Þ;
where n = 1, 2, . . . , N; Eq. (A17) can be rewritten in an alternative form:
where:
1
Eqs. (A15) and (A16) allow us to rewrite Eq. (A14) in the form: ( ) n1 X mðtj Þ þ mðtj1 Þ Xðtn ; sj Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðt n ; sj Þ Dt mðtn Þð1 þ g n Þ þ mðtn1 Þg n þ 2 t n sj j¼1
tj
ðA14Þ
mðtn1 Þg n
ðA18Þ
Note that gn is a continuous function of Dt, and gn ? 0 when Dt ? +0. Hence, 1 + gn – 0 for a sufficiently small timestep Dt. For n = 1 the sum in Eq. (A18) is equal to zero and m(t0) is a known constant (see above). This allows us to calculate m (t1) in an explicit form from Eq. (A18). As soon as the value of m(t1) is found we can use Eq. (A18) to find m(t2) etc. At the nth timestep, Eq. (A18) is used for calculating m(tn), based on the values of m(t0), m(t1), . . . , m(tn1) calculated at the previous timesteps. At this stage all terms in the sum Pn1 j¼1 are already known. In the limiting case when Rd(t) = Rd(s) = Rd = const, functions X(t, s) and x(t, s) are simplified to:
1 pffiffiffiffi
Xðt; sÞ ¼ pffiffiffiffi jHðtÞ;
p
ðA19Þ
" # ( ) pffiffiffiffi 1 Rd R2d : xðt; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jHðtÞ exp jðt sÞ pðt sÞ jðt sÞ ðA110Þ
where n = 0, 1, . . . , N. Note that
mðt0 Þ ¼ mð0Þ ¼ 2l0 ð0Þ
Appendix B
is the known constant derived in C. The first (n 1) terms in the sum in Eq. (A14) can be approximated as:
Z
tj
The integrals in Eqs. (28) and (40) have the same type of integrable singularity as the integral in Eq. (22). The following analysis will focus on the latter equation which will enable us to simplify the notation. Let us rewrite this equation as:
Xðt n ; sÞ ffi þ xðt n ; sÞ ds mðsÞ pffiffiffiffiffiffiffiffiffiffiffiffi tn s t j1 (
)
mðtj Þ þ mðtj1 Þ Xðtn ; sj Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðt n ; sj Þ Dt; 2 tn sj
ðA15Þ
v ðR; ^tÞ ¼
where j ¼ 1; 2; . . . ; n 1; sj ¼ t j This approximation is valid since all functions in the integrand are continuous, and we look for the solution in the class of continuous functions (m(t) should be continuous for t P 0). In this approximation the known functions are taken at the points s = sj (middle of the time interval [tj1, tj]), while the unknown function is taken as an arithmetic mean of its values at the times tj1 and tj. The in the sum in Eq. (A14) has an integrable singularplast ffiffiffiffiffiffiffiffiffiffiterm ffi ity 1= t s when s ? t 0 (recall that functions X(t, s) and x(t, s) are continuous for s 2 [0, t]). This allows us to approximate this term as: Dt . 2
tn Xðt n ; sÞ ffi þ xðt n ; sÞ ds mðsÞ pffiffiffiffiffiffiffiffiffiffiffiffi tn s t n1 Z tn mðtn Þ þ mðtn1 Þ ds pffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðt n ; sn ÞDt Xðtn ; sn Þ 2 tn s t n1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mðtn Þ þ mðtn1 Þ 2Xðtn ; sn Þ t n tn1 þ xðt n ; sn ÞDt ¼ 2 ¼ ðmðt n Þ þ mðtn1 ÞÞg n ;
B.1. Numerical calculation of the improper integrals in Eqs. (28) and (40)
Z
N Z X j¼1
tj
mðsÞGð^t; s; RÞds;
ðA21Þ
t j1
where ^t ¼ t N ; t n ¼ nDt; n ¼ 0; 1; 2; . . . ; N; Dt ¼ ^t=N. In all integrals we can replace m(s) with the average values over the corresponding time interval (m(tj1) + m(tj))/2. Moreover, in all integrals, except the last one, we can replace Gð^t; s; RÞ with Gð^t; sj ; RÞ, where sj = (tj1 + tj)/2. As a result, Eq. (A21) can be presented in a more explicit form:
v ðR; ^tÞ ¼
N1 X mðtj1 Þ þ mðtj Þ ^ Gðt; sj ; RÞDt 2 j¼1 Z mðtN1 Þ þ mðtN Þ tN ^ Gðt; s; RÞds: þ 2 t N1
ðA22Þ
Firstly we assume that an a priori chosen R is not equal to Rd ð^tÞ. In this case Gð^t; s; RÞ, as defined by Eq. (23), approaches 0, when s ! ^t 0. Hence the singularity in the integrand is not present and the last timestep can be treated as in all the previous timesteps. This allows us to simplify Eq. (A22) to:
ðA16Þ
v ðR; ^tÞ ¼
N X mðtj1 Þ þ mðtj Þ ^ Gðt; sj ; RÞDt: 2 j¼1
ðA23Þ
1288
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
In the case when R ¼ Rd ð^tÞ the first exponent in Eq. (23) tends to 1 when s ! ^t 0. This leads to a singularity ð^t sÞ1=2 in the integrand in Eq. (A22). As a result, the integral in this equation can be presented as:
Z
tN
Gð^t; sN ; RÞds ¼
t N1
" # pffiffiffiffi ( j ðR Rd ðsN ÞÞ2 pffiffiffiffi exp 2 p 4jð^t sN Þ " # ) Z tN ds ðR þ Rd ðsN ÞÞ2 Dt pffiffiffiffiffiffiffiffiffiffiffi exp pffiffiffiffiffiffiffiffiffiffiffiffiffi ^t sN 4jð^t sN Þ ^t s t N1 " # pffiffiffiffi ( j pffiffiffi ðR Rd ðsN ÞÞ2 2 exp ¼ pffiffiffiffiffiffiffi 2jDt 2p " #) ðR þ Rd ðsN ÞÞ2 pffiffiffiffiffiffi Dt : ðA24Þ exp 2jDt
The latter equation allows us to simplify the equation for v ðR; ^tÞ for R ¼ Rd ð^tÞ to: pffiffiffiffi N 1 X mðtj1 Þ þ mðtj Þ ^ ðmðt N1 Þ þ mðt N ÞÞ j pffiffiffiffiffiffiffi v ðR; ^tÞ ¼ Gðt; sj ; RÞDt þ 2 2 2p j¼1 ( " # " #) 2 pffiffiffi ðR Rd ðsN ÞÞ ðR þ Rd ðsN ÞÞ2 pffiffiffiffiffiffi 2 exp Dt : exp 2jDt 2jDt ðA25Þ
Appendix C C.1. Derivation of the expression for l0(0) Having substituted Eq. (37) into Eq. (36) and integrating by parts we obtain:
U 0R ðt; RÞR¼R
d ðtÞ
¼
Z
Reff
0
@G1 ðt; R; fÞ ðfT d0 ðfÞÞ df @f R¼Rd ðtÞ
¼ ðfT do ðfÞÞ½G2 ðt; Rd ðtÞ; fÞjf¼Reff Z Reff ½G2 ðt; R; fÞjR¼Rd ðtÞ ðfT d0 ðfÞÞ0f df; 0
ðA31Þ
where
G2 ðt; R; fÞ ¼
1
½G0 ðt; RðtÞ fÞ þ G0 ðt; RðtÞ þ fÞ " ! !# 1 ðR fÞ2 ðR þ fÞ2 ¼ pffiffiffiffiffiffiffiffiffi exp þ exp : 4jt 4jt 2 pjt
j
One can see that the first term on the right hand side of Eq. (A31) approaches zero when t ? +0 (there is no singularity at t = 0) since Reff Rd0 > 0. Hence, remembering Eq. (33), we obtain:
U 0R ðt; RÞR¼R
d ðtÞ
! ðfT do ðfÞÞ0f f¼R ; d0
ðA32Þ
when t ? +0. Following the same approach and remembering Eq. (33), we obtain:
Uðt; RÞjR¼Rd ðtÞ ! ðfT do ðfÞÞjf¼Rd0 ¼ Rd0 T d0 ðRd0 Þ;
ðA33Þ
when t ? +0. Remembering Eqs. (A32) and (A33), we obtain Eq. (39) as a limiting case of Eq. (35) when t ? +0. References [1] B. Abramzon, W.A. Sirignano, Droplet vaporization model for spray combustion calculations, Int. J. Heat Mass Transfer 32 (1989) 1605–1668. [2] S.S. Sazhin, Advanced models of fuel droplet heating and evaporation, Prog. Energy Combust. Sci. 32 (2006) 162–214. [3] S.S. Sazhin, P.A. Krutitskii, I.G. Gusev, M.R. Heikal, Transient heating of an evaporating droplet, Int. J. Heat Mass Transfer 53 (2010) 2826–2836. [4] S.S. Sazhin, I.N. Shishkova, I.G. Gusev, A. Elwardany, P.A. Krutitskii, M.R. Heikal, Fuel droplet heating and evaporation: new hydrodynamic and kinetic models, in: Proceedings of the 14th International Heat Transfer Conferences, Washington, paper IHTC14-22320, 8–13 August 2010. [5] B. Abramzon, S.S. Sazhin, Convective vaporization of fuel droplets with thermal radiation absorption, Fuel 85 (2006) 32–46. [6] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, 1986. [7] E.M. Kartashov, Analytical Methods in the Thermal Conductivity Theory of Solids, Vyshaya Shkola, 2001 (in Russian). [8] A.N. Tikhonov, A.A. Samarsky, Equations of Mathematical Physics, Nauka Publishing House, 1972 (in Russian). [9] S.S. Sazhin, A. Elwardany, P.A. Krutitskii, G. Castanet, F. Lemoine, E.M. Sazhina, M.R. Heikal, A simplified model for bi-component droplet heating and evaporation, Int. J. Heat Mass Transfer 53 (2010) 4495–4505. [10] T. Kristyadi, V. Deprédurand, G. Castanet, F. Lemoine, S.S. Sazhin, A. Elwardany, E.M. Sazhina, M.R. Heikal, Monodisperse monocomponent fuel droplet heating and evaporation, Fuel 89 (2010) 3995–4001.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Transient heating of an evaporating droplet with presumed time evolution of its radius S.S. Sazhin a,⇑, P.A. Krutitskii b, I.G. Gusev a, M.R. Heikal a a b
Sir Harry Ricardo Laboratories, Centre for Automotive Engineering, School of Computing, Engineering and Mathematics, University of Brighton, Brighton BN2 4GJ, UK Keldysh Institute for Applied Mathematics, Department 4, Miusskaya Sq. 4, Moscow 125047, Russia
a r t i c l e
i n f o
Article history: Received 8 April 2010 Accepted 14 October 2010 Available online 1 December 2010 Keywords: Droplets Diesel fuel Heating Evaporation Moving boundary Analytical solution Stefan problem
a b s t r a c t New solutions to the heat conduction equation, describing transient heating of an evaporating droplet, are suggested, assuming that the time evolution of droplet radius Rd(t) is known. The initial droplet temperature is assumed to be constant or allowed to change with the distance from the droplet centre. Since Rd(t) depends on the time evolution of the droplet temperature, an iterative process is required. Firstly, the time evolution of Rd(t) is obtained using the conventional approach, when it remains constant during the timestep, but changes from one timestep to another. Then these values of Rd(t) are used in the new solutions to obtain updated values of the time evolution of the distribution of temperatures inside the droplet and on its surface. These new values of droplet temperature are used to update the function Rd(t). This process continues until convergence is achieved, which typically takes place after about 15 iterations. The results of these calculations are compared with the results obtained using the previously suggested approach when the droplet radius was assumed to be a linear function of time during individual timesteps for typical Diesel engine-like conditions. For sufficiently small timesteps the time evolutions of droplet temperatures and radii predicted by both approaches coincided. This suggests that both approaches are correct and valid. Similarly to the case when droplet radius is assumed to be a linear function of time during the timestep, the new solutions predict lower droplet temperature and slower evaporation when the effects of the reduction of Rd are taken into account. It is shown that in the case of constant droplet initial temperature, models both taking and not taking into account the changes in the initial droplet temperature with the distance from the droplet centre predict the same results. This indicates that both models are correct. Crown Copyright Ó 2010 Published by Elsevier Ltd. All rights reserved.
1. Introduction The conventional approach to the modelling of droplet heating and evaporation is based on the assumption that droplet radius remains constant during each timestep, but changes from one timestep to another due to the evaporation process [1,2]. In our previous paper [3], we took into account the time evolution of droplet radius Rd(t) assuming that Rd(t) is a linear function of time during the timestep. It was shown that this approach not only allows the reduction of the number of timesteps required for the calculation, but also leads to the prediction of lower droplet temperatures and longer evaporation times compared with the conventional approach. These new and in practice important results encouraged us to investigate this problem further and relax our assumption that Rd(t) is a linear function of time. The results of
⇑ Corresponding author. Tel.: +44 0 1273642677; fax: +44 0 1273642309. E-mail address: [email protected] (S.S. Sazhin).
our investigation are presented in this paper. Some preliminary results and their analysis were presented in [4]. Basic equations and approximations are presented and discussed in the next section. In Section 3, the solution based on the assumption that Rd(t) is a linear function of time, obtained in [3], is reproduced. The new solution, based on the assumption that Rd(t) is an arbitrary (but twice continually differentiable) function of time, but the initial droplet temperature Td0 does not depend on the distance from the droplet centre R, is derived in Section 4. In Section 5, the most general solution for arbitrary Rd(t) and Td0(R) is presented. The application of the new solutions to calculation of Rd(t) and T(R, t) of heated and evaporated droplets is described in Section 6. The results of calculation of droplet temperatures and radii for Diesel engine-like conditions using the conventional approach (Rd = const during the timestep), the previously suggested approach (Rd is the linear function of t during the timestep) and the new approach are presented, compared and discussed in Section 7. The main results of the paper are summarised in Section 8.
0017-9310/$ - see front matter Crown Copyright Ó 2010 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2010.10.018
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
1279
Nomenclature BM c fn G h h0 H(t) k L M(t) qn R Rd t T u U(t, R)
v vn W Y
Greek symbols evaporation rate defined by Eq. (9) Hn function defined by Eq. (11) j thermal diffusivity k, kn eigenvalues defined by Eq. (15) l(t) function defined by Eq. (5) l0(t) function defined by Eq. (12) m(t) function defined by Eq. (22) q density x(t, s) function defined by Eq. (A.13) X(t, s) function defined by Eq. (A.12)
Spalding mass transfer number specific heat capacity function introduced in Eq. (11) function introduced in Eq. (22) convection heat transfer coefficient R0 ðtÞR ðtÞ hðtÞ Rd ðtÞ 1 d 2jd kl function introduced in Eq. (5) thermal conductivity specific heat of evaporation function introduced in Eq. (3) function introduced in Eq. (16) distance from the droplet centre droplet radius time temperature TR function defined by Eq. (30) u R Td0 or function defined by Eq. (34) function defined by Eq. (13) function defined by Eq. (17) mass fraction
a
Subscripts d droplet e evaporation g gas l liquid n timesteps in the integration p constant pressure s surface v vapour 0 initial
for t 2 (0, te), R 2 (0, Rd(t)) with the boundary conditions
2. Basic equations and approximations
Let us assume that a stationary evaporating droplet is immersed into a homogeneous hot gas at temperature Tg. The droplet is heated by convection with the convection heat transfer coefficient h(t) depending on time t and droplet radius Rd(t), and cooled down due to evaporation. Rd(t) is a twice continuously differentiable function of time in the range 0 6 t 6 te, where te is the evaporation time. Both Rd(t) and h(t) are assumed to be known. Effects of thermal radiation are ignored, which is justified for typical Diesel fuel droplets and gas temperatures about or less than 1000 K, when the effect of radiation on droplet evaporation time is typically less than 1% [5]. The changes in the droplet temperatures are described by the heat conduction equation in the form [6,7]:
@T @ 2 T 2 @T ¼j þ @t @R2 R @R
! ð1Þ
for 0 < t < te, 0 < R < Rd(t), where j = kl/(cl ql) = const, kl is the liquid thermal conductivity, cl is the liquid specific heat capacity, ql is the liquid density, R is the distance from the centre of the droplet. Eq. (1) needs to be solved subject to the boundary condition:
@T kl ¼ hT g þ ql LR_ d ðtÞ; þ hT @R R¼Rd ðtÞ
ð2Þ
T is finite and continuous at R ? +0, Ts T(Rd(t), t) is the droplet’s surface temperature, L is the latent heat of evaporation. We took into account that R_ d ðtÞ dRd =dt 6 0. Eq. (2) is the energy balance condition at R = Rd(t). The initial condition is taken in the form T(t = 0) = Td0(R), where 0 6 R 6 Rd0 = Rd(t = 0). Let us rewrite boundary condition (2) in the form
@T h h q ¼ T g þ l LR_ d ðtÞ MðtÞ; þ T @R kl kl kl R¼Rd ðtÞ
ð3Þ
and introduce the new variable u = TR. Using this new variable, we can rewrite Eq. (1) as
@u @2u ¼j 2 @t @R
ð4Þ
@u þ HðtÞu ¼ lðtÞ; @R R¼Rd ðtÞ
ujR¼0 ¼ 0;
ð5Þ ð6Þ
where
HðtÞ ¼
hðtÞ 1 ; kl Rd ðtÞ
lðtÞ ¼ MðtÞRd ðtÞ
for t 2 (0, te). The initial condition is
uðRÞjt¼0 ¼ RT d0 ðRÞ
ð7Þ
for R 2 (0, Rd0). 3. Solution for the case when Rd(t) is a linear function Let us assume that Rd(t) is the linear function of t during the time step
Rd ðtÞ ¼ Rd0 ð1 þ atÞ;
ð8Þ
where Rd0 is the initial droplet radius. Remembering the physical background to the problem, the expression for a can be presented as [2]:
a¼
kg lnð1 þ BM Þ
ql cpg R2d0
;
ð9Þ
where BM = Yvs/(1 Yvs) is the Spalding mass transfer number, Yvs is the mass fraction of fuel vapour near the droplet surface, kg and cpg are gas thermal conductivity and specific heat capacity respectively. For Rd(t) defined by (8) the solution to Eq. (1), subject to the above boundary and initial conditions, was obtained in the form [3]:
" # 1 aRd0 R2 TðRÞ ¼ pffiffiffiffiffiffiffiffiffiffiffi exp 4jRd ðtÞ R Rd ðtÞ " # 1 X R l ðtÞ R ; Hn ðtÞ sin kn þ 0 Rd ðtÞ 1 þ h0 Rd ðtÞ n¼1
ð10Þ
1280
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
where
"
Hn ðtÞ ¼ Hn ð0Þ exp " exp
k2n R2d0
j a
k2n R2d0
j a
1 1 1 þ at
#
1 1 1 þ at 1 þ as
pffiffiffiffiffiffiffiffiffiffiffi
l0 ðtÞ l~ ðtÞ Rd ðtÞ exp
where
þ fn
# ds;
Z 0
t
Gðt; s; RÞ ¼ G0 ðt s; R Rd ðsÞÞ G0 ðt s; R þ Rd ðsÞÞ; pffiffiffiffi j x2 : G0 ðt; xÞ ¼ pffiffiffiffiffiffi exp 4jt 2 pt
dl0 ðsÞ ds ð11Þ
0 Rd ðtÞRd ðtÞ ; 4j
ð12Þ
l~ ðtÞ ¼ MðtÞR2d ðtÞ; MðtÞ ¼ khl T g þ qkll LR_ d ðtÞ; f n ¼ kvsinn kk2nk2 , n
v n ðnÞ ¼ sin kn n
ðn ¼ 1; 2; . . .Þ ! 1 sin 2kn 1 h0 2 ; ¼ 1 1þ 2 kv n k ¼ 2 2 2kn h0 þ k2n
ð13Þ ð14Þ
Expression G(t, s, R) can be presented in an alternative form: ( " # " #) pffiffiffiffi j ðR Rd ðsÞÞ2 ðR þ Rd ðsÞÞ2 exp : Gðt; s; RÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 4jðt sÞ 4jðt sÞ 2 pðt sÞ ð23Þ
Note that G(t, s, R = 0) = 0. m(t) is an unknown continuous function to be found later from one of the boundary conditions. Function v(R, t) is known as a single layer heat potential and it has the following properties for any continuous function m(t) [7,8]: (1) (2) (3) (4)
kn are positive solutions to the equation
k cos k þ h0 sin k ¼ 0;
ð15Þ
presented in ascending order,
h0 ¼
hðtÞ R0 ðtÞRd ðtÞ Rd ðtÞ 1 d kl 2j
Z t @ v ðR; tÞ mðtÞ @Gðt; s; RÞ þ ¼ m ð s Þ ds: @R R!Rd ðtÞ0 2 @R 0 R¼Rd ðtÞ
is assumed to be constant during the time step,
Hn ð0Þ ¼ qn þ l0 ð0Þfn ; qn ¼
Z
1 2
kv n k
W 0 ðnÞv n ðnÞdn:
3=2
W 0 ðnÞ ¼ Rd0 nT 0 ðnRd0 Þ exp
R0d ð0ÞRd0 2 n : 4j
ð17Þ
4. Solution for the case of arbitrary Rd(t) but Td0(R) = const The analysis of this section is based on the assumption that Td0(R) = Td0 = const. In this case we can introduce the new variable v = u R Td0 and rearrange Eq. (4) as:
@v @2v ¼j 2 @t @R
ð18Þ
@Gðt; s; RÞ 1 ¼ pffiffiffiffiffiffiffi @R 4 pjðt sÞ3=2 R¼Rd ðtÞ (
# ðRd ðtÞ Rd ðsÞÞ2 ðRd ðtÞ Rd ðsÞÞ exp 4jðt sÞ " #) ðRd ðtÞ þ Rd ðsÞÞ2 : ðRd ðtÞ þ Rd ðsÞÞ exp 4jðt sÞ
@v ¼ l0 ðtÞ; þ HðtÞv @R R¼Rd ðtÞ
ð19Þ
v jR¼0 ¼ 0
ð20Þ
for t 2 (0, te) and the initial condition
Since Rd(t) is a continuously differentiable function, we obtain in the limit s ? t 0:
@Gðt; s; RÞ 1 p ffiffiffiffiffiffiffiffiffiffi ffi : / O @R ts R¼Rd ðtÞ
mðtÞ
þ
Z
t
mðsÞ
0
Z t @Gðt; s; RÞ d s þ HðtÞ mðsÞGðt; s; Rd ðtÞÞds @R 0 R¼Rd ðtÞ
¼ l0 ðtÞ; or
mðtÞ 2
hðtÞ l0 ðtÞ ¼ T d0 HðtÞRd ðtÞT d0 þ lðtÞ ¼ Rd ðtÞT d0 þ lðtÞ: kl We look for the solution to the problem (18)–(21) in the form:
Z
0
t
mðsÞGðt; s; RÞds;
ð22Þ
ð26Þ
It follows from this equation, that there is an improper integral on the right hand side of Eq. (24). In view of Eqs. (24) and (22) we can rewrite the boundary condition (19) as:
þ
Z 0
t
mðsÞ
( ) @Gðt; s; RÞ þ HðtÞGðt; s ; R ðtÞÞ ds ¼ l0 ðtÞ; d @R R¼Rd ðtÞ
ð21Þ
for R 2 (0, Rd(t)), where
"
ð25Þ
2
for t 2 (0, te) and R 2 (0, Rd(t)) with the boundary conditions
v jt¼0 ¼ 0
This means that for any continuous function m(t) the potential satisfies Eq. (18) and boundary and initial conditions (20) and (21). Choice of the function m(t), satisfying integral Eq. (24), should be made in such a way that the remaining boundary condition (19) is satisfied as well. From Eq. (23) it follows that
v(R, t)
The first step in the derivation of this equation is based on the reduction of the problem (1) and (2) to the problem (4)–(7) with the corresponding initial conditions. Note that the condition h0 = const can be violated for t in the immediate vicinity of te for finite timesteps. This imposes a rather severe restriction on the choice of the timestep just before the droplet fully evaporates. Condition (8) allows us to use the above solution for individual timesteps but not for the whole period of droplet heating and evaporation. In what follows this assumption will be relaxed and we consider arbitrary Rd(t) provided that this function is twice continuously differentiable (natural condition remembering the physical background to the problem).
v ðR; tÞ ¼
ð24Þ
ð16Þ
1
0
It satisfies Eq. (18) for 0 < t < te and 0 < R < Rd(t); It satisfies Conditions (20) and (21); It is continuous at R ? Rd 0; For the derivative @ v (R, t)/@R the following limiting formula is valid:
ð27Þ where G(t, s, R) and its derivative with respect to R are defined by Eqs. (23) and (25). As follows from the definition of G(t, s, Rd(t)) (see Eq. (23)), in the limit s ? t 0 this function has the singularity:
1 Gðt; s; Rd ðtÞÞ / O pffiffiffiffiffiffiffiffiffiffiffi ts
1281
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
(cf. Eq. (26)). Therefore, the integral in Eq. (27) is defined as an improper integral. Eq. (27) is an integral equation of Volterra type. It has a unique continuous solution. A scheme for its numerical solution is described in Appendix A. This solution is then substituted into Eq. (22). The final distribution of temperature inside the droplet can be calculated from the following expression: pffiffiffiffi Z t j pffiffiffiffi 2R p 0 ( " # " #) 2 2 m ð sÞ ðR Rd ðsÞÞ ðR þ Rd ðsÞÞ pffiffiffiffiffiffiffiffiffiffiffi exp exp ds: 4jðt sÞ 4jðt sÞ ts
Tðt; RÞ ¼ T d0 þ
ð28Þ
Details of the numerical calculation of the integral on the right hand side of Eq. (28) are given in Appendix B.
Having substituted Eq. (34) into Eq. (4) and boundary and initial conditions (5)–(7), we obtain problem (18)–(21) for v(t, R) in which
l0 ðtÞ ¼ U 0R ðt; RÞ þ HðtÞUðt; RÞ R¼Rd ðtÞ þ lðtÞ:
ð35Þ
The solution of the latter problem is similar to the one discussed in Section 4. The expression for l0(t) contains
U 0R ðt; RÞR¼R
d ðtÞ
¼
Z
Reff 0
@G1 ðt; R; fÞ ðfT do ðfÞÞ df; @R R¼Rd ðtÞ
ð36Þ
where
@G1 ðt; R; fÞ 1 @ G ¼ ð ðt; R fÞ G ðt; R þ fÞ Þ 0 0 @R j @R R¼Rd ðtÞ R¼Rd ðtÞ 1 @ ¼ ðG0 ðt; R fÞ G0 ðt; R þ fÞÞ ; j @f R¼Rd ðtÞ ð37Þ
5. Solution for the case of arbitrary Rd(t) and Td0(R) Let us assume that an arbitrary twice continuously differentiable function Td0(R) is defined for 0 6 R 6 Rd0. This definition is extended for the whole range 0 6 R < 1:
8 > < T d0 ðRÞ when 0 6 R 6 Rd0 ; T d0 ðRÞ ¼ T out ðRÞ when Rd0 < R 6 Reff ; > : 0 when R > Reff ;
ð29Þ
G0 is the same as in Eq. (22). Remembering the latter equation we ðt;R;fÞ can rewrite the expression for @G1@R in a more explicit form:
( R¼Rd ðtÞ " # @G1 ðt; R; fÞ 1 ðR fÞ2 ¼ ðR fÞ exp p ffiffiffi ffi @R 4jt 4 pðjtÞ3=2 R¼Rd ðtÞ " #) ðR þ fÞ2 : ðR þ fÞ exp 4jt R¼Rd ðtÞ
where
Hence, we obtain an explicit expression for l0(t) in the form:
h io 1n Rd0 T d0 ðRd0 Þ þ ðR Rd0 Þ ðRT d0 ðRÞÞ0R R¼R T out ðRÞ ¼ ; d0 R
l0 ðtÞ ¼ pffiffiffiffi 3=2 4 pðjtÞ
1
Reff is the effective outer radius such that Reff > Rd0. Function Td0(R) defined by Eq. (29) is continuously differentiable in the range 0 6 R 6 Reff. Let us now introduce a new function U(t, R) defined as:
Uðt; RÞ ¼
Z
Reff
ðfT do ðfÞÞG1 ðt; R; fÞdf;
ð30Þ
0
where
G1 ðt; R; fÞ ¼
1
j
ð31Þ
(1) It satisfies Eq. (18) for 0 < t < te and 0 < R < 1; (2) It satisfies the boundary Condition (20) for 0 < t < te; (3) It satisfies the initial condition
RT d0 ðRÞ when 0 6 R 6 Reff ; 0
when R > Reff :
ð32Þ
The latter relation follows from the property of the delta-function:
lim
adelta !1
adelta pffiffiffiffi exp a2delta x2 ¼ dðxÞ: p
ð33Þ
We look for the solution to Eq. (4) in the form:
uðt; RÞ ¼ Uðt; RÞ þ v ðt; RÞ:
0
ðRd ðtÞ fÞ2 4jt
#
"
ð34Þ
i
l0 ð0Þ ¼ ðfT do ðfÞÞ0f f¼Rd0 þ Hð0ÞRd0 T d0 ðRd0 Þ þ lð0Þ:
ð38Þ
ð39Þ
Combining Eqs. (22) and (34) we can present the final solution to our problem in the form:
Tðt; RÞ ¼
Note that G1(t, R = 0) = 0. Function U(t, R) has the following properties [7,8]:
"
ðfT d0 ðfÞÞ ðRd ðtÞ fÞ exp
#) ðRd ðtÞ þ fÞ2 df ðRd ðtÞ þ fÞ exp 4jt ( " # Z Reff HðtÞ ðRd ðtÞ fÞ2 pffiffiffiffiffiffiffiffiffi ðfT d0 ðfÞÞ exp 4jt 2 pjt 0 " #) ðRd ðtÞ þ fÞ2 df þ MðtÞRd ðtÞ: exp 4jt
h
G0(t, x) is the same as in Eq. (22). Remembering the latter equation, the expression G1(t, R) can be presented in an alternative form:
Uðt; RÞjt¼þ0 ¼
(
Reff
In the limit t ? 0+ the expression for l0(t) is simplified to (see Appendix C):
½G0 ðt; R fÞ G0 ðt; R þ fÞ;
( " # " #) 1 ðR fÞ2 ðR þ fÞ2 exp : G1 ðt; RÞ ¼ pffiffiffiffiffiffiffiffiffi exp 4jt 4jt 2 pt j
Z
" ( " # pffiffiffiffi Z t 1 j mðsÞ ðR Rd ðsÞÞ2 pffiffiffiffiffiffiffiffiffiffiffi exp Uðt; RÞ þ pffiffiffiffi R 4jðt sÞ 2 p 0 ts " #) # ðR þ Rd ðsÞÞ2 ds ; exp 4jðt sÞ
ð40Þ
where m(s) is the solution to Eq. (27) with l0(t) given by Eq. (38), U(t, R) is given by Eq. (30). Note that taking into account the initial distribution of temperature along R is absolutely essential when the solution is applied to individual timesteps. In the solution described in the last two sections, however, the same formulae describe the time evolution of droplet temperatures during the whole period of their evaporation. It is anticipated that in this case the effect of the initial distribution of droplet temperatures is not important in most practical applications. Hence, the solution described in Section 4 is expected to be more practically important than the solution described in this section.
1282
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Note that although the analysis presented so far refers to stationary droplets, it can be generalised in a straightforward way to the case of moving droplets, based on the effective thermal conductivity model [1,2]. 6. Implementation of the new solutions into a numerical code In the solutions presented in the last two sections it was assumed that Rd(t) is known. From the point of view of the physical background to the problem, however, Rd(t) depends on the time evolution of the droplet temperature T(R, t), which is the solution to the problem. Hence, an iterative process is required. Firstly, the time evolution of droplet radius Rd(t) is obtained using the conventional approach, when it remains constant during the timestep, but changes from one timestep to another due to the evaporation process. Then these values of Rd(t) are used in the new solutions to obtain updated values of the time evolution of the distribution of temperatures inside the droplet and on its surface T(R, t). These new values of droplet temperature are used to update the function Rd(t). This process is continued until convergence is achieved, which typically takes place after about 15 iterations. 7. Results In Figs. 1–7 we compared the results of calculations of droplet surface temperatures and radii, taking into account the effects of evaporation, using the new integral solution for arbitrary Rd(t) but constant Td0 (Eq. (28)), the previously reported solution, based on the linear approximation of Rd(t) (Eq. (10)), and the conventional approach based on the assumption that droplet radius does not change during the timestep. Droplets are assumed to be those of n-dodecane (Mf = 170 kg/kmole), and ambient air is assumed to
be at the pressure of p = 30 atm = 3000 kPa (typical values for Diesel engine-like conditions). The results of calculations for Rd0 = 5 lm, and ambient air temperature 1000 K are shown in Fig. 1. One thousand timesteps were used in the conventional approach and the approach based on Solution (10). In the integral solution based on Eq. (28) the integral (A21) was approximated as the sum of 100 terms, and up to 15 iterations were used. At the first iteration the time evolution of the droplet radius was assumed to be the same as predicted by the conventional model. As follows from Fig. 1, the results predicted by the integral solution (28) and linear solution (10) practically coincide, which suggests that both approaches are correct and valid. Both these solutions predict lower droplet temperatures and longer evaporation times in agreement with the results reported in [3]. Note that deviations between the predictions of the integral and linear solutions were observed in the immediate vicinity of the time when the droplet completely evaporates. There were obvious numerical problems when we approached this time due to the fact that the time derivative of Rd becomes infinitely large. In practice the extrapolation, based on the assumption that the second derivative of Rd(t) is constant, was used for these times. This leads to small deviations between the predicted evaporation times. In the case shown in Fig. 1, the evaporation times predicted by the conventional model, linear solution, and integral solution were 0.595 ms, 0.622 ms and 0.628 ms respectively. That means that the difference between the evaporation times predicted by the linear and integral solutions was less than 1% and can be safely ignored in most practical applications (this error can be reduced further if required). The same comment applies to other cases considered below. The effect of the choice of the number of iterations on the prediction of the integral solution is illustrated in Fig. 2 for the same
Fig. 1. The plots of Ts versus time (a) and Rd versus time (b) using the conventional model (thick solid), integral model based on Eq. (28) (dashed) and linear model (thin solid) for a stationary n-dodecane (Mf = 170 kg/kmole) droplet with an initial radius 5 lm, evaporating in ambient air at a pressure of p = 30 atm = 3000 kPa and temperature 1000 K.
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Fig. 2. The same as Fig. 1 but for different numbers of iterations in the integral solution.
Fig. 3. The same as Fig. 1 but for a droplet with initial radius 50 lm.
1283
1284
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Fig. 4. The same as Figs. 1 and 3 but for a droplet with initial radius 100 lm.
Fig. 5. The same as Fig. 1 but for an ambient gas temperature equal to 800 K.
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
1285
Fig. 6. The same as Figs. 1 and 5 but for an ambient gas temperature equal to 1200 K.
Fig. 7. The plots of T versus n = R/Rd for a stationary n-dodecane (Mf = 170 kg/ kmole) droplet with initial radius 5 lm, evaporating in ambient air at a pressure of p = 30 atm = 3000 kPa and temperature 1000 K. The moments of time are indicated near the curves. The calculations were performed based on Eq. (28) (derived for constant initial distribution of temperatures inside droplets) and (40) (derived for the general case of the arbitrary distribution of the initial temperature inside the droplet).
case as shown in Fig. 1. This effect is shown only for the times when the deviation between the results predicted by the linear and integral solutions is maximal. For the first iteration, the time evolution of droplet radius is the same as predicted by the conventional model. The deviation of the corresponding droplet temperatures predicted by the integral and linear solutions appears to be quite noticeable. For the fifth iteration the droplet surface temperatures predicted by the integral and linear solutions
practically coincide up to t 0.45 ms. The corresponding plots of Rd(t), predicted by the integral solution, turned out to be closer to those predicted by the linear solution than those predicted by the conventional model. The closeness between the plots predicted by the linear and integral solutions improved as the number of iterations increased. However, even for the 15th iteration the deviation between the results remains visible, although not important for practical applications (cf. Fig. 1). For higher iterations the results are practically indistinguishable from those predicted by the 15th iteration. Interestingly, odd iterations predicted smaller Rd(t) and even iterations predicted larger Rd(t) compared with those predicted by the linear solution. At the qualitative level this could be related to the fact that a faster evaporation rate, assumed for the first iteration (conventional model), leads to a lower droplet surface temperature. At the second iteration, this lower droplet surface temperature leads to a slower evaporation rate etc. As to the computational efficiency of the new integral model, we note that for a PC Xeon 3000 Hz (the calculations were processed on one kernel only) with 2 GB RAM, the conventional approach requires 3586 s to calculate 1191 steps. Once these calculations have been completed, the integral model requires 453 s to calculate 15 iterations. This makes this model potentially suitable for incorporation into computational fluid dynamics (CFD) codes. The results, similar to those shown in Fig. 1, but for droplets with initial radii 50 and 100 lm are shown in Figs. 3 and 4 respectively. As can be seen from these figures, the plots of droplet surface temperatures and radii are largely unaffected by the initial droplet radii. This agrees with similar results reported in [3] (see Figs. 4–6 in that paper). The results, similar to those shown in Fig. 1, but for the gas temperatures 800 and 1200 K are shown in Figs. 5 and 6 respectively.
1286
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Comparing Figs. 1 and 5 one can see that a decrease in gas temperature by 20% leads to an increase in the evaporation time predicted by all models by more than 50%. The results predicted by the integral and linear solutions for gas temperature equal to 800 K, are much closer to those predicted by the conventional model than is the case for gas temperature equal to 1200 K. In the case of droplet surface temperature, they are hardly distinguishable, and we had to use an insert figure to zoom the difference between them. As in the case shown in Fig. 1, the results predicted by the integral and linear models for gas temperatures 800 and 1200 K practically coincide. Comparing Figs. 1 and 6 one can see that an increase in gas temperature by 20% leads to a decrease in the evaporation time predicted by all models of more than 25%. The deviation between the results predicted by the integral and linear models and those predicted by the conventional model is noticeably larger in the case shown in Fig. 6 than in the case shown in Figs. 1 and 5, in agreement with similar results reported in [3] (see Fig. 7 in that paper). Also, for the case shown in Fig. 6, the deviation between the results predicted by the integral and linear models is much more visible compared with the previous plots, especially at the final stage of droplet evaporation. The plots of T versus n = R/Rd for T(R, t = 0) Td0(R) = 300 K, Rd0 = 5 lm, Tg = 1000 K and various moments of time (indicated near the curves) are shown in Fig. 7. The calculations were performed based on Eq. (28) (the case of constant Td0(R)) and Eq. (40) (the general case of arbitrary Td0(R)). As can be seen from this figure, the predictions of Eqs. (28) and (40) coincide. This shows the correctness of both approaches to the problem. In agreement with the earlier models described in [2] and [3], Fig. 7 shows that the gradient of temperature inside droplets cannot be ignored at least at the initial stage of droplet heating and evaporation. The effect of non-constant initial distribution of droplet temperature on the time and space evolution of this distribution is illustrated in Fig. 8. Two cases are compared in this figure. In both cases, the initial droplet radii are assumed to be equal to 5 lm, and gas temperature is assumed to be constant and equal to Tg = 1000 K. In the first case the initial distribution of temperature was assumed to be independent on R (or n) and equal to 300 K, as in the case shown in Fig. 7. In this case the analysis was based on Eq. (28). In the second case, the initial distribution of droplet temperature was approximated as
T d0 ðRÞ ¼ 300 þ 10ðR=Rd0 Þ2 ¼ 300 þ 10ðnÞ2 ;
ð41Þ
and the analysis was based on Eq. (40).
Comparing the plots referring to both cases, shown in Fig. 8, one can see that these plots visibly converge with time. This can be related to the fact that increased droplet surface temperature in the general case leads to decreased convective heating of droplets. Hence the droplet surface temperature increases more slowly in the general case than in the case of constant initial temperature inside droplets. We appreciate that the errors associated with the conventional assumption that the droplet radii remain constant during the timestep can be comparable with or even larger than those associated with other effects, including uncertainties in gas temperature measurements, convection heat transfer coefficient approximations and effect of interactions between droplets in realistic sprays. The importance of the latter effect is discussed in our recent papers [9,10], but its analysis lies beyond the scope of this paper. 8. Conclusions Two new solutions to the heat conduction equation, describing transient heating of an evaporating droplet, are suggested, assuming that the time evolution of droplet radius Rd(t) is known. The initial droplet temperature is assumed to be constant or allowed to change with the distance from the droplet centre. The results turned out to be the simplest in the first case and the main focus of our analysis is upon these. Since Rd(t) depends on the time evolution of the droplet temperature, an iterative process is required. Firstly, the time evolution of Rd(t) is obtained using the conventional approach, when it remains constant during the timestep, but changes from one timestep to another. The droplet surface temperature in this case is obtained from the analytical solution to the heat conduction equation inside the droplet. It is assumed that this droplet is heated by convection from the ambient gas, and its radius remains constant during the timestep. Then these values of Rd(t) are used in the new solutions to obtain updated values of time evolution of the distribution of temperatures inside the droplet and on its surface. These new values of droplet temperature are used to update the function Rd(t). This process continues until convergence is achieved, which typically takes place after about 15 iterations. The results of the calculations of droplet surface temperature, using this approach, are compared with the results obtained using the previously suggested approach when the droplet radius was assumed to be a linear function of time during individual timesteps for typical Diesel engine-like conditions. For sufficiently small timesteps the time evolutions of droplet surface temperatures and radii, predicted by both approaches coincide. This suggests that both approaches are correct and valid. Similarly to the case when droplet radius is assumed to be a linear function of time during the timestep, the new solution predicts lower droplet temperatures and slower evaporation when the effects of the reduction of Rd are taken into account. It is shown that in the case of constant droplet initial temperature, models both taking and not taking into account the changes in initial droplet temperature with the distance from the droplet centre, predict the same results. This suggests that both models are correct. It is shown that the temperatures predicted by the models based on the assumption of constant initial droplet temperature, and the one taking into account the increase in this temperature with the distance from the droplet centre, tend to converge with time. Acknowledgements
Fig. 8. The same as Fig. 7 but for the general solution (Eq. (40) applied to the case when the initial distribution of temperature inside the droplet is given by Eq. (41)).
The authors are grateful to the European Regional Development Fund Franco-British INTERREG IVA (Project C5, Reference 4005) and the Royal Society (joint project with Russia, Reference JP 071790) for the financial support of this project.
1287
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
Appendix A
where
A1. Numerical solution of the integral Eq. (27)
pffiffiffiffiffiffi xðt n ; sn ÞDt : g n ¼ Xðt n ; sn Þ Dt þ 2
Remembering Eqs. (23) and (25) we can rewrite Eq. (27) as:
mðtÞ þ
Z
t
0
1
mðsÞ pffiffiffiffiffiffiffiffiffiffiffi Xðt; sÞ þ xðt; sÞ ds ¼ 2l0 ðtÞ; ts
ðA11Þ
Rd ðtÞ Rd ðsÞ pffiffiffiffi þ jHðtÞ ts p 2 j " # ðRd ðtÞ ðRd ðsÞÞ2 ; exp 4jðt sÞ 1
Xðt; sÞ ¼ pffiffiffiffi pffiffiffiffi
Rd ðtÞ þ Rd ðsÞ pffiffiffiffi jHðtÞ ts " # ðRd ðtÞ þ ðRd ðsÞÞ2 : exp 4jðt sÞ 1
ðA12Þ
1
ðA13Þ
Functions X(t,psffiffiffiffiffiffiffiffiffiffi ) and ffi x(t, s) are continuous for s 2 [0, t]. Hence, the singularity 1= t s of the kernel in Eq. (A11) is presented in an explicit form. We look for the solution of Eq. (A11) for t 2 ½0; ^t, where ^t 6 t e is an arbitrary, but fixed positive constant. Let Dt ¼ ^t=N and tn = nDt, where N is the total number of timesteps, n = 0, 1, . . . , N is the number of the current timestep. Note that t0 = 0 and t N ¼ ^t. Discretisation of Eq. (A11) gives:
mðtn Þ þ
Xðt n ; sÞ ffi þ xðt n ; sÞ ds ¼ 2l0 ðtn Þ; mðsÞ pffiffiffiffiffiffiffiffiffiffiffiffi tn s t j1
j¼1
ðA17Þ
n1 1 X mðtj Þ þ mðtj1 Þ 1 þ g n j¼1 2 1 þ gn ( ) Xðt n ; sj Þ 2l ðtn Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðtn ; sj Þ Dt þ 0 : 1 þ gn t n sj
mðtn Þ ¼
xðt; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pðt sÞ 2 j
n Z X
¼ 2l0 ðt n Þ;
where n = 1, 2, . . . , N; Eq. (A17) can be rewritten in an alternative form:
where:
1
Eqs. (A15) and (A16) allow us to rewrite Eq. (A14) in the form: ( ) n1 X mðtj Þ þ mðtj1 Þ Xðtn ; sj Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðt n ; sj Þ Dt mðtn Þð1 þ g n Þ þ mðtn1 Þg n þ 2 t n sj j¼1
tj
ðA14Þ
mðtn1 Þg n
ðA18Þ
Note that gn is a continuous function of Dt, and gn ? 0 when Dt ? +0. Hence, 1 + gn – 0 for a sufficiently small timestep Dt. For n = 1 the sum in Eq. (A18) is equal to zero and m(t0) is a known constant (see above). This allows us to calculate m (t1) in an explicit form from Eq. (A18). As soon as the value of m(t1) is found we can use Eq. (A18) to find m(t2) etc. At the nth timestep, Eq. (A18) is used for calculating m(tn), based on the values of m(t0), m(t1), . . . , m(tn1) calculated at the previous timesteps. At this stage all terms in the sum Pn1 j¼1 are already known. In the limiting case when Rd(t) = Rd(s) = Rd = const, functions X(t, s) and x(t, s) are simplified to:
1 pffiffiffiffi
Xðt; sÞ ¼ pffiffiffiffi jHðtÞ;
p
ðA19Þ
" # ( ) pffiffiffiffi 1 Rd R2d : xðt; sÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jHðtÞ exp jðt sÞ pðt sÞ jðt sÞ ðA110Þ
where n = 0, 1, . . . , N. Note that
mðt0 Þ ¼ mð0Þ ¼ 2l0 ð0Þ
Appendix B
is the known constant derived in C. The first (n 1) terms in the sum in Eq. (A14) can be approximated as:
Z
tj
The integrals in Eqs. (28) and (40) have the same type of integrable singularity as the integral in Eq. (22). The following analysis will focus on the latter equation which will enable us to simplify the notation. Let us rewrite this equation as:
Xðt n ; sÞ ffi þ xðt n ; sÞ ds mðsÞ pffiffiffiffiffiffiffiffiffiffiffiffi tn s t j1 (
)
mðtj Þ þ mðtj1 Þ Xðtn ; sj Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðt n ; sj Þ Dt; 2 tn sj
ðA15Þ
v ðR; ^tÞ ¼
where j ¼ 1; 2; . . . ; n 1; sj ¼ t j This approximation is valid since all functions in the integrand are continuous, and we look for the solution in the class of continuous functions (m(t) should be continuous for t P 0). In this approximation the known functions are taken at the points s = sj (middle of the time interval [tj1, tj]), while the unknown function is taken as an arithmetic mean of its values at the times tj1 and tj. The in the sum in Eq. (A14) has an integrable singularplast ffiffiffiffiffiffiffiffiffiffiterm ffi ity 1= t s when s ? t 0 (recall that functions X(t, s) and x(t, s) are continuous for s 2 [0, t]). This allows us to approximate this term as: Dt . 2
tn Xðt n ; sÞ ffi þ xðt n ; sÞ ds mðsÞ pffiffiffiffiffiffiffiffiffiffiffiffi tn s t n1 Z tn mðtn Þ þ mðtn1 Þ ds pffiffiffiffiffiffiffiffiffiffiffiffiffi þ xðt n ; sn ÞDt Xðtn ; sn Þ 2 tn s t n1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mðtn Þ þ mðtn1 Þ 2Xðtn ; sn Þ t n tn1 þ xðt n ; sn ÞDt ¼ 2 ¼ ðmðt n Þ þ mðtn1 ÞÞg n ;
B.1. Numerical calculation of the improper integrals in Eqs. (28) and (40)
Z
N Z X j¼1
tj
mðsÞGð^t; s; RÞds;
ðA21Þ
t j1
where ^t ¼ t N ; t n ¼ nDt; n ¼ 0; 1; 2; . . . ; N; Dt ¼ ^t=N. In all integrals we can replace m(s) with the average values over the corresponding time interval (m(tj1) + m(tj))/2. Moreover, in all integrals, except the last one, we can replace Gð^t; s; RÞ with Gð^t; sj ; RÞ, where sj = (tj1 + tj)/2. As a result, Eq. (A21) can be presented in a more explicit form:
v ðR; ^tÞ ¼
N1 X mðtj1 Þ þ mðtj Þ ^ Gðt; sj ; RÞDt 2 j¼1 Z mðtN1 Þ þ mðtN Þ tN ^ Gðt; s; RÞds: þ 2 t N1
ðA22Þ
Firstly we assume that an a priori chosen R is not equal to Rd ð^tÞ. In this case Gð^t; s; RÞ, as defined by Eq. (23), approaches 0, when s ! ^t 0. Hence the singularity in the integrand is not present and the last timestep can be treated as in all the previous timesteps. This allows us to simplify Eq. (A22) to:
ðA16Þ
v ðR; ^tÞ ¼
N X mðtj1 Þ þ mðtj Þ ^ Gðt; sj ; RÞDt: 2 j¼1
ðA23Þ
1288
S.S. Sazhin et al. / International Journal of Heat and Mass Transfer 54 (2011) 1278–1288
In the case when R ¼ Rd ð^tÞ the first exponent in Eq. (23) tends to 1 when s ! ^t 0. This leads to a singularity ð^t sÞ1=2 in the integrand in Eq. (A22). As a result, the integral in this equation can be presented as:
Z
tN
Gð^t; sN ; RÞds ¼
t N1
" # pffiffiffiffi ( j ðR Rd ðsN ÞÞ2 pffiffiffiffi exp 2 p 4jð^t sN Þ " # ) Z tN ds ðR þ Rd ðsN ÞÞ2 Dt pffiffiffiffiffiffiffiffiffiffiffi exp pffiffiffiffiffiffiffiffiffiffiffiffiffi ^t sN 4jð^t sN Þ ^t s t N1 " # pffiffiffiffi ( j pffiffiffi ðR Rd ðsN ÞÞ2 2 exp ¼ pffiffiffiffiffiffiffi 2jDt 2p " #) ðR þ Rd ðsN ÞÞ2 pffiffiffiffiffiffi Dt : ðA24Þ exp 2jDt
The latter equation allows us to simplify the equation for v ðR; ^tÞ for R ¼ Rd ð^tÞ to: pffiffiffiffi N 1 X mðtj1 Þ þ mðtj Þ ^ ðmðt N1 Þ þ mðt N ÞÞ j pffiffiffiffiffiffiffi v ðR; ^tÞ ¼ Gðt; sj ; RÞDt þ 2 2 2p j¼1 ( " # " #) 2 pffiffiffi ðR Rd ðsN ÞÞ ðR þ Rd ðsN ÞÞ2 pffiffiffiffiffiffi 2 exp Dt : exp 2jDt 2jDt ðA25Þ
Appendix C C.1. Derivation of the expression for l0(0) Having substituted Eq. (37) into Eq. (36) and integrating by parts we obtain:
U 0R ðt; RÞR¼R
d ðtÞ
¼
Z
Reff
0
@G1 ðt; R; fÞ ðfT d0 ðfÞÞ df @f R¼Rd ðtÞ
¼ ðfT do ðfÞÞ½G2 ðt; Rd ðtÞ; fÞjf¼Reff Z Reff ½G2 ðt; R; fÞjR¼Rd ðtÞ ðfT d0 ðfÞÞ0f df; 0
ðA31Þ
where
G2 ðt; R; fÞ ¼
1
½G0 ðt; RðtÞ fÞ þ G0 ðt; RðtÞ þ fÞ " ! !# 1 ðR fÞ2 ðR þ fÞ2 ¼ pffiffiffiffiffiffiffiffiffi exp þ exp : 4jt 4jt 2 pjt
j
One can see that the first term on the right hand side of Eq. (A31) approaches zero when t ? +0 (there is no singularity at t = 0) since Reff Rd0 > 0. Hence, remembering Eq. (33), we obtain:
U 0R ðt; RÞR¼R
d ðtÞ
! ðfT do ðfÞÞ0f f¼R ; d0
ðA32Þ
when t ? +0. Following the same approach and remembering Eq. (33), we obtain:
Uðt; RÞjR¼Rd ðtÞ ! ðfT do ðfÞÞjf¼Rd0 ¼ Rd0 T d0 ðRd0 Þ;
ðA33Þ
when t ? +0. Remembering Eqs. (A32) and (A33), we obtain Eq. (39) as a limiting case of Eq. (35) when t ? +0. References [1] B. Abramzon, W.A. Sirignano, Droplet vaporization model for spray combustion calculations, Int. J. Heat Mass Transfer 32 (1989) 1605–1668. [2] S.S. Sazhin, Advanced models of fuel droplet heating and evaporation, Prog. Energy Combust. Sci. 32 (2006) 162–214. [3] S.S. Sazhin, P.A. Krutitskii, I.G. Gusev, M.R. Heikal, Transient heating of an evaporating droplet, Int. J. Heat Mass Transfer 53 (2010) 2826–2836. [4] S.S. Sazhin, I.N. Shishkova, I.G. Gusev, A. Elwardany, P.A. Krutitskii, M.R. Heikal, Fuel droplet heating and evaporation: new hydrodynamic and kinetic models, in: Proceedings of the 14th International Heat Transfer Conferences, Washington, paper IHTC14-22320, 8–13 August 2010. [5] B. Abramzon, S.S. Sazhin, Convective vaporization of fuel droplets with thermal radiation absorption, Fuel 85 (2006) 32–46. [6] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, 1986. [7] E.M. Kartashov, Analytical Methods in the Thermal Conductivity Theory of Solids, Vyshaya Shkola, 2001 (in Russian). [8] A.N. Tikhonov, A.A. Samarsky, Equations of Mathematical Physics, Nauka Publishing House, 1972 (in Russian). [9] S.S. Sazhin, A. Elwardany, P.A. Krutitskii, G. Castanet, F. Lemoine, E.M. Sazhina, M.R. Heikal, A simplified model for bi-component droplet heating and evaporation, Int. J. Heat Mass Transfer 53 (2010) 4495–4505. [10] T. Kristyadi, V. Deprédurand, G. Castanet, F. Lemoine, S.S. Sazhin, A. Elwardany, E.M. Sazhina, M.R. Heikal, Monodisperse monocomponent fuel droplet heating and evaporation, Fuel 89 (2010) 3995–4001.