(a) Governing equations

Thermal diffusion in a single-phase, one-dimensional medium unaffected by distributed heat sources is governed by the equation

where X is distance in the direction of heat flow, θ the temperature above an arbitrary datum (e.g. the transition temperature), t time, U the medium velocity in the X direction and k, ρ, c _p are the medium thermal conductivity, density, and specific heat capacity respectively. When X is measured from a point fixed relative to a condensed phase i growing by phase transition from another phase j the above equation applied to each phase gives

for phases i and j respectively (see Fig. 1).

Fig. 1.

Transition in a one-dimensional domain.

Consider first the growing phase. To convert the variables to quantities of order unity, it is necessary to introduce new (normalized) variables x = X/X _c, ϕ = θ/θ _c, and τ = t/t _c, where the subscript c denotes a (positive) characteristic or reference quantity: it is not essential that the characteristic quantities be known at this point but merely that they exist. The first equation above may then be re-written as

where

is a Fourier number and k _c is a reference value of conductivity: the subscript i has been omitted for clarity. It is obvious that difficulties arising from variable properties in this phase are most closely associated with k because it is under the influence of the operator ∂/∂X. Fortunately, a simple transformation will remove the difficulty, at least from that point in the equation, for the important circumstances when k depends upon either temperature or position. If k = k(ϕ), the transformation with respect to ϕ.

will suffice. On the other hand, if k = k(x), the transformation with respect to x,

is required. In either case Equation (2) will then take the form

in which it is understood that ϕ and x are defined according to the dependency of k, and F ₀ is the corresponding Fourier number: note that the latter may or may not contain the variables ϕ or x. Thus the temperature distribution in a growing phase with temperature or spatially dependent properties may be determined from a “constant property” form given by Equation (3).

A similar procedure may now be followed for the diminishing phase, j, the only additional complication arising from the inclusion of bulk motion, i.e. advection. The velocity U _j in the positive X direction may be determined from the density difference ρ _i −ρ _j and the rate at which the excess (or deficit) volume appears. If E(t) is the location of the interface, then

and hence the normalized form of the second of Equations (1) is

where ξ(τ) = E(t)/X _c. The characteristic quantities inferred in this equation refer to phase j: the subscripts have been omitted except on the density. It might be expected that advection would be negligible when ρ _i ≈ρ _j and Equation (4) confirms this quantitatively. The magnitude of the error introduced by neglecting advection is evidently (ρ _i −ρ _j )/ρ _j , since this term alone gives the relative magnitude of the second and first terms on the right-hand side of the equation; from the definition of ϕ, ξ, x and τ, the derivatives may be taken to be of unit order.

(b) Boundary and interface conditions

A general treatment of the plane phase-change problem requires solutions of the system of equations (3) and (4). Before considering these solutions let us first consider the various boundary conditions which they may have to satisfy. Applying the first law of thermodynamics at an i−j interface we obtain

where L is the latent heat withdrawn in forming phase i from phase j, p the system pressure (assumed the same in each phase) and J _Q is the heat transfer rate to phase i, supplied externally. Kinetic and potential energy have been ignored. The second term on the left-hand side is the work required to introduce (or remove) the excess volume arising from the fact that ρ _i ≠ ρ _j . This “p−v” work is evidently negligible if p(ρ _i −ρ _j )/ρ _i ρ _j L ≪ 1. For an ice–water system this requires that p ≪ 10⁸ bars and since this pressure is rarely achieved the above equation may be simplified accordingly. It is interesting to note that during sublimation this term reduces to p/ρ _v L, where ρ _v is the vapour density, and may not be negligible.

The heat supply rate, J _Q may take several forms depending upon the nature of phase j and the environment. For two condensed phases

, in which case Equation (5), in normalized form, reduces to

where F _0i = κ _i tc _i /X _ci and Ste _i = c _pi θ _ci /L, will be called the Stefan number. Note the presence of characteristic quantities from both phases. From the definitions of ϕ _i , ϕ _j , x _i and x_j the derivatives on the right-hand side may be assumed to be of order unity and therefore the importance of the temperature distribution in the diminishing phase is seen to be governed by the magnitude of the non-dimensional group k _j X _ci θ _cj /k _i X _cj θ _ci . This is evidently small when

• a. the conductivity of the diminishing phase is much less than that of the forming phase, e.g. water freezing,

• b. the characteristic (or actual) depth of the forming phase is small compared with that of the diminishing phase,

• c. the departure from the transition temperature in the diminishing phase is small compared with the corresponding departure in the forming phase.

For transition in an ice-water system at the edge of a “semi-infinite” domain in which temperatures in the diminishing phase never greatly depart from the transition temperature it is reasonable to expect that k _j X _ci θ _cj /k _i X _cj θ _ci ≪ 1 is a condition which will frequently hold good, thus permitting the reduction of the above equation to the approximate form

each quantity referring to the forming phase at the interface.

At a moving interface, Equation (5) or its reduced form must be satisfied along with the fact that the temperature must be the transition temperature, here assumed fixed. (A transition temperature range may be accommodated by treating the problem as one of variable properties without latent heat (Reference Carslaw and JaegerCarslaw and Jaeger, 1959, p. 289; Reference MartynovMartynov, 1956, English translation p. 133).) If the interface is stationary, i.e. phase j is not being transformed to phase i at that point, the left-hand side of Equation (5) vanishes and the boundary condition assumes the explicit (Neumann) form

In this case the temperature at the interface need not be the transition temperature and in general it will be time-dependent. This dependency provides an alternative (Dirichlet) form of the boundary condition. The broadest single form arises when J _Q is given by

where h_Q is the heat transfer coefficient, σ the Stefan–Boltzmann constant, ϵ the surface emittance and the subscripts s and ∞ respectively refer to conditions at the interface and a substantial distance away from it: included in this expression are convection, emitted radiation and incident radiation F(t). Using normalized variables this may be combined with Equation (7) to give

where Φ_s = T _s/T _c, Nu is the Nusselt number h _Q X _ci /k _i and

a nondimensional “radiation number”: θ _c would now refer to the overall system, e.g. air temperature minus freezing temperature. If the environment is at a higher temperature than the interface the convective contribution is a heat addition and the positive sign in front of the Nusselt number is taken: the reverse situation requires the minus sign. The left-hand side of the equation would, of course, change sign if J_Q was directed into the interface from the left, instead of as shown in Figure 1. Equation (9) is a fairly complete statement of the heat transfer rate at a surface though it is occasionally extended to accommodate a mass transfer. For example, at the water–water vapour interface of a layer of water lying over melting ice, Equation (5) applied to the interface gives

where w denotes the water phase. The growth term can be re-written as Lh _M(C _∞−C _s), where h _M is the mass transfer coefficient and C _∞, C _S are the water-vapour concentrations at infinity and the interface respectively. This is equivalent to adding a term to Equation (7) and involving the vapour phase only as a boundary condition in the ice–water system. Clearly, this is permissible only if the velocity of the ice–water interface is much greater than that of the water–water vapour interface. With this understood, Equation (9) may be written in the extended approximate form

where Sh is the Sherwood number h _M X _ci /D, D is the mass diffusion coefficient, Q _R = DLC _∞/k _i θ _ci is a measure of the relative importance of latent and sensible heat effects and c _s = C _s/C _∞

In using the above equations it is assumed that Nu and Sh are available from a theoretical analysis of the convective system or from experimental data: that is, it is assumed that we know the relations

and

where Re is the Reynolds number, Sc = ν/D is the Schmidt number (the ratio of kinematic viscosity y to the mass diffusion coefficient D) and P _r = ν/κ the Prandtl number (the ratio of kinematic viscosity to thermal diffusivity κ). It is worth noting that for water vapour in air, Pr ≈ 0.7 and Sc ≈ 0.6 and hence for such a system the functions ϒ ₁, ϒ ₂ are almost identical. The resulting simplification will be particularly valuable for complex conditions, e.g. when the Reynolds number is time-dependent.

(c) Asymptotic solutions

The main advantage of writing the interface conditions in normalized form is to be found in the discussion of reasonable approximations and the range of their validity. Following Equation (5) it was demonstrated that “p − v”work would normally have a negligible effect on a moving interface and that conduction in a diminishing condensed phase could be ignored when

. Similarly, for a stationary interface, the importance of the convective transport of heat and mass is entirely dependent upon the magnitudes of Nu and Sh and Q _R in Equation (10). Since the dependent and independent variables in this equation are normalized it is reasonable to expect that the significance of convection can be measured by comparing Nu and the product Sh Q _R, with unity, neglecting the effects if they are very much less than this. The same argument applies to the role of G and f(τ) in determining the significance of radiation. The important feature of ordering terms in this way is that it is quantitative and enables approximations to be justified (or denied) in a simple, explicit manner a priori.

Another important consequence of ordering the moving interface Equation (5) is to be found in its special and reduced forms Equations (5a) and (6). By definition, the order of magnitude of the term in square brackets on the right-hand side of Equation (5a) and of its reduced form, the partial derivative in Equation (6), is unity. Now from either of these equations it is clear that each side of the expression is essential to the description of conditions at the interface and therefore both sides must be retained and must be equal in magnitude. Consistent with the normalization process, this implies that each term is of unit order and therefore

or

For a constant-property medium, this provides a convenient definition for t _c (or X _c) given by

from which it follows that Ste = 1/F ₀ and therefore the Equation (3), for conduction in the forming phase, becomes

The above equation, and the corresponding form of Equation (4), each reveal the essential role of the Stefan number in the determination of the temperature distribution. When sensible and latent heat effects are about equal, i.e. Ste = 0(1), the complete conduction equation must be solved. However, in many situations common to ice–water systems latent heat effects dominate, i.e. Ste ≪ 1, which suggests solutions for the temperature in the form

Substituting this into the conduction equation above leads to the infinite set of equations

and

for which the solutions are

where the prime denotes differentiation with respect to τ.

In this paper discussion will be limited to the asymptotic situation Ste → 0. This amounts to the neglect of all the terms in the above expansion for temperature except ϕ ₀(x, τ) given by Equation (11). Since ϕ ₁, ϕ ₂, etc. exhibit the curvature in the temperature profile it follows that when Ste = 0 the profile must be linear. Prima facie, a straight-line profile appears almost unreasonable, but closer examination shows that it is completely consistent with the implication that latent heat effects dominate the sensible heat effects. A very large latent heat implies a layer whose growth is so slow that it can be treated as a slab of virtually constant thickness. A negligible specific heat, on the other hand, suppresses the transient aspects of conduction within a slab. Taken together, these physical characteristics agree exactly with the mathematical consequences of Ste = 0. The error in the solution when Ste is small but non-zero can be estimated from the expansion for ϕ ₀(x, τ), which indicates that it will be of the order of the largest neglected term, i.e. Ste (since ϕ ₁ = 0(1)). It will be demonstrated later that typical errors confirm this expectation. Strictly, the above discussion applies to a constant-property medium only but the arguments for a variable-property medium are identical and lead to the same conclusions. The only difference in the solution would be in the definitions of ϕ or x used in developing Equation (3).

The use of asymptotic solutions subject to boundary conditions typified by Equations (5a) and (10) covers a very wide range of problems. Our discussion will be limited to those which occur at the edge of a “semi-infinite” domain and in the absence of mass transfer. This will remove very little of the generality of the discussion (since the convective heat- and mass-transfer systems are analogous) without masking the essential simplicity of the approximate solutions. For this type of discussion we may adopt the approximate form (Equation (6)) for the moving interface and since this implies the neglect of conduction in the diminishing phase, Equation (4) may now be disregarded. In short, the problem is reduced to the solution of Equation (3) subject to the interface condition, Equation (5), and boundary conditions typified by Equation (9).

The asymptotic solution of Equation (3), given by Equation (11), may be substituted into the interface condition giving

or

for a constant-property medium. From this it can be seen that arbitrary boundary conditions, as they appear though a(τ), may be handled without much difficulty. In general, the right-hand side of this equation will depend upon ξ and τ and hence we may take dξ/dτ = s(ξ, τ), say, as the most general form of the interface equation for asymptotic conditions. In the most general circumstances s(ξ, τ) will be a non-linear implicit function and therefore the equation cannot be integrated directly, in which case a numerical or graphical solution, e.g. the method of isoclines, must be sought. In many other circumstances the form of s(ξ, τ) permits integration with comparative ease and leads to simple closed-form or series solutions. In this paper the latter circumstances will be used to assess the validity of the asymptotic solutions.