- Email: [email protected]
A Reduced Model for Fast Simulation of Crystal Growth
7b-045
Copyright © 1996 IFAC 13th Triennial World Congress, San Francisco, USA
A REDUCED MODEL FOR FAST SIMULATION OF CRYSTAL GROWTH E. Olivari, N. Van den Bogaert, V. Wertz and F. Dupret Unite de Mecanique Appliquee. CESAME. Universite catholique de Louvain, Al'. C. Lemaftre, 4. 8-1348 Louvain-la-Neuve. Belgium. (e-mail: [email protected])
Abstract : A reduced dynamic simulation model is presented. which predicts the global heat transfer, the shape of the solidification front. and the evolution of the crystal radius, for single crystals grown by the Czochralski process. This model is robust, accurate and fast, which allows for efficiently using its predictive capability during the growth process. Keywords: crystal growth. modelling. model based regulation. numerical simulation
I. INTRODUCTION
Process control is a key issue in engineering sciences, which is often solved by building up a differential model of the system in order to capture its governing transients (Clarke, 1994; Morari and Zafiriou, 1989). A major obstacle to applying this method in bulk crystal growth lies in the fact that most relevant quantities. which characterise crystal quality. cannot be measured during growth. in view of the very high temperatures required to melt most materials, and of the extreme sensitivity of the crystal lattice with respect to the environment of the liquid-solid interface (Hurle et al.. 1994). This forbids any use of sensors in this region. Nevertheless, for several years now. detailed simulation models are available, which can accurately describe the dynamical evolution of the growth (Brown. 1988; Dupret and Van den Bogaert, 1994). However, these models, which integrate PDEs using finite element methods, are extremely demanding in computing time and cannot be used when the crystal is grown. To circumvent this difficulty, a reduced model is described in this paper, which has been developed in order to predict the evolution of the temperature field in the overall furnace. together with the crystal shape. while working much faster than detailed models. During growth. the reduced model is continually acting as an observer, which approximates accurately the state of the system in order to facilitate the evaluation of non-measurable variables.
Global time-dependent models can predict the whole growth process from, as sole input. the geometry of the furnace. the material properties of its constituents. and the evolution of the processing parameters. Our method to accelerate the solution of such expensive detailed systems has consisted in (i) strongly reducing the number of unknowns characterising the thermal environment of the crystal and the melt. and (ii) performing as many complete simulations as possible before actual growth. Hence. the reduced model is based on using a projection technique. in order to decrease the system size, together with a data bank, in order to recover the resulh of detailed simulations performed beforehand. The difficulty of simulating bulk crystal growth can be related to the combined non-linear effects of radiative transfer in the furnace - since most crystals are grown at high temperatures. which involves a significant emission in the infra-red range - , and convection in the melt, as induced by the high temperature gradients prevailing in the liquid during growth. In the sequel. we will focus on the Czochralski (Cz) growth of semi-conductors, although similar methods can be elaborated for other processes and materials. In Cz growth (Fig. I). radiative exchanges between the constituents play a major role. and only global models (taking radiation into account) can provide significant results. Our reduced model is based on the same principles
6313
6314
can be "radiative" (the furnace enclosures). "twodimensional" (e.g .. the solid bodies), or "one-dimensional" (the thin shells). The melt-crystal is a particular twodimensional macro-element, whose shape is a priori unknown, and where latent heat is released at the solidliquid interface. The set of interfaces connecting the macroelements constitutes the skeleton of the global domain. At any new time step. each macro-element includes a set of unknowns. which are the new nodal temperatures, represented by the vector T. the incremental generalised nodal heat fluxes on the skeleton. represented by the vector Q+. and the geometrical unknowns. represented by nodal positions X. The nodal temperatures can be split into two parts. namely T+. which represents the temperatures on the skeleton. and T*. which represents all other temperatures. The equations of the detailed model will now be described for each kind of macro-element. In the absence of melt convection. the equations governing heat transfer in the melt-crystal macro-element can be put in the following compact form : A +(X) T+ + Q+ = E+(X)
(I)
T* = C*(X) + A *(X) T+ .
(2)
where eq. (I) relates the unknowns located on the skeleton, while eq. (2) expresses the internal temperatures as functions of the nodal temperatures on the skeleton. Note that matrices and column vectors in these equations depend on X, which signifies the dependence of heat transfer upon the geometry (position of solid-liquid interface, for example). The solidification heat release. proportional to the pull rate and to the time evolution of the solid-liquid interface position, is included in matrices C* and E+. To complete the description of this macro-element, geometrical equations must be considered :
'f (X) = 0 .
(3)
These latter include the discrete form of the Young-Laplace equation (which characterises the meniscus shape), the contact angle constancy condition - which is imposed between crystal wall and melt surface - , the time evolution of the interface radius, and the mass conservation. Finally. the crystal-melt interface position is constrained to be the melting isotherm : (4)
where the vector Tint is a subset of T and represents the nodal temperatures along the solidification front, while T m denotes the thermodynamic melting temperature. When convection in the melt is taken into account, the equations governing the crystal-melt macro-element become much more non-linear. A method for
approximating convection in dynamic simulations is developed in Dupret and Van den Bogaert (1994). Briefly, coupled global quasi-steady solutions are calculated at several stages of growth. The velocity field is stored in the form of the Stokes stream function, and subsequently used in dynamic simulations by interpolating between these precalculated results. Hence. only the level of the melt surface with respect to the crucible (and possibly the crystal and crucible rotation rates) play a role when evaluating the flow field. This means that eqs. (I) and (2) keep the same form, while the dependence of the matrices with respect to geometrical unknowns is more complex. Heat transfer in the other two-dimensional macro-elements can be summarised as follows: A + T+ + Q + = E +(W) ,
(5)
T* = C*(W) + A* T+.
(6)
which form is similar to that of the melt-crystal element, except that no geometrical coupling exists and that the heater power W is considered. For one-dimensional macroelements, an equation similar to (5) is obtained. Radiative exchanges relate the tluxes to the fourth powers of nodal temperatures along the boundary of the enclosures: Q+ = B(X) F(X) T+4
.
(7)
This mode of heat transfer is dependent upon the viewed and hidden parts of the furnace, which are represented in F. Hence, this matrix depends on geometrical unknowns, such as the crystal radius and the melt surface shape. Eqs. (I) to (7) fully describe the Cz system, as represented in the detailed model. The total number of nodal temperature unknowns on the skeleton is N. Equations are non-linear, due to the fourth power of the temperatures and due to the geometrical dependence of various matrices and vectors. The numerical technique for solving this system consists in using a decoupled scheme, in which eqs. (l) (2) (3) (4) on the one hand, and eqs. (1) (5) (7) on the other hand, are solved alternatively (Dupret et al., 1990; Ryckmans et a!., 1990). The Newton-Raphson method is used to solve each subsystem. Eq. (6) is solved subsequently. Although a decoupled scheme is used, the system size is so large that a huge amount of computer time is necessary for solving the problem. The cost is higher when performing time-dependent calculations with an unknown crystal radius, since two or three radiative matrices F - corresponding to different positions of the a priori unknown tri-junction - must be computed at every time step. The objective of the reduced model is to build up a system of equations whose number of unknowns is much lower than that of the detailed model. The idea is to provide a robust, efficient and accurate time-dependent prediction, so
6315
6316
6317
6318
Copyright © 1996 IFAC 13th Triennial World Congress, San Francisco, USA
A REDUCED MODEL FOR FAST SIMULATION OF CRYSTAL GROWTH E. Olivari, N. Van den Bogaert, V. Wertz and F. Dupret Unite de Mecanique Appliquee. CESAME. Universite catholique de Louvain, Al'. C. Lemaftre, 4. 8-1348 Louvain-la-Neuve. Belgium. (e-mail: [email protected])
Abstract : A reduced dynamic simulation model is presented. which predicts the global heat transfer, the shape of the solidification front. and the evolution of the crystal radius, for single crystals grown by the Czochralski process. This model is robust, accurate and fast, which allows for efficiently using its predictive capability during the growth process. Keywords: crystal growth. modelling. model based regulation. numerical simulation
I. INTRODUCTION
Process control is a key issue in engineering sciences, which is often solved by building up a differential model of the system in order to capture its governing transients (Clarke, 1994; Morari and Zafiriou, 1989). A major obstacle to applying this method in bulk crystal growth lies in the fact that most relevant quantities. which characterise crystal quality. cannot be measured during growth. in view of the very high temperatures required to melt most materials, and of the extreme sensitivity of the crystal lattice with respect to the environment of the liquid-solid interface (Hurle et al.. 1994). This forbids any use of sensors in this region. Nevertheless, for several years now. detailed simulation models are available, which can accurately describe the dynamical evolution of the growth (Brown. 1988; Dupret and Van den Bogaert, 1994). However, these models, which integrate PDEs using finite element methods, are extremely demanding in computing time and cannot be used when the crystal is grown. To circumvent this difficulty, a reduced model is described in this paper, which has been developed in order to predict the evolution of the temperature field in the overall furnace. together with the crystal shape. while working much faster than detailed models. During growth. the reduced model is continually acting as an observer, which approximates accurately the state of the system in order to facilitate the evaluation of non-measurable variables.
Global time-dependent models can predict the whole growth process from, as sole input. the geometry of the furnace. the material properties of its constituents. and the evolution of the processing parameters. Our method to accelerate the solution of such expensive detailed systems has consisted in (i) strongly reducing the number of unknowns characterising the thermal environment of the crystal and the melt. and (ii) performing as many complete simulations as possible before actual growth. Hence. the reduced model is based on using a projection technique. in order to decrease the system size, together with a data bank, in order to recover the resulh of detailed simulations performed beforehand. The difficulty of simulating bulk crystal growth can be related to the combined non-linear effects of radiative transfer in the furnace - since most crystals are grown at high temperatures. which involves a significant emission in the infra-red range - , and convection in the melt, as induced by the high temperature gradients prevailing in the liquid during growth. In the sequel. we will focus on the Czochralski (Cz) growth of semi-conductors, although similar methods can be elaborated for other processes and materials. In Cz growth (Fig. I). radiative exchanges between the constituents play a major role. and only global models (taking radiation into account) can provide significant results. Our reduced model is based on the same principles
6313
6314
can be "radiative" (the furnace enclosures). "twodimensional" (e.g .. the solid bodies), or "one-dimensional" (the thin shells). The melt-crystal is a particular twodimensional macro-element, whose shape is a priori unknown, and where latent heat is released at the solidliquid interface. The set of interfaces connecting the macroelements constitutes the skeleton of the global domain. At any new time step. each macro-element includes a set of unknowns. which are the new nodal temperatures, represented by the vector T. the incremental generalised nodal heat fluxes on the skeleton. represented by the vector Q+. and the geometrical unknowns. represented by nodal positions X. The nodal temperatures can be split into two parts. namely T+. which represents the temperatures on the skeleton. and T*. which represents all other temperatures. The equations of the detailed model will now be described for each kind of macro-element. In the absence of melt convection. the equations governing heat transfer in the melt-crystal macro-element can be put in the following compact form : A +(X) T+ + Q+ = E+(X)
(I)
T* = C*(X) + A *(X) T+ .
(2)
where eq. (I) relates the unknowns located on the skeleton, while eq. (2) expresses the internal temperatures as functions of the nodal temperatures on the skeleton. Note that matrices and column vectors in these equations depend on X, which signifies the dependence of heat transfer upon the geometry (position of solid-liquid interface, for example). The solidification heat release. proportional to the pull rate and to the time evolution of the solid-liquid interface position, is included in matrices C* and E+. To complete the description of this macro-element, geometrical equations must be considered :
'f (X) = 0 .
(3)
These latter include the discrete form of the Young-Laplace equation (which characterises the meniscus shape), the contact angle constancy condition - which is imposed between crystal wall and melt surface - , the time evolution of the interface radius, and the mass conservation. Finally. the crystal-melt interface position is constrained to be the melting isotherm : (4)
where the vector Tint is a subset of T and represents the nodal temperatures along the solidification front, while T m denotes the thermodynamic melting temperature. When convection in the melt is taken into account, the equations governing the crystal-melt macro-element become much more non-linear. A method for
approximating convection in dynamic simulations is developed in Dupret and Van den Bogaert (1994). Briefly, coupled global quasi-steady solutions are calculated at several stages of growth. The velocity field is stored in the form of the Stokes stream function, and subsequently used in dynamic simulations by interpolating between these precalculated results. Hence. only the level of the melt surface with respect to the crucible (and possibly the crystal and crucible rotation rates) play a role when evaluating the flow field. This means that eqs. (I) and (2) keep the same form, while the dependence of the matrices with respect to geometrical unknowns is more complex. Heat transfer in the other two-dimensional macro-elements can be summarised as follows: A + T+ + Q + = E +(W) ,
(5)
T* = C*(W) + A* T+.
(6)
which form is similar to that of the melt-crystal element, except that no geometrical coupling exists and that the heater power W is considered. For one-dimensional macroelements, an equation similar to (5) is obtained. Radiative exchanges relate the tluxes to the fourth powers of nodal temperatures along the boundary of the enclosures: Q+ = B(X) F(X) T+4
.
(7)
This mode of heat transfer is dependent upon the viewed and hidden parts of the furnace, which are represented in F. Hence, this matrix depends on geometrical unknowns, such as the crystal radius and the melt surface shape. Eqs. (I) to (7) fully describe the Cz system, as represented in the detailed model. The total number of nodal temperature unknowns on the skeleton is N. Equations are non-linear, due to the fourth power of the temperatures and due to the geometrical dependence of various matrices and vectors. The numerical technique for solving this system consists in using a decoupled scheme, in which eqs. (l) (2) (3) (4) on the one hand, and eqs. (1) (5) (7) on the other hand, are solved alternatively (Dupret et al., 1990; Ryckmans et a!., 1990). The Newton-Raphson method is used to solve each subsystem. Eq. (6) is solved subsequently. Although a decoupled scheme is used, the system size is so large that a huge amount of computer time is necessary for solving the problem. The cost is higher when performing time-dependent calculations with an unknown crystal radius, since two or three radiative matrices F - corresponding to different positions of the a priori unknown tri-junction - must be computed at every time step. The objective of the reduced model is to build up a system of equations whose number of unknowns is much lower than that of the detailed model. The idea is to provide a robust, efficient and accurate time-dependent prediction, so
6315
6316
6317
6318