E. M.

Sparrow
S. V. Patankar Analwsis of Melting in the Presence
S. Ftamadhyani
Department of Mechanical Engineering
University ol Minnesota, Minneapolis, Minn. of Natural Confection in the Melt
Region
An analysis of multidimensional melting is performed which takes account of natural
convection induced by temperature differences in the liquid melt. Consideration is given
to the melt region created by a heated vertical tube embedded in a solid which is at its fu-
sion temperature. Solutions were obtained by an implicit finite-difference scheme tai-
lored to take account of the movement of the liquid-solid interface as melting progresses.
The results differed decisively from those corresponding to a conventional pure-conduc-
tion model of the melting problem. The calculated heat transfer rate at the tube wall de-
creased at early times and attained a minimum, then increased and achieved a maximum,
and subsequently decreased. This is in contrast to the pure conduction solution whereby
the heat transfer rate decreases monotonically with time. The thickness of the melt region
was found to vary along the length of the tube, with the greatest thickness near the top.
This contrasts with the uniform thickness predicted by the conduction solution. These
findings indicate that natural convection effects, although unaccounted for in most phase
change analyses, are of importance and have to be considered.

Introduction literature. Therefore, many of the available conduction-based phase

change solutions do not fully reflect physical reality.
Melting and freezing processes are encountered in a wide range of The accounting of natural convection exposes differences between
technologies, encompassing such diverse applications as casting of freezing and melting problems whose conduction-based solutions are
metals and glass, freeze drying of foodstuffs, and storage of thermal interchangeable. In this connection, attention may be directed to the
energy. The number and diversity of melting/freezing applications counterpart problems of freezing a liquid which is at its fusion tem-
has motivated an extensive literature. As will be discussed shortly, perature and of melting a solid that is also at the fusion temperature.
an examination of this literature shows that, in the main, solid/liquid In the freezing problem, the temperature of the unfrozen liquid re-
phase change problems are analyzed as heat conduction problems. mains at the fusion value, thereby precluding the possibility of natural
Conduction can be accepted as the sole transport mechanism in sit- convection. On the other hand, the liquid produced in the melting
uations where the temperature of the liquid is uniform and equal to problem is necessarily nonuniform in temperature, so that natural
the fusion temperature.1 However, when the temperature of the liquid convection motions may occur. Since the most commonly analyzed
is nonuniform, consideration has to be given to the possible role of melting/freezing problems encompass media that are initially at the
natural convection. fusion temperature, it is worthy of note that the solutions of the
The presence of temperature variations in the liquid will, in most melting and freezing problems are not interchangeable due to natural
cases, activate natural convection motions. Under these conditions, convection.
the heat conduction problem in the solid has to be solved simulta- It is also relevant to note that there are some instances where the
neously with the natural convection flow and heat transfer problem presence of temperature variations in the liquid does not lead to
in the liquid. Historically, the possible existence of natural convection natural convection motions. One such instance is stable stratification,
in the liquid phase has largely gone unnoticed in the heat transfer wherein liquid of lesser density is situated above liquid of greater
density. Even when heavier fluid is situated above lighter fluid, nat-
ural convection does not occur until an instability threshold is reached.
In many phase change problems, natural convection is relegated to
an ineffectual role due to the high resistance offered to fluid motion
1
by passages of small dimensions. For instance, such effects are oper-
Excluding situations where radiative transfer is involved. ative in low-permeability fibrous or granular materials.
Contributed by the Heat Transfer Division for publication in the JOURNAL
OF HEAT TRANSFER. Manuscript received by the Heat Transfer Division The foregoing general discussion of the role of natural convection
June 27,1977. in solid/liquid phase change problems is intended to set the stage for

S20 / VOL 99, NOVEMBER 1977 Transactions of the ASME

Copyright © 1977 by ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 07/24/2013 Terms of Use: http://asme.org/terms

 the analysis and the results to be presented in this paper. The moti-
vating application for the study is the storage of thermal energy by
the melting of a salt. Qualitative laboratory experiments performed
in connection with the storage work revealed a decisive influence of
natural convection. In view of those findings, it appeared mandatory
to take account of natural convection in the analysis and solution of
the melting problem, and that dictum was followed here.
The physical situation that was analyzed is pictured schematically
in Fig. 1. A vertical tube of height H and outer radius r 0 is embedded
in a solid which is at its fusion temperature Tsat. Initially, no liquid
is present, and the solid and the tube are in direct contact. At time t
= 0, the temperature of the outer surface of the tube is raised to a value
Tw > Tsat and melting begins. Although conduction plays a dominant
role at early times when the tube is surrounded by a thin melt layer,
natural convection ultimately takes over. As a consequence, the in-
Fig. 1 Schematic diagram of the melting problem
terface between the solid and liquid is not vertical (in contrast with
the vertical interface for pure conduction), but is sloped as indicated
schematically in the figure. Under the influence of buoyancy, there
is a recirculating flow within the liquid with an upflow along the
heated tube and a downflow along the interface.
contain velocity terms, the heat transfer processes are not affected
The analysis and solutions were aimed at two general objectives. by buoyancy, and conduction remains as the sole transport mecha-
One was to obtain information such as heat transfer coefficients and nism. In [5] and [6], the freezing or melting of metals and alloys was
interface positions that are applicable to design. The second was to analyzed by solving the one-dimensional heat conduction equation
determine streamline patterns and temperature distributions in order for the solid; the heat transfer processes in the liquid were accounted
to obtain insights into the flow and energy transport processes. for via empirical natural convection heat transfer coefficients. Ex-
As formulated, the analysis involves the two spatial coordinates periments intended to clarify the interaction of natural convection
r, z as well as the time t. The complexity of the two-dimensional un- and phase change in one-dimensional geometries are reported in [7]
steady flow and heat transfer problem calls for the use of a numerical and [8].
solution procedure. The main parameters that were varied for the It thus appears that unsteady, multidimensional coupled natural
solutions were the Rayleigh number and the length-to-radius ratio convection and melting has not yet been analyzed from first principles
of the tube. The effects of two other parameters, the Stefan number in the published literature. In particular, the complex recirculating
and the Prandtl number, were also examined. natural convection flows which can occur in the liquid region have not
Attention will now be turned to previously published work dealing been analytically investigated, nor have their heat transfer ramifi-
with problems of melting and freezing. In the heat transfer literature, cations.
it has been standard to analyze such problems by assuming that
conduction is the sole transport mechanism, and the same approach The Analysis
has been followed in the applied mathematics literature. The effects The starting point of the analysis is the conservation equations for
of natural convection have been taken into account in a few papers. axisymmetric, unsteady flow and heat transfer in the melt layer. The
In [l], 2 natural convection was included in an analysis of the steady first task in specializing these equations is the formulation of the
state phase change problem associated with continuous casting. buoyancy term via the vertical pressure gradient and body force that
Natural convection effects have been investigated in the problem of appear in the z -momentum equation, that is
freezing or melting of water in horizontal layers, and a survey of this
work is available in [2]. Experiments on the melting of a sphere of ice -dp/dz - pg (1)
situated in a container of liquid water are described in [3].
The density appearing in the foregoing expression can be related to
In contrast to the foregoing, there appears to be a greater awareness the temperature via the Boussinesq model, according to which density
in other literatures, particularly metallurgy, of the role of natural variations are considered only insofar as they contribute to buoyancy,
convection in phase change problems. This awareness was exhibited but are otherwise neglected. If Tsat is the fusion temperature and psat
in [4] in connection with a one-dimensional analysis of freezing in a is the corresponding liquid density then, according to the Boussinesq
vertical slot. However, since the relevant energy equation does not model

P - Psat - PPsat(T - Tsat) (2)

'• Numbers in brackets designate References at end of paper. Furthermore, if a reduced pressure p * is defined as

.Nomenclature.

cp = specific heat at constant pressure Ri = dimensionless interface coordinate, u = axial velocity

Fo = Fourier number, at/r02 nlra v = radial velocity
g = acceleration of gravity r = radial coordinate z = axial coordinate
hsi = latent heat of fusion a = thermal diffusivity
r; = interface coordinate
k = thermal conductivity (3 = thermal expansion coefficient
P = dimensionless pressure, equation (12) ro = tube radius
f = dimensionless axial coordinate, z/ro
Pr = Prandtl number Ste = Stefan number, cp{Tw - Tsat)/hsL rj = transformed radial coordinate,
p = pressure T = temperature (r - r 0 )/(n - r 0 )
p* = reduced pressure, equation (3) TSat = fusion temperature 6 = dimensionless temperature,
Q = heat transfer rate at tube surface Tw = tube wall temperature (T — Tsat)/(TW — Tsat)
R = dimensionless radial coordinate, r/ro t = time elapsed from the start of heating v = kinematic viscosity
Ra = Rayleigh number, U = dimensionless axial velocity, u/(u/ro) p = density
(gP(Tu - Tmt) r 0 3 /« 2 )Pr V = dimensionless radial velocity, u/(a/ro) T = dimensionless time, equation (13)

Journal of Heat Transfer NOVEMBER 1977, VOL 99 / 521

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 07/24/2013 Terms of Use: http://asme.org/terms

 P* = P + PsatgZ (3) in which hsL is the latent heat for solid-liquid phase change. The
motivation for using a/ro rather than v/ro as a reference velocity and
then (1) becomes
for involving the Stefan number in the dimensionless time variable
-dp*/dz+gpPsat(T-Tsat) (4) will soon become apparent.
If the conservation equations were to be transformed in strict ad-
the second term of which represents the buoyancy force. herence with (10) and (11), the result would be a set of equations no-
Next, to further adapt the conservation equations, the constant table for length and complexity. In the main, the terms which elongate
property assumption is introduced and the subscript sat is deleted and complicate the transformed equations are spawned by the
from psat in (4). With the foregoing, the equations representing con- right-most member of equation (10), which involves the slope cfa/cfe
servation of mass, z -momentum, /'-momentum, and energy can be of the interface. The inclusion of these terms would add significantly
written as to the complexity and execution time of numerical solutions. Fur-
du 1 d (rv) thermore, it is highly likely that they would adversely affect the
— + — =0 (5) convergence of an iterative solution scheme.
dz r dr
Inasmuch as the present research represents the first analysis of
du du2 1 dlruv) ldp* multidimensional, unsteady natural convection phase change, it does
—+ + -= - + g0(T - Tmt) + vV2u (6)
dt dz r dr p dz not seem unreasonable to adopt a somewhat simplified model to en-
to d(uu) 1 d{ru2) 1 dp* able solutions to be obtained with reasonable computational effort.
--tL~+v(V2v-(v/r2)) (7) It will, therefore, be assumed that the interface radius varies slowly
dt dz r dr p dr
with z, so that the terms involving drjdz and d2rt/dz2 can be ne-
glected. This assumption is of the same nature as the local similarity
(8)
dt dz r dr model commonly used in the analysis of boundary layer flow and heat
transfer.
in which V2 is the Laplace operator in r, z coordinates. These equa-
tions, along with an interface energy balance and other boundary With the aforementioned assumption, the transformed version of
conditions, are to be solved in a domain whose boundaries do not all the conservation equations can be written as
fall along coordinate lines in a cylindrical coordinate system. In par- dU 1 d(RV)
ticular, the radius r; of the liquid-solid interface varies with z as well 0 (15)
df ' R(Rt - 1) dV
as with time as melting proceeds. This state of affairs precludes an
analytical solution and requires that, even in a numerical solution, 1
To dU
M dU2 1 d(RUV)l
— Ste — + +
approximations be- made in order to obtain results with a realistic PrL dr 3f R(Ri-l) dv J
computational effort. dP
For the computational model to be employed here, it will be as- + Ra9 + V2LT (16)
sumed that there is a small time lag between the heat delivery to the
dt
interface and the resulting interface motion. Specifically, during a 1 d(RV2)l
small time interval, the fluid flow and heat transfer in the liquid region PrL
Pr L dr 3f R(Ri - 1) dr, J
will be solved assuming that the interface is stationary. The heat dP -
transferred to the interface during that interval will then be used to +V2V-V/R2 (17)
compute the small finite displacement associated with the melting. dr\
With the new liquid domain which corresponds to the displaced in- „ de d(ue) 1 d(RVe) 2
terface, the solution for the next time interval is performed as just Ste — + — + ve (18)
described. dr 3f R(Rt - 1) dr,
The quantities Ri and R respectively represent the dimensionless
To transform the solution domain into a more tractable shape, new
radius of the interface, which depends only on z during a given time
coordinates f, r, are introduced. Consistent with the aforementioned
interval, and the dimensionless radial coordinate
model, the interface radius r,- is independent of time during any given
computational time interval. As a consequence, r; = ry(z) during each Ri = n/r0, R = r/r0 = r,(Ri - 1) + 1 (19)
interval. The new coordinates are defined as 2
Furthermore, V is the transformed Laplace operator
f = z/r0, (r - r0)/(n - r0) (9) -„ d2 I d
(20)
so that
V2 = — +
dt1 R(Ri-l)2dri
,
(»s)
V dr; d It may be noted that because of the choice of the dimensionless
(10)
dz ro<V (n ~r0) dz dr, variables, the inertia terms in equations (16) and (17) are multiplied
by (1/Pr) and the storage terms (i.e., d/<9r terms) contain Ste as a
a _ . i d .
(11) factor. Aside from the explicit presence of Pr and Ste, the fluid
dr" fa - r0) dr, properties that make up Pr appear implicitly in Ra, while Ste is in-
The key result of the transformation is that the range of r, extends cluded in the dimensionless time variable T. For typical salt storage
from zero to one at all f during any time interval. It may be noted that applications, both (1/Pr) and Ste « 1. The ramifications of these facts
the f, r, coordinates are nonorthogonal, as are the transformed coor- will emerge when the results are presented.
dinates in conventional boundary layer analyses. The four conservation equations (15)-(18) contain five unknowns.
In addition, to reduce the number of parameters that have to be These include the velocity components U and V, the pressure P, the
specified for the numerical solutions, dimensionless variables and temperature 6, and the interface radius R,. The fifth equation arises
groups are introduced as from the energy balance which applies at each point on the liquid-solid
interface. This balance states that the heat delivered by the fluid to
v the interface during a small time interval At equals the latent heat
U--
a/ro'
v-- P=- 2
p(a/r0) Pr absorbed by the melting material. The heat delivered AQ,- can be
T-Tsi written as
0= (12)
-* w 1 sat AQ,- = 2wridz(k(dT/dz)(dn/dz) - k(dT/dr))At (21)
T = (at/r02)(cp(Tw - Tsat)lhSL) = (Fo)(Ste) (13)
where the expression in brackets is a representative interfacial heat
g0(Tw - Tsat)r0s „ ^ c„(Tw - Tsat) flux for the time interval. Since the interface is a surface of constant
Ra = -Pr, (14)
hsL temperature (T = Tsat), it follows that

522 / VOL 99, NOVEMBER 1977 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 07/24/2013 Terms of Use: http://asme.org/terms

 dT/dz = -(dn/dz)(dTldr) (22)point in the f,j/ grid corresponds to different physical locations at time
TJ and time T,-+1. Since the present calculations were performed with
so that (21) becomes
a fixed array of grid points in the 1,-q plane, the field values at these
AQ; = -k(dT/dr)2irndz(l + (drJdz)2)At (23) grid points correspond to one set of physical locations at time T,- and
2
to another set of locations at T,+I.
For gradually or even moderately sloping interfaces, (drt/dz) « 1. Re-evaluation of the field variables at the grid points has to be
The volume swept out in response to the heat transfer AQ; is performed to take cognizance of this change of relationship between
2 2 2
jrA(r,' )rfz if (drt/dz) is neglected compared with one, where A(/"; ) the physical and transformed planes. The re-evaluation was carried
= nj+12 — rtj2 and j and 0 + 1) refer to the beginning and end of a out by interpolation at all grid points except those which correspond
time interval. With this and with the expression for AQ,-, the interface to the newly created melt region between ft;(f,TJ+1) and ft;(f,T,). In
energy balance is the latter, the velocities were assumed to be zero and the temperature
~k(dT/dr)2nAt = PhSLA(ri ) 2
(24) taken to be Tmt (i.e., 8 = 0). The interpolation for the velocities was
performed by first evaluating the stream function and then differ-
or, in dimensionless terms, entiating it at the locations of interest. This approach was employed
to ensure that both overall and local mass conservation were fulfilled.
To ensure energy conservation, the temperature interpolation was
carried out via an enthalpy integral, the subsequent differentiation
The boundary conditions will now be discussed. For the velocity, of which yielded the local value of 8.
both components are zero on all of the surfaces which bound the liq- To avoid potential computational difficulties at T = 0 (e.g., those
uid.3 The temperature is prescribed as being uniform on the tube associated with l/(ft; — 1) -»• °°), the solution was started by assuming
surface (T = Tw), while at the interface the temperature is Tsat', the existence of a very thin melt region concentric to the tube. The
correspondingly, 8 = 1 at r) = 0 and 8 = 0 at rj = 1. The upper and lowervalue of T corresponding to the assumed melt layer thickness was
bounding surfaces, z = H and z - 0 respectively, are assumed to be determined from the Stefan solution. Auxiliary computations were
adiabatic. performed to ensure that the assumed layer thickness was sufficiently
If the complete set of governing equations and boundary conditions small so as not to affect the subsequent results.
is examined, it is seen that there are four prescribable parameters: Ra, At each time step, it is necessary to solve a set of simultaneous al-
H/ro, Pr, and Ste. The relative importance of these parameters will gebraic equations which represent the implicit finite-difference
become evident from the numerical results. discretization of equations (15)-(18). This two-dimensional elliptic
problem was dealt with using a numerical scheme which is an adap-
Numerical Solutions tation of that of [9], In actuality, the elliptic finite-difference scheme
The solutions were carried out by means of an implicit finite-dif- described in [9] is employed there as one of the building blocks in a
ference procedure. The most challenging problem that had to be faced procedure for dealing with three-dimensional parabolic flows. From
in devising a solution methodology was the movement of the interface. the elliptic scheme of [9], which is set forth there for Cartesian coor-
After examination of various possibilities, it appeared convenient and dinates, the discretization, grid structure, and solution method have
effective to adopt the already outlined procedure which allows for a been utilized here in constructing a computer code tailored to the f,t;
small time lag between the convective flow and heat transfer in the coordinates of the present problem.
liquid and the melting and movement of the interface. The finite-difference grid encompassed 12 X 14 nodal points for
To elucidate the approach, suppose that at an instant of time TJ, H/r0 = 4 and 12 X 20 nodal points for H/r0 = 10, respectively in the
all the field variables (i.e., U, V, P, and 6) are known, as is the di- •q and f directions. The points were uniformly distributed in ij. On the
mensionless interface positionft;(f, TJ). Then, for a small time interval other hand, the distribution along the f direction was skewed so as
AT = (T/+I — TJ), the interface is regarded as fixed. For this fixed- to have a greater concentration of points near the top of the domain
boundary domain, the implicit finite-difference form of equations in order to accommodate the rapid turning of the recirculating flow
(15)-(18) is solved, subject to their boundary conditions, to yield the in that region. The time duration of the solutions extended to T ~ 0.20
distributions of the field variables at time T ; +I. From the distribution to 0.25, which corresponded to a maximum value offt,-of about 3. The
of 6, the values of the derivative dSld-q along the interface, which are time step AT required for accuracy and stability depended on the
proportional to the local heat flux, can be determined. Then, by Rayleigh number, the length-radius ratio, and the value of T itself. In
solving the interface energy balance (25), the position of the interface general, smaller values of AT were needed at higher Ra and H/r0. The
at time T,-+I can be evaluated; that is overall range of AT used in the computations extended from 0.001 to
0.01.

Results and Discussion

where, on the right-hand side, ft; =ft;(f,T/).
In selecting values of the governing parameters to be used as input
With the determination of ft; from equation (26), all quantities, i.e.,
to the numerical solutions, guidance was taken from thermal storage,
the field variables andft,-,are known at time T/+I. An identical pro-
although the results are in no way limited to that application. The first
cedure can then be employed to proceed to time TJ+2, and so forth.
estimates of Ra, Pr, and Ste were based on the sodium nitrate-sodium
The effective implementation of this approach requires that AT be
hydroxide eutectic which had been used in the laboratory experiments
chosen sufficiently small to avoid errors associated with the lag be-
mentioned in the Introduction. These estimates gave Ra ~ 7 X 105,
tween the heat delivery to the interface and the resultant interface
Pr ~ 7, and Ste ~ 0.15. For a parametric study, Rayleigh numbers
displacement.
between 7 X 104 and 7 X 106 appeared to be relevant.
Another feature of the solution method is the re-evaluation of the
The main calculations were, therefore, carried out for Rayleigh
field variables subsequent to the determination of ft; from equation
numbers of 7 X 104,7 X 106, and 7 X 106, with Pr = 7 and Ste = 0.15.
(26). The need for the re-evaluation stems from the fact that the re-
For each of these cases, results were obtained for two length-radius
lationship between a given radial position r and the transformed
ratios of the heated tube, H/r0 = 4 and 10. In recognition of the higher
coordinate i] changes as r; changes (e.g., equation (9)). Thus, a fixed
Prandtl numbers and lower Stefan numbers that are characteristic
of the liquid phase of certain phase change materials, check runs were
made for Pr = 70 and for Ste j= 0.05.
3
At the interface, the normal component of the velocity is zero when the density The heat transfer results will be presented first, followed succes-
change associated with the phase change is neglected, as in the present analy- sively by representative temperature field and flow field information
sis. and by the time evolution of the liquid-solid interface.

Journal of Heat Transfer NOVEMBER 1977, VOL 99 / 523

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 07/24/2013 Terms of Use: http://asme.org/terms

 If the rate of heat transfer at the tube surface is denoted by Q, then at higher Rayleigh numbers and for shorter heating tubes. However,
Q/H is the heat transfer per unit length of tube. In Fig. 2, the heat owing to the different paces of the various time histories, there is some
transfer per unit length, represented in dimensionless terms by crossing of the curves which locally invalidates these general
(Q/H)/k{Tw - Tsat), is plotted as a function of time T (= Fo Ste) for trends.
all of the cases that were investigated. The solid and dashed lines Fig. 2 also shows that when the Rayleigh number is used to char-
appearing in the figure correspond respectively to length-radius ratios acterize the natural convection effects, the results do not depend on
H/ro of 4 and 10. The upper, middle, and lower pairs of curves re- the Prandtl number for the range investigated (Pr = 7-70). Fur-
spectively portray the results for Ra = 7 X 106,7 X 106, and 7 X 104. thermore, it is evident from equations (16)-(18) that if Pr has no effect
For all of these cases, calculations were performed with Pr = 7 and for this range, it will not have an effect for any Pr > 7. On the other
Ste = 0.15. However, as noted in the figure, the uppermost curve hand, had the Grashof number (Ra/Pr) been used, then a Prandtl
corresponds to calculations made with both Pr = 7 and Pr = 70. In number dependence would have been encountered.
addition, one of the middle curves represents the solutions for both It is also seen from Fig. 2 and equations (16)-(18) that the in-
Ste = 0.15 and 0.05. volvement of the Stefan number in the dimensionless time variable
Inspection of the figure reveals a complex variation of the heat T eliminates the dependence of the results on Ste in the range from
transfer with time which requires careful consideration. However, one 0 to 0.15. In view of these findings, both the Prandtl number and the
fact is immediately evident, namely, that the heat transfer results are Stefan number can be eliminated from the list of active parameters
drastically different from the one for pure conduction phase change, in the aforementioned ranges. At higher values of Ste, a parametric
which is characterized by a monotonic decrease with time. Clearly, dependence of the results on the Stefan number may be antici-
another transport mode, natural convection, is playing a first order pated.
role. Attention will now be turned to illustrating representative types
Within the time scale of the figure, a greater number of events in of temperature distributions that exist in the melt layer at various
the time history of the heat transfer are in evidence at higher Rayleigh stages in its time history. Fig. 3 has been prepared for this purpose.
numbers. From the curves for Ra = 7 X 106, it is seen that the heat The figure consists of three panels. In each panel, radial temperature
transfer first decreases and attains a minimum, then increases and profiles spanning the melt layer are plotted at four axial stations, z/H
achieves a maximum, and subsequently decreases. At lower Rayleigh = 0.15,0.45,0.70, and 0.93. The Rayleigh number and the time T that
numbers, the occurrence of these events is shifted to larger and larger characterize the profiles in each panel are indicated therein, and H/ro
times, so that the maximum and the subsequent decrease are not seen = 4.
in the figure. The results in the left-hand panel are representative of early times,
These events in the time history of the heat transfer can be tied to when the influence of natural convection is small. The dashed lines
certain physical processes which will be mentioned briefly here and plotted in the figure represent the logarithmic temperature distri-
illustrated later via the temperature and fluid flow results. At early bution for steady state, one-dimensional radial conduction. The solid
times, the natural convection, although present, is very weak, so that and dashed lines are nearly coincident, thereby not only confirming
conduction is the dominant transport mode and, as a consequence, the dominance of the conduction mode, but also indicating a negligible
the heat transfer decreases with time. As the melt layer thickens, effect of transient energy accumulation in the fluid.
natural convection grows stronger, thereby arresting the decrease in The middle panel is intended to illustrate a temperature field which
heat transfer. The further strengthening of the natural convection is strongly affected by natural convection, but where the boundary
causes the heat transfer to increase. During this period, heat is layer regime has not yet been established. In the panel at the right,
transported by the recirculating natural convection flow as well as the boundary layer regime is clearly in evidence. The temperature is
directly across the melt layer. uniform across most of the thickness of the melt layer, with sharp
When the melt layer becomes sufficiently thick, boundary layers
are respectively established adjacent to the tube wall and adjacent
to the interface. As a result, heat is no longer transported directly
across the melt layer, but is carried from the tube wall to the interface 1,0
only by the recirculating flow. In addition, the growing resistance - \
encountered by the recirculating flow along the top and bottom walls
0.5 Z/H = 0.93 \
(primarily the former) tends to impede the fluid motion. It is believed
that these factors are responsible for the decrease in heat transfer at
large times that is evident from the curves for Ra = 7 X 106. 0 „ Ro=7xl04 ^ Ro=7xl0 4 Ra= 7x10
In general, the heat transfer per unit length appears to be greater 1.0 T =0.015 T =0.21 T =0.20

0.5 -
0.70 ~ ~ \

0 -
.Pr- 7 and 70 T-T 5 0 |
/ Ro = 7xl0 6 1.0
* \ H/r0 T„-T s a ,

l0
^<o->
0.5 -
"i 0.45~""\ "A
<-" 6 0 0
k
Ste = 0.15 and 0.05 1.0
CONDUCTION
7x!0 5
40

_ _ _ _ _
0.5
X. \0.I5
0
~>
•20
, 1 1 1 1 1 1
I 1
0.5 0.5 0.5

V V

Fig. 2 Timewise variation of the tube heat transfer rate Fig. 3 Representative temperature profiles (H/r0 = 4)

524 / VOL 99, NOVEMBER 1977 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 07/24/2013 Terms of Use: http://asme.org/terms

 gradients adjacent to the tube wall and the liquid-solid interface. conduction dominates, a natural convection flow definitely exists. The
The nature of the velocity field and the natural convection recir- flow pattern that typifies most of the time history of the melt layer
culation can be well illustrated by streamline maps, and three rep- is illustrated by the middle diagram. There is a single recirculation
resentative flow patterns are shown in Fig. 4. Each diagram is an r,z zone that is centered in the upper portion of the melt region, where
plane encompassing the melt region. For all cases, the general char- the melt layer thickness is'greatest. The three-eddy array shown in
acteristic of the flow pattern is an upflow adjacent to the heated tube the right-hand diagram was encountered only for the longer melt layer
and a downflow adjacent to the liquid-solid interface. (i.e., H/r0 - 10) and only at times corresponding to the neighborhood
The diagram at the left is representative of early times, and the of the heat transfer maximum (Fig. 2).
trace of the melt layer in the r,z plane is nearly rectangular. Although The streamline diagrams indicate that the radial extent of the melt
layer varies along the length of the heating tube, with the smallest
layer thickness at the bottom and the greatest layer thickness at the
top. This is in contrast to the melt layer of uniform thickness that
r/r0 r/r0 r/rn would be predicted by an analysis in which conduction is the only
mode of transport.
I.I 1.0 2.0 1.5 1.0
i.o A more complete presentation of melt layer thicknesses and in-
C7\ terface shapes is provided by Figs. 5 and 6, respectively for H/r0 = 4
and 10. Each figure contains two panels, corresponding respectively
to Ra = 7 X 104 and 7 X 106. Each panel is an r,z plane in which the
heating tube encompasses the region 0 < r/ro < 1 and the phase
change material occupies the space r/r0 > 1. The sequence of curves
represent the successive positions of the interface at the times indi-
cated in the figure.
The nonuniformity of the melt layer thickness that was in evidence
in the streamline diagrams is confirmed by Figs. 5 and 6. In addition,
these figures indicate that the interface shape and its evolution are
affected markedly by the Rayleigh number and, to a lesser extent, by
the length-radius ratio. At high Rayleigh numbers, the interface slopes
0.5
gradually outward, with the slope being more or less uniform. On the
other hand, at low Rayleigh numbers, the interface tends to lie parallel
to the tube for a substantial length starting at the bottom and pro-
ceeding upward. Near the top of the region, the interfaces flare out.
From these results, it is evident that at lower Rayleigh numbers,
significant natural convection motions are confined to the uppermost
portion of the melt region; at higher Rayleigh numbers, significant
flows also occur in the lower portion of the melt. In the limit as Ray-
leigh number approaches zero, the interface will be vertical since
natural convection vanishes and conduction is the sole transport
mechanism.
Figs. 5 and 6 enable an assessment to be made of the conditions
under which the interface slope dri/dz is small, as assumed in the
(a) (b) (c) analysis. At high Rayleigh numbers, and especially at large H/r0, the
Ra = 7 x l 0 4 Ra = 7 x l 0 6 Ra = 7 x l 0 6 slopes are small to moderate. For lower Rayleigh numbers, the in-
T = 0.0I5 T= 0 . 0 4 4 T = 0.040 terface slope is small along most of the length, but may be quite large
H/r 0 = 4 H/r 0 = 4 H/r 0 = IO near the very top of the melt region. It is here that some deviations
of the computed results from reality are to be expected, but the
Fig. 4 Representative flow patterns qualitative trends should not be affected.

Fig. 5 Melt layer thicknesses and interface shapes, H/r0 = 4

Journal of Heat Transfer NOVEMBER 1977, VOL 99 / 525

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 07/24/2013 Terms of Use: http://asme.org/terms

 The calculated heat transfer rates at the tube wall were found to differ
markedly from the monotonic timewise decrease that is characteristic
of pure conduction phase change. Rather, the heat transfer decreased
at early times and attained a minimum, then increased and achieved
a maximum, and subsequently decreased. The times at which these
events occurred were markedly affected by the Rayleigh number. At
lower Rayleigh numbers, the pace of the timewise heat transfer
variations was much slower than at higher Rayleigh numbers. The
events in the time history can be tied to physical processes which were
illustrated by temperature and velocity field information presented
in the paper.
Another indication of the marked departures from the pure con-
duction model can be found in the shape of the melt region. For con-
duction, the thickness of the melt region is uniform along the entire
length of the heating tube. On the other hand, in the presence of
natural convection, the thickness of the melt layer varies along the
length, with the smallest layer thickness at the bottom and the
greatest layer thickness at the top.
These findings indicate that natural convection effects, although
unaccounted for in most phase change analyses, are of first order
importance and have to be considered.

Acknowledgment
This research was performed under the auspices of ERDA contract
No. E(ll-l)-2595.

References
1 Kroeger, P. G., and Ostrach, S., "The Solution of a Two-Dimensional
Freezing Problem Including Convection Effects in the Liquid Region," In
ternational Journal of Heat and Mass Transfer, Vol. 17, 1974, pp. 1191-
1207.
2 Seki, N., Fukusako, S., and Sugawara, M., "A Criterion of Onset of Free
Convection in a Horizontal Melted Water Layer with a Free Surface," JOUR-
NAL OF HEAT TRANSFER, TRANS. ASME, Series C, Vol. 99, 1977, pp.
92-98.
3 Vanier, C. R., and Tien, C, "Free Convection Melting of Ice Spheres,"
Fig. 6 Melt layer thicknesses and Interface shapes, H/r0 = 1" AIChE Journal, Vol. 16,1970, pp. 76-82.
4 Szekeley, J., and Stanek, V., "Natural Convection Transients and Their
Effects on Unidirectional Solidification," Metallurgical Trans., Vol. 1,1970,
pp. 2243-2251.
5 Szekeley, J., and Chhabra, P. S., "The Effect of Natural Convection on
Concluding Remarks the Shape and Movement of the Melt-Solid Interface in the Controlled Solid-
In contrast to the conventional conduction-based treatment of ification of Lead," Metallurgical Trans., Vol. 1,1970, pp. 1195-1203.
solid-liquid phase change, the present analysis has taken account of 6 Chiesa, F. M., and Guthrie, R. I. L., "Natural Convective Heat Transfer
Rates During the Solidification and Melting of Metal and Alloy Systems,"
natural convection induced by temperature differences in the liquid JOURNAL OF HEAT TRANSFER, TRANS. ASME, Series C, Vol. 96,1974,
melt. The analysis was performed for melting of a solid due to heat pp. 377-384.
supplied by a vertical finite-length tube embedded in it. Initially, no 7 Cole, G. S., "Temperature Measurement and Fluid Flow Distributions
liquid is present, and the solid, which is at its fusion temperature, is Ahead of Solid-Liquid Interface," Trans. AIME, Vol. 239, 1967, pp. 1287-
1295.
in contact with the tube. At subsequent times, the tube temperature 8 Heertjes, P. M., Jongenelen, J. A., and deLeeuw den Bouter, J. A., "The
is maintained at a uniform value that exceeds the fusion temperature, Effects of a Moving Boundary on Heat Transfer by Free Convection," Chemical
and the progressive melting of the solid gives rise to a melt region Engineering Science, Vol. 25,1970, pp. 1881-1890.
which grows larger as time passes. 9 Patankar, S. V., and Spalding, D. B., "A Calculation Procedure for Heat,
Mass, and Momentum Transfer in Three-Dimensional Parabolic Flows," In-
Solutions were obtained by an implicit finite-difference procedure ternational Journal of Heat and Mass Transfer, Vol. 15, 1972, pp. 1787-
tailored to take account of the movement of the solid-liquid interface. 1806.

526 / VOL 99, NOVEMBER 1977 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 07/24/2013 Terms of Use: http://asme.org/terms