et for ViscElMat with Law2_ScGeom_ViscElPhys_Basic, bug or feature?

Asked by Daniel Kiracofe on 2020-05-27

I had an issue with ViscElMat, started to write a question, ended up solving it myself. So now my question is really: is this a bug or a feature?

If you have a material like this
ViscElMat(young=200e9,poisson=0.7,frictionAngle=0, en = 0.5 , density=rho , label='ve')

and use InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], [Ip2_ViscElMat_ViscElMat_ViscElPhys()],[Law2_ScGeom_ViscElPhys_Basic()] ), in O.engines the resulting behavior is completely elastic in the normal direction. The coefficient of restitution is completely ignored. Same behavior regardless of the value of en. (A full working example is below).

Now if you do:

ViscElMat(young=200e9,poisson=0.7,frictionAngle=0, en = 0.5 , et=1.0, density=rho , label='ve')

Then you get the expected viscoelastic behavior, with differing response as en is changed. But you also get the exact same behavior with this

ViscElMat(young=200e9,poisson=0.7,frictionAngle=0, en = 0.5 , et=0.0, density=rho , label='ve')

or even this

ViscElMat(young=200e9,poisson=0.7,frictionAngle=0, en = 0.5 , et=1000.0, density=rho , label='ve')

The value of et is completely ignored, as long as its specified. So if the value is ignored, why force the user to specify it? Why not just ignore it all the time? Is this a bug, or was this done intentionally for some reason?

FWIW, The fact that et is ignored is mentioned in the documentation for Law2_ScGeom_ViscElPhys_Basic, but the fact that you get elastic behavior if et is omitted is not.

Using yade version 20200511-3819~5bf8512~buster1

################################################

import matplotlib.pyplot as pyplot
from yade import qt, plot
qt.View() #open the controlling and visualization interfaces

box_x = 0.05
box_y = 0.05
box_z = 0.05

particle_dia = 0.005

young2 = 200e9
rho = 8230

mn = Vector3(0, box_y, 0)
mx = Vector3(box_x,2*box_y, box_z) #corners of the initial packing
thick = 2*particle_dia # the thickness of the walls

global ballIds
#first create a very tall loose pack
# sp=pack.SpherePack()
bigmx = (mx[0], 3 * mx[1], mx[2])

restitution = 0.8

#this always gives elastic behavior. ball rebounds more or less exactly to starting height
O.materials.append(ViscElMat(young=young2,poisson=0.7,frictionAngle=0, en = restitution , density=rho , label='ve'))
#this give viscoelastic behavior. ball rebounds to 64% of initial height as expected. et is ignored, as long as it is specificed.
#O.materials.append(ViscElMat(young=young2,poisson=0.7,frictionAngle=0, en = restitution , et=1000, density=rho , label='ve'))

ball = sphere( (mx[0]/2, 2*mx[1] , mx[2]/2 ) , particle_dia/2, material='ve' )
ballIds = O.bodies.append(ball)

walls=utils.aabbWalls([mn,bigmx],thickness=thick,oversizeFactor=1.5,material='ve')
wallIds=O.bodies.append(walls)

#turn on gravity and let it settle
O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()] ),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
                [Law2_ScGeom_ViscElPhys_Basic()] ),
 NewtonIntegrator(gravity=(0,-9.81,0),damping=0.0),
        PyRunner(command='addPlotData()', iterPeriod=100),
]
O.dt=.05*PWaveTimeStep()

def addPlotData():
        plot.addData(time=O.time , pos = (O.bodies[ballIds].state.pos[1]-box_y)/0.15 )

plot.plots={'time':('pos')}
plot.plot()

#O.run(-1, True)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
2020-05-28
Last query:
2020-05-28
Last reply:
2020-05-28
Robert Caulk (rcaulk) said : #1

Hello,

Thanks a lot for the clear description and the minimal working example :-).

>FWIW, The fact that et is ignored is mentioned in the documentation for Law2_ScGeom_ViscElPhys_Basic, but the fact that you get elastic behavior if et is omitted is not.

Can you please link to where it says this? I cannot see it.

I admit I have not used this law but looking at the source code, et is used to compute cs and ks [1], which are both used to compute the torque [2].

Possible reasons you see no effect:
Your frictionAngle in your script is zero which might mean you never reach [3].
It is not immediately apparent that you would activate any shear in your example script.

Cheers,

Robert

[1]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticPM.cpp#L279
[2]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticPM.cpp#L66
[3]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticPM.cpp#L164

Daniel Kiracofe (kiracodl) said : #2

Robert, the documentation I'm referring to this here: https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Law2_ScGeom_ViscElPhys_Basic

There are four different cases. I'm trying to do the fourth one:

"If Yound modulus (young), poisson’s ratio (poisson), normal and tangential restitution coefficient (en,et) are precised, the equivalent stiffnesses are evaluated as previously: Kn=2kn1kn2kn1+kn2, knx=Exdx, Ks=2(ks1ks2)/(ks1+ks2), ksx=vknx. The damping constant is computed at each contact in order to fulfill the normal restitution coefficient en=(en1en2)/(en1+en2). This is achieved resolving numerically equation 21 of [Schwager2007] (There is in fact a mistake in the article from equation 18 to 19, so that there is a change in sign). Be careful in this configuration the tangential restitution coefficient is set to 1 (no tangential damping). This formulation imposes directly the normal restitution coefficient of the collisions instead of the damping constant."

The relevant line in the code is https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticPM.cpp#L263 and the following 3 lines.

if et and en are specified, cn is calculated from en, and cs is set to zero (so et is ignored). But if either et or en is omitted, then both cs and cn are not set to anything (so default zero).

so why not just change line 263 to be:

} else if (isfinite(mat1->en) ) {

> so why not just change line 263 to be:
> } else if (isfinite(mat1->en) ) {

I have no precise idea now (don't know that code in much details) but a suggestion is usually a good thing.
If I get it correctly it would make behavior in scenario 4 more consistent, would it change anything for cases 1-3?
Bruno

Daniel Kiracofe (kiracodl) said : #4

Bruno: I don't see any changes to cases 1-3. Those are handled by the other branches of the if-then-else tree starting here https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ViscoelasticPM.cpp#L215 Only the last case would be affected.

Daniel Kiracofe (kiracodl) said : #6

Yes exactly. Thank you!

Daniel Kiracofe (kiracodl) said : #7

Thanks Bruno Chareyre, that solved my question.