New question #692918 on Yade: https://answers.launchpad.net/yade/+question/692918
Dear Dr. Robert Caulk, I'm trying to model three-point bending test in following paper. "Modeling acoustic emissions inheterogeneous rocks duringtensile fracture with the DiscreteElement Method" https://opengeomechanics.centre-mersenne.org/article/OGEO_2020__2__A2_0.pdf After applying the model parameters and ran a simulation to obtained the results as shown in figure 7a of the paper. The outcome is different from the reported one in the paper. I understand that maximum load ratio is obtained by the current load divided by the maximum load during the test. If I compare the vertical deflections in my simulation and the outcome from the paper, they are different.( at the peak load paper has deflection value of 0.2mm while my simulation shows 19 mm) please see the comparison here, https://www.dropbox.com/sh/dzqwt14z3jljpi5/AADLoCXg-2iOMCBYXLVsPyiua?dl=0 Could you please help me to resolve this? I have followed following steps to model the problem 1. random dense pack algorithm is used to create the particle packing with radius range mentioned in the paper. 2. apply cpm material to the packing, used interaction detection factor value of 1.5 3. attached three facetcylinder with diameter of 14mm ( top middle one for applying load, remaining two at the corner of bottom to simulate the roller support) please see the particle assembly and supports here, https://www.dropbox.com/s/5l09ujh20tz1z9b/assembly.PNG?dl=0 4. Then, constant velocity was imposed on the top middle facetcylinder to apply loading. The Minimum working example as follows, # -*- coding: utf-8 -*- # Code for three point bending test with crack mouth opening displacement and deflection from yade import pack ############################################ ### DEFINING VARIABLES AND MATERIALS ### ############################################ from yade import pack, plot,qt,ymport # Damping: damp=0.7 #Material parameters young = 30e9 poisson = .3 frictionAngle = 0.33 sigmaT = 40e6 epsCrackOnset = 3e-4 relDuctility = 30 density = 2600 # stress path ''' # 3 point bending test sample preparation parameters L= length of the sample W= width of the sample H= height of the sample C=clearance from two ends(canteliver length A= width of the crack opening B= height of the crack opening R=radius of the rollers/piston(radius of facet clyinder) all mm s (10^-3) scale Length along x axis, height along y axis and width along z axis ''' L=240e-3 #length of the sample W=40e-3 #width of the sample H=80e-3 #height of the sample R=7e-3 #radius of the rollers/piston(radius of facet clyinder) all mm s (10^-3) scale C=(12e-3+R) A=3e-3 B=50e-3 psdmean=1.5e-3 rmin=1.5*psdmean prepared=0 v=.01*10000 # applied velocity on piston concMat=O.materials.append(CpmMat(young=young,poisson=poisson,density=density,damLaw=1,frictionAngle=frictionAngle, sigmaT=sigmaT,epsCrackOnset=epsCrackOnset,neverDamage=False,label='concrete',relDuctility=relDuctility,)) frictMat = O.materials.append(CpmMat(young=100*young,poisson=poisson,frictionAngle=0,sigmaT=0,epsCrackOnset=epsCrackOnset,relDuctility=1.0001,neverDamage=True,density=3000,label='walls')) if prepared<0.5: pred=pack.inAlignedBox((0,0,0),(L,H,W)) SpherePack=pack.randomDensePack(pred,spheresInCell=500,radius=psdmean,rRelFuzz=0.25,color=(0,0,1),material=concMat,returnSpherePack=True) sand=SpherePack.toSimulation() # assign all particles to sand object else: O.bodies.append(ymport.textExt('tpb_sample.txt',format='x_y_z_r',material=concMat)) print ('I have previously prepared packing!!') for b in O.bodies: if isinstance(b.shape,Sphere): b.shape.color = (0.2,0.6,0.6)# give a colour #piston for apply load piston=O.bodies.append(geom.facetCylinder(center=(.5*L,H+R-0.3e-3,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat)) #piston radius is 7mm #rollers for support left_roller=O.bodies.append(geom.facetCylinder(center=(C,-R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat)) right_roller=O.bodies.append(geom.facetCylinder(center=(L-C,-R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat)) newton=NewtonIntegrator(damping=damp) #interaction detection factor enlargeFactor=1.5 #colour the particles for b in O.bodies: if isinstance(b.shape,Sphere): b.shape.color = (0.2,0.6,0.6) print len(O.bodies) O.engines=[ ForceResetter(), InsertionSortCollider([ Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'), Bo1_Facet_Aabb() ]), InteractionLoop( [ Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ig2ss'), Ig2_Facet_Sphere_ScGeom() ], [ Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=O.iter+5), Ip2_FrictMat_CpmMat_FrictPhys(), Ip2_FrictMat_FrictMat_FrictPhys(), ], [ Law2_ScGeom_CpmPhys_Cpm(), Law2_ScGeom_FrictPhys_CundallStrack(), ], ), ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf) #TranslationEngine(ids=[piston.id],translationAxis=[0,0,1],velocity=1), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3), PyRunner(iterPeriod=1,command='apply()',label='apply'), PyRunner(iterPeriod=1,command='addPlotData()',label='graph'), newton, ] #Display spheres with 2 colors for seeing rotations better Gl1_Sphere.stripes=0 qtr=yade.qt.Renderer() qtr.bgColor=(1,1,1) yade.qt.Controller(), yade.qt.View() #run 1 step and set interaction detection factor to 1 O.step() bo1s.aabbEnlargeFactor=ig2ss.interactionDetectionFactor=1. cylinderIds=[] for i in range(0,len(piston)): cylinderIds.append(piston[i]) #print piston[i] #print cylinderIds ''' #for i in range(0,len(left_roller)): #O.bodies[left_roller[i]].state.blockedDOFs='y' #for i in range(0,len(right_roller)): #O.bodies[right_roller[i]].state.blockedDOFs='y' #for i in range(0,len(piston)): #O.bodies[piston[i]].state.blockedDOFs='XYZxyz' ''' #function to apply load def apply(): for i in range(0,len(piston)): O.bodies[piston[i]].state.vel[1]=-v def addPlotData(): f_2=utils.sumForces(ids=cylinderIds,direction =(0,1,0))# measure total force induced on the cylinder delta_z=v*O.time*1000 #beam deflection plot.addData(t = O.time,f_2 = f_2,delta_z=delta_z) plot.plots={'delta_z':('f_2',)} plot.plot() from yade import qt qt.Controller() qt.View() -- 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
