Abstract: Based on the similarity of integration area of random inclusions, similar sub-domain boundary element method scheme has been presented for the first time in this paper. Then, the 2D solid with randomly distributed inclusions can be reduced to a multiply connected domain of the matrix with inner boundary conditions. In this way, computational efficiency is enhanced significantly comparing to the conventional multi-domain approach of the FEM or BEM. As numerical examples, plenty of numerical computation for the solid with randomly distributed circular or elliptic inclusions is performed using the similar sub-domain BEM scheme. Furthermore, the interface between the matrix and inclusion can be not only ideal interface, but also the interface with interphase layers. The above-mentioned computation provides reliable numerical simulation methods for the investigation of the macroscopically effective properties of the corresponding fiber-reinforced composites.