wall friction

Asked by William

Hi,
I have set the friction angle of the left, right and top walls to 0, but why are the left and right walls still subjected to vertical force and the top wall still subjected to horizontal force in the second script.
I have two scripts and they should run in sequence.

The first script:

#########################################
### Defining parameters and variables ###
#########################################

#Material constants
Density = 3000
FrictionAngle = 50
PoissonRatio = 0.5
Young = 300e6
Damp = 0.5
AvgRadius = 0.00006
N_particles = 10000

#Wall constants
WDensity = 0
WFrictionAngle = 0.0
WPoissonRatio = 0.5
WYoung = 50e9

#Packing variables
mn = Vector3(0,0,0)
mx = Vector3(0.1,0.020,0.001)

#Confining variables
ConfPress1 = -1e5 #pre-compression
ConfPress = -1.0e5

#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT

########################################
#import necessary packages
from yade import pack,plot,os,timing
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

########################################
### Sample Preparing ###################
########################################

#Create materials for spheres and plates
SphereMat = O.materials.append(FrictMat(young = Young, poisson = PoissonRatio, frictionAngle = radians(FrictionAngle), density = Density))
WallMat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle), density = WDensity))

#Create walls for packing
wallIds = O.bodies.append(aabbWalls([mn,mx],thickness=0.000001,material=WallMat))

#Use SpherePack object to generate a random loose particles packing
#O.periodic = True
#O.cell.setBox(8,3,8)

sp = pack.SpherePack()

#psdSizes,psdCumm=[0.00142,0.00151,0.00160,0.00169,0.00179,0.00188,0.00197,0.00207,0.00216,0.00225,0.00235,0.00244,0.00253,0.00263,0.00272,0.00281],[0,0.067,0.133,0.2,0.267,0.333,0.4,0.467,0.533,0.6,0.667,0.733,0.8,0.867,0.933,1]

sp.makeCloud(Vector3(0,0,0.0005),Vector3(0.1,0.020,0.0005),num=N_particles,rMean=AvgRadius,rRelFuzz=0.33,distributeMass=False)

sp.toSimulation(material = SphereMat)

O.usesTimeStepper=True
O.trackEnerty=True
#################################
#####Defining triaxil engines####
#################################

###first step: compression#######
triax1=TriaxialStressController(
    wall_back_activated = False, #for 2d simulation
    wall_front_activated = False,
    thickness = 0.000001,
    maxMultiplier=1.+1.5e5/Young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+4e3/Young,
    #maxMultiplier = 1.002,
    internalCompaction = True, # If true the confining pressure is generated by growing particles
    #max_vel = 1.5,
    stressMask = 7,
    computeStressStrainInterval = 10,

    goal1 = ConfPress1,
    goal2 = ConfPress1,
    #goal3 = ConfPress,
)

newton=NewtonIntegrator(damping=Damp)

###engine
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_FrictPhys()],
     [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax1,
    newton,
    PyRunner(realPeriod=10,command='checkUnbalanced()',label='check'),
]
# Simulation stop conditions defination
def checkUnbalanced():
    unb=unbalancedForce()
    mStress = (triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.0
    s1 = triax1.stress(triax1.wall_right_id)[0]
    s2 = triax1.stress(triax1.wall_top_id)[1]
    print(unb,s2,s1)
    if unb<0.01 and abs(ConfPress1-s2)/abs(ConfPress1)<0.01:
       O.pause()
       O.save('initial.yade.bz2')

The second script:

#Wall constants
WDensity = 0
WFrictionAngle = 0.0
WPoissonRatio = 0.5
WYoung = 50e9
Mat = O.materials.append(FrictMat(young = WYoung, poisson = WPoissonRatio, frictionAngle = radians(WFrictionAngle), density = WDensity))
FrictionAngle = 35

Damp = 0.25

#Confining variables
force = -10

#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT

#restart
O.load('initial.yade.bz2')

########################################
#import necessary packages
from yade import pack,plot,os,timing

bodylist = []
for i in range(0,165):
    bodylist.append(O.bodies.append(sphere([-0.060+0.001*i,-0.0005,0.0005],0.0005,material=Mat)))
bottomwallclump = O.bodies.clump(bodylist)

for b in range(6):
    O.bodies.erase(b)

id_box1 = O.bodies.append(box((-0.0005,0.01,0.0005),(0.0005,0.01,0.0005),fixed=True,material=Mat))# left wall
id_box2 = O.bodies.append(box((0.1005,0.01,0.0005),(0.0005,0.01,0.0005),fixed=True,material=Mat))# right wall
id_box3 = O.bodies.append(box((0.050,0.0205,0.0005),(0.050,0.0005,0.0005),fixed=False,material=Mat,color=(1,0,1)))# top wall
#id_box4 = O.bodies.append(box((51,-0.5,0.5),(70,0.5,0.5),fixed=True,material=Mat))# bottom wall

print('bottom',O.bodies[-4].state.pos)
print('top',O.bodies[-1].state.pos)
print('right',O.bodies[-2].state.pos)
print('left',O.bodies[-3].state.pos)

topwall = O.bodies[id_box3]
topwall.state.blockedDOFs = 'xzXYZ'
O.bodies[-4].state.blockedDOFs = 'xyzXYZ'
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.shape.color=(0.,0.,1.)

for i in range(-164,-1):
    O.bodies[i].shape.color=(1.,0.,1.)

O.forces.setPermF(O.bodies[-1].id,(0,force,0))

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_FrictPhys()],
     [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    NewtonIntegrator(damping=Damp),
    PyRunner(realPeriod=2,command='endCheck()',label='check')
# PyRunner(realPeriod=10,command='checkUnbalanced()',label='check'),
]

# Simulation stop conditions defination
def endCheck():
    unb=unbalancedForce()
    topstress = (10-O.forces.f(O.bodies[-1].id)[1])/0.0001
    print(unb,topstress)
    if unb<0.001 and (abs(topstress)-1e5)/1e5 < 0.0001:
        if O.iter<60000:
            return
        else:
            O.pause()
            setContactFriction(radians(FrictionAngle))
            O.save('Compression.yade.bz2')

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

> ... in the second script.

it depends on the definition of "in ... script" and how you run the script(s) and when you investigate the forces.
If e.g. you investigate the forces after
setContactFriction(radians(FrictionAngle))
then friction s not 0 any more.

Cheers
Jan

Revision history for this message
William (qfxx-123) said :
#2

Hi, Jan.

The program stopped before I changed the friction angle. That's to say even if i delete 'setContactFriction(radians(FrictionAngle))', the left and right walls still subjected to vertical force and the top wall still subjected to horizontal force.

Xujin

Revision history for this message
Jan Stránský (honzik) said :
#3

> the top wall still subjected to horizontal force
>
> id_box3 = O.bodies.append(box(...))# top wall
> O.forces.setPermF(O.bodies[-1].id,(0,force,0))

you yourself prescribed horizontal force

Why using O.bodies[-1].id and not much more readable topwall.id or id_box3 ?
Similarly for other occurrences of O.bodies[-N]

> the left and right walls still subjected to vertical force

difficult to say without a MWE [1], M=minimal, both from code AND number of particles, running time etc. Please try to reduce your script.

You can check yourself the interactions of the side walls and the interactions properties (frictionAngle, contact point, forces, ...)
E.g. it is possible that the vertical force is not due to friction, but simply because some particles pushes the wall from bottom or from top.

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Can you help with this problem?

Provide an answer of your own, or ask William for more information if necessary.

To post a message you must log in.