We study a class of shear-free, homogeneous but anisotropic cosmological models with imperfect matter sources in the context of f(R) gravity. We show that the anisotropic stresses are related to the electric part of the Weyl tensor in such a way that they balance each other. We also show that within the class of orthogonal f(R) models, small perturbations of shear are damped, and that the electric part of the Weyl tensor and the anisotropic stress tensor decay with the expansion as well as the heat flux of the curvature fluid. Specializing in locally rotationally symmetric spacetimes in orthonormal frames, we examine the late-time behaviour of the de Sitter universe in f(R) gravity. For the Starobinsky model of f(R), we study the evolutionary behavior of the Universe by numerically integrating the Friedmann equation, where the initial conditions for the expansion, acceleration and jerk parameters are taken from observational data.