Freezing and Sublimation under various
Boundary Conditions for Cylinders and Spheres.
Rahul Basua,
a Department of Mechanical Engineering,
Sambram Insiitute of Technology, VTU MS Palya, Jalahalli E, Bangalore 560097. email: [email protected]
. Abstract
A computer simulation using Wolfram was performed on the equations describing Moving Boundary
Problems in Porous substrates, as related to solidification and Self sublimation under various
boundary conditions. The effect of sources and sinks was also included, and exact solutions were
derived for the phase front velocity. Only the thermal effects are considered, neglecting the effect of
convective flow. Very little work has appeared dealing with exact solutions, due to the complexity of
solving the Moving Boundary Problem which is highly Nonlinear, and depends on solution of
transcendental equations.. It was found that the geometry as well as the type of boundary conditions
affected the phase front velocity, and specific results obtained earlier, have been verified.
keywords: Diffusion, MBP (Moving Boundary Problem ),Phase transformation, Porous Media, Stefan Problem,.
1. Introduction
The majority of the literature has dealt with stationary boundaries in connection with
solidification and sublimation problems. The literature in heat conduction is replete with
adaptations from text books where standard boundary conditions are applied and various
methods of solution are displayed with numerical schemes and variations, relying on the
computer. However, before the computer arrived on the scene, much of the calculation had
to be done by hand, and in case of the moving boundary problem described here it was a
difficult task more so because of the inherent non linearity of the MBP. Not only was the
theoretical solution extremely difficult to obtain, the subsequent solution was even more
intractable because of the transcendental equation appearing for the unknown moving
boundary.The location of the moving boundary was also to be obtained as part of the
solution. According to Carslaw and Jaeger [1], only a few solutions are available for
various geometries for moving heat sources. The classic Moving Boundary Problem was
known as the Stefan Problem, originating with a study of the Arctic Polar ice over a 100
years ago [2]. As of today, people are still struggling with this problem, and no exact
solution has been obtained for many standard cases.
The Problem of freezing, melting and sublimation with a moving interface has been
described also in Ingersoll [3], Crank [4], Ozisik [5], Incropera [6] and several others.
Various approaches to tackle the problem are transformation of coordinates and variables,
variable parameter scaling, use of Greens Functions, Perturbation analysis, and Enthalpy
transformations.
Many applications of the work apart from formation of polar ice can be cited today.
Among these are laser surface modification, geology, formation of ice on underground
cables, ablation of space craft shields, food preservation, heat treatment of steel and
cryosurgery. Modification of surface properties by point and line sources is used in various
applications like etching, rapid prototyping, and IC manufacture, supplementing and
replacing traditional treatments like nitriding, carburizing and surface work hardening.
1.1. Mathematical Formulation of the General Problem
I consider first the problem of solidification or melting of a cylinder and a sphere
without any mass diffusion terms. A solution is given by Paterson [7] for a cylindrical
and a spherical boundary propagating from a line and a point source. The initial
temperature is referring to the melting point as zero. A boundary of fusion r = R(t)
propagates into the medium. Denoting the solid region as 1 and liquid region as 2,
assume a solution
1 = -A Ei(- r2/4 Kt ) –B (1)
2 = -C Ei (-r2/4Kt) –D (2)
R(t) = a(t) , where a line heat source exists given by Q =q(t)0.5
, (3)
considering a point heat source in an infinite medium. The solutions in the two
regions are
1 = A[ (K1t)0.5
/r exp( -r2/K1t) – (
erfcrt
Similarly for region 2
Applying the Boundary conditions, when t =0, R =0, and after some simple algebra
the energy balance at the interface is given by:
K1|grad 1| - K2|grad 2| = L V = L 1t/|grad 1| (5)
2. Simulations using WOLFRAM for the Sphere
Setting particular values for the physical parameters, solutions can be obtained by
solving the transcendental equation. The solution is a single positive root as the LHS
decreases monotonically ( as seen in the attached graph).
FIG 1. Solution for the sphere, liquid solid equilibrium,
Where we are using the values Water: k2 =0.00144 Cal/cm sec K, K2 = 0.00144
cm2/sec,Ice: k1 = 0.0053 Cal/cm.sec.K, K1 = 0.0155 cm2/sec, L = 73.6, Q is taken as 2.38
cal/cm.sec, =-2
The solution is obtained as 0.09763.
Consider now the boundary with a zero flux condition imposed at r =0. By transforming to the rectilinear case, it is seen to be equivalent to an insulated boundary. In terms of the foregoing analysis, this can be arrived at by putting the constant C in the general solution
to be zero, since d/dt =C (K1t)/r2 exp ( - r2/4K1t), hence for zero flux, C is zero.
The interface position can be obtained by neglecting the first term in the thermal balance equation. The solution for this case gives , -.01607 , which is less than the prior case with the heat sink.
X =-0.016017
For a zero heat flux at the boundary, there is a singularity of the eror function derivatives at x=0. Thus, an adiabatic condition ( dT/dx =0) at the boundary of a rectilinear slab problem or origin in case of the sphere , would lead to a flat profile for the temperature. That is the general solution given in [ ] cannot be satisfied at the boundary with one temperature and the phase interface at the melting temperature (taken as 0 in reduced variables). Several authors have attempted the derivative boundary condition including, [10],[11] . An exact solution has eluded all who have attempted the derivative condition because of the difficulty posed by this property of the error function. A similar problem was analysed by McCue [ 10 ], where the problem is set up and stated that no exact solution exists. However, as the problem is stated, the conditions of constant temperature and adiabatic wall cannot lead to compatibility of the solution at the wall or the interface, since a flat or constant temperature is the only solution possible. The wall temperature cannot be different from the interface temperature. If it is different there would be a jump in temperature somewhere along the inner phase, most likely at the moving boundary. However the condition at the moving boundary assumes both phases have the same transition temperature, so this condition cannot be satisfied. The addition of sources and sinks can lead to a thermal profile which is flat at the boundary, and is given in Carslaw Jaeger [1], for the adiabatic wall.
2.1 SIMULATION FOR THE CYLINDRICAL CASE: Using a similar procedure, the same equation is solved with a slight change in the terms, replacing SphErf with Ei and the cubic functional dependence on the rhs is seen to become quadratic. A solution for the cylinder with the line sink was given by Paterson [7 ]. The solution is recomputed with Wolfram, and found to give x= 0.0669178, as compared to x= 0.06637 earlier [ 7].
When the heat sink is removed, the solution is imaginary. If the sign of L is changed, a real root is found as 0. ( Note Ei(x) is defined for x>0 and x<0, at x=0 the function has a branch cut and can be defined as bracketed between two limits, and thus defined as infinite.
1. . [ from Abramowitz
and Stegun[12]
Graphical representation of the root for the cylinder:
, x=0
Fig.2 Simulation for the melting cylinder with zero heat sources, giving zero solution.
3. SUBLIMATING CYLINDER AND SPHERE, WITH AND WITHOUT HEAT
SOURCE.
3.1 Sublimating cylinder In the subsequent analysis and results, a cylinder and sphere are solved for direct
phase change from solid to vapour, termed as ablation or sublimation. In such a case the
net latent heat is the sum of the intermediate phase transition heats. The parameters for
water vapour are taken as density = 1.694 kg/m3 at 1 bar, diffusivity 2.338 x 10
-5 m
2/s,
conductivity varying from 0.16 to 0.248 W/m/deg K, L = 540 cal/gm.
With the same values for the heat source, the thermal balance equation is set up for
ice and vapour, using the Latent heat value for sublimation from solid to gas as
(540+73.6_ or 613.6 cal/gm.
For the cylinder results are depicted in the figure below, with a value for the
sublimation eigenvalue below each figure .
Fig 4 Simulation for sublimating cylinder with heat sink
3.2 SELF FREEZING ( no heat sink)
The figure 7 shows the solution when the cylinder sublimates drawing energy from the
solid phase behind the phase front, without passing through the intermediate liquid
phase. The equation is directly obtained from the previous thermal balance when q =0
(no imposed heat source). The thermal energy for the phase front motion is taken from
the latent heat of sublimation
Fig 5 Simulation for self freezing Cylinder with zero heat sink
Solution
----------------------------------------------------------------------
3.2. Sublimating Sphere with and without heat source:
Similar calculations were done for the sphere, replacing the Ei function with the
Spherical error function. Results are appended below.
Figure 6. Spherical simulation with heat source for ablating or sublimating sphere,
solution
For the case where no heat source or sink is present, (Self freezing or self
sublimating), the solution is obtained from Wolfram as in Figure 7:
Figure 7: Self sublimating sphere, no heat source
For opposite sign of latent heat, we apply Wolfram solver:
4. RESULTS
The numerical results can be tabulated in Tables 1and 2 in terms of coordinate
geometry, and presence of source or sinks.
TABLE 1: water ice system
Geometry Heat sink adiabatic wall(self)
Cylindrical 0.0669 0.000
Spherical 0.09763 -.01607
TABLE 2: Ice vapour system
Geometry Heat sink Adiabatic(self)
Cylindrical 0.0766 0.000
Spherical 0.0847 -0.0002
(note these values are very sensitive to round offs and errors from division of terms.
The negative value implies the interface moves in the opposite direction)
In addition, the values for the sublimating freezing cases with change of sign for L give completely different results from those for the positive L, due to the highly nonlinear transcendental equation which must be solved. The above values are some of the physically plausible values which were derived.
The cylinder has a markedly different behavior from the sphere when a zero flux is
applied at one boundary. In both liquid solid and solid vapour transitions, the phase
velocity is zero because of the property of the Ei function.
4. DISCUSSION
The examples calculated here are for the water ice vapour systems to illustrate the
solution methodology, for one particular set of Boundary conditions treated by Paterson
in 1952, before electronic computers were common. The example given for the cylinder
with a line heat sink has been verified to 3 decimal places. Other geometries such as the
sphere are also computed. In order to tackle other boundary conditions like the Dirichlet
and convective, most published references have tackled these by using numerical
methods. Pedroso Domoto [8 ] applied perturbation methods to the sphere. More recent
ones are those since 1980’s by Soward,[ 9 ] continuing that of Pedroso Domoto. The
problem tackled by McCue[10 ] for the adiabatic boundary with a zero heat flux, is
found to be physically inconsistent in that the error function cannot have a zero derivative
at one end and continue as a curve rising to or falling to the phase transition temperature.
It has to be a constant to be consistent with the theoretical solution of Paterson, which
satisfies the diffusion equation at all boundary conditions. The only other possibility is
for the temperature to vary linearly with x, which is a steady state solution with no
moving boundary (at infinite t). This is not a valid solution for the transient MBP
however.
5. CONCLUSION
Using the theoretical solutions for the cylinder and sphere, the moving boundary problem
was solved for the moving interface under various conditions of heat sink and adiabatic
boundary. It was found that for the cylinder the boundary does not move in the adiabatic
(zero heat source) case, whereas in general, the spherical boundary propagates faster than
the cylindrical case both for heat sink and adiabatic ( zero heat source) conditions.
6. References [1] Carslaw H S and Jaeger J C: Conduction of Heat in Solids, Clarendon, Oxford (1959)
[2] Stefan J:AnnPhys Chem, NF42, pp269-86, (1891)
[3] Ingersoll LR, Zobel O J, Ingersoll AC, Heat Conduction With Engineering, Geological, and Other
Applications, Oxford IBH, Calcutta, (1969)
[4] Crank J C. : The Mathematics of Diffusion, Clarendon, Oxford,(1975)
[5] Ozisik MN: Heat Conduction, J.Wiley, NYC, (1980)
[6] Incropera FP, Dewitt DP, Bergman TL, Lavine AS: Fundamentals of Heat and Mass Transfer, J.Wiley,
NDelhi, (2007)
[7] Paterson S: Propagation of a Boundary of Fusion, Glasgow Math Assn Proc., 1, pp42-47, (1952)
[8] Pedroso RL and Domoto GA: Perturbation Solution, J Heat Transfer, 95,pp 42-46, (1973)
[9] Soward AM: A Unified Approach to Stefans Problem for Spheres and Cylinders, Proc Roy Soc Lond A,
373, pp131-147 ( 1980)
[10] McCue S W, King J R , Riley DS: Extinction Behaviour for two dimensional inward solidification
problems, Proc Roy Soc (A), 459, pp 977-999 (2003)
[11] Stewartson K and Waechter RT, On Stefans Problem for Spheres, Proc Roy Soc Lond A, 348,pp415-
426, (1976) [12]. Abramowitz, Milton; Irene Stegun: Handbook of Mathematical Functions with Formulas, Graphs, and
Mathematical Tables. Abramowitz and Stegun. New York: Dover. ., Chapter 5, p229, (1964)
Nomenclature
A,B,C,D constants
a(t),R(t) interface position with time
K diffusivity
k conductivity
L latent heat of transformation (melting or vapourization)
Q heat sink/source strength
s, r coordinates in rectilinear and cylindrical, or spherical coordinates
t time
v velocity
density
non dimensional temperature
erf error function
grad gradient function, Ei exponential integral function