yade-users team mailing list archive

[Li Zeng <question693437@xxxxxxxxxxxxxxxxxxxxx>]

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.