Non-local effects and size-dependent properties in Stefan problems with Newton cooling

TL;DR

This paper investigates how size-dependent thermal conductivity and Newton cooling influence the growth of a solid in Stefan problems, revealing that non-local effects diminish with smaller Biot numbers through numerical and asymptotic analysis.

Contribution

It introduces a modified Stefan problem model incorporating size-dependent conductivity and analyzes the impact of small Biot numbers on solidification.

Findings

Non-local effects are less significant at small Biot numbers.

Size-dependent properties influence solidification dynamics.

Asymptotic methods effectively analyze the problem.

Abstract

We model the growth of a one-dimensional solid by considering a modified Fourier law with a size-dependent effective thermal conductivity and a Newton cooling condition at the interface between the solid and the cold environment. In the limit of a large Biot number, this condition becomes the commonly used fixed-temperature condition. It is shown that in practise the size of this non-dimensional number is very small. We study the effect of a small Biot number on the solidification process with numerical and asymptotic solution methods. The study indicates that non-local effects become less important as the Biot number decreases.

Figures

Tables

Table 1. Table 1 : Thermophysical parameters for Silicon at 1000 K [ 26 , 27 ] .

$k^{\ast }$ [W/m$\cdot$K]	$c^{\ast }$ [J/kg$\cdot$K]	${\rho }^{\ast }$ [kg/m${}^{\text{3}}$]	$L_m^{\ast }$ [kJ/kg]	$T_{\text{f}}^{\ast }$ [K]	${\ell }^{\ast }$ [$\mu$m]
43.67	864.89	2296	1787	1687	12.84

Equations

$$c^{\ast }{\rho }^{\ast }\frac{\partial T^{\ast }}{\partial t^{\ast }}=-\frac{\partial q^{\ast }}{\partial x^{\ast }}.$$

$$q^{\ast }=-k^{\ast }\frac{\partial T^{\ast }}{\partial x^{\ast }}.$$

$$T^{\ast }(s^{\ast },t^{\ast })=T_{\text{f}}^{\ast }.$$

$$q^{\ast }(0,t^{\ast })=h^{\ast }(T_{\text{e}}^{\ast }-T^{\ast }(0,t^{\ast })),$$

$${\rho }^{\ast }{L}_m^{\ast }\frac{\mathrm{d}{s}^{\ast }}{\mathrm{d}{t}^{\ast }}=-q^{\ast }(s^{\ast },t^{\ast }),$$

$$s^{\ast }(0)=0.$$

$$k_{\text{eff}}^{\ast }(L^{\ast })=2k^{\ast }{\left(\frac{L^{\ast }}{{\ell }^{\ast }}\right)}^2\left(\sqrt{1+{\left(\frac{{\ell }^{\ast }}{L^{\ast }}\right)}^2}-1\right).$$

$$q^{\ast }=-k_{\text{eff}}^{\ast }(s^{\ast })\frac{\partial T^{\ast }}{\partial x^{\ast }},$$

$$c^{\ast }{\rho }^{\ast }\frac{\partial T^{\ast }}{\partial t^{\ast }}=k_{\text{eff}}^{\ast }(s^{\ast })\frac{{\partial }^2{T}^{\ast }}{\partial x^{\ast 2}},$$

$$-k^{*}_{\text{eff}}(s^{*})\frac{\partial T^{*}}{\partial x^{*}}\bigg{|}_{x^{*}=0}=h^{*}(T^{*}_{\text{e}}-T^{*}(0,t^{*})).$$

$$\rho^{*}L^{*}_{m}\frac{\mathrm{d}s^{*}}{\mathrm{d}t^{*}}=k^{*}_{\text{eff}}(s^{*})\frac{\partial T^{*}}{\partial x^{*}}\bigg{|}_{x^{*}=s^{*}}.$$

$$T_t=f(s)T_{xx},$$

$$f(s)T_x=\text{Bi}(1+T),$$

$$T=0,$$

$$\beta s_t=f(s)T_x,$$

$$s=0,$$

$$s{u}_t=\xi s_t{u}_{\xi }+F(s)u_{\xi \xi },$$

$$F(s)u_{\xi }=\text{Bi}(1+u),$$

$$u=0,$$

$$\beta s_t=F(s)u_{\xi },$$

$$\frac{1}{2}\beta \lambda p\left(t^{p-1}+\lambda t^{2p-1}\right)=u_{\xi },$$

$$s\approx \lambda t,s_t\approx \lambda ,F\approx 2(1-\lambda t),$$

$$\lambda \xi u_{\xi }+2u_{\xi \xi }=0,$$

$$2u_{\xi }(0)=\text{Bi}(1+u(0)),$$

$$u(1)=0.$$

$$u(\xi )=\text{Bi}\frac{\text{erf}\left(\sqrt{\lambda }\xi /2\right)-\text{erf}\left(\sqrt{\lambda }/2\right)}{2\sqrt{\lambda /\pi }+\text{Bi}\text{erf}\left(\sqrt{\lambda }/2\right)},$$

$$\lambda +\frac{\sqrt{\pi }}{2}\text{Bi}\sqrt{\lambda }\text{erf}\left(\frac{\sqrt{\lambda }}{2}\right)=\frac{\text{Bi}}{\beta }{e}^{-\lambda /4}.$$

$$\text{erf}\left(\frac{\sqrt{\lambda }}{2}\right)\approx \frac{\sqrt{\lambda }}{\sqrt{\pi }}-\frac{\lambda \sqrt{\lambda }}{12\sqrt{\pi }},\exp \left(\frac{\lambda }{4}\right)\approx 1+\frac{\lambda }{4},$$

$$\left(1+\frac{\text{Bi}}{2}\right)\lambda +\frac{1}{4}\left(1+\frac{\text{Bi}}{3}\right){\lambda }^2+O({\lambda }^3)=\frac{\text{Bi}}{\beta }.$$

$$\lambda =\delta {\lambda }_0+{\delta }^2{\lambda }_1+O({\delta }^3).$$

$$\left(1+\frac{\text{Bi}}{2}\right){\lambda }_0=1,\left(1+\frac{\text{Bi}}{2}\right){\lambda }_1+\frac{1}{4}\left(1+\frac{\text{Bi}}{3}\right){\lambda }_0^2=0,$$

$$\lambda =\frac{2}{2+\text{Bi}}\frac{\text{Bi}}{\beta }-\frac{2(\text{Bi}+3)}{3(2+\text{Bi}{)}^3}\frac{{\text{Bi}}^2}{{\beta }^2}+O\left(\frac{{\text{Bi}}^3}{{\beta }^3}\right).$$

$$u(\xi )\approx \frac{\text{Bi}(\xi -1)}{2+\text{Bi}},$$

$$\frac{\partial u}{\partial t}\approx \frac{u_i^{n+1}-u_i^n}{\Delta t},\frac{\partial u}{\partial \xi }\approx \frac{u_{i+1}^{n+1}-u_{i-1}^{n+1}}{2\Delta \xi },\frac{{\partial }^{2}u}{\partial {\xi }^2}\approx \frac{u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1}}{(\Delta \xi {)}^2},$$

$$A_i^n{u}_{i-1}^{n+1}+B_i^n{u}_i^{n+1}+C_i^n{u}_{i+1}^{n+1}=s^n{u}_i^n,1\le i\le N-1,$$

Full text

1 Introduction

It is widely accepted that heat conduction at small length scales differs from the classical description [1, 2, 3]. There exist a wide range of theoretical models which extend Fourier’s law to account for non-local effects which become dominant on length scales comparable to the phonon mean free path (MFP). The models may be classified into micro-, meso- or macroscopic models depending on whether the aim is focused on describing the behaviour of the individual heat carriers, the evolution of their distribution or related macroscopic quantities, such as the temperature or the heat flux. Examples for micro- and mesoscopic models are molecular dynamics, Monte Carlo simulations [4], the Boltzmann transport equation (BTE) [5] or the equation of phonon radiative transfer [6]. The equations involved in these approaches often require excessive computational effort and it is therefore more sensible to aim for a macroscopic equation which captures the physics of the problem. Popular macroscopic models are the Maxwell–Cattaneo law [7], the thermomass model [8, 9] or the Guyer–Krumhansl (GK) equation [10, 11] and the framework of phonon hydrodynamics developed from it [9, 12, 13, 14].

Based on the BTE and on the framework of extended irreversible thermodynamics [12], Alvarez and Jou [15] derived a model for heat flow with a size-dependent effective thermal conductivity (ETC). This was able to capture non-local effects and showed good agreement with experimental data. Hennessy et al. [16] show that such an ETC can also be derived when the GK equation is included into the formulation of the Stefan problem. Recently, Font [17] proposed a simpler non-classical formulation of the Stefan problem which includes an effective Fourier law where a size-dependent ETC replaces the bulk thermal conductivity. He showed that, even when the temperature at the boundary is instantly set to the temperature of the cold environment, the problem of an initially infinite solidification rate is avoided by incorporating a size-dependent ETC.

The popular fixed-temperature condition is a specific limit of a more general boundary condition known as the Newton cooling condition. In this paper we aim to generalise the model proposed by Font by considering cooling conditions at the interface between the material and its environment. For large Biot (or Nusselt) numbers, this condition converges to the fixed-temperature condition. Formulations of the Stefan problem where Newton cooling conditions are applied, can be found in the framework of nanoparticle melting [18, 19, 20] or nanowire melting [21]. Hennessy et al. [20] recently proposed an extended form of cooling condition accounting for memory effects. However, the case of small Biot number is not considered in these studies. In Sec. 2.4 we explain the importance of this limit.

2 Mathematical model

We consider a liquid bath, initially at the phase change temperature $T^{*}_{\text{f}}$, that occupies the space $x^{*}\geq 0$. The ∗ notation refers to dimensional quantities. Due to a low external temperature $T^{*}_{\text{e}}$ applied at $x^{*}=0$, the liquid starts to undergo a solidification process and a solid starts to grow into the liquid occupying the space $0\leq x^{*}\leq s^{*}(t^{*})$; see Fig. 1. For simplicity, the specific heat $c^{*}$, the density $\rho^{*}$ and the bulk thermal conductivity $k^{*}$ are assumed to be constant, and no supercooling will be considered.

2.1 Classical formulation

The process is driven by the heat flow in the solid phase. The liquid remains at the phase change temperature $T^{*}_{\text{f}}$ for all time. The temperature $T^{*}$ and heat flux $q^{*}$ in the solid are related via conservation of energy, which in one spatial dimension takes the form

[TABLE]

In the classical formulation, the heat flux is determined by the gradient of the temperature,

[TABLE]

At the solid-liquid interface, the temperature of the solid is equal to the phase change temperature,

[TABLE]

The heat exchange at the interface between the material and the environment is governed by a Newton cooling condition

[TABLE]

where $h^{*}$ is the heat transfer coefficient. Contrary to the fixed temperature condition $T^{*}=T^{*}_{\text{e}}$, which corresponds to the limit $h^{*}\to\infty$, the cooling condition does not lead to an infinite initial solidification rate [18, 21, 22, 23].

The evolution of the interface is determined by an energy balance which is known as the Stefan condition In its simplest form, this is

[TABLE]

where $L^{*}_{m}$ is the latent heat of fusion.

Finally, at $t^{*}=0$ we assume that no solid has been produced,

[TABLE]

2.2 Effective Fourier Law

For a system of characteristic size $L^{*}$, Alvarez and Jou [15] propose substituting the bulk thermal conductivity $k^{*}$ by an ETC of the form

[TABLE]

where $\ell^{*}$ is the phonon mean free path. The ratio $\ell^{*}/L^{*}$ is the Knudsen number and it determines the dominance of non-local effects in heat transport. In the limit $L^{*}\ll\ell^{*}$ (Kn$\gg 1$), Eq. (7) reduces to $k^{*}\approx 2k^{*}_{0}L^{*}/\ell^{*}$, which is in accordance with experimental data [24]. Conversely, in the limit $\ell^{*}\gg L^{*}$ (Kn$\ll 1$) the classical limit $k^{*}_{\text{eff}}=k^{*}$ is retrieved. In the current problem the system length-scale is determined by the size of the solid region, hence $k^{*}_{\text{eff}}(L^{*})=k^{*}_{\text{eff}}(s^{*})$ and we obtain

[TABLE]

which we assume to govern heat conduction through the solid. Upon combining Eqs. (1) and (8) we obtain

[TABLE]

The temperature at the solid-liquid interface is determined by Eq. (3), whereas at $x^{*}=0$ Eq. (4) becomes

[TABLE]

The Stefan condition may now be written as

[TABLE]

2.3 Dimensionless formulation

The temperature scale is given by the temperature jump at $x^{*}=0$, that is, $\Delta T^{*}=T^{*}_{\text{f}}-T^{*}_{\text{e}}$. The natural length scale of the problem is the mean free path $\ell^{*}$ and hence we choose the corresponding diffusive time scale $c^{*}\rho^{*}\ell^{*2}/k^{*}$, which balances the terms in Eq. (1). Consequently, we introduce the dimensionless quantities $x=x^{*}/\ell^{*}$, $s=s^{*}/\ell^{*}$, $t=k^{*}t^{*}/(c^{*}\rho^{*}\ell^{*2})$ and $T=(T^{*}-T^{*}_{\text{f}})/\Delta T^{*}$. Upon writing derivatives as indices, the dimensionless equations are

[TABLE]

where $\beta=L^{*}_{m}/(c^{*}\Delta T)$ and $\text{Bi}=h^{*}\ell^{*}/k^{*}$ are the Stefan number and the Biot number respectively. In melting problems, the latter is also called Nusselt number because the material in contact with the environment is in liquid state. Finally, $f(s)=2s\left(\sqrt{s^{2}+1}-s\right)$ is the non-dimensional form of $k^{*}_{\text{eff}}$.

Note, due to the way in which the problem has been scaled, the mean free path has been eliminated from the ETC and is only present in the Biot number. In fact, the dimensionless position of the interface $s(t)$ can now be understood as an effective Knudsen number depending on time, since it corresponds to the ratio of the size of the growing solid to the mean free path.

The effect of the Stefan number on the evolution of the interface in formulations with a fixed temperature has been studied by many authors; see the books by Gupta [22] or Alexiades and Solomon [23], for example. However, when Newton cooling conditions are applied, the behaviour of the system when varying Bi has to be studied as well.

2.4 Parameter estimation

In this paper we will use the thermophysical parameters of silicon, since it is material that is widely used in theoretical studies [15, 17, 24, 25, 26]. The values of the relevant physical quantities are given in Table 1.

The Stefan number can be parametrized in terms of the temperature change as $\beta=\mathcal{T}^{*}/\Delta T$ and where $\mathcal{T}^{*}=L^{*}_{m}/c^{*}$. Using the values given in Table 1 we find $\mathcal{T}^{*}\approx 2754$ K. Therefore, even for a temperature drop of 100 K we still obtain $\beta\approx 27.5$. For other materials such has tin, lead or gold, we find that $\mathcal{T}^{*}$ is of the order of hundreds of Kelvin [18, 20, 28], and hence we can still expect large Stefan numbers for temperature changes of tens of Kelvin. However, since the liquid is assumed to be initially at the phase change temperature, only a small temperature change is needed to drive the solidification process and thus we expect $\beta\gg 1$.

The Biot number can be expressed as $\text{Bi}=h^{*}/\mathcal{H}^{*}$, where for silicon $\mathcal{H}^{*}=k^{*}/\ell^{*}\approx 3.4\times 10^{9}$ W/m${}^{2}\cdot$K. Determining the heat transfer coefficient $h^{*}$ is complicated, since its value depends on the environment and the material. Nonetheless, there is a maximal value $h^{*}_{\text{max}}$ beyond which the material would simply vaporise. Ribera and Myers [18] provide the expression $h_{\text{max}}^{*}=c^{*}\sqrt{\rho^{*}B^{*}/3}$, where $B^{*}$ is the bulk modulus. For silicon we have $B^{*}\approx 10^{11}$ kg/m$\cdot$s2 [29], therefore $h^{*}_{\text{max}}\approx 7.6\times 10^{9}$ W/m${}^{\text{2}}\cdot$K and hence $\text{Bi}_{\text{max}}\approx 2.2$. However, the order of magnitude of $h^{*}_{\text{max}}$ is extremely large and is never reached in practical situations [30], thus $h^{*}\ll h^{*}_{\text{max}}$ generally.

Hence, throughout this study we will assume $\beta\gg 1$ and $\text{Bi}\ll 1$. If necessary, the relative size of $\beta$ to Bi will be discussed during the analysis.

3 Solution methods

There exist few analytical solutions to practical Stefan problems. However, in the case of a constant thermal conductivity, i.e. $f(s)=1$, and a fixed temperature condition at $x=0$, it is possible to obtain an exact solution, termed the Neumann solution [22, 23, 31].

Such a solution does not exist with a Newton cooling condition, not even in the case of constant thermal conductivity, and thus we need to explore numerical and approximate methods to solve the problem with $f\neq 1$. Furthermore, the expected sizes of the Biot and Stefan numbers allow us to perform an asymptotic analysis to distinguish different time regimes of interest and obtain analytical expressions in most of them.

3.1 Numerical solution

A usual approach to obtain numerical solutions of Stefan problems is to first introduce an alternative variable which transforms the moving domain into the unit interval. The numerical scheme to solve the resulting problem consists of discretizing explicitly for the temperature and implicitly for the interface position and the interface speed [17, 18, 28, 32]. In addition, since the solid does not exist initially, a small time analysis must be performed to obtain a valid initial condition for the numerical scheme.

3.1.1 Boundary fixing transformation

We define the new variable $\xi=x/s$ and rewrite the temperature as $T(x,t)=u(\xi,t)$, which transforms Eqs. (12a)–(12d) into

[TABLE]

where we have defined $F(s)=2(\sqrt{1+s^{2}}-s)$ for simplicity. This formulation breaks down as $t\to 0$ and hence the initial condition for $s$ must be substituted by an approximate value at a small time $t_{0}$.

3.1.2 Small-time solution

Since the solid phase does not exist initially, a small-time analysis has to be carried out to study the initial dynamics of the system. To facilitate the analysis, we assume the remaining dimensionless numbers are order one in magnitude. In this way, the small-time behaviour applies to more limits that may arise in other situations where, for instance, the Stefan number is expected to be small [33]. An analysis of the results for asymptotic limits of these parameters can be performed afterwards. Let now $t\ll 1$ and assume that the solid-liquid interface is approximated by an expression of the form $s\approx\lambda t^{p}$ (hence $s_{t}\approx p\lambda t^{p-1}$), where $\lambda$ and $p$ are constants to be determined. Neglecting terms of order $t^{2p}$ gives $F\approx 2(1-\lambda t^{p})$ and thus Eq. (13d) can be written as

[TABLE]

where we have used $(1-z)^{-1}\approx 1+z$ for $z\ll 1$ and $u_{\xi}$ is evaluated at $\xi=1$. The only possible way to balance with the r.h.s. of Eq. (14) is obtained by setting $p=1$, hence

[TABLE]

provided $t\ll 1$. In particular, Eq. (15) predicts that the solid initially grows at a finite rate $\lambda$ independent of the choice of boundary condition at $x=0$. The fixed boundary condition is known to predict an infinite phase change rate [18, 21, 22, 23, 28, 32], so the size-dependent ETC produces a more realistic growth rate.

Substituting Eq. (15) into Eqs. (13a)–(13c) and taking the limit $t\to 0$ yields a second order boundary value problem for $u$

[TABLE]

The solution to (16) is

[TABLE]

where $\text{erf}(z)=2\pi^{-1/2}\int_{0}^{z}\exp(-t^{2})dt$ is the error function. Finally, $\lambda$ is determined by substituting Eqs. (15) and (17) into Eq. (14) and taking the limit $t\to 0$, which yields an equation for $\lambda$

[TABLE]

For a fixed value of Bi it is trivial to find $\beta(\lambda)$, which can be used to plot the exact solution on the $(\beta,\lambda)$-plane. The same argument is valid for $\text{Bi}(\lambda)$ for a fixed $\beta$. However, finding $\lambda$ for fixed values of Bi and $\beta$ requires the use of numerical methods.

Recall, Eq. (18) has been obtained under the assumption $\text{Bi},\beta=O(1)$. Since we are interested in the limits $\text{Bi},\beta^{-1}\ll 1$, we can perform an asymptotic analysis on Eq. (18) to find an approximate solution in this specific limit.

In the limit $\text{Bi},\text{Bi}/\beta\to 0$, Eq. (18) yields $\lambda\to 0$, which indicates $\lambda=O(\delta)$ for some $\delta=\delta(\text{Bi},\beta)\ll 1$. We can make the approximations

[TABLE]

which reduces Eq. (18) to

[TABLE]

Note, the r.h.s of Eq. (20) represents the driving force of the solidification process, therefore a balance with the l.h.s is required. Since $\lambda\ll 1$, this balance is only achieved if we choose $\lambda=O(\text{Bi}/\beta)$. Then we seek a solution for $\lambda$ as a power series

[TABLE]

where $\delta=\text{Bi}/\beta$. Introducing this expansion into Eq. (18) yields a subproblem for each power of $\delta$. The first- and second-order subproblems are

[TABLE]

which gives $\lambda_{0}=2/(2+\text{Bi})$ and $\lambda_{1}=-(3+\text{Bi})\lambda_{0}^{3}/12$ and hence

[TABLE]

In Fig. 2 we compare this approximation to the exact solution of Eq. (18) and observe an excellent agreement across the entire range $\text{Bi},\beta^{-1}\in(10^{-3},1)$.

Consider now the first-order approximation. Substituting into Eq. (17) and using Eq. (19) we can reduce $u$ to

[TABLE]

which may be used to approximate the temperature profile at a small time $t_{0}$ for $\text{Bi},\beta^{-1}\ll 1$.

3.1.3 Numerical scheme

The problem in the fixed domain is now discretized on a spatial grid of $N+1$ points of the form $\xi_{i}=\frac{i}{N}$ and a time grid with points of the form $t_{n}=t_{0}+n\Delta t$. In interior points of the domain, derivatives are approximated using second order central differences:

[TABLE]

where $u^{n}_{i}=u(\xi_{i},t_{n})$ and $\Delta\xi=N^{-1}$. Upon substituting these expression into Eq. (13a) and writing $s=s^{n}$ and $\mathrm{d}s/\mathrm{d}t=s_{t}^{n}$ we obtain $N-1$ algebraic equations of the form

[TABLE]

with coefficients

[TABLE]

Using Eq. (13c) we obtain

[TABLE]

whereas the Newton condition at $\xi=0$ becomes, after using second order forward differences,

[TABLE]

where $\alpha^{n}=F(s^{n})/(\text{Bi}\Delta\xi)$. To update the position of the moving boundary we discretize Eq. (13d) explicitly for $s$ to find $s^{n+1}$,

[TABLE]

where we have used second order backward differences to discretize $u_{\xi}$ at $\xi=1$. To avoid stability issues, we discretize time as a logarithmically spaced grid. In this way, the time steps are smaller at the beginning of the process, where the growth rate is expected to be larger. Numerical experimentation shows that no unstabilities appear for the considered range of parameters.

3.2 Asymptotic solution

In [17], an asymptotic analysis for the formulation with a fixed-temperature condition is performed under the assumption of a large Stefan number. In the case of the Newton cooling condition, the asymptotic analysis depends also on Bi and potentially on the relative size of Bi to $\beta$. In order to simplify the analysis we will focus only on three time regimes, which correspond to $s\ll 1$, $s=O(1)$ and $s\gg 1$.

The first time regime is given by $t=O(\varepsilon)$, where $\varepsilon\ll 1$ is an artificial parameter, and it has already been partially studied during the small-time analysis. It describes the initial stage of the solidification, where the solid has not grown much yet, therefore $x,s=O(\epsilon_{1})$ for some small $\epsilon_{1}\ll 1$. In addition, the small Biot number indicates that the influx of heat is very small and therefore the temperature is not expected to differ much from its initial value, $T=O(\epsilon_{2})$ for $\epsilon_{2}\ll 1$. Upon defining the scaled variables $t=\varepsilon\hat{t}$, $x=\epsilon_{1}\hat{x}$, $s=\epsilon_{1}\hat{s}$ (thus $f\approx 2\epsilon_{1}\hat{s}$) and $T=\epsilon_{2}\hat{T}$, by balancing terms in the Newton condition we find $\epsilon_{2}=\text{Bi}$. To ensure that solidification occurs at $\hat{t}>0$ we need to balance both sides of Eq. (12d), which yields $\epsilon_{1}=\text{Bi}\beta^{-1}\varepsilon\ll\varepsilon$. At the leading order, the system in the new variables reads

[TABLE]

After applying the boundary conditions (31b) and (31c), the solution to Eq. (31a) is

[TABLE]

and therefore the position of the interface is determined by

[TABLE]

In the original non-dimensional variables we find

[TABLE]

for $t\ll 1$. Notice that this temperature profile coincides with the small-time profile given in Eq. (24) when terms of order $\text{Bi}^{2}$ are neglected.

The second time regime of interest is when $s=O(1)$, which invalidates the approximation $f\approx 2s$. From Eq. (34) we find $t=O(\text{Bi}^{-1}\beta)$ and $T=O(\text{Bi})$. Upon defining the new variables $x=x^{\prime}$, $s=s^{\prime}$, $t=\text{Bi}^{-1}\beta t^{\prime}$ and $T=\text{Bi}T^{\prime}$, the leading order problem for $T^{\prime}$ becomes

[TABLE]

whose solution is

[TABLE]

The solid-liquid interface is therefore determined by

[TABLE]

subject to the matching condition

[TABLE]

In the original dimensional variables the solution in this time regime is therefore

[TABLE]

for $t=O(\text{Bi}^{-1}\beta)$, which captures the previous time regime as well due to the behaviour of $f$ for small values of $s$.

The third time regime of interest is when $s\gg 1$ and therefore $f\approx 1$. The correct balance is obtained by choosing $t=O(\text{Bi}^{-2}\beta)$, $x,s=O(\text{Bi}^{-1})$ and $T=O(1)$, which yields, upon writing $t=\text{Bi}^{-2}\beta\tilde{t}$, $x=\text{Bi}^{-1}\tilde{x}$, $s=\text{Bi}^{-1}\tilde{s}$ and $T=\tilde{T}$,

[TABLE]

The solution to (40) is given by

[TABLE]

and thus the interface is determined by

[TABLE]

which can be integrated to give

[TABLE]

where $C$ is a constant of integration to be determined to match the previous time regime. Using Eq. (39) we find $\tilde{s}\sim\tilde{t}$ for $\tilde{t}=O(\text{Bi})$, which gives $C=0$ at leading order. In the original dimensionless variables the temperature profile is

[TABLE]

To obtain $s$ we must invert Eq. (44b) numerically. Equivalently, it can be calculated integrating

[TABLE]

which captures all the previous time regimes and can thus be solved using the initial condition $s(0)=0$.

We can consider a fourth time regime, where the Newton condition converges to the fixed temperature condition, as it can be observed by assuming $\tilde{x}\gg 1$ in Eq. (40b). However, observe that Eq. (45) already captures this regime, since $s\to\infty$ yields $f\to 1$ and $s_{t}\propto s^{-1}$, from where the classical behaviour $s\sim\sqrt{t}$ is recovered.

4 Results and discussion

In Fig. 3 we show the numerical and asymptotic solutions of the non-classical formulation and compare them against the classical solution for constant thermal conductivity. This is done for different choices of Bi and $\beta$.

Firstly, we observe that the asymptotic solution is in good agreement with the numerical solution. Only in Fig. 3(a) do some discrepancies appear for larger times. This is due to the ratio $\text{Bi}/\beta$, whose value in the case of Fig. 3(a) is at least one order of magnitude larger than in the other cases. Recall that the asymptotic solution is based on the assumption $\text{Bi}/\beta\ll 1$ and hence we should account for higher-order terms to increase the accuracy of the asymptotic solution.

Secondly, we observe that the speed of the solid-liquid interface decreases as we decrease the Biot number and the usual square-root profile transforms into a linear one. Furthermore, as we decrease Bi both classical and non-classical formulations tend to the same solution and non-local effects described by the size-dependent ETC seem to become less important. For $\text{Bi}=0.01$, which we do not show here, both the classical and non-classical profiles are indistinguishable for $\beta\in(10,100)$. In fact, our asymptotic analysis indicates that we expect a linear profile until $t=O(\beta/\text{Bi})$ with a slope $\text{Bi}/\beta$. The same slope is obtained in the small-time analysis of the classical formulation (see A). Physically, the convergence of both solutions to the same profile can be understood by recalling that small Biot numbers correspond to a situation where heat flow is limited by the heat exchange with the environment rather than heat conduction through the bulk [34], which implies that the form of thermal conductivity employed should not matter in this case.

In Fig. 3(c) it can also be observed that large deviations between the non-classical and classical formulations appear for large values of $\beta$. This may be understood by recalling that the Stefan number corresponds to a slow solidification, since $s_{t}\propto\beta^{-1}$. Hence, non-local effects remain for a longer time period due to the small growth rate, which causes these larger deviations between both formulations.

To illustrate the general case of an arbitrary Biot number, in Fig. 4 we show the evolution of the solid-liquid interface for different values of the Stefan and Biot numbers according to the classical and non-classical formulations. For $\beta=10$, differences between both formulations are considerably less than for the larger value $\beta=100$. In Fig. 5 we have plotted the evolution in time of the absolute difference of both formulations for different values of $\beta$ and Bi. We observe the existence of two regimes: in the first regime, due to the presence of non-local effects, differences between the classical and modified formulations increase, whereas in the second time regime they decrease as non-local effects disappear. By observing Figs. 4(a) we see that this change in behaviour occurs when $s\approx 1$. In the case of $\beta=100$ the second time regime enters later due to the fact that solidification is slower and hence the presence of non-local effects is significant for a longer period of time.

Hence, these results suggest that non-local effects become less important if we decrease the Biot number, which corresponds to a poor heat conduction through the solid, or by decreasing the Stefan number, which corresponds to fast solidification. In the first case, the choice of thermal conductivity does not alter the evolution of $s(t)$ significantly. In the latter, non-local effects disappear due to the fast growth of the solid phase.

The effect of the Biot number on the temperature can be observed in Fig. 7, which shows the evolution of the temperature profiles according the non-classical and classical formulations for different values of Bi and for $\beta=10$. For early times the size of the solid is small, which leads to a small ETC and hence changes in temperature in the non-classical formulation occur mainly near $x=0$. This happens independent of the value of Bi, as it can be seen by comparing panels (a), (c) and (e). In addition, it can seem that diffusion is enhanced as Bi decreases, in contradiction to the fact that conduction becomes worse in this limit. We must recall that the spatial coordinate used in Fig. 7 is the alternative variable $\xi=x/s$, and hence the growth rate of the solid is not shown here. Hence, the temperature profiles must be observed together with the growth of the solid for these values of Bi, shown in Figs. 4(a). Another indicator for a worse conduction for small values of Bi is the range of temperatures shown in the color bars of the corresponding plots.

In the classical formulation, diffusion through the solid occurs normally and changes in temperature appear gradually through the crystal. Similarly as in the case of the non-classical formulation, panels (b), (d) and (f) of Fig. 7 must be observed by having in mind the corresponding growth rates of the solid-liquid in, shown in Fig. 4(b).

We can also see that the initial differences between the temperature profiles for a fixed Bi disappear as time increases, which is caused by the fact that the growth of the solid leads to an increase of the ETC, which finally converges to the classical conductivity.

5 Conclusions

In this paper we have formulated a mathematical model for a one-dimensional solidification process with a size-dependent thermal conductivity and a Newton cooling condition at the fixed boundary. Numerical and asymptotic solution methods have been proposed and compared, showing good agreement for the range of parameters considered. In contrast to previous studies, we have studied the limit of a small Biot number, which seems more reasonable after a proper analysis. We have observed that non-local effects tend to disappear as the Biot number decreases, which seems surprising since the initial solidification rate is proportional to Bi and thus non-local effects would have been expected to keep present for longer periods of time. This result has been attributed to the poor conduction in both the classical and the modified formulations due to the Biot number being small, since this implies that the heat exchange with the environment at the fixed boundary dominates over the heat conduction through the growing solid. Similarly, as we decrease the Stefan number the presence of non-local effects is reduced to an initial stage of the process that is overcome rapidly due to the higher speed of the solid-liquid interface.

Acknowledgements

M. C. acknowledges that the research leading to these results has been funded by “La Caixa” Foundation and by the CERCA programme of the Generalitat de Catalunya. The author thanks T. G. Myers and M. G. Hennessy for a critical reading of the manuscript during its preparation.

Appendix A Small-time behaviour of the classical Stefan problem with Newton cooling conditions

As in the small-time analysis for the non-classical case, we assume $\text{Bi},\beta=O(1)$. After applying the transformation $T(x,t)=u(\xi,t)$, the classical problem with Newton conditions takes the form

[TABLE]

Let $t\ll 1$. Assuming only $s=O(t^{p})$ requires $p=0.5$ to balance terms in the Stefan condition, which is unphysical since we recover the initially infinite solidification rate which should be avoided due to the Newton condition. Furthermore, only scaling $s$ transforms the Newton condition into $u_{\xi}\approx 0$ and therefore no heat exchange would occur at the boundary to drive the solidification process. In fact, balance between both sides of Eq. (46b) can only happen if $u=O(t^{p})$. Assuming $u,s=O(t^{p})$ in Eq. (46d) implies that $p=1$ is required for balancing terms. Hence, for an arbitrary $\varepsilon\ll 1$ we define the scaled variables $u=\varepsilon\bar{u}$, $s=\varepsilon\bar{s}$, $t=\varepsilon\bar{t}$ and $\xi=\bar{\xi}$. The Stefan condition becomes

[TABLE]

where $\bar{u}_{\bar{\xi}}$ is evaluated at $\bar{\xi}=1$. Upon neglecting terms of order $\varepsilon$, the problem for $\bar{u}$ is

[TABLE]

which gives $\bar{u}(\bar{\xi})=\bar{s}\text{Bi}(\bar{\xi}-1)$. Substituting this expression into Eq. (47) yields the initial value problem

[TABLE]

whose solution is $\bar{s}(\bar{t})=\text{Bi}\beta^{-1}\bar{t}$. In the original dimensionless variables, the small-time behaviour of the classical problem is therefore

[TABLE]