Analysis of Melt Flow Mixing in Czochralski Crystal Growth Process

Article pubs.acs.org/IECR

Analysis of Melt Flow Mixing in Czochralski Crystal Growth Process Mojtaba Izadi and Stevan Dubljevic* Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta 26G 2V4, Canada ABSTRACT: The prime example of a relevant industrial process in which fluid flow dynamics impacts desired process specifications is given by the Czochralski (CZ) crystal growth process. Melt flow in the crucible induces a variety of transport phenomena with profound effects on the mass, momentum, and heat transfer during the process. To provide novel insight into the transport phenomena associated with the CZ process melt flow, the approach of discerning the mixing templates by identifying Lagrangian coherent structures (LCS) is explored. In particular, the LCS have been extracted to identify transport features in the melt during the crystal growth. The finite-element method has been employed to simulate realistic axisymmetric melt flow accounting for physically relevant geometries (including the melt−crystal and melt−ambient interfaces) for identification of the LCS. These structures are known to represent dynamically driven material surfaces, and they provide insight into mixing properties of the melt flow that influence the characteristics of the grown crystal.

1. INTRODUCTION The Czochralski (CZ) crystal growth process is the most important method for manufacturing large-scale silicon crystals for use in the semiconductor and electronics industries.1 In this method, a crystal rod is pulled out vertically at a rate of a few centimeters per hour from the surface of a heated pool of melt contained in a crucible, as shown schematically in Figure 1. The

well as radiation. Sackinger et al. improved the TCM to generate an algorithm, called hydrodynamic TCM (HTCM), that includes convection in the melt.5 A comprehensive introduction and review of other relevant research prior to 1997 is given by Prasad et al.6 Liu and Kakimoto developed a three-dimensional global simulation of the CZ process using a mixed two-/threedimensional discretization technique7 in which the melt and crystal are discretized in three dimensions, whereas other domains including pulling rod, crucible, pedestal, heater, and heat insulators are discretized in two dimensions. Turbulent momentum and heat transport in an idealized CZ process were investigated by Wagner and Friedrich by means of direct numerical simulation.8 They focused on the influence of crystal and crucible rotation on the flow structure and development of temperature fluctuations. Three-dimensional unsteady analysis of the melt turbulent convection was presented by Kalaev et al., who included the three-dimensional calculation of the melt−crystal interface and incorporated two-dimensional defect evolution in the crystal.9 A hybrid finite-volume/finite-element method was used by Lipchin and Brown to model the heat transfer, turbulent convection, and dopant transport in industrial-scale Czochralski crystal growth systems.10 Abbasoglu and Sezai carried out a series of unsteady three-dimensional numerical simulations to examine the melt flow during the CZ process considering the effect of convection driven by surface tension on the free surface of the melt.11 It has been demonstrated that, during the formation of a crystal from its melt, mixing properties of the melt flow affect transfer phenomena.12 If complete mixing in the melt does not occur, a gradient of impurity concentrations exists in the melt near the crystal interface.13 This will affect the crystal characteristics such as compositional uniformity of the grown crystal.

Figure 1. Czochralski crystal growth process setup.5

thermal field in the grown crystal, the melt−crystal interface shape, and the melt flow in the crucible during crystal growth have significant effects on the crystalline structure and formation of microdefects, which influence the quality and properties of the grown crystal. Over the past few decades, many studies on the dynamics and characteristics of the CZ crystal growth process have been performed using various models. The proper account of model geometries (melt−crystal and melt−ambient interfaces) is directly linked to the quality of the grown crystal, and it also serves as the most important parameter in the optimization and regulation of the CZ growth process. The so-called thermal-capillary model (TCM) developed by Derby and Brown2,3 and Derby et al.4 is realized through a finite-element method for analysis of the thermal characteristics of the CZ process. This model assumes axisymmetric heat conduction within the crystal and melt domains and includes the calculation of the moving geometric boundaries as © 2012 American Chemical Society

Received: Revised: Accepted: Published: 8675

December 18, 2011 February 17, 2012 March 31, 2012 April 1, 2012 dx.doi.org/10.1021/ie2029656 | Ind. Eng. Chem. Res. 2012, 51, 8675−8683

Industrial & Engineering Chemistry Research

Article

(see Figure 2), are defined as a measure of fluid elements being stretched in forward and backward time, respectively. The ridges of forward FTLE fields reveal repelling LCS, and the ridges of backward FTLE fields show attracting LCS. Whereas the ridges of the FTLE field have been found to be the accurate definition for LCS in many applications, Haller presented a more accurate definition for LCS as “locally the strongest repelling or attracting material surface” and developed a variational theory to identify observable LCS in an objective manner.25 If the stable manifolds of the hyperbolic point shown in Figure 2 coincide with the unstable manifold of another hyperbolic point, this heteroclinic connection is often called separatrix; it separates dynamically distinct regions in the phase space. It is well-known that, in time-independent systems, heteroclinic connections are typically broken by imposing periodic perturbations. Time-periodic systems can be viewed by considering their evolution at fixed intervals of time period (through a Poincaré section), resulting in time-independent systems. The stagnation points of the unperturbed system often remain fixed in a Poincaré section; however, the heteroclinic connection will typically break and transversely intersect. The intersection of these manifolds creates regions called lobes. Because all manifolds are invariant, meaning that fluid particles do not cross them, the fluid that is trapped in lobes is confined to remain in them as time evolves. Therefore, the motion of lobes in terms of entrainment and detrainment helps explain the transport and mixing processes.16 Heteroclinic and lobe-like structures are present for fully unsteady flows as well, but the available analytic techniques cannot be used to obtain the detailed lobe dynamic structure in general time-dependent flows. Such an analysis is important, for example, to be able to quantify transport rates, especially for engineering or biological applications. However, relying on techniques for locating hyperbolic manifolds, a computational study using LCS can reveal lobe dynamics for general aperiodic flows, without the need for a perturbative assumption or periodicity or the use of Poincaré sections.26,27 In this article, to study the mixing template of the melt flow in the CZ process, a simplified linear version of the finiteelement hydrodynamic thermal capillary model (HTCM)5 is used to determine fluid flow dynamics. The HTCM is a detailed model of governing equations in the form of Navier− Stokes, continuity, and heat equations with the Boussinesq approximation. This approach accounts for geometric features of boundaries and interface shapes, and the free surface problem is formulated to determine geometry with high accuracy. Also, the quasi-steady-state (QSS) assumption is incorporated in the HTCM because of the slow pulling rate (a few centimeters per hour) in comparison to the CZ process dynamics, namely, the evolution of the temperature and melt velocity fields, and the process is reduced to a series of steady-state solutions of governing equations. Using the velocity field data, Lagrangian coherent structures are computed in an objective manner by following the procedure given in ref 25. Also, forward and backward FTLE fields for the fluid flow system are computed to identify repelling and attracting LCS as defined in ref 23. Visually, we observe a coincidence between the LCS obtained from the two definitions. Finally, the distribution of impurities is discussed by analyzing the mixing patterns revealed by LCS.

To substantiate the influence of the melt flow on the grown crystal, Carruthers and Nassau14 and Kobayashi and Arizumi15 studied flow patterns and stagnation surfaces that prevent complete mixing of the melt in the crucible. In this work, we use another approach to study bulk mixing of the melt flow. In particular, fluid flow problems can be considered within the framework of dynamical systems.16 The behavior of a flow problem posed as a dynamical system is usually formulated in Lagrangian framework as investigation of the (particle) trajectories in phase space, whereas the fluid flow behavior at fixed points in space is described by the Eulerian approach. It is well-known that special flow patterns exist within the fluid domain that determine the governing transport structure and mechanics of the flow. These patterns are referred to as coherent structures, and when captured in terms of quantities derived from particle trajectories, they are called Lagrangian coherent structures (LCS). Recently, the identification of these structures has become an important tool for examining transport phenomena in fluid mechanics, because they describe mixing templates that govern the transport in a fluid system. (See the works of Peng and Dabiri,17 Shadden et al.,18 and Lekien and Ross19 for different applications.) To illustrate this concept, consider the phase portrait of a dynamical system near a hyperbolic point as shown in Figure 2,

Figure 2. A hyperbolic point constitutes the intersection of stable and unstable manifolds. Points x1 and x2, initially located close on either side of the stable manifold, will move apart after time τ > 0. Points x1 and x3, initially located on either side of the unstable manifold, will come closer after the same time interval.

which can be seen as a two-dimensional flow. Two particles initially located on either side of the stable manifold will diverge when advected by the flow after time interval τ > 0. This stable manifold is considered as a repelling LCS along which particles tend to move apart from each other. On the other hand, the unstable manifold can be considered as an attracting LCS along which particles located at either side tend to come close to each other. In this way, the features from dynamical systems theory, such as hyperbolic steady state or saddle point, are related to quantifiable material lines in the flow. In the contributions, Haller and Yuan20 and Haller21,22 utilized the framework of dynamical systems theory to demonstrate the existence of LCS as finite-time invariant manifolds in extended phase space for general time-dependent flows. They found LCS to be codimension one material surfaces; for two-dimensional flows, LCS are one-dimensional curves, and for three-dimensional flows, they are two-dimensional surfaces. Shadden et al.23 and Lekien et al.24 explored the definition of the LCS as “ridges of the finite-time Lyapunov exponent (FTLE) field” and provided the necessary computational procedure for the identification of FTLE fields in two- and n-dimensional systems. In their works, the forward and backward FTLEs, that is, those for τ > 0 and τ < 0 8676

dx.doi.org/10.1021/ie2029656 | Ind. Eng. Chem. Res. 2012, 51, 8675−8683


Article

where ∇ is the dimensionless differential operator; Gr and Pr are the Grashof and Prandtl numbers, respectively, given in Table 1; σ is the dimensionless stress tensor; and T and v are

2. FORMULATION The physical configuration of the mathematical model for the Czochralski crystal growth is shown schematically in Figure 3.

Table 1. Dimensionless Groups dimensionless group

symbol

definition

Biot number capillary number Grashof number Prandtl number

Bi Ca Gr Pr

heRc/ks μα/γRc gβTmpRc3/v2 ν/α

the dimensionless temperature and velocity, respectively. The stress tensor is defined as σ = −pI + σd (3) where p and σd are the dimensionless pressure and deviatoric stress, respectively, and I is the identity tensor. σd for a Newtonian melt is given by σd = Pr[∇v + (∇v)′] Figure 3. Schematic geometry of the Czochralski crystal growth process.

where the prime symbol (′) denotes the transpose. The radius of the crucible, Rc, and the melting-point temperature, Tmp, are the length and temperature scales, respectively, and the pressure and stress components are scaled by μα/Rc2. Other process parameters are introduced in Tables 2 and 3.

Different physical regions of the system are crystal and melt domains (Ds and Dm, respectively) with different thermophysical properties. To obtain melt flow dynamics in the CZ crystal growth process, the HTCM finite-element approach was developed in ref 5 considering the following basic assumptions:5 (1) The domain is axisymmetric, and the cylindrical coordinate system (r, z) is centered at the bottom of crucible, as shown in Figure 3. (2) Except for cases with the melt flow with a high rotational Reynolds number (i.e., rotation of crucible or crystal) or under the influence of external fields, such as a magnetic field, it is not necessary to consider turbulent and inertia effects.7 (3) The radius of crystal R is constant and known a priori. (4) The quasi-steady-state (QSS) approximation is used to derive the HTCM. In this approximation, the effects of batchwise decrease of the melt volume is neglected,5 and the operation of the system is simulated by a sequence of steady-state calculations at different values of the melt volume. Derby and Brown studied the validity of the QSS approximation and found it to be accurate for long crystals in systems with good radius control.28 (5) Convection heat transfer is negligible because of the low Prandtl number (on the order of 0.01) of silicon melt.4 2.1. Governing Equations. The main purpose of this work is the extraction of Lagrangian coherent structures in the CZ melt flow for which the patterns of the fluid flow are the most important issue. Here, we intend to simulate melt flow that adequately matches those given in the literature and not to include the most detailed model. To obtain a linear model, we neglect inertial effects (the nonlinear term in the momentum equation) as a first approximation of the fluid flow field, so that the momentum and continuity equations can be written in dimensionless form using the Boussinesq approximation as ∇·σ + GrPr 2(T − 1)ez = 0

(1)

∇·v = 0

(2)

(4)

Table 2. Thermophysical Properties of Silicon property

symbol

magnitude5

thermal conductivity of the melt thermal conductivity of the crystal melting-point temperature thermal diffusivity in the melt coefficient of thermal expansion of the melta surface tension emissivity viscosity of the melt kinematic viscosity of the melt contact angle of the melt−crystal interface

km [W/(m K)] ks [W/(m K)] Tmp (K) α (m2/s) β (1/K)

64 22 1683 2.6 × 10−5 10−5

γ (N/m) ε μ (N s/m2) ν (m2/s) ϕ (deg)

0.72 0.3 7 × 10−4 2.89 × 10−7 11

a

Thermal expansivity is chosen from the range given in ref 5 to regenerate results corresponding to one of the reported cases.

Table 3. Operating Parameters of the Czochralski Crystal Growth Process parameter convective heat transfer coefficient5 total heat transfer coefficient equivalent heat-transfer coefficient of radiation radius of crucible5 ambient temperature5 crucible wall temperature melt volume

symbol

magnitude

h [W/(m2 K)] he [W/(m2 K)] hr[W/(m2 K)]

10 273 263

Rc (m) Ta (K) Th Vmelt

0.073 1430

Also, neglecting the convective term due to the upward velocity of the crystal in the QSS, the energy equations for the crystal and melt are given by

∇·K i∇T = 0

(5)

where Ki, i = s, m, is the dimensionless thermal conductivity scaled by the conductivity of the crystal, ks. In the HTCM, the 8677



Article

the boundary H1(r) should be known. Instead of using eq 6, the kinematic condition

momentum and energy equations (eqs 1 and 5) are coupled because of the buoyancy forces appearing from Boussinesq approximation. 2.2. Boundary Conditions. The boundary conditions of eqs 1 and 2 are associated with the specification of the velocity or stress on the melt boundary, ∂Dm (see Figure 3). On the crucible walls, ∂D3 and ∂D4, and the melt−crystal interface, ∂D0, no-slip condition is applied to the velocity components. The appropriate boundary condition for the melt−ambient interface, ∂D1, involves the balance of forces at the surface, as well as a kinematic condition stating that there is no velocity component normal to the boundary. Neglecting the forces from the ambient and the Marangoni effect (i.e., flow driven by a surface tension gradient), the balance of forces in terms of stress components yields29

1 2/ Ca

σnn =

is imposed as the boundary condition for the momentum equation;5,29 then, eqs 6−8 are solved numerically with the boundary conditions dH1 = −cot(ϕ), dr

dH1 = 0, dr

(7)

2/ =

dr 2

⎡ ⎢⎣1 +

3/2 dH1 2 ⎤ dr

( ) ⎥⎦

1 + r ⎡ ⎢⎣1 +

1/2 dH1 2 ⎤ dr

(8)

σSBε(T̅ 4 − Ta 4) T̅ − Ta

x(t0) = x 0

(14)

(15)

The flow map encodes the Lagrangian paths of the particles used to capture the local stretching in the flow. Consider a point located at x0 ∈ D at time t0. When advected by the flow to Ftt00+τ after a time interval τ, the symmetric Cauchy−Green finite-time deformation tensor Ctt00+τ is given by

(∇Ftt00+ τ )′∇Ftt00+ τ

(9)

(16)

with real positive eigenvalues sorted as λ1(x0,t0,τ) ≤ ··· ≤ λn(x0,t0,τ) and associated eigenvectors ξ1(x0,t0,τ), ..., ξn(x0,t0,τ). For linearized flow, if an infinitesimal perturbation to the trajectory is given by δx0, maximum stretching of this perturbation after the time interval τ occurs when δx0 is aligned with ξn(x0,t0,τ) and

where Bi is the Biot number given in Table 1, n is the outward unit vector normal to the boundary, and Ta is the ambient temperature. To model radiative heat transfer, an equivalent linearized heat transfer coefficient, hr, is associated with radiation as hr =

(13)

Ftt0 : D → D: x 0 ↦ Ftt0(x 0) = x(t ; x 0, t0)

The boundary conditions of eqs 5 are associated with the specification of temperature or heat flux on the melt and crystal boundaries, ∂Dm and ∂Ds, respectively. On the crucible side wall, ∂D4, temperature Th is specified. The crucible bottom, ∂D3, is assumed to be isolated. Thermal boundary conditions are applied to the exposed surfaces of the system, ∂D1 and ∂D2 ⎛ T ⎞ K i n ·∇T = Bi⎜⎜T − a ⎟⎟ Tmp ⎠ ⎝

(12)

with the dot denoting differentiation with respect to time t, x ∈ D ⊂ n , and the smooth velocity field v(x,t) defined over fluid domain D over the time interval t = [t0, t0 + τ]. The solution of system 14 is the trajectory x(t;x0,t0) of a particle at time t that was initially (i.e., at t = t0) located at x0. Fixing the initial time t0 and the final time t, the flow map Ftt0 is defined as the map that takes the point from its initial position x0 in the fluid domain at time t0 to its final position x in the fluid domain at time t

dH1 dr

( ) ⎥⎦

r = Rc

at

ẋ = v(x , t ),

where σnn and σnt are normal and tangent stresses, respectively, at any point on the interface; Ca is the capillary number given in Table 1; and 2/ is the dimensionless axisymmetric mean curvature, given by d2H1

at r = R

where ϕ is the contact angle of the melt−crystal interface (see Figure 3). 2.4. Lagrangian Coherent Structures. In this section, we present the formulation of the LCS identification for a general time-dependent fluid flow problem. The fluid flow evolution can be considered as a dynamical system of the form

(6)

σnt = 0

(11)

v·n = 0

max ||δx(x 0, t0 , τ )|| = δx 0

(10)

λn(x 0, t0 , τ ) ||δx 0||

(17)

where δx 0 is aligned with ξn(x0,t0,τ). Equation 17 can be written as

where σSB is the Stefan−Boltzmann constant and T̅ is the mean temperature along the boundary, which is assumed to be the average of ambient (Ta) and heating (Th) temperatures. The Biot number is based on the total radiative and convective heat transfer, he = h + hr. Finally, on ∂D5, symmetry boundary conditions are applied for both the momentum and energy equations. 2.3. Shapes of Interfaces. The shapes of the melt−crystal and melt−ambient interfaces are determined by the melt dynamics, heat transfer, and capillarity. The solidification interface, H0(r), is determined along the melting-point isotherm in the whole domain, so that T = 1. Equations 6 and 8 characterize the shape of the melt−ambient interface. To apply eq 6 as a boundary condition to the momentum equation, the shape of

max ||δx(x 0, t0 , τ )|| = exp⎡⎣Λ tt00+ τ(x 0)|τ |⎤⎦||ξn(x 0, t0 , τ )|| δx 0

(18)

where Λ tt00+ τ(x 0) =

1 ln |τ |

λ n (x 0 , t 0 , τ )

(19)

Λtt00+τ(x0) represents the (largest) finite-time Lyapunov exponent (FTLE) over integration time τ. In general, FTLE measures how much particles initially located at x0 separate after time interval τ. 8678



Article

Figure 4. Finite-element mesh for four time instances during the process: (a) H2 = 0.1 m, Vmelt = 7 × 10−4 m3, and DOF = 13621; (b) H2 = 0.135 m, Vmelt = 5.241 × 10−4 m3, and DOF = 11221; (c) H2 = 0.17 m, Vmelt = 3.501 × 10−4 m3, and DOF = 10119; and (d) H2 = 0.19 m, Vmelt = 2.494 × 10−4 m3, and DOF = 9017. (DOF is the number of degrees of freedom of the finite-element analysis including unknown pressures, velocities, and temperatures at corresponding nodes.)

The regions of maximum separation or attraction of fluid elements are known as local maximizing curves of the FTLE field.20 Shadden et al. defined Lagrangian coherent structures as ridges of the FTLE field and proposed two mathematical definitions for the ridge.23 The ridges of a smooth function is a set of curves whose points are local maxima of the function in at least one dimension. They estimated that the flux across an LCS is small and negligible for a well-defined LCS. In other words, these ridges can be seen as transport barriers that can be detected in the fluid flow by the algorithm proposed by Senatore and Ross.30 In this definition, the ridges of forward (τ > 0) and backward (τ < 0) FTLE fields reveal repelling and attracting LCS, respectively. However, it is shown that attracting LCS can be extracted from local minima of the forward FTLE field,31 which allows both repelling and attracting LCS to be obtained from a single numerical run. Although excellent agreement of FTLE-based LCS with stable and unstable manifolds has been found in physical flows, Haller provided counterexamples for LCS defined as FTLE ridges.25 Then, through a variational theory, he specified computable sufficient and necessary criteria for the extraction of LCS in terms of invariants of the Cauchy−Green strain tensor field. Applying these criteria to two-dimensional flows, an LCS is a one-dimensional material surface that is approximately normal to the eigenvector field associated with the largest eigenvalues of the Cauchy−Green strain tensor field. With these conditions, the (hyperbolic) LCS is the strongest repelling or attracting surface among all nearby surfaces. Repelling LCS at time t0 are obtained for τ > 0, and attracting LCS at t0 are obtained for τ < 0 (see eq 16).

conditions. For the melt domain Dm, velocity components are approximated by nine-node quadratic Lagrangian interpolation functions, and pressure is approximated by a four-node bilinear approximation that is continuous across element boundaries. Temperature is interpolated by nine-node quadratic Lagrangian interpolation functions for both domains. To solve the precedent equations using the finite-element method, an approximate position is assumed for the unknown boundaries. Solving the finite-element model, the boundaries are modified as stated in section 2.3, and this procedure continues until the shapes of the interfaces converge with the desired accuracy. As two conditions, the volume of the melt, Vmelt, and the crystal radius are kept constant with predefined values. When the crystal grows, because of the decrease in the surface area of heat transfer from the crucible wall to the melt as well as the increase in the surface area of heat transfer from the crystal to the ambient, the heating temperature, Th, increases to maintain a constant crystal radius. The radius of the crystal is defined by the radial coordinate of the triple point (crystal−melt−ambient meeting point) at which the melt− crystal [H0(r)] and melt−ambient [H1(r)] interfaces meet. In the QSS approximation, the process evolution is considered as a series of time steps (snapshots) for which the steady governing equations are solved. From the current time step to the next one, only the geometry of the problem changes, preserving the total crystal and melt mass, and the steady finiteelement problem is solved. The computational algorithm is given as follows: For each time step calculate crystal height and melt volume initialize boundaries modify heating temperature for crystal radius to be constant Do while boundaries are not converged finite-element analysis modify boundaries preserving melt volume End Do End For

3. NUMERICAL MODEL The governing equations are solved using the finite-element method to obtain velocity data for the CZ melt flow dynamics. Subsequently, the Lagrangian coherent structures are identified for another computational grid associated with the computation of the deformation tensor field. 3.1. Finite-Element Model. A finite-element model is developed to solve eqs 1, 2, and 5 with appropriate boundary 8679



Article

Figure 5. Evolution of streamlines and velocity vectors of the melt flow during the process in which one convective cell transforms into a two-cell configuration.

Figure 6. Temperature contours (K) during the CZ crystal growth process.

Figure 7. Evolution of the (top) normalized forward FTLE field and (bottom) repelling LCS of the melt flow for τ = 5, 10, 25, and 60 s from left to right, respectively.

3.2. Computation of Lagrangian Coherent Structures. To compute the flow map and deformation tensor, a discrete uniform grid of computational points over the melt domain is defined. Usually, this grid should be fine enough to account for the high FTLE grow rate near the LCS.23 For each point in the grid, the flow map ϕtt00+τ should be obtained, that is, the position of each point when advected by the flow after the time interval τ. To this end, the trajectory of each point is determined from the provided velocity data evolution calculated by the HTCM finite-element analysis using a fourth-order Runge−Kutta integration algorithm. Note that common finite-element interpolation functions are used for spatial interpolation of discrete velocity data. Finally, the gradient of the flow map is computed by finite

differencing, for example, for point x0i,j located at [ri,j(t0), zi,j(t0)], this gradient is given by ∇Ftt00+ τ (x 0i , j) ⎡ ri + 1, j(t0 + τ ) ⎢ ri + 1, j(t0) ⎢ =⎢ ⎢ zi + 1, j(t0 + τ ) ⎢ ri + 1, j(t0) ⎢⎣

ri , j + 1(t0 + τ ) − ri , j − 1(t0 + τ ) ⎤ ⎥ − ri − 1, j(t0) zi , j + 1(t0) − zi , j − 1(t0) ⎥ ⎥ − zi − 1, j(t0 + τ ) zi , j + 1(t0 + τ ) − zi , j − 1(t0 + τ ) ⎥ ⎥ − ri − 1, j(t0) zi , j + 1(t0) − zi , j − 1(t0) ⎥⎦ − ri − 1, j(t0 + τ )

(20)

Once the two-dimensional flow map is obtained, computation of the deformation tensor and corresponding eigenvalues and eigenvectors is straightforward. For all computations, finite 8680



Article

Figure 8. Evolution of the (top) normalized backward FTLE field and (bottom) attracting LCS of the melt flow for τ = 5, 10, 25, and 60 s from left to right, respectively.

Figure 9. Intersections of attracting and repelling LCS define strongly hyperbolic points in the flow. Different regions of the melt flow are attracted to and repelled from the attracting and repelling LCS, respectively. The time step from one panel to the next is 6 s.

the figure, one convective cell in the melt changes into two cells as the crystal grows and the volume of the melt decreases. Temperature contours are shown in Figure 6. In what follows, Lagrangian coherent structures are identified for the melt flow using the procedure given in ref 25, accompanied by the computation of forward and backward FTLE fields. Note that the time interval τ for computation of the finite time deformation tensor is chosen long enough so that important dynamics have sufficient time to induce separation. Increasing the integration time tends to sharpen the ridges of the FTLE field and extend their length in the flow domain. To represent the observable evolution of LCS for representative time instances, the FTLE fields are scaled to have values in the interval [0, 1]. Figure 7 shows the evolution of the forward FTLE field and the corresponding repelling LCS during crystal growth. There is one notable LCS located below the crystal−melt interface that surrounds the triple point and is associated with the melt undergoing phase transformation. This LCS does not move, and its shape does not change significantly. Moreover, because

differencing is used to perform differentiation, and all variables are linearly interpolated across computational points.

4. RESULTS The thermophysical properties and operating parameters used in the modeling of silicon crystal growth are listed in Tables 2 and 3. For the conditions presented here, the radius of the crystal is assumed to be constant (R = 0.035 m), for which the heating temperature Th should be changed during the process. The finite-element mesh and its specifications for four time instances during the CZ crystal growth process are shown in Figure 4, and each mesh is refined until no significant change in velocity and temperature fields is manifested. These instances were selected to represent the time evolution of the process, and for all of them, the boundaries at the solid−melt and melt− ambient interfaces converged after five iterations, as shown in Figure 4. Figure 5 shows the evolution of the melt flow pattern during the process, for which streamlines were obtained using the Stokes stream function of the axisymmetric flow. As depicted in 8681



■

a closed curve (pocket) is made with this LCS and domain boundaries, it is viewed as the transport barrier in the flow. This fact reveals that impurities cannot flow out of this pocket and that they will be incorporated into the grown crystal. Another LCS is formed in the middle of the melt domain as an open ring at the beginning of the process. As the volume of the melt decreases and two convective cells appear, this LCS is expanded and prevents the silicon melt from flowing freely in the crucible during the process, which leads to different mixing properties and movement of impurities. Therefore, the physical properties of the grown crystal when there is one convective cell in the melt is different from those when there are two convective cells. Figure 8 shows the evolution of the backward FTLE field and the corresponding attracting LCS in the melt during the CZ crystal growth process. Although the attracting LCS below the crystal is not closed, its pattern prevents the free flow of the melt. Also, the same discussion can be made for the LCS formed in the middle of the crucible, when there are two convective cells. A visual comparison of the FTLE fields and corresponding LCS in Figures 7 and 8 shows that the local maxima of the FTLE fields can identify the same LCS as obtained from the objective procedure. All of these repelling and attracting LCS satisfy conditions given in ref 25. One property of the LCS is shown in Figure 9: the intersections of attracting and repelling LCS define strongly hyperbolic points in the flow. In this figure, both repelling and attracting LCS for the last representative case (right panels of Figures 7 and 8) located in the middle of the crucible are overlaid (compare with Figure 2), and a sequence of snapshots of particle trajectories in forward time is depicted. Note how different regions of the melt flow with different colors are attracted and repelled from the identified LCS, whereas no region flows across the LCS to be mixed in other regions.

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The computational part of this work used the numerical servers as infrastructure and resources of AICT (Academic Information and Communication Technologies) of the University of Alberta. Also, the authors thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the article. Financial support by Natural Science and Engineering Research Council of Canada (NSERC) Discovery Grant 386508-2011 is gratefully acknowledged.

■

REFERENCES

(1) Yu, P. Y.; Cardona, M. Fundamentals of Semiconductors: Physics and Materials Properties; Springer Verlag: Berlin, 2010. (2) Derby, J. J.; Brown, R. A. Thermal-capillary analysis of Czochralski and liquid encapsulated Czochralski crystal growth: I.I: Simulation. J. Cryst. Growth 1986, 74, 605−624. (3) Derby, J. J.; Brown, R. A. Thermal-capillary analysis of Czochralski and liquid encapsulated Czochralski crystal growth: II. Processing strategies. J. Cryst. Growth 1986, 75, 227−240. (4) Derby, J. J.; Atherton, L. J.; Thomas, P. D.; Brown, R. A. Finiteelement methods for analysis of the dynamics and control of Czochralski crystal growth. J. Sci. Comput. 1987, 2, 297−343. (5) Sackinger, P. A.; Brown, R. A.; Derby, J. J. A finite element method for analysis of fluid flow, heat transfer and free interfaces in Czochralski crystal growth. Int. J. Numer. Methods Fluids 1989, 9, 453− 492. (6) Prasad, V.; Zhang, H.; Anselmo, A. P. In Transport Phenomena in Czochralski Crystal Growth Processes; Hartnett, J. P., Irvine, T. F., Cho, Y. I., Greene, G. A., Eds.; Advances in Heat Transfer; Elsevier: New York, 1997; Vol. 30; pp 313−435. (7) Liu, L.; Kakimoto, K. Partly three-dimensional global modeling of a silicon Czochralski furnace: I. Principles, formulation and implementation of the model. Int. J. Heat Mass Transfer 2005, 48, 4481−4491. (8) Wagner, C.; Friedrich, R. Direct numerical simulation of momentum and heat transport in idealized Czochralski crystal growth configurations. Int. J. Heat Fluid Flow 2004, 25, 431−443. (9) Kalaev, V. V.; Lukanin, D. P.; Zabelin, V. A.; Makarov, Y. N.; Virbulis, J.; Dornberger, E.; von Ammon, W. Prediction of bulk defects in CZ-Si crystals using 3D unsteady calculations of melt convection. Mater. Sci. Semicond. Process. 2002, 5, 369−373. (10) Lipchin, A.; Brown, R. A. Hybrid finite-volume/finite-element simulation of heat transfer and melt turbulence in Czochralski crystal growth of silicon. J. Cryst. Growth 2000, 216, 192−203. (11) Abbasoglu, S.; Sezai, I. Three-dimensional modelling of melt flow and segregation during Czochralski growth of GexSi1− single crystals. Int. J. Thermal Sci. 2007, 46, 561−572. (12) Fedoseyev, A. I.; Alexander, J. I. D. Investigation of vibrational control of convective flows in Bridgman melt growth configurations. J. Cryst. Growth 2000, 211, 34−42. (13) Pimputkar, S. M.; Ostrach, S. Convective effects in crystals grown from melt. J. Cryst. Growth 1981, 55, 614−646. (14) Carruthers, J.; Nassau, K. Nonmixing cells due to crucible rotation during Czochralski crystal growth. J. Appl. Phys. 1968, 39, 5205−5214. (15) Kobayashi, N.; Arizumi, T. Computational analysis of the flow in a crucible. J. Cryst. Growth 1975, 30, 177−184. (16) Wiggins, S. Chaotic Transport in Dynamical Systems; Springer Verlag: Berlin, 1992.

5. CONCLUSIONS Lagrangian coherent structures (LCS) have been identified for the melt flow of the Czochralski crystal growth process to study the bulk mixing in the crucible. A linear finite-element method based on the hydrodynamic thermal capillary model (HTCM) is developed to solve for the governing equations and to find free boundary shapes. The velocity fields are used to compute the finite-time deformation tensor using particle trajectories to identify Lagrangian coherent structures during the process. The LCS are identified by two relevant methodologies, namely, FTLE ridges23 and the objective approach,25 and their correspondence is provided. The presence of LCS during the process demonstrates the existence of barriers that prevent complete mixing in the silicon melt, as they are identified in the regions close to the melt− solid interface. Furthermore, the incomplete mixing of the melt induces a nonuniform distribution of impurities in the grown crystal and, consequently, affects the physical properties of the grown crystal that are of significant importance in the semiconductor industry. We have provided the foundation for the Lagrangian framework of dealing with incomplete mixing, and we plan to investigate control methodologies to address mixing in the CZ crystal growth process. 8682



Article

(17) Peng, J.; Dabiri, J. O. Transport of inertial particles by Lagrangian coherent structures: application to predator-prey interaction in jellyfish feeding. J. Fluid Mech. 2009, 623, 75−84. (18) Shadden, S. C.; Astorino, M.; Gerbeau, J. F. Computational analysis of an aortic valve jet with Lagrangian coherent structures. Chaos 2010, 20, 017512. (19) Lekien, F.; Ross, S. D. The computation of finite-time Lyapunov exponents on unstructured meshes and for non-Euclidean manifolds. Chaos 2010, 20, 017505. (20) Haller, G.; Yuan, G. Lagrangian coherent structures and mixing in two-dimensional turbulence. Physica D: Nonlinear Phenom. 2000, 147, 352−370. (21) Haller, G. Lagrangian coherent structures from approximate velocity data. Phys. Fluids 2002, 14, 1851−1861. (22) Haller, G. Distinguished material surfaces and coherent structures in three-dimensional fluid flows. Physica D: Nonlinear Phenom. 2001, 149, 248−277. (23) Shadden, S. C.; Lekien, F.; Marsden, J. E. Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows. Physica D: Nonlinear Phenom. 2005, 212, 271−304. (24) Lekien, F.; Shadden, S. C.; Marsden, J. E. Lagrangian coherent structures in n-dimensional systems. J. Math. Phys. 2007, 48, 065404. (25) Haller, G. A Variational Theory of Hyperbolic Lagrangian Coherent Structures. Physica D: Nonlinear Phenom. 2011, 240, 574− 598. (26) Shadden, S. C.; Dabiri, J. O.; Marsden, J. E. Lagrangian analysis of fluid transport in empirical vortex ring flows. Phys. Fluids 2006, 18, 047105. (27) Du Toit, P.; Mezić, I.; Marsden, J. E. Coupled oscillator models with no scale separation. Physica D: Nonlinear Phenom. 2009, 238, 490−501. (28) Derby, J. J.; Brown, R. A. On the quasi-steady-state assumption in modeling Czochralski crystal growth. J. Cryst. Growth 1988, 87, 251−260. (29) Reddy, J. N.; Gartling, D. K. The Finite Element Method in Heat Transfer and Fluid Dynamics; CRC Press: Boca Raton, FL, 2010. (30) Senatore, C.; Ross, S. D. Detection and characterization of transport barriers in complex flows via ridge extraction of the finite time Lyapunov exponent field. Int. J. Numer. Methods Eng. 2011, 86, 1163−1174. (31) Haller, G.; Sapsis, T. Lagrangian coherent structures and the smallest finite-time Lyapunov exponent. Chaos 2011, 21, 023115.

8683


Analysis of Melt Flow Mixing in Czochralski Crystal Growth Process

Recommend Documents