APPROXIMATE

METHOD

FOR

TEMPERATURE

FIELD

WITH ACCOUNT

OF THERMAL

V.

V.

Salomatov

CALCULATING

IN A G R O W I N G

a n d A.

THE

CRYSTAL

RADIATION D.

Gorbunov

UDC 536.212.2

An engineering method is described for solving the heat-conduction problem in a growing single c r y s t a l for the case in which heat is lost at the surface through radiation. An integral substitution is used to t r a n s f e r the nonlinearity in the boundary condition into the basic equation for the p r o c e s s ; the nonlinear complex in the t r a n s f o r m e d equation can be evaluated. The approximate solution is c o m p a r e d with the results of a c o m p u t e r numerical calculation. Because of the increasingly stringent requirements on c r y s t a l quality, considerable attention is being devoted at p r e s e n t to the development of analytic methods for calculating the transient t e m p e r a t u r e fields existing during crystallization from a hot melt. The physical p r o p e r t i e s and the quality of the c r y s t a l s grown are known to depend to a large extent on the thermal conditions. Knowledge of the t e m p e r a t u r e fields in c r y s t a l s during their growth will allow a determination of the basic c r i t e r i a affecting the growth conditions and the s t r u c t u r e and will allow a study of optimum methods for controlling crystallization. The problem with which we deal here - heat conduction in a growing c r y s t a l - belongs to the class of m a t h e m a t i c a l - p h y s i c s problems involving an important nonlinearity in the boundary conditions (the heat dissipation is proportional to T 4 and there is a moving boundary). Such problems have usually been solved by various approximate analytic and numerical methods for linear conditions of convective heat t r a n s f e r . Let us examine the basic studies which have appeared on this subject. 1.

Studies

with

Convective

Boundary

Conditions

The s t e a d y - s t a t e two-dimensional problem was solved in [1-4], and a model of a constant heat input through the crystallization front was proposed [4]. Some studies [2, 3] have taken the finite c r y s t a l size into account. The nonsteady-state problem has been solved only in [5], and the growth rate was taken into account, but here a semiinfinite cylinder was a s s u m e d . 2.

Boundary

Radiation

Only the simplest, p r i m a r i l y s t e a d y - s t a t e , problems of this type have been solved in closed form for both finite and infinite regions. Interestingly, a closed solution can be found when there is a certain t e m perature dependence for the t r a n s f e r coefficients; in [6], for example, the t h e r m a l conductivity was p r o portional to l / T , while in [7, 8] it was proportional to T 3. The radiation law was linearized through the use of a convective coefficient in the approximate solutions in [3]. A numerical method was used in [9] to solve the two-dimensional s t e a d y - s t a t e problem involving both radiation and convection and for the case of a variable t h e r m a l conductivity. The v a r i e t y of models which have been used led to a refinement [10] of the mathematical formulation of the heat-conduction problem in a growing c r y s t a l . S. M. Kirov T o m s k Polytechnical Institute. Translated from Izvestiya Vysshikh Uchebnykh Zavedenii Fizika, No. 2, pp. 89--95, F e b r u a r y , 1971. Original article submitted F e b r u a r y 16, 1970.

9 1973 Consultants Bureau, a division of Plenum Publishing Corporation, 227 West 17th Street, New York, N. Y. 10011. All rights reserved. This article cannot be reproduced for any purpose whatsoever without permission of the publisher. A copy of this article is available from the publisher for $15.00.

212

This b r i e f review has shown that it is not possible, on the b a s i s of the existing solutions, to propose a method for the calculation of the t e m p e r a t u r e field which s a t i s f i e s the engineering r e q u i r e m e n t s both of s i m p l i c i t y and of sufficient a c c u r a c y . Below we witl analyze the following m a t h e m a t i c a l model, a s s u m i n g a semiinfinite cylindrical ingot: the heat-conduction equation

aO(R,Z, Fo) 6Fo

O2O(R,Z, Po) + 1 OO(R,Z, Fo) ~.O'-'O(R,Z, OR2 R OR Og~

where RE [0, 1], ZE [0, o~], FoE [0, oo], and OE (|

Fo)

p e a O ( R , Z, Fo),

(1)

c)Z

1], with the initial condition

O (R, Z, 0) = ~ (R, Z)

(2)

and the boundary conditions O0 (1, Z, Fo)

OR dO (0, Z, Fo)

OR

--

Sk [O' (1, Z, Fo) - - O'm],

O,

OO (R, ~ , Fo)

az

=

(3)

O,

(4)

o (R, o, Fo) = ~ ( O , Fo).

(5)

H e r e R = r / R 0, Z = z / R 0, S k = a b T ~ R 0 / h , P e = VRo/a, | 0, | = T m / T 0 , To, a n d T m a r e t h e c r y s t a l lization t e m p e r a t u r e and the t e m p e r a t u r e of the surrounding medium, and ~b(R, Z) and ~b(R, Fo) a r e any differentiable functions. Approximate

Analytic

Method

of Calculation

T o reduce b o u n d a r y - v a l u e p r o b l e m (1)-(5) to a function equation, we use s u c c e s s i v e l y finite Hankel integral t r a n s f o r m a t i o n s with r e s p e c t to coordinate R with the kernel 1

g(~, z, Fo) = {O (~, z, Fo)}. = S O (t~, Z, Fo) R~o(~R) dR, 0

(6)

9

where ~z a r e the roots of II(/D = O, and we use a Laplace t i m e t r a n s f o r m a t i o n with the kernel co

O, z , s) = {{e (R, z , Fo)}~}~ = ~ (~, Z, Fo) e-sFodFo.

(7)

0

Solving the resulting s y s t e m of equations for the t r a n s f o r m s , and c a r r y i n g out the i n v e r s e t r a n s f o r mations, we find y,

ooI".~(P.,z)- dZ} -F2VpZ+p2+s~zy~(~,Z, ~yt(v.,Z, s)

z , s)

~ Io(,R) I

-

(.,

y,,(r~,Z,s)

dz _

[ cz-. r

- Sk

bt2(,,7, ; J +VP'~P+

s.:~y, O, z, s) j , - '

i8)

where 7(#, Z) = {~b(R, Z)}H is the t r a n s f o r m e d initial condition; ~{/~, s) = {{goOR, FO)}H}L is a twice d i d ferentiable boundary condition; Yl,2 = e ( P ~ ~ Z ; indicates the i n v e r s e Laplace t r a n s f o r m a t i o n .

P = P e / 2 , ~--(Z, s) = {[|

Z, F o ) - |

and L -1

Equation (8) f o r m a l l y gives the t e m p e r a t u r e at any point of the growing c r y s t a l in t e r m s of the s u r face t e m p e r a t u r e . To find 0(1, Z, Fo) we s e t R = 1 in (8); then f r o m the resulting e x p r e s s i o n , which is a nonlinear V o l t e r r a integral equation of the second kind, we can find the s u r f a c e t e m p e r a t u r e by the method of s u c c e s s i v e a p p r o x i m a t i o n s . The convergence of this solution was d e m o n s t r a t e d in [11] for the onedimensional mode l. Analyzing (8), we conclude that for s m a l l values of Sk, c o r r e s p o n d i n g to the usual conditions holding during c r y s t a l growth, we can neglect the second t e r m in (8). Then, for the m a x i m u m and simple values ~(R, Z) = 1 and ~(R, Fo) = 1, we find 213

T A B L E 1. C o m p a r i s o n of the R e l a t i v e T e m p e r a t u r e s C a l c u l a t e d f r o m Eqs. (11) and (13) and T h o s e Found in the C o m p u t e r Calculation for Sk = 0.1, Pe = 0, and | = 0

Fo• R

I

0,4 C0mp. Calc.

0,9969 0,9692 0,9690 0,9429 D,9515 3,9278 3,918~ 3,8987 3,9121 9,8929 0,9035 0,885~ 0,903~ 0,8856

0,1 0,5 1,0 5,0 10,0 50.0 100,0

1,0 ~% Comp. Calc.

0,9981 0,9654 0,9717 0,9390 0,9517 0,9216 0,8986 0,8749 0,8835 0.8616 0,8775 0,8561 10 8775 0,8561

--0,1 +0,3 --0,3 +0,4 --0,2 -}-0,7 ,},2,1 -}-2,6 +3,1 +3,5 +2,9 +3,2 +2,9 +3,2

0.9959 0,9624 0,9416 0,9111 0,9018 0,8763 0,8212 0,8038 0,8046 0,7885 0,7906 0,7748 0,7906 0,7748

0,9974 0,9618 0,9471 0,9112 0,9061 0,8757 9,8017 9,7828 0,7738 0,7574 0,7625 0,7470 0,7625 0,7470

5,0 2,0 I ~% Comp. Calc. ~% Comp. Calc. --0,1 0,9959 +0,1 I 0,9621 --0, 5 I 0,9304 --0,9 ] O,8998 --0, 51 O,8701 +0, I I 0,8465 -}-2,31 0,7222 +2,5 0,7111 +3,8 0,6874 -}-3, ~ I O,6784 -}-3,5 I 0,6668 -}-3,7 I 0,6588 +3,5 ] 0,6668 +3,7 ] o, 6588

0,9973 0,9617 0,9367 0,9016

I[--0, 1 I +0,1 --0,7 --0,1

0,9959 0 9973 0,9621 0,9617 0,9296 0,9359 0,8989 0,9008 0,8779 --0,9 0,8626[0,8712 0,8504 --0.5 9 8394 0,8444 0,7168 -}-}-0,7 0,6268 0,6454 0,7O42 +l,0 9,6201 0,6370 0.6724 -}-2,1 0.5331 0,5296 0,5590 0,6626 -}-2,3 0,5543 0,6538 ,},1,8 0, 867 o,5o97 0,6451 ,},2,1 0,4853 0,5054 0,6538 ,},1,8 o,4867 i0,5097 0,6451 i ,},2,1 0,4853[0,5064

10,0 ~% Comp. Calc. +0,11

--0,

,.

~%

.

--0,7 --0,2 --1,0 --0.6 --2,3 --2, 1 -3,3 --3,2 4,7104007 I0,4133 --3,2 -4,5 I0,4006 I0,4118 --3,1 -4,510,4006 I0,4,32 --3,2 --4,5 i 0,4005 [0,4117 --3,1

--0,7 I 0,9296 0,9359 --0,2 ] 0,8989 0,9008 --1,010,8626 0,8712 --0,6]0,8394 ]0,8444 -2,9]0,623010,6379 --2,7 0,6163 10,6299 --4,810,5124 0,5296 --4,610,5093 ]0,5257

O ( R, Z, Fo) = ~'~I~(~R---~)[ epZ ( E (~,, Z, VF-~) -- e-,~'voE (a~.,Z, V'F-~)) (9) whe re

cr

p 2,

a.z=p,E(~x,z; V-F-o)=

1[ e - ~ z - e r f C /~ 2 vz~ ~

which c a n be used a s a f i r s t a p p r o x i m a t i o n for c a l c u l a t i n g the t e m p e r a t u r e Engineering

Method

field in the i n t e r v a l 0 < Sk < 0.05.

of Calculation

I n t r o d u c i n g the i n t e g r a l a n a l o g of the t e m p e r a t u r e ,

~(R,Z, Fo)=

C.exp

(10)

04(R,Z, Fo)_Om4 ,

w h e r e C = const, we can c o n v e r t s y s t e m (1)-(5) to

O~(R,Z, Fo) OFo

O~',I(R,Z,Fo) +I O~(R,Z, Fo) + O%~(R,Z, Fo) p e 0 ~ ( R , Z , F o) ~-f(R, Z, Fo), OR" R OR OZ2 OZ (R, Z, 0) = C. exp ( ~ dO/iO 4 (R, Z, 0) - One]) = ~o (R, Z), 0~(1, Z, Fo) _. oR 0~ (0, Z, Fo) - - 0 ,

OR

Sk~(1, Z, Fo),

(1,) (2')

(3')

O~(R' co' F~ = 0,

(4')

OZ

(5')

~q(R, O, Fo) = C. exp ( ~ dO / [O 4 (R, 0, F o ) - Ore4]), where 403 (R, Z, F~ - 1

[ (0~ (RO_~, F~

+

is the n o n l i n e a r c o m p l e x and w h e r e we have 0 -< ~(R, Z, Fo) - ~70. The n o n l i n e a r t e r m fiR, Z, Fo) m a k e s it difficult to find a g e n e r a l c l o s e d solution of s y s t e m (1')-(5'); we will i n s t e a d c o n s t r u c t an a p p r o x i m a t e solution. Setting fiR, Z, Fo) = 0 with | Z, 0) = | 0, Fo) = 1, and a p p l y i n g (6) and (7) to the r e s u l t i n g s y s t e m of equations, we find a f i r s t a p p r o x i m a t i o n for ~?(R, Z, Fo): t

214

~, (R, Z, Fo) = ~o Y, AJo p.

fiR) [epz (E (~,, Z, ]fP-~) -- e-~~176(~,,, Z, V 1%)) +

e-~ 'Fo] ,

(11)

w h e r e A # = 2Sk/Io(/~)[Sk 2 + #2] i s the i n i t i a l t h e r m a l a m p l i t u d e , and # a r e the r o o t s of the e q u a t i o n p i t ( p )

= Ski0(#). E q u a t i o n (11) with (10), w h i c h c a n be r e p r e s e n t e d g r a p h i c a l l y w i t h o u t a n y s e r i o u s l o s s of a c c u r a c y , c a n be u s e d to d e t e r m i n e the d e s i r e d t e m p e r a t u r e | Z, Fo) in the f i r s t a p p r o x i m a t i o n . T a b l e 1 s h o w s the d i m e n s i o n l e s s t e m p e r a t u r e s a t the c e n t e r and s u r f a c e of the g r o w i n g c r y s t a l c a l c u l a t e d by a n u m e r i c a l m e t h o d on a n M - 2 0 c o m p u t e r and c a l c u l a t e d in t h e f i r s t a p p r o x i m a t i o n b y o u r p r o c e d u r e . The c a l c u l a t i o n s w e r e c a r r i e d out f o r Sk = 0.1, P e = 0, and | = 0. A n a l y s i s of the t e m p e r a t u r e field s h o w s t h a t with Sk -- 0.1; 0 < P e -< 0.15; and 0 --- O m -< 0.4 the e r r o r s a r e e v e n s m a l l e r t h a n t h o s e q u o t e d in T a b l e 1. T h i s s i t u a t i o n a r i s e s b e c a u s e i n c r e a s e s in Pe and O m a r e to s o m e e x t e n t e q u i v a l e n t to a d e c r e a s e in Sk. T o c a l c u l a t e the t e m p e r a t u r e s for h i g h e r v a l u e s of Sk, we c o n s t r u c t the s e c o n d a p p r o x i m a t i o n , s o l v i n g s y s t e m ( 1 ' ) - ( 5 ' ) , t a k i n g a c c o u n t of the c o m p l e x f(R, Z, Fo) c a l c u l a t e d f r o m the f i r s t a p p r o x i m a t i o n . We f i r s t d e t e r m i n e the f u n c t i o n a l d e p e n d e n c e

0 (R, Z, Fo) = Expanding 5 d@/[|

Z, Y o ) - |

O[t

(R, Z, Fo)].

from (10) in a series for the selected value C = exp (~/4|

(R,Z, Fo)=exp

- - ~ m~ 3 e a ( R , Z , F o ) F 707 (R, Z, Fo) . . . .

and r e t a i n i n g o n l y the f i r s t t e r m in the s e r i e s - a l e g i t i m a t e p r o c e d u r e f o r s m a l l | an explicit expression for | Z, Fo)

(up to 0.4) - we find

0 ( R , Z, Fo) = [ - - 3 1 n ~(R, Z, Fo)] - ' a .

(13)

U s i n g (13) and (1t), we find A

(R, z, Fo) = f(~,)

=

4[-31n,h (te.Z, Fo)l-~ -- I [ (Or. (R,Z, Fo)) a _I_(O~, (R,Z, Fo))'2 ] [--3 In -q, (R, Z, Fo)]-'o -- (@r~ OR oa _

(1r

and then the s e c o n d a p p r o x i m a t i o n :

(R,Z,

Fo) = ~q (R, Z, Fo)

-FEA~.Io(I~R)

y~ (~, Z, s) 21/p~+ ~"+s

Yl (P, Z, s) 2 1/t)~ -? F2 + s ,

I. z s)

I

(15)

Y l ( ~ , Z , s) ] c - "

where ~t(/~, Z s) = {{ft(R, Z FO)}H}L. A n a l o g o u s l y , we c a n w r i t e the n - t h a p p r o x i m a t i o n :

~(R,Z, Eo)=~_,~R,Z, Fo) + 2A,~.Io(,~R){ 9

o

,~

y2(~,Z,s) 2]/-p~ +?2~_

s

2Vp+/'4-s~ y,(~,z,~)'j~_,

The convergence of this approximation procedure method was demonstrated in [12] for f = f(Z, Fo). C o m p a r i s o n of t h e s e r e s u l t s with the c o m p u t e r r e s u l t s s h o w s t h a t the f i r s t a p p r o x i m a t i o n is good w i t h i n 10% f o r Sk --- 0.25 and t h e s e c o n d is good w i t h i n the s a m e e r r o r for 0.25 --< Sk -< 1.5. A n a l y s i s of the c o m p u t e r r e s u l t s a l s o s h o w s the following: 1. The e r r o r i n t r o d u c e d b y n e g l e c t i n g the r a t e of c r y s t a l g r o w t h d o e s not e x c e e d 2% f o r P e -< 0.15 (in p r a c t i c e the l i m i t i n g value) for SI~ < 0.1.

215

2. F o r s m a l l values of the radiation c r i t e r i o n (Sk -< 0.15), c o r r e s p o n d i n g to r e a l situations, we can a s s u m e the t e m p e r a t u r e field of the c r y s t a l to be uniform within 4% in the interval 0 < Pe -< 0.2. CONCLUSIONS Even in the f i r s t approximation, this method gives quite a c c u r a t e l y the t e m p e r a t u r e field in a growing single c r y s t a l for the p a r a m e t e r ranges Sk -< 0.1, Pe -< 0.15, and | ~< 0.4, which c o v e r the working r a n g e s for c r y s t a l l i z a t i o n f r o m m e l t s . The r e s u l t s of a n u m e r i c a l calculation have been used to evaluate the e r r o r s a r i s i n g due to the possible neglect of radial gradients and the motion of the growing ingot. LITERATURE 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.

216

CITED

E. Billig, P r o c . Roy. Soc., 235A, No. 1200 (1956). G . P . Boikov and V. D. Kuchin, Izv. VUZ. Fiz., No. 2 (1959). Ya. T s ' u i and F. T s ' o u , T e p l o p e r e d a c h a , No. 3, 116 (1963). W . R . Wilcox and R. Duty, ibid., 88, No. 1 (1966). A . A . Uglov, in: P r o c e e d i n g s of the Moscow Institute of Railroad T r a n s p o r t a t i o n E n g i n e e r s [in Russian], No. 189 (1965!. E. Billig, P r o c . Roy. Soc., 229A, 346 (1955). V . V . Salomatov, Izv. VUZ. Fiz., No. 5 (1965). V . V . Salomatov, ibid., No. 1 (1966). K. Akiyama and J . Yamaguchi, J . Appl. Phys., 33, No. 5 (1962). L . A . Goryainov, in: P r o c e e d i n g s of the Moscow Institute of R a i l r o a d T r a n s p o r t a t i o n E n g i n e e r s [in Russian], No. 254 (1967). V . V . Salomatov, I n z h . - F i z . Zh., 17, No. 2 (1969). V . V . Salomatov and A. D. Gorbunov, Izv. VUZ. Energetika, No. 3 (1969). -

-

j