yade-users team mailing list archive

[Zoheir Khademian <question698778@xxxxxxxxxxxxxxxxxxxxx>]

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

I am trying to define a cavity at the center of a sphere pack. I remove some spheres at the center before making the cells at the center part of the cavity [1]. Then, I assign a cavity flux using [2]. However, when I run the model, I dont get any outflow. I monitor outflow both at the model boundary and cavity boundaries [3]. 

[1] https://yade-dev.gitlab.io/trunk/yade.wrapper.html#yade.wrapper.FlowEngine.imposeCavity 
[2] https://yade-dev.gitlab.io/trunk/yade.wrapper.html#yade.wrapper.FlowEngineT.cavityFlux
[3] https://yade-dev.gitlab.io/trunk/yade.wrapper.html#yade.wrapper.FlowEngineT.getCavityFlux


Here is the MWE: 


from yade import pack, ymport, plot, utils, export, timing
import numpy as np

young=5e6

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)


O.materials.append(FrictMat(young=young*100,poisson=0.5,frictionAngle=0,density=2600e10,label='walls'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600e10,label='spheres'))


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

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0015,rRelFuzz=0.333,num=200,seed=11) 
sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
 FlowEngine(dead=1,label="flow",multithread=False),
 NewtonIntegrator(damping=0.5)
]

O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1

tri_pressure = 1000
triax.goal1=triax.goal2=triax.goal3=-tri_pressure
triax.stressMask=7
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<0.1 and abs(-tri_pressure-triax.meanStress)/tri_pressure<0.001:
    break

triax.internalCompaction=False
minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
dz=maxZ-minZ
dy=maxY-minY
dx=maxX-minX
print ("minX:",minX,"maxX:",maxX,"minY:",minY,"maxY:",maxY,"minZ:",minZ,"maxZ:",maxZ)
CavityList=[]
for b in O.bodies:
	if isinstance(b.shape,Sphere):
		if np.linalg.norm(b.state.pos-(maxX-dx/2.,maxY-dy/2.,maxZ-dz/2.))<dz/4.: 
			CavityList.append(b.state.pos)
			O.bodies.erase(b.id)

for b in O.bodies:
	if isinstance(b.shape, Sphere): 
		b.dynamic=False # mechanically static
flow.dead=0

flow.debug=False
# add flow
flow.permeabilityMap = False
flow.pZero = 0
flow.meshUpdateInterval=5
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=-1e-5
flow.viscosity= 0.001
flow.decoupleForces =  False
flow.bndCondIsPressure=[0,0,0,0,0,1]
flow.bndCondValue=[0,0,0,0,0,0]


# defining cavity
for b in CavityList:
	ini=flow.imposeCavity(b)
	#print("ini",ini)
flow.controlCavityPressure=1
flow.cavityFlux=-1e-5
#flow.isCavity=1

flow.emulateAction()
O.dt=1e-8
O.dynDt=False


from yade import plot

def history():
 p1=flow.getBoundaryFlux(flow.wallIds[flow.xmin])
 p2=flow.getBoundaryFlux(flow.wallIds[flow.xmax])
 p3=flow.getBoundaryFlux(flow.wallIds[flow.ymin])
 p4=flow.getBoundaryFlux(flow.wallIds[flow.ymax])
 p5=flow.getBoundaryFlux(flow.wallIds[flow.zmin])
 p6=flow.getBoundaryFlux(flow.wallIds[flow.zmax])
 Outflow=(p1+p2+p3+p4+p5+p6)##
 CavityFlow=flow.getCavityFlux
 print('Outflow',Outflow, 'CavityFlow',CavityFlow)

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]

O.run(1000)

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