PFV problems with sedimentation in fluid

And some more problems are here:
(1) Can I use facetBox to build box other than aabbWalls([mn,mx],thickness=0,material='walls')\
(2) What does the triax.stressMask means from 0 to 7, I am not sure about it.
(3) If I want the water level can go down gradually after its leakage or go out form the bottom, is it possible to simulate with PFV?
(4) How can I simulate the problems like "4 Immersed granular deposition" [Catalano2014a] Catalano, E., Chareyre, B., Barthélémy, E. (2014), Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects.
Thank you!

hi,

just a short answer to question (2):

stressMask is not good explained in documentation, have a look at examples/triax-tutorial for a better documentation

https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py#L87

hope it helps,

christian

Hi,

I think it is an error in the boundary conditions. Check probably here: flow.bndCondIsPressure

Bests
Donia

Hi, Jakob, Donia.
Thank you!

It is true, the boundary conditions seems not working when flow.bndCondIsPressure is set as zero flux boundary in six sides.
And flow.bndCondIsPressure=[0,0,1,1,0,0] works.

But the FlowEngine has more things I am not sure:
(1) The default viscosity of fluid is 1, not 0.001Pa.s for water, why? Is the default Values of unit not International unit ?
(2) There is no setting for fluid density, What the default value of it?
(3) There is no buoyant force, when a sedimentation simulation, how to calculate it for different spheres.

>(1) The default viscosity of fluid is 1, not 0.001Pa.s for water, why?

Some fluids are not water.

>(2) There is no setting for fluid density, What the default value of it?

Tthe current model is assuming Stokes flow, which is non-inertial flow.
Yet, you are right that buoyancy depends on fluid density ....

>(3) There is no buoyant force, when a sedimentation simulation

... so, you should define the material like this:
material.density = (solidDensity-fluidDensity) #submeged density, independent of fluid flow.

Hi Bruno,
Thank you!

>>>(1) The default viscosity of fluid is 1, not 0.001Pa.s for water, why?
>>Some fluids are not water.

So if the viscosity of fluid is 0.001Pa.s , I can set the viscosity value to 0.001, is it true?

And I have test the problems similar to the sedimentation in "4 Immersed granular deposition" [Catalano2014a] Catalano, E., Chareyre, B., Barthélémy, E. (2014), Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects.

It works for some time, but it breaks down 35s , could you please give some advice for the following scripts, thank you!
## ---------------------------------start

# basic simulation showing sphere falling ball gravity,
# bouncing against another sphere representing the support
from yade import pack
# DATA COMPONENTS
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,1.5,0.15) # 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,0.03,0.333,num_spheres,False, 0.61,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 max_vel = 0.005,
 internalCompaction=0, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=0.1,gravity=(0,-9.8,0))

# simulation loop -- see presentation for the explanation
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.wall_bottom_activated=False
triax.wall_top_activated=False
triax.wall_left_activated=False
triax.wall_right_activated=False
triax.wall_back_activated=False
triax.wall_front_activated=False

# Blocked certain degress of freedom to make 2D-Model in plane-XY
##for k in O.bodies:
## if isinstance(k.shape, Sphere): k.state.blockedDOFs='zXY'

# set timestep to a fraction of the critical timestep
# the fraction is very small, so that the simulation is not too fast
# and the motion can be observed
O.dt=.5e-4*utils.PWaveTimeStep()

O.run(1,1)
print "###------ ###"

# save the simulation, so that it can be reloaded later, for experimentation
O.saveTmp()

###---------------------------------------------------------------------------------------------------add by lc

####A. Check bulk modulus of the dry material from load/unload cycles
###triax.stressMask=2
###triax.goal1=triax.goal3=0

###triax.internalCompaction=False
###triax.wall_bottom_activated=False
####load
###triax.goal2=11000; O.run(2000,1)
####unload
###triax.goal2=10000; O.run(2000,1)
####load
###triax.goal2=11000; O.run(2000,1)
###e22=triax.strain[1]
####unload
###triax.goal2=10000; O.run(2000,1)

###e22=e22-triax.strain[1]
###modulus = 1000./abs(e22)

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

O.run(1,1)

Qin = flow.getBoundaryFlux(2)
Qout = flow.getBoundaryFlux(3)
permeability = abs(Qin)/O.dt
print "Qin=",Qin," Qout=",Qout," permeability=",permeability

#newton.damping=0

zeroTime=O.time

from yade import plot

## a function saving variables
def history():
 plot.addData(P1 = flow.averageSlicePressure(0.000),P2 = flow.averageSlicePressure(0.3),P3 = flow.averageSlicePressure(0.6),P4 = flow.averageSlicePressure(0.9),P5 = flow.averageSlicePressure(1.2),P6 = flow.averageSlicePressure(1.5),t=O.time)
   #plot.addData(Qin = flow.getBoundaryFlux(5),Qout = flow.getBoundaryFlux(4),t=O.time,permeability = abs(flow.getBoundaryFlux(5))/2/0.75/O.dt)
   #plot.addData(e22=triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
##make nice animations:
#O.engines=O.engines+[PyRunner(iterPeriod=200,command='flow.saveVtk()')]

from yade import plot
plot.plots={'t':('P1','P2','P3','P4','P5','P6')}
plot.plot()
O.saveTmp()
O.timingEnabled=1
from yade import timing
print "starting oedometer simulation"
O.run(200,1)
timing.stats()

## Make more steps to see the convergence to the stationnary solution

##----------------------------------------------------------------------------------------------------------------------------------------------end

> So if the viscosity of fluid is 0.001Pa.s , I can set the viscosity value to 0.001, is it true?

If you define pressure in Pa and time in s, yes. Otherwise, you need to be consistent.
I'm afraid I will not have time to test your script in the near future. Donia (above) is already waiting for my help for another PFV problem ;)

Something I missed:
> It is true, the boundary conditions seems not working when flow.bndCondIsPressure is set as zero flux boundary in six sides.

Correct. Incompressible fluid flow with only Neuman-type (no-flux) BCs is ill-posed. It is enough to impose an arbitrary value of pressure in one point though. No need to change the boundary condition.

Hi Bruno,
Ok, thank you.

Do you have any ideas about the following message, what do they mean? What are the most possibility for "CHOLMOD warning: matrix not positive definite"? And what are LDLt?
--------------
32 : Vh==NULL!! id=32 Point=0.641734 -0.368851 0.125888 rad=0.0394611
59 : Vh==NULL!! id=59 Point=0.785606 -0.130624 0.0825895 rad=0.0389221
......
1005 : Vh==NULL!! id=1005 Point=1.41177e+07 -3.31642e+07 -1.11211e+07 rad=0.0200399
infinite K1!
infinite K1!
CHOLMOD warning: matrix not positive definite
something went wrong in Cholesky factorization, use LDLt as fallback this time
AREA <= 0!!
AREA <= 0!!
9 : Vh==NULL!! id=9 Point=1.19656e+09 -1.65309e+10 6.17722e+09 rad=0.0399201
12 : Vh==NULL!! id=12 Point=3.15611e+10 -2.81586e+10 -9.68463e+09 rad=0.0398603
......
Floating point exception (core dumped)

>>>Correct. Incompressible fluid flow with only Neuman-type (no-flux) BCs is ill-posed. It is enough to impose an arbitrary value of pressure in one point though. No need to change the boundary condition.

Do you mean keep the flow.bndCondIsPressure=[0,0,0,0,0,0], flow.bndCondValue=[0,0,0,0,0,0], but how to impose an arbitrary value of pressure in one point?

Le 28/04/2014 15:11, liucheng83 a écrit :
> Question #247705 on Yade changed:
> https://answers.launchpad.net/yade/+question/247705
>
> liucheng83 posted a new comment:
>>>> Correct. Incompressible fluid flow with only Neuman-type (no-flux)
> BCs is ill-posed. It is enough to impose an arbitrary value of pressure
> in one point though. No need to change the boundary condition.
>
> Do you mean keep the flow.bndCondIsPressure=[0,0,0,1,0,0],
> flow.bndCondValue=[0,0,0,0,0,0], but how to impose an arbitrary value of
> pressure in one point?
>
Hi,

You can impose pressure to a point by:

flow.imposePressure((xpos,ypos,zpos.),value)

Cheers,
Timos

Hi Timos,

You are right, I have tested it. Thank you!

Cheng

Hi,

How to use the flow.pressureProfile(), it cannot be ploted or exported.

flow.pressureProfile(wallUpY, wallDownY) returns a file in which you have the values of pressure at different Y positions as described and calculated by the function: measurePressureProfile(wallUpY,wallDownY)

Bests
Donia

Hi,

The problem is that the balls went out of the wall when the sedimentation reach the end, what may be the reason?

Hi Donia,

I have found the file named Pressure_profile, and it really generated a better pic, thank you!

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