In a dense granular system, particles interact in networks containing many particles and interaction times are long compared with the particle binary collision time. In these systems, the streaming part of the granular stress is negligible. We only consider the collisional stress in this paper. The average behavior of particle contacts is studied. By following the statistical method developed recently by the authors [Zhang and Rauenzahn, J. Rheol. 41, 1275 (1997)], we derive an evolution equation for the collisional stress. This equation provides guidance to collateral numerical simulations, which show that the probability distribution of particle contact times is exponential for long contact times. This can be explained by network interactions in a dense granular system. In general, the relaxation of the collisional stress is a combined effect of the decay of the contact time probability and the relaxation of collisional forces among particles. In the numerical simulations, the normal force between a pair of particles is modeled as parallel connect of a spring and a dashpot. In this case, the relaxation of the force magnitude conditionally averaged given a specific contact time is negligible, and the major contribution to the stress relaxation is from the exponential decay of the contact time probability. We also note that the probability decay rate is proportional to the imposed strain rate. Consequently, in a simple shear flow with a constant particle volume fraction, as the shear rate approaches zero, the shear stress approaches a finite value. This value is the yield stress for that particle volume fraction. Hence, the evolution equation of the collisional stress predicts viscoplasticity of dense granular systems.