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.

_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to     : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Reply via email to