Dynamic three-dimensional simulation of facet formation and segregation in Bridgman crystal growth

(1)

Journal of Crystal Growth 303 (2007) 287–296

Dynamic three-dimensional simulation of facet formation and

segregation in Bridgman crystal growth

C.W. Lan

, C.J. Chen

Department of Chemical Engineering, National Taiwan University, Taipei 10617, Taiwan, ROC Available online 30 December 2006

Abstract

Facet formation is an important kinetic phenomenon during crystal growth and its simulation remains a non-trivial task in crystal-growth modeling. Numerical models for facet formation have been developed, and its dynamic three-dimensional (3D) simulations, coupled heat ﬂow and segregation in a Bridgman growth process are presented. The growth of yttrium aluminum garnet (YAG) crystals in the [1 1 1] direction having {1 1 1}, {2 1 1}, and {1 1 0} facets are taken as examples. Two numerical schemes for the interface tracking based on the geometric and kinetic models are considered, and the calculated results are found in good agreement. By using different segregation coefﬁcients at the facets and the rough surface, both radial and axial abnormal segregations, such as the dark core, due to facet formation are further simulated. Multiple facets are considered as well.

PACS: 44.25.+f; 47.27.Te; 81.10.Fq; 02.60.C6; 02.70.Fj

Keywords: A1. Facet; A1. Kinetics; A1. Segregation; B2. Bridgman method

1. Introduction

Facet formation is an important kinetic phenomenon during single-crystal growth. Especially, the segregation and twinning phenomena are associated with faceting

[1–3]. Facets often appear at low-index crystallographic planes having a lower surface energy. For oxides, such as yttrium aluminum garnet (YAG) and gadolinium gallium garnet (GGG), having the growth in the /1 1 1S direction, faceting at {2 1 1} crystallographic planes is typical [4,5]. Other facets in {1 1 1} and {1 1 0} were also observed. Faceting can affect the crystal quality in several ways. The most common ones are the stress concentration at the singular edges and the severe segregation due to the different segregation coefﬁcients in the facets and the rough surface[1]. The dark core in the growth of Bi12Si20O

(BSO) is a typical example of the ion segregation due to faceting[6].

The modeling of facet formation dates back to the early 1970s [7,8] using a simple geometric approach, assuming the thermal gradient at the interface is known. The ﬁrst attempt to simulate the faceting coupled with the heat transfer in oxide growth was made in Brandon’s group

[9,10] using a two-dimensional (2D) finite element ap-proach. Recently, Weinstein and Brandon [11,12] also developed algorithms for faceting, taking into account the coupling of different kinetic mechanisms and extended them to 3D simulation as well[13]. Nevertheless, melt flow and segregation were not considered in their models because the faceting, except {1 1 1}, breaks axisymmetry and a 3D simulation is required. Lan and Tu[14]presented for the first time the 3D simulation of {2 1 1} facet formation in a Bridgman YAG growth, where a pseudo steady-state model was used [14]. Melt convection, heat transfer, and dopant segregation were simulated by using a finite volume method. The abnormal segregation due to the different segregation coefficients at the facets was presented.

In viewing the numerical modeling of the facet forma-tion, a simple geometric approach, similar to the one www.elsevier.com/locate/jcrysgro

Corresponding author. Tel./fax: +886 2 2363 3917. E-mail address:[email protected] (C.W. Lan).

(2)

proposed in the classic model [7,8], was used by Lan and Tu[14]. Ma et al.[15]also used a similar approach for the 2D growth of YAG crystals; however, dynamic formation of the facets was further considered. On the other hand, Weinstein and Brandon [11,12] used a kinetic model approach that moves the interface shape according to the local kinetic models related to crystallographic orientation and undercooling. Different mechanisms for facet growth could be considered.

In this paper, the model proposed by Lan and Tu[14]is extended to time-dependent simulation having different facets, so that the dynamic facet formation can be simulated. Radial and axial dopant segregations, as well as dark core formation, during crystal growth can also be simulated. Besides, multiple facets are further considered. For comparison purposes, the kinetic model is also implemented and the facet formation through different mechanisms is illustrated. In the next section, the numerical models are brieﬂy described. Section 3 is devoted to results and discussion, followed by conclusions in Section 4. 2. Mathematical model and numerical solution

A generic Bridgman crystal growth conﬁguration is shown in Fig. 1a. Due to the faceting at the interface (or imperfect growth environment), the problem is 3D and

thus described by Cartesian coordinates (x, y, z) with the origin (ﬁxed and stationary) at the sample bottom. The furnace or ambient thermal distribution Ta(x, t) is assumed

to be known. Both ampoule and the furnace thermal proﬁle are allowed to move. However, the ampoule translation speed Uais used only for pseudo-steady calculations, while

the translation of the thermal proﬁle at the speed Uhis used

for time-dependent growth simulation. Again, if the pseudo-steady state is assumed, the magnitude of the ampoule pulling speed Ua is assumed to be the same as

the axial growth rate. The flow, temperature, and dopant fields, as well as the melt/crystal interface, are also represented in the Cartesian coordinate. Furthermore, it is also assumed that the dopant concentration is so low that its effects on the flow and liquidus temperature can be neglected. Solid-state diffusion is also not considered and the melt is assumed to be opaque [16]. To represent the governing equations in dimensionless forms, the dimen-sionless variables are defined by scaling length with the crystal diameter Dc, time with Dc2/am, velocity with am/Dc,

pressure with rmam/Dc2, temperature with the melting point

Tm, and dopant concentration with its initial concentration

C0in the melt, where amis the thermal diffusivity and rm

the melt density. For the convenience of discussion, the variables used in the equations are all dimensionless unless otherwise stated. Based on the Boussinesq approximation,

(3)

the governing equations in dimensionless form for the transport processes in the melt during crystal growth can be described by the conservation laws for the mass, momentum, energy, and dopant as follows[17]:

r v ¼ 0, (1) qv qtþv r v ¼ rP þ Prr 2_{v Pr Ra} TðT 1ÞeG, (2) qT qt þv rT ¼ r 2_T, ₍₃₎ qC qt þv rC ¼ Pr Scr 2_C. ₍₄₎

In the above equations, v, t, P, T, and C are the dimensionless velocity, time, pressure, temperature, and dopant concentration, respectively. Also, eG¼ cos (g)ex

sin (g)eyis the unit vectors in the gravitational direction and

g is set to zero in this report; exand eyare the unit vectors

in the x and y directions, respectively. The associated dimensionless thermal Rayleigh number is deﬁned by RaT¼

b_Tg₀TmD3c

nmam

,

where bT is the thermal expansion coefﬁcient, g0 the

gravitational acceleration, and nmthe kinematic viscosity.

In addition, Pr ¼ nm/am is the Prandtl number and

Sc ¼ nm/D the Schmidt number, where D the dopant

diffusivity in the melt.

In the crystal (c) and the ampoule (a), only heat transfer needs to be considered, and their formulations and related discussion could be found elsewhere[14]. In short, for most oxides, except being highly doped, the crystal is semitran-sparent. Accordingly, internal radiation cannot be ignored. Detailed calculation of internal radiation is complicated and computing-extensive. Even with the P1approximation,

the computational effort is not trivial [18]. However, to reduce the computational cost and to concentrate the discussion on the facet formation, the internal radiation inside the crystal is treated by the Rosseland diffusion approximation [19] using the so-called effective thermal conductivity.

The thermal boundary condition at the melt/crystal interface is set by the heat ﬂux balance:

Qj_mQj_cþg_cSt vcþ

qhc

qt

exn ¼ 0, (5)

Q|m and Q|care the dimensionless total heat ﬂuxes at the

melt and the crystal sides, respectively. Also, gc is the

density ratio of the crystal to the melt and vc the crystal

moving speed; the crystal moving speed is equal to the ampoule translation speed Ua. The Stefan number

St ¼ DH/CpmTm scales the heat of fusion DH released

during solidiﬁcation to the sensible heat in the melt; Cpmis

the speciﬁc heat of the melt and n the unit normal vector at the interface pointing to the melt. The mass ﬂux balance at the interface can be written as the following, where the solid-state diffusion is neglected:

n rCjmþ Sc Prg_cð1 KÞC vcþ qhc qt ðexnÞ ¼ 0, (6)

where K is the segregation coefﬁcient. Since we allow different segregation coefﬁcients at the faceted interface and the rough surface, K ¼ Kf is set inside the facets;

typical K is used for the rough surface. For simplicity, the furnace temperature Ta(x, t) is assumed to be a hyperbolic

tangent function and the crystal growth can be carried out by moving the thermal proﬁle upwards at the speed Uh(as

illustrated inFig. 1a). The stationary condition at t ¼ 0 is deﬁned at the middle of the system, i.e., Ta(L/2, 0) ¼ Tm.

Other boundary conditions are similar to the ones reported before [14].

The treatment of temperature or growth speed at the melt/crystal interface considering the kinetic effect needs a special attention. During facet formation, the interface temperature is not at the melting point Tm, so that the

Stefan condition (Eq. (5)) alone is not adequate for the interface speed (or position) and temperature. Therefore, an additional equation, relating the normal interface speed (Vn) to the undercooling (DT), is needed. The simplest

kinetic model (dimensional) for this at the interface can be written as

Vn¼bðy; DTÞ DT, (7)

where b is the kinetic coefﬁcient, which can be a function of the misorientation angle y and the interfacial undercooling DT. The misorientation angle is the angle between the normal vectors of the interface and the low-index plane, where the low-index plane eventually develops into a facet. For a rough surface, b ¼ broughcan be regarded as a

constant [20]. However, for facet formation, different mechanisms may be involved. Usually, a new facet plane can be formed through 2D nucleation or step motion initiated from defects, such as the dislocations and twins. For step motion, the kinetic coefﬁcient can be represented by bSM¼BSM|sin (y)|, where BSM is the

step kinetic coefﬁcient; y is usually less than 11 for facet growth [21]. The step (spiral) growth can also take place through dislocations. In such a case, the kinetic coefﬁcient can be written as bDG¼BDG cos (y)DT; BDG

is a kinetic constant. On the other hand, without the defects to initialize the growth, 2D nucleation has to kick in and the kinetic coefﬁcient can be written as b2DN¼B2DNcos (y) exp (A2DN/DT)[22,23].

Eqs. (5) and (7) can be coupled to obtain the interface temperature and speed (or interface position); the normal interface velocity Vn¼(qhc/qt)(exn). Nevertheless, the

coupled solution with heat ﬂow calculations requires a robust numerical scheme. A simple way to do this is the geometric model approach [7,8,14,15]. During iterations, once the interface is obtained from the isotherm at Tm, DT

(4)

is calculated from the kinetic equation (Eq. (7)). Then, the interface is cut by the given faceting plane, such that Eq. (7) is satisﬁed at the facet center having the largest under-cooling. This approach is rather simple from the numerical implementation point of view. However, the facet size needs to be larger than the size of grid cell. Otherwise, this can cause divergence during iterations. For transient simulation, this is indeed a problem. If the initial time step is too small, the supercooling is not large enough to form a facet that is large enough to cover a grid cell, the scheme may fail.

Another numerical scheme, the kinetic model approach, is similar to the one proposed by Weinstein and Brandon

[11–13]. It uses Eq. (7) directly by marching the new interface position (hcnew) explicitly from the old one (hcold):

hnew_c ¼hold_c þDtbðyold; DToldÞDTold=ðexnÞ, (8)

where Dt is the time step; both the rough and facet surfaces are treated by the same equation, but having different kinetic coefﬁcients. This approach allows a small under-cooling to be considered, so that growth mechanisms having a small undercooling could be easily incorporated without the restriction of the grid size. In fact, the facet formation through different mechanisms has been simu-lated successfully by Weinstein and Brandon [11–13]. Nevertheless, to precisely represent the facet, the cell size still needs to be large enough. Also, to insure numerical

stability, the time step size (Dt) during the explicit time integration using Eq. (8) cannot be too large.

In reality, since a facet should not grow faster than a rough surface, the interface movement needs to be selected from the smaller velocity, i.e., Vn¼min[Vn(rough), Vn(facet)].

However, once the facet is formed, the fastest growth mode dominates in the facet, i.e., Vn(facet)¼max[Vn(DG), Vn(2DN),

Vn(SM)]. Therefore, the selection of the kinetic coefﬁcient

can be formulated as follows[11–13]:

bðy; DT Þ ¼ min½b_rough; max½b_DG; b_2DN; b_SM. (9) Based on this selection rule, Eq. (7) can be properly used for simulation; if the interface speed is known, the undercooling can be calculated. For example, as the growth starts by moving the heating proﬁle upwards, the undercooling DT is small, as are b2DN and bDG.

As a result, the step growth is dominant in the facet formation. Accordingly, the kinetic coefﬁcient b(y, DT) ¼ min[brough, BSM| sin (y)|] is used to move the interface.

Because at the beginning there is no facet, the interface moves as a rough surface except at the faceted direction having a small y; also bSMpbroughneeds to be satisﬁed. As

the growth proceeds further, because the facet grows at a lower speed, the facet size continues to increase. Mean-while, the undercooling continues to build up inside the facet until its speed catches the rough one or the next growth mode kicks in.

Fig. 2. Simulated facets and the undercooling on the interface using the pseudo-steady-state model; from left to right are {1 1 1}, {1 1 0}, and {2 1 1} facets, respectively. The mesh is also shown on the interface. DT ¼ 2 K and Ua¼ 1 cm/h.

(5)

Once the kinetic information is known, to implement the facet calculation is straightforward by settling a geometric relationship for the crystallographic planes. Take the growth in the [1 1 1] direction having {2 1 1} facets as an example. We first align the [1 1 1] direction with the x-axis. By a simple rotation transformation, ½¯1 0 1 is in the z-direction, while ½¯1 2 1 is in the y-direction. Furthermore, the plane normal of the three {2 1 1} facets can be easily obtained. For other facet planes, they can be treated in a similar way. In addition, because the interface is defined by the nodes, the vertices of the mesh are then interpolated from the node values. Both nodes and vertices are used for remeshing. This procedure is robust for interface move-ment, and no special treatment is needed after finding the interface shape.

The above governing equations and their associated boundary conditions are solved by an efﬁcient ﬁnite volume method (FVM) scheme with multigrid acceleration

[24] for the moving boundary problem. Three meshes are considered in our simulations. In the first level (mesh M1), there are 16 21 30 (in the radial, angular, and axial directions, respectively) finite volumes in the melt and 16 21 15 in the crystal, and 5 21 45 in the ampoule. The second level (mesh M2) doubles the finite volumes in each direction, i.e., 32 42 60 cells in the melt. The third mesh (M3) doubles the grid in the radial and angular directions only, but not involved in the multigrid calcula-tions. A sample of mesh M2 is shown in Fig. 1b for reference. The calculations are performed in personal computers (Pentium/2 GHz CPU with 1 G RAM), and one calculation takes about 4 h of CPU time. For most of the cases in this report, the results based on these three meshes are about the same. Therefore, the results presented

Fig. 3. (a) Evolution of the interface shape having {2 1 1} facets and (b) the interface moving speeds at different positions with and without {2 1 1} facets. The dashed line in (a) is the interface shape at t ¼ 1200 s without facet formation.

Fig. 4. (a) Evolution of the interface moving speed at different facets and (b) evolution of the maximum undercooling at different facets.

(6)

here are based on mesh M2. A detailed description of the numerical method can be found elsewhere[24].

3. Results and discussion

For comparison purposes, we again take the growth of YAG in a molybdenum crucible [13] as the examples for study; the growth direction is /1 1 1S. Most of the simulation parameters are the same as the ones used previously[14](the Rosseland mean absorption coefﬁcient aR¼5 cm

1

), except the kinetic parameters; brough¼4

102cm/(s K), BSM¼2 10 3

cm/(s K), A2DN¼149 K,

B2DN¼8.13 1013cm/(s K2), and BDG¼0.3–2 104cm/

(s K2)[11–13]. These kinetic data are used in the simulation unless otherwise stated.

Before showing the results from time-dependent simula-tion, we ﬁrst illustrate three as-grown facets and the mesh on the interface based on the pseudo steady-state simula-tion [14]. The results are shown in Fig. 2 as having a maximum undercooling of 2 K; Ua¼ 1 cm/h. The

inter-face shapes on the x–z plane are further presented for

comparison. As shown, with the same undercooling, the {1 1 1} facets are larger than {2 1 1} facets. The {1 1 0} facets are the smallest. This is consistent with the prediction by the classical model[7], where the facet size b can be represented by

b ¼pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi8 DTR=G, (10) where R is the radius of curvature of the melting point isotherm and G the local thermal gradient. The radius of the isotherm curvature near the centerline, where the {1 1 1} facet is located, is the largest, so that the facet there is larger. On the other hand, near the ampoule, due to the smaller radius of curvature of the isotherm, the {1 1 0} facets are smaller. The thermal gradient G is about the same along the interface. With the basic geometric under-standing of the three facets, dynamic simulations are further presented. In the dynamic simulation, the kinetic coefficients are as just mentioned and the translation speed of the thermal profile is set at Uh¼5 mm/h.

The formation of the {2 1 1} facet at y ¼ 0 is shown in

Fig. 3a; the time interval is 200 s. Indeed, at the beginning,

Fig. 5. Effect of the segregation coefficient at the facet (Kf) on the flow and dopant fields: (a) Kf¼0.05, (b) Kf¼0.1, and (c) Kf¼0.2; at the rough surface,

(7)

the growth speed at the facet is the slowest. Its speed increases as the undercooling increases as a result of the heater movement. After 1200 s, a steady-state growth is reached and the interface shape remains about the same. The calculated interface shape (dashed-line) without facet formation is also plotted in Fig. 3a. As shown, the difference with the one having facet formation is not much except in the faceted region. In fact, this is consistent with that reported in Ref.[14]. The effect of facets on the overall heat transfer is not signiﬁcant. The dynamic interface speeds at different positions with and without {2 1 1} facets are shown in Fig. 3b. Again, the interface speed at the facet center is slower than others. The moving speed of the rough surface at the centerline is also not affected much by the existence of the facets. Without the facets, the interface speeds at different positions are also quite close.

The comparison of the time evolutions of the interface moving speeds at different facets is shown in Fig. 4a. As shown, the {1 1 1} facet, which develops into a bigger facet, moves slower than the others, while the {1 1 0} facets move faster. The comparison of the evolution of the largest undercooling at different facets is shown in Fig. 4b. As shown, the undercooling at the {1 1 1} facet, which moves slower, also builds up slower; this can be simply predicted from the kinetic equation, i.e., Eq. (7), using the dislocation growth model.

As observed in the experiments, the segregation is affected by facet formation due to different segregation coefﬁcients at the facets (Kf) and the rough surface (K)[1].

To illustrate this,Fig. 5shows the effect of Kfon the ﬂow

and segregations at t ¼ 4000 s, while the segregation coefﬁcient at the rough surface remains at 0.1. The Kf

value is set to 0.05, 0.1, and 0.2 in Fig. 5a–c, respectively.

Fig. 6. Dopant segregations in the grown crystal: (a) {2 1 1}, (b) {1 1 1}, and (c) {1 1 0} facets; Kf¼0.2 and K ¼ 0.1. The left are the segregations on the x–z

(8)

Fig. 7. (a) The evolution of the kinetic coefﬁcient b at the interface using the kinetic model approach, (b) comparison of the interface shapes at different growth periods simulated by two numerical models, and (c) comparison of dopant segregations at different growth periods simulated by two numerical models.

(9)

As shown, the calculated flow and dopant fields in the melt are not altered much by varying the segregation coefficients on the facets. Nevertheless, near the interface, due to the different segregation values, the dopant distribution is significantly affected. The dopant concentration in the facets having a lower Kf (0.05) is lower, while the melt

concentration is higher (Fig. 5a). On the other hand, for the higher Kf value (0.2), the dopant inside the facets is

enriched. For YAG, due to its large viscosity, the convection in the melt is weak. Nevertheless, due to the large Sc value (1860), the dopant field is significantly affected by the melt convection. The thermal convection near the interface pushes the dopant from the center to the periphery of the interface. The zigzag dopant concentration fields at the facet/rough boundary are due to the finite grid size. The zigzag becomes smaller as the mesh is further refined. Nevertheless, we have found that mesh M2 is adequate, in terms of accuracy, for dopant calculations. The result obtained from mesh M3 gives a smoother boundary, but the dopant distribution is not changed much.

Because dynamic simulation is carried out, both radial and axial segregations in the grown crystal can be calculated. Fig. 6a shows the concentration ﬁelds on the x–z plane at y ¼ 0.3 cm (left) and on the y–z plane at x ¼ 4.5 cm (right) for Kf¼0.2 (the same as Fig. 5c),

respectively. As shown, the cores develop clearly due to facet formation. In the simulation, because the transient period is short, being about 500 s, the cores also develop quickly. Similar simulations are also carried out for {1 1 1} and {1 1 0} facets, and the dopant ﬁelds on the x–z and x–y

planes are shown inFigs. 6b and c, respectively. Again, as shown, the dopant segregations are affected by the facets signiﬁcantly. Because {1 1 0} facets appear near the ampoule wall where the dopant concentration in the melt is the highest, as a result of thermal convection, the dopant concentration in the {1 1 0} facets is also much higher than that in the {2 1 1} and {1 1 1} facets.

The kinetic model approach is further used for similar simulations. As mentioned previously, the main advantage of using this approach is that the mechanisms having low supercooling at the beginning can be adopted. In addition, inside a facet, different mechanisms can be incorporated for facet formation as well. However, to insure the numerical stability, the time step needs to be much smaller than that used in the geometric model (Dt ¼ 40 s). We have found that Dt ¼ 0.1 s is good enough for most of the simulations for the kinetic model approach using the explicit scheme described by Eq. (8). Fig. 7a shows the evolution of the kinetic coefﬁcient b for the same growth condition, except the kinetics of the rough surface and step motion mechanisms are considered. The selection of b is based on Eq. (9). At the rough surface, brough¼2

102cm/s (s K) and the value is ﬁxed. However, inside the facet, the kinetic value changes with the misorientaion angle y and the undercooling DT dynamically. In Fig. 7a, the step motion is dominant at the initial growth and near the rim of the facets, where the undercooling is smaller. On the other hand, the dislocation growth is dominant at the facet center due to the larger undercooling, especially after the initial stage. Once the dislocation growth dominates and the undercooling is established, the kinetic coefﬁcient

(10)

approaches to a constant.Fig. 7bshows the comparison of the interface shapes on the x–y plane from two different numerical models; in the geometric model, only the dislocation growth mechanism is considered. As shown, they are in reasonable agreement. The comparison of the dopant ﬁelds at the interface (crystal side) along the z-axis (y ¼ 0.2 cm) from two numerical models is shown in

Fig. 7c. As shown, the difference in the transient period is larger due to the different facet size. However, once the facets develop, both models give similar results.

Multiple facets can be easily considered for both numerical approaches. Fig. 8 shows the evolution of the dopant distribution at the interface considering {2 1 1} and {1 1 0} facets; the kinetic model approach is used for simulation. As shown, the highest dopant concentration appears at the edge of {1 1 0} facets, and this is due to the effect of thermal convection that pushes the dopant from the center to the rim of the interface, as well as the larger segregation coefﬁcient (Kf¼0.2) considered in the

simula-tion.

4. Conclusions

In this study, we have presented the dynamic 3D simulation of facet formation for Bridgman YAG crystal growth. The coupling of the facet formation with the heat flow and dopant fields is considered. In addition to {2 1 1} facets, the formation of {1 1 1} and {1 1 0} facets is also illustrated. Two numerical schemes, based on the geometric and kinetic model approaches, are adopted and compared. If we do not pay much attention to the low undercooling, the geometric approach is quite efficient for facet simula-tion because the time step for integrasimula-tion can be much larger. However, the kinetic model approach allows the incorporation of the low-undercooling mechanisms, which are important at the initial stage of growth. Different

growth mechanisms can be considered in the facet formation as well. Accordingly, a more realistic simulation can be carried out. Through the 3D dynamic simulation of facet formation, both radial and axial segregations in the crystal can be predicted, and this should be useful for the control of crystal uniformity.

References

[1] K.F. Hulme, J.B. Mullin, J. Electron. Control 3 (1957) 160. [2] D.T.J. Hurle, J. Crystal Growth 147 (1995) 239.

[3] J. Amon, F. Dumke, G. Mu¨ller, J. Crystal Growth 187 (1998) 1. [4] B. Cockayne, J. Crystal Growth 42 (1977) 413.

[5] A.G. Petrosyan, Kh.S. Bagdasarov, J. Crystal Growth 34 (1976) 110. [6] X. Xu, J. Liao, B. Shen, P. Chao, X. Chen, C. He, J. Crystal Growth

133 (1993) 267.

[7] J.C. Brice, J. Crystal Growth 6 (1970) 205. [8] V.V. Voronkov, Krystallograﬁya 17 (1972) 909.

[9] Y. Liu, A. Virozub, S. Brandon, J. Crystal Growth 205 (1999) 333.

[10] A. Virozub, S. Brandon, Model. Sim. Mater. Sci. Eng. 10 (2002) 57.

[11] O. Weinstein, S. Brandon, J. Crystal Growth 270 (2004) 232. [12] O. Weinstein, S. Brandon, J. Crystal Growth 268 (2004) 299. [13] O. Weinstein, S. Brandon, J. Crystal Growth 284 (2005) 235. [14] C.W. Lan, C.Y. Tu, J. Crystal Growth 233 (2001) 523.

[15] Y. Ma, L.L. Zhang, D.J. Larson Jr., J. Crystal Growth 266 (2004) 257.

[16] Brandon, J.J. Derby, J. Crystal Growth 121 (1992) 473. [17] C.W. Lan, I.F. Lee, B.C. Yeh, J. Crystal Growth 254 (2003) 503. [18] C.W. Lan, Int. J. Heat Mass Transfer 46 (2003) 1629.

[19] D. Vizman, I. Nicoara, G. Mu¨ller, J. Crystal Growth 212 (2000) 334. [20] D. Turnbull, Thermodynamics in Metallurgy, American Society of

Metals, Metals Park, OH, 1949.

[21] V.V. Voronkov, Crystals grown, properties, and applications, in: A.A. Chernov (Ed.), Modern Theory of Crystal Growth I, Vol. 9, Springer, New York, 1983.

[22] M. Volmer, Kinetik der Phasenbildung, Steinkopf, Dresden and Leipzig, 1939.

[23] R. Becker, W. Do¨ring, Ann. Phys. 24 (1935) 719. [24] C.W. Lan, M.C. Liang, J. Comp. Phys. 152 (1999) 55.