Servo control of flexible membrane using ServoPIDController

Hello,

> don't know how to use the ServoPIDController to achieve servo control of the flexible membrane with confining pressure.

To achieve confining pressure, probably the easiest approach is to assign corresponding center-oriented forces to the grid nodes.

For the start, I **personally** would not use servo control, I **personally** would start with a simpler approach - prescribe force to the nodes directly using O.forces.setPermF.

If you insist on servo control, just use ServoPIDController for each grid node.
Something like (not tested):
###
confiningStress = ... # N/m2

def node2servo(node):
    x,y,z = node.state.pos
    area = ... # cylinder surface / number of nodes ?
    f = confiningStress * area
    return ServoPIDController(ids=[node.id],axis=axis,iterPeriod=...,target=f,...)

nodes = [b for b in O.bodies if isinstance(b.shape,GridNode)]
servos = [node2servo(node) for node in nodes]
O.engines = O.engines + servos
###

Cheers
Jan

Hi! Jan,

Thanks for your reply. It is really helpful. But I still have two questions about this function:

1) return ServoPIDController(ids=[node.id],axis=axis,iterPeriod=...,target=f,...)
for the 'axis' mentioned in PIDController, should it be defined individually for each node like this?

confiningStress=5e4
cen_x = x_cyl # x coordinate of the centre of the crosssection
cen_y = y_cyl # y coordinate of the centre of the crosssection
def node2servo(node):
    x,y,z = node.state.pos
    area = 2* pi*r_cyl*h_cyl/len(nodes) # cylinder surface / number of nodes ?
    f = confiningStress * area
    axis = vector(x-cen_x, y-cen_y, z)
    return ServoPIDController(ids=[node.id],axis=axis,iterPeriod=100,target=f)

2) If I want to monitor the average stress on the flexible membrane, what should I do? My idea is that iterate all nodes and sum up the force in the direction of the connection line between the centre of the crosssection and the node, and then divide it by the area of the flexible membrane. If so, how to get the force in the direction of the connection line between the centre of the crosssection and the node for an individual node?

Cheer,
Kyle

> for the 'axis' mentioned in PIDController, should it be defined individually for each node

yes

> like this?

in general yes, with a few problems:

> axis = vector(x-cen_x, y-cen_y, z)

z should be 0

> area = 2* pi*r_cyl*h_cyl/len(nodes) # cylinder surface / number of nodes ?

The question mark in my comment is meant seriously.
For regular grid it is probably OK.
Also consider that the "boundary" nodes (bottom and top) maybe should have half force then "inner" nodes?

2) ...

the main idea sounds OK.

stress is backward what is used
stress = forceMagnitude / area

forceMagnitude = forceVector.dot(directionUnitVector)
forceVector = O.forces.f(body.id)

the direction vector may be as first iteration equal to axis from servo controller
x,y,z = node.state.pos
directionUnitVector = Vector3(x,y,0).normalized()

Cheers
Jan

Hi! Jan

Thanks for your reply. Your answer is rigorous and helpful.

> area = 2* pi*r_cyl*h_cyl/len(nodes) # cylinder surface / number of nodes ? The question mark in my comment is meant seriously. For regular grid it is probably OK. Also consider that the "boundary" nodes (bottom and top) maybe should have half force then "inner" nodes?

Yes. That's a really critical problem that needed to be clarified. What confused me currently is that I found the number of nodes is not consistent with what I expected. For example, in the beginning, we defined the number of facets by the following code [1]:
nw=40 # number of facets along the perimeter
nh=25 # number of facets along the height
nodes = [b for b in O.bodies if isinstance(b.shape, GridNode)]
print(len(nodes))

If we print the number of nodes, we will get 8000 where my initial suppose is 40*26 = 1040. That indicates that the ID of nodes forming Pfacet are independent of each other. In other words, more than one node can be created at the same position. So, I am wondering if maybe it is correct to use ' confining Stress * cylinder surface / number of nodes' to represent the average force applied on each node? We don't need to classify the top and bottom nodes anymore as the force applied on nodes is the same. It is identical to obtaining the node force by calculating the force applied on each Pfacet and distributing it into three nodes.

[1] https://answers.launchpad.net/yade/+question/681410

> the number of nodes is not consistent with what I expected

You have to question whether expectation or realization is wrong.
Here the expectation seems correct, so the there is a problem in realization.

> more than one node can be created at the same position

Of course then can have same position, but they should not have same position.
The purpose of flexible membrane (or using PFacets in general) is that the PFacets are sharing the nodes.

> it is correct to use ' confining Stress * cylinder surface / number of nodes' to represent the average force applied on each node?

Yes, it is correct to use the formula for "average force applied on each node".
However, if you have more nodes in one position, the result is meaningless..
So first you have to fix the usage of nodes.

> for r in range(nw):
> for h in range(nh):
> ...
> V1=(O.bodies.append(gridNode(v1,rNode,...)))
> V2=(O.bodies.append(gridNode(v2,rNode,...)))
> V3=(O.bodies.append(gridNode(v3,rNode,...)))
> V4=(O.bodies.append(gridNode(v4,rNode,...)))

The problem is here.
You need first to create nodes, one (not four) per each nw-nh cycle.
That way your nodes would not be duplicated.
Then you use these unique nodes to create gridConnections and PFacets.

Cheers
Jan

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

> axis = Vector3(x-x_cyl, y-y_cyl, 0)

probably it should be opposite, Vector3(x_cyl-x, y_cyl-y, 0) to point inside the cylinder (not outside)?
Currently the membrane is stretched rather then compressed (tested with maxVelocity=100,iterPeriod=1)

And maybe different parameters kP, kI, kD then default 1.0?

I don't know, then the behavior is also strange, the membrane is pulled inward once there is node-sphere interaction...

Cheers
Jan

Thanks for your reply, Jan!

> probably it should be opposite, Vector3(x_cyl-x, y_cyl-y, 0) to point inside the cylinder (not outside)?
Currently the membrane is stretched rather then compressed (tested with maxVelocity=100,iterPeriod=1)
> That's true. We have to reverse the direction vector.

>And maybe different parameters kP, kI, kD then default 1.0?
> Yes. But the difficulty is that it's hard to determine proper proportional parameters for kP, kI, and kD.

> I don't know, then the behavior is also strange, the membrane is pulled inward once there is node-sphere interaction...
> Yes. I also encounter this problem. I try to reduce the maximum velocity of servo control but still failed, only delaying the time of membrane destruction.

Actually, the PFacet is demonstrated to be effective in modelling flexible membranes successfully [1]. The problem is that there are some critical parameters such as velocity are not mentioned in this paper, which increases the difficulty for descendants to replicate the work.

[1] https://doi.org/10.1061/(ASCE)GM.1943-5622.000134

Thanks for your reply, Jan!

I'd also like to ask whether the function for calculating the stress on the loading plate and confining stress is correct or not:

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)

Many thanks
Kyle

> Yes. But the difficulty is that it's hard to determine proper proportional parameters for kP, kI, and kD.

I have no personal experience.
I suggest reading the wikipedia article and then testing the parameters on a simple (e.g. one node-sphere interaction or a few PFacets - a few spheres) example

> I'd also like to ask whether the function for calculating the stress on the loading plate and confining stress is correct or not:

yes, it looks OK

> dirV = Vector3(x-x_cyl,y-y_cyl,0).normalized()

maybe also reverse the axis, but is just changed sign..

Cheers
Jan