Dynamical Analysis of Melt Flow in the Bridgman Process - Industrial

Oct 17, 2014 - The Bridgman crystal growth process represents a prime example of a relevant industrial ... Industrial & Engineering Chemistry Research...
1 downloads 0 Views 3MB Size
Article pubs.acs.org/IECR

Dynamical Analysis of Melt Flow in the Bridgman Process Mojtaba Izadi,† Youssef Belhamadia,‡ and Stevan Dubljevic*,† †

Department of Chemical and Materials Engineering and ‡Department of Biomedical Engineering, University of Alberta, Edmonton, Alberta T6G 2R3, Canada ABSTRACT: The Bridgman crystal growth process represents a prime example of a relevant industrial process in which fluid flow and mixing features impact desired process specifications. In particular, we analyze the bulk mixing properties of the melt flow in a vertical Bridgman process through the identification of Lagrangian coherent structures (LCS). The Navier−Stokes, continuum, and heat transport equations are formulated in a mixed finite-element universal integration scheme to simulate the multiphase problem using the finite-element method. The highly resolved time-varying accurate numerical solution of the velocity field is utilized to identify the transport barriers in terms of finite time Lyapunov exponent (FTLE) ridges whose movements govern the mixing in the melt flow. The ridge structures are dynamically driven materials surfaces which provide the insight into the efficiency of mixing in two distinct operating regimes of the Bridgman crystal growth process.

1. INTRODUCTION The dynamical analysis applied to fluid flow systems has attracted the interest of fluid dynamicists and engineers for its relevance in understanding fundamental fluid flow features and for providing practitioners with a new insight into a large number of industrial processes (chemical, pharmaceutical, petroleum, food). In particular, in the last two decades, methods and quantities derived from dynamical system theory emerged as Lagrangian or kinematic theory of transport and mixing, and spread out as methodology to address features of flow properties, laboratory equipment, static and stirred mixers, and industrial-scale processes.1 One of the important crystal manufacturing processes for the production and/or purification of grown crystal is the Bridgman crystal growth process. The Bridgman process is one of the common methods used for the production of melt-grown crystals,2 where the material is introduced in a crucible placed in a furnace. In the classical vertical Bridgman method (see Figure 1), the upper zone of the furnace has temperatures above the melting point of the material and the lower zone with a temperature below melting point separated by an adiabatic zone as a baffle between the

two. The crucible is raised into the upper zone and solid material is melted, possibly with a remaining nonmelted singlecrystal seed for the initiation of crystallization at the crucible’s bottom in the lower zone. After the temperature stabilizes, the crucible is lowered slowly into the lower zone, so that the crystal grows from the melt. In contrast, in an inverted Bridgman process, the crucible is heated from the lower zone and cooled from the upper zone of the furnace (see Figure 1). In the melt-grown crystallization processes, thermal convection of the melt flow plays an important role by affecting the mixing properties, as well as the heat- and mass-transfer phenomena.3 The consequence of incomplete mixing in the melt is that a gradient of impurities concentration exists in the melt near the grown-crystal/solid-melt interface.4 This affects the crystal characteristics, such as compositional uniformity of the grown crystal. In order to substantiate the influence of the melt flow on the grown crystal, Carruthers and Nassau,5 and also Kobayashi and Arizumi,6 studied flow patterns and stagnation surfaces that prevent complete mixing of the melt in the crucible. The Lagrangian coherent structures (LCS) approach, which is conceptually formulated from the overlapping areas of nonlinear dynamics and fluid dynamics, is utilized in an analysis of the bulk mixing of the melt flow in the Bridgman crystal growth process. In other words, a novel way of understanding transport in complex fluid flows is given by identification of the LCS structures that define material transport across the boundaries of the LCS structures over a selected period of time. These structures are distinguished material lines, i.e., continuous smooth curves of fluid elements advected by the flow in a two-dimensional (2-D) space, that play the dominant role in attracting and repelling neighboring fluid elements. Transport by material advection suggests the kinematic or Lagrangian perspective of the fluid flow description where Received: July 30, 2014 Revised: October 15, 2014 Accepted: October 17, 2014

Figure 1. Schematic configuration of the classical (left) and inverted (right) vertical Bridgman process. Th and Tc are heating and cooling temperatures, respectively. © XXXX American Chemical Society

A

dx.doi.org/10.1021/ie503030z | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Lagrangian, i.e., the flux of the material across them is small. In addition to finding LCSs from FTLE fields, Haller presented a variational theory to identify LCSs that result in Lagrangian material lines.12 Recently, the identification of these structures has become an important tool for examining the transport phenomena in fluid mechanics, since they describe mixing templates that govern the transport in a fluid system (see the works of Dabiri,13 Shadden et al.,14 Lekien and Ross,15 and Izadi and Dubljevic16 for different applications). Although the identification of the LCS structures defines objective indexes (FTLE fields) that quantify the quality of the mixing process within each of the invariant regions associated with the mixing flow, it does not provide any information about how these region-dependent indexes can be merged together to construct some global measure of mixing in the entire flow domain. In other words, the pure kinematic approach to the analysis of the mixing templates in the Bridgman crystal growth process lacks the presence of a diffusive transport mechanism, and this insight should be included in future research. In this paper, we address the analysis of mixing features in the melt flow in the vertical Bridgman crystal growth process through the kinematic approach of identification of Lagrangian coherent structures. In order to study the mixing templates of the melt flow in the Bridgman process, a velocity field is obtained based on the mixed finite-element formulation of the set of nonlinear governing equations, which include the Navier−Stokes, continuity, and energy equations. Using the obtained velocity field data, forward and backward FTLE fields for the fluid flow system are computed to identify repelling and attracting LCSs in the Bridgman crystal growth process. The simulation and associated data reveal the relation between the impurities distributions in two important operational settings of the Bridgman crystal growth process. In particular, the two process regimesone characterized with a low Rayleigh number with two formed stable vortexes and less efficient mixing, and another one characterized with a high Rayleigh number and aperiodic fluid flow and more efficient mixingare presented from theoretical and practical perspectives.

individual fluid elements are tracked as they are advected by the flow through the flow domain. By contrast, the Eulerian point of view considers the properties of a flow field at each fixed point in space and time. A detailed understanding of Lagrangian transport already exists for time-independent flows such as the steady saddlepoint flow in Figure 2. Fluid parcels close to the saddle point

Figure 2. A (hyperbolic) saddle point in a time-independent twodimensional (2-D) flow, which constitutes the intersection of stable (red) and unstable (blue) manifolds. A fluid element close to the saddle point on the two sides of the stable manifold (repelling LCS) is drawn away along the unstable manifold (attracting LCS).

are drawn away from the stable manifold, that is the line of flowing material serving as a repulsive transport barrier, toward a transverse unstable manifold, that is the material line constituting an attractive transport barrier. The unstable manifold acts as a core organizing structure in the vicinity of the saddle point, attracting all nearby fluid particles, which then stretch out to adopt its shape. Also, one manifold (for τ > 0) looks like the other with the time reversed (i.e., τ < 0). In addition to steady flows, prominent material lines are also known to exist in periodic and quasi-periodic flows, where they serve as skeletons of observed tracer patterns. However, the identification of dynamical skeletons for material patterns in flows with complex spatial and temporal structure is an ongoing challenge. Assuming that the velocity field is observed for times ranging over a finite interval, the LCSs of the flow during that interval are the material lines that repel or attract nearby fluid trajectories at the highest local rate relative to other material lines nearby. Overall, the repelling and attracting LCSs play similar roles to the stable and unstable manifolds, respectively, of the saddle point in Figure 2, where the repelling LCSs direct particles to different domains of the attracting LCSs. The framework of the dynamical systems theory to demonstrate the existence of LCS as finite-time invariant manifolds in the extended phase space for general timedependent flows is studied in the work of Haller and Yuan7 and Haller.8 A pioneering insight into Lagrangian features of the velocity field data in these contributions is the landscape of the finite time Lyapunov exponent (FTLE) field. The Lyapunov exponent is a scalar measure of how much initially adjacent particles from a given location have separated. Regions of high separation are locally of the strongest diverging flow. Performing the same procedure in backward time, one identifies regions with high backward-time FTLE values with strongest local convergence in forward time. The complex patterns of FTLE distributions for physical flow fields were connected to LCSs in ref 9, where the ridges of maxima in the FTLE field are defined as repelling LCSs in forward time and attracting LCSs in backward time. This definition was later explored by Shadden et al.10 and Lekien et al.11 in 2-D and ndimensional systems, showing that the ridges are almost

2. FORMULATION Kaenton et al. study the transition flow regimes for the solidification of gallium in an inverted vertical Bridgman configuration.17 Depending on the Rayleigh number (Ra), steady symmetric flow with two counter rotating cells, steady asymmetric flow with one dominant cell, symmetric periodic motion with one or two frequencies, and asymmetric quasiperiodic and nonperiodic regimes are reported for the 2-D Cartesian melt flow simulation. 2.1. Governing Equations. Following the same geometry as that given in ref 17 and shown schematically in Figure 3, a velocity field in the melt is obtained in this work for pure gallium whose thermophysical properties in solid and liquid states are the same. The fluid is assumed to be Newtonian incompressible, and the Boussinesq approximation is considered. The governing equationsincluding the Navier−Stokes, continuity, and energy equationscan be written in a dimensionless form as

B

g ∂v + v∇·v = Pr ∇·(∇v) − ∇p + Av + RaPr T ∂t g

(1)

∇·v = 0

(2) dx.doi.org/10.1021/ie503030z | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

⎧ 0 if T < Tf ϕ = F (T ) = ⎨ ⎩1 if T > Tf ⎪



(4)

In applications, phase change is not always instantaneous and may occur in a small temperature range [Tf − ϵ, Tf + ϵ]. The relationship for F(T) can thus be replaced by a regularized function Fϵ(T). 2.2. Finite Element Model. Solving the system defined by eqs 1−4 is a challenging task. In this paper, we adopted the numerical method used in Belhamadia et al.,18 which consists of using a mixed formulation for Navier−Stokes (see eqs 1 and 2). The finite-element discretization is based on quadratic polynomials (written by using a hierarchical basis) for the velocity and linear (continuous) polynomials for the pressure. Using the test functions (u,q), the variational formulation can be written as ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩

Figure 3. Geometry of the inverted Bridgman process: W and L are the width and height of the crucible and LΔT is the length of the isolated baffle.

1 ∂ϕ ∂T − ∇·(∇T ) = 0 + v ·∇T + ∂t Ste ∂t

(3)

where v is the velocity, p the pressure, and T the temperature. The Boussinesq approximation, in common use, is used to describe the melt’s convective motion. In the classical Boussinesq model, the buoyancy term varies linearly with the temperature and is given by RaPr(g/∥g∥)T. The dimensionless formulation is obtained by using W as the scale factor for the length, W2/α for the time, α/W for the velocity, and ρα2/W2 for the pressure. Pr, Ra, and Ste are dimensionless parameters. Pr is the Prandtl number, which is defined as Pr =

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩

v α

gβ(Th − Tc)W αv



∫Ω q∇·v dΩ = 0









∫Ω ⎝ ∂∂Tt uT + (v·∇T )uT ⎠ dΩ + ∫Ω ⎝ Ste1 ∂∂ϕt uT ⎠ dΩ + ∫ (∇T ·∇uT ) dΩ = 0 Ω ⎜







∫Ω (ϕ − Fϵ(T ))uϕ dΩ = 0

In all our numerical simulations, a fully implicit backward second-order scheme is used for time discretization. For instance, given approximate solutions Tn−1, Tn, and Tn+1 at times tn−1, tn, and tn+1, respectively, the time derivative at time tn+1 is approximated by

3

Ste is the Stefan number, which is defined as Ste =





A mixed formulation is also employed for the energy (see eqs 3 and 4). The finite-element discretization is based on quadratic polynomials for both temperature T and phase change variable ϕ. Using the test functions (uT,uϕ), the corresponding variational formulation can be written as

Ra is the Rayleigh number, which is defined as Ra =



∫Ω ⎝ ∂∂vt u + (v∇·v)u⎠ dΩ + ∫ (Pr ∇v∇u − p∇u − A(φ)vu) dΩ Ω g Tu dΩ = ∫ RaPr Ω g

∂T n + 1 3T n + 1 − 4T n + T n − 1 (t ) ≃ ∂t 2Δt

cp(Th − Tc)W 3 L

Therefore, the overall algorithm for solving the proposed model is described as follows. (1) Starting from (Tn,ϕn), the solutions (vn+1,pn+1) are obtained by solving the system

In these dimensionless parameters, β is the thermal expansion coefficient, α the thermal diffusivity, ν the kinematic viscosity, and L the latent heat. The specific heat (cp) and the density (ρ) are assumed to be equal in the liquid and solid phases. A is defined as

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

⎡ (1 − ϕ)2 ⎤ A = − a⎢ 3 ⎥ ⎣ ϕ +b ⎦

where b = 10−6 is a constant introduced to avoid division by zero and a = 1015. In the fluid domain (Ωl), A is zero and has no influence. However, in the totally solid region (Ωs), A has a large value and forces the velocity to be zero. The function ϕ is a phase change variable that satisfies a simple algebraic equation of the form C



∫Ω ⎜⎝ 3v

n+1

− 4v n + v n − 1 ⎞ u ⎟ dΩ 2Δt ⎠

∫Ω ((vn+1·∇vn+1)u + Pr∇vn+1∇u) dΩ + ∫ ( −p(n + 1) ∇v − A(φn)v(n + 1)u) dΩ Ω g (n) = ∫ RaPr T u dΩ Ω g +

∫Ω q∇·vn+1 dΩ = 0 dx.doi.org/10.1021/ie503030z | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

ξ1(x0,t0,τ), ξ2(x0,t0,τ). The deformation tensor is positive definite; therefore, it has positive eigenvalues. For the linearized flow, if an infinitesimal perturbation to the trajectory is defined by δx0, after the time interval τ, maximum stretching of this perturbation occurs when δx0 is aligned with ξ2(x0, t0, τ) and

and the Newton’s method is used to solve for the zeros of a large nonlinear system of algebraic equations. (2) Starting from (vn+1,pn+1), the solutions (Tn+1,ϕn+1) are obtained by solving the system ⎧ ⎛ (n + 1) − 4T (n) + T (n − 1) ⎞ ⎪ ⎜⎜ 3T uT ⎟⎟ ⎪ Ω⎝ 2Δt ⎠ ⎪ (n + 1) (n + 1) (v )u T dΩ ·∇T ⎪ + Ω ⎪ ⎪ ⎛ 1 3φ(n + 1) − 4ϕ(n) + ϕ(n − 1) ⎞ ⎨ uT ⎟⎟ dΩ + ⎪ Ω ⎜⎜ Ste 2Δt ⎝ ⎠ ⎪ ⎪ + (∇T (n + 1)·∇uT ) dΩ = 0 Ω ⎪ ⎪ (n + 1) ⎪ (ϕ − Fϵ(T (n + 1)))uϕ dΩ = 0 ⎩ Ω



max δx(x0 , t0 , τ ) = δx0





t0+ τ

max δx(x0 , t0 , τ ) = e Λt0

Λ tt00+ τ(x0) =

0

(5)

Φ(t ) =

(6)

The flow map, in essence, encodes the Lagrangian or kinematic paths of the particles used to capture the local stretching in the flow by which the particles are advected. By considering at time t0 a point located at x0 ∈ Ωl, which is advected by the flow to ϕtt00+τ after a time interval τ, the symmetric finite-time Green deformation tensor is obtained and defined by the following: Ctt00+ τ = (∇ϕtt0+ τ )′∇ϕtt0+ τ 0

0

(9)

1 ln |τ |

λ 2 (x 0 , t 0 , τ )

(10)

where the Λtt00+τ(x0) represents the (largest) finite-time Lyapunov exponent (FTLE) over integration time τ. Generally, from the kinematic point of view, the FTLE scalar field describes how many particles initially located at x0 separate after time interval τ by being advected by flow. The Lagrangian coherent structures (LCSs) are defined as ridges of the scalar FTLE field.10 This implies that the intuitive notion of the ridge suffices to identify LCS; however, Haller provided counterexamples for the LCSs defined as FTLE ridges.12 In refs 12 and 19, variational theory is presented based on the definition of the LCS as “locally the strongest repelling or attracting material surface.” In this theory, the mathematically given sufficient and necessary conditions for extraction of LCS, in terms of eigenvalues and eigenvectors of the Green strain tensor field, are presented. Application of these criteria to 2-D flows defines an LCS as an one-dimensional (1-D) material surface that is approximately normal to the eigenvector field associated with the largest eigenvalues of the Green strain tensor field. Therefore, the (hyperbolic) LCS is the strongest repelling or attracting surface among all nearby material surfaces. To differentiate among repelling and attracting ones, the repelling LCS at time t0 are obtained for going forward in time (τ > 0), while attracting LCS at t0 are obtained for going backward in time (τ < 0). In addition, Farazmand and Haller provided a computational algorithm to capture the LCS in a 2D flow.20 The identification of the LCS is a useful tool to provide assessment of transport fluxes across the LCS structures. In particular, integration of fluxes over the LCS, which is given by

with x denoting the position in the fluid domain, given as x ∈ Ω1 ⊂ n and the smooth velocity field v(x,t) being defined over the fluid domain Ωl. The system given by eq 5 is evolved over a time interval of t = [t0,t0 + τ]. Solution of the system described by eq 5 is the trajectory x(t;x0, t0) of a particle at time t initially (i.e., at t = t0), located at x0. By fixing the initial time t0 and the final time t of the dynamical system evolution, the flow map ϕtt0 is defined as the map which 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: 0

ξ2(x0 , t0 , τ )

where

with Newton’s method being used to obtain the solution. (3) Return to step (1). 2.3. Lagrangian Coherent Structures. The standard formulation of the Lagrangian Coherent Structures (LCS) identification methodology usually assumes a general timedependent fluid flow problem. Namely, the stationary fluid flow systems reveal time-invariant LCS structures, while periodic or fully time-varying flows with possible transfer to chaotic flows depict more interesting features of the LCS evolution. In the following discussion, we describe the method of extracting the LCS using the method of finite-time Lyapunov exponents (FTLE). In particular, since the FTLE method is a kinematic approach, we must compute a flow map that maps the particles forward in time along their trajectories. Differentiation of the flow map will provide a measure of separation of two nearby particles advected by the flow, so that we identify the LCS to be a material surface in the flow on which separation is locally maximal. Let the fluid flow system be considered as a dynamical system of the form

φtt : Ω l → Ω l : x0 → φtt (x0) = x(t ; x0 , t0)

(2x0) |τ |

δx0



x(t0) = x0

(8)

where δx0 is aligned with ξ2(x0,t0,τ). The aforementioned expression for the maximal stretch perturbation of the fluid element (eq 8) can be written as



dx = v(x , t ), dt

λ 2(x0 , t0 , τ ) δx0

∫LCS ddLt ds

(11)

where the integral is taken over the length of the LCS, where L(x,t) is the LCS material line (surface), which is obtained as the level set L(x,t) = 0 and provides insight into quantitative transport features. In addition, as it will be demonstrated in the Results and Discussion section, the intersection of repelling and attracting LCS reveals a homoclinic tangle, and the time evolution of the LCS describes the continuous motion of the lobes defined by these intersections. In other words, some regions of the fluid are transported by the mechanism of lobe dynamics.

(7)

with two real positive eigenvalues in the case of 2-D fluid flow λ1(x0, t0, τ) ≤ λ2(x0, t0, τ) and associated eigenvectors D

dx.doi.org/10.1021/ie503030z | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Rayleigh number (Ra = 4.5 × 104), yields a relatively uniform planar (straight line) interface between solid−liquid and formation of two similar vortexes coexisting within the melt domain. In contrast, the more thermally intensive process given by the higher Rayleigh number (Ra = 106) is characterized by the more intensive heat transfer, and yields completely timevarying aperiodic flow with the parabolic type of solid/melt line interface. From the crystal growth practitioner view, the higher Ra value would be more desired; however, the problem is that the thermal gradients across the same height of the solidified crystal are very large and may impact the grown crystal quality. In other words, the material point at the solid/melt interface and at x = 0.5 and corresponding solidified material points along the horizontal direction and at the same height generate very large temperature gradients, which might be detrimental for the grown crystal. In this sense, there is a tradeoff between the mixing and the quality of produced grown crystal. Along the same line, when the processing crystal operational conditions favored a less-intensive process, Figure 5 depicts a sequence of the attracting (blue) and repelling (red) LCS for Ra = 4.5 × 104 for which flow regime is characterized by symmetric partition of two vortexes. There is an LCS at the symmetry line, as is expected, at which repelling and attracting LCS material lines almost coincide. Another attracting LCS is formed at the lower part of the crucible, and the repelling LCS formed below the liquid/solid interface and expanded toward the sides and bottom of the crucible; then, it moves downward (see Figure 5). Since there is no melt flow across the LCS, because of transport barriers, the movement of the LCS shown in Figure 5 governs the mixing of the melt in a symmetric, lessintensive thermal regime. Figure 6 shows a sequence of attracting LCSs (denoted by blue lines, which are obtained by integrating the flow map backward in time) and repelling LCSs (denoted by red lines, which are obtained by integrating the flow map forward in time) for Ra = 106, for which the flow regime is asymmetric and fully time-varying. Since the magnitude of velocity vectors is larger than that in the previous case (Ra = 4.5 × 104), the integration time (τ) is smaller. In general, the value for τ is selected in such a way that fluid element separations reveal the desired dynamics of the fluid system.10 For this flow regime, mixing is better accomplished, because of the fact that, initially, the melt flow pattern is not symmetric and the middle LCS does not exist. In this case, the dynamic lobe is formed and mixing takes a set of particles by the flow from one side of the crucible to the other one (see Figure 7), and, later, more intersections of the lobe set occur, which characterizes observable LCSs using the degree and extent of the melt mixing.

3. RESULTS AND DISCUSSION In the general formulation, the nondimensional parameters that define the problem are given in Table 1 with the associated Table 1. Nondimensional Parameters parameter

magnitude

Prandtl number, Pr Stefan number, Ste L/W LΔT/W

0.01 1 2 0.25

values used in this work. Theses parameters represent pure gallium, whose properties in solid and liquid phases are the same, and have been taken from ref 17, where the solidification of this material in an inverted Bridgman apparatus have been investigated. The computational domain is the rectangle [0,1] × [0,2]. A 64 × 128 grid is used for the simulation. Since we use quadratic discretization, this grid leads to 74 691 and 66 306 degrees of freedom for the Navier−Stokes equation and the continuity and energy formulations, respectively. Initial conditions for the simulations include stationary fluid with temperature Th and solid with temperature Tc with a flat interface and Ra is applied at t = 0, leading to evolution of the solution with a nondimensional time step of Δt = 6.8 × 10−4. In order to compute the flow map and deformation tensor, a discrete uniform grid of computational points over the melt domain is defined. A rule of thumb is that the computational grid for the FTLE field should be fine enough to account for the high FTLE growth rate near the LCS.10 Therefore, the point fluid element at each grid point in the computational domain is advected by the flow map over the time interval τ resulting in the flow map ϕtt00+τ. The fluid element trajectory is determined from the velocity data (provided by the solution of the finite-element model) using a fourth-order Runge−Kutta explicit integration algorithm. For spatial interpolation of discrete velocity data standard linear finite element interpolation functions are used, while the gradients of the flow map were computed by finite differencing. Having the 2-D flow map, computation of the deformation tensor and corresponding eigenvalues is straightforward. Figure 4 shows representatives of FTLE fields for different Rayleigh numbers. As can be noticed, the numerical simulation reveals two different processing configurations of the Bridgman growth process. The first one, which is less temperature and heat-transfer intensive and is characterized by the lower

4. CONCLUSION In this paper, we demonstrate a novel approach for understanding the interplay between geometric features and dynamics of mixing in the melt flow in the vertical Bridgman process. In particular, the finite-element discretization of a mixed formulation for the Navier−Stokes equation, and continuity and energy balances, provide highly resolved velocity data. These data are used to identify Lagrangian coherent structures (LCSs) as ridges of the finite-time Lyapunov exponent (FTLE) field, representing time-dependent transport barriers to the flow. The presence of attracting and repelling LCSs in the melt during the Bridgman process shows the existence of material surfaces that govern the mixing template

Figure 4. Representative backward FTLE field for Ra = 4.5 × 104 and τ = −150Δt (left) and forward FTLE field for Ra = 106 and τ = 20Δt (right). E

dx.doi.org/10.1021/ie503030z | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Figure 5. A sequence of attracting LCS (blue lines) and repelling LCS (red lines) for Ra = 4.5 × 104 and τ = ±150Δt shown at t = 800, 925, ..., 1300Δt.

Figure 6. A sequence of attracting LCS (denoted by blue lines) and repelling LCS (denoted by red lines) for Ra = 106 and τ = ±20Δt shown at t = 1101, 1104, ..., 1337Δt.



ACKNOWLEDGMENTS 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), through Discovery Grant No. 386508-2011, is gratefully acknowledged.

■ Figure 7. A sequence of attracting LCS (denoted by blue lines) and repelling LCS (denoted by red lines) for Ra = 106 and τ = ±20Δt, which are enclosed by the red and green lobes, which evolve from one side of the domain into another with constant intersections, with attracting and repelling LCSs.

of the melt. Incomplete mixing of the melt results in nonuniform distribution of the material in the manufactured crystal. This fact highly affects the physical properties of the product, which have significant importance in the semiconductor industry.



REFERENCES

(1) Ottino, J. M. The Kinematics of Mixing: Stretching, Chaos, and Transport; Cambridge University Press: Oxford, U.K., 1989. (2) Kasap, S.; Capper, P. Springer Handbook of Electronics and Photonics Materials; Springer: New York, 2006. (3) 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. (4) Pimputkar, S. M.; Ostrach, S. Convective effects in crystals grown from melt. J. Cryst. Growth 1981, 55, 614−646. (5) Carruthers, J.; Nassau, K. Nonmixing cells due to crucible rotation during Czochralski crystal growth. J. Appl. Phys. 1968, 39, 5205−5214. (6) Kobayashi, N.; Arizumi, T. Computational analysis of the flow in a crucible. J. Cryst. Growth 1975, 30, 177−184. (7) Haller, G.; Yuan, G. Lagrangian coherent structures and mixing in two-dimensional turbulence. Phys. D (Amsterdam, Neth.) 2000, 147, 352−370. (8) Haller, G. Lagrangian coherent structures from approximate velocity data. Phys. Fluids 2002, 14, 1851−1861. (9) Haller, G. Distinguished material surfaces and coherent structures in three-dimensional fluid flows. Phys. D (Amsterdam, Neth.) 2001, 149, 248−277. (10) Shadden, S. C.; Lekien, F.; Marsden, J. E. Definition and properties of Lagrangian coherent structures from finite-time

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. F

dx.doi.org/10.1021/ie503030z | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Lyapunov exponents in two-dimensional aperiodic flows. Phys. D (Amsterdam, Neth.) 2005, 212, 271−304. (11) Lekien, F.; Shadden, S. C.; Marsden, J. E. Lagrangian coherent structures in n-dimensional systems. J. Math. Phys. 2007, 48, 065404. (12) Haller, G. A Variational Theory of Hyperbolic Lagrangian Coherent Structures. Phys. D (Amsterdam, Neth.) 2011, 240, 574−598. (13) 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. (14) Shadden, S. C.; Astorino, M.; Gerbeau, J. F. Computational analysis of an aortic valve jet with Lagrangian coherent structures. Chaos 2010, 20, 017512. (15) Lekien, F.; Ross, S. D. The computation of finite-time Lyapunov exponents on unstructured meshes and for non-Euclidean manifolds. Chaos 2010, 20, 017505. (16) Izadi, M.; Dubljevic, S. Analysis of Melt Flow Mixing in Czochralski Crystal Growth Process. Ind. Eng. Chem. Res. 2012, 51, 8675−8683. (17) Kaenton, J.; Semma, E.; Timchenko, V.; El Ganaoui, M.; Leonardi, E.; de Vahl Davis, G. Effects of anisotropy and solid/liquid thermal conductivity ratio on flow instabilities during inverted Bridgman growth. Int. J. Heat Mass Transfer 2004, 47, 3403−3413. (18) Belhamadia, Y.; Kane, A.; Fortin, A. An Enhanced Mathematical Model for Phase Change Problems with Natural Convection. Int. J. Numer. Anal. Model., Ser. B 2012, 3, 192−206. (19) Farazmand, M.; Haller, G. Erratum and addendum to “A variational theory of hyperbolic Lagrangian coherent structures” [Physica D 240 (2011) 574−598]. Phys. D (Amsterdam, Neth.) 2012, 241, 439−441. (20) Farazmand, M.; Haller, G. Computing Lagrangian coherent structures from their variational theory. Chaos 2012, 22, 013128.

G

dx.doi.org/10.1021/ie503030z | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX