yade-users team mailing list archive

[Amiya Prakash Das <question641624@xxxxxxxxxxxxxxxxxxxxx>]

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.