New question #699779 on Yade: https://answers.launchpad.net/yade/+question/699779
Hi everybody I want to carry out triaxial simulation on sand with considering a non-linear contact law (Hertz-Mindlin). I mention the code in below which is working very well with Cohesion moment law but in terms of Hertz- MIndlin there are some problems in the case of creating a pack such as an unacceptable unbalanced (almost 0.67) force for desire void ratio which is around 0.83. Do you know what should I do to well-create a pack with this non-linear contact law ? Also, I need to get to know the Hertz-Mindlin equations that are using on YADE, do you know the references of this contact law? # encoding: utf-8 import matplotlib; matplotlib.rc('axes',grid=True) import matplotlib.pyplot as plt from yade import pack import pylab from numpy import * import numpy as np from math import * utils.readParamsFromTable(seed=1,num_spheres=10000,compFricDegree = 25.0) from yade.params import table seed=table.seed num_spheres=table.num_spheres# number of spheres compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process) confiningS=-1e5 rate=5 psdSizes,psdCumm=[0.00012,0.00014,0.00018,0.00021,],[0.0030,0.039,0.21,1] psdSizesArray=np.array(psdSizes) psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter sp=pack.SpherePack() mn,mx=Vector3(0,0,0),Vector3(0.0001,0.0001,0.0001) sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed) sp.psd(bins=10,mass=True) ## create material #0, which will be used as default O.materials.append(FrictMat(young=5e9,poisson=.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres')) O.materials.append(FrictMat(young=5e9,poisson=.3,frictionAngle=0,density=0,label='frictionless')) ## create walls around the packing walls=aabbWalls((mn,mx),thickness=0,material='frictionless') wallIds=O.bodies.append(walls) O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp]) triax=TriaxialStressController( internalCompaction=False, goal1=confiningS, goal2=confiningS, goal3=confiningS, max_vel=10, label="triax" ) newton=NewtonIntegrator(damping=0) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_MindlinPhys()], [Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=False,includeMoment=False)] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'), triax, newton ] #s11=[];e22=[] while 1: O.run(1000,True) unb=unbalancedForce() print ('unb=', unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1),'e=', triax.porosity/(1-triax.porosity), 'Size=', (triax.width)) if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01 and (triax.porosity)/(1-triax.porosity)<0.8399999999 : break print ("Particle.Distribution.is.stable") ######################################################## ##Turn into triaxil test ######################################################## finalFricDegree = 18.0 triax.internalCompaction=False setContactFriction(radians(finalFricDegree)) while 1: O.run(1000,True) unb=unbalancedForce() if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01: break #ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2] #si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2] triax.depth0=triax.depth triax.height0=triax.height triax.width0=triax.width triax.stressMask=5 triax.goal1=triax.goal3=confiningS triax.goal2=-rate triax.wall_bottom_activated=0 #O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True VTK.dead=1 #cohesiveLaw.always_use_moment_law=False #recorder.dead=0 Exaxil=[5.98e-04,6.11e-04,0.00125,0.0013,0.00249,0.0031,0.00491,0.00729,0.010,0.017,0.025] Exdev=[8.69,16.30,51.09,81.52,92.40,110.88,142.40,170.67,196.77,218.54,229.44] ee=[];e22=[];s22=[] file=open('Results'+'.txt',"w") while 1: O.run(1000,True) ee.append((triax.porosity)/(1-triax.porosity)) s22.append(((-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3) e22.append((-triax.strain[1])) print('ax=',(-triax.strain[1]), 'deviatoric=',(((-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3), 'ee=', (triax.porosity)/(1-triax.porosity), 's11=',((-triax.stress(triax.wall_top_id)[1])*1e-3), 's33=',(((-triax.stress(triax.wall_right_id)[0])-(triax.stress(triax.wall_front_id)[2]))/2.0)*1e-3) if abs(-triax.strain[1])>0.1 : break print(" Triaxil test is finished ") file.write(str(s22)+" "+str(e22)+"\n") file.write(str(s22)+" "+str(ee)+"\n") file.close() plt.figure() plt.xlabel('Axial strain') plt.ylabel('deviatoric') plt.grid(True) plt.plot(Exaxil,Exdev,linestyle='-',marker="*", markersize=8,color='g',linewidth=2) plt.plot(e22,s22,linestyle='-',marker="^", markersize=8,color='r',linewidth=2) plt.legend(('Exprimental data','simulation data'),loc='upper right', shadow=True) plt.show() -- 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