flow.savePoreNetwork()

Hi,

The answer is in your question. savePoreNetwork() is not a member of FlowEngine (but of TwoPhaseFlowEngine for instance, as you mention it). So, if "flow" stands for a FlowEngine instance (impossible to know without a Minimal Working Example), flow.savePoreNetwork() can not work.

You have to define a TwoPhaseFlowEngine "flow" guy in your YADE simulation. Making it work for e.g. just one iteration would possibly not hurt, even though you wanna simulate saturated media.

Dear Jérôme,

Thank you so much for your quick response and valuable information.

Would you please let me know how I can define a TwoPhaseFlowEngine "flow" based on your recommendation? Where should I use this?

I attached the codes for your consideration:

TwoPhaseFlowInjection.py [2] ‘TwoPhaseFlowEngine’:

##
import os
import numpy as np
import pandas as pd
from yade import pack, ymport
from yade import utils, plot, timing

num_spheres=1000# number of spheres
young=1.e6
compFricDegree = 3 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(1.,1.,0.4) # corners of the initial packing
graindensity=2600
errors=0
toleranceWarning =1.e-11
toleranceCritical=1.e-6

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=graindensity,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

O.engines = [
  TwoPhaseFlowEngine(dead=0,label="flow"),
]

press=1000.

flow.dead=0
flow.meshUpdateInterval=-1
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=0.1

flow.bndCondIsWaterReservoir=[0,0,1,0,0,0]

flow.bndCondIsPressure=[0,0,1,0,0,0]
flow.bndCondValue=[0,0,press,0,0,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
flow.iniVoidVolumes=True

flow.surfaceTension = 0.0
flow.drainageFirst=False
flow.isDrainageActivated=False
flow.isImbibitionActivated=True
flow.isCellLabelActivated=True
flow.initialization()
cs=flow.getClusters()
c0=cs[1]

#flow.getCellVoidVolume(4)
#flow.getCellInSphereRadius(100)
flow.savePoreNetwork() #extract info about the pore

##

oedometer.py PFV model [1] 'FlowEngine' :

from __future__ import print_function
from builtins import range
from yade import pack
from yade import export

num_spheres=1000# number of spheres
young=1e6
compFricDegree = 3 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

psdSizes,psdCumm=[.01,0.03,0.04,.05,.06,.08,.12],[0.,0.1,0.3,0.3,.3,.7,1.]
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

yade.export.text('x.txt', mask=-1)

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 max_vel = 0.005,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=0.2)

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()],label="iloop"
 ),
 FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

triax.goal1=triax.goal2=triax.goal3=-10000

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

setContactFriction(radians(finalFricDegree))

## ______________ Oedometer section _________________

#A. Check bulk modulus of the dry material from load/unload cycles
triax.stressMask=2
triax.goal1=triax.goal3=0

triax.internalCompaction=False
triax.wall_bottom_activated=False
#load
triax.goal2=-11000; O.run(2000,1)
#unload
triax.goal2=-10000; O.run(2000,1)
#load
triax.goal2=-11000; O.run(2000,1)
e22=triax.strain[1]
#unload
triax.goal2=-10000; O.run(2000,1)

e22=e22-triax.strain[1]
modulus = 1000./abs(e22)

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=10
flow.bndCondIsPressure=[0,0,1,1,0,0]
flow.bndCondValue=[0,0,1,0,0,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
O.dt=0.1e-3
O.dynDt=False

O.run(1,1)
Qin = flow.getBoundaryFlux(2)
Qout = flow.getBoundaryFlux(3)
permeability = abs(Qin)/1.e-4 #size is one, we compute K=V/∇H
print("Qin=",Qin," Qout=",Qout," permeability=",permeability)

#C. now the oedometer test, drained at the top, impermeable at the bottom plate
flow.bndCondIsPressure=[0,0,0,1,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
flow.updateTriangulation=True #force remeshing to reflect new BC immediately
newton.damping=0

#we want the theoretical value from Terzaghi's solution
#keep in mind that we are not in an homogeneous material and the small strain
#assumption is not verified => we don't expect perfect match
#there can be also an overshoot of pressure in the very beginning due to dynamic effects
Cv=permeability*modulus/1e4
zeroTime=O.time
zeroe22 = - triax.strain[1]
dryFraction=0.05 #the top layer is affected by drainage on a certain depth, we account for it here
drye22 = 1000/modulus*dryFraction
wetHeight=1*(1-dryFraction)

def consolidation(Tv): #see your soil mechanics handbook...
 U=1
 for k in range(50):
  M=pi/2*(2*k+1)
  U=U-2/M**2*exp(-M**2*Tv)
 return U

triax.goal2=-11000

from yade import plot

## a function saving variables
def history():
 plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000)
 #plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
##make nice animations:
#O.engines=O.engines+[PyRunner(iterPeriod=200,command='flow.saveVtk()')]

from yade import plot
plot.plots={'t':(('e22','b--'),('e22_theory','b-'),None,('s22','g--'),('p','g-'))}
plot.plot()
O.saveTmp()
O.timingEnabled=1
from yade import timing
print("starting oedometer simulation")
O.run(200,1)
timing.stats()

print("\nPress ▶ (the start button) to see graph.\n")

###

Many thanks in advance.

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/FluidCouplingPFV/oedometer.py
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/FluidCouplingPFV/TwoPhaseFlowInjection.py

This question was expired because it remained in the 'Open' state without activity for the last 15 days.