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: raulbasu@gmail.com

. 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