In this chapter, a Computational Grain is developed for direct numerical modelling of composites with nanoscale inclusions considering both interface stretching and bending effects, using a large number of CGs in a representative volume element. The CGs developed in this chapter are by far the first and the only numerical tool for direct numerical modelling of nanocomposites with a large number of inclusions with Steigman-Ogden matrix/inclusion interfaces. By using a new boundary-type multi-field variational principle together with Papkovich-Neuber potentials, the stiffness matrices of CGs can be directly evaluated and assembled. Together with the parallel algorithms, it is found that very efficient simulations of nanocomposites can be realized, for example, a SERVE containing 10,000 nano inclusions only takes fifty minutes on a sixteen-core computer. The influence of spatial distributions of the nano inclusions on the overall properties of nanocomposites is also investigated in this chapter. We also study the influence of interface bending resistance parameters on the effective modulus of nanocomposites. Numerical results show that interface bending resistance parameters affect the shear modulus of nanocomposites but their effect on the bulk modulus is negligible.