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