This paper develops a three-dimensional (3D) model for a heterogeneous half-space with inclusions distributed periodically beneath its surface subject to elastohydrodynamic lubrication (EHL) line-contact applied by a cylindrical loading body. The model takes into account the interactions between the loading body, the fluid lubricant and the heterogeneous half-space. In the absence of subsurface inclusions, the surface contact pressure distribution, the half-space surface deformation and the lubricant film thickness profile are obtained through solving a unified Reynolds equation system. The inclusions are homogenized according to Eshelby’s equivalent inclusion method (EIM) with unknown eigenstrains to be determined. The disturbed half-space surface deformations induced by the subsurface inclusions or eigenstrains are iteratively introduced into the lubricant film thickness until the surface deformation finally converges. Both time-independent smooth surface contact and time-dependent rough surface contact are considered for the lubricated contact problem.