New question #707087 on Yade: https://answers.launchpad.net/yade/+question/707087
I was conducting a triaxial test on the argillite, which contained a matrix and two inclusions, to investigate its small strain stiffness. I tried to use spherical grains and JCF model to simulate this kind of rock. And I seperated the particle into three groups and assigned them different material properties (Young's modulous and Poisson's ratio). But the results were not desired. Because the small strain stiffness of the sample is almost equivalent to the sample containg only one material, which is confusing. Theoretically, the change of materials in the sample will affect the properties of the rock. How can I realize this goal? Could you give me some advice? from yade import ymport, utils , plot O.load('isocompact200kpa2.yade.bz2') O.engines[6].dead=True #---------------- ENGINES ARE DEFINED HERE try: os.mkdir('post') except: pass PACKING='121_1k.spheres' OUT=PACKING+'_1MPa_r0.02' iterMax=50000 # maximum number of iterations saveVTK=iterMax/5 # Vtk files record interval YOUNG2=100e9 YOUNG3=50e9 ALPHA2=0.06 ALPHA3=0.27 bodies2=[b for b in O.bodies if 5.0<b.id<178.0] bodies3=[b for b in O.bodies if 177.0<b.id<349.0] for b in bodies2: b.mat.young=YOUNG2 b.mat.poisson=ALPHA2 for b in bodies3: b.mat.young=YOUNG3 b.mat.poisson=ALPHA3 confinement=-7e6 strainRate=-0.0001 triax=O.engines[3] #### triaxial Engine triax.stressMask=5 triax.goal1=confinement triax.goal2=strainRate triax.goal3=confinement triax.max_vel=1 #### simulation is defined here (DEM loop, interaction law, servo control, recording, etc...) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=-1.,label='Saabb')]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=-1.,label='SSgeom'),Ig2_Box_Sphere_ScGeom()], [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')], [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key='data',label='interactionLaw')] ), triax, GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5,defaultDt=0.1*utils.PWaveTimeStep()), NewtonIntegrator(damping=0.5,label="newton"), PyRunner(command='addPlotData()',iterPeriod=10,label='record'), VTKRecorder(iterPeriod=int(5000),initRun=True,fileName='./post/p1-',recorders=['spheres','jcfpm','cracks'],Key='data',label='vtk') ] #### custom recording functions #tensCks=shearCks=0 # if you want to plot during simulation #---------------- SIMULATION STARTS HERE #### manage interaction detection factor during the first timestep and then set default interaction range (intRadius=1) def addPlotData(): global ini_e22a global ini_e22b e22=-triax.strain[1] unb = unbalancedForce() mStress = -(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3.0 volume = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*(O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2]) Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if isinstance(b.shape,Sphere))/volume plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2], ev=-triax.strain[0]-triax.strain[1]-triax.strain[2], s11=-triax.stress(triax.wall_right_id)[0], s22=-triax.stress(triax.wall_top_id)[1], s33=-triax.stress(triax.wall_front_id)[2], ub=unbalancedForce(), dstress=(-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]), disx=triax.width, disy=triax.height, disz=triax.depth, i=O.iter, porosity=Porosity, cnumber=avgNumInteractions(), tc=interactionLaw.nbTensCracks, sc=interactionLaw.nbShearCracks, te=interactionLaw.totalTensCracksE, se=interactionLaw.totalShearCracksE, ) #if (e22 - ini_e22a) > 0.001 and e22 < 0.05: # O.save('save'+'{:.4f}'.format(e22)+'yade.bz2') # ini_e22a = e22 #if (e22 - ini_e22b) > 0.00005: print ('unbalanced force: %f, mean stress: %f, s11: %f, s22: %f, s33: %f, coordination number: %f, porosity: %f' %(unb,mStress,-triax.stress(triax.wall_right_id)[0],-triax.stress(triax.wall_top_id)[1],-triax.stress(triax.wall_front_id)[2],avgNumInteractions(),Porosity)) #Output() #ini_e22b = e22 plot.saveDataTxt('loadinglog3.txt.bz2') plot.plots={'e22':('s22',None,'ev')} plot.plot() #O.run(iterMax) -- 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