In a crystal growth system, the crystal quality is greatly affected by the coupling properties between unsteady melt flow and thermal transfer. In this paper, an improved lattice Bolzmann method is proposed. This incompressible axisymmetric model based method transforms the fluid equations of cylindrical coordinate into those of the two-dimensional Cartesian coordinate and constructs the evolutionary relationship of the external force terms, such as rotational inertia force and the thermal buoyancy. In the unsteady melt, the temperature distribution and the rotational angular velocity are determined based on the D2Q4 model and the velocity of axisymmetric swirling fluid is calculated based on the D2Q9 model. The mirror bounce format is adopted as the boundary conditions of the free surface and the axis symmetry. For the remaining boundary conditions, the non-equilibrium extrapolation format is used. In the simulation, 12 sets of flow function results are obtained by choosing different sets of Grashof number and Reynolds number. By comparing with the finite crystal growth results, the effectiveness of the proposed method can be shown. Furthermore, by studying the convection shape and the temperature distribution of the melt under coupling between high Grashof number and high Reynolds number, it can be concluded that the thermal coupling properties and flow in the unsteady melt relate to Grashof number and Reynolds number. By adjusting the high Reynolds number generated by the crystal and crucible rotation, the strength of the forced convection in the melt can be changed. Therefore, the natural convection in the melt can be suppressed effectively and the temperature distribution results can be improved significantly. In addition, it is worth mentioning that the findings in this paper can be straightforwardly extended to the silicon single crystal growth experiment by turning the dimensionless crystal rotation Reynolds number and crucible rotation Reynolds number into the actual rotation speed.