Question #693437 on Yade changed:
https://answers.launchpad.net/yade/+question/693437
Status: Needs information => Open
Li Zeng gave more information on the question:
Hi Jan,
below is new one.
from __future__ import print_function
from yade import plot,pack
"""
A simple grain test
A sphere with radius 2.5e-3
Load in z direction by z-perpendicslar walls.
"""
# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(noTableOk=True, # unknownOk=True,
weibullCutOffMax=10,
weibullCutOffMin=0.1,
xSectionShape = 2,
xSectionScale = 0,
momentRadiusFactor=4,
young = 4e10,
cohesion = 3.5e7,
poisson = .4,
density = 2706,
frictionAngle=0.7,
iterper=1000,
intRadius =1.329,
velocity = 1e-2,
specimenRadius = .024,
sphereRadius = 1.65e-3,
rRelFuzz=0.33333,
)
from yade.params.table import *
identifier='shape-velocity
'+str(velocity)+'-weibullshape'+str(xSectionShape)+'-momentRadiusFactor'+str(momentRadiusFactor)
outputDir='out_'+identifier
if not os.path.exists(outputDir):
os.mkdir(outputDir)
if not os.path.exists('txt'):
os.mkdir('txt')
output = './'+outputDir+'/'+identifier
# material
JCFmat = O.materials.append(JCFpmMat(young=young, cohesion = cohesion,
density=density, frictionAngle = frictionAngle, poisson=poisson,
jointTensileStrength=3.5e7, jointNormalStiffness = 4e10,
jointShearStiffness = 1.6e10, jointFrictionAngle= radians(30),
label='JCFpmMat'))
# spheres
sp=pack.randomDensePack(
pack.inSphere((0,0,0),specimenRadius),
radius = sphereRadius, rRelFuzz = rRelFuzz, returnSpherePack = True
)
sp.toSimulation()
# walls
wallMat =
O.materials.append(FrictMat(young=100e10,poisson=.05,frictionAngle=.25,density=70000,
label='frictionlessWalls'))
zMin,zMax=[pt[2] for pt in aabbExtrema()]
wallIDs = O.bodies.append([wall((0,0,z),2) for z in (zMin,zMax)])
walls = wallMin,wallMax = [O.bodies[i] for i in wallIDs]
v = velocity
wallMin.state.vel = (0,0,0)
wallMax.state.vel = (0,0,-v)
# engines
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Wall_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(
cohesiveTresholdIteration=10, label='jcf',
xSectionWeibullShapeParameter=xSectionShape,weibullCutOffMin=weibullCutOffMin,
weibullCutOffMax=weibullCutOffMax
)],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,Key=identifier,
label='interactionLaw',
neverErase=True,recordCracks=True,recordMoments=True,momentRadiusFactor=momentRadiusFactor),
Law2_ScGeom_FrictPhys_CundallStrack()],
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.4,
defaultDt=0.1*utils.PWaveTimeStep()),
VTKRecorder(dead=0,iterPeriod=iterper*2,initRun=True,fileName=(output+'-'),recorders=['jcfpm','cracks','facets','moments','intr'],Key=identifier,label='vtk'),
NewtonIntegrator(damping=0.7),
PyRunner(iterPeriod=100,command='addPlotData()',initRun=True),
PyRunner(iterPeriod=100,command='stopIfDamaged()'),
]
# stop condition
def stopIfDamaged():
if O.iter < 5000: # do nothing at the beginning
return
fMax = max(plot.data["f"])
f = plot.data["f"][-1]
if f/fMax < 0.6:
print("Damaged, stopping.")
print("ft = ",max(plot.data["f"]))
O.pause()
# plot stuff
def addPlotData():
# forces of walls. f1 is "down", f2 is "up" (f1 needs to be negated for
evlauation)
f1,f2 = [O.forces.f(i)[2] for i in wallIDs]
f1 *= -1
# average force
f = .5*(f1+f2)
# displacement
wall = O.bodies[wallIDs[1]]
dspl = -1* wall.state.displ()[2]
# store values
plot.addData(
t = O.time,
i = O.iter,
dspl = dspl,
f1 = f1,
f2 = f2,
f = f,
)
plot.saveDataTxt('txt/data'+identifier+'.txt',vars= ('i', 'dspl','f'))
# plot dspl on x axis, loadon y1 axis and f,f1,f2 in y2 axis
plot.plots={'dspl':('f',)}
O.dt = 0
O.step() # to create initial contacts
# now reset the interaction radius and go ahead
# time step
O.dt = 0.004 * utils.PWaveTimeStep()
# run simulation
plot.plot()
O.run()
# when running with yade-batch, the script must not finish until the simulation
is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()
regards,
Li
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
