With the rapid development of nuclear industry in China, an increasing number of nuclear power plants (NPPs) are unavoidable to be constructed upon non-lithological foundations. Therefore, it is essential to optimize the current evaluation methods for non-lithological foundation reinforcement plans. By using the Octree discretization method combining with the scaled boundary polyhedron finite element method (SBFEM), a refined calculation model of a Generation III+ NPP is constructed. Due to the highly nonlinear dynamic characteristics of the system involved, this paper adopts a generalized plasticity model for the soil, a plastic damage model for the piles, and a nonlinear overburden layer-based seismic input method. Subsequently, large scale plastic damage analysis in detail to the NPPs under beyond-design basis earthquake inputs is performed. The results show that the pile foundation damage characteristics and distribution acquired from simulations agree well with the data from actual earthquake disasters. The most severe damage to the pile foundation occurred at the interface between the end of the pile and the base mat. The inclusion of plastic damage showed that the floor response spectra of the nuclear island experiences significant changes in three dimensions. Therefore, the effect of pile damage to the dynamic response of the power plants should be considered for strong earthquake excitation. The computational and analysis methods proposed in this paper can effectively reflect the dynamic nonlinear response of the nuclear power plants pile foundation-soil interaction.