In this article, we present a $T$-matrix method for numerical computation of second-harmonic generation from clusters of arbitrarily distributed spherical particles made of centrosymmetric optical materials. The electromagnetic fields at the fundamental and second-harmonic (SH) frequencies are expanded in series of vector spherical wave functions, and the single sphere $T$-matrix entries are computed by imposing field boundary conditions at the surface of the particles. Different from previous approaches, we compute the SH fields by taking into account both local surface and nonlocal bulk polarization sources, which allows one to accurately describe the generation of SH in arbitrary clusters of spherical particles. Our numerical method can be used to efficiently analyze clusters of spherical particles made of various optical materials, including metallic, dielectric, semiconductor, and polaritonic materials.