yade-users team mailing list archive

[Luis Barbosa <question700119@xxxxxxxxxxxxxxxxxxxxx>]

Question #700119 on Yade changed:
https://answers.launchpad.net/yade/+question/700119

Luis Barbosa gave more information on the question:
Hi all. Just adding another option I found. But it raises other
questions.

Starting the simulation with imbibition till it reach the desired
saturation and then start drainage.

But, in this case I am facing issues.

1. From getCellInvVoidVolume it seems that the void volume increases over time, but I would expect it to decrease due to increment of water in the cells.
2. From getCellThresholdSaturation and getCellMergedID I getting zero value.
3. As a option to (2) I used getCellSaturation and calculated the saturation of the sample. It worked, but it increases linearly with pressure in the W-phase. I would expect something non-linear because in my simulations I have different pore sizes.

Here is a very simple script:

#!/usr/bin/python
# -*- encoding=utf-8 -*-
#*************************************************************************

from yade import pack
from yade import bodiesHandling
from yade import export
from yade import utils
from yade import ymport

############################################
###   DEFINING VARIABLES AND MATERIALS   ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
    num_spheres=3000,# number of spheres
    compFricDegree = 1, # contact friction during the confining phase
    key='_triax_base_', # put you simulation's name here
    unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key

targetPorosity = 0.39 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 1 # contact friction during the deviatoric loading
rate=0 # loading rate (strain rate)
damp=0.1 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
#2e4+70e4medio 1e4+70e4bom 1e4+60e4bom 3e4+90e4+w3,1,-1-the best
young=10e5 # contact stiffness200e4
young2=10e5
youngcoat=90e4
bondstr=10e8#2e7
bondstr2=0.8e5
bondstrcoat=5e4
mn,mx=Vector3(0,0,0),Vector3(0.006,0.006,0.006) # corners of the matrix
mncoat,mxcoat=Vector3(0,0,0),Vector3(0.006,0.0005,0.006)# corners of the coating
mnbox,mxbox=Vector3(0,0,0),Vector3(0.006,0.006,0.006)

## create materials for spheres and plates
mat=O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,tensileStrength=bondstr,cohesion=bondstr,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres'))
O.materials.append(JCFpmMat(type=1,young=20e7,poisson=0.3,frictionAngle=radians(0),density=2600,tensileStrength=0,cohesion=0,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=0,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='walls'))
O.materials.append(JCFpmMat(type=1,young=youngcoat,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,tensileStrength=bondstrcoat,cohesion=bondstrcoat,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstrcoat,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spherescoat'))

## create walls around the packing
walls=aabbWalls([mnbox,mxbox],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp1=pack.SpherePack()
#sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=[0.0002,0.0012,0.001205], psdCumm=[0.4,0.6,1.],num=500,distributeMass=True,seed=1) #"seed" make the "random" generation always the same
sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=0.0006,num=num_spheres,seed=1)

O.bodies.append([sphere(center,rad,material='spheres') for center,rad in
sp])


############################
###   DEFINING ENGINES   ###
############################

triax=TriaxialStressController(
    maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
    thickness = 0,
    stressMask = 0,
    internalCompaction=False, # If true the confining pressure is generated by growing particles
    wall_front_activated=True,
    wall_back_activated=True,
    wall_top_activated=True,
    wall_bottom_activated=True,
    wall_left_activated=True,
    wall_right_activated=True,
    goal1=-20,
    goal2=-20,
    goal3=-20,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[

	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=-1,xSectionWeibullShapeParameter=1.5, xSectionWeibullScaleParameter=2,label='Physicspheres')],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=False)]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
	triax,
	newton

]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

###################################################
###   REACHING A SPECIFIED POROSITY PRECISELY   ###
###################################################

import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
    # we decrease friction value and apply it to all the bodies and contacts
#    compFricDegree = 0.95*compFricDegree
    setContactFriction(radians(compFricDegree))
    print ("\r Friction: ",compFricDegree," porosity:",triax.porosity,
    sys.stdout.flush())
    # while we run steps, triax will tend to grow particles as the packing
    # keeps shrinking as a consequence of decreasing friction. Consequently
    # porosity will decrease
    O.run(500,1)

print ("###    Compacted state saved      ###")
print(triax.stress(3)[1])

##############################
###   DEVIATORIC LOADING   ###
##############################

setContactFriction(radians(finalFricDegree))


#==========Change Cohesion==============
for b in O.bodies:
	b.mat.cohesion=bondstr2
	b.mat.young=young2

for i in O.interactions:
	i.phys.cohesion=bondstr2
	i.phys.young=young2
#=======================================

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 2
triax.wall_bottom_activated=0
#now goal2 is the target strain rate
triax.goal1=rate#triax.stress(1)[0]
triax.goal3=rate#triax.stress(5)[2]
triax.goal2=triax.stress(3)[1]

#####################################################
###    Example of how to record and plot data     ###
#####################################################

#from yade import plot
from yade import plot
O.run(1000,True)

#strain is logarithmic strain or true strain which is ls=(ln1+e) where e=dl/L (strain) 
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]

## a function saving variables

def history():
    plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
            ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
            s11=-triax.stress(triax.wall_right_id)[0]-si0,
            s22=-triax.stress(triax.wall_top_id)[1]-si1,
            s33=-triax.stress(triax.wall_front_id)[2]-si2,
            pc=unsat.bndCondValue[2],
            sw=vvtsat,#unsat.getCellThresholdSaturation(),
            z1=O.bodies[3].state.pos[1],
            i=O.iter)

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

plot.plots={'pc':('sw')}
plot.plot()

#######################################################
##     Drainage Test under oedometer conditions     ###
#######################################################
##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()

unsat.bndCondIsWaterReservoir=[0,0,1,0,0,0]
unsat.bndCondIsPressure=[0,0,1,0,0,0]
unsat.iniVoidVolumes=True
unsat.surfaceTension = 0.0728
unsat.drainageFirst=False
unsat.isDrainageActivated=False
unsat.isImbibitionActivated=True
unsat.initialization()

voidvol=0.0
voidvoltot=0.0
vvt=0.0
nvoids=unsat.nCells()
initialvol=[0.0] * (nvoids)
vv=[0.0] * (nvoids)

#==================================
ts=O.dt
pgstep= 4500000000*ts #30Pa/s
print (pgstep)
pgmax= 600000#9316 #Pa
print ()
for pg in arange(1,pgmax,pgstep):
  unsat.bndCondValue=[0,0,pg,0,0,0]
  print(unsat.getCellMergedID(1))

  for ii in range(nvoids):
    initialvol[ii]=1./unsat.getCellInvVoidVolume(ii)#
    voidvoltot+=initialvol[ii]
    vv[ii]=unsat.getCellSaturation(ii)
    vvt+=vv[ii]
    vvtsat=vvt/nvoids
  unsat.computeCapillaryForce()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.