An inclusion refers to localized eigenstrains appearing in such processes as thermal expansion and plastic deformation. In view of micromechanics, the existence of inclusions may significantly influence the mechanical properties of the engineering materials. A micromechanical model is proposed to determine the variation of the strain energies in the presence of the near-surface inclusions. The corresponding inclusion problem in a half-space is usually difficult to be solved analytically. In this work, the strain energy is evaluated numerically via the method of images, which superposes the counterpart solutions in full-space and eliminates the tractions on the boundary surface of the half-space. The validity of the present work is confirmed by comparing with the published results and the finite element method (FEM).