← Back to team overview

yade-users team mailing list archive

Re: [Question #701670]: Wrong permeability values during compaction

 

Question #701670 on Yade changed:
https://answers.launchpad.net/yade/+question/701670

    Status: Answered => Open

Soheil Safari is still having a problem:
Dear Bruno,

Thank you very much for your kind reply.

I fixed the wall number problem and now it works properly.

But I still get weird values for the permeability.

I do not know what the issue is. Because the permeability should have a
really small number like the oedometer.py example. I could gain the
small value just in the first iteration and after that it increased
surprisingly.

# iter		permeability		    porosity		                   time
1000	0.01918070551830009	0.40189088549692115	2.554930150811923
3000	564.9408214756033	0.3889401788504288	2.754930150812345
5000	371.23771416004183	0.3754329114171155	2.954930150812767
7000	418.8671964612367	0.3616597792495562	3.1549301508131893
9000	502.25383359539313	0.34705109690924985	3.3549301508136113
11000	613.7106945556807	0.33216990357515797	3.5549301508140334
13000	679.398330939237	        0.31659321147545794   3.7549301508144555
15000	582.7989699519044	0.29873633244629394	3.9549301508148775
17000	601.5434724622095	0.2823106604884174	4.154930150814612
19000	627.2095474981151	0.26721453576039933	4.354930150814146

Here is my code:

###

from __future__ import print_function
from builtins import range
from yade import pack, plot
from yade import export

num_spheres=400 # number of spheres
young=1e6
compFricDegree = 3 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(1,1,2) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=1e99,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0.1,material='walls')
wallIds=O.bodies.append(walls)
top = walls[5]
bottom = walls[4]
leftb = walls[3]
rightb = walls[1]
rightf = walls[2]

psdSizes,psdCumm=[.01,0.03,0.04,.05,.06,.08,.12],[0.,0.1,0.3,0.3,.3,.7,1.]
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

yade.qt.views()

triax=TriaxialStressController(
	maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
	finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
	thickness = 0,
	stressMask = 7,
	max_vel = 0.005,
	internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=0.2)


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()],label="iloop"
	),
	FlowEngine(dead=1, label="flow"),  #introduced as a dead engine for the moment, see 2nd section
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	newton
]

triax.goal1=triax.goal2=triax.goal3=-10000

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

setContactFriction(radians(finalFricDegree))


triax.dead = True # (!)
top.state.vel = (0, 0, -0.3)
bottom.state.vel = leftb.state.vel = rightb.state.vel = rightf.state.vel =(0,0,0) # zero velocity of all other walls

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead = 0
flow.defTolerance = 0.3
flow.meshUpdateInterval = 200
flow.useSolver = 3
flow.permeabilityFactor = 1
flow.viscosity = 10
flow.bndCondIsPressure = [0, 0, 1, 1, 0, 0]
flow.bndCondValue = [0, 0, 1, 0, 0, 0]
flow.boundaryUseMaxMin = [0, 0, 0, 0, 0, 0]
O.dt = 0.1e-3
O.dynDt = False


O.engines=O.engines+[PyRunner(iterPeriod=2000,command='flow.saveVtk()')]


# export spheres function every 5000 iterations
O.engines += [PyRunner(command="exportSpheres()", iterPeriod=2000, initRun=True)]


def exportSpheres():
    i = O.iter
    fName = f"yade-spheres-{i:06d}.txt"
    export.text(fName)


def get_permeability():
    Qin = flow.getBoundaryFlux(5)
    Qout = flow.getBoundaryFlux(4)
    permeability = abs(Qin) / 1.e-4  #size is one, we compute K=V/∇H
    print("Qin=", Qin, " Qout=", Qout, " permeability=", permeability)
    return permeability

# porosity
O.engines += [PyRunner(iterPeriod=2000,command="plotAddData()",initRun=True)]  # runs  every 5000 iterations

def plotAddData():
    porosity = utils.porosity() 
    permeability = get_permeability()
    plot.addData(iter=O.iter, time=O.time, porosity=porosity, permeability=permeability)
    plot.saveDataTxt("result.txt")
  
# run
O.stopAtIter = 20100
O.run()

###

Many thank in advance
Best regards,
Soheil

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