Quasilinearization solution of the Stefan problem

12%

Shorter Communications

diameter of the cohunn height of the column diameter of the glass spheres flow rate void volume of the investigated system void fraction volume of the injection zone volume of the measuring zone mean residence time in the investigated zone mean residence time in the injection zone mean residence time in the measuring zone

6cm 50 cm imm 1.31 I/hr 564 ml 39.4% 94.5 ml 90.6 ml 30 min 5.02 min 4.82 min

Numerical integration was executed to calculate F,,,(t) by use of a quadratic interpolation polynomial and numerical differentiation to calculate f,&f) by use of a second-order smoothing curve, using intervals of h = 0.01 t/T2in both cases. The influence of the other terms of Eq. (5) is shown in Fig. 4, too. The term 7, . n. fL,(f) sometimes can be neglected, if the values of 7, and ~~are sufficiently small. Then no numerical differentiation is necessary. Technische Universitiit Berlin Znstitutfiir Technische Chemie 1 Berlin 12 Strasse des I7 Juni 122-128, Germany

E. OTTO? W. GESTRICH

REFERENCES

____)L

t [Mid

Fig. 4. Calculation of the RTD cumulative density function of the packed column (dashed line) due to Eq. (5). tNow with CIBA-GEIGY Marienberg GmbH, D.-6140, Bensheim. Germanv.

Chemical

Engineering

Science, 1974, Vol. 29, pp. 1296-1298.

Pergamon

Quasilinearization

Press.

[l] C. Brie, M. Guivarch and J.-L. Petit, Proc. 1st Znt. Conf. Calorimetry and Thermodynamics p. 73 Warschau, Aug.-Sept. 1969. [2] P. C. Gravelle, Advan. Catal. 1972 22 191. [3] 0. Levenspiel and .I. C. R. Turner, Chem. Engng Sci. 1970 25 1605. [4] 0. Levenspiel, B. W. Lai and C. Y. Chatlynne, Chem. Engng Sci. 1970 25 1611. [5] J. C. R. Turner, Chem. Engng Sci. 1971 26 549. [6] D. K. Schmalzer and H. E. Hoelscher, A.I.Ch.E. Jl 1971 17 241. [7] C. G. Sinclair and K. J. McNaughton, Chem. Engng Sci. 1965 20 261. [8] P. ialoudik, Brit. Chem. Engng 1969 14 657. [9] W. Gestrich and E. Otto, Chem. Zng. Techn. 1971 43 1241. [lo] W. Gestrich and E. Otto, Chem. Zng. Techn. 1973 45 577.

Printed

in Great Britain

solution of the Stefan problem

(Received 25 May, 1973; accepted 11 October 1973) There is a large class of important problems which involve the solution of either the diffusion or heat conduction equation with a moving boundary. These problems are referred to generically as the Stefan problem. The presence of a moving boundary imposes an additional nonlinearity on the nonlinear differential equations. This topic has been reviewed by Rubinsteinrl] and Bankoff[2]. In this paper we consider a version of a Stefan-type problem, “the melting problem” in a semi-infinite solid with temperature-dependent thermal properties. Consider a solid below its melting point, T,. At some arbitrary time one surface is brought to and maintained at

a fixed temperature, T, > T,,,. Heat is conducted from the surface through the solid resulting in two phases. The Iiquid phase has a temperature profile bounded by T,, and T, whereas the moving boundary is maintained at T,,,. The solid has a temperature profile bounded by T, and T,. The problem is to determine the respective temperature profiles as well as the position of the moving boundary X(t) as a function of time. By employing the standard Boltzmann and Zener transformations the dynamic equations are transformed into two point boundary-value problems. We employ an additional error function transformation in order to reduce the domain of the resulting

1297

Shorter Communications system of equations from- 1 to 1 and then apply the quasilinearization technique to solve the resulting nonlinear ordinary differential equations. Hamill and Bankoff 131have previously considered the plane melting problem and have given a time-independent similarity solution for it. We defer to this paper for a statement of the problem. After Boltzmann and Zener transformation of the liquid and solid phase heat conduction equations and accompanying boundary conditions the problem reduces to a consideration of the following cases. A. CASE

a, = I,

3

a, = I + IOU,, a*= t+ IO”.

00 00

where I is an arbitrary constant, erf is the error function, and y =x/X(t). By this additional transformation the semi-infinite domain of the system of equations is reduced to the domain - 1 to + 1. This error function transformation substantially decreases the size of the mesh, the necessary computer storage space, and the required computation time. Additionally, we are able to overcome the difficulty of determining the finite distance, yO, at which the reduced temperature vanishes. As we need not know the exact value of y0 it can be estimated. The value of 1 is chosen to make erf f (yO- 1) - 1 (cf. Ref. 4). By use of the quasilinearization technique[4] the equations and boundary conditions may then be linearized and solved by the superposition method[5]. Sets of particular and homogeneous solutions were obtained by numerical integration of the resulting linearized iterative differential equation (complete details may be obtained from the authors). The integration proceeds outward in both directions from z = 0 and the Runge-Kutta fourth-order method was employed for this purpose.

2.0

AND

2m

29 NO. 5-P

Fig. 1. Temperature-distance curve for a = I+ au + bu*. E=l, m=l, K=l, e=O, 0=0.3, T-CT,,,, and p,= p*= 1.

o.,.l 0.6

@aa,

’

I-02u,

i

3

0.0

a*.----

I-O 5th ,&--_ l-0 7u,

IO

2.0

25

LX= curve for Temperature-distance Fig. 2. l/(l+au+bU*). E=l, m=l, K=l, c=O, B= 0.3, T,
@ o=O. a

0~2.5,

b=O b=O

DISCUSSION

If the thermal parameters are independent of temperature the exact solutions for Case 1 and Case 2 are both known[6,7]. We have compared the quasilinearization solution with the exact solutions for these two cases. The difference between them was found to be less than 1 X 10e5 over the entire range of reduced temperature. To illustrate the use of the quasilinearization scheme with the nonlinear Stefan’s problem (temperature!ependent thermal parameters) we considered the following representative cases for which the thermal conductivities ki depend on temperature (k, = kiOoi(ui): ai = 1 + a,U# +b,UF, l/(l+a,Ui+biUiz), and exp(aiU,), where ai and bi are arbitrary constants. Numerical results for these cases are shown in Figs. 1-3. The results are plotted as

CES VOL.

3.0

x

2

initial temperature of the solid, T,, is equal to the transition temperature, T,,. The resulting system of equations and boundary conditions are considerably simplified over Case 1 and the domain of the equations goes from 0 to 1. Hence the error function transformation need not be employed for this case. However, we define the variable .z as z = 1 - y and then use the quasilinearization technique to linearize the resulting system and boundary conditions. The superposition method was also used in a manner similar to that of Case 1. The

RESULTS

I+5u,,a*=I+5u.

IO

Define: z = erf I(y - 1)

NUMERICAL

I

a,=

1

initial temperature of the solid T,, is less than the transition temperature, T,,. The

B. CASE

a*=

a

0

Fig. 3. Temperature-distance curve for a = m=l,A=l,e=O,f3=O,T,=T,,andp,=l.

1+

au + bu2.

Shorter Communications

1298

NOTATION

constant, coefficient of polynomial in U constant, coefficient of polynomial in U specific heat index, 1 = liquid phase, 2 = solid phase thermal conductivity parameter in error function order of iteration temperature of liquid phase temperature of solid phase melting point of solid constant surface temperature uniform temperature of solid phase time reduced temperature, defined in Eq. (9) position of moving boundary Zener’s variable, defined in Eq. (9) distance at which reduced temperature vanishes

OE.0 Q E = -0.5 @E

= -1.0

> 04

02

0.0 00

IO

20

3.0

x 2a

Fig. 4. Temperature-distance curve for various e(= lE = 1, m = 1, K = 1, 0 = 0.3, T, < T,,,, a = e”‘l, a = es?, and p, = 1.

(&p,)).

Uu. x/(2%‘/h2t) in a manner similar to Ref. 4 in which the dependence of reduced concentration on the analogous Boltzmann variable x/2vD0t) was displayed. In order to ascertain the effect of the relative density of the two phases on the U vs. 7 profile we display this dependence for the case a = eO” in Fig. 4. It appears that the effect is small for large differences in density. In Figs 1,2 the curves are given for a specified value of tl = (T, - T-)&T, - T_). For two phases the position of the moving boundary y = X(t)/Zvh,t where A2= k2Jp2Cp2,, depends on 0 in contrast to a one-phase pure liquid system. y is given by the inflection point in the curve. Liquid phase exists for data points above the inflection point and solid phase below. Although we have ascertained the temperature profiles for single-phase liquid systems having various functional forms for the thermal conductivity we illustrate the type of result obtained in Fig. 3 for the case cr = 1 + au + bu*. The moving boundary y is located at U = 0 for these cases. The convergence rate of this method is discussed in detail in Ref. [4]. Depending on the functional form of the temperature dependence only four to nine iterations are required to obtain five-digit accuracy by Runge-Kutta shooting. The basic advantage of this technique is, in fact, that a recursive solution of the linearized equations has the property of quadratic convergence. The quasilinearization scheme appears to be a very effective tool for solving nonlinear Stefan-type problems.

Chemical

Engineering

Science, 1974, Vol. 29, pp. 1298-1300.

Pergamon

Press.

Y 8 A P

,‘,‘r function, defined in Eq. (16) functional form of temperature dependence mal conductivity parameter in Zener variable constant reduced temperature constant density

School of Chemical Engineering and Materials Science University of Oklahoma Norman, Oklahoma 73069, U.S.A.

of ther-

CHING-SE

Department of Chemistry RUSSELL Texas A&M University College Station, Texas 77843, U.S.A.

YEN

D. LARSEN

REFERENCES

Rubinstein L. I., The Stefan Problem, American Mathematical Society, Providence, R. I. 1971. [2] Bankoff S. G., in Advances in Chemical Engineering, (Eds. T. B. Drew et al.), Vol. 5 p. 75. Academic Press, New York, N.Y. 1949. [3] Hamill T. D. and Bankoff S. G., Ind. Engng Chem.

[I]

Fundls 1964 3 177. [4] Yen C. S. and Larsen R. D., J. Comput. Phys. 1972 10

5.54. [5] Lee E. S., Quasilinearization and Invariant Imbedding. Academic Press, New York, N.Y. 1968. [6] Zener, C., J. Appl. Phys. 1949 20 950. [7] Kehoe, P., Chem. Engng Sci. 1971 26 809.

Printed

Inlet and outlet residence-time

in Great Britain

distributions

(Received 4 October 1973)

In reactor design problems we are often concerned with residence-time distribution of material leaving the reactor, since a knowledge of this distribution may tell us much

about the performance of the reactor. In general, we look at the material coming out of the equipment at any instant and wish to know at what times, stretching backwards to