getSaturation - 2PFV drainage

Hi all,

I am running this script [1], where I have replaced particles by aggregates.
The issue is that I am getting negative values of saturation (getSaturation) in the end of the drainage.

Does this have to do with pack boundaries?

[1]

#!/usr/bin/python
# -*- encoding=utf-8 -*-
#*************************************************************************
# Copyright (C) 2010 by Bruno Chareyre *
# bruno.chareyre_at_grenoble-inp.fr *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************/

from yade import pack
from yade import bodiesHandling
from yade import export
from yade import utils
from yade import ymport

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
    num_spheres=3000,# number of spheres
    compFricDegree = 1, # contact friction during the confining phase
    key='_triax_base_', # put you simulation's name here
    unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.37 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=0 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=110e4 # contact stiffness
bondstr=1.5e4
mn,mx=Vector3(0,0,0),Vector3(0.01,0.01,0.01) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,tensileStrength=bondstr,cohesion=bondstr,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres'))
O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(0),density=2600,tensileStrength=0,cohesion=0,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=0,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=0.0006,num=num_spheres,seed=1) #"seed" make the "random" generation always the same rRelFuzz=1,num=num_spheres,
#O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

for p in sp:
 a=p[0]
# print (a)
 f1=O.bodies.append(ymport.textExt("aggcoating2.txt", format='x_y_z_r', shift=a-Vector3(0,0,0.0006), scale=1.0,material='spheres',color=(0,1,1)))

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
    ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
    ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
    maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
    thickness = 0,
    ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
    ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
    ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
    ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
    ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
    stressMask = 0,
    internalCompaction=False, # If true the confining pressure is generated by growing particles
    wall_front_activated=True,
    wall_back_activated=True,
    wall_top_activated=True,
    wall_bottom_activated=True,
    wall_left_activated=True,
    wall_right_activated=True,
    goal1=-20,
    goal2=-20,
    goal3=-20,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[

 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=-1, xSectionWeibullShapeParameter=3, xSectionWeibullScaleParameter=1)],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=False)]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
# VTKRecorder(iterPeriod=2000,recorders=['all','cracks'],fileName="/home/user/Área de Trabalho/Paper_Biopore/SandBox_Laser/Video/Videotest"),
# newton
# NewtonIntegrator(damping=0.4,gravity=[0,0,0]),
 newton

]

#Display spheres with 2 colors for seeing rotations better
#Gl1_Sphere.stripes=0
#if nRead==0: yade.qt.Controller(), yade.qt.View()

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

## We will reach a prescribed value of porosity with the REFD algorithm
## (see http://dx.doi.org/10.2516/ogst/2012032 and
## http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)

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
# compFricDegree = 0.95*compFricDegree
    setContactFriction(radians(compFricDegree))
    print ("\r Friction: ",compFricDegree," porosity:",triax.porosity,
    sys.stdout.flush())
    # while we run steps, triax will tend to grow particles as the packing
    # keeps shrinking as a consequence of decreasing friction. Consequently
    # porosity will decrease
    O.run(500,1)

#O.save('compactedState'+key+'.yade.gz')
print ("### Compacted state saved ###")
print(triax.stress(3)[1])

##############################
### DEVIATORIC LOADING ###
##############################

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
#triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 2
triax.wall_bottom_activated=0
#now goal2 is the target strain rate
triax.goal1=rate#triax.stress(1)[0]
triax.goal3=rate#triax.stress(5)[2]
triax.goal2=triax.stress(3)[1]

##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
#O.saveTmp()

#####################################################
### Example of how to record and plot data ###
#####################################################

#from yade import plot
from yade import plot
O.run(1000,True)

#strain is logarithmic strain or true strain which is ls=(ln1+e) where e=dl/L (strain)
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]

## a function saving variables
def history():
    plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            s11=-triax.stress(triax.wall_right_id)[0]-si0,
            s22=-triax.stress(triax.wall_top_id)[1]-si1,
            s33=-triax.stress(triax.wall_front_id)[2]-si2,
            pc=-unsat.bndCondValue[2],
            sw=unsat.getSaturation(False),
            i=O.iter)

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

#plot.plots={'pc':('sw',None,'e22')}
#plot.plot()

#######################################################
## Drainage Test under oedometer conditions ###
#######################################################
##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.

##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,-1e8,0,0,0]
unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
unsat.initialization()
unsat.surfaceTension = 0.0728

#start invasion, the data of normalized pc-sw-strain will be written into pcSwStrain.txt
f1=open('Cells1.txt',"w")
f2=open('Cells2.txt',"w")
f3=open('Cells3.txt',"w")
f4=open('Cells4.txt',"w")
f5=open('SwPc.txt',"w")

ts=O.dt
pgstep= 10000000*ts
pgmax= 9316

for pg in arange(1.0e-8,pgmax,pgstep):
  unsat.bndCondValue=[0,0,(-1.0)*pg,0,0,0]

#for pg in arange(1.0e-5,10,0.01):
# #increase gaz pressure at the top boundary
# unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  #compute the evolution of interfaces
  unsat.invasion()
  #save the phases distribution in vtk format, to be displayed by paraview
# unsat.savePhaseVtk("/home/user/Área de Trabalho/Paper_Biopore/SandBox_Laser/Video/Pressurefield")
  #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))

  if pg==0.60001:
    cels=unsat.nCells()
    celsW1 = [0.0]*cels
    celsV1 = [0.0]*cels
    celsBar1 = [0.0]*cels
    for ii in range(cels):
      celsW1=unsat.getCellIsWRes(ii)
      celsV1=unsat.getCellVolume(ii)
# celsBar1=unsat.getCellBarycenter(ii)
      celsBar1=unsat.getCellCenter(ii)
      f1.write(str(celsW1)+" "+str(celsV1)+" "+str(celsBar1)+"\n")
    f1.close()

  if pg==2.30001:
    cels2=unsat.nCells()
    celsW2 = [0.0]*cels2
    celsV2 = [0.0]*cels2
    celsBar2 = [0.0]*cels2
    for jj in range(cels2):
      celsW2=unsat.getCellIsWRes(jj)
      celsV2=unsat.getCellVolume(jj)
      celsBar2=unsat.getCellCenter(jj)
      f2.write(str(celsW2)+" "+str(celsV2)+" "+str(celsBar2)+"\n")
    f2.close()

  if pg==3.10001:
    cels3=unsat.nCells()
    celsW3 = [0.0]*cels3
    celsV3 = [0.0]*cels3
    celsBar3 = [0.0]*cels3
    for gg in range(cels3):
      celsW3=unsat.getCellIsWRes(gg)
      celsV3=unsat.getCellVolume(gg)
      celsBar3=unsat.getCellCenter(gg)
      f3.write(str(celsW3)+" "+str(celsV3)+" "+str(celsBar3)+"\n")
    f3.close()

  if pg==9.90001:
    cels4=unsat.nCells()
    celsW4 = [0.0]*cels4
    celsV4 = [0.0]*cels4
    celsBar4 = [0.0]*cels4
    for hh in range(cels4):
      celsW4=unsat.getCellIsWRes(hh)
      celsV4=unsat.getCellVolume(hh)
      celsBar4=unsat.getCellCenter(hh)
      f4.write(str(celsW4)+" "+str(celsV4)+" "+str(celsBar4)+"\n")
    f4.close()

  #reac
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    if unb<0.001:
      break
  f5.write(str(pg)+" "+str(unsat.getSaturation(False))+" "+str(triax.strain[1])+"\n")
f5.close()

plot.plots={'pc':('sw',None,'e22')}
plot.plot(noShow=True).savefig('Fig2.png')

exit()

Aggregate - aggcoating2.txt
#format x_y_z_r
3.76558e-05 -8.41309e-05 0.000908592 0.0001
-0.00024319 -0.000125352 0.000435266 0.0001
0.000294296 -0.000222202 0.000612602 0.0001
0.000418583 -0.000219287 0.000751958 0.0001
-0.000226268 2.43074e-05 0.000783996 0.0001
0.000160139 -0.000294018 0.000350308 0.0001
0.000433878 0.000103233 0.000418752 0.0001
-0.000158174 -0.000227146 0.000681728 0.0001
-0.000194443 0.000125398 0.000639618 0.0001
-0.000251319 -0.00015249 0.000787068 0.0001
-0.000250173 9.64227e-05 0.000907745 0.0001
-0.00030227 0.000287221 0.000455403 0.0001
0.000485042 -3.575e-05 0.000508324 0.0001
0.000302121 0.000388348 0.000537567 0.0001
0.000326139 3.14592e-05 0.000816209 0.0001
0.000122141 -0.000162135 0.000864344 0.0001
-3.96352e-05 0.000354408 0.000527824 0.0001
-0.000201977 1.54357e-05 0.000227817 0.0001
-7.4623e-05 -3.98969e-05 0.00100381 0.0001
-0.000147642 0.000100428 0.00104245 0.0001
-0.000256531 -0.000378955 0.000586158 0.0001
0.000308172 8.11204e-05 0.00030273 0.0001
8.47225e-07 -0.000213415 0.00101332 0.0001
-0.000132391 -9.4304e-05 0.000806951 0.0001
0.000174136 0.000319067 0.000820194 0.0001
0.000326592 8.61222e-05 0.00050047 0.0001
-0.000253401 -0.000262311 0.000386586 0.0001
0.000202625 1.75135e-05 0.000953477 0.0001
-0.000130475 3.49049e-05 0.00051105 0.0001
0.000156978 0.000378562 0.000465153 0.0001
0.000112672 -9.94457e-05 0.000153187 0.0001
5.01333e-05 8.32437e-05 0.00103545 0.0001
0.000147004 0.000241594 0.000687315 0.0001
0.000148048 -2.97394e-06 0.000723499 0.0001
0.000247982 -0.000258142 0.000833374 0.0001
0.000261184 -0.00033164 0.000495643 0.0001
-0.000290397 -0.000181989 0.000588953 0.0001
-9.28162e-05 0.000422466 0.000672301 0.0001
1.25775e-05 -0.000269665 0.000298369 0.0001
-0.000119542 -0.000291999 0.000840892 0.0001
-0.000111854 -0.000372176 0.000388634 0.0001
-2.60069e-05 0.000176844 0.000966913 0.0001
2.60352e-05 0.00027016 0.000406357 0.0001
-0.000250068 0.000313828 0.000723481 0.0001
3.08219e-05 0.000169685 0.00025657 0.0001
0.000151379 -0.000350054 0.000568136 0.0001
0.00020777 -5.67146e-05 0.000843063 0.0001
-8.30843e-05 0.000277365 0.000270681 0.0001
-9.0589e-05 -0.000208896 0.000385832 0.0001
0.00025288 4.20635e-06 0.000399811 0.0001
-0.00035799 -0.000264158 0.000514701 0.0001
-0.000351678 -0.000106289 0.000683854 0.0001
2.42899e-05 8.92898e-05 0.000420469 0.0001
-0.000428749 5.3352e-05 0.000448192 0.0001
0.000399031 0.000243883 0.000569854 0.0001
0.000176455 0.000413723 0.000659182 0.0001
-0.000355453 0.00017135 0.000545837 0.0001
-0.000194565 -0.000209788 0.000967668 0.0001
4.11024e-05 -0.000182733 0.00048278 0.0001
0.000277921 -0.000285809 0.00030946 0.0001
-0.000363441 0.000104945 0.000729084 0.0001
-0.000367004 -2.7924e-05 0.000838584 0.0001
0.000272903 -0.000336054 0.000688 0.0001
6.99681e-05 0.000222928 0.000539046 0.0001
0.000163231 -4.72276e-05 0.00045087 0.0001
0.000119068 -0.000147811 0.000329868 0.0001
7.52657e-05 0.00012236 0.000692622 0.0001
3.25153e-05 0.00031401 0.000720474 0.0001
-0.000115277 -0.000472376 0.000544238 0.0001
-6.3985e-05 3.7724e-05 0.000694756 0.0001
-3.96344e-05 0.000146132 0.000556236 0.0001
0.000437613 7.38593e-05 0.000608536 0.0001
-5.42288e-05 -0.000305212 0.000545797 0.0001
0.000287961 -3.86083e-05 0.0005703 0.0001
-6.72336e-06 -0.000359846 0.000673851 0.0001
-0.000409214 0.000168504 0.000373636 0.0001
0.000146369 -8.9648e-05 0.000612921 0.0001
0.000382267 -0.000113681 0.00057611 0.0001
-0.000140111 0.000370419 0.00041522 0.0001
0.000153788 0.000229417 0.000331849 0.0001
-3.32316e-05 -0.000418565 0.000833671 0.0001
-0.000171179 -0.00028338 0.000528477 0.0001
-0.000166367 -6.22616e-05 0.000644039 0.0001
-0.000257597 0.000106314 0.00036711 0.0001
-0.000402476 -0.000141045 0.000379145 0.0001
0.000438376 0.000170724 0.000711756 0.0001
1.36135e-05 -6.6382e-05 0.000755 0.0001
0.000195611 0.000135453 0.000815171 0.0001
0.000151216 -0.000469842 0.000641459 0.0001
0.000158198 4.34732e-05 0.000262611 0.0001
0.00038393 -9.57436e-05 0.00074336 0.0001
3.08201e-05 0.000384211 0.000349408 0.0001
-7.32696e-05 0.000144854 0.000775437 0.0001
0.000265532 0.000151693 0.000641517 0.0001
-0.000461932 5.01701e-05 0.000570563 0.0001
0.000286433 0.000210237 0.000443172 0.0001
0.00021887 0.000189744 0.00022012 0.0001
-7.96387e-05 -0.000167518 0.000899771 0.0001
2.77186e-05 0.000448184 0.00078581 0.0001
-6.95824e-05 -9.4496e-05 0.000330986 0.0001
0.000126438 -0.000213294 0.000731194 0.0001
1.96227e-05 -1.29656e-05 0.00058299 0.0001
0.000295241 4.49876e-05 0.00069004 0.0001
-6.59752e-05 0.000248671 0.000669664 0.0001
-0.000178501 0.000269147 0.000573993 0.0001
-0.000234289 0.000224858 0.000312635 0.0001
-0.000303332 1.95424e-05 0.000632745 0.0001
0.0002551 -0.000143107 0.000720083 0.0001
-0.000349276 0.000221435 0.000833498 0.0001
0.000111761 -0.000247921 0.000609618 0.0001
-0.000297051 -0.000286166 0.000715954 0.0001
-0.000231734 -0.000130819 0.000275358 0.0001
0.000153383 0.000127635 0.000409071 0.0001
0.000313707 0.000199103 0.000819924 0.0001
-8.08601e-05 0.000117801 0.000342111 0.0001
-0.000221782 0.000388779 0.000576177 0.0001
0.00029968 0.000293505 0.000685586 0.0001
-0.000141029 -4.74788e-06 0.000365925 0.0001
-0.000205918 0.000197864 0.000775098 0.0001
0.000311767 -0.000184511 0.000431775 0.0001
-0.00031711 -0.000199012 0.000879434 0.0001
-0.000418758 -8.82513e-05 0.000517954 0.0001
-8.6727e-05 -6.1194e-05 0.000159197 0.0001
-0.000141739 -0.000253983 0.000238511 0.0001
-0.000160239 0.00039061 0.000768023 0.0001
-0.000333443 0.000230133 0.000641908 0.0001
0.000185753 -0.000167971 0.00100093 0.0001
0.000214614 0.000350336 0.000352933 0.0001
-1.35542e-05 -0.000239018 0.000766767 0.0001
3.3344e-05 -0.000427444 0.000529541 0.0001
-0.000164105 0.000227587 0.000909519 0.0001
0.000322746 -0.000133394 0.000897852 0.0001
-0.000141818 -0.000385113 0.000692781 0.0001
5.48302e-05 -1.58745e-05 0.000303377 0.0001
-0.000163875 0.000136352 0.000191111 0.0001
0.000242334 -0.000111171 0.000271347 0.0001
-0.000227601 -5.7618e-05 0.000959154 0.0001
-0.000284654 -2.08033e-05 0.000507157 0.0001
0.000117134 -0.000366719 0.000762292 0.0001
-0.000331363 -6.50477e-06 0.000343825 0.0001
0.000173963 0.000187624 0.000952326 0.0001
-0.000123856 -0.000128282 0.000516013 0.0001
9.89772e-05 -3.565e-05 0.00104006 0.0001
1.71099e-05 -4.82906e-05 0.000462197 0.0001
4.79151e-05 0.000337391 0.000934679 0.0001
-2.73732e-05 -0.000158971 0.000637287 0.0001
6.75922e-05 -0.000294129 0.000905222 0.0001
-8.95706e-05 4.79077e-05 0.000883271 0.0001
-0.000476521 -1.67566e-05 0.000702634 0.0001
-7.82808e-05 0.000320102 0.000843521 0.0001
-3.63076e-06 -0.000162705 0.000206596 0.0001
6.26698e-05 0.000384733 0.000587514 0.0001
4.4008e-05 0.000210102 0.000845065 0.0001
-2.19028e-05 4.43025e-05 0.00018558 0.0001
6.29598e-05 5.68316e-05 0.000864066 0.0001
0.000115839 5.154e-05 0.000124301 0.0001
0.000207907 0.000266118 0.000533277 0.0001
0.000153701 8.40376e-05 0.000558781 0.0001
-0.000120721 0.000219598 0.000441192 0.0001
-0.000240225 0.000128607 0.000483583 0.0001
0.000473712 -2.25817e-05 0.000692164 0.0001
5.69195e-05 -0.000309473 0.000428843 0.0001
0.00019845 -0.000181459 0.000477741 0.0001
0.000374005 -7.0736e-05 0.000393175 0.0001
-0.000226909 -0.00034777 0.000833712 0.0001