Shear Flow of non-spherical particles

Dear lammps users,

I’m new to lammps and trying to simulate shear flow of non-spherical particles (dumbblle shape) using fix rigid and granular package.

The particles are create from two overlapping spheres (see attached figure) and initially configured in a 25x25x25 lattice (using with a vol fraction of 0.6

I am using the following script for my simulation

units si
atom_style sphere
boundary p p p
newton off
comm_modify mode single vel yes

#Read the Cluster of compound particles and add mol property
fix addMolecule all property/atom mol ghost yes
read_data dumbells.dat fix addMolecule NULL Molecules

#Relax the particles
pair_style gran/hooke/history 2000 NULL 50 NULL 1 0
pair_coeff * *

velocity all create 1e15 1234 loop geom

neighbor 0.003 bin
neigh_modify delay 0 every 1 check yes

timestep 1e-5
fix 1 all rigid molecule

run 30000

#Remove Previous Settings
unfix 1
uncompute 1
uncompute 2

reset_timestep 0

#New Seetings
change_box all triclinic remap

pair_style gran/hooke/history 2000 NULL 50 NULL 0.5 1
pair_coeff * *

timestep 1e-5

fix 1 all rigid molecule

fix shearRun all deform 1 xy erate ${shearRate} remap v units box

run ${NSTEPS}

I expect to get a linear velocity profile but I get an S shape profile. This is independant of the time step size. I have attched velocity profile results for three different shear rates. At very low rates ({shearRate}~1) the profile is almost linear but this quickly changes to an S shape by slightly increasing the shear rate ({shearRate}>2). The simulation also seems to have reached steady state since I tested it with two different ${NSTEPS} with larger nstep twice the smaller one.

I appriciate any suggestions from the developpers or the community.

Thank you