Use of ServoPIDController for simulating cylindrical triaxial test

Asked by Ruidong LI

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 encountered with some problems. They are summarized as follows:

1) When implementing the MWE, the top plate is supposed to move, but actually does not. And the cylindrical walls are going to be destroyed. What causes this to happen and how to solve this?
2) I would like to apply confining pressure on cylindrical walls. How to achieve this?
3) For the ServoPIDController, do 'TranslationEngine', 'CombinedKinematicEngine' and 'ServoPIDController' need to be used together?

Many thanks.

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

# 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(),
   ]
  ),
     TranslationEngine(translationAxis=[0, 0, 1], velocity=-0.1, ids=top_plate, dead=False, label='translat'),
        CombinedKinematicEngine(ids=top_plate, label='combEngine', dead=True) +
        ServoPIDController(axis=[0, 0, 1], maxVelocity=0.1, iterPeriod=100, ids=top_plate, target=5.0e5, 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'),
  NewtonIntegrator(gravity=(0,-9.81,0),damping=0.1,label='newton')
]

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

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/PIDController.py

Question information

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

Hello,

> 1) ... plate is supposed to move, but actually does not ...

It actually does move.

Instead of
O.run()
I used
###
top0 = O.bodies[top_plate[0]]
print(top0.state.displ()) # Vector3(0,0,0)
O.run(2000,True)
print(top0.state.displ()) # Vector3(0,0,-1.999999998947288304e-06)
###
with non-zero result.

> 3) For the ServoPIDController, do 'TranslationEngine', 'CombinedKinematicEngine' and 'ServoPIDController' need to be used together?

No.
In the example it is combined with rotation, but for linear motion, you can use just ServoPIDController

> 1) ... And the cylindrical walls are going to be destroyed. What causes this to happen and how to solve this?
> 2) I would like to apply confining pressure on cylindrical walls. How to achieve this?

please focus on one problem per one question ([2], point 5).

Cheers
Jan

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

Revision history for this message
Ruidong LI (kyle2000) said :
#2

Hi, Jan

Thanks for your reply.
> It actually does move
> Yes, I can capture the movement After adding the 'Bo1_Facet_Aabb()' in the 'InsertionSortCollider' part of the engine.

> In the example it is combined with rotation, but for linear motion, you can use just ServoPIDController
OK. Thx.

> please focus on one problem per one question ([2], point 5).
> OK. I will open one new question and close this one. Many thanks.

Best,
Kyle