yade-users team mailing list archive

[Ruidong LI <question707278@xxxxxxxxxxxxxxxxxxxxx>]

Question #707278 on Yade changed:
https://answers.launchpad.net/yade/+question/707278

    Status: Answered => Open

Ruidong LI is still having a problem:
Thanks for your reply, Jan!

Your answer solved my problem. I rewrote the code for generating the
flexible membrane. Moreover, I used your suggested servo-control mode
for the top, bottom, and cylindrical walls. I tried to measure whether
the confining stress was reached or not. The result indicated that the
confining stress in the top and bottom direction reached the target
value while that in the cylindrical direction ly reached the half of
target value. I don't know where the problem is. The assignment of
cylindrical stress or measure?

The updated MWE is presented below:

from __future__ import print_function
from yade import pack, ymport, qt
import gts, os.path, locale, plot, random
import numpy as np
from yade.gridpfacet import *
import math

locale.setlocale(
        locale.LC_ALL, 'en_US.UTF-8'
) #gts is locale-dependend. If, for example, german locale is used, gts.read()-function does not import floats normally

############################################
###                  1                   ###
### DEFINING VARIABLES AND MATERIALS     ###
############################################

# The following lines are used parameter definitions
readParamsFromTable(

        # 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

)
from yade.params.table import *

# create materials for spheres and walls
wallMat = O.materials.append(FrictMat(young = young_w, poisson = poi_w, frictionAngle = friAngle_w, density = den_w, label = 'walls'))
sphereMat = O.materials.append(FrictMat(young = young_g, poisson = poi_g, frictionAngle = friAngle_g, density = den_g, label = 'spheres'))

###############################################
###                    2                    ###
###   IMPORTING GRAINS AND CREATING CLUMPS  ###
###############################################

# spheres
pred = pack.inCylinder((x_cyl, y_cyl, z_cyl), (x_cyl, y_cyl, h_cyl), r_cyl)
sp = SpherePack()
sp = pack.randomDensePack(pred, spheresInCell=2000, radius=0.005, returnSpherePack=True)
spheres = sp.toSimulation(color=(0, 1, 1), material=sphereMat)

# assign unique color for each sphere
currentSphereId = O.bodies[0].id
color = randomColor()
for b in O.bodies:
    if b.id != currentSphereId:#change color each time clumpId changes
        currentSphereId = b.id
        color = randomColor()
    if isinstance(b.shape,Sphere):#colorize spheres
        b.shape.color = color

# 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)])
h0_cyl = min([b.state.pos[2] - b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)])
idTplate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h_cyl),radius=r_cyl*1.2,material=wallMat,height=0,segmentsNumber=40,color=(1,1,1),wire=True))
idBplate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h0_cyl),radius=r_cyl*1.2,material=wallMat,height=0,segmentsNumber=40,color=(1,1,1)))

#########################
###          3        ###
###   DEFINE ENGINES  ###
#########################

#target confining stress
confiningStress = 5e4
global wAz
wAz = pi * (1.2*r_cyl) * (1.2*r_cyl)
axialforce = confiningStress * wAz

# Define engine
O.engines = [
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),
     Bo1_PFacet_Aabb(),
     Bo1_Facet_Aabb()],
     sortThenCollide=True),
  InteractionLoop(
   [
    Ig2_GridNode_GridNode_GridNodeGeom6D(),
    Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
    Ig2_Sphere_Sphere_ScGeom6D(),
    Ig2_Sphere_PFacet_ScGridCoGeom(),
    Ig2_Facet_Sphere_ScGeom(),
   ],
   [
          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(),
   ]
  ),
  NewtonIntegrator(gravity=(0,0,0),damping=0.7,label='newton'),
  ServoPIDController(axis=[0, 0, 1], maxVelocity=10, iterPeriod=100, ids=idTplate, target=axialforce, kP=1.0, kI=1.0, kD=1.0),
  ServoPIDController(axis=[0, 0, -1], maxVelocity=10, iterPeriod=100, ids=idBplate, target=axialforce, kP=1.0, kI=1.0, kD=1.0),
  PyRunner(command='addPlotData()', iterPeriod=100, label='graph'),
]

################################
###           3              ###
### CREATE FLEXIBLE MEMBRANE ###
################################

## Generate flexible membrane
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' ))

global height
pfacets=[]
width=2*r_cyl #diameter of cylinder
height=h_cyl #height of cylinder
nw=40 # number of facets along perimeter
nh=25 # number of facets 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]
rCyl2 = 0.5*width / cos(pi/float(nw))
vCenter = Vector3(x_cyl, y_cyl, 0)

# create gridnodes 
nodesIds=[]
for r in range(nw):
 for h in range(nh+1):
    v = Vector3(rCyl2*cos(2*pi*r/float(nw)),rCyl2*sin(2*pi*r/float(nw)), height*h/float(nh))+vCenter
    V = (O.bodies.append(gridNode(v,rNode,wire=False,fixed=False, material='gridNodeMat',color=color2)) )
    # append node ids to seperate matrix for later use
    nodesIds.append(V)
# create grid connection
nodeIDmin = min(nodesIds)
cylIdsV=[]
cylIdsH=[]
cylIdsI=[]
for r in range(nw):
   for h in range(nh):
      if r == nw-1:
         #determine grid node
         V1 = nodeIDmin+(r+0)*(nh+1)+h+0
         V2 = nodeIDmin+(r+0)*(nh+1)+h+1
         V3 = nodeIDmin+(0+0)*(nh+1)+h+1
         V4 = nodeIDmin+(0+0)*(nh+1)+h+0
      else:
         #determine grid node
         V1 = nodeIDmin+(r+0)*(nh+1)+h+0
         V2 = nodeIDmin+(r+0)*(nh+1)+h+1
         V3 = nodeIDmin+(r+1)*(nh+1)+h+1
         V4 = nodeIDmin+(r+1)*(nh+1)+h+0
      #create grid connection 
      cylIdsV.append(O.bodies.append(gridConnection(V1,V2,rNode,material='gridNodeMat',color=color3)))
      cylIdsH.append(O.bodies.append(gridConnection(V1,V4,rNode,material='gridNodeMat',color=color3)))
      cylIdsI.append(O.bodies.append(gridConnection(V1,V3,rNode,material='gridNodeMat',color=color3)))
      if h == nh-1:
         cylIdsH.append(O.bodies.append(gridConnection(V2,V3,rNode,material='gridNodeMat',color=color3)))
# create PFacet
for r in range(nw):
   for h in range(nh):
      if r == nw-1:
         #determine grid node
         V1 = nodeIDmin+(r+0)*(nh+1)+h+0
         V2 = nodeIDmin+(r+0)*(nh+1)+h+1
         V3 = nodeIDmin+(0+0)*(nh+1)+h+1
         V4 = nodeIDmin+(0+0)*(nh+1)+h+0
      else:
         #determine grid node
         V1 = nodeIDmin+(r+0)*(nh+1)+h+0
         V2 = nodeIDmin+(r+0)*(nh+1)+h+1
         V3 = nodeIDmin+(r+1)*(nh+1)+h+1
         V4 = nodeIDmin+(r+1)*(nh+1)+h+0
      #create Pfacets
      pfacets.append(O.bodies.append(pfacet(V1,V3,V2,wire=True,material='pFacetMat',color=color3)))
      pfacets.append(O.bodies.append(pfacet(V1,V4,V3,wire=True,material='pFacetMat',color=color4)))

wAc = pi * 2 * r_cyl * h_cyl
nAPFacet = wAc/(nw*nh)
def node2servo(node):
    global x_cyl
    global y_cyl
    global nAPFacet
    global height
    global confiningStress
    x, y, z = node.state.pos
    if z == 0 or z == height:
        area = nAPFacet * 0.5
    else:
        area = nAPFacet
    f = confiningStress * area
    axis = Vector3(x-x_cyl, y-y_cyl, 0)
    return ServoPIDController(ids=[node.id],axis=axis,iterPeriod=100,target=f,kP=1.0, kI=1.0, kD=1.0)
    
nodes = [b for b in O.bodies if isinstance(b.shape,GridNode)]
servos = [node2servo(node) for node in nodes]
O.engines = O.engines + servos

def addPlotData():
 # axial stress
 global wAz
 global wAc
 global x_cyl
 global y_cyl
 fwt = Vector3(0, 0, 0)
 for i in idTplate:
  fwt += O.forces.f(i)
 fwb = Vector3(0, 0, 0)
 for j in idBplate:
  fwb += O.forces.f(j)
 # cylindrical stress
 fcyl = 0
 for n in nodes:
  fnode = O.forces.f(n.id)
  x,y,z = n.state.pos
  dirV = Vector3(x-x_cyl,y-y_cyl,0).normalized()
  fnodeMagn = fnode.dot(dirV)
  fcyl += fnodeMagn
 plot.addData(z=O.iter, swt=fwt[2]/wAz, swb=fwb[2]/wAz, swc = fcyl/wAc)

qt.View()
r = qt.Renderer()
r.bgColor = 1, 1, 1

plot.plots = {'z': ('swt', 'swb', 'swc')}
plot.plot()
O.run()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.