This study presents model extensions for a lattice Boltzmann (LB) approach to thermal axisymmetric flow including swirl or rotation. An incompressible axisymmetric lattice Boltzmann D2Q9 model was applied to solve the axial and radial velocities through inserting source terms into the two-dimensional lattice Boltzmann equation. The equations governing azimuthal (or swirling) velocity and the temperature were also solved by the LBM. It is found that this scheme is much more stable and consistent compared to previous hybrid schemes. It provides a significant advantage in simulation of melt flows with high Reynolds number and high Grashof number. The present scheme was validated by comparing the LB results with benchmark solutions for melt flow in Czochralski crystal growth. Unsteady flows with high Grashof numbers were studied in detail. The critical Grashof number for the onset of the oscillation is found to be about 2.5 x 10(6). The oscillation amplitude psi(max) is proportional to (Gr - Gr(c))(0.5) for 2.5 x 10(6) < Gr < 6 x 10(6). The frequencies and flow patterns of the unsteady flows are also analyzed. The distributions of the mean quantities of the temperature and rms of temperature at Grashof number as high as 6 x 10(7) is found to be similar to those obtained by 3D simulations. (C) 2014 Elsevier Ltd. All rights reserved.