# Shear loading

I made a simple model to understand what’s going on. I have two spheres:

sphere(

sphere(

The following is the boundary condition:

O.bodies[

O.bodies[

O.bodies[

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:/

[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.

# PARTICLES

O.bodies.append([

# 0

sphere(

# 1

sphere(

])

nbSpheres = 2;

# BOUNDARY CONDITION

O.bodies[

O.bodies[

O.bodies[

# FUNCTIONAL COMPONENTS

# 1

def addPlotData1():

# 2

def addPlotData2():

plot.saveDataT

#PRINTING

from yade import plot

from pprint import pprint

# SIMULATION LOOP

O.engines=[

PyRunner(

PyRunner(

ForceReset

InsertionS

InteractionLoop(

[Ig2_

[Ip2_

[Law2_

),

GlobalStif

NewtonIntegrat

]

#TIME STEP

O.dt=0.

#PLOTTING

plot.plots=

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

Hello,

> uniaxial force

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

1) because Ig2_Sphere_

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/

> 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

[3] https:/

[4] https:/

[5] https:/

mrhappy (mrhappy) said : | #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(

Thanks!

Jan Stránský (honzik) said : | #3 |

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

Jan

[6] https:/

mrhappy (mrhappy) said : | #4 |

Thank you so much for clearifying!