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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
