In the present study a phase field model for high surface energy anisotropy is developed to model the morphologies of voids in UO2. In order to precisely account for the high anisotropy, an alternative forward-backward strategy based on a staggered grid with an averaged interface normal scheme is used in the numerical procedure. A variety of equilibrium void shapes are reproduced with respect to a constant volume condition. The facet areas and facet energies are calculated. The simulations show excellent agreement with the analytic predictions obtained through Wulff constructions. For the void shapes with high Miller index facets, it is discovered that a slight decrease in total surface area will result in a substantial increase in the total surface energy.