New question #662105 on Yade:
https://answers.launchpad.net/yade/+question/662105
Hello there,
I have been struggling with this for a while now and could not figure out what
is wrong. My simulation (shear test with periodic boundary conditions) explodes
as soon as the flow engine is activated (see the last 2 lines) and I don't know
if it is due to the technical problems that I am facing with the flow engine or
if there is something wrong with my script. Could someone please try it and see
if it works on his side?
Many thanks
Luc
###
from yade import pack
sp=pack.SpherePack()
O.periodic=True
# dimensions of sample
RADIUS=0.05
length=10*(2*RADIUS)
height=length
width=length
thickness=RADIUS/10.
# microproperties
DENS=2600
E=1e8
P=0.5
compFRIC=10.
FRIC=30.
TENS=0.
COH=0.
# loading conditions
SN=5e6
RATE=0.1
####
O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)
O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=DENS,young=E,poisson=P,frictionAngle=radians(compFRIC),normalCohesion=TENS,shearCohesion=COH,label='sphereMat'))
upBox =
utils.box(center=(0,2*height+thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox =
utils.box(center=(0,height-thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])
sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.2,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for
s in sp])
effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol
#### engines
flow=PeriodicFlowEngine(
isActivated=0
,useSolver=3
### what should we define in the next 2 lines?
,defTolerance=-1 # with this value, triangulation is done according to
meshUpdateInterval
,meshUpdateInterval=1000 # this value by default
,duplicateThreshold=0.5 # how do we define this?
,boundaryUseMaxMin=[0,0,0,0,0,0] # what is this for? should we have
[0,0,1,1,0,0] because of the top and bottom walls?
,wallIds=[-1,-1,1,0,-1,-1] # top wall id=0, bottom wall id=0
,wallThickness=thickness
,bndCondIsPressure=[1,1,0,0,1,1] # drained in the horizontal plane
,bndCondValue=[0,0,0,0,0,0]
,permeabilityFactor=1.e-3
,viscosity=1.
,label='flowEng'
)
O.engines=[
ForceResetter()
,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
,InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
)
,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-1.e5/volRatio,-1.e5/volRatio,-1.e5/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='compactionDone()',label='triax')
,flow
,NewtonIntegrator(damping=0.2,label='newton')
]
def compactionDone():
print 'Changing contact properties now'
tt=TriaxialCompressionEngine()
tt.setContactProperties(FRIC)
print 'DONE! iter=',O.iter,'| isotropic confinement done:
stresses=',volRatio*triax.stress,
triax.dead=True
O.pause()
#### Initialization
print 'SAMPLE PREPARATION!'
O.run(1000000,1)
#### Applying normal stress
print 'NORMAL LOADING! iter=',O.iter
stage=0
stiff=fnPlaten=currentSN=0.
def servo():
global stage,stiff,fnPlaten,currentSN
if stage==0:
currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
unbF=unbalancedForce()
boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
O.bodies[0].state.vel[1]=boundaryVel
if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
stage+=1
fnPlaten=O.forces.f(0)[1]
print 'Normal stress =',currentSN,' | unbF=',unbF
for i in O.interactions.withBody(O.bodies[0].id):
stiff+=i.phys.kn
print 'DONE! iter=',O.iter
O.pause()
if stage==1:
fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
boundaryVel=copysign(min(RATE,abs(0.35*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
# why 0.35?
O.bodies[0].state.vel[1]=boundaryVel
O.engines =
O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]
O.run(1000000,1)
print 'Normal stress (platen) =
',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) =
',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]
#### preparing for shearing
print 'Gluing spheres to boundary walls'
## gluing particles in contact with the walls
for i in O.interactions:
if i.isReal:
if isinstance(O.bodies[i.id1].shape,Box):
O.bodies[i.id2].shape.color[0]=1
if isinstance(O.bodies[i.id2].shape,Box):
O.bodies[i.id1].shape.color[0]=1
for i in O.interactions:
if i.isReal and ( O.bodies[i.id1].shape.color[0]==1 and
O.bodies[i.id2].shape.color[0]==1 ):
O.bodies[i.id1].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
O.bodies[i.id2].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
O.bodies[i.id1].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
O.bodies[i.id2].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
i.phys.initCohesion=True
#### shearing
print 'SHEARING! iter=',O.iter
flow.isActivated=1 # sample explodes as soon as fluid is activated???
O.bodies[0].state.vel[0]=RATE
O.run(10)
--
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 : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
