boundary thickness is 0.001 by default in FlowEngine

[Christian Jakob]

Hi,

I want to get the pore pressure of a saturated soil sample in an undrained triaxial test when the sample volume is decreasing.

...
TriaxialStressController(thickness=0,stressMask=7,internalCompaction=False,label='triax')
...

 triax.goal1 = 35000
 triax.goal2 = 35000
 triax.goal3 = 70000
 triax.wall_back_activated = False
...
 flow.viscosity = 1.3e-3
 flow.fluidBulkModulus = 2e9
 flow.bndCondValue=[0,0,0,0,0,0]
 flow.bndCondIsPressure=[0,0,0,0,0,0]
...

For one of my simulations a decrease of ca. 1,8 % in porosity (=-dV/V) is obtained. When I calculate corresponding pressure increase I get

dp = K*(-dV/V) = 2e9*0.018 = 36000000 Pa

But the output from PFV shows maximum average pressures around of 16 Pa (not MPa!). As long as I use SI units the output of flow.averagePressure() should be Pa, right?

So is there a problem in the code or did I miss something in the compressibility option of PFV?

regards,

christian

###additional info:

I double-checked my script and did not see anything wrong. So I think it is on C++ side. So I was going on testing PFV method:

I also tested a drained situation with incompressible fluid. The porosity is decreasing due to a triggered collapse of the grain structure, but the pore pressure gets negative (see attachment).
So I think it is a bug in pressure calculation of PFV method.

[Bruno Chareyre (bruno-chareyre)]

Hi Christian,
Could you please show how to reproduce this in a simple case?
I didn't find any results like what you report before.

[Christian Jakob (jakob-ifgt)]

Hi Bruno,

I am already working on a small script for reproducing, but it will take a while ...

[Christian Jakob (jakob-ifgt)]

i send you a "small" version of my script. step 1 took around 8 minutes to prepare the sample.
the final calc. is unfortunately running for ca. 2 hours. i hope this is ok for testing ...

what i see is, that for small increase or decrease of porosity, the pressure is ok (like in the beginning of the calc.):

p: 0.0785083522065 /poro: 48.1453194129
p: 0.0861023388763 /poro: 48.1453043476
p: 0.0870515283585 /poro: 48.1452806796
p: 0.124337357905 /poro: 48.1452662981
p: 0.134103911243 /poro: 48.1452555768
p: 0.134944854374 /poro: 48.1452507787
p: 0.13938079395 /poro: 48.1452458346
p: 0.141438612959 /poro: 48.145242923
p: 0.152935808325 /poro: 48.1452409334

...

but when the grain structure collapses, pressure calculation does not work correctly anymore:

p: 11.42792547 /poro: 48.1451353243
p: 11.4313053052 /poro: 48.1451358473
p: 11.4343726311 /poro: 48.1451363719
p: 11.4371363766 /poro: 48.1451368869
p: 11.4392032959 /poro: 48.1451373207
p: 11.4424940715 /poro: 48.145137852
p: 8.11778479291 /poro: 48.1448508288
p: 9.8711045028 /poro: 48.1447914659
p: 10.7528636056 /poro: 48.1450854371
p: 10.5477578602 /poro: 48.1447098899
p: 8.11642447213 /poro: 48.1465865647
p: 6.99041462318 /poro: 48.1486620985
p: 8.99533593814 /poro: 48.1506954255
p: 7.83218662093 /poro: 48.1531095185
p: 8.69203259997 /poro: 48.154684897
p: 11.6359525072 /poro: 48.1548101512
p: 12.2310464532 /poro: 48.1519503998
p: 11.8209039334 /poro: 48.1508466024
p: 10.41266522 /poro: 48.1496083854
p: 13.8113509815 /poro: 48.1501534546
p: 11.7659257219 /poro: 48.1449743123
p: 11.3078736187 /poro: 48.1416428342
p: 11.7379254629 /poro: 48.1358210037
p: 7.75968666335 /poro: 48.1314928616
p: 9.11740236569 /poro: 48.1293477162
p: 9.30441402016 /poro: 48.1290937942
p: 10.6806525624 /poro: 48.1307995448
p: 10.256471075 /poro: 48.1334583699
p: 11.1890787209 /poro: 48.1358109485

...

p: 27.3255613877 /poro: 48.1124053232
p: 18.8932169787 /poro: 48.0921313593
p: 25.0517632071 /poro: 48.0918358101
p: 30.6684972931 /poro: 48.0860760522
p: 24.5449112901 /poro: 48.0819263772
p: 22.9203772051 /poro: 48.081861069
p: 23.4897854262 /poro: 48.0809934653
p: 30.0214400989 /poro: 48.0765479153
p: 28.6346694789 /poro: 48.0705626027
p: 30.475391 /poro: 48.0628890719
p: 27.1475228671 /poro: 48.0583210558
p: 29.2405442094 /poro: 48.0497216434
p: 25.9130183935 /poro: 48.0443911735
p: 20.737201863 /poro: 48.0378116206
p: 27.6899778636 /poro: 48.035147687
p: 27.2237453694 /poro: 48.0277675295
p: 29.6125202476 /poro: 48.0245254687
p: 27.7152319951 /poro: 48.018681696
p: 27.2864857718 /poro: 48.0154804492
p: 24.0510331304 /poro: 48.0148798888
p: 23.9897576147 /poro: 48.0158389365
p: 23.9914264679 /poro: 48.0165658403
p: 26.6901780408 /poro: 48.0155119429
p: 24.9508433898 /poro: 48.0144790611
p: 28.3914458976 /poro: 48.0112270994
p: 26.2457685113 /poro: 48.0033967489
p: 28.4702987162 /poro: 47.9955400221
p: 30.7188956685 /poro: 47.9883818643
p: 28.8575853606 /poro: 47.9814830293

...

p: 106.320892489 /poro: 46.75595998
p: 106.306309536 /poro: 46.7559601475
p: 106.294913322 /poro: 46.7559601918
p: 106.285210065 /poro: 46.7559603014
p: 106.273826458 /poro: 46.7559605049
p: 106.263787039 /poro: 46.7559605891
p: 106.2516...

[Bruno Chareyre (bruno-chareyre)]

Didn't try the script yet, but one problem is that the boundaries of your system are defined using min/max coordinates of the spheres, it does not reflect the movement of the walls at all. Thus the total volume change is wrong.

You need to assign the walls as boundary conditions, just like in the oedometer example:
flow.boundaryUseMaxMin=[0,0,0,0,0,0]

And if the numbering of the walls is not the same as in triaxial scripts (0 to 5), you need to list them explicitely:
flow.wallIds=[x-,x+,y-,y+,z-,z+]

[Christian Jakob (jakob-ifgt)]

> You need to assign the walls as boundary conditions, just like in the oedometer example:
> flow.boundaryUseMaxMin=[0,0,0,0,0,0]

I did so, but this messages occur:

...
Starting calculation, number of steps: 50000000
6 : Vh==NULL!! id=6 Point=0.000512399 0.00124487 0.00125213 rad=3.15073e-05
7 : Vh==NULL!! id=7 Point=0.000587947 0.000954619 0.000983527 rad=0.000124813
8 : Vh==NULL!! id=8 Point=0.00146642 0.000434536 0.00102062 rad=7.286e-05
9 : Vh==NULL!! id=9 Point=0.00151061 0.0011991 0.000915917 rad=9.27561e-05
10 : Vh==NULL!! id=10 Point=0.00160129 0.00132019 0.00131859 rad=0.000116441
...
279 : Vh==NULL!! id=279 Point=-nan -nan -nan rad=5.56929e-05
280 : Vh==NULL!! id=280 Point=-nan -nan -nan rad=8.59516e-05
281 : Vh==NULL!! id=281 Point=-nan -nan -nan rad=9.2951e-05
282 : Vh==NULL!! id=282 Point=-nan -nan -nan rad=0.000119065
CHOLMOD error: invalid xtype
CHOLMOD error: argument missing
Segmentation fault

hm ...

[Bruno Chareyre (bruno-chareyre)]

> I did so, but this messages occur:
>
> ...
> Starting calculation, number of steps: 50000000
> 6 : Vh==NULL!! id=6 Point=0.000512399 0.00124487 0.00125213 rad=3.15073e-05

It means that the geometry is weird and triangulation can't be
constructed. "nan" positions does not sound good either.
Does it happen at the first iteration after activating the fluid?
It could be the triax engine not working well when combined with fluid
flow, thus moving boundaries the wrong way. I would try to assign the
movement of the walls by myslef to see if it works correctly in that case.
Also, did you inspect the positions of particles and walls? Why are some
coordinates "nan"?
Isn't it possible to reduce the number of particles to make the script
faster. Waiting 2h to see a bug is scary.
Bruno

[Christian Jakob (jakob-ifgt)]

> Does it happen at the first iteration after activating the fluid?

yes

> Isn't it possible to reduce the number of particles to make the script
> faster. Waiting 2h to see a bug is scary.

yes, can be done by reducing dimension (dim vector) in the script.

> Why are some coordinates "nan"?

i hoped that you can tell me why :(

[Bruno Chareyre (bruno-chareyre)]

>i hoped that you can tell me why :(

See #8 again. My guess is TriaxEngine can't work properly because it needs an estimate of the boundary stiffness (usually cumulated stiffness of contacts) and the way it is obtained is not appropriate in your case (if undrained then the fluid stiffness should be also included in the global stiffness).

Try imposing velocities to the boundaries of the undrained sample and see if the results are consistent in terms of pressure changes. If so, you can control loading with your own python function or find a tuning parameter of TriaxEngine to make it more stable.

[Christian Jakob (jakob-ifgt)]

> Try imposing velocities to the boundaries of the undrained sample and see if the
> results are consistent in terms of pressure changes.

I rewrote the script. It is now taking 8 minutes for preparing and 2 minutes for final calc. The script saves the model after preparation. So you can set load = True at the beginning of the script to load prepared sample.

changes:
Triax engine is deactivated after preparation steps and velocity of top boundary is imposed, so that the sample is compressed. In the undrained case I get the following result:

Running script pfv-bug-cell-pressure-too-low.py
set fluid bulk modulus to: 2.000000e+09
Starting calculation, number of steps: 100000
p: -0.0 /poro: 48.147175
p: 77.3254702635 /poro: 47.5619234375
p: 128.265866485 /poro: 46.971240625
p: 141.818936358 /poro: 46.3730234375
p: 190.076801101 /poro: 45.7726046875
p: 256.087922324 /poro: 45.163634375
p: 286.638303292 /poro: 44.5447875
p: 336.686906526 /poro: 43.9271546875
p: 370.445026634 /poro: 43.3041359375
p: -246.733682291 /poro: 42.6770140625

If I uncomment line 145: flow.boundaryUseMaxMin=[0,0,0,0,0,0] triangulation is not working correctly...

Running script pfv-bug-cell-pressure-too-low.py
set fluid bulk modulus to: 2.000000e+09
Starting calculation, number of steps: 100000
6 : Vh==NULL!! id=6 Point=0.000512399 0.00124487 0.00125213 rad=3.15073e-05
7 : Vh==NULL!! id=7 Point=0.000587947 0.000954619 0.000983527 rad=0.000124813
8 : Vh==NULL!! id=8 Point=0.00146642 0.000434536 0.00102062 rad=7.286e-05
9 : Vh==NULL!! id=9 Point=0.00151061 0.0011991 0.000915917 rad=9.27561e-05
...
...
279 : Vh==NULL!! id=279 Point=0.000526466 0.000269253 0.00137925 rad=5.56929e-05
280 : Vh==NULL!! id=280 Point=0.000373563 0.000806535 0.000942971 rad=8.59516e-05
281 : Vh==NULL!! id=281 Point=0.000918819 0.00104285 0.00064448 rad=9.2951e-05
282 : Vh==NULL!! id=282 Point=0.000928043 0.0006173 0.000905637 rad=0.000119065
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
AREA <= 0!!
p: nan /poro: 48.147175
...

It seems, that all cells get zero volume with each triangulation step.

c

[Bruno Chareyre (bruno-chareyre)]

> If I uncomment line 145: flow.boundaryUseMaxMin=[0,0,0,0,0,0]
> triangulation is not working correctly...

This is the real problem. When using max/min the pressure is weird but
it is not a bug, it is expected behavior. Basically you can use that
only when the boundaries are fixed or they have imposed pressure.

Thanks for the simplified script. I've seen the problem. You define
walls with thickness 0 while the default in FlowEngine is 0.001.
Thus the positions of flow boundaries are wrong and it breaks triangulation.
Assigning the correct value solves the problem:
flow.wallThickness=0

Clearly, having 0.001 thickness by default is stupid. It should be zero,
this is the bug. I will fix that.

Bruno

[Christian Jakob (jakob-ifgt)]

Thanks a lot, Bruno!

[Bruno Chareyre (bruno-chareyre)]

Fixed
https://bugs.launchpad.net/yade/+bug/1362090

[Christian Jakob (jakob-ifgt)]

fixed in

https://github.com/yade/trunk/commit/5c4632b49976ea59c870daf734486db4f401060f