# 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

import numpy as np

##################################################################################################################
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

Ft[i] = step_inc * np.floor(i/cyc)

# 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 [-]
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
]
)

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=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.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
O.engines[-1].force = Vector3(0,Ft[O.iter],-Fn) # apply normal forces (negative z-axis)
O.engines[-1].ids =  # associate external loads with the free sphere

# user-defined function to track values of interest
displacement = O.bodies.state.pos, # store the tangential displacement of sphere
Ff = -1*O.interactions[0,1].phys.shearForce # store the developed friction force
)
print (O.engines[-1].force, -1*O.interactions[0,1].phys.shearForce , O.interactions[0,1].phys.normalForce , O.bodies.state.vel) # 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) # appending the shear force [N] to the data stream
SphereDisp = np.append(SphereDisp, O.bodies.state.pos) # 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))

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

## Question information

Language:
English Edit question
Status:
Expired
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Malek Jaran (malek-utwente) said on 2023-03-14: #1

 Revision history for this message Launchpad Janitor (janitor) said on 2023-03-26: #2