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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
