Purpose: The purpose of this paper is to preset a spherical element subdivision method for the numerical evaluation of nearly singular integrals in three-dimensional (3D) boundary element method (BEM). Design/methodology/approach: In this method, the source point is first projected to the tangent plane of the element. Then two cases are considered: the projection point is either inside or outside the element. In both cases, the element is subdivided into a number of patches using a sequence of spheres with decreasing radius. Findings: With the proposed method, the patches obtained are automatically refined as they approach the projection point and each patch of the integration element is “good” in shape and size for standard Gaussian quadrature. Therefore, all kinds of nearly singular boundary integrals on elements of any shape and size with arbitrary source point location related to the element can be evaluated accurately and efficiently. Originality/value: Numerical examples for planar and slender elements with various relative location of the source point are presented. The results demonstrate that our method has much better accuracy, efficiency and stability than conventional methods.