yade-users team mailing list archive

[Luc Scholtès <question662105@xxxxxxxxxxxxxxxxxxxxx>]

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.