Yade

contact-normal in Yade

Asked by Fu zuoguang

Dear Prof. Chareyre and all users,

Fabric analysis is one of the most common methods for detecting those micro-level natures of such granular materials. Recently I have tried to do it in 2D-simulation task in Yade-1.14-daily version with the simplest isotropic condition. But unfortunately, I’m lost in getting the same results as that in Ref.[1]. All the details of my simulation are in Fig.1

Fig dir: http://tinypic.com/view.php?pic=2dwhukl&s=9#.V2kZvbEs-VU

After the consolidation process finished, I listed all the directional data I want to use as the following mode:

#################################################################
Body_Numbers, Body_Cur_P_x, Body_Cur_P_y, Body_Cur_P_z,
           1, -1.28536e-01, 0.00000e+00, 0.00000e+00,
           2, 1.28528e-01, 0.00000e+00, 0.00000e+00,
           3, 0.00000e+00, -1.28642e-01, 0.00000e+00,
           4, 0.00000e+00, 1.28632e-01, 0.00000e+00,
           5, 0.00000e+00, 0.00000e+00, -5.00000e-01,
           6, 0.00000e+00, 0.00000e+00, 5.00000e-01,
           7, 1.20924e-01, 6.60472e-02, 0.00000e+00,
           8, -1.22713e-01, 1.02525e-01, 0.00000e+00,
           9, 4.91498e-02, 1.84187e-02, 0.00000e+00,
          10, -4.19541e-02, -3.18193e-02, 0.00000e+00,
           .....
###########################################

So，I firstly summarized the contact-normal as the vector variable and transformed it from Cartesian to polar axis (using format 1-1) but got a strange result as shown in Fig.2. All the contact-normals are not symmetrically distributed around the origin of principle axes, but most of them are located in the area of (math.pi to 3/2math.pi).

Then, I respectively used the particle position and the contact-point position to execute the same process as above for ensuring that I got a correct output of those data, and plotted them in the same style as in Fig.3 and Fig.4. Now, the statistical results of the two group data show me the same trend. They are distributed axisymmetrically around the principle axes. And in the same quadrant, the number of data located in diagonal area is more than that in others, which can meet the deviatoric-shear condition.

Why is only the contact normal showing the non-symmetric distribution here?
[1] Are there anything wrong in my understanding of contact normal distribution?
[2] Or Yade has its own pattern in defining the contact normal?

Seeking your help!

Ref.
[1] Distribution of directional data and fabric tensors. 1984_Int.J.Eng.Sci

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Fu zuoguang
Solved:
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#1

>[1] Are there anything wrong in my understanding of contact normal distribution?

Looking at fig.2, your maths are wrong for sure. I can't tell where exactly.
Even figs.3 and 4 are weird, the extrema are like shifted by pi/4 from what is expected with vertical+horizontal walls.

>[2] Or Yade has its own pattern in defining the contact normal?

Not at all.

Bruno

Revision history for this message
Fu zuoguang (zgfu1985) said : 		#2

Dear Prof. Chareyre，

Follow your suggestion, I made a simple procedure (3th floor) with no strict controlling in its running for testing my problems again.

In this test, after the simulation process over, I did the following things:
(1) I listed all the contact normal data in O.interactions;
(2) I divided the total Cartesian system into 4 quadrants and then recorded respectively the number of contact normal located in each area with ignoring all the sphere-wall contacts.

From the final report, it is for sure that the unexpected results still exist now, maybe

(1) I exactly got a wrong output of contact normal data from Yade;
(2) The statistical pattern I used is already inappropriate;
(3) There is something incorrect in my understanding of contact normal and the concept of this analysis pattern is originally wrong.

Seeking your help!

Revision history for this message
Fu zuoguang (zgfu1985) said : 		#3

# encoding: utf-8
# the script demonstrates a simple case of triaxial simulation using TriaxialCompressionEngine. More elaborated examples can be found in the triax-tutorial folder
from yade import pack

sp=pack.SpherePack()
## corners of the initial packing
mn,mx=Vector3(0,0,0.5),Vector3(1,1,0.5)

## box between mn and mx, avg radius ± ½(20%), 2k spheres
sp.makeCloud(minCorner=mn,maxCorner=mx,rRelFuzz=.2,num=2000)

## create material #0, which will be used as default
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=radians(30),density=2600,label='spheres'))
O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=0,density=0,label='frictionless'))

## copy spheres from the packing into the scene
## use default material, don't care about that for now
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

## create walls around the packing
walls=aabbWalls(thickness=1e-10,material='frictionless')
wallIds=O.bodies.append(walls)

# block the dof in z-direction.
for b in O.bodies:
 if isinstance(b.shape, Sphere):b.state.blockedDOFs = 'zXY'

# define the engine.
triax=TriaxialStressController(
    wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],wall_left_id=wallIds[0],
 wall_right_id=wallIds[1],wall_back_id=wallIds[4],wall_front_id=wallIds[5],

    # keep all the walls from moving except the back and front walls.
    wall_back_activated = 0, wall_front_activated = 0,
    thickness = 1e-9, internalCompaction = 0,

    # stress controlling condition.
    stressMask = 7,
        goal1 = -100000,goal2 = -100000,goal3 = -0.0,

 label="triax", )

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 # you can add TriaxialStateRecorder and such here…
 NewtonIntegrator(damping=.4)
]

# run this simulation.
O.run(200000, True)
################################## running simulation process is over!

# record all the contact normals in the final state.
contact_normal_list=[[i.geom.normal[0],i.geom.normal[1]] for i in O.interactions]

# divide the whole Cartesian system into 4 quadrants as: (1)x>0 and y>0; (2)x<0 and y>0 ...
# then summarize respectively the number of contact normals located in each area with
# ignoring cursoryly the sphere-wall contacts.

quadrant_1_num, quadrant_2_num, quadrant_3_num, quadrant_4_num, =0, 0, 0, 0

for i in xrange(len(contact_normal_list)):
    if (contact_normal_list[i][0]>0 and contact_normal_list[i][1]>0): quadrant_1_num=quadrant_1_num+1
    if (contact_normal_list[i][0]<0 and contact_normal_list[i][1]>0): quadrant_2_num=quadrant_2_num+1
    if (contact_normal_list[i][0]<0 and contact_normal_list[i][1]<0): quadrant_3_num=quadrant_3_num+1
    if (contact_normal_list[i][0]>0 and contact_normal_list[i][1]<0): quadrant_4_num=quadrant_4_num+1
    else:pass

print 'the total number of contacts is: ', len(contact_normal_list)
print 'the number of sphere-sphere contacts is: ',(quadrant_1_num+quadrant_2_num+quadrant_3_num+quadrant_4_num)
print 'the number of contacts in 1 quadrant is: ', quadrant_1_num
print 'the number of contacts in 2 quadrant is: ', quadrant_2_num
print 'the number of contacts in 3 quadrant is: ', quadrant_3_num
print 'the number of contacts in 4 quadrant is: ', quadrant_4_num

Revision history for this message
Jan Stránský (honzik) said : 		#4

Hello,

there really must be something wrong in your calculations.

first of all, as I understand it, plotting normal is meaningful only for
values (0,pi), the other half may be mirrored if needed. From contact
direction point of view, vectors (1,0,0) and (-1,0,0) are identical..

with your script, I got perfect results.. either with
quadrant_X_num, quadrant_1_num+quadrant_3_num is almost exactly equal to
quadrant_2_num+quadrant_4_num (see above why there is the sum :-)

Also putting the directions into more "bins", I got e.g.
##############################################
def contactDirs(nBins=12):
normals = []
for i in O.interactions: # use only sphere-sphere contacts
if not isinstance(O.bodies[i.id1].shape,Sphere):
continue
if not isinstance(O.bodies[i.id2].shape,Sphere):
continue
normals.append(i.geom.normal)
# until now could be done in one line, but would look ugly in email
bins = [0 for i in range(nBins)]
for n in normals:
angle = atan2(n[1],n[0]) # atan2 returns value in range (-pi,pi)
if angle<0:
angle += pi
i = int(angle/pi * nBins)
bins[i] += 1
return bins
##############################################

Yade [1]: contactDirs(6)
 -> [1]: [475, 483, 484, 533, 486, 518]

Yade [2]: contactDirs(12)
 -> [2]: [268, 207, 264, 219, 249, 235, 262, 271, 255, 231, 235, 283]

which seems to be reasonable

cheers
Jan

2016-06-22 10:13 GMT+02:00 Fu zuoguang <<email address hidden>
>:

> Question #295521 on Yade changed:
> https://answers.launchpad.net/yade/+question/295521
>
> Fu zuoguang gave more information on the question:
> # encoding: utf-8
> # the script demonstrates a simple case of triaxial simulation using
> TriaxialCompressionEngine. More elaborated examples can be found in the
> triax-tutorial folder
> from yade import pack
>
> sp=pack.SpherePack()
> ## corners of the initial packing
> mn,mx=Vector3(0,0,0.5),Vector3(1,1,0.5)
>
> ## box between mn and mx, avg radius ± ½(20%), 2k spheres
> sp.makeCloud(minCorner=mn,maxCorner=mx,rRelFuzz=.2,num=2000)
>
> ## create material #0, which will be used as default
>
> O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=radians(30),density=2600,label='spheres'))
>
> O.materials.append(FrictMat(young=15e7,poisson=.4,frictionAngle=0,density=0,label='frictionless'))
>
> ## copy spheres from the packing into the scene
> ## use default material, don't care about that for now
> O.bodies.append([sphere(center,rad,material='spheres') for center,rad in
> sp])
>
> ## create walls around the packing
> walls=aabbWalls(thickness=1e-10,material='frictionless')
> wallIds=O.bodies.append(walls)
>
> # block the dof in z-direction.
> for b in O.bodies:
> if isinstance(b.shape, Sphere):b.state.blockedDOFs = 'zXY'
>
> # define the engine.
> triax=TriaxialStressController(
>
> wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],wall_left_id=wallIds[0],
>
> wall_right_id=wallIds[1],wall_back_id=wallIds[4],wall_front_id=wallIds[5],
>
> # keep all the walls from moving except the back and front walls.
> wall_back_activated = 0, wall_front_activated = 0,
> thickness = 1e-9, internalCompaction = 0,
>
> # stress controlling condition.
> stressMask = 7,
> goal1 = -100000,goal2 = -100000,goal3 = -0.0,
>
> label="triax", )
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_ScGeom_FrictPhys_CundallStrack()]
> ),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
> triax,
> # you can add TriaxialStateRecorder and such here…
> NewtonIntegrator(damping=.4)
> ]
>
> # run this simulation.
> O.run(200000, True)
> ################################## running simulation process is over!
>
> # record all the contact normals in the final state.
> contact_normal_list=[[i.geom.normal[0],i.geom.normal[1]] for i in
> O.interactions]
>
> # divide the whole Cartesian system into 4 quadrants as: (1)x>0 and y>0;
> (2)x<0 and y>0 ...
> # then summarize respectively the number of contact normals located in
> each area with
> # ignoring cursoryly the sphere-wall contacts.
>
> quadrant_1_num, quadrant_2_num, quadrant_3_num, quadrant_4_num, =0, 0,
> 0, 0
>
> for i in xrange(len(contact_normal_list)):
> if (contact_normal_list[i][0]>0 and contact_normal_list[i][1]>0):
> quadrant_1_num=quadrant_1_num+1
> if (contact_normal_list[i][0]<0 and contact_normal_list[i][1]>0):
> quadrant_2_num=quadrant_2_num+1
> if (contact_normal_list[i][0]<0 and contact_normal_list[i][1]<0):
> quadrant_3_num=quadrant_3_num+1
> if (contact_normal_list[i][0]>0 and contact_normal_list[i][1]<0):
> quadrant_4_num=quadrant_4_num+1
> else:pass
>
> print 'the total number of contacts is: ', len(contact_normal_list)
> print 'the number of sphere-sphere contacts is:
> ',(quadrant_1_num+quadrant_2_num+quadrant_3_num+quadrant_4_num)
> print 'the number of contacts in 1 quadrant is: ', quadrant_1_num
> print 'the number of contacts in 2 quadrant is: ', quadrant_2_num
> print 'the number of contacts in 3 quadrant is: ', quadrant_3_num
> print 'the number of contacts in 4 quadrant is: ', quadrant_4_num
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Fu zuoguang (zgfu1985) said : 		#5

Dear Jan，

Thanks for your suggestion and e.g. I am afraid that I had a wrong understanding of contact normal analysis and plotting.

As my original standpoint, all the contact normal data obtained by (i.geom.normal) should be uniformly located in the circle edge (Fig.1) in Cartesian system.

Fig.1 dir：http://tinypic.com/view.php?pic=66gfmf&s=9

They can be used to represent the unit vector with the start point O.
So, in my test codes uploaded yesterday, I used (++), (-+), (--) and (+-) to define respectively 4 quadrants and considered that the number of data in each group may almost the same in isotropic condition. But I was wrong, they are not the same at all. As reported by many tested cases in Yade, the number of data in 3rd quadrant is the largest and that in 4th quadrant is larger than the others. I think it is caused by Yade calculation system, and another DEM platform may have the same case.

Your opinion is very useful as the solution. You skillfully transformed all the vectors which are originally in 1st and 2nd quadrant, respectively to 3 and 4 quadrant BY (angle+pi). As shown in Fig.1, vector 1 is to be equivalent to 3 and 2 to 4. Thus, all the thetas are just in (0, pi), which can perfectly meet your point as

‘plotting normal is meaningful only for values (0,pi)’.

And your dealing with
'quadrant_1_num+quadrant_3_num= quadrant_2_num+quadrant_4_num' is reasonable.

I hope that I had no wrong in getting your points. But I had not seen all the details above in previous papers and simultaneously I am afraid that this method may get fabric tensor calculation into trouble.

Did you have some related papers which have recorded these details? If necessary, please tell me that.

cheers
zuoguang

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#6

>I am afraid that this method may get fabric tensor calculation into trouble.

Fabric tensor is based on squared normal, so the problem disappears: n and -n gives the same fabric.

Bruno

Revision history for this message
Fu zuoguang (zgfu1985) said : 		#7

Dear Prof. Chareyre and all users,

Thanks for Yade-2016, because the fabric tensor function in it now can be used for non-periodic cases. I tested this calculating process in two ways, one was using 'utils.fabricTensor()' directly and the other was using my own python codes just like

#########
def fabric_tensorcalculation(data_list):
        fabric_list = np.array([(i,i).dyadic_product() for i in data_list])
        total_fabric_tensor = sum(fabric_list) / float(len(fabric_list))
        return total_fabric_tensor
#########

They are almost the same.

my concern is exactly unnecessary, n and -n can provide the same single fabric tensor.