yade-users team mailing list archive

[Xiaolong ZHAO <question692626@xxxxxxxxxxxxxxxxxxxxx>]

Question #692626 on Yade changed:
https://answers.launchpad.net/yade/+question/692626

    Status: Needs information => Open

Xiaolong ZHAO gave more information on the question:
Helle, Prof. Robert,
Sorry for my late reply and thank you for your response.  I  have checked the script and deleted all the unrelated items in engines.  Your question is good, and I have run it also for CohFrict and FrictMat. The same problem also exists for these two contact models. 
I think it is just a problem of  'pack.regularHexa' and it is not related to contact models.
Now, I can simplify the question.  I want to generate a spherical HCP packing with size L, and the elementary ball radius is R. If L/R is kept in the times of 2, such as 1/2, 2, 4, 8, 16, the fabric of the packing is strictly the same, but for other cases, such as L/R=3,5,6,7,..., the CN and number of interaction is slightly different.  And this difference will cause a pronounced difference in the following simulation.
I think it is caused by the error of float number of L and R. I once defined the value of R in a difference of 0.0001, which caused the CN  and number of interactions significantly different.
Below is the simplified script for JCFPM.
-------------------------------------------------------------
from yade import utils,pack,plot,timing
from math import *
import numpy as np
readParamsFromTable(LL=1.1e-2,RR=4e-4)   #LL=3.0e-2,RR=1e-3   LL=1.1e-2,RR=4e-4
enlarge_factor=1.0
young=3.3e9
poisson=.3
frictionAngle=radians(30)
tensileStrength=22.5e6
cohesion=225e6
gamma_int=1.0
from yade.params import table

rayon_s=float(table.RR*enlarge_factor)*1.0    #elementary sphere
rayon_g=float(table.LL/2.0*enlarge_factor)    #0.05  #large sphere
 
O.materials.append(JCFpmMat(type=1,young=young,poisson=poisson,frictionAngle=frictionAngle,density=2680,
                   tensileStrength=tensileStrength,cohesion=cohesion,
                   label='rock'))
sp=pack.regularHexa(yade._packPredicates.inSphere(center=(0,0,0),radius=rayon_g),gap=0,radius=rayon_s,material='rock')

O.bodies.append(sp)
count=0 
for i in O.bodies:
    if isinstance(O.bodies[i.id].shape,Sphere):  
        count=count+1
print('ball number=',count)

O.engines=[ForceResetter(),
           InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=gamma_int,label='is1aabb')]),
           InteractionLoop([Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=gamma_int,label='iss2sc'),
            ],
           [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1),            
            ],
           [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=False,smoothJoint=False)]),
            NewtonIntegrator(damping=0.8),
          ]
O.dt=0
O.step() #for CoheFrictM to release stress within agglomerate; for JCF model to create initial contacts
CN=avgNumInteractions(skipFree=True)
nbCohIntrs=0
for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
        if i.phys.isCohesive:
            nbCohIntrs+=1
print ("nbCohIntrs=",nbCohIntrs)
print('CN=',CN)
is1aabb.aabbEnlargeFactor=1.
iss2sc.interactionDetectionFactor=1.
O.dt=utils.PWaveTimeStep()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.