Checklist for model preparation

For the melt flow computation the following procedure can be suggested.

  1. First the Reynolds numbers of the crystal and crucible rotations and the Grashof number based on the crucible radius should be computed. For computation of the characteristic numbers see section Computing Characteristic Numbers. The square root from the Grashof number and rotational Re numbers can be used for estimation of the melt flow nature. For laminar stationary flow a converged numerical solution can be obtained with natural material parameters of the melt. Special care should be taken on the mesh refinement in the boundary layers under the crystal and on the crucible wall.

    If characteristic numbers exceed the laminary and stationary mode, then is clear that the 2D modeling supplies only an approximation of the realistic 3D melt flow. In 2D two possibilities are being used. First, the melt viscosity can be increased. A factor 10 is normally suffcient for the medium size Cz melts in order to get the converged stationary solution.

    The second possibility is an application of the k-epsilon turbulence model to the melt flow. The validity of the assumtions behind of the turbulence model is conditional in case of the complex Czochralski melt flow with rotation and/or by small Prandtl numbers of the semiconductor melts. Very fine numerical mesh and a strong underrelaxation is required for the 2D computation of the turbulent melt flow. Generally the solutions of the turbulence model are not much more reliable as results of the laminar computations with increased viscosity, but the stability and computational time of the turbulence modeling is much more critical.

    Once the melt flow mode and applied flow model are clarified, one should prepare the numerical mesh, material data and boundary conditions.

  2. The mesh density in the melt domain should be sufficient for resolution of the thin boundary layers, which may be especially difficult to resolve at the crystal melt interface. For estimation of the boundary layer thickness see the first table row Boundary layer thickness(m) in section Structured mesh Properties.

    The careful mesh preparation is required also on the unstructured mesh side. First, the temperature in the nodes of the structured mesh at the interface between the structured and unstructured mesh is interpolated linearly between neighbouring nodes of the unstructured mesh. Therefore the unstructured mesh density defines the field resolution along the boundary of the hybrid mesh. It is also the case for interface separating two structured regions with different materials. The thermal degrees of freedom are again the nodes of the unstructured mesh at the region boundary hiddened by the structured mesh. If this boundary is radiating, the radiative heat fluxes are assembled also in the nodes of the unstructured mesh.

    Second, the edges of the structured mesh match exactly with the edges of the unstructured mesh at the phase boundary crystal-melt. Therefore first the suffficient mesh density should be created for the unstructured mesh there and then the generation of the structured mesh can proceed.

    A number of edges along the phase boundary should be sufficient for the computational stage while the phase boundary became strongly curved either towards the crystal or towards the melt. The requirement of the high spatial resolution under the phase boundary leads to narrow control volumes along the phase boundary. If the number of control volumes (edges) along the phase boundary is poor, then one obtains control volumes with high aspect ratio. The line connecting central nodes of the neighboring control volumes will leave the area of them by some curvature radius of the curved phase boundary because of the angle between neighboring edges. This creates large interpolation errors and can lead to the numerical instability. One counteracts to this problem increasing the number of edges along the phase boundary still in the stage of the unstructured mesh generation.

  3. The material parameters which are important for execution of the melt flow are

    • dynamic viscosity
    • thermal conductivity
    • density
    • specific heat capacity
    • thermal expansion coefficient
    • Marangoni coefficient

    If the Lorentz force of the electromagnetic field is considered, then additionally electric conductivity of the melt and of all other auxilliary materials should be entered.

  4. If the computation should begin the first time and no temperature distribution was still obtained, then control the value of Start Residual in the Other dialog box of the Computation -> Numerical Parameters -> Forward tab or deactivate the checkbox Track interface in the Other dialog box of the Computation -> Numerical Parameters -> Forward tab. The phase tracking should begin after an initial temperature distribution will be obtained. If the start of the phase tracking is intended immediately, control the underrelaxation of the phase tracking in the Start Residual dialog window. Its value should not exceed 0.01, otherwise jumps in the iterative deformation of the numerical mesh may occur and the mesh will be corrupted.

    Set reduced value of the underrelaxation for enthalpy equation residual in the same dialog window in the Forward relaxation factor dialog field and select the GSSV solver. These measures increase the stability of the iteration at the beginning of the computation. Later the underrelaxation can be set back to unity.

  5. Control the start value of the heater powers. If the regulated heaters are used, then control additionally the thereshold value of the enthalpy equation for start of the inverse modelling, the underrelaxation of the inverse model and used control points. The parameters of the inverse modeling are accessible in the dialog tab Computation -> Numerical Parameters -> Inverse. The parameter Controller relaxation factor should be set equal to small value of about 0.1. The underrelaxation damps strong variation of the heating power during the inverse computation. The Start controller at parameter in the same dialog tab prescribes the residual of the enthalpy equation (or the same, of the forward simulation) which should be achieved before the inverse modelling starts.

    If the system is heated inductivelly, then it is suggested to use only Fixed power type of heater and prescribe the electric current. Later the value of electric currents can be adjusted manually on the temperature in the regulation points. The manual adjustments of electric currents is preferred, because the inverse modelling in combination with inductive heating is working not reliably.

    For the power controlled computation check parameters of the control points. One useful control point should be placed on the triple point location crystal-melt-gas in the Cz configuration or on the crystal-melt-liquid encapsulant location in the LEC or VCz configuration. Important for behaviour of the inverse modeling is the parameter toler: in the dialog window Control Points. Smaller value in Kelvin causes higher acccuracy for absolute temperature but requires more iterations in order to solve the inverse problem. The range between 1 and 3 Kelvin is advised for the control point coinciding with the triple point of Cz configuration. The tolerances of another control points in the system should be selected with account of the positioning and measurement precision of the thermoelements in the Cz puller.

    The heater powers regulation can be activated immediately or later, after an initial solution for temperature is obtained. For switching off of the inverse modelling just change a type of the controlled heaters to the Fixed power in the Heaters dialog window.

  6. Control settings for the specific latent heat of the solidification and the crystal pulling velocity in the dialog window Computation -> Proceess Parameters -> General. Corresponding values are entered into the Growth rate [m/s]: and Latent heat [J/m³] dialog fields. The volumetric latent heat means the same that the phase transition heat parameter of the two phase material in the enthalpy method multiplied with the crystal density.
  7. Control settings of the boundary conditions for the radial, axial and azimuthal velocities at the melt interface and if available at the interface of the liquid encapsulant. For settings of boundary conditions see the firsdt part of this tutorial.
  8. Select the melt (and liquid encapsulant) region and mark them for convection computation using the dialog window Physical Phenomena. If the turbulent convection should be computed in the melt, in the first step do not activate the turbulence in the melt region, but increase the melt viscosity and try to obtain the convrged solution for the laminary flow. This solution can serve as a start point for the modeling with turbulence.