bonded particle interaction

Dear Developers,

I have a trouble in understanding the rotational bonded particle interaction.

If two same-radius particles are bonded in the X-Y plane, and one of them is set a certain angular velocity in Z direction, then I expect that the two particles will have the same angular acceleration. However, the simulation gives different angular accelerations.

If the two particles interact via RotFric model, the simulation produces the same angular accelerations.

The simple code and geometry files is as follows.

Script 1
two_bonded_particles.py
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
from math import *
import numpy

#create a simulation container object:
sim = LsmMpi (numWorkerProcesses=1, mpiDimList=[1,1,1])
maxRadius = 1.001
minRadius = 1.001
sim.initNeighbourSearch (
 particleType = "RotSphere",
 gridSpacing = 2.500*maxRadius,
 verletDist = 0.200*minRadius
)

#specify the number of timesteps and timestep increment:

sim.setTimeStepSize(1.0000e-3)

#enforce 2d computation
sim.force2dComputations(True)

#specify the spatial domain for simulation:
domain = BoundingBox ( Vec3 (-200.0,-200.0,0.0), Vec3 (200.0,200.0,0.0) )
sim.setSpatialDomain (domain)

sim.readGeometry("2.geo")

density = 1.0
#Set particle density
sim.setParticleDensity (
 tag = 1,
 mask = -1,
 Density = density #depends on the porosity
)

sim.setParticleDensity (
 tag = 2,
 mask = -1,
 Density = density #depends on the porosity
)

elasticContactModulous = 10.000
contactStrength = 3000000000.000

#create rotational elastic-brittle bonds between particles:
cement = sim.createInteractionGroup (
 RotBondPrms(
  name="cement",
  normalK=elasticContactModulous,
  shearK=1.0*elasticContactModulous,
  torsionK=0.0,
  bendingK=2.0*elasticContactModulous,
  normalBrkForce = contactStrength,
  shearBrkForce = 2.0*contactStrength,
  torsionBrkForce = 0.0,
  bendingBrkForce = 20000.0*contactStrength/3.0,
  tag = 1, ##if tag = 0, there is non bond; if tag = 1, there is bond.
  scaling = False
 )
)

#create frictional interaction between particles:
sim.createInteractionGroup (
 RotFrictionPrms(
  name = "friction",
  normalK = elasticContactModulous,
  dynamicMu = 100.0,
  staticMu = 100.0,
  shearK = 1.0*elasticContactModulous,
  scaling = False,
  rigid=True
 )
)

##create an exclusion between bonded and frictional interactions:
sim.createExclusion(
 interactionName1 = "cement",
 interactionName2 = "friction"
)

##set the damping as zero in order to check whether the moments on connecting particles are same.
sim.createInteractionGroup (
  LocalDampingPrms(
     name="damping1",
     viscosity=0.0
  )
)

sim.createInteractionGroup (
  RotLocalDampingPrms(
     name="damping2",
     viscosity=0.0
  )
)

##fix particle non translational
#sim.setParticleNonTranslational(1)
#sim.setParticleNonTranslational(2)

sim.runTimeStep()
particleList=sim.getParticleList()
for pp in particleList:
 print pp.getId(),pp.getAngularAcceleration()[2]

##set particle angular velocity
sim.setParticleAngularVelocity(0,Vec3(0.0,0.0,0.01))

sim.runTimeStep()
particleList=sim.getParticleList()
for pp in particleList:
 print pp.getId(),pp.getAngularAcceleration()[2]

sim.runTimeStep()
particleList=sim.getParticleList()
for pp in particleList:
 print pp.getId(),pp.getAngularAcceleration()[2]

sim.exit()

##Script 2
2.geo
LSMGeometry 1.2
BoundingBox -10.0 -10.0 0.0 10.0 10.0 0.0
PeriodicBoundaries 0 0 0
Dimension 2D
BeginParticles
Simple
2
0.0 1.0 0.0 1.001 0 2
0.0 -1.0 0.0 1.001 1 1
EndParticles
BeginConnect
1
0 1 1
EndConnect

Best regards,
Jiadun