Bug in Calculating getSaturation in TwoPhaseFlow Engine

Asked by Saeed on 2019-09-22

hi,
I had a problem using "TwoPhaseFlow Engine" in 'Yadedaily 2018.02b-290bf6a54e~xenial' and 'Yade-2019.01a'.
The value of Saturation Degree 'getSaturation(False)' was shown properly while doing drainage,but after the dry completed, I close the pressure related to water in bottom of box by 'bndCondIsPressure=[0,0,0,1,0,0]' . I provide Oedometer conditions for the loading.

According to the formula:
Porosity=VolumeOfVoids/BoxVolume
VoidVolumes=Porosity*BoxVolume

The volume of Voids containing the air and Water.
We have blocked the way out water by 'bndCondIsPressure=[0,0,0,1,0,0]' While Loading, So Volume of Water while loading is Constant.
And only the volume of air in the ٰVoid decreases by increasing loading.

According to the formula:
DegreeOfSaturation=VolumeOfWater/VolumeOfVoids

while Oedometer loading:
VolumeOfWater > Constant
VolumeOfVoids > Decresed

So DegreeOfSaturation must be Increased that The same is true with the laboratory results and the literature
but in TwoPhaseFlow Engine Degree of Saturation No changed while loading!!!
-----------------------------------
CODE:

import matplotlib; matplotlib.rc('axes',grid=True)
from yade import pack
import pylab
from numpy import *

compFricDegree = 30.0
confiningS=-1e5
psdSizes,psdCumm=[.599,0.6,0.849,0.85,1.0,1.40],[0.,0.35,0.35,0.70,.70,1.]
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(10,10,10)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=1000,seed=1)

O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=0,density=0,label='frictionless'))

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=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)

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,
 NewtonIntegrator(damping=0.4,label="newton"),
 TwoPhaseFlowEngine(label="unsat")
]

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01:
    break

##################################
## Drainage Test under oedometer conditions #
##################################
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0

meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True
unsat.initialization()
unsat.surfaceTension = 10

for pg in arange(1.0e-5,5.0,0.1):
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  unsat.invasion()
  unsat.computeCapillaryForce()
  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
  print 'Drainage:::','PC:',-unsat.bndCondValue[2],' Sw:',unsat.getSaturation(False),' VoidVolume:',triax.porosity*triax.boxVolume

print '----------LOADING--------------'
## Oedometer conditions
unsat.bndCondIsPressure=[0,0,0,1,0,0]
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=False
loadingMatrix=[-500e3,-1000e3,-1500e3,-2000e3,-2500e3,-3000e3]
for i in arange(0,len(loadingMatrix),1):
 triax.goal2=loadingMatrix[i]
 O.run(2000,True)
 print 'Load:',loadingMatrix[i],' Sw:',unsat.getSaturation(False),' VoidVolume:',triax.porosity*triax.boxVolume
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-------------------------------------------------------------------------
end of code
-------------------------------------------------------------------------

The results I see is wrong:

Load: -500000.0 Sw: 0.268726543553 VoidVolume: 436.590507582
Load: -1000000.0 Sw: 0.268726543553 VoidVolume: 433.836864791
Load: -1500000.0 Sw: 0.268726543553 VoidVolume: 431.014178921
Load: -2000000.0 Sw: 0.268726543553 VoidVolume: 428.158278656
Load: -2500000.0 Sw: 0.268726543553 VoidVolume: 425.283191863
Load: -3000000.0 Sw: 0.268726543553 VoidVolume: 422.394919453
---------------------------------------------------------

It is observed that due to the closure of the water outlet, the volume of Voids decreased but the amount of saturation remained constant during loading, which is incorrect.
It seems that when loading and changing the sample size, the effect of the saturation change is not calculated.

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-09-22
Last reply:
2019-09-27

This question was reopened

Saeed (saeed0668191) said : #3

hi there,
I am adding some codes to the loading section also made no difference and degree of saturation not changed!

loadingMatrix=[-500e3,-1000e3,-1500e3,-2000e3,-2500e3,-3000e3]
for i in arange(0,len(loadingMatrix),1):
 triax.goal2=loadingMatrix[i]
 O.run(2000,True)
 unsat.invasion()
 unsat.computeCapillaryForce()
 for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
 print 'Load:',loadingMatrix[i],' Sw:',unsat.getSaturation(False),' VoidVolume:',triax.porosity*triax.boxVolume

Saeed (saeed0668191) said : #4

Hello to all
Can't anyone help me?

Robert Caulk (rcaulk) said : #5

The getsaturation function uses porevolumes which depend on “computePoreBodyVolume()” and that is only called when you run initialization(). Two options are to re run initialization() each time you want to compute saturation or we need to wrap computePoreBodyVolume(). I’ve never used this engine so i don’t know the consequences of re running initialization() multiple times. Report back, if it doesn’t work I will wrap computePoreVolume() for you so you can call it in python before computing saturation.

[1]https://gitlab.com/yade-dev/trunk/blob/master/pkg/pfv/TwoPhaseFlowEngine.cpp#L3228

Can you help with this problem?

Provide an answer of your own, or ask Saeed for more information if necessary.

To post a message you must log in.