The simulation of the crystal growth by the Czochralski method (Cz) is associated with serious difficulties. The most critical issue by the Cz modelling is the numerical method for modelling of the melt flow. The combination of the forces acting on the melt contains the buoyancy, centrifugal, Coriolis, viscose shear stress, Marangoni flow and sometimes also the Lorentz force provided the electromagnetic field is applied. The relations between different forces and dimensionless properties of the melt flow can be expressed by means of the characteristic numbers.

- Prandtl number - ratio between the kinematic viscosity and thermal diffusivity, is a material property
- Reynolds number - ratio between the inertia and viscose shear forces
- Peclet number - ratio between the convective and conductive heat fluxes, is a product of the Reynolds and Prandtl number
- Grashof number - ratio between the buoyancy and viscose shear forces
- Rayleigh number - a product of the Grashof number with the Prandtl number
- Marangoni number - ratio between the surface shear force due to the Marangoni effect to the viscose forces.
- Hartmann number - criterion of the strength of the stationary magnetic field influence on the melt flow
- Magnetic Reynolds number - criterion of the strength of the velocity field influence on the magnetic field
- Magnetic Taylor number - ratio between the Lorentz force of the rotating magnetic field to viscose forces.

The melt flow mode depends on the characteristic numbers and geometric aspect ratio, or the same on the melt size, on the material properties of the melt, on the direction of the heat flux, on the temperature difference and on the applied rotation rates of the crystal and of the melt.

The most simple case for the numerical modeling takes place if the characteristic numbers of the melt flow are inside of the area of the stationary and laminar flow. Such melt flows correspond to very small sizes of the crucible and/or very small temperature gradients and rotation rates of the crystal and the melt.

For the stationary and laminar melt flow a relatively simple system of transport equations can be formulated and the solution can be obtained for the velocity field, heat transport within the melt and shape of the equilibrium phase boundary crystal-melt. This solution is thermally fully coupled with the global 2D model. Also the fluid flow in the media other than melt can be computed simultineously and the interaction effect between the fluid flow in melt and in gas or in melt and in the liquid encapsulant at their shared interface can be accounted for.

Unfortunately most of the practical Cz applications cannot satisfy to the parameters of the laminar and stationary flow regime because larger and larger crystal sizes are produced in the industrial scale. The sizes of the melt volume are increased correspondingly. A good example of teh crystal scale up is the silicon Czochralski where the standard crystal cross-section and weight were rapidly increased until over 300 mm diameter and 300 kG. The melt flow is in any case very complex by the large scale Cz. It is 3D and almost always turbulent.

The time-dependent melt flow in the Cz-configuration is frequently 3D. Therefore a rigorous 3D time-dependent formulation is required for the quantitative description of the most practical cases of the Cz melt flow.

Current approaches are dealing with different kinds of compromises between the computational cost and accuracy of prediction. Application of the direct numerical simulation (DNS) are most expensive. They require weeks of parallel computing and can be applied today in the best case to the medium size Cz melt flow by Grashof numbers about 10^9. The turbulent melt flow at such characteristic numbers occur by the melts diameters about 20 cm and more. Another possibility to compute the turbulent Cz melt flow is to apply the Large Eddy Simulation technique (LES). LES allows to reduce the linear mesh size by at least factor 2 in comparison with DNS. The time step may be increased in LES according to the CFL criterion too. The resulting reduction of the computational cost was about factor 20 in the known LES applications. Also in this case of the LES application for the medium scale Cz flow remains computationally intensive but it can be executed under certain conditions within several days on the parallel cluster or on the multicore workstation.

Another possibility to compute the turbulent Cz melt flow is an application of the quasi-DNS. This method means the straight forward application of the DNS method with coarser numerical mesh than required for the best resolution of the Kolmogorov length scales of the turbulent motion. The smallest flow structures couldn't be resolved in the quasi-DNS method. The cut-off of the fine flow structures is also the principle of the LES method. The effect of the fine scale on the resolved flow is modelled in the LES by spatial filtering. The effect of the coarse mesh on the resolved flow in quasi-DNS is expected to be analogious to the LES due to the artificial viscosity introduced into the discretization by the coarse mesh. The coarse mesh effect of the quasi-DNS should be still studied in detail in comparison with the LES performance.

The typical aim of the Cz modeling is to describe the effect of the melt flow on the timely averaged heat and species transfer and the shape of the phase boundary crystal melt. The lost effect of the fine flow structure on the coarse mesh is mainly its contribution to the turbulent mixing. This contribution can be considered as an effect on the effective transport coefficients for the momentum, heat and species transport. With the progressing mesh refinement the amount of the turbulent energy contained in the cut-off fine scale is reduced and its contribution to the considered averaged thermal result of the melt flow should decrease.

The practical solution which can be applied for 2D modeling of Cz systems with parameters exceeding the stationary and laminary limit are the laminar model with increased value of the melt dynamic viscosity and application of the standard k-epsilon turbulence model. One should be aware that crude simplifications of any accessible turbulence model in the 2D limit may conflict with the realistic complex 3D dynamics of the large scale Cz melt flow.