ILASS – Europe 2010, 23rd Annual Conference on Liquid Atomization and Spray Systems, Brno, Czech Republic, September 2010 An analytical solution to the spherically symmetric species diffusion equation: application to modelling of heating and evaporation of bi-component droplets a∗ S.S. Sazhin , P.A. Krutitskiib , A. Elwardanya , G. Castanetc , F. Lemoinec , M.R. Heikala a Sir Harry Ricardo Laboratories, School of Environment and Technology, University of Brighton, Brighton, BN2 4GJ, UK b Keldysh Institute for Applied Mathematics, Department 4, Miusskaya Sq. 4, Moscow 125047, Russia c LEMTA, Nancy-Université, 2, Avenue de la Forêt de Haye, BP 160, 54504 Vandoeuvre-lès-Nancy, France Abstract An analytical solution to the spherically symmetric species diffusion equation is derived. This solution includes an additional term, exponentially increasing with time, in the series. This term does not lead to an unphysical result due to the fact that the solution of the species mass diffusion equation is valid for mass fractions in the range from 0 to 1. Based on this solution, a numerical code for modelling bi-component droplet heating and evaporation is developed and applied to the analysis of the observed average temperature of droplets in a monodisperse spray. The code takes into account all key processes, which take place during droplet heating and evaporation, including the distribution of temperature, the diffusion of liquid species inside the droplet and the effects of the non-unity activity coefficient (ideal and non-ideal models). The predicted time evolution of the average temperatures is reasonably close to the measured ones. The predicted non-zero gradient of the mass fraction of species inside the droplet shows the limitations of widely-used models, based on the assumptions of either zero or infinitely large species diffusivities inside droplets. Introduction In [1] a model for the heating and evaporation of monocomponent droplets, based on the analytical solution of the heat conduction equation inside droplets, was suggested. The practical application of this model, however, is limited by the fact that most real-life fuels are multi-component. The models of multi-component droplet heating and evaporation could be subdivided into two main groups: those based on the analysis of individual components [2], applicable in the case when a small number of components needs to be taken into account, and those based on the probabilistic analysis of a large number of components (see review [3]). In the second family of models a number of additional simplifying assumptions were used, including the assumption that species inside droplets mix infinitely quickly. The focus of this paper is on the extension of the model, developed in [1] for monocomponent droplets, to the case of multi-component droplets. Only the first group of multi-component droplet heating and evaporation models is considered, which allows us to perform the deterministic analysis of individual species. Moreover, at this stage only the simplest case of bi-component droplets is considered, as in [2], which allows us to gain a better understanding of the underlying physics of the processes involved. The new model is based on the analytical solution of the spherically symmetric species diffusion equation, which has not yet been reported to the best of our knowledge. It is much simpler than the previously suggested models (e.g. [2]) which makes it potentially attractive for implementation into computational fluid dynamics (CFD) codes. Basic Equations and Approximations Assuming that the processes inside droplets are spherically symmetric (no convection), equations for mass fractions of liquid species Yli ≡ Yli (t, R) inside them can be presented in the following form ∂Yli = Dl ∂t  ∂ 2 Yli 2 ∂Yli + 2 ∂R R ∂R  , (1) where i = 1, 2, Dl is the liquid mass diffusivity. Equations (1) are solved with the following boundary conditions α(ǫi − Ylis ) = −Dl ∂Yli ∂R , (2) R=Rd −0 ∗ Corresponding author: S.Sazhin@brighton.ac.uk 1 ILASS – Europe 2010 An analytical solution to the spherically symmetric species diffusion equation and the initial conditions Yli (t = 0) = Yli0 (R), where Ylis = Ylis (t) are liquid components’ mass fractions at the droplet’s surface, α = |ṁd |/(4πρl Rd2 ), ṁd is the droplet evaporation rate. We are interested only in a solution which is continuously differentiable twice in the whole domain. This implies that Yli should be bounded for 0 ≤ R < Rd . Moreover, the physical meaning of Yli , as the mass fraction, implies that 0 ≤ Yli ≤ 1. Although Equations (1) with boundary conditions (1) look similar to the equation considered in [1], the solution suggested in [1] cannot be used for the above problem. The new analytical solution to this equation, subject to the above boundary and initial conditions, is presented in the next section. This solution was incorporated into a zero dimensional code used for modelling the heating and evaporation of bi-component droplets, along with a number of other models, taking into account the distribution of temperature inside the droplet and the effects of the non-unity activity coefficient (ideal and non-ideal models). The effects of recirculation in the moving droplets on heat and mass diffusion are taken into account using the effective thermal conductivity and the effective diffusivity models [3]. Analytical Solution In this section, the solution of Equation (1) (Yli (t, R)) for t ≥ 0 and 0 ≤ R < Rd is presented. Remembering the physical background of the problem, we will look for a solution which is continuously differentiable twice in the whole domain. Let us rewrite the boundary condition (2) in the form  ∂Yli α Yli − ∂R Dl  =− R=Rd αǫi (t) . Dl (3) The initial condition is Yli (t = 0) = Yli0 (R). We look for the solution to Equation (1) in the form Yli (t, R) = y(t, R) + ǫ(t), (4) where the subscript i at ǫ is omitted to simplify the notation. Having substituted (4) into Equation (1) and the boundary condition (3) we can rewrite this equation and the corresponding boundary and initial conditions in the form ∂y = Dl ∂t   ∂2y 2 ∂y + ∂R2 R ∂R ∂y α y − ∂R Dl   − dǫ(t) , dt (5) = 0, (6) R=Rd y|t=0 = Yli0 (R) − ǫ(0) ≡ Yli0 (R) − ǫ0 . (7) Introduction of the new variable u(t, R) = y(t, R)R allows us to rewrite Equation (5) and the corresponding boundary and initial conditions in the form ∂2u dǫ(t) ∂u = Dl −R , 2 ∂t ∂R dt u|R=0 =  h0 ∂u u + ∂R Rd  (8) = 0, (9) R=Rd u|t=0 = R (Yli0 (R) − ǫ0 ) , (10)   αRd where h0 = − 1 + . Note that the change of the variable from y to u leads to the need for a second Dl boundary condition at R = 0. Our assumption that the solution is continuously differentiable twice implies that y 2 ILASS – Europe 2010 An analytical solution to the spherically symmetric species diffusion equation is finite everywhere in the domain 0 ≤ R < Rd . Hence, the boundary condition (9) at R = 0. The solution to the problem (8) - (10) for h0 > −1 was earlier reported in [1]. Here the focus will be on the case h0 < −1, which is directly relevant to the problem of diffusion of species inside droplets. We look for the solution to Equation (8) in the form u= ∞ X Θn (t)vn (R), (11) n=0 where vn (R) is the full set of non-trivial solutions to the equation ∂2v + pv = 0, ∂R2 (12) subject to the boundary conditions v|R=0 =  h0 ∂v v + ∂R Rd  = 0. (13) R=Rd Equation (12) with boundary conditions (13) is the well known Sturm-Liouville problem. Our first task is to find eigenvalues p for this problem. The cases p = 0, p < 0 and p > 0 will be considered separately. For p = 0 the general solution to Equation (12) can be presented as v = A + BR. The condition v|R=0 = 0 implies that A = 0. The boundary condition at R = Rd leads to the following equation B(1 + h0 ) = 0. Since h0 ∈ (−∞, −1), the latter equation is satisfied only when B = 0. This leads to the trivial solution v = 0 which is disregarded in our analysis. Hence, Equation (12) has no non-trivial solutions for p = 0. Assuming that p = −λ2 < 0 we write the general solution to (12) as  R v(R) = A cosh λ Rd   R + B sinh λ Rd  , (14) where A and B are arbitrary constants. The boundary condition at R = 0 (see (13)) implies that A = 0. The boundary condition at R = Rd leads to the following equation B (λ cosh λ + h0 sinh λ) = 0. Rd (15) B in this equation is not equal to zero as we do not consider the trivial solution v = 0. Hence, Equation (15) can be re-written as tanh λ = − λ . h0 (16) It is easy to show that for h0 ∈ (−∞, −1) Equation (16) has three solutions λ = 0 and λ ± λ0 , where λ0 ∈ (0, +∞), for h0 < −1 and it has no solutions for h0 > −1. For h0 = −1 it has one solution λ = 0. The solution λ = 0 leads to the trivial solution v = 0, which is disregarded in our analysis. The solutions λ = ±λ0 lead to Solutions (14) (eigen functions) which differ only by the sign of B. Since the value of the coefficient B is determined by the normalisation condition only (see below), the solution λ = −λ0 can be disregarded. Hence, we can conclude that the solution of Equation (16) gives only one eigenvalue λ = λ0 > 0 and the corresponding eigen function   R , v0 (R) = sinh λ0 Rd (17) where the normalisation leading to B = 1 has been chosen. 3 ILASS – Europe 2010 An analytical solution to the spherically symmetric species diffusion equation Assuming that p = λ2 > 0 we write the general solution to (12) as     R R + B sin λ , v(R) = A cos λ Rd Rd (18) where A and B are arbitrary constants. The boundary condition at R = 0 (see (13)) implies that A = 0. The boundary condition at R = Rd leads to the following equation B (λ cos λ + h0 sin λ) = 0. Rd (19) B in this equation is not equal to zero as we do not consider the trivial solution v = 0. Hence, tan λ = − λ . h0 (20) As in the case p < 0 we disregard the solutions to this equation corresponding to zero and negative λ. A countable set of positive solutions to this equation (positive eigenvalues) λn are arranged in ascending order: 0 < λ1 < λ2 < λ3 < ... The corresponding eigen function can be presented as   R , vn (R) = sin λn Rd (21) where the normalisation leading to B = 1 has been chosen as in the case p < 0. The orthogonality of functions vn can be proven by direct integration of their products over R from 0 to Rd . Note that R can be presented as a Fourier series with respect to functions vn R= ∞ X Qn vn (R), (22) n=0 where  2 Rd 1 (1 + h0 ) sinh λ0 ||v0 ||2 λ0 Qn =  2  Rd 1   (1 + h0 ) sin λn ||vn ||2 λn     − when n=0 when n≥1 (23) Having substituted Expressions (11) and (22) into Equation (8), we rewrite the latter equation in the form ∞ X Θ′n (t)vn (R) = Dl n=0 ∞ X Θn (t)vn′′ (R) − ǫ′ (t) ∞ X Qn vn (R), (24) n=0 n=0 dΘn d2 vn dǫ where Θ′n = ; vn′′ (R) = ; ǫ′ (t) = ≡ ǫ′ . Since the expansion in the series with respect to dt dR2 dt vn (Fourier series) is unique, Equation (24) is satisfied only when it is satisfied for each term in this expansion. Remembering that  2  2 λn λ0 ′′ ′′ v0 and vn = − vn (n ≥ 1), v0 = Rd Rd we can see that this implies that Θ′0 (t) = Dl  λ0 Rd 2 Θ0 (t) − ǫ′ Q0 , Θ′n (t) = −Dl when n ≥ 1. 4  λn Rd 2 Θn (t) − ǫ′ Qn , ILASS – Europe 2010 An analytical solution to the spherically symmetric species diffusion equation These equations need to be solved subject to the initial conditions for Θn (t) (n ≥ 0). To find these initial conditions we substitute (11) into the initial condition (10) and expand RYli0 (R) into a Fourier series with respect to vn . Remembering that the expansion with respect to vn is unique, this leads us to the following equation for Θn (0) Θn (0) = qin − ǫ(0)Qn , (25) where Qn is defined by (23),0 qin = 1 Z ||vn ||2 Rd RYli0 (R)vn (R)dR, (26) 0 n ≥ 0. The solutions to equations for Θn subject to the initial condition (25) can be presented in the form " Θ0 (t) = exp Dl  " Θn (t) = exp −Dl λ0 Rd  2 # Z t [qi0 − ǫ(0)Q0 ] − Q0 # "   2 dǫ(τ ) λ0 (t − τ ) dτ, exp Dl dτ Rd t 0 λn Rd 2 # Z t [qin − ǫ(0)Qn ] − Qn t 0 # "  2 dǫ(τ ) λn (t − τ ) dτ, exp −Dl dτ Rd (27) (28) where n ≥ 1. Having substituted (17), (21), (27) and (28), into (11), we can present the final solution to Equation (1), satisfying boundary condition (2) and the corresponding initial condition in the form: # # (" "   # "     Z t 2 2 dǫi (τ ) R λ0 λ0 1 t [qi0 − ǫi (0)Q0 ] −Q0 (t − τ ) dτ sinh λ0 exp Dl exp Dl Yli = ǫi + R Rd dτ Rd Rd 0 + ∞ X n=1 " " exp −Dl  λn Rd 2 # Z t [qin − ǫi (0)Qn ] −Qn 0 t # # "  )  2 dǫi (τ ) R λn , (t − τ ) dτ sin λn exp −Dl dτ Rd Rd (29) where Qn , qin , λ0 and λn (n ≥ 1) are defined by Equations (23), (26), (16) and (20) respectively; the subscript i at ǫ has been restored. Note that Expression (29) contains the term which exponentially increases with time. This, however, will not lead to an unphysical solution to Equation (1), since this equation is valid only for 0 ≤ Yli ≤ 1. Once Yli reaches one of its limiting values it will remain equal to this value. Results The previously obtained analytical solution of the transient heat conduction equation has been incorporated in the numerical code alongside the original analytical solution to the species diffusion equation inside droplets, described in the previous section. The code takes into account all key processes, which take place during droplet heating and evaporation, including the distribution of temperature, the diffusion of liquid species inside the droplet and the effects of the non-unity activity coefficient (ideal and non-ideal models). The predicted time evolutions of surface, average and central droplet temperatures have been compared with the results of direct measurements of droplet average temperatures for the case of various mixtures of ethanol and acetone. The ambient gas temperature and droplet velocities were measured simultaneously. This allowed us to ignore the effects of droplets on the ambient gas, and consider only the effects of gas on droplets. It has been shown that the code underpredicts the average droplet temperature at the final stages of droplet heating and evaporation for the case of acetone and acetone rich droplets. There is a general agreement between the predicted and observed average temperatures in the case of pure ethanol and 50% ethanol – 50% acetone mixture droplets. In the case of 75% ethanol – 25% acetone mixture droplets, the predicted average droplet temperature is several degrees higher compared with the observed one. It has been shown that the temperatures predicted by the numerical code and the earlier reported 5 ILASS – Europe 2010 An analytical solution to the spherically symmetric species diffusion equation vortex model are reasonably close. Also, the temperatures predicted by the ideal and non-ideal models differ by not more that several degrees. This can justify the application of the code with the activity coefficient equal to 1 for the interpretation of the time evolution of temperatures measured with similar errors. The predicted radial distributions of temperature and species mass fractions inside droplets have been shown to be consistent with the physical nature of the phenomenon. Conclusions An analytical solution to the spherically symmetric species diffusion equation is presented and discussed. This solution includes an additional term, exponentially increasing with time, in the series. This term does not lead to an unphysical result due to the fact that the solution of the species mass diffusion equation is valid for mass fractions in the range from 0 to 1. This imposes the restrictions on the range of its validity. This solution has been implemented into the numerical code alongside the previously obtained analytical solution of the transient heat conduction equation. The effects of the non-unity activity coefficient (ideal and non-ideal models) are taken into account. This code has been applied to the analysis of the observed average temperature of droplets in a monodisperse spray. The predicted time evolution of the average temperatures is reasonably close to the measured ones, especially in the case of pure ethanol and 50% ethanol – 50% acetone mixture droplets. The temperatures predicted by the ideal and non-ideal models differ by not more than several degrees. The predicted radial distributions of temperature and species mass fractions inside droplets are consistent with the physical nature of the process under consideration. The non-zero gradient of the mass fraction of ethanol inside the droplet shows the limitations of widely-used models, based on the assumptions of either zero or infinitely large species diffusivities inside droplets Acknowledgments The authors are grateful to the European Regional Development Fund Franco-British INTERREG IVA (Project C5, Reference 4005) and the Royal Society (UK) (Project JP71790) for the financial support of this project. Nomenclature D diffusion coefficient [m2 / s] h0 parameter introduced in Eq. (6) m mass [kg] p eigen values introduced in Eq. (12) qin parameter defined by Eq. (26) R distance from the droplet centre [m] Rd droplet radius [m] t time [s] u yR [m] vn eigen function y variable introduced in Eq. (4) Y mass fraction α droplet evaporation rate [m/s] ǫ evaporation rate of species Θn function introduced in Eq. (11) λn eigen values ρ density [kg·s-3 ] Subscripts d droplet i species l liquid 0 initial or zero-order References [1] Sazhin, S.S., Krutitskii, P.A., Abdelghaffar, W.A., Mikhalovsky, S.V., Meikle, S.T. and Heikal, M.R. Int J Heat Mass Transfer 47:3327-3340 (2004). [2] Maqua, C., Castanet, G. and Lemoine, F. Fuel 87:2932-2942 (2008). [3] Sazhin, S.S. Prog. in Energy Combustion Sci. 32:162-214 (2006). 6