How to define appropriate k, r and R values for irregular polyhedra in Potential Particles?

Asked by weijie on 2020-11-05

Dear all,

When I use Potential Particles, I don’t know how to define the appropriate k, r, and R values[1] for the polyhedron to display its image correctly.

Thanks in advance.
Jie

[1]https://gitlab.com/yade-dev/trunk/-/blob/master/examples/PotentialParticles/cubePPscaled.py#L89

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Vasileios Angelidakis
Solved:
2020-11-12
Last query:
2020-11-12
Last reply:
2020-11-10

This question was reopened

Hi Jie,

Just to make sure, are you facing a purely visualisation issue, e.g. the particle surface looks grainy/step-like, or are you asking which k,r,R values might work for your particle from a modelling perspective?

If it is the former, please have another look in the documentation of the Potential Particles [1-2].
If it is the latter, please provide a MWE and I will try to help.

All the best,
Vasileios

[1] http://yade-dem.org/doc/potentialparticles.html#visualization
[2] http://yade-dem.org/doc/potentialparticles.html#axis-aligned-bounding-box

weijie (amandajoe) said : #2

Hi Vasileios, and thank you again.

I consider which k, r, R value is more suitable for particles from the perspective of modeling.This is my MWE[1].I found that the particles created with this are very different from the particles of Potential Blocks with the same parameters, and it seems that the volume and Principal inertia cannot be directly calculated with b.shape.volume, b.shape.inertia.

Best regards,
Jie

[1]
distanceToCentre= 0.025
b=Body()
b.mask=1
b.aspherical=True
wire=False
color=Vector3(random.random(),random.random(),random.random())
highlight=False
aa= [-0.7295950392572228, 0.46966065591120043, 0.978755366731517, 0.4159312753165438, -0.7851930775036051, 0.32622712056941044, 0.2487629499110137, -0.32802297814039366, 0.22757387395233566, -0.34064170979743924, 0.6535371560905189, -0.534030672328059, -0.9483691932498606]
bb= [-0.5797330029669138, -0.8677701055770507, 0.0908583947065325, -0.4698271127991969, 0.37706716796452805, 0.04463765109291187, 0.9516154740960145, -0.9017571077715157, 0.6925524728306885, -0.10267768395748879, 0.3532767701547774, 0.8086775588116711, -0.0540657356796434]
cc= [-0.3627681407761932, -0.1624620329672607, -0.18380066432307973, 0.7786293458972144, -0.49121500575982796, -0.9442369119611344, 0.18040228438849618, 0.28148720112193015, -0.6845299148104668, 0.9345696971138125, 0.6693912975817099, 0.24668167116191958, 0.3125264301143626]
dd= [0.01697034586028, 0.01248528317756, 0.01173016874627, 0.01854063787194, 0.01405551559515, 0.01349450583878 , 0.01047944143112, 0.01105292067252 , 0.02109918761652, 0.01609025523102 , 0.01787444301714,0.01542370609022 , 0.01180116564426 ]
r=min(dd)/2
b.shape=PotentialParticle(k=0.1, r=r, R=0.5*R, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=color, wire=wire, highlight=highlight, minAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), minAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), AabbMinMax=True, id=count)
V=(2*distanceToCentre)**3
geomInert=(1./6.)*V*(2*distanceToCentre)**2
utils._commonBodySetup(b, V, Vector3(geomInert,geomInert,geomInert), material='frictionless', pos=[0,0,0], fixed=False)
b.state.pos =[0,0,0]
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random())
O.bodies.append(b)

Hi Jie,

Indeed, the volume, inertia and initial Aabb are not calculated automatically for the Potential Particles. You have to do something manually. E.g. if your particles are fairly angular, you can create first a Potential Block with the same a,b,c,d and use their inertial characteristics for the PP. (This is a hack, but if the PP is kind of similar to the PB, it can lead to a somewhat realistic representation of its volume/inertia/bounds etc).

Also, the principal orientations are not calculated manually, so a PP has to be defined to its principal inertial system atm, to get correct angular momenta and rotations.

Adding to this, since the PPs are represented by a smooth implicit function, their visualisation is subjected to resolution issues, based on the sampling density (using the sizeX, sizeY, sizeZ parameters in the Gl functor) and the size of the drawing space considered (minAabb, maxAabb).

Regarding the parameters k,r,R, they control how rounded the edges and faces of the particles will be. Having a look in Eq. 2 [1], you can see that "k" controls a spherical term, determining how curved the faces are. Also, since the first term of Eq. 2 is normalised by dividing with "r", while the second term by dividing with "R", the ratio of r/R affects the overall shape.

In general, small values of k and r result in more angular shapes, although we do not suggest very small values for these, since in such cases we might have a bad calculation of contact normals. Please find a script below, showing a PP and a PB using your a,b,c,d. To have an empirical automatic rule, you can try to associate r and R with the largest "d" values somehow, e.g. R=max(dd) and r=0.2*R.

Best,
Vasileios

[1] https://yade-dem.org/doc/potentialparticles.html#equation-ppformulanormalized

distanceToCentre= 0.025
R=distanceToCentre/2.

O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(0.0),density=2000,label='frictionless'))
count=0

# Potential Particle
b1=Body()
b1.aspherical=True
aa= [-0.7295950392572228, 0.46966065591120043, 0.978755366731517, 0.4159312753165438, -0.7851930775036051, 0.32622712056941044, 0.2487629499110137, -0.32802297814039366, 0.22757387395233566, -0.34064170979743924, 0.6535371560905189, -0.534030672328059, -0.9483691932498606]
bb= [-0.5797330029669138, -0.8677701055770507, 0.0908583947065325, -0.4698271127991969, 0.37706716796452805, 0.04463765109291187, 0.9516154740960145, -0.9017571077715157, 0.6925524728306885, -0.10267768395748879, 0.3532767701547774, 0.8086775588116711, -0.0540657356796434]
cc= [-0.3627681407761932, -0.1624620329672607, -0.18380066432307973, 0.7786293458972144, -0.49121500575982796, -0.9442369119611344, 0.18040228438849618, 0.28148720112193015, -0.6845299148104668, 0.9345696971138125, 0.6693912975817099, 0.24668167116191958, 0.3125264301143626]
dd= [0.01697034586028, 0.01248528317756, 0.01173016874627, 0.01854063787194, 0.01405551559515, 0.01349450583878 , 0.01047944143112, 0.01105292067252 , 0.02109918761652, 0.01609025523102 , 0.01787444301714,0.01542370609022 , 0.01180116564426 ]
r=min(dd)/5.
b1.shape=PotentialParticle(k=0.1, r=r, R=0.5*R, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=[0,1,0], wire=False, minAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabb=sqrt(3)*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), maxAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), minAabbRotated=1.02*Vector3(distanceToCentre,distanceToCentre,distanceToCentre), AabbMinMax=True, id=count)
V=(2*distanceToCentre)**3 # FIXME: Not true
geomInert=(1./6.)*V*(2*distanceToCentre)**2 # FIXME: Not true
utils._commonBodySetup(b1, V, Vector3(geomInert,geomInert,geomInert), material='frictionless', pos=[0,0,0], fixed=False)
b1.state.pos =[0,0,0]
#b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random())
O.bodies.append(b1)

# Potential Block with same a,b,c,d
b2=Body()
b2.aspherical=True
b2.shape=PotentialBlock(r=r, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=[1,0,0], wire=True, AabbMinMax=True)
utils._commonBodySetup(b2, b2.shape.volume, b2.shape.inertia, material='frictionless', pos=[0,0,0], fixed=False)
b2.state.pos =[0,0,0]
b2.state.ori = b2.shape.orientation
O.bodies.append(b2)

from yade import qt
v=qt.View()
v.sceneRadius=0.25

v.eyePosition=Vector3(0.06420229007600809779,0.07990533579869238401,0.08296260554396843456)
v.upVector=Vector3(-0.240421083299934335,0.7853946365571540245,-0.5703972015816188845)
v.viewDir=Vector3(-0.4868619564481135864,-0.6059420633056463723,-0.6291261012550753984)

# I use the values below to achieve a good visualisation resolution in Qt.
Gl1_PotentialParticle.sizeX=100
Gl1_PotentialParticle.sizeY=100
Gl1_PotentialParticle.sizeZ=100

weijie (amandajoe) said : #4

Thanks Vasileios Angelidakis, that solved my question.