Owing to its excellent corrosion, radiation and high temperature resistance, ZrO ₂has been considered as a strong candidate material for inert fuel for the incineration of actinides. In this paper, a combination of thermal spike model and molecular dynamics is used to simulate the phase transition process of ZrO ₂in the nuclear radiation environment. Based on the thermal spike model, two coupled diffusion equations are established with considering the multiple physical process of energy deposition and transmission after the implantation of swift heavy ions into target material. The space-time evolution characteristics of ZrO ₂lattice temperature are obtained by solving the coupled diffusion equations numerically. Then the phase transformation of ZrO ₂form monoclinic phase to tetragonal phase under the thermal spike is investigated on an atomic scale by means of molecular dynamics. It is found that a cylindrical track with a radius of 7 nm is generated in the center of ZrO ₂after the implantation of swift heavy ion with an electronic energy loss of 30 keV·nm ^–1. The lattice melts immediately in the center of track, accompanied with the coordination number of Zr decreasing from 7 to 4–6. Then at about 2 ps, the melting zone gradually turns cool and recrystallized. And in the center of the melting zone, voids begin to form and are surrounded by a highly disordered amorphous region. Meanwhile, tetragonal phase of ZrO ₂, whose coordination number of Zr is 8, is formed at the periphery of the amorphous region, which is also confirmed by the XRD calculation results. As energy transfers from track center to the surround, the tetragonal region gradually develops into the whole system, accompanied with the increase of voids size. The simulation results indicate that the irradiation of ZrO ₂with swift heavy ions can lead to a transformation from the monoclinic to the tetragonal phase when the deposited electronic energy loss exceeds an effective threshold ~21 keV·nm ^–1, greater than the experimental value (12 keV·nm ^–1), which was mainly due to the large difference between the simulated and measured incident ion fluences and the accuracy of the force field used in the molecular dynamics.