MindlinDeresiewitz Contact Law Issue

Hello!

I am trying to compare behavior of MindlinDeresiewitz contact model with HertzMindlin's.

I was able to successfully run HertzMindlin contact law with no issues. However, when applying "Law2_ScGeom_MindlinPhys_MindlinDeresiewitz()" law, no interaction is detected between the sphere and the wall. In other words, upon inspecting the scene, there is no penetrationDept -> kn=0, and ks=0.

Below is the script to run and test the two contact models... I kindly ask you to comment and uncomment the Law2 functors to see the difference between the two contact models. It will automatically show the pre-sliding behavior when HertzMindlin constitutive law is used while no results can be seen when you activate the MindlinDeresiewitz law:

##########################################################################
# Script for validating the MindlinDeresiewitz contact model
# Simulating a sliding sphere with known coefficient of friction
# The script is used to evaluate the instantaneous compliances and dependence on both the current and previous state loading sequences

from yade import plot
import numpy as np

##################################################################################################################
### Section to define the thrust loading function
cyc = 10000 # number of cycles with constant force application (sufficiently large to reach steady state)
step_inc = 20 # amount of force increment [N]
Ftmax = 280 # maximum tangential force to reach during loading phase[N]
Ftmin = 200 # minimum tangential force to reach during unloading phase [N]

N = (Ftmax-Ftmin)/step_inc + 1
Ft = np.zeros(int(abs(Ftmax*cyc/step_inc + cyc))) # pre-allocating memory to store tangential force values

for i in range(0,len(Ft)): # loop to define the loading phase
    Ft[i] = step_inc * np.floor(i/cyc)

Ft = np.append(Ft, np.flip(Ft[int(-N*cyc):-cyc])) # appending the UNLOADING phase to the previous loading phase

Ft = np.append(Ft, np.flip(Ft[int(-N*cyc):-cyc])) # appending the RELOADING phase to the previous unloading phase

# Note: the maximum thrust force is below the maximum static friction that causes the body to accelerate
##################################################################################################################

nu1 = 0.3 # Poisson's ratio of the first material (infinite wall)
nu2 = 0.3 # Poisson's ratio of the second material (sphere)
E1 = 2.13e11 # Young's Modulus (Pa) of the first material (infinite wall)
E2 = 2.13e11 # Young's Modulus (Pa) of the second material (sphere)
mu = 0.3 # coeffcient of static/dynamic friction [-]
Rad1 = 0.001 # defining radius of the sphere [m]
Fn = 1e3 # defining the constant normal force [N]
targetiter = len(Ft) # target iteration before quiting the simulation

Mat1 = O.materials.append(FrictMat(young=E1, poisson=nu1, frictionAngle = np.arctan(mu), label='Mat1')) # assembling material properties
Mat2 = O.materials.append(FrictMat(young=E2, poisson=nu2, frictionAngle = np.arctan(mu), label='Mat2')) # assembling material properties

# Bodies to add in the simulation (wall and sphere)
O.bodies.append(
        [
                # add an infinite wall to the assembly (with z-axis being its normal)
                utils.wall( Vector3(0,0,0), axis = 2, sense=0, material = Mat1),
                # add a sphere tagent to the constructed wall
                utils.sphere((0, 0, Rad1), radius=Rad1, material = Mat2)
        ]
)

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb()], label="collider"),
    InteractionLoop(
            [Ig2_Wall_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_MindlinPhys()],
            [Law2_ScGeom_MindlinPhys_Mindlin()],
            #[Law2_ScGeom_MindlinPhys_MindlinDeresiewitz()],
    ),
        NewtonIntegrator(damping=0.1), # gravity is deactivated
        PyRunner(iterPeriod=1, command='loading_fun()', label='external_load'), # function to apply the external loads on the sphere
        PyRunner(iterPeriod=cyc, command='addPlotData()', label='datastream'), # plots results in real-time
        PyRunner(iterPeriod=cyc, command='PARSE()', label='postprocessing'), # storing data of concern (for post-processing)
        PyRunner(iterPeriod= targetiter-1, command='SAVING()', label='tab_delimited') # saves outputs of concern to a neutral file format
]

O.bodies[1].state.blockedDOFs = 'XYZ' # suppress rotations of the sphere
O.engines += [ForceEngine()]
O.dt = utils.PWaveTimeStep()
FrictionForce, SphereDisp, data = [], [], [] # preserving memory for variables to be parsed
plot.plots = {'displacement': ('Ff',)}
plot.plot()
O.stopAtIter = targetiter
O.run()

# user-defined function to apply external loads on the free sphere
def loading_fun():
        O.engines[-1].force = Vector3(0,Ft[O.iter],-Fn) # apply normal forces (negative z-axis)
        O.engines[-1].ids = [1] # associate external loads with the free sphere

# user-defined function to track values of interest
def addPlotData():
        plot.addData(
               displacement = O.bodies[1].state.pos[1], # store the tangential displacement of sphere
               Ff = -1*O.interactions[0,1].phys.shearForce[1] # store the developed friction force
                )
        print (O.engines[-1].force[1], -1*O.interactions[0,1].phys.shearForce[1] , O.interactions[0,1].phys.normalForce[2] , O.bodies[1].state.vel[1]) # inspecting loads on sphere and its velocity to check how close is the analysis to quasi-static condiiton

# user-defined function to parse data and save in tab-delimited txt file (postprocessing)
def PARSE():
        global FrictionForce # updating variable to global kind to allow storage
        global SphereDisp # updating variable to global kind to allow storage
        global data # updating variable to global kind to allow storage
        FrictionForce = np.append(FrictionForce , -1*O.interactions[0,1].phys.shearForce[1]) # appending the shear force [N] to the data stream
        SphereDisp = np.append(SphereDisp, O.bodies[1].state.pos[1]) # exracting sliding displacemennt [m]
        data = np.transpose(np.array([FrictionForce, SphereDisp]))

def SAVING():
        comment = '''The following results are for a sphere with %s mm radius being pressed against an elastic half-space under constant normal load...
Sphere Material Properties: Young's Modulus = %s GPa, Poisson's ratio = %s
Half-space Material Properties: Young's Modulus = %s GPa, Poisson's ratio = %s
Coefficient of static friction = %s
Frictional Force [N]\tDisplacement[m]'''%(str(Rad1*1000), str(E2/1e9), str(nu2), str(E1/1e9), str(nu1), str(mu))
        np.savetxt(r'Presliding_Behavior_MindlinDeresiewitz.txt', data, delimiter='\t', header=comment)

###########################################################################