- Email: [email protected]
Three-dimensional modelling of radial segregation due to weak convection
ARTICLE IN PRESS
Journal of Crystal Growth 269 (2004) 454–463
Three-dimensional modelling of radial segregation due to weak convection . Carlbergb Minh Do-Quanga,*, Gustav Amberga, Torbjorn a b
Department of Mechanics, Royal Institute of Technology, Osquars Backe 18, 100 44 Stockholm, Sweden Department of Engineering, Physics and Mathematics, Mid Sweden University, 851 70 Sundsvall, Sweden Received 13 November 2003; accepted 17 May 2004 Available online 6 July 2004 Communicated by K.W. Benz
Abstract A comprehensive three-dimensional, time-dependent model of heat, momentum and solute transfer during solidification is carried out to illustrate the influence of weak convection, caused by surface tension forces, on radial dopant segregation occurring in crystal growth under micro-gravity conditions. 3D adaptive finite element method is used in order to simulate the motion and deformation of the solidification interface. The geometry studied is a Bridgman configuration with a partly coated surface. The small slots in the coating gives a free surface in a controlled way, and is varied in order to alter the Marangoni flow. In this study, A comparison is made between the numerical results and the experimental results. A good agreement has been observed for the effective distribution coefficient keff and for the radial segregation Dc0 : The radial dopant segregation is affected by weak convection. r 2004 Elsevier B.V. All rights reserved. PACS: 64.75.+g; 81.10.Mx; 44.25.+f; 02.60.Cb; 02.70.Dh Keywords: A1. 3D numerical simulation; A1. Adaptive finite element; A1. Automated code generation; A1. Radial segregation; A2. Floating zone technique; A2. Microgravity conditions
1. Introduction Several theoretical studies have shown that weak convection, with flow rates of the same *Corresponding Author. Tel.:+46-8-790-6871; fax:+46-8796-9850. E-mail addresses: [email protected] (M. Do-Quang), [email protected] (G. Amberg), [email protected] (T. Carlberg).
order of magnitude as the growth rate, can give strong radial segregation [1–4]. Such convection can be expected to occur in space experiments if parameters like temperature gradients, sample dimensions and growth rates are unfavourably chosen. Some experimental indications have been obtained in Soviet experiments referred to in Ref. [4], and similar effects have been observed in samples with weak Marangoni convection [5]. There it was
0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.05.076
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
455
shown that very strong radial inhomogeneities occurred if the flow rates were of the same order of magnitude as the growth rate. These results can, however, only be regarded as qualitative verifications of the phenomenon. Recently, experiments were done where thermocapillarity was used to obtain controlled weak flows [6]. In the present paper results from numerical simulations of the experimental set up is presented. The aim is to obtain increased knowledge of the coupling between weak flow and radial segregation, both in this particular case and in a more general manner. The project is the first phase of a research program, which will continue with experiments in the ISS, during which the effect of residual gravity on the radial segregation will be investigated.
2. Experimental results Details of the experimental technique, ground based work and analysis of thermal field has been published in Ref. [7]. Preliminary results and more extensively treated results have been published in Refs. [8,6], respectively. The aim of the experiments were to study effects of weak convection on radial segregation of tin in antimony samples. The experiments were performed with coated samples, which had slots in the coating of different widths, and the free surfaces thus formed gave thermocapillary flow varying with the slot size. The slots were positioned on one side of the samples to give flow vortices extending over the whole sample diameter and thus an asymmetric dopant distribution. The sample diameters were 7 mm and a length of 15 mm were remelted and solidified during an experimental run. During a flight in the STS-108 in December 2001 seven experiments were performed. They had the following percentages of free surface: z ¼ 0; 10%, 15%, 25%, 40%, 55% and 75%. The chemical composition of the samples were analysed by a wet chemical technique, i.e. dissolution followed by plasma atomic emission spectroscopy (ICP-AES). In Fig. 1 axial concentration profiles for three of the samples are shown together with the composi-
Fig. 1. Axial concentration profiles for three space experiments and for a rod before remelting.
tion in a rod before processing. The data points on position 1 mm correspond to the starting material concentration, and are actually measured in a slice in the not melted part close to the remelt boundary. From these curves approximate values of effective distribution coefficients can be obtained from the start of solidification by dividing the concentration at the beginning of growth (extrapolated) with the content of the starting material. A typical appearance for all three curves is a relatively strongly increased concentration along the sample length to a position around 7– 9 mm: At a sample length of about 10 mm the power of the furnaces was turned off, i.e. it is only the first 10 mm that are solidified in a controlled way. The decreased concentration towards the end reveals a build up of solute in the centre of the samples. In Fig. 2 radial concentration profiles for the same three space experiments are shown together with a ground reference sample. The slot side is to the left in the figure. The space samples have a clearly stronger radial segregation than the reference sample. For 10% free surface the concentration at the slot side is high, i.e. maximum concentration is on the slot side. There is a weak tendency that at increasing fractions of free
ARTICLE IN PRESS 456
M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
L
αL
z
Free surface
Coated surface
x Thermocouple Fig. 2. Radial concentration profiles.
y
surface, the area of high concentration is shifted to the side opposite of the slots.
3. Mathematical modelling In the experiments with partially coated surfaces, the surface consisted of several slots of free surface and this is used in order to control the different levels of thermo-capillary convection. The dimensions of the slots were much smaller than the size of the molten zone and occupied one quarter of the periphery of the molten zone, Fig. 3. The temperature is fully controlled by a feedback system in the experiments so that the solidification interface moved up with a constant speed (Us ¼ 1:5 mm min1 ) during the cooling down stage, and passed each and every of the free slots, one by one. Therefore, accounting for solidification of a binary alloy, the movement of the solid– liquid interface, and the evolution of latent heat near the interface, the numerical method must consider time-dependent, three-dimensional flow in the entire zone. A new version of femLego [9], a three-dimensional parallel adaptive finite element method, therefore, was developed within the scope of this project. This tool was based on the automated code generation [10].
Fig. 3. Geometry of the floating zone configuration.
3.1. Governing equations In space, the only coupling between the energy equation and the momentum equations is through the boundary conditions at the free surface. Therefore, the governing transport dimensionless equations where the melt is treated as an incompressible Newtonian liquid read as r ðw~ u Þ ¼ 0;
ð1Þ
q~ u 1 n þ ð~ u rÞ~ u ¼ rp þ w~ u; r2 ~ uþ qt ReTC HðwÞ
ð2Þ
qy 1 1 qw þ ð~ u rÞy ¼ ; ðr2 yÞ þ qt ReTC Pr Ste qt
ð3Þ
where, ~ u is the fluid velocity, p is the pressure and y is the temperature. The extra term in the Navier– Stokes equation, Eq. (2), represents the porous medium model at the solid–liquid interface. w is the fraction of volume occupied by melt, w has a value of 0 in the solid region, 1 in the completely molten region. HðwÞ is the permeability of the mush [11]. HðwÞ is infinite in the molten region and very rapidly decreases to very small values (say, 108 ) in the solid region, so that the velocity ~ u become effectively zero in the solid region.
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
The dimensionless parameters occurring in Eqs. (1)–(3) are the thermo-capillary Reynolds number ReTC that based on the thermo-capillary velocity U 0 ¼ gDTm1 ; the Prandtl number Pr and the Stefan number. They are defined as gDTL ; ð4Þ ReTC ¼ mn n Pr ¼ ; a
ð5Þ
CpDT ; ð6Þ L where, g is a characteristic value of the coefficient of surface tension, m and n are dynamic and kinematic viscosities, a is the thermal diffusivity. In this study, the thermal conductivity k of melt and solid are the same and shown in Table 1. The fraction of volume that is occupied by solid can be computed, according to Amberg [11], qw ¼ k y yliq : ð7Þ qt Here, k is a numerical parameter and chosen k ¼ 2=Ste in order to satisfy the numerical stability problem. This equation allows a simple way of determination of the solid fraction w at the new time level by using the discretized form of the lefthand side: If w tends to increase beyond 1 or decrease below 0; the new value of w is taken to be 1 or 0; respectively. This means that, y is then allowed to differ from the nondimensionless of the Ste ¼
457
melting temperature yliq ; which is, of course, quite correct in the solid (w ¼ 0) or the liquid (w ¼ 1) region. The numerical parameter k is a small constant, used to amplify any deviation of the temperature from yliq and results in a rapid melting or freezing that restores y to yliq : The redistribution of dopant in the base material is obtained from the mass transfer equation. qcm 1 rcm : þ ð~ u rÞcl ¼ r ð8Þ Pe qt Here, cm is the net composition of the solid–liquid mixture, cl is the composition of the liquid. Pe ¼ U 0 L=D is the Peclet number with D is the molecular diffusivity of the dopant. This diffusion can be treated as partially solidified diffusion (Ds ) with a solute diffusion of the dopant (Dl ) as follows: D ¼ Ds ð1 wÞ þ Dl w:
ð9Þ
3.2. Boundary conditions Top surface. The temperature of the top surface is a function in time and follows the temperature history in the experiment. The dopant concentration at the top is a no-flux condition. qc ¼ 0: qz
ð10Þ
The velocities follow: uz ¼ 0
Table 1 Physical properties of Sn/Sb Physical property
Symbol
Value Sn/Sb
Density Liquidus temperature Kinetic viscosity Dynamic viscosity Thermal conductivity Specific heat Temp. coef. of surface tension Latent heat of fusion Solute diffusivity Solidified diffusivity Segregation coefficient of Sn in Sb
r ðkg m3 ) Tliq ðKÞ n ðm2 s1 Þ m ðN sm2 Þ k ðW m1 K1 Þ Cp ðJ kg1 K1 ) dg=dT ðN m1 KÞ
6:684 903:9 1:9 107 1:2 103 21:3 260 5:0 105
L ðJ kg1 Þ Dl ðm2 s1 Þ Ds ðm2 s1 Þ k0
1:63 105 2:0 109 5 1015 0:3
and
qux quy ¼ ¼ 0: qz qz
ð11Þ
Bottom surface. The temperature of the bottom surface is also a function in time and following the experiment temperature. The top and the bottom temperature are shown in Fig. 4. The bottom position shall never be melted. Therefore, the velocities can be zero. The dopant concentration satisfies n% rc ¼ 0:
ð12Þ
Free surface. The gradient of surface tension g in tangential direction i along the free surface in the slots must be balanced by shear stress, the latter being proportional to the normal flux of the
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
458
1100 1000 900 Temperature˚C
800
Top surface
700 600
Bottom surface
500 400 300
Heating up stage
200
Cooling down stage
100 0 0
100 200 300 400 500 600 700 800 time (s)
Fig. 4. The experiment temperature at the top and the bottom surface of GAS experiment.
tangential velocity ~ ui m
q~ u i qg qy : ¼ qn qy qi
ð13Þ
Here, m is the dynamic viscosity of the fluid and qg=qy is the coefficient of surface tension, n is the normal to the surface directed outward from the liquid region. The boundary condition for the dopant concentration as follows: n% rc ¼ 0:
ð14Þ
Coated surface. At the entirely coated surface the velocity obeys a no-slip boundary condition. n% rc ¼ 0;
~ u ¼ 0:
and Quanrtapelle [12]. The Poisson equation for the pressure is solved by using the well-known Conjugate-Gradient method. The convective term in Eqs. (2) and (3) and the dopant mass transfer equation (8) are calculated implicitly by using the GMRES method. Here, the streamline-upwind concept was used for the test functions in order to increase the stability [13]. In this study, the adaptive mesh was used in order to solve the diffusion problem at the solidification interface. The minimum of mesh size was limited at h ¼ 0:01 mm: Therefore, we would have at least 8 elements placed in the boundary layer due to diffusion, d0d ¼ 0:08 mm [14]. Generally, a mesh with about 1:2 million elements and 200 000 nodes was used in the present study. During each adaptive time step, there have been about 2000 elements that would be merged and 2000 new elements that would be created by the refinement stuff. All computations were performed on a PC-cluster with 16 CPUs, AMD Athlon MP 2400+ with 1GB memory for each. This cluster is a home made product of Mechanics Department of Royal Institute of Technology. The longest necessary computation for the whole cooling stage took about one week. A symbolic computational tool called femLego has been applied to create a parallel code that can solve the present time dependent, three-dimension mathematical problem.
ð15Þ
3.3. Numerical method In this study, a three-dimensional, parallel finite element method has been used in conjunction with an adaptive mesh refinement and coarsening. The equations have been discretised using a finite element approach on an unstructured grid. The type of elements which are used for both velocity and pressure variables are piecewise linear tetrahedral elements. The restrictions imposed by the Babuska–Brezzi condition are avoided by adding a pressure stabilisation term. The pressure and velocity solutions are split by using a fractional step algorithm that was introduced by Guermond
4. Results and discussion The problem studied in this paper is the numerical simulation of the space experiment (STS-108) that was performed in a Get Away Special (GAS) module in December 2001. Both experiment and numerical method focus on the radial segregation of tin in antimony samples due to weak convection. The physical properties of Sn/ Sb are listed in Table 1. The evolution of the axial dopant concentration distribution in a vertical slice at j ¼ 0 for case z ¼ 55% is shown in Fig. 5. The solidification interface was moving with a constant speed Usol E0:025 mm/s, exactly as it was programmed in the experiment [7]. That it to say, although the solidification interface was deformed
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
459
Fig. 5. Evolution of the concentration distributation at t ¼ 440 s (A), 580 s (B), 720 s (C) and 860 s (D). 55% free window.
Fig. 6. Velocity field at the vertical slice for the different free window, (A) z ¼ 10%; (B) z ¼ 25% and (C) z ¼ 55%:
by the convection, the speed of the interface was constant as expected. In this figure, local segregation close to the free windows appear as result of a complicated flow near the free surface. The other local segregation appears in the opposite side of the slots due to the recirculating fluid. This recirculating fluid together with the main flow field, swept dopant from the diffusion layer into the bulk liquid and induced radial segregation across the radial axis. That is the reason why the highest dopant concentration did not appear in the opposite side of the slots. They appeared instead in the middle of the sample (at r ¼ 0:5–2:5 mm).
A comparison of the flow field in a vertical slice for the different free window surface: z ¼ 10%; 25% and 55% is shown in Fig. 6. In this figure, the fluid flow is accelerated at the free window surface and decelerated at the coated surface. The overall flow fields are similar to the surface velocity solution for the mixed boundary, see [5]. In case of strong convection flow, the shape of the solidification interface is bent. Specially, at the early stage of cooling, when the solidification interface is placed around z ¼ 3 mm; the boundary at the partly coated side was not the well mixed region any more. It became the boundary layer of the ‘‘flat
ARTICLE IN PRESS 460
M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
Fig. 7. The concentration distributation at t ¼ 860 s with the different free window z ¼ 0% (A); z ¼ 10% (B); z ¼ 25% (C) and z ¼ 55% (D). 1.05 Effective distribution coefficient
plate’’ with an input flow from up to down. When the flow was as strong as in the case of 55% free window, Fig. 6C, the area that is close to the ‘‘flat plate’’ appears as a recirculating fluid. This is the reason why the solidification interface in Fig. 5A is more blend than it is in the cases (B)–(D). Fig. 7 shows a dopant field computed for the segregation coefficient (ko ¼ 0:3) of tin in antimony at different free window surface fraction: z ¼ 0%; 10%; 25% and 55%: Obviously, the maximum dopant concentration in the crystal was shifted from the slot side to the other side, as the convection increases. The radial segregation can be observed in Fig. 7. It is totally influenced by the convection that was controlled by the fraction of free window surface. A similar shift of maximum dopant concentration was also observed in the experimental results, Fig. 2. However, a direct quantitative agreement with the computed ditributions is not obtained. There are several possible reasons to weaknesses in the experimental results, e.g. the free surfaces at the slots might have some contamination influencing the surface tension coefficient. The measurements are also made with a method that cannot reveal all the small scale variations of concentration variations that the computations predict. The experimental results in Fig. 1, showing increased concentration along the axial direction can also qualitatively get some support from the calculated concentration
Numerical simulation Experiment
1
0.95
0.9
0.85
0.8 0
10
20
30
40
50
60
70
80
Free surface area [%]
Fig. 8. Comparison of the effective distribution coefficient between numerical and experiment data.
fields shown in Fig. 7, but also here the experimental data are more pronounced than the calculated. The effective distribution coefficients keff was approximated as follows [2]: keff ¼
cS ; cinit
ð16Þ
where cS was the average concentration across the crystal and cinit is the initial concentration of the material, in dimensionless form cinit ¼ 1: The numerical results of the effective distribution
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
461
Fig. 9. The dopant profile at z ¼ 10 mm; (z ¼ 0%, 10%, 25% and 55%).
coefficients together with the treatment data from experiment are shown in Fig. 8. A good agreement between the numerical results and the experimental results has been observed. Both results have keff ranging from 0:83 to 1: Such data indicates that a very low convection level occurred in the space samples and it is obviously that keff decreases with increasing free surface windows fraction.
Analyzing the numerical data in the same way as the experimental, a horizontal slice with 3 mm thickness is produced to obtain the radial segregation profiles. The results of four different z; the percentage of free window surfaces over the partially coated surface, are shown in Fig. 9. The maximum dopant concentration for the Sn-droped Sb looks like an island and shifts from the slot side
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
462
1.05
Table 2 Comparison of the radial segregation due to the weak convection
1 0.95
Sn Conc.
0.9 0.85 0.8
zð%Þ
cmin
cmax
cavg
Dc0
10 25 55
0.8160 0.6399 0.5890
1.0146 1.0973 1.1830
0.9547 0.9470 0.8196
0.2080 0.4830 0.7248
0.75 0.7 0.65 0.6 0.55
00% free window 10% free window 25% free window 55% free window
-3
-2
-1
0 1 radius [mm]
2
3
Fig. 10. Variation of dopant concentration in the radial axis.
(j ¼ 0) to the entirely coated surface side (j ¼ p), as the convection increases. The islands appear as the levels of the thermo-capillary convection were not high enough to distribute the dopant over whole diameter of the sample. The maximum concentration island for 10% free window surface appears near the centre of slice cut and shifts far to the slot side in case of z ¼ 25%: In those two cases, the gradient in the dopant concentration is much sharper between the partially coated surface and the island of maximum concentration. It shows that, the influence of convective dopant transport is much more active, compressing the concentration lines. When the convection is strong enough, case z ¼ 55%; a recirculating flow appeared and drove the dopant from the entirely coated surface side to the centre. Therefore, the maximum dopant concentration was compressed in the radial axis direction and transported to the periphery at j ¼ p=2 and 3p=2: In this case the gradient in the dopant concentration was more or less the same between two sides of the island. The variation of dopant concentration along a diameter is shown in Fig. 10. The radial segregation defined as Dc0 ¼
cmax cmin ; cavg
ð17Þ
where, cmax ; cmin and cavg are maximum, minimum and average dopant concentration across the
horizontal slice of the sample. See Table 2. Obviously, the results listed in Table 2 show the influence of convection flow to the radial segregation. By increasing convection, the radial segregation will also be increased. At this point, the conclusion from the experimental result is different than the numerical result. In the experiment, cavg is dependent on the content of the starting material and is totally different in the various cases. Meanwhile, in the numerical method, the starting material remained constant for all cases.
5. Conclusions Melt growth experiments of semiconductor material using the gradient freeze technique were simulated. A corresponding, the space experiments were directly modelled. The moving solidification interface was captured using an adaptive mesh method with very fine resolution at the interface. When the convection was strong, the shape of the solidification interface was deformed and some recirculating flow regions appeared at the periphery. The radial segregation in the solidified material was obtained and a clear quantitative relation between the free surface area and the radial segregation was revealed. The dopant concentration results show a slight shift of the maximal concentration in the radial profiles.
References [1] S.A. Nikitin, V.L. Polezhayev, A.L. Fedyushkin, J. Crystal Growth 52 (1981) 471. [2] C.L. Chang, R.A. Brown, J. Crystal Growth 63 (1983) 343. [3] J.I.D. Alexander, J. Quazzini, F. Rosenberger, J. Crystal Growth 97 (1989) 285.
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463 [4] V.I. Polezhayev, S.A. Nikitin, Oxford, 1989, p. 237. [5] M. Levenstam, G. Amberg, E. Tillberg, T. Carlberg, J. Crystal Growth 104 (1990) 641. [6] T. Carlberg, in: International Symposium on Physical Sciences in Space, 2004, submitted for publication. [7] T. Carlberg, C. Winkler, G. Amberg, in: Proceedings of the First International Symposium on Microgravity Research & Applications in Physical Sciences & Biotechnology, Sorrento, 2000. . [8] P. Holm, K. Loth, B. Larsson, T. Carlberg, Proc. of the 16th ESA Symp. on European Rocket and Balloon
[9] [10] [11] [12] [13] [14]
463
Programmes and Related Research, St. Gallen, Switzerland, 2003, ESA SP-530, p. 113–117. M. Do-Quang, I. Loginova, G. Amberg, in preparation. . G. Amberg, R. Tonhardt, C. Winkler, Math. Comput. Simulation 49 (1999) 149. G. Amberg, Int. J. Heat Mass Transfer, 34 (1) (1991) 17. J.L. Guermond, L. Quanrtapelle, Int. J. Numer. Meth. Fluids. 26 (1998) 1039. O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, Vol. 3, Butterworth Heinemann, Oxford, 2000. C. Winkler, G. Amberg, T. Carlberg, J. Crystal Growth 210 (2000) 573.
Journal of Crystal Growth 269 (2004) 454–463
Three-dimensional modelling of radial segregation due to weak convection . Carlbergb Minh Do-Quanga,*, Gustav Amberga, Torbjorn a b
Department of Mechanics, Royal Institute of Technology, Osquars Backe 18, 100 44 Stockholm, Sweden Department of Engineering, Physics and Mathematics, Mid Sweden University, 851 70 Sundsvall, Sweden Received 13 November 2003; accepted 17 May 2004 Available online 6 July 2004 Communicated by K.W. Benz
Abstract A comprehensive three-dimensional, time-dependent model of heat, momentum and solute transfer during solidification is carried out to illustrate the influence of weak convection, caused by surface tension forces, on radial dopant segregation occurring in crystal growth under micro-gravity conditions. 3D adaptive finite element method is used in order to simulate the motion and deformation of the solidification interface. The geometry studied is a Bridgman configuration with a partly coated surface. The small slots in the coating gives a free surface in a controlled way, and is varied in order to alter the Marangoni flow. In this study, A comparison is made between the numerical results and the experimental results. A good agreement has been observed for the effective distribution coefficient keff and for the radial segregation Dc0 : The radial dopant segregation is affected by weak convection. r 2004 Elsevier B.V. All rights reserved. PACS: 64.75.+g; 81.10.Mx; 44.25.+f; 02.60.Cb; 02.70.Dh Keywords: A1. 3D numerical simulation; A1. Adaptive finite element; A1. Automated code generation; A1. Radial segregation; A2. Floating zone technique; A2. Microgravity conditions
1. Introduction Several theoretical studies have shown that weak convection, with flow rates of the same *Corresponding Author. Tel.:+46-8-790-6871; fax:+46-8796-9850. E-mail addresses: [email protected] (M. Do-Quang), [email protected] (G. Amberg), [email protected] (T. Carlberg).
order of magnitude as the growth rate, can give strong radial segregation [1–4]. Such convection can be expected to occur in space experiments if parameters like temperature gradients, sample dimensions and growth rates are unfavourably chosen. Some experimental indications have been obtained in Soviet experiments referred to in Ref. [4], and similar effects have been observed in samples with weak Marangoni convection [5]. There it was
0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.05.076
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
455
shown that very strong radial inhomogeneities occurred if the flow rates were of the same order of magnitude as the growth rate. These results can, however, only be regarded as qualitative verifications of the phenomenon. Recently, experiments were done where thermocapillarity was used to obtain controlled weak flows [6]. In the present paper results from numerical simulations of the experimental set up is presented. The aim is to obtain increased knowledge of the coupling between weak flow and radial segregation, both in this particular case and in a more general manner. The project is the first phase of a research program, which will continue with experiments in the ISS, during which the effect of residual gravity on the radial segregation will be investigated.
2. Experimental results Details of the experimental technique, ground based work and analysis of thermal field has been published in Ref. [7]. Preliminary results and more extensively treated results have been published in Refs. [8,6], respectively. The aim of the experiments were to study effects of weak convection on radial segregation of tin in antimony samples. The experiments were performed with coated samples, which had slots in the coating of different widths, and the free surfaces thus formed gave thermocapillary flow varying with the slot size. The slots were positioned on one side of the samples to give flow vortices extending over the whole sample diameter and thus an asymmetric dopant distribution. The sample diameters were 7 mm and a length of 15 mm were remelted and solidified during an experimental run. During a flight in the STS-108 in December 2001 seven experiments were performed. They had the following percentages of free surface: z ¼ 0; 10%, 15%, 25%, 40%, 55% and 75%. The chemical composition of the samples were analysed by a wet chemical technique, i.e. dissolution followed by plasma atomic emission spectroscopy (ICP-AES). In Fig. 1 axial concentration profiles for three of the samples are shown together with the composi-
Fig. 1. Axial concentration profiles for three space experiments and for a rod before remelting.
tion in a rod before processing. The data points on position 1 mm correspond to the starting material concentration, and are actually measured in a slice in the not melted part close to the remelt boundary. From these curves approximate values of effective distribution coefficients can be obtained from the start of solidification by dividing the concentration at the beginning of growth (extrapolated) with the content of the starting material. A typical appearance for all three curves is a relatively strongly increased concentration along the sample length to a position around 7– 9 mm: At a sample length of about 10 mm the power of the furnaces was turned off, i.e. it is only the first 10 mm that are solidified in a controlled way. The decreased concentration towards the end reveals a build up of solute in the centre of the samples. In Fig. 2 radial concentration profiles for the same three space experiments are shown together with a ground reference sample. The slot side is to the left in the figure. The space samples have a clearly stronger radial segregation than the reference sample. For 10% free surface the concentration at the slot side is high, i.e. maximum concentration is on the slot side. There is a weak tendency that at increasing fractions of free
ARTICLE IN PRESS 456
M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
L
αL
z
Free surface
Coated surface
x Thermocouple Fig. 2. Radial concentration profiles.
y
surface, the area of high concentration is shifted to the side opposite of the slots.
3. Mathematical modelling In the experiments with partially coated surfaces, the surface consisted of several slots of free surface and this is used in order to control the different levels of thermo-capillary convection. The dimensions of the slots were much smaller than the size of the molten zone and occupied one quarter of the periphery of the molten zone, Fig. 3. The temperature is fully controlled by a feedback system in the experiments so that the solidification interface moved up with a constant speed (Us ¼ 1:5 mm min1 ) during the cooling down stage, and passed each and every of the free slots, one by one. Therefore, accounting for solidification of a binary alloy, the movement of the solid– liquid interface, and the evolution of latent heat near the interface, the numerical method must consider time-dependent, three-dimensional flow in the entire zone. A new version of femLego [9], a three-dimensional parallel adaptive finite element method, therefore, was developed within the scope of this project. This tool was based on the automated code generation [10].
Fig. 3. Geometry of the floating zone configuration.
3.1. Governing equations In space, the only coupling between the energy equation and the momentum equations is through the boundary conditions at the free surface. Therefore, the governing transport dimensionless equations where the melt is treated as an incompressible Newtonian liquid read as r ðw~ u Þ ¼ 0;
ð1Þ
q~ u 1 n þ ð~ u rÞ~ u ¼ rp þ w~ u; r2 ~ uþ qt ReTC HðwÞ
ð2Þ
qy 1 1 qw þ ð~ u rÞy ¼ ; ðr2 yÞ þ qt ReTC Pr Ste qt
ð3Þ
where, ~ u is the fluid velocity, p is the pressure and y is the temperature. The extra term in the Navier– Stokes equation, Eq. (2), represents the porous medium model at the solid–liquid interface. w is the fraction of volume occupied by melt, w has a value of 0 in the solid region, 1 in the completely molten region. HðwÞ is the permeability of the mush [11]. HðwÞ is infinite in the molten region and very rapidly decreases to very small values (say, 108 ) in the solid region, so that the velocity ~ u become effectively zero in the solid region.
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
The dimensionless parameters occurring in Eqs. (1)–(3) are the thermo-capillary Reynolds number ReTC that based on the thermo-capillary velocity U 0 ¼ gDTm1 ; the Prandtl number Pr and the Stefan number. They are defined as gDTL ; ð4Þ ReTC ¼ mn n Pr ¼ ; a
ð5Þ
CpDT ; ð6Þ L where, g is a characteristic value of the coefficient of surface tension, m and n are dynamic and kinematic viscosities, a is the thermal diffusivity. In this study, the thermal conductivity k of melt and solid are the same and shown in Table 1. The fraction of volume that is occupied by solid can be computed, according to Amberg [11], qw ¼ k y yliq : ð7Þ qt Here, k is a numerical parameter and chosen k ¼ 2=Ste in order to satisfy the numerical stability problem. This equation allows a simple way of determination of the solid fraction w at the new time level by using the discretized form of the lefthand side: If w tends to increase beyond 1 or decrease below 0; the new value of w is taken to be 1 or 0; respectively. This means that, y is then allowed to differ from the nondimensionless of the Ste ¼
457
melting temperature yliq ; which is, of course, quite correct in the solid (w ¼ 0) or the liquid (w ¼ 1) region. The numerical parameter k is a small constant, used to amplify any deviation of the temperature from yliq and results in a rapid melting or freezing that restores y to yliq : The redistribution of dopant in the base material is obtained from the mass transfer equation. qcm 1 rcm : þ ð~ u rÞcl ¼ r ð8Þ Pe qt Here, cm is the net composition of the solid–liquid mixture, cl is the composition of the liquid. Pe ¼ U 0 L=D is the Peclet number with D is the molecular diffusivity of the dopant. This diffusion can be treated as partially solidified diffusion (Ds ) with a solute diffusion of the dopant (Dl ) as follows: D ¼ Ds ð1 wÞ þ Dl w:
ð9Þ
3.2. Boundary conditions Top surface. The temperature of the top surface is a function in time and follows the temperature history in the experiment. The dopant concentration at the top is a no-flux condition. qc ¼ 0: qz
ð10Þ
The velocities follow: uz ¼ 0
Table 1 Physical properties of Sn/Sb Physical property
Symbol
Value Sn/Sb
Density Liquidus temperature Kinetic viscosity Dynamic viscosity Thermal conductivity Specific heat Temp. coef. of surface tension Latent heat of fusion Solute diffusivity Solidified diffusivity Segregation coefficient of Sn in Sb
r ðkg m3 ) Tliq ðKÞ n ðm2 s1 Þ m ðN sm2 Þ k ðW m1 K1 Þ Cp ðJ kg1 K1 ) dg=dT ðN m1 KÞ
6:684 903:9 1:9 107 1:2 103 21:3 260 5:0 105
L ðJ kg1 Þ Dl ðm2 s1 Þ Ds ðm2 s1 Þ k0
1:63 105 2:0 109 5 1015 0:3
and
qux quy ¼ ¼ 0: qz qz
ð11Þ
Bottom surface. The temperature of the bottom surface is also a function in time and following the experiment temperature. The top and the bottom temperature are shown in Fig. 4. The bottom position shall never be melted. Therefore, the velocities can be zero. The dopant concentration satisfies n% rc ¼ 0:
ð12Þ
Free surface. The gradient of surface tension g in tangential direction i along the free surface in the slots must be balanced by shear stress, the latter being proportional to the normal flux of the
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
458
1100 1000 900 Temperature˚C
800
Top surface
700 600
Bottom surface
500 400 300
Heating up stage
200
Cooling down stage
100 0 0
100 200 300 400 500 600 700 800 time (s)
Fig. 4. The experiment temperature at the top and the bottom surface of GAS experiment.
tangential velocity ~ ui m
q~ u i qg qy : ¼ qn qy qi
ð13Þ
Here, m is the dynamic viscosity of the fluid and qg=qy is the coefficient of surface tension, n is the normal to the surface directed outward from the liquid region. The boundary condition for the dopant concentration as follows: n% rc ¼ 0:
ð14Þ
Coated surface. At the entirely coated surface the velocity obeys a no-slip boundary condition. n% rc ¼ 0;
~ u ¼ 0:
and Quanrtapelle [12]. The Poisson equation for the pressure is solved by using the well-known Conjugate-Gradient method. The convective term in Eqs. (2) and (3) and the dopant mass transfer equation (8) are calculated implicitly by using the GMRES method. Here, the streamline-upwind concept was used for the test functions in order to increase the stability [13]. In this study, the adaptive mesh was used in order to solve the diffusion problem at the solidification interface. The minimum of mesh size was limited at h ¼ 0:01 mm: Therefore, we would have at least 8 elements placed in the boundary layer due to diffusion, d0d ¼ 0:08 mm [14]. Generally, a mesh with about 1:2 million elements and 200 000 nodes was used in the present study. During each adaptive time step, there have been about 2000 elements that would be merged and 2000 new elements that would be created by the refinement stuff. All computations were performed on a PC-cluster with 16 CPUs, AMD Athlon MP 2400+ with 1GB memory for each. This cluster is a home made product of Mechanics Department of Royal Institute of Technology. The longest necessary computation for the whole cooling stage took about one week. A symbolic computational tool called femLego has been applied to create a parallel code that can solve the present time dependent, three-dimension mathematical problem.
ð15Þ
3.3. Numerical method In this study, a three-dimensional, parallel finite element method has been used in conjunction with an adaptive mesh refinement and coarsening. The equations have been discretised using a finite element approach on an unstructured grid. The type of elements which are used for both velocity and pressure variables are piecewise linear tetrahedral elements. The restrictions imposed by the Babuska–Brezzi condition are avoided by adding a pressure stabilisation term. The pressure and velocity solutions are split by using a fractional step algorithm that was introduced by Guermond
4. Results and discussion The problem studied in this paper is the numerical simulation of the space experiment (STS-108) that was performed in a Get Away Special (GAS) module in December 2001. Both experiment and numerical method focus on the radial segregation of tin in antimony samples due to weak convection. The physical properties of Sn/ Sb are listed in Table 1. The evolution of the axial dopant concentration distribution in a vertical slice at j ¼ 0 for case z ¼ 55% is shown in Fig. 5. The solidification interface was moving with a constant speed Usol E0:025 mm/s, exactly as it was programmed in the experiment [7]. That it to say, although the solidification interface was deformed
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
459
Fig. 5. Evolution of the concentration distributation at t ¼ 440 s (A), 580 s (B), 720 s (C) and 860 s (D). 55% free window.
Fig. 6. Velocity field at the vertical slice for the different free window, (A) z ¼ 10%; (B) z ¼ 25% and (C) z ¼ 55%:
by the convection, the speed of the interface was constant as expected. In this figure, local segregation close to the free windows appear as result of a complicated flow near the free surface. The other local segregation appears in the opposite side of the slots due to the recirculating fluid. This recirculating fluid together with the main flow field, swept dopant from the diffusion layer into the bulk liquid and induced radial segregation across the radial axis. That is the reason why the highest dopant concentration did not appear in the opposite side of the slots. They appeared instead in the middle of the sample (at r ¼ 0:5–2:5 mm).
A comparison of the flow field in a vertical slice for the different free window surface: z ¼ 10%; 25% and 55% is shown in Fig. 6. In this figure, the fluid flow is accelerated at the free window surface and decelerated at the coated surface. The overall flow fields are similar to the surface velocity solution for the mixed boundary, see [5]. In case of strong convection flow, the shape of the solidification interface is bent. Specially, at the early stage of cooling, when the solidification interface is placed around z ¼ 3 mm; the boundary at the partly coated side was not the well mixed region any more. It became the boundary layer of the ‘‘flat
ARTICLE IN PRESS 460
M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
Fig. 7. The concentration distributation at t ¼ 860 s with the different free window z ¼ 0% (A); z ¼ 10% (B); z ¼ 25% (C) and z ¼ 55% (D). 1.05 Effective distribution coefficient
plate’’ with an input flow from up to down. When the flow was as strong as in the case of 55% free window, Fig. 6C, the area that is close to the ‘‘flat plate’’ appears as a recirculating fluid. This is the reason why the solidification interface in Fig. 5A is more blend than it is in the cases (B)–(D). Fig. 7 shows a dopant field computed for the segregation coefficient (ko ¼ 0:3) of tin in antimony at different free window surface fraction: z ¼ 0%; 10%; 25% and 55%: Obviously, the maximum dopant concentration in the crystal was shifted from the slot side to the other side, as the convection increases. The radial segregation can be observed in Fig. 7. It is totally influenced by the convection that was controlled by the fraction of free window surface. A similar shift of maximum dopant concentration was also observed in the experimental results, Fig. 2. However, a direct quantitative agreement with the computed ditributions is not obtained. There are several possible reasons to weaknesses in the experimental results, e.g. the free surfaces at the slots might have some contamination influencing the surface tension coefficient. The measurements are also made with a method that cannot reveal all the small scale variations of concentration variations that the computations predict. The experimental results in Fig. 1, showing increased concentration along the axial direction can also qualitatively get some support from the calculated concentration
Numerical simulation Experiment
1
0.95
0.9
0.85
0.8 0
10
20
30
40
50
60
70
80
Free surface area [%]
Fig. 8. Comparison of the effective distribution coefficient between numerical and experiment data.
fields shown in Fig. 7, but also here the experimental data are more pronounced than the calculated. The effective distribution coefficients keff was approximated as follows [2]: keff ¼
cS ; cinit
ð16Þ
where cS was the average concentration across the crystal and cinit is the initial concentration of the material, in dimensionless form cinit ¼ 1: The numerical results of the effective distribution
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
461
Fig. 9. The dopant profile at z ¼ 10 mm; (z ¼ 0%, 10%, 25% and 55%).
coefficients together with the treatment data from experiment are shown in Fig. 8. A good agreement between the numerical results and the experimental results has been observed. Both results have keff ranging from 0:83 to 1: Such data indicates that a very low convection level occurred in the space samples and it is obviously that keff decreases with increasing free surface windows fraction.
Analyzing the numerical data in the same way as the experimental, a horizontal slice with 3 mm thickness is produced to obtain the radial segregation profiles. The results of four different z; the percentage of free window surfaces over the partially coated surface, are shown in Fig. 9. The maximum dopant concentration for the Sn-droped Sb looks like an island and shifts from the slot side
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463
462
1.05
Table 2 Comparison of the radial segregation due to the weak convection
1 0.95
Sn Conc.
0.9 0.85 0.8
zð%Þ
cmin
cmax
cavg
Dc0
10 25 55
0.8160 0.6399 0.5890
1.0146 1.0973 1.1830
0.9547 0.9470 0.8196
0.2080 0.4830 0.7248
0.75 0.7 0.65 0.6 0.55
00% free window 10% free window 25% free window 55% free window
-3
-2
-1
0 1 radius [mm]
2
3
Fig. 10. Variation of dopant concentration in the radial axis.
(j ¼ 0) to the entirely coated surface side (j ¼ p), as the convection increases. The islands appear as the levels of the thermo-capillary convection were not high enough to distribute the dopant over whole diameter of the sample. The maximum concentration island for 10% free window surface appears near the centre of slice cut and shifts far to the slot side in case of z ¼ 25%: In those two cases, the gradient in the dopant concentration is much sharper between the partially coated surface and the island of maximum concentration. It shows that, the influence of convective dopant transport is much more active, compressing the concentration lines. When the convection is strong enough, case z ¼ 55%; a recirculating flow appeared and drove the dopant from the entirely coated surface side to the centre. Therefore, the maximum dopant concentration was compressed in the radial axis direction and transported to the periphery at j ¼ p=2 and 3p=2: In this case the gradient in the dopant concentration was more or less the same between two sides of the island. The variation of dopant concentration along a diameter is shown in Fig. 10. The radial segregation defined as Dc0 ¼
cmax cmin ; cavg
ð17Þ
where, cmax ; cmin and cavg are maximum, minimum and average dopant concentration across the
horizontal slice of the sample. See Table 2. Obviously, the results listed in Table 2 show the influence of convection flow to the radial segregation. By increasing convection, the radial segregation will also be increased. At this point, the conclusion from the experimental result is different than the numerical result. In the experiment, cavg is dependent on the content of the starting material and is totally different in the various cases. Meanwhile, in the numerical method, the starting material remained constant for all cases.
5. Conclusions Melt growth experiments of semiconductor material using the gradient freeze technique were simulated. A corresponding, the space experiments were directly modelled. The moving solidification interface was captured using an adaptive mesh method with very fine resolution at the interface. When the convection was strong, the shape of the solidification interface was deformed and some recirculating flow regions appeared at the periphery. The radial segregation in the solidified material was obtained and a clear quantitative relation between the free surface area and the radial segregation was revealed. The dopant concentration results show a slight shift of the maximal concentration in the radial profiles.
References [1] S.A. Nikitin, V.L. Polezhayev, A.L. Fedyushkin, J. Crystal Growth 52 (1981) 471. [2] C.L. Chang, R.A. Brown, J. Crystal Growth 63 (1983) 343. [3] J.I.D. Alexander, J. Quazzini, F. Rosenberger, J. Crystal Growth 97 (1989) 285.
ARTICLE IN PRESS M. Do-Quang et al. / Journal of Crystal Growth 269 (2004) 454–463 [4] V.I. Polezhayev, S.A. Nikitin, Oxford, 1989, p. 237. [5] M. Levenstam, G. Amberg, E. Tillberg, T. Carlberg, J. Crystal Growth 104 (1990) 641. [6] T. Carlberg, in: International Symposium on Physical Sciences in Space, 2004, submitted for publication. [7] T. Carlberg, C. Winkler, G. Amberg, in: Proceedings of the First International Symposium on Microgravity Research & Applications in Physical Sciences & Biotechnology, Sorrento, 2000. . [8] P. Holm, K. Loth, B. Larsson, T. Carlberg, Proc. of the 16th ESA Symp. on European Rocket and Balloon
[9] [10] [11] [12] [13] [14]
463
Programmes and Related Research, St. Gallen, Switzerland, 2003, ESA SP-530, p. 113–117. M. Do-Quang, I. Loginova, G. Amberg, in preparation. . G. Amberg, R. Tonhardt, C. Winkler, Math. Comput. Simulation 49 (1999) 149. G. Amberg, Int. J. Heat Mass Transfer, 34 (1) (1991) 17. J.L. Guermond, L. Quanrtapelle, Int. J. Numer. Meth. Fluids. 26 (1998) 1039. O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, Vol. 3, Butterworth Heinemann, Oxford, 2000. C. Winkler, G. Amberg, T. Carlberg, J. Crystal Growth 210 (2000) 573.