Adjusting the friction angle between dynamic components (spheres) and non-dynamic components (facets)

Asked by Paul Schröder

Hello dear YADE community,

I would like to load particles (spheres) in a simulation. These should have the following properties:
1. within the bulk of particles a certain friction value should be set in the form of an angle of repose (e.g. frictionAngle=0.65; i.e. approx. 37.5° for the material "Granodiorith").
2. furthermore, the friction coefficient μ ( approx. 0.75) between the particle material and the simulation surfaces made of facets should also be set.

I have already set the value for 1. via the "FrictMat" with "frictionAngle". For the 2nd value, I have made a conversion for the angle of repose according to [1] with the arctan and also set via the "FrictMat" with "frictionAngle". To check this friction value I wrote the script below. However, this shows me that I have not set the friction value required in 2.: μ = 0.087
In the script itself, a normal force is first applied to a particle by "Platte1" and "Platte2". Then "Platte3" moves towards the particle, which counteracts the frictional force between "Platte1", "Platte2" and the particle. The friction force is measured by the force reached by "plate3". The mean value of 10 of these frictional forces in relation to the normal force is then used to determine the friction value.

My main question is how to set the coefficient of friction for the facet material to 0.75. By just trying it out myself, I also come up with the following questions:
Have I reached this friction value requested by me and an error in this script?
Are there any other ways to set friction properties for the materials, or what does the command "setContactFriction" do if it does not affect the non-dynamic materials? This parameter was also mentioned in a similar question [2] in this context.
Does the type of calculation algorithm [3] have something to do with the fact that I get a different result? What must be done so that I still achieve the required setting?
How can I set a friction value to non-dynamic bodies like in 2.

With kind regards, Paul

[1] It is an external link which says: μ = tan(frictionAngle), if necessary just ask for link.
[2] https://answers.launchpad.net/yade/+question/700843
[3] https://yade-dem.org/doc/formulation.html?highlight=frictphys

########################################################################
# YADE-Version: 2022.01a
########################################################################
## Importieren zusaetzlicher Module
from yade import geom, pack, plot, qt
from math import pi, atan
from statistics import mean

## Definieren der verwendeten Parameter
dP = 10e-3 # Partikeldurchmesser

DIM = 2* dP
dPlt = 0.01*dP
Ueberlappung = 5e-2*dP # realisiert Anpresskraft
Platte3vel = 5e-3 # wirkt Reibkraft entgegen
Platte1vel = 5e-3
Platte2vel = -Platte1vel
Versuchszahl = 10
Normalkraft = 50

# Definieren der Materialparameter
EingabegroesseReibungStahl = 0.75
print("angestrebter Reibwert:", EingabegroesseReibungStahl)

idVersuchsgut = O.materials.append(FrictMat(young=70e9,poisson=.25, density=2700, frictionAngle=.65, label='Versuchsgut'))
idStempelPlatte = O.materials.append(FrictMat(young=.25e9,poisson=.28, density=7850, frictionAngle=atan(EingabegroesseReibungStahl), label='Platte'))

## Generieren der Simulationskomponenten
Simulationskoerper = 0

PartikelIDs = O.bodies.append(utils.sphere((0,0,0),radius=dP/2, material='Versuchsgut'))
Partikelmenge0 = Simulationskoerper
Partikelmenge1 = len(O.bodies)
Partikelmenge = len(O.bodies)-Simulationskoerper
Simulationskoerper += Partikelmenge

Platte1IDs = O.bodies.append(geom.facetBox((0,dP/2+dPlt+1e-2*dP,0), (DIM, dPlt, DIM), color=(0, 1, 0), wire=False, material='Platte'))
Platte1groesse0 = Simulationskoerper
Platte1groesse1 = len(O.bodies)
Platte1groesse = len(O.bodies)-Simulationskoerper
Simulationskoerper += Platte1groesse

Platte2IDs = O.bodies.append(geom.facetBox((0,-dP/2-dPlt-1e-2*dP,0), (DIM, dPlt, DIM), color=(0, 1, 0), wire=False, material='Platte'))
Platte2groesse0 = Simulationskoerper
Platte2groesse1 = len(O.bodies)
Platte2groesse = len(O.bodies)-Simulationskoerper
Simulationskoerper += Platte2groesse

Platte3IDs = O.bodies.append(geom.facetBox((-dP/2-dPlt-5e-3*dP,0,0), (dPlt, dP, DIM), color=(0, 1, 0), wire=False, material='Platte'))
Platte3groesse0 = Simulationskoerper
Platte3groesse1 = len(O.bodies)
Platte3groesse = len(O.bodies)-Simulationskoerper
Simulationskoerper += Platte3groesse

## Definieren des Ausgangszustandes von Laufvariablen
Platte1kraft = 0
Platte2kraft = 0
Platte3kraft = 0
Partikelkraft = 0
maxPartikelvel = 0
Haftreibwert = []
Bewegung = False
Iteration = False

Phase0 = "Anfang"
Phase01 = "Normalkraft"
Phase1 = "Annaeherung"
Phase2 = "Entfernung"
Phase3 = "Ende"
Iterationsphase = Phase0

## Definieren der Engines
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(),]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()],),
    NewtonIntegrator(damping=.3),
    PyRunner(command='Statemaschine()',realPeriod=1),]
O.dt=0.2*PWaveTimeStep()

# Steuereung des Simulationsablaufes: Statemaschine
def Statemaschine():
    global Iterationsphase, Platte1kraft, Platte2kraft, Platte3kraft, Partikelkraft, maxPartikelvel, Bewegung, Iteration, Haftreibwert

    aktPlatte1kraft = 0
    for i in range(Platte1groesse0,Platte1groesse1):
        aktPlatte1kraft += O.forces.f(i)[1]
    deltaPlatte1kraft = aktPlatte1kraft - Platte1kraft
    Platte1kraft = aktPlatte1kraft
    #print("Platte1kraft",Platte1kraft)

    aktPlatte2kraft = 0
    for i in range(Platte2groesse0,Platte2groesse1):
        aktPlatte2kraft += O.forces.f(i)[1]
    deltaPlatte2kraft = aktPlatte2kraft - Platte2kraft
    Platte2kraft = aktPlatte2kraft
    #print("Platte2kraft",Platte2kraft)

    aktPlatte3kraft = 0
    for i in range(Platte3groesse0,Platte3groesse1):
        aktPlatte3kraft += O.forces.f(i)[0]
    deltaPlatte3kraft = aktPlatte3kraft - Platte3kraft
    Platte3kraft = aktPlatte3kraft
    #print("Platte3kraft",Platte3kraft)

    aktPartikelkraft = 0
    for i in range(Partikelmenge0,Partikelmenge1):
        aktPartikelkraft += O.forces.f(i)[0]
    deltaPartikelkraft = aktPartikelkraft - Partikelkraft
    Partikelkraft = aktPartikelkraft
    #print("Partikelkraft",Partikelkraft)

    maxPartikelvel = 0
    for i in range(Partikelmenge0,Partikelmenge1):
        Partikelvel = O.bodies[i].state.vel[0]
        if Partikelvel > maxPartikelvel:
            maxPartikelvel = Partikelvel
    #print("maxPartikelvel",maxPartikelvel)
    #print("len(Haftreibwert)",len(Haftreibwert))

    Krit0 = (abs(Platte3kraft)>0)
    Krit1 = (abs(Partikelkraft)>1e-10)
    Krit2 = (abs(maxPartikelvel)>1e-10)
    Krit3 = (len(Haftreibwert)>=Versuchszahl-1)
    Krit4 = ((abs(Platte1kraft-deltaPlatte1kraft)+abs(Platte2kraft-deltaPlatte2kraft))/2>Normalkraft)

    #print("Krit0",Krit0)
    #print("Krit1",Krit1)
    #print("Krit2",Krit2)
    #print("Krit3",Krit3)
    #print("Iterationsphase",Iterationsphase)

    if Iterationsphase == Phase0:
        Iteration = True
        Iterationsphase = Phase01
        print("Partikel in Ruhe")
        O.engines += [TranslationEngine(translationAxis=[0, 1, 0], velocity=-Platte1vel, ids=Platte1IDs)]
        O.engines += [TranslationEngine(translationAxis=[0, 1, 0], velocity=-Platte2vel, ids=Platte2IDs)]
    elif Iterationsphase == Phase01:
        if Krit4:
            Iterationsphase = Phase1
            print("Normalkraft:", round((abs(Platte1kraft)+abs(Platte2kraft))/2,2))
            O.engines += [TranslationEngine(translationAxis=[0, 1, 0], velocity=0, ids=Platte1IDs)]
            O.engines += [TranslationEngine(translationAxis=[0, 1, 0], velocity=0, ids=Platte2IDs)]
            O.engines += [TranslationEngine(translationAxis=[1, 0, 0], velocity=Platte3vel, ids=Platte3IDs)]
    elif Iterationsphase == Phase1:
        if Krit2:
            Iterationsphase = Phase2
            #print("Tagentialkraft1:", round(abs(Platte3kraft),5))
            #print("Tagentialkraft2:", round(abs(Platte3kraft-deltaPlatte3kraft),5))
            Wert = 2*abs(Platte3kraft)/(abs(Platte1kraft)+abs(Platte2kraft))
            Haftreibwert.append(Wert)
            #print(abs(Platte1kraft-deltaPlatte1kraft)+abs(Platte2kraft-deltaPlatte2kraft))
            print("Haftreibwert", Wert)
            O.engines += [TranslationEngine(translationAxis=[1, 0, 0], velocity=0, ids=Platte3IDs)]
            if Krit3:
                Iterationsphase = Phase3
            else:
                print("Partikel in Bewegung -> Entlastung")
                O.engines += [TranslationEngine(translationAxis=[1, 0, 0], velocity=-Platte3vel, ids=Platte3IDs)]
    elif Iterationsphase == Phase2:
        if not(Krit1 or Krit2 or Krit0):
            Iterationsphase = Phase1
            print("Partikel in Ruhe -> Belastung")
            O.engines += [TranslationEngine(translationAxis=[1, 0, 0], velocity=Platte3vel, ids=Platte3IDs)]
    elif Iterationsphase == Phase3:
        print("Mittelwert aus",Versuchszahl,"Messwerten ergibt:",mean(Haftreibwert))
        O.engines += [TranslationEngine(translationAxis=[1, 0, 0], velocity=0, ids=Platte3IDs)]
        O.pause()
        Iterationsphase = Phase2

    plot.addData(t_in_s=O.time,Platte1kraft_in_N=Platte1kraft,Platte2kraft_in_N=Platte2kraft,Platte3kraft_in_N=Platte3kraft,Partikelkraft_in_N=Partikelkraft,groesste_Partikelgeschwindigkeit_in_m_je_s=maxPartikelvel)

    return Iterationsphase, Platte1kraft, Platte2kraft, Platte3kraft, Partikelkraft, maxPartikelvel, Bewegung, Iteration, Haftreibwert

## Data mining
plot.plots={'t_in_s': (('Platte1kraft_in_N','r-'),('Platte3kraft_in_N','c--'),('Partikelkraft_in_N','b--'), None, ('groesste_Partikelgeschwindigkeit_in_m_je_s','g--'))}
plot.plot(subPlots = False)

qt.Controller()
qt.View()
O.run()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Best Jan Stránský (honzik) said :
#1

Hello,

> frictionAngle=0.65; i.e. approx. 37.5°
> the friction coefficient μ ( approx. 0.75)

μ=0.75 means frictionAngle=atan(0.75)=36.
Almost the same as case 1.
So the question is if to bother with two values at all :-)

> To check this friction value I wrote the script.

There are two aspects.
1) What value is actually assign
2) How it behaves in actual simulation

Your script seems to do the option 2, where many factors may influence the result (e.g. dynamic behavior of the sphere).
I suggest first the method 1, i.e. take the interaction and check the actual value, i.e.:
###
intr = O.interactions[0,1] # for example bodies 0 and 1
tanfa = intr.phys.tangensOfFrictionAngle # [4]
###

> My main question is how to set the coefficient of friction for the facet material to 0.75.

As you do in the script

How to apply it to different material combinations, see below.

> Have I reached this friction value requested by me and an error in this script?

As mentioned above, the "real" simulation may influence the results.
For friction measurement, I suggest making all particles with prescribed motion such that forces are just measured but do not influence the particles.

> ...

In Yade design, material does not care on what bodies or shapes it is assigned. So there is nothing like "non-dynamic material"

setContactFriction [5] is an auxiliary function to change friction angle on all
The documentation sentence (and implementation) "The friction for non-dynamic bodies is not modified" is a bit misleading, it actually may or may not be modified, according to how the materials are shared among particles.

To use different values of certain quantity for different material combinations, specify MatchMaker for frictAngle parameter of Ip2_FrictMat_FrictMat_FrictPhys

Cheers
Jan

[4] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FrictPhys.tangensOfFrictionAngle
[5] https://yade-dem.org/doc/yade.utils.html#yade._utils.setContactFriction
[6] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Ip2_FrictMat_FrictMat_FrictPhys.frictAngle

Revision history for this message
Paul Schröder (paul5678) said :
#2

Thank you for your answer Jan,
as far as I can see now, the core question is answered with it. However, I would still be interested in the further, whether I can also make a distinction in static and sliding friction? Is such a distinction possible with a standard material model?

If someone would like to learn from this question in the future, I summarize my knowledge here:
> 2) How it behaves in actual simulation
In order to examine the particle contact in the script belonging to the question, I have inserted the code suggested by Jan accordingly at my example within the function "Statemaschine()" after the definition of "Krit4":
###
    if Krit4:
        intr = O.interactions[Partikelmenge0+0,Platte1groesse0+4]
        tanfa = intr.phys.tangensOfFrictionAngle
        print("actually assigned value:", tanfa)
###

Then I also noticed that further settings the smaller friction value of both contact partners is taken over. To change that, I have applied the suggested "MatchMaker" according to [7]. The appropriate formulation, which is supplemented into the definition of the engines, I adapted for it as follows. I have added here the friction value leading steel-steel, even if this is not needed in my simulation.
###
        [Ip2_FrictMat_FrictMat_FrictPhys(frictAngle=MatchMaker(matches=((idVersuchsgut, idStempelPlatte, atan(0.75)),(idVersuchsgut, idVersuchsgut, .65),(idStempelPlatte, idStempelPlatte, atan(0.12)))))],
###

[7] https://answers.launchpad.net/yade/+question/244498

Revision history for this message
Paul Schröder (paul5678) said :
#3
Revision history for this message
Jan Stránský (honzik) said :
#4

> whether I can also make a distinction in static and sliding friction?

in general (e.g. implementing it from scratch) yes.

> Is such a distinction possible with a standard material model?

I am not sure, maybe you can open a new specific question o this topic (as this might be "deep" in the thread for people able to answer to actually answer).

Cheers
Jan

Revision history for this message
Paul Schröder (paul5678) said :
#5

Thank you for your feedback, Jan. I'll be happy to follow your advice and open a new question on the subject [8]. You have already been of good help to me.
With best wishes, Paul

[8] https://answers.launchpad.net/yade/+question/707430