The main objective of this thesis aims at developing a new generation of software products, in order to obtain a fully automatic simulator predicting the entire Czochralski process while handling correctly the switches between the different growth stages. First of all, new efficient, robust and high-quality automatic mesh generation algorithms with enough flexibility for any complex geometry were implemented, including a 1D mesh generator by global grade-adaptive method, a 2D initial triangulation algorithm by improved sweep line technique and an automatic 2D shape-quality unstructured mesh generator by modified incremental Delaunay refinement technique. Secondly, a Finite Element Navier-Stokes solver based on unstructured meshes was developed and validated. Enhanced turbulence models based on the classical mixing-length or k-l model, together with a generic transformation method to avoid negative k when solving the turbulent kinetic energy equation by the Newton-Raphson iterative method were introduced and implemented. Moreover, laminar and turbulent mathematical models governing the gas convection, thermal distribution and oxygen concentration were developed, and Finite Element numerical methods to solve these governing equations on unstructured meshes were implemented, while appropriate numerical approaches to capture the wall shear stress exerted by the gas flow and experienced by the silicon melt were investigated. Finally, a series of numerical experiments devoted to investigate the industrial Czochralski crystal growth process under various growth conditions are presented based on all the developments implemented. Comparisons of the simulation results with literature and available experimental observations are also presented, and conclusions are drawn based on these simulation results and observations.