yade-users team mailing list archive

[chanaka Udaya <question692918@xxxxxxxxxxxxxxxxxxxxx>]

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

chanaka Udaya posted a new comment:
Dear Dr. Robert Caluk,

Many thanks for sharing the code.

I was able to run the simulation using the provided script with JCFPM
model.(result-https://www.dropbox.com/s/7tx32am1b33pbwb/JCFPM-fail-
matched.PNG?dl=0  )

Later, I have used CPM model with same particle packing  file
(240x80mmBeam_1.5mmrad.sphres)provided by you. (result-
https://www.dropbox.com/s/cbfcut6up5dvrcp/cpm-same%20packing.PNG?dl=0 )

It worked !

Then, I tried to use randomdensepackfunction to generate the packing
instead of your packing file.

This is the line I changed
-------------------------------------------------
numSpheres=10000
mn,mx=Vector3(-0.04, -.140,-0.020),Vector3(0.04, .140, 0.020)
pred = pack.inAlignedBox(mn, mx) 
sp = pack.randomDensePack(pred, radius=0.0015, spheresInCell=numSpheres, rRelFuzz=0.25, material='JCFmat', returnSpherePack=True)
sp.toSimulation()
----------------------------------------------

Then I tried to run the program again with above changes and after some time program stopped with the error "Core dumped-Aborted")
what could be the reason.

The working CPM model script for attached image as follows,

-------------------------------------------------------------------------------------------------------------------------------

from yade import ymport, utils, pack, export
from yade import plot
from pylab import *
import os
import math
import numpy as np
import numpy.linalg as la
import time
timeStr = time.strftime('%m-%d-%Y')

readParamsFromTable(noTableOk=True,
	conf = 0,
	weibullCutOffMax=10,
	weibullCutOffMin=0.1,
	xSectionShape = 4,
	xSectionScale = 0,

)

from yade.params.table import *

mn,mx=Vector3(-0.04, -.140,-0.020),Vector3(0.04, .140, 0.020)

young = 30e9
poisson = 0.3
targetPorosity = 0.55
density = 5000
relDuctility=30
finalFricDegree = 19 
intRadius= 1.329
sigmaT = 10e6
epsCrackOnset=(10e6/30e9)
iterper=1000
cohesion=40e6
dispVel = -0.02 # m/s
highFric = 45

identifier='shape-'+str(xSectionShape)+'-'+timeStr
outputDir='out_'+identifier
if not os.path.exists(outputDir):
	os.mkdir(outputDir)
if not os.path.exists('txt'):
	os.mkdir('txt')
	
output = './'+outputDir+'/'+identifier

wallMat =
O.materials.append(FrictMat(young=80e9,poisson=.45,frictionAngle=radians(highFric),density=7000,label='frictionlessWalls'))


JCFmat=O.materials.append(CpmMat(young=young, poisson=poisson, density=density,frictionAngle=radians(finalFricDegree),
sigmaT=cohesion,epsCrackOnset=epsCrackOnset,neverDamage=False,label='JCFmat',relDuctility=0,))
#### preprocessing to get dimensions
sp = O.bodies.append(ymport.textExt('240x80mmBeam_1.5mmrad.spheres', 'x_y_z_r',color=(0,0.2,0.7), material='JCFmat'))

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

r=X/15.

O.reset()

wallMat =
O.materials.append(FrictMat(young=80e9,poisson=.45,frictionAngle=radians(highFric),density=7000,label='frictionlessWalls'))

# Rigid blocks
block1 = O.bodies.append(geom.facetCylinder(center=(xinf-0.971*r,yinf+2*r,0),radius=r,height=Z,orientation=Quaternion((0, 0, 0), 90),segmentsNumber=20,wire=False,material=wallMat)) 

block2 =
O.bodies.append(geom.facetCylinder(center=(xinf-0.915*r,ysup-2*r,0),radius=r,height=Z,orientation=Quaternion((0,
0, 0), 90),segmentsNumber=20,wire=False,material=wallMat))

# Loading piston
piston=O.bodies.append(geom.facetCylinder(center=(xsup+0.943*r,0,0),radius=r,height=Z,dynamic=False,orientation=Quaternion((1,0, 0),0),segmentsNumber=20,wire=False,material=wallMat))

p0=O.bodies[piston[0]].state.pos[1]


JCFmat=O.materials.append(CpmMat(young=young, poisson=poisson, density=density,frictionAngle=radians(finalFricDegree),
sigmaT=cohesion,epsCrackOnset=epsCrackOnset,neverDamage=False,label='JCFmat',relDuctility=0,))
# Specimen
beam = O.bodies.append(ymport.textExt('240x80mmBeam_1.5mmrad.spheres', 'x_y_z_r',color=(0,0.2,0.7), material='JCFmat'))

AEfile = 'AEcounts_'+identifier+'.txt'

f = open('txt/'+AEfile, 'w')
f.write('time iteration AEcount P deflection pistonDisp\n')
f.close

# remove rigid block DOFs
for i in block1:
	O.bodies[i].state.blockedDOFs='xyzXYZ'	

for i in block2:
	O.bodies[i].state.blockedDOFs='xyzXYZ'

# set engine list
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='Saabb'),Bo1_Facet_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='SSgeom'),Ig2_Facet_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=10, label='jcf')],
  [Law2_ScGeom_CpmPhys_Cpm(label='interactionLaw'),Law2_ScGeom_FrictPhys_CundallStrack()]
 ),

        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.4, defaultDt=0.1*utils.PWaveTimeStep()),
	TranslationEngine(ids=piston,translationAxis=(1,0,0), velocity=dispVel, label='pistonEngine'),
	TranslationEngine(ids=block1,translationAxis=(1,0,0), velocity=0, label='block1Engine'),
	TranslationEngine(ids=block2,translationAxis=(1,0,0), velocity=0, label='block2Engine'),
        #VTKRecorder(dead=0,iterPeriod=iterper*2,initRun=True,fileName=(output+'-'),recorders=['jcfpm','cracks','facets','moments','intr'],Key=identifier,label='vtk'),
	NewtonIntegrator(damping=0.2)
]

# options for AE model
#interactionLaw.momentRadiusFactor=3
#interactionLaw.clusterMoments=True
#interactionLaw.neverErase=True
#interactionLaw.useIncrementalForm=True
#interactionLaw.useStrainEnergy = True

O.dt = 0.004 * utils.PWaveTimeStep()

# collect data for active plotting and post processing
from yade import plot
O.engines=O.engines[0:7]+[PyRunner(dead=0,iterPeriod=int(iterper/4),command='stressStrainHist()',label='dataCollector')]+O.engines[7:8]

O.engines=O.engines[0:8]+[PyRunner(dead=0,iterPeriod=iterper,command='stopifDamaged()',label='damageCheck')]+O.engines[8:9]

# engage blocks to ensure force balance and remove dynamics of system
O.engines=O.engines[0:9]+[PyRunner(dead=0, iterPeriod=1,command='ensureBlockEngagement()',label='blockEngagement')]+O.engines[9:10]

def ensureBlockEngagement():
	Pblock1, Pblock2 = checkBlockEngagement()
	P = getPistonForce()
	maxVel = 0.01
	if Pblock1 != P/2.:
		multiplier = (P/2. - Pblock1)/(P/2.)
		block1Engine.velocity = multiplier*maxVel
	if Pblock2 != P/2.:
		multiplier = (P/2. - Pblock2)/(P/2.)
		block2Engine.velocity = multiplier*maxVel

def getPistonForce():
	P=0
	for i in piston:
		P += la.norm(O.forces.f(i))
	return P

def getPistonDisp():
	pistonDisp = O.bodies[piston[0]].state.displ().norm()
	block1Disp = O.bodies[block1[0]].state.displ().norm()
	block2Disp = O.bodies[block2[0]].state.displ().norm()
	return pistonDisp + (block1Disp + block2Disp)/2.

def checkBlockEngagement():
	Pblock1 = Pblock2 = 0
	for i in block1:
		Pblock1+= la.norm(O.forces.f(i))
	for i in block2:
		Pblock2+= la.norm(O.forces.f(i))
	return Pblock1, Pblock2


def stressStrainHist():
	plot.addData(
		i=O.iter,
		time = O.time,
		disp = -O.time*dispVel,
		P = getPistonForce(),
		pistonVel = pistonEngine.velocity,
		PistonDisp = getPistonDisp(),
		Pblock1 = checkBlockEngagement()[0],
		Pblock2 = checkBlockEngagement()[1])
	plot.saveDataTxt('txt/data'+identifier+'.txt',vars= ('i', 'disp','P', 'pistonVel','PistonDisp'))
	#momentFile = 'moments_'+identifier+'.txt'

	#AEcount=0
	#if os.path.isfile(momentFile):
		#AEcount = sum(1 for line in open(momentFile))

	#f = open('txt/'+AEfile, 'a')
	#P, pistonDisp = plot.data['P'],plot.data['PistonDisp']
	#f.write('%g %g %g %g %g\n' % (O.time, O.iter, AEcount, P[-1], pistonDisp[-1]))
	#f.close

def stopifDamaged():
	P=plot.data['P']
	if O.iter > 10000:
		if  max(P)>2000 and P[-1] < 0.5*max(P):
			print('failure reached')
			sigma_t = 3*max(P)*(2*(ysup-2*r))/(2*Z*X**2)
			print("Tensile strength=",sigma_t/1e6)
			O.pause()

plot.plots={'PistonDisp':('P'), 'i':('Pblock1','Pblock2')}
plot.plot(subPlots=True)

O.step()
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

O.run()

waitIfBatch()
----------------------------------------------------------------------------------------------------------

I got a bit confusion about following lines,


# Rigid blocks
block1 = O.bodies.append(geom.facetCylinder(center=(xinf-0.971*r,yinf+2*r,0),radius=r,height=Z,orientation=Quaternion((0, 0, 0), 90),segmentsNumber=20,wire=False,material=wallMat)) 

block2 =
O.bodies.append(geom.facetCylinder(center=(xinf-0.915*r,ysup-2*r,0),radius=r,height=Z,orientation=Quaternion((0,
0, 0), 90),segmentsNumber=20,wire=False,material=wallMat))

# Loading piston
piston=O.bodies.append(geom.facetCylinder(center=(xsup+0.943*r,0,0),radius=r,height=Z,dynamic=False,orientation=Quaternion((1,0, 0),0),segmentsNumber=20,wire=False,material=wallMat))

 from there, I can see

xinf-0.971*r
xinf-0.915*r
xsup+0.943*r

should this be changed to

xinf-1*r
xinf-1*r
xsup+1*r

this or, I would like to know why there are different
coefficients(0.971, 0.915,0.943)?

Could you please help me on this?

May I know how the packing is generated also?

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