We simulated two particle-based fluid models, namely multiparticle collision dynamics and dissipative particle dynamics, under shear using reverse nonequilibrium simulations (RNES). In cubic periodic simulation boxes, the expected shear flow profile for a Newtonian fluid developed, consistent with the fluid viscosities. However, unexpected secondary flows along the shear gradient formed when the simulation box was elongated in the flow direction. The standard shear flow profile was obtained when the simulation box was longer in the shear-gradient dimension than the flow dimension, while the secondary flows were always present when the flow dimension was at least 25% larger than the shear-gradient dimension. The secondary flows satisfy the boundary conditions imposed by the RNES and give a total flow field with a lower rate of viscous dissipation than the corresponding unidirectional flows. This work highlights a previously unappreciated limitation of RNES for generating shear flow in simulation boxes that are elongated in the flow dimension, an important consideration when applying RNES to complex fluids like polymer solutions.