yade-users team mailing list archive

[Ziyu Wang <question702316@xxxxxxxxxxxxxxxxxxxxx>]

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

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

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