  We numerically investigate viscoelastic phase separation in polymer solutions
under shear using a time-dependent Ginzburg-Landau model. The gross variables
in our model are the polymer volume fraction and a conformation tensor. The
latter represents chain deformations and relaxes slowly on the rheological time
giving rise to a large viscoelastic stress. The polymer and the solvent obey
two-fluid dynamics in which the viscoelastic stress acts asymmetrically on the
polymer and, as a result, the stress and the diffusion are dynamically coupled.
Below the coexistence curve, interfaces appear with increasing the quench depth
and the solvent regions act as a lubricant. In these cases the composition
heterogeneity causes more enhanced viscoelastic heterogeneity and the
macroscopic stress is decreased at fixed applied shear rate. We find steady
two-phase states composed of the polymer-rich and solvent-rich regions, where
the characteristic domain size is inversely proportional to the average shear
stress for various shear rates. The deviatoric stress components exhibit large
temporal fluctuations. The normal stress difference can take negative values
transiently at weak shear.
