About KinemCNDEngine and KinemCNLEngine for simple shear test

Hello,

I'm currently working on a simulation involving simple shear under constant volume and constant normal force conditions. I'm using the kinemCNDEngine for this purpose, but I'm facing some challenges with the implementation. Despite referencing the documentation and forum discussions, I haven't been able to achieve the desired behavior in my simulation. I'm hoping someone could provide some insights or guidance to help me resolve these issues.

In the code, I'm attempting to use the kinemCNDEngine to simulate simple shear. However, after compaction, I'm not observing any wall movement to induce simple shear. Interestingly, when I check the values of f0 and y0, I find that they are both set to 0.0. I suspect that I might be missing something crucial in the implementation, but I can't pinpoint the issue.

I would greatly appreciate any assistance or suggestions that could help me get this simulation on track.

Thanks in advance for the help.

from yade import pack

compFricDegree=33
targetPorosity = 0.40
compFricDegree = 33
finalFricDegree = 33
rate = -0.01
damp = 0.2
stabilityThreshold = 0.001
young = 25e6
mn, mx = Vector3(0, 0, 0), Vector3(15, 15, 15)

O.materials.append(FrictMat(young=young, poisson=0.4, frictionAngle=radians(compFricDegree), density=2650, label='spheres'))
O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=0, density=0, label='walls'))

walls = aabbWalls([mn, mx], thickness=0, material='walls')
wallIds = O.bodies.append(walls)

sp = pack.SpherePack()

psdSizes,psdCumm=[0.1,0.6,0.752,0.8,0.837,0.874,0.9,0.915,0.945,0.967,1.0,1.046,1.112,1.2,1.6],[0,0.01,0.025,0.0550,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.972,1.0] # particle size distribution

sp.makeCloud(mn, mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,porosity=0.46,seed=1)

O.bodies.append([sphere(center, rad, material='spheres') for center, rad in sp])

simple_shear= KinemCNDEngine(
 dead = True,

)

triax = TriaxialStressController(
        maxMultiplier=1. + 2e4 / young,
        finalMaxMultiplier=1. + 2e3 / young,
        thickness=0,

        stressMask=7,
        internalCompaction=True,
)

def checkUnbalanced():
 if O.iter < 5000:
  return
 if unbalancedForce() > .01:
  return
 triax.dead=True
 simple_shear.dead= False
 print("Shearing")
 checker.dead = True

newton = NewtonIntegrator(damping=damp)

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=100, timestepSafetyCoefficient=0.8),
        triax,
        newton,
        simple_shear,
        PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]

triax.goal1 = triax.goal2 = triax.goal3 = -50000

from yade import plot

def history():
            plot.addData(

         e11=-triax.strain[0],
         e22=-triax.strain[1],
         e33=-triax.strain[2],
         ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
                s11=-triax.stress(triax.wall_right_id)[0],
                s22=-triax.stress(triax.wall_top_id)[1],
                s33=-triax.stress(triax.wall_front_id)[2],
                i=O.iter,
                #u_0=-flow.getPorePressure((5,0.0015,5)),
                #u_25=-flow.getPorePressure((5, 3, 5)),
                #u_50=-flow.getPorePressure((5, 5, 5)),
                #u_avr=flow.averagePressure(),
                #Qin = flow.getBoundaryFlux(2),
  #Qout = flow.getBoundaryFlux(3),

                )

O.run(100,True)