- Email: [email protected]
melt interface during Czochralski silicon crystal growth
International Journal of Heat and Mass Transfer 142 (2019) 118463
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
The influence mechanism of melt flow instability on the temperature fluctuation on the crystal/melt interface during Czochralski silicon crystal growth Junling Ding, Lijun Liu ⇑ Key Laboratory of Thermo-Fluid Science and Engineering, Ministry of Education, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
a r t i c l e
i n f o
Article history: Received 4 April 2019 Received in revised form 26 June 2019 Accepted 20 July 2019
Keywords: Temperature fluctuation Flow instability Silicon Interface Crystal growth
a b s t r a c t The temperature fluctuation on the crystal/melt interface during Czochralski silicon crystal growth is firstly investigated with LES method based on the relationship between kinetic undercooling and growth rate. The effects of crystal and crucible rotation on the temperature fluctuation are analyzed. The results show that when crystal counter-rotates with crucible, the temperature fluctuations on the interface and that in the melt close to crystal are more intense compared with co-rotation. The characteristic frequencies of temperature fluctuations on the interface and that in the melt are identical for counter-rotation. However, this consistency of the characteristic frequencies is not obvious for co-rotation. In particular, the temperature fluctuations on the interface are about 0.4 K and 0.7 K for co-rotation and counterrotation, respectively. Importantly, the influence mechanism of melt flow instability on temperature fluctuation is clarified by the analyses of cross-correlation and local convective heat transport. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction The Czochralski (CZ) crystal growth method is the dominant technology for growing high quality bulk monocrystalline silicon for electronic applications. In the past decade, the melt volume and crucible diameter gradually increase for growing large-size crystal, reducing the cost of growing crystal. The diameter of silicon wafer has increased from 200 mm to 300 mm, or even the 450 mm [1,2]. Due to the unique crystal growth environment of CZ method, there are many forces in the melt, such as the thermal buoyancy and thermocapillary forces, centrifugal and Coriolis forces caused by the crystal and crucible rotations, the shear force influenced by the argon flow on the melt free surface [3,4]. With the increase of crucible size, these forces interact violently in the melt, causing the melt flow to turbulence characterized with velocity and temperature fluctuations [5]. There are several studies on the melt turbulence in a real ellipsoidal industrial crucible [5–8], mainly about the temperature fluctuation in the melt. Based on the experimental results with respect to 150 mm diameter crucible for growing 75 mm diameter crystal, Hoshikawa et al. [6] reported that the magnitude of temperature fluctuation at the central axis of melt region increases first and
⇑ Corresponding author. E-mail address: [email protected] (L. Liu). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118463 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.
then decreases with the distance from crystal/melt interface. The maximum temperature fluctuation is about 5.5 K at 20 mm from the interface. Vizman et al. [7] observed that the maximum temperature fluctuation near the crystal/melt interface is about 9 K for 104 mm crystal growth using 360 mm crucible. The similar results also have been found by Cen et al. [8] through threedimensional transient numerical simulations. Recently, Liu et al. [5] numerically investigated the instantaneous behaviors and statistical features of the melt turbulence in 750 mm crucible for growing 300 mm silicon crystal. The study showed that the maximal temperature fluctuation near the crystallization zone with an approximate value of 7 K. During the crystal growth, strong oscillations of temperature and melt flow in the crucible, especially the fluctuation of local temperature at the solidification front, directly leading the unstable temperature and temperature gradient on the crystal/melt interface. The unstable temperature at the interface causes the fluctuation of undercooling at the interface, which strongly influences the variation of crystal growth rate according to the kinetics of crystal growth [9]. Based on the so-called Voronkov criterion often used to evaluate dominant intrinsic point defect in the crystal [10,11], the variations of temperature gradient and growth rate on the crystal/melt interface are closely related to the formation of micro-defects in crystal. Therefore, it is important to study the characteristic of temperature fluctuation on the crystal/melt interface and investigate the influence mechanism of melt flow instabil-
2
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
ity on the interface temperature fluctuation. However, as mentioned above, all these works [5–8] only focused on the temperature fluctuation in melt, and to the authors’ knowledge, there are no relevant reports yet about the temperature fluctuation on the crystal/melt interface, which has a close correlation to the microdefect formation and impurity concentration in grown crystal [12]. In addition to the works mentioned above, most researchers in the crystal growth simulation consider that the crystal/melt interface is isothermal equal to the crystal melting temperature [13–16], which is the Dirichlet boundary condition and not suitable for solving temperature oscillation on the interface. To solve this problem, in our present study, the relationship between kinetic undercooling and growth rate is considered to develop a numerical model that can calculate the unstable melt flow and fluctuant temperature on the crystal/melt interface simultaneously. During the development of numerical model, the turbulent large eddy simulation (LES) technique is used to study the instability of melt flow. Based on the relationship between the release of latent heat and energy conservation at the crystal/ melt interface, the temperature fluctuation in melt region is coupled with the temperature distribution in crystal. Then, considering the relationship between undercooling and crystal growth rate, the temperature fluctuation on the crystal/melt interface is obtained. In addition, the rotations of crystal and crucible as the most easily controlled factors during CZ crystal growth, has been widely studied on its effects on the melt flow and temperature oscillation [17–19]. Therefore, the influences of crystal and crucible rotation on the temperature fluctuation in the melt and that on the crystal/melt interface are also considered. This paper is organized as follows: In Section 2, the computational model and corresponding grids are presented. Section 3 contains the numerical method including the LES governing equations and boundary conditions. The calculation of initial fields for LES and the independence verification of time step are shown in Section 4. Next, in Section 5, the results obtained in our work are discussed in detail. Finally, conclusions are presented. 2. Model description Fig. 1 shows the illustrations of the configuration and computational grids. For the simulation, the current popular 300 mm silicon wafer growing in 750 mm industrial crucible is considered. The depth of melt is 300 mm. And the crystal length is 600 mm. And the shape of meniscus is calculated according to the Ref. [20]. The crystal rotation rate xcry = 8 rpm, while two different crucible rotation rates xcru = 6 rpm and 6 rpm are performed. The
negative value indicates the clockwise rotation. Due to the geometric complexity of the crucible and crystal/melt interface, the computational grids are generated in body-fitted curvilinear coordinate system [21]. Fig. 1(b) shows the three-dimensional structured grids used for the computation of melt flow and heat transfer in the melt and crystal. The near-wall grid is refined for better solution of velocity and temperature fluctuation in the melt. 3. Numerical method 3.1. LES governing equations The melt flow and thermal fields in crystal and melt are calculated by solving the 3D time-dependent equations of continuity, momentum and energy conservation. The silicon melt is considered as incompressible Newtonian fluid with Boussinesq assumption. These equations are filtered based on the implicit top-hat filter in finite-volume solution methodology [22]. The filtered equations in Cartesian coordinates for the resolved scales of melt turbulence are written as follows:
@ ui ¼0 @xi
ð1Þ
q
@ ui @ ui @p @ 2 ui ¼ þ leff qgbT ðT T 0 Þ þ q uj @t @xj @xi @xj @xj ! qð2x ui þx ðx r ÞÞ
ð2Þ
@T @T @2 T ¼ aeff þ uj @t @xj @xj @xj
ð3Þ
where the ‘‘overline” denotes filtered variables. x is the angular velocity of reference frame, bT is the thermal expansion coefficient of silicon melt, T0 is the reference temperature and g is the gravitational acceleration. The effective eddy viscosity and effective thermal diffusivity are defined as leff = l + lSGS and aeff = a + lSGS/rSGS respectively. And rSGS is SGS turbulent Prandtl number with a value of 0.9. The dynamic eddy-viscosity model is given as:
lSGS ¼ qC S D2 S
ð4Þ
ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi @u where S ¼ 2 Sij Sij , and Sij ¼ 12 @@xuji þ @xij is the strain rate tensor.
D ¼ ðD1 D2 D3 Þ1=3 is the grid-filter width, and D1, D2, D3 are the grid
Fig. 1. The illustration of computational model. (a) Geometry and boundary conditions of the computational model (2D). (b) Sketch of partial computational grids (3D). Here, a few grids are displayed for better visibility.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
spacing in each space direction respectively. The model coefficient CS is a local time-dependent variable obtained by the approach proposed by Germano [23] and modified by Lily [24]. As known from Ghosal [25], the filtering operation and the differential operation do not commute in non-uniform grids, leading to the occurrence of numerical errors, which always exists when solving the complex geometry model. Owing to the mesh is uniform in computational space, we solve the governing equations in computational space to eliminate the numerical errors. And the SGS eddy viscosity lSGS is calculated with the dynamic SGS model for body-fitted grids described in our previous paper [5], which can satisfactorily predict turbulent fluctuation in the melt with much shorter CPU time and much smaller computer memory during the CZ crystal growth. Besides, in the regions near the crucible wall and crystal/melt interface, lSGS was estimated with a three-layer log-law wall function [26,27]. The LES governing equations are solved by using the finite volume method on non-orthogonal curvilinear grids. SIMPLEC algorithm is adopted to couple the velocity and pressure fields. The second-order central differencing scheme is applied to diffusion term, and QUICK scheme is used for convective term by deferred correction method [28]. For the time term, second-order implicit backward differencing scheme is adopted. The present numerical code for LES calculation has been validated by comparison with DNS results of Raufeisen et al. [29], which can be seen in our previous work [5]. The physical properties of silicon used in this work are provided in Table 1. 3.2. Boundary conditions No-slip velocity boundary conditions are applied at the crucible wall and crystal/melt interface. At the melt free surface, Marangoni convection is considered, which is induced by the surface tension caused by temperature gradients at the free surface. So, the shear stress is balanced with surface tension as follows [5]:
sgg ¼ l
@ug ¼0 @g
ð5Þ
sgn ¼ l
@un @r dr @T ¼ ¼ @g @n dT @n
ð6Þ
sgf ¼ l
@uf @r dr @T ¼ ¼ @g @f dT @f
ð7Þ
where n and f are local coordinates tangential to the melt surface and g is orthogonal to both n and f. The temperature coefficient of surface tension dr/dT is a constant (see Table 1). At the melt free surface and crystal surface, the heat loss is considered by Stefan–Boltzmann equation:
Table 1 The physical properties of silicon used in the work. Parameters
Symbols
Values
Density, crystal Density, melt Thermal conductivity, crystal Thermal conductivity, melt Heat capacity, crytsal Heat capacity, melt Dynamic viscosity Thermal expansion coefficient Latent heat of phase change Surface tension coefficient
qc qm
dr dT
2330 kg/m3 2530 kg/m3 22 W/(mK) 64 W/(mK) 1059 J/(kgK) 1000 J/(kgK) 7 104 kg/(ms) 1.4 104 K1 1.411 106 J/kg 1.0 104 N/(km)
Tm
0.58 0.3 1685 K
Emissivity factor, crystal Emissivity factor, melt Melting temperature
kc km Cp-c Cp-m
l
bT DH
ec em
keff
@T ¼ rSB eðT 4 T 4env Þ @n
3
ð8Þ
where keff is the effective thermal conductivity. rSB is the StefanBoltzmann constant. The ambient temperature Tenv of melt free surface is 1620 K, and the ambient temperature Tenv of crystal surface decreases linearly with the increase of crystal height z (see Fig. 1 (a)), from 1620 K to 500 K, which have been proved to be effective for simulating the crystal heat dissipation compared with global simulation [30]. According to the experimental data [31], a parabolic function relationship between the temperature profile Tcru at crucible wall and the crucible height is considered (see Fig. 1(a)). The minimum temperature is located on the bottom and the upper edge of the crucible. And the maximum temperature is located on the middle of the crucible wall. At the crystal/melt interface, the heat flux is balanced considering latent heat release, which satisfies the following equation:
@T @T k k ¼ qc DHV n @n c @n m
ð9Þ
where the subscript c and subscript m denote values of crystal and melt side. DH is the latent heat of crystallization. And Vn is the local normal growth rate, which is related to the undercooling on the interface. The undercooling on the interface is necessary for crystal growth, which is the difference between theoretical crystallization temperature and actual crystallization temperature. Fujiwara et al. [9] have deeply studied the relationship between growth rate and undercooling by conducting silicon growth experiments in recent years. They focus on the low undercooled melt, which is similar to the industrial growth process of bulk silicon crystal when using CZ or Floating zone (FZ) methods. They found the growth rate linearly increased with the undercooling degree:
V n ¼ 7:5 106 T m T f m=s
ð10Þ
where Tm is the melting temperature of silicon, and Tf is the local temperature on the crystal/melt interface. This relationship is adopted in this work. 4. Preparation for LES computation 4.1. Initial field and crystal/melt interface shape In the initial stage of LES calculation, reasonable temperature and velocity fields are necessary. In this work, the initial temperature and velocity fields used for LES are calculated by the RANS model presented in our previous work [13]. Besides, the crystal/ melt interface is a coupling boundary between melt and crystal. In order to obtain temperature fluctuation on the interface, different from previous studies [5–8], in the current study, temperature on the interface is not fixed, which is related to the heat flux near the interface and the release of latent heat. And the growth rate is obtained by Eq. (10), according to the local temperature Tf on the interface. The initial isothermal interface related to crystal growth rate is required for the reasonable calculation of temperature fluctuation on the interface by LES. Owing to the heat generated by heater (not shown in this work, see Ref. [13]) is transferred to melt through crucible wall, thus the crucible wall is considered as the heat source in present computational model. Therefore, the temperature Tcru on crucible wall need to be adjusted to ensure that the crystal/melt interface is isothermal. In the process of achieving an isothermal interface, crystal pulling velocity in the iso-diameter growth stage is selected as 0.8 mm/ min, thus temperature on the interface should be 1683.23 K based on the relationship between growth rate and undercooling (see
4
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
formula (10)). Besides, considering the corresponding temperature Tcru on crucible wall is not easy to known when the interface temperature is 1683.23 K, the standard form of PID controller is adopted to adjust the temperature Tcru conveniently, which as follows:
Z 1 t deðtÞ DT cru ðtÞ ¼ K p eðtÞ þ eðtÞdt þ sd dt si 0
ð11Þ
where DTcru(t)is the change of temperature on crucible wall, then the manipulated variable Tcru is adjusted to Tcru = Tcru(0) + DTcru(t) with Tcru(0) is previous temperature on crucible wall.
e(t) = Tref T tri (t) is difference between the reference temperature
Tref (1683.23 K) and current mean temperature T tri ðtÞ on the trijunction edge considering the circumferential temperature is not perfectly uniform in the 3D model. Kp is the gain, si and sd are the integral time and the derivate time, respectively. Due to the control process is a timeless problem, iteration step is used as the sampling interval in this calculation. As a starting point, we first present the distribution of monitoring points used in the calculation as shown in Fig. 2, as a basis for further analysis in subsequent sections. The situation of xcry = 8 rpm and xcru = 6 rpm is selected to illustrate the calculation process. The calculation algorithm involves two steps: (1) Using PID controller adjust the temperature on crucible wall
to make T tri close to Tref. This stage provides reasonable temperature and velocity fields for the next moving grid stage, ensuring the mesh movement can be carried out smoothly. Fig. 3 shows the convergence process of temperature at
monitoring points and T tri . It is obvious that the tempera-
tures in the melt and on the interface remain unchanged
approximately after 12,000 iteration steps, and T tri eventually stabilizes to the target temperature Tref. This illuminates that PID controller is an effective and convenient method for
controlling the crucible wall temperature to ensure T tri approach Tref. At the end of this stage, the crystal/melt interface is a flat non-isothermal interface as shown in Fig. 4. And the minimum temperature is located at the tri-junction edge. (2) Moving mesh to ensure the crystal/melt surface is isothermal and the temperatures in melt and crystal are converged.
During the mesh movement, T tri is also controlled by adjusting crucible wall temperature using PID controller. The location of the crystal/melt interface is determined in such a way that it is updated iteratively in the iterations according to
(Tf T tri ), until (Tf T tri ) approaches zero. As can be seen in
Fig. 5, with the grids moving, the temperature T tri and crucible wall temperature Tcru tend to stabilize gradually. After
two hundred mesh move steps, the temperature T tri reach to the reference temperature Tref 1683.23 K. In addition, the maximum temperature difference at the crystal/melt interface decrease with moving mesh. After three hundred move steps, the maximum temperature difference is less than 0.2 K. Thus, it is considered that the crystal/melt interface is approximate isothermal surface at this moment. At the end of this stage, the crystal/melt interface is an isothermal curved surface as shown in Fig. 6. 4.2. Selection of time step Considering the computational accuracy and numerical cost comprehensively, the grid size of 50 50 85 in the three coordinate directions of curvilinear coordinate system is used. Owing to the temperature fluctuation near the crystal/melt interface is the focus of our attention, the grids near the interface is specially refined. In order to obtain a solution independent of time step,
Fig. 2. Distribution of monitoring points on the crystal/melt interface. P1-P4 are close to the center region of interface, P5-P8 are near the edge. In addition to the above eight monitoring points, there are also eight monitoring points P1_B-P8_B in the melt located at depth of 20 mm from the interface corresponding to the eight points on the interface respectively, which are not shown for better visibility.
Fig. 3. Varieties of temperature at monitoring points with iteration steps.
Fig. 4. Temperature field in the melt at the end of step one.
Fig. 5. Process of obtaining the isothermal crystal/melt interface.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
5
Fig. 6. Temperature field in the melt at the end of step two. Fig. 8. Time-averaged fields of melt flow in the meridian plane when xcry = 8 rpm and xcru = 6 rpm. The left part: temperature distribution. The interval is 0.5 K. The right part: stream function distribution. The interval is 0.05 kg/s.
Fig. 9. Time-averaged fields of melt flow in the meridian plane when xcry = 8 rpm and xcru = 6 rpm. The left part: temperature distribution. The interval is 0.5 K. The right part: stream function distribution. The interval is 0.05 kg/s.
Fig. 7. Spectral analyses of temperature fluctuation at point P2_B for different time steps. (a) 0.02 s time step. (b) 0.01 s time step. (c) 0.005 s time step.
the effects of time steps on the characteristics of temperature fluctuation near interface are investigated. Fig. 7 shows the power spectra density (PSD) distributions of temperature fluctuation at point P2_B in melt for different time steps 0.02 s, 0.01 s, 0.005 s, respectively. A total of 60,000 values in a time range of 600 s are used for the spectral analysis by fast Fourier transformed method, which is also used in Section 5. Obviously, when the time steps are selected as 0.01 s and 0.005 s, there is an identical characteristic frequency of the temperature fluctuation. However, the identical characteristic frequency disappears with time step 0.02 s, indicating that it is not suitable for the analysis of characteristics of temperature fluctuation in the melt. Considering the computing efficiency, thus the time step 0.01 s is selected in the subsequent calculations. 5. Results and discussion 5.1. Time-averaged field of melt flow To analyze the statistical behavior of the melt flow, the velocity and temperature fields are averaged in time and in the circumferential direction. Figs. 8 and 9 display the distributions of statistically averaged temperature and stream function for different crystal and crucible rotation conditions. Both for counter-rotation and co-rotation, it is obvious that three convection rolls exist in the melt, the Taylor-Proudman vortex, the buoyant vortex and a secondary vortex. The counterclockwise buoyant vortex is driven by thermal buoyancy, which is induced by the axial temperature
difference in the melt, holds the great area of melt region. The buoyant melt flow upward along the crucible wall, then reach the melt free surface, causing the clockwise secondary vortex around the corner. In the center melt region below the crystal, the buoyant vortex is damped by the clockwise Taylor-Proudman vortex, which is induced by the rotation of crystal and crucible. The difference is that when the crystal counter-rotates with crucible, the Taylor-Proudman vortex is larger than that in corotation. The larger Taylor-Proudman vortex, the more complex melt flow under the crystal. With regard to the time-averaged temperature, it is found that the trend of temperature distribution in the melt remains unchanged for different crystal and crucible rotation conditions. However, by comparing the distribution of isotherms, the melt temperature for counter-rotation is higher than that for co-rotation, especially near the crystal/melt interface. This indicates that more heat is transferred to the crystal/melt interface when crystal counter-rotates with crucible. 5.2. Characteristic of temperature fluctuation In this section, we firstly study the characteristics of temperature fluctuation on the crystal/melt interface and that in the melt near interface. In order to display the information of temperature fluctuation more clearly, temperature fluctuations at points P1P8 are presented by boxplots as shown in Fig. 10. The boxplot is a flexible exploratory-data-analysis tool for graphically representing the key values from numerical data in detail. From the boxplot, information about the location, spread, skewness, and longtailedness of observations can be quickly known [32]. In the paper, the whiskers are determined by the maximum and minimum temperatures at monitoring points. The square in the box represents the mean value of temperature fluctuations. The horizontal line in the box expresses the median of temperature data. The upper edge of the box represents the 75th percentile, and the lower edge is the 25th percentile. For details of the box diagram, one can see Ref. [32]. In summary, the height of the box or the distance between
6
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Fig. 10. Temperature fluctuations at monitoring points P1-P8 on the interface for different crystal and crucible rotation conditions. (a) xcry = 8 rpm, xcru = 6 rpm. (b) xcry = 8 rpm, xcru = 6 rpm.
the ends of whiskers represents the range of temperature fluctuation approximatively. From Fig. 10, it can be noticed that when crystal counter-rotates with crucible, both the height of the box and the distance between the ends of whiskers are longer compared to that when crystal and crucible rotate in the same direction. This observation illustrates that temperature oscillation on the crystal/melt interface becomes more intense when crystal and crucible rotate in the opposite direction. Besides, the changes of cut line at the median which represents location of box are smaller when crystal counter rotate with crucible, especially in the central region of interface. This reveals that temperature fluctuation on the interface is more uniform for counter-rotation. Considering Fig. 10(a) and (b) synthetically, based on the height of box, the conservative temperature fluctuation on the crystal/melt interface is about 0.4 K for corotation and 0.7 K for counter-rotation. Fig. 11 shows the temperature fluctuations at points P1_B-P8_B in the melt. Similarly, the temperature oscillation in the melt near interface is more intense for counter-rotation compared with that obtained for co-rotation. Combining the results of time-averaged melt flow in Section 5.1, the probable explanation for this observation is as follows: when the crystal counter-rotates with the crucible, the Taylor-Proudman vortex always appears in the region below the crystal/melt interface (see Fig. 5 in Ref. [33]). In the Taylor-Proudman vortex, the melt near crucible wall spirals up to the region near interface, resulting in the intense mixing of melt below crystal, then giving rise to the larger temperature fluctuation near interface. However, the Taylor column vanishes when the rotation rate of crystal and crucible is same both in magnitude and direction [33]. This explanation is also proved by the analysis
of local heat transfer in the melt in Section 5.4. Carefully, the rotation of crucible is slightly higher than that of crystal in our present simulation. Thus, a small Taylor-Proudman vortex exist in our simulation, as shown in Fig. 9, which also been observed in Ref. [33]. From Fig. 11, it can also be found that the ranges of temperature fluctuation at all monitoring points are approximately identical for counter-rotation of crystal and crucible, while the temperature fluctuation in the central region is weaker than that near the crystal edge for co-rotation. Comprehensively, based on the height of box, the temperature fluctuations are about 2.8 K and 2.0 K for counter-rotation and co-rotation, respectively. In Refs. [7,31], experiments about the observation of temperature fluctuation near interface for counter-rotation of crystal and crucible have been performed. The difference between the maximum and minimum temperatures is 9.2 K in Ref. [7], while it is 10 K in Ref. [31]. According to the distance between the ends of whiskers in Fig. 11(a), the range of temperature fluctuation in the melt below crystal in our numerical results are in good agreement with those obtained by experiments. 5.3. Relation between temperature fluctuation on the interface and that in the melt The relationship between the temperature fluctuation on the interface and that in the melt close to crystal is analyzed in this section. The cross correlation as a general method of estimating the relevance of two time series is adopted in the work. Fig. 12 illustrates the maximum cross-correlation coefficients of temperature fluctuations at monitoring points P1-P8 distributed on the interface and points P1_B-P8_B in the melt for different crystal
Fig. 11. Temperature fluctuations at monitoring points P1_B-P8_B in the melt for different crystal and crucible rotation conditions. (a) xcry = 8 rpm, xcru = 6 rpm. (b) xcry = 8 rpm, xcru = 6 rpm.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
7
Fig. 13. The cross-correlations of temperature fluctuations at monitoring points P2 and P6 on the interface with their corresponding points P2_B and P6_B in the melt respectively when xcry = 8 rpm, xcru = 6 rpm. (a) P2 and P2_B. (b) P6 and P6_B.
Fig. 12. The maximum cross-correlation coefficients of temperature fluctuations on the interface with that in the melt for different crystal and crucible rotation conditions. (a) xcry = 8 rpm, xcru = 6 rpm. (b) xcry = 8 rpm, xcru = 6 rpm.
and crucible rotation conditions. It is obvious that there is a strong relationship between the temperature variation at the points distributed on the interface and that just below these points. And the correlation is slightly enhanced for counter-rotation compared with co-rotation. In the calculation process, a similar characteristics of temperature fluctuations at points in the center region of interface and melt close to crystal has been found. And the similarity also exists for the points near the edge. Thus, for the sake of clarity, the points P2 (P2_B) and P6 (P6_B) are selected to represent the center region and edge area of the interface (melt) in the latter discussion, respectively. Fig. 13 illustrates the cross correlation series with a maximum time delay of 600 s when the crystal counter-rotates with crucible in detail. There is a strong correlation at the time delay of about 0.61 s at points P2 and its corresponding point P2_B. And the maximum correlation is achieved at the time delay of 0.54 s at points P6 and P6_B. Besides, when crystal rotate with crucible in the same direction, the time delays related to the maximum correlation are 0.41 s and 0.2 s for P2&P2_B and P6&P6_B respectively, as shown in Fig. 14. The positive value of time delay represents that the temperature variation on the crystal/melt interface lags behind that in the melt. According to our results, it is found that the values of time delays at all points P1-P8 with their corresponding points P1_B-P8_B for counter-rotation and that at the points P1-P4 with their corresponding points P1_B-P4_B for co-rotation are positive. On the contrary, the values of time delays are negative in the region close to fringe of crystal when crystal co-rotates with crucible. The reasons for this will be explained in the next section. The frequency analyses of temperature fluctuations on the crystal/melt interface and in the melt are performed, in order to obtain the characteristics of temperature fluctuation, as shown in Figs. 15 and 16. For the convenience of comparison, the PSD is normalized by its maximum value, which is denoted as PSDN. It can be seen
Fig. 14. The cross-correlations of temperature fluctuations at monitoring points P2 and P6 on the interface with their corresponding points P2_B and P6_B in the melt respectively when xcry = 8 rpm, xcru = 6 rpm. (a) P2 and P2_B. (b) P6 and P6_B.
that there are good agreements for dominant frequencies at P2 and P6 with their corresponding points in the melt when crystal counter-rotates with crucible. However, the consistency of dominant frequencies is not obvious for co-rotation. And no frequencies consistent with the crystal or crucible rotational frequency are found in two figures, which are 0.133 Hz and 0.1 Hz respectively. For the counter-rotation, the dominant frequency at P2 is 0.0136 Hz, one-tenth of crystal rotational frequency, which is probably influenced by crystal rotation. Different from the corotation, PSDN is strong in the lower frequency range for counterrotation, indicating that larger periodic flows occur in the melt below the crystal/melt interface when crystal counter-rotates with crucible. The reason for this is that when crystal counter-rotates with crucible, the melt flow in the crucible is dominated by the buoyancy, which is a large periodic circulation passing crucible wall, melt free surface and bottom of crystal successively [5,8]. While when crystal co-rotates with crucible, the melt flow is mainly affected by rotations of crystal and crucible, which is characterized by small period [33].
8
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Fig. 15. Frequency analyses of temperature fluctuations at points P2, P6 on the interface and the corresponding points P2_B, P6_B in the melt when xcry = 8 rpm, xcru = 6 rpm.
Fig. 16. Frequency analyses of temperature fluctuations at points P2, P6 on the interface and the corresponding points P2_B, P6_B in the melt when xcry = 8 rpm, xcru = 6 rpm.
5.4. Analysis of heat transfer in the melt The instantaneous local convective heat flux proposed by Shang et al. [34,35] is adopted to study the heat transfer in the melt near crystal/melt interface, according to local velocity U(r,t) and local temperature T(r,t). The definition of instantaneous local convective heat flux is as follows:
jðr; tÞ ¼
C pm qm U ðr; t Þ½T ðr; tÞ T 0 km DT=L
ð12Þ
where T0 is the mean temperature of bulk melt, DT is the temperature difference between crucible and crystal/melt interface, L is the depth of melt. The local convective heat flux is normalized by conductive heat flux kmDT/L. Besides, Shang et al. proved that the diffusive contribution can be ignored compared to the convective heat flux [34]. According to the Eq. (12), it is known that the warm fluctuation (T(r, t) T0) > 0) produces positive flux when they move
upward to the crystal/melt interface (UZ > 0), and the falling (UZ < 0) cold plume (T(r, t) T0) < 0) also contributes to a positive flux. The positive flux at point P2_B or P6_B in the melt means that the heat is transferred to crystal/melt interface.
5.4.1. Crystal counter-rotates with crucible Fig. 17 shows the time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuations at points P2_B and P6_B in the melt close to crystal. At point P2_B, as shown in Fig. 17(a), the fluctuation of circumferential velocity is slightly skewed toward the positive direction, which is directly influenced by rotations of crystal and crucible. At point P6_B, the region close to fringe of crystal, the skew is further strengthened as shown in Fig. 17(b). While, the radial and axial velocity fluctuate around zero value both at points P2_B and P6_B. It is inferred that the effects of crystal and crucible rotation on melt flow in the radial and axial directions are small at the loca-
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
9
Fig. 17. Time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuation in the melt when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
tion 2 cm from interface in our present work, and the melt has a free movement in both directions. Fig. 18 displays the histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions at points P2_B and P6_B, respectively. To display the histograms of j at different directions in a same graph, the histograms are normalized by their maximum value H0. When the thermal plumes exists, the temperature fluctuation is related to the velocity fluctuation, producing positive fluctuations [34]. In Fig. 18(a), it is seen that the H(jUT)/H0 and H(jUR)/H0 are symmetric and mean values of jUT and jUR are zero, indicating that no such correlation exists between the temperature and circumferential, radial velocity fluctuations. Duo to the falling cold plums and rising warm plumes produce positive fluxes in the vertical direction as mentioned above, the H(jUZ)/H0 becomes highly skewed toward the positive direction, even though velocity fluctuations are symmetric as shown in Fig. 17(a). The similar principle can be used to explain the symmetry of circumferential heat flux, although the circumferential velocity is slightly skewed toward the positive direction. Besides, the same shape of H(jUT)/H0 and H(jUR)/H0 suggests that there are same statistical characteristics of instantaneous local heat flux in the two directions. In Fig. 18(b), for point P6_B, the region close to fringe of crystal, because of the obvious circumferential motion, the circumferential heat flux becomes highly skewed toward the positive
direction with a long decay tail, indicating that the heat flux fluctuations produced by thermal plumes contain some large but rare events [34]. Besides, the radial heat flux is slightly skewed toward the negative direction, approximately same as the negative flux fluctuations in circumferential direction. Similar to P2_B, H(jUZ)/ H0 of P6_B also becomes skewed toward the positive direction with the same shape of P2_B. In summary, when crystal counter-rotates with crucible, the heat transfer is dominated by axial local heat flux in the center region of melt near crystal. While in the region close to the fringe of crystal, the heat transfer is complex with the asymmetries of local heat flux in three directions.
5.4.2. Crystal co-rotates with crucible Fig. 19 shows the time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuations at points P2_B and P6_B in the melt close to crystal when crystal co-rotates with crucible. Similar to the situation for counter-rotation, the fluctuations of radial and axial velocity components are symmetric with zero mean values both at P2_B and P6_B. The difference is that the fluctuation of circumferential velocity is skewed toward the negative direction when crystal co-rotates with crucible. From Figs. 17 and 19, it is found that the direction of circumferential velocity at the location 2 cm from
Fig. 18. Histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
10
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Fig. 19. Time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuation in the melt when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
interface is the same as that of crucible rotation, indicating that the melt flow in crucible is mainly influenced by crucible rotation. Fig. 20 illustrates the histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions at P2_B and P6_B, respectively. Similar to the situation for counter-rotation, at P2_B, the fluctuations of the circumferential and radial heat flux are symmetric with zero mean values, and the H(jUZ)/H0 also becomes slightly skewed toward the positive direction. By comparing Figs. 18(a) and 20(a), it is clear that the values of instantaneous local heat flux j at P2_B for counter-rotation are larger than that for co-rotation, which is the reason for the occurrence of more violent temperature oscillations when crystal counter-rotates with crucible, especially in the central area of the interface, as displayed in Fig. 11. For point P6_B, duo to the reversal of the direction of crucible rotation for co-rotation compared with counter-rotation, the H(jUT)/H0 becomes highly skewed toward the negative direction. And the H(jUR)/H0 is skewed toward the negative direction. It is important to note that the axial heat flux is also slightly skewed toward the negative direction, which is clearly different from the situation in P2_B for co-rotation and that in P2_B and P6_B for counter-rotation. This is the real cause of the only negative delay time of P6_B when crystal co-rotates with crucible, while which is positive in other three conditions.
6. Conclusions In this paper, temperature fluctuations on the crystal/melt interface during CZ silicon crystal growth are investigated by LES with the curvilinear dynamic Smagorinsky subgrid-scale model, considering the relationship between kinetic undercooling and growth rate. And the influence mechanism of melt flow instability near interface on the temperature fluctuations are clarified in detail. The following conclusions are obtained from this study: (1) When crystal counter-rotates with crucible, temperature fluctuations on the crystal/melt interface and that in the melt close to crystal are more intense compared with corotation. And the temperature fluctuations on the crystal/ melt interface are about 0.4 K for co-rotation and 0.7 K for counter-rotation, respectively. (2) The temperature variation at points distributed on the interface has a strong relationship with that just below these points. And the temperature variation on the crystal/melt interface lags behind that in the melt for the counterrotation and in the central region when co-rotation occurs. However, the situation is exactly opposite in the region close to fringe of crystal when crystal co-rotates with crucible.
Fig. 20. Histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Besides, the characteristic frequencies of temperature fluctuations on the interface and that in the melt are identical for counter-rotation. However, this consistency of dominant frequencies is not obvious for co-rotation. (3) For both counter-rotation and co-rotation, the heat transfer is dominated by axial local heat flux in the center region of melt close to crystal. The difference is that larger local heat flux exists for counter-rotation compared with co-rotation. And in the region close to fringe of crystal, the heat transfer is complex with the asymmetries of local heat flux in three directions. Besides, the axial heat transfer in the region close to the fringe of crystal for co-rotation is significantly different from the situations of counter-rotation and that in the central region for co-rotation. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 51676154). References [1] Z. Lu, S. Kimbel, Growth of 450mm diameter semiconductor grade silicon crystals, J. Cryst. Growth 318 (2001) 193–195. [2] E. Kamiyama, J. Vanhellemont, K. Sueoka, K. Araki, K. Izunome, Thermal stress induced void formation during 450 mm defect free silicon crystal growth and implications for wafer inspection, Appl. Phys. Lett. 102 (2013) 082108. [3] T. Shen, C.M. Wu, Y.R. Li, Experimental investigation on the effect of crystal and crucible rotation on thermocapillary convection in a Czochralski configuration, Int. J. Therm. Sci. 104 (2016) 20–28. [4] A. Muiznieks, A. Krauze, B. Nacke, Convective phenomena in large melts including magnetic fields, J. Cryst. Growth 303 (2007) 211–220. [5] L.J. Liu, X. Liu, Y. Wang, Large-eddy simulation of melt turbulence in a 300-mm Cz–Si crystal growth, Int. J. Heat Mass Transf. 55 (2012) 53–60. [6] H. Hirata, K. Hoshikawa, Three-dimensional numerical analyses of the effects of a cusp magnetic field on the flows, oxygen transport and heat transfer in a Czochralski silicon melt, J. Cryst. Growth 125 (1992) 181–207. [7] D. Vizman, O. Grabner, G. Muller, 3D numerical simulation and experimental investigations of melt flow in an Si Czochralski melt under the influence of a cusp-magnetic field, J. Cryst. Growth 236 (2002) 545–550. [8] X.R. Cen, Y.S. Li, J.M. Zhan, Three dimensional simulation of melt flow in Czochralski crystal growth with steady magnetic fields, J. Cryst. Growth 340 (2012) 135–141. [9] K. Fujiwara, K. Maeda, N. Usami, G. Sazaki, Y. Nose, A. Nomura, T. Shishido, K. Nakajima, In situ observation of Si faceted dendrite growth from low-degreeof-undercooling melts, Acta Mater. 56 (2008) 2663–2668. [10] V.V. Voronkov, The mechanism of swirl defects formation in silicon, J. Cryst. Growth 59 (1982) 625–643. [11] J. Vanhellemont, Intrinsic point defect incorporation in silicon single crystals grown from a melt, revisited, J. Appl. Phys. 110 (2011) 063519. [12] K. Kakimoto, T. Shyo, M. Eguchi, Correlation between temperature and impurity concentration fluctuations in silicon crystals grown by the Czochralski method, J. Cryst. Growth 151 (1995) 187–191.
11
[13] L.J. Liu, K. Kakimoto, Partly three-dimensional global modeling of a silicon Czochralski furnace. I. Principles, formulation and implementation of the model, Int. J. Heat Mass Transf. 48 (2005) 4481–4491. [14] A. Chatterjee, V. Prasad, A full 3-dimensional adaptive finite volume scheme for transport and phase-change processes, part I: formulation and validation, Numer. Heat Tranf. A-Appl. 37 (2000) 801–821. [15] L.J. Liu, Q.H. Yu, X.F. Qi, W.H. Zhao, G.X. Zhong, Controlling solidification front shape and thermal stress in growing quasi-single-crystal silicon ingots: process design for seeded directional solidification, Appl. Therm. Eng. 91 (2015) 225–233. [16] V. Kumar, G. Biswas, G. Brenner, F. Durst, Effect of thermocapillary convection in an industrial Czochralski crucible: numerical simulation, Int. J. Heat Mass Transf. 46 (2003) 1641–1652. [17] S.S. Son, K.W. Yi, Characteristics of thermal fluctuation in a low Pr number melt at a large crucible for Czochralski crystal growth method, J. Cryst. Growth 275 (2005) E259–E264. [18] C. Wagner, R. Friedrich, Direct numerical simulation of momentum and heat transport in idealized Czochralski crystal growth configurations, Int. J. Heat Fluid Flow 25 (2004) 431–443. [19] T. Shen, C.M. Wu, L. Zhang, Y.R. Li, Experimental investigation on effects of crystal and crucible rotation on thermal convection in a model Czochralski configuration, J. Cryst. Growth 438 (2016) 55–62. [20] T. Duffar, Crystal Growth Processes Based on Capillarity: Czochralski, Floating Zone, Shaping and Crucible Techniques, Wiley, Wiltshire, 2010, pp. 479–481. [21] J.F. Thompson, Z.U.A. Warsi, C.W. Mastin, Boundary-fitted coordinate systems for numerical solution of partial differential equations-a review, J. Comput. Phys. 47 (1982) 1–108. [22] E.M.J. Komen, L.H. Camilo, A. Shams, B.J. Geurts, B. Koren, A quantification method for numerical dissipation in quasi-DNS and under-resolved DNS, and effects of numerical dissipation in quasi-DNS and under-resolved DNS of turbulent channel flows, J. Comput. Phys. 345 (2017) 565–595. [23] M. Germano, U. Piomelli, P. Moin, W. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A 3 (1991) 1760–1765. [24] D.K. Lilly, A proposed modification of the Germano subgrid-scale closure method, Phys. Fluids A 4 (1992) 633–635. [25] S. Ghosal, P. Moin, The basic equations for the large eddy simulation of turbulent flows in complex geometry, J. Comput. Phys. 118 (1995) 24–37. [26] I.Y. Evstratov, V.V. Kalaev, A.I. Zhmakin, Y.N. Makarov, A.G. Abramov, N.G. Ivanov, E.M. Smirnov, E. Dornberger, J. Virbulis, E. Tomzig, W.V. Ammon, Modeling analysis of unsteady three-dimensional turbulent melt flow during Czochralski growth of Si crystals, J. Cryst. Growth 230 (2001) 22–29. [27] L. Temmerman, M.A. Leschziner, C.P. Mellen, J. Frohlich, Investigation of wallfunction approximations and subgrid-scale models in large eddy simulation of separated flow in a channel with streamwise periodic constrictions, Int. J. Heat Fluid Flow 24 (2003) 157–180. [28] T. Hayase, J.A.C. Humphrey, R. Greif, A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures, J. Comput. Phys. 98 (1992) 108–118. [29] A. Raufeisen, M. Breuer, T. Botsch, A. Delgado, DNS of rotating buoyancy- and surface tension-driven flow, Int. J. Heat Mass Transf. 51 (2008) 6219–6234. [30] J.L. Ding, L.J. Liu, Real-time prediction of crystal/melt interface shape during Czochralski crystal growth, Crystengcomm 20 (2018) 6925–6931. [31] S. Enger, O. Grabner, G. Muller, M. Breuer, F. Durst, Comparison of measurements and numerical simulations of melt convection in Czochralski crystal growth of silicon, J. Cryst. Growth 230 (2001) 135–142. [32] Y. Benjamini, Opening the box of a boxplot, Am. Stat. 42 (1988) 257–262. [33] J.R. Carruthers, K. Nassau, Nonmixing cells due to crucible rotation during Czochralski crystal growth, J. Appl. Phys. 39 (1968) 5205–5214. [34] X.D. Shang, X.L. Qiu, P. Tong, K.Q. Xia, Measurements of the local convective heat flux in turbulent Rayleigh-Benard convection, Phys. Rev. E 70 (2004) 026308. [35] X.D. Shang, P. Tong, K.Q. Xia, Scaling of the local convective heat flux in turbulent Rayleigh-Benard convection, Phys. Rev. Lett. 100 (2008) 244503.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
The influence mechanism of melt flow instability on the temperature fluctuation on the crystal/melt interface during Czochralski silicon crystal growth Junling Ding, Lijun Liu ⇑ Key Laboratory of Thermo-Fluid Science and Engineering, Ministry of Education, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
a r t i c l e
i n f o
Article history: Received 4 April 2019 Received in revised form 26 June 2019 Accepted 20 July 2019
Keywords: Temperature fluctuation Flow instability Silicon Interface Crystal growth
a b s t r a c t The temperature fluctuation on the crystal/melt interface during Czochralski silicon crystal growth is firstly investigated with LES method based on the relationship between kinetic undercooling and growth rate. The effects of crystal and crucible rotation on the temperature fluctuation are analyzed. The results show that when crystal counter-rotates with crucible, the temperature fluctuations on the interface and that in the melt close to crystal are more intense compared with co-rotation. The characteristic frequencies of temperature fluctuations on the interface and that in the melt are identical for counter-rotation. However, this consistency of the characteristic frequencies is not obvious for co-rotation. In particular, the temperature fluctuations on the interface are about 0.4 K and 0.7 K for co-rotation and counterrotation, respectively. Importantly, the influence mechanism of melt flow instability on temperature fluctuation is clarified by the analyses of cross-correlation and local convective heat transport. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction The Czochralski (CZ) crystal growth method is the dominant technology for growing high quality bulk monocrystalline silicon for electronic applications. In the past decade, the melt volume and crucible diameter gradually increase for growing large-size crystal, reducing the cost of growing crystal. The diameter of silicon wafer has increased from 200 mm to 300 mm, or even the 450 mm [1,2]. Due to the unique crystal growth environment of CZ method, there are many forces in the melt, such as the thermal buoyancy and thermocapillary forces, centrifugal and Coriolis forces caused by the crystal and crucible rotations, the shear force influenced by the argon flow on the melt free surface [3,4]. With the increase of crucible size, these forces interact violently in the melt, causing the melt flow to turbulence characterized with velocity and temperature fluctuations [5]. There are several studies on the melt turbulence in a real ellipsoidal industrial crucible [5–8], mainly about the temperature fluctuation in the melt. Based on the experimental results with respect to 150 mm diameter crucible for growing 75 mm diameter crystal, Hoshikawa et al. [6] reported that the magnitude of temperature fluctuation at the central axis of melt region increases first and
⇑ Corresponding author. E-mail address: [email protected] (L. Liu). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118463 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.
then decreases with the distance from crystal/melt interface. The maximum temperature fluctuation is about 5.5 K at 20 mm from the interface. Vizman et al. [7] observed that the maximum temperature fluctuation near the crystal/melt interface is about 9 K for 104 mm crystal growth using 360 mm crucible. The similar results also have been found by Cen et al. [8] through threedimensional transient numerical simulations. Recently, Liu et al. [5] numerically investigated the instantaneous behaviors and statistical features of the melt turbulence in 750 mm crucible for growing 300 mm silicon crystal. The study showed that the maximal temperature fluctuation near the crystallization zone with an approximate value of 7 K. During the crystal growth, strong oscillations of temperature and melt flow in the crucible, especially the fluctuation of local temperature at the solidification front, directly leading the unstable temperature and temperature gradient on the crystal/melt interface. The unstable temperature at the interface causes the fluctuation of undercooling at the interface, which strongly influences the variation of crystal growth rate according to the kinetics of crystal growth [9]. Based on the so-called Voronkov criterion often used to evaluate dominant intrinsic point defect in the crystal [10,11], the variations of temperature gradient and growth rate on the crystal/melt interface are closely related to the formation of micro-defects in crystal. Therefore, it is important to study the characteristic of temperature fluctuation on the crystal/melt interface and investigate the influence mechanism of melt flow instabil-
2
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
ity on the interface temperature fluctuation. However, as mentioned above, all these works [5–8] only focused on the temperature fluctuation in melt, and to the authors’ knowledge, there are no relevant reports yet about the temperature fluctuation on the crystal/melt interface, which has a close correlation to the microdefect formation and impurity concentration in grown crystal [12]. In addition to the works mentioned above, most researchers in the crystal growth simulation consider that the crystal/melt interface is isothermal equal to the crystal melting temperature [13–16], which is the Dirichlet boundary condition and not suitable for solving temperature oscillation on the interface. To solve this problem, in our present study, the relationship between kinetic undercooling and growth rate is considered to develop a numerical model that can calculate the unstable melt flow and fluctuant temperature on the crystal/melt interface simultaneously. During the development of numerical model, the turbulent large eddy simulation (LES) technique is used to study the instability of melt flow. Based on the relationship between the release of latent heat and energy conservation at the crystal/ melt interface, the temperature fluctuation in melt region is coupled with the temperature distribution in crystal. Then, considering the relationship between undercooling and crystal growth rate, the temperature fluctuation on the crystal/melt interface is obtained. In addition, the rotations of crystal and crucible as the most easily controlled factors during CZ crystal growth, has been widely studied on its effects on the melt flow and temperature oscillation [17–19]. Therefore, the influences of crystal and crucible rotation on the temperature fluctuation in the melt and that on the crystal/melt interface are also considered. This paper is organized as follows: In Section 2, the computational model and corresponding grids are presented. Section 3 contains the numerical method including the LES governing equations and boundary conditions. The calculation of initial fields for LES and the independence verification of time step are shown in Section 4. Next, in Section 5, the results obtained in our work are discussed in detail. Finally, conclusions are presented. 2. Model description Fig. 1 shows the illustrations of the configuration and computational grids. For the simulation, the current popular 300 mm silicon wafer growing in 750 mm industrial crucible is considered. The depth of melt is 300 mm. And the crystal length is 600 mm. And the shape of meniscus is calculated according to the Ref. [20]. The crystal rotation rate xcry = 8 rpm, while two different crucible rotation rates xcru = 6 rpm and 6 rpm are performed. The
negative value indicates the clockwise rotation. Due to the geometric complexity of the crucible and crystal/melt interface, the computational grids are generated in body-fitted curvilinear coordinate system [21]. Fig. 1(b) shows the three-dimensional structured grids used for the computation of melt flow and heat transfer in the melt and crystal. The near-wall grid is refined for better solution of velocity and temperature fluctuation in the melt. 3. Numerical method 3.1. LES governing equations The melt flow and thermal fields in crystal and melt are calculated by solving the 3D time-dependent equations of continuity, momentum and energy conservation. The silicon melt is considered as incompressible Newtonian fluid with Boussinesq assumption. These equations are filtered based on the implicit top-hat filter in finite-volume solution methodology [22]. The filtered equations in Cartesian coordinates for the resolved scales of melt turbulence are written as follows:
@ ui ¼0 @xi
ð1Þ
q
@ ui @ ui @p @ 2 ui ¼ þ leff qgbT ðT T 0 Þ þ q uj @t @xj @xi @xj @xj ! qð2x ui þx ðx r ÞÞ
ð2Þ
@T @T @2 T ¼ aeff þ uj @t @xj @xj @xj
ð3Þ
where the ‘‘overline” denotes filtered variables. x is the angular velocity of reference frame, bT is the thermal expansion coefficient of silicon melt, T0 is the reference temperature and g is the gravitational acceleration. The effective eddy viscosity and effective thermal diffusivity are defined as leff = l + lSGS and aeff = a + lSGS/rSGS respectively. And rSGS is SGS turbulent Prandtl number with a value of 0.9. The dynamic eddy-viscosity model is given as:
lSGS ¼ qC S D2 S
ð4Þ
ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi @u where S ¼ 2 Sij Sij , and Sij ¼ 12 @@xuji þ @xij is the strain rate tensor.
D ¼ ðD1 D2 D3 Þ1=3 is the grid-filter width, and D1, D2, D3 are the grid
Fig. 1. The illustration of computational model. (a) Geometry and boundary conditions of the computational model (2D). (b) Sketch of partial computational grids (3D). Here, a few grids are displayed for better visibility.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
spacing in each space direction respectively. The model coefficient CS is a local time-dependent variable obtained by the approach proposed by Germano [23] and modified by Lily [24]. As known from Ghosal [25], the filtering operation and the differential operation do not commute in non-uniform grids, leading to the occurrence of numerical errors, which always exists when solving the complex geometry model. Owing to the mesh is uniform in computational space, we solve the governing equations in computational space to eliminate the numerical errors. And the SGS eddy viscosity lSGS is calculated with the dynamic SGS model for body-fitted grids described in our previous paper [5], which can satisfactorily predict turbulent fluctuation in the melt with much shorter CPU time and much smaller computer memory during the CZ crystal growth. Besides, in the regions near the crucible wall and crystal/melt interface, lSGS was estimated with a three-layer log-law wall function [26,27]. The LES governing equations are solved by using the finite volume method on non-orthogonal curvilinear grids. SIMPLEC algorithm is adopted to couple the velocity and pressure fields. The second-order central differencing scheme is applied to diffusion term, and QUICK scheme is used for convective term by deferred correction method [28]. For the time term, second-order implicit backward differencing scheme is adopted. The present numerical code for LES calculation has been validated by comparison with DNS results of Raufeisen et al. [29], which can be seen in our previous work [5]. The physical properties of silicon used in this work are provided in Table 1. 3.2. Boundary conditions No-slip velocity boundary conditions are applied at the crucible wall and crystal/melt interface. At the melt free surface, Marangoni convection is considered, which is induced by the surface tension caused by temperature gradients at the free surface. So, the shear stress is balanced with surface tension as follows [5]:
sgg ¼ l
@ug ¼0 @g
ð5Þ
sgn ¼ l
@un @r dr @T ¼ ¼ @g @n dT @n
ð6Þ
sgf ¼ l
@uf @r dr @T ¼ ¼ @g @f dT @f
ð7Þ
where n and f are local coordinates tangential to the melt surface and g is orthogonal to both n and f. The temperature coefficient of surface tension dr/dT is a constant (see Table 1). At the melt free surface and crystal surface, the heat loss is considered by Stefan–Boltzmann equation:
Table 1 The physical properties of silicon used in the work. Parameters
Symbols
Values
Density, crystal Density, melt Thermal conductivity, crystal Thermal conductivity, melt Heat capacity, crytsal Heat capacity, melt Dynamic viscosity Thermal expansion coefficient Latent heat of phase change Surface tension coefficient
qc qm
dr dT
2330 kg/m3 2530 kg/m3 22 W/(mK) 64 W/(mK) 1059 J/(kgK) 1000 J/(kgK) 7 104 kg/(ms) 1.4 104 K1 1.411 106 J/kg 1.0 104 N/(km)
Tm
0.58 0.3 1685 K
Emissivity factor, crystal Emissivity factor, melt Melting temperature
kc km Cp-c Cp-m
l
bT DH
ec em
keff
@T ¼ rSB eðT 4 T 4env Þ @n
3
ð8Þ
where keff is the effective thermal conductivity. rSB is the StefanBoltzmann constant. The ambient temperature Tenv of melt free surface is 1620 K, and the ambient temperature Tenv of crystal surface decreases linearly with the increase of crystal height z (see Fig. 1 (a)), from 1620 K to 500 K, which have been proved to be effective for simulating the crystal heat dissipation compared with global simulation [30]. According to the experimental data [31], a parabolic function relationship between the temperature profile Tcru at crucible wall and the crucible height is considered (see Fig. 1(a)). The minimum temperature is located on the bottom and the upper edge of the crucible. And the maximum temperature is located on the middle of the crucible wall. At the crystal/melt interface, the heat flux is balanced considering latent heat release, which satisfies the following equation:
@T @T k k ¼ qc DHV n @n c @n m
ð9Þ
where the subscript c and subscript m denote values of crystal and melt side. DH is the latent heat of crystallization. And Vn is the local normal growth rate, which is related to the undercooling on the interface. The undercooling on the interface is necessary for crystal growth, which is the difference between theoretical crystallization temperature and actual crystallization temperature. Fujiwara et al. [9] have deeply studied the relationship between growth rate and undercooling by conducting silicon growth experiments in recent years. They focus on the low undercooled melt, which is similar to the industrial growth process of bulk silicon crystal when using CZ or Floating zone (FZ) methods. They found the growth rate linearly increased with the undercooling degree:
V n ¼ 7:5 106 T m T f m=s
ð10Þ
where Tm is the melting temperature of silicon, and Tf is the local temperature on the crystal/melt interface. This relationship is adopted in this work. 4. Preparation for LES computation 4.1. Initial field and crystal/melt interface shape In the initial stage of LES calculation, reasonable temperature and velocity fields are necessary. In this work, the initial temperature and velocity fields used for LES are calculated by the RANS model presented in our previous work [13]. Besides, the crystal/ melt interface is a coupling boundary between melt and crystal. In order to obtain temperature fluctuation on the interface, different from previous studies [5–8], in the current study, temperature on the interface is not fixed, which is related to the heat flux near the interface and the release of latent heat. And the growth rate is obtained by Eq. (10), according to the local temperature Tf on the interface. The initial isothermal interface related to crystal growth rate is required for the reasonable calculation of temperature fluctuation on the interface by LES. Owing to the heat generated by heater (not shown in this work, see Ref. [13]) is transferred to melt through crucible wall, thus the crucible wall is considered as the heat source in present computational model. Therefore, the temperature Tcru on crucible wall need to be adjusted to ensure that the crystal/melt interface is isothermal. In the process of achieving an isothermal interface, crystal pulling velocity in the iso-diameter growth stage is selected as 0.8 mm/ min, thus temperature on the interface should be 1683.23 K based on the relationship between growth rate and undercooling (see
4
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
formula (10)). Besides, considering the corresponding temperature Tcru on crucible wall is not easy to known when the interface temperature is 1683.23 K, the standard form of PID controller is adopted to adjust the temperature Tcru conveniently, which as follows:
Z 1 t deðtÞ DT cru ðtÞ ¼ K p eðtÞ þ eðtÞdt þ sd dt si 0
ð11Þ
where DTcru(t)is the change of temperature on crucible wall, then the manipulated variable Tcru is adjusted to Tcru = Tcru(0) + DTcru(t) with Tcru(0) is previous temperature on crucible wall.
e(t) = Tref T tri (t) is difference between the reference temperature
Tref (1683.23 K) and current mean temperature T tri ðtÞ on the trijunction edge considering the circumferential temperature is not perfectly uniform in the 3D model. Kp is the gain, si and sd are the integral time and the derivate time, respectively. Due to the control process is a timeless problem, iteration step is used as the sampling interval in this calculation. As a starting point, we first present the distribution of monitoring points used in the calculation as shown in Fig. 2, as a basis for further analysis in subsequent sections. The situation of xcry = 8 rpm and xcru = 6 rpm is selected to illustrate the calculation process. The calculation algorithm involves two steps: (1) Using PID controller adjust the temperature on crucible wall
to make T tri close to Tref. This stage provides reasonable temperature and velocity fields for the next moving grid stage, ensuring the mesh movement can be carried out smoothly. Fig. 3 shows the convergence process of temperature at
monitoring points and T tri . It is obvious that the tempera-
tures in the melt and on the interface remain unchanged
approximately after 12,000 iteration steps, and T tri eventually stabilizes to the target temperature Tref. This illuminates that PID controller is an effective and convenient method for
controlling the crucible wall temperature to ensure T tri approach Tref. At the end of this stage, the crystal/melt interface is a flat non-isothermal interface as shown in Fig. 4. And the minimum temperature is located at the tri-junction edge. (2) Moving mesh to ensure the crystal/melt surface is isothermal and the temperatures in melt and crystal are converged.
During the mesh movement, T tri is also controlled by adjusting crucible wall temperature using PID controller. The location of the crystal/melt interface is determined in such a way that it is updated iteratively in the iterations according to
(Tf T tri ), until (Tf T tri ) approaches zero. As can be seen in
Fig. 5, with the grids moving, the temperature T tri and crucible wall temperature Tcru tend to stabilize gradually. After
two hundred mesh move steps, the temperature T tri reach to the reference temperature Tref 1683.23 K. In addition, the maximum temperature difference at the crystal/melt interface decrease with moving mesh. After three hundred move steps, the maximum temperature difference is less than 0.2 K. Thus, it is considered that the crystal/melt interface is approximate isothermal surface at this moment. At the end of this stage, the crystal/melt interface is an isothermal curved surface as shown in Fig. 6. 4.2. Selection of time step Considering the computational accuracy and numerical cost comprehensively, the grid size of 50 50 85 in the three coordinate directions of curvilinear coordinate system is used. Owing to the temperature fluctuation near the crystal/melt interface is the focus of our attention, the grids near the interface is specially refined. In order to obtain a solution independent of time step,
Fig. 2. Distribution of monitoring points on the crystal/melt interface. P1-P4 are close to the center region of interface, P5-P8 are near the edge. In addition to the above eight monitoring points, there are also eight monitoring points P1_B-P8_B in the melt located at depth of 20 mm from the interface corresponding to the eight points on the interface respectively, which are not shown for better visibility.
Fig. 3. Varieties of temperature at monitoring points with iteration steps.
Fig. 4. Temperature field in the melt at the end of step one.
Fig. 5. Process of obtaining the isothermal crystal/melt interface.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
5
Fig. 6. Temperature field in the melt at the end of step two. Fig. 8. Time-averaged fields of melt flow in the meridian plane when xcry = 8 rpm and xcru = 6 rpm. The left part: temperature distribution. The interval is 0.5 K. The right part: stream function distribution. The interval is 0.05 kg/s.
Fig. 9. Time-averaged fields of melt flow in the meridian plane when xcry = 8 rpm and xcru = 6 rpm. The left part: temperature distribution. The interval is 0.5 K. The right part: stream function distribution. The interval is 0.05 kg/s.
Fig. 7. Spectral analyses of temperature fluctuation at point P2_B for different time steps. (a) 0.02 s time step. (b) 0.01 s time step. (c) 0.005 s time step.
the effects of time steps on the characteristics of temperature fluctuation near interface are investigated. Fig. 7 shows the power spectra density (PSD) distributions of temperature fluctuation at point P2_B in melt for different time steps 0.02 s, 0.01 s, 0.005 s, respectively. A total of 60,000 values in a time range of 600 s are used for the spectral analysis by fast Fourier transformed method, which is also used in Section 5. Obviously, when the time steps are selected as 0.01 s and 0.005 s, there is an identical characteristic frequency of the temperature fluctuation. However, the identical characteristic frequency disappears with time step 0.02 s, indicating that it is not suitable for the analysis of characteristics of temperature fluctuation in the melt. Considering the computing efficiency, thus the time step 0.01 s is selected in the subsequent calculations. 5. Results and discussion 5.1. Time-averaged field of melt flow To analyze the statistical behavior of the melt flow, the velocity and temperature fields are averaged in time and in the circumferential direction. Figs. 8 and 9 display the distributions of statistically averaged temperature and stream function for different crystal and crucible rotation conditions. Both for counter-rotation and co-rotation, it is obvious that three convection rolls exist in the melt, the Taylor-Proudman vortex, the buoyant vortex and a secondary vortex. The counterclockwise buoyant vortex is driven by thermal buoyancy, which is induced by the axial temperature
difference in the melt, holds the great area of melt region. The buoyant melt flow upward along the crucible wall, then reach the melt free surface, causing the clockwise secondary vortex around the corner. In the center melt region below the crystal, the buoyant vortex is damped by the clockwise Taylor-Proudman vortex, which is induced by the rotation of crystal and crucible. The difference is that when the crystal counter-rotates with crucible, the Taylor-Proudman vortex is larger than that in corotation. The larger Taylor-Proudman vortex, the more complex melt flow under the crystal. With regard to the time-averaged temperature, it is found that the trend of temperature distribution in the melt remains unchanged for different crystal and crucible rotation conditions. However, by comparing the distribution of isotherms, the melt temperature for counter-rotation is higher than that for co-rotation, especially near the crystal/melt interface. This indicates that more heat is transferred to the crystal/melt interface when crystal counter-rotates with crucible. 5.2. Characteristic of temperature fluctuation In this section, we firstly study the characteristics of temperature fluctuation on the crystal/melt interface and that in the melt near interface. In order to display the information of temperature fluctuation more clearly, temperature fluctuations at points P1P8 are presented by boxplots as shown in Fig. 10. The boxplot is a flexible exploratory-data-analysis tool for graphically representing the key values from numerical data in detail. From the boxplot, information about the location, spread, skewness, and longtailedness of observations can be quickly known [32]. In the paper, the whiskers are determined by the maximum and minimum temperatures at monitoring points. The square in the box represents the mean value of temperature fluctuations. The horizontal line in the box expresses the median of temperature data. The upper edge of the box represents the 75th percentile, and the lower edge is the 25th percentile. For details of the box diagram, one can see Ref. [32]. In summary, the height of the box or the distance between
6
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Fig. 10. Temperature fluctuations at monitoring points P1-P8 on the interface for different crystal and crucible rotation conditions. (a) xcry = 8 rpm, xcru = 6 rpm. (b) xcry = 8 rpm, xcru = 6 rpm.
the ends of whiskers represents the range of temperature fluctuation approximatively. From Fig. 10, it can be noticed that when crystal counter-rotates with crucible, both the height of the box and the distance between the ends of whiskers are longer compared to that when crystal and crucible rotate in the same direction. This observation illustrates that temperature oscillation on the crystal/melt interface becomes more intense when crystal and crucible rotate in the opposite direction. Besides, the changes of cut line at the median which represents location of box are smaller when crystal counter rotate with crucible, especially in the central region of interface. This reveals that temperature fluctuation on the interface is more uniform for counter-rotation. Considering Fig. 10(a) and (b) synthetically, based on the height of box, the conservative temperature fluctuation on the crystal/melt interface is about 0.4 K for corotation and 0.7 K for counter-rotation. Fig. 11 shows the temperature fluctuations at points P1_B-P8_B in the melt. Similarly, the temperature oscillation in the melt near interface is more intense for counter-rotation compared with that obtained for co-rotation. Combining the results of time-averaged melt flow in Section 5.1, the probable explanation for this observation is as follows: when the crystal counter-rotates with the crucible, the Taylor-Proudman vortex always appears in the region below the crystal/melt interface (see Fig. 5 in Ref. [33]). In the Taylor-Proudman vortex, the melt near crucible wall spirals up to the region near interface, resulting in the intense mixing of melt below crystal, then giving rise to the larger temperature fluctuation near interface. However, the Taylor column vanishes when the rotation rate of crystal and crucible is same both in magnitude and direction [33]. This explanation is also proved by the analysis
of local heat transfer in the melt in Section 5.4. Carefully, the rotation of crucible is slightly higher than that of crystal in our present simulation. Thus, a small Taylor-Proudman vortex exist in our simulation, as shown in Fig. 9, which also been observed in Ref. [33]. From Fig. 11, it can also be found that the ranges of temperature fluctuation at all monitoring points are approximately identical for counter-rotation of crystal and crucible, while the temperature fluctuation in the central region is weaker than that near the crystal edge for co-rotation. Comprehensively, based on the height of box, the temperature fluctuations are about 2.8 K and 2.0 K for counter-rotation and co-rotation, respectively. In Refs. [7,31], experiments about the observation of temperature fluctuation near interface for counter-rotation of crystal and crucible have been performed. The difference between the maximum and minimum temperatures is 9.2 K in Ref. [7], while it is 10 K in Ref. [31]. According to the distance between the ends of whiskers in Fig. 11(a), the range of temperature fluctuation in the melt below crystal in our numerical results are in good agreement with those obtained by experiments. 5.3. Relation between temperature fluctuation on the interface and that in the melt The relationship between the temperature fluctuation on the interface and that in the melt close to crystal is analyzed in this section. The cross correlation as a general method of estimating the relevance of two time series is adopted in the work. Fig. 12 illustrates the maximum cross-correlation coefficients of temperature fluctuations at monitoring points P1-P8 distributed on the interface and points P1_B-P8_B in the melt for different crystal
Fig. 11. Temperature fluctuations at monitoring points P1_B-P8_B in the melt for different crystal and crucible rotation conditions. (a) xcry = 8 rpm, xcru = 6 rpm. (b) xcry = 8 rpm, xcru = 6 rpm.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
7
Fig. 13. The cross-correlations of temperature fluctuations at monitoring points P2 and P6 on the interface with their corresponding points P2_B and P6_B in the melt respectively when xcry = 8 rpm, xcru = 6 rpm. (a) P2 and P2_B. (b) P6 and P6_B.
Fig. 12. The maximum cross-correlation coefficients of temperature fluctuations on the interface with that in the melt for different crystal and crucible rotation conditions. (a) xcry = 8 rpm, xcru = 6 rpm. (b) xcry = 8 rpm, xcru = 6 rpm.
and crucible rotation conditions. It is obvious that there is a strong relationship between the temperature variation at the points distributed on the interface and that just below these points. And the correlation is slightly enhanced for counter-rotation compared with co-rotation. In the calculation process, a similar characteristics of temperature fluctuations at points in the center region of interface and melt close to crystal has been found. And the similarity also exists for the points near the edge. Thus, for the sake of clarity, the points P2 (P2_B) and P6 (P6_B) are selected to represent the center region and edge area of the interface (melt) in the latter discussion, respectively. Fig. 13 illustrates the cross correlation series with a maximum time delay of 600 s when the crystal counter-rotates with crucible in detail. There is a strong correlation at the time delay of about 0.61 s at points P2 and its corresponding point P2_B. And the maximum correlation is achieved at the time delay of 0.54 s at points P6 and P6_B. Besides, when crystal rotate with crucible in the same direction, the time delays related to the maximum correlation are 0.41 s and 0.2 s for P2&P2_B and P6&P6_B respectively, as shown in Fig. 14. The positive value of time delay represents that the temperature variation on the crystal/melt interface lags behind that in the melt. According to our results, it is found that the values of time delays at all points P1-P8 with their corresponding points P1_B-P8_B for counter-rotation and that at the points P1-P4 with their corresponding points P1_B-P4_B for co-rotation are positive. On the contrary, the values of time delays are negative in the region close to fringe of crystal when crystal co-rotates with crucible. The reasons for this will be explained in the next section. The frequency analyses of temperature fluctuations on the crystal/melt interface and in the melt are performed, in order to obtain the characteristics of temperature fluctuation, as shown in Figs. 15 and 16. For the convenience of comparison, the PSD is normalized by its maximum value, which is denoted as PSDN. It can be seen
Fig. 14. The cross-correlations of temperature fluctuations at monitoring points P2 and P6 on the interface with their corresponding points P2_B and P6_B in the melt respectively when xcry = 8 rpm, xcru = 6 rpm. (a) P2 and P2_B. (b) P6 and P6_B.
that there are good agreements for dominant frequencies at P2 and P6 with their corresponding points in the melt when crystal counter-rotates with crucible. However, the consistency of dominant frequencies is not obvious for co-rotation. And no frequencies consistent with the crystal or crucible rotational frequency are found in two figures, which are 0.133 Hz and 0.1 Hz respectively. For the counter-rotation, the dominant frequency at P2 is 0.0136 Hz, one-tenth of crystal rotational frequency, which is probably influenced by crystal rotation. Different from the corotation, PSDN is strong in the lower frequency range for counterrotation, indicating that larger periodic flows occur in the melt below the crystal/melt interface when crystal counter-rotates with crucible. The reason for this is that when crystal counter-rotates with crucible, the melt flow in the crucible is dominated by the buoyancy, which is a large periodic circulation passing crucible wall, melt free surface and bottom of crystal successively [5,8]. While when crystal co-rotates with crucible, the melt flow is mainly affected by rotations of crystal and crucible, which is characterized by small period [33].
8
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Fig. 15. Frequency analyses of temperature fluctuations at points P2, P6 on the interface and the corresponding points P2_B, P6_B in the melt when xcry = 8 rpm, xcru = 6 rpm.
Fig. 16. Frequency analyses of temperature fluctuations at points P2, P6 on the interface and the corresponding points P2_B, P6_B in the melt when xcry = 8 rpm, xcru = 6 rpm.
5.4. Analysis of heat transfer in the melt The instantaneous local convective heat flux proposed by Shang et al. [34,35] is adopted to study the heat transfer in the melt near crystal/melt interface, according to local velocity U(r,t) and local temperature T(r,t). The definition of instantaneous local convective heat flux is as follows:
jðr; tÞ ¼
C pm qm U ðr; t Þ½T ðr; tÞ T 0 km DT=L
ð12Þ
where T0 is the mean temperature of bulk melt, DT is the temperature difference between crucible and crystal/melt interface, L is the depth of melt. The local convective heat flux is normalized by conductive heat flux kmDT/L. Besides, Shang et al. proved that the diffusive contribution can be ignored compared to the convective heat flux [34]. According to the Eq. (12), it is known that the warm fluctuation (T(r, t) T0) > 0) produces positive flux when they move
upward to the crystal/melt interface (UZ > 0), and the falling (UZ < 0) cold plume (T(r, t) T0) < 0) also contributes to a positive flux. The positive flux at point P2_B or P6_B in the melt means that the heat is transferred to crystal/melt interface.
5.4.1. Crystal counter-rotates with crucible Fig. 17 shows the time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuations at points P2_B and P6_B in the melt close to crystal. At point P2_B, as shown in Fig. 17(a), the fluctuation of circumferential velocity is slightly skewed toward the positive direction, which is directly influenced by rotations of crystal and crucible. At point P6_B, the region close to fringe of crystal, the skew is further strengthened as shown in Fig. 17(b). While, the radial and axial velocity fluctuate around zero value both at points P2_B and P6_B. It is inferred that the effects of crystal and crucible rotation on melt flow in the radial and axial directions are small at the loca-
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
9
Fig. 17. Time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuation in the melt when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
tion 2 cm from interface in our present work, and the melt has a free movement in both directions. Fig. 18 displays the histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions at points P2_B and P6_B, respectively. To display the histograms of j at different directions in a same graph, the histograms are normalized by their maximum value H0. When the thermal plumes exists, the temperature fluctuation is related to the velocity fluctuation, producing positive fluctuations [34]. In Fig. 18(a), it is seen that the H(jUT)/H0 and H(jUR)/H0 are symmetric and mean values of jUT and jUR are zero, indicating that no such correlation exists between the temperature and circumferential, radial velocity fluctuations. Duo to the falling cold plums and rising warm plumes produce positive fluxes in the vertical direction as mentioned above, the H(jUZ)/H0 becomes highly skewed toward the positive direction, even though velocity fluctuations are symmetric as shown in Fig. 17(a). The similar principle can be used to explain the symmetry of circumferential heat flux, although the circumferential velocity is slightly skewed toward the positive direction. Besides, the same shape of H(jUT)/H0 and H(jUR)/H0 suggests that there are same statistical characteristics of instantaneous local heat flux in the two directions. In Fig. 18(b), for point P6_B, the region close to fringe of crystal, because of the obvious circumferential motion, the circumferential heat flux becomes highly skewed toward the positive
direction with a long decay tail, indicating that the heat flux fluctuations produced by thermal plumes contain some large but rare events [34]. Besides, the radial heat flux is slightly skewed toward the negative direction, approximately same as the negative flux fluctuations in circumferential direction. Similar to P2_B, H(jUZ)/ H0 of P6_B also becomes skewed toward the positive direction with the same shape of P2_B. In summary, when crystal counter-rotates with crucible, the heat transfer is dominated by axial local heat flux in the center region of melt near crystal. While in the region close to the fringe of crystal, the heat transfer is complex with the asymmetries of local heat flux in three directions.
5.4.2. Crystal co-rotates with crucible Fig. 19 shows the time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuations at points P2_B and P6_B in the melt close to crystal when crystal co-rotates with crucible. Similar to the situation for counter-rotation, the fluctuations of radial and axial velocity components are symmetric with zero mean values both at P2_B and P6_B. The difference is that the fluctuation of circumferential velocity is skewed toward the negative direction when crystal co-rotates with crucible. From Figs. 17 and 19, it is found that the direction of circumferential velocity at the location 2 cm from
Fig. 18. Histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
10
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Fig. 19. Time series of the circumferential velocity (UT), radial velocity (UR), axial velocity (UZ), and temperature fluctuation in the melt when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
interface is the same as that of crucible rotation, indicating that the melt flow in crucible is mainly influenced by crucible rotation. Fig. 20 illustrates the histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions at P2_B and P6_B, respectively. Similar to the situation for counter-rotation, at P2_B, the fluctuations of the circumferential and radial heat flux are symmetric with zero mean values, and the H(jUZ)/H0 also becomes slightly skewed toward the positive direction. By comparing Figs. 18(a) and 20(a), it is clear that the values of instantaneous local heat flux j at P2_B for counter-rotation are larger than that for co-rotation, which is the reason for the occurrence of more violent temperature oscillations when crystal counter-rotates with crucible, especially in the central area of the interface, as displayed in Fig. 11. For point P6_B, duo to the reversal of the direction of crucible rotation for co-rotation compared with counter-rotation, the H(jUT)/H0 becomes highly skewed toward the negative direction. And the H(jUR)/H0 is skewed toward the negative direction. It is important to note that the axial heat flux is also slightly skewed toward the negative direction, which is clearly different from the situation in P2_B for co-rotation and that in P2_B and P6_B for counter-rotation. This is the real cause of the only negative delay time of P6_B when crystal co-rotates with crucible, while which is positive in other three conditions.
6. Conclusions In this paper, temperature fluctuations on the crystal/melt interface during CZ silicon crystal growth are investigated by LES with the curvilinear dynamic Smagorinsky subgrid-scale model, considering the relationship between kinetic undercooling and growth rate. And the influence mechanism of melt flow instability near interface on the temperature fluctuations are clarified in detail. The following conclusions are obtained from this study: (1) When crystal counter-rotates with crucible, temperature fluctuations on the crystal/melt interface and that in the melt close to crystal are more intense compared with corotation. And the temperature fluctuations on the crystal/ melt interface are about 0.4 K for co-rotation and 0.7 K for counter-rotation, respectively. (2) The temperature variation at points distributed on the interface has a strong relationship with that just below these points. And the temperature variation on the crystal/melt interface lags behind that in the melt for the counterrotation and in the central region when co-rotation occurs. However, the situation is exactly opposite in the region close to fringe of crystal when crystal co-rotates with crucible.
Fig. 20. Histograms of the instantaneous local heat flux j in the circumferential, radial and axial directions when xcry = 8 rpm, xcru = 6 rpm. (a) P2_B point. (b) P6_B point.
J. Ding, L. Liu / International Journal of Heat and Mass Transfer 142 (2019) 118463
Besides, the characteristic frequencies of temperature fluctuations on the interface and that in the melt are identical for counter-rotation. However, this consistency of dominant frequencies is not obvious for co-rotation. (3) For both counter-rotation and co-rotation, the heat transfer is dominated by axial local heat flux in the center region of melt close to crystal. The difference is that larger local heat flux exists for counter-rotation compared with co-rotation. And in the region close to fringe of crystal, the heat transfer is complex with the asymmetries of local heat flux in three directions. Besides, the axial heat transfer in the region close to the fringe of crystal for co-rotation is significantly different from the situations of counter-rotation and that in the central region for co-rotation. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 51676154). References [1] Z. Lu, S. Kimbel, Growth of 450mm diameter semiconductor grade silicon crystals, J. Cryst. Growth 318 (2001) 193–195. [2] E. Kamiyama, J. Vanhellemont, K. Sueoka, K. Araki, K. Izunome, Thermal stress induced void formation during 450 mm defect free silicon crystal growth and implications for wafer inspection, Appl. Phys. Lett. 102 (2013) 082108. [3] T. Shen, C.M. Wu, Y.R. Li, Experimental investigation on the effect of crystal and crucible rotation on thermocapillary convection in a Czochralski configuration, Int. J. Therm. Sci. 104 (2016) 20–28. [4] A. Muiznieks, A. Krauze, B. Nacke, Convective phenomena in large melts including magnetic fields, J. Cryst. Growth 303 (2007) 211–220. [5] L.J. Liu, X. Liu, Y. Wang, Large-eddy simulation of melt turbulence in a 300-mm Cz–Si crystal growth, Int. J. Heat Mass Transf. 55 (2012) 53–60. [6] H. Hirata, K. Hoshikawa, Three-dimensional numerical analyses of the effects of a cusp magnetic field on the flows, oxygen transport and heat transfer in a Czochralski silicon melt, J. Cryst. Growth 125 (1992) 181–207. [7] D. Vizman, O. Grabner, G. Muller, 3D numerical simulation and experimental investigations of melt flow in an Si Czochralski melt under the influence of a cusp-magnetic field, J. Cryst. Growth 236 (2002) 545–550. [8] X.R. Cen, Y.S. Li, J.M. Zhan, Three dimensional simulation of melt flow in Czochralski crystal growth with steady magnetic fields, J. Cryst. Growth 340 (2012) 135–141. [9] K. Fujiwara, K. Maeda, N. Usami, G. Sazaki, Y. Nose, A. Nomura, T. Shishido, K. Nakajima, In situ observation of Si faceted dendrite growth from low-degreeof-undercooling melts, Acta Mater. 56 (2008) 2663–2668. [10] V.V. Voronkov, The mechanism of swirl defects formation in silicon, J. Cryst. Growth 59 (1982) 625–643. [11] J. Vanhellemont, Intrinsic point defect incorporation in silicon single crystals grown from a melt, revisited, J. Appl. Phys. 110 (2011) 063519. [12] K. Kakimoto, T. Shyo, M. Eguchi, Correlation between temperature and impurity concentration fluctuations in silicon crystals grown by the Czochralski method, J. Cryst. Growth 151 (1995) 187–191.
11
[13] L.J. Liu, K. Kakimoto, Partly three-dimensional global modeling of a silicon Czochralski furnace. I. Principles, formulation and implementation of the model, Int. J. Heat Mass Transf. 48 (2005) 4481–4491. [14] A. Chatterjee, V. Prasad, A full 3-dimensional adaptive finite volume scheme for transport and phase-change processes, part I: formulation and validation, Numer. Heat Tranf. A-Appl. 37 (2000) 801–821. [15] L.J. Liu, Q.H. Yu, X.F. Qi, W.H. Zhao, G.X. Zhong, Controlling solidification front shape and thermal stress in growing quasi-single-crystal silicon ingots: process design for seeded directional solidification, Appl. Therm. Eng. 91 (2015) 225–233. [16] V. Kumar, G. Biswas, G. Brenner, F. Durst, Effect of thermocapillary convection in an industrial Czochralski crucible: numerical simulation, Int. J. Heat Mass Transf. 46 (2003) 1641–1652. [17] S.S. Son, K.W. Yi, Characteristics of thermal fluctuation in a low Pr number melt at a large crucible for Czochralski crystal growth method, J. Cryst. Growth 275 (2005) E259–E264. [18] C. Wagner, R. Friedrich, Direct numerical simulation of momentum and heat transport in idealized Czochralski crystal growth configurations, Int. J. Heat Fluid Flow 25 (2004) 431–443. [19] T. Shen, C.M. Wu, L. Zhang, Y.R. Li, Experimental investigation on effects of crystal and crucible rotation on thermal convection in a model Czochralski configuration, J. Cryst. Growth 438 (2016) 55–62. [20] T. Duffar, Crystal Growth Processes Based on Capillarity: Czochralski, Floating Zone, Shaping and Crucible Techniques, Wiley, Wiltshire, 2010, pp. 479–481. [21] J.F. Thompson, Z.U.A. Warsi, C.W. Mastin, Boundary-fitted coordinate systems for numerical solution of partial differential equations-a review, J. Comput. Phys. 47 (1982) 1–108. [22] E.M.J. Komen, L.H. Camilo, A. Shams, B.J. Geurts, B. Koren, A quantification method for numerical dissipation in quasi-DNS and under-resolved DNS, and effects of numerical dissipation in quasi-DNS and under-resolved DNS of turbulent channel flows, J. Comput. Phys. 345 (2017) 565–595. [23] M. Germano, U. Piomelli, P. Moin, W. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A 3 (1991) 1760–1765. [24] D.K. Lilly, A proposed modification of the Germano subgrid-scale closure method, Phys. Fluids A 4 (1992) 633–635. [25] S. Ghosal, P. Moin, The basic equations for the large eddy simulation of turbulent flows in complex geometry, J. Comput. Phys. 118 (1995) 24–37. [26] I.Y. Evstratov, V.V. Kalaev, A.I. Zhmakin, Y.N. Makarov, A.G. Abramov, N.G. Ivanov, E.M. Smirnov, E. Dornberger, J. Virbulis, E. Tomzig, W.V. Ammon, Modeling analysis of unsteady three-dimensional turbulent melt flow during Czochralski growth of Si crystals, J. Cryst. Growth 230 (2001) 22–29. [27] L. Temmerman, M.A. Leschziner, C.P. Mellen, J. Frohlich, Investigation of wallfunction approximations and subgrid-scale models in large eddy simulation of separated flow in a channel with streamwise periodic constrictions, Int. J. Heat Fluid Flow 24 (2003) 157–180. [28] T. Hayase, J.A.C. Humphrey, R. Greif, A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures, J. Comput. Phys. 98 (1992) 108–118. [29] A. Raufeisen, M. Breuer, T. Botsch, A. Delgado, DNS of rotating buoyancy- and surface tension-driven flow, Int. J. Heat Mass Transf. 51 (2008) 6219–6234. [30] J.L. Ding, L.J. Liu, Real-time prediction of crystal/melt interface shape during Czochralski crystal growth, Crystengcomm 20 (2018) 6925–6931. [31] S. Enger, O. Grabner, G. Muller, M. Breuer, F. Durst, Comparison of measurements and numerical simulations of melt convection in Czochralski crystal growth of silicon, J. Cryst. Growth 230 (2001) 135–142. [32] Y. Benjamini, Opening the box of a boxplot, Am. Stat. 42 (1988) 257–262. [33] J.R. Carruthers, K. Nassau, Nonmixing cells due to crucible rotation during Czochralski crystal growth, J. Appl. Phys. 39 (1968) 5205–5214. [34] X.D. Shang, X.L. Qiu, P. Tong, K.Q. Xia, Measurements of the local convective heat flux in turbulent Rayleigh-Benard convection, Phys. Rev. E 70 (2004) 026308. [35] X.D. Shang, P. Tong, K.Q. Xia, Scaling of the local convective heat flux in turbulent Rayleigh-Benard convection, Phys. Rev. Lett. 100 (2008) 244503.