wet phase distribution

Asked by Luis Barbosa on 2021-03-05

Hi all,

When using yade.wrapper.TwoPhaseFlowEngine is it possible to quantify the volume distribution of wet phase (not only the total volume) over the sample?

Cheers,
Luis

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
2021-03-16
Last query:
2021-03-16
Last reply:
2021-03-08
Luis Barbosa (luis-pires-b) said : #2

Thanks Robert.

I have just two more questions:

1. getCellIsWRes is giving if there is wet phase or not. Wehrever if it is trapped or not. Is that right?

2. I want to spatially locate the wet cells, so what is the difference between getCellCenter and getCell?

Kind regards,
Luis

Hi, documentation might help here:

In [4]: TwoPhaseFlowEngine.getCell?
Docstring:
getCell( (TwoPhaseFlowEngineT)arg1, (float)X, (float)Y, (float)Z) -> int :
    get id of the cell containing (X,Y,Z).
Type: function

In [5]: TwoPhaseFlowEngine.getCellCenter?
Docstring:
getCellCenter( (TwoPhaseFlowEngineT)arg1, (int)id) -> Vector3 :
    get voronoi center of cell 'id'.

They do nearly the opposite.

Bruno

Luis Barbosa (luis-pires-b) said : #4

Thank you all.
When I run the code [1], getCellBarycenter giving some strange values (coordinates totally outside my sample):

For instance, these values 125107.48 or -125097.48:
True 0.01797170859509369 Vector3(125107.4828372159245,1.62543863864436311,3.652322448466255977)
True 0.011656083836606964 Vector3(9.624066355278376506,1.530915196593882976,4.766563274628489388)
True 0.020528870575860336 Vector3(7.413446355544288657,-125097.4827497912484,7.885475000566313675)
True 0.01815466910920833 Vector3(125106.2605586981226,-125096.2810642897675,6.157849618659231083)

is it some mistake in the usage?

[1]
import matplotlib; matplotlib.rc('axes',grid=True)
from yade import pack
import pylab
from numpy import *

utils.readParamsFromTable(seed=1,compFricDegree = 15.0)
from yade.params import table

seed=table.seed
#num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process),num_spheres=1000
confiningS=-1e5

## creat a packing with a specific particle side distribution (PSD)
#psdSizes,psdCumm=[.599,0.6,0.849,0.85,1.0,1.40],[0.,0.35,0.35,0.70,.70,1.] psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,
sp=pack.SpherePack()
sp1=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(10,2,10)
mnn,mxx=Vector3(0,2,0),Vector3(10,10,10)
mnc,mxc=Vector3(0,0,0),Vector3(10,10,10)
sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=0.15,seed=seed)
sp1.makeCloud(minCorner=mnn,maxCorner=mxx,rMean=0.4,seed=seed)

## create material #0, which will be used as default
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=radians(compFricDegree),density=2600,label='spheres1'))
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=0,density=0,label='frictionless'))

## create walls around the packing
walls=aabbWalls((mnc,mxc),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
O.bodies.append([utils.sphere(center,rad,material='spheres1') for center,rad in sp1])

triax=TriaxialStressController(
 internalCompaction=True,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)

newton=NewtonIntegrator(damping=0.4)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
# VTKRecorder(iterPeriod=100,recorders=['all'],fileName="/home/user/Área de Trabalho/PVF/vtk/Spheres"),
 newton
]

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
    break

#############################
## REACH NEW EQU. STATE ###
#############################
finalFricDegree = 30 # contact friction during the deviatoric loading

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False
# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.001 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
    break

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
O.save('1kPacking.yade') #save the packing, which can be reloaded later.

O.run(1000,True)
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]

from yade import plot
O.engines=O.engines+[PyRunner(iterPeriod=20,command='history()',dead=1,label='recorder')]

def history():
 plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
  s11=-triax.stress(0)[0]-si0,
  s22=-triax.stress(2)[1]-si1,
  s33=-triax.stress(4)[2]-si2,
  pc=-unsat.bndCondValue[2],
  sw=unsat.getSaturation(False),
  i=O.iter
  )

plot.plots={'pc':('sw',None,'e22')}
plot.plot()

#######################################################
## Drainage Test under oedometer conditions ###
#######################################################
##oedometer conditions
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0
recorder.dead=0

##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.

##set boundary conditions, the drainage is controlled by decreasing W-phase pressure and keeping NW-phase pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
unsat.initialization()
unsat.surfaceTension = 10

#start invasion, the data of normalized pc-sw-strain will be written into pcSwStrain.txt
f1=open('Cells1.txt',"w")
f5=open('SwPc.txt',"w")
for pg in arange(1.0e-5,15.0,0.1):
  #increase gaz pressure at the top boundary
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  #compute the evolution of interfaces
  unsat.invasion()
  #save the phases distribution in vtk format, to be displayed by paraview
# unsat.savePhaseVtk("/home/user/Área de Trabalho/PVF/vtk/Pressurefield")
  #compute and apply the capillary forces on each particle
  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))

  if pg==0.60001:
    cels=unsat.nCells()
    celsW1 = [0.0]*cels
    celsV1 = [0.0]*cels
    celsBar1 = [0.0]*cels
    for ii in range(cels):
      celsW1=unsat.getCellIsWRes(ii)
      celsV1=unsat.getCellVolume(ii)
      celsBar1=unsat.getCellBarycenter(ii)
      f1.write(str(celsW1)+" "+str(celsV1)+" "+str(celsBar1)+"\n")
    f1.close()

  #reac
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    if unb<0.001:
      break
  f5.write(str(pg)+" "+str(unsat.getSaturation(False))+" "+str(triax.strain[1]-ei1)+"\n")
f5.close()

Best Robert Caulk (rcaulk) said : #5

>When I run the code [1], getCellBarycenter giving some strange values (coordinates totally outside my sample):

Well, nCells includes all cells in the triangulation, including cells associated with the boundary. The boundary would by definition not be inside your specimen. You can filter those out using spatial coordinates if you only want cells "inside your specimen."

Luis Barbosa (luis-pires-b) said : #6

Thank!

Luis Barbosa (luis-pires-b) said : #7

Thanks Robert Caulk, that solved my question.