When the individual particles in an otherwise quiescent suspension of freely suspended spherical particles are acted upon by external couples, the resulting suspension-scale fluid motion is characterized by a non-symmetric state of stress. Viewed at the interstitial scale (i.e. microscopic scale), this coupling between translational and rotational particle motions is a manifestation of particle–particle hydrodynamic interactions and vanishes with the volume fraction $\phi$ of suspended spheres. The antisymmetric portion of the stress is quantified by the suspension-scale vortex viscosity $\mu_v$, different from the suspension's shear viscosity $\mu$. Numerical boundary element method (BEM) simulations of such force-free suspensions of spheres uniformly dispersed in incompressible Newtonian liquids of viscosity $\mu_0$ are performed for circumstances in which external couples (of any specified suspension-scale position-dependence) are applied individually to each of the suspended particles in order to cause them to rotate in otherwise quiescent fluids. In the absence of external forces acting on either the spheres or boundaries, such rotations indirectly, through interparticle coupling, cause translational motions of the individual spheres which, owing to the no-slip boundary condition, drag neighbouring fluid along with them. In turn, this combined particle–interstitial fluid movement is manifested as a suspension-scale velocity field, generated exclusively by the action of external couples. Use of this scheme to create suspension-scale particle-phase spin fields $\mbox{\boldmath \Omega}$ and concomitant velocity fields $\bm v$ enables both the vortex and shear viscosities of suspensions to be determined as functions of $\phi$ in disordered systems. This scheme is shown, inter alia, to confirm the constitutive equation, ${\textsfbi T}^{\hspace*{1pt}a} \,{=}\, 2\mu_v\bm \mbox{\boldmath \varepsilon}\,{\bm\cdot}\,[(1/2) {\bm \nabla}\,{\times}\, \bm v - \mbox{\boldmath \Omega}],$ proposed in the continuum mechanics literature for the linear relation between the antisymmetric stress ${\textsfbi T}^{\hspace*{1pt}a}$ and the disparity existing between the particle-phase spin rate $\mbox{\boldmath \Omega}$ and half the suspension's vorticity, ${\bm \nabla} \,{\times}\, \bm v$ (with the third-rank pseudotensor $\bm \mbox{\boldmath \varepsilon}$ the permutation triadic). Our dynamically based BEM simulations confirm the previous computations of the Prosperetti et al. group for the dependence of the vortex viscosity upon the solids volume fraction in concentrated disordered suspensions, obtained by a rather different simulation scheme. Moreover, our dynamically based rheological calculations are confirmed by our semi-independent, energetically based, calculations that equate the rates of working (equivalently, the energy dissipation rates) at the respective interstitial and suspension scales. As an incidental by-product, the same BEM simulation results also verify the suspension-scale Newtonian constitutive equation, ${\textsfbi T}^{\hspace*{1pt}s}\,{=}\, \mu[{\bm \nabla}\bm {v} + ({\bm \nabla}\bm {v})^{\dag}]$, as well as the functional dependence of the shear viscosity $\mu$ upon $\phi$ found in the literature.