Stressing Over Stress Tensor Symmetry in the Navier-Stokes Equations How do we know that the stress tensor must be symmetric in the Navier-Stokes equation? Here are some papers that discuss this issue beyond the usual derivations:
Behavior of a Vorticity Influenced Asymmetric Stress Tensor In Fluid Flow http://www.dtic.mil/dtic/tr/fulltext/u2/a181244.pdf 
Vorticity and Stress Tensor https://www.rgnpublications.com/journals/index.php/jims/article/view/1061 
Diffusion Stratification in a Rotating Linearly Stratified Fluid https://link.springer.com/article/10.1007/BF02698703
First, consider a shear flow v = (Cy,0,0) where C is a constant with the dimensions of 1/time. Diffusion of momentum at the molecular dynamic level implies that there is a shear stress in the x-direction along planes that have their normals in the y-direction. If the stress-tensor is symmetric, then there must be exactly the same magnitude of shear-stress in the y-direction along planes that have their normals in the x-direction, even though there is no momentum in the y-direction to diffuse to account for such a stress in the y-direction. If such an additional stress exists, by what molecular dynamic process is it brought into existence and how do we model it mathematically? 
Note that, at the level of the Navier-Stokes equations, the term that is used to symmetrize the stress tensor only makes a non-vanishing contribution if the flow either has density gradients and/or viscosity gradients – e.g. the flow is (roughly) either compressible and/or has temperature gradients. Thus a simple flow that might illustrate the presence of a non-zero antisymmetric (aka rotational) viscosity might be: an incompressible viscous flow of uniform density in a cylinder of radius R with full slip boundary conditions with a large and constantly maintained temperature gradient in the radial direction that is initially rotating uniformly (i.e. same angular speed everywhere) as fast as possible without inducing turbulence. 
Physically, one would expect the higher temperature of the outer bands (of constant radius) would initially cause a net diffusion of angular momentum radially inward, create a viscous torque density, and induce angular acceleration of the inner bands – even without an initial rate of deformation of the volume elements. One would expect angular momentum to be conserved globally since there are no external torques acting on the system – not because there are not equal and opposite torques between interacting fluid elements within the system - (i.e. the stress tensor does not necessarily need to be symmetric for angular momentum to be conserved globally).
Mathematically, in cylindrical polar coordinates, the theta component of the Navier-Stokes equation for the covariant velocity produces a viscous torque density term that is in fact proportional to the gradient of the rotational viscosity times the vorticity. Thus, if the rotational viscosity does not vanish, the initial uniform rotation of the fluid would induce a viscous torque on the inner bands that would cause them to undergo angular acceleration – consistent with what one might expect physically.
I’m not sure how to think about the final steady state that this system evolves to other than to suspect that: it conserves angular momentum globally, it dissipates as much angular kinetic energy as possible under this constraint, and all accelerations and force densities vanish. If one plugs in an ansatz such that the final distribution of angular speed is proportional to the radius raised to some power, one can show that the angular speed vanishes everywhere except at the origin (where it is infinite), that angular momentum is still conserved globally, that all the angular kinetic energy dissipates into heat, and that all force densities vanish.
There is some talk in the literature of angular speeds possibly blowing up at the origin for some Navier-Stokes flows. Needless to say, taking the limit as r goes to zero doesn’t make much sense for a physical fluid as it would be less than the inter-particle spacing (which implies a rotational inertia of the fluid that does not vanish and an angular speed that is consequently not infinite). There is also some talk of Coriolis forces and increases of angular speeds being quite noticeable for rotating flows with density (and presumably temperature) gradients.
So how do we know with certainty that the rotational viscosity is identically zero? And if it is identically zero, what molecular dynamic mechanism accounts for the presence of orthogonal shear forces in a pure shear flow?
 A: There are a lot of questions here, so I'll only answer two of them:
1. Why is the stress tensor symmetric in the Navier-Stokes equations?
It's because we've defined it that way for this model/these equations. Remember the definition of the stress tensor in the Navier-Stokes model:
$$\bar{\bar{\sigma}} = \lambda (\nabla\cdot\vec{u})\bar{\bar{I}} + \mu (\nabla\vec{u} + \nabla\vec{u}^T)$$
Both individual components of this expression are obviously symmetric, and therefore the stress tensor in this model must be. You may say that you do not believe this model for the stress tensor is the correct model, in which case you are saying the Navier-Stokes equations do not apply.
2. When is the stress tensor symmetric in general?
The symmetry of the stress tensor follows directly from conservation of angular momentum if all angular momentum in the fluid is as a result of moments derived from applied forces: see here for a derivation. If the fluid which you are dealing with is undergoing volumetric torques, or if it is capable of transmitting surface couples, then the stress tensor is not symmetric and those liquids are referred to (confusingly) as polar fluids. These liquids are usually considered exotic.
A: It looks like the method of control volumes can be used to generate additional physical equations from the Navier-Stokes equation (in conservation form). If you take the cross product of both sides with the radial vector and pull the cross product through the divergence operators, additional terms must be subtracted out as per the chain rule. Taking integrals of both sides over some control volume and using the divergence theorem yields that the time rate of change of the angular momentum within the control volume is equal to (negative of) the rate of transported loss of angular momentum through its surface plus the torque generated on the surface - which must from an equation that holds on its own by physical grounds (and is also equivalent to a continuity diffusion equation for angular momentum with a torque source)  - PLUS the volume integral of the term shed off by the chain rule (now integration by parts) which must therefore be zero. This vanishing latter term is just the volume integral of the stress tensor contracted with the Levi-Civita symbol, and so the stress tensor is arguably symmetric.
If one takes the dot product of both sides of the Navier-Stokes equations with the velocity vector and follows the analogous steps, one gets the enthalpy- continuity kinetic energy-diffusion equation with a pressure source plus a second continuity-diffusion equation for internal kinetic energy with mechanical and dissipative sources.
As far as deriving the value of the shear stress in the y-direction perpendicular to the shear flow above (which has no momentum in the y-direction) directly from molecular dynamic considerations (as opposed to deriving it as a consequence of the conservation of angular momentum), it is still not clear to me how to visualize the causal mechanism. If one puts a set of blocks in a row on (the x-axis of) a table and exerts a shear force on them in the x-direction, they might rotate except that a shear force would be generated in the y-direction on each of their front faces as they collide with the back face of the block in front of them... Is this the right way to visualize the molecular dynamics that cause the stress tensor to be symmetric? 
Regarding the fluid beginning in rigid body rotation with a radial temperature gradient, it looks like the viscous torque term from the gradient of the anti-symmetric/rotational viscosity times the vorticity (twice the initial uniform angular rotation rate) does indeed lead to an increase in angular speed everywhere (except at the boundary) generating angular momentum "out of nowhere" during the first time step dt. Physically, it still seems to me that the (temperature enhanced) rate of inward diffusion of tangential velocity should induce angular acceleration of the inner bands during the first time step away from rigid body rotation - but it is not clear to me how to model this. 
Anyway - thank you for your help!
