In this study, three-dimensional Treffz method is developed for micromechanical modeling of heterogeneous materials with nanoscale inclusions, considering the interface stress effect. Spherical harmonics are used as the Papkovich-Neuber potentials to derive the Trefftz trial displacement fields in the matrix material as well as in the nano-inclusions. Characteristic lengths are used to scale the Trefftz trial functions, to avoid ill-conditioning of the derived system of linear equations. The compatibility between the independently assumed fields are enforced by collocation at a finite number of points, along the interface between the matrix and inclusions. The stress jump across the matrix/inclusion interface is described by the generalized Young–Laplace equation for solids, which is also enforced by the collocation method. Numerical results by the proposed Trefftz method are consistent with available analytical solutions in the literature, demonstrating the high accuracy of the present method. The effective elastic constants of solids containing nano-inclusions show significant size-dependency, in contrast to those for composites without considering interface effects. Interactions of multiple nano-inclusions are also studied, which can be used as benchmark solutions in future studies.