Three-dimensional numerical simulation of temperature and flow fields in a Czochralski growth of germanium

1. Introduction

The transport processes (heat and mass transfer including melt and gas convection) problem in Czochralski crystal growth technique represents one of the biggest challenges of crystal growth modelling [6,9,24,34]. Crowly [42] developed a mathematical model of the heat transfer in the holm region in the germanium Czochralski furnace including two moving boundaries, the phase change surface and the air-liquid meniscus applying the enthalpy method. Dupret etal. [43] implemented a simulation of the temperature field and the position of the crystallization front during growth of germanium crystal using radiative heat transfer. Bogaert and Dupret [44,45] computed the time-dependent Czochralski growth of germanium crystal including heat exchange by conduction and diffuse grey radiation, in and between all the setup constituents, together with all the transient effects induced by the growth of crystal without melt convection. Bykova [46] simulated a two-dimensional time-dependent Ge crystal growth by the AHP method under microgravity conditions and showed that axial microaccelerations will have no impact on the forced melt flow along the crystallization front. Abbasoglu [47] performed a transient 3D numerical simulation to examine the roll of crystal and crucible rotations on the fluid flow and the radial segregation of silicon during the growth of Ge_x Si_1−x crystals by the Czochralski technique using microgravity conditions. Recently, Honarmandnia etal. [48,49] have carried out a 2D global simulation of an RF Czochralski apparatus for different stages of germanium crystal growth in order to analyse of thermal field, the convexity of the crystal–melt interface and thermo-elastic stresses of the grown crystal. Besides, an asymmetrical configuration leads to asymmetric and 3D thermal fields in the growth setup, which in turn, changes the crystal quality, Figure 1. For this reason, a 3D calculation is required and has to be performed.

It is usual to provide a heat shield assembly disposed above the molten source material (Si and Ge) and surrounding the ingot as it is pulled upward from the melt to shield the grown crystal against heat radiated from the crucible and the heater surrounding the crucible. Heat shield assemblies are typically constructed of graphite. Because of a relatively high emissivity of graphite, a conventional heat shield has a high ability to emit radiant energy from its surface, and so it radiates a substantial amount of heat toward the grown crystal, thereby inhibiting cooling of it. Consequently, the crystal pulling rate from the melt can be increased in crystal puller setups. On the other hand, the heat shield is refereed sometimes to a gas flow guide shield. Because it guides the inert inlet gas (usually Argon) toward the ingot to aid the cooling of the single crystal (reducing the possibility of microdefect nucleation). Another important effect of the guided gas flow is that the melt free surface is blown by it in such a manner that the gas flow takes the impurity away from the free surface regularly. Therefore, the obtained grown crystal will contain a reduced concentration of impurities and so a crystal with improved quality is achieved.

The goal of this work is to apply an appropriate 3D numerical approach of a resistance heated CZ furnace for a Ge crystal before the seeding process (i.e., including only melt and gas), and to analyze the structure of fluid flow, temperature field and global heat transfer using finite volume method (CFD FLUENT Package [50]). Our attention is specially focused on the region where the crystal is grown. It is worth to note that the obtained numerical results are interesting and quite important for the growers, which present detail information about the nonsymmetrical thermal field in the growth setup. During the growth process, from seeding to cooling crystal, a real-time control and monitoring of the thermal field is crucial task because it directly effects on the shape and quality of the grown crystal. For example, a nonsymmetrical thermal field can change the cylindrical shape (circular cross section) of the crystal to a non-cylindrical shape [51-53]. Therefore, improvement and reduction of the non-cylindrical symmetric conditions of the growth setup, and precisely controlled the growth procedure are possible using the simulation results presented here.

2. Model description

2.1. Governing equations

Figure 1 represents the schematic diagram of the CZ growth, which is used in this study. Our model covers the following key assumptions: (1) The system is not time dependent (i.e., steady state). (2) Melt and gas are incompressible Newtonian fluids keeping the Boussinesq approximation. (3) The flow is laminar. (4) Viscous dissipation is not included. (5) The surface of melt is flat, i.e., there is not any meniscus at the crucible wall. (5) The thermo-physical properties of the fluids are constant except for the density difference in the buoyancy force term. (6) Fixed temperature is set on the chamber walls.

Fig 1. Schematic model of the resisting heated Czochralski system for Ge growth.

Thus, the basic equations are;

(a) Fluid flow in the melt and gas:

                   (1)

which is the Navier-Stokes equation with Boussinesq approximation.

(b) Continuity equation:

                                                                    (2)

c) Energy equation:

   (convection and conduction in the melt and gas)                                                     (3)

   (conduction in the solid parts)              (4)

where the subscripts f and s denote fluids (gas and melt), solids (crucible, heat shield, insulation and chamber), respectively. The boundary conditions are:

a) at the solid-fluid interfaces:

                                                                   (5)

b) at the melt-gas interface:

                            (6)

                              (7)

                                                                       (8)

where the subscripts m and g denote melt and gas, respectively. The Equation (7) represents the thermal Marangoni phenomena at this interface.

c) at the outer surfaces of the chamber (the water cooled wall):

                                                              (9)

d) at the solid-gas interfaces:

                          (10)

which denotes the gas cooling combined with surface to surface heat radiation exchange from those surfaces.

e) at the gas inlet:

                                                           (11)

                                             (12)

f) at the gas outlet:

                                           (13)

    In the above equations, = (u, v, w) is the fluid velocity vector in the cylindrical coordinate system (r,φ z), p the pressure, T the temperature,  the z-directional unit vector, g the acceleration due to gravity, β the thermal expansion coefficient, µ the dynamic viscosity, α the thermal diffusivity, γ the surface tension, ε the emissivity, k the thermal conductivity, ρ the density,  the unit normal vector,  the unit tangential vector along the meridional direction and σ the Stefan-Boltzmann constant.

2.2. Numerical method

In the next step of the calculation procedure, the current problem of great interest is to obtain a very high-order accurate and efficient numerical solution for the system of governing equations in the computation domain. It is as follows:

    1. The governing equations were solved employing the finite-volume based CFD FLUENT Package with three-dimensional double precision second-order discretization. The finite volume method is a common way used in several CFD codes and has some advantages in memory usage and solution speed, especially for large problems (like our problem).

    2. Since the domain of interest has a complex geometrical shape, an unstructured grid of more than 5×10⁵ finite control volume was employed, Figure 3.

    3. The convergence is checked by monitoring residuals of momentum, continuity and energy equations which were set at 10^-5, 10^-5 and 10^-6, respectively.

2.3. The calculation conditions

All parts of the considered CZ setup have cylindrical symmetry except for gas ducts and outlets. Placing of two gas ducts in the bottom insulation and the same outlets in the chamber breaks the 2D rotational symmetry of the system. This means that the CZ setup has still a symmetry axis (z) and two perpendicular planes of symmetry (zx and zy), Figure 2. The gas ducts and outlets are located in the zx plane.

     The temperature at the central point of the melt-gas interface (i.e., the position of the seed crystal) is 1210 K, i.e., 20 K above the melting point of Ge (real condition before the seeding). It is not a boundary condition but it must be a calculation result by adaption of the heater power. The inlet gas flow rate is set to be 1 lit/min. The thermophysical properties used for our simulations are selected from the FLUENT materials database.

     In order to analyse and obtain detailed information about the effects of the asymmetrical geometry and the gas flow on the thermal conditions during the Ge growth, we have considered two computational cases, the configuration contains (a) only gas (i.e., without any melt and crystal), and (b) gas and melt, confirming to the real condition just before the Ge seed contacts the melt. We will discuss the flow and temperature fields in these two casesspecifically.

Fig 2. Two perpendicular planes of symmetry zx and zy in the growth setup.

Fig 3. The applied unstructured grid used for the calculation.

Fig 4. Temperature field (right hand side) and the vector velocity (left hand side) in the symmetry plane zx for Case a - configuration containing only gas.

3. Results and discussion

3.1. Case a: Configuration contains only gas

Figures 4 and 5 show the temperature field (right hand side) and flow arrows (left hand side) in the setup for both perpendicular planes of symmetry zx and zy. In both planes, a counter clockwise gas vortex caused by free convection (buoyancy driven) with average speed  exists and occupies the upper entire volume of the configuration (i.e., independent of the asymmetrical conditions). This flow carries the heat from the upper side of the heat shield to the middle space (location of the seed holder) and so prevents the inlet gas flow in that part. For this reason, the inlet gas flow is completely deflected from the central part toward the chamber side wall. Then it is directed toward the central part along the upper side of the heat shield, downward into the crucible, passing the area between the heat shield and the crucible, and then again downward to the space below the crucible. In that part of the system, the gas flow is guiding into two ducts placed in the bottom insulation (zx plane) and speeding it up (average speed ), and then exits the setup via the chamber outlets with average speed . For this reason, the behaviour of gas flow is completely three-dimensional and non-cylindrical symmetric in that region.  

Fig 5. Temperature field (right hand side) and the flow arrows (left hand side) in the symmetry plane zy for Case a - configuration containing only gas.

Fig 6. A three-dimensional view of the temperature field in the symmetry planes zx and zy for Case a - configuration containing only gas.

The temperature gradient in the setup is considerably affected by the orientation of gas convection, and so it is completely three-dimensional too, Figure 6. In the growth setup, the role of gas flow is quite important as other mechanisms of the heat transport phenomena (such as conduction and radiation). Therefore, the isotherms have a curvature commensurate to the structure and intensity of gas flow. The isotherms are pushed along the gas flow (proportional to its velocity) in the upper part of the system and deflected toward the bottom ducts and outlets. The temperature maximum in the setup is and is located at the upper part of the heating elements.

Fig 7. Temperature field (right hand side) and the flow arrows (left hand side) in the symmetry plane zx for Case b - configuration containing melt and gas.

Fig 8. Temperature field (right hand side) and the flow arrows (left hand side) in the symmetry plane zy for Case b - configuration containing melt and gas.

Fig 9. Temperature field (right hand side) and the flow arrows (left hand side) in the symmetry planes zx and zy in the melt for Case b.

Fig 10. A three-dimensional view of the temperature field in the symmetry planes zx and zy  for Case b - configuration containing melt and gas.

3.2 Case b: Configuration contains melt and gas

Figures 7 and 8 represent the temperature field (right hand side) and flow arrows (left hand side) in the growth system for both perpendicular planes of symmetry (zx and zy). The orientation of argon flow is completely similar to Case a, which is totally three-dimensional as has been already described. The gas average speed is which is a little higher than Case a.

     In the molten material, a counterclockwise eddy that is produced by the free convection overlapped with the thermal Marangoni flow happens and occupies its entire volume with an average speed , Figure 9. This vortex has expanded from the crucible wall to the centerline. The melt flow is started along the hot crucible sidewall (buoyant flow) and the melt-gas interface where the temperature gradient via the surface tension produces the thermal Marangoni convection [40,47].

     The temperature field in the growth setup is distinctly influenced by the gas flow similar to Case a. The influence of gas flow is more essential than conduction and radiation in the gas part. Therefore, the isotherm shapes have a deviation commensurate to the gas convection approach and flow arrangement. The gas isotherms are moved downward along the inner surface of the heat shield, upward adjacent to inner crucible sides, deflected downward in the space between the crucible and heater and deformed toward the gas ducts and outlets, Figure 10. The temperature maximum in the furnace is and is located at the lower part of the heating zone contrary to Case a. In the Ge melt which its location is in the lower portion of the crucible sidewall and close to its bottom.

     Figure 11 illustrates the temperature variation of the melt surface and its sensitivity to the gas flow above the melt surface. It shows a non-cylindrical symmetric structure, which arises from the different gas flow in the zx and zy symmetry planes. In the zx plane, the argon flow rate more accelerates and becomes stronger than zy plane and so it can cool the melt surface more efficiently in that direction. Therefore, the temperature maximum of the melt surface is located in the zy symmetry plane with more compared to the zx plane. It should be noted here that the non-symmetrical conditions of thermal field in the growth setup is quite similar in both cases considered here which is arises directly from different argon flow in the zx and zy symmetry planes. 

     It has been mentioned that any non-symmetrical condition for setup thermal field is not desirable and can lead to crucial problems related to the crystal quality. Consequently, we have to prevent or reduce these non-symmetrical thermal conditions of the growth furnace by special considerations during the system design, for example, using four symmetrical gas ducts and outlets in the presented growth setup. 

4. Conclusions and outlook

We have presented and exhibited the results of a three-dimensional global numerical calculation of the flow structure and temperature field for a non-cylindrical Ge crystal puller containing only gas (Case a), and melt and gas (Case b). From these simulations, we can conclude:

1- The presented numerical solution is only the first result of the considered non-cylindrical Ge Czochralski configuration including seed and crystal. Therefore, the obtained results are useful for the seeding process because the detailed information about the thermal field at the melt surface is a crucial factor for a successful seeding process (before, during after touching the seed to the melt).

Fig 11. Temperature variation of the melt free surface, which shows a non-rotationally symmetric distribution.

2- Although the temperature gradient is symmetrical around the center of melt surface but it is completely non-cylindrical symmetric close to the crucible wall which arises from the three-dimensional orientation of argon flow above it. The temperature difference along the y-direction is  more than x-direction, Figure 11. Therefore, this condition can influence on the control of crystal diameter and its uniformity. For this reason, further calculations including the seed and grown crystal have to be performed.