Rotary cement kilns are complex systems that involve occurrence of several simultaneous physical phenomena in both the bed and gas region. In this paper, a comprehensive 3D mathematical model for rotary cement kiln including heavy fuel oil combustion, hot flow, zone heat flux model for clinker formation and rotating wall was developed. The developed CFD model was able to predict key parameters like burning zone position, the maximum gas temperature with and without clinker formation, the hot spot zone, shell temperature and species concentration. The comparison of predicted data and experimental measurements under industrial operation were quite satisfactory. The burning zone temperature is about 1900 K and is located between 5 and 20 m along the rotary kiln. These results indicated that the present model may be used to improve the understanding of the process, the optimization and diagnosis of the operational problems.