Analytical method for the second-order homogenization of two-phase composites within Mindlin–Toupin strain gradient elasticity theory is proposed. Direct approach and self-consistent approximation are used to reduce the homogenization problem to the problem of determination of averaged Cauchy stresses, double stresses and static moments of Cauchy stresses inside the inclusions under prescribed quadratic boundary conditions. The ellipsoidal shape of inclusions and orthotropic properties of phases are assumed. Extended equivalent inclusion method with linear eigenstrain is proposed to derive the explicit relations between the Eshelby-like tensors and corresponding concentration tensors that are used to define the averaged field variables inside the inclusions. Obtained analytical solutions allow to evaluate the effective classical and gradient elastic moduli of composite materials accounting for the phases properties, volume fraction, shape and size of inclusions. Presented solution for the effective gradient moduli covers the full range of volume fraction and correctly predicts the absence of gradient effects for the homogeneous classical Cauchy medium. Examples of calculations for the composites with spherical inclusions are given. Micro-scale definition of the well-known simplified strain gradient elasticity theory is provided. Namely, it is shown that this theory is the phenomenological continuum model for the composites with isotropic matrix and with the small volume fraction of stiff isotropic spherical inclusions, which have the same Poisson’s ratio to those one of matrix phase.