Heterogeneous materials are commonly found in nature (e.g. soil, rock, blood) andengineering applications (e.g. paints, lubricants, sintered materials). Understanding the microstructure-property relations is crucial to predicting how these materials fail and behave during processing. In this thesis, we propose a generalization of classical homogenization estimates foranisotropic linear two-phase composites. These generalized estimates are consistent with known bounds for the effective properties and satisfy certain duality properties relating the primal and dual formulations of the homogenization problem. Particularly for composites with large phase contrast, the new estimates are expected to have improved predictive capa- bilities as they can capture – at least qualitatively – the onset of clustering of the inclusion phase at non-dilute volume fractions. The new estimates also have important implications for the phase averages of the stress and strain-rate fields as well as the field fluctuations. In addition, for the special case of monodisperse suspensions of rigid spheres in a Newtonian fluid, the new estimates for the shear viscosity are validated against experimental results. The generalized estimates are then used to generate corresponding estimates for theeffective behavior of nonlinear viscoplastic composites. To that end, we employ the variational linear comparison methods, which make use of optimally chosen linear comparison composite media, allowing the direct conversion of estimates for linear composites into corresponding estimates for nonlinear composites. For suspensions of rigid spheres in a Herschel-Bulkley fluid, variational estimates for the effective yield and flow stress are ob- tained and validated against experimental results. In addition, second-order linear com- parison homogenization estimates for nonlinear viscoplastic composites with infinite phase contrast are generated and the dependence of the effective behavior on the macroscopic loading conditions is investigated. Special attention is devoted to the method’s ability to describe the highly localized anisotropic deformation patterns that emerge in composites with a finite yield stress. Under pure shear loading conditions the dependence of the effec- tive yield stress is found to be non-analytic for dilute inclusion volume fractions which is shown to be connected to the emergence of shear bands.