Heat transfer to an infinite solid through a moving boundary

Int. Comm. Heat Mass Transfer, Vol. 24, No. 2, pp. 181-189, 1997 Copyright © 1997 Elsevier Science Ltd Printed in the USA. All rights reserved 0735-1933/97 $17.00 + .00

Pergamon

PII S0735-1933(97)00004-3

H E A T T R A N S F E R TO A N I N F I N I T E S O L I D THROUGH A MOVING BOUNDARY

C.H. Leung, J.E.J. Staggs and A.C. McIntosh Department of Fuel and Energy University of Leeds, Leeds. LS2 9JT ENGLAND

J. Brindley School of Applied Mathematical Studies University' of Leeds, Leeds. LS2 9JT ENGLAND

(Communicated by P.J. Heggs)

ABSTRACT The objective of this paper is to investig~ate the classical Stefan problem of a one dimensional heat transfer proolem with a moving boundary. To solve this problem we use a known perturbation scheme to obtain an approximate solution when the ratio of sensible heat to the latent heat is small. We compare a numerical solution of this problem with the asymptotic and investigate the variation of the asymptotic regression rate with key, parameters. This method is valid for small values of the Stefan parameter and the motivation of this model is that it can be applied to a similar problem involving char formation. © 1997 Elsevier Science Ltd

Basic Physical Model and Mathematical Formulation

We consider a one dimensional semi-infinite solid subjected to constant heat flux F on its top surface. The solid heats up until the top surface reaches the material's ablation temperature. The material then vaporizes and a phase change develops. The temperature in the solid obeys the one dimension heat conduction equation,

T, = ~ T ~ . 181

(3)

182

C.H. Leung et al.

Vol. 24, No. 2

On the moving boundary, some of the applied heat flux will vaporize the material aud lhe rest is conducted into the solid. Moreover, the moving boundary is always al the ablation temperature, Tp. Suppose y = ,> is the position of the movillg boundary with respect to a stationary reference frame. Then, at y = s, the temperature of the material is the ablation temperature,

T=G.

(2)

At the boundary, the heat input is equal to the sum of the energy used for vaporization and the heat that is conducted into the solid, thus F=pL~

ds

- ,~ 7~.

(3)

Note that the second term is negative indicating that heat travels in the direction of negative temperature gradient. The boundary condition as g + ~c is given by, T --+ 0 as 9 ---> vc.

(4)

In this work we have made an assumption that the ambient temperature is zero. This is not crucial but simplifies the analysis. (Another way to look at the problem is to have Tp equal to the difference between the ablation temperature and the ambient temperature. ) Dimensionless Form The dimensionless form of the problem is obtained in a similar way to that defined in references [1], [2], by first introducing an evaporation controlled limit for the speed of the moving front. This limit arises when the temperature distribution ahead of the evaporating boundary approaches a steady state which migrates at a constant velocity with the moving boundary. All velocities are scale against this limiting velocity v . For a constant heat flux F applied on the surface for time at, an amount 5s of material is evaporated. By conservation of energy, p(Lv + c~T~) as = m t Rearranging the above equation and letting at -+ 0, we obtain ds

F

d-~-Ilimit - P (Lv + cp Tp) - constant - v,

(5)

Vol. 24, No. 2

HEAT TRANSFER TO AN INFINITE SOLID

183

where L, is the Latent heat of ablation and cv is the specific heal capacity. \¥e set the following dimensionless variables, T

e

y

G

vl

¢=7

T

s

/ '

(6)

(--7'

where l is the length scale defined as, l -

v

-

tp

(7)

'

where e = cv Tv/Lv is the Stefan parameter. After the nOnldimensionalisation, the problem becomes,

(9,=(9~.

(<¢
(S)

with b o u n d a r y conditions, =0,

and

( 9 = 1,

at

¢=(,

(9)

and @--+0 as ~ - * oc.

(10)

Transformation Since the material is bounded by a moving boundary of unknown location, instead of front tracking, a front fixing method is used to 'immobilize' the boundary as well as 'fixing' the domain. A transformation is used to map an adaptive grid onto the material domain so that the domain always lies between 0 and 1 in the new coordinate. The transformation is

(II)

z = tanh[a(~ - ((r))], where a is a scaling factor in this coordinate transformation.

After the transformation

equation (8) becomes

®~ = Ozz

+ O~ \dCZ + dr] dr] '

o<

z

< i,

(12)

with b o u n d a r y conditions dr-1

(17-e)

dz[,, - f

and @ = 1

~dzl,~

at

z=O,

(13)

and (9--~0

at

z=

1.

(14)

184

C.H. Leung eta[.

Vol. 24, No. 2

Numerical Scheme

A second order central difference method is used to approximate the spatial derivatives, a forward difference in time is used, and the Crank-Nicolson method [3] is applied to the difference equations. The advection term in equation (12) does not, in general, affect 1he stability of the numerical scheme [4]. Moreover, the scheme can be shown to converge by successive runs with a finer mesh and a reduced time step . The algorithm used in solving the PDE, iterates the Stefan condition (13) tor a suitable time step A t , with a fixed step change in the moving boundary position A( (that is 'A( constant). The temperature distribution in the solid is found by solving equation (12) by controlling the change in ( instead of time and using the approximation that d(/dr ~ A ( / A r provided A{ and A r are sufficiently small. The Stefan condition (13) is re-written as A~ = A~

d-~

(1

dz I,, -

d-dT[~d~ I,~

at

z = 0,

15)

where A t ( = "r2 - T1) is the finite change in time and A( is the finite change in the position of the moving boundary over the elapsed time A~- .

Results and Discussion

Pre-Ablation Region (d~/d'r = O) In the initial stages, the solid is at ambient temperature and a constant heat flux is applied to the material top surface. We solve equation (8) (with ~ = 0) for the temperature distribution in a semi-infinite solid with a fixed heat flux at one end and a constant temperature at the other end. The fixed heat flux can be obtained by using the boundary condition - 0 0 / 0 ~

= (1 + e)/¢ in place of (9). The solution to equation (8) is trivial and it

can be shown [5], that when the surface temperature reaches unity, the normalised form of the solution is,

O=exp ( - \ {~(l---+e)'~) 2)

~(l+e)erfc(~(l+e)~

i

\

)

(16)

Vol. 24, No. 2

HEAT TRANSFER TO AN 1NFINITE SOLID

185

The numerical solution used in this paper and the analyt;ical solution above are in good agreement. Post Ablation In the post ablation region, we are interested in the variation of the regression rate with time and also with the Stefan parameter e. A perturbation solution can be derived on the basis that e is small [1]. The uniformly valid asymptotic expansion of the Stefan condition (9) is then, d~

2[ 1

+((~erfc(7w2)--(rrT)-;e-~ )][arcsm ( ( 1 1

which has errors of 0(~) when

1

]-

1

lr

.

T = O(e 2) and of O(e 2) when

¼rre2/r)½)]

w = 0(i)

.

In F i g . l , the normalized regression rate against normalized time is plotted for different values of A~, with ~ = 0.2.

This figure shows that for large 5~, the solutions are less

accurate. It also shows that the solutions oscillate for large A~. However, by halving ZX~ successively, we obtain a convergent solution . Fig.2 compares the p e r t u r b a t i o n solution (17) with the numerical scheme. The two solutions are very close to each other when e is small. When e = 0.05 they almost overlap, whereas for e = 0.2, the numerical solution is about 4% higher than the p e r t u r b a t i o n solution, and e = 0.75, the numerical solution is much greater than the p e r t u r b a t i o n solution . The value of the scaling factor r~ has an effect on the behaviour of the post ablation solution. The effect of rr on the error of the post ablation result is analysed (see Fig.3). Our numerical results suggest that the regression rate tends to the limit

d~/dr =

1 as t i m e increases. We define a measure of the total difference between

the regression rate over a t i m e domain and the limiting regression rate as Ilvl[=

(1 -

)2dw

(18)

Note that the regression rate approaches unity when the time approaches infinity. Therefore the limit of the integral in equation (18) was chosen arbitrarily to correspond roughly to the case e = 0.05 where the regression rate is reasonably close to 1 at r = 0.,5. Since the solution converges, improved spatial discretization ( that is, increased number of mesh points ), improves it, and the difference between two solutions of different spatial discretization is a

186

C.H. Leung et al.

Vol. 24, No. 2

l O1 ~" 0.8

= 0.6

"it 0.4

0.2

I::11

•

I:.'// Illl//

O.Ol " 0.0

-" ---- --

//

" I 0.1

t 0.2

~ 0.3

I 0.4

dimensionless

I 0.5

A*

= o.o2s : 0.0125

A~ = 0.0015625

I 0.6

t 0.7

I 0.8

I 0.9

time ('r)

FIG. 1 Comparison in the accuracy of regression rate-time relationship with different A~ values

q u a n t i t a t i v e measure of the convergence. This can be made specific by calculating the value of

II v II for

two different meshes and taking the difference between them. It should be noted

that this integral only measures the difference in solutions up to the time level of r = 0.5, instead of the complete solution. For this reason, we chose a particular solution where the regression rate had reached a value reasonably close to 1 within a finite time span . By repeating this process for a range of a values, we obtain a range of II v II values for the same e. Repeating the above process for different values of e results in a different convergence curve for different e. In the analysis, the mesh sizes chosen were 400 and 450, at which converged solutions were observed. The result is plotted in Fig.3 for different values of ~. In the figure, the absolute difference between norm 1 (11 v II evaluated with mesh size 400) and norm 2 (11 v II evaluated with mesh size 450) is calculated and is plotted on a logarithmic scale against the scaling factor ~ . It is interesting to note that, for each value of e, there is an optimum cr such that the

Vol. 24, No. 2

HEAT TRANSFER TO AN INFINITE SOLID

187

1.0-

~0.8-

.~ 0 . 6 -

g ~0.4:

p e r t (~ = 0.75)

.[

pert (a = 0.05) hum (Z = 0,05) pert (~: = 0.2)

40.2 -

h u m (e = 0.2)

0.0 0.0

I

I

I

I

0.1

0.2

0.3

0.4

0.5

0.6

0.7

O.S

I

I

0.9

1.0

dimensionless time (~-)

FIG. 2 Comparison between numerical and perturbation solution at different e values

rate of convergence of the solution is maximised. For example, from Fig.3,when e = 0.1, the o p t i m u m a is about 1. Also notice that there is a tendency for the optimum value of ~r to decrease as the Stefan parameter e increases.

Summary and Conclusions

T h e numerical solution of the pre-ablation region and the exact solution are in excellent agreement. The order of the average error varies from less than 0.1 to just above 0.0001 depending on the transformation factor er . In the post-ablation region the numerical solution agrees with the perturbation solution for small values of the perturbation variable e. The difference between the two solutions depends on e. With e = 0.2, the difference is about 4%. The numerical solution is assumed to be valid for the whole domain of e o11 the basis of the agreement with the perturbation solution at small e . Provided the constant regression at each time level is sufficiently small, the numerical

188

C.H. Leung et al.

Vol. 24, No. 2

1E-05-

....... ¢=0.I ~~ =0.2

7E-06-

4E-06 -

\\\ \

0.3

2E-06-

2

IE-067E-07 -

4E-07-

2E-07 0.2

J 0.3

I 0.4

t 0.5

I L I 0.6 0.7 0.8 scaling factor (o)

I 0.9

I 1.0

I 1.1

I 1.2

FIG. 3 Variation of norm (ll ~ II) against a for different ~ values

scheme will be well behaved, (that is to say, it converges to a correct value and is stable). Moreover, the rate of convergence can be optimised by choosing a suitable value of the transformation factor a . Further development will involve simulating a more realistic heat transfer problem of a charring material. There has been much interest in seeking the behaviour of the regression rate when a layer of another material, for example, char, is formed during the ablation process. The phase change means t h a t heat input to the original material is shielded from a layer of second material which forms as a product of ablation. Hence, in practice, the original m a t e r i a l receives only part of the constant external heat input. The long term goal is to develop a model which can be used to predict such realistic behaviour under char forming conditions . Nomenclature

Dimensional variables T

t e m p e r a t u r e (K)

Tp

material ablation t e m p e r a t u r e (K)

Vol. 24, No. 2

HEAT TRANSFER TO AN INFINITE SOLID

t

time (s)

y

spatial coordinate

v

limiting velocity (velocity scale) (m/s)

l

length scale (m)

s

position of the moving boundary (m)

F

constant uniform heat flux (W/m 2)

cp

specific heat capacity (J/kgK)

Lv

latent of ablation (J/kg)

z

transformed spatial coordinate (dimensionless)

c~

thermal diffusivity (m2/s)

A

thermal conductivity (W/mK)

p

density (kg/m a)

189

(111_)

Dimensionless variables 0 temperature spatial coordinate r

time

(

position of the moving boundary

e

Stefan parameter

a

scaling factor in corrdinate transformation References

1.

J.G. Andrews and D.R. Atthey, J. Inst. Maths Applics, 1.5, p.59 (1975)

2.

H.G. Landau, Quarterly of Appl Maths, 8, p.81 (1950)

3.

G.D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Clarendon Press, Oxford, 3rd ed. (1985)

4.

R.D. Richtmyer and K.W. Morton, Difference Methods for Initial value Problems, Interscience, New York, 2nd ed. (1967)

5.

H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, 1st ed. (1947) Received August 26, 1996