yade-users team mailing list archive

[Leonard <question698594@xxxxxxxxxxxxxxxxxxxxx>]

Question #698594 on Yade changed:
https://answers.launchpad.net/yade/+question/698594

    Status: Answered => Open

Leonard is still having a problem:
Hello Jan and Luc,
Thanks very much for your reply.

Here I made two MWE for rigid wall case and periodic boundary case
respectively. It may take around 8 min for each of the simulation, or
you can also find the results of stress-strain and volumetric strain
evolution at [1].

>same strain rate (I think)
Because in rigid wall which uses TriaxialStressController, I can specify the loading rate by triax.goal. But as far as I understand, it seems that for PeriTriaxController, it can only specify the loading rate by maxStrainRate. So I made triax.strain=1 (rigid case) and triax.maxStrainRate = (1., 1., 1.) in periodic case. Do you know how to specify the strain rate in PeriTriaxController?

>the two cases use ... the same of sample size
Now I know they are not the same. Just explain here why I though they were same before. If you run the periodic MWE below, it prints out the cell size just before the deviatoric loading, this size is almost the same as the size in rigid case.

>initial porostiy
Sorry for not clearly describe it. Yes, the initial porosity I mentioned is the porosity at the begining of the deviatoric loading. In the following two MWE, they are both 0.358 getting by function of porosity(). But note that, in rigid boundary case, I also print the porosity at the same state by using triax.porosity, which retures 0.33.

>coordination number
Thanks Luc, that is a good point to mention. In the two MWE, I print the coordination number (CN) just before the deviatoric loading. And the CN of rigid boundary is 7.28, while the CN in periodic is 6.3. Which means there are more contacts in rigid wall case, indicating a denser state in general. Then my question is why they have same inputtings (porosity, particle number...) but they have different CN? And another question is how can I make the same sample under the two boundary conditions? It would be conflict if I control them as the same from porosity and from CN.

################################### rigidWall case ##########################
from __future__ import division
from yade import pack, plot
import numpy as np
num_spheres=7000# number of spheres
targetPorosity = 0.33 
compFricDegree = 30 
finalFricDegree = 35 
damp=0.6 
stabilityThreshold=0.001 
confinement=100e3
mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)  #determine the size of the sample
young=5e6
## create materials for spheres and plates
MatWall=O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(30),\
										 density=2650.0,label='sand'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])

Gl1_Sphere.quality=3
triax=TriaxialStressController(
	maxMultiplier=1.+2e4/young,
	finalMaxMultiplier=1.+2e3/young,
	thickness = 0,
	stressMask = 7,
	internalCompaction=True, 
)
newton=NewtonIntegrator(damping=damp)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
	newton,
	PyRunner(iterPeriod=100,command='history()',label='recorder'),
	PyRunner(command='stopIfDamaged()',iterPeriod=200,label="endSimulation"),
]

recorder.dead=True
endSimulation.dead=True
Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

def getCN():
	# get coordination number
	CN=0
	for i in O.bodies:
		if isinstance(i.shape,Sphere):
			CN +=len(i.intrs())
	Z = CN/ 7000.0
	return Z

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
	  break

print "###    state 1 Isotropic completed     ###"

import sys
while triax.porosity - targetPorosity>0.0001:
	compFricDegree = 0.95*compFricDegree
	setContactFriction(radians(compFricDegree))
	print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
	sys.stdout.flush()
	O.run(500,1)

print "###  state 2  Reach target porosity completed      ###"
print "particle radii is: ", O.bodies[2000].shape.radius ## here get the particle radius which will be used as input in makecloud in periodic boundary case.
print 'triax.porosity at 100 kPa',triax.porosity,'porosity()',porosity()
print ("coordination number at 100 kPa is",getCN())
topStress_initial=-triax.stress(triax.wall_top_id)[1]
bottomStress_initial=-triax.stress(triax.wall_bottom_id)[1]

# start deviatoric loading
triax.internalCompaction=False
triax.wall_bottom_activated=False
recorder.dead=False
endSimulation.dead=False
triax.stressMask=5
triax.goal2=-1
triax.goal1=-confinement
triax.goal3=-confinement
setContactFriction(radians(finalFricDegree))

def history():
	plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
			ev=triax.strain[0]+triax.strain[1]+triax.strain[2],
			s11=-triax.stress(triax.wall_right_id)[0],
			s22=-triax.stress(triax.wall_top_id)[1],
			s33=-triax.stress(triax.wall_front_id)[2],
			q=-triax.stress(triax.wall_top_id)[1]-100000,)

def stopIfDamaged():
	if -triax.strain[1]>0.4:
		O.pause()
		plot.saveDataTxt('results_rigidWall.txt')

plot.plots={'e22':('q',),'e22 ':('ev')}
plot.plot()
###################################### end of the 1st MWE #############

#################### periodic case #######################3
from __future__ import print_function
sigmaIso = -100e3
from yade import pack, qt, plot

O.periodic = True
O.cell.setBox(.09,.18,.09)
young=5e6
compFricDegree=23 # using this value can make the porosity at 100 kpa be 0.358, which is the same as that of rigid wall
finalFricDegree = 35
MatSand = O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(30),density=2650.0,label='sand'))
sp = pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.1,0.2,0.1)
sp.makeCloud(mn,mx,0.0025,0,7000,False, 0.95,seed=1)
sp.toSimulation(material="sand")
# here set contact friction angle to reach a low porosity at target confining pressure.
setContactFriction(radians(compFricDegree))
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        PeriTriaxController(
                label='triax',
                goal=(sigmaIso, sigmaIso, sigmaIso),
                stressMask=7,
                dynCell=True,
                maxStrainRate=(10, 10, 10),
                maxUnbalanced=.1,
                relStressTol=1e-3,
                doneHook='compactionFinished()'
        ),
        NewtonIntegrator(damping=.6),
        PyRunner(command='addPlotData()', iterPeriod=100,label="recorder"),
]
O.dt = .5 * PWaveTimeStep()
recorder.dead=True

def addPlotData():
	plot.addData(
	        unbalanced=unbalancedForce(),
	        i=O.iter,
	        sxx=-triax.stress[0],
	        syy=-triax.stress[1],q=-triax.stress[1]-100000,
	        szz=-triax.stress[2],
	        exx=-triax.strain[0],ev=triax.strain[0]+triax.strain[1]+triax.strain[2],
	        eyy=-triax.strain[1],
	        ezz=-triax.strain[2],

        )

O.trackEnergy = True

plot.plots = {
        'eyy': ('q',),
        'eyy  ': ('ev'),}
plot.plot()

def compactionFinished():
	print ("porosity at 100 kPa is",porosity())
	print ("coordination number at 100 kPa is",getCN())
	print ("O.cell.size is", O.cell.size)
	setContactFriction(radians(finalFricDegree))
	O.pause()
	O.cell.trsf = Matrix3.Identity
	triax.goal = (sigmaIso,-0.4, sigmaIso)
	triax.stressMask = 5
	triax.maxStrainRate = (1., 1., 1.)
	triax.doneHook = 'triaxFinished()'
	triax.maxUnbalanced = 10
	O.step() # if I didnot use O.step(), the initial eyy will not start from 0.
	recorder.dead=False
def triaxFinished():
	print('Finished')
	O.pause()
	plot.saveDataTxt('results_periodic.txt')
def getCN():
	CN=0
	for i in O.bodies:
		if isinstance(i.shape,Sphere):
			CN +=len(i.intrs())
	Z = CN/ 7000.0
	return Z
##################################################

Thanks
Leonard
[1]https://we.tl/t-xlkvIPB9zx
ps: both of the MWE are modified from YADE trunk.

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