A one-phase Stefan problem with size-dependent thermal conductivity
aa r X i v : . [ c ond - m a t . o t h e r] M a r A one-phase Stefan problem with size-dependentthermal conductivity
F. FontDepartment of Physics, Universitat Polit`ecnica de Catalunya, Av.Dr. Mara˜non 44-50, 08028 Barcelona, SpainMarch 21, 2018
Abstract
In this paper a one-phase Stefan problem with size-dependent thermalconductivity is analysed. Approximate solutions to the problem are foundvia perturbation and numerical methods, and compared to the Neumannsolution for the equivalent Stefan problem with constant conductivity. Wefind that the size-dependant thermal conductivity, relevant in the contextof solidification at the nanoscale, slows down the solidification process. Asmall time asymptotic analysis reveals that the position of the solidificationfront in this regime behaves linearly with time, in contrast to the Neumannsolution characterized by a square root of time proportionality. This hasan important physical consequence, namely the speed of the front predictedby size-dependant conductivity model is finite while the Neumann solutionpredicts an infinite and, thus, unrealistic speed as t → The Stefan problem, describing the phase change of a material, is one of the mostpopular problems in the moving boundary problem literature. Typically, it re-quires solving heat equations for the temperature in the two phases ( e.g. solid andliquid), while the position of the front separating them, the moving boundary, isdetermined from an energy balance referred to as the Stefan condition. The Stefanproblem has been studied in great detail since Lam´e and Clapeyron formulated it inthe 19th century [1]. There are several reference books that the reader may consult1or a comprehensive background on the classical problem [2, 3, 4, 5, 6, 7]. The prob-lem only admits an exact solution in the Cartesian one-dimensional case, referredto as the Neumann solution (see, for instance, [6]). Solutions in other geometries orhigher dimensions are usually found numerically or via asymptotic/perturbationtechniques (see, for example, [8] for an investigation of the classical two-phase Ste-fan problem in a sphere or an analysis of the solidification of a liquid half-space[9]). A complete literature review of the classical Stefan problem is not the purposeof this study and the reader is referred to the new edition of the book by Gupta[7] for an up to date bibliography on the problem.The classical problem has been modified in different ways to introduce newphysical phenomena such as supercooling or curvature dependent phase changetemperature [10, 11, 12, 13, 14, 15, 16, 17]. Modifications are usually linked witha thermophysical property of the material which changes with, for example, thegeometry of the system ( e.g. curvature-induced melting point depression [12, 13,14, 15, 16]), the speed of the moving boundary ( e.g. supercooling [10, 11, 18]), orthe temperature itself [17]. Motivated by recent experimental studies on Siliconnanofilms and nanowires showing that the thermal conductivity decreases as thesize of the physical system decreases [19, 20], in this work we will consider theeffect of size-dependent thermal conductivity on the solidification process of a one-dimensional slab. In Jou et al. [21] an analytical expression for the dependence ofthe thermal conductivity of a solid on the size of the physical system is derivedfrom the Extended Irreversible Thermodynamics theory [22]. Assuming that allthe phonon mean-free paths and relaxation times are equal, their expression forthe thermal conductivity takes the form k = 2 k L l r l L − ! , (1)where k represents the bulk thermal conductivity of the solid, L the size of thesolid, l the phonon mean free path [21]. In a follow-up paper, expression (1) istested against experimental data showing good agreement [23]. The main goal ofthis paper is to study the effect of the size-dependent thermal conductivity on asolidification process by introducing (1) in the formulation of the one-phase Stefanproblem.The paper is organized as follows. In the next section we formulate the Stefanproblem with size-dependent thermal conductivity and discuss how the Neumannsolution can be retrieved. In section 3 we provide a perturbation solution basedon large Stefan number. In section 4 we present the numerical strategy to solvethe problem and analyse the small time limit, which is needed to initialize thenumerical scheme. In section 5 we discuss our results and in section 6 we draw ourconclusions. 2 (kg/m ) c (J /kg · K) k (W/m · K) T f K ∆ H (J/kg) l (m)2320 770 120 1687 1926 × × − Table 1: Thermophysical properties of Silicon.
Consider a liquid initially at its equilibrium freezing temperature, T f , occupyingthe space x ≥
0. Suddenly, the temperature is lowered to T c < T f on the edge x = 0 and the liquid starts to solidify. The newly created solid phase will start togrow occupying the space 0 < x < s ( t ), where s ( t ) represents the position of thesolidification front as well as the size of the solid phase. The temperature of thesolid phase, T ( x, t ), is described by ρc ∂T∂t = k ( s ( t )) ∂ T∂z on 0 < z < s ( t ) , (2)where ρ is the density, c the specific heat, and k ( s ( t )) the thermal conductivitywhich depends on the size of the solid phase and is obtained by setting L = s ( t )in (1). That is, k ( s ( t )) = 2 k s ( t ) l s l s ( t ) − ! . (3)Note, in the current study, we consider the one-phase approximation, i.e. we as-sume the liquid to be at the equilibrium freezing temperature, so no equation isneeded for the temperature of the liquid phase. The temperature of the solid issubject to the boundary conditions T (0 , t ) = T c , T ( s ( t ) , t ) = T f . (4)Finally, at the solidification front we have the Stefan condition ρ ∆ H dsdt = k ( s ( t )) ∂T∂z (cid:12)(cid:12)(cid:12)(cid:12) z = s ( t ) , (5)where ∆ H is the latent heat of fusion. To obtain physically meaningful results toour problem we use the parameter values for Silicon in Table 1. Introducing the nondimensional variables T ∗ = T − T f T f − T c , x ∗ = xl , t ∗ = kl ρc t , s ∗ = sl ,
3n (2)-(5) and removing the star notation, we obtain the dimensionless model ∂T∂t = 2 s (cid:16) √ s + 1 − s (cid:17) ∂ T∂z on 0 < z < s ( t ) , (6a) T (0 , t ) = − , (6b) T ( s ( t ) , t ) = 0 , (6c) β dsdt = 2 s (cid:16) √ s + 1 − s (cid:17) ∂T∂z (cid:12)(cid:12)(cid:12)(cid:12) z = s , (6d) s (0) = 0 , (6e)where β = ∆ H/ρc ( T f − T c ) is the Stefan number (ratio of latent heat to sensibleheat). The term 2 s (cid:0) √ s + 1 − s (cid:1) in (6a)-(6e) represents the thermal conductivity(3) in dimensionless form, which takes the values 0 for s = 0 and 1 for s → ∞ . The classical one-phase Stefan problem is retrieved by setting 2 s ( √ s + 1 − s ) to1 in (6a) and (6e), i.e. by considering the bulk value of the thermal conductivity.In this case, the problem (6a)-(6e) is reduced to an initial value problem via thesimilarity transformation η = x/ √ t . The resulting system has an exact solution,whose temperature and position of the moving boundary are given by T ( x, t ) = − (cid:16) x √ t (cid:17) erf( λ ) , s ( t ) = 2 λ √ t . (7)where erf ( z ) = √ π R z e − y dy is the error function. The constant λ is found bysolving the transcendental equation β √ πλ erf( λ ) e λ = 1 . (8)System (7)-(8) is referred to as the Neumann solution [3]. The complexity introduced by the size-dependent thermal conductivity in the gov-erning equation prevents an exact similarity solution as in the classical Stefanproblem. To make analytical progress a perturbation solution in the limit of largeStefan number is developed. Physically, a large Stefan number corresponds toslow solidification (as can be seen from β ∝ / ( T f − T c ), so large β implies smalltemperature drop). This is consistent with the physical properties of the referencematerial in Table 1, for which β is always larger than 1 (even T c = 0 gives β = 1 . t = β ˆ t and defining the small parameter δ = β − , theproblem becomes δ ∂T∂ ˆ t = 2 s (cid:16) √ s + 1 − s (cid:17) ∂ T∂z on 0 < z < s ( t ) , (9a) T (0 , ˆ t ) = − , (9b) T ( s (ˆ t ) , ˆ t ) = 0 , (9c) dsd ˆ t = 2 s (cid:16) √ s + 1 − s (cid:17) ∂T∂z (cid:12)(cid:12)(cid:12)(cid:12) z = s , (9d) s (0) = 0 , (9e)suggesting an expansion in the form T ( x, ˆ t ) = T + δ T + O ( δ ). We obtain theleading and first order problems O (1) : 0 = ∂ T ∂z , T (1 , ˆ t ) = − , T ( s (ˆ t ) , ˆ t ) = 0 , O ( δ ) : ∂T ∂ ˆ t = 2 s (cid:16) √ s + 1 − s (cid:17) ∂ T ∂z , T (1 , ˆ t ) = 0 , T ( s (ˆ t ) , ˆ t ) = 0 , with solutions T = − xs , T = s ˆ t s (cid:0) √ s + 1 − s (cid:1) x ( s − x ) . These lead to the temperature profile T ( x, ˆ t ) = − xs + 1 β s ˆ t s (cid:0) √ s + 1 − s (cid:1) x ( s − x ) + O (cid:0) δ (cid:1) . (10)Inserting T ≈ T + δ T into (9d) we obtain the following expression for the speedof the moving boundary s ˆ t = 6 β (1 + 3 β ) (cid:16) √ s + 1 − s (cid:17) , (11)which can be readily integrated to give s + s √ s + 1 + arcsin( s ) = 12 β (1 + 3 β ) ˆ t , (12)where the initial condition (9e) was applied.5 Boundary immobilisation and numerical solu-tion
A typical problem when seeking numerical solutions to moving boundary problemsis how to deal with the discretization of a domain whose size changes in time.A way to overcome this difficulty consists in defining a new space variable thattransforms the variable domain into a fixed domain. This method is referred toas the boundary immobilization method. To immobilise the boundary s ( t ) wemap the space variable x to the unit domain via the Landau-type transformation ξ = x/s ( t ). Problem (6a)-(6e) then becomes s ∂T∂t = s t ξ ∂T∂ξ + 2 (cid:16) √ s + 1 − s (cid:17) ∂ T∂ξ on 0 < ξ < T (0 , t ) = − , (13b) T (1 , t ) = 0 , (13c) βs t = 2 (cid:16) √ s + 1 − s (cid:17) ∂T∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =1 , (13d) s (0) = 0 . (13e)A Backward Euler semi-implicit finite difference scheme is used on (13a), dis-cretising implicitly for T ( ξ, t ) and explicitly for s ( t ) and s t ( t ), and using secondorder central differences in space [24]. The semi-implicit scheme allows equation(13a) (containing time dependent coefficients) to be formulated, after discretising,as a matrix linear system which can be solved by inverting the matrix of the sys-tem at each time step. The position of the solidification front is found via (13d),using a backward difference for s t and a one-sided second order difference for thepartial derivative. A recurring complication associated with the numerical solution of Stefan problemsis how to initiate the computation in a region which initially has zero thickness.A typical approach to tackle this issue is to find a small time asymptotic approx-imation for T and s and use this as the initial condition in the numerical scheme[25].To study the small time limit we assume s ( t ) to be a power of time thatmatches the initial condition s (0) = 0. Hence, we set s = λ t p where λ and p are constants. We rescale time as t = ετ , where ǫ ≪
1, to obtain the expressions s = λ ε p τ p and s t = pλ ε p − τ p − for the position and velocity of the solidification6ront, respectively. Substituting s and s t into (13d) leads to β pλ ε p − τ p − = (cid:16) √ ε p τ p + 1 − λ ε p τ p (cid:17) ∂T∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =1 , which, after multiplying the right hand side by ( √ ε p τ p + 1+ λ ε p τ p )( √ ε p τ p + 1+ λ ε p τ p ) − , rearranging and noting that √ ε p τ p + 1 ≈ ε p τ p ≪
1, gives β pλ ε p − τ p − (1 + λ ε p τ p ) ≈ ∂T∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =1 . (14)Now we need to choose p . We realize that, in order for the front to move, theleft hand side of (14) must be O (1), thereby balancing the right hand side. Thebalance is only satisfied for p = 1, so we obtain s = λ t , (15)for t →
0. Substituting (15) into (13a) and taking the limit t → ξ λ ∂T∂ξ + ∂ T∂ξ , with solution T = − (cid:16) √ λ ξ (cid:17) erf (cid:16) √ λ (cid:17) . (16)Finally, the constant λ is found by substituting (15) and (16) into (13d) andtaking the limit t →
0, leading to the transcendental equation β √ π p λ erf (cid:18) p λ (cid:19) e λ / = 1 . (17)From a physical standpoint, expression (15) states that the speed of the solidifica-tion front, s t , is constant as t →
0, in contrast to the Neumann solution (8) whichgives s t → ∞ as t → ξ = x/s ( t ) into (16) and using (15) we obtainan initial size and temperature profile for the solid phase at some small time t >
0, which we use to initialise our numerical scheme. We note that by defining √ λ / λ the transcendental equation (17) becomes equal to that in (8).7 Results and discussion
In Figure 1 we present the perturbation and numerical solution to the Stefanproblem with size-dependent thermal conductivity for two values of the Stefannumber, β = 20 and β = 5, along with the corresponding Neumann solution tothe classical problem. The perturbation and numerical solutions to the problemwith variable thermal conductivity show excellent agreement, even for a relativelysmall Stefan number β = 5, thereby validating the accuracy of the numericalsolution.Figure 1: Temperature profiles at different times and position of the solidificationfront for the Stefan problem with size-dependent thermal conductivity and thestandard Stefan problem. The circles, the solid line and the dash-dotted line rep-resent the perturbation solution, the numerical solution and the Neumann solution,respectively. Panels (a)-(b) correspond to β = 20 and (c)-(d) to β = 5.By comparing the solutions of the variable conductivity problem with the Neu-mann solution we can assert that the overall effect of the size-dependent ther-mal conductivity on the solidification process is to delay the propagation of thetemperature and the solidification front. Indeed, the delay is caused by a smallvalue of the thermal conductivity, represented by the term 2 s ( √ s + 1 − s ), in the8arly stages of the solidification process when s ∼ O (1). However, the value of2 s ( √ s + 1 − s ) quickly converges to the maximum value of 1 (for instance, s ≈ s ( √ s + 1 − s ) ≈ .
95) and the solution to the Stefan problem withsize-dependent thermal conductivity converges to the Neumann solution. This isfurther illustrated in Figure 2.In Figure 2 we show the position ((a) panel) and velocity ((b) panel) of thesolidification front spanning several orders of magnitude in time for the case β =5. The solution of the Stefan problem with size-dependent thermal conductivityrepresented by the solid line (for clarity we only show the numerical solution),clearly shows different qualitative behaviour for small and large times. For smalltimes, the position of the front behaves linearly with time (according to eq. (15))which is represented by the dashed line which has a slope 1. For large times, theposition is well approximated by the Neumann solution, which is proportional to √ t (according to eq. (8)), and so gives a slope of 0.5. The fact that s ∝ t as t → s t ∝ / √ t , leading to an infinite speed for the solidification front at the beginningof the process.Figure 2: Position and velocity of the solidification front as a function of timespanning several orders of magnitude for the case β = 5. The solid line, dashedline and the dash-dotted line represent the numerical, the small time asymptoticand the Neumann solutions, respectively. In this paper we analysed a one-phase Stefan problem subject to size-dependentthermal conductivity. The results indicate that the size-dependent thermal con-9uctivity induces a delay in the propagation of the solidification front, which iscaused by the reduced value of the conductivity when the solid phase is small. Asmall time limit analysis revealed that the speed of the solidification front is con-stant as t →
0, whereas the Neumann solution to the classical problem predictsan unrealistic infinite velocity as time goes to zero. This result is very relevant asit gives a physically realistic mathematical description of how solidification phe-nomena initiates. When the size of the solid phase is approximately order 1, thequalitative behaviour of the solution switches to the standard √ t proportionalityand the solution to the problem tends to the classical Neumann solution. Acknowledgements
This work was partially funded by the Mathematics for Industry Network (MI-NET, COST Action TD1409). The author would like to thank Vincent Creganfor a critical reading of the manuscript.
References [1] G. Lame, B. Clapeyron., M´emoire sur la solidification par refroidissementd’un globe solide, Annales Chimie Physique 47 (1831) 250–256.[2] L. Rubinstein, The Stefan Problem, 1st Edition, Translations of MathematicalMonographs Vol. 27, (Ame. Math. Soc.), 1971.[3] J. Crank, Free and Moving Boundary Problems, 1st Edition, Oxford Univer-sity Press, 1996.[4] J. M. Hill, One-Dimensional Stefan Problems: An Introduction, 1st Edition,Longman Scientific & Technical, 1987.[5] V. Alexiades, A. Solomon, Mathematical Modelling of Freezing and MeltingProcesses, 1st Edition, Hemisphere Publishing Corporation, 1993.[6] H. Carslaw, J. Jaeger, Conduction of Heat in Solids, 1st Edition, OxfordUniversity Press, 1973.[7] S. Gupta, The Classical Stefan Problem: Basic Concepts, Modelling andAnalysis with Quasi-Analytical Solutions and Methods, 2nd Edition, Else-vier, 2017. 108] S. W. McCue, B. Wu, J. M. Hill, Classical two-phase Stefan problem forspheres, Proceedings of the Royal Society of London A: Mathematical, Phys-ical and Engineering Sciences 464 (2096) (2008) 2055–2076.[9] A. M. Wallman, J. R. King, D. S. Riley, Asymptotic and numerical solutionsfor the two-dimensional solidification of a liquid half-space, Proceedings of theRoyal Society of London A: Mathematical, Physical and Engineering Sciences453 (1962) (1997) 1397–1410.[10] X. Weiqing, The Stefan problem with a kinetic condition at the free boundary,SIAM Journal on Mathematical Analysis 21 (2) (1990) 362–373.[11] J. D. Evans, J. R. King, Asymptotic results for the Stefan problem with kineticundercooling, The Quarterly Journal of Mechanics and Applied Mathematics53 (3) (2000) 449–473.[12] B. Wu, S. W. McCue, P. Tillman, J. M. Hill, Single phase limit for meltingnanoparticles, Applied Mathematical Modelling 33 (5) (2009) 2349–2367.[13] F. Font, T. Myers, Spherically symmetric nanoparticle melting with a variablephase change temperature, Journal of Nanoparticle Research 15 (12) (2013)2086.[14] F. Font, T. G. Myers, S. L. Mitchell, A mathematical model for nanoparticlemelting with density change, Microfluidics and Nanofluidics 18 (2) (2015)233–243.[15] J. M. Back, S. W. McCue, T. J. Moroney, Including nonequilibrium inter-face kinetics in a continuum model for melting nanoscaled particles, ScientificReports 4 (2014) 7066.[16] F. I. Dragomirescu, K. Eisenschmidt, C. Rohde, B. Weigand, Perturbation so-lutions for the finite radially symmetric Stefan problem, International Journalof Thermal Sciences 104 (2016) 386–395.[17] J. E. Sunderland, S. H. Cho, Phase change problems with temperature-dependent thermal conductivity, Journal of Heat Transfer 96 (2) (1974) 214–217.[18] F. Font, S. Mitchell, T. Myers, One-dimensional solidification of supercooledmelts, International Journal of Heat and Mass Transfer 62 (2013) 411–421.[19] D. Li, Y. Wu, P. Kim, L. Shi, P. Yang, A. Majumdar, Thermal conductivityof individual silicon nanowires, Applied Physics Letters 83 (14) (2003) 2934–2936. 1120] W. Liu, M. Asheghi, Phononboundary scattering in ultrathin single-crystalsilicon layers, Applied Physics Letters 84 (19) (2004) 3819–3821.[21] D. Jou, J. Casas-Vazquez, G. Lebon, M. Grmela, A phenomenological scalingapproach for heat transport in nano-systems, Applied Mathematics Letters18 (8) (2005) 963–967.[22] D. Jou, G. Lebon, J. Casas-V´azquez, Extended Irreversible Thermodynamics,4th Edition, Springer Netherlands, 2010.[23] F. X. Alvarez, D. Jou, Memory and nonlocal effects in heat transport: Fromdiffusive to ballistic regimes, Applied Physics Letters 90 (8) (2007) 083109.[24] J. Caldwell, Y. Y. Kwan, Numerical methods for one-dimensional stefan problems,Communications in Numerical Methods in Engineering 20 (7) (2004) 535–545. doi:10.1002/cnm.691 .URL http://dx.doi.org/10.1002/cnm.691http://dx.doi.org/10.1002/cnm.691