Question #707148 on Yade changed:
https://answers.launchpad.net/yade/+question/707148
Status: Needs information => Open
Ruidong LI gave more information on the question:
Hi! Jan
Thanks for your reply.
> The loading is simply such that a constant velocity is prescribed
> (s.state.vel = ...)
But which velocity value should I define to achieve quasi-equilibrium
compression? Can you show me the code? I am sorry for my stupidity.
>Therefore I ask for a MWE and better description of "are destroyed" (anyway
>useful even with the code).
Actually, what I mean is that when you run the 'triax' py without any change,
the triangle facets will break. And my current problem is not with the 'triax'
py.
My problem is that I just finished the part about generating flexible membranes
using Pfacet and don't know what to do next. What I want to achieve is
isotropic compression with a confining pressure equal to 50 kPa. And then,
maintain the confining pressure in the x and y directions and apply a
displacement-controlled movement of the top wall.
my MWE is presented as below
################################
### 1 ###
### LOAD INITIAL PACKING ###
################################
dir_inipacking = '/home/kyle/Yade/install/bin/lentils/xml/iniPacking.xml'
O.load(dir_inipacking)
# Define engine
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),
Bo1_PFacet_Aabb(),],
sortThenCollide=True),
InteractionLoop(
[
Ig2_GridNode_GridNode_GridNodeGeom6D(),
Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
Ig2_Sphere_Sphere_ScGeom6D(),
Ig2_Sphere_PFacet_ScGridCoGeom(),
],
[
Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
Ip2_FrictMat_FrictMat_FrictPhys()
],
[
Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
Law2_ScGeom_FrictPhys_CundallStrack(),
Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),
]
),
GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8,label='ts'),
NewtonIntegrator(gravity=(0,-9.81,0),damping=0.1,label='newton')
]
# encoding: utf-8
from yade import qt
from yade.gridpfacet import *
import math
############################################
### 2 ###
### DEFINING VARIABLES AND MATERIALS ###
############################################
# The following lines are used parameter definitions
readParamsFromTable(
# directory
dir_vtk = '/home/kyle/Yade/install/bin/lentils/vtk/20', # directory of
storing vtk files
dir_txt = '/home/kyle/Yade/install/bin/lentils/txt/20/Num20_servo.txt',
# directory of storing txt files
dir_xml = '/home/kyle/Yade/install/bin/lentils/xml/20', # directory of
storing xml files
# material parameters
young_w = 5e9, # Young's modulus of plates and walls
young_g = 1.8e6, # Young's modulus of grains [1]
den_w = 7874, # density of plates and walls
den_g = 980, # density of grains
poi_w = 0.3, # possion's ratio of plates and walls
poi_g = 0.25, # possion's ratio of grains
friAngle_w = 0.5,#radians(15), # friction angle of plates and walls
friAngle_g = 0.5,#radians(29), # friction angle of grains
# Parameters of cylindrical walls
x_cyl = 0.0547, # x-coordinate of center of cylindrical walls
y_cyl = 0.0535, # y-coordinate of center of cylindrical walls
z_cyl = 0, # z-coordinate of center of cylindrical walls
r_cyl = 0.0358, # radius of the cylinder
h_cyl = 0.14, # height of the cylinder
# confining pressure
preStress = -50000,
# axial strain rate
strainRate=-100
)
from yade.params.table import *
########################################
### 3 ###
### CREATE TOP AND BOTTOM PLATES ###
########################################
# erase rigid cylindrical facets
for b in O.bodies:
if isinstance(b.shape,Facet): O.bodies.erase(b.id)
# create top and bottom plates
h_cyl = max([b.state.pos[2] + b.shape.radius for b in O.bodies if
isinstance(b.shape, Sphere)])
top_plate =
O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h_cyl),radius=r_cyl*1.2,height=0,segmentsNumber=40,color=(1,1,1)))
bottom_plate =
O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,0),radius=r_cyl*1.2,height=0,segmentsNumber=40,color=(1,1,1)))
'''
vel = strainRate * (h - rParticle * 2 * bcCoeff)
for s in bot:
s.shape.color = (1, 0, 0)
s.state.blockedDOFs = 'xyzXYZ'
s.state.vel = (0, 0, -vel)
for s in top:
s.shape.color = Vector3(0, 1, 0)
s.state.blockedDOFs = 'xyzXYZ'
s.state.vel = (0, 0, vel)
'''
######################################
### 4 ###
### GENERATE FLEXIBLE MEMBRANE ###
######################################
# Materials
O.materials.append(CohFrictMat(young=1e7,poisson=1,density=2650,frictionAngle=radians(30),normalCohesion=3e7,shearCohesion=3e7,momentRotationLaw=True,label='gridNodeMat'))
O.materials.append( FrictMat(
young=1e7,poisson=0.1,density=2650,frictionAngle=radians(30),label='pFacetMat'
))
# Generate flexible membrane
nodesIds=[]
cylIds=[]
pfacets=[]
width=2*r_cyl #diameter of cylinder
height=h_cyl #height of cylinder
nw=40 # no of nodes along perimeter
nh=25 # no of nodes along height
rNode=width/100 #radius of grid node
color1=[255./255.,102./255.,0./255.]
color2=[0,0,0]
color3=[248/255,248/255,168/255]
color4=[156/255,160/255,98/255]
#color3=[165/255,245/255,251/255]
#color4=[101/255,153/255,157/255]
rCyl2 = 0.5*width / cos(pi/float(nw))
vCenter = Vector3(x_cyl, y_cyl, 0)
for r in range(nw):
for h in range(nh):
v1 =
Vector3(rCyl2*cos(2*pi*(r+0)/float(nw)),rCyl2*sin(2*pi*(r+0)/float(nw)),
height*(h+0)/float(nh))+vCenter
v2 =
Vector3(rCyl2*cos(2*pi*(r+1)/float(nw)),rCyl2*sin(2*pi*(r+1)/float(nw)),
height*(h+0)/float(nh))+vCenter
v3 =
Vector3(rCyl2*cos(2*pi*(r+1)/float(nw)),rCyl2*sin(2*pi*(r+1)/float(nw)),
height*(h+1)/float(nh))+vCenter
v4 =
Vector3(rCyl2*cos(2*pi*(r+0)/float(nw)),rCyl2*sin(2*pi*(r+0)/float(nw)),
height*(h+1)/float(nh))+vCenter
V1=(O.bodies.append(gridNode(v1,rNode,wire=False,fixed=False,
material='gridNodeMat',color=color2)) )
V2=(O.bodies.append(gridNode(v2,rNode,wire=False,fixed=False,
material='gridNodeMat',color=color2)) )
V3=(O.bodies.append(gridNode(v3,rNode,wire=False,fixed=False,
material='gridNodeMat',color=color2)) )
V4=(O.bodies.append(gridNode(v4,rNode,wire=False,fixed=False,
material='gridNodeMat',color=color2)) )
# append node ids to seperate matrix for later use
nodesIds.append(V1)
nodesIds.append(V2)
nodesIds.append(V3)
nodesIds.append(V4)
#create grid connection
cylIds.append(O.bodies.append(gridConnection(V1,V2,rNode,material='gridNodeMat',color=color3)))
cylIds.append(O.bodies.append(gridConnection(V2,V3,rNode,material='gridNodeMat',color=color3)))
cylIds.append(O.bodies.append(gridConnection(V1,V3,rNode,material='gridNodeMat',color=color3)))
cylIds.append(O.bodies.append(gridConnection(V3,V4,rNode,material='gridNodeMat',color=color3)))
cylIds.append(O.bodies.append(gridConnection(V4,V1,rNode,material='gridNodeMat',color=color3)))
#create Pfacets
pfacets.append(O.bodies.append(pfacet(V1,V2,V3,wire=True,material='pFacetMat',color=color3)))
pfacets.append(O.bodies.append(pfacet(V1,V3,V4,wire=True,material='pFacetMat',color=color4)))
# calculate void ratio
def myClumpVoidRatio(VolTot, density):
mass=sum([ b.state.mass for b in O.bodies if b.isClump])
VolSolid = mass/density
poro = (VolTot-VolSolid)/VolSolid
return poro
h_cyl = max([b.state.pos[2] + b.shape.radius for b in O.bodies if
isinstance(b.shape, Sphere)])
Vol = math.pi * h_cyl * r_cyl * r_cyl
poro_ini = myClumpVoidRatio(Vol, den_g)
print('Initial void ratio:', str(poro_ini))
--
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 : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
