Periodic flow and permeability

Asked by Jibril Coulibaly on 2019-11-08

I am trying to perform the seemingly simple computation of the permeability of a periodic assembly using the PFV module. From that previous question (https://answers.launchpad.net/yade/+question/662105), I gather that I need to define a macroscopic pressure gradient, as well as to fix the pressure in a given cell.

Trying to do so, Yade complains "setCellPImposed" and "setCellPressure" being out of range with max value 0 and returns wrong values for my boundary flows. Following is my input script, I would very much appreciate your help on figuring this out:

###############################################

O.periodic=True
O.cell.hSize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
sp=pack.SpherePack()
O.bodies.append([sphere(s[0],s[1]) for s in sp])

O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
NewtonIntegrator(damping=0.),
]

O.dt=0.1e-8
O.dynDt=False

flow.defTolerance=0.3
flow.meshUpdateInterval=300
flow.useSolver=0 # Solver >0 not supported for PBC
flow.permeabilityFactor=1
flow.viscosity=1
flow.setCellPImposed(0,True)
flow.setCellPressure(0,0)

O.run(1,1)
Qin = flow.getBoundaryFlux(0)
Qout = flow.getBoundaryFlux(1)
print("Qin=",Qin," Qout=",Qout)
###############################################

Thank you,
Jibril B. Coulibaly

Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2020-02-03
2020-02-04
 Robert Caulk (rcaulk) said on 2019-11-09: #1

Hello Jibril,

I'm not sure you'll find solace in my answer, but a quick review of your post lead me to these two points:

1st point: In the MWE you present, you are indexing cells before a triangulation is created. FlowEngine builds the triangulation during the first time step that it is active. So it would be necessary to make a step (O.run(1,1)) before assigning pressures. Alternatively, you could use flow.emulateAction()[1] to create the triangulation without running a mechanical timestep.

2nd point:
># Solver >0 not supported for PBC

Should not be true for yadedaily versions and compiled trunk versions.

Cheers,

Robert

 Jibril Coulibaly (jibril-coulibaly) said on 2019-11-09: #2

Hello Robert,

I am running Yade 2018.02b that I installed running "sudo apt-get install yade", that must be why I do not have access to "solver>0".

Regarding the building of the triangulation, moving the "O.run(1,1)" before assigning pressures indeed stops the error. However, it still yields zero boundary flux, so I must be missing a piece of the setup still.

I tried running a second time using "O.run(1,1)" after assigning the pressures, hoping it would run the system with the properly assigned pressures but I get a segmentation fault instead.

Jibril B. Coulibaly

 Launchpad Janitor (janitor) said on 2019-11-25: #3

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

 Robert Caulk (rcaulk) said on 2019-11-25: #4

 Bruno Chareyre (bruno-chareyre) said on 2019-11-25: #5

Hi,
A periodic problem has no boundaries and therefore "flow.getBoundaryFlux(1)" is meaningless (Body #1 is not a boundary anyway).
What you want is "flow.averageVelocity()".

In order to impose pressure at (x,y,z) coordinates you need flow.imposePressure(). Better avoid playing with mesh elements directly since they change all the time.
Cheers
Bruno

 Jibril Coulibaly (jibril-coulibaly) said on 2020-02-03: #6

Hi,

Thank you very much rcaulk and bruno-chareyre for your answers. Using flow.imposePressure() and flow.averageVelocity() provide sensible results.

I have a couple remaining questions:

- (1) : If I try to run the script in a loop to obtain permeability in all 3 directions, I get sensible results on the first run, then I get zero or nan on subsequent runs. I did not get these issues when using walls (see script below)

- (2) : Running a similar script, but reading the particles coordinates from an external file instead, I obtain the error message: "GS did not converge in 20k iterations (maybe because the reference pressure is zero?)".
I tried to figure it out from the code but I am not sure what this message is warning me about.

Jibril B. Coulibaly

###########

O.periodic=True
O.cell.hSize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
pos_center = Vector3(0.05,0.05,0.05)
sp=pack.SpherePack()
O.bodies.append([sphere(s[0],s[1]) for s in sp])

O.engines=[
ForceResetter(),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
NewtonIntegrator(damping=0.),
]

O.dt=0.1e-8
O.dynDt=False

flow.defTolerance=0.3
flow.meshUpdateInterval=300
flow.useSolver=0 # Solver >0 not supported for PBC
flow.permeabilityFactor=1
flow.viscosity=1
flow.imposePressure(pos_center,1.0)

for i in [0,1,2]:
if i==0:
elif i==1:
elif i==2:

flow.updateTriangulation=True
flow.updateTriangulation=True
O.run(1,1)
aveVel = flow.averageVelocity()
print("average velocity=(",aveVel[0],aveVel[1],aveVel[2],")")

#######################

 Bruno Chareyre (bruno-chareyre) said on 2020-02-04: #7

Hi,
Your script reveals that writing "flow.updateTriangulation=True" repeatedly, at each iteration, leads to problems with periodic solver. It needs to be investigated. Most likely some internals are not updated correctly after assigning gradP multiple times.

In fact changing the boundary conditions repeatedly at one point in time (i.e. in 0 iteration) is not something supposed to happen in a concrete simulation and probably it was not anticipated.

Nevertheless, you can always use three engines instead of one, so you are sure there is no such issue. Below is a shorter, working, version of your script. Results with useSolver=3 (default):

average velocity= Vector3(-4.537650575699425e-6,-4.7346547328556296e-8,-5.300733352662439e-8)
average velocity= Vector3(-4.734654732786302e-8,-4.190262426966041e-6,5.93967279201477e-8)
average velocity= Vector3(-5.300733352643585e-8,5.9396727919603324e-8,-4.479154269170494e-6)

Bruno

O.periodic=True
O.cell.hSize=Matrix3(0.1,0,0, 0,0.1,0, 0,0,0.1)
pos_center = Vector3(0.05,0.05,0.05)
sp=pack.SpherePack()
O.bodies.append([sphere(s[0],s[1]) for s in sp])

for p in [Vector3(1,0,0),Vector3(0,1,0),Vector3(0,0,1)]:
flow=PeriodicFlowEngine()
flow.imposePressure(pos_center,0)
flow.emulateAction()
flow.saveVtk()
print("average velocity=",flow.averageVelocity())

 Bruno Chareyre (bruno-chareyre) said on 2020-02-04: #8

The problem is reported here: https://gitlab.com/yade-dev/trunk/issues/141