# Shear loading

Asked by mrhappy on 2020-01-23

I made a simple model to understand what’s going on. I have two spheres:
sphere(center=(0,0,0),radius=0.5,fixed=True,material='mat1'),
sphere((0,0,1),radius=0.5,material='mat1'),
The following is the boundary condition:
O.bodies[1].state.blockedDOFs = "xyzXYZ"
O.bodies[1].state.vel = (0.0001,0,0)
O.bodies[1].state.angVel = (0,0,0)

I am applying a velocity in the x direction on sphere #1 which should cause a shear reaction force to occur. The following are my questions:
1. When I track the forces on sphere #1, I notice that there is a force in the x-direction (which is expected) and a force in the z-direction (uniaxial). Where is the uniaxial force coming from if I only applied a shear boundary condition? The contrary example to this would be to apply a uniaxial loading condition which would result in only a uniaxial force.
2. I have read [1] and [2] thoroughly. I understand that the tangential plane rotates as the spheres shear and that the shear force is calculated incrementally. Is deltaUs in [2] the change in tangential displacement with respect to time or is it from the initial configuration?
3. I tried to do a simple DEM model on MATLAB where I had two 2D spheres and I applied the same loading condition as my YADE model while calculating the uniaxial and shear forces on sphere#1. I could not get the same results as YADE predicts. Is there a document that explains clearly the steps. I read [1], but I want something with more detail.
4. My code is running well, so I assume that the my code is correct.

[1] https://yade-dem.org/doc/formulation.html#fig-shear-2d
[2] A DEM model for soft and hard rocks: Role of grain interlocking on strength

Thank you soo much for your time.

MY CODE:
# MATERIAL PROPERTIES
intR=20 # allows near neighbour interaction (can be adjusted for every packing)
DENS=2500 # Density
YOUNG=1800
FRICT=7
ALPHA=0.1
TENS=100
COH=100
iterMax = 100

O.materials.append(JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH,label='mat1'))

# PARTICLES
O.bodies.append([
# 0
sphere(center=(0,0,0),radius=0.5,fixed=True,material='mat1'),
# 1
sphere((0,0,1),radius=0.5,material='mat1')
])

nbSpheres = 2;

# BOUNDARY CONDITION
O.bodies[1].state.blockedDOFs = "xyzXYZ"
O.bodies[1].state.vel = (0.0001,0,0)
O.bodies[1].state.angVel = (0,0,0)
# FUNCTIONAL COMPONENTS
# 1
def addPlotData1():
plot.addData(t=O.time,DX1 = O.bodies[1].state.pos[0],DY1 = O.bodies[1].state.pos[1],DZ1 = O.bodies[1].state.pos[2], FX1 = O.forces.f(1)[0], FY1 = O.forces.f(1)[1], FZ1 = O.forces.f(1)[2])

# 2
def addPlotData2():
plot.saveDataTxt('/home/nabid/Desktop/SCRIPTS/dataFile.txt',vars=('t','DX1','DY1','DZ1','FX1','FY1','FZ1'))

#PRINTING
from yade import plot
from pprint import pprint

# SIMULATION LOOP
O.engines=[
PyRunner(command='addPlotData1()',iterPeriod=1),
PyRunner(command='addPlotData2()',iterPeriod=1),
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR)],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
NewtonIntegrator(damping=0.1),
]

#TIME STEP
O.dt=0.5e-4*PWaveTimeStep()

#PLOTTING
plot.plots={'DX1':('FX1'),'DX1':('FZ1')}
plot.plot()

# SAVE SIMULATION

## Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
mrhappy
Solved:
2020-01-30
Last query:
2020-01-30
Last reply:
2020-01-23
 Jan Stránský (honzik) said on 2020-01-23: #1

Hello,

> uniaxial force

usually it is called normal (force, direction, component...) rather than uniaxial

1) because Ig2_Sphere_Sphere_ScGeom computes normal force according to penetration depth, which is
penetrationDepth = radius1 + radius2 - distanceOfCenters
and since distance of centers changes, penetrationDepth and thus normal force changes

3+4) your Matlab model is probably "more simple" (does not take into account the actual centers' distance, only its normal projection), basically based on different assumptions/formulation. Therefore you get different results.

> I read [1], but I want something with more detail.

read not only a short part, but rather the whole section Kinematic variables [3]. IMO there is enough details. I am not sure if some papers give more details.
Anyway, the source code is the best source how to investigate what really is computed :-) normal component [4] and shear components [5] (not only the referenced line, but also the code around).

cheers
Jan

 mrhappy (mrhappy) said on 2020-01-28: #2

Hi,

I read [1] very carefully. I dont understand one of the formula (maybe there is a typo). Under the Contact Model (exmaple) section, the formula with an if statement for calcuating FT. If the if statement is true, then FT = FtT*norm(FN)*tan(angle)/FtT. Now FtT (in the denominator) is a vector. How can you divide with a vector? Can someone clearfy this for me?

Thanks!

 Jan Stránský (honzik) said on 2020-01-28: #3

it is a typo, should be norm of the trial shear force (as in the source code [6])
Jan

 mrhappy (mrhappy) said on 2020-01-30: #4

Thank you so much for clearifying!

To post a message you must log in.