New question #641624 on Yade:
https://answers.launchpad.net/yade/+question/641624
Hi
I am using this two phase flow engine to simulate the behaviour of unsaturated
soils. So, i am running a simple script in which , for the first objective i am
desaturating my sample to say some 80% and then applying an axial loading. What
i am interested is to simulate something very similar to 1 d loading of
unsaturated soil, under constant suction.
I reached my first objective, but for the next step i.e., loading...my program
just stalls without doing anything (i neither get any kind of error message nor
my prog terminates just the virtual time stops and the real time is still
running). I appreciate some further help in this regard.
import matplotlib; matplotlib.rc('axes',grid=True)
from yade import pack, qt
import pylab
from numpy import *
targetPorosity=0.40
friction=0.6
angle=atan(friction)
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.5,0.5,0.5)
sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=0.0005,num=10000,seed=1)
## create material #0, which will be used as default
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=angle,density=2600,label='spheres'))
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=0,density=0,label='frictionless'))
## create walls around the packing
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in
sp])
triax=TriaxialStressController(
internalCompaction=True,
goal1=-10000,
goal2=-10000,
goal3=-10000,
max_vel=10,
label="triax"
)
newton=NewtonIntegrator(damping=0.4)
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()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
newton
]
while 1:
O.run(1000,True)
unb=unbalancedForce()
print triax.meanStress
if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
break
print "out"
import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
## we decrease friction value and apply it to all the bodies and
contacts
friction = 0.95*friction
setContactFriction(friction)
print "\r Friction: ",friction," porosity:",triax.porosity,
sys.stdout.flush()
O.run(500,1)
#############################
## REACH NEW EQU. STATE ###
#############################
triax.internalCompaction=False
setContactFriction(0.5)
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.001 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
break
print "out of loop"
triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
O.save('1kPacking.yade') #save the packing, which can be reloaded later.
print "package saved"
O.run(1000,True)
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]
from yade import plot
O.engines=O.engines+[PyRunner(iterPeriod=20,command='history()',dead=1,label='recorder')]
def history():
plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1,
e33=-triax.strain[2]-ei2,
s11=-triax.stress(0)[0]-si0,
s22=-triax.stress(2)[1]-si1,
s33=-triax.stress(4)[2]-si2,
pc=-unsat.bndCondValue[2],
sw=unsat.getSaturation(True),
i=O.iter
)
plot.plots={'pc':('sw',None,'e22')}
plot.plot()
#######################################################
## Drainage Test under oedometer conditions ###
#######################################################
##oedometer conditions
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0
recorder.dead=0
##Instantiate a two-phase engine
unsat=UnsaturatedEngine()
##set boundary conditions, the drainage is controlled by decreasing W-phase
pressure and keeping NW-phase pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-50000,50000,0,0]
unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
unsat.initialization()
unsat.surfaceTension = 74
##start invasion, the data of normalized pc-sw-strain will be written into
pcSwStrain.txt
file=open('pcSwStrain.txt',"w")
for pg in arange(-50000,-25000,1000):
#increase gaz pressure at the top boundary
unsat.bndCondValue=[0,0,(-1.0)*pg,50000,0,0]
#compute the evolution of interfaces
unsat.invasion()
#save the phases distribution in vtk format, to be displayed by paraview
unsat.savePhaseVtk("./")
#compute and apply the capillary forces on each particle
unsat.computeCapillaryForce()
for b in O.bodies:
O.forces.setPermF(b.id, unsat.fluidForce(b.id))
#reac
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.01:
break
file.write(str(pg)+" "+str(unsat.getSaturation(True))+"
"+str(triax.strain[1]-ei1)+"
"+str(unsat.bndCondValue[3]-unsat.bndCondValue[2])+"\n")
file.close()
triax_iso=TriaxialStressController(
stressMask=2,
wall_bottom_activated=0,
internalCompaction=False,
goal1=0.,
goal2=-500000,
goal3=0.,
max_vel=0.001,
label="triax_iso"
)
pc_ini = (unsat.bndCondValue[3]-unsat.bndCondValue[2])
print pc_ini
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()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax_iso,
PyRunner(command='capillary()',iterPeriod=1),
TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
NewtonIntegrator(damping=0.4),
PyRunner(command='addPlotData()',iterPeriod=1000),
PyRunner(command='saveAddData()',iterPeriod=1000)
]
def capillary():
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndValue=[0,0,-pc_ini,0,0,0]
unsat.computeCapillaryForce()
Cap = '%f_cap.txt' %(O.iter)
f2=open(Cap, 'w')
for b in O.bodies:
O.forces.setPermF(b.id, unsat.fluidForce(b.id))
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.01:
break
f2.write(str(O.iter)+" "+str(pg)+" "+str(unsat.getSaturation(True))+"
"+str(triax.strain[1]-ei1)+" "+str(-unsat.bndCondValue[2])+"\n")
def addPlotData():
contact_stress=utils.getStress()
plot.addData(
i=O.iter,
S11=-triax_iso.stress(0)[0],S22=-triax_iso.stress(2)[1],S33=-triax_iso.stress(4)[2],
E11=-triax_iso.strain[0],E22=-triax_iso.strain[1],E33=-triax_iso.strain[2],
poro=utils.porosity(),
stress=-triax_iso.meanStress,stressxx=contact_stress[0][0],
stressxy=contact_stress[0][1], stressxz=contact_stress[0][2],
stressyx=contact_stress[1][0], stressyy=contact_stress[1][1],
stressyz=contact_stress[1][2], stresszx=contact_stress[2][0],
stresszy=contact_stress[2][1], stresszz=contact_stress[2][2],
pc_odeo=-unsat.bndCondValue[2], sat=unsat.getSaturation(False)
)
name = '%f.txt' %(O.iter)
f=open(name, 'w')
for i in O.interactions:
if i.id1>5:
par1 = O.bodies[i.id1].state.pos
par2 = O.bodies[i.id2].state.pos
fn = i.phys.normalForce.norm()
f.write(str(O.iter)+' '+str(i.id1)+' '+str(i.id2)+'
'+str(par1[0])+' '+str(par1[1])+' '+str(par1[2])+' '+str(par2[0])+'
'+str(par2[1])+' '+str(par2[2])+' '+str(fn)+'\n')
#set time step and run simulation:
O.dt=0.5*PWaveTimeStep()
plot.plots={'i ':('S11','S22','S33'),' i':('E11','E22','E33'),
'stress':('poro')}
#show the plot
plot.plot()
def saveAddData():
#save the plot
plot.saveDataTxt('Data_stress.txt',vars=('i','stressxx','stressxy','stressxz','stressyx','stressyy','stressyz','stresszx','stresszy','stresszz','pc_odeo','sat','poro'))
Thanks
Amiya
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
