# yade-users team mailing list archive

## [Question #699838]: SWRC curve for fine-grained soil

```New question #699838 on Yade:

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
import pylab
from numpy import *
import numpy as np
from math import *

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(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)

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),
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 ###
#############################

#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)

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.run(1000,True)
ei0=-triax.strain;ei1=-triax.strain;ei2=-triax.strain
si0=-triax.stress(0);si1=-triax.stress(2);si2=-triax.stress(4)

def history():
s11=-triax.stress(0)-si0,
s22=-triax.stress(2)-si1,
s33=-triax.stress(4)-si2,
pc=-unsat.bndCondValue,
sw=unsat.getSaturation(False),

i=O.iter
)

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

#######################################################
##     Drainage Test under Triaxial conditions     ###
#######################################################
##oedometer conditions
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)
triax.goal2=goalTop
triax.wall_bottom_activated=0
triax.wall_left_activated=0
triax.wall_right_activated=0

##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
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)

#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*1e-3))
if unb<0.1:
break

--