
{Page 83}
Faruk Civan
School of Petroleum and Geological Engineering, University of Oklahoma, Norman, OK 73019
C.M. Sliepcevich
School of Chemical Engineering and Materials Science, University of Oklahoma, Norman, OK 73019
Of the various approaches for modeling problems in heat conduction with phase change (the classical Stefan problem) the apparent heat capacity formulation still appears to be one of the preferred methods, judging from the recent literature. The principal advantage of this approach is that temperature is the primary dependent variable that derives directly from the solution. However, it is well known that this formulation suffers from a singularity problem for a phase change that occurs at a fixed temperature. This difficulty can be circumvented by assuming that the phase change occurs over a small temperature interval [Hashemi and Sliepcevich (1)]. Another alternative would be to utilize an enthalpy formulation [Civan and Sliepcevich (2,3)]. However, the purpose of this presentation is to quantify the error introduced in the apparent heat capacity formulation as a consequence of assuming a finite temperature interval for the phase change.
To investigate the effect of this phasechange temperature interval the twophase Stefan problem for a onedimensional system undergoing a phase transition at a fixed temperature will be solved for the case of constant physical properties in each phase. The results will be compared with the analytical solution given by Neumann in Carslaw and Jaeger (4). For this purpose, a simple and accurate computational method developed in previous work by the authors for semiinfinite media will be used. Normalized variables, along with the Boltzmann similarity variable and coordinate mapping transformation eliminate computational difficulties and problems that frequently arise with conventional techniques, such as the finite difference method [Goodrich (5)]. For the present purpose, a typical example of a pure incompressible substance undergoing cooling and freezing (or conversely, heating and melting) will be solved. (This same methodology can be applied to other processes involving discontinuities, such as shock phenomena and interfacial mass transport.)
Consider the energy equation for a twophase system of a pure substance:
(¶ / ¶t) (r_{1}h_{1}f_{1} + r_{2}h_{2}f_{2}) + Ñ•[r_{1}h_{1}V_{1} + r_{2}h_{2}V_{2} + (q_{1}f_{1} + q_{2}f_{2})]  Dp / Dt = 0  (1) 
in which subscripts distinguish the phases, r and h are the density and enthalpy, and f is the volume fraction (volume of one phase divided by the volume occupied by both phases), respectively. V and q are the volumetric and heat fluxes associated with each of the phases, respectively. D/Dt represents a substantive derivative (¶ / ¶t) + V • Ñ); r is pressure and t is time. By definition
f_{1} + f_{2} = 1  (2) 
Eliminating f_{2} by means of eq 2 and rearranging permits eq 1 to be written as
(¶ / ¶t) [r_{2}h_{2} + (r_{1}h_{1}  r_{2}h_{2})f_{1}] + Ñ•{(r_{1}h_{1}V_{1} + r_{2}h_{2}V_{2}) + [q_{2} + (q_{1}  q_{2})f_{1} ]}  Dp / Dt = 0  (3) 
Since the phases are at thermal equilibrium during the phase transition, they are simultaneously at the same temperature. On the assumption that all heat transferred is by conduction, then Fourier's equation applies:
{Page 84}
q = k•ÑT  (4) 
in which k is the thermal conductivity tensor and T is temperature. Eq 3 can be rearranged as
{¶(r_{2}h_{2}) / ¶T + [¶(r_{1}h_{1}) / ¶T  ¶(r_{2}h_{2}) / ¶T]f_{1} + (r_{1}h_{1}  r_{2}h_{2}) (df_{1} / dT)}¶T / ¶t + Ñ•{(r_{1}h_{1}V_{1} + r_{2}h_{2}V_{2})  [k_{2} + (k_{1}  k_{2})f_{1}]•ÑT}  Dp / Dt = 0  (5) 
The energy equation in the form of eq 5 is called an apparent heat capacity formulation in which the quantities enclosed in the braces preceding ¶t / ¶t represent an "apparent" volumetric heat capacity. Since the volumetric latent heat effect of phase change, (r_{1}h_{1}  r_{2}h_{2}), is included in the apparent volumetric heat capacity definition, the apparent heat capacity formulation allows for a continuous treatment of a system involving phase transfer. In contrast the Stefan formulation requires solutions of separate, singlephase equations for each of the phases; these solutions are later matched along the phase change interface boundary for compatibility and consistency.
If the phase transition takes place instantaneously at a fixed temperature, then a mathematical function such as
f_{1} = U(T  T_{t})  (6) 
is representative of the volumetric fraction of the phase 1 as shown in Fig. 1(a). U is a step function whose value is zero when T<T_{t} but equals one otherwise. Its derivative, i.e., the variation of the phase 1 fraction with temperature sketched in Fig. 1 (b), is
df_{1} / dT = d(T  T_{t})  (7) 
in which d (T  T_{t} ) is the Dirac delta function whose value is infinity at the transition temperature, T_{t} , but is zero at all other temperatures. To alleviate this singularity the Dirac delta function can be approximated by a uniform distribution over a finite temperature interval as sketched in Fig. 2(a). Thus eq 7 is replaced by
df_{1} / dT = 1 / (2DT_{t})  (8) 
Alternatively the normal distribution function sketched in Fig. 2(b) can be used as a smooth variation [Hashemi and Sliepcevich (1)]
df_{1} / dT = (Îp^{ 1/2}) exp[ Î^{2} ( 1 / (2DT_{t})T  T_{t})^{2}]  (9) 
in which Î is so chosen that erf (Î D T_{t})= 1.0 l, where DT_{t} is onehalf of the assumed phase change interval and l is a sufficiently small positive number. Consequently, the integrals of eqs 8 and 9 yield the linear and the error function approximations, respectively, for the phase 1 fraction as sketched in Figs. 3(a) and 3(b). With conventional finite difference and finite element methods the phase 1 fraction
{Page 85}
derived from eq 8 or 9 by integration should be used to avoid the numerical instabilities arising from the jump in the values of the volumetric fraction of phase 1 from zero to one as shown in Fig. 1(a). However, as demonstrated in the following application for onedimensional transientstate heat transfer, the jump property represented by eq 6 does not create any problems since the numerical solution is carried out in terms of the Boltzmann variable, which combines the time and the distance variables into one variable.
To demonstrate the effect of the magnitude of the small temperature interval introduced by means of the approximations given by eq 8 or 9, consider an infinitely long onedimensional medium (xdirection) containing a material in a thermodynamic state represented by phase 1, the temperature of which is uniform. A transientstate heat transfer process is initiated by contacting the leading surface of this material with a sink of infinite extent whose temperature is such that it results in the formation of phase 2 as shown in Fig. 4. In other words, the material undergoes a transition from phase 1 to phase 2 with the interface progressing gradually with time from the leading surface to infinite distance.
For simplification we assume that the
{Page 86}
pressure is constant, the phases are isotropic and homogeneous, the phase densities are equal, and there is no motion in either of the phases. Accordingly, eq 5 reduces to
rc (¶T / ¶t) = (¶ / ¶x)(k¶T / ¶x  (10) 
in which
r = r_{1} = r_{2}  (11) 
c = c_{2}+(c_{1}  c_{2})f_{1} + Ldf_{1} / dT  (12) 
k = k_{2} + (k_{1}  k_{2})f_{1}  (13) 
In eq 12 c denotes the specific heat capacity and L the specific latent heat for change from phase 2 to phase 1, which are determined according to the specific enthalpies of phases 1 and 2, respectively,
h_{1} = L +T_{t}òT c_{1}dT for T > T_{t}  (14) 
and
h_{2} = T_{t}òT c_{2}dT for T £ T_{t}  (15) 
For comparison with an analytical solution consider the initial and boundary conditions as
T = T_{i} , t = 0,0£x <¥  (16) 
T = T_{s} , t >0 , x = 0  (17) 
T = T_{i} , t >0 , x®¥  (18) 
To achieve the optimum accuracy in the numerical solution of eq 10 apply the Boltzmann similarity transformation
y = x / t^{1/2}  (19) 
and the transformation to map the semiinfinite domain to a unitsize domain
z = 1  exp(y / b)  (20) 
as demonstrated in the earlier publications by Civan and Sliepcevich (2,3). In eq 20 b is an arbitrary scaling factor. In this study it is assumed that
b = 2^{1/2}(k / rc)_{01/2}  (21) 
in which (k / rc)_{0} is the maximum value of the heat diffusivity in the range (T_{i} T_{s}) . In addition a normalized temperature is defined by
u = T_{i}òT kdT / T_{i}òT _{s}kdT  (22) 
As a result, eqs 10, 16, 17, and 18 are transformed to
ln(1 z) (du / dz) = y(d / dz)[(du / dz)]  (23) 
u = 1 , z = 0  (24) 
u = 0 , z = 1  (25) 
In eq 23 y is the diffusivity ratio given by
y = k / rc) / (k / yc)_{0}  (26) 
The solution to eq 23 is accomplished by separating it into two firstorder, ordinary differential equations:
du / dz = v / (1  z)  (27) 
u = 1 , z = 0  (28) 
u = 0 , z = 1  (29) 
and
dv / dz = (1 / y) ln(1  z) (du / dz)  (30) 
The numerical solution has been obtained by solving eqs 27 and 30 simultaneously by a variablestep RungeKuttaFehlberg four (five) method (Fehlberg (6)] using double precision and by requiring that truncation errors are less than 10^{12} in using an IBM 3081 computer at the University of Oklahoma.
For numerical calculations the data from Goodrich (5) has been considered in which the thermal conductivities of the frozen and unfrozen phases are given as 2.25 and 1.75 W/m • K, the volumetric heat capacities as 1.5 x 10^{6} and 2.5 X 10^{6} J/m^{3}•K, respectively. The transition temperature is 0 °C and the volumetric latent heat is 1 x
{Page 87}
10^{8} J/m^{3}. Computations were carried out for cases in which the initial and leading surface temperatures, T_{i} and T_{s}, were considered to be 2 and 2 °C, 2 and 10 °C, and 10 and 10 °C, respectively.
Calculations were made by approximating the Dirac delta function by a normal distribution and by a uniform distribution over the temperature intervals of ±0.01, 0.1, 0.5, 1.0, and 1.5 °C. Since the results with the uniform distributions, eq 8, were about twice as greatly in error as with the normal distribution, eq 9, only the results with the normal distribution are shown in Figs. 57. Furthermore, these figures are expressed in terms of the normalized variables in order to magnify the errors.
It is evident that the accuracy of the apparent heat capacity method is sensitive to the magnitude of the assumed temperature interval that is arbitrarily selected to approximate the Dirac delta function or to distribute the latent heat. In fact for DT_{t} <0.1 °C, values were essentially identical to those obtained from Neumann's analytical solution as shown in Figs. 5 and 6; even for DT_{t}< 1.5 °C (Fig. 7) the differences between the numerical and the analytical results were of the order of 10^{3}.
It is evident that the error introduced by approximating a singularity (Dirac delta function) by a gradual change over a small temperature interval is more pronounced in the unfrozen region, which is effectively of infinite thickness, than in the thin, frozen region near the leading surface; the latter conforms closely to the exact solution. Another observation is that the effect of distributing the latent heat over a temperature interval diminishes as the ratio of the
{Page 88}
temperature interval to the overall temperature variation of the system becomes smaller. Since the apparent heat capacity method approaches the exact analytical solution as the assumed temperature interval approaches zero, it is evident that the apparent heat capacity method can produce accurate solutions.
It appears from these numerical solutions that the results obtained by approximating the phase change at a fixed temperature by a gradual change over a small temperature interval should be acceptable if
(2DT_{t}) /  (T_{i}  T_{s}  <0.1  (31) 
Italic Symbols
b 
scaling factor 
c 
heat capacity (J/kg•K) 
h 
specific enthalpy (J/kg) 
k 
thermal conductivity (W/m•K)(scalar) 
L 
latent heat (J/kg) 
p 
pressure (N/m^{2}) 
t 
time(s) 
T 
temperature (°C) 
u 
normalized temperature 
v 
a function resulting from separating the secondorder ordinary differential equation into two firstorder ordinary differential equations 
x 
distance (m) 
y 
similarity variable 
z 
transform variable 
Bold Symbols
k 
thermal conductivity (W/m•K)(tensor) 
q 
heat flux(J/M^{2}•s) (vector) 
V 
volume flux (m^{3}/m^{2}•) (vector) 
Greek Symbols
r 
density (kg/m^{3}) 
l 
a sufficiently small number 
Î 
a parameter 
f 
volume fraction (m^{3}/m^{3}) 
y 
diffusivity ratio 
Mathematical Symbols
d 
ordinary differentiation operator 
¶ 
partial differential operator 
d 
Dirac delta function 
D/Dt 
substantive derivative (¶ / ¶t + V•Ñ) 
exp 
exponential base, e 
U 
unit step function 
D 
finite difference operator 
Ñ 
divergence operator (vector) 
¥ 
infinity 
Subscripts
i 
initial 
a 
leading surface 
t 
phase trasition point 
1,2 
phases of pure substance 
1. H.T. Hashemi and C.M. Sliepcevich, Chem. Eng. Progr., Symp. Ser., 63, No. 79, 3441 (1967).
2. F. Civan and C.M. Sliepcevich, Int. J. Heat Mass Transfer 27: 14281430 (1984).
3. F. Civan and C.M. Sliepcevich, Water Resources Res. 21: 407410 (1985).
4. H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, Clarendon, Oxford, 1969.
5. L.E. Goodrich, Int. J. Heat Mass Transfer 21: 615621 (1978).
6. E. Fehlberg, LowOrder Classical RungeKutta Formulas with Stepsize Control and Their Application to Some Heat Transfer Problems, NASA TR R315, NASA, Huntsville, AL, July 1969.