yade-users team mailing list archive

[Diego Abarca <question693568@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Diego Abarca is still having a problem:
Thank you for your answer and for your time Jan,

I implemented your suggestions and rewrote my code so
PeriTriaxController is dead during the shear phase, however stresses
still rapidly fall to zero. In the modified MWE I changed the output to
verify that there are no significant changes in volume during shear and
from the results it appears that volume is indeed constant throughout
the phase.  I didn't know that PeriTriaxController overwrites
O.cell.velGrad, so that was a changed that I had to make anyways,
however it appears that there's yet another problem in my implementation
that drives stresses to zero in the shear phase.

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 `yuessk'
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 18:31:19.944231
Cloud created, number of particles =  480
Compacting to target porosity:  2020-10-21 18:31:19.962880
[[ ^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 18:31:43.531963
Starting shear simulation:  2020-10-21 18:33:52.318887
Current iteration =  1510831 , s11 =  -99810.48919734705 , s22 =  -99810.48919734705 , s33 =  -99810.48919734705 , vol =  4.181077915190452e-08 , dvol =  6.617444900424222e-24
Current iteration =  1810831 , s11 =  -101110.82061300645 , s22 =  -101110.82061300645 , s33 =  -101110.82061300645 , vol =  4.181077915189381e-08 , dvol =  -1.070702584888639e-20
Current iteration =  2110831 , s11 =  -85591.4815078114 , s22 =  -85591.4815078114 , s33 =  -85591.4815078114 , vol =  4.181077915188469e-08 , dvol =  -1.9825864921670967e-20
Current iteration =  2410831 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , vol =  4.1810779149430266e-08 , dvol =  -2.4742494133788155e-18
Current iteration =  2710831 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , vol =  4.181077914504756e-08 , dvol =  -6.856956701150176e-18
Current iteration =  3010831 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , vol =  4.1810779139635944e-08 , dvol =  -1.2268570791819095e-17
Current iteration =  3310831 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , vol =  4.1810779131895016e-08 , dvol =  -2.000949901765774e-17
Current iteration =  3610831 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , vol =  4.181077912360294e-08 , dvol =  -2.8301574376918016e-17
Current iteration =  3910831 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , vol =  4.181077911514276e-08 , dvol =  -3.676175268898377e-17
Current iteration =  4210831 , s11 =  0.0 , s22 =  0.0 , s33 =  0.0 , vol =  4.1810779112342816e-08 , dvol =  -3.9561699505471663e-17
Finished simulation phase:  2020-10-21 18:38:15.582393

Here is the modified MWE:

# -*- 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]
vol0 = 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
    global vol0
    # Save current strain
    [e110, e220, e330] = triax.strain
    # Save current volume
    vol0 = O.cell.volume
    # 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.dead = 1
    goalChecker.dead = 0    
    print("Starting shear simulation: ", datetime.datetime.now())
      
def showValues():    
    # This function is just for the MWE
    print("Current iteration = ", O.iter,
          ", s11 = ", utils.getStress()[1,1],  
          ", s22 = ", utils.getStress()[1,1],
          ", s33 = ", utils.getStress()[1,1],
          ", vol = ", O.cell.volume,
          ", dvol = ", O.cell.volume-vol0)
          
        
def checkGoals():
    global phase
    global batch
    if phase == 3:
        # These values are just for the MWE
        chk1 = log(O.cell.trsf[0,0]) - e110 > 0.025
        chk2 = log(O.cell.trsf[1,1]) - e220 < -0.05
        chk3 = log(O.cell.trsf[2,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()

Thanks again for your time, I hope you have a great day.

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