A novel numerical method called the Thermal Discrete Dipole Approximation (T-DDA) isproposed for modeling near-field radiative heat transfer in three-dimensional arbitrary geometries. The T-DDA is conceptually similar to the Discrete Dipole Approximation, except that the incident field originates from thermal oscillations of dipoles. The T-DDA is described in details in the paper, and the method is tested against exact results of radiative conductance between two spheres separated by a sub-wavelength vacuum gap. For all cases considered, the results calculated from the T-DDA are in good agreement with those from the analytical solution. When considering frequency-independent dielectric functions, it is observed that the number of sub-volumes required for convergence increases as the sphere permittivity increases.Additionally, simulations performed for two silica spheres of 0.5 µm-diameter show that the resonant modes are predicted accurately via the T-DDA. For separation gaps of 0.5 µm and 0.2 µm, the relative differences between the T-DDA and the exact results are 0.35% and 6.4%, respectively, when 552 sub-volumes are used to discretize a sphere. This work suggests that the T-DDA is a robust numerical approach that can be employed for solving a wide variety of near-field thermal radiation problems in three-dimensional geometries.