Elastic potential energy of polyhedra

Asked by Vasileios Angelidakis on 2019-09-06

Hi,

Following the discussion in [1], the normal stiffness "kn" in the Law2_PolyhedraGeom_PolyhedraPhys_Volumetric can have different units depending on the value of "volumePower"; e.g. for the default value volumePower=1, kn is in [N/m^3].

When calculating the elastic potential energy, we currently use [2,3,4]:
scene->energy->add(0.5*(normalForce.squaredNorm()/phys->kn+shearForce.squaredNorm()/phys->ks),
     "elastPotential",elastPotentialIx,true);
where "kn" should be in [N/m] for the units to give energy values [Joules=N*m].

I think the energy calculation in the normal direction should change either as 0.5*(normalForce.squaredNorm()/(phys->kn*phys->area),
where area the projection of the overlapping volume perpendicularly to the normal direction (this would work only if volumePower=1), or more safely, use the equivalentPenetrationDepth to begin with (which works regardless of the units of kn). Btw, looking at the calculation of the equivalentPenetrationDepth [5,6], is the exact "area" calculated or is an estimation made?

Cheers,
Vasileios

[1] https://answers.launchpad.net/yade/+question/681383
[2] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/Polyhedra.cpp#L517
[3] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/Polyhedra.cpp#L584
[4] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/Polyhedra.cpp#L589
[5] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/Polyhedra_Ig2.cpp#L121
[6] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/Polyhedra_Ig2.cpp#L219

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Vasileios Angelidakis
Solved:
2019-09-11
Last query:
2019-09-11
Last reply:
2019-09-06
Jan Stránský (honzik) said : #1

Hi Vasileios,

good point. The equations [2,3,4] are wrong. In general, elastic energy is
int F(u) du
in case F=k*u, int F du = int k*u du = 0.5*k*u^2 = 0.5*F^2/k
But it is not the case of polyhedrons interaction, where the function F(u) is actually unknown.
But probably some estimation could be done, e.g. F=k*V and assume V being proportional to u^3.
Or the approaches mentioned by you (today is too late for me to investigate them :-)

> Btw, looking at the calculation of the equivalentPenetrationDepth [5,6], is the exact "area" calculated or is an estimation made?

estimation [7]

cheers
Jan

[7] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/Polyhedra_Ig2.cpp#L116

Hi Jan,

Thanks for taking a look at this.
I think assuming that V is proportional to u^3 using a fixed estimation can be risky, since it depends highly on the geometry of the overlap region. :p

I like the idea of an equivalentPenetrationdepth, since it simplifies our problem in that the elastic potential work is F*equivalentPenetrationDepth (the integral of the force along the distance it acts on).

Currently, this distance is estimated by dividing the penetrationVolume by the area perpendicularly to the contact normal direction.
I think it would be more straightforward if we calculated two opposite points of the overlap volume, along the normal direction, passing from the contact point, and calculate the equivalentPenetrationDepth as the distance of these points.

Calculating the equivalentPenetrationDepth deterministically can give us an exact calculation for the "area" as well (=penetrationVolume/equivalentPenetrationDepth), in case you have any other uses for this parameter.

All the best,
Vasileios

Jan Stránský (honzik) said : #3

> I think assuming that V is proportional to u^3 using a fixed estimation can be risky, since it depends highly on the geometry of the overlap region. :p
> the function F(u) is actually unknown.

this is what I meant :-)

> two opposite points of the overlap volume, along the normal direction, passing from the contact point, and calculate the equivalentPenetrationDepth as the distance of these points.

then it would be needed to distinguish equivalentPenetrationDepth and this "estimatedPEnetrationDepth"

cheers
Jan

Good point, thanks Jan!

Regards,
Vasileios