Yade

Stresses drop to zero in periodic triaxial shear simulation

Asked by Diego Abarca

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

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Diego Abarca
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said : 		#1

Hello,

In your code
###
def prepareShear(rate,confinement):
    ...
    triax.stressMask = 7
    triax.goal[1] = -1e3*confinement
    triax.goal[0] = triax.goal[2] = -triax.goal[1]/2
    O.cell.velGrad=Matrix3(rate/2,0,0,0,-rate,0,0,0,rate/2)
   ...
###

> triax.stressMask = 7
> triax.goal[1] = -1e3*confinement

you set stress-controll, so strain/volume is NOT kept constant (print also strain/volume, not only stress) but is the result of simulation (trying to reach defined stress).

> triax.goal[0] = triax.goal[2] = -triax.goal[1]/2

you require **tensile stress** in two directions. As the triax engine tries to reach this, it enlarges the cell, probably causing volume enlargement and "particle separation". I suspect this also causes the final zero stress.

> O.cell.velGrad=Matrix3(rate/2,0,0,0,-rate,0,0,0,rate/2)

this is overwritten by PeriTriaxController.

Proposed solution:
- set triax.dead = True (not to overwrite cell deformation)
- set O.cell.velGrad as it is now (without triax, the volume should be kept constant)

cheers
Jan

Revision history for this message
Diego Abarca (diabag) said : 		#2

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.

Revision history for this message
Jan Stránský (honzik) said : 		#3

> ", s11 = ", utils.getStress()[1,1],
> ", s22 = ", utils.getStress()[1,1],
> ", s33 = ", utils.getStress()[1,1],

use correct indices :-)

I have tried your script, with the same result.
Try to print number of interactions - the number is also decreasing.

cheers
Jan

Revision history for this message
Diego Abarca (diabag) said : 		#4

Thanks Jan, I modified the showValues function so it includes your comments:

def showValues():
    # This function is just for the MWE
    print("Current iteration = ", O.iter,
          ", s11 = ", utils.getStress()[0,0],
          ", s22 = ", utils.getStress()[1,1],
          ", s33 = ", utils.getStress()[2,2],
          ", dvol = ", O.cell.volume-vol0,
          ", coorNum = ", utils.avgNumInteractions(skipFree=False))

And these are the results I get:

Start: 2020-10-22 10:16:24.653917
Cloud created, number of particles = 480
Compacting to target porosity: 2020-10-22 10:16:24.674490
Target porosity passed, porosity = 0.49986744040494263
Starting isocompaction simulation: 2020-10-22 10:16:47.926089
Starting shear simulation: 2020-10-22 10:19:05.654546
Current iteration = 1563876 , s11 = -99919.50024856672 , s22 = -100201.48450320616 , s33 = -100092.81615722248 , dvol = -6.617444900424222e-24 , coorNum = 2.5125
Current iteration = 1863876 , s11 = 0.0007813920667472948 , s22 = -0.018439345804731364 , s33 = -0.0008659837088851886 , dvol = -2.0696058926076752e-19 , coorNum = 0.004166666666666667
Current iteration = 2163876 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , dvol = -3.751303782597383e-18 , coorNum = 0.0
Current iteration = 2463876 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , dvol = -8.511132637799019e-18 , coorNum = 0.0
Current iteration = 2763876 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , dvol = -1.3271298982690577e-17 , coorNum = 0.0
Current iteration = 3063876 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , dvol = -1.7382916022671158e-17 , coorNum = 0.0
Current iteration = 3363876 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , dvol = -2.150842307949773e-17 , coorNum = 0.0
Current iteration = 3663876 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , dvol = -2.686728951686537e-17 , coorNum = 0.0
Current iteration = 3963876 , s11 = 0.0 , s22 = 0.0 , s33 = 0.0 , dvol = -3.423459680362606e-17 , coorNum = 0.0
Finished simulation phase: 2020-10-22 10:22:12.170592

So indeed, the number of interactions drops to zero at the same time that stresses drop to zero, while volume stays "constant".
dvol goes from around e-24 to e-19, and I would think that is realistically like going from zero to zero (volume is of the order of 1e-8), so I still don't understand why the particles in my simulation have the space to loosen up so dramatically. Is there another way to measure volume of the cell that would show that the cell is actually growing?

Many thanks.

Revision history for this message
Jan Stránský (honzik) said : 		#5

> Is there another way to measure volume of the cell that would show that the cell is actually growing?

there are many ways, but basically equivalent to what O.cell.volume does [1,2,3]. If O.cell.volume is not changing, then the volume of the cell is not changing..

The shape is changing. Maybe it has similar effect as vibration compaction - by changing the configuration, you let particles to relax, find "more stress-free" positions.

Also I don't know Ip2_FrictMat_FrictMat_MindlinPhys() and Law2_ScGeom_MindlinPhys_Mindlin(), maybe they can have effect on this?

cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Cell.hpp#L250
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Cell.hpp#L177
[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Cell.hSize

Revision history for this message
Diego Abarca (diabag) said : 		#6

Thank you again for you time Jan,

Indeed, by all measures volume appears to stay constant. What you say about shape changing may be the issue, but it feels weird that stresses drop to zero just because of particles relaxing. I would expect particles to reaccomodate and stresses to drop under constant volume conditions in some cases, but particles floating is not something I had in mind. I guess I'll have to find another way to impose constant volume.

About Ip2_FrictMat_FrictMat_MindlinPhys() and Law2_ScGeom_MindlinPhys_Mindlin(), they originate from the contact force model "Hertz-Mindlin no micro-slip solution" that was implemented (I think) in [1]. I don't think they were the problem as I also ran my code using Ip2_FrictMat_FrictMat_FrictPhys() and Law2_ScGeom_FrictPhys_CundallStrack() and the issues persisted.

Anyways, thanks again for all the help Jan.

[1] http://www2.eng.ox.ac.uk/civil/publications/theses/modenese_pdf , section 4.2.1.

Revision history for this message
Jérôme Duriez (jduriez) said : 		#7

> it feels weird that stresses drop to zero just because of particles relaxing

Does relaxing not mean to feel less or no stress ? .. ;-)
I think it is the same in (granular) mechanics !

Revision history for this message
Diego Abarca (diabag) said : 		#8

Thank you for your answer Jérôme,

> I think it is the same in (granular) mechanics !

Indeed, from granular mechanics I know that for the shear phase of an undrained triaxial test there is a family of combinations of initial stress and void index where it would be expected that stresses drop. However, in laboratory tests with real soils under these conditions it would usually happen at larger vertical strains than the strain the model is experiencing when stresses drop to zero, the final stresses of the real soil would drop to small but nonzero values, and the effects of this drop would be smaller when starting from a lower void index.

That was what I was expecting to happen when imposing an undrained condition by maintaining constant volume. I know this can be done in YADE as it has been used to simulate soil liquefaction [1], but there are many things I could be doing wrong as I don't know the specifics of how this was achieved.

Thank you for your time Jérôme, I hope you have a nice day.

[1] http://lbezone.ust.hk/bib/b1781031, page 46 under figure 4.1.

Revision history for this message
Jan Stránský (honzik) said : 		#9

Try some other (even "non-realistic") parameters - confinement, friction, porosity (is 0.5 realistic value?), ...

You can also try some "shaking" phases before iso-confinement. A real sample undergoes at least gravity, possibly also some other effects.

> ... I don't know the specifics of how this was achieved.

have you tried to contact the author?

cheers
Jan

Revision history for this message
Diego Abarca (diabag) said : 		#10

Thanks again for your answers Jan,

I'm happy to report that using an initial porosity of 0.37 I'm getting non zero stresses during the whole phase. While I'm still not sure if the results I'm getting are what I expected, the particles are no longer floating. This means that all along particles were not close enough to guarantee that they don't have space to move freely, so I suppose I have to reframe what I think is a reasonable porosity when working with spherical particles.

This has been a problem for weeks for me, so thank you again for all your help and time.

I hope you have a nice weekend.