In this paper, the Gaussian quadrature formula is used to approximately calculate two-dimensional hypersingular integrals of the form ∬Sf(x,y)r3dS on the circular domain. The two-dimensional hypersingular integrals are changed into iterated integrals by polar transformation. The inner layer of the integral has singularity and Gaussian quadrature formula is constructed for this part. The outer layer without singularity is calculated by Gauss–Legendre quadrature formula of Riemann integral. The corresponding theoretical analysis is given. Furthermore, this type of two-dimensional hypersingular integral is extended to the form ∬Sf(x,y)rp+1dS(p>2) and an error estimate is also derived. Several numerical examples are presented to verify the effectiveness of the proposed algorithm, no matter the singularity is located at the center or near the boundary of the circular domain.