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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
