Next Article in Journal
Flexible and Structural Coloured Composite Films from Cellulose Nanocrystals/Hydroxypropyl Cellulose Lyotropic Suspensions
Next Article in Special Issue
Side-Chain Liquid Crystal Co-Polymers for Angular Photochromic Anti-Counterfeiting Powder and Fiber
Previous Article in Journal
Synthesis, X-ray Single Crystal, Conformational Analysis and Cholinesterase Inhibitory Activity of a New Spiropyrrolidine Scaffold Tethered Benzo[b]Thiophene Analogue
Previous Article in Special Issue
Interaction and Polarization Energy Relationships in σ-Hole and π-Hole Bonding
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

3D Numerical Analysis of the Asymmetric Three-Phase Line of Floating Zone for Silicon Crystal Growth

Research Institute for Applied Mechanics, Kyushu University, 6-1 Kasuga-koen, Kasuga, Fukuoka 816-8580, Japan
*
Author to whom correspondence should be addressed.
Crystals 2020, 10(2), 121; https://doi.org/10.3390/cryst10020121
Submission received: 24 January 2020 / Revised: 13 February 2020 / Accepted: 14 February 2020 / Published: 15 February 2020

Abstract

:
A numerical simulation has been carried out to study the asymmetric heat transfer, fluid flow, and three-phase line to explain the phenomenon of the spillage of the melt in floating zone (FZ) silicon growth. A three-dimensional high-frequency electromagnetic (EM) field is coupled with the heat transfer in the melt and crystal calculation domains. The current density along the three-phase line is investigated to demonstrate the inhomogeneous heating along the three-phase line. The asymmetric heating is found to affect the flow pattern and temperature distribution of the melt. The three-dimensional solid–liquid interface results show that, below the current supplies, the interface is deflected due to strong heating below the current supplies. The calculated asymmetric three-phase line shows a similar trend as the experimentally observed results. The results indicate that the re-melting and spillage phenomenon could occur below the current supplies.

1. Introduction

With the expansion of the power device market, highly cost effective floating zone (FZ) silicon has been developed. A needle-eye inductor is used to provide high-frequency electromagnetic (EM) heating in FZ silicon with a diameter larger than 100 mm. Because the needle-eye inductor is a one-turn coil, asymmetric heating becomes more obvious when the diameter of silicon is larger according to the previous study [1]. The asymmetric heating can not only affect the melting speed of the feed rod but also affect the shape of the solid–liquid interface between the melt and the single crystal. Inhomogeneous heat generation at the three-phase line could lead to spillage of the melt or to the formation of a bulge at the grown single crystal.
According to a previous study [1], the FZ experiment for silicon crystal was conducted. During the process, the inhomogeneous three-phase line and spillage phenomena have been observed. The three-phase line is deflected downward under the main slit of the inductor. This phenomenon needs to be analyzed by numerical studies. The first 2D numerical study on the FZ was published by Mühlbauer et al. [2]. The 2D global numerical studies were later conducted and were found to provide a good fit with experimental results [3,4,5,6,7]. However, the physical phenomena along the azimuthal direction cannot be analyzed with a 2D numerical model. Therefore, 3D numerical simulation results for the EM field, heat transfer, and impurities have been conducted [8,9,10,11,12]. The asymmetric three-phase line analysis by numerical simulation has not yet been reported. The asymmetric three-phase line, therefore, requires investigation. In a previous study, we proposed a 3D global model to calculate the asymmetric interface [13]. The proposed model could calculate the shape of crystallization interface in 3D by using enthalpy method to track the crystallization interface. The asymmetric 3D deflection of the interface is useful to evaluate the operating conditions in the experiment. Because serious deflection of crystallization interface periphery causes spillage of melt. In the present paper, on the basis of our previous model, the EM field model is improved using the impedance boundary condition offered by the commercial software COMSOL [1]. The three-phase line was calculated to investigate the experimental phenomena.

2. Numerical Models

Induction heating in the FZ is typically conducted at a high frequency. The typical frequency used in FZ is 3 MHz. The skin depth in the silicon melt is 0.26 mm. The dimension of the skin layer is far too small compared with the calculation domain. Therefore, an impedance boundary condition should be imposed to model such a high-frequency EM field. Figure 1 shows the model constructed with COMSOL.
The alternating current frequency is high, so the EM field is calculated on the basis of the time-harmonic assumption:
× ( μ 0 1 × A ) + σ V + j ω σ A = J e ,
where A is the magnetic vector potential, σ is electrical conductivity, μ 0 is the permeability, V is the electric potential, j is imaginary unit, ω is the angular frequency, and J e is generated current density. Table 1 shows physical properties and calculation conditions.
The hydrodynamic and temperature fields are calculated using the 3D model (Figure 2) with the finite volume method (FVM). The model was constructed using the open-source mesh software SALOME, and the calculation was performed with FzGlobalFOAM [13], which was developed on the basis of the open-source library OpenFOAM.
In case of steady-state calculation, the melt flow is solved in the following equations:
· ( ρ U ) = 0 ,
· ( ρ U U ) = · ( μ U ) p + ρ g ,
where U is the velocity vector in the melt, ρ is the melt density, μ is viscosity,   g is the gravitational acceleration constant, and p is the pressure. The temperature field is calculated by the following equation:
ρ U · ( h ) = · ( a h ) ,  
where h is enthalpy and a is thermal diffusivity. The thermal diffusivities differ at different temperatures in the melt and in the crystal as a result of the thermal conductivity of silicon being temperature dependent. The latent heat Q L is calculated as follows:
Q L = ρ V g L ,
where V g is the local crystal growth rate [14] and L is the latent heat of silicon. The surface current density distribution from the EM model is coupled with the heat generation at the free surface of the melt in the heat transfer model in the following equation [15]:
Q E M = J s 2 δ σ ,
where Q E M and J s are the surface power density and the current density at the melt free surface, respectively, and skin depth δ is calculated by the following equation:
δ = 1 π μ 0 f σ ,
where f is frequency.
The rotation of crystal is considered at the crystallization interface between crystal and melt. The crystallization interface is calculated according to the temperature. The shapes of melting interface and free surface are determined using the data from previous study [13]. The phase boundary of 3D free surface is imposed as a fixed wall. For the temperature field, first-type boundary condition is applied to the melting interface. The temperature at melting interface is melting point. Marangoni and EM forces are imposed as second-type boundary condition at free surface [15]:
F s u r f a c e = F E M + F M a = 1 4 μ 0 δ s ( J s 2 ) + γ T s ( T ) ,
where γ is surface tension. Our previous study gives more explanations of the model in detail [13]. Table 2 shows the process parameters and physical properties.
In the calculations, we assume that the shape of free surface is not changing with the temperature and flow fields. The geometry of high-frequency EM model is not updated with fluid flow and heat transfer model. The high-frequency EM model only provides surface current result. The surface current result is used to calculate the EM force F E M and EM heat generation Q E M . For the calculation of fluid flow, EM force F E M is imposed as a fixed gradient boundary condition at free surface of melt. For the calculation of heat transfer, Q E M is imposed as a fixed gradient boundary condition at free surface of melt. EM force F E M and EM heat generation Q E M are updated iteratively if the mean temperature of triple points is not equal to the melting point.

3. Results and Discussion

3.1. Current Density Analysis

The high-frequency current density is obtained from the EM calculation (Figure 3). The current density is low below the side slits. The current density is high below the tips of side slits. Compared with the current density below the three side slits, that in the vicinity of the main slit is lower. However, near the three-phase line, the current density below the current supplies is higher. Figure 4 shows the surface current density distribution along the three-phase line. The surface current density distribution below the main slit is steep and high. Except for the main slit, the surface current density below the three side slits is also higher than in the other areas because of the strong magnetic field below the tips of the side slits. This asymmetric EM field at the free surface is used to calculate the EM heat generation in the heat transfer model.

3.2. Melt Flow Analysis

3D heat transfer and fluid flow calculation were carried out, and the temperature distribution is shown in Figure 5. The pattern of temperature distribution is similar with that of current density distribution (Figure 3). The pattern is asymmetric due to the rotation of crystal. The results in the melt are shown in Figure 6. In the cross section paralleled with current supplies (X plane in Figure 6a), the maximum temperature occurs below the side slit tip. The temperature is lower below the main slit than below the side slit. The position of the solid–liquid interface is higher in the left side (below the main slit). In the plane perpendicular to the main slit (Y plane in Figure 6b), the temperature distribution is more symmetric than that of the X plane. The temperature is not strictly symmetric because the crystal rotation shifts the temperature distribution. The flow field is shown in Figure 7. The maximum velocity occurs at the inner triple point (ITP) between silicon melt, polycrystalline feed rod and gas because the temperature gradient at the inner triple point is very large. The Marangoni force at the ITP is strong. According to the asymmetric temperature distribution in Figure 6a, the velocity is asymmetric in X plane (Figure 7a). The velocity at ITP is higher in the right part than that in the left part due to larger temperature gradient below side slit than that below main slit. Ratnieks et al. also mentioned this phenomenon [8]. Figure 7b shows the flow field in the Y plane. The presence of current supplies does not strongly affect the symmetric flow pattern in the Y plane. The melt converges and flows downward due to strong EM force below the three side slits. However, below the main slit, there is no similar flow pattern driven by EM force because the EM force is weak compared to the side slits. In the actual growth process, with the rotation of the crystal, the flow pattern repeatedly varies from symmetric to asymmetric. The flow separation point, which is the origin of accumulated dopant, also changes repeatedly. These calculation results can explain the asymmetric distribution of the measured radial resistivity distribution in the as-grown crystal [4].

3.3. Three-Phase Line Analysis

The solid–liquid interface is affected by the asymmetric heating. Figure 8 shows the 3D shape of the solid–liquid interface. The color legend indicates the vertical height from the bottom of the interface. The symmetric interface center demonstrates that the asymmetric heating at the melt free surface does not strongly influence the center of the interface. However, the interface periphery is obviously inhomogeneous because the inhomogeneous heat is generated at the three-phase line (Figure 4).
The three-phase line is calculated and compared with experimentally observed results [1] in Figure 9. Below the current supplies, the position is low because of the high current density and high heat power below the current supplies. The experiment was conducted in IKZ (Leibniz Institute for Crystal Growth). From the experimental observation (Figure 9a), the three-phase line is deflected below the current supplies. The deflection is regarded as the reason for spillage down the silicon melt. In the calculation results (Figure 9b), along the rotation direction, the three-phase line descends and ascends below the current supplies. In the experimental observation, the deflection is smaller than description in Figure 9a and still visible to the naked eye [1].
Through the qualitatively comparison, the calculated deflection shows the similar trend with the experimental observation. Because of the inhomogeneity of the interface along the azimuthal positions, the local crystal growth rate is inhomogeneous at the three-phase line. Along the crystal rotation direction, the local crystal growth rate is lower than the pulling rate when the three-phase line descends and the local crystal growth rate is higher than the pulling rate when the three-phase line ascends. The latent heat, which is dependent on the local crystal growth rate, is considered using the method suggested in a previous study [14]. The inhomogeneous local growth rate causes the asymmetric three-phase line deflection below the current supplies. Moreover, when the local growth rate is lower than the pulling rate, the re-melting phenomenon occurs.
Figure 10 shows the three-phase line position along the azimuthal direction. The three-phase line is asymmetric because the rotation shifts the fluid and temperature field. The three-phase line is deflected downward in the vicinity of the main slit and side slits as a result of high temperature in the vicinity of the slits (Figure 5). The position of three-phase line is higher between −90°and 90° because the temperature is lower in the vicinity of the main slit. Below the main slit, the defection is about 1 mm. It is demonstrated that even though the crystal is rotated at 5 RPM, the deflection of three-phase line still exists.

4. Conclusions

Three-dimensional numerical simulation of the EM field and heat transfer were carried out to analyze the influence of asymmetric heating on the asymmetric temperature distribution and three-phase line deflection. The induction heating at the free surface of the melt is asymmetric because of the asymmetric inductor design. Using the three-dimensional induction heating calculation result, the asymmetric results of temperature distribution and flow field were calculated. The calculation results of the three-phase line show a similar trend as the experimental results. We confirmed that the spillage phenomenon can be caused by asymmetric heat generation at the three-phase line. According to the calculation result, the position of maximum heat generation is close to the three-phase line. In order to reduce the risk of spillage of melt, the position of the maximum heat generation should be designed farther from the three-phase line. For example, decreasing the length of the side slits in high-frequency inductor is beneficial to prevent spillage phenomenon.

Author Contributions

Conceptualization, H.H. and Y.M.; Funding acquisition, K.K.; Project administration, K.K.; Software, X.L. and S.N.; Supervision, K.K.; Writing—original draft, X.-F.H.; Writing—review & editing, X.-F.H. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partly supported by the New Energy and Industrial Technology Development Organization (NEDO) under the Ministry of Economy, Trade and Industry.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Menzel, R. Growth Conditions for Large Diameter FZ Si Single Crystals. Ph.D. Thesis, Technical University of Berlin, Berlin, Germany, 2013. [Google Scholar]
  2. Mühlbauer, A.; Erdmann, W.; Keller, W. Electrodynamic convection in silicon floating zones. J. Cryst. Growth 1983, 64, 529–545. [Google Scholar] [CrossRef]
  3. Mühlbauer, A.; Muiznieks, A.; Virbulis, J.; Lüdge, A.; Riemann, H. Interface shape, heat transfer and fluid flow in the floating zone growth of large silicon crystals with the needle-eye technique. J. Cryst. Growth 1995, 151, 66–79. [Google Scholar] [CrossRef]
  4. Mühlbauer, A.; Muiznieks, A.; Raming, G.; Riemann, H.; Lüdge, A. Numerical modelling of the microscopic inhomogeneities during FZ silicon growth. J. Cryst. Growth 1999, 198, 107–113. [Google Scholar] [CrossRef]
  5. Wünscher, M.; Lüdge, A.; Riemann, H. Growth angle and melt meniscus of the RF-heated floating zone in silicon crystal growth. J. Cryst. Growth 2011, 314, 43–47. [Google Scholar] [CrossRef] [Green Version]
  6. Sabanskis, A.; Virbulis, J. Simulation of the influence of gas flow on melt convection and phase boundaries in FZ silicon single crystal growth. J. Cryst. Growth 2015, 417, 51–57. [Google Scholar] [CrossRef]
  7. Raming, G.; Muiznieks, A.; Mühlbauer, A. Numerical investigation of the influence of EM-fields on fluid motion and resistivity distribution during floating-zone growth of large silicon single crystals. J. Cryst. Growth 2001, 230, 108–117. [Google Scholar] [CrossRef]
  8. Ratnieks, G.; Muiznieks, A.; Mühlbauer, A.; Raming, G. Numerical 3D study of FZ growth: Dependence on growth parameters and melt instability. J. Cryst. Growth 2001, 230, 48–56. [Google Scholar] [CrossRef]
  9. Rudevics, A.; Muiznieks, A.; Ratnieks, G.; Riemann, H. 3D modeling of the molten zone shape created by an asymmetric HF EM field during the FZ crystal growth process. Magnetohydrodynamics 2005, 41, 123–146. [Google Scholar]
  10. Sabanskis, A.; Surovovs, K.; Virbulis, J. 3D modeling of doping from the atmosphere in floating zone silicon crystal growth. J. Cryst. Growth 2017, 457, 65–71. [Google Scholar] [CrossRef]
  11. Plāte, M.; Krauze, A.; Virbulis, J. Three-dimensional modelling of thermal stress in floating zone silicon crystal growth. In Proceedings of the IOP Conference Series: Materials Science and Engineering, Vanderbijlpark, South Africa, 23–26 October 2018; Volume 355, p. 012005. [Google Scholar]
  12. Plāte, M.; Krauze, A.; Virbulis, J. Mathematical modelling of the feed rod shape in floating zone silicon crystal growth. J. Cryst. Growth 2017, 457, 85–91. [Google Scholar] [CrossRef]
  13. Han, X.-F.; Liu, X.; Nakano, S.; Harada, H.; Miyamura, Y.; Kakimoto, K. 3D Global Heat Transfer Model on Floating Zone for Silicon Single Crystal Growth. Cryst. Res. Technol. 2018, 53, 1700246. [Google Scholar] [CrossRef]
  14. Kakimoto, K.; Liu, L. Analysis of local segregation of impurities at a silicon melt–crystal interface during crystal growth in transverse magnetic field-applied Czochralski method. J. Cryst. Growth 2009, 311, 2313–2316. [Google Scholar] [CrossRef]
  15. Ratnieks, G.; Muiznieks, A.; Buligins, L.; Raming, G.; Mühlbauer, A.; Lüdge, A.; Riemann, H. Influence of the three dimensionality of the HF electromagnetic field on resistivity variations in Si single crystals during FZ growth. J. Cryst. Growth 2000, 216, 204–219. [Google Scholar] [CrossRef]
Figure 1. Illustration of simulation model for calculating high-frequency EM field in the FZ method: (a) polycrystalline feed rod, single crystal, and inductor; (b) high frequency inductor with four slits. The single crystal diameter is 100 mm. The model is constructed with unstructured grids.
Figure 1. Illustration of simulation model for calculating high-frequency EM field in the FZ method: (a) polycrystalline feed rod, single crystal, and inductor; (b) high frequency inductor with four slits. The single crystal diameter is 100 mm. The model is constructed with unstructured grids.
Crystals 10 00121 g001
Figure 2. 3D finite volume model for temperature and flow calculation in the melt and crystal. The length of the crystal is 0.3 m. The diameter of the crystal is assumed as 100 mm (4 inch). The hexahedron grids were constructed with the open-source software SALOME.
Figure 2. 3D finite volume model for temperature and flow calculation in the melt and crystal. The length of the crystal is 0.3 m. The diameter of the crystal is assumed as 100 mm (4 inch). The hexahedron grids were constructed with the open-source software SALOME.
Crystals 10 00121 g002
Figure 3. Calculated surface current density distribution at the free surface of the melt, as determined from the electromagnetic (EM) model.
Figure 3. Calculated surface current density distribution at the free surface of the melt, as determined from the electromagnetic (EM) model.
Crystals 10 00121 g003
Figure 4. Surface current density distribution along the three-phase line. Azimuthal angle of the main slit is regarded as 0°. The angles of the three side slits are −90°, 90°, and 180°, respectively.
Figure 4. Surface current density distribution along the three-phase line. Azimuthal angle of the main slit is regarded as 0°. The angles of the three side slits are −90°, 90°, and 180°, respectively.
Crystals 10 00121 g004
Figure 5. Temperature distribution at the free surface of melt. Below the three side slits, the temperature is high because of high power density.
Figure 5. Temperature distribution at the free surface of melt. Below the three side slits, the temperature is high because of high power density.
Crystals 10 00121 g005
Figure 6. Temperature distribution in the melt: (a) cross section in the X plane (in the vicinity of current supplies); (b) cross section in the Y plane (perpendicular to the main slit).
Figure 6. Temperature distribution in the melt: (a) cross section in the X plane (in the vicinity of current supplies); (b) cross section in the Y plane (perpendicular to the main slit).
Crystals 10 00121 g006
Figure 7. Flow field in the melt: (a) cross section in the X plane (in the vicinity of the current supplies); (b) cross section in the Y plane (perpendicular to the main slit). The flow separation point is indicated, where the two major vortices meet at the solid–liquid interface.
Figure 7. Flow field in the melt: (a) cross section in the X plane (in the vicinity of the current supplies); (b) cross section in the Y plane (perpendicular to the main slit). The flow separation point is indicated, where the two major vortices meet at the solid–liquid interface.
Crystals 10 00121 g007
Figure 8. Top view of 3D interface shape. The color shows the vertical position distribution of the solid–liquid interface. The position of the bottom of the interface is regarded as 0 mm. The deflection of the interface is 16 mm.
Figure 8. Top view of 3D interface shape. The color shows the vertical position distribution of the solid–liquid interface. The position of the bottom of the interface is regarded as 0 mm. The deflection of the interface is 16 mm.
Crystals 10 00121 g008
Figure 9. Comparison of the shape of the three-phase line between (a) experimental observation results [1] and (b) calculation results. The deflection of the three-phase line from the experimental observation result is smaller than that shown in the schematic. The calculated three-phase line exhibits a shape similar to that experimentally observed.
Figure 9. Comparison of the shape of the three-phase line between (a) experimental observation results [1] and (b) calculation results. The deflection of the three-phase line from the experimental observation result is smaller than that shown in the schematic. The calculated three-phase line exhibits a shape similar to that experimentally observed.
Crystals 10 00121 g009
Figure 10. Three-phase line along the azimuthal direction. The rotation direction is from right to left. Azimuthal angle of the main slit is regarded as 0°. The angles of the three side slits are −90°, 90°, and 180°, respectively.
Figure 10. Three-phase line along the azimuthal direction. The rotation direction is from right to left. Azimuthal angle of the main slit is regarded as 0°. The angles of the three side slits are −90°, 90°, and 180°, respectively.
Crystals 10 00121 g010
Table 1. Physical properties and calculation conditions in the high-frequency EM model.
Table 1. Physical properties and calculation conditions in the high-frequency EM model.
PropertiesValue and Unit
High-frequency inductor current1000 A
Current frequency3 MHz
Skin depth of melt0.26 mm
Skin depth of crystal1.3 mm
Skin depth of high frequency inductor0.037 mm
Electrical conductivity of melt1.2 × 10 6 S/m
Electrical conductivity of crystal5 × 10 4 S/m
Electrical conductivity of high frequency inductor6 × 10 7 S/m
Table 2. Process parameters and physical properties.
Table 2. Process parameters and physical properties.
PropertiesValue and Unit
Melt thermal conductivity67 W/(m∙K)
Crystal thermal conductivity 98 − (9.43 × 10 2 ) T − (2.89 × 10 5 ) T2 W/(m∙K)
Heat capacity of silicon900 J/(kg∙K)
Melt density2420 kg/m3
Crystal density2330 kg/m3
Marangoni coefficient of melt−1.0 × 10 4 N/ (m∙K)
Latent heat of silicon1.8 × 10 6 J/kg
Emissivity of melt0.27
Emissivity of crystal0.46
Crystal pulling rate3 mm/min
Crystal rotation speed5 RPM
Feed rod rotation speed−20 RPM
Single crystal diameter100 mm

Share and Cite

MDPI and ACS Style

Han, X.-F.; Liu, X.; Nakano, S.; Harada, H.; Miyamura, Y.; Kakimoto, K. 3D Numerical Analysis of the Asymmetric Three-Phase Line of Floating Zone for Silicon Crystal Growth. Crystals 2020, 10, 121. https://doi.org/10.3390/cryst10020121

AMA Style

Han X-F, Liu X, Nakano S, Harada H, Miyamura Y, Kakimoto K. 3D Numerical Analysis of the Asymmetric Three-Phase Line of Floating Zone for Silicon Crystal Growth. Crystals. 2020; 10(2):121. https://doi.org/10.3390/cryst10020121

Chicago/Turabian Style

Han, Xue-Feng, Xin Liu, Satoshi Nakano, Hirofumi Harada, Yoshiji Miyamura, and Koichi Kakimoto. 2020. "3D Numerical Analysis of the Asymmetric Three-Phase Line of Floating Zone for Silicon Crystal Growth" Crystals 10, no. 2: 121. https://doi.org/10.3390/cryst10020121

APA Style

Han, X. -F., Liu, X., Nakano, S., Harada, H., Miyamura, Y., & Kakimoto, K. (2020). 3D Numerical Analysis of the Asymmetric Three-Phase Line of Floating Zone for Silicon Crystal Growth. Crystals, 10(2), 121. https://doi.org/10.3390/cryst10020121

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop