New question #701924 on Yade:
https://answers.launchpad.net/yade/+question/701924
Hi,
Due to the needs of laboratory tests, I need to use cylindrical specimens for
triaxial compression numerical simulation.
Since TriaxialStressController does not seem to be applicable to cylinder
specimens, I refer to the example in [1], but [1] does not mention measuring
transverse strain. Can anyone help?
My code as follows:
----------------------------------------
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys
rate=-0.2
damp=0.4
stabilityThreshold=0.001
key='_triax_base_'
young=100e9
name='JCFPM_triax'
compFricDegree=30
poisson=0.4
name='cylinder'
preStress=-20e6
strainRate = -1
OUT=str(preStress)+'_JCFPM_triax'
nw=24
nh=15
rParticle=0.000731723
bcCoeff = 5
width = 0.025
height = 0.05
mn,mx=Vector3(-0.025,-0.025,0),Vector3(0.025,0.025,0.1)
O.materials.append(JCFpmMat(type=1,density=2640,young=young,poisson=poisson,tensileStrength=15e6,cohesion=20e6,frictionAngle=radians(25),label='sphere'))
O.materials.append(JCFpmMat(type=1,label='wall',young=young,poisson=poisson,frictionAngle=radians(25)))
spheres=O.bodies.append(ymport.text('packing-'+name+'.spheres',scale=1,shift=Vector3(0,0,0),material='sphere',color=(0.6,0.5,0.15)))
bot = [O.bodies[s] for s in spheres if
O.bodies[s].state.pos[2]<rParticle*bcCoeff]
top = [O.bodies[s] for s in spheres if
O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
vel = strainRate*(height-rParticle*2*bcCoeff)
for s in bot:
s.shape.color = (1,0,0)
s.state.blockedDOFs = 'xyzXYZ'
s.state.vel = (0,0,-vel)
for s in top:
s.shape.color = Vector3(0,1,0)
s.state.blockedDOFs = 'xyzXYZ'
s.state.vel = (0,0,vel)
facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
for h in range(nh):
v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),
rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),
rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),
rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),
rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
f1 = facet((v1,v2,v3),color=(0,0,1),material='wall')
f2 = facet((v1,v3,v4),color=(0,0,1),material='wall')
facets.extend((f1,f2))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
f.state.mass = mass
f.state.blockedDOFs = 'XYZz'
plot.plots = { 'e':('s',), }
def addForces():
for f in facets:
n = f.shape.normal
a = f.shape.area
O.forces.addF(f.id,preStress*a*n)
def stopIfDamaged(maxEps=0.001):
extremum = max(abs(s) for s in plot.data['s'])
s = abs(plot.data['s'][-1])
e = abs(plot.data['e'][-1])
if O.iter < 1000 or s > .5*extremum and e < maxEps:
return
if abs(s)/abs(extremum) < 0.8 :
print('Simulation finished')
presentcohesive_count = 0
for i in O.interactions:
if hasattr(i.phys, 'isCohesive'):
if i.phys.isCohesive == True:
presentcohesive_count+=1
print('the number of cohesive bond now
is:',presentcohesive_count)
print('Max stress and strain:',extremum,e)
O.pause()
O.dt=.5*utils.PWaveTimeStep()
newton=NewtonIntegrator(damping=damp)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ss2sc'),Ig2_Facet_Sphere_ScGeom()],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT+'_Crack',label='interactionLaw')]
),
PyRunner(iterPeriod=1,command="addForces()"),
newton,
PyRunner(command='plotAddData()',iterPeriod=10),
PyRunner(iterPeriod=50,command='stopIfDamaged()'),
VTKRecorder(iterPeriod=2000,initRun=True,fileName='triax/JFFPM-',recorders=['spheres','jcfpm','cracks'],Key=OUT+'_Crack',label='vtk',dead=0),
]
def plotAddData():
f1 = sum(O.forces.f(b.id)[2] for b in top)
f2 = sum(O.forces.f(b.id)[2] for b in bot)
f = .5*(f2-f1)
s = f/(pi*.25*width*width)
e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) /
(height-rParticle*2*bcCoeff)
plot.addData(
i = O.iter,
s = -s,
e = -e,
tc = interactionLaw.nbTensCracks,
sc = interactionLaw.nbShearCracks,
)
plot.saveDataTxt(OUT)
O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1
cohesiveCount = 0
for i in O.interactions:
if hasattr(i.phys, 'isCohesive'):
if i.phys.isCohesive == True:
cohesiveCount+=1
print('the origin total number of cohesive bond is:',cohesiveCount)
plot.plot()
O.run()
----------------------------------------------
Thanks for kindly help!
[1]https://gitlab.com/yade-dev/trunk/-/blob/master/examples/concrete/triax.py
--
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
