Heavy damage in materials in the limit of high-dose irradiation is important both for fundamental materials science and nuclear industrial applications. Inspired by a recent study (Phys. Rev. Mater. 4, 023605 (2020)), we computationally investigate heavy radiation damage in alpha zirconium (Zr) on the atomic scale with molecular statics calculation of Frenkel pair accumulation (FPA), which contains repeated random insertions of Frenkel pairs into the sample and energy minimizations between the insertions. This method allows us to constantly introduce point defects into a Zr crystal to reach a high dose of several displacement per atom (dpa), and meanwhile, to avoid thermally-activated migration and recovery of defects. We observe an asymptotic approach to a steady damage level, which features a high point defect content (∼3%), nanoscale dislocation lines and loops, and significant anisotropic deformation strains. We also find that small interstitial clusters play an important role in the initial evolution of the radiation damage. Finally we simulate the damage produced by repeated displacement cascades in Zr at 0.1 K using molecular dynamics (MD) and compare it with that by FPA, obtaining a scaling relationship which can quantitatively describe the difference between the radiation damage caused by the two methods.