In traditional kinetic theory for a granular flow, it is usually assumed that particle interactions are instantaneous and binary. For a dense granular system, these assumptions are usually invalid. In this paper, we use an ensemble averaging technique to examine the effects of finite particle interaction time and multiparticle collisions. The main objectives of this paper are to develop a method and to provide a tool to study dense granular materials. As an example, we study flows of granular particles coated with thin layers of resin. To model particle elasticity and resin viscosity, the force between a pair of particles is approximated by a serial connection of a linear spring and a dashpot. Subsequently, a viscoelastic model is developed from the averaging method. In order to determine coefficients in the constitutive model, direct numerical simulations are performed. When the particle concentration is relatively low, the shear stress is quadratically proportional to the shear rate, in agreement with kinetic theories. At a high particle concentration, the shear stress depends linearly on the rate of strain. The transition between this quadratic and linear dependence is similar to a phase transition. In a dense system, when the shear rate exceeds a critical value, shear band formation is also observed.