After isotropic consolidation, I check the contact normal direction and contact normal force distribution which should be evenly distributed in all directions, and I found they are not uniform. I don't know what's wrong, can you help me, thank you! 1. The code of isotropic consolidation is as follows: ##______________ First section, generate sample_________ from __future__ import print_function from yade import pack, qt, plot from math import * import matplotlib; matplotlib.rc('axes',grid=True) import pylab nRead=readParamsFromTable( ## model parameters num_spheres=40000, targetPorosity= .33, confiningPressure=-50000, ## material parameters compFricDegree=25,#contact friction during the confining phase finalFricDegree=30,#contact friction during the deviatoric loading young=2e8, poisson=.2, density=2600, alphaKr=7.5, alphaKtw=0, competaRoll=.1, finaletaRoll=.1, etaTwist=0, normalCohesion=0, shearCohesion=0, ## fluid parameters fluidDensity=1000, dynamicViscosity=.001, ## control parameters damp=0, stabilityThreshold=.001, ## output specifications filename='suffusion', unknowOk=True ) from yade.params.table import * mn,mx=Vector3(0,0,0),Vector3(0.03,0.03,0.03) psdSizes=[0.00042,0.0005,0.00208,0.0024] psdCumm=[0.0,0.35,0.35,1.0] # create materials for spheres #shear strength is the sum of friction and adhesion, so the momentRotationLaw=True O.materials.append(CohFrictMat(alphaKr=alphaKr,alphaKtw=alphaKtw,density=density,etaRoll=competaRoll,etaTwist=etaTwist,frictionAngle=radians(compFricDegree),momentRotationLaw=True,normalCohesion=normalCohesion,poisson=poisson,shearCohesion=shearCohesion,young=young,label='spheres')) O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls')) walls=aabbWalls([mn,mx],thickness=0,material='walls') wallIds=O.bodies.append(walls) # generate particles packing sp=pack.SpherePack() sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=1) sp.toSimulation(material='spheres') O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='contact',setCohesionNow=False,setCohesionOnNewContacts=False)], [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)], ), #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), TriaxialStressController(label='triax', # specify target values and whether they are strains or stresses goal1=confiningPressure,goal2=confiningPressure,goal3=confiningPressure, stressMask=7, # type of servo-control, external walls compaction internalCompaction=False, thickness=0, ), NewtonIntegrator(damping=0) ] qt.View() O.dt=0.5*PWaveTimeStep() import sys while True: O.run(1000,True) unb=unbalancedForce() print('unbalanced force:',unb,' mean stress: ',triax.meanStress) if unbtargetPorosity: compFricDegree=0.95*compFricDegree setContactFriction(radians(compFricDegree)) print('Friction:',compFricDegree,'porosity:',triax.porosity) sys.stdout.flush() while True: O.run(500,True) unb=unbalancedForce() if unbmaxX or b.state.pos[0]maxY or b.state.pos[1]maxZ or b.state.pos[2]0 else utils.aabbExtrema()) avg=sum([nums[i]*counts[i] for i in range(len(nums))])/(1.*sum(counts)) pylab.bar(nums,counts,fc=[1,1,0],alpha=.7,align='center') pylab.xlabel('Interactions per body (avg. %g)'%avg) pylab.axvline(x=avg,linewidth=3,color='r') pylab.ylabel('Body count') if noShow: return pylab.gcf() else: pylab.savefig('fig1.jpg') pylab.ion() pylab.show() plotdirections(sphSph=True) 3. The code of contact normal force distribution is as follows: import pylab,math O.load('compactedStatesuffusion3.yade.gz') nBins=20 def contact_normal_force_distribution(nBins=20): '''calculate the contact normal force distribution param nBins: the number of histograms return: * ''' contact_normal_force_XY=[0 for i in range(nBins)] contact_normal_force_XZ=[0 for i in range(nBins)] contact_normal_force_YZ=[0 for i in range(nBins)] contact_nXY=[0 for i in range(nBins)] contact_nXZ=[0 for i in range(nBins)] contact_nYZ=[0 for i in range(nBins)] # use only sphere-sphere contacts for i in O.interactions: if not isinstance(O.bodies[i.id1].shape,Sphere): continue if not isinstance(O.bodies[i.id2].shape,Sphere): continue contact_normal=i.geom.normal fn=i.phys.normalForce.norm() angleXY=atan2(contact_normal[1],contact_normal[0]) angleXZ=atan2(contact_normal[2],contact_normal[0]) angleYZ=atan2(contact_normal[2],contact_normal[1]) if angleXY<0: angleXY=+math.pi if angleXZ<0: angleXZ=+math.pi if angleYZ<0: angleYZ=+math.pi nXY=int(angleXY/math.pi*nBins) nXZ=int(angleXY/math.pi*nBins) nYZ=int(angleXY/math.pi*nBins) print(nXY,nXZ,nYZ) if nXY==nBins: nXY=nBins-1 if nXZ==nBins: nXZ=nBins-1 if nYZ==nBins: nYZ=nBins-1 contact_normal_force_XY[nXY]+=fn contact_normal_force_XZ[nXZ]+=fn contact_normal_force_YZ[nYZ]+=fn contact_nXY[nXY]+=1 contact_nXZ[nXZ]+=1 contact_nYZ[nYZ]+=1 for i in range(nBins): contact_normal_force_XY[i]=contact_normal_force_XY[i]/contact_nXY[i] contact_normal_force_XZ[i]=contact_normal_force_XZ[i]/contact_nXZ[i] contact_normal_force_YZ[i]=contact_normal_force_YZ[i]/contact_nYZ[i] return contact_normal_force_XY, contact_normal_force_XZ, contact_normal_force_YZ, contact_nXY, contact_nXZ, contact_nYZ contact_normal_force_XY, contact_normal_force_XZ, contact_normal_force_YZ, contact_nXY, contact_nXZ, contact_nYZ = contact_normal_force_distribution(nBins=nBins) theta1=[0 for i in range(nBins)] for i in range(nBins): theta1[i]=i*math.pi/nBins+math.pi/(2*nBins) theta2=[x+math.pi for x in theta1] theta=theta1+theta2 contact_normal_force_XY+=contact_normal_force_XY contact_normal_force_XZ+=contact_normal_force_XZ contact_normal_force_YZ+=contact_normal_force_YZ contact_normal_force=list() contact_normal_force.append(contact_normal_force_XY) contact_normal_force.append(contact_normal_force_XZ) contact_normal_force.append(contact_normal_force_YZ) for axis in [0,1,2]: fc=[0,0,0] fc[axis]=1. subp=pylab.subplot(220+axis+1,polar=True) pylab.bar(theta,contact_normal_force[axis],width=math.pi/nBins,fc=fc,alpha=.7,label=['xy','xz','yz'][axis]) pylab.text(.5,.25,['xy','xz','yz'][axis],horizontalalignment='center',verticalalignment='center',transform=subp.transAxes,fontsize='xx-large') pylab.savefig('Fig_contact_normal_force_Distribution'+str(axis)+'.jpg') pylab.ion() pylab.show() The figure links of contact normal distribution and contact normal force distribution are as follows: 1. https://i.loli.net/2021/09/02/Y3nIWrmTx2FLfUc.jpg 2. https://i.loli.net/2021/09/02/oTCZmLUqjkaVB79.jpg

When	Who	Action	New State
2021-09-03	Jan Stránský	Request for more information	Needs information
2021-09-03	Zhicheng Gao	Give more information	Open
2021-09-03	Jan Stránský	Answer	Answered
2021-09-05	Zhicheng Gao	Reopen	Open
2021-09-06	Jan Stránský	Comment	Open
2021-09-12	Zhicheng Gao	Confirm	Solved

Yade

History of question #698605