# contact normal is not spherical

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

jsonscript (jsonscript) said : | #2 |

Hello，

Triaxial compression test with sand..

I meaurse the contact normal by i.geom.normal in O.interactions.

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:

PeriTriaxCo

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

jsonscript (jsonscript) said : | #4 |

Hello.

some codes are listed here:

O.engines=[

ForceResetter(),

InsertionSortC

InteractionLoop(

[Ig2_

[Ip2_

[Law2_

),

GlobalStiffnes

# PeriTriaxContro

NewtonIntegrat

PyRunner(

]

After consolidation, I meaurse the contact normal by i.geom.normal in O.interactions.

The rose diagram is listed here:

http://

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)

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-

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

cheers

Jan

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:/

[3] https:/

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.

sp=pack.

radius=5e-3

num=sp.

O.bodies.

O.engines=[

ForceResetter(),

InsertionSortC

InteractionLoop(

[Ig2_

[Ip2_

[Law2_

),

PeriTriaxContr

NewtonIntegrat

]

def triaxDone():

print ('Here we are: stress'

for i in O.interactions:

n=i.geom.normal

f=open(

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=PWaveTimeS

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:

# Function to test isotropy.

# For random direction counts normals "close" to this direction

def countSimilarNor

# 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 = [countSimilarNo

directions,counts = zip(*data)

# some statistics

cmax = max(counts)

import statistics

cmean = statistics.

cstddev = statistics.

print("Isotropy fulfilled with",cstddev/

# plot

dvals = [d*c/cmax for (d,c) in zip(directions,

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(

ax.auto_

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:

> for j=pi/18:

> A(hang,

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,

% 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