Metallic nanoparticles (NPs) support localized surface plasmon resonances (LSPRs), which enable to concentrate sunlight at the active layer of solar cells. However, full-wave modeling of the plasmonic solar cells faces great challenges in terms of huge computational workload and bad matrix condition. It is tremendously difficult to accurately and efficiently simulate near-field multiple scattering effects from plasmonic NPs embedded into solar cells. In this work, a preconditioned volume integral equation (VIE) is proposed to model plasmonic organic solar cells (OSCs). The diagonal block preconditioner is applied to different material domains of the device structure. As a result, better convergence and higher computing efficiency are achieved. Moreover, the calculation is further accelerated by two-dimensional periodic Green's functions. Using the proposed method, the dependences of optical absorption on the wavelengths and incident angles are investigated. Angular responses of the plasmonic OSCs show the super-Lambertian absorption on the plasmon resonance but near-Lambertian absorption off the plasmon resonance. The volumetric method of moments and explored physical understanding are of great help to investigate the optical responses of OSCs.