A three-dimensional numerical model was established based on the volume of fluid method. The thermal-fluid process in the laser cleaning thermal barrier coatings is analyzed in detail. The effect of scanning speed on the cleaning process was investigated by tracing the migration position of the air-coating phase interface, with consideration of gravity and vapor impact pressure. It was demonstrated that phase change of the zirconia coating is dominated by the vaporization-evaporation process when it is ablated by a millisecond laser; the ablation process of the C263 alloy substrate is a solidification-melting process. The bulge at the edge of the gully will affect the morphology of the gully by influencing the temperature distribution.