Problem with imposing flux in a cavity in flow engine

Asked by Zoheir Khademian

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)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said :
#1

Hey Zoheir,

Thanks for the good question.

So when I was creating this functionality, I came to the conclusion that it was actually more effective to keep the bodies in the cavity, but "unbound" them so that they do not mechanically interact with the rest of the specimen. This allows the triangulation to more properly triangulate the cavity space, so we don't end up with possibly overlapping cell centers or cell centers outside of cells.

So to do this I would recommend *not* deleting any bodies. Instead do something like:

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)
   b.bounded = False # this body will not find interactions with other bodies. Careful if the specimen itself needs displacement/rotation as these unbounded bodies will not displace/rotate with it. If it is the case, we need to assign the same motion to these bodies as we do to the rigid specimen.

and then flow.imposeCavity is meant for the cells, not the bodies (cells are the ones participating in the flow problem. Better to iterate through your cell list and find the cells that match your spatial criterion:

for i in range(0,flow.nCells()):
  coords = flow.getCellCenter(i)
  if np.linalg.norm(coords[0] - maxX-dx/2., coords[1] - maxY-dy/2., coords[2] - maxZ-dz/2.))<dz/4.:
    flow.imposeCavity((coords[0], coords[1], coords[2]))

then I would simply assign the total flux to one of your cavity cells:

flow.imposeFlux((maxX-dx/2,maxY-dy/2.,maxZ-dz/2.),-1e-5)

Make sure to keep flow.controlCavityPressure=False unless you want all the cells to have an equivalent imposed pressure.

You can then use flow.getCavityFlux() to get the total flux across the boundary of this assigned cavity.

Please let me know if you are able to get the results you expect. If you have any issues don't hesitate to continue the thread - this is a good feature to ensure we have working and documented.

Thanks for your help,

Robert

Revision history for this message
Robert Caulk (rcaulk) said :
#2

By the way, if you are using ThermalEngine, it should unbound cavity bodies automatically. So in this case you would only need to define the cells of the cavity and let yade take care of the rest. No need to loop on bodies and unbound those.

Revision history for this message
Zoheir Khademian (zkhademian) said :
#3

Thanks Robert for the detailed information. I followed your comments and made changes to the MWE. I injected 1e-4 m3/s into the cavity but the output of flow.getCavityFlux() [1] is about 8e-18 m3/s. Would you please let me know what I am missing? By the way, does flow.imposeFlux apply volumetric flux or mass flux?

[1] https://yade-dev.gitlab.io/trunk/yade.wrapper.html#yade.wrapper.FlowEngineT.getCavityFlux

Here is the updated 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)
   b.bounded = False #
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=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,0]

flow.emulateAction()
# defining cavity
for i in range(0,flow.nCells()):
  coords = flow.getCellCenter(i)
  if np.linalg.norm((coords-(maxX-dx/2.,maxY-dy/2.,maxZ-dz/2.)))<dz/4.:
    flow.imposeCavity((coords[0], coords[1], coords[2]))

flow.controlCavityPressure=0
O.dt=1e-8
O.dynDt=False
flow.imposeFlux((maxX-dx/2,maxY-dy/2.,maxZ-dz/2.),-1e-4)

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)

Revision history for this message
Best Robert Caulk (rcaulk) said :
#4

Hey Zoheir,

>> I injected 1e-4 m3/s into the cavity but the output of flow.getCavityFlux() [1] is about 8e-18 m3/s. Would you please let me know what I am missing?

Before I dive into the script myself (I won't have time until the weekend), I will propose some possible reasons. You are using a compressible fluid, and quite a small time step with very few time steps before checking the cavity flux. So technically, the governing equations would allow for an unequal fluid flux exiting the cavity to that which is entering the cavity. However that contrast in fluxes seems quite stark and so it is very possible something else is wrong.

Have you saved the VTK output to verify that the cavity is defined as you expect it to be? This will also allow you to look at the pressure gradients and verify that they are also as you expect them to be.

You can also compute the flux yourself by iterating on the cavity cells. This would verify whether or not the function "getCavityFlux()" is correct. (btw it is useful to store the list of cavityCells by their ID so you can just iterate on those directly eachtime you want to make a calculation such as this).

totalVolume = 0
v = np.array([0,0,0])
cavityCells = []
for i in range(0,flow.nCells()):
 coords = flow.getCellCenter(i)
 if flow.getCellPImposed(i): continue
 if np.linalg.norm((coords-(maxX-dx/2.,maxY-dy/2.,maxZ-dz/2.)))<dz/4.:
  cavityCells.append(i)
  velocityVector = np.array(flow.getCellVelocity((coords[0],coords[1],coords[2])))
  cellVol = flow.getCellVolume((coords[0],coords[1],coords[2]))
  v = v + cellVol*velocityVector
  totalVolume += cellVol

q = np.linalg.norm(v)/totalVolume

>> By the way, does flow.imposeFlux apply volumetric flux or mass flux?
In FlowEngine, it is a volumetric flux.

Cheers,

Robert

Revision history for this message
Robert Caulk (rcaulk) said :
#5

Also, I notice that you define the cavity cells *after* emulating action in flow and setting flow.dead=0. Please do this *before* to avoid any kind of undefined behavior in the solver.

Revision history for this message
Zoheir Khademian (zkhademian) said :
#6

Thanks so much Robert. After checking the VTK output for pressure, I noticed the cell I was injecting to was not registered as a cavity. So, simply applying the flux to one of the cavity cells in the cavity list solved the problem. Now, the influx and outflux are the same.

One more question though:
In your flux calculation above, the unit of q the flux seams to be m/s. Do we need to have the area of the cavity boundary to obtain the flux in m3/s? if so, what is the easiest way to get that area?

Thanks again for the help!

Revision history for this message
Zoheir Khademian (zkhademian) said :
#7
Revision history for this message
Robert Caulk (rcaulk) said :
#8

Hey Zoheir,

Glad you figured out the problem.

Re q: I believe the units are m^3/s. We are

v = v + cellVol*velocityVector # multiplying the velocity vector (m/s) by the volume (m^3/s), m^4/s
totalVolume += cellVol # summing volumes of all cavity cells to get total cavity volume (m^3)

q = np.linalg.norm(v)/totalVolume # m^4/s divided by m^3 -> m^3/s

Robert