# yade-users team mailing list archive

## [Question #698605]: The distribution of contact normal direction and contact normal force are not uniform after isotropic consolidation

```New question #698605 on Yade:

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

## model parameters
num_spheres=40000,
targetPorosity= .33,
confiningPressure=-50000,
## material parameters
compFricDegree=25,#contact friction during the confining phase
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
)

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(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
# 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 unb<stabilityThreshold and abs((confiningPressure-triax.meanStress)/confiningPressure)<0.001:
break

# reaching a specified porosity precisely
while triax.porosity>targetPorosity:
compFricDegree=0.95*compFricDegree
print('Friction:',compFricDegree,'porosity:',triax.porosity)
sys.stdout.flush()
while True:
O.run(500,True)
unb=unbalancedForce()
if unb<stabilityThreshold and abs((confiningPressure-triax.meanStress)/confiningPressure)<0.001:
break

# check the position of spheres
maxX=O.bodies[1].state.pos[0]
minX=O.bodies[0].state.pos[0]
maxY=O.bodies[3].state.pos[1]
minY=O.bodies[2].state.pos[1]
maxZ=O.bodies[5].state.pos[2]
minZ=O.bodies[4].state.pos[2]
ball_out_of_wall=list()
for b in O.bodies:
if b.state.pos[0]>maxX or b.state.pos[0]<minX:
ball_out_of_wall.append(b.id)
print('the walls or balls are moving too fast, the time step is too big')
elif b.state.pos[1]>maxY or b.state.pos[1]<minY:
ball_out_of_wall.append(b.id)
print('the walls or balls are moving too fast, the time step is too big')
if b.state.pos[2]>maxZ or b.state.pos[2]<minZ:
ball_out_of_wall.append(b.id)
print('the walls or balls are moving too fast, the time step is too big')
print(ball_out_of_wall)

# change the contact parameters to the final calibration value
for b in O.bodies:
b.mat.etaRoll=finaletaRoll
for i in O.interactions:
i.phys.etaRoll=finaletaRoll

print('Compacted state saved')
binsSizes, binsProc, binsSumCum=psd(bins=20,mass=True)
pylab.semilogx(binsSizes, binsProc,label='Mass PSD of (free) %d random spheres'%len(binsSizes))
pylab.show()

2. The code of  contact normal distribution is as follows:
"""Plot 3 histograms for distribution of interaction directions, in yz,xz and xy planes and
(optional but default) histogram of number of interactions per body. If sphSph only sphere-sphere interactions are considered for the 3 directions histograms.

:returns: If *noShow* is ``False``, displays the figure and returns nothing. If *noShow*, the figure object is returned without being displayed (works the same way as :yref:`yade.plot.plot`).
"""
import pylab,math
for axis in [0,1,2]:
fc=[0,0,0]
fc[axis]=1.
subp=pylab.subplot(220+axis+1,polar=True)
theta1=d[0]
theta2=[x+math.pi for x in theta1]
theta=theta1+theta2
# 1.1 makes small gaps between values (but the column is a bit decentered)
#pylab.title(['yz','xz','xy'][axis]+' plane')
pylab.text(.5,.25,['yz','xz','xy'][axis],horizontalalignment='center',verticalalignment='center',transform=subp.transAxes,fontsize='xx-large')
if numHist:
pylab.subplot(224,polar=False)
nums,counts=utils.bodyNumInteractionsHistogram(aabb if len(aabb)>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

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

--