contact normal is not spherical

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

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!

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

> 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

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!

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

Hello,

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

Thanks!

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

> 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

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

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

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