New question #706413 on Yade:
https://answers.launchpad.net/yade/+question/706413
Thanks to Jean, I was able to come this far.
I am doing an analysis on a direct(simple) shear test.
I have created a program like below.
Why do "some" of the particles exit from the top of the facet when the plate is
used to compress the particles?
Is there a solution to this problem?
##############################################
from __future__ import print_function
from yade import pack,plot,polyhedra_utils,geom
from yade import export,qt
from yade.gridpfacet import *
import math
import numpy as np
import os
###############################################
#parameter (パラメーター in japanese)
###############################################
frictangle = np.radians(25)#45do
density = 3400.0
young = 3e8
gravity = (0.0, 0.0, 0)
velocity_cyl=2.5
velocity_cyl_idsCyl2=2.5
velocity_cyl_Box=2.5
###############################################
#Create sphere (モデルの作成 in japanese)
###############################################
mat_sp = CohFrictMat(alphaKr=.2,alphaKtw=.2,young = young, poisson =
0.35,frictionAngle=(frictangle),density=density)
O.materials.append(mat_sp)
pred=pack.inCylinder(centerBottom=(0,0,0.0),centerTop=(0,0,.02),radius=.03)
sp0=pack.randomDensePack(pred,spheresInCell=3000,radius=0.0004)
O.bodies.append(sp0)
###############################################
#Boundary creation
###############################################
idsCyl1 = O.bodies.append(geom.facetCylinder((0, 0, 0.0075), radius=.03,
height=.015, segmentsNumber=30, wallMask=6))
idsCyl2 = O.bodies.append(geom.facetCylinder((0, 0, 0.015), radius=.03,
height=.01, segmentsNumber=30, wallMask=5))
#idsCyl3 = O.bodies.append(geom.facetCylinder((0, 0, .0225), radius=0.10,
height=.005, segmentsNumber=100, wallMask=2))
radius1 = .03
radius2 = .15
fixBoxIds=[]
moveBoxIds=[]
upper_plate=box((0,0,-0.001),(0.06,0.06,0.001))
upper_plateID=O.bodies.append(upper_plate)
O.bodies[upper_plateID].state.blockedDOFs = 'xyXY' #xy:xyの方向角の固定、XY：XY変位の固定
upper_plate=box((0,0,-0.001),(0.06,0.06,0.001))
upper_plateID=O.bodies.append(upper_plate)
O.bodies[upper_plateID].state.blockedDOFs = 'xyXY' #xy:xyの方向角の固定、XY：XY変位の固定
#O.forces.addF(upper_plateID,(0,0,50e3),permanent=True)
O.forces.setPermF(upper_plateID,(0,0,100e3))
for i in range(0,365,5):
r = math.radians(i)
x1 = radius1 * math.cos(r)
y1 = radius1 * math.sin(r)
x2 = radius2 * math.cos(r-math.radians(5))
y2 = radius2 * math.sin(r-math.radians(5))
x3 = radius2 * math.cos(r+math.radians(5))
y3 = radius2 * math.sin(r+math.radians(5))
x4 = radius1 * math.cos(r-math.radians(5))
y4 = radius1 * math.sin(r-math.radians(5))
x5 = radius1 * math.cos(r+math.radians(5))
y5 = radius1 * math.sin(r+math.radians(5))
x6 = radius2 * math.cos(r)
y6 = radius2 * math.sin(r)
fixBoxIds.append(O.bodies.append(facet([(x1,y1,.015),(x2,y2,.015),(x3,y3,.015)])))
fixBoxIds.append(O.bodies.append(facet([(x4,y4,.015),(x5,y5,.015),(x6,y6,.015)])))
moveBoxIds.append(O.bodies.append(facet([(x1,y1,.015),(x2,y2,.015),(x3,y3,.015)])))
moveBoxIds.append(O.bodies.append(facet([(x4,y4,.015),(x5,y5,.015),(x6,y6,.015)])))
###############################################
#engines （制御 in japanese）
###############################################
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
# interaction loop
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(),
###TranslationEngine(translationAxis=[0,0,1],velocity=-1,ids=idsCyl3),
##TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_idsCyl2,ids=idsCyl2),
##TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_Box,ids=moveBoxIds),
#TranslationEngine(ids=idsCyl2,label="translation"),
#TranslationEngine(ids=moveBoxIds,label="translation2"),
#PyRunner(iterPeriod=1,command="setTranslation()"),
#PyRunner(iterPeriod=1,command="setTranslation2()"),
PyRunner(command='checkStress()', iterPeriod=10),
PyRunner(command='plotPlotData()',iterPeriod=1),
#PyRunner(command='checkStressNO()', iterPeriod=10),
]
###############################################
#Definition of distortion
###############################################
def checkStress1():
if O.bodies[idsCyl2[0]].state.displ()[0] > 0.015:
O.pause()
def checkStress2():
if O.bodies[idsCyl2[0]].state.displ()[0] > 0.015 and
O.bodies[idsCyl2[0]].state.displ()[1] > 0 :
O.pause()
###############################################
#Circle Control（円の制御 in japanese）
###############################################
def setTranslation():
x = O.bodies[idsCyl2[0]].state.displ()[0]
y = O.bodies[idsCyl2[0]].state.displ()[1]
angle = atan2(y,x)
a = Vector3(-sin(angle),cos(angle),0)
v = 50
translation.translationAxis = a
translation.velocity = v
def setTranslation2():
x = O.bodies[idsCyl2[0]].state.displ()[0]
y = O.bodies[idsCyl2[0]].state.displ()[1]
angle = atan2(y,x)
a = Vector3(-sin(angle),cos(angle),0)
v = 50
translation2.translationAxis = a
translation2.velocity = v
###############################################
#output（出力 in japanese）
###############################################
sp_num = len([b for b in O.bodies if isinstance(b.shape, Sphere)]) # 粒子の総数をカウント
print("number of spheres = ",sp_num)
def checkStress():
cylDisplacement_x = O.bodies[idsCyl2[0]].state.displ()[0]
cylDisplacement_y = O.bodies[idsCyl2[0]].state.displ()[1]
print('{:1.04f}'.format(cylDisplacement_x)," ",
'{:1.04f}'.format(cylDisplacement_y))
O.dt=0.5*PWaveTimeStep()
def plotPlotData():
pipe_force_x = sum([O.forces.f(i)[0] for i in idsCyl2])
pipe_force_y = sum([O.forces.f(i)[1] for i in idsCyl2])
cylDisplacementx = O.bodies[idsCyl2[0]].state.displ()[0]
cylDisplacementy = O.bodies[idsCyl2[0]].state.displ()[1]
cylDisplacementz = O.bodies[idsCyl2[0]].state.displ()[2]
Dz=upper_plate.state.displ()[2]
plot.addData(
i=O.iter,
disp=O.time*velocity_cyl,
force_x=pipe_force_x,
force_y=pipe_force_y,
Dz=Dz,
cylDisplacementx=cylDisplacementx,
cylDisplacementy=cylDisplacementy,
cylDisplacementz=cylDisplacementz,
)
plot.saveDataTxt('pipe_force.out',vars=('disp','force_x','force_y',
'Dz',
'cylDisplacementx','cylDisplacementx','cylDisplacementx'
))
plot.plots={
'disp':('force_x','force_y'),
'i':('Dz'),
'cylDisplacementx':('cylDisplacementy',)
}
plot.plot()
#plot.saveDataTxt(O.tags['id'] + '.txt')
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
