yade-users team mailing list archive

[Fedor <question706740@xxxxxxxxxxxxxxxxxxxxx>]

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

Fedor posted a new comment:
Hello Bruno, thank you for reply. 
Script is attached below. I want to see as in theory the result of increasing of pore pressure up to a constant constant value. 



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 = 100
vtk_dir = 'vtk_output'
vtk_recorders = ['all']
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 = 0.002
density= 1e12
sigmaIso=-1e5
young=np.abs(k*sigmaIso*R*2)
targetVoid=0.7
frictionAngle=radians(30)
alphaKr=2
alphaKtw=2
etaRoll=.15

#SETTINGS
# removed O.periodic = True
# removed O.cell.hSize = Matrix3(dimcell , 0, 0, 0, dimcell , 0, 0, 0, dimcell )


sp = pack.SpherePack()


diameter=[0.000075,0.000106,0.00015,0.00025,0.00030,0.000425,0.000850]

part=[0.0037,0.0263,0.1319,0.7136,0.8681,0.9939,1]


sp.makeCloud((0, 0, 0), (dimcell, dimcell, dimcell), psdSizes=diameter,psdCumm=part,distributeMass=True,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,
        max_vel=0.05,
        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.1 * PWaveTimeStep()
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)


def addPlotData():
    plot.addData(
        i=O.iter,
        Ezz=-triax.strain[2],
        Eyy=triax.strain[1],
        Exx=triax.strain[0],
        e=utils.porosity()/(1-utils.porosity()),
        q=abs(utils.getStress()[2,2]-utils.getStress()[0,0]),
        V=triax.volumetricStrain,
        p=flow.getPorePressure((0.0015,0.0015,0.0015)),
        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.defTolerance = 0.8
    flow.meshUpdateInterval = 800
    #flow.useSolver = 3
    flow.permeabilityFactor = 1
    flow.viscosity = 0.00298
    flow.bndCondIsPressure = [0, 0, 0, 0, 0, 0]
    flow.fluidBulkModulus=0.00001  
    flow.bndCondValue = [0, 0, 0, 0, 0, 0]
    flow.boundaryUseMaxMin = [0, 0, 0, 0, 0, 0]

    #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.0001
    triax.max_vel=0.0001
    triax_finished.dead = False
    setContactFriction(radians(17))
    O.dt = 0.1e-3
      

    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.01) <= 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 = {
        'Eyy ': ( 'q'),
        ' Exx': ('V'),
        ' i ':('e'),
        'Ezz': ('p')
              }
plot.plot()

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

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