Survey of quantitative analyses of the effects of capillary shaping on crystal growth

74

Journal of Crystal Growth 82 (1987) 74—80 North-Holland, Amsterdam

SURVEY OF QUANTITATIVE ANALYSES OF THE EFFECTS OF CAPILLARY SHAPING ON CRYSTAL GROWTH V.A. TATARCHENKO Institute of Solid State Physics. Academy of Sciences of the USSR, Chernogoloi’ka, Moscow District 142432, USSR Received 22 July 1986

Factors influencing the dynamical stability of the meniscus for capillary shaping techniques of crystal growth are reviewed using the method of Lyapunov, and the application of the analysis is demonstrated for special cases involving the Czochralski and Stepanov processes.

1. Introduction The methods of crystallization from the melt without any contact of the side faces of a growing crystal with the container walls are now becoming increasingly widespread. The shape of the growing crystal is mainly determined by a meniscus crosssection formed in the vicinity of the solid—liquid interface. Since the meniscus is formed by capillary forces, these methods are termed the capillary shaping techniques. It is interesting to carry out a comparative estimation of the applicability of the different capillary shaping techniques for growing shaped crystals. In this work, we propose to analyze the dynamical stability of crystallization processes according to Lyapunov. This analysis allows us, for the conventional crystallization techniques involving capillary shaping (the Czochralski, Stepanov, Verneuil and floating zone techniques), to find the steady-state modes for growing crystals with a controlled shape of the cross-section. We shall state a general approach and demonstrate some of its applications. The general scheme involves coupled solutions of the equation of melt flow (the Navier—Stokes equation) with the boundary condition at the free meniscus surface (the Laplace equation), the continuity equation (the mass conservation law), the heat conduction equation (the energy conservation

0022-0248/87/$03.50

law), the diffusion equation (the impurity mass conservation), and the use of the condition of the growth angle constancy, specific for crystallization from the melt. This system of equations is general for all the crystallization techniques involved, the singularities of each technique being characterized by the system of boundary conditions and real values of the parameters included in the equations. Each of the crystallization techniques considered here may be characterized by a finite number of the main variables. In the Czochralski, Stepanov and Kyropoulos techniques, these are the crystal size and the position of the crystallization front with respect to the free surface of the melt or the shaper edges; in the floating zone technique, they are supplemented by the melt zone volume and the position of the melting front. These variables (x,), in the general case, are functions of the control parameters of the crystallization process, and they fit the following system of differential equations derived on the basis of the above-mentioned conservation laws: =

i~ (t,

x

1~...x~,c),

where I is time, n is which depends on the and the cross-sectional tamed, and c represents

© Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

i = 1, 2

n~

(1)

the number of variables, method of crystallization shape of the crystals oball of the control parame-

V.A. Tatarchenko

/

Effects ofcapillary shaping on crystal growth

ters of crystallization, and thermophysical and other constants of the crystallized material, We are interested in the steady-state values of the variables x,0, i.e., the growing crystals with a constant cross-section at a steady-state position of the crystallization front and other parameters of the system. These values of x,0 satisfy the system of equations (1) with a vanishing left side:

75

complete solution of the problem, as stated, is mathematically a very difficult proposition. It is therefore of interest to estimate the effect of some factors on memscus formation. In crystal growth, three types of forces participate in this formation: (a) the forces of inertia related to the melt flow; (b) capillarity forces; and (c) gravity forces. An estimate of the effects of these forces is given by the 2L/y characteristic numbers: Weber number We of characterizing the comparative effect pV inertia forces and capillarity forces; Froude number Fr V/ ~/~Efor inertia and gravity forces; and Bond number Bo pgL2/-y for gravity and capillarity forces. Here p is the material density, L is a characteristic linear dimension, V is a characteristic velocity, -y is the surface tension coefficient, g is the gravity acceleration. It appears that, for crystallization rates up to 10_102 cm/s, the effect of inertia forces as com=

x~0, c) 0. However, only steady-state solutions of eq. (2) physically realizable. According to Lyapunov, solutions of equation (1) are stable if they stable for the corresponding linearized system =

=

dt

~

~

—

X~).

(2) are the are

(3)

k=1 aXk

All partial derivatives in (3) are evaluated at x, x, 0. Eq. (3), in turn, is stable when all the roots ~ of the characteristic equation =

det( 8Xk \

—

~

k)

(4)

have negative real parts (~jk is Kronecker’s symbol). The determination of the non-stationary functions j, which depend explicitly on time, requires the solution of the equation and, usually, entails considerable difficulties. The difficulties can be overcome by using the quasi-stationary approach. Mullins and Sekerka [1], for instance, used this approach to analyse the morphological stability of crystallization front perturbations and obtained results in agreement with experimental data. However, in each specific case, the quasi-stationary approach calls for substantiation.

2. The conditions of formation of the molten column (capillarity problem)

=

=

pared with capillarity and gravity forces can be ignored (fig. 1), i.e. instead of the solution of the Navier—Stokes equation, it is sufficient to solve the equation. In addition to the flowLaplace caused capillarity by crystallization, convective flows with velocities substantially higher than the crystallization rates may appear in the liquid column. The effect of these forces on the meniscus formation and on the processes of heat transfer in the liquid phase will be discussed later. Elementary geometric manipulations demonstrate (fig. 2) that, for constant growth angle, the

We 100

1

Ba 0.01

0.1

$1

10 1 ~10~100

~Bo~

0.1 2

In capillary shaping, the crystal cross-section coincides with that of the molten column at the crystallization surface. The molten colunm shape can be found by solving the Navier—Stokes equation with the Laplace capillarity equation serving as a boundary condition at the free surface. A

0.01 We

______________________

Fig. 1. Schematic representation of the effects of inertia (1), capillary (2) and gravity (3) forces on the molten column formation [2}.We, Fr and Bo are the Weber, Froude and Bond numbers, respectively, and Fr2 =WeBo.

76

V.A. Tatarchenko

/

Effects of capillary shaping on crystal growth

in the liquid and solid phases, ~ is the latent heat of melting per unit volume, J~ V— dh/dt is the crystallization rate, h is the crystallization front

Z

=

1

position, and V is the pulling rate. From eq. (6), 1

V

dh/dt

\

~

R~ ~

2 --

—

-

-

--

-

—

- --

--

r --

Fig. 2. Change of the growing crystal diameter at a0 ~ a1: (1) crystal, (2) melt.

crystal radius R depends on time according to =

—

l’~tan(a0

—

cxi),

(5)

where V~is the crystallization rate, ~ is the angle between the meniscus and the horizontalthelinetangent at theto three-phase intersection point, and a 1 ~ In the case of steady-state growth, a0 ae. The angle a0 as a function of the crystal size, crystallization position and other factors is found, together with the meniscus shape, by solving the Navier—Stokes equation with the corresponding system of boundary conditions, depending on the specific crystallization technique. =

v5G5).

—

(7)

The quantities GL and G5, being the functions of the crystal size, position of the crystallization front, and other parameters of the process, have to be found by solving the heat conductivity equation for the melt and crystal, respectively. It was mentioned above that in the general case this requires that non-stationary equations be solved. The quasi-stationary approach used here implies that the temperature distribution at any moment satisfies the stationary heat conductivity equation for the instantaneous position of the crystallization front. This approach is valid if the characteristic time of the temperature relaxation is much smaller than that of front relaxation, after a perturbation, back to the steady-state position. Both GL and G5 are then found by solving a one-dimensional con.

.

.

vective equation taking into account heat exchange at the lateral surface [4]: 2T, V dT, ~ d I( T i~) 0, (8) ~

—

—

—

=

—

=

3. Thermal conditions of crystallization The condition of the heat balance at the crystallization front enables one to find the interface velocity 2~fr’.. ~ PLGLPSGS—

V— ~(vLGL

~

‘~

dR/dt

=

Here VL, v~ the are heat conductivities of the liquid and the solid phases, GL and G 5 are the temperature gradients at the crystallization front

where i L,S (L for the melt and S for the crystal), T, is the crystal temperature, ,c~is the thermal diffusivity coefficient, z is the coordinate, ~u, is the coefficient of heat exchange with the ambient medium, I is the ratio of the perimeter of the crystal cross-section to its area, Tm is the ambient temperature. The second term in (8) takes into account the melt flow and the crystal motion with the rate V, the third term the heat exchange across the side surface with the ambient medium. Several reasons are decisive in selecting this equation for a description of heat transfer in a crystal—melt system. In contrast to three-dimensional equations, eq. (8) has simple solutions which furnish a good fit to the true temperature distribution for small Biot =

.

S

S

‘

numbers (Bi (1s/v)R ~ 1). Nevertheless, eq. (8) gives a qualitatively correct description of GL and =

V.A. Tatarchenko

/

Effects of capillary shaping on crystalgrowth

G~as functions of the relevant parameters even for not too small Biot numbers (Bi ~ 1). The temperature at the interface is equal to 1~,and if account is taken of the surface kinetics, we may assume a relationship between the growth rate and supercooling. The remaining boundary conditions depend on the specific growth tecimique and conditions. In growing single crystal of optically transparent materials, the radiation flux is taken into account by introducing an effective thermal conductivity coefficient proportional to T3. If the heat transfer from the lateral surface of the crystal proceeds by radiation according to the Stefan—Boltzmann law, we arrive at an equation of the type of (8) with T replaced by T4. The result obtained can therefore be applied to this case if T is replaced by T”4. When opaque materials with high melting points are crystallized, linearization of radiation heat transfer and introduction of an effective heat exchange coefficient for the lateral surface enable one to use eq. (8) up to 2000°C.

4. The effect of surface tension gradient on formation of the meniscus and heat transfer processes in the liquid phase Normally, crystallization is conducted with a temperature gradient in the liquid phase. This may result in convection and, in particular, free convection involving gravity and characterized by the Rayleigh number Ra=

g/3z~tTL3 ,

coefficient. The number L2 ~ g$ A obtained as the ratio of the Rayleigh and Marangoni numbers, characterizes the relative contributions of gravitational and capillarity forces to the formation of convective flow. This number is small for small characteristic dimensions of the system when capillanty forced predominate. It was proved experimentally [6] that it is the Marangoni convection which develops in a zone with characteristic dimension L 5 mm. The characteristic dimensions in the capillary shaping techniques are normally still smaller; therefore we discuss the effects of only this convection. If the melt is superheated, the surface tension is increased at the crystallization front and, therefore, the melt is transported toward the front at the free meniscus surface, and away from it at the centre because of the pressure gradient generated. This pressure gradient is directed against that produced by the gravity field, and hence increases the convexity of the molten columns. The pattern is reversed in supercooled melts. It was mentioned, however, that flow affects the meniscus shape appreciably beginning at flow rates of 10_102 cm/s. The actually observed flow rates in ref. [6] are about 1 cm/s, therefore the effect of the above factors on the meniscus shape will be insignificant. Convection in the liquid phase results in melt stirring and so increases the heat influx to the crystallization front from a superheated melt. The magnitude of these flow rates can be estimated from the boundary condition at the meniscus free surface: -~

dV =

and Marangoni convection involving the surface tension gradient and characterized by the Marangoni number Ma

=

8yz~TL ~ ~

Here ,6 is the bulk expansion coefficient, z2iT is the characteristic temperature difference, L is a characteristic linear dimension, i~ is the dynamic viscosity coefficient and y is the surface tension

77

—

~grad

T,

and thus V

~

~

8T h’q S

where R is the crystal radius (the charactenstic transverse dimension of the crystal), and h is the meniscus height. The Marangoni number gives an estimate of the ratio of the heat flux transferred by the melt flow at the flow rate V to that transferred by diffusion. The estimates given by

78

/

V.A. Tatarchenko

Effects of capillary shaping on crystal growth

the above formula exceed the observed values by an order of magnitude [6]. This discrepancy is caused by friction on the surface of the crystal. It was also mentioned in the cited reference that the convective contribution to heat transfer becomes important if zones are longer than 1 mm. Note that the convective heat flux to the crystallization front increases as R rises and h diminishes; this is clear from the estimate of the flow rates. These changes operate in the same direction as the corresponding changes in the heat flux transported by diffusion. In other words, the thermal coefficients, which are of interest in an analysis of stability, change their magnitudes but do not alter their signs. Summarizing, the effect of convection on heat transfer can be taken into account by introducing formally the coefficients of the effective heat conductance and heat exchange from the lateral surface.

5. Stability of the shape of the crystallization front The crystallization front was assumed plane in the preceding considerations, and only its coordinate was analyzed for stability. In the general case, an arbitrary perturbation of the front can be treated as a superposition of the shape and coordinate perturbations. Formally, this situation could be treated by expanding the perturbation in a series of an appropriate set of orthogonal functions, with the first term taking into account the change of the crystallization front coordinate and the remaining terms representing shape distortion at a constant front coordinate. All the terms of the expansion are regarded as independent, since this S

S

S

S

analysis of stability is restricted to a linearized problem. So far, only the perturbation of the first term has been considered. Perturbations of the front shape at constant front coordinate (perturbation in all the terms in the expansion excepting the first one) do not change the capillarity problem as discussed above, so that the shape stability is determined by the thermal effects only; this was investigated by Mullins and Sekerka [I] and by others. It was found that, with impurity effects neglected, a plane front is stable in growth from a superheated melt, while the instability ap-

pearing in the case of growth from a supercooled metal can be compensated for, at low supercooling, by the Gibbs—Thompson effect (the dependence of the crystallization temperature on the surface curvature) and the effect of kinetics at the crystallization front. Qualitatively, the stability behavior is reproduced if the initial front shape is not planar. The stability of the crystallization front is determined by a joint action of the capillanty and thermal conditions of the crystallization process. Obviously, the magnitudes of the thermal coefficients must change but the signs will not be altered. The stability of the front shape must be analyzed in a manner similar to that of ref. [1].

6. Czochralski and Stepanov techniques In Czochralski and Stepanov growth, two parameters can be varied independently in the first approximation: the crystal size and the position of the crystallization front (a two-degrees-offreedom system). The linearized equations (5) and (7) are conveniently written as ~ =ARR 6R +ARh ~h, (9) ~ A~ 81? where A

=

—

~

~

=

(10)

A

BR

RR

A

6h,

+ Ahh

=

=

A

~

z BR

—

V-~—-~

Rh

‘

=

z

hh

‘

S

S

‘

‘

The necessary and sufficient conditions of stability of the steady-state solutions R0 and h0, according to Gurvitz s criterion, are: .

-

.

.

ARR + Ahh <0, A A A A hh

Rh

(11) (12)

> ~

—

RR

,

hR

-

The results of an analysis of A,k coefficient signs are tabulated (see table 1). The process kinetics (crystallization temperature as a function of growth rate) were taken into account, and showed that only small changes were caused by a decrease of Ahh. If the melt is slightly supercooled with respect to the equilibrium temperature, this effect may

V.A. Tatarchenko

/

Effects of capillary shaping on crystal growth

79

Table 1 Coefficients A

5 which determine the crystallization stability in the Czochralski and Stepanov techniques

Crystallization Conditions

ARR

Czochralski growth smalldiameter crystal

o large diameter crystal R>a or a plate

0.

o

C’)

>0

boundary catching condition

>

ARh

r0

~o

~.~0Ro singlevalued miiiisci

boundary catching condition

_________

<0

<0 —

—

<0 r0-R0>m a0>a1

—

two>0 valued ~5 menisci atR0r0

<0

boundary wetting condition

<0

<0

0 r0- R0< m

a°
lead to an actual superheating of the melt, that is, to a transition to the thermal stability region. The results obtained enable one to derive quantitative estimates which confirm applicability of the quasi-stationary approach in treating thermal phenomena when stability is analyzed. It was mentioned above that this approach is justified if the characteristic time of the crystallization front relaxation to the steady state, ‘r,~,is greater than the corresponding relaxation time of temperature ~ From the solution of eq. (8), 1

T~f~ AhhI

_________

P(Tm~

)•

(13)

~T

can be estimated as

proportional to the characteristic relaxation time) is higher for materials with higher thermal conductivity A and growth angle 4~and with smaller latent heat of melting The stability can be increased by raising the melt superheating and by forced air cooling of the crystal in the crystallization zone. Not only is the stability increased, but the range of stable growth is widened as well in these cases. It is interesting to analyse how the type of stability changes with the changing size of the .~.

3 2

As follows from 2/sc, so thatref. [4], TrT h

(‘5

T

(14)

1

where CT is thermal capacitance. It is also of interest to estimate the effect is various thermophysical and other material constants on the crystallization stability. We find that, with other conditions being identical, the stability (taken to be

0

1T

pCT(Tm—

~

_

102_103

>>

1,

0

0.2

0.4

0.6

0.8

1

R0/r0 Fig

3

Stability index, S (s

I),

as a function of crystal

radius/shaper radius ratio, R0/ro, for three materials: (1) Al, (2) Si, (3) Al203, grown under identical conditions.

80

V.A. Tatarchenko

/

Effects of capillary shaping on crystal growth

crystals grown. For instance, when thin rods are grown by the Stepanov technique from a superheated melt, an increased rod diameter (this can be achieved, for example, by diminishing the pulling rate) changes the stability type in the following sequence: “saddle” type stability —s “unstable node” “unstable focus” stable focus “stable node”. The characteristic relaxation times in a steady-state mode are, for different materials, in the range of 1—10 s. When the stability boundary is approached, the relaxation times increase. In principle, an optimal R0/r0 ratio can be determined on the basis of such calculations. Fig. 3 plots the stability index (S) as a function of —*

—

this ratio at a fixed thermal regime for three materials (Al, Si, A12O3).

References [1] W.W. Mullins and R.F. Sekerka, J. Appl. Phys. 35 (1964)

—~

[2] A.V. Lykov, Teploobmen (Energiya, Moscow, 1975). [3] V.V. Voronkov, Fiz. Tverd. Tela 5 (1963) 571. [41G. Karlsrow and D. Eger, Teploprovodnost Tverdykh Tel (Nauka, Moscow, 1964) p. 385. [5] VA. Tatarchenko, J. Crystal Growth 37 (1977) 272. [6] D. Schwabe, A. Scharmann, F. Preisser and R. Oeder, J. Crystal Growth 43 (1978) 305,