Cylindrical triaxial test using Pfacet

Hello,

1) Yes
2) triax.py
3) How do I achieve servo control?
4) [1]
    - what does "are destroyed" mean?
    - please provide a MWE
5) triax.py

Of course, triax.py is just a primitive approach - a cylindrical surface divided in triangles with fixed force.
It will not work if the surface is significantly distorted and it does not contain any "membrane effect".

You can use other shapes for the surface, maybe a bit overlapped?
Or PFacets to also simulate a membrane effect.
Or ...

Cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ServoPIDController

Dear Jan,

Thanks for your reply.

1) I'll check. Many thanks.
2) I'll check. Many thanks.
3) How do I achieve servo control? I am still not clear about this. I guess there is no displacement-based servo control in the code 'triax.py'.
4) This means the cylindrical surface is divided into triangles. Maybe you can try to run the code 'triax.py' and you will understand what I am talking about.
5) I'll check. Many thanks.

Cheers,
Kyle

> How do I achieve servo control? I am still not clear about this

e.g. using ServoPIDController [1]

> I guess there is no displacement-based servo control in the code 'triax.py'

Yes.
The loading is simply such that a constant velocity is prescribed (s.state.vel = ...)

>>> But the problem is that when implementing the code, facets generated using the command 'Pfacet' are destroyed
>> what does "are destroyed" mean? please provide a MWE
> Maybe you can try to run the code 'triax.py' and you will understand what I am talking about.

No, I will not understand, since there are no PFacets in current triax.py.
Therefore I ask for a MWE and better description of "are destroyed" (anyway useful even with the code).

Cheers
Jan

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

>> 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?

Define such velocity which and whose effect on the sample satisfy your definition of "quasi-equilibrium compression".

> Can you show me the code?

[2]

> I am sorry for my stupidity.

stupidity is not the right word here.
"Lack of experience" fits better

>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.

the same, what does "will break" mean?
That they do not touch each other?

> 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.

I do not use PFacet myself, so I cannot help here..
But thanks for description, maybe you can open a new more specific question (here the description is maybe too "deep" in the thread)

Cheers
Jan

[2] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/concrete/triax.py#L73

Hi,

> "My problem is that I just finished the part about generating flexible membranes using Pfacet and don't know what to do next."

I did not read all the details but if you are at this step I suppose "what to do next" is to assign forces to the nodes of the membranes to reflect external pressure. It certainly needs to account for some area and local orientation.
I never tried such thing myself so I'm not 100% sure what to expect. It should work but it probably needs to update the external forces after large deformations.

Cheers
Bruno

Hi! Jan

Thanks for your reply.
>the same, what does "will break" mean? That they do not touch each other?
The particles went out of the cylinder and flew all around [1].

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

Cheers,
Kyle

Can someone help me with servo-control part? Many thanks.

Have you tried ServoPIDController [1]?
If yes, why it did not help?
Cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ServoPIDController

Hi, Jan.

>Have you tried ServoPIDController?
I try to use it but can't find any tutorials. In PFC [1], the velocity value of the wall is determined by multiplying a factor by the stress difference (target stress minus current stress). This value is updated regularly. I am not sure whether the principle of the ServoPIDController is the same as that in PFC.

Looking forward to your reply.
Kyle

[1]https://docs.itascacg.com/pfc700/pfc/pfcmodule/doc/manual/wall_manual/wall_commands/cmd_wall_servo.html?highlight=servo#command:wall.servo

> I am not sure whether the principle of the ServoPIDController is the same as that in PFC

yes, it looks pretty similar [4]

Cheers
Jan

[4] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/common/KinematicEngines.cpp#L201

Hi, Jan.

>yes, it looks pretty similar [4]
Could you please provide an example of using ServoPIDController? Because there is no tutorial about how to use it. I am inclined to think that this command's explanation is not easy to understand without an example. Many thanks.

Looking forward to your reply.
Kyle

There is an example [5].
you can play with iterPeriod (try 10 or 1 instead of 1000) and/or other parameters.
The referenced source code [4] may be useful, as well as documentation itself [6] and links therein.
Cheers
Jan

[5] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/PIDController.py
[6] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ServoPIDController

You may also consider creating a new question aiming to the servo control, as here it might be already too "deep" in the thread.
Also, it is a separate problem and as such should have its own question [7]
Cheers
Jan

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