Topology optimization has become a useful tool for designing Phononic Crystals (PnCs) to achieve desired elastic wave propagation properties. However, topological design of three-dimensional (3D) PnCs remains a significant challenge due to extremely high computational cost required for the repeated solution of the nested double-loop problem of design optimization and eigenvalue analysis to determine the band structure during the optimization iteration process. This paper presents a more efficient topology optimization method based on the concept of successive iteration of analysis and design. The method employs a sequential approximation of the eigenpairs through an inverse iteration-like procedure to avoid solving the expensive eigenvalue problems in each design iteration. To further improve the efficiency of the design iteration, the eigenmodes of the intermediate PnCs designs are first computed on a relatively coarse finite element mesh. Then they are mapped onto the fine mesh and used as the initial trial eigenmodes in the eigenpair analysis and design sensitivity analysis. The high efficiency of the proposed method, as compared with conventional nested double-loop optimization approaches, is demonstrated by 3D numerical examples with over one million degrees of freedom. It is shown that this approach can achieve simultaneous convergence of the eigenvectors along with the design evolution of the unit cell, therefore substantially reduces the computational burden for high-resolution topological design of 3D PnCs.