yade-users team mailing list archive

[Leonard <question690248@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi,
I am wondering that how to reach target confining pressure in a triaxial simulation.

I followed Bruno's example[1], using radius expansion method. The target confining pressure is 100kPa.

After "Isotropic state saved" completed, the mean stress is OK but the stress of the top wall could be very different from 100kPa sometimes. In addition, after "REACHING A SPECIFIED POROSITY PRECISELY" completed, the stresses changed again.

Do you have any suggestions to reach target confining pressure?

Another question is that when I run the same script twice with the same parameters, I got different results with the stress of the top wall. Sometimes I got a reasonable stress of the top wall (i.e., the stress is close to target stress 100kPa), but when I run it again, I got different results. How to solve this problem? You may run the following MWE twice to see the difference if you are interested, it costs around 2 minutes.
####
from __future__ import division
from yade import pack, plot
import numpy as np

num_spheres=7000
targetPorosity = 0.44
compFricDegree = 25 
finalFricDegree = 18
rate=-0.01
damp=0.4
stabilityThreshold=0.001
young=1e8
eta=0.14
confinement=100e3
KrKt=7.0


mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

MatWall=O.materials.append(FrictMat(young=young*10,poisson=0.3,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(CohFrictMat(isCohesive=True,young=young,alphaKr=KrKt,alphaKtw=KrKt,\
                                         poisson=0.3,frictionAngle=radians(30),etaRoll=eta,etaTwist=eta,\
                                         density=2650.0,normalCohesion=1e6, shearCohesion=1e6,\
                                         momentRotationLaw=True,label='sand'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

sp.makeCloud(mn,mx,-1,0.3,num_spheres,False, 0.95,seed=1)

O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])

Gl1_Sphere.quality=3

###   DEFINING ENGINES   ###

triax=TriaxialStressController(
    maxMultiplier=1.+2e7/young, 
    finalMaxMultiplier=1.+2e5/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_CohFrictMat_CohFrictMat_CohFrictPhys(),
		 Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
	newton,
]

Gl1_Sphere.stripes=0

triax.goal1=triax.goal2=triax.goal3=-confinement

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  state completed     ###"

#REACHING A SPECIFIED POROSITY PRECISELY   ###

import sys
while triax.porosity -targetPorosity > 0.000001:
	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 'top', -triax.stress(triax.wall_top_id)[1], 'right',\
-triax.stress(triax.wall_right_id)[0],'front',-triax.stress(triax.wall_front_id)[2]

topStress_initial=-triax.stress(triax.wall_top_id)[1]

f = open("outputTopStress.txt", "a+")
f.write('porosity{}_topWallStress{} \n' .format(triax.porosity,topStress_initial))
f.close()
###

Thanks,
Leonard
[1]https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py



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