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