how to get the change of permeability during triaxial compression in flowengine

Asked by Ziyu Wang

Hi,
My purpose is to study the change law of permeability of cylinder sample during triaxial compression, which largely refers to some methods in [1].But different from [1], I want to monitor the real-time change of permeability, that is, add it to the plot.addData().

Referring to the method in [1], I defined a function as the record of permeability. However, this leads to expensive calculation costs (in other words, because the permeability is calculated every time, the simulation runs very slowly and can hardly run)

Is there any way to solve this problem, or is there a better way to achieve my goal of recording permeability changes?

Thanks!

############MWE#############
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys
damp=0.4
young=66.2e9
name='JCFPM_triax'
poisson=0.522
name='cylinder'
preStress=-30e6
strainRate = -0.1
OUT=str(preStress)+'_JCFPM_triax'
r1=0.002
r2=0.01
nw=24
nh=15
rParticle=0.000731723
bcCoeff = 5
width = 0.025
height = 0.05

allx,ally,allz,r=np.genfromtxt('packing-cylinder.spheres', unpack=True)
mnx=min(allx)*1.01 #same as aabbExtrema, get the xyz of the lowest corner, multiply by factor to reduce it to let the wall surround all the spheres.
mny=min(ally)*1.01
mnz=min(allz)*0.99
mxx=max(allx)*1.01 #same as aabbExtrema, get the xyz of the top corner, multiply by factor to increase it to let the wall surround all the spheres.
mxy=max(ally)*1.01
mxz=max(allz)*1.01
Sy=2e6

O.materials.append(JCFpmMat(type=1,density=2640,young=young,poisson=poisson,tensileStrength=8e6,cohesion=40e6,frictionAngle=radians(25),label='sphere'))
O.materials.append(JCFpmMat(type=1,label='wall1',young=young,poisson=poisson,frictionAngle=radians(25)))
O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='wall2'))

mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0,material='wall2')
wallIds=O.bodies.append(walls)

bottom,top=Vector3(0,0,0),Vector3(0,0,0.05)
sp=pack.SpherePack()
pred=pack.inCylinder(bottom,top,radius)
sp=pack.randomDensePack(pred,radius=0.0005,rRelFuzz=0,returnSpherePack=True,memoizeDb='/tmp/triax.sqlite')
spheres=sp.toSimulation(color=(0.9,0.8,0.6),material='sphere')

bot = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]<rParticle*bcCoeff]
top = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
vel = strainRate*(height-rParticle*2*bcCoeff)
for s in bot:
 s.shape.color = (1,0,0)
 s.state.blockedDOFs = 'xyzXYZ'
 #s.state.vel = (0,0,-vel)
for s in top:
 s.shape.color = Vector3(0,1,0)
 s.state.blockedDOFs = 'xyzXYZ'
 s.state.vel = (0,0,vel)

facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
 for h in range(nh):
  v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
  v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
  v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
  v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
  f1 = facet((v1,v2,v3),color=(0,0,1),material='wall1')
  f2 = facet((v1,v3,v4),color=(0,0,1),material='wall1')
  facets.extend((f1,f2))

O.bodies.append(facets)
mass = O.bodies[7].state.mass
for f in facets:
 f.state.mass = mass
 f.state.blockedDOFs = 'XYZz'

def addForces():
 for f in facets:
  n = f.shape.normal
  a = f.shape.area
  O.forces.addF(f.id,preStress*a*n)
O.dt=.5*utils.PWaveTimeStep()
newton=NewtonIntegrator(damping=damp)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ss2sc'),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT+'_Crack',label='interactionLaw')]
 ),
 PyRunner(iterPeriod=1,command="addForces()"),
 FlowEngine(dead=1,label="flow"),
 newton,
 PyRunner(command='plotAddData()',iterPeriod=10),
]

def plotAddData():
 perme = permeability()
 f1 = sum(O.forces.f(b.id)[2] for b in top)
 f2 = sum(O.forces.f(b.id)[2] for b in bot)
 f = .5*(f2-f1)
 s = f/(pi*.25*width*width)
 e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / (height-rParticle*2*bcCoeff)
 plot.addData(
  i = O.iter,
  s = -s,
  e = -e,
  tc = interactionLaw.nbTensCracks,
  sc = interactionLaw.nbShearCracks,
  permeab = perme
 )
 plot.saveDataTxt(OUT)
O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1

flow.debug=False
flow.permeabilityMap = False
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=1
flow.viscosity= 0.001
flow.decoupleForces = False
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,Sy]
flow.dead=0
flow.emulateAction()

#################blockcells################
numPoints = 100
xs = np.linspace(1.05*mnx,1.05*mxx,numPoints) #np.linspace is to divide the distance between 2 points
ys = np.linspace(1.05*mny,1.05*mxy,numPoints) #if the factor for mx change to less than 1, code will not work properly.
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)

cellsHit = [] #create array
cylcenterx=(mxx+mnx)/2
cylcentery=(mxy+mny)/2
for x,y,z in itertools.product(xs, ys, zs):
 cellId = flow.getCell(x,y,z) #getCell return cell id for a point position xyz
 if cellId in cellsHit: continue #this prevents counting a cell multiple times
 if np.sqrt((x-cylcenterx)**2+(y-cylcentery)**2) > width:
  cellsHit.append(cellId)

for i in cellsHit:
 flow.blockCell(i,blockPressure=True)
 flow.setCellPressure(i,0)
O.run(1,1)

flow.meshUpdateInterval=-1 #these two lines to prevent remeshing after 1000 run which unblock all cells in cellsHit
flow.defTolerance=-1

numPoints1 = 5
xr = np.linspace(1.05*mnx,1.05*mxx,numPoints) #np.linspace is to divide the distance between 2 points
yr = np.linspace(1.05*mny,1.05*mxy,numPoints)
zr = np.linspace(0.95*mnz,0.95*mxz,numPoints)
width1=width*0.98
def permeability():
 totalVolume = 0
 v = np.array([0,0,0])
 cellsHit2=[]
 for x,y,z in itertools.product(xr,yr,zr):
  cellId2=flow.getCell(x,y,z)
  if cellId2 in cellsHit2:
   continue
  if np.sqrt((x-cylcenterx)**2+(y-cylcentery)**2) < width1:
   cellsHit2.append(cellId2)
   velocityVector = np.array(flow.getCellVelocity((x,y,z)))
   velMag = np.linalg.norm(velocityVector)
   cellVol = flow.getCellVolume((x,y,z,))
   v = v + cellVol*velocityVector
   totalVolume += cellVol
 Cylinder_Vel = np.linalg.norm(v)/totalVolume
 return Cylinder_Vel

plot.plots = { 'e':('s',), 'elateral':('s'),}
plot.plot()
O.run()
################################
[1]https://answers.launchpad.net/yade/+question/689330

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

Hello,

The method you are using here is inherently expensive. Why not simply use a gradient across the specimen and get the flow.getBoundaryFlux() function?

You should be putting your custom permeability() function into a PyRunner(). Inside the PyRunner, you indicate the iterPeriod, which lets you decide how often that function is executed. I hope it is obvious that this is one way to reduce computational effort.

Cheers,

Robert

Revision history for this message
Ziyu Wang (ziyuwang1) said :
#2

Thanks Robert for your suggestion.
>Why not simply use a gradient across the specimen and get the flow.getBoundaryFlux() function?

I don't know the meaning of ‘gradient’ very well.If I set the flow condition as follows, the value of flow.getboundaryflux() I get seems to remain unchanged..
###
flow.permeabilityFactor=1
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,Sy]
###
(About cylinder,is boundaryUseMaxMin=[0,0,0,0,0,0] suitable?)

>You should be putting your custom permeability() function into a PyRunner().

This is a good idea. I can reduce the calculation time by increasing iterperiod, but if I use this method, how should I record the calculated permeabilities?(such as how to add it into plot.addData()?)

Cheers,

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

Calling a function like this in your addData:

def getPermeability():

 # set new bcs for perm check
 flow.bndCondIsPressure = [0, 0, 0, 0, 1, 1]
 flow.bndCondValue = [0, 0, 0, 0, deltaP, 0]

 # force update the triangulation
 flow.updateTriangulation = True

 flow.emulateAction()
 flow.emulateAction()

 Qin = flow.getBoundaryFlux(4)
 Qout = flow.getBoundaryFlux(5)

 perm = abs(Qin) * flow.viscosity * L / (A * deltaP)

 print("Q", Qin, 'perm', perm, "A", A, "L",L)

        # reset boundary conditions before returning to sim
 flow.bndCondIsPressure = [0, 0, 0, 1, 0, 0]
 flow.bndCondValue = [0, 0, 0, 0, 0, 0]
 flow.updateTriangulation = True

 return perm

Revision history for this message
Ziyu Wang (ziyuwang1) said :
#4

Hi,Robert.

Thanks for the function you provide,but I don't particularly understand this method.

In my understanding, the flow pressure (flow.bndcondvalue) should be continuously applied to the sample, but the method you provide seems to be applied only during measurement. Is there a mistake in my understanding of the method of flowengine?

Cheers,

Revision history for this message
Ziyu Wang (ziyuwang1) said :
#5

Hi,
I have known how to measure the permeabilty about my question,but still not understand the bndcondvalue you provided..I will open an new thread asking about the bndcondValue.

Thanks again Robert!

Revision history for this message
Ziyu Wang (ziyuwang1) said :
#6

Thanks Robert Caulk, that solved my question.