Mathematical modelling and control system design of Czochralski and liquid encapsulated Czochralski processes: the basic low order mathematical model

j. . . . . . . .

ELSEVIER

CRYSTAL GROWT

H

Journal of Crystal Growth 154 (1995) 172-188

Mathematical modelling and control system design of Czochralski and liquid encapsulated Czochralski processes: the basic low order mathematical model G.A. Satunkin Centralnaja 18-55, 142432 Chernogolot,ka, Moscow District, Russian Federation

Received 16 May 1994; manuscript received in final form 31 March 1995

Abstract

The low order model (LOM) taking into account a simple form of curvature of the crystallization front is proposed for Czochralski (Cz) and liquid encapsulated Czochralski (LEC) methods. The scope of the LOM application is determined.

1. I n t r o d u c t i o n

The application of a modern automatic control theory to the growing processes is largely suppressed by the absence of generally accepted mathematical models of the controlled object. An interesting " l u m p e d " model of Czochralski (Cz) and liquid encapsulated Czochralski (LEC) methods are given by Gevelber and coworkers [1,2], which allows one, in principle, to propose the controller of the complex structure, when there are sufficient system p a r a m e t e r sensors. The TV system of precise temperature measurements, presented by Wargo and Witt [3], permits the advanced control of the thermal environment of the crystal in a puller with several independent heaters. However, the main heat action on the growing crystal is realized now in Cz and LEC methods with the aid of heating of the melt by a single heater. The output weighing control is the source of information about the current state of the

system. When the behaviour of impurities during the crystallization is not taken into account, the simplest mathematical low order model (LOM) for Cz and LEC methods has to be built by use of three conservation laws: (i) the heat balance at the crystallization front, (ii) the mass conservation, and (iii) the constancy of the growth angle e 0 with linearised response approach. This approach is used by Hurle et al. [4], and is nearest in spirit to this article. The simplest linear model of Cz and LEC methods is proposed here and the scope of this model is determined. For more clarity, the initial data from Hurle and coworkers [4,5] are used in calculations of the transfer processes in this article. The proposed model was tested in design of digital diameter control systems for the growth of LiNbO 3, and was successfully applied for industrial production of this crystal, and the elaboration of special weighing sensors [6]. In the following publications we will discuss the task of design of optimal diameter controllers including the state controllers, the ex-

0022-0248/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0022-0248(95)00050-X

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

173

(3) The melt temperature response at the meniscus base T0(0,t) (Fig. 1) is the main perturbation of the stationary growth. (4) We neglect temperature expansion and evaporation of the melt and incapsulant.

perimental determination of the controlled object parameters and filtering of the measurements noise.

2. Adopted approximations (1) The melt/solid and solid/gas interfaces are considered in an isotropic approximation. The front coincides with the crystallization isotherm except for a very small region 10-100 mK [7] near the three-junctions line (Fig. 1). Any perturbation of stationary growth leads to the local curvature change, and the macroscopic front curvature R[-t does not change, at least during the considered time of the transfer process. (2) The crystallization front Z = S ( r , t ) has the simple form (convex or concave, Fig. 1) with the small macroscopic curvature S'r ---> O. The temperature gradients in solid G s = T~z and liquid G l = T{~ phases make the main contribution to the heat balance at the front -,~sVnTs = - l~lVnT ' + p y c L

3. Thermal analysis It is sufficient to examine the thermal conditions in the meniscus and in a crystal part near the crystallization front for investigation of the transfer processes (Fig. 1). It is necessary to point that the localization of the current crystal radius variation 3 r ( t ) f o r b i d s o n e to use a simple differentiation on r the solution of the stationary thermal problem. It is shown below that the value of the temperature gradient variation at the crystal front G~,, given by Hurle et al. [4], is a factor of four smaller than found in the present work. Namely, the localization of the crystal radius variation allows one to propose the pathway of heat task solution, based on the appearance of a small parameter.

(1)

and Eq. (1) is always valid.

~/

L Zo

C~s ~

;-"

2

R

O'

,

r--

I:z Fig. 1. Scheme for low order model of Cz and LEC methods and the symbols, adopted in the text. Analysis of the heat transfer is performed for the region defined by the dashed line. The front curvature R~-1 (convex or concave) does not change during transfer processes. The curvature near the triphase line is determined by the G i b b s - T o m p s o n effect and is always in equilibrium, so the sum of the three vectors of the surface free energy is equal to zero ~so + ~SL + YLG = 0.

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

174

Let the crystal part with length l 0 near the crystallization front has grown u n d e r the disturbance of stationary growth conditions. We shall insert the dimensionless variables u s = (T s T a ) / ( T c - Ta), u I = (T I - T a ) / ( T c - Ta), x = r / r o , y = z / l o, ~ = z / h o, v = V J V o using as scale units r 0, the crystal radius, h 0, the meniscus height u n d e r stationary crystallization with constant pulling rate V0, Tc = const, and T~ are the crystallization and the ambient temperatures. In general, if the quasi-stationary heat conditions are fulfilled, it is necessary to solve the heat-transfer equation for each phase i = s,1. Using the inserted dimensionless variables, this equation for the solid phase with the coordinates at the crystallization front may be rewritten as: 1~ Ous

vo am

- - - - + - - - - V

XsZs ats* ~E2 s

Xs lo

Oy

1 0 / x O U s ] + _ _a2us • x ~xk ~-x] 0y 2

'

(2)

and for the melt with the coordinates at the meniscus base, as h 2 Ou l

V0 +----V--

XlZl 0tl*

Xl h0

O#

1 o {xOUl~ o2m I ) --. lxTx Tx + a;

= E2 -

(3)

H e r e e s = l o / r o and E! = h / r o are small parameters and h ( t ~ ) ~ h o is a current meniscus height, ~s and ~'l are the characteristic time scales of thermal relaxation processes in the solid and liquid phases. W e shall discuss the thermal relaxation process on the time scale of the crystallization transfer process. In other words, as we wish to study the transfer processes during crystal growth, we chose the solidification time scale as more suitable for our calculations [8]. This time r is determined by the dimensionless variables from Eq. (1):

U'sy - ( A l l o / Z s h o ) u ' l = z ( V o / l o ) v ,

(4)

as r = 1 2 L / x C ( T c Ta), where X = A J p s C is a thermal diffusivity and C the heat capacity of the solid phase. Simple evaluation shows that thermal

processes can be discussed in the quasi-stationary approximation. If we consider z s = ~'l = z in Eqs. (2) and (3), so 12/xs'r = C ( T c - T a ) / L = ETs and h~/xjz--evshoXJloX I = e I. T h e contribution of the t i m e - d e p e n d e n t c o m p o n e n t s must be determined by comparison of the coefficients of the right and left parts of Eqs. (2) and (3). The quasi-stationary condition in the solid phase is correct at ETs < e2s or l o > rOeTsl/2. Since ETsl/2 < 0.2 for (To - Ta) < 5 ° C and for the thermal constants of Ge, Si, and GaAs, the evaluation shows that the quasi-stationary a p p r o a c h is correct after growth l 0 > 1 - 2 m m for the crystal radius r 0 = 1 5 cm (see Fig. 4a). The axial rate of melt motion in the meniscus is equal to the crystallization rate Vc in the simplest approximation. More strictly, the p a r a m e t e r 1 2 V o v / X J o = esPe~(t) ~ EiPel(t), where Pei(t) = r o V c ( t ) / X i is a Peclet number, d e p e n d i n g upon the current time. The order of Pei(t), which determines the convective c o m p o n e n t of heattransfer, may be presented as Pe i = O(1) at the " b e g i n n i n g " of the transfer processes when [~V~(t)/Vco[ > 1 and Pei = es(1) is less at the "ending" of the transfer processes when [~Vc(t)/V~o[ < 1 (see Figs. 4 and 7). U n d e r small deviations from the initial state, the new t e m p e r a t u r e gradients G s and G l with accuracy up to second order in the deviations are presented a s G i = Gio + G~r~r + G~hSh + O(E), i = s,l. T h e current velocity of crystallization is V~ = Vc0 + ~Vc, where Gio is the gradient and Vc0 = V0 - H 0 the crystallization rate in the initial stationary state, and /40 = - p s r 2 V o / ( P l R 2 psr 2) is the initial fall rate of the bulk melt in a cylindrical crucible with radius R. In general, the crystallization rate is given by Vc = Vo dS(r,t)/dt = V o - S~ - S'rr;. Since S ( r , t ) = +H(t)+h(r,t) and S ' ~ 0 , we get V ~ = V g / 4 ( t ) - h ( t ) = V~0 + 3 V ~ ( t ) o r 8Vc = ~ V 0 - 8 H 8/~ for the crystallization rate variation. T h e heat

balance equation (1) after the substitutions is described by the expression:

I~Sr + Ih~h = PsL~Vc.

(5)

H e r e lj = AjG'U - hsG~j , j = r,h. Analysis of transfer processes in the controlled object requires explicit functions for G~j.

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

3.1. The thermal conditions in the melt

=

a2uj

1 a [xaUl~

aul / ) - e ( y ~ - =0. ~-x

-~- O(E~). The width of the boundary layer can be evaluated by the values O(el) or O(ho). So the temperature distribution in the meniscus is one-dimensional with small correction, appreciable near the lateral meniscus surface only. One can express from Eq. (10) the temperature gradient in the melt as Gl(r,z) = (T c - T ( O ) ) / h + O(e 0 so the term G'lr = 0 and EltO0(p,~" ) "J- ~~(.Ol(p,~') -[- "'"

To determine function G~ and G~h in Eq. (5), we shall take into consideration some geometrical peculiarities of the meniscus. For crystals with initial radius r 0 >> h 0 --* 1.4a it is possible to consider the meniscus as a thin round disc with a vertical lateral surface. H e r e a (20rig/pig) 1/2 is the capillary constant, which does not exceed some m m for the most materials. Considering elPe I = e17 in Eq. (3), where the order of the p a r a m e t e r y = O(1), we have:

-e2 o~~ + ' x ~

175

(6)

G]h = -- (T c - T(O))/h~.

(11)

The thermocouple measurements of the linear temperature distribution in a Si melt and the "central" part of the meniscus were done in Ref. [12].

The boundary conditions for Eq. (6) are u,]~= 1 = 1,

(7)

u,l¢=o=(To- T a ) / ( L u',~ + Biu,k_ , = 0,

Ta) =Uo,

(8) (9)

where T 0 = T(0) is the temperature at the meniscus base. Eqs. (7) and (8) determine the dimensionless t e m p e r a t u r e at the top and bottom of the meniscus. Eq. (9) determines the "intensity" of meniscus heat-exchange with the surroundings, including possible radiative heat-exchange with Biot number Bi aro/A~, were a is the effective heat-transfer coefficient of a meniscus surface with the surroundings [9]. The "intensity level" of the heat-exchange is determined by the Biot coefficient. We shall assume the worse case, when Bi = O(1) corresponds to intense heat-exchange between the melt and surroundings. Using the approach of Refs. [10,11], the solution of Eqs. (6)-(9) can be represented as w ( x , ( ) = tOo(X,() + elwl(x,~) + ... + O ( e ( ) . After substitution into Eq. (6) we find that the " m a i n " t e m p e r a t u r e distribution in the meniscus is linear" =

~0(~) = (1 - u0)~" + Uo.

(10)

The approach to the analysis of the temperature distribution near the lateral meniscus surface requires the coordinate transformation p = ( 1 x ) / a l, considering ( as unchanged. For the internal solution in this case one can write w ( p , O =

3.2. The thermal conditions in the crystal The thermal problem for the solid phase, in dimensionless variables, can be written as

a2u~ @2 •

au~ Oy

1 a (xaU~ ) = s x ax k ~ O,

es,Y...T - + e 2 - __

U y=0

s 0
Here Pe~ = 7

r

Usn + Bius]r(y ) = 0. depends, in general, on the time

t/.r, and U's, is a temperature gradient along the external normal ~ to the lateral crystal surface. The boundary conditions on the upper crystal end at y - o ~ are not important in this case, because the task consists in the determination of the temperature gradient function G~j(0), and not the temperature itself. Assumption of small macroscopic curvature of the crystallization front allows one to represent the task as one-dimensional by averaging of temperature u~(x,y) over the cross section of the crystal S(y)=~rr2(y). Presence of local curvature near the three-junction line does not influence upon integral conditions of the crystal heat transfer. After input of the new variable Ts(Y) =

2~(r(y)Us( X'y )x "10

~

dx,

integration of the heat transfer equation termwise, considering Us(X,y)[r(y ) = Ts(y), gives the

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

176 one-dimensional equation:

d (r2(y)dTs(Y))

d--y-

~

-

esYr2(Y)

eTa(y) dy

- 2BiE2r(y)T~(y) = 0.

(12)

Taking into account the slope of the lateral surface of the perturbed part of the crystal r'z = tg/3, the boundary condition can be written as U'sx = -2 t t - B i u s + e s rr(Y)Usy, where x = r ( y ) and Bi = aro/AsCOS ~. Eq. (12) with the corresponding boundary condition describes the quasi-stationary state (on the scale of time r ) of the perturbed crystal for the current crystallization time. Denoting as Ts, G~0 and ?0 the current values of averaged temperature, gradient and crystal radius, satisfying Eq. (12), we seek a variation of the temperature gradient at the crystallization front 8G s as a function of the crystal radius variation ~r. After linearization of the Eq. (12) with due account of the order Peclet number P%(t/z) = Yo + ~y(t/~') in consequence of the variable crystallization rate V~(t) = Vc0 + gVc(t), one gets the following equation: 2Bi _ ~G~, - e~y0~Gs + e2~-o ~T~ 2[

-

2Bi_

]

?o (G---s°~r)'r+Esy°G---~°~r+es~oT~r]'

(13)

value Gs(0) may be expressed from the initial gradient in the crystal Gs0 = Gs0. It is interesting that the convective component esY0~G~ may make a more substantial contribution to the heat-exchange during transfer processes in comparison with the lateral crystal surface heat-exchange (e22Bi/?0)~T~. The exact solution of Eq. (13) demands to take account of the unknown shape of the lateral perturbed surface. For purpose of analysis it is enough to use solution (14). Note here that the values br(t), ~h(t), ~Vc(t) and so on may be determined as a solution of the set of LOM equations (see Eqs. (33) and (39)). With the lower reliability these conclusions may be propagated on non-stationary growth of a crystal with axial symmetry, for example, on the growth of the initial " c o n e " with / 3 0 ) = 0. The functions G~o(t) and T~(t) are slowly changing and may be considered as any constant during the long time of the transfer process. It allows one to use the approximation of the "frozen" coefficients [13] in design of feedback control system. Finally, taking into consideration variations of the melt temperature at the meniscus base ~T(t) and the pulling rate variations ~Vo(t), after corresponding substitutions into the Eq. (5) we get one of the equations of the linear model in the common case of the non-stationary initial state as

~h(t) = A h r ( t ) ~ r ( t ) + A h h ( t ) ~ h ( t ) -- ~l£I(t) + Ahr( t )~T( t ) + ~Vo( t ),

(15)

with boundary conditions ~T~(0)= ~T~() = 0. If we assume e~ ~ 0, Eq. (13) is simplified and its solution is

where the heat coefficients of the model will be written as

Gs(0 ) = - (2G-~0/P0)3r(0) .

Ahr(t) = _ (2XsGso(t))/o~Lr(t),

(14)

With second order accuracy Eq. (13) becomes _ , + esYo~so~r ] ~G£y - esY0gG~ = - ~2 [(G~0gr)y and its solution is equal to Eq. (14) too. This means that the value of the gradient Gs(0) variation is independent of the heat-exchange conditions to second order accuracy. Physically, this is due to locally variation of the crystal geometry near the front crystallization. Taking in to account the quasi-stationary approximation, the

Ahh(t) = A l ( T c - T ( t ) ) / P s L h 2 ( t ) , Ahr ( t ) = AJpsLh( t ).

(16)

4. The capillary conditions Let us take into consideration the slope angle of the crystal surface at the three-junction line ~(t) = a(t) - e o (Fig. 1), where e0 is the constant growth angle [14-17]. The equation of the crystal

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188 d i a m e t e r change is k(t) = Vc(t)tg/30). If the initial n o n p e r t u r b e d state is " c o n e " growth with /30 = 0 ) = 0 , linearization of ~(t) to a second o r d e r accuracy gives

crystal mass is: t 2

m(t)=TrPsfor ( t l V c ( t ) d t + p s O ( r ( t ) ) .

t*(t)

+Arh(t)~h(t ) +Artt(t)~tgI(t)

=

(17)

-Y-p,g2(r).

T h e s e coefficients are:

Arr( t ) = Vc( t ) /3'r/COS2/3,

Z rh(t) : Vc(t)/3h/COS2/3, ArH(t)=--tg/3,

A~v(t ) = tg/3.

(18)

. /3r = ah' =-(h'Jh'~). For H e r e /3h' = O ~. h / h. a . and the crystal size r 0 >> a the Tsivinskii [18] approximations of h =f(r,a) are appropriate, so we have

a2sino( . . . . acoso . 4r

1

4r

PI(Vg ~ O(r))

= 7rp2r(t)h(t) + 7rp,ar2(t)cos a(t)

+ Arv~Vo( t ) .

h',-

2rcoso) B

asina

H e r e 12(r(t)) is the crystal volume b o u n d e d by the crystallization front (Fig. 1). In the simplest case the front has the shape of a sphere s e g m e n t with curvature radius R 1 >> r. T h e u p p e r signs in Eqs. (20) and (21) c o r r e s p o n d to the convex to the melt crystallization front, and the lower ones to the concave front. So in Eq. (19) we have

~M( t ) = 7rp,R2~tgf( t ), 8th(t) = zrpsr[2Vc~r(t ) + rSVc(t)] + ps.Q'r~t:

+ o(22), 8/2(t) =/x'rSt:(t) + iz'hSh(t ) T-plO'rS~ ,

t

(21)

B

and a 2 COS a [ h'~4r 2 1

(201

T h e meniscus mass is d e t e r m i n e d by the relation:

~i(t) = A ~ ( t ) S r ( t ) +Arh(t)Sh(t )

Ark(t) =--tg/3,

177

a cos a 4r B ,

)

where B = (1 - sin a + a 2 COS 2 a/16r2) -1/2. It is also possible to adopt H u r l e ' s solution [19]. T h e current values are Vc(0)= Vc0 = V0 - H 0 , /3(0)= Eo, h(t)= h o, r(O)= r o in the simplest stationary initial state when /3(0)= 0 and so ArhArtf =Arv = 0. In general /3 4:0 and V c ( t ) = V¢(/3(t)), see Ref. [20].

(22)

where /x'r = 7rpl(2rh + a 2 cos a - a2r s i n ( a ) a ' r) and /z'h = 7rpl(r 2 - aZr sin(a)a~). Substituting Eqs. (22) into Eq. (19) one gets the last linear equation of the Czochralski L O M :

815I( t) =Al_lr( t)Sr( t) + Al4~( t)8~( t) +AHh(t)ah(t ) + A u v ( t ) S V o.

(23)

H e r e the model coefficients are:

AHr(t ) = -2psVcr/A , AHh(t)=(TrPsr2--~h)/'rr A, AH~(t)

5. T h e m a s s b a l a n c e

We shall denote current values of the melt mass in a crucible, the pulling crystal and meniscus mass as M(t), m(t) and /x(t), respectively. T h e mass balance M(t) + m(t) + / x ( t ) = const, in the conservative Cz and L E C processes allows one to write 8m(t) + 8tz(t) + 8M(t) = 0 or 8rh(t) + 8~(t) = -Sh~/(t).

(19)

According to J o h a n s e n ' s analysis [21], taking into account the crystallization front curvature, the

:

[ - - ~ r ~- ( P s - P,)O'r]/erA,

At_iv(t) = -p2sr/A,

(24)

where A = p l R 2 - p s re and by use of the first approximation R 1 = const., we have .(2'r = rrr3(R 2 _ r 2) 1/2.

6. T h e o u t p u t e q u a t i o n

It is easy to d e t e r m i n e the form of the output equation for the Cz and L E C m e t h o d s with the

G.A. Satunkin / Journal of Crystal Growth 154 (1995) 172-188

178

ideal sensitivity of the weighing sensor assumption. The weighing control was discussed in more detail in Ref. [6]. It is obvious that the rate of the weight change rather than the weight

mf under encapsulant layer with density pf, can be expressed as f a ( t ) =gpf(mf(t)/ps +ix(t)/p,) = Eam(t) + F a ( t ) ,

F(t)/g = rrpsfotr2(t)Vc(t ) dt + Ix(t) + (ps - p , ) a ( r ) 2f t •

= -1rplR JoH(t) dt, should give information about the current object state. Again the plus corresponds to the convex and the minus to the concave crystallization front. The output equation Y(t) = 8F(t)/g is found by subtraction in g-m d(F0 + 8F)/dt of the program mass rate Fo/g = ~rpsrgVco = - r r p l R 2 H 0 and by substituting in the resulting equation the corresponding components from Eqs. (15), (17) and (23). Finally we get: y = -TrplR2~I2I = Cll~r + Cl2gh + d l l g T + d12~V, (25)

Y(t) = C X ( t ) + DU(t),

(26)

where C = (c H c~2) = -'n'plR2(AH~ +AHi.all + AHha21 AH~al2 +AHha22) and D = (d H d 1 2 ) = - ~ p l R Z ( A H i b H + AHhb21 AHibl2 + AHhb22 + AHV). It is principally important, that the output equation (25) directly includes the input disturbances ~T(t) and gV(t). In Refs. [1,2,4,22-24] etc. the output equation is derived from mass conservation and growth angle constancy laws only. The heat balance condition during crystallization was not taken into consideration and so the possibility of the system response determination on the heat perturbation is lost. The weighing signal in the LEC method will be determined by the equation: t 2

F(t) =Trpsgfor Vcdt +gIx + g ( p s - p , ) ~ I - Fa 2

for the initial stationary growth IX = const., r(t) = r 0 , V~=Vc0, E a = c o n s t . and F ( t ) = F 0 = 7rpsgr2Vco =-rrplgR2121o as in ordinary Cz growth. The derivative of the weight signal gives information about the state variables. Under small deviations of the growth condition we get with second order accuracy: ~k~(t) = g ( ~ r h ( t ) + ~/2(t)) +_g(ps--pl)dJ')(t) -- ~il0a( t ),

t

= -~pl R gfoHdt-Ea. The varying Archimed buoyancy force (A-force), affecting on the crystal part, including meniscus

(28)

aFa(t) = 81~.dm(t) + gila(t). Assuming Pf/Ps = pf/Pl = r/ and taking into account the relation gFa(t) = gTl(grh(t) + gt2(t)) = -g~7gA)I(t), the output equation (28) for the time interval before

were

the crystal appears from beneath encapsulant may be written as

Y(t) = - g ( 1 - r/)~3;/(t).

or in the matrix form

(27)

(29)

In general, the variation of gFam(t) in current time t, when the crystal appears from beneath the encapsulant may be written as:

Fam( to) = g~7~mf( t )

(30) In Eq. (30) the first integral describes a full increment of crystal mass variation during the crystallization time, the second one a mass variation of a crystal part, appearing from the encapsulant in current time. Here we suppose as a first approach the constancy of the time delay r d when the crystal, appearing from the encapsulant, has a simple cylindrical form. It is important that the meniscus has not a "memory" about the preliminary states: the meniscus geometry is determined by the current crystal diameter at the crystallization front. So the second component gFa, of full variation of the A-force is determined by variation of the meniscus volume in current time: gFa(t) =

179

G.A. Satunkin/Journalof CrystalGrowth154 (1995)172-188 gr/8#(t)-T-pfS~Q(t). After differentiating 8Fam(t) and 8F~(t) and substituting the results in Eq. (28) we get a final form of the output equation:

all=Arr+ArhAhr, a21 =

Ahr -- AHr -ArrAHi--AhrArHAHt 1 +AHh

8/~(t) = - g ( 1 - r t ) S M ( t ) + grtSrh(t - rd). (31) Taking into account the balance mass condition - 8 M ( t - %) = 8rh(t - r d) + 8/2(t - r d) _ g.(Ps pl)5~J(t - rd).one can rewrite Eq. (31) as 5 F ( t ) = - g ( 1 - r l ) s m ( t ) - g~78aTl(t - %) - grlSfz(t - r d) 7+g(p s - p l ) 8 ~ ( t - %). The output equation in matrix form for the LEC method is:

al2=Arh +ArhAhh, Ahh --ArhAHi --AhhArnAtti a22 =

[B(t) = ( bll b12 / b21 b22 ] ArhAhT AhT(1 - A r f I A H ~ ) 1 +AHh

=C1X(t ) +nlU(t )

(32)

The matrices C(t) and D ( t ) for LEC are equal to the matrix coefficients of conventional Cz providing for -~-pRl(1 - rt) 2. The matrices C 2 and D 2 have the components c2 H =~g~T&ro(2Vco-

roArh),

C212 --

~grlpsr2Ahh,

d211 =

--Trg~Tpsr2AhT and d212 = 0. In contrast to common Cz, the output equation for LEC includes the delaying components Y 2 ( t - rd), which "cipher" information about the current deviation on the crystallization front.

The set of LOM equations (15), (17) and (23) are equivalent for Cz and LEC. In a standard matrix form we have:

X ( t) = &( t ) X ( t) + • ( t ) U ( t ) .

(33)

Here X ( t ) = ( x l ( t ) , x 2 ( t ) ) T = ( S r , 8h) T is the vector of independent state variables and U ( t ) = (ul(t), Ue(t)) T = (ST, 8V0) T is the input vector. The components a~j and b~j of the matrices A(t) and B(t) are: f~(t) = ( a21 all with

a22 a12)

(34)

0 1 -AHv

(35)

1 +AHh

The values of parameters A i j ( t ) are clear from Eqs. (16), (18) and (24). The LOM is simplified in the autonomous case (aij, b,j are time independent)

&(t)

[a21 a= lx2(t)

[0 01(, t, )

+ b21 7. The linear model in continuous and in discrete time

'

and

Y ( t ) = Y , ( t ) + Y 2 ( t - %)

+ C 2 X ( t - T d ) +D2V(t-~'d).

1 +AHh

b=

lu2(t)

"

The LOM described by Eqs. (33) may be called a "non-inertia" model. The input components can be described by arbitrary functions, but in accordance with its physical nature. The temperature perturbation ut(t) may be caused, for example, by the heat flux changes in the melt or by heating power change. In the latter case it is necessary to determine the inertia of the melt or heating power change. In the latter case it is necessary to determine the inertia of the melt temperature response to instantaneous heat power changes d N ( t ) (T-control channel). As a rule, the real objects have appreciable inertia of temperature response which can be approximated mathematically by experimental results [25,26]. This approach complements the LOM and will be discussed later. As a rule, it is possible to ignore the input inertia of the second component u2(t),

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

180

pulling rate disturbances (in of the V-control channel (in The Laplace transform functions of the controlled

open state) or input closed state). yields the transfer object in the form X(s) = (sI - A)'IBU(s) + (sI - A)-IX(0), the output signal for weighing control Y(s)= [ C ( I s A)-IB + ff)]U(s) + C(Is - A)-IX(0) or in expanded form

~r( s) = WrT( S)~T( s ) -b Wrv( s)~Vo( s), h(s) = Wh (s) T(s ) + Whv(S) Vo(S ), s 14( s) = s) + w . v ( s) Vo( s).

are new matrices of the discrete LOM, Td is the digitization time. For the real roots of the equation D(s) = 0 we have: h 2-aii

AI - - a i i

fii = exp(h iTd) - -

A2-A

exp(AeTa)--, 1

A2 -A 1

air [exp(h2Ta) - exp(hlTd)], fir - h2 _ hi i,j = 1,2;i = j .

(40)

In the case of the complex roots AL2 = a _+ ifl we get (37)

In the autonomous case (and if X ( 0 ) = 0) one gets:

f~i

= exp(aTd) [sin ( fiTj)( aii- a)/[3 + cos(/3Ta)],

fir = exp( aTa)air sin(flTd)//3"

(41)

WrT( S) =ArhAhT/D( s), Wrv ( S) = Arh (1 -- A n v ) / D

8. Discussion

( s),

Why(s ) = (SAhT - ArrAhT ) /D( s ), Whv( S ) = Is(1 - A n y ) - A r r ( a - A n v ) ] / D ( s), W/4T(S) = [s2AHhAhT + S(ArhAhTAH~--ArrAhTAHh)

+ ArhAhTAHr]/SD(S), Wttv( S ) = [s2(AHh + AI-IV ) + s( ArhAl4 e - ArrAHh - A r r A I 4 v - - A h h A H v ) + A r r A hh A Hv -- Arh A hrA HV

+ ArhAHr]/sD(s),

(38)

where D(s) = s 2 - (all + a22)s + (alia22 a12a21 ) = S 2 -- (spur A)s + det/% is the characteristic equation of the controlled object. The procedure for transformation of the linear model in discrete time is standard [13]. By rearranging, we get the original equation in the canonical form

X ( k + 1) =A-X(k) + B U ( k ) , Y(k) = CX(k) + DU(k). Here

[ii1

./T= exp(ATa) = [f2x

B= (.h-I)A-1B,

f22 '

C = C, D = D

(39)

Let us pay attention to some peculiarities of the proposed model. (1) In spite of the fact that three conservation laws form the basis of the linear model, but only two state variables x~ = 5r(t) and x 2 = 5h(t) are independent. The output equation (26) for weighing control practically describes a dynamics of the third variable ~ / 4 ( t ) = - Y / ~ p l R2 and obviously includes the system response to the input perturbations. The LOM of the common Cz is presented in Fig. 2. It is followed from the Fig. 2, that both the subsystem of the meniscus x2(t) and the subsystem of the crystal x~(t) are fully observable and the disturbance of the meniscus subsystem triggers the change of crystal diameter. The completely observability for the LEC method required the additional condition ~'d = const. [6], but, in general, Za(t) depends on the time. (2) It should be noted that the stability of the weighing control in approach of the ideal weighing sensor is fully determined by the position of the roots of the general characteristic equation D(s) = 0. The weighing sensors of the secondorder (for example, self compensated balance [6]) have their own dynamics and in this case autooscillations and sensor instability are possible. The analysis of the closed and open states stability of the controlled object will be carried out in detail in the next paper.

G.A. Satunkin / Journal of Crystal Growth 154 (1995) 172-188

181

Fig. 2. Principal decomposition scheme of the controlled object for the convenient Cz method, defined by Eqs. (26) and (33). The scheme serves to extract the meniscus, crystal and weighing control subsystems. The state variables x = [~r(t),~h(t)] are completely observable under weighing control. The scalar output signal Y(t) is the sum of the response of the meniscus and crystal subsystems and includes the direct connection with the input disturbances u = [~T(t),~Vo(t)].

(3) When the heat conditions are considered in this LOM, the local variation ~r(t) is taken into account. Physically this is reduced to the heat flux constancy through the crystallization front J = -Trr2,~sGso = const, due to a negligible change in

the heat surrounding of the growing crystal [26]. The variations ~J(r) lead directly to relationship (14). The melt temperature changes are considered as the principle ones, as in Refs. [4,8]. However, Hurle et al. [4] determined G~r- Gso/2ro,

Table 1 Thermo-physical properties of g e r m a n i u m and operating conditions of experiments described in Refs. [4,5] utilised in calculations of transfer processes Property Thermal conductivity (solid) Thermal conductivity (liquid) Latent heat (per unit volume) Surface tension Laplace constant Growth angle Melting point (MP) Density at MP (solid) Density at MP (liquid) Specific heat Pulling speed Crystal radius Crucible radius T e m p e r a t u r e gradient in crystal at interface Meniscus base temperature Mean mass of the melt during growth period analyzed in Ref. [4] Interface displacement (convex) Time period of truncated square wave variation 5T = 5:2 ° C

Symbol

Value

As

0.06 cal cm - 1 K - 1 0.12calcm 1K-1 601 cal cm -3 616 erg c m - 1 4.77 m m 15° 1209 K 5.26 g cm i 5.51 g c m -1 0.74 cal g - 1 K - 1 4cm h -l

A1

ps L OrLG a E0

L Ps Pl Cp

vo F0

1 cm

R

2 cm 100 K cm 1 1224 K 150 g 2 mm 200 s

6so 7"0 M hfr

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

182

which is a factor of four smaller than in this work. That is an important difference because the stability and transfer processes are highly sensitive to the gradient value Gs0 and to the melt "overheating" value ( T 0 ( 0 , t ) - T c ) . Note here that practically all LOM parameters, which are necessary for the quantitative calculations, are available and well known in advance. Only three parameters gives rise to some uncertainty: the gradient at the crystal front Gs0, the temperature of the meniscus base TO and the flexure hfr (or the curvature) of the front (see also Refs. [25,26]). (4) The scope of the applicability of the small curvature approximation can be evaluated for the simple shape of the crystal front, for a convex front hfr = r0(1 - sin 0 ) / c o s 0, where 0 is the slope of the front to the vertical. The maximal flexure hfr under stationary growth when /3 = 0 ° is determined by the angle 0 = y (Fig. 1) This angle y like e 0 is the physical constant of the material and may be calculated from the known growth and welting angles for materials with known surface tension [17,25]. It is supposed in the LOM, that the maximal front flexure under

2.

0

stationary growth does not exceed the meniscus height hfr = r0(1 - sin y ) / c o s Y < h0 < 1.4a. Otherwise the flexure is too large and the meniscus base temperature loses physical sense. For instance for Ge y = 50 ° [17], the maximal flexibility (hfr/r o) = 0.364 during stationary growth. Physically large front flexibility is possible, but only when the crystal expands (/3 > 0 °). These simple estimations cannot be done for a concave front. (5) Let us use the LOM for some calculations of the transfer processes in Cz growth of the Ge crystal. As initial date it is convenient to use the experimental data and the parameters from Hurl and coworkers [4,5]. These data are presented in Table 1. Additionally, the crystallization front flexure was taken to be hfr -- 2 mm. Fig. 3 shows the results of the calculations by Eqs. (39) of the transfer processes caused by the non-inertia step wise increase (curves 1-5) and decrease (6-8) of the melt temperature. Stationary growth, corresponding to the point A(ro,h o) on the curve of the stationary meniscus height h0(r0,e0), occurs for Ge with e 0 = 15°. The general stability of the system (in open state) allows it to transfer to new

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A(ro.ho)

h(ro.~o)/a

O. 5-

o.o

. . . . . . . . . O.

I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.5

1.0

!

5

2

0

i ......... 2.5

i 3.0

r/a Fig. 3. Transfer processes during step-wise temperature increase and decrease of the meniscus base. The parameters for Ge are presented in Table 1. The initial stationary crystal growth with r 0 = 10 mm and flexure of the convex front hfr = 2 mm corresponds to point A(r,h). The stationary meniscus heights are calculated by use of the Tsivinskii approximation h(r,a) = [a2(1 - sin a) + a4cos 2 o~/16r2] 1/2 - a2cos a / 4 r [18]. The maximal heights are taken from Mika and Uelhoff [28]. The values 8T = 30, 24, 16, 8, 2°C corresponds to curves (1-5); gT = - 2 , - 1 8 , - 1 6 ° C corresponds to curves (6-8). The meniscus break should be waited at gT = + 24° C after 56-60 s (curve 2)(see also Fig. 5).

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

for the crystal to break away from the melt as a result of the attainment of the maximal meniscus height Hmax. This height Hma~ is the numerical solution of a Jacobi equation and is presented in

stationary states and in the Cz-method in the first instance it is determined mainly by "overheating" of the melt and the gradient along the crystal. There is the critical gT = 24 ° C when it is possible

"~,

- . ~ - " :zJ--

o.

~

-/-

"~

-

~•-

E

.~

~

-8

. . . . . ~. . . . . . Me(t)

:

/

--;

i . . . . . . . . . . . . . . . . . . . . . .

'

!

/D \ ~ 5 : ~",,t, '-

~:' ......

' d. o ~-~l'~t)

L

~

-~ mm~

_,.

~ .............

~ .............

:

i

1 0',0_ 0

!, _

_

1 5:0.

Y(t)

!

~

::

183

2_ :

-

~..~s

(a) !

~. . . . . . . . . .

Timei(x2sec) 0

2 0;0.

_~ ...........

0

! 2 5:0-

!,

,!

L.............

i

0

8T=+2 C

-1.5-

.r,4

8.0 ................

; ...............

: .............

; ............

; ........

-,b(i

~

::

:: 5 T = - 1 6 C

::

::

' '!

~I\\ w(')

:

::

, ::

::

Lo

q/k

:,

~,-

q]N.

:i

b

~/ -,

o3!

,_ ---. ~ - - ~: ~

~...-f

.........

i ~ ( x z s e o )

:

',

,,

,,

L ............

~

,,

8h(') I

. . . . . . . . . . . .

" ............

:: . . . . . . . .

".

. . . . . .

:

Fig. 4. Transfer processes are calculated with Eqs. (33) and (39) for the positive (a) and negative (b,c) changes of the meniscus base temperature. The characteristic roots of equation D ( s ) = 0 (see Eq. (38)) are A1 = -0.00644, A2 = -0.2176 and all = 1.0, a l a = 1 . 1 9 8 × 1 0 - e , a 2 1 = 4 . 6 9 4 X 1 0 -2 , a 2 2 = 9 . 4 4 3 × 1 0 1, bl I = _ 7 . 7 0 7 × 1 0 S, bl 2=_1.749×10 -2 , b 2 1 = 1 . 2 3 9 × 1 0 5, b2 a = 2.812, c n = 2.232 × 10 -2, Cl: = -3.096 × 10 -2, d l l = 4.060 × 10 6, dl 2 = 3.091. When the disturbance increases to ~ T = + 4 ° C , the calculated crystallization rate becomes negative ~Vcr < 0 ((a), dashed line). It is considered that inside segment ED the crystallization rate Vet(t)= 0 (i.e. remelting is absent). This means the crystal does not change its radius t:(t)= 0 in the time interval t E - t D and the meniscus height increase at a rate V0 - / t ( t E) to the moment when Vcr(t) > 0 again ((b), dashed line). Taking into account the assumption, mentioned above, that the calculated profile curve ~r(t)ldr--16oc coincides with the curve A E G (c). The output signal Y(t) is initially opposite in sign to the crystal radius change; this is possible at the scalar input for materials with density ratio p~/p~ < 1.

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

184

O. 0 .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

(c) :dT=-24C: 6.

0 ¸

:dT=-16C! ,~4.

O.

f D !i

A(ro,ho);; O. 0

G

i'

2. 0 ¸

O. 0

1.0

'

i

i

i

2.0

3.0

4.0

5.0

7

6.0

Lcr(mm) Fig. 4 (continued).

the tables in Refs. [27,28] by Mika and Uelhoff. A more general arrangement and the solution of the meniscus stability task is given in Ref. [29]. (6) It is clear that the input perturbation u ( t ) = [SVo(t),gT(t)] may be arbitrary. The responses of

the growth system, that ~Vo(t) - ~[t(t) - ~I:I(t) may be easy calculated, For a specified input infinite large variations

are ~r(t), 5h(t), 8Vc(t) = are predetermined and using Eqs. (33) and (39). u(t) the LOM excludes 8Vc(k) = 8Vo(k) (a21 + -

1.5

..................

~_-<

. . . . . .

!

,,

i ..........

"4--

D

-

. . . . . . . . .

:. . . . . . . . . . . . . . . . . . . . . . . .

!

]

i

!

i .........

~____- _" _--_,~ . . . . .

!

v

..........

0.5

~1

0.o(

5

l; 0

i --0.

1 5

" ........

roho) i

2~ 0 ~

!

~

i

0

a; 5

--i

!

4; 0

E

fi

Rmax/a;

r(t)/a

Fig. 5. Trajectories o f transfer processes for temperature disturbances with allowance for the calculated segment ED, corresponding to ~Vcr < 0 (see also Fig. 3). The line BC is determined by a chosen flexure o f the front. The values ~T = 30, 24, 16, 8, 2 ° C corresponds to curves (1-5); 5T = - 2, - 18, - 16, - 24, - 30 ° C corresponds to curves (6-10). The curve o f the limit break height Hmax/Rmax is taken from Ref. [28]. The negative meniscus heights (curves 9 and 10) cannot formally be considered in the scope of the proposed model.

G.A. Satunkin / Journal of Crystal Growth 154 (1995) 172-188

c11)~r(k) +(b22 +

185

under positive S T > 0. These calculations show the boundaries of the proposed LOM reliability. Remelting of the crystal is not described by this model when it actually takes place in the transfer processes. However, if one takes into account that inside the time interval ED the crystallization rate is equal to zero, one may consider the transfer processes as continuous in the scope of the LOM (dashed line in Fig. 4b), but in this case a bend in point E on the crystal surface r(l) appears (Fig. 4c). Using this approximation, one can build the allowable trajectories of the transfer processes during the step wise input temperature perturbations, which are presented in Fig. 5 (and added to Fig. 3). Here ED are the intervals, on which Vet(tED) = 0 is taken and in addition the responses on the large negative temperature jumps are shown (curves 9 and 10), when the output from the limits of allowable application of the LOM is possible (h < 0 and h < hfr). (8) For the chosen growth parameters, breakage of the meniscus is expected at St = 24°C after 60 s, when the maximal meniscus height has been achieved (curve 2 in Fig. 5). Note here, that

- (a22 + C l 2 ) ~ h ( k ) + (b21 + dlt)~T(k) d12)6Vo(k) for any time t = kTd. In this

regard the results of Wilde et al. [24] and Johansen [30] are interesting. If the crystal geometry (or the disturbances 8r(t)) is given by some arbitrary function, then the crystallization rate determination according to two conservation laws results in divergence! (7) The time-dependent transfer processes for the scalar input u = [ul,0] are given in Fig. 4a and Fig. 4b. The temperature rise decreases the crystal diameter and increases the meniscus height. Note here that the output signal is initially positive Y(t) > 0; that is natural for materials such as Ge with pl > Ps- When the step in ~T is + 2 ° C, the crystallization r a t e Vcr > 0, but as soon as 6T = + 4° C "running" of he crystallization front to the cold zone produces a negative crystallization rate (region ED in Fig. 4a). For a step wise decrease of the melt temperature ~T = 16° C, there is an interval of time ED (Fig. 4b) for which Vcr > 0. In this case it is difficult to speak about possible remelting of the crystal during cooling of the melt, although this is the case -

1.,5

1

,,~ 1. o

"~

1

:

1

:

'"

!o

~"

H.~Rr~

L

~o.5

...............

i

3

'

.

~

!

!

:

,

-e .......................

:

1 B

O.

o

,

0.5

,

,

i

,

,

,

,

,

I 1

0

r

,

,

,

,

,

f

i

, 1

5

,

,

,

.

.

,

,

.

,

p

i 2 . 0

.

r

.

i

.

.

.

,

,

i

C

,

,

,

i 2 . 6

,

,

,

,

,

,

,

,

, 0

0

Rm~x/~: r(t)/~ Fig. 6. Trajectories of transfer processes for pulling rate disturbances with allowance for the segment ED, corresponding to 6Vcr < 0 (see also Figs. 3 and 5). This means in this case, that the crystal does not change its radius t:(t)= 0 and the meniscus height increases at a rate Vo + 6V0 - / t ( t E) to the time tD, when Vcr(t) > 0 again. The values ~Vo = ( - 0 . 5 , - 0.9)V0 corresponds to curves (1) and (2); 8Vo = (1, 3, 7, 9, 11, 15)Vo corresponds to the curves (3)-(8). The break of the crystal from the melt at the step 8V0 = 7V0 occurs after 90-100 s.

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

186

the results of similar calculations can be used for the experimental identification of the process parameters. (9) The calculation results in the case of rate perturbation u(t)= [0,~V 0] are quoted in Fig. 6, 7a and 7b. Naturally there is a limit of pulling rate increase, above which m e n i s c u s breakage

v

1

takes place. In this case the breakage happens at ~ V 0 = 7 V 0 after 90 s. In Fig. 7a t i m e - d e p e n d e n t transfer p r o c e s s e s are p r e s e n t e d at the positive step ~ V 0 = 2 V 0. During the transfer process the crystallization rate remains positive. At ~V 0 = 3V 0 on curve V~(t) (shown by a dashed line) s e g m e n t E D , w h e r e V~(t)< 0, appears. W h e n the steps of

.. :

:

,

,

(t) / / z

i

/

....

i

~. _

v o(t)

. . . . .

,~ ...............

. . . . .

:: .-

o

,

:

.'- . . . . . . . . . . . .

(a)',

: .......

:

i

i

',__

:

:

,,- . . . . . . . . . . . .

,~ . . . . . . . . . .

1ram I ~Vc=2V

o

I --3

0

O[

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

\

v o~

,,,.,

v

,t=

o~ f-~

ar(O

: O.

.

0

.

.

.

.

.

.

1510. 0

~:

.

.

.

.

20~.

_

av~(o

.

.

.

.

0 .

2510. .

.

.

~ - 1 . O


5Ve="O.SVo

;>

--2.

O

. . . . . . . . . . . . . .

L . . . . . . . . . . .

• . . . . . . . . . . . .

Time

i

L .

.

.

.

.

.

.

.

.

(×2see) L . . . . . . . . . . .

,

Fig. 7. T r a n s f e r p r o c e s s e s for positive (a) and negative (b) pulling rate changes. W h e n the d i s t u r b a n c e i n c r e a s e s to ~Vo = 3Vo, the c a l c u l a t e d c r y s t a l l i z a t i o n rate b e c o m e s negative 3Vcr < 0 (dashed line). Inside s e g m e n t E D it is assumed V~r(t) = 0 (the remelting is absents) and J:(t) = 0. The meniscus height increases at a rate Vo + 8V0 - / ~ ( t E) to the time tD, w h e n Vcr(t) > 0. The o u t p u t signal Y(t) is initially opposite in sign to 3r(t). The crystallization rate variation 3V~r(t) is shown at ~Vo = - 0 . 5 V o (b). For the crystal growth conditions and physical parameters, the crystallization rate Vor stays positive w h e n - Vo < ~V0 < 0.

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188 :3.

0 .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

187 .

.

.

.

.

.

.

~T=+-2C ~

2,

O

4

1 . O-

..........

60.0

1. 0 .

.

.

.

-\-~

: . . . \. . .\ 1

":-- \

/ I f * ( \. . . . . . . .

-'< \ 1 , ' ~ - ~

. . . . .1 7 - " . . .\.-.'.i .i

(~rramp -2.

0

.

.

.

.

.

.

"

~. . . . . . . . . . . . . . . . . . . . . . . . . . .

L. . . . . . . . . .

Time

(xlOsee)

Fig. 8. The calculation of transfer processes for truncated square wave variation 5T0 = _+2°C with time period 200 s. The growth parameters are given in Table 1 and agree with data of Hurle and coworkers [4,5]. The change 5r(t) at ~T0 = _+4° C is shown by a small dashed line. In both cases, at the initial heating of the melt, an apparent expansion of the crystal appears. Additional introduction of the constant ramp 3.4 × 10-5 cm s 1 leads to stabilization of crystal diameter (curve ~ r r a m p ) a f t e r 700 s (see Ref. [41). pulling rate changes are negative, the condition Vc(t) > 0 is always fulfilled (Fig. 7b) (for the chosen growth parameters). Note that the weighing signal Y(t) behaves " a n o m a l o u s " (Figs. 7a and 7b) and this is d e t e r m i n e d first of all by the fact that pl > Ps, but not by e 0 > 0. T h e possibilities of vector input u = [ul,u 2] allows one to make an " a n o m a l o u s " situation for materials with the " n o r m a l " densities relation p¿ Ps. T h e o u t p u t consist in a sum of response of meniscus and crystal subsystems that will be considered in following publications. (10) Using data of Hurle and coworkers [4,5], we consider the response to the input ~T = _+2 ° C of a truncated square wave with a period of 200 s. The calculation results are given in Fig. 8. T h e crystal radius variation u n d e r 3 T = _+4°C is shown by the dashed line. If in the first half period of the oscillation the t e m p e r a t u r e increases (corresponding to the firs decrease on the curve ~r(t), see Ref. [4], Figs. 4 and 8a), then the long-time crystal expansion is marked. However, at the end of the calculated time the transfer process finishes and a p p a r e n t expansion of the crystal disappears. This effect, to which attention

was paid by Hurle et al. [4], does not d e m a n d nonlinear effects in the scope of the L O M . Characteristically, if at the first oscillation the crystal is first cooled, narrowing is observed instead of expansion. Probably, the small load mass (150 g) did not permit observation of the end of the transfer process in the experiments. However, this fact allows one to neglect the thermal inertia of the melt t e m p e r a t u r e response to the change of heater power in the calculations. (11) For completeness of the comparisons, in Fig. 8 the calculated crystal radius response during modulation and at a constant increase of the pulling rate (with ramp rate 3.4 × 10 -5 cm s 1, see Fig. 5 of Ref. [4]) is presented. Really, this is a vector input that allows one to overcome the crystal "expansion" after 700 s from the start of the transfer processes.

References [1] M.A. Gevelber and G. Stephanopoulos, J. Crystal Growth 84 (1987) 647. [2] M.A. Gevelber, G. Stephanopoulos and M.J. Wargo, J. Crystal Growth 91 (1988) 199.

188

G.A. Satunkin /Journal of Crystal Growth 154 (1995) 172-188

[3] M.J. Wargo and A.F. Witt, J. Crystal Growth 116 (1992) 213. [4] D.T.J. Hurle, G.C. Joyce, M. Ghassempoory, A.B. Crowley and E.J. Stern, J. Crystal Growth 100 (1990) 11. [5] D.T.J. Hurle, G.C. Joyce, G.C. Wilson, M. Ghassempoory and C. Morgan, J. Crystal Growth 74 (1986) 480. [6] G.A. Satunkin and A.G. Leonov, J. Crystal Growth 102 (1990) 592. [7] V.V. Voronkov, Kristallographia 23 (1978) 249. [8] A.B. Crowley, IMA. J, Appl. Math. 30 (1983) 173. [9] H.S. Carslaw and J.C. Jager, Conduction of Heat in Solids (Clarendon, Oxford, 1959) ch 8. [10] A.H. Nayfen, Introduction to perturbation techniques (Wiley, New York, 1981) p. 535. [11] I.E. Zino and E.A. Tropp, Assymptotic methods in the task of heat transform and termoelasicity (LGU, St. Petersburg, 1978) p. 210. [12] E. Kuroda and H. Kozuke, J. Crystal Growth 63 (1983) 276. [13] B.C. Kuo, Digital control systems (Holt, Rinehart and Winston, New York, 1980) p. 445. [14] V.V. Voronkov, Phys. Status Solidi (USSR) 5 (1963) 571. [15] W. Bardsley, F.C. Frank, G.W. Green and D.T.J. Hurle, J. Crystal Growth 23 (1974) 341. [16] T. Surek, Scripta Met. 10 (1976) 425. [17] G.A. Satunkin, V.A. Tatarchenko and V.I. Shaitanov, J. Crystal Growth 50 (1980) 133.

[18] S.V. Tsivinskii, Inzh.-Fiz. Zh. Akad. Nauk Belorussk. SSR 5 (1962) 59. [19] D.T.J. Hurle, J. Crystal Growth 63 (1983) 13. [20] G.A. Satunkin, S.N. Rossolenko, V.A. Kurlov, B.S. Red'kin, V.A. Tatarchenko and A.M. Avrutik, Cryst. Res. Technol. 21 (1986) 1257. [21] T.H. Johansen, J. Crystal Growth 118 (1992) 353. [22] W. Bardsley, D.T.J. Hurle and G.C. Joyce, J. Crystal Growth 40 (1977) 13. [23] V. Nikolov, K. lliev and P. Peshev, J. Crystal Growth 112 (1991) 451. [24] J.P. Wilde, L. Hessenlink and R.S. Feigelson, J. Crystal Growth 113 (1990) 337. [25] G.A. Satunkin, B.S. Red'kin, V.A. Kurlov, S.N. Rossolenko, V~A. Tatartarchenko and Yu.A. Tuflin, Cryst. Res. Technol. 21 (1986) 995. [26] G.A. Satunkin and S.N. Rossolenko, Cryst. Res. Technol. 21 (1986) 1125. [27] K. Mika and W. Uelhoff. J. Crystal Growth 30 (1975) 9. [28] W. Uelhoff and K. Mika, Berichte Institute fur Festkorperforschung, Kernforschungsanlage, Julich, May 1975, pp. 1-135. [29] A.V. Zhdanov, G.A. Satunkin, And R.P. Ponomareva, J. Colloid Interface Sci. 104 (1985) 334. [30] T.H. Johansen, J. Crystal Growth 126 (1993) 347.