The Stefan problem with variable thermophysical properties and phase change temperature

TL;DR

This paper develops a generalized Stefan problem model incorporating variable thermophysical properties and size or velocity-dependent phase change temperatures, with applications to nanoparticle melting and dendrite formation, validated through numerical simulations.

Contribution

It introduces a novel formulation of the Stefan problem that accounts for property variations and dynamic phase change temperatures, addressing limitations in existing models.

Findings

The model accurately predicts nanoparticle melting behavior.

Finite difference simulations validate the theoretical approach.

Highlighting errors in previous literature on phase change modeling.

Abstract

In this paper we formulate a Stefan problem appropriate when the thermophysical properties are distinct in each phase and the phase-change temperature is size or velocity dependent. Thermophysical properties invariably take different values in different material phases but this is often ignored for mathematical simplicity. Size and velocity dependent phase change temperatures are often found at very short length scales, such as nanoparticle melting or dendrite formation; velocity dependence occurs in the solidification of supercooled melts. To illustrate the method we show how the governing equations may be applied to a standard one-dimensional problem and also the melting of a spherically symmetric nanoparticle. Errors which have propagated through the literature are highlighted. By writing the system in non-dimensional form we are able to study the large Stefan number formulation and an energy-conserving one-phase reduction. The results from the various simplifications and assumptions are compared with those from a finite difference numerical scheme. Finally, we briefly discuss the failure of Fourier's law at very small length and time-scales and provide an alternative formulation which takes into account the finite time of travel of heat carriers (phonons) and the mean free distance between collisions.

Tables

Table 1. Table 1: Values of the physical parameters for tin and gold.

	Symbol	Tin	Gold
Density of liquid kg/m3	${\rho }_l$	6980	17300
Density of solid kg/m3	${\rho }_s$	7180	19300
Thermal conductivity of liquid (bulk) W/(m$\cdot$K)	$k_l$	30	106
Thermal conductivity of solid (bulk) W/(m$\cdot$K)	$k_s$	67	317
Specific heat of liquid J/(kg$\cdot$K)	$c_l$	268	163
Specific heat of solid J/(kg$\cdot$K)	$c_s$	230	129
Melt temperature (bulk) K	$T_m$	505	1337
Latent heat (bulk) J/kg	$L_m$	58500	63700
Solid-liquid surface energy N/m	${\sigma }_{sl}$	0.064	0.27
Vapour-liquid surface energy N/m	${\sigma }_{la}$	0.589	1.15
Maximum heat transfer coefficient W/(m${}_{}{}^2\cdot$K)	$h_{\text{max}}$	$4.7\times {10}^9$	$4.9\times {10}^9$

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0\text{leavevmode},$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{v}\right)+\nabla \cdot \left(\rho \mathbf{v}\otimes \mathbf{v}+p\underline{\mathbf{I}}\right)=0,$$

$$\frac{\partial }{\partial t}\left(\frac{1}{2}\rho v^2\right)+\nabla \cdot \left(\frac{1}{2}\rho v^2\mathbf{v}+p\mathbf{v}\right)-p\nabla \cdot \mathbf{v}=0,$$

$$\frac{\partial }{\partial t}\left(\rho \left[u+\frac{1}{2}{v}^2\right]\right)+\nabla \cdot \left(\rho \left[u+\frac{1}{2}{v}^2\right]\mathbf{v}+\mathbf{q}+p\mathbf{v}\right)=0.$$

$$\nabla \cdot \mathbf{v}=0.$$

$$\rho \left(\frac{\partial u}{\partial t}+\mathbf{v}\cdot \nabla u\right)=-\nabla \cdot \mathbf{q}.$$

$$\frac{\partial F}{\partial t}+\nabla \cdot \mathbf{G}=0,$$

$$[F{]}_{-}^{+}{v}_n=[\mathbf{G}\cdot \hat{\mathbf{n}}{]}_{-}^{+},$$

$$({\rho }_l-{\rho }_s)v_n=\left({\rho }_l{\mathbf{v}}_l-{\rho }_s{\mathbf{v}}_s\right)\cdot \hat{\mathbf{n}}$$

$${\rho }_l\left(u_l+\frac{v_l^2}{2}\right)(v_n-{\mathbf{v}}_l\cdot \hat{\mathbf{n}})-{\rho }_s\left(u_s+\frac{v_s^2}{2}\right)(v_n-{\mathbf{v}}_s\cdot \hat{\mathbf{n}})=\left[({\mathbf{q}}_l-{\mathbf{q}}_s)+(p_l{\mathbf{v}}_l-p_s{\mathbf{v}}_s)\right]\cdot \hat{\mathbf{n}}$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{v}\cdot \hat{\mathbf{n}}\right)+\nabla \cdot \left[\left(\rho \mathbf{v}\otimes \mathbf{v}+p\underline{\mathbf{I}}\right)\cdot \hat{\mathbf{n}}\right]=0$$

$$\left[\left({\rho }_l{\mathbf{v}}_l-{\rho }_s{\mathbf{v}}_s\right)\cdot \hat{\mathbf{n}}\right]v_n$$

$$({\rho }_l-{\rho }_s)s_t$$

$$p_s-p_l={\rho }_l{v}_l(v_l-s_t)$$

$${\rho }_l\left(u_l+\frac{v_l^2}{2}\right)(s_t-v_l)-{\rho }_s{u}_s{s}_t$$

$$\underline{\mathbf{K}}\sim ϵ(\nabla \varphi \otimes \nabla \varphi )$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{v}\right)+\nabla \cdot \left(\rho \mathbf{v}\otimes \mathbf{v}+p\underline{\mathbf{I}}\right)+\sigma \kappa \delta (r-s)\hat{\mathbf{n}}=0,$$

$$\frac{\partial }{\partial t}\left(\frac{1}{2}\rho v^2\right)+\nabla \cdot \left(\frac{1}{2}\rho v^2\mathbf{v}+p\mathbf{v}\right)-p\nabla \cdot \mathbf{v}+\sigma \kappa \delta (r-s)v_n=0,$$

$$\frac{\partial }{\partial t}\left(\rho \left[u+\frac{1}{2}{v}^2\right]\right)+\nabla \cdot \left(\rho \left[u+\frac{1}{2}{v}^2\right]\mathbf{v}+\mathbf{q}+p\mathbf{v}\right)+\sigma \kappa \delta (r-s)v_n=0.$$

$$p_s-p_l$$

$${\rho }_l\left(u_l+\frac{v_l^2}{2}\right)(s_t-v_l)-{\rho }_s{u}_s{s}_t$$

$$\left[{\rho }_s\left(u_l-u_s+\frac{v_l^2}{2}\right)-{\sigma }_{sl}\kappa -\left(1-\frac{{\rho }_s}{{\rho }_l}\right)p_l\right]s_t$$

$$u_s$$

$$u_l$$

$$\mathrm{d}u=c\mathrm{d}T.$$

$$\mathbf{q}=-k\nabla T$$

$${\rho }_s{c}_s\frac{\partial T_s}{\partial t}$$

$${\rho }_l{c}_l\left(\frac{\partial T_l}{\partial t}+{\mathbf{v}}_{\mathbf{l}}\cdot \nabla T_l\right)$$

$${\rho }_s\left[(c_l-c_s)(T_I-T_m)+L_m-\frac{1}{2}\left(1-{\left(\frac{{\rho }_s}{{\rho }_l}\right)}^2\right)s_t^2\right]s_t$$

$$v_l(t)=\left(1-\frac{{\rho }_s}{{\rho }_l}\right)\frac{\mathrm{d}s}{\mathrm{d}t}.$$

$$\frac{\mathrm{d}s}{\mathrm{d}t}=C_0{\mathrm{e}}^{-C_1/T_I}(T_m-T_I),$$

$$T_I=T_m-\Phi s_t\text{leavevmode}.$$

$$({\rho }_l-{\rho }_s)\frac{\mathrm{d}R}{\mathrm{d}t}={\rho }_l{v}_l(R,t)\text{leavevmode}.$$

$$v_l(r,t)=\left(1-\frac{{\rho }_s}{{\rho }_l}\right)\frac{R^2}{r^2}\frac{\mathrm{d}R}{\mathrm{d}t}.$$

$$R_b^3=\frac{{\rho }_s}{{\rho }_l}{R}_0^3+\left(1-\frac{{\rho }_s}{{\rho }_l}\right)R^3\text{leavevmode},$$

$$\displaystyle\frac{\partial T_{s}}{\partial r}\bigg{|}_{r=0}=0.$$

Full text

1 Introduction

The classical Stefan problem is described by two heat equations, one in each of the two material phases. The domain over which each equation holds is determined by an energy balance, the Stefan condition, which equates the energy released during the phase change with that conducted away through either phase. Thermophysical properties are assumed constant (often taking the same value in each phase) and the melt temperature is also constant. Whilst exact solutions do exist for certain configurations and boundary conditions, there is possibly just one physically useful one: the Neumann solution. This involves a semi-infinite material with constant thermophysical properties changing phase due to an instantaneous switch in the boundary temperature, which subsequently stays at the same, constant value. The classic problem, and Neumann solution, are described in many texts, see [1, 15, 29] for example.

In practice an instantaneous switch to a new temperature is impossible to achieve and it is difficult to maintain a constant boundary temperature throughout a phase-change process. Thermophysical properties are temperature dependent (admittedly sometimes only weakly so), but they always differ between phases (often quite significantly). There also exist situations where properties are size or velocity dependent. This may occur, for example, when the curvature is high, such as in dendritic growth or nanoscale phase change, or with the solidification of supercooled materials.

The size dependence of material properties at small length-scales has been demonstrated experimentally, by molecular dynamics simulations, and also theoretically. One of the main reasons for this dependence is the large ratio of surface to interior atoms. Interior atoms are free to bond with other atoms on all sides. Surface atoms are exposed to the environment and do not have the ability to bond on all sides, consequently they behave differently to interior atoms. At the macro-scale the very high ratio of interior to surface atoms means that the interior dominates the material’s behaviour. However, for sufficiently small samples the ratio decreases until surface properties become important. At this stage size dependence must be considered in any material description. This argument is equally valid whenever the surface curvature is sufficiently high that the local surface to bulk ratio is large.

Properties that vary with size include the catalytic, ferromagnetic, and mechanical characteristics, colour, surface tension, phase change temperature, and latent heat [23, 39]. For example, colloidal gold with particles smaller than 40 nm has a colour ranging from clear to yellow, between 40–100 nm it has a reddish hue, between 100–120 nm the colour becomes purple. In terms of phase change, experiments have shown a decrease in the melt temperature of approximately 500 K below the bulk value for gold nanoparticles with radii of the order of 1 nm [8] and 70 K and 200 K in tin and lead nanoparticles, respectively [14]. Molecular dynamics simulations have shown a decrease of 800 K (approximately 60%) for gold nanoparticles of radius 0.8 nm [45]. It has also been suggested that latent heat and surface tension decrease with size. For tin decreases in the latent heat of the order 70% below the bulk value were reported by [33] although subsequent refinements to the technique showed even greater reductions [31]. Sun and Simon [49] concluded that the melt temperature is well approximated by the Gibbs–Thomson relation while the latent heat decrease is much greater than theoretically predicted. Ribera and Myers [43] review theoretical predictions for latent heat reduction and conclude that none of the models examined match the experimental data. Decreases in surface tension are usually less significant than those observed in the latent heat and melt temperature, typically of the order 15% below the bulk value for $R=5$ nm [50]. With supercooling the melt temperature $T_{I}$ varies nonlinearly with the fluid velocity. In [19] examples are provided for copper, where the front velocity varies between 0 and 2.9 m/s as the degree of supercooling varies between 0 and 1000 K.

There has been particular focus on spherical nanoparticle melting, since this provides a simple framework to investigate the effect of size-dependent properties. The theoretical investigations on spherical nanoparticle phase change [4, 5, 20, 35, 52, 53] have the same density in both phases but include melting point depression. The mathematical model in these papers involves standard heat equations and a Stefan condition with an ‘effective’ latent heat of the form $L_{m}+\Delta c(T_{I}-T_{m})$, where $T_{m}$ is the bulk melt temperature, $\Delta c$ is the jump in specific heat between the liquid and solid phases and $T_{I}(t)$ the temperature at the phase change interface. This expression is also employed in studies of the solidification of supercooled materials [1, 40]. Melting point depression and constant (but different) densities in each phase are included in the model of [21]. The change in density forces fluid motion, so introducing advection into the liquid heat equation and kinetic energy into the Stefan condition. The generation of kinetic energy means there is less energy available for the phase change which results in significantly longer melt times. The model of [21] is an extension of the Cartesian version of [1] which, due to a lack of curvature effects, does not include melting point depression. The thesis of Back [3, §7.1-7.4] discusses a number of variations to the standard Stefan problem including the advection/kinetic energy terms. A model is also discussed including a size-dependent latent heat, following the formula of [33], which is released at the bulk melt temperature. In [39] the governing equations include an effective latent heat which is the sum of the size-dependent latent heat, the kinetic energy, and the energy required to make new surface. Surface energy is also included in the model of [15]. In [43] it is shown that the formula of [33] underestimates the value of latent heat near the bulk value (where the data is most reliable). They also demonstrate that other formulations, such as those of [46, 54], also provide a poor match with experimental data. They go on to propose an exponential fit to the data for tin nanoparticles which provides excellent agreement in the limit of large particles and only shows a noticeable discrepancy for $R<8$ nm. This fit is used in the model developed in [39] to demonstrate significantly faster melt times than those with a constant latent heat. However, the well-known liquid skin model [33] indicates that the experimental calculation of latent heat does not coincide with the actual latent heat release. This is because the liquid skin is thought to form spontaneously, with a resultant decrease in surface energy and also an energy cost in forming the skin. Thus, the experimentally measured latent heat should be interpreted as an effective latent heat rather than the actual value for the nanoparticle.

The Stefan condition, which is simply an energy balance, is derived in different ways. Obviously one can examine a small change around the front and let the time taken for this change to tend to zero, to provide the differential form. Davis [15] applies this approach on an arc of surface, taking into account the latent heat release and energy required for creating new surface, hence introducing surface tension. Fedorov and Shulkin [17] simply work with the definition of latent heat in terms of internal energy in the current and bulk configurations to arrive at the effective latent heat $L_{m}+\Delta c(T_{I}-T_{m})$. This then replaces the standard latent heat in the Stefan condition. The same expression is obtained by [1] through the conservation equations defined in [6]. If density change between phases is included the conservation route also naturally leads to the inclusion of kinetic energy (with or without curvature) in [1]. Gupta [26] defines a different effective latent heat $L_{m}+[(\rho_{l}/\rho_{s})c_{l}-c_{s}]T_{m}$ but neglects kinetic energy. A model including curvature and hence surface energy is derived in [39], they also go on to show that the term $\Delta c(T_{I}-T_{m})$ plays a very minor role, a result confirmed in [43].

A common issue with studies of phase change is the applied boundary conditions. The simplest boundary condition to impose in thermal and Stefan problems is the fixed-temperature condition, where the material, which is initially at a constant temperature, has the boundary temperature instantaneously changed to a different, constant value. Clearly this is unphysical and chosen primarily for mathematical simplicity. The work described in [21] showed that, when including density variation, the initial infinite melt velocity which results from this condition had a significant effect on melt times, and that this effect was more noticeable with smaller particles. However, even up to the macroscale a discrepancy of approximately 15% was found between the single- and two-density models for gold particles. In [43] it is shown that when a Newton cooling boundary condition is employed this discrepancy disappears, even for very high values of the heat transfer coefficient, $h$. Note, they also point out that there is a limit on $h$, beyond which the material would simply vaporise: the mathematical limit $h\rightarrow\infty$, which reduces Newton cooling to the fixed-temperature condition, is not physically possible. Their calculations use the maximum value for $h$. Since it is orders of magnitude larger than typical values the results are perhaps rather forgiving in that they are relatively close to the fixed-temperature results. In reality much larger differences from the fixed-temperature boundary condition should be observed.

In the following we will discuss phase change with certain variable properties and Newton’s law of cooling at the boundary. The Newton condition is more realistic than a fixed-temperature condition, which may be reproduced by introducing an infinite heat transfer coefficient. When a specific configuration is required we will focus on the melting of a spherically symmetric nanoparticle. The melting of a spherically symmetric nanoparticle is the ideal basic physical situation for developing models of phase change with size-dependent parameters: the symmetry reduces the problem to one dimension, while the geometry requires the introduction of curvature-driven effects. This means that the mathematical model can remain relatively simple while still elucidating the key physical features.

In the following sections we will closely examine the derivation of the Stefan problem via conservation equations and show that they lead to the omission of crucial terms, such as kinetic energy and momentum, in the Stefan condition. These terms have also been missed by researchers following other routes. The correct formulation will then be given. A main motivation behind this work is the large variety of models and specifically the Stefan condition being applied in recent models. We will demonstrate that what should be a seemingly straightforward derivation can in fact be problematic resulting in missing terms in many of the previously studied models. Even though we will often talk in the context of nanoscale melting the model will be based on the continuum assumption. The validity of this assumption in the context of phase change has been discussed in detail in [20, 41], with the basic conclusion that it holds down to around between 2 and 5 nm, depending on the material. The one-phase Stefan problem is a popular simplification, which also has its own significance in the context of ablation, where solid is ablated and subsequently carried away by an external flow [36]. With size-dependent parameters the standard form has been shown to lose energy [16, 42, 40]. We will also discuss the one-phase reduction for the current formulation. Finally, we briefly mention other forms of heat equation specifically designed for very small or very short time scales, which may subsequently be coupled to the new form of Stefan condition.

2 Governing equations for phase change

Phase change involves the energy-driven transformation of mass from one form into another. It is therefore natural to formulate the governing equations in terms of conservation laws and the thermodynamics of a two-phase mixture. The appropriate conservation equations are provided in [6]. In the specific context of phase change they are applied in [1, 39].

2.1 Formulation of the bulk equations

The standard governing equations for phase change may be derived from the following conservation laws:

[TABLE]

These represent conservation of mass, momentum, mechanical energy and total energy, see [6, §3.1, 3.2, 3.3, 10.1], [1, §2.3E]. The various parameters are density $\rho$, velocity $\mathbf{v}$ (where $v=|\mathbf{v}|$), pressure $p$, internal energy per unit mass $u$, and the conductive heat flux $\mathbf{q}$. Equation (1a) is a standard mass balance. Equation (1b) states that the change in momentum is balanced by the convection of momentum and the effect of pressure. Note $\otimes$ is the tensor product and $\mathbf{\underline{I}}$ the unit tensor. Since we will be working in three-dimensional space we may define $\mathbf{v}\otimes\mathbf{v}=\mathbf{v}\mathbf{v}^{T}$, which is a symmetric 3$\times$3 matrix. Equation (1b) therefore differs from the rest in that it is a vector conservation law rather than a scalar form. We will deal with this later. Equation (1c) states that the change in kinetic energy is balanced by the input of kinetic energy through the bulk flow, the work done by pressure, $p$, and the rate of reversible conversion to internal energy. Equation (1d) states that the rate of change of total energy depends on energy flow through convection, conduction and the work done by pressure.

In writing down these equations we have assumed that gravity and viscous effects are negligible. To simplify the discussion in the following we will work in the context of a solid to liquid transition but the model will, of course, be valid for other transitions. We assume that the solid and liquid phases are incompressible, having constant densities $\rho_{s}$ and $\rho_{l}$ that, in general, will not be equal. For the equations to hold throughout the material requires an assumption that there is a two-phase mixture or at least a small region, the interface, where there is a transition from solid to liquid. In applying the Rankine–Hugoniot condition below we impose a sharp-interface limit, i.e. the width of the interface tends to zero.

In regions of constant density mass conservation reduces to

[TABLE]

Consequently the final term of equation (1c) is zero. Subtracting equation (1c) from equation (1d) leads to the equation for conservation of internal energy

[TABLE]

The standard heat equation is often derived from this balance by expressing $u$ and $\mathbf{q}$ in terms of the temperature. These are discussed later.

2.2 Jump conditions at the interface

The density and internal energy are discontinuous across a sharp solid-liquid interface. Conservation of mass and energy links the jump in these two quantities to the motion of the interface. For a scalar conservation law written in divergence form,

[TABLE]

the Rankine–Hugoniot condition states that the jump in $F$ across a discontinuity moving with normal velocity $v_{n}$ is given by

[TABLE]

where $\hat{\mathbf{n}}$ is the unit vector normal to the discontinuity, $v_{n}$ represents the normal velocity of the interface, and $[\cdot]_{-}^{+}$ denotes the jump. In the following examples the normal vector points into the liquid so that the jump is measured from the solid into the liquid, i.e., $[\cdot]_{-}^{+}=[\cdot]_{s}^{l}$. In one-dimensional Cartesian problems the interface is generally denoted by $s(t)$, in one-dimensional cylindrical and radial problems by $R(t)$. In the following we will use either depending on the context.

Equation (5) may be applied immediately to the mass and energy balances to give

[TABLE]

where all variables are evaluated at the solid-liquid interface. Equation (6b) is an extension of the classical Stefan condition accounting for kinetic energy and pressure-volume work. This form of Stefan condition is derived through conservation laws in [1, §2.3E], [40].

The momentum balance is a vector relation, so we may not immediately apply the scalar form of Rankine–Hugoniot. A correct form may be obtained by dotting with the unit vector. If we define $\mathbf{\underline{R}}=\mathbf{v}\otimes\mathbf{v}$ then, provided $\hat{\mathbf{n}}$ is constant and $\mathbf{\underline{R}}$ symmetric, $(\nabla\cdot\mathbf{\underline{R}})\cdot\hat{\mathbf{n}}=\nabla\cdot(\mathbf{\underline{R}}\cdot\hat{\mathbf{n}})$. The appropriate scalar conservation form is then

[TABLE]

and Rankine–Hugoniot may now be applied

[TABLE]

The second relation is obtained after applying $\mathbf{\underline{I}}\cdot\hat{\mathbf{n}}=\hat{\mathbf{n}}$ and $[(\mathbf{v}\otimes\mathbf{v})\cdot\hat{\mathbf{n}}]\cdot\hat{\mathbf{n}}=(\mathbf{v}\cdot\hat{\mathbf{n}})^{2}$.

In [1] a number of cases of solidification and melting of a semi-infinite one-dimensional bar are described. Following the convention for one-dimensional Cartesian analyses, the interface velocity is written $v_{n}=s_{t}$ and $\mathbf{v}_{l}\cdot\hat{\mathbf{n}}=v_{l}$ is the liquid velocity in the $x$-direction. Here we choose an example where the solid is fixed such that $\mathbf{v}_{s}=\mathbf{0}$. The one-dimensional solid and liquid conductive fluxes are given by $q_{s}=\mathbf{q}_{s}\cdot\hat{\mathbf{n}}$ and $q_{l}=\mathbf{q}_{l}\cdot\hat{\mathbf{n}}$, respectively. The jump conditions (6a)–(6b) and (8) for mass, momentum and energy then reduce to

[TABLE]

and they hold at $x=s(t)$. Hence a form of the jump conditions given in [1, Ch. 2, eqs 25-27] is easily retrieved.

The momentum condition (9b) highlights a problem with this derivation form. It demonstrates that if the density varies between the phases then the resultant change in velocity causes a change in the pressure. This is fine when dealing with a flat interface. In the cases of nanowire and nanosphere melting examined in [5, 35, 40] the melting proceeds from the outside; following standard notation the interface velocity is denoted $v_{n}=R_{t}$ and $\mathbf{v}_{l}\cdot\hat{\mathbf{n}}=v_{l}$ is the (purely) radial velocity of the liquid, $\mathbf{v}_{s}\cdot\hat{\mathbf{n}}=0$ (since the central solid region is fixed). Consequently, the jump conditions which arise for both of these cases will be identical to equations (9a)–(9c) (with $s_{t}$ replaced by $R_{t}$). However, for curved interfaces it is well known that there is also a pressure jump caused by the surface tension $p_{s}-p_{l}=\sigma_{sl}\kappa$, where $\kappa\propto 1/R$ is the curvature and $\sigma_{sl}$ the surface tension between the solid and liquid phases: equation (9b) misses this contribution. The problem can be traced to the requirement that the conservation equations hold throughout the whole domain and so in fact represent a single material with varying physical properties, rather than two separate materials. The phase-change interface is the region where properties change and, for a sharp interface, this requires letting its thickness tend to zero, but this is the limit of a model which contains no interface and so there is no mechanism to account for surface tension.

To correctly account for surface tension in this model a Korteweg stress tensor of the form

[TABLE]

could be introduced into the conservation equations, where $\epsilon$ is related to surface tension and $\phi$ is the volume fraction of one of the phases. However, this then requires that a form be specified for $\phi$ which correctly reduces to the expected surface tension effect. A simpler ‘fix’ is to add a source term into the conservation equations which accounts for surface tension at the interface:

[TABLE]

If an outer surface exists then a second source term may be required (but this will not impact on the derivation of the Stefan condition below).

The Rankine–Hugoniot equations may also be derived by simply integrating the governing equations across a finite region which includes the jump and then letting the thickness of this region tend to zero. Noting that $\int_{-a}^{a}\delta(x)f(x)\,\mathrm{d}x=f(0)$ for any $a>0$ and that the velocity $v_{n}$ may be written $s_{t}$ in this configuration, the momentum and total energy equations now become

[TABLE]

where all quantities are defined at $r=s(t)$. The pressure jump now correctly accounts for surface tension, which also affects the energy balance. The addition of surface energy has no effect on condition (6a). The final equation may be rearranged after invoking equation (9a) to give a form of Stefan condition for a stationary solid phase

[TABLE]

again all quantities are defined at $r=s(t)$. Obviously this may be easily modified for different liquid and solid velocity combinations.

2.3 Internal Energy and Heat Flux

The Stefan problem is most commonly written in terms of temperature rather than internal energy and heat flux. The internal energy per unit mass (or specific internal energy) depends on the heat energy and work done against an external pressure. The heat energy is related to temperature by the heat capacity, $c$, which is defined as the ratio of heat energy entering a system to the temperature change induced by that energy $c=q/\Delta T$ (if $c$ varies with temperature then we define $q=\int c\,\mathrm{d}T$). It is standard to discuss heat capacities measured either at constant pressure, $c_{p}$, or constant volume, $c_{v}$. For solids and liquids, which are virtually incompressible, $c_{p}\simeq c_{v}$, so in the following we make no distinction between the two and simply write $c$. To convert from internal energy to temperature we follow [1, 6] and write

[TABLE]

The bulk phase change temperature $T_{m}$ is measured at a specified reference pressure, denoted $p_{\text{ref}}$. It is an arbitrary but standard choice to use the bulk value. This choice means that in most Stefan problems, when the temperature reaches the bulk phase change value (for the specified pressure) the enthalpy, $h_{s}=c_{s}(T_{s}-T_{m})$, of the solid is zero and the jump between the liquid and solid enthalpies is the latent heat. In [1] a reference internal energy is added to the problem $u_{\text{ref}}=p_{\text{ref}}/\rho_{s}$. However, given that it is a reference state it cancels throughout the calculations hence, for simplicity we will omit this term.

In [6] a differential form for the internal energy is provided, $\mathrm{d}u=c_{v}\mathrm{d}T$. For incompressible materials the volume is constant in each phase and we may employ this relation in deriving the heat equation. Note, in [1] a slightly different route is followed to obtain $\mathrm{d}u=c_{p}\mathrm{d}T$ (recall subscripts $v$, $p$ denote fixed volume and fixed pressure) despite the fact that pressure may vary in the phases. However, since the phases are incompressible and so the specific heats are equivalent the result is equivalent and we may write

[TABLE]

The standard form of heat flux is based on Fourier’s law

[TABLE]

where $k$ is the thermal conductivity. However, there are situations where this does not hold. These are described in more detail in §8.

The heat equations may be obtained from (3) by noting that $\mathrm{d}u=c\,\mathrm{d}T$ and using the appropriate law for the heat flux, in this case (16). Since the material properties are constant within each phase we obtain

[TABLE]

$s(t)$ is replaced by $R(t)$ in cylindrical and spherical problems. These equations may also be obtained by substituting for $u$, via equations (14a, 14b), into equation (3). This leads to terms of the form $\mathrm{d}(p/\rho)/\mathrm{d}t$. In [6, ch. 11.2] these are shown to be zero for incompressible phases. Using definitions (14a, 14b) the Stefan condition (13) may now be written, with the help of (12a), as

[TABLE]

where $T_{I}$ is the interface temperature (not necessarily equal to $T_{m}$). Equations (17)–(19) define the two phase Stefan problem when there is a density difference between the phases and a stationary solid.

In many studies kinetic energy is neglected in the Stefan condition, either because it is small compared to the latent heat or because density variation is not considered. However, it is worth pointing out that when applying a fixed-temperature boundary condition, $T_{s}(0,t)=T_{0}$, where $T_{0}$ is constant, then the initial solidification velocity is infinite and kinetic energy dominates. Also in the final stages of melting a cylinder or spherical material the melt rate tends to infinity. Hence there are times when kinetic energy is dominant but it is omitted from virtually all analyses, see [21]. Another important point to note is that equation (19) does not match the form provided in the textbook [1]. We discuss this now since their equations have been taken up by subsequent authors, see [21, 43, 44] for example. The issue may be traced to the solution of their Problem 24, §2.3. Their equation (48a) defines the internal energy in the solid at the interface as the reference value, this requires choosing $p_{s}(s)=p_{\text{ref}}$. To derive equation (48b) they neglect a term in the interface value of the internal energy of the liquid, to do this requires $p_{l}(s)=p_{\text{ref}}$, i.e. $p_{l}(s)=p_{s}(s)$. Yet their equation (33) defines the pressure jump $p_{l}(s)-p_{s}(s)\neq 0$. The goal of their equation (37) is to determine $p_{s}(s)-p_{\text{ref}}$ which then introduces an erroneous kinetic energy term, so affecting the kinetic energy in the final equation and thus all expressions of subsequent authors who take this as their base model.

3 Application to standard configurations

3.1 One-dimensional solidification of a semi-infinite bar

To illustrate the method we now focus on the standard problem of the solidification of a semi-infinite one-dimensional bar driven by cooling at $x=0$. The solid region occupies $x\in[0,s(t))$, the liquid occupies $x>s(t)$, where $s(0)=0$. In this scenario it makes sense to keep the solid material fixed and allow the fluid to move due to the density change, so that $\mathbf{v}_{s}=\mathbf{0}$ and then we may apply equations (17)–(19) without any amendment for a moving solid. The interface velocity is written $v_{n}=s_{t}$.

The liquid energy equation (18) requires knowledge of the velocity $v_{l}(x,t)=\mathbf{v}_{l}\cdot\hat{\mathbf{n}}$, where $\hat{\mathbf{n}}=\hat{\mathbf{x}}$. In general this may be obtained via the mass conservation equation (2), which here reduces to $\partial v_{l}/\partial x=0$ and hence $v_{l}(x,t)=v_{l}(t)$. Since the liquid velocity is independent of $x$ we only require its value at a single position to determine the value throughout the liquid. It may be calculated from the jump condition at $x=s(t)$, equation (6a),

[TABLE]

This in turn depends on the interface velocity, $s_{t}$, which is determined by the Stefan condition (19). The problem formulation is completed by specifying appropriate initial and boundary conditions. Here an initial temperature must be specified in the liquid, the solid does not exist at $t=0$ so requires no initial condition, while the phase change front satisfies, $s(0)=0$. In the case of melting we would specify an initial solid temperature and no liquid region. At $x=0$ any of the standard conditions may be applied. In the far-field the liquid temperature approaches the initial value. At the moving boundary the temperature depends on the situation. Most commonly the interface temperature is set to the bulk melt temperature $T_{I}=T_{m}$, so simplifying the Stefan condition. However, even with a flat interface $T_{I}$ may vary. For example when dealing with the solidification of supercooled materials [1, §2.4F], [42] there exists a nonlinear relation between the degree of supercooling and the front velocity [2, 19]

[TABLE]

where $C_{i}$ are constants. This equation determines the phase change temperature. A linearised version, valid for small supercooling, is more standard in mathematical studies [1, 19]

[TABLE]

The variation of $T_{I}$ with curvature is discussed in the following section.

With the above definitions the solidification of a semi-infinite one-dimensional bar may now be carried out by appropriate numerical or approximate methods.

3.2 Spherically symmetric nanoparticle melting

A popular problem in the literature is that of spherically symmetric nanoparticle melting. We now examine this since it permits a discussion of the model reduction in spherical co-ordinates and also introduces size-dependent parameters. The physical situation is depicted in Figure 1.

The nanoparticle is envisioned as being surrounded by a hot ambient environment with temperature $T_{a}$. As the solid melts, two interfaces evolve; these are the solid-liquid and the liquid-air. Due to the density difference between the solid and liquid, the liquid-air interface will grow in order to accommodate the change in volume that occurs during melting. The solid and liquid occupy the domains given by $0<r<R(t)$ and $R(t)<r<R_{b}(t)$. In this geometry, the velocity of the solid-liquid interface is denoted $v_{n}=\mathrm{d}R/\mathrm{d}t$ and the unit normal vector to the interface is $\hat{\mathbf{n}}=\hat{\mathbf{r}}$.

The heat equations are defined by the spherically symmetric versions of (17)–(18). As before an expression for the liquid velocity is obtained first by integrating the mass conservation equation (2) which indicates $v_{l}(r,t)=A(t)/r^{2}$. The function $A(t)$ is determined by the jump condition (6a) at $r=R(t)$

[TABLE]

Hence we find the relation between the liquid velocity and the rate at which the front moves

[TABLE]

This specifies the liquid velocity in equation (18).

At the outer boundary the liquid velocity is simply $v_{l}(R_{b},t)=\mathrm{d}R_{b}/\mathrm{d}t$. Noting that if we set $r=R_{b}$, $v_{l}(R_{b},t)=\mathrm{d}R_{b}/\mathrm{d}t$ in (24) and integrate then the following relation is obtained

[TABLE]

where we have imposed $R_{b}(0)=R(0)=R_{0}$. This determines the position of the outer boundary.

The interface velocity is defined by the Stefan condition (19).

The model of nanoparticle melting is completed by imposing initial and boundary conditions. Boundedness of the temperature at the origin requires

[TABLE]

At the solid-liquid interface, the solid and liquid temperatures are equal to the phase change temperature,

[TABLE]

At a curved interface (in the absence of supercooling) the temperature may be approximated by the Gibbs-Thomson relation

[TABLE]

with curvature $\kappa$. This may be viewed as a limit of the generalised Gibbs-Thomson which accounts for energy contributions due to changes in pressure and specific heat between the phases, see [1, 20]. If $\kappa$ is sufficiently small then the relation reduces to $T_{I}=T_{m}$.

Since this is a finite domain a thermal condition is also required at the outer boundary. Thermal energy is gained from the surrounding material, proportional to the difference in temperatures across the boundary,

[TABLE]

where $T_{a}$ is the temperature of the surrounding material and the constant of proportionality $h$ is the heat transfer coefficient. However energy is used up by the resistance to the change in particle size carried out on an environment which has constant pressure $p_{a}$

[TABLE]

Just inside the particle the change in energy is specified by

[TABLE]

(equation (31) is derived in the following section). Conservation of energy requires these three expressions to balance, leading to

[TABLE]

at $r=R_{b}(t)$. The Young-Laplace equation states that the pressure jump across the interface is given by $p_{l}-p_{a}=\sigma_{la}\kappa$, this removes all terms involving $v_{l}$. Then, despite the pressure work and creation of new surface the standard Newton cooling condition is retrieved

[TABLE]

Letting $h\rightarrow\infty$ leads to the standard fixed-temperature boundary condition $T_{l}(R_{b})=T_{a}$.

The initial radius of the nanoparticle is denoted as $R(0)=R_{b}(0)=R_{0}$. An initial solid temperature must be specified, in the literature this is often taken as the melt temperature $T_{s}(r,0)=T_{I}(R_{0})$. In this way melting starts at $t=0$ and there is no need to investigate an initial heating phase. Due to the lack of a liquid phase at $t=0$, an initial condition for the liquid temperature is not required. However, a small-time analysis can be used to calculate the liquid temperature for arbitrarily small times, which may then be applied in a numerical solution; see Appendix A.

4 Energy conservation

The derivation of the previous governing equations was based on conservation equations and consequently energy is naturally conserved, provided exchange at an outer boundary is correctly accounted for. However, a common simplification of the Stefan problem is the one-phase reduction. It has been demonstrated that in the presence of melting point depression the standard reduction does not conserve energy. We will now detail the energy conservation argument to show that energy loss or gain only occurs at the boundary. The argument will then inform the discussion of the one-phase reduction of §5.2 and will also clarify the expression for the boundary condition (31).

To verify that energy is conserved, we note that the governing equations have been derived from conservation laws, even the Stefan condition denotes conservation across the interface, so energy conservation should be both automatic and obvious. We may see this explicitly by considering the total energy in the spherically symmetric system. If we use the standard form of total energy conservation, equation (1d) (which does not account for interfaces), then the surface energy must be added at both the exterior and melt front. The total energy may then be written

[TABLE]

where the terms in the integrals obviously coincide with those in the total energy expression of equation (1d) (after taking into account the fact that the solid velocity is zero). This equation may also be derived, without explicitly adding the surface energy terms via equation (11c) after including a second source term to account for the outer boundary.

The rate of change of energy must balance that passing through the boundary. Taking the derivative of (34) with respect to time gives

[TABLE]

We may use equation (1d), with $\mathbf{v}_{s}=\mathbf{0},\mathbf{v}_{l}=v_{l}\hat{\mathbf{r}}$, to remove the time derivatives inside the integrals

[TABLE]

Now we express $v_{l}(R,t)$ in terms of $\mathrm{d}R/\mathrm{d}t$, via (23), and $v_{l}(R_{b},t)=\mathrm{d}R_{b}/\mathrm{d}t$. In spherical co-ordinates we may write $\kappa r^{2}=2r$ then after some algebra the above expression may be reduced to

[TABLE]

The terms in the square brackets evaluated at $r=R$ are simply the Stefan condition (13) and may therefore be set to zero leaving

[TABLE]

So, energy only leaves or enters the system at the outer boundary. The change is a result of the heat flux there, the work done against the ambient pressure to change the size of the particle and also the creation of new surface. This balance is the appropriate form required in the derivation of the Newton cooling boundary condition of the previous section.

5 Model reductions

5.1 Non-dimensionalisation

The relative importance of the terms in the Stefan problem can be assessed by casting the model into dimensionless form. We define

[TABLE]

where the time, length and temperature scales are $\tau,L,\Delta T$ respectively. These choices are not unique and should be based on the physical situation. A standard time-scale comes from writing the liquid heat equation in non-dimensional form. Taking the spherical version (and immediately dropping the hat notation) we obtain

[TABLE]

and consequently set $\tau=\rho_{l}c_{l}L^{2}/k_{l}=L^{2}/\alpha_{l}$, where $\alpha_{l}$ is the thermal diffusivity of the liquid, and the velocity scale $V=L/\tau$. An obvious length-scale would be the initial particle radius, $L=R_{0}$, while, since the phase change is driven by the external temperature, we choose $\Delta T=T_{a}-T_{m}$. With these scalings we obtain

[TABLE]

where $\alpha=\alpha_{s}/\alpha_{l}$. With the flux described by Fourier’s law the Stefan condition becomes

[TABLE]

where $\beta=L_{m}/(c_{l}\Delta T)$ is the Stefan number (or inverse Stefan number, both forms are common throughout the literature), $\gamma=(c_{l}-c_{s})\Delta T/L_{m}$, $\Gamma=2T_{m}\sigma_{sl}/(\rho_{s}L_{m}R_{0}\Delta T)$, $\rho=\rho_{s}/\rho_{l}$ and $\delta=R_{0}^{2}\left(1-\rho^{2}\right)/(L_{m}\tau^{2})$. If the Stefan number is large equation (42) indicates $R_{t}=\mathcal{O}(1/\beta)\ll 1$ and the phase change is slow. To focus on the evolution of the particle a new time variable $t^{*}=\hat{t}/\beta=t/(\beta\tau)$ may be used.

The dimensionless position of the outer boundary is given by

[TABLE]

and the velocity of the liquid is

[TABLE]

The boundedness condition at the origin remains unchanged. The Gibbs-Thomson equation becomes $T_{I}=-\Gamma/R(t)$ and

[TABLE]

The Newton condition may be written as

[TABLE]

where $Nu=hR_{0}/k_{l}$ is the Nusselt number. The initial conditions are

[TABLE]

The non-dimensional numbers provide information about the physical process as well as identifying dominant and negligible terms. In Table 1 typical parameter values are presented for tin and gold, these are taken from [20, 21, 43]. The value for $h_{\text{max}}$ for gold is calculated from $h_{\text{{max}}}=c_{s}\sqrt{\rho_{s}B/3}$, where $B=220\times 10^{9}$ Pa is the bulk modulus, see [43]. Four numbers are independent of the applied temperature difference: using the values for tin shows

[TABLE]

With a 10 nm sphere this indicates a heat conduction times-scale on the order $10^{-12}$ s. The term $\gamma\Gamma$ represents the effect of the melting point variation on the Stefan condition. Only in the very final stages of melting will this play a significant role: when $R={\cal O}(1)$ nm we expect this to cause a variation of the order 10%. This finding is in keeping with that reported in [39]. The kinetic energy term $\delta$ is order unity and so may play an important role, in keeping with the findings of [21]. This is especially true when combined with a fixed-temperature boundary condition. The Nusselt number $Nu$ is also of order unity, but this is when employing the maximum possible heat transfer coefficient. In general $h\ll 10^{9}$ W/(m2$\cdot$K) and we would expect the Nusselt number to be negligible, that is, it appears sufficient to apply a zero flux condition $\partial T_{l}/\partial r=0$ at the outer boundary. However, this is not the whole story, clearly to start the melting process some input of heat is required, i.e. a non-zero flux. So either the Nusselt number term must be retained or a rescaling applied at least for small times. Once melting begins the fact that the melt temperature of the solid decreases with decreasing radius as does the surface area to be melted means that the energy required to melt the front also decreases. This suggests a runaway process which, once started, will only increase in speed. This may be seen in all published results for the radius of a melting nanoparticle where $R_{t}\rightarrow-\infty$ as $R\rightarrow 0$, see [4, 5, 3, 20, 21, 35, 43, 52] for example. Similarly the analysis of [28], using the Maxwell-Cattaneo heat equation, shows ‘supersonic melting’, where the speed of the melt front is faster than the speed of heat propagation.

The remaining non-dimensional numbers depend on the temperature scale, when $\Delta T=10$ K:

[TABLE]

The melting time-scale is $\beta\tau$, which is of order $10^{-11}$ s, i.e. 10 picoseconds. The Stefan number $\beta$ may be considered large, indicating that a pseudo-steady solution will provide good accuracy. This is discussed in Section 5.3. The difference between the bulk melt temperature and the initial temperature is represented by $\Gamma$ which is here of order unity.

For gold we find similar values and hence the same conclusions: the main difference coming through a larger kinetic energy contribution due to a greater density difference between solid and liquid phases.

5.2 One-phase models

A common simplification in the study of Stefan problems is the one-phase reduction, whereby one of the phases is neglected. In the absence of supercooling or size-dependent properties this reduction may be achieved by simply setting the temperature in the neglected phase to the bulk melt temperature. The heat equation is then automatically satisfied while the corresponding conduction term in the Stefan condition is zero. However, when supercooling or size dependence occurs this reduction does not conserve energy, see [16, 42].

Taking into account both surface-tension effects and supercooling, an energy-conserving formulation based on negligible conductivity in the solid phase, i.e. $k_{s}\ll k_{l}$, is derived in [16]. They highlight a number of past papers where incorrect reductions are made, for example by naively setting $k_{s}=0$. Wu et al [52] support the assertions of [16] going on to state that if the initial temperature is different to the phase change temperature then the one-phase limit may only be derived under the assumption $k_{s}\ll k_{l}$. Due to the way heat is conducted through a material the thermal conductivity of a solid is nearly always greater than in its liquid phase: it is hard to conceive of a material where $k_{s}\ll k_{l}$. In [42] an asymptotic reduction is derived which is valid for the more physically sensible case $k_{s}>k_{l}$. Consequently the energy-conserving one-phase reduction in the limit $k_{s}\gg k_{l}$ is provided in [42].

In the context of nanoparticle melting the issue of one-phase reductions is discussed in more detail in [40] where it is shown that previous problems arise due to an inconsistency in assigning temperatures. Since most authors start from the standard governing equations they are often unaware that continuity of temperature at the interface has already been imposed. If the one-phase reduction is then made by setting $T_{s}$ to some value not equal, at all times, to the temperature imposed during the derivation then the inconsistency manifests itself in a loss of energy. It is concluded that energy-conserving forms may be written down immediately, without resorting to asymptotic reductions, provided consistency is maintained.

A one-phase model can be systematically derived from the non-dimensional system (41)–(42) after assuming $\rho c=(\rho_{s}c_{s})/(\rho_{l}c_{l})={\cal O}(1)$ and $k=k_{s}/k_{l}\gg 1$. Expanding in terms of $1/k\ll 1$ the solid temperature in a melting spherical nanoparticle is given to first order by

[TABLE]

where boundedness at $r=0$ has been applied. The leading-order term $T_{I}$ reflects the fact that the solid is a good conductor and hence rapidly transmits the boundary temperature. The Stefan condition states that the melt front is driven by the heat flux in the solid and liquid. With Fourier’s law the heat flux is proportional to the temperature gradient, which is zero to leading order in the solid since $T_{I}$ is a function of time. However, in equation (42) the solid temperature gradient is multiplied by $k$ meaning that the first-order term enters the leading-order balance and must be accounted for. The appropriate Stefan condition is then

[TABLE]

The final term in the square brackets derives from the energetic contribution of the solid to the melting process.

Given that the above analysis involves a varying solid temperature the term one-phase reduction may not be strictly appropriate. Taking only the leading order temperature $T_{s}=T_{I}(R(t))$ would be more appropriate for a true one-phase reduction, but this removes the $\Gamma$ term (i.e. the fourth term) from the Stefan condition which is at leading order, rather than a small correction. Further, the leading-order reduction for the temperature only satisfies the heat equation in the limit $k\rightarrow\infty$, so the form given by (48)–(49) is preferable.

The main message of [40] is that energy conservation with the one-phase reduction is largely an issue of consistency, with problems arising due to the use of accepted forms for the heat equations and Stefan condition without taking into account the assumptions applied during their derivation. With the current derivation energy conservation is not an issue. The Stefan condition (13) does not specify a form for the temperature hence the energy balance (38) holds, independent of the form of temperature. Depending on the choice of temperature an inconsistency can arise when trying to satisfy the heat equation. The choice of a constant temperature $T_{s}=T_{m}$ satisfies the heat equation, but not the boundary condition at $r=R$. The choice of $T_{s}$ specified by (48) is accurate to ${\cal O}(1/k^{2})$ while $T_{s}=T_{I}(R(t))$ only satisfies the heat equation in the limit $k\rightarrow\infty$. All consistent choices can conserve energy provided the Stefan condition is used in the form (13). A common error in a number of papers is the use of the form (19), which is based on setting the interface temperature $T_{s}=T_{l}=T_{I}(t)$: the subsequent choice $T_{s}=T_{m}$ then leads to loss of energy conservation.

5.3 Limit of large Stefan number

The Stefan number characterises the relative time scales of melting and thermal diffusion. The standard large Stefan number expansion, which corresponds to slow melting, requires a re-scaling of (dimensionless) time such that $t=\beta t^{*}$ and reduces the governing equations to

[TABLE]

where now $v_{l}=(1-\rho)(R^{2}/r^{2})R_{t}$ and the star notation has been dropped. The Stefan condition becomes

[TABLE]

The coefficient $1/\beta$ in the heat equations indicates that the system is now pseudo-steady, that is the heat equations at leading order assume their steady form indicating that the growth is slow compared to heat conduction. Time enters through the position of the moving boundary, specified by equation (51).

The leading-order temperatures are

[TABLE]

where $A(t)$, $B(t)$ are

[TABLE]

The kinetic energy term in the Stefan condition now has a factor $1/\beta^{2}$ indicating it may be neglected. However, it should be treated with care. In the standard problem where the boundary temperature is fixed $T_{l}(R,t)=T_{a}$ the initial front velocity $R_{t}$ is infinite. When the particle size tends to zero, with all standard boundary conditions $R_{t}\rightarrow\infty$. In both cases kinetic energy may be dominant, indicating the need for re-scaling. Even with a Newton boundary condition the reduction may not be trivial, as discussed below.

Provided $R_{t}\neq{\cal O}(\beta)$ substituting the temperatures (52) into the Stefan condition (51) yields a nonlinear differential equation for $R$. At leading order this is

[TABLE]

subject to the initial condition $R(0)=1$. Note, the limit $Nu\rightarrow 0$ indicates $A\rightarrow 0$ and hence $R_{t}\rightarrow 0$ as discussed in the previous section. With a fixed-temperature boundary condition the initial front velocity is infinite, in the examples studied later we use the highest value of heat transfer possible so leading to the possibility of high velocities. Here we note that $A(0)=-Nu(1+\Gamma)$, which from (54) implies that the initial speed of the interface is given by $R_{t}(0)=-\rho^{-1}Nu(1+\Gamma)(1-\gamma\Gamma)^{-1}$. The role of interface kinetic energy can then be estimated from the magnitude of

[TABLE]

The quadratic factors of the Nusselt number $Nu$ and $1+\Gamma$ that appear in the numerator of (55) imply that heat transfer and melting point depression, both of which are related to the size of the nanoparticle, can amplify the effect of kinetic energy. In the case of a 10 nm gold nanoparticle subject to a heating of $\Delta T=10$ K, we find that $\beta\simeq 39$, $\delta\simeq-54$, $Nu\simeq 0.46$, and $\Gamma\simeq 5.9$, leading to $K\simeq 0.15$. This value of $K$ is not particularly small, especially compared to $\beta^{-1}\simeq 0.03$, the small number used in the asymptotic expansion that led to (54). Consequently we anticipate inaccuracy in the expansion (which we will see in the numerical results). The problem could be dealt with via a more thorough analysis involving re-scaling near the boundary and then matching to an outer expansion. However, the issue really arises since we have chosen such a high value of heat transfer coefficient, in general this high velocity should not be an issue. Consequently we take a much simpler solution, by leaving kinetic energy in the Stefan condition and writing

[TABLE]

The kinetic energy term enters plays a significant role when $K$ is large, as $K$ reduces its effect will become less noticeable. The solutions for the temperature remain the same and are given by (52)–(53). Numerically solving equation (56) is relatively straightforward but requires dealing with the cubic equation for $\mathrm{d}R/\mathrm{d}t$ at each time step (rather than the simpler linear equation when kinetic energy is neglected).

6 Numerical solution

There are numerous published methods for the numerical solution of Stefan problems, such as finite difference, finite element, level set, enthaply, integral and heat balance methods [9, 30, 37, 38]. Our goal here is to demonstrate the form of solution and hence we will simply follow the finite difference scheme of [21].

For the spherical problem, the change of variable $\mathrm{u}(r,t)=rT_{l}(r,t)$ and $\mathrm{v}(r,t)=rT_{s}(r,t)$ is made, which effectively transforms the problem geometry to Cartesian coordinates. The heat equations (41) become

[TABLE]

The boundary conditions transform to

[TABLE]

and the Stefan condition (42) is

[TABLE]

We now immobilise the free boundaries by introducing space-like variables defined as $\zeta=r/R(t)$ and $\eta=(r-R(t))/(R_{b}(t)-R(t))$ in the heat equations for the solid and liquid, (57a) and (57b), respectively, which leads to

[TABLE]

These equations hold in fixed domains given by $0<\zeta<1$ and $0<\eta<1$. Furthermore, the liquid velocity in the new coordinate system is

[TABLE]

The boundary conditions become

[TABLE]

and the Stefan condition is

[TABLE]

Numerically solving the model requires an initial condition for the liquid temperature. This can be obtained from a small-time analysis of the governing equations, as detailed in Appendix A. Following the analysis of Appendix A we impose at some time $t_{*}\ll 1$, the following initial conditions:

[TABLE]

where $C$ is the solution of the cubic equation (76).

Spatial derivatives in the heat equations (60) are discretised using standard second-order centered finite differences. To avoid the use of ghost points, the first-order space derivatives with respect to $\eta$ and $\zeta$ in the Stefan condition (65) are approximated by second-order forwards and backwards difference formulae, respectively. The spatial coordinates $\zeta$ and $\eta$ are discretised into a grid of $I$ nodes. Semi-implicit time stepping of the heat equations is performed, whereby all of the terms involving $R$, $R_{b}$, and their derivatives are treated explicitly, and all terms involving the temperature are treated implicitly. Therefore, at each time step, we solve an $2I\times 2I$ linear system, given by the discretisation of (60), to determine the new temperature profiles. We then substitute the temperatures into the Stefan condition (65) and solve the resulting cubic equation for $\mathrm{d}R/\mathrm{d}t$ to obtain the radius $R$ at the next time step using an explicit update. Further details may be found in [21].

6.1 Transformed one-phase models

The governing equation for the liquid is given by (60b) with boundary and initial conditions (62) and (66), respectively. The problem is closed via the Stefan condition

[TABLE]

which is accurate to ${\cal O}(1/k^{2})$. If the effect of the solid is to be neglected completely then the final term in the square brackets disappears.

7 Results and discussion

The model proposed in Section 3.2 as well as its reductions will now be used to study the melting of spherical gold nanoparticles. Gold has been chosen because it is one of the most commonly used nanomaterials and, as can be seen from Table 1, it has a relatively large contrast between the density and thermal conductivity of each phase. The governing equations will be solved using the numerical scheme described in Section 6 and the parameter values provided in Table 1. This includes the maximum value of heat transfer coefficient which although unrealistically high leads to results in line with those previously published in the literature. We first focus on the role of density variation during melting. Then, the validity of the large Stefan number reduction is explored. Finally, we compare the two one-phase models ($k\gg 1$ and $k\to\infty$) to the full two-phase model.

7.1 Impact of density variation during melting

Font et al. [21] demonstrated that including density change in the model for the melting of a gold nanoparticle will significantly slow the melting process. This is primarily due to the addition of kinetic energy to the Stefan condition. Melting times more than double those of the constant density model were predicted. This study applied a fixed-temperature boundary condition, which results in an initial infinite melt rate, so it not so surprising that kinetic energy played an important role. Ribera and Myers [43] subsequently studied the melting of tin subject to a Newton condition (33). In this case the difference between variable-density and constant-density models was only a few percent. This change may be attributed to two factors, firstly the Newton condition does not lead to an infinite melt rate, secondly tin has a much smaller density contrast than gold. However, it clearly demonstrated that the effect of density change should not be as severe as suggested by Font et al. [21]. Unfortunately, as discussed in §2.3 the models used in both [21, 43] were derived from that of the textbook [1] which incorrectly represents the kinetic energy. As shown by equation (19) the term $(\rho_{s}/2)(1-\rho_{s}/\rho_{l})^{2}R_{t}^{3}$ used in the previous studies should be replaced with $(\rho_{s}/2)[1-(\rho_{s}/\rho_{l})^{2}]R_{t}^{3}$. By comparing the coefficients of $R_{t}^{3}$ in these terms (with the parameter values of Table 1) we find that kinetic energy in the present model is roughly 18 times stronger in the case of gold, and 70 times stronger in the case of tin, indicating that the role of density variation needs to be re-assessed.

Figure 2 displays the numerical solution for the variation of the radius with time with (solid lines) and without (dashed lines) density variation in the model. Initial radii of 10 nm and 100 nm were considered for various imposed temperature differences $\Delta T$. As stated, with density variation kinetic energy removes some of the energy available for melting which then acts to slow down the process. From the figure it is clear that for the 10 nm particle the differences in predictions are large, the melting time increases by 72%, 80%, and 140% for temperature differences of $\Delta T=1$ K, $10$ K, and $100$ K, respectively. With the 100 nm particles the effect is not so dramatic with increases of 14%, 15%, and 30% for $\Delta T=1$ K, $10$ K, and $100$ K.

Increasing the initial radius, up to 10 $\mu$m, the melting time settles at approximately 15% higher than the fixed-density model for all tested values of $\Delta T$. From this we see, in line with the conclusions of [21], that density plays an important role in the melting process well beyond the nanoscale and so should not be neglected in any practical study. Decreasing the heat transfer coefficient leads to reductions in the discrepancy: if we take $h=10^{-5}\,h_{\text{max}}$ then the difference in melt time is around 8% for large particles. Note, the effect of kinetic energy is reflected in the size of $\delta\propto 1/R^{2}$. For a 10 $\mu$m gold nanoparticle $\delta\approx 10^{-6}$, which is obviously negligible. With a small heat transfer coefficient there is no initial high velocity, so the differences cannot be due to kinetic energy. Density change also affects the position of the outer boundary, consequently we conclude that the difference in melt time for large particles must be due primarily to the fact that the fixed-density model gives errors in the boundary position, which in turn affects the domain over which the external heat must travel.

7.2 Validation of reduced models for large Stefan number

The reduction of the model in the limit of large Stefan number was described in Section 5.3. This led to a single differential equation for the position of the melt front, equation (54). It was then shown that with the present (maximal) value of heat transfer coefficient, the initial size of the kinetic energy term is $K\approx 0.15$. Thus, neglecting kinetic energy may lead to errors on the order of 15% from the very beginning of the melt process. Consequently, a second differential equation which accounts for kinetic energy, given by (56), was presented. As the Nusselt number decreases this equation reduces to the standard large Stefan number form.

In Figure 3 (a)–(b), we compare the numerical solution of the full model with the predictions of the large Stefan number reduction, with and without kinetic energy: symbols denote the numerical solution, solid lines the standard reduction, and dashed lines the reduction including the kinetic energy term. As before, we consider gold nanoparticles with initial radii of 10 nm and 100 nm. The $\Delta T$ values of 1 K, 10 K, and 100 K correspond to Stefan numbers $\beta=390$, $39$, and $3.9$, respectively. With such high values of $\beta=390,39$ we would expect almost exact agreement between all solutions, yet in the case of a 10 nm particle, the deviation between the full numerical and simple reduction is substantial, with errors in the melting times of 30%, 33%, and 53% for $\Delta T=$ 1 K, 10 K, and 100 K. Retaining kinetic energy results in much closer agreement with errors in the melting time of 4.1%, 3.6%, and 0.46%. With a 100 nm particle the value of $K$ is significantly lower (due to the reduction in $\delta$ by a factor of 100) and then both large Stefan reductions show excellent agreement with the numerical solution.

The importance of density change can be observed by comparing Figures 2 (b) and 3 (b). Neglecting the density change leads to melting-time errors between 14 and 30% in Figure 2 (b), whereas in Figure 3 (b) the results are almost identical. Figure 3 (b) demonstrates that for larger particles kinetic energy does not play such an important role and thus confirms the previous conclusion that errors, for sufficiently large particles, are a consequence of the fact that the constant-density model does not capture the growth of the outer boundary at $r=R_{b}(t)$. Practically this means that for sufficiently large material samples we could use an equation of the form (54), i.e. neglecting kinetic energy, but retaining the density difference in the term representing the motion of the outer boundary.

7.3 Comparison of one- and two-phase models

Two approaches for performing a one-phase reduction of the full two-phase model are discussed in Section 5.2. The first approach is based on systematically exploiting the large ratio of solid to liquid thermal conductivity, $k\gg 1$, and led to the modified Stefan condition given by (49). The second approach is based on simply neglecting the solid phase altogether and setting $T_{s}=T_{I}(R(t))$, which is technically only valid in the limit $k\to\infty$. The accuracy of the reduced one-phase models is now assessed through comparison with the full two-phase model. Figure 4 shows numerical simulations of the one- and two-phase models for various values of $\Delta T$ for a fixed value of $R_{0}=10$ nm. The ‘asymptotic’ one-phase model based on $k\gg 1$, shown as solid lines, is exceptionally accurate; in all cases the solutions perfectly coincide with those of the two-phase model (shown as symbols). The high degree of accuracy is remarkable given that the ratio of thermal conductivity for gold is $k\sim 3.0$, which is not particularly large. Even the one-phase model that treats the solid temperature as uniform (i.e. assumes $k\to\infty$, shown as dashed lines) is accurate, and overestimates the melting time by only 2%. This overestimate occurs because the conduction of heat from the bulk of the solid to the interface, which has a lower temperture than the bulk due to size effects, is neglected from this version of the one-phase model. If the initial radius of the nanoparticle is increased to 100 nm, then both one-phase models are in perfect agreement with the two-phase model.

8 Non-Fourier heat transfer

The results presented so far have been in non-dimensional form and so hide the fact that melting times for the 10 nm particle are on the order of picoseconds. This is comparable to the thermal relaxation time for many materials. Practically this means that Fourier’s law may no longer accurately describe the transfer of heat that occurs with phase change at very small length (and also time) scales. The breakdown of Fourier’s law, at both small length and time scales, has been predicted theoretically, via molecular dynamics, and observed experimentally, see [24, 55]. Consequently we will now briefly mention alternative heat conduction models.

There are many different models that aim to correctly describe heat transfer in situations where Fourier’s law is invalid, they can be categorized into micro-, meso- and macroscopic approaches. Micro- and mesoscopic models focus on each of the heat carriers or on their distribution, whereas the macroscopic models aim to describe the system in terms of quantities such as the temperature or the heat flux. The computational cost of approaches belonging to the first two categories, such as molecular dynamics [55], the Boltzman transport equation (BTE) [7] or the equation of phonon radiative transfer (EPRT) [34] is high and therefore only very small systems are usually considered. An approach that blends well with the PDE formulation used so far is a macroscopic model of the form

[TABLE]

This is called the Guyer–Krumhansl or GK equation. The first new parameter introduced here is $\tau_{R}$ which represents the thermal resistive relaxation time, i.e., the mean time between collisions among the heat carriers (phonons) which do not conserve their momentum; the second parameter $\ell$ is the phonon mean free path, i.e., the mean distance that a heat carrier travels between collisions. The first term in (68a) is said to add memory to the system and is dominant on time scales comparable to $\tau_{R}$, the second-order derivatives of the flux describe non-local effects that dominate heat transport in systems with length scales comparable to $\ell$. In the limit $\ell\rightarrow 0$ a simpler equation, the Maxwell–Cattaneo law, is retrieved [13, 51]. The Maxwell–Cattaneo equation gives rise to the hyperbolic or relativistic heat equation, which describes heat propagation with finite speed. Although the Guyer–Krumhansl equation was originally derived for situations with extremely low temperatures, recent studies have shown that it may be accurate at much higher temperatures [25]. It has been used to successfully predict the size dependence of the thermal conductivity in nanowires for temperatures ranging from 1 to 350 K [12, 11], a feature that cannot be captured using Fourier’s law. The effect of using the Guyer-Krumhansl equation on phase-change models has been studied in detail by Hennessy et al. [27]. The Maxwell–Cattaneo equation is used in the framework of melting nanoparticles in [28].

The use of a non-Fourier conduction law in a phase-change model can result in a temperature jump across the interface formed between the two phases [22, 47, 48]. When Fourier’s law is replaced by the Maxwell–Cattaneo equation, this temperature jump is needed to ensure that the speed of the interface is less than the finite speed at which thermal energy can be delivered to it. For example in the one-dimensional experimental setup described in [27] the temperature jump is given by

[TABLE]

where $\mathcal{V}=(\alpha/\tau_{R})^{1/2}$ is the speed of hyperbolic heat transport and $Q=-3\ell(\partial q/\partial x)$ is the flux of the flux [32].

In Fig. 5 we present typical solutions for the position of a solidification front predicted via the Guyer–Krumhansl equation using the finite-difference method described in [10]: we refer to this paper for more details on the scheme. In this case a one-phase situation is assumed, allowing us to write

[TABLE]

where $\gamma$ and $\eta$ are the non-dimensional relaxation time and the mean free path. For the simulation the boundary temperature was fixed $T(0,t)=-1$ and the temperature at the interface constant $T(s,t)=0$, initially $s=1$ with $T=q=0$. The position of the solid-liquid interface is plotted for typical values of $\gamma$ and $\eta$ as well as the classical case ($\gamma=\eta=0$). When the relaxation time is high, $\gamma=10$, it takes longer for heat to be transmitted from the boundary to the front and so solidification is slower than in the classical case. If non-local effects were neglected $\eta=0$ (which holds for the Maxwell–Catteneo model) this would be true for any relaxation time $\gamma>0$. However, from the curve with $\eta^{2}=10$ we see that the diffusion of heat flux can make the solidification occur more rapidly, despite the fact $\gamma>0$. The figure makes it clear that non-Fourier effects can play a significant role at small length and time scales, however for sufficiently large times (and hence sufficiently large length-scale $s(t)$) we expect non-Fourier effects to disappear, this is apparent as all solutions coincide for large times.

9 Conclusions

The governing equations describing the phase change of a material with distinct thermophysical properties in each phase were derived from conservation laws. This is a natural approach for a two-phase mixture, with a finite region over which properties change. As the width of this region approaches zero the sharp-interface limit, used in many Stefan problems, is retrieved. The one drawback to this approach is that being the limit of a mixture model there is no mechanism to account for surface tension. This could be dealt with formally by introducing an appropriate Korteweg stress tensor, here we simply added a source term at the interface. From the conservation laws the standard heat equations are obtained. Applying the Rankine–Hugoniot condition to the conservation laws leads to relations between the material and interface velocities, the pressure and energy jump. The latter expression provides the Stefan condition.

An important point to note is that the obtained Stefan condition differs to those found in the literature. Various expressions, with minor differences to that derived here, have been provided for the effective latent heat and the contribution of the difference between interface and bulk melt temperature. A particular error concerns the kinetic energy contribution which has subsequently propagated into subsequent studies. The Stefan condition derived here clarifies the various energetic contributions and is written in such a way that it may be readily incorporated into non-Fourier based studies as well.

To illustrate the use of the governing equations, the system was first applied to the standard problem of solidification of a semi-infinite bar and subsequently to the popular problem of a melting spherical nanoparticle. The nanoparticle case is particularly useful, since it permits the study of size dependent parameters in a simple one-dimensional setting. Using this vehicle, we were able to determine the significant effect of density change on melting. At small length-scales it is very large and the effect continues into the macroscale with errors of the order 10% observed. It was shown that for small particles kinetic energy plays an important role, while at larger scales it is primarily the error in the position of the outer boundary that leads to incorrect melting times. Overall, the difference between the fixed- and changing-density model depends on the relative density of the two phases, the heat transfer at the boundary, and the initial particle size.

Although the examples were primarily focused on the case of spherically symmetric nanoparticle melting the formulation is such that it can be readily applied to other practical geometries where size-dependent material properties play an important role, such as the case of cylindrical nanowire melting [18]. The evolution of micro- and macro-scopic systems can also be influenced by size-dependent material properties, a well-known example of which is the stabilisation of solid dendritic structures forming in a supercooled melt [15] by the change in melt temperature with curvature.

Finally we briefly discussed the problem of heat flow and phase change in situations where Fourier’s law no longer holds. This is becoming increasingly important as advances are made in nanotechnology. This field is currently wide open, with many opportunities for mathematical advancement into exciting new and practically interesting problems.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 707658. MC acknowledges that the research leading to these results has received funding from ‘la Caixa’ Foundation. TM acknowledges financial support from the Ministerio de Ciencia e Innovación Grant No. MTM2017-82317-P. The authors have been partially funded by the CERCA Programme of the Generalitat de Catalunya.

Appendix A Small-time solution

We may focus on the small-time behaviour by introducing an artificial small parameter $\epsilon\ll 1$ and then re-scaling time $t=\epsilon\bar{t}$. The radius $R$ is obviously close to its initial value, $R(t)=1-\epsilon\bar{R}(\bar{t})$. An expansion of (43) in terms of $\epsilon$ shows that $R_{b}=1+\epsilon(\rho-1)\bar{R}+O(\epsilon^{2})$. The Gibbs–Thomson relation can be expanded as $T_{m}=-\Gamma/R=-\Gamma-\epsilon\Gamma\bar{R}+O(\epsilon^{2})$. Upon writing $T_{s}(r,t)=\bar{T}_{s}(r,\bar{t})$, the leading-order problem for the temperature of the solid becomes $\partial\bar{T}_{s}/\partial\bar{t}=0$, with boundary and initial conditions given by $\partial\bar{T}_{s}/\partial r=0$ at $r=0$, $\bar{T}_{s}(1,\bar{t})=-\Gamma$, and $\bar{T}_{s}(r,0)=-\Gamma$. The appropriate solution is then $\bar{T}_{s}(r,t)=-\Gamma$ and is constant.

Since the liquid layer is initially very thin we must scale the radial co-ordinate in accordance with the variation of $R$, hence $r=1-\epsilon\bar{r}$. The liquid temperature is written as $T_{l}(r,t)=-\Gamma+\epsilon\bar{T}_{l}(\bar{r},\bar{t})$, which ensures that there is a balance of terms in the Newton condition. The leading-order problem for the liquid temperature is then

[TABLE]

subject to

[TABLE]

This has solution

[TABLE]

Since the leading order solid temperature is constant for small times it does not enter the Stefan condition and the melting is driven by the conduction of heat from the exterior through the liquid. The leading-order Stefan condition can be written as

[TABLE]

At small times the solution typically takes a power law form $\bar{R}=C\bar{t}^{\alpha}$. Here we find $\alpha=1$ and $C$ satisfies the cubic equation

[TABLE]

The small-time solution for the temperature and positions of the interfaces, in terms of the original dimensionless variables, can be written as

[TABLE]

where $C$ satisfies (76).