Large scale hydrodynamic simulations are presented for the rheology of sheared suspensions of elastic spheres. Suspensions dynamics are studied for constant shear rate and constant stress. For dilute to moderately concentrated conditions up to volume fraction φ = .50, the particles remain disordered and suspension behavior resembles hard sphere suspensions. At higher volume fractions in the range .50 < φ < .60, the suspensions may undergo a disorder-order transition with hexagonal packing patterns. Three dimensions probability distributions g(x, y, z) are presented to characterize the suspension microstructure, while simulation videos are shown to illustrate the dynamics of structure evolution. The phase boundary is a function of the volume fraction, dimensionless elastic modulus and the strength of lubrication forces as characterized by a minimum effective separation gap δ__min for hydrodynamic force calculation. The calculation of the effective lubrication strength requires solution of a time dependent differential equation for the coupled elastic deformation-viscous film drainage quations for each particle contact. Constant shear rate simulations show large fluctuations in instantaneous viscosity at high volume fractions, φ > .50. Elastic spheres present a useful model system for investigating the approach to discontinuous shear thickening and jamming phenomena in concentrated suspensions. While direct simulation of hard sphere suspensions above φ = .50 is difficult, the use of elastic spheres with increasingly large elastic modulus provides a useful model system to study hydroclusters and the rapid shear thickening regime at large φ. These systems also represent an interesting intermediate regime between rigid hard spheres and concentrated suspensions of deformable droplets.