← Back to team overview

yade-users team mailing list archive

[Question #699603]: getSaturation - 2PFV drainage

 

New question #699603 on Yade:
https://answers.launchpad.net/yade/+question/699603

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





-- 
You received this question notification because your team yade-users is
an answer contact for Yade.