PROCEEDINGS OF THE OKLAHOMA ACADEMY OF SCIENCE Volume 67—1987

{Page 83}

Limitation in the Apparent Heat Capacity Formulation for Heat Transfer With Phase Change

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

INTRODUCTION

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 two-phase Stefan problem for a one-dimensional 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 semi-infinite 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.)

FORMULATION

Consider the energy equation for a two-phase system of a pure substance:

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

Eliminating f₂ by means of eq 2 and rearranging permits eq 1 to be written as

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}

in which k is the thermal conductivity tensor and T is temperature. Eq 3 can be rearranged as

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₁h₁ - r₂h₂), 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, single-phase 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

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

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

Alternatively the normal distribution function sketched in Fig. 2(b) can be used as a smooth variation [Hashemi and Sliepcevich (1)]

in which Î is so chosen that erf (Î D T_t)= 1.0 -l, where DT_t is one-half 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 one-dimensional transient-state 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.

APPLICATION

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 one-dimensional medium (x-direction) containing a material in a thermodynamic state represented by phase 1, the temperature of which is uniform. A transient-state 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

in which

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,

and

For comparison with an analytical solution consider the initial and boundary conditions as

To achieve the optimum accuracy in the numerical solution of eq 10 apply the Boltzmann similarity transformation

and the transformation to map the semi-infinite domain to a unit-size domain

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

in which (k / rc)₀ is the maximum value of the heat diffusivity in the range (T_i T_s) . In addition a normalized temperature is defined by

As a result, eqs 10, 16, 17, and 18 are transformed to

In eq 23 y is the diffusivity ratio given by

The solution to eq 23 is accomplished by separating it into two first-order, ordinary differential equations:

and

The numerical solution has been obtained by solving eqs 27 and 30 simultaneously by a variable-step Runge-Kutta-Fehlberg 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.

RESULTS and Discussion

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⁶ and 2.5 X 10⁶ J/m³•K, respectively. The transition temperature is 0 °C and the volumetric latent heat is 1 x

{Page 87}

10⁸ J/m³. 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. 5-7. 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

SYMBOLISM

REFERENCES

1.   H.T. Hashemi and C.M. Sliepcevich, Chem. Eng. Progr., Symp. Ser., 63, No. 79, 34-41 (1967).

2.   F. Civan and C.M. Sliepcevich, Int. J. Heat Mass Transfer 27: 1428-1430 (1984).

3.   F. Civan and C.M. Sliepcevich, Water Resources Res. 21: 407-410 (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: 615-621 (1978).

6.   E. Fehlberg, Low-Order Classical Runge-Kutta Formulas with Stepsize Control and Their Application to Some Heat Transfer Problems, NASA TR R-315, NASA, Huntsville, AL, July 1969.