I made a simple model to understand what’s going on. I have two spheres:
The following is the boundary condition:
O.bodies.state.blockedDOFs = "xyzXYZ"
O.bodies.state.vel = (0.0001,0,0)
O.bodies.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  and  thoroughly. I understand that the tangential plane rotates as the spheres shear and that the shear force is calculated incrementally. Is deltaUs in  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 , but I want something with more detail.
4. My code is running well, so I assume that the my code is correct.

 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

# PARTICLES
O.bodies.append([
# 0
# 1
])

nbSpheres = 2;

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

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

#PRINTING
from pprint import pprint

# SIMULATION LOOP
O.engines=[
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:
Assignee:
No assignee Edit question
Solved by:
mrhappy
Solved:
2020-01-30
Last query:
2020-01-30
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
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 , but I want something with more detail.

read not only a short part, but rather the whole section Kinematic variables . 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  and shear components  (not only the referenced line, but also the code around).

cheers
Jan

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

Hi,

I read  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 )
Jan

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

Thank you so much for clearifying!