yade-users team mailing list archive

[ytang <question693734@xxxxxxxxxxxxxxxxxxxxx>]

New question #693734 on Yade:
https://answers.launchpad.net/yade/+question/693734

Hi All,

I'm trying to use this constitutive model (Law2_ScGeom_MindlinPhys_Mindlin()) to simulate the triaxial tests.
here is my code:
###################################
############################
### DEFINING PARAMETERS  ###
############################
from yade import pack, plot
targetPorosity=0.3682
compFricDegree=45.0
finalFricDegree=19.5
rate=-0.5
damp=0.7
stabilityThreshold=0.001
young=4e8
mn,mx=Vector3(0,0,0),Vector3(4e-3,8e-3,4e-3)
##################################################################
###################### material for sphere and wall ##############
##################################################################
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2648,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
psdSizes=[0.12e-3,0.14e-3,0.16e-3,0.18e-3,0.19e-3,0.20e-3,0.22e-3,0.24e-3,0.29e-3,0.38e-3]
psdCumm=[0,11.1,22.2,33.3,44.4,55.5,66.6,77.7,88.8,100]
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=11526,distributeMass=1,seed=1)
sp.toSimulation(material='spheres')
O.dt=.5*PWaveTimeStep()
############################
###   DEFINING ENGINES   ###
############################
triax=TriaxialStressController(
	maxMultiplier=1.0+2e5/young,
	finalMaxMultiplier=1.0+2e4/young,
	thickness=0,
	stressMask=7,
	internalCompaction=False
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_FrictMat_FrictMat_MindlinPhys( betan =0.2 , betas = 0.2, krot = 0.2 , eta = 0.5  )],
			[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom_MindlinPhys_Mindlin()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
	triax,
	#VTKRecorder(fileName='george-triaxial-',recorders=['all'],iterPeriod=1000),
	#TriaxialStateRecorder(iterPeriod=2000,file='WallStresses-george'+table.key),
	newton
]
########################################
####   APPLYING CONFINING PRESSURE   ###
########################################
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
	O.run(3000,True)
	unb=unbalancedForce()
	print 'unbalanced force: ',unb,'mean stress: ',triax.meanStress
	if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.001:
		break
####################################################
####   REACHING A SPECIFIED POROSITY PRECISELY   ###
####################################################
##import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
	compFricDegree = 0.95*compFricDegree
	setContactFriction(radians(compFricDegree))
	print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
	sys.stdout.flush()
	O.run(500,1)
#set new friction angle
triax.wall_bottom_activated=True
triax.wall_top_activated=True
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True
setContactFriction(radians(finalFricDegree))
while 1:
	O.run(3000,True)
	unb=unbalancedForce()
	print 'unbalanced force: ',unb,'mean stress: ',triax.meanStress
	if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.001:
		break
print "porosity= ",triax.porosity
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
    O.run(1000, True)
    #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
    unb=unbalancedForce()
    print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
    if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.001:
        break
#O.save('compactedState'+key+'.yade.gz')
print "###    Compacted state saved      ###"
print "porosity= ",triax.porosity
re11=-triax.strain[0]
re22=-triax.strain[1]
re33=-triax.strain[2]
################################
#####   DEVIATORIC LOADING   ###
################################
triax.stressMask = 5
newton=NewtonIntegrator(damping=0.0)
#now goal2 is the target strain rate
triax.goal2=rate
# we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-100000
triax.goal3=-100000
triax.strainRate2=1
triax.strainRate1=triax.strainRate3=1000.0
######################################################
####    Example of how to record and plot data     ###
######################################################
def history():
  	plot.addData(e11=-triax.strain[0]-re11, e22=-triax.strain[1]-re22, e33=-triax.strain[2]-re33,
  		    ev=-triax.strain[0]-re11-triax.strain[1]-re22-triax.strain[2]-re33,
		    s11=-triax.stress(triax.wall_right_id)[0],
		    s22=-triax.stress(triax.wall_top_id)[1],
		    s33=-triax.stress(triax.wall_front_id)[2],
	            i=O.iter)  
	plot.saveDataTxt('calibration-george.txt')
if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1000,command='history()',label='recorder')]+O.engines[5:7]
  #O.engines=O.engines[0:5]+[PyRunner(iterPeriod=2000,command='savedata()',label='recorder_1')]+O.engines[5:7]
  ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
  O.engines[4]=PyRunner(iterPeriod=1000,command='history()',label='recorder')
  #O.engines[4]=PyRunner(iterPeriod=2000,command='savedata()',label='recorder_1')
O.run(1000000,True)
####################################################
I'm using High-performance computing with 22 CPUs.
########################################
There are no errors.
#######################
The total number is around 10000. 
The unbalanced force did not reach the threshold (= 0.001) after one day.
I want to know if this situation is reasonable? (because when I use this constitutive model: Law2_ScGeom6D_CohFrictPhys_CohesionMoment,  the unbalanced force reaches the threshold just 7 hours ).
#############################

best,
Yong








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