## Numerical simulation of the three-phase boundary in the Czochralski process

**Dr Alexander Raufeisen (MSc)**

**VT supervisor: Prof. Tilman Botsch**

The aim of this thesis was to develop a numerical method for fully transient three-dimensional simulations of Czochralski (Cz) crystal growth processes. Up to now, the 3D numerical simulations were limited to the quasi-steady state for the crystallization, i.e. the shape and position of the melt-crystal interface and the crystal were fixed. Only the flow and heat transfer were allowed to change in time. This assumption is only valid for well-controlled crystal growth processes for silicon during the bulk growth. However, in the starting and ending phase of the Cz process, the geometry of the crystal and the melt-crystal interface are subject to rapid changes due to the unstable conditions. Therefore, there was a need for a numerical method capable of fully time-dependent simulations of the Cz process to be able to gain knowledge of the transient behaviour of the crystallization procedure in the early and late stages. This method should include all known phenomena occurring in the process, starting from the turbulent flow and heat transfer in the melt, the heat conduction in the crystal, the thermal radiation from the free surface of the melt and the crystal, via the species transport in the melt, up to the movement of the crystallization interface, the melt free surface and the three-phase boundary (TPB), which is the contact line between the free surface of the melt, the crystal, and the surrounding atmosphere.

Such a method was developed, implemented, and applied in different simulations in the present work. As the fundament, the multi-physics CFD code FASTEST-3D was utilized. It is based on the finite-volume method using curvilinear coordinates. Presently, second-order accurate discretizations are used for all terms, and implicit three-point backward time marching (also second-order accurate) is employed. A sophisticated moving grid algorithm using the ALE method is implemented.

For the purpose of the present work, the code was modified and extended. In order to be able to compute all phenomena listed above simultaneously in one simulation, the coupling procedure available in the code had to be changed as well as the moving grid algorithm.

For the sake of this unified approach, the existing implementation of the crystallization interface tracking using the Stefan condition had to be adapted. The method for the computation of the free surface movement and, finally, the algorithm for the simulation of the dynamic TPB were added and integrated into the solution procedure. For the free surface algorithm, major changes in the calculation of the pressure had to be applied and a formulation for the computation of the curvature had to be found. New boundary conditions were introduced and an innovative way to move the free surface boundary together with the TPB was installed.

For the improvement of the grid quality after the different procedures for grid movement, several smoothing and node redistribution methods were implemented. All new or modified implementations were thoroughly tested and successfully validated with results from the literature. Furthermore, simulations were conducted to investigate the physical phenomena occurring in Cz crystal growth separately and thus to determine their effect on the crystal growth process.

First of all, a direct numerical simulation (DNS) of the turbulent flow and heat transfer of liquid silicon in an idealized cylindrical crucible was carried out in order to gain knowledge about the flow structure and to obtain reference data for the later comparison with LES computations. The cylinder with 340 mm diameter was rotating at a rate of 5 rpm and heated from the side and bottom. The temperature profile ranging from about 1702 to 1720 K was adapted from experimental data. A counterrotating disc (100 mm diameter, 20 rpm) at the top center of the melt symbolized the crystal. With these boundary conditions, the characteristic dimensionless numbers of this case are approximately Pr = 1.3×10-2, Re = 4.7×104, Gr = 2.7×109, Ra = 2.8×107, Ma = 2.8×104 and Ro = 5.3×10-1. The high Reynolds and Rayleigh number suggest that the flow is fully turbulent. The grid resolution was chosen to contain 8.4 million cells to capture almost all turbulent scales. Grid refinement towards the walls led to y+ values below 1 nearly everywhere. The resulting flow field showed a three-dimensional time-dependent turbulent structure consisting of thermal plumes rising from the bottom to the top of the melt. Large vortical structures formed around the center of the crucible and remained stable throughout the whole simulation moving around very slowly, yet rotating faster than the crucible due to the Coriolis acceleration. The averaged flow field showed a Benard cell-like structure typical for flow configurations heated from below and with a free surface like Rayleigh-Benard-Marangoni convection. The turbulent kinetic energy and the dissipation rate were highest at the corner of the crystal. Below the crystal, a small shear layer formed due to the rotation. Furthermore, the analysis of the anisotropy invariant map revealed that the flow is predominantly anisotropic. These results were published in Raufeisen, Breuer, Kumar, Botsch & Durst (2007) and Raufeisen et al. (2007, 2008a, 2008b). Since DNS is too expensive for practical applications, in the next step the large-eddy simulation technique was tested for Cz simulations.

In order to study the effects of grid resolution, SGS model (Smagorinsky’s, dynamic or no model), and discretization scheme (second-order accurate CDS or HLPA, a combination of second-order and upwind) on the result, several LES computations of the same setup with different parameters were carried out. Naturally, the simulation with the highest resolution (about 1 million CVs), showed the best results. Only small deviations from the DNS reference data were recognized. The influence of the SGS model was still low at this resolution. Reducing the number of cells stepwise by a factor of about 2, the quality of the results deteriorated as expected, until with a grid of less than 130,000 CVs, strong qualitative deviations occurred. With the reduction of the resolution, the impact of the SGS model became larger. Whereas no substantial difference between Smagorinsky’s and the dynamic model could be found, the advantage over simulations without SGS modelling were significant. Throughout the simulations, the usage of the second-order accurate CDS for the convective terms delivered results far superior to those using the HLPA discretization. Thus, it can be concluded that the grid resolution plays the most important role for the quality of the results of LES computations, followed by the discretization method. Acceptable results can still be obtained using only about 1/32 of the amount of CVs of a DNS, and using about 1/8 of the DNS resolution leads to results very close to DNS data. The speedup of LES computations using this range of grid resolutions compared to DNS lies between 20 and 200, whereas the choice of the SGS model or discretization scheme does not lead to a significant reduction of computing time. Thus, with the choice of the grid resolution, the decision between accuracy and computational effort has to be made. These results were published in Raufeisen, Breuer, Kumar, Botsch & Durst (2007) and Raufeisen et al. (2007, 2008a, 2009b).

The LES technique was directly applied to a realistic case from the literature. However, the focus was not directed towards the flow and heat transfer, but to the oxygen distribution in the melt. A new formulation including not only the dissolution of oxygen from the crucible and the incorporation into the crystal, but also considering the evaporation at the free surface, was employed. This is found to be more realistic. For validation, two test cases were computed. The first case was originally published by Togawa et al. (1995) and consisted of a crucible of 198 mm diameter with a crystal of 78 mm diameter, modelled only as a fixed disc on top of the melt. The results of the present 3D LES computations with the new formulation for the oxygen evaporation showed better agreement with the experimental data than the published results of the 2D axisymmetric simulation carried out by Togawa et al. (1995).

In addition, another test case adapted from Togawa et al. (1996) consisting of a Cz crucible of 178 mm diameter with a full-size model of the crystal (82 mm diameter) and a moving phase interface was simulated both with and without rotation of the crucible and the crystal, respectively. The comparison of the present 3D simulation results with the available 2D simulation data by (Jana et al. 2006) shows only little differences. The oxygen concentration at the free surface is of about the same value; however, the distribution is more homogeneous in the present 3D simulation, which is more realistic. However, these results, which were published in Raufeisen et al. (2006, 2007), suggest as a conclusion that the accurate computation of the flow and heat transfer (3D LES instead of 2D axisymmetric laminar) has a greater impact on the prediction of the oxygen distribution than the model of the oxygen distribution itself.

Finally, several simulations applying the algorithm for the combined movement of the free surface, the crystallization interface, and the TPB were conducted. In preliminary test simulations, the free surface was still kept fixed and the TPB was restricted to move only in horizontal direction. The shape and position of the melt-crystal interface was dynamically determined by the Stefan condition. The geometry of the crucible and crystal and the boundary conditions were taken from Togawa et al. (1996) – the same as in the second oxygen simulation case described above. First, a high pulling velocity was applied to the crystal. The simulation result showed a concave crystallization front and a decreasing crystal diameter, which represents the realistic behaviour. In the second simulation, a lower pulling velocity was applied and the temperatures at the crucible were increased. With these boundary conditions, it was expected that the crystal diameter increases and the shape of the phase interface becomes convex. This behaviour was qualitatively correctly reproduced by the simulation. However, as the TPB was still fixed in vertical direction and the crystal actually tended to grow downwards fast, the calculated growth rates at the border of the crystal-melt interface were higher than the actual growth. This error indicated that there was a need for fully dynamic computations of the TPB. These results were published in Raufeisen et al. (2008c).

The simulations with all possible degrees of freedom for the TPB in connection with a dynamic free surface and phase interface were carried out first using a small crucible model of Liu & Kakimoto (2005a, 2005b). This setup consists of a cylindrical crucible with 64 mm diameter and 28 mm height with a crystal of 32 mm diameter and 28.5 mm height. The crystal was initially located 1.5 mm above the melt level, such that a meniscus formed at the contact line between the crystal and the free surface.

Several simulations were conducted varying different parameters to investigate their influences on the growth process (see also Hirt (2010)). However, the simulations could only be run for a short period of time (about 7 s real time or 7,000 time steps, respectively) owing to the distortions of the grid and the highly sensitive moving free surface. Although grid smoothing techniques were applied, the distortions could not be completely eliminated. Thus, numerical oscillations in the pressure field were generated which could not be damped effectively such that the solution diverged. Nevertheless, some interesting effects could be investigated.

At first, the crystal pulling rate was varied between 0.3 and 3 mm/min. The results only showed a significant difference in the height of the crystal and thus the TPB caused by the different pulling speed. The shape of the phase interface was almost the same in all three simulations. This is due to the short simulation time that did not allow the slow temperature field, which has the largest influence on the growth process, to change and thus lead to a different crystal shape.

In a second study, the rotation rates of the crystal and the crucible, respectively, were varied from 1 to 5 rpm for the crucible and between -5 and -10 rpm for the crystal (counterrotation). From the results, it can be concluded that the variation of the crystal rotation did not lead to significant changes in the crystal growth process, because the overall flow structure in the melt almost stayed the same. The shear flow induced by the counterrotating crystal has a spatially rather limited influence. On the contrary, the acceleration of the crucible rotation totally altered the flow field. The faster rotation causes stronger centrifugal and Coriolis forces which can now overcome the forces due to the thermal buoyancy, such that the direction of the large convection rolls is inverted. Thus, in the center of the crucible, the flow is directed upwards instead of downwards and, consequently, at the free surface the melt is driven away from the crystal instead of being driven towards it. This leads to lower temperatures at the crystal border and below the phase interface, such that the crystal growth is accelerated.

In the third parameter study, the temperature boundary conditions at the crucible were changed. Compared to the initial case, the temperatures were lowered by 10 K or increased by 15 K, respectively. In the simulation results, the effect of this variation can be clearly recognized. In case of the lower temperatures at the crucible, the flow structure changes, similar to the case with higher crucible rotation. Due to the weaker thermal buoyancy, the centrifugal and Coriolis forces induced by the rotation gain influence, such that the flow in the center below the crystal is directed downwards. Therefore, the crystal growth is accelerated. However, clearly the lower temperatures in general promote the solidification mostly, such that also a significant increase of the crystal diameter can be noticed. The opposite is the case regarding the simulation with higher temperatures prescribed at the crucible. The flow structure remains almost the same, yet the flow velocities become higher due to the stronger thermal convection. Thus, more hot melt also with higher temperatures is transported towards the crystal, such that the crystallization process is hampered and the crystal diameter is decreased.

From this parametric study with a small Cz crucible, it can be concluded that the algorithm for the movement of the TPB is able to reproduce the real crystal growth process qualitatively correctly despite the stability issues mentioned above.

Furthermore, fully transient simulations of a large Cz crucible were conducted in order to prove the capability of the dynamic algorithm to generate realistic results even from large-scale fully turbulent Cz setups.

The original geometry of an industrial crucible with 340 mm diameter and 100 mm depth was used with a 100 mm diameter crystal, which measured about 200 mm in length. The temperature boundary conditions at the crucible were taken from an experimental measurement ranging from ca. 1702 K to 1720 K. The crucible was rotating at 5 rpm and the crystal at 20 rpm in the other direction. A numerical grid consisting of ca. 600,000 CVs was designed to meet the resolution criteria for good LES results, which were presented above.

One simulation with no crystal pulling was carried out. Due to the numerical oscillations, only 6,600 time steps (6.6 s real time) could be computed. However, the desired results were obtained: Owing to the lack of pulling, a high crystallization rate developed and the crystal significantly grew in diameter. The inhomogeneity of the temperature distribution in the melt below the crystal could also be found analogously in the crystal growth rate, i.e. in a region at the border of the crystal to which hot melt was transported, the crystallization was hindered, such that the crystal diameter did not increase as much as anywhere else and a “bump” could be recognized in the TPB and the crystal, respectively.

Another simulation with a high crystal pulling rate of 5 mm/min was conducted in order to reproduce the effect of crystal necking, i.e. a diameter decrease. As this simulation had to be started with the crystal high above the melt level, the numerical instabilities due to the grid skewness were even stronger than in the previous case. Thus, the initial conditions had to be “idealized”, i.e. a steady temperature field was assumed instead of a pre-computed turbulent field. Although the simulation diverged after 5,000 time steps (5 s real time), the expected behaviour could be recognized qualitatively in the results, i.e. the crystal diameter decreased significantly.

Therefore, it can be stated that the algorithm for the fully transient simulation of the Cz process, including the movement of the crystallization interface, the free surface, and the three-phase boundary, and all other physical phenomena, was implemented. For the first time, fully transient simulations of even large industrial Cz crucibles were conducted and showed qualitatively correct results (see Raufeisen et al. (2008d, 2009a, 2009c, 2009d, 2010)). However, numerical oscillations were present in all simulations, such that only a few seconds of process time could be computed. This issue should be addressed in future work to be able to simulate longer periods of the complete crystal growth processes.