negative volume for an ordinary pore

Asked by Zhang Wenyue

Hi,

I'm working on DEM-PFV coupling and get some problems.

1)When I run the simulation, I get such warning:
negative volume for an ordinary pore (temp warning, should still be safe)

2)I start the simulation through the controller and inspect it in the 3D view, I find that the small particles are not moving under the fluid force.
The FlowEngine is added before the newton integrator, why doesn't the fluid force acting on the particles?

Can anyone give me some suggestion? Thanks a lot.
I'm running Yade 2018.02b, Ubuntu 18.04.

My code is shown as follow:

====================================================
from yade import pack

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(2,2,6) # 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)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.3,rRelFuzz=0,num=100,periodic=False)
sp.makeCloud(mn,mx,rMean=0.03,rRelFuzz=0,num=900,periodic=False)
sp.toSimulation(material='spheres')

triax=TriaxialStressController(
 thickness = 0,
 stressMask = 7,
 max_vel = 0.005,
 internalCompaction=0, # If true the confining pressure is generated by growing particles
    wall_bottom_activated=False
)

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 _________________

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

#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,0,0,1,1]
flow.bndCondValue=[0,0,0,0,10,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
O.dt=0.1e-3
O.dynDt=False

newton.damping=0
======================================================

Best,
Zhang Wenyue

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
Last reply:

This question was reopened

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

Hi, could you describe in which way you modified the example script to make it disfunctional?

Revision history for this message
Zhang Wenyue (zhangwenyue) said :
#2

Sorry, the script should be like this.
I find that the small particles did move slowly, which means the engine works. But the warning of negative volume for an ordinary pore still appear.

============================================

from yade import pack

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(2,2,6) # 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)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.3,rRelFuzz=0,num=50,periodic=False)
sp.makeCloud(mn,mx,rMean=0.03,rRelFuzz=0,num=900,periodic=False)
sp.toSimulation(material='spheres')

triax=TriaxialStressController(
 thickness = 0,
 stressMask = 7,
 max_vel = 0.005,
 internalCompaction=0, # If true the confining pressure is generated by growing particles
    wall_bottom_activated=False
)

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 _________________

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

#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,0,0,1,1]
flow.bndCondValue=[0,0,0,0,100,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
O.dt=0.1e-3
O.dynDt=False

newton.damping=0

Revision history for this message
Zhang Wenyue (zhangwenyue) said :
#3

I have noticed a similar question:
https://answers.launchpad.net/yade/+question/634175
I don't understand the answer in that question. Can anyone give me more help?

I also got the following error:
-----------------------------------------
......
1048 : Vh==NULL!! id=1048 Point=21.9765 -0.29241 -8.71654 rad=0.03
1049 : Vh==NULL!! id=1049 Point=-nan -nan -nan rad=0.03
1050 : Vh==NULL!! id=1050 Point=17.5632 13.8312 -3.46853 rad=0.03
1051 : Vh==NULL!! id=1051 Point=47969.9 27856.3 55384.9 rad=0.03
1052 : Vh==NULL!! id=1052 Point=1.07385 -0.399759 -1.86005 rad=0.03
1053 : Vh==NULL!! id=1053 Point=-nan -nan -nan rad=0.03
1054 : Vh==NULL!! id=1054 Point=-136772 298537 61692.5 rad=0.03
1055 : Vh==NULL!! id=1055 Point=-nan -nan -nan rad=0.03
Vh==NULL!! id=1054 Point=-13046.5 28479.2 5889.76 rad=0.03
CHOLMOD error: invalid xtype. file: ../Cholesky/cholmod_analyze.c line: 431
CHOLMOD error: argument missing. file: ../Cholesky/cholmod_factorize.c line: 121
Segmentation fault (core dumped)
--------------------------------------------

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

Hi, could you describe in which way you modified the example script to make it disfunctional?

Revision history for this message
Robert Caulk (rcaulk) said :
#5

Hello,

You have two problems that require two solutions. One of your problems is an interesting phenomenon that needs to be documented, so I went ahead and answered Bruno's question for you ;-)

<< First problem >> negative volume for ordinary pore:

Example reply to Bruno: I notice that when I add the small particles to the script, I get this error. But when I only use the large particles, the script is functional.

Solution: Reduce the time step. Smaller particles require smaller time steps. Try an order of magnitude less:

O.dt=0.1e-4

<< Second problem >> Vh==NULL!!:

Example reply to Bruno: When I use stressMask=7 (triaxialStress goals) with the Z-direction pressure gradient (flow.bndCondValue=[0,0,0,0,100,0], I notice the packing drifts in the Z-direction until the particles become uncontained and fly away, resulting in the Vh==NULL!! error.

Explanation: Your imposed pressure gradient is pushing on the particles in the Z-direction, those particles are pushing on your walls in the Z-direction. You commanded your walls to maintain constant stress, so they move to accommodate the changing stress.

Solution: The problem can be corrected by locking one of the walls [2] associated with the pressure boundary condition:

triax.wall_back_activated=False

Cheers,

Robert

[1]https://answers.launchpad.net/yade/+question/634175
[2]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TriaxialStressController.wall_back_activated

Revision history for this message
Zhang Wenyue (zhangwenyue) said :
#6

Thanks Robert Caulk, that solved my question.

Revision history for this message
Zhang Wenyue (zhangwenyue) said :
#7

Dear Robert,

I try O.dt=0.1e-4. At first it performs well, but when I run long enough, the warning of negative volume still appears. I wonder whether there is a method which can adjust dt automatically. This time the program didn't crash. Can I just ignore the warning? Will it have any negative effect on the simulation? Or is there any method can adjust the timestep automatically so that the warning won't show up?

Thanks

Revision history for this message
Best Robert Caulk (rcaulk) said :
#8

>> I wonder whether there is a method which can adjust dt automatically.

Probably! If you determine such a method I encourage you to share with the community.

>> Can I just ignore the warning? Will it have any negative effect on the simulation?

Here is the origin of the warning [1], looks like the author is using a triple product to evaluate each cell's volume [2] so the pseudoscalar output needs to be signed [3]. In other words, the volume calc may yield a negative value despite being perfectly safe. Hence the "should still be safe". However, in your case you are likely seeing this as a result of the large particle size mismatch in your simulation and the fast-moving particles. I think this combination results in very acute tetrahedral angles and if these particles move faster than the triangulation updates, they end up easily inverting themselves, resulting in a negative triple product. You are still getting the correct volume of a pore represented by the locations of those particles, but the forces are probably pretty inaccurate since the triangulation is outdated. If you are scared of this you could try updating the triangulation more frequently with meshUpdateInterval.

[1]https://github.com/yade/trunk/blob/master/pkg/pfv/FlowEngine.ipp.in#L629
[2]https://github.com/yade/trunk/blob/master/pkg/pfv/FlowEngine.ipp.in#L628
[3]https://github.com/yade/trunk/blob/master/pkg/pfv/FlowEngine.ipp.in#L630

Revision history for this message
Zhang Wenyue (zhangwenyue) said :
#9

Thanks Robert Caulk, that solved my question.