Materials Science and Engineering A292 (2000) 219 – 223 www.elsevier.com/locate/msea
Boundary element analysis of the heat transfer in Bridgman growth process of semi-transparent crystals Wen-Qiang Lu Department of Physics, Di6ision of Thermal Science, The Graduate School at Beijing, Uni6ersity of Science and Technology of China, Academia Sinica, PO Box 3908, Beijing 100039, People’s Republic of China
Abstract The dual reciprocity boundary element method (DRBEM) has been developed to solve unsteady heat transfer problems with phase change moving interface and non-linear thermophysical properties. Numerical simulation of the heat transfer in Bridgman growth process of semi-transparent crystals has been made. Some results of transient temperature fields and time marching solid-liquid phase-change moving interface in the crystal growth are presented in this paper. © 2000 Elsevier Science S.A. All rights reserved. Keywords: Semi-transparent crystal growth; Phase change moving interface; Nonlinear thermophysical properties; DRBEM
1. Introduction During crystal growth process, the shapes of the solid–liquid interface and crystal-pulling rate have important effects on the defect structure and stress in the final crystal. As is now well known, a concave solid-liquid interface usually produces a polycrystalline ingot. Recently, in order to deeply understand the essence of these phenomena, the numerical simulations of thermal process during crystallization process have been made [1 –6]. Since the role of internal radiation transfer within the solid and the melt phase during vertical Bridgman growth process of semi-transparent crystals can not be neglected, the diffusion approximation for radiation had been induced in papers [7,8]. A review of calculating the unsteady process with phase-change moving interface by finite difference methods had been made [9]. The methods are separated into both kinds: (Fig. 6) (1) the strong numerical methods are to directly employ the physical formulation of the process as the Stefan condition, calculating and locating phase change moving boundary and finding temperature profiles at each time step; (2) the weak numerical methods are to reformulate the problem in such a way that the Stefan condition is implicitly incorporated in a new form of equation, in which explicit attention to the nature of E-mail address:
[email protected] (W.-Q. Lu).
the moving boundary is avoided, for example, the enthalpy method, the heat integration method, the effective capacity method, the source based method, the apparent capacity method, and so on. Comparing method (1) with method (2), the former methods can correctly capture phase change moving interface, however the complicated schemes and much more computational time are required. For example, when calculate the phase change moving interface by strong finite difference methods, it is necessary to repeatedly calculate the body-fitted coordinates systems, therefore, it is very difficult for solving transient three-dimensional complex shape phase change moving interface problems. The success of finite element method (FEM) lies in their ability to handle complex geometries, however much more storage and computational requirements are necessary since the calculations are made in total domain. In papers [7,8], FEM is used to calculate this problem by some simplified assumptions. For example, the release of latent heat of solidification was considered to be negligible in paper [7]; this eliminates the need for moving boundary calculations. In paper [8], a simplified approximation of steady process had been made. The boundary element method (BEM) has provided potential advantages of storage and computational requirements over FDM and FEM. And the discretization of phase change moving interface can be directly made in BEM, therefore BEM is very well suited to the modeling of these kinds of problems.
0921-5093/00/$ - see front matter © 2000 Elsevier Science S.A. All rights reserved. PII: S0921-5093(00)00998-9
W.-Q. Lu / Materials Science and Engineering A292 (2000) 219–223
220
Recently, DRBEM is developed as an effective pure boundary integral method [10]. The method was successfully used to model heat conduction problems, however, the modeling problems of the heat conduction with phase change and nonlinear thermophysical properties have been received in a few papers [11,12]. In this paper, DRBEM has been developed to solve unsteady heat transfer problems with phase change moving interface and non-linear thermophysical properties. As an example, numerical simulation for Bridgman growth process of semi-transparent crystals is made. Since the method is a pure boundary integral method, the calculation and storage requirements have been greatly reduced. Some valuable results for the furnace design are presented in the paper.
In some kinds of semi-transparent crystal growth process, for example, Bridgman crystal growth process of slowly pulling rate in microgravity environment, the natural convection in the melt is negligible. The heat transfer in the crystal, melt, ampoule and surroundings is by conduction and radiation. In order to consider the energy transport by radiation in the crystal and melt, the diffusion approximation [7,8] has been used. In this model, the medium behaves like a conductor that has a thermal conductivity dependent on temperature as follows: k= k mol +k rad = k mol +
16n 2sT 3 3Ro
(1)
Where QL and U o are the latent heat and the velocity of the phase change, respectively.
3. The dual reciprocity boundary element method It is well known that the nonlinear right terms of Eq. (2) can be transformed into linear operator form by Kirchhoff’s transform [10]. The basic idea is to construct a new variable U= U(T)=
&
T
kl dT To
(U = Cl 92U (t
(T =9[kl 9T] (t
(2)
Where, rl, cl and kl are the density, the specific heat and the thermal conductivity in the lth domain including the crystal (cry), melt (mel) and ampoule (a). Along the outer surfaces of the ampoule (Gi ) the boundary conditions are assumed as the follow: ka(n i 9T)= hi [Ta − Tf]+ os[T 4a −T 4f ]
(3)
Where n i is the outward pointing normal along the considered surface Gi. hi is the heat transfer coefficient. o is the surface emissivity. Tf is the variable furnace temperature profile. The temperature fields on the interface between two adjacent domain (cry, mel, a) are continuous. On phase change interface, the Stefan condition is used:
(5)
Where Cl = kl /(rlcl ) is thermal diffusivity. To is an arbitrary temperature reference value. Further inducing another new variable t= ((tl /(t =Cl ), Eq. (5) are transformed into the follows [10]: (U = 92U (tl
(6)
Eq. (6) can be transformed into pure boundary integral equations by the dual reciprocity theorem [10] as follows: wiUi +
&
rad
Where k and k are the molecular and the radiative conductivity, respectively. n is the refractive index. s is the Stefan– Boltzmann constant and Ro is the Rosseland mean absorption coefficient. T is the temperature. Employing the above approximation assumption, the control equation can be written as transient heat conduction equations: rlcl
(4)
so that the right terms of Eq. (2) become Laplace operator form in the new variable.
2. Physical model
mol
kcry(n cry · 9Tcry)− kmel(n cry · 9Tmel)= rmelQLU o · n cry
N+L
G
q*UdG−
& &
= % aj (tj ) wiT. ij + j=1
T*qdG
G
G
q*T. j dG−
& G
T*qˆj dG
(7)
Discretizing the above equations and time differential term, the following algebraic equations are obtained:
1 C. + uu H U M + 1 − uq GQM + 1 Dt
=
n
1 C. − (1−uu )H U M + (1− uq )GQM Dt
(8)
Where, wi = g/2p, g is internal angle at point i. N and L are the numbers of boundary and internal nodes, respectively. C. = − (HT. − GQ. )F − 1. H and G are coefficient matrices. The elements of matrices U, Q, T. , Q. and F are U, q=(U/(n, T. , qˆ = (T. /(n and f, respectively. For two-dimensional problem: f= 1+r, T. =r 2/ 4+ r 3/9, qˆ = ((r/(n)(1/2 + r/3). Where r is the distance between node i and node j. uu and uq are parameters that position the values of U and q between time levels m and m+ 1, respectively. The Green function: T*= ln(1/r) for two-dimensional. y is Cartesian coordinates in the generating plane. q*=(T*/(n. Dtj = Cj Dt is the step value of the modified time variable at node j.
W.-Q. Lu / Materials Science and Engineering A292 (2000) 219–223
The corresponding boundary conditions (3) along the outer surfaces of the ampoule (Gi ) in the transform space can be transformed into: (U dU q= =n i · 9U = (n i · 9T) = ka(n i · 9T) (n dT = hi [Ta −Tf]+os[T 4a −T 4f ]
(9)
(10)
4. Iterative schemes At per time step, the convergence for relaxation iterative calculation of the moving interface by the Stefan condition (10) must be achieved. At per given moving interface, the iterative calculation of the nonlinear Eq. (8) must be convergent by the use of Newton–Raphson scheme. Eq. (8) are written in abbreviated form as L(c)= v(c), where c is the unknown. Newton–Raphson scheme:
Fig. 1. Computational domain.
Fig. 2. Variable furnace temperature profile.
P(c m)c m + 1 = P(c m)c m − L. (c m), L. (c)= L(c)− v(c),
P(c)= (L. (c)/(c,
+1 = From Eq. (8), ones easily obtain: when c m j m+1 m+1 m+1 q , when c j =Uj , P ij = Cij /Dtj +uuHij − −1 +1 +1 +1 /(Um ) + Cij ( U m − Um / uq Gij ( ( q m j j j j ) ( ( ( Dtj ) m+1 (U j ). Where, m+1 j
The corresponding condition of phase change interface (4) in the transform space is transformed into: qcry − qmel = rmelQLU o · n cry
221
((Dtj ) − 1 rj cj dkj kj dcj kj drj =− 3 − − m+1 m+1 m+1 +1 (U j k j Dt dT j cj dT j rj dT m j
when assuming the density and the specific heat as constant, ((Dtj ) − 1 rj cj 16n 2sT 2j =− 3 m+1 (U j k j Dt Ro +1 +1 For given heat flux, (q m /(U m = 0; on the outer j j surfaces of the ampoule (Gi ), +1 hj 4osT 3j (q m j = + +1 (U m kj kj j
5. Numerical results Bridgman growth of semi-transparent crystals is numerically simulated (see Fig. 1). The thermophysical properties of the materials (CaF2) are chosen as the follows. The density: rcry = rmel = 3200 kg/m3. The thermal conductivity: k mol W/m/K, k mol mel = 0.6 cry = 6 W/m/K. The specific heat: ccry =cmel = 880 J kg − 1 K − 1. Rosseland absorption coefficient: Rocry = 0.3 cm − 1, Romel = 3.0 cm − 1. The refractive index:n = 1.44. The melting temperature: Tm =1653K. The latent heat:QL = 318 kJ kg − 1. The thermophysical properties of ampoule are chosen as the follows. The density:ra = 1900 kg m − 3. The specific heat: ca = 1400 J kg − 1 K − 1. The emissivity: o= 0.81. The thermal conductivity:ka = 50 W m − 1 K − 1. Ampoule wall thickness: l =1 mm. Furnace temperature gradient: Grad = 20 K cm − 1. Gradient zone length:lgra =20 cm. The rate of change of the furnace temperature profile:Vf = 0.5 cm h − 1 (see Fig. 2). Initial bottom zone temperature:To = 1600 K. Crystal half width: lcry = 1.0 cm. Ambient heat transfer coefficients: h1 = 0.04 W cm − 2 K − 1, h2 = 0.07 W cm − 2 K − 1, h3 = 0.05 W cm − 2 K − 1. Figs. 3 and 4 show the initial time marching case of crystal-melt interface and the temperature distribution at t=3 h, respectively. Initially, the curvature of the interface gradually decreases, then gradually increases. The temperature distributions are dependent on applied furnace temperature gradient and thermal process in the ampoule. During time marching, the solid–melt interfaces are all almost convex to the melt since the conductivity of solid phase is greatly larger than of melt. However, at
222
W.-Q. Lu / Materials Science and Engineering A292 (2000) 219–223
quality in the nearly top of grown crystals, with a dendrite, gas bubbles and voids. The effects of different furnace temperature gradients (Grad = 20 K cm − 1,Grad = 30 K cm − 1) are numerically simulated. Fig. 6 shows the effect of the furnace temperature gradients on the temperature distributions at central axis at t= 10 h. The temperature gradients of the melt and crystal obviously increase as increasing the furnace temperature gradient. During initial time marching, the changes of the furnace temperature gradient can affect the shape of the melt–crystal interface (see Fig. 7). This is comparable to the steady results of the lager aspect ratio in paper [8]. Fig. 3. Initial time marching case of the phase change moving interface.
6. Conclusions DRBEM can better used to simulate the transient thermal process in Bridgman growth of semi-transparent crystals since it has pure boundary integral character. The furnace temperature gradient affects the
Fig. 4. The temperature fields at t =3 h.
Fig. 6. The temperature distribution on y-axis in different furnace temperature gradients at t= 10 h.
Fig. 5. The variable case of phase interface at the later period of the solidification process.
the later period of the solidification process, the interface gradually becomes concave (see Fig. 5) due to boundary effect. These results are comparable to numerical results in paper [7]. And it can also be an explanation of some experimental results of the poor
Fig. 7. The effect of different furnace temperature gradients on the shape of phase interface at 3 h.
W.-Q. Lu / Materials Science and Engineering A292 (2000) 219–223
temperature distributions of the crystal and melt. The interface shape depends on thermal process in the ampoule and boundary effects.
Acknowledgements The project is supported by National Natural Science Foundation of China.
References [1] P.M. Adornato, R.A. Brown, J. Cryst. Growth 80 (1987) 155. [2] M.C. Liang, C.W. Lan, J. Cryst. Growth 167 (1996) 320. [3] D.H. Kim, R.A. Brown, J. Cryst. Growth 114 (1991) 441.
.
223
[4] K. Hoai, K. Sonnenberg, H. Wenzl, J. Cryst. Growth 137 (1994) 59. [5] H. Ouyang, W. Shyy, Int. J. Heat Mass Transfer 39 (1996) 2039. [6] F. Dupret, N. Van den Bogaert, Modeling Bridgman and Czochralski Growth, in: D.T.J. Hurle (Ed.), Handbook of Crystal Growth, vol. 2b, North-Holland, New York, 1994. [7] F. Barvinschi, I. Nicoara, J.L. Santailer, et al., Modelling Simul. Mater. Sci. Eng. 6 (1998) 691. [8] D. Vizman, I. Nicoara, D. Nicoara, J. Cryst. Growth 169 (1996) 161. [9] H. H.U., S.A. Argyropoulos, Modelling Simul. Mater. Sci. Eng. 4 (1996) 371. [10] P.W. Partridge, C.A. Brebbia, L.C. Wrobel, The Dual Reciprocity Boundary Element Method, Elsevier Applied Science, London, 1992, p. 175. [11] L.C. Wrobel, C.A. Brebbia, Comput. Methods Appl. Mech. Eng. 65 (1987) 147. [12] B. Sarler, G. Kuhn, Eng. Anal. Bound. Elem. 21 (1998) 53.