← Back to team overview

yade-users team mailing list archive

[Question #693568]: Stresses drop to zero in periodic triaxial shear simulation

 

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

Hello everyone,

I'm trying to simulate a periodic undrained triaxial test by mantaining constant volume  instead of explicitly simulating water (as a sidenote, I never could get periodicFlowEngine to work correctly for this).  My simulation is divided in three phases, for the first one a cloud of spheres is isotropically compacted to a target porosity, then for the second phase  the packing is isotropically compacted again to a target stress, and in the third one a velGrad is imposed in the cell  so there is compression  a certain rate in one direction and there  is extension at half the compression rate in the other two directions. The idea of this final phase is to impose a constant volume condition as to simulate an undrained behaviour.

The problem is that even though at the end  of the isoCompaction phase the stress is equal to a target stress (100000 for expample) and volume in the shear phase is indeed constant, after a number of  iterations stresses drop to zero and start to oscillate around that value. This is the result I get from my MWE that is at the end of this post:

Welcome to Yade 20201018-4279~61d08ab~bionic1 
Using python version: 3.6.9 (default, Oct  8 2020, 12:12:24) 
[GCC 8.4.0]
TCP python prompt on localhost:9000, auth cookie `acusye'
XMLRPC info provider on http://localhost:21000
qt5ct: using qt5ct plugin
Running script /home/geotesis/DEM/Yade/PeriodicModel_noWater_MWE.py
Start:  2020-10-21 16:09:53.985841
Cloud created, number of particles =  480
Compacting to target porosity:  2020-10-21 16:09:54.003525
[[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]]

Target porosity passed, porosity =  0.49986744040494263
Starting isocompaction simulation:  2020-10-21 16:10:16.938606
Current iteration =  1534900 , s11 =  -100157.84795970531 , s22 =  -100145.52886621605 , s33 =  -100199.74978550954 , e =  0.6462475530919375
Starting shear simulation:  2020-10-21 16:12:28.666616
Current iteration =  1534901 , s11 =  -100154.30036085597 , s22 =  -100148.62230105326 , s33 =  -100200.53394271889 , e =  0.6462475530919375
Current iteration =  1834901 , s11 =  -55393.437466646115 , s22 =  -65305.28482971604 , s33 =  -51735.14392560943 , e =  0.6462749983880726
Current iteration =  2134901 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , e =  0.6462749983315953
Current iteration =  2434901 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , e =  0.6462749981327709
Current iteration =  2734901 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , e =  0.6462749980919655
Current iteration =  3034901 , s11 =  -0.11991207951660474 , s22 =  -0.22419400597402317 , s33 =  -0.026624167711505737 , e =  0.6462749980486365
Current iteration =  3334901 , s11 =  -0.047911199264998226 , s22 =  -0.00023733397986893054 , s33 =  -0.013137608101300035 , e =  0.6462749979832385
Current iteration =  3634901 , s11 =  -0.018439277204848392 , s22 =  -0.026673245415610117 , s33 =  -0.049101414208267816 , e =  0.6462749979320861
Current iteration =  3934901 , s11 =  -0.00045092480566355026 , s22 =  -0.0041975100217730744 , s33 =  -0.07083510893219645 , e =  0.6462749978859719
Current iteration =  4234901 , s11 =  -10.470723417543306 , s22 =  -5.878484222842999 , s33 =  -2.042145920074655 , e =  0.6462749977701003
Current iteration =  4534901 , s11 =  -0.014448912716733427 , s22 =  -0.07253827703395103 , s33 =  -0.013568831616481059 , e =  0.6462749976951889
Current iteration =  4834901 , s11 =  -0.004931408116028016 , s22 =  -0.03483611806022092 , s33 =  -0.0062370493144523085 , e =  0.6462749976127734
Current iteration =  5134901 , s11 =  0.00012374680500377676 , s22 =  -0.0009835150202924818 , s33 =  -0.0049291412638700184 , e =  0.6462749975535781
Current iteration =  5434901 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , e =  0.6462749974492507
Finished simulation phase:  2020-10-21 16:18:16.761407

As you can see, after 600000 iterations stresses in all directions drop from 100000 to zero. When I visually check what's happening in the simulation, It's quite obvious that particles start to separate in this phase and start to float around until the target strains are reached.  I know that there's something I'm not doing or  that I'm doing wrong, but I can't wrap my head around the fact that volume stays constant but particles have space to float around at the same time. This also happens if the isoCompaction goal is higher, though it takes more iterations.

As a summary, my CU TX periodic simulation has zero stresses in shear phase, even though it started from isotropic confinement. Expected behaviour would be that stresses change from initial confinement (s11 and s33 equal, and s22 different), wihout all of them reaching zero, while volume stays constant.

# -*- coding: utf-8 -*-
######################
###   INITIALIZE   ###
######################

from yade import pack
from yade import plot

## Time of simulation
import time
import datetime
print ("Start: ",datetime.datetime.now())
start = time.time()

## These are default parameters for batch execution in full code
confinement=100
e0=1.0
rate=0.01
batch = 0

#######################################
###   BASIC SIMULATION PARAMETERS   ###
#######################################

### Set Periodic cell
O.periodic=True
x1 = 0.006
y1 = x1
z1 = x1
mn,mx = Vector3(0,0,0),Vector3(x1,y1,z1) # corners of the initial packing
O.cell.setBox((x1,y1,z1))

### PSD in [mm]

# Silica Sand #7
psdSizes = [0.075,0.106,0.250,0.425,0.850]
psdCumm = [3,7,88,100,100]

## From [mm] to [m] and scale in homotetic way
psdScale = 3 # 1 is around 10000 particles, 2 around 1500, 3 around 500
psdSizes = [x*(10**-3)*psdScale for x in psdSizes]
psdCumm = [x*(10**-2) for x in psdCumm]

## Material of spheres
young=70e9
compFricDegree = 3 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
O.materials.append(FrictMat(young=young,
                            poisson=0.17,
                            frictionAngle=radians(compFricDegree),
                            density=2600,
                            label="spheres"))

###################
###   ENGINES   ###        
###################

sigmaIso = -1e8
triax = 	PeriTriaxController(
    # specify target values and whether they are strains or stresses
    goal = (sigmaIso,sigmaIso,sigmaIso), stressMask = 7,
    # type of servo-control
    dynCell = True, mass = 0.5, maxStrainRate = (1,1,1),
    # wait until the unbalanced force goes below this value
    maxUnbalanced = 0.01, relStressTol = 0.002,
    doneHook="changePhase()"
)

O.engines=[
    ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb()]),
    InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys()],
        [Law2_ScGeom_MindlinPhys_Mindlin()]
	),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,label="GSTS"),
    triax,    
    NewtonIntegrator(damping = 0.2),
    PyRunner(command = "checkPorosity()",
             iterPeriod = 1000,
             label = "checkPorosityPy"),  
    PyRunner(iterPeriod=10000,command="checkGoals()",
         label="goalChecker",dead=1) 
]

### Functions ###
phase = 1
[e11, e22, e33] = [0, 0, 0]
[e110, e220, e330] = [0, 0, 0]

def changePhase():
    # Phase 1: Isotropic compaction to reach target porosity
    # Phase 2: Isotropic compaction to reach confinement
    # Phase 3: Triax shear
    global phase
    global saveL
    global confinement
    global rate
    global key
    global sigmaIso
    if phase == 1:
        # End targetPorosity phase
        print ("Target porosity passed, porosity = ", utils.porosity())
        checkPorosityPy.dead = 1
        
        # Begin isoCompaction phase
        prepareIsoCompaction(confinement)
        
    elif phase == 2:             
        ## Start Shear phase
        phase = 3
        # Add engine to show that stresses go to zero
        O.engines = O.engines+[PyRunner(iterPeriod=300000,
                                command="showValues()",
                                label="values")]       
        prepareShear(rate,confinement)

    elif phase == 3:
        stopSimulation()            
        
    else:
        print("There's something wrong with value of phase, ", phase)
        O.pause()
        
def checkPorosity():              
    if utils.porosity()<targetPorosity:
        changePhase()
        
def prepareIsoCompaction(confinement):
    global phase
    # Begin isoCompaction phase
    phase = 2
    setContactFriction(radians(finalFricDegree))            
    O.cell.setBox(O.cell.size)           
    triax.maxStrainRate = (0.1,0.1,0.1)
    sigmaIso = -1e3*confinement # kPa          
    triax.goal[0]=triax.goal[1]=triax.goal[2]=sigmaIso
    print("Starting isocompaction simulation: ", datetime.datetime.now()) 
   
def prepareShear(rate,confinement):
    global phase
    global e110,e220,e330
    # Save current strain
    [e110, e220, e330] = triax.strain
    # So things move in the correct direction
    triax.stressMask = 7
    triax.goal[1] = -1e3*confinement
    triax.goal[0] = triax.goal[2] = -triax.goal[1]/2
    # Impose strainRate (du/dx,du/dy,du/dz,dv/dx,dv/dy,dv/dz,dw/dx,dw/dy,dw/dz)
    O.cell.velGrad=Matrix3(rate/2,0,0,0,-rate,0,0,0,rate/2)
    triax.maxStrainRate = [rate/2,rate,rate/2]
    
    goalChecker.dead = 0    
    
    showValues()
    print("Starting shear simulation: ", datetime.datetime.now())
      
def showValues():    
    # This function is just for the MWE
    print("Current iteration = ", O.iter,
          ", s11 = ", triax.stress[0],  
          ", s22 = ", triax.stress[1],
          ", s33 = ", triax.stress[2],
          ", e = ", utils.porosity()/(1-utils.porosity()))
          
        
def checkGoals():
    global phase
    global batch
    if phase == 3:
        # These values are just for the MWE
        chk1 = triax.strain[0] - e110 > 0.025
        chk2 = triax.strain[1] - e220 < -0.05
        chk3 = triax.strain[2] - e330 > 0.025
        if chk1 & chk2 & chk3:
            stopSimulation()
    else:
        print("There's something wrong with value of phase, ", phase)
        O.pause()            

def stopSimulation():
    print("Finished simulation phase: ", datetime.datetime.now())
    O.pause()

#################################
###   START SIMULATION   ###
#################################

### Create initial cloud with desired PSD (distributed by mass)
sp = SpherePack()
sp.makeCloud(mn,mx,
             psdSizes=psdSizes,
             psdCumm=psdCumm,
             periodic=True,
             distributeMass=True,
             seed=1)
print ("Cloud created, number of particles = ", len(sp))

## Send spheres to simulation
sp.toSimulation(material="spheres")

### Reach target porosity ###
print ("Compacting to target porosity: ", datetime.datetime.now())

## Porosity of cloud and targetPorosity
targetVoidRatio = e0
targetPorosity = (targetVoidRatio)/(1+targetVoidRatio)

### Start simulation ###             
O.dt = utils.PWaveTimeStep()
O.run()

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