- Email: [email protected]
Combined global 2D–local 3D modeling of the industrial Czochralski silicon crystal growth process
Journal of Crystal Growth 368 (2013) 72–80
Contents lists available at SciVerse ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Combined global 2D–local 3D modeling of the industrial Czochralski silicon crystal growth process T. Jung a, J. Seebeck a, J. Friedrich a,b,n a b
Fraunhofer Institut IISB, Schottkystr. 10, 91058 Erlangen, Germany Fraunhofer THM, Am St.-Niclas-Schacht 13, 09599 Freiberg, Germany
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 November 2012 Received in revised form 11 January 2013 Accepted 17 January 2013 Communicated by J.J. Derby Available online 29 January 2013
A global, axisymmetric thermal model of a Czochralski furnace is coupled to an external, local, 3D, timedependent flow model of the melt via the inclusion of turbulent heat fluxes, extracted from the 3D melt model, into the 2D furnace model. Boundary conditions of the 3D model are updated using results from the 2D model. In the 3D model the boundary layers are resolved by aggressive mesh refinement towards the walls, and the Large Eddy Simulation approach is used to model the turbulent flow in the melt volume on a relatively coarse mesh to minimize calculation times. It is shown that by using this approach it is possible to reproduce fairly good results from Direct Numerical Simulations obtained on much finer meshes, as well as experimental results for interface shape and oxygen concentration in the case of growth of silicon crystals with 210 mm diameter for photovoltaics by the Czochralski method. & 2013 Elsevier B.V. All rights reserved.
Keywords: A1. Computer simulation A1. Heat transfer A1. Fluid flow A1. Mass transfer A2. Czochralski method B2. Semiconducting silicon
1. Introduction Modeling of temperature field and oxygen distribution for the growth of large silicon crystals using the Czochralski method has actually advanced to a point, where it can be successfully applied to optimize the puller and the crystal growth process, recent examples of the state of the art can be found in e.g. Refs. [1–9]. Based on the results published in literature so far the following statements can be made:
(1) Accurate numerical results especially for the growth of the cylindrical crystal body can be obtained for industrial CZ silicon crystal growth when the global, quasi-stationary, axisymmetric model of the whole puller is combined with a local three-dimensional model of the melt–crystal–crucible region. (2) In the local 3D model convection and advection in the melt need to be calculated using a full time-dependent approach, steady-state RANS modeling seems not to be adequate in order to obtain accurate results [3,10]. However, it is well
n Corresponding author at: Fraunhofer Institut IISB, Schottkystr. 10, 91058 Erlangen, Germany. Tel.: þ49913176 1269; fax: þ 49913176 1280. E-mail address: [email protected] (J. Friedrich).
0022-0248/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jcrysgro.2013.01.026
established to use a Large Eddy Simulation (LES) turbulence model to reduce the number of mesh cells [11–13]. (3) In the local 3D model near-wall-treatment remains difficult, possibilities are the direct resolution of the boundary layers, the use of wall functions, or hybrid models using RANS models near the wall. There seems not to be yet a general agreement about which approach to pursue in Czochralski growth, a recent more general overview can be found e.g. in Ref. [14]. The coupling of 2D and 3D models is not trivial, especially when results shall be achieved in reasonable computation times on computer hardware available today in industry. The examples given in Refs. [1–9] use monolithic applications containing both the 2D and the 3D model. Such closed applications, however, do not offer generally both the best solution for the 2D and the 3D part of the model. This gives rise to the question if one could use independent tools and devise a tailored coupling strategy. But, e.g. Derby et al. [15] have shown that a simple iterative exchange of boundary conditions between two such models might not converge at all in some situations, and generally only very slowly. A different possibility to couple the effect of thermal convection in the melt to the global temperature field in the furnace has been proposed by Fainberg et al. [16]. Here, the melt flow and heat transfer are computed in a time-dependent, 3D model. Temperature and velocities are averaged azimuthally and in time, and a turbulent heat flux is calculated from the observed fluctuations of
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
temperature and velocity. These average values plus the turbulent heat fluxes are then imported into a Reynolds averaged form of the temperature equation inside the 2D model. It has already been shown in Ref. [16] that this approach can be faster than a simple exchange of boundary conditions. In Ref. [17] this approach has been applied to the modeling of growth of oxide crystals by the Czochralski method, but no comparisons to experimental results or discussion of grid dependence of the results have been done. In the present work, a general description of the governing equations and of the numerical procedure of the coupled global 2D–local 3D LES model is given. The combined model is implemented by using CrysMAS [18] for the global model and OpenFOAM [19] for the local 3D model. By comparison to data from Direct Numerical Simulation from Ref. [20] it will be shown that convective heat transport in the melt volume can be predicted on a relatively coarse mesh using the local 3D LES model, when boundary layers are resolved by an aggressive refinement of the mesh in a direction perpendicular to the wall, without compromising numerical stability, and without adding additional modeling parameters or assumptions as required e.g. by wall functions. Then, the coupled 2D–3D model is used to compute the heat and oxygen transfer during growth of silicon crystals with 210 mm diameter in an industrial Czochralski puller. The interface shapes and the oxygen concentration obtained from the 2D–3D model will be compared to the experimental data. Finally, the results obtained with the 2D–3D model will be discussed with respect to convergence behavior, accuracy and computational effort.
2. Description of the combined 2D–3D model 2.1. Mathematical description 2.1.1. 2D model The basics of the present 2D–3D-coupling scheme are presented in Refs. [16,17]. The starting point is the global 2D model, which computes the temperature field in the whole furnace including the shape of the solid–liquid interface by adjusting the heating power in such way, that the temperature at the prescribed position of triple point matches the melting point [21]. The solidification of liquid silicon in 2D is modeled by using an interface tracking algorithm, which considers the latent heat sources according to the given growth velocity which is assumed to be equal to the pulling speed Vp. In the melt region the influence of the unsteady 3D melt motion is taken into account using the Reynold’s averaged form of the enthalpy equation: ! rUCU! u UrT g ¼ rUðlg rT g ÞrU a ð1Þ
r—density of liquid silicon; C—heat capacity of liquid silicon; ! u —temporally and azimuthally averaged velocity field; Tg—temperature field (axisymmetric); lg—effective thermal ! conductivity; a —turbulent heat fluxes. ! The temporally and azimuthally averaged velocity field u , ! effective thermal conductivity lg and turbulent heat fluxes a are obtained from the local 3D model as it will be explained further below. 2.1.2. Data transfer from 2D to 3D From the global 2D temperature field the relevant boundary conditions along the periphery of the local 3D model are extracted, i.e. all data transferred from 2D to 3D are axisymmetric:
(i) geometry of the melt domain including the shape of solid– liquid interface;
73
(ii) at the crucible–melt interface the temperature distribution Tc (r, z); (iii) at the solid–liquid interface the melting isotherm T0; (iv) at the free melt surface the heat flux Q2d (r, z). Q2d (r, z) considers the heat flux which is absorbed at the free melt surface from the surrounding environment through radiation and conduction. Q2d enters the temperature boundary condition at the free melt surface of the local 3D model in the following way: @T leff ! ¼ esT 4 Q 2d ð2Þ @n ! where n —surface normal vector (pointed outwards from melt), e—emissivity of the melt, s—Stefan–Boltzmann constant, Q2d— absorbed heat flux.
2.1.3. 3D model Using the geometrical information about the melt region including the present shape of the solid–liquid interface obtained from the global 2D model, the geometry and mesh of the local 3D computational domain is generated. The transient 3D calculation of the convective heat transport is performed with fixed geometry and mesh. The boundary conditions for the velocity field of the local 3D model are a fixed angular velocity at the solid–liquid interface, and at the crucible wall to take into account crystal and crucible rotation, as well as Marangoni shear stress at the free melt surface. In order to consider the turbulent behavior of the melt the LES approach is applied using the Smagorinsky sub-grid model with van Driest wall damping [22]. The resulting governing equations for mass, momentum and energy of the local 3D model are as follows: !
rU U ¼ 0
ð3Þ
! ! @ U ! ! ! þ U Ur U ¼ rP þ rU 2meff S þ rbðT 0 T Þ g rU @t
ð4Þ
@T ! þ U Ur T ¼ rUðleff rT Þ ð5Þ @t ! where U —velocity, P—pressure, meff—effective viscosity of the melt, S—strain rate tensor, which is defined in vector form as ! ! follows S ¼ ðr U þ ðr U ÞT Þ=2, b—thermal expansion coefficient, ! T0—reference temperature (melting point of the melt), g —gravity, C—heat capacity of the melt, leff—effective thermal conductivity of the melt. In the formulation of the governing equations (Eqs. (3)–(5)) an effective melt viscosity meff and an effective thermal conductivity of the melt leff account for the influence of the sub-grid scale turbulence according to the following relations: pffiffiffiffiffiffiffiffiffiffiffiffi mt ¼ rUðC s dÞ2 U 2USUS ð6Þ
rUC
meff ¼ m þ mt leff ¼ CU
m Pr
ð7Þ
þ
mt Prt
ð8Þ
with m—laminar viscosity, mt—turbulent viscosity, calculated according to Smagorinsky–Lilly model with constant Cs ¼0.167 and d—cubic-root of the mesh cell volume. The laminar Prandtl number is Pr ¼0.011479 and the influence of the sub-grid scale turbulence on heat transfer was modeled by using the turbulent Prandtl number Prt ¼0.9 as e.g. in Ref. [11].
74
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
When these relations (Eqs. (11)–(13)) are introduced into the enthalpy equation (Eq. (1)) and the averaging operation (Eq. (9)) is applied one gets from the local 3D model the following equation for the axisymmetric temperature field at stationary state in the melt: !
!0
rUCU/ U Sr/TS ¼ rUð/leff Sr/TSÞ þ rUðrUCU/ U T 0 S/l0 rT 0 SÞ ð14Þ The local 3D (Eq. (14)) and the global 2D model (Eq. (1)) must yield the same temperature fields, i.e. oT4 ¼ Tg. Therefore, the following definitions of the variables of the 2D model using the corresponding values of the 3D model can be made which results in a closed mathematical formulation of the 2D/3D problem: ! ! u ¼ /U S
ð15Þ
lg ¼ /leff S
ð16Þ
!0 ! 0 a ¼ ðrUCU/ U T 0 S/l rT 0 SÞ
ð17Þ !0 0 The temporally and azimuthally averaged values / U T S and 0 /l rT 0 S in (Eq. (17)) are evaluated in local 3D model according to the following relations: !0 ! ! / U rT 0 S ¼ / U rTS/ U SU/rTS 0
/l rT 0 S ¼ /leff rTS/leff SU/rTS
Fig. 1. Flow chart for the coupled global 2D and local 3D model.
2.1.4. Data transfer from 3D to 2D In order to transfer the results obtained from the local 3D model to the global 2D model, the values and cross correlation terms of temperature, velocity and thermal conductivity fields computed in the local 3D model are averaged in time and in the azimuthal direction. The averaging operation for scalar (V) and ! vector values ( V ) is defined according to the following equations: Z 2p Z t0 þ t 1 /VS ¼ VUdtUdj ð9Þ 2pUt 0 t0 "Z # 2p Z t 0 þ t ! 1 ! ! V U e r UdtUdj U e r 2pUt 0 t0 ""Z # #! 2p Z t 0 þ t ! ! ! þ V U e z UdtUdj U e z
ð18Þ ð19Þ
2.1.5. Model for the oxygen transport in the 2D–3D model In addition to the convective heat transport the oxygen transport is calculated by the coupled 2D–3D approach. For that purpose the equation for oxygen transport in the melt, which is analogous to Eq. (5), is solved in the local 3D model. Thereby, the oxygen concentration in the melt is coupled to the transport of SiO in the gas phase of the axisymmetric global model via a flux balance and balancing of the SiO vapor pressure against the azimuthally averaged oxygen concentration at the free melt surface, using the same boundary conditions and material properties as given in Ref. [23]. In the current implementation this coupling of the oxygen concentration in the melt to the SiO concentration in the gas is based on an explicit iterative exchange of the boundary conditions at the free melt surface.
! /V S ¼
0
2.2. Implementation and iteration procedure ð10Þ
t0
where t—time interval for averaging, t0—start instant for the ! ! averaging in time, e r , e z —radial and axial base vectors of the cylindrical reference frame, j—azimuthal angle. From this averaging procedure the distribution of the con! ! vective velocity u , the turbulent heat fluxes a and the effective thermal conductivity lg are determined which are needed in Eq. (1) of the global 2D model. ! Temperature T, velocity U and thermal conductivity leff of the local 3D model can be formally expressed as a sum of a mean ! temporally and azimuthally averaged value oV 4 or o V 4 and ! its corresponding value of the instantaneous fluctuations V’ or V ’ respectively: T ¼ /TS þ T 0
ð11Þ
! ! !0 U ¼ / U Sþ U
ð12Þ
leff ¼ /leff Sþ l0
ð13Þ
Fig. 1 shows a flow chart of the coupled 2D–3D calculations. The coupling is initialized by generating the 2D mesh of the whole CZ furnace in CrysMAS [18]. After the global 2D simulation, the corresponding axisymmetric shape of the melt is used to generate the 3D mesh using Gambit [24]. For the 1st iteration the influence of the unsteady 3D melt motion is not considered directly, instead an initial guess for the effective thermal conductivity lg in the melt is made. During the further iterations the converged axisymmetric solution of the global 2D model is used to update the 3D geometry of the melt region, based on the actual shape of the solid–liquid interface, and to extract the boundary conditions for the local 3D model. For that purpose corresponding interfaces are developed and used to exchange the data between CrysMAS and OpenFOAM. Furthermore, routines are implemented in the local 3D OpenFOAM model to actualize the 3D mesh according to the current shape of the solid–liquid interface. As the changes of the interface shape during the iterations can be considered as relatively small a stretching approach for the mesh points is used instead of the regeneration of the whole 3D mesh.
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
75
Fig. 2. Time- and azimuthally averaged temperature T (a), stream function (b), radial (c) and axial component and (d) of turbulent heat fluxes. Left: this work, LES model with 108,000 cells and aggressive wall refinement; right: DNS model with 8 million cells (from Ref. [7]).
with surface element ds:
Table 1 Used material properties of silicon. Material Property Kinematic viscosity n Thermal expansion b Melting point Tm Thermal conductivity (solid) ls Thermal conductivity (liquid) ll Emissivity (solid) es Emissivity (liquid) el Latent heat L Density r Heat capacity C
Value 3 10 7 m2/s 1.4 10 4 K 1 1683 K 96017/T1.149 W/(m K) 60.26 W/(m K) 0.7 0.3 1.8 106 J/kg 2520 kg/m3 915 J/(kg K)
Next, the transient 3D flow calculation is carried out in OpenFOAM until the system on average reaches a stable state. The solver used is based on the ‘‘buoyantBoussinesqPimpleFoam’’ solver provided with OpenFoam [19]. All calculations are using central difference interpolation schemes, in a reference system rotating with the crucible. Then, after the stable state is reached, the timely and azimuthally averaged fields (Eqs. (15)–(17)) are transferred from the local 3D OpenFOAM to the global 2D CrysMAS model. After interpolation of the 3D results onto the mesh of the 2D model the ! averaged velocities u are corrected by addition of the gradient of !0 a scalar field p to yield a corrected velocity field u in order to ensure conservation of mass and energy: !0 ! u ¼ u þ rp Z
!0 u Uds ¼ 0-
Z
rpUds ¼
Z
! u Uds
ð20Þ
This yields an equation system for p which is solved using the finite volume method implemented in CrysMAS, i.e. forming an equation system by summing up the boundary fluxes. The correction step also includes forcing normal components of averaged velocity and turbulent heat fluxes at axis of symmetry and outer walls to zero. No efforts have been made to enforce zero radial gradients of velocity at the axis of symmetry. On one hand, this results difficult with a simple collocated, equal order scheme on an unstructured triangular mesh for velocities and the pressure-like scalar field. On the other hand, the results shown later indicate that after long enough averaging times this results naturally to a sufficient degree. !0 With the corrected velocity field u the energy equation (Eq. (1)) is solved in the global 2D model. In order to ensure numerical stability of the coupled system, the convective term of Eq. (1) is ramped up during the first few iterations. The update of the boundary conditions for the 3D flow model from the 2D model is done using an under-relaxation factor of typical 0.25. This coupling procedure described above is repeated until the 2D/3D thermal fields and the shape of the phase interface are converged.
3. Test cases 3.1. Comparison to benchmark data from Direct Numerical Simulation As time averaged data from time dependent calculations are required for several 2D–3D iterations the flow calculation procedure in the local 3D model must be as efficient as possible while
76
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
maintaining a reasonable degree of precision. The main problem is posed by the necessity to resolve or model the boundary layers for temperature, momentum, and oxygen. In this work large efforts were made to really resolve the boundary layers down to the laminar region, corresponding to a dimensionless wall distance y þ of the cell nearest to the wall of less than 5 according to von-Karman’s law [25]. This was achieved by a radical refinement of the mesh perpendicular to the wall, and accepting aspect ratios of near-wall-cells up to almost 1000. While the average cell size in the volume of the melt is 5 mm, it goes down to 3 mm at the walls, details are shown later for the example of the realistic Czochralski crucible. In order to investigate the validity of this approach the results obtained from the local 3D LES model with radical mesh refinement close to the walls are compared to benchmark data from a Direct Numerical Simulation (DNS) with approximately 8,000,000 cells for a model Czochralski configuration published by Raufeisen et al. [20]. Fig. 2 shows the time and azimuthally averaged temperature, stream function, and turbulent heat fluxes obtained from the present 3D model for a mesh within total 107,000 cells and the corresponding DNS data [20]. The intention was not to reproduce these results exactly, but to find a suitable compromise between precision and calculation time. The results do not match exactly the DNS data, but all important features are reproduced, and the absolute numbers match to within about 20%. It is believed that this agreement is good enough for an engineering model, considering the extremely low number of cells resulting in calculation speeds high enough to calculate around 600 s of real-time in typically 24 h on the used AMD Opteron 2352 cluster with time steps between 0.02 and 0.05 s. Test runs with smaller time steps did not show any significant difference. 3.2. Test for a realistic Cz crucible without 2D–3D-iteration Next, a simplified case consisting in 2D and 3D only of a 65 kg heavy silicon melt with fixed solid–liquid interface in a crucible with a diameter of 580 mm is considered in order to test the implementation of the coupled 2D–3D model. The thermal boundary conditions for the reduced 2D and 3D models were obtained from a typical, global 2D CZ furnace model containing the crucible with this melt. The crystal and crucible rotation rates were ocrys ¼ 10 rpm and ocruc ¼10 rpm, respectively. The material properties of silicon used in the simulations are summarized in Table 1. The local 3D hydrodynamic LES simulation is carried out on a mesh with 108,000 elements for more than 1000 s real time while the mesh in the 2D simulation has around 33,000 elements in the melt region. Fig. 3 shows the temporal behavior of the average temperature, flow velocity and turbulent heat flux at a monitor point in the center of the melt. It is obvious that after an initial transient time of around 200 s the system reaches a stable state, i.e. there are only minor variations of the averaged field values. Therefore, a total real time of 600 s is considered to be sufficient in the local 3D model because the system is stable for at least 400 s which gives enough data for the averaging procedure to calculate the timely and azimuthally averaged variables. Fig. 4 shows the comparison of the temporally and azimuthally averaged temperature results from the 2D and the 3D model. The average temperature distribution on the right side is obtained from the local transient 3D LES simulation. The corresponding distribution on the left side is calculated in the quasi-stationary 2D model according to Eq. (1) using the averaged values from the 3D model. It is obvious that there is a large discrepancy between the 2D and 3D ! results when turbulent heat fluxes a are not considered in the 2D ! model (Fig. 4a). When the turbulent heat fluxes a are taken into account while solving Eq. (1) in the 2D model, a quantitative agreement between temperature fields of the 2D and 3D models is
Fig. 3. Time averaged temperature o T4 (a), absolute value of velocity oU 4 (b) and turbulent heat flux o a4 and (c) at a monitor point in the melt versus computed real time.
found (Fig. 4b). This result highlights the importance of the turbulent ! heat fluxes a when exchanging data between the 2D and 3D models. Almost no further improvement (Fig. 4c) is achieved by including the terms resulting from the fluctuations of turbulent conductivity (Eq. (17)). The deviation is still less than 2 K. The remaining small discrepancies between the 2D and 3D model are most likely results of limited averaging times as well as of errors from the interpolation from 3D to 2D mesh, and possibly of the consequent flow correction step described in Eq. (20). However, this existing deviation is considered to be negligible compared to all other assumptions and errors which are made when setting up a global model for a Czochralski puller. 3.3. Mesh dependence of the 3D model The dependence of results of the 3D model on the mesh refinement close to the walls of the melt domain has been analyzed for a 24 in. crucible containing 90 kg of silicon. The thermal boundary conditions were extracted from a global 2D furnace model, and certain crystal and crucible rotation rate were used. The mesh refinement at the walls was done using the ‘‘refineWallLayer’’ tool provided with OpenFoam [19], which allows the refinement of selected cells in specified directions. In order to estimate the effect of the mesh refinement an initial mesh with 79,000 cells was generated and different levels of the wall refinement were set. Fig. 5 shows the resulting maximal and averaged dimensionless wall distances yþ at the growth interface as a function of the averaged wall distance of the center of the cells at the boundary for the different levels of the wall
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
77
Fig. 4. Temperature field obtained from global 2D model (left) and averaged temperature field obtained from local 3D model (right); (a) without considering turbulent heat fluxes in 2D; (b) with turbulent heat fluxes in 2D, but without turbulent heat conductivity (i.e. lg ¼ llam); and (c) with turbulent heat fluxes and turbulent heat conductivity (i.e. lg ¼ o leff 4) in 2D.
Fig. 5. Average and maximal dimensionless wall distance y þ at growth interface versus real distance y of the cell center to the wall at the center of the interface for the coarse volume mesh with initial 79,000 cells and different levels of wall refinement (lines) and for the fine mesh with in total 800,000 cells (symbols).
until yþ is everywhere less than 5. In this case the total number of cells is about 800,000. In this case the resulting dimensionless wall distance yþ is slightly smaller than for the coarse mesh with 108,000 cells and with radical wall refinement only. For the coarse volume mesh with the different refinement levels at the walls as well as for the fine volume mesh with the final refinement level the convective heat and oxygen transport in the silicon melt was computed with the local 3D time-dependent LES model assuming that the oxygen concentration is zero at the free melt surface. Fig. 6 shows the timely averaged, integrated heat flux across the growth interface as well as the maximal and averaged oxygen concentration at the interface for the different meshes. It can be seen, that when the dimensionless wall distance yþ becomes less than 5, the change of the heat flux at the interface from the coarse to the fine volume mesh is only about 10%, and the change of the oxygen concentration only about 15%. Therefore, the overall mesh dependent error on the coarse mesh is assumed to be below 25%.
4. Application of the coupled 2D–3D model to an industrial Czochralski process refinement of the initial mesh with 79,000 cells. It is obvious that a cell size of less than 10 mm is needed in order to fulfill the criterion yþ o5. In addition to the wall refinement of the coarse mesh, a fine volume mesh was generated where the initial mesh was refined by dividing each cell into 8 cells, and again refining the wall layers
The combined 2D–3D model has been applied for an industrial Czochralski puller used to grow silicon crystals with diameter of approximately 210 mm for photovoltaic applications from a silica crucible with 24in. diameter. In this case the Grashof number Gr ¼ gUR3 UDTUb =n2 , is¼2.2 1010 with crucible radius R¼12 in.,
78
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
DT¼50 K as a typical maximal temperature difference from crystal to crucible and material properties n and b as given in Table 1. 4.1. Comparison of simulated solid–liquid interface shape with experimental data Experimental data had been provided where the pulling speed Vp was varied during the growth of the crystal body while keeping the other process parameters such as crystal and crucible rotation constant. Lateral Photovoltage Scanning technique [26] has been used to visualize the interface shape on longitudinal cuts. The 2D–3D model was applied to two cases: case A corresponds to a growth stage with a melt weight of around 180 kg at a pulling speed VA p while for case B the melt weight is 130 kg at a pulling speed VBp 4VA p , the actual values being confidential. For both cases, the numerical meshes consist of 40,000 elements in the melt region of the global 2D model and 150,000 elements in the local 3D model. Fig. 7 shows the solid–liquid interface shape for both cases together with the calculated phase boundaries as bold red lines. In view that
there has not been any model adjustment at all, the input consisted merely in a CAD drawing of the puller as well as of process and material parameters, a good correlation between experiment and simulation was obtained. The convergence behavior of the 2D–3D-coupling scheme is shown in Fig. 8 for a case, where the growth rate Vp has been reduced by 50%, starting from a converged solution. Heater power, interface position, maximal crucible temperature and the average temperature difference in the melt between the 3D and 2D model are basically converged after five 2D–3D iterations. 4.2. Comparison of simulated oxygen concentrations to experimental data Apart from the results shown before, the measured interstitial oxygen concentrations had been available for a different growth run with fixed pulling rate Vp in the same puller as in the section above. Calculations including oxygen transport have been performed for three different lengths of the cylindrical crystal body of 2 cm, 102 cm,
Fig. 6. Integrated heat flux Q over growth interface (a) and oxygen concentration [O] at center of growth interface (b) as function of mesh refinement, expressed by the averaged dimensionless wall distance yþ at the interface for the coarse volume mesh with initial 79,000 cells and different levels of wall refinement and for the fine mesh with in total 800,000 cells. The maximum values and the values averaged over the whole area of the interface are shown.
Fig. 7. Comparison of the calculated interface shapes (bold lines) with experimental data (background) for an industrial CZ system at two growth stages; top — low B growth rate VA p , 180 kg melt weight, bottom — high growth rate Vp, 130 kg melt weight. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
79
Fig. 8. Convergence behavior of the 2D–3D-coupling scheme: heater power (a), axial position of solid–liquid interface at the symmetry axis (b), maximum crucible temperature (c), and average temperature deviation in the melt between 2D and 3D model (d).
and 168 cm. These calculations have been done applying two kinds of boundary conditions for the oxygen concentration at the free melt surface. In case O1 the oxygen concentration at the free melt surface was set to zero, which means that there is no coupling of the oxygen transport between the local 3D model of the melt and the global axisymmetric model of the gas atmosphere. In case O2 the oxygen transport in the melt was coupled to the transport in the argon stream, as described in Section 2. The meshes for the three growth stages consist again of typically 100,000 to 150,000 cells with radical refinement at the walls until yþ becomes less than 5. The results of the simulations for both cases O1 and O2 are shown in Fig. 9. The oxygen concentration [O] measured at the center of wafers cut at the three different crystal lengths decreases from the top to the tail end of the crystal. This tendency can only be reproduced for case O2 when the oxygen transport in the melt is fully coupled to the SiO transport in the argon stream. Without coupling the oxygen concentration increases at the tail end again. Taking into account that there has not been any model adjustment at all, the measured and computed values agree again reasonably well.
5. Discussion The examples which were shown in Section 4 highlight that the presented combination of global 2D and local 3D model is able to deliver accurate results for industrial silicon Czochralski crystal growth. The solution of the global 2D problem is not time critical. Therefore, the overall computation time is strongly determined by the speed of the local, transient 3D model and how many iterations are needed to get a converged solution. Taking the example from above one iteration which comprises the computation of in total 600 s real time with 150,000 elements in the local 3D model and one quasi-stationary solution of the global 2D model takes around 1 day, where the cost of the included 2D calculation is negligible. For a converged solution five iterations are typically needed, i.e. the
Fig. 9. Comparison of measured oxygen concentration at the crystals axis versus grown length of the cylindrical crystal body with calculated values for case O1 (without coupling to the gas atmosphere) and case O2 (with coupling to the gas atmosphere).
overall computation time for one set of process parameters is five days. It is not possible to compare this computation time directly with other methods used for combined 2D–3D models because there were no data published. A significant reduction of the time for solving the problem can be expected either when for special cases the local 3D simulations can be significantly accelerated by using for example the approach published in Ref. [27] for the calculation of the influence of horizontal magnetic fields on the melt flow. But, it is also believed that the performance of the implementation of local 3D model on the basis of the popular open source framework OpenFOAM [19] used in present work can be improved further by employing improved parallelization on computer clusters and platforms with high graphical processing power or by utilizing more efficient calculation algorithms.
80
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
6. Summary and conclusions A previously presented idea for modeling of Czochralski processes by using a coupled global 2D–local 3D model was applied to describe the transport phenomena occurring in large industrial scale crystal growth systems. Using a coarse mesh in the melt volume, but a radically refined mesh at the walls, the turbulent 3D motion in the silicon melt was modeled with the LES approach by using the axisymmetric phase interface and boundary conditions from global 2D model. The influence of the unsteady 3D melt flow was accounted by using averaged 3D fields in the temperature equation for the melt region of the global 2D model. The combined model was implemented in CrysMAS and OpenFOAM codes. The model was used to simulate the industrial Czochralski process for growing silicon crystals with 210 mm from a crucible with 2400 diameter. A good quantitative correlation between the measured and simulated shape of the solid–liquid interface for two different pulling speeds was found. The computed and measured oxygen concentration in the crystal for a different growth run agreed also fairly well, when the oxygen transport in the melt is coupled to the SiO transport in the argon stream.
[2] [3] [4] [5]
[6] [7] [8] [9] [10] [11] [12] [13]
[14] [15] [16] [17]
Acknowledgment The authors would like to acknowledge our former co-worker Dr. Andis Rudevics who has made significant contributions to the implementation of the coupling procedure as well as Bosch Solar Energy AG for providing the experimental data. This work, which was funded by the German Ministry of Education and Research (BMBF) was carried out within the ‘‘CZSil’’ project (FKZ03SF0379C) which is part of the German cluster of excellence ‘‘Solarvalley Mitteldeutschland’’.
[18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
References [1] W. Su, R. Zuo, K. Mazaev, V. Kalaev, Journal of Crystal Growth 312 (2010) 495–501.
B. Gao, K. Kakimoto, Journal of Crystal Growth 312 (2010) 2972–2976. L. Liu, K. Kakimoto, Journal of Heat and Mass Transfer (2005) 4481–4491. L. Liu, K. Kakimoto, Journal of Crystal Growth 275 (2005) 1521–1526. V.V. Kalaev, D.P. Lukanin, V.A. Zabelin, N. Makarov, J. Virbulis, E. Dornberger, W. von Ammon, Materials Science in Semiconductor Processing 5 (2003) 369–373. O.V. Smirnova, V.V. Kalaev, N. Makarov, Ch. Frank-Rotsch, M. Neubert, P. Rudolph, Journal of Crystal Growth 266 (2004) 67–73. O.V. Smirnova, V.V. Kalaev, Yu.N. Makarov, N.V. Abrosimov, H. Riemann, Journal of Crystal Growth 266 (2004) 74–80. O.V. Smirnova, V.V. Kalaev, Yu.N. Makarov, N.V. Abrosimov, H. Riemann, Journal of Crystal Growth 287 (2006) 281–286. O.V. Smirnova, V.V. Kalaev, Yu.N. Makarov, N.V. Abrosimov, H. Riemann, V.N. Kurlov, Journal of Crystal Growth 303 (2007) 141–145. ¨ A. Seidl, G. Muller, E. Dornberger, E. Tomzig, B. Rexer, W.v. Ammon, Proceedings of the Electrochemical Society 98-1 (1998) 417–428. D.P. Lukanin, V.V. Kalaev, Yu.N. Makarov, T. Wetzel, J. Virbulis, W. von Ammon, Journal of Crystal Growth 266 (2004) 20–27. A. Krauze, N. Jekabsons, A. Muiznieks, A. Sabanskis, U. Lacis, Journal of Crystal Growth 312 (2010) 3225–3234. I.Yu. Evstratov, V.V. Kalaev, A.I. Zhmakin, Yu.N. Makarov, A.G. Abramov, N.G. Ivanov, E.M. Smirnov, E. Dornberger, J. Virbulis, E. Tomzig, W. von Ammon, Journal of Crystal Growth 230 (2001) 22–29. J. Larsson, S. Kawai, Annual Research Briefs, Center for Turbulence Research, 2010 pp. 39–46. J.J. Derby, L. Lun, A. Yeckel, Journal of Crystal Growth 303 (2007) 114–123. ¨ J. Fainberg, D. Vizman, J. Friedrich, G. Muller, Journal of Crystal Growth 303 (2007) 124–134. T. Tsukada, K.I. Sugioka, C.J. Jing, M. Kobayashi, Journal of Crystal Growth 312 (2010) 997–1004. ¨ M. Kurz, A. Pusztai, G. Muller, Journal of Crystal Growth 198–199 (1999) 101–106. OpenFOAM (R), /www.openfoam.comS. A. Raufeisen, M. Breuer, T. Botsch, A. Delgado, International Journal of Heat and Mass Transfer 51 (2008) 6219–6234. M. Kurz, Ph.D. Thesis, University Erlangen-Nueremberg, 1998. E.R. van Driest, Journal of the Aeronautical Sciences 23 (1956) 1007–1011. A.D. Smirnov, V.V. Kalaev, Journal of Crystal Growth 310 (2008) 2970–2976. Gambit (R), /www.ansys.comS. T. von Karman, Nachrichten von der Gesellschaft der Wissenschaften zu ¨ Gottingen Fachgruppe 1 (Mathematik) 5 (1930) 58–76. ¨ ¨ N.V. Abrosimov, A. Ludge, H. Riemann, W. Schroder, Journal of Crystal Growth 237 (2002) 356–360. Y. Collet, O. Magotte, N. Van den Bogaert, R. Rolinsky, F. Loix, M. Jacot, V. Regnier, J. Marchal, F. Dupret, Journal of Crystal Growth 360 (2012) 18–24.
Contents lists available at SciVerse ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Combined global 2D–local 3D modeling of the industrial Czochralski silicon crystal growth process T. Jung a, J. Seebeck a, J. Friedrich a,b,n a b
Fraunhofer Institut IISB, Schottkystr. 10, 91058 Erlangen, Germany Fraunhofer THM, Am St.-Niclas-Schacht 13, 09599 Freiberg, Germany
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 November 2012 Received in revised form 11 January 2013 Accepted 17 January 2013 Communicated by J.J. Derby Available online 29 January 2013
A global, axisymmetric thermal model of a Czochralski furnace is coupled to an external, local, 3D, timedependent flow model of the melt via the inclusion of turbulent heat fluxes, extracted from the 3D melt model, into the 2D furnace model. Boundary conditions of the 3D model are updated using results from the 2D model. In the 3D model the boundary layers are resolved by aggressive mesh refinement towards the walls, and the Large Eddy Simulation approach is used to model the turbulent flow in the melt volume on a relatively coarse mesh to minimize calculation times. It is shown that by using this approach it is possible to reproduce fairly good results from Direct Numerical Simulations obtained on much finer meshes, as well as experimental results for interface shape and oxygen concentration in the case of growth of silicon crystals with 210 mm diameter for photovoltaics by the Czochralski method. & 2013 Elsevier B.V. All rights reserved.
Keywords: A1. Computer simulation A1. Heat transfer A1. Fluid flow A1. Mass transfer A2. Czochralski method B2. Semiconducting silicon
1. Introduction Modeling of temperature field and oxygen distribution for the growth of large silicon crystals using the Czochralski method has actually advanced to a point, where it can be successfully applied to optimize the puller and the crystal growth process, recent examples of the state of the art can be found in e.g. Refs. [1–9]. Based on the results published in literature so far the following statements can be made:
(1) Accurate numerical results especially for the growth of the cylindrical crystal body can be obtained for industrial CZ silicon crystal growth when the global, quasi-stationary, axisymmetric model of the whole puller is combined with a local three-dimensional model of the melt–crystal–crucible region. (2) In the local 3D model convection and advection in the melt need to be calculated using a full time-dependent approach, steady-state RANS modeling seems not to be adequate in order to obtain accurate results [3,10]. However, it is well
n Corresponding author at: Fraunhofer Institut IISB, Schottkystr. 10, 91058 Erlangen, Germany. Tel.: þ49913176 1269; fax: þ 49913176 1280. E-mail address: [email protected] (J. Friedrich).
0022-0248/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jcrysgro.2013.01.026
established to use a Large Eddy Simulation (LES) turbulence model to reduce the number of mesh cells [11–13]. (3) In the local 3D model near-wall-treatment remains difficult, possibilities are the direct resolution of the boundary layers, the use of wall functions, or hybrid models using RANS models near the wall. There seems not to be yet a general agreement about which approach to pursue in Czochralski growth, a recent more general overview can be found e.g. in Ref. [14]. The coupling of 2D and 3D models is not trivial, especially when results shall be achieved in reasonable computation times on computer hardware available today in industry. The examples given in Refs. [1–9] use monolithic applications containing both the 2D and the 3D model. Such closed applications, however, do not offer generally both the best solution for the 2D and the 3D part of the model. This gives rise to the question if one could use independent tools and devise a tailored coupling strategy. But, e.g. Derby et al. [15] have shown that a simple iterative exchange of boundary conditions between two such models might not converge at all in some situations, and generally only very slowly. A different possibility to couple the effect of thermal convection in the melt to the global temperature field in the furnace has been proposed by Fainberg et al. [16]. Here, the melt flow and heat transfer are computed in a time-dependent, 3D model. Temperature and velocities are averaged azimuthally and in time, and a turbulent heat flux is calculated from the observed fluctuations of
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
temperature and velocity. These average values plus the turbulent heat fluxes are then imported into a Reynolds averaged form of the temperature equation inside the 2D model. It has already been shown in Ref. [16] that this approach can be faster than a simple exchange of boundary conditions. In Ref. [17] this approach has been applied to the modeling of growth of oxide crystals by the Czochralski method, but no comparisons to experimental results or discussion of grid dependence of the results have been done. In the present work, a general description of the governing equations and of the numerical procedure of the coupled global 2D–local 3D LES model is given. The combined model is implemented by using CrysMAS [18] for the global model and OpenFOAM [19] for the local 3D model. By comparison to data from Direct Numerical Simulation from Ref. [20] it will be shown that convective heat transport in the melt volume can be predicted on a relatively coarse mesh using the local 3D LES model, when boundary layers are resolved by an aggressive refinement of the mesh in a direction perpendicular to the wall, without compromising numerical stability, and without adding additional modeling parameters or assumptions as required e.g. by wall functions. Then, the coupled 2D–3D model is used to compute the heat and oxygen transfer during growth of silicon crystals with 210 mm diameter in an industrial Czochralski puller. The interface shapes and the oxygen concentration obtained from the 2D–3D model will be compared to the experimental data. Finally, the results obtained with the 2D–3D model will be discussed with respect to convergence behavior, accuracy and computational effort.
2. Description of the combined 2D–3D model 2.1. Mathematical description 2.1.1. 2D model The basics of the present 2D–3D-coupling scheme are presented in Refs. [16,17]. The starting point is the global 2D model, which computes the temperature field in the whole furnace including the shape of the solid–liquid interface by adjusting the heating power in such way, that the temperature at the prescribed position of triple point matches the melting point [21]. The solidification of liquid silicon in 2D is modeled by using an interface tracking algorithm, which considers the latent heat sources according to the given growth velocity which is assumed to be equal to the pulling speed Vp. In the melt region the influence of the unsteady 3D melt motion is taken into account using the Reynold’s averaged form of the enthalpy equation: ! rUCU! u UrT g ¼ rUðlg rT g ÞrU a ð1Þ
r—density of liquid silicon; C—heat capacity of liquid silicon; ! u —temporally and azimuthally averaged velocity field; Tg—temperature field (axisymmetric); lg—effective thermal ! conductivity; a —turbulent heat fluxes. ! The temporally and azimuthally averaged velocity field u , ! effective thermal conductivity lg and turbulent heat fluxes a are obtained from the local 3D model as it will be explained further below. 2.1.2. Data transfer from 2D to 3D From the global 2D temperature field the relevant boundary conditions along the periphery of the local 3D model are extracted, i.e. all data transferred from 2D to 3D are axisymmetric:
(i) geometry of the melt domain including the shape of solid– liquid interface;
73
(ii) at the crucible–melt interface the temperature distribution Tc (r, z); (iii) at the solid–liquid interface the melting isotherm T0; (iv) at the free melt surface the heat flux Q2d (r, z). Q2d (r, z) considers the heat flux which is absorbed at the free melt surface from the surrounding environment through radiation and conduction. Q2d enters the temperature boundary condition at the free melt surface of the local 3D model in the following way: @T leff ! ¼ esT 4 Q 2d ð2Þ @n ! where n —surface normal vector (pointed outwards from melt), e—emissivity of the melt, s—Stefan–Boltzmann constant, Q2d— absorbed heat flux.
2.1.3. 3D model Using the geometrical information about the melt region including the present shape of the solid–liquid interface obtained from the global 2D model, the geometry and mesh of the local 3D computational domain is generated. The transient 3D calculation of the convective heat transport is performed with fixed geometry and mesh. The boundary conditions for the velocity field of the local 3D model are a fixed angular velocity at the solid–liquid interface, and at the crucible wall to take into account crystal and crucible rotation, as well as Marangoni shear stress at the free melt surface. In order to consider the turbulent behavior of the melt the LES approach is applied using the Smagorinsky sub-grid model with van Driest wall damping [22]. The resulting governing equations for mass, momentum and energy of the local 3D model are as follows: !
rU U ¼ 0
ð3Þ
! ! @ U ! ! ! þ U Ur U ¼ rP þ rU 2meff S þ rbðT 0 T Þ g rU @t
ð4Þ
@T ! þ U Ur T ¼ rUðleff rT Þ ð5Þ @t ! where U —velocity, P—pressure, meff—effective viscosity of the melt, S—strain rate tensor, which is defined in vector form as ! ! follows S ¼ ðr U þ ðr U ÞT Þ=2, b—thermal expansion coefficient, ! T0—reference temperature (melting point of the melt), g —gravity, C—heat capacity of the melt, leff—effective thermal conductivity of the melt. In the formulation of the governing equations (Eqs. (3)–(5)) an effective melt viscosity meff and an effective thermal conductivity of the melt leff account for the influence of the sub-grid scale turbulence according to the following relations: pffiffiffiffiffiffiffiffiffiffiffiffi mt ¼ rUðC s dÞ2 U 2USUS ð6Þ
rUC
meff ¼ m þ mt leff ¼ CU
m Pr
ð7Þ
þ
mt Prt
ð8Þ
with m—laminar viscosity, mt—turbulent viscosity, calculated according to Smagorinsky–Lilly model with constant Cs ¼0.167 and d—cubic-root of the mesh cell volume. The laminar Prandtl number is Pr ¼0.011479 and the influence of the sub-grid scale turbulence on heat transfer was modeled by using the turbulent Prandtl number Prt ¼0.9 as e.g. in Ref. [11].
74
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
When these relations (Eqs. (11)–(13)) are introduced into the enthalpy equation (Eq. (1)) and the averaging operation (Eq. (9)) is applied one gets from the local 3D model the following equation for the axisymmetric temperature field at stationary state in the melt: !
!0
rUCU/ U Sr/TS ¼ rUð/leff Sr/TSÞ þ rUðrUCU/ U T 0 S/l0 rT 0 SÞ ð14Þ The local 3D (Eq. (14)) and the global 2D model (Eq. (1)) must yield the same temperature fields, i.e. oT4 ¼ Tg. Therefore, the following definitions of the variables of the 2D model using the corresponding values of the 3D model can be made which results in a closed mathematical formulation of the 2D/3D problem: ! ! u ¼ /U S
ð15Þ
lg ¼ /leff S
ð16Þ
!0 ! 0 a ¼ ðrUCU/ U T 0 S/l rT 0 SÞ
ð17Þ !0 0 The temporally and azimuthally averaged values / U T S and 0 /l rT 0 S in (Eq. (17)) are evaluated in local 3D model according to the following relations: !0 ! ! / U rT 0 S ¼ / U rTS/ U SU/rTS 0
/l rT 0 S ¼ /leff rTS/leff SU/rTS
Fig. 1. Flow chart for the coupled global 2D and local 3D model.
2.1.4. Data transfer from 3D to 2D In order to transfer the results obtained from the local 3D model to the global 2D model, the values and cross correlation terms of temperature, velocity and thermal conductivity fields computed in the local 3D model are averaged in time and in the azimuthal direction. The averaging operation for scalar (V) and ! vector values ( V ) is defined according to the following equations: Z 2p Z t0 þ t 1 /VS ¼ VUdtUdj ð9Þ 2pUt 0 t0 "Z # 2p Z t 0 þ t ! 1 ! ! V U e r UdtUdj U e r 2pUt 0 t0 ""Z # #! 2p Z t 0 þ t ! ! ! þ V U e z UdtUdj U e z
ð18Þ ð19Þ
2.1.5. Model for the oxygen transport in the 2D–3D model In addition to the convective heat transport the oxygen transport is calculated by the coupled 2D–3D approach. For that purpose the equation for oxygen transport in the melt, which is analogous to Eq. (5), is solved in the local 3D model. Thereby, the oxygen concentration in the melt is coupled to the transport of SiO in the gas phase of the axisymmetric global model via a flux balance and balancing of the SiO vapor pressure against the azimuthally averaged oxygen concentration at the free melt surface, using the same boundary conditions and material properties as given in Ref. [23]. In the current implementation this coupling of the oxygen concentration in the melt to the SiO concentration in the gas is based on an explicit iterative exchange of the boundary conditions at the free melt surface.
! /V S ¼
0
2.2. Implementation and iteration procedure ð10Þ
t0
where t—time interval for averaging, t0—start instant for the ! ! averaging in time, e r , e z —radial and axial base vectors of the cylindrical reference frame, j—azimuthal angle. From this averaging procedure the distribution of the con! ! vective velocity u , the turbulent heat fluxes a and the effective thermal conductivity lg are determined which are needed in Eq. (1) of the global 2D model. ! Temperature T, velocity U and thermal conductivity leff of the local 3D model can be formally expressed as a sum of a mean ! temporally and azimuthally averaged value oV 4 or o V 4 and ! its corresponding value of the instantaneous fluctuations V’ or V ’ respectively: T ¼ /TS þ T 0
ð11Þ
! ! !0 U ¼ / U Sþ U
ð12Þ
leff ¼ /leff Sþ l0
ð13Þ
Fig. 1 shows a flow chart of the coupled 2D–3D calculations. The coupling is initialized by generating the 2D mesh of the whole CZ furnace in CrysMAS [18]. After the global 2D simulation, the corresponding axisymmetric shape of the melt is used to generate the 3D mesh using Gambit [24]. For the 1st iteration the influence of the unsteady 3D melt motion is not considered directly, instead an initial guess for the effective thermal conductivity lg in the melt is made. During the further iterations the converged axisymmetric solution of the global 2D model is used to update the 3D geometry of the melt region, based on the actual shape of the solid–liquid interface, and to extract the boundary conditions for the local 3D model. For that purpose corresponding interfaces are developed and used to exchange the data between CrysMAS and OpenFOAM. Furthermore, routines are implemented in the local 3D OpenFOAM model to actualize the 3D mesh according to the current shape of the solid–liquid interface. As the changes of the interface shape during the iterations can be considered as relatively small a stretching approach for the mesh points is used instead of the regeneration of the whole 3D mesh.
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
75
Fig. 2. Time- and azimuthally averaged temperature T (a), stream function (b), radial (c) and axial component and (d) of turbulent heat fluxes. Left: this work, LES model with 108,000 cells and aggressive wall refinement; right: DNS model with 8 million cells (from Ref. [7]).
with surface element ds:
Table 1 Used material properties of silicon. Material Property Kinematic viscosity n Thermal expansion b Melting point Tm Thermal conductivity (solid) ls Thermal conductivity (liquid) ll Emissivity (solid) es Emissivity (liquid) el Latent heat L Density r Heat capacity C
Value 3 10 7 m2/s 1.4 10 4 K 1 1683 K 96017/T1.149 W/(m K) 60.26 W/(m K) 0.7 0.3 1.8 106 J/kg 2520 kg/m3 915 J/(kg K)
Next, the transient 3D flow calculation is carried out in OpenFOAM until the system on average reaches a stable state. The solver used is based on the ‘‘buoyantBoussinesqPimpleFoam’’ solver provided with OpenFoam [19]. All calculations are using central difference interpolation schemes, in a reference system rotating with the crucible. Then, after the stable state is reached, the timely and azimuthally averaged fields (Eqs. (15)–(17)) are transferred from the local 3D OpenFOAM to the global 2D CrysMAS model. After interpolation of the 3D results onto the mesh of the 2D model the ! averaged velocities u are corrected by addition of the gradient of !0 a scalar field p to yield a corrected velocity field u in order to ensure conservation of mass and energy: !0 ! u ¼ u þ rp Z
!0 u Uds ¼ 0-
Z
rpUds ¼
Z
! u Uds
ð20Þ
This yields an equation system for p which is solved using the finite volume method implemented in CrysMAS, i.e. forming an equation system by summing up the boundary fluxes. The correction step also includes forcing normal components of averaged velocity and turbulent heat fluxes at axis of symmetry and outer walls to zero. No efforts have been made to enforce zero radial gradients of velocity at the axis of symmetry. On one hand, this results difficult with a simple collocated, equal order scheme on an unstructured triangular mesh for velocities and the pressure-like scalar field. On the other hand, the results shown later indicate that after long enough averaging times this results naturally to a sufficient degree. !0 With the corrected velocity field u the energy equation (Eq. (1)) is solved in the global 2D model. In order to ensure numerical stability of the coupled system, the convective term of Eq. (1) is ramped up during the first few iterations. The update of the boundary conditions for the 3D flow model from the 2D model is done using an under-relaxation factor of typical 0.25. This coupling procedure described above is repeated until the 2D/3D thermal fields and the shape of the phase interface are converged.
3. Test cases 3.1. Comparison to benchmark data from Direct Numerical Simulation As time averaged data from time dependent calculations are required for several 2D–3D iterations the flow calculation procedure in the local 3D model must be as efficient as possible while
76
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
maintaining a reasonable degree of precision. The main problem is posed by the necessity to resolve or model the boundary layers for temperature, momentum, and oxygen. In this work large efforts were made to really resolve the boundary layers down to the laminar region, corresponding to a dimensionless wall distance y þ of the cell nearest to the wall of less than 5 according to von-Karman’s law [25]. This was achieved by a radical refinement of the mesh perpendicular to the wall, and accepting aspect ratios of near-wall-cells up to almost 1000. While the average cell size in the volume of the melt is 5 mm, it goes down to 3 mm at the walls, details are shown later for the example of the realistic Czochralski crucible. In order to investigate the validity of this approach the results obtained from the local 3D LES model with radical mesh refinement close to the walls are compared to benchmark data from a Direct Numerical Simulation (DNS) with approximately 8,000,000 cells for a model Czochralski configuration published by Raufeisen et al. [20]. Fig. 2 shows the time and azimuthally averaged temperature, stream function, and turbulent heat fluxes obtained from the present 3D model for a mesh within total 107,000 cells and the corresponding DNS data [20]. The intention was not to reproduce these results exactly, but to find a suitable compromise between precision and calculation time. The results do not match exactly the DNS data, but all important features are reproduced, and the absolute numbers match to within about 20%. It is believed that this agreement is good enough for an engineering model, considering the extremely low number of cells resulting in calculation speeds high enough to calculate around 600 s of real-time in typically 24 h on the used AMD Opteron 2352 cluster with time steps between 0.02 and 0.05 s. Test runs with smaller time steps did not show any significant difference. 3.2. Test for a realistic Cz crucible without 2D–3D-iteration Next, a simplified case consisting in 2D and 3D only of a 65 kg heavy silicon melt with fixed solid–liquid interface in a crucible with a diameter of 580 mm is considered in order to test the implementation of the coupled 2D–3D model. The thermal boundary conditions for the reduced 2D and 3D models were obtained from a typical, global 2D CZ furnace model containing the crucible with this melt. The crystal and crucible rotation rates were ocrys ¼ 10 rpm and ocruc ¼10 rpm, respectively. The material properties of silicon used in the simulations are summarized in Table 1. The local 3D hydrodynamic LES simulation is carried out on a mesh with 108,000 elements for more than 1000 s real time while the mesh in the 2D simulation has around 33,000 elements in the melt region. Fig. 3 shows the temporal behavior of the average temperature, flow velocity and turbulent heat flux at a monitor point in the center of the melt. It is obvious that after an initial transient time of around 200 s the system reaches a stable state, i.e. there are only minor variations of the averaged field values. Therefore, a total real time of 600 s is considered to be sufficient in the local 3D model because the system is stable for at least 400 s which gives enough data for the averaging procedure to calculate the timely and azimuthally averaged variables. Fig. 4 shows the comparison of the temporally and azimuthally averaged temperature results from the 2D and the 3D model. The average temperature distribution on the right side is obtained from the local transient 3D LES simulation. The corresponding distribution on the left side is calculated in the quasi-stationary 2D model according to Eq. (1) using the averaged values from the 3D model. It is obvious that there is a large discrepancy between the 2D and 3D ! results when turbulent heat fluxes a are not considered in the 2D ! model (Fig. 4a). When the turbulent heat fluxes a are taken into account while solving Eq. (1) in the 2D model, a quantitative agreement between temperature fields of the 2D and 3D models is
Fig. 3. Time averaged temperature o T4 (a), absolute value of velocity oU 4 (b) and turbulent heat flux o a4 and (c) at a monitor point in the melt versus computed real time.
found (Fig. 4b). This result highlights the importance of the turbulent ! heat fluxes a when exchanging data between the 2D and 3D models. Almost no further improvement (Fig. 4c) is achieved by including the terms resulting from the fluctuations of turbulent conductivity (Eq. (17)). The deviation is still less than 2 K. The remaining small discrepancies between the 2D and 3D model are most likely results of limited averaging times as well as of errors from the interpolation from 3D to 2D mesh, and possibly of the consequent flow correction step described in Eq. (20). However, this existing deviation is considered to be negligible compared to all other assumptions and errors which are made when setting up a global model for a Czochralski puller. 3.3. Mesh dependence of the 3D model The dependence of results of the 3D model on the mesh refinement close to the walls of the melt domain has been analyzed for a 24 in. crucible containing 90 kg of silicon. The thermal boundary conditions were extracted from a global 2D furnace model, and certain crystal and crucible rotation rate were used. The mesh refinement at the walls was done using the ‘‘refineWallLayer’’ tool provided with OpenFoam [19], which allows the refinement of selected cells in specified directions. In order to estimate the effect of the mesh refinement an initial mesh with 79,000 cells was generated and different levels of the wall refinement were set. Fig. 5 shows the resulting maximal and averaged dimensionless wall distances yþ at the growth interface as a function of the averaged wall distance of the center of the cells at the boundary for the different levels of the wall
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
77
Fig. 4. Temperature field obtained from global 2D model (left) and averaged temperature field obtained from local 3D model (right); (a) without considering turbulent heat fluxes in 2D; (b) with turbulent heat fluxes in 2D, but without turbulent heat conductivity (i.e. lg ¼ llam); and (c) with turbulent heat fluxes and turbulent heat conductivity (i.e. lg ¼ o leff 4) in 2D.
Fig. 5. Average and maximal dimensionless wall distance y þ at growth interface versus real distance y of the cell center to the wall at the center of the interface for the coarse volume mesh with initial 79,000 cells and different levels of wall refinement (lines) and for the fine mesh with in total 800,000 cells (symbols).
until yþ is everywhere less than 5. In this case the total number of cells is about 800,000. In this case the resulting dimensionless wall distance yþ is slightly smaller than for the coarse mesh with 108,000 cells and with radical wall refinement only. For the coarse volume mesh with the different refinement levels at the walls as well as for the fine volume mesh with the final refinement level the convective heat and oxygen transport in the silicon melt was computed with the local 3D time-dependent LES model assuming that the oxygen concentration is zero at the free melt surface. Fig. 6 shows the timely averaged, integrated heat flux across the growth interface as well as the maximal and averaged oxygen concentration at the interface for the different meshes. It can be seen, that when the dimensionless wall distance yþ becomes less than 5, the change of the heat flux at the interface from the coarse to the fine volume mesh is only about 10%, and the change of the oxygen concentration only about 15%. Therefore, the overall mesh dependent error on the coarse mesh is assumed to be below 25%.
4. Application of the coupled 2D–3D model to an industrial Czochralski process refinement of the initial mesh with 79,000 cells. It is obvious that a cell size of less than 10 mm is needed in order to fulfill the criterion yþ o5. In addition to the wall refinement of the coarse mesh, a fine volume mesh was generated where the initial mesh was refined by dividing each cell into 8 cells, and again refining the wall layers
The combined 2D–3D model has been applied for an industrial Czochralski puller used to grow silicon crystals with diameter of approximately 210 mm for photovoltaic applications from a silica crucible with 24in. diameter. In this case the Grashof number Gr ¼ gUR3 UDTUb =n2 , is¼2.2 1010 with crucible radius R¼12 in.,
78
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
DT¼50 K as a typical maximal temperature difference from crystal to crucible and material properties n and b as given in Table 1. 4.1. Comparison of simulated solid–liquid interface shape with experimental data Experimental data had been provided where the pulling speed Vp was varied during the growth of the crystal body while keeping the other process parameters such as crystal and crucible rotation constant. Lateral Photovoltage Scanning technique [26] has been used to visualize the interface shape on longitudinal cuts. The 2D–3D model was applied to two cases: case A corresponds to a growth stage with a melt weight of around 180 kg at a pulling speed VA p while for case B the melt weight is 130 kg at a pulling speed VBp 4VA p , the actual values being confidential. For both cases, the numerical meshes consist of 40,000 elements in the melt region of the global 2D model and 150,000 elements in the local 3D model. Fig. 7 shows the solid–liquid interface shape for both cases together with the calculated phase boundaries as bold red lines. In view that
there has not been any model adjustment at all, the input consisted merely in a CAD drawing of the puller as well as of process and material parameters, a good correlation between experiment and simulation was obtained. The convergence behavior of the 2D–3D-coupling scheme is shown in Fig. 8 for a case, where the growth rate Vp has been reduced by 50%, starting from a converged solution. Heater power, interface position, maximal crucible temperature and the average temperature difference in the melt between the 3D and 2D model are basically converged after five 2D–3D iterations. 4.2. Comparison of simulated oxygen concentrations to experimental data Apart from the results shown before, the measured interstitial oxygen concentrations had been available for a different growth run with fixed pulling rate Vp in the same puller as in the section above. Calculations including oxygen transport have been performed for three different lengths of the cylindrical crystal body of 2 cm, 102 cm,
Fig. 6. Integrated heat flux Q over growth interface (a) and oxygen concentration [O] at center of growth interface (b) as function of mesh refinement, expressed by the averaged dimensionless wall distance yþ at the interface for the coarse volume mesh with initial 79,000 cells and different levels of wall refinement and for the fine mesh with in total 800,000 cells. The maximum values and the values averaged over the whole area of the interface are shown.
Fig. 7. Comparison of the calculated interface shapes (bold lines) with experimental data (background) for an industrial CZ system at two growth stages; top — low B growth rate VA p , 180 kg melt weight, bottom — high growth rate Vp, 130 kg melt weight. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
79
Fig. 8. Convergence behavior of the 2D–3D-coupling scheme: heater power (a), axial position of solid–liquid interface at the symmetry axis (b), maximum crucible temperature (c), and average temperature deviation in the melt between 2D and 3D model (d).
and 168 cm. These calculations have been done applying two kinds of boundary conditions for the oxygen concentration at the free melt surface. In case O1 the oxygen concentration at the free melt surface was set to zero, which means that there is no coupling of the oxygen transport between the local 3D model of the melt and the global axisymmetric model of the gas atmosphere. In case O2 the oxygen transport in the melt was coupled to the transport in the argon stream, as described in Section 2. The meshes for the three growth stages consist again of typically 100,000 to 150,000 cells with radical refinement at the walls until yþ becomes less than 5. The results of the simulations for both cases O1 and O2 are shown in Fig. 9. The oxygen concentration [O] measured at the center of wafers cut at the three different crystal lengths decreases from the top to the tail end of the crystal. This tendency can only be reproduced for case O2 when the oxygen transport in the melt is fully coupled to the SiO transport in the argon stream. Without coupling the oxygen concentration increases at the tail end again. Taking into account that there has not been any model adjustment at all, the measured and computed values agree again reasonably well.
5. Discussion The examples which were shown in Section 4 highlight that the presented combination of global 2D and local 3D model is able to deliver accurate results for industrial silicon Czochralski crystal growth. The solution of the global 2D problem is not time critical. Therefore, the overall computation time is strongly determined by the speed of the local, transient 3D model and how many iterations are needed to get a converged solution. Taking the example from above one iteration which comprises the computation of in total 600 s real time with 150,000 elements in the local 3D model and one quasi-stationary solution of the global 2D model takes around 1 day, where the cost of the included 2D calculation is negligible. For a converged solution five iterations are typically needed, i.e. the
Fig. 9. Comparison of measured oxygen concentration at the crystals axis versus grown length of the cylindrical crystal body with calculated values for case O1 (without coupling to the gas atmosphere) and case O2 (with coupling to the gas atmosphere).
overall computation time for one set of process parameters is five days. It is not possible to compare this computation time directly with other methods used for combined 2D–3D models because there were no data published. A significant reduction of the time for solving the problem can be expected either when for special cases the local 3D simulations can be significantly accelerated by using for example the approach published in Ref. [27] for the calculation of the influence of horizontal magnetic fields on the melt flow. But, it is also believed that the performance of the implementation of local 3D model on the basis of the popular open source framework OpenFOAM [19] used in present work can be improved further by employing improved parallelization on computer clusters and platforms with high graphical processing power or by utilizing more efficient calculation algorithms.
80
T. Jung et al. / Journal of Crystal Growth 368 (2013) 72–80
6. Summary and conclusions A previously presented idea for modeling of Czochralski processes by using a coupled global 2D–local 3D model was applied to describe the transport phenomena occurring in large industrial scale crystal growth systems. Using a coarse mesh in the melt volume, but a radically refined mesh at the walls, the turbulent 3D motion in the silicon melt was modeled with the LES approach by using the axisymmetric phase interface and boundary conditions from global 2D model. The influence of the unsteady 3D melt flow was accounted by using averaged 3D fields in the temperature equation for the melt region of the global 2D model. The combined model was implemented in CrysMAS and OpenFOAM codes. The model was used to simulate the industrial Czochralski process for growing silicon crystals with 210 mm from a crucible with 2400 diameter. A good quantitative correlation between the measured and simulated shape of the solid–liquid interface for two different pulling speeds was found. The computed and measured oxygen concentration in the crystal for a different growth run agreed also fairly well, when the oxygen transport in the melt is coupled to the SiO transport in the argon stream.
[2] [3] [4] [5]
[6] [7] [8] [9] [10] [11] [12] [13]
[14] [15] [16] [17]
Acknowledgment The authors would like to acknowledge our former co-worker Dr. Andis Rudevics who has made significant contributions to the implementation of the coupling procedure as well as Bosch Solar Energy AG for providing the experimental data. This work, which was funded by the German Ministry of Education and Research (BMBF) was carried out within the ‘‘CZSil’’ project (FKZ03SF0379C) which is part of the German cluster of excellence ‘‘Solarvalley Mitteldeutschland’’.
[18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
References [1] W. Su, R. Zuo, K. Mazaev, V. Kalaev, Journal of Crystal Growth 312 (2010) 495–501.
B. Gao, K. Kakimoto, Journal of Crystal Growth 312 (2010) 2972–2976. L. Liu, K. Kakimoto, Journal of Heat and Mass Transfer (2005) 4481–4491. L. Liu, K. Kakimoto, Journal of Crystal Growth 275 (2005) 1521–1526. V.V. Kalaev, D.P. Lukanin, V.A. Zabelin, N. Makarov, J. Virbulis, E. Dornberger, W. von Ammon, Materials Science in Semiconductor Processing 5 (2003) 369–373. O.V. Smirnova, V.V. Kalaev, N. Makarov, Ch. Frank-Rotsch, M. Neubert, P. Rudolph, Journal of Crystal Growth 266 (2004) 67–73. O.V. Smirnova, V.V. Kalaev, Yu.N. Makarov, N.V. Abrosimov, H. Riemann, Journal of Crystal Growth 266 (2004) 74–80. O.V. Smirnova, V.V. Kalaev, Yu.N. Makarov, N.V. Abrosimov, H. Riemann, Journal of Crystal Growth 287 (2006) 281–286. O.V. Smirnova, V.V. Kalaev, Yu.N. Makarov, N.V. Abrosimov, H. Riemann, V.N. Kurlov, Journal of Crystal Growth 303 (2007) 141–145. ¨ A. Seidl, G. Muller, E. Dornberger, E. Tomzig, B. Rexer, W.v. Ammon, Proceedings of the Electrochemical Society 98-1 (1998) 417–428. D.P. Lukanin, V.V. Kalaev, Yu.N. Makarov, T. Wetzel, J. Virbulis, W. von Ammon, Journal of Crystal Growth 266 (2004) 20–27. A. Krauze, N. Jekabsons, A. Muiznieks, A. Sabanskis, U. Lacis, Journal of Crystal Growth 312 (2010) 3225–3234. I.Yu. Evstratov, V.V. Kalaev, A.I. Zhmakin, Yu.N. Makarov, A.G. Abramov, N.G. Ivanov, E.M. Smirnov, E. Dornberger, J. Virbulis, E. Tomzig, W. von Ammon, Journal of Crystal Growth 230 (2001) 22–29. J. Larsson, S. Kawai, Annual Research Briefs, Center for Turbulence Research, 2010 pp. 39–46. J.J. Derby, L. Lun, A. Yeckel, Journal of Crystal Growth 303 (2007) 114–123. ¨ J. Fainberg, D. Vizman, J. Friedrich, G. Muller, Journal of Crystal Growth 303 (2007) 124–134. T. Tsukada, K.I. Sugioka, C.J. Jing, M. Kobayashi, Journal of Crystal Growth 312 (2010) 997–1004. ¨ M. Kurz, A. Pusztai, G. Muller, Journal of Crystal Growth 198–199 (1999) 101–106. OpenFOAM (R), /www.openfoam.comS. A. Raufeisen, M. Breuer, T. Botsch, A. Delgado, International Journal of Heat and Mass Transfer 51 (2008) 6219–6234. M. Kurz, Ph.D. Thesis, University Erlangen-Nueremberg, 1998. E.R. van Driest, Journal of the Aeronautical Sciences 23 (1956) 1007–1011. A.D. Smirnov, V.V. Kalaev, Journal of Crystal Growth 310 (2008) 2970–2976. Gambit (R), /www.ansys.comS. T. von Karman, Nachrichten von der Gesellschaft der Wissenschaften zu ¨ Gottingen Fachgruppe 1 (Mathematik) 5 (1930) 58–76. ¨ ¨ N.V. Abrosimov, A. Ludge, H. Riemann, W. Schroder, Journal of Crystal Growth 237 (2002) 356–360. Y. Collet, O. Magotte, N. Van den Bogaert, R. Rolinsky, F. Loix, M. Jacot, V. Regnier, J. Marchal, F. Dupret, Journal of Crystal Growth 360 (2012) 18–24.