yade-users team mailing list archive

[Ruidong LI <question707278@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi! I am trying to model a cylindrical triaxial test based on Pfacet and ServoPIDController. I would like to model a triaxial compression test with confining pressure equal to 50kPa. The bottom and top loading plates are created using the 'facet' element. The top one is assumed to move slowly in the vertical direction while the bottom one is fixed. The cylindrical walls are generated using 'Pfacet' elements. Referring to [1], I wrote the current MWE but don't know how to use the ServoPIDController to achieve servo control of the flexible membrane with confining pressure.

The 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)])
top_plate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h_cyl),radius=r_cyl*1.2,material=wallMat,height=0,segmentsNumber=10,color=(1,1,1),wire=True))
bottom_plate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,0),radius=r_cyl*1.2,material=wallMat,height=0,segmentsNumber=40,color=(1,1,1)))
#O.bodies[top_plate].state.vel = (0, 0, -0.01)
#O.bodies[bottom_plate].state.vel = (0, 0, 0)
#O.bodies[bottom_plate].state.blockedDOFs = 'xyzXYZ'

# 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,-9.81,0),damping=0.1,label='newton'),
        TranslationEngine(translationAxis=[0, 0, 1], velocity=-2, ids=top_plate, dead=False, label='translat'),
        CombinedKinematicEngine(ids=top_plate, label='combEngine', dead=True) +
        ServoPIDController(axis=[0, 0, 1], maxVelocity=2, iterPeriod=100, ids=top_plate, target=1.0e7, kP=1.0, kI=1.0, kD=1.0),
        PyRunner(command='addPlotData()', iterPeriod=1000, label='graph'),
        PyRunner(command='switchTranslationEngine()', iterPeriod=45000, nDo=2, label='switchEng'),
]

## 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' ))

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]
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)))

from yade import qt
qt.View()
r = qt.Renderer()
r.bgColor = 1, 1, 1

def addPlotData():
	fMove = Vector3(0, 0, 0)
	for i in top_plate:
		fMove += O.forces.f(i)
	plot.addData(z=O.iter, pMove=fMove[2], pFest=fMove[2])

def switchTranslationEngine():
	print("Switch from TranslationEngine engine to ServoPIDController")
	translat.dead = True
	combEngine.dead = False

plot.plots = {'z': ('pMove', 'pFest')}
plot.plot()
O.run()




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