Two phase flow

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