New question #686273 on Yade: https://answers.launchpad.net/yade/+question/686273
Hello All, I want to simulate a triaxial test of sample of granular material with real shape. I am using clumps instead of spheres. I want to plot micro strain of sampel in paraview. Is it possible to get micro strain when the particles are clumps? Because my code does not work. Thank you very much for your help in advance. Best Regards, Alireza ########my code######### O.reset() from yade import utils, plot from yade import pack, qt from datetime import datetime #==================================COKE AGGREGATE CLUMP TEMPLATES================================================================ radz1=[0.355155e-3,0.505113e-3,0.397713e-3,0.465286e-3,0.484395e-3,0.394534e-3,0.493151e-3,0.487328e-3,0.613619e-3,0.413455e-3] poz1= [[1.13418e-3,-0.703895e-3,-1.20338e-3],[-0.390408e-3,0.476061e-3,-0.150612e-3],[-0.556545e-3,0.451341e-3,1.1495e-3],[-0.633942e-3,0.498253e-3,0.348231e-3],[0.0256934e-3,0.388855e-3,-0.733445e-3],[-0.218563e-3,0.504478e-3,1.54117e-3],[-0.319601e-3,0.104778e-3,0.742895e-3],[0.650678e-3,-0.76675e-3,-0.289908e-3],[0.0113115e-3,-0.207684e-3,0.00255944e-3],[0.594902e-3,-0.301473e-3,-0.878654e-3]] template1= [] template1.append(clumpTemplate(relRadii=radz1,relPositions=poz1)) radz2=[0.330164e-3,0.504115e-3,0.399587e-3,0.614205e-3,0.466444e-3,0.495302e-3,0.394324e-3,0.486898e-3,0.444037e-3,0.489396e-3] poz2= [[1.16386e-3,-0.70706e-3,-1.25222e-3],[-0.39596e-3,0.482385e-3,-0.144097e-3],[-0.554928e-3,0.448475e-3,1.14685e-3],[-0.00361766e-3,-0.198211e-3, 0.00106559e-3],[-0.633362e-3,0.490424e-3,0.357391e-3],[0.0434148e-3,0.367924e-3,-0.736319e-3],[-0.218749e-3,0.504703e-3,1.54165e-3],[-0.311777e-3, 0.101954e-3,0.750235e-3],[0.621565e-3,-0.779387e-3,-0.179029e-3],[0.711114e-3,-0.503202e-3,-0.795602e-3]] template2= [] template2.append(clumpTemplate(relRadii=radz2,relPositions=poz2)) radz3=[0.959101e-3,0.774782e-3,0.924682e-3,0.882986e-3,0.572207e-3,0.875338e-3,0.499054e-3,0.582586e-3,1.03184e-3,0.561428e-3] poz3=[[-0.742455e-3,-0.539002e-3,0.030705e-3],[0.800775e-3,1.00778e-3,0.366095e-3],[-0.217508e-3,-0.423924e-3,-0.671114e-3],[0.65016e-3,0.741405e-3,0.305558e-3],[-0.770717e-3,0.423006e-3,0.33259e-3],[0.20475e-3,-0.707669e-3,-0.332745e-3],[-0.0356616e-3,1.34122e-3,0.703809e-3],[0.64889e-3,0.340976e-3,-0.456228e-3],[0.258817e-3,0.13357e-3,0.178886e-3],[0.35466e-3,1.36934e-3,0.583507e-3]] template3= [] template3.append(clumpTemplate(relRadii=radz3,relPositions=poz3)) radz4=[1.18842e-3,1.0381e-3,0.664711e-3,0.531753e-3,0.853808e-3,0.789023e-3,0.717061e-3,1.0081e-3,0.967001e-3,0.535189e-3] poz4=[[-0.0354242e-3,-0.245642e-3,-0.0514299e-3],[0.444228e-3,-0.342353e-3,-0.327639e-3],[-0.815227e-3,1.41647e-3,-0.246128e-3],[1.38493e-3,-0.405625e-3,-0.69207e-3],[0.608272e-3,0.0171296e-3,0.863712e-3],[-0.191982e-3,-1.05752e-3,-0.341411e-3],[-1.08258e-3,-0.336001e-3,-0.482846e-3],[-0.396508e-3,0.47029e-3,0.0820364e-3],[0.486211e-3,0.345623e-3,0.492724e-3],[-0.0502768e-3,0.977767e-3,0.467401e-3]] template4= [] template4.append(clumpTemplate(relRadii=radz4,relPositions=poz4)) #================= define the materials ======================= O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=1.95e7, density=1532.2, poisson=0.3, frictionAngle= 0.0, fragile=False, label='aggregate-814')) O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=4e9, density=1523.6, poisson=0.3, frictionAngle= 0.0, fragile=False, label='wall')) #========= creating walls ====================== #fc=yade.geom.facetCylinder((0,0,0.05),12.7e-3, 0.1, orientation=Quaternion((0, 0, 1), 0), segmentsNumber=100, wallMask=7, angleRange=None, closeGap=False, radiusTopInner=-1, radiusBottomInner=-1,material='wall') #O.bodies.append(fc) walls=aabbWalls([(-0.05,-0.05,-0.05),(0.05,0.05,0.05)],thickness=0.0003,oversizeFactor=1.0,material='wall') wallIds=O.bodies.append(walls) ############################ ### DEFINING ENGINES ### ############################ O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]), InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()], [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()] ), # TriaxialStateRecorder(iterPeriod=1000,file='WallStresses'), NewtonIntegrator(damping=0.4,gravity=[0,0,-10]) ,PyRunner(command='calm()',iterPeriod=10,label='calmEngine') ] O.dt=1e-7 #====================================================Generating the aggregates=================================================== coke=((1.875e-3,100),(0.9475e-3,367),(0.50125e-3,500),(0.50125e-3,500)) nums=['t','t','t','t'] temps=[template1,template2,template3,template4] mats=['aggregate-814','aggregate-814','aggregate-814', 'aggregate-814'] for i in range(len(nums)): nums[i]=pack.SpherePack() nums[i].makeCloud((-0.045,-0.045,-0.045),(0.045,0.045,0.045),rMean=coke[i][0],rRelFuzz=0.0,num=coke[i][1]) O.bodies.append([utils.sphere(c,r,material=mats[i]) for c,r in nums[i]]) O.bodies.replaceByClumps(temps[i],[1.0],discretization=5) O.step() print 'number of interactions is ',len([i for i in O.interactions]) print '' print '' print '' ############################ ### Display settings ### ############################ for x in range(len(O.bodies)): if (O.bodies[x]): if isinstance(O.bodies[x].shape,Sphere): if O.bodies[x].isStandalone: O.bodies[x].shape.color=(0.2,0.3,0.6) else: O.bodies[x].shape.color=(0.7,0.63,0.0) else: O.bodies[x].shape.color=(0.2,0.3,0.1) qtr=qt.Renderer() qtr.bgColor=(1,1,1) O.engines=O.engines+[PyRunner(command='check_it()',iterPeriod=500,label='checkEngine')] def check_it(): print 'number of interactions is ',len([i for i in O.interactions]) pp=max(i.geom.penetrationDepth for i in O.interactions) print 'max pen is ',pp print 'unbalanced force is: ',unbalancedForce() print '' #O.run(5000,True) #=============================================== #=============== Compaction ==================== #=============================================== triax=TriaxialStressController( maxMultiplier=1.000, finalMaxMultiplier=1.000, thickness = 0, stressMask = 7, internalCompaction=False, ) O.engines=O.engines+[triax] triax.goal1=-1.0e5 triax.goal2=-1.0e5 triax.goal3=-1.0e5 triax.wall_back_activated=True sigmaIso=-1.0e5 O.dt=2e-6 calmEngine.dead=True checkEngine.dead=True while 1: O.run(100, True) unb=unbalancedForce() meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3 print 'mean Stress:',triax.meanStress,'porosity:', triax.porosity,'meanS:',unb if triax.porosity<0.385 and abs(sigmaIso-triax.meanStress)/abs(sigmaIso)<0.01: print 'Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:',-triax.strain[2] break print "### Isotropic state saved ###" triax.depth0=triax.depth triax.height0=triax.height triax.width0=triax.width O.save('RVE-sizeDis-solid-Isoe5-Isopart.yade') #==================================================================== #===================== DEVIATORIC LOADING ======================= #==================================================================== triax.stressMask = 3 triax.goal1 = sigmaIso triax.goal2 = sigmaIso triax.goal3 = -0.1 #==================================================================== #==================== Micro Strain =========================== #==================================================================== TW=TesselationWrapper() TW.computeVolumes() TW.volume(20) TW.setState(0) step=0 while 1: step +=1 O.run(1000,True) TW.setState(1) TW.defToVtk("strain_"+str(step)+".vtk") O.engines=O.engines+[PyRunner(iterPeriod=20000,command='finishIt()',label='calmEngine')] def finishIt(): if abs(-triax.strain[2]) > 0.5: # makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000) print O.iter print 'time is ', O.time print 'porosity:',triax.porosity print 'strain 1:',-triax.strain[0], 'strain 2:',-triax.strain[1], 'strain 3:',-triax.strain[2] O.save('RVE-sizeDis-solid-Isoe5-Devpart.yade') O.pause() #================================================================= #======================= data collecting ========================= #================================================================= if 1: O.engines=O.engines[0:5]+[PyRunner(iterPeriod=100,command='history()',label='recorder')]+O.engines[5:7] else: O.engines[4]=PyRunner(iterPeriod=100,command='history()',label='recorder') def history(): plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2], ev=(triax.strain[0]+triax.strain[1]+triax.strain[2]), sxx=-triax.stress(triax.wall_right_id)[0], syy=-triax.stress(triax.wall_top_id)[1], szz=-triax.stress(triax.wall_front_id)[2], mass=sum(b.state.mass for b in O.bodies if isinstance(b.shape,Sphere)), q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]), p=triax.meanStress, R=triax.stress(triax.wall_front_id)[2]/sigmaIso, por=triax.porosity, i=O.iter), plot.saveDataTxt('RVE-sizeDis-solid-Isoe5-Devpart.txt') -- 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