Yade

Poor pore pressure transmission in FlowEngine

Asked by Chu

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:
Last reply:
Revision history for this message
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

Revision history for this message
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
>
>
>

Revision history for this message
Chu (chu-1) 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

Revision history for this message
Chu (chu-1) 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.

Revision history for this message
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

Revision history for this message
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?

Revision history for this message
Chu (chu-1) 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

Revision history for this message
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)

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#10

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

Revision history for this message
Chu (chu-1) 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

Revision history for this message
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
>
>
>

Revision history for this message
Chu (chu-1) 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

Revision history for this message
Launchpad Janitor (janitor) said : 		#14

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

Revision history for this message
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.