This paper presents a combined numerical-theoretical study of the macroscopic behavior and local field distributions in a special class of two-dimensional periodic composites with viscoplastic phases. The emphasis is on strongly nonlinear materials containing pores or rigid inclusions. Full-field numerical simulations are carried out using a fast Fourier transform algorithm [Moulinec, H., Suquet, P., 1994. A fast numerical method for computing the linear and nonlinear properties of composites. C. R. Acad. Sci. Paris II 318, 1417–1423.], while the theoretical results are obtained by means of the ‘second-order’ nonlinear homogenization method [Ponte Castañeda, P., 2002. Second-order homogenization estimates for nonlinear composites incorporating field fluctuations. I. Theory. J. Mech. Phys. Solids 50, 737–757]. The effect of nonlinearity and inclusion concentration is investigated in the context of power-law (with strain-rate sensitivity m) behavior for the matrix phase under in-plane shear loadings. Overall, the ‘second-order’ estimates are found to be in good agreement with the numerical simulations, with the best agreement for the rigidly reinforced materials. For the porous systems, as the nonlinearity increases ( m decreases), the strain field is found to localize along shear bands passing through the voids (the strain fluctuations becoming unbounded) and the effective stress exhibits a singular behavior in the dilute limit. More specifically, for small porosities and fixed nonlinearity m > 0 , the effective stress decreases linearly with increasing porosity. However, for ideally plastic behavior ( m = 0 ) , the dependence on porosity becomes non-analytic. On the other hand, for rigidly-reinforced composites, the strain field adopts a tile pattern with bounded strain fluctuations, and no singular behavior is observed (to leading order) in the dilute limit.