yade-users team mailing list archive

[Fedor <question707197@xxxxxxxxxxxxxxxxxxxxx>]

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

Summary changed to:
Undrained triaxial (NOT CONSTANT VOLUME)

Description changed to:
Please help! 
I have been struggled with flow engine since March and I experience depression or something. 
I set up the triaxial undrained simulation and faced some difficulties.
First thing, pore pressure responce is wierd. I could not get the increase of pore pressure up to max value. I change viscosity, time step, everyting , compressibility of fluid (diffrent bulk modules). The pore pressure responce is 1e-10 or something and goes to negative value. Please give me some suggestions how to approach a adequete results (which is corresponding to soil mechanics theory). 
  
from __future__ import print_function

import time
import datetime
import os

from yade import pack, plot, export
import numpy as np

# VTK RECORDER PARAMS
vtk_output = True
vtk_iter = 1000
vtk_dir = 'vtk_output'
#vtk_recorders = ['all']
vtk_recorders = ['spheres',
                 'facets',
                 'boxes',
                 'id',
                 'mass',
                 'clumpId',
                 'colors',
                 'mask',
                 'materialId',
                 'coordNumber',
                 'velocity',
                 'stress',
                 'force',
                 'bstresses']
if vtk_output == True:
    vtk_dir += '/' + datetime.datetime.now().strftime("%Y-%m-%d_%H.%M.%S") + '/'
    if not os.path.exists(vtk_dir):
        os.makedirs(vtk_dir)

#FIXED PARAMETERS
k=3E7
poisson=0.2
R=1e-4
dimcell = 2
density= 1e12
sigmaIso=-1e5
young=np.abs(k*sigmaIso*R*2)
targetVoid=0.77
frictionAngle=radians(30)
alphaKr=2
alphaKtw=2
etaRoll=.15


sp = pack.SpherePack()

sp.makeCloud((0,0,0),(dimcell, dimcell, dimcell),rMean=.1,seed=1)

pp = O.materials.append(CohFrictMat(
    young=young,
    poisson=poisson,
    frictionAngle=radians(30),
    density=density,
    isCohesive=False,
    momentRotationLaw=True,
    etaRoll=etaRoll,label='spheres'
    ))

walls = aabbWalls(((0, 0, 0), (dimcell, dimcell, dimcell)),
                  thickness=0,
                  material='spheres')
wallIds = O.bodies.append(walls)
sp.toSimulation(material='spheres')

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()]),
    InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(), Ig2_Box_Sphere_ScGeom6D()],
                    [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                    [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
                        useIncrementalForm=True,
                        always_use_moment_law=True),
                     Law2_ScGeom_FrictPhys_CundallStrack()]),
    TriaxialStressController(
        goal1=sigmaIso,
        goal2=sigmaIso,
        goal3=sigmaIso,
        thickness=0,
        stressMask=7,
        internalCompaction = False,
        label='triax',
    ),
    FlowEngine(dead=True,
               label="flow",
               ),
    NewtonIntegrator(damping=.2),
    PyRunner(command='addPlotData()',
             iterPeriod=500,
             label='plotter',
             dead=True),
    PyRunner(command='testHook()',
             iterPeriod=10,
             label='test_hook',
             dead=False),
    PyRunner(command='triaxFinished()',
             iterPeriod=10,
             label='triax_finished',
             dead=True),
    VTKRecorder(fileName=vtk_dir+'3d-vtk-',
                recorders=vtk_recorders,
                iterPeriod=vtk_iter,
                dead=not vtk_output)
]

O.dt = 0.001 
O.dynDt = False
print('time step',O.dt)

def triaxDone():
    u=utils.porosity()
    ev=u/(1-u)
    frictionAngle = radians(30)
    while ev > targetVoid:
        frictionAngle *= 0.99
        setContactFriction(frictionAngle)
        O.run(3000, True)
        u=utils.porosity()
        ev=u/(1-u)
        print ("Iteration: ", O.iter, "\tFriction angle: ", frictionAngle,"\tVoid:",ev)

startExx = startEyy = startEzz = None
def initPlotValues():
    global startExx
    global startEyy
    global startEzz

    startEzz=-triax.strain[2]
    startEyy=triax.strain[1]
    startExx=triax.strain[0]

def addPlotData():
    global startExx
    global startEyy
    global startEzz

    if startExx is None:
        initPlotValues()

    plot.addData(
        i=O.iter,
        Ezz=triax.strain[2] - startEzz,
        Eyy=triax.strain[1] - startEyy,
        Exx=triax.strain[0] - startExx,
        e=utils.porosity()/(1-utils.porosity()),
        q=abs(utils.getStress()[2,2]-utils.getStress()[0,0]),
        V=triax.volumetricStrain,
        p=flow.getPorePressure((1,1,1)),
        t=O.time,
    )

def testHook():
    if abs(sigmaIso)*0.99 <= abs(triax.meanStress) <= abs(sigmaIso)*1.01:
        O.pause()
        test_hook.dead = True
        print("FIRST STAGE DONE")

def Deviatoric():
    flow.dead = False
    flow.meshUpdateInterval = 5000
    #flow.useSolver = 3
    #flow.permeabilityFactor = 1
    flow.viscosity = 0.00298
    flow.bndCondIsPressure = [0, 0, 0, 0, 0, 0]
    flow.fluidBulkModulus=0.0001
    flow.bndCondValue = [0, 0, 0, 0, 0, 0]
    flow.boundaryUseMaxMin = [0, 0, 0, 0, 0, 0]
    O.dt = 0.001
    #flow.imposePressure(Vector3(triax.width/2,triax.height/2,triax.depth/2),100.001)
    triax.wall_bottom_activated = True
    triax.wall_top_activated = True
    triax.wall_back_activated = True
    triax.wall_front_activated = True
    triax.wall_left_activated = True
    triax.wall_right_activated = True
    triax.stressMask = 3
    triax.goal1 = sigmaIso
    triax.goal2 = sigmaIso
    triax.goal3 = -0.01
    triax_finished.dead = False
    setContactFriction(radians(17))
    
    O.run(1,1)
    flow.updateTriangulation=True #force remeshing to reflect new BC immediately
    newton.damping=0

    plotter.dead = False

def triaxFinished():
    if abs(triax.strain[2])*0.99 <= abs(0.05) <= abs(triax.strain[2])*1.01:
        O.pause()
        triax_finished.dead = True
        fname = "Drained1000"
        print(abs(utils.getStress()[2,2]-utils.getStress()[0,0]))
        plot.saveDataTxt(fname +'.data.bz2')
        #import sys
        #sys.exit(0)

plot.plots = {
        'Ezz ': ( 'q'),
        ' Ezz': ('V'),
        ' i ':('e'),
        'Ezz': ('p')
              }
plot.plot()

#O.run(-1, True)
#triaxDone()
#O.save('1.gz')
O.load('1.gz')
Deviatoric()
O.run(-1, False)

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