About simulating constant volume condition in a triaxial test

Hi,

I would like to simulate a constant volume condition in a triaxial test (i.e., undrained triaxial test).

In general, we can simulate the constant volume by setting triax controller like:

triax.stressMask = 0 ## all strain control
triax.goal2=rate
triax.goal1=-0.5*rate
triax.goal3=-0.5*rate

For a constant volume condition, volumetric strain should be zero theoretically, while there could be some fluctuation in volumetric strain since it is a numerical simulation. It is OK if I only have a small amount of volumetric strain during the simulation (e.g. below 1e-5). However, I always get accumulation in volumetric strain to a relative large value (e.g. up to 1e-3), which doesn't satisfy the constant volume condition.

Do you have any ideas about simulating a constant volume condition more precisely? For example, what parameters we can adjust to help us simulate constant volume more accurately?

A MWE below modified from Bruno's triaxial script shows the accumulation in volumetric strain during undrained loading.
####################
from yade import pack

nRead = readParamsFromTable(
        num_spheres=1000,
        compFricDegree=30,
        key='_triax_base_',
        unknownOk=True
)
from yade.params import table

num_spheres = table.num_spheres
key = table.key
targetPorosity = 0.43
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate = -0.01
damp = 0.2
stabilityThreshold = 0.01
young = 5e6
mn, mx = Vector3(0, 0, 0), Vector3(1, 1, 1)

O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=radians(compFricDegree), density=2600, label='spheres'))
O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=0, density=0, label='walls'))

walls = aabbWalls([mn, mx], thickness=0, material='walls')
wallIds = O.bodies.append(walls)

sp = pack.SpherePack()

sp.makeCloud(mn, mx, -1, 0.3333, num_spheres, False, 0.95, seed=1)
O.bodies.append([sphere(center, rad, material='spheres') for center, rad in sp])

triax = TriaxialStressController(
        maxMultiplier=1. + 2e4 / young,
        finalMaxMultiplier=1. + 2e3 / young,
        thickness=0,

        stressMask=7,
        internalCompaction=True,
)

newton = NewtonIntegrator(damping=damp)

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
        GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=100, timestepSafetyCoefficient=0.8),
        triax,
        TriaxialStateRecorder(iterPeriod=100, file='WallStresses' + table.key),
        newton
]

Gl1_Sphere.stripes = 0
if nRead == 0:
 yade.qt.Controller(), yade.qt.View()

triax.goal1 = triax.goal2 = triax.goal3 = -10000

while 1:
        O.run(1000, True)
        unb=unbalancedForce()
        print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
        if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
                break

print "### Isotropic state saved ###"

import sys
while triax.porosity>targetPorosity:
        compFricDegree = 0.95*compFricDegree
        setContactFriction(radians(compFricDegree))
        print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
        sys.stdout.flush()
        O.run(500,1)

triax.internalCompaction=False

setContactFriction(radians(finalFricDegree))

triax.stressMask = 0
triax.goal2=rate
triax.goal1=-0.5*rate
triax.goal3=-0.5*rate

newton.damping=0.6

from yade import plot

def history():
        plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
                s11=-triax.stress(triax.wall_right_id)[0],s22=-triax.stress(triax.wall_top_id)[1],s33=-triax.stress(triax.wall_front_id)[2],i=O.iter)

O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

O.run(100,True)

plot.plots={'e22':('s11',None,'ev')}

plot.plot()
##############

Thanks
Leonard