2D simulation with clump using PeriTriaxController

Asked by nie jiayan on 2019-03-19

Hello, everyone. Now I want to do 2D DEM simulation with clumps in periodic boundary condition. But some confusions exits. Why when I run the simulation, the interactions between the constituent balls in clumps are also active. And the PeriTriaxController() failed. Anyone can help me? Many thanks in advance.
The following are my codes.

from yade import pack, qt, plot
import time
import numpy as np

############################################
### Read external file information ###
############################################

filename1='pebblex.txt'
filename2='pebbley.txt'
filename3='pebbler.txt'
filename4='clumpn.txt'

pbx=np.loadtxt(filename1)
pby=np.loadtxt(filename2)
pbr=np.loadtxt(filename3)
cmn=np.loadtxt(filename4)

cn=cmn.size-1

## create materials for initial consolidation
O.materials.append(FrictMat(density=2650,young=6.0e8,poisson=0.8,frictionAngle=0.01))

## process the external file information to get the clump configuration list

for ii in range(0,cn):
    aa=int(np.sum(cmn[:ii+1]))
    bb=int(np.sum(cmn[:ii+2]))
    bodyList=[]
    for jj in range(aa,bb):
        bodyList.append(O.bodies.append(sphere([pbx[jj],pby[jj],0.05],pbr[jj])))
    O.bodies.clump(bodyList,discretization=10)

############################################
### DEFINING VARIABLES AND BLOCK DEGREES ###
############################################

ts = time.time()
pressure = -1e5
size=0.30

## setup the periodic boundary
O.periodic=True
O.cell.hSize = Matrix3(size,0,0, 0,size,0, 0,0,0.1)

for p in O.bodies:
        p.state.blockedDOFs = 'zXY'

print(len(O.bodies))
############################
### DEFINING ENGINES ###
############################

triax=PeriTriaxController(
  # specify target values and whether they are strains or stresses
  goal=(pressure,pressure,0),stressMask=3,
  # type of servo-control
  dynCell=True,maxStrainRate=(0.5,0.5,0),
  # wait until the unbalanced force goes below this value
  maxUnbalanced=0.001,relStressTol=1e-3,
                # call this function when goal is reached and the packing is stable
  doneHook='consolidationFinished()'
 )

newton=NewtonIntegrator(damping=.1)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),],
  [Ip2_FrictMat_FrictMat_FrictPhys(),],
  [Law2_ScGeom_FrictPhys_CundallStrack(),]
 ),
 triax,
 newton,
        PyRunner(command='addPlotData()',iterPeriod=100),
]

O.dt=.5*PWaveTimeStep()

def addPlotData():

   pp = utils.porosity() #this is the porosity of the cell.
   ee = pp / (1-pp) #this is the void ratio of the 3D cell.
   a = [i for i in O.bodies if isinstance(i.shape, Clump) and i.intrs()==[]]
   plot.addData(unbalanced=unbalancedForce(),i=O.iter,
                e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
  s11=-triax.stress[0],
  s22=-triax.stress[1],
  s33=-triax.stress[2],
                void = ee, floatPt=len(a),
   )

# define what to plot//since they have the same x (i), the latter i should have a space in it 'i '
plot.plots={'i':('unbalanced',),'i ':('s11','s22','s33'),' i':('e11','e22','e33'),
   ' i ':('void',),'i ':('floatPt',),
   #each 'i' should have different number or location of space. Or the latter will overwrite the former one.
}
# show the plot
plot.plot()

yade.qt.Controller(), yade.qt.View()

def consolidationFinished():
        # set the current cell configuration to be the reference one
 O.cell.trsf=Matrix3.Identity
 print 'consolidationFinished'
        print('time: '+str((time.time()-ts)/60.)+'min')
        O.save('dense_packing_100kPa.yade.gz')
 O.pause()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
nie jiayan
Solved:
2019-03-22
Last query:
2019-03-22
Last reply:
nie jiayan (jiayann) said : #1

Dear all, I rechecked my codes and found that using O.bodies.clump() method to generate clumps, then the interactions between spheres within the clumps also have normal forces. So when I use PeriTriaxController() servo, the initial stress sigma_x sigma_y, also exits. So I wonder this is unreasonable. I have scanned some similar questions in launchpad, and someone said that there is no interactions between those spheres within clump.
 Ps, I also run the manual example with 06-periodic-triaxial-test.py , it is ok. I am afraid if some details I miss.
Thanks for your help. Regards, Jiayan.

Robert Caulk (rcaulk) said : #2

Tips to get improved replies:

>>And the PeriTriaxController() failed

Describe the "failure", include copies of the error message etc.

>> I have scanned some similar questions in launchpad, and someone said that there is no interactions between those spheres within clump

Reference the questions by link

nie jiayan (jiayann) said : #3

Hello, Robert,

    Thanks for your reply and tips.

    Actually, for the first one, I have generated the clumps successfully within the periodic cell, which do not contact with each other. so when I run the sample using triax=PeriTriaxController(). it is reasonable that values of triax.stress increase from zero to the target value, while it is not that in my simulation. It begins from a large value, and that is strange, so I checked the contact force of the overlapping spheres within clump, I found that exists. So I do not figure it out. Why is that in my code?

    As for the second one, I reference this link that there is no interaction within the clumps.
    https://answers.launchpad.net/yade/+question/403339

    Regards,
    Jiayan

nie jiayan (jiayann) said : #4

After rebuilding the source code, the problem has been solved.