contact normal is not spherical

Asked by jsonscript on 2020-06-27

Hello everyone:

   The distribution of contact normal after the isotropic consolidation was found to be not spherical using the PeriTriaxController engine...
    May I ask what is the solution? Is there a problem in the process of preparing sample?

Thank you very much!!

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
jsonscript
Solved:
2020-07-05
Last query:
2020-07-05
Last reply:
2020-07-05
Jan Stránský (honzik) said : #1

Hello,
please provide a MWE [1] illustrating result, how you measured "contact normal" and what exactly is "spherical".
cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

jsonscript (jsonscript) said : #2

Hello,
   Triaxial compression test with sand..
   I meaurse the contact normal by i.geom.normal in O.interactions..then draw the rose diagram..
   The distribution diagram is not a sphere, there is little contact in the vertical direction, and the middle is recessed...

    Here is part of my code:
    PeriTriaxController(dynCell=True,maxUnbalanced=2e-3,relStressTol=2e-5,goal=[-1e5,-1e5,-1e5],stressMask=7,globUpdate=5,maxStrainRate=[0.1,0.1,0.1],doneHook='triaxDone()',label='triax'),

Thanks!

Jan Stránský (honzik) said : #3

Hello,
please provide a MWE [1] illustrating result, how you measured "contact normal".
To not-just-guess (out of maaaany possible reasons), we need your code.
cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

jsonscript (jsonscript) said : #4

Hello.
   some codes are listed here:
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=False)]
 ),

 GlobalStiffnessTimeStepper(timeStepUpdateInterval=100,label='ts',targetDt=1e-2),
 # PeriTriaxController(dynCell=True,maxUnbalanced=2e-3,relStressTol=2e-4,goal=[-1e5,-1e5,-1e5],stressMask=7,globUpdate=5,maxStrainRate=[1,1,1],doneHook='triaxDone()',label='triax'),
 NewtonIntegrator(damping=.2,label="newton"),
 PyRunner(iterPeriod=100,command='defData()'),
]

         After consolidation, I meaurse the contact normal by i.geom.normal in O.interactions..then draw the rose diagram..
        The rose diagram is listed here:
http://chuantu.xyz/t6/739/1593503924x-1404793492.jpg

Jan Stránský (honzik) said : #5

> some codes are listed here:

instant replay:
please read [1] and provide a MWE* illustrating result and how you measured "contact normal".
To not-just-guess (out of maaaany possible reasons), we need your code.
cheers
Jan

*M = minimal
*W = working (complete is a necessary condition)

[1] https://www.yade-dem.org/wiki/Howtoask

jsonscript (jsonscript) said : #6

hello.
   Sorry, I just want to know how to prepare uniform samples to finish triaxial test..
   Because my sample after consolidation is always anisotropic..

Thanks!

Jan Stránský (honzik) said : #7

Hello,

please read [1] and provide a MWE illustrating your result and how you measured "contact normal".
Is it a problem?

Of course it is possible to give you some code to prepare uniform samples, but it can be very different from your needs and therefore useless.
It is also possible that your sample is already isotropic and you just have some mistake in the isotropy-evaluation.

So a MWE provided by you is the most reasonable way to start a reasonable discussion..

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

jsonscript (jsonscript) said : #8

Hello,

Peritriaxcontroller we use is to produce isotropic stress state..
but they cannot guarantee us that we would get a contact normal isotropy.

Thanks!

Jan Stránský (honzik) said : #9

Yes, PeriTriaxController cannot guarantee a contact normal isotropy.

Which strongly depends on "contact normal isotropy" definition.
The contacts are random and discrete with finite number of contacts, so the isotropy in the strict sense is simply impossible.

On the other hand, the contact normal distribution should be "very close" to isotropy (depending on definition and evaluation).
Mechanical properties should be "very close" to isotropy, too.

If it is not like that in your case and you want a help, please provide a MWE.

cheers
Jan

Jan Stránský (honzik) said : #10

> I just want to know how to prepare uniform samples to finish triaxial test..

Out of maaaany options, see [2] figure 2.3 created by [3].

But as I said, it is possible (and very likely IMO) that your packing already is isotropic and there is a mistake in your evaluation.

cheers
Jan

[2] https://github.com/stranskyjan/phd-thesis/blob/master/pdfs/PhD_thesis_Stransky_2018_v2.pdf
[3] https://github.com/stranskyjan/phd-thesis/blob/master/codes/scripts/randomPeriPack/packingIsotropy.sh

jsonscript (jsonscript) said : #11

Hello Jan,

I realise I wasn't clear enough in my previous comments. I apologise for that.

Here is the MWE script (Runtime 2-3 mins) :

#-----------------------------------------------#
from yade import pack
from yade import plot

O.periodic=True
O.cell.setBox(.3,.3,.3)
sp=pack.SpherePack()
radius=5e-3
num=sp.makeCloud(minCorner=(0,0,0),maxCorner=(.3,.3,.3),rMean=radius,rRelFuzz=0.5,num=10000,periodic=True)

O.bodies.append([sphere(s[0],s[1]) for s in sp])

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PeriTriaxController(dynCell=True,maxUnbalanced=2e-3,relStressTol=2e-4,goal=[-1e5,-1e5,-1e5],stressMask=7,globUpdate=5,maxStrainRate=[1.,1.,1.],doneHook='triaxDone()',label='triax'),
 NewtonIntegrator(damping=.2)
]

def triaxDone():
  print ('Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff)
  for i in O.interactions:
    n=i.geom.normal
    f=open("normals.txt","a")
    print ("%s %s %0.9f %0.9f %0.9f" % (i.id1, i.id2, n[0], n[1], n[2]), file=f)
  f.close()
  O.pause()

O.dt=PWaveTimeStep()
O.run()
O.wait()
#----------------------------------------------------------#

I tried to consolidate a specimen isotropically with 100kPa and later checked for the isotropy in contact-normal distribution.
From the contact-normals recorded in the text file "normals.txt", I plotted the spherical histogram using Matlab. It seems to me, the distribution is not even near to isotropic.

I couldn't attach my histogram plot here, respecting "No-external link" rule of this forum. But, I see the distribution is not isotropic.

Am I doing something wrong here? if not, is there some other preparation method to generate isotropic sample w.r.t. both stress and contact normals?

Jan Stránský (honzik) said : #12

Hello,

thanks for the code. I have tried it and the results are IMHO nicely isotropic.
Evaluation:
I chose 200 random unit vectors and for these vectors I counted normals that are inclined from the vector maximally by given angle.
From these counts I computed relative standard deviation, whose value is:
1.1 % for angle 30 deg
3.2 % for angle 15 deg
11 % for angle 5 deg
Plotting the values gives relatively nice semi-sphere.
Of course it does not give the same number for every direction, but I did not find any systematic problem ("there is little contact in the vertical direction, and the middle is recessed")
###
# Extract normals from file
with open("normals.txt") as f:
   lines = f.readlines()
lines = [l.split()[2:] for l in lines]
lines = [[float(v) for v in l] for l in lines]
normals = [Vector3(*l) for l in lines]
print("Number of normals:",len(normals))

# Function to test isotropy.
# For random direction counts normals "close" to this direction
def countSimilarNormals(limAngle=pi/6):
   # random unit vector with positive z coordinate
   from random import gauss
   v = [gauss(0, 1) for i in range(3)]
   mag = sqrt(sum(pow(x,2) for x in v))
   v = Vector3(*[x/mag for x in v])
   if v[2] < 0: v *= -1
   # counting "similar" normals
   ret = 0
   for n in normals:
      # compute angle
      dot = v.dot(n)
      dot = abs(dot)
      angle = acos(dot)
      # if it is within the limit, count it
      if angle <= limAngle:
         ret += 1
   return v,ret

# count similar normals for a few random directions
data = [countSimilarNormals(pi/6/2/3) for i in range(200)]
directions,counts = zip(*data)
# some statistics
cmax = max(counts)
import statistics
cmean = statistics.mean(counts)
cstddev = statistics.stdev(counts)
print("Isotropy fulfilled with",cstddev/cmean*100,"% relative standard deviation")

# plot
dvals = [d*c/cmax for (d,c) in zip(directions,counts)]
x,y,z = [[d[i] for d in dvals] for i in (0,1,2)]
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes(projection="3d")
ax.auto_scale_xyz((-1,1),(-1,1),(0,1))
ax.scatter3D(x,y,z)
plt.show()
# sorry, I don't know how to set aspect ratio, so the result is displayed as semi-ellipsoid
###

> I plotted the spherical histogram using Matlab.
> Am I doing something wrong here?

please provide the code. As discussed, this stage may contain a mistake.

> I couldn't attach my histogram plot here

yes, thanks for respecting the rule, but you can add the data (i.e. a few lines of e.g. angle-number pairs) or the code producing the plot or the data.

cheers
Jan

Jan Stránský (honzik) said : #14

Thanks for the code

> for i=pi/18:2*(pi/36):(2*pi)
> for j=pi/18:(pi/18):(pi)
> A(hang,5)=length(find(PHI(:)>ic&PHI(:)<i&Theta(:)>jc&Theta(:)<j)); %Distribution to different bins

I don't think uniform separate distribution of theta and phi are correct here.
Such division has higher "area" around equator and lower "area" near poles (just have a look at the "globus earth model" with meridians and parallels).
This implies, that there is higher probability of normals to fall into equator bin than to the polar bin.
I.e., you should have bins with more normals near equator and less normals near poles.
Is this what you are observing?

cheers
Jan

jsonscript (jsonscript) said : #15

Hello Jan,

>Is this what you are observing?
Yes, your explanation clarifies my query. In my code, I just tried to extend the 2D Rose diagram idea to spherical distribution plot. I haven't thought about the area distribution factor. It makes sense now. I have been looking at the wrong distribution plot.

So, a sphere divided with equal areas, should give a more uniform distribution. I would try potting that.

Thank you so much. This answers my question.

Jan Stránský (honzik) said : #16

Just for future reference, this is exactly why I insisted on providing the code.

I think you can leave the distribution of the points, but change the condition a bit, something like:
%%% not tested, just an idea
ia = 0.5*(i+ic); % average i
ja = 0.5*(j+jc); % average j
d = pi/18; % or something like that
A(hang,5)=length(find(PHI(:)>ia-d&PHI(:)<ia+d&Theta(:)>ja-d&Theta(:)<ja+d));
% now the selecting region should be "square", equal for all directions
% I am not sure about the limit cases near 0, pi, 2*pi...
%%%

cheers
Jan