yade-users team mailing list archive

[Huang peilun <question692704@xxxxxxxxxxxxxxxxxxxxx>]

Question #692704 on Yade changed:
https://answers.launchpad.net/yade/+question/692704

Huang peilun posted a new comment:
Hi Robert,

Thank you for answering my question. I'm sorry that I didn't read the
posting guidelines carefully. I'll follow the guidelines in the future.

As for the issue, I tried decreasing the time step and turning 
O.dynDt=False.
The result is almost the same. I noticed that the pressure field can be calculated by the following formula:

GP=Ev+Qq+Qp

where G is the conductivity matrix, P is the column vector containing
all values of pressure, and E is the matrix defining the rates of volume
change. Qq and Qp are flux vectors reflecting boundary conditions,
respectively source terms (imposed fluxes) in Qq and imposed pressures
in Qp.

I checked the initial velocity of my sample, it appears that all
particles have a high velocity, thus in the second iteration gives a
'unreasonable' pressure field. I tried to slow down the particles, and
the result is more and more 'reasonable'. Therefore, I think the way I
generate the sample is improper. The code that I use to generate the
sample is as followed:

####################################################
####################################################

from yade import pack

SoilMat=CohFrictMat(young=1e7,poisson=0.3,density=2650,frictionAngle=0.7,alphaKr=50,alphaKtw=50,momentRotationLaw=True)
O.materials.append((SoilMat))
O.materials.append(CohFrictMat(young=1e7,poisson=0.5,frictionAngle=0,density=0,label='walls',momentRotationLaw=True))

O.bodies.append(geom.facetBox((0,0,0.8),(.25,.25,0.8),wallMask=31,material='walls'))

sp=pack.SpherePack()
sp.makeCloud((-.25,-.25,0),(.25,.25,1.6),rMean=.0175,rRelFuzz=.3)
sp.toSimulation(material=SoilMat)

#Engine of the simulation
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
		[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
		[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
	),
	#introduced as a dead engine for the moment
	NewtonIntegrator(gravity=(0,0,-9.81)),
	#Run the specific controlled motion.
	PyRunner(command='start()',realPeriod=1,label='checker'),
]

O.dt=.5*PWaveTimeStep()

def start():
	newton.damping=0.3
	print('\n##### Simulation begins #####')
	print('Stage: Gravitity desposition.')
	print('Time step is %ss.' % O.dt)
	checker.command='checkUnbalanced()'

def checkUnbalanced():	
	print('%s:UnbalancedForce=%s' % (O.iter,unbalancedForce()))
	if max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])>0.7: return
	if unbalancedForce()>0.01: return 
	Initial_height=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
	print('The initial height of the simulation is %s' % Initial_height)	
	print('Stage: Compression.')
	checker.command='Compression()'

def Compression():
	O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
	global plate
	plate=O.bodies[-1]
	plate.state.vel=(0,0,-.2)
	print('Wall successfully created! Initial veloctity is %s, initial position is %s.' % (plate.state.vel[2],plate.state.pos[2]))
	global Volume	
	Volume=0
	for b in O.bodies:
		if isinstance(b.shape,Sphere):
			Volume=Volume+b.shape.radius*b.shape.radius*b.shape.radius*4*np.pi/3
	print('The volume of the pack is %s' % Volume)
	SoilMat.frictionAngle=0
	checker.command='Compression_S2()'

def Compression_S2():
	void_ratio=plate.state.pos[2]*.5*.5/Volume-1
	print('%s: Void ratio=%s' % (O.iter,void_ratio))
	if void_ratio>0.6: return
	SoilMat.frictionAngle=0.6
	plate.state.vel=(0,0,.2)
	print('Wait for the particles to settle down after compression.')
	checker.command='Density_Checker()'

def Density_Checker():
	void_ratio=plate.state.pos[2]*.5*.5/Volume-1
	print('%s: Void ratio=%s, Wall Position=%s' % (O.iter,void_ratio,plate.state.pos[2]))
	if abs(O.forces.f(plate.id)[2])>0: return
	if unbalancedForce()>0.01: return 
	print(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]))
	checker.command='Record()'

def Record():
	print('The sample is generated successfully!')
	O.pause()

####################################################
####################################################

I think the way I generate the sample is not stable. And also, to spend
less time, the material is different between when I generate the sample
and when I do the seepage test.

I'm now trying to use the radius expansion method to generate the sample and use the same material in both test.
I hope this time it will work so that I can post my result as the end of this question.

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