Yade

Poor pore pressure transmission in FlowEngine

Asked by Chu on 2019-06-19

Hello,

I use FlowEngine to simulate Fluid-Structure Interaction.I generate a dense sample,but I get a very strange pore pressure field.[1]When I set small Viscosity or large CondValue,I get correct pore pressure field.[2]But I must use small CondValue.Could you give me some suggestions?

Thanks for your any suggestion.

[1]https://s2.ax1x.com/2019/06/19/VXZA7F.png
[2]https://s2.ax1x.com/2019/06/19/VXm8Wd.png

#############MWE###################
from yade import pack

psdSizes,psdCumm=[0.000075,0.00015,0.00118,0.00236],[0,0.35,0.35,1.]
num_spheres=50000# number of spheres
young=1e8
compFricDegree = 0 # initial contact friction during the confining phase
finalFricDegree = 30# contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(0.01,0.01,0.01) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.2,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=1e10,poisson=0.2,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,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

triax=TriaxialStressController(
 thickness = 0,
 stressMask = 7,
 internalCompaction=False, # 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
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

triax.goal1=triax.goal2=triax.goal3=-200000

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

setContactFriction(radians(finalFricDegree))
newton=NewtonIntegrator(damping=0)

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
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 newton
]

flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=0.001
flow.bndCondIsPressure=[1,1,0,0,0,0]
flow.bndCondValue=[0,10,0,0,0,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
flow.updateTriangulation=True

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-07-01
Last reply:
2019-07-16
Robert Caulk (rcaulk) said : #1

Hello,

Have you tried decreasing the time step and or running the simulation for a longer duration? Decreasing permeability requires smaller time step to maintain numerical stability and increases the time required to reach the final pressure field.

Robert

Chareyre (bruno-chareyre-9) said : #2

Hi, maybe the first picture has pressure dominated by pressure feedback
from grains motion. What happens if you block positions and set null
velocities?
Bruno

Le lun. 24 juin 2019 14:22, Robert Caulk <
<email address hidden>> a écrit :

> Question #681486 on Yade changed:
> https://answers.launchpad.net/yade/+question/681486
>
> Status: Open => Answered
>
> Robert Caulk proposed the following answer:
> Hello,
>
> Have you tried decreasing the time step and or running the simulation
> for a longer duration? Decreasing permeability requires smaller time
> step to maintain numerical stability and increases the time required to
> reach the final pressure field.
>
> Robert
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
>
>

Chu (arcoubuntu) said : #3

Hi Robert,
>Have you tried decreasing the time step and or running the simulation for a longer duration?
 Yes,I tried.Decreasing the time step seems no effect ,but running the simulation for a longer duration (30k iters) looks better[1].Is it means that the FlowEngine works correctly?

In addition,I tried using flow.useSlover=0,and I got correct pore pressure field.But it caculate very slowly,and GS did not converge in 20k iterations.
[1]https://s2.ax1x.com/2019/06/25/ZVRuon.png

Chu (arcoubuntu) said : #4

Hi Bruno,
>What happens if you block positions and set null velocities?
 I added the following code,and it seems no effect.

###########################
for i in O.bodies:
    if isinstance(i.shape,Sphere):
        i.state.blockedDOFs='xyzXYZ'
        i.vel=(0,0,0);i.angVel=(0,0,0)
###########################

As I said in #3,it looks better to run the simulation for a longer duration.Is it means that the FlowEngine works correctly?Whether it will cause the calculation of fluid force is incorrect?
In addition,I tried using flow.useSlover=0,and I got correct pore pressure field.But it caculate very slowly,and GS did not converge in 20k iterations.

Bruno Chareyre (bruno-chareyre) said : #5

> Does it mean that the FlowEngine works correctly?

Until now it's unclear to me why you think something is wrong. So I would assume everything is correct.
Bruno

Robert Caulk (rcaulk) said : #6

>>In addition,I tried using flow.useSlover=0,and I got correct pore pressure field.But it caculate very slowly,and GS did not converge in 20k iterations.

To clarify: you run the same simulation for the same number of time steps using useSolver=0 and useSolver=3 (or 4), and you get entirely different pressure fields?

Chu (arcoubuntu) said : #8

Thanks Bruno and Robert.

I apologize for my mistake in #7.

>Until now it's unclear to me why you think something is wrong.

I consider the pore pressure field should be gradient distribution like [1],but actually it looks irregular[2].

>To clarify: you run the same simulation for the same number of time steps using useSolver=0 and useSolver=3 (or 4), and you get entirely different pressure fields?

I made a mistake here.Actually I find that I run the simulation for the same number of time steps and I got same results.But for one step,I got the pore pressure field[1].It looks well.But for two steps,I got the pore pressure field[2].They're entirely different.I haven't tried to run the simulation for more steps.Why they're so different?

[1]https://s2.ax1x.com/2019/06/27/ZufDNF.png
[2]https://s2.ax1x.com/2019/06/27/Zuf639.png

Bruno Chareyre (bruno-chareyre) said : #9

Please check what happens after two steps while blocking all positions (not just spheres).

for i in O.bodies:
        i.state.blockedDOFs='xyzXYZ'
        i.vel=(0,0,0)

Bruno Chareyre (bruno-chareyre) said : #10

With #2 in mind I still don't see the results as necessarily wrong.

Chu (arcoubuntu) said : #11

Hi Bruno,

I tried what you said,and I get same results.It is worth mentioning that I changed the range of the color map between 0 and 10[1] in Paraview before.Because the boundaryValue is 10.But original range is between -446 and 1168[2].Maybe it is useful.

[1]https://s2.ax1x.com/2019/06/28/ZuLUTs.png
[2]https://s2.ax1x.com/2019/06/28/ZuqoZj.png

Chareyre (bruno-chareyre-9) said : #12

You need to turn triax off to, because it's moving the boundaries.
Basically you superimpose an external (total) stress of 200kPa and a fluid
pressure of 0.01kPa. No surprise if the response is dominated by the former.
Bruno

Le jeu. 27 juin 2019 19:17, Chu <email address hidden> a
écrit :

> Question #681486 on Yade changed:
> https://answers.launchpad.net/yade/+question/681486
>
> Status: Answered => Open
>
> Chu is still having a problem:
> Hi Bruno,
>
> I tried what you said,and I get same results.It is worth mentioning that
> I changed the range of the color map between 0 and 10[1] in Paraview
> before.Because the boundaryValue is 10.But original range is between
> -446 and 1168[2].Maybe it is useful.
>
> [1]https://s2.ax1x.com/2019/06/28/ZuLUTs.png
> [2]https://s2.ax1x.com/2019/06/28/ZuqoZj.png
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
>
>

Chu (arcoubuntu) said : #13

Thanks Bruno,

I added triax.dead=1 to the script,but the pore pressure field seems no better[1].
I also find that if I use equal particle size(all particle size = 0.2mm)[2] or bigger particle size[3](psdSizes,psdCumm=[0.0002,0.0003,0.0015,0.002],[0,0.35,0.35,1.]),without triax.dead=1,the result seems better.

I use very small particle size so that the size of model is very small.So boundary pressure is 0.01kPa which results in a proper hydraulic gradient.

Thanks again for your suggestions.

[1]https://s2.ax1x.com/2019/07/01/Z88Uje.png
[2]https://s2.ax1x.com/2019/07/01/Z88dnH.png
[3]https://s2.ax1x.com/2019/07/01/Z88wBd.png

Launchpad Janitor (janitor) said : #14

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

Bruno Chareyre (bruno-chareyre) said : #15

In my opinion there is something moving in your problem - something you perhaps do not suspect - and the pressure field responds to that. I suggest to freeze all positions, run a few steps, then only after that you start solving the flow (after making sure that all positions are strictly constant across multiple timesteps).

Also I would not recommend to redefine O.engines in the middle of a simulation with same lables, since you may have unexpected behavior (and its really purposeless).

If it still doesn't work please send a working script produciing exactly the vtk fileds you are showing.