# Modeling presence of water( a degree of saturation) in consolidation test

I have been trying to work on the Oedometric consolidation test. So far ( with help from Jan :) ) I have been able to work out the test on a dry granular media . The current trial script is as follows:

import math,time
from math import sin,cos

"""creates cylinder made of boxes. Axis is parallel with z+"""
center = Vector3(center)
angles = [i*2*pi/nSegments for i in range(nSegments)]
ret = []
for i,pt in enumerate(pts):
es = (.5*thick,l,.5*height)
ori = Quaternion((0,0,1),i*2*pi/nSegments)
ret.append(box(pt,es,ori,**kw))
return ret

surf = cylSurf((1,2,3),4,5,nSegments=12,thick=.1,wire=True)
O.bodies.append(surf)
for b in surf: b.state.blockedDOFs='xyzXYZ'
pred=pack.inCylinder((1,2,3),(1,2,8),4)
O.bodies.append(sp)
lower_plateID=O.bodies.append(box((1,2,3),(4,4,0.0002)))
O.bodies[lower_plateID].state.blockedDOFs = 'xyzXYZ'
upper_plate=box((1,2,8),(4,4,0.0002))
#upper_plate.state.mass= 1.5
upper_plateID=O.bodies.append(upper_plate)
O.dt= 0.3*PWaveTimeStep()

Dz=upper_plate.state.displ()[2]
t= O.time

plot.plots= {'Tx':('Dz')}
plot.plot()

I am wiondering if there is a possibility to include and thus model the presence of water ( degree of saturation) to the dry granular media. I have been going through the documentation but so far the only relevant thing I found was the use of flowEngines &. there are no examples for me to look at and understand the application.

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Jérôme Duriez (jduriez) said on 2017-10-05: #1

Hello,

FlowEngine is for saturated cases. For unsaturated i.e. partially saturated cases, you should give a look to

- the Capillary model: see [1] and related classes, and the wiki [2]. This is for pendular regime only (see the wiki)
- or the more general UnsaturatedEngine (probably not included by default in any YADE version), see [3]

 Revision history for this message Swapnil (swapspace) said on 2017-10-05: #2

Okay Jerome, thanks. going through it :) Will come up with my questions soon ;)

 Revision history for this message Swapnil (swapspace) said on 2017-10-10: #3

Jerome,
I have been thinking of first trying the Full Saturation condition of granular media here.

I referred to the Oedometer.py script utilisong the particle-fluid coupling and building up my case from the scenario depicted there.

the following is the script:

num_spheres=1000# number of spheres
young=1e6
compFricDegree = 3 # initial contact friction during the confining phase
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

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,-1,0.3333,num_spheres,False, 0.95,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,
max_vel = 0.005,
internalCompaction=True, # 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"
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
newton
]

triax.goal1=triax.goal2=0
triax.goal3=-10000 ---> Is this correct?

# full saturation condition

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3 -----?????
flow.permeabilityFactor=1
flow.viscosity=1
flow.bndCondIsPressure=[0,0,1,1,0,0] ---??? format meaning
flow.bndCondValue=[0,0,1,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) ---??? Number denote which boundary
Qout = flow.getBoundaryFlux(3)
permeability = abs(Qin)/1.e-4 #size is one, we compute K=V/∇H
print "Qin=",Qin," Qout=",Qout," permeability=",permeability

#C. now the oedometer test, drained at the top, impermeable at the bottom plate
flow.bndCondIsPressure=[0,0,0,1,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
newton.damping=0

I have gone through the documentation but there are not enough examples to catch the theme right.

My current doubts with this script are:

1) How to apply stress only in Z-direction. StressMask needs to hold what value? 2 or 3 ?

2) I did not understand: flow.useSolver=3.

3) In flow.bndCondIsPressure=[0,0,1,1,0,0], flow.bndCondValue=[0,0,1,0,0,0], flow.boundaryUseMaxMin=[0,0,0,0,0,0], how to interpret the sequence 0,0,1,1,,.... Could you give a couple of examples to illustrate.

4) Qin = flow.getBoundaryFlux(2) ---??? Number denote which boundary
Qout = flow.getBoundaryFlux(3)

for example if we are testing 1D consolidation along Z .axis ( as in my case) what should I use?

5) Is it not important to provide the density of the fluid we are using. for example water ..

5) I want to record the change in thickness of sample along the Z-direction. How do I accomplish that. Do I have to use strain 'ezz' or is there another direct way?

 Revision history for this message Jérôme Duriez (jduriez) said on 2017-10-11: #4

Hello,

5) (the 2d one) Yes, I guess the zz strain ezz actually is a good measure of the change in z-length ;-)

For other questions, I'm not familiar with FlowEngine myself, I would advice (as usual) to search through the mailing list in addition to this doc, and to open another (specific) question if this does not solve your questions.

Jérôme

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

2) The short answer is "keep it =3". If =0 you will use a less efficient strategy for solving. Other numbers will not work. Misleading, I agree.

3) The 6 values are for the six boundaries, sorted as in x-,x+,y-,y+,z-,z+. I think the example script oedometer.py illustrates that very well, and it answers directly your sub-question on "1D consolidation along Z .axis".

4) see 2)

5-1) ---> "Is it not important to provide the density of the fluid we are using. for example water .."
No it is not. A viscous flow (assumption of this coupling) depends on viscosity only.
You need to keep in mind that "pressure" in the model corresponds to excess pore pressure, i.e. something in addition to the hydrostatic pressure. The latter depends on density as you suggest but it causes no flow, it is not reflected in FlowEngines's pressure..

5-2) Please check the example script again.

Bruno

 Revision history for this message Swapnil (swapspace) said on 2017-10-12: #6

Thanks Jerome & Bruno :)

My current script is ( for 1D consolidation along z axis):

import math

num_spheres=1000# number of spheres
young=1e6
compFricDegree = 3 # initial contact friction during the confining phase
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

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,-1,0.3333,num_spheres,False, 0.95,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,
max_vel = 0.005,
internalCompaction=True, # 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"
),
triax,
newton
]

triax.goal1=triax.goal2=0
triax.goal3= -10000

# full saturation condition

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3 #-----?????
flow.permeabilityFactor=1
flow.viscosity=10
flow.bndCondIsPressure=[0,0,0,0,0,1]
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(4)
Qout = flow.getBoundaryFlux(5)
permeability=abs(Qin)/0.0001
print "Qin=",Qin," Qout=",Qout," permeability=",permeability

#C. now the oedometer test, drained at the top, impermeable at the bottom plate
flow.bndCondIsPressure=[0,0,0,0,0,1]
flow.bndCondValue=[0,0,0,0,0,0]
newton.damping=0

Dz=-triax.strain[2]
t= O.time

plot.plots= {'Tx':('Dz')}
plot.plot()

I have the following questions now:
1) StressMask : I am still trying to get acquainted with the concept applied there. But I found this useful note in one of the scripts:
### switch stress/strain control using a bitmask. What is a bitmask, huh?!
## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.

This helped me get a better idea ;)

2)

a) For a 1D (z axis flow) we would need: triax.goal1=triax.goal2=0
triax.goal3= -10000 .

Is it correct?

b) For flow drained at the top and bottom undrained it should be:

flow.bndCondIsPressure=[0,0,0,0,0,1]
flow.bndCondValue=[0,0,0,0,0,0]
I think this is how the bit stream works here?

3) Plot:
Is the addPlot data section to measure the change in thickness fine? I am not sure where I am going wrong. All I see is the Dz value constantly 0 with time.

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

I can recognize my prose. Thank you for calling it useful. :)
Any suggestion on how/where it should be explained is welcome, since it is a very frequent question.

2.a) yes
2.b) yes
3) the code looks ok, I did not try the script. Most likely the top/bottom plates are not moving (should they?) or the engine is idle hence not updating the strain values.

Bruno

 Revision history for this message Swapnil (swapspace) said on 2017-10-13: #8

Thanks Bruno :)

1) . I would recommend to have it written along with the explanation of StressMask in the section for Triaxial stress controller. That would be the best spot for people to understand the usage clearly.

I actually had a question in mind but forgot it mid-way. Will write as soon as I remember.

3) In a standard 1D oedometer test, the top plate is loaded and the sample is drained/undrained (from top/bottom).

a) I can see the spheres growing in size. I am not sure if they are doing so with the expulsion of water from the pores ( because I cannot visualize the fluid phase in the simulation. I guess Yade does not include that yet).
Going by the concept of soil mechanics, should the volume of the solid particles not remain constant? Why do the spheres have to change their volume. It should only be the total volume of the sample that would change with the consolidation occuring. Am I getting something wrong :) ?

b) Could you please try running the script. I am unable to figure out the issue with the plot. It isn't updating.

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

> I can see the spheres growing in size.
They grow in size because you request so (internalCompaction=True). In the original script you started with this is turned off after the sample preparation, but you removed it together with "Part A" section of that script, so they keep changing size.

>with the expulsion of water from the pores
No, although it would not be difficult to include.

>should the volume of the solid particles not remain constant?
It should indeed.

> the plot [...] isn't updating.
addPlotData() is not inserted in the engines loop, hence not recording. You need a PyRunner (as found in any plotting example).

Bruno

 Revision history for this message Swapnil (swapspace) said on 2017-10-19: #10

1) So it means that the volume of the spheres remains constant although they appear to grow. clear.

When internal compaction takes place it either results in the expulsion of air (voids) or just acts as a method to apply a certain initial compaction without changing the porosity. Is it?

2) But then, in theory the voids change with the expulsion of water (after the full saturation condition as in this example). How do we include that part in the code?

3) The engine I am using has the Pyrunner with the addPlot data.

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

1) > it either results in the expulsion of air

As I told you it is not expulsing any fluid when growing because the change of size is not reflected in the fluid model (the situation is in fact more complex than that, this is the short answer). Why would it be expulsing "air" now?? This is a one-phase flow model, remember.

> or just acts as a method to apply a certain initial compaction without changing the porosity.

Growing particles is obviously changing the porosity. What do you mean?

2) I am a bit lost. What do you want to include and why?

3) Lost again. In #6 you show a script which is not plotting and you ask why (and you suggest me to test it). That script was simply not recording, hence not plotting. Now you say you actual script is different and is really recording. Why then was I supposed to test the script in #6?

I would suggest:
1/ Please be consistent in your questions, it will save time for us all.
2/ Please study the example scripts seriously, since in this thread and in the one on triaxial test it seems you are faciing problems which are addressed in the examples.

Bruno

 Revision history for this message Swapnil (swapspace) said on 2017-10-20: #12

Okay. sorry again if I caused confusion.

Yes, it appears I missed some segment again like the triaxial test problem in another thread.

thanks for reminding. Will try it out in "full" :)

 Revision history for this message winingercorene (winingercorene) said on 2021-07-26: #19

Well. sorry again if i caused confusion. Yes, it seems that I lost some segment again, like the triax test problem in another thread. You can see my web: https://businessdollars.us/

 Revision history for this message jackob yuosf (jackobyuosf) said on 2021-09-05: #20

For other questions, I am not familiar with FlowEngine, I would recommend (as usual) searching the mailing list in addition to this document, and opening another (specific) question if this does not solve your questions. You can see my blog here: https://justiry.com/