yade-users team mailing list archive

[Paul Schröder <question707367@xxxxxxxxxxxxxxxxxxxxx>]

New question #707367 on Yade:
https://answers.launchpad.net/yade/+question/707367

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()


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