# The distribution of contact normal direction and contact normal force are not uniform after isotropic consolidation

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

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Zhicheng Gao
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2021-09-03: #1

Hello,

Please provide a MWE. M = minimal, i.e. also make the computation time minimal if possible (which is this case, see below).
I am running your code for a while (more then 1 hour) with veeeeery slow convergence.
If you provide the final packing and a script to load it, make one iteration to create interactions and O.save for the plot scripts, it is sufficient.

> 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.

What is your definition of "even distribution" and "uniform"? They cannot be perfect from the discrete nature of the simulation.

> def plotdirections(...

any difference from utils.plotDirections? what difference? why?

Why are you using different code for forces and directions? Aren't they essentially the same?
If you have code for forces, you can directly use it for directions (?)

The picture of force distribution seems overall ok to me (as discussed, you cannot expect perfect circle).

The force distributions are exactly the same for all 3 sections, should not be like this (!?)

> contact_normal_force_XY[nXY]+=fn
> contact_nXY[nXY]+=1

I would cut out "almost perpendicular" forces, or used projection of the force to the plane.
This way, a perpendicular force is assigned its full strength to a random bin.

> contact_normal_force_XY[i]=contact_normal_force_XY[i]/contact_nXY[i]

Why do you count the numbers and use an average force for each bin? Isn't total force of the bin more meaningful?
Imagine you have 10 same forces. 9 forces in one bin and 1 forces is in another bin. With this approach you would get the same height for both bins. Is this what you intend?

Cheers
Jan

 Revision history for this message Zhicheng Gao (zhichenggao) said on 2021-09-03: #2

Dear Jan,

>What is your definition of "even distribution" and "uniform"?
By uniform, I mean that the height of the polar histogram is almost the same.

>any difference from utils.plotDirections? what difference? why?
Since the function utils.plotDirections can only display the contact distribution in a semicircle, I made the contact distribution to be displayed in the entire circle by symmetrical processing.

>Why do you count the numbers and use an average force for each bin? Isn't total force of the bin more meaningful?

I did that to calculate the contact normal force distribution defined by Rothenburg[1]

Usually, in the literature, the contact direction distribution and contact normal force distribution after isostatic consolidation are almost a circle. For example, like this picture： https://i.loli.net/2021/09/03/pwxedUM7cLyrGPK.jpg. This picture is from the article [2]
[2]https://www.sciencedirect.com/science/article/pii/S0266352X2030255X

I want to know if there is a problem in my isostatic consolidation process, which leads to uneven contact and contact normal force distribution. Thank you！

Best wishes

 Revision history for this message Jan Stránský (honzik) said on 2021-09-03: #3

> By uniform, I mean that the height of the polar histogram is almost the same.

just shift: how do you define "almost the same"?

> Since the function utils.plotDirections can only display the contact distribution in a semicircle, I made the contact distribution to be displayed in the entire circle by symmetrical processing.

since the other semicircle is symmetric, is it worth to change a "standard" function? what are the results of the original utils.plotDirections? same? better? worse?

> I did that to calculate the contact normal force distribution defined by Rothenburg[1]

I do not have access, can you provide the formula / method description?

> I want to know if there is a problem in my isostatic consolidation process, which leads to uneven contact and contact normal force distribution.

I would first suspect the plotting code. As discussed:
- the same forces distribution for the three planes is very suspicious
- the "randomness of bins of perpendicular forces" does not look very rigorous
- different computation for direction and forces, but it is intrinsically the same. For directions, are sphere-wall interactions included or not?
- I don't like the bin-wise averaging
- do you have some testing sample for which you know the result? Such test would strongly influence the code "confidence"
- ...

Cheers
Jan

 Revision history for this message Zhicheng Gao (zhichenggao) said on 2021-09-05: #4

>how do you define "almost the same"?
The difference is not more than 0.5 times the average value, but it can be seen from the figure that the maximum value is three times the average value

>since the other semicircle is symmetric, is it worth to change a "standard" function? what are the results of the original utils.plotDirections? same? better? worse?

The result of the original utils.plotDirections of course is the same

>I do not have access, can you provide the formula / method description?
Despite the apparent randomness in variation of contact forces, regular trends emerge when normal and tangential components of interparticle forces are averaged over groups of contacts with similar orientations.

 Revision history for this message Jan Stránský (honzik) said on 2021-09-06: #5

> Despite the apparent randomness in variation of contact forces, regular trends emerge when normal and tangential components of interparticle forces are averaged over groups of contacts with similar orientations.

ok. Just have in mind the case that in one bin you have one force 9 N, in another bin 3 forces 3 N. Depending on what you really want, plain sum or average gives (of course) totally different results

Anyway, in your case, you assign full force to the bins, no matter their relation to the plane.
So if the force is perpendicular to the plane, you still assign its full value, even if it sounds logical not to include it at all.
Consider assigning just the projection of the force to the plane.
Or just select "not too much perpendicular" interactions.
Or ...
There are really many options of the evaluation.

I suggest:
- using (trying at least) the same method for both forces and direction
- investigate why the 3 forces plot is exactly the same
- investigate if wall interactions were used for direction plot (they are used in "standard" utils.plotDirections)
- take into account force direction w.r.t. plane
- test the method on artificial input with known results

For further reasonable discussion, please provide a MWE (final packing and plot evaluation, no simulation).

Cheers
Jan

 Revision history for this message Zhicheng Gao (zhichenggao) said on 2021-09-12: #6

Dear Jan,
I am very grateful for your patient answers, and I am very sorry for not responding to you in time. I have taken some exams these days. I will try your suggestion.

Best wishes
Gao