how to get the change of permeability during triaxial compression in flowengine
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!
#######
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(
r1=0.002
r2=0.01
nw=24
nh=15
rParticle=
bcCoeff = 5
width = 0.025
height = 0.05
allx,ally,
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.
O.materials.
O.materials.
mn,mx=Vector3(
walls=aabbWalls
wallIds=
bottom,
sp=pack.
pred=pack.
sp=pack.
spheres=
bot = [O.bodies[s] for s in spheres if O.bodies[
top = [O.bodies[s] for s in spheres if O.bodies[
vel = strainRate*
for s in bot:
s.shape.color = (1,0,0)
s.state.
#s.state.vel = (0,0,-vel)
for s in top:
s.shape.color = Vector3(0,1,0)
s.state.
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(
v2 = Vector3( rCyl2*cos(
v3 = Vector3( rCyl2*cos(
v4 = Vector3( rCyl2*cos(
f1 = facet((
f2 = facet((
facets.
O.bodies.
mass = O.bodies[
for f in facets:
f.state.mass = mass
f.state.
def addForces():
for f in facets:
n = f.shape.normal
a = f.shape.area
O.forces.
O.dt=.5*
newton=
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
PyRunner(
FlowEngine(
newton,
PyRunner(
]
def plotAddData():
perme = permeability()
f1 = sum(O.forces.
f2 = sum(O.forces.
f = .5*(f2-f1)
s = f/(pi*.
e = (top[0]
plot.addData(
i = O.iter,
s = -s,
e = -e,
tc = interactionLaw.
sc = interactionLaw.
permeab = perme
)
plot.saveDataT
O.step()
ss2sc.interacti
is2aabb.
flow.debug=False
flow.permeabili
flow.fluidBulkM
flow.useSolver=4
flow.permeabili
flow.viscosity= 0.001
flow.decoupleForces = False
flow.boundaryUs
flow.bndCondIsP
flow.bndCondVal
flow.dead=0
flow.emulateAct
#######
numPoints = 100
xs = np.linspace(
ys = np.linspace(
zs = np.linspace(
cellsHit = [] #create array
cylcenterx=
cylcentery=
for x,y,z in itertools.
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(
cellsHit.
for i in cellsHit:
flow.blockCell
flow.setCellPr
O.run(1,1)
flow.meshUpdate
flow.defToleran
numPoints1 = 5
xr = np.linspace(
yr = np.linspace(
zr = np.linspace(
width1=width*0.98
def permeability():
totalVolume = 0
v = np.array([0,0,0])
cellsHit2=[]
for x,y,z in itertools.
cellId2=
if cellId2 in cellsHit2:
continue
if np.sqrt(
cellsHit2.
velocityVector = np.array(
velMag = np.linalg.
cellVol = flow.getCellVol
v = v + cellVol*
totalVolume += cellVol
Cylinder_Vel = np.linalg.
return Cylinder_Vel
plot.plots = { 'e':('s',), 'elateral':('s'),}
plot.plot()
O.run()
#######
[1]https:/
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: