This paper describes a new approach to the high-fidelity simulation of axisymmetric free-surface flows. A boundary-fitted grid is coupled with a new projection method for the solution of the Navier–Stokes equations with second-order accuracy in space and time. Two variants of this new method are developed by adapting two existing algorithms, suitable for prescribed velocity boundary conditions, to the case of prescribed normal and tangential stresses at the free boundary. A normal-mode analysis for a fixed-boundary problem confirms the second-order accuracy of the algorithms. The approach is validated by comparison with a Rayleigh–Plesset solution for an oscillating spherical bubble, with an analysis of shape oscillations, and with existing results for the buoyant rise of a deforming bubble for Reynolds numbers up to 200 and Weber numbers up to 12. In addition to the simulation of axisymmetric free-surface flows of intrinsic interest, the present approach is suitable for the validation of genuinely three-dimensional calculations.