yade-users team mailing list archive

[Hossein <question699838@xxxxxxxxxxxxxxxxxxxxx>]

New question #699838 on Yade:
https://answers.launchpad.net/yade/+question/699838

Hi everybody 
I'm trying to achieve the SWRC curve base on the experiment SWRC. first of all, The  soil is fine-grained and I encountered an error 
(rmin>rmax) and then I changed the type of pressure method entering and this error solved. Since, the soil is fine-grained a large 
number of suction need  (around 900 kPa ) during the simulation while, the water content (omega) should decrease slightly. My
real problem is that after 50 kPa in suction rate, water content decreased sharply and become almost zero (0.001) however, after 
900 kPa must become 0.001.
I mention my code in below and I guess the main problem is in terms of pore network which is not well molded since the PSD couldn't match with the real PSD due to the very small particles.

###################################
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
from yade import pack
import pylab
from numpy import *
import numpy as np
from math import *

utils.readParamsFromTable(seed=1,num_spheres=20000,compFricDegree =35.0)
from yade.params import table

seed=table.seed
num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
confiningS=-1e5


## creat a packing with a specific particle side distribution (PSD)
psdSizes,psdCumm=[0.000016,0.000025,0.000036,0.000074,0.00014,0.00017,0.00024,0.00029,0.00043,0.00060,0.00083,0.0014,0.002],[0.081,0.089,0.094,0.12,0.20,.23,.29,.38,.42,.57,.71,.88,1]

sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.01,0.01,0.01)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed)
sp.psd(bins=17,mass=True)

#key='_etaRoll_'+str(etaRoll)+'_rate_'+str(rate/2)


## create material #0, which will be used as default
O.materials.append(CohFrictMat(young=2e8,poisson=.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=2e8,poisson=.3,frictionAngle=0,density=0,label='frictionless'))

## create walls around the packing
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

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

triax=TriaxialStressController(
 internalCompaction=False,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)
newton=NewtonIntegrator(damping=0.4)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
		#Finally, two different contact laws for sphere-box and sphere-sphere
		[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
			useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
			always_use_moment_law=False,  #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
			label='cohesiveLaw')],
  ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'),
 triax,
 newton
]

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  print ('unb=', unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1),'e=', (triax.porosity)/(1-triax.porosity), 'Size=', (triax.width),)
  if unb<0.1 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.05 and (triax.porosity)/(1-triax.porosity)<0.80999999999 :
    break

print ("Particle.Distribution.is.stable")
#############################
##   REACH NEW EQU. STATE ###
#############################
finalFricDegree = 30 # contact friction during the deviatoric loading

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False
# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.1 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.05:
    break

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
#O.save('1kPacking.yade') #save the packing, which can be reloaded later.

O.run(1000,True)
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]

from yade import plot
O.engines=O.engines+[PyRunner(iterPeriod=20,command='history()',dead=1,label='recorder')]

def history():
	plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
		s11=-triax.stress(0)[0]-si0,
		s22=-triax.stress(2)[1]-si1,
		s33=-triax.stress(4)[2]-si2,
		pc=-unsat.bndCondValue[2],
		sw=unsat.getSaturation(False),
		
		i=O.iter
		)

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


#######################################################
##     Drainage Test under Triaxial conditions     ###
#######################################################
##oedometer conditions
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0
triax.wall_left_activated=0
triax.wall_right_activated=0
recorder.dead=0



##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.
unsat.entryPressureMethod=2
print ('meanDiameter=',meanDiameter)

##set boundary conditions, the drainage is controlled by decreasing W-phase pressure and keeping NW-phase pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
unsat.initialization()
#unsat.keepTriangulation=True
unsat.surfaceTension =0.0728


PgEnd=10

print ('e=', triax.porosity/(1-triax.porosity))

##start invasion, the data of normalized pc-sw-strain will be written into pcSwStrain.txt
file=open('pcSwStrain'+'.txt',"w")
for pg in arange(1.0e-5,PgEnd,0.1):
  #increase gaz pressure at the top boundary
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  #compute the evolution of interfaces
  unsat.invasion()
  #save the phases distribution in vtk format, to be displayed by paraview
  unsat.savePhaseVtk("./Phase1")
  unsat.saveVtk("./Pressure1")   #Save pressure field in vtk format. Specify a folder name for output. The cells adjacent to the bounding spheres are generated conditionally based on withBoundaries (not compatible with periodic boundaries)
  
  Watercontent=((unsat.getSaturation(False)*(triax.porosity/(1-triax.porosity))/270)
  waterContatnMatrix.append(Watercontent)
  suctionMatrix.append(-unsat.bndCondValue[2])
 
 #timing.reset()
scripts use 0.5 or even 0.1*PWaveTimeStep()
#ts.active=False
  #compute and apply the capillary forces on each particle
  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
  #reac
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    print ('e=', triax.porosity/(1-triax.porosity), 'suction=' , (-unsat.bndCondValue[2]*1e-3))
    if unb<0.1:
     break
     

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