Bevezetés

Kezdetben a pályaszerkezeteket szinte a forgalom nagyságától és a talaj teherbírásától függetlenül, közel azonos anyagból és azonos vastagságban építették meg. Később a forgalom rohamos növekedésének hatására szükségessé vált, hogy az utak pályaszerkezetének vastagsága mind a forgalom nagyságával, mind a földmű minőségével arányos legyen. Ebben az időszakban dolgozták ki a gyakorlati megfigyeléseken, nagyminta kísérleteken és elméleti elgondolásokon alapuló empirikus és szemiempirikus méretezési módszereket.

Az igények és az igénybevételek növekedése miatt csakhamar előtérbe kerültek a mechanikai alapokon nyugvó eljárások is. Ezt támasztja alá az Európai Unió által finanszírozott AMADEUS kutatás 1998-1999-ben folyt zárójelentése, mely kimutatta, hogy bár Európában – és világszerte másutt is – többféle megközelítéssel élnek az útburkolatok méretezésére, alapvető közös jellemzőjük mégis az, hogy a pályaszerkezetben a tengelysúlyok áthaladásának hatására keletkező erők, feszültségek és alakváltozások számításán alapulnak. Ezek a mechanikai alapú méretezési módszerek a pályaszerkezetet rugalmasságtani alapon számítható szerkezetként fogják fel. A réteges szerkezetek méretezésére kidolgozott és alkalmazott mechanikai eljárások:

  1. az egyenérték vastagságon (Odemark),
  2. a rétegenkénti analitikus számításokon,
  3. a véges elemek módszerén (VEM),
  4. és a diszkrét elemek módszerén (DEM) alapulnak.

Az egyenérték-vastagság elvén alapuló módszerek estében a réteges szerkezeteket egy végtelen féltérré alakítják át, melyre már érvényesek és alkalmazhatók lesznek a Boussinesq-féle állapotegyenletek. A rétegenkénti analitikus modellek általában Burmister elgondolásain alapulnak, melynek segítségével a pályaszerkezet bármely pontjában számíthatóvá válnak „közelítően” az ébredő erők, feszültségek és alakváltozások.

Jelenleg a mérnöki gyakorlatban leggyakrabban használt numerikus eljárás, a végeselem-módszer (Finite Element Method, FEM). A végeselem-módszer egy numerikus eljárás parciális differenciálegyenletek közelítő megoldására. A véges elemek módszerének lényege a közelítő eljárásoknál a geometriai és a matematikai függvénytér finitizálásával együtt járó sajátos bázisfüggvény-megválasztási technika [Bojtár és Gáspár; 2003]. A keresett függvény leírása a tér kisebb, már matematikailag jól kezelhető részekre bontásával elemi függvényekből összeállítható. Az elemi függvények értelmezési tartományai nem fednek át és összegük a globális függvény értelmezési tartományát adja. A modellezéshez leginkább használatos végeselem a három kontrollpont alkotta térbeli háromszög. A függvényt e háromszögön belül kell bizonyos peremfeltételek mellett megkomponálni [Czimber; 2002].

A numerikus eljárások között ma még inkább csak kísérleti jelleggel használt módszer a szerkezetet különálló elemekből felépítő, úgynevezett diszkrét elemes módszer (DEM). Ez a vizsgálat alapjaiban különbözik az imént ismertetett FEM technikától. Az anyagot itt nem kontinuumnak tételezzük fel, hanem apró szemcsék együtteseként modellezzük. A makroszintű anyagi jellemzőket az egyes szemcsék közti kapcsolatok mechanikai jellemzői, illetve a halmaz geometriája szabják meg. A külső mechanikai hatások eredményeképpen az egyes szemcsék eltolódhatnak, illetve elfordulhatnak. Az elmozdulások Newton II. törvénye alapján számíthatók. Ennél a modellezési technikánál a megoldás értelmezése is igen összetett feladat, hiszen az eredményt nem a klasszikus elméletekben megszokott folytonos függvények formájában kapjuk meg, hanem az egyes szemcsék elmozdulásait és közöttük átadódó erőket ismerjük meg csupán. A fő eltérés a FEM és a DEM között az, hogy míg a végeselem-módszer, makrójellemzők alkalmazásával kontinuumnak (esetünkben homogénnek) tekinti a tartomány viselkedését, addig a diszkrételem-módszer alulról építkezik, a vizsgált tartományt különálló szemcsék együtteseként kezeli [Nasztanovics és Füstös; 2000].

Most következő cikkünkben a hajlékony útpályaszerkezek DEM alapú modellezésének lehetőségeit kívánjuk bemutatni. Az elméleti háttér ismertetésekor, tanulmányunk nagymértékben támaszkodik Bagi [2008] munkájára, mely mindez idáig az egyetlen magyar nyelven megjelent összefoglaló írás a témáról.

A diszkrét elemes modell alapjai

A diszkrét elemes modell (DEM) különálló elemekből és az elemek érintkezésével létrejövő kapcsolatokból áll. Bagi [2008] szerint egy numerikus eljárást – némileg pontosítva Cundall és Hart 1992-es definícióját – akkor tekinthetünk diszkrét elemes modellnek, ha egymástól egyértelműen elkülönülő elemekből épül fel; az elemek önálló elmozdulási szabadságfokokkal rendelkeznek oly módon, hogy a modell képes követni az elemek véges nagyságú eltolódásait és elfordulásait (esetleg deformációit is); az elemek között új kapcsolatok jöhetnek létre és meglévő kapcsolatok szűnhetnek meg. Végtelen merev elemekkel dolgozó diszkrét elemes modell futása hat fő szakaszra bontható [Jakob és Konietzky; 2012]:

1. A geometriai modell megalkotása: a kezdeti és a határfeltételek rögzítése, részecskék anyagtulajdonságainak beállítása és a szemcsehalmaz generálása (ismert az elemek kezdeti helyzete).

2. A részecskék és a határfelületek között kialakuló kapcsolatok meghatározása.

3. A kapcsolatoknál fellépő erők (\(F\)) és nyomatékok (\(M\)) számítása.

4. A részecskék gyorsulás (\(\ddot{u}\), \(\dot{\omega}\)), sebesség (\(\dot{u}\),\(\omega\)) és végül az elmozdulás növekmények (\(u\)) és elfordulások (\(\varphi\)) meghatározása.

5. Az összes részecske új helyzetének (pozíciójának) számítása.

6. A 2-es és 5-ös pont közötti fázisokat adott időközönként (időlépésenként) addig ismételjük, míg el nem érjük a leállási feltételt.

A számítás lényege, hogy a rendszer adott tehernövekményből az elmozdulásnövekményeket kis lépések sorozataként határozza meg. A diszkrét elemes modell a terhelési folyamatot tehát kis elmozdulás-növekmények sorozataként állítják elő.

Fontos megemlítenünk, hogy Jakob és Konietzky [2012] megfogalmazása nem általában a diszkrételemes modellekre, hanem azok egy speciális fajtájára, a végtelen merev elemeket alkalmazó modellekre vonatkozik. Deformálható elemek esetén másképp működnek a modellek, ezekkel azonban cikkünkben nem foglalkozunk.

A bemutatott eljárás úttörője Peter A. Cundall aki az Imperial College doktoranduszaként a 70-es évek legelején mutatta be első numerikus modelljét [Bagi; 2008]. A továbbiakban az imént bemutatott eljárást részletiben mutatjuk be.

A geometriai modell felvétele

A diszkrét elemes modellezés első lépése a rendszer geometriai modelljének elkészítése. A rendszert alkotó részecskéket a legtöbb esetben merev testként kezeljük, deformációt nem szenvednek. Ilyenkor az elemek közötti kapcsolatok viselkedésébe „sűrítjük” az anyagi tulajdonságokat. A másik lehetőség az, hogy az elemek deformálhatóak, meg kell adni a feszültség-alakváltozás összefüggést. A leggyakoribb modellezési eljárás az, hogy az elemek egyszerű belső végeselemes hálóval rendelkeznek, így a rájuk ható külső erőkből számíthatók az elemek deformációi. A részecskék alakja elvileg bármilyen lehet, de a számítások szempontjából a gömb a leghatékonyabb és legpraktikusabb választás. Természetesen léteznek olyan eljárások ahol a részecskéket poligon, ellipszis vagy éppen körökből összetett elemtípussal (clump) reprezentálják.

A geometriai modell felállításának egyik legfontosabb és legnehezebb része az egymással érintkező elemekből álló halmazok generálása. Számítógép segítségével szabályos halmazokat könnyű ugyan gyorsan generálni, de a szabályos geometriájú modellek megbízhatósága több szempontból is megkérdőjelezhető. Valós problémák elemzésénél ugyanis kulcs fontosságú, hogy a modellt alkotó szemcsék eloszlása jól egyezzen a vizsgált jelenségnél tapasztalttal. A véletlenszerű elrendezések generálásának hatékony lehetőségeit ezért számos matematikus és mérnök kutatja jelenleg is [Bagi; 2008].

A PFC szoftverekben első lépésként lehetőség van az ún. SSI eljárás alkalmazására egymással nem érintkező elemek véletlenszerű generálására, majd a felhasználó – választása szerint – valamely dinamikus halmaztömörítési eljárással csökkentheti a kívánt értékre a porozitást.

Érintkezések felismerése (ütközésdetektálás)

Az ütközésdetektálás azokat az időpontokat és helyzeteket azonosítja, amikor két részecske ütközik egymással. Az ütközésdetektálást három fő részre lehet felosztani: ütközésdetektálás, ütközés-meghatározás és válasz az ütközésre [Nagy; 2011]. Az ütközésdetektálás sebessége kulcsfontosságú a DEM rendszerek számára, mivel a részecskék számának növekedésével a számítási idő is négyzetesen növekszik. Ha a modelltér \(n\) mozgó és \(m\) statikus részecskét tartalmaz, akkor egy naiv megközelítés [Nagy; 2011]:

\[nm+\left(\frac{n}{2}\right)\]

ütközésvizsgálatot hajtana végre minden egyes lépés után. A kifejezés első tagja a statikus és dinamikus részecskék tesztjeinek a számát, míg a második tag dinamikus részecskék egymással történő tesztjeinek a számát adja meg. Az \(m\) és \(n\) növekedésével az elvégzendő tesztek száma nagy mértékben megnő! Ezért mindenképpen szükséges optimalizálni a kapcsolatok felismerését a DEM rendszerek számára. Ennek az első lépése az, hogy a nem lehetséges kapcsolatokat kizárjuk a vizsgálatból, hogy minél kevesebb részecske-párost kelljen vizsgálni. Ezután már csak a megmaradt kapcsolatokkal kell részletesebben foglalkoznunk. Ha találunk valódi kapcsolatokat a részecskék között, akkor alkalmaznunk kell a kapcsolati törvényt (lásd később) és a segítségével számolni a kapcsolati erőket és nyomatékokat. A lehetséges kapcsolatok meghatározására különféle algoritmusokat dolgoztak ki.

A modelltér felosztása

Az egyik leggyakrabban használt struktúra a szabályos térháló, implementálása és bejárása igen egyszerű. Amennyiben a modellterünk síkbeli, helyezzünk rá egy szabályos, (ax, ay) méretű, koordinátatengelyekkel párhuzamos oldalú téglalaprácsot. Térbeli esetben egy (ax, ay, az) kiterjedésű téglatesthálót használunk, ami a teljes modellteret befoglaló doboz felosztásával kapunk. Az előfeldolgozás alatt minden cellára meghatározzuk az adott cellába eső részecskéket (1. ábra). Ehhez azt kell megvizsgálni, hogy az adott részecskének van-­e közös része a cellával. Minden cellához informatikai oldalról egy láncolt lista tartozik. A láncolt listákban azoknak a részecskéknek az azonosítóját tároljuk el melyeknek van közös része a cellával. Ütközésdetektáláskor meg kell határoznunk, hogy a vizsgált részecske melyik cellába esik, ezután a tényleges ütközésdetektálást már csak a cellában lévő részecskékre kell elvégezni. Természetesen arra is van lehetőség, hogy az adott cella szomszédjainak hatásait is figyelembe vegyük, ez főleg molekuladinamikai számítások esetén fontos (Cell-list). Az eljárás hátránya, hogy nem veszi figyelembe a részecskék egymáshoz viszonyított helyzetét, azaz ugyanolyan aprólékosan oszt fel olyan térrészeket is, ahol nincs is részecske és másik végletként, megjelenhetnek „túlzsúfolt” cellák is, amelyekhez rengeteg részecske tartozik. Ezért elmondható, hogy bár lekérdezések szempontjából gyors és egyszerű ez az adatszerkezet, a gyakorlatban tárolás szempontjából már nem praktikus. Az optimális cellaméretet a DEM szoftverek heurisztikus módszerek segítségével becslik meg a kezdeti geometria felvétele után.

Szabályos térháló alapú felosztás (Cell-list)
1_cell-list

Egy másik lehetőség a részecskék nyilvántartására a Verlet-lista vagy táblázat. Vizsgáljunk meg egy \(N\) részecskéből álló rendszert. Képzeletben húzzunk minden részecske köré egy \(r_{c}\) és egy \(r_{l}=r_{c}+d\) sugarú kört, ahol \(d\) a részecske átmérőjét jelöli (2. ábra). Az \(r_{l}\) és \(d\) alapján megbecsülhetjük azon részecskék \(N_{max}\) maximális számát, amelyekkel egy adott részecske egy adott időpillanatban kölcsönhatásban állhat [Varga; 2002]:

\[N_{max}\sim\frac{r_{l}^{2}}{d^{2}}\]

A részecskék nyilvántartása a Verlet-lista segítségével
2_verlet-list

Létre kell hoznunk egy listát amiben az i indexű szemcse \(r_{l}\) környezetébe eső részecskék indexeit eltároljuk. Futásidőt azzal nyerünk, hogy a szomszédsági viszonyokat tároló listát nem kell minden iterációs lépés után frissíteni, mert az egymást követő időlépésekben nem biztos, hogy a szemcse kiment az \(r_{l}\) sugarú körből vagy oda új szemcse érkezett. A frissítések sűrűsége függ a részecskék között előforduló legnagyobb \(v_{max}\) sebességtől és a \(\Delta t\) iterációs lépés hosszától. Frissítést elegendő minden \(n\)-edik iterációs lépésben elvégezni [Kun; 2011]:

\[n=\frac{d}{2\Delta{t}\dot{}v_{max}}\]

ahol \(n\) számértéke akár 10 – 20 is lehet, ami szignifikánsan redukálja a program futásidejét. Még ilyen egyszerűsítések mellet is a DEM programok a számítási idő 90 %-át a részecskék közötti kapcsolatok detektálására fordítják [O’Sullivan; 2011].

Befoglalókeretek alkalmazása

A modelltér felosztása lehetővé teszi, hogy egy adott szemcse környezetét gyorsan lekérdezzük. Ez azonban még nem jelenti azt, hogy minden szomszédjával kapcsolatba is lép az adott szemcse. Az igen költséges ütközésszámítást tovább gyorsíthatjuk, ha bevezetjük a befoglaló keretek fogalmát. A befoglaló keret egy egyszerű geometriájú objektum, tipikusan gömb (kör) vagy téglatest (téglalap), amely egy-egy bonyolultabb objektumot teljes egészében tartalmaz. Az ütközésdetektálás alatt először a befoglaló kereteket próbáljuk elmetszeni egymással. Ha nincs metszéspont, akkor nyilván a befoglalt objektummal (szemcsével) sem lehet metszéspont, így a bonyolultabb számítást megtakaríthatjuk [Szirmay-Kalos, Antal és Csonka; 2003]. Amennyiben a befoglaló alakzatok egymásba hatolnak, akkor folytatni kell a vizsgálatot. Az egyik objektumot összevetjük a másik befoglaló alakzatával, és ha itt is ütközés mutatkozik, akkor magával az objektummal (3. ábra). A koncepciót kör alakú részecskék és koordinátatengelyekkel párhuzamos élű befoglaló dobozok (Axis-Aligned Bounding Box) esetén mutatjuk be (1. algoritmus).

algoritmus1

A befoglaló keretek (AABB) koncepciója
3_boundingbox

Az AABB előállítása egyszerű, a csúcspontok maximális és minimális x, y, z koordinátái alkotják a téglatest két szemközti csúcsának koordinátáit. Két AABB akkor hatol egymásba, ha valamennyi alábbi egyenlőtlenség fennáll: \(x_{min}^{1}\le x_{max}^{2}\), \(x_{max}^{1}\ge x_{min}^{2}\), \(y_{min}^{1}\le y_{max}^{2}\), \(y_{max}^{1}\ge y_{min}^{2}\), \(z_{min}^{1}\le z_{max}^{2}\), \(z_{max}^{1}\ge z_{min}^{2}\).

Tényleges ütközésszámítás

A részecskék befoglaló dobozainak közös metszete szükséges, de nem elégséges feltétele az ütközésnek. A pontos ütközésszámítást kör alakú szemcsék esetén mutatjuk be (2. algoritmus). A számítás a Pitagorasz-tételen alapul. Ütközésről akkor beszélünk, ha a két szemcse sugarának összege kisebb vagy egyenlő, mint a két szemcse között értelmezett távolság (4. ábra). Poligonok ütközésszámításáról részletesen olvashatunk Cozic [2006] internetes munkájában.

algoritmus2

Kör alakú szemcsék pontos ütközésszámítása
4_collison

A kapcsolatok leírása

Merev elemek esetén a kapcsolatok anyagjellemzőinek megfelelő megválasztásával kell biztosítani, hogy a halmaz egésze valósághűen viselkedjen. A legegyszerűbb modellekben a kapcsolatok pontszerűen kicsinyek, tehát koncentrált erőt (és esetleg koncentrált nyomatékot) adhatnak át egymásnak az érintkező elemek. A kapcsolat mechanikai anyagmodellje egyrészt azt írja le, hogy ezek az erők ill. nyomatékok hogyan függenek a két elem relatív elmozdulásától, másrészt teherbírási korlátokat fogalmazhat meg, illetve azt, hogy képlékeny állapotba kerülése után a további terhelés esetén hogyan viselkedik a kapcsolat [Bagi; 2008]. A kapcsolati törvényt akkor kell alkalmazni, amikor két szemcse egymással érintkezik (ütközik). A kapcsolati törvényt (viselkedést) a reológia segítségével fogalmazhatjuk meg. A reológia három alapvető modell segítségével írja le az anyagok tulajdonságait:

  • elasztikus (Hooke-test)
  • viszkózus (Newton-test)
  • plasztikus (Saint-Venant-test)

Elasztikus az a test vagy kapcsolat, amelynek relatív alakváltozása arányos a testre ható mechanikai feszültség értékével („Amilyen a nyúlás, olyan az erő”). Jellemző sajátossága a rugalmassági modulusz vagy Young-modulus. Ez anyagi tulajdonság, amely a feszültség és a relatív alakváltozás értékét kapcsolja össze (5.a ábra). A Hooke-test esetében reológiai hatás nincsen, a testre ható terhelés hatására keletkező alakváltozás az időtől független [Dulácska, Fekete és Varga; 1982]. Jelképe a helyettesítő vázlatokon: a rugó.

Reológiai alapmodellek
5_rheology

Viszkózus az a test vagy kapcsolat, amelyen állandó \(\tau\) nyírófeszültség (csúsztatófeszültség) állandó \(\frac{\mbox{d}\epsilon}{\mbox{d}t}\) sebességű folyási jelenséget hoz létre. Newton törvénye értelmében: \(\tau=\eta\frac{\mbox{d}\epsilon}{\mbox{d}t}\), ahol \(\eta\) a dinamikai viszkozitási együttható (5.b ábra). A Newton-testnél a pillanatnyi, igen rövid ideig tartó terhelés hatására nincsen alakváltozás. Tartós terhelés esetén lép csak fel alakváltozás, mely a terhelési időintervallum és a terhelés nagyságától függ [Dulácska, Fekete és Varga; 1982]. Jelképe a helyettesítő vázlatokon: a dugattyúhoz hasonló ábra.

Ha a plasztikus testre vagy kapcsolatra a \(\tau_{0}\) határfeszültségnél kisebb hat, alakváltozás nem jön létre (5.c ábra). A határfeszültséget elérve az alakváltozás minden határon túl növekszik. Fontos megjegyezni, hogy a jelenség nem írható le az idővel összefüggésben, ezért a sebesség fogalma itt nem értelmezhető. Az alakváltozás létrejötte csupán egy feltételtől függ: a feszültség az ok, és az alakváltozás az okozat. A hétköznapi életben ilyen jelenség a súrlódás. Szintén plasztikus jelenség a törés, a Mohr-féle törési elmélet szerint. Jelképe a helyettesítő vázlatokon: két egymásban súrlódásosan mozgó test.

Lineárisan rugalmas kapcsolat

A legegyszerűbb mechanikai modell, amit lényegében minden DEM programban megtalálunk egy lehetséges opcióként, a lineárisan rugalmas kapcsolat [Cundall és Strack; 1979]. A részecskék normálmerevségét \(K_{n}\), míg nyírómerevségét \(K_{s}\) rugóállandók írják le. Itt a normálerő csak nyomás lehet (konvenció szerint a húzóerők pozitívak), és nagysága az érintkezést alkotó két anyagi pont relatív eltolódási komponensével (\(u_{n}\)) arányos (6. ábra):

\[F_{n}=K_{n}\dot{}u_{n}\]
ahol \(F_{n} \le{}0\).

A lineárisan rugalmas kapcsolat
6_particle-model1

Ha a normálerő pozitívvá válna (azaz a két elem eltávolodik egymástól), a kapcsolat megszűnik [Bagi; 2008]. A kapcsolati nyíróerő az érintőirányú relatív eltolódással (\(u_{s}\)) arányos:

\[F_{s}=K_{s}\dot{}u_{s}\]

Abban az esetben, ha rúgóállandókat nem a kapcsolathoz hanem a részecskékhez rendeljük, mint anyagjellemző, akkor az egyenértékű normál- és nyírómerevségeket a következő összefüggésel becsülhetjük [Ardic; 2006]:

\[K_{n}=\frac{k_{n1}\dot{}k_{n2}}{k_{n1}+k_{n2}}\]

és

\[K_{s}=\frac{k_{s1}\cdot k_{s2}}{k_{s1}+k_{s2}}\]

ahol \(k_{ni}\) és \(k_{si}\) az \(i\)-edik részecskéhez tartozó saját normál- és nyírómerevség paraméter. A kapcsolatok merevségét a  \(\lambda=\frac{K_{s}}{K_{n}}\) arány jellemzi.

A lineárisan rugalmas Coulomb-súrlódásos kapcsolat

A szilárd testek között fellépő súrlódást száraz súrlódásnak vagy Coulomb-súrlódásnak nevezzük. A Coulomb-féle súrlódás az érintkezési felületen fellépő jelenségeknek csupán egy lehetséges (nagyon egyszerű) modellje, számos más modell is létezik. Hatását egy Hooke-test és egy Saint-Venant-test sorba kötésével vehetjük figyelembe. A súrlódásról fontos megjegyezni, hogy az mindössze reakcióerő, ami azt jelenti, hogy csak akkor lép fel, ha egy aktív erő a testet el akarja mozdítani vagy már elmozdította, és ilyenkor mindig a pillanatnyi elmozdulással ellentétes irányban hat. Ha egy nyugalomban lévő testre nem hat a támaszkodó felületével párhuzamos irányú erő, akkor ott nem is ébred súrlódás. A nyíróerő nagyságának a Coulomb-feltétel szab határt (7. ábra):

\[F_{s}^{max}=F_{n}\dot{}\tan\left(\varphi\right)=F_{n}\dot{}\mu\]

ahol \(\varphi\) a rendszerrel jellemző súrlódási szög, \(\mu\) pedig a dimenzió nélküli súrlódási tényező (beszélhetünk statikus és dinamikus súrlódási tényezőről).

A lineárisan rugalmas Coulomb-súrlódásos kapcsolat
7_particle-model2

A súrlódási tényező az érintkező felületek anyagminőségétől függő empirikus mennyiség. Ha a nyíróerő eléri ezt az értéket, a további relatíveltolódás-növekmény már állandó erő mellett (zérus nyírási merevséggel) történik, mindaddig, amíg a relatív eltolódás iránya meg nem változik [Bagi; 2008]. A szimuláció alatt minden időlépésre számoljuk az érintkező részecskék között fellépő normál (\(u_{n}\)) és érintőirányú (\(u_{s}\)) relatív elmozdulásokat, majd ezek segítségével a normál- és nyíróerőket. A nyíróerők új értékét \(F_{s}^{new}\) az előző időlépés alatt számolt nyíróerők \(F_{s}^{old}\) és az újonnan számolt \(\Delta F_{s}\) érték összegéből képezzük [Jakob és Konietzky; 2012]:
\[F_{s}^{new}=F_{s}^{old}+\Delta F_{s}\]
ahol
\[\Delta F_{s}=-K_{s}\dot{}\Delta u_{s}\]

A teljes kapcsolati erő, ami a részecskéket éri:

\[F_{c}=F_{n}+F_{s}\]

A bemutatott kapcsolati törvény jól alkalmazható triaxiális vizsgálatok szimulációjára és számos kutatási eredmény igazolta (közelítő volta ellenére) helyességét [Scholtes; 2009].

A Hertz-Mindlin féle kapcsolat

A részecskék merevségét leíró (\(k_{n}\)és (\(k_{s}\)) rugóállandókat a Hertz-Mindlin-féle közelítés segítségével becsülhetjük meg a szemcsék nyírási modulusa (\(G_{1}\), \(G_{2}\)), Poisson-tényezője (\(\nu_{1}\), \(\nu_{2}\)), sugaraik (\(R_{1}\), \(R_{1}\)) és normál irányú relatív elmozdulásuk (\(u_{n}\)) segítségével [Arévalo-Mendoza et al., 2014]:

\[k_{n}=\left(\frac{2\bar{G}\sqrt{2\bar{R}}}{3\left(1-\bar{\nu}\right)}\right)\sqrt{u_{n}}\\]
ahol
\[\bar{R}=\frac{2R_{1}R_{2}}{R_{1}+R_{2}},\]
\[\quad\bar{G}=\frac{G_{1}+G_{2}}{2},\]
\[\quad\bar{\nu}=\frac{\nu_{1}+\nu_{2}}{2}\]

A kapcsolat nyírómerevségét (amely a nyomóerő nagyságától függ) Mindlin és Deresiewicz [1953] és Deresiewicz [1958] munkája alapján becsülhetjük [Van Baars; 1996]:

\[k_{s}=\left(\frac{2\left(3\bar{G}^{2}\left(1-\bar{\nu}\right)\bar{R}\right)^{\frac{1}{3}}}{2-\bar{\nu}}\right)\cdot\left|F_{n}\right|^{\frac{1}{3}}\]

Fel kell hívnunk arra a figyelmet, hogy az idézett cikkekben a levezetések csak igen speciális esetekre vonatkoznak (azonos méretű, tökéletesen rugalmas-képlékeny gömbök kapcsolata), így hiába bonyolultak, mégis csak igen durva közelítésnek tekinthetők [Bagi; 2008].

A kapcsolati modell viszkózus csillapítása

Amikor két érintkező felület elmozdul egymáson a deformációhoz szükséges munka egy része nem a rugalmas alakváltozást fedezi, hanem hővé alakul. A valóságos folyamatok helyes leírása miatt csökkentenünk kell a rendszer mozgási energiáját, azaz foglalkoznunk kell az energiadisszipácó kérdésével (visszafordíthatatlanul hőenergiává alakuló mozgási energia).

A kapcsolati modellt ki kell egészíteni a viszkózus csillapítással, melyet az jellemez, hogy a csillapító erő arányos a tömeg sebességével (8. ábra). Ezt a csillapítást azért nevezik viszkózusnak, mert a folyadék csillapítását modellezi. A \(c\) arányossági tényezőt (vagyis \(\eta\) dinamikai viszkozitási együtthatót) csillapítási tényezőnek nevezik, az összefüggés általános esetre:

\[D_{c}=-cv=-c\dot{x}=-c\frac{dx}{dt}\]

ahol \(D_{c}\) a csillapítási erő (damping force).

A kapcsolati modell viszkózus csillapítása
8_particle-model3

A kapcsolati viszkózus csillapítás lényege az, hogy a kapcsolati erő mindegyik komponenséhez hozzáadódik egy-egy olyan erő (ill. nyomaték), amelyeknek nagysága a kapcsolati relatív elmozdulás-sebesség megfelelő komponenseinek nagyságával arányos, iránya pedig olyan, hogy épp lassítsa a relatív sebességet (Bagi, 2008). A normálerőhöz hozzáadódó komponens nagysága például:

\[F_{n}=k_{n}u_{n}-c_{n}\dot{u}_{n}\]

A nyíró erő számításánál is hozzá kell adnunk a viszkózus csillapítási tényező hatását (\(c_{s}\cdot\Delta\dot{u}_{s}\)):

\[\Delta F_{s}=-k_{s}\cdot\Delta u_{s}-c_{s}\cdot\Delta\dot{u}_{s}\]

ahol \(c_{n}\) és \(c_{s}\) a normál és nyírási csillapítás tényezők. Ha a csillapítás eléggé kicsi, a rendszer rezegni fog, de idővel megszűnik a rezgése. Ez az eset a gyakorlat szempontjából a legfontosabb. Ha a csillapítást addig növeljük, míg a rendszer éppen megszűnik rezegni, akkor elértük a kritikus csillapítást. A kritikus csillapításhoz tartozó csillapítási tényező értéke egytömegű lengőrendszer esetén:

\[c^{crit}=2\sqrt{k\cdot{} m}\]

ahol \(k\) arányossági együttható a rugómerevség, egysége erő/elmozdulás (N/m).

A rendszer csillapításának mértékéül a csillapítási viszonyt (\(\beta\)) vezették be (ezt csillapítási faktornak is hívják). Ez a csillapítási viszony a tényleges csillapítási tényező és a kritikus csillapítási tényező hányadosa. Képlete az egytömegű lengőrendszer esetén:

\[\beta=\frac{c}{2\sqrt{k\cdot{} m}}\]

Az összefüggések segítségével számíthatjuk a \(c_{n}\) és \(c_{s}\) értékét:

\[c_{n}=\beta_{n}c_{n}^{crit}=2\beta_{n}\sqrt{k_{n}\cdot{}\bar{m}}\]
\[c_{s}=\beta_{s}c_{s}^{crit}=2\beta_{s}\sqrt{k_{s}\cdot{}\bar{m}}\]

ahol \(\bar{m}\) a rendszer hatékony tömege:

\[\bar{m}=\dfrac{m_{1}\cdot{} m_{2}}{m_{1}+m_{2}}\]

A \(\beta\) csillapítási viszonyt az \(e\) ütközési szám (coefficient of restitution) alapján becsülhetjük [Schafer et al.; 1996]:

\[\beta=\frac{\ln\left(e\right)}{\sqrt{\ln^{2}\left(e\right)+\pi^{2}}}\\]
ahol
\[e=\frac{v_{a}}{v_{b}}\]

és \(v_{a}\) az ütközés utáni, \(v_{b}\) pedig az ütközés előtti relatív sebessége a részecskéknek.

A mozgásegyenlet

A klasszikus mechanikában a merev test a véges nagyságú szilárd test idealizált modellje, amelynél az alakváltozást elhanyagolják. Más szóval a merev test bármely két pontjának távolsága időben állandó, függetlenül az esetleg rá ható erőhatásoktól. A merev test idealizációja szigorúan csak a klasszikus mechanika szerint használható. A diszkrét elemek helyzetének és mozgásának leírására háromdimenziós, derékszögű Descartes-koordinátarendszert fogunk alkalmazni, amelyet az (x, y, z) tengelyek alkotnak. A merev testnek hat szabadságfoka van: derékszögű koordináta-rendszerben az x, y és z tengely irányába történő elmozdulás (transzláció) és az x, y és z tengely körüli elfordulás (rotáció).

Minden részecskének kijelölünk egy referenciapontot, amely elvileg az elem bármely pontja lehet, de a számításokat egyszerűsíti, ha egybeesik az elem tömegközéppontjával. Merev test modell esetén a részecskék deformálhatóságát a pontszerűen kialakuló kapcsolatokba sűrítjük.

Az \(\vec{x}_{1}\) és \(\vec{x}_{2}\) vektorok a részecskék referenciapontjára mutatnak, az \(\vec{n}\) és \(\vec{t}\) pedig a normál és érintő irányú egységvektorokat jelölik (9. ábra). A következő összefüggések adódnak:

\[d_{a}=\left|\vec{x}_{2}-\vec{x}_{1}\right|\:\]
\[n=\dfrac{\vec{x}_{2}-\vec{x}_{1}}{d_{a}}\:\]
\[u_{n}=R_{1}+R_{2}-d_{a}\]

Két részecske érintkezése
9_particle-model-main

A két szemcse között kialakuló kapcsolati pont \(\vec{x}_{c}\) vektorát a következő képen számíthatjuk:

\[\vec{x}_{c}=\vec{x}_{1}+(R_{1}-1/2\cdot{}u_{n})\cdot\vec{n}\]

A részecskék között kialakuló kapcsolati erő \(F_{c}\) normál és nyíró irányú összetevőjét a \(K_{n}\) és \(K_{s}\) merevségi értékek határozzák meg. A normál komponens a részecskék között kialakuló kapcsolat alatt közvetlenül számítható. A nyírási komponens ilyenkor még nulla, az egyes időlépések alatt fog növekedni vagy csökkenni az értéke:

\[\Delta F_{s}=-K_{s}\cdot\Delta u_{s}\]
és
\[F_{n}=K_{n}\cdot{}u_{n}\]

A nyírási elmozdulások összege egy időlépés alatt \(\Delta u_{s}\) amit a fejezet végén az (47) összefüggéssel számoljuk. A két szemcse között kialakuló kapcsolati erőt a következő képlet adja meg:

\[\vec{F}_{c}=F_{n}\cdot{}\vec{n}+F_{s}\cdot{}\vec{t}\]

Az egy részecskére ható teljes erő \(\vec{F}\left(\hat{=}F\right)\), a szomszédos részecskéről átadódó kapcsolati erőkből és a nehézségi erőből tevődik össze \(F_{g}=m\cdot\vec{g}\):

\[F =\sum_{c}\vec{F}_{c}+\vec{F}_{g}\]

A teljes erő az i-edik irányban \((i\in\left\{ 1,\,2,\,3\right\})\) \(F_{i}\), egy részecskére kifejtett hatását Newton II. törvénye határozza meg. Az \(F_{i}\) értékét a szemcse tömegének (\(m\)) és az öt ért gyorsulások \(\ddot{u}_{i}\) és \(g_{i}\) összegének szorzataként határozzuk meg:

\[F_{i}=m\cdot{}\left(\ddot{u}_{i}+g_{i}\right)\]

A fenti képlet átrendezésével kapjuk:

\[\ddot{u}_{i}=\frac{F_{i}}{m}-g_{i}\]

Az idő (\(t\)) szerinti első integráltja a gyorsulásnak a sebességet, második integráltja pedig az elmozdulást adja meg:

\[\dot{u}_{i}=\int\ddot{u}_{i}dt\\]
és
\[\qquad u_{i}=\int\dot{u}_{i}dt\]

A forgó mozgás leírására vezessük be az \(\vec{r}_{c}\) vektor, mely a szemcse referenciapontjából a ütközési (kapcsolati) pontba mutat. Az \(\vec{M}\) forgató nyomatékot a következő összefüggés adja:

\[\vec{M}=\sum_{c}\left(\vec{r}_{c}\times\vec{F}_{c}\right)\]

ahol \(\vec{r}_{c}=\vec{x}_{c}-\vec{x}\) és \(\vec{x}=\vec{x}_{1}\) az 1-es szemcse esetén.

Az \(M_{i}\) forgató nyomaték az i-edik irányban a tehetetlenségi nyomaték (\(\Theta\)) és a szöggyorsulás (\(\beta=\dot{\omega}=\ddot{\varphi}\)) szorzataként számítható:

\[M_{i}=\Theta\cdot\beta_{i}=\Theta\cdot\dot{\omega}_{i}\]

A tehetetlenségi nyomaték, a tömeggel analóg mennyiség forgómozgásnál. Vagyis a tehetetlenségi nyomaték a forgást végző merev test forgási tehetetlensége. Gömb alakú részecskék tehetetlenségi nyomatéka:

\[\Theta=\frac{2}{5}\cdot{}mR^{2}\]

Körlapra merőleges tengelyre (korong):

\[\Theta=\frac{1}{2}\cdot{}mR^{2}\]

Az (40) összefüggés átrendezésével és idő szerinti integrálásával kapjuk a szögsebességet (\(\omega\)):

\[\omega_{i}=\int\frac{5M_{i}}{2mR^{2}}dt\]

A test érintőirányú (tangenciális) sebességét \(v_{t}\) (kerületi sebességét) a következőképpen számíthatjuk általános esetben:

\[v_{t}=\frac{ds}{dt}=r\cdot{}\frac{d\varphi}{dt}=r\cdot{}\omega\]

ahol az \(r\) a kör sugarát jelöli és \(s=r\cdot{}\varphi\) a körmozgást végző test útfüggvénye. A két összetalálkozott szemcse sebessége az ütközés után:

\[\vec{v}_{1}^{‘}=\vec{v}_{1}+\vec{\omega}_{1}\times\left(\vec{x}_{c}-\vec{x}_{1}\right)\]

\[\vec{v}_{2}^{‘}=\vec{v}_{2}+\vec{\omega}_{2}\times\left(\vec{x}_{c}-\vec{x}_{2}\right)\]

ahol \(\vec{v}_{1}\hat{=}\dot{u}_{i1}\) és\(\vec{v}_{2}\hat{=}\dot{u}_{i2}\) a részecskék haladó sebessége, \(\vec{\omega}_{1}\) és \(\vec{\omega}_{2}\) pedig szögsebessége. A részecskék ütközés utáni sebességkülönbsége (\(\vec{v}_{rel}\)):

\[\vec{v}_{rel}=\vec{v}_{2}^{‘}-\vec{v}_{1}^{‘}=\left(\vec{v}_{2}+\vec{\omega}_{2}\times\left(\vec{x}_{c}-\vec{x}_{2}\right)\right)-\left(\vec{v}_{1}+\vec{\omega}_{1}\times\left(\vec{x}_{c}-\vec{x}_{1}\right)\right)\]

Az ütközés utáni normálirányú sebességkülönbség pedig (\(v_{n}\)):

\[v_{n}=\left|\vec{v}_{rel}\cdot{}\vec{n}\right|\]

Az érintő irányú elmozdulásnövekményt (\(\Delta u_{s}\)) minden egyes időlépés alatt a nyírósebesség (\(v_{s}\)) és a \(\Delta t\) időlépés szorzataként számítjuk:

\[\Delta u_{s}=v_{s}\cdot\Delta t\]

A nyírósebességet pedig a \(v_{rel}=\left|\vec{v}_{rel}\right|\) és a normálirányú sebességkülönbség különbségéből kapjuk:

\[v_{s}=v_{rel}-v_{n}\]

A mozgásegyenlet matematikai formáját tekintve egy közönséges differenciálegyenlet. A numerikus megoldási eljárások az úgynevezett véges differencia módszerre épülnek. Ezek alapgondolata az, hogy az időt, mint független változót egyenközű \(\Delta t\) lépésekkel diszkretizáljuk és a megoldásfüggvényt a csonkolt Taylor sorával közelítjük [Kun; 2011]. Ezzel szemben a DEM szoftverekben a numerikus megoldás stabilitásának biztosítása érdekében nem célszerű egyenközű \(\Delta t\) eljárást alkalmazni (lásd 3.4 pont). A PFC-ben is változó időlépésekkel történik az időintegrálás. A DEM modellekben alkalmazott eljárásokról részletesen olvashatunk Bagi [2008] és Kun [2011] munkáiban. A bemutatott eljárás dinamikai modelljére a PFC program ismertetésénél térünk ki részletesebben.

Diszkrét elemes modellezés PFC2D/3D-szoftverrel

Az elméleti alapok ismertetése után, most röviden kitérnénk az egyik legnépszerűbb diszkrételemes modellezést támogató szoftvercsomag bemutatására. Az Itasca Consulting Group Inc. által kifejlesztett PFC (Particle Flow Code) programok a gömbszerű szemcsék mozgásait és kölcsönhatásait modellezik két vagy három dimenzióban Cundall és Strack 1979-ben publikált diszkrét elemes módszere [Cundall, Strack; 1979] alapján. A PFC programok fő jellemzőit Mészöly [1999], Tóth [2004], Fischer és Horvát [2010] munkái alapján foglaltuk össze a teljes igénye nélkül. Akik további részletekre kíváncsiak, azoknak javasoljuk a felsorolt munkák részletes tanulmányozását. A PFC program a szemcsehalmazok modellezésénél a következő feltételezésekkel él [Mészöly; 1999]:

  • A részecskék merev testek.
  • Minden szemcse kör (2D) vagy gömb (3D) alakú, de a szemcsék egymáshoz rögzítésével tetszőleges alakú merev testek hozhatók létre.
  • A kapcsolatok egy pontban jönnek létre.
  • A kapcsolati erők csak az érintkezési pontokon adódnak át.
  • Puha kapcsolatot feltételez, a részecskék átfedhetik egymást.
  • Két szemcse között ébredő kapcsolati erő ezen átfedés mértékével arányos, amit az úgynevezett „büntető merevsége” jellemez. Ily módon a szemcsék alakja ugyan nem változik, de mégis megjelenik a kapcsolati deformációknak megfelelő viselkedés.
  • A kapcsolati deformáció a szemcse méretéhez képest mindig kicsi.
  • Kötés modell esetén húzásnak, nyomásnak, nyírásnak és hajlításnak is ellenálló kötések létesülhetnek a kapcsolati pontokon a szemcsék között.

A PFC2D alkalmazás 2D-os modellezésre alkalmas: a mozgástörvényekben csak két-két (síkban fekvő) elmozdulás és erő, illetve egy-egy (síkra merőleges) elfordulás és nyomaték komponens szerepel. A kör alakú szemcséket a program vagy adott sugarú gömbökként, vagy – a síkra merőleges irányba kiterjedő – adott szélességű hengerekként kezeli, és ennek megfelelően számítja tehetetlenségi jellemzőiket [Tóth; 2004].

A PFC programok két elemet ismernek: a szemcsét (partice) és a falat (wall). Mind a szemcsék, mind a falak végtelen merev elemek, a köztük (szemcse-fal) és egymás közt (szemcse-szemcse) fellépő érintkezések, kapcsolatok tulajdonságai határozzák meg az anyag viselkedését. Kapcsolat fal-fal között nem jöhet létre. A szemcsék hivatottak szimbolizálni a szemcsés anyagban (vagy egyéb diszkrét elemekkel modellezhető anyagokban) a „szemeket”, a falak pedig határoló funkciót látnak el alapesetben, de a perem- és kezdeti feltételek megadásánál is fontos szerepet játszanak [Fischer, Horvát; 2010].

A fal mozgása megadható és segítségével korlátlan erőt lehet kifejteni, mivel a program nem ír fel rá Newton-törvényt. Ahogy a falaknak úgy a szemcséknek is lehet előírt sebessége, mely a számítás alatt nem változik. Két szemcsét a felhasználó által megadott kapcsolati erő köthet össze. Minden elemhez megadható a normál- és nyírási merevség.

A mozgásegyenletek

A rendszer dinamikai viselkedése egy időléptető algoritmuson alapul, ahol az egyes időlépéseken belül a sebességek és gyorsulások állandóak. A feltételezés szerint az időlépés olyan rövidre választott, hogy a hatás csak a szomszédos részecskékre jut el, így minden pillanatban a következő állapot kiszámításához csak a szomszédos részecskékkel való kölcsönhatást kell figyelembe venni [Mészöly; 1999].

A szemcserendszer deformációja a szemcsék elcsúszásából és elfordulásából származik, nem pedig az egyes szemcsék deformációjából. Az energiaelnyelés pedig csúszási súrlódás által történik. Elképzelhető az is, hogy a csúszási modell nem működik az adott modellben, így szükség van az energia elnyeléséhez helyi csillapítások alkalmazására. A PFC programban alkalmazott csillapítás konstans, csak gyorsuló mozgást csillapít és frekvenciától függetlenül csillapítja a sajátrezgéseket [Mészöly; 1999].

A számítás Newton II. törvényének és a kapcsolati erő-elmozdulás törvénynek az alkalmazása között alternál (10. ábra): Newton II. törvénye alapján kiszámíthatók a különálló szemcsék elmozdulásai az épp aktuális kapcsolati és térfogati erők hatására, míg a kapcsolati erő-elmozdulás törvény segítségével módosíthatók a relatív elmozdulások miatt megváltozott kapcsolati erők [Tóth; 2004].

Számítási ciklus a PFC programban
10_dem-flow

A PFC alkalmazás a sebességeket (\(\dot{x_{i}}\) és \(\omega_{i}\)) átlagolt \(t\pm\Delta t/2\) időközönként, az \(x_{i}\), \(\ddot{x}_{i}\),\(\dot{\omega}_{i}\), \(F_{i}\) és \(M_{i}\) értékeket pedig \(t+\Delta t\) intervallumra számítja. A (37), (38) és (43) egyenletek numerikus megoldása a centrális differenciák módszerével az alábbi formában írható fel:

\[\ddot{u}_{i}^{(t)} =\frac{1}{\Delta t}\left(\dot{u}_{i}^{(t+\Delta t/2)}-\dot{u}_{i}^{(t-\Delta t/2)}\right)\]

\[\dot{\omega}_{i}^{(t)} =\frac{1}{\Delta t}\left(\omega_{i}^{(t+\Delta t/2)}-\omega_{i}^{(t-\Delta t/2)}\right)\]

A mozgástörvény egy ciklusa [Mészöly; 1999]:

1. Adott a szemcsék tömegközéppontjainak sebessége \(\dot{u}_{i}^{(t-\Delta t/2)}\) és szögsebessége \(\omega_{i}^{(t-\Delta t/2)}\) a \(t-\Delta t/2\) pillanatban, a szemcsék tömegközéppontjának helyzete \(x_{i}^{(t)}\) és a rá ható eredő erő \(F_{i}^{(t)}\) és nyomaték \(M_{i}^{(t)}\) a \(t\) időpillanatban.

2. Az (49), (50), (36) és (40) egyenletek segítségével a szemcsék tömegközéppontjának sebessége és szögsebessége számítható a \(t+\Delta t/2\) időpillanatra:

\[\dot{u}_{i}^{(t+\Delta t/2)} =\dot{u}_{i}^{(t-\Delta t/2)}+\left(\frac{F_{i}^{(t)}}{m}+g_{i}\right)\Delta t\]
\[\omega_{i}^{(t+\Delta t/2)} =\omega_{i}^{(t-\Delta t/2)}+\left(\frac{5M_{i}^{(t)}}{2mR^{2}}\right)\Delta t\]

3. Ismeretükben a szemcse \(t+\Delta t\)-beli helyzete meghatározható:

\[x_{i}^{\left(t+\Delta t\right)}=x_{i}^{(t)}+\dot{u}_{i}^{\left(t+\Delta t/2\right)}\Delta t\]

4. A számítás újrakezdődik az (27), (28), (29),összefüggésektől, az erőket és elmozdulásokat az új pozíció segítségével becsüljük a következő időlépésre.

A számítás további részleteit az érdeklődő olvasó megtalálja a PFC program kézikönyvében [Itasca; 2008].

Kapcsolati modellek

A PFC-ben különféle kapcsolati modelleket és kötéseket lehet alkalmazni a szemcsék között, amik döntően befolyásolják a halmaz egészének viselkedését [Tóth; 2004].

A kapcsolati modell (Contact Model)

Ahogy azt már említettük, mind a szemcsék, mind a falelemek végtelen merevek, ezért a köztük fellépő érintkezések tulajdonságai határozzák meg a modell viselkedését [Fischer, Horvát; 2010]. A PFC program két kapcsolati modellt alkalmazhat: lineáris rugómodell vagy egyszerűsített Hertz-Mindlin modellt. Két olyan szemcse között, melyekre nem azonos modellt alkalmazunk, nem jöhet létre kapcsolat, mert viselkedésük bizonytalan lenne. A Hertz-modell összeférhetetlen minden kötési modellel, mivel ott nincsenek értelmezve a vonzóerők [Mészöly; 1999].

A súrlódásos modell (Slip Model)

Elcsúszás akkor következik be egy kapcsolatnál, ha nincs érintkezési kötés és a kapcsolati erő nyíró komponensének nagysága a maximálisan megengedhető nyíróerőt 0,1%-nál jobban megközelíti:

\[F_{s}>0,99\cdot{}F_{s}^{max}\]

Az érintkezési kötés és a csúszási modell közül mindig csak az egyik lehet aktív. A kapcsolati modellel viszont működhet párhuzamosan. Egy kapcsolatban az egymással érintkező két szemcse súrlódási tényezője közül a kisebbet veszi figyelembe.

A kötési modell (Bond Model)

Annak érdekében, hogy ne csupán belső kötés nélküli szemcsés anyagokat lehessen modellezni a PFC programban, olyan lehetőség is rendelkezésére áll a felhasználónak, hogy a szimulációban normálirányú igénybevételt, de akár nyomatékbíró kapcsolatokat is létre tud hozni [Fischer, Horvát; 2010].

Részecskék összekapcsolására két kapcsolati kötési modellt használ a program. Az egyik az érintkezési-, a másik a párhuzamos kötési modell. Az előbbi csak erőátadásra képes és végtelen kis területen hat (egy pontban), utóbbi véges méretű területen erők és nyomatékok átadására is alkalmas. Mindkét esetben a kapcsolat felszakad, ha a hozzájuk tartozó erőt túllépjük, vagy a kötési erők értékét nullára állítjuk [Mészöly; 1999].

Az érintkezési kötési modell (Contact-Bond Model) egy rugópárként fogható fel. A csúszási modell ilyenkor nem működik! Két paraméterrel, az érintkezéses kapcsolat normál és nyíró komponensével lehet megadni. A nyíróerőt a nyíró kötési erő, a húzóerőket pedig a normál kötési erő korlátozza. Megszűnik a kötés, ha a húzóerő eléri vagy túllépi a kapcsolat normál erejét, vagy ha a nyíróerő túllépi a kapcsolat nyíró kötési erejét. Ez a kapcsolati erőket nem változtatja meg. A kapcsolat végtelen kis területen, egy pontban jön létre, így nyomaték átadására nem alkalmas.

A párhuzamosan kötött kapcsolatok (Parallel-Bond Modell) opciójával vagyunk képesek modellezni az aszfalt vagy cementes kötésű betonszerkezetek (11. ábra). Ezek nem csak erő, de nyomaték felvételére is képes kapcsolatot biztosítanak a szemcsék között. Úgy kell ezt a típusú kapcsolatot elképzelni, mintha rugalmas ragasztóval ragasztották volna össze a szemcséket. Párhuzamosan kötött kapcsolatoknál is ugyanazok az elvek érvényesek az érintkezések kezelésénél, mint az érintkezéses kapcsolatoknál [Fischer, Horvát; 2010].

A párhuzamos kötési modell
11_parallelbound

Csillapítások

A különféle csillapítások célja az, hogy a rendszer mozgási energiáját csökkentsék, részben az egyensúlyhoz való gyors konvergálás érdekében, részben pedig a valóságban is lejátszódó energiadisszipáció modellezése céljából [Bagi; 2008]. Alapvetően kétféle csillapítási eljárást alkalmaznak a BALL-típusú modellek (és így a PFC is), az egyik a lokális a másik pedig a kapcsolati viszkózus csillapítás. Az utóbbi lényegét részletesen bemutattuk a 2.3.4. pontban, ezért most csak a lokális csillapításra térünk ki.

A lokális csillapítás lényege, hogy minden elem mozgásegyenletében az elemekre ható eredőt módosítjuk úgy, hogy hozzáadunk egy csillapítóerőt, amely épp az elem sebességével ellentétes irányú, nagysága pedig a megfelelő erőkomponens adott \(\alpha\)-szorosa [Bagi; 2008]. Az i-edik elemre értelmezett lokális csillapítások általános egyenlete:

\[F_{i}+F_{d_{i}} =m\cdot{}\ddot{u}_{i}\]
\[M_{i}+M_{d_{i}} =\Theta\cdot{}\ddot{\varphi}_{i}\]

ahol \(F_{d_{i}}\) és \(M_{d_{i}}\) lokális csillapítási erők és nyomatékok, számításuk a következő képpen történhet 2D-ben [Ardic; 2006]:

\[F_{x,d_{i}} =-\alpha\cdot{}\textrm{sign}\left(\dot{u}_{x,i}\right)\cdot{}\left|F_{x,i}\right|\]
\[F_{y,d_{i}} =-\alpha\cdot{}\textrm{sign}\left(\dot{u}_{y,i}\right)\cdot{}\left|F_{y,i}\right|\]
\[M_{x,d_{i}} =-\alpha\cdot{}\textrm{sign}\left(\omega_{x,i}\right)\cdot{}\left|M_{x,i}\right|\]
\[M_{y,d_{i}} =-\alpha\cdot{}\textrm{sign}\left(\omega_{y,i}\right)\cdot{}\left|M_{y,i}\right|\]

ahol \(\alpha\) tényező dimenzió nélküli csillapítási tényező. Az \(\alpha\) tényezőt a felhasználó választja meg, alapértelmezett értéke a PFC programokban 0,70. Ez a fajta csillapítás azoknak az elemeknek csökkenti leginkább a gyorsulását, amelyek legkevésbé vannak egyensúlyi állapotban [Bagi; 2008].

Az időintegrálás

Mivel a PFC szoftver explicit időintegrálásos modellel dolgozik, a \(\Delta t\) alatt keletkező elmozdulás-növekményeket a számítás törvényszerűen túlbecsüli – mivel nincsen merevségi mátrix – a pontos megoldáshoz képest (annál jobban, minél nagyobb a \(\Delta t\) intervallumhossz). A kialakuló túlzott elmozdulás-növekmények miatt pedig túl nagy belső erők keletkeznek az elmozdulások irányába, és ezek az erők a következő számítási lépésben mintegy „visszalökik” a rendszert. Ez a jelenség egyfajta oszcilláló mozgást alakít ki, a pontos megoldást jelentő folyamat körül [Bagi; 2008]. A felvázolt kedvezőtlen numerikus jelenséget csökkenthetjük, ha korlátozzuk a felvehető időlépés hosszát (critical time step). A PFC programban egy egyszerű módszerrel becsülik meg a megengedhető időlépés nagyságát, minden számítási ciklusra. A számítás az egy szabadságfokú (ideális lengő) rendszeren (Single Degree-of-Freedom, SDOF) alapul. Az egytömegű ideális lengőrendszer összefüggése a következő:

\[F=m\ddot{x}=-kx\]

A tömegre ható erők összege az alábbi közönséges differenciálegyenletre vezet:

\[m\ddot{x}+kx=0\]

A mozgásegyenlet alapján levezethető kritikus időlépés a következő alakú lesz:

\[\Delta t_{crit} =\frac{T}{\pi}\]
\[T =2\pi\sqrt{\frac{m}{k}}\]

ahol \(m\), \(k\) és \(T\) a rendszer tömege, rúgómerevsége és periódus ideje. A 12. ábrán tömegközéppontok és rugók végtelen sorozatát láthatjuk. Az egyenértékű rendszer kritikus időlépése a fentiek alapján [Behzad; 2012]:

\[\Delta t_{crit}=2\sqrt{\frac{m}{4k}}=\sqrt{\frac{m}{k}}\]

Sorba kötött rugórendszer
12_springsystem

A bemutatott rendszerhez képest a DEM szimulációk a PFC2D és a PFC3D programokban sokkal bonyolultabb szemcse és rugó rendszerekből tevődnek össze, továbbá minden szemcsének és rugónak más és más tömege és merevsége van. Így a \(\Delta t_{crit}\) értéket minden egyes szemcsére külön-külön kell meghatározni szabadságfokonként. Az egész rendszer kritikus időlépése (az időlépés megengedhető maximális hossza) pedig mindezek minimuma lesz:

\[\Delta t_{crit} \le{}\underset{\left(p\right)}{\min}\left\{ \min\left(\sqrt{\frac{m^{p}}{k_{trans}^{p}}},\sqrt{\frac{\Theta^{p}}{k_{rot}^{p}}}\right)\right\}\]

Ebben a kifejezésben \(k_{trans}^{p}\) és \(k_{rot}^{p}\) a \(p\) szemcse legnagyobb eltolódási ill. elfordulási merevsége (ezek a \(p\) elem kapcsolatainak merevségeiből számíthatók), \(\Theta_{p}\) pedig a legnagyobb elfordulási merevség irányához tartozó forgási tehetetlenség [Bagi; 2008].

Mint már említettük az eljárás iterációs jellegű, így a dinamikus folyamatok fokozatosan csengnek le, és – ha az egyáltalán kialakulhat – közelítik az egyensúlyi állapotot. Abszolút értelemben azonban nem beszélhetünk nyugalmi helyzetről egy összetett modell esetén, mert abban mindig marad némi vibrálás. Ezért célszerű egy olyan határértéket definiálni, amit elérve már egyensúlyi állapotként fogadjuk el az adott helyzetet [Tóth; 2004].

A szemcsék alkotta halmaz mozgási állapotot az úgynevezett „átlagos kiegyensúlyozatlan erő” jellemzi, aminek értékét a program maga frissíti lépésről lépésre. Mozgások ugyanis csak a kiegyensúlyozatlan erők hatására jöhetnek létre. Bármely modellnél felvehető egy megfelelő küszöbérték, amit elérvén az „átlagos kiegyensúlyozatlan erő” értéke már nem történnek számottevő mozgások a halmazban, így leáll a futtatás, a vizsgálat szempontjából egyensúlyi helyzet állt be.

Mértékegységek

A PFC programban nincsenek előre megkötött mértékegységek, így azokat önkényesen kell definiálnunk. Szabadon megválasztható az idő (pl. sec), a tömeg (pl. g vagy kg), és a hossz (pl. m vagy cm) mértékegysége, és ezekből már egyértelműen következik, hogy mi lesz az erő, a merevségek, a gravitációs gyorsulás stb. mértékegysége. A mértékegységek tisztázása, egyértelmű megadása azért fontos, mert így nyílhat lehetőség a modell paraméter identifikációjára, melynek köszönhetően a modellezni kívánt szerkezet tényleges adatait vihetjük be a rendszerbe, és így a futtatások eredménye is valóságos értékek lesznek (Tóth, 2004).

Munkakörnyezet

A PFC program használata nem egyszerű, mivel nem grafikus felületen keresztül, hanem parancssorból lehet irányítani. Ennek az az oka, hogy a felhasználóknak teljes mértékben szabad kezet ad a szoftver a szimulációs modell felépítésében. Az időlépéses algoritmuson, valamint a beépített kapcsolati modelleken túl valamennyi folyamatot, az összes szemcse helyét, kapcsolatát, jellemzőit, a gravitációt, stb. a felhasználónak kell definiálnia. Ehhez nagy segítséget jelent a programba beépített úgynevezett „FISH” programozási nyelv, mellyel különféle, a PFC-ben futtatható algoritmusokat lehet megadni [Tóth; 2004]. Egy egyszerű programot mutat be az 13. ábra.

Munka a PFC2D programban
13_pfc2d

Hajlékony útpályaszerkezetek modellezése

A hajlékony útpályaszerkezetek modellezésére mára már számos eljárást kidolgoztak. A legtöbb módszer az útpályaszerkezetet egy rugalmas végtelen féltéren nyugvó többrétegű rendszerként veszi számításba. Az egyik legrégebbi és legelfogadottabb szoftver – ami ezt az elgondolást követi – a SHELL Kutatóközpont által kifejlesztett BISAR (BItumen Stress Analysis in Roads). A programmal feszültséget, megnyúlást és elmozdulást lehet számolni egy függőleges erővel terhelt rugalmas többrétegű rendszerben. A rétegeket a rétegvastagság, a rugalmassági modulus, a Poisson-féle tényező valamint a rétegek határán értelmezett tapadás jellemzi. Az egész rendszert legalul egy végtelen rugalmas féltér támasztja alá. A kérdés az, hogy az analitikus számítások és a DEM modell eredményei hogyan viszonyulnak egymáshoz?

Dondi és mtsai [2007] tanulmányukban pontosan erre próbáltak meg válaszolni. Elkészítették egy hajlékony (aszfalt) útpályaszerkezet DEM modelljét a PFC3D program segítségével. A most következő megállapítások és eredmények az idézet közleményből származnak. A PFC3D programban a modellezett pályaszerkezet rész szélességét és hosszát egyaránt 2 m-re vették fel, a szerkezet teljes vastagsága pedig 58 cm-re adódott. A pályaszerkezet rétegrendje felülről lefele haladva (14. ábra):

  • aszfalt réteg (23 cm),
  • szemcsés alapréteg (30 cm),
  • és földmű (5 cm).
A hajlékony útpályaszerkezet modellje a PFC programban [Dondi és mtsai; 2007]
14_dempave

A numerikus modell paramétereit a laboratóriumi vizsgálatokból vezették le. A rétegeket alkotó szemcsehalmazokat a szemeloszlási vizsgálat alapján generálták. A generált szemcsehalmazok szemeloszlási görbéjét a laboreredmények felskálázása útján nyerték (15. ábra). A szemcsehalmazok felskálázására azért volt szükség, hogy a rétegek szemcseszáma még kezelhető legyen a számítógépek számára, ezzel csökkentve a számítási időt. Az aszfaltréteget és az alapréteget így is 16800 és 11100 darab szemcse építette fel. A földműt mint legalsó réteget, azonos átmerőjű (2,5 cm) szemcsékből állították össze.

Az aszfaltkeverék és a szemcsés réteg szemeloszlási görbéje [Dondi és mtsai; 2007]
15_sieve_analisys

Az aszfaltréteg „rugalmas-viszkózus-képlékeny” anyag. Időfüggő merevsége (az i irányban), a négyparaméteres Burgers-féle modellel írható le (16. ábra):

\[k_{i}=\left[\frac{1}{K_{0}^{i}}+\frac{t}{C_{\infty}^{i}}+\frac{1}{K_{1}^{i}}\left(1-e^{-t/\tau^{i}}\right)\right]^{-1}\]

ahol,
\(t\) terhelési idő,
\(\tau^{i}\) relaxációs idő, τi=Ci1/Ki1,
\(K^{i}_{0}\) a Hooke-test rugalmassági modulusa,
\(C^{i}_{\infty}\) a Newton-test viszkozitása,
\(K^{i}_{1}\) a Voigt-Kelvin test rugalmassági modulusa,
\(C^{i}_{1}\) a Voigt-Kelvin test viszkozitása,
Burgers-féle modell.
16_burgersmodel

A Burgers-féle modellről magyar nyelven Gömze és Kovács [2005] munkájában olvashatunk. Aszfaltkeverékek diszkrét elemes modellezéséről pedig Linbing [2010] könyvében találunk részletesebb ismertetőt. A kapcsolati modell paramétereit az 1. táblázat foglalja össze.

Az aszfaltréteg modellparaméterei, T=20° [Dondi és mtsai; 2007]
Normál merevség kn 6,4*108 MN/m Nyíró merevség ks 6,4*107 MN/m
K0n 1,0*109 MN/m K0s 1,0*108 MN/m
K1n 1,0*108 MN/m K1s 1,0*107 MN/m
Cn 5,0*107 MNs/m Cs 5,0*106 MNs/m
C1n 5,0*106 MNs/m C1s 5,0*105 MNs/m

A modell paramétereit önkényesen úgy vették fel, hogy a laboratóriumi vizsgálattal nyert eredmények (görbék alakja és nagysága) egyezzenek a szimulációból származókkal [Collop és mtsai; 2004].

A szemcsés rétegek (alsó alap és földmű) viselkedését rugalmas kapcsolati modellel írták le a laboratóriumi vizsgálatok alapján (2. táblázat). A szemcsék között fellépő súrlódási tényezőt (\(\mu=0,8\)) a súrlódási szög (\(\varphi=35^{\circ}\)) alapján becsülték meg. Azért, hogy a számítási időt lecsökkentsék, a földmű anyagának sűrűségét a valódinál nagyobb értékkel vették figyelembe. Ez megegyezik azzal, mintha a földmű szintjét egy állandó függőleges nyomásértékkel terhelnénk.

A szemcsés rétegek modellparaméterei [Dondi és mtsai; 2007]
Réteg kn, [MN/m] ks, [MN/m]
Alsó alap 3,6*1010 3,1*109
Földmű 1,6*1010 2,9*109

A pályaszerkezetet két kör alakú – szemcsékből felépített – tárcsán keresztül terhelték. A tárcsák közötti távolság 10 cm volt (ikerabroncs). Mindegyik tárcsát 30 kN erő (\(N\)) terhelte, sugaruk pedig 0,109 m-re adódott az alábbi képlet alapján (17. ábra):

\[R=\left[N/\left(\pi p\right)\right]^{0,5}\]

ahol a \(p\) keréknyomás 8 bar volt.

A szemcsékből álló terhelési felületek [Dondi és mtsai; 2007]
17_dempave-patch

A tanulmányban bemutatott DEM modell segítségével lehetőség nyílik arra, hogy a pályaszerkezet viselkedését dinamikus terhelés hatására vizsgáljuk. A idő figyelembevételévé a pályaszerkezetek leromlási folyamata is jobban megérthető a terhelés nagyságának, sebességének és ismétlődési számának függvényében (pl. keréknyomképződés).

A diszkrét elemes modell felületén átadott terhelés hatására – a szemcsék között – kialakul a kapcsolati erők hálózata. Az időben vizsgálva a folyamatot megfigyelhetjük, hogy a legnagyobb erők a burkolat felületén – az ikerabroncs alatt – jelentkeznek, majd sugárirányban terjedve oszlanak el az egyes rétegekben a mélységgel arányosan. A rendszerben kialakuló elmozdulásokat és feszültségeket a z-x szimmetria síkon a két terhelési felület középvonalában vizsgálták részletesen (18. ábra, 19. ábra). Az aa sík az ikerabroncs és az aszfaltréteg, a bb sík az aszfaltréteg és a szemcsés alapréteg, a cc réteg pedig az alapréteg és a földmű között van értelmezve. A DEM modellel kapott eredményeket a már említett BISAR program segítségével ellenőrizték. A felépített szerkezet BISAR-ban alkalmazott modellparamétereit az 3. táblázat mutatja be. A legutolsó réteg – a végtelen féltér – nagyon magas Young-modulussal és Poisson-féle tényezővel szerepe, hogy a BISAR modell jól egyezzen a DEM modell peremfeltételeivel. A DEM modell esetében az elmozdulásokat és feszültségeket a réteghatároknál nem lehet pontosan meghatározni, ezért a határ feletti és alatti rétegek eredményeit együttesen kell értelmezni.

A BISAR programban használt rétegparaméterei [Dondi és mtsai; 2007]
Réteg Vastagság, [m] Young-modulus, [MPa] Poisson-féle tényező, [-]
Aszfalt 0,23 2200 0,35
Alsó alap 0,30 800 0,45
Földmű 0,05 150 0,47
Végtelen féltér 1012 0,50

A szimuláció lefuttatása alapján elmondható, hogy a DEM és a BISAR modellek hasonló eredményeket szolgáltattak. Az elmozdulások és feszültségek eloszlása mindkét számítási módszernél azonos trendet mutattak. A legjobb egyezést a bb síkon tapasztalták, szemben a cc síkkal. Ezt a szemcsék beékelődéséből származó hatásnak (interlocking effect) tulajdonították. Az alakkal zárás mértéke függ a szemcsék nagyságától és eloszlásától. Mivel a földmű réteg azonos nagyságú szemcsékből áll, ez csökkenti a beékelődés mértékét és így a rétegek együttdolgozását is.

Függőleges elmozdulások a réteghatároknál [Dondi és mtsai; 2007]
18_displacements
Nyomófeszültségek eloszlása a réteghatároknál [Dondi és mtsai; 2007]
19_normalstress-horizontal
Feszültségek a terhelések tengelyében [Dondi és mtsai; 2007]
20_stress_vertical

A pályaszerkezet viselkedését a terhelési felületek tengelyében a z-y síkban is elemezték. A DEM modellből levezett \(\sigma_{zz}\), \(\sigma_{xx}\) és \(\tau_{zx}\) feszültségértékeket az 20. ábra mutatja be. Láthatjuk, hogy ebben az esetben is jó egyezést találunk a BISAR programmal. A tanulmány alapján a következő főbb megállapításokat tehetjük:

  • A pályaszerkezeti rétegeket alkotó szemcsék felskálázásával nagymértékben csökkenthető a számítási idő anélkül, hogy a numerikus eredmények romlanának.
  • Összehasonlítva a BISAR és DEM rendszereket, megállapítható, hogy az adatok jó egyezést mutatnak, így mind kvalitatív és kvantitatív módon is kiértékelhetők az eredmények. Míg a BISAR az elmozdulások és feszültségek „átlagértékével” számol, addig a DEM „csúcsértékekkel” dolgozik, ami jobban megfelel valóságnak.
  • A DEM alkalmas valós jelenségek és így útpályaszerkezetek modellezésére is, lehetőséget biztosít általános következtetések levonására. Viszont komoly nehézséget okoz a megfelelő modellálandók meghatározása. A mikroszkopikus és makroszkopikus paraméterek összevetése a szakirodalomban fellelt és laboratóriumi mérésekkel igazolt eredményekkel nem egyszerű feladat. A probléma csak intenzív kísérleti és kutatási tevékenységgel oldható meg.
  • A DEM módszerek hátránya, hogy nagy számítási kapacitást igényelnek, hosszadalmas folyamat a modellek értékelése és ellenőrzése.

Összefoglalás

A bemutatott tanulmánnyal az elsődleges célunk az volt, hogy a diszkrét elemes modellezés elvét (annak végtelen merev elemeket alkalmazó speciális változatát) az érdeklődő mérnökök számára összefoglaljuk. Az alapelvek ismeretében reményeink szerint bárki képes lesz saját egyszerű merevelemes (BALL-típusú) DEM modelljét megalkotni, akár saját programozás útján is. Ez utóbbi célt komolyan gondolva, elkezdtük egy ingyenes és nyílt forráskódú DEM rendszer alapjait lefektetni „DEMeter” néven, melynek fejlesztésébe bárki belefolyhat, annak megoldásait átveheti. A fejlesztés még kezdeti fázisban van, de amint elkészül az első működékepés változat, azt az Útügyi Lapok oldalain közre fogjuk adni.

Ismertetőnkben kitértünk a hajlékony útpályaszerkezek DEM alapú modellezésére a PFC szoftver felhasználásával. Igaz csak felvázoltunk egy lehetséges megvalósítást, mégis abban bízunk, hogy egyre több kutató mérnök fogja a rugalmas útburkolatok viselkedését ezzel a módszerrel tanulmányozni. Röviden bemutattuk a PFC szoftvert felépítését és működési elvét, remélve azt, hogy segíti az érdeklődőket az elindulásban. Természetesen cikkünk nem térhetett ki mindenre, mivel a téma jelenleg is intenzíven kutatott. Azoknak, akik mélyebb ismeretekre szeretnének szert tenni, azt ajánljuk, hogy a cikkünkhöz felhasznált irodalmakat tanulmányozzák át. Kiemeljük ezek közül Mészöly [1999], Nasztanovics és Füstös [2000], Tóth [2004], Bagi [2008], Fischer és Horvát [2010] magyar nyelvű munkáit.

Köszönetnyilvánítás

Primusz Péter „Szemcserendszerek diszkrételemes modellezése a PFC szoftverrel” publikációt megalapozó kutatása a TÁMOP-4.2.4.A/2-11/1-2012-0001 azonosító számú Nemzeti Kiválóság Program — Hazai hallgatói, illetve kutatói személyi támogatást biztosító rendszer kidolgozása és működtetése konvergencia program című kiemelt projekt keretében zajlott.

Hivatkozások

Ardic, Ö. (2006): Analysis of Bearing Capacity Using Discrete Element Method. Msc Thesis. Middle East Technical University, Ankara, http://etd.lib.metu.edu.tr/upload/12607866/index.pdf (Utolsó letöltés: 2014.12.28.)

Arévalo-Mendoza, G., Ramos-Cañón, A., Prada-Sarmiento, L. (2014): Análisis de confiabilidad en un modelo de descarga de silos de almacenamiento mediante el Método de Elementos Discretos DEM. (Reliability analysis in an unloading model of silo storage by means of the Discrete Element Method DEM), Obras y Proyectos 15, 21-30

Bagi Katalin (2007): A diszkrét elemek módszere, BME Tartószerkezetek Mechanikája Tanszék, Budapest, ISBN 978-963-420-929-4, http://www.epito.bme.hu/me/dolgozok/feltoltesek/kbagi/dem.pdf (Utolsó letöltés: 2014.12.28.)

Behzad M. (2012): Discrete element method applied to the vibration process of coke particles, Master of Science (M.Sc.), Université Laval, p.79 http://www.theses.ulaval.ca/2012/29386/

Bojtár I. és Gáspár Zs. (2003): Végeselemmódszer építőmérnököknek, Terc, Budapest, 339 o.

Collop A.C., Mcdowell G.R. és Lee Y. (2004): Use of the distinct element method to model the deformation behavior of an idealised asphalt mixture, International Journal of Pavement Engineering, n. 5, pp. 1-7

Cozic, L. (2006): 2D Polygon Collision Detection, http://www.codeproject.com/Articles/15573/D-Polygon-Collision-Detection (Utolsó letöltés: 2014.10.26.)

Cundall, P.A. és Strack, O.D.L. (1979): A discrete numerical model for granular assemblies. Geotechnique, Vol. 29 (1), pp. 47-65

Cundall, P.A. és Hart, D.H. (1992): Numerical modelling of discontinua. Journal of Engineering Computations, Vol. 9, pp. 101-113

Czimber K. (2002): Geoinformatika. Elektronikus jegyzet, Sopron, p. 101.

Deresiewicz, H. (1958): Mechanics of granular matter, Advances in Applied Mechanics, volume 5, pp. 233-306

Dondi G., Bragaglia M. és Vignali V. (2007): Flexible pavement simulation with distinct particle element method, 4th International Siiv Congress – Palermo (ITALY), 12-14 september

Dulácska E., Fekete S. és Varga L. (1982): Az altalaj és az építmény kölcsönhatása, Akadémiai Kiadó, Budapest, p. 330

Fischer Sz. és Horvát F. (2010): A georács erősítésű vasúti zúzottkő ágyazat diszkrét elemes modellezési lehetőségei, Közlekedésépítési szemle, (60. évf.) 8. sz. 20-29. old.

Gömze A. L. és Kovács Á. (2005): Aszfaltkeverékek reológiai tulajdonságainak vizsgálata In: Építőanyag : a finomkerámia-, üveg-, cement-, mész-, beton-, tégla- és cserép-, kő- és kavics-, tűzállóanyag-, szigetelőanyag-iparágak szakmai lapja, ISSN 0013-970x , (57. évf.) 2. sz. 34-38. old.

Itasca (2008): PFC3D Version 4.0 Theory and Background. Manual.

Jakob, C. és Konietzky, H. (2012): Particle Methods, An Overview, Technical University Bergakademie Freiberg, p. 23

Kun F. (2011): Számítógépes modellezés ́es szimuláció, Debreceni Egyetem, Elméleti Fizikai Tanszék, Debrecen

Linbing W. (2010): Mechanics of Asphalt, Microstructure and Micromechanics, The McGraw-Hill Companies, Inc, p.464, ISBN: 978-0-07-164097-8

Mészöly T. (1999): Silóvizsgálat mikroszerkezeti modell segítségével, TDK, BME, Tartószerkezetek Mechanikája Tanszék, Budapest

Mindlin, R.D. és Deresiewicz, H. (1953): Elastic Spheres in Contact under Varying Oblique Forces, ASME Journal of Applied Mechanics, Vol. 20, pp. 327-344

Nagy A. (2011): Fejlett grafikai megoldások, egyetemi tananyag, ISBN 978-963-279-516-4

Nasztanovics F., Füstös A. (2000): Lyukkal gyengített tárcsa feszültségeloszlásának numerikus vizsgálata. TDK, BME, Tartószerkezetek Mechanikája Tanszék, Budapest, http://store.naszta.hu/articles/2000_01.pdf (Utolsó letöltés: 2014.12.28.)

O’Sullivan, C. (2011): Particulate Discrete Element Modelling: a Geomechanics Perspective. Spon Press, Taylor&Francis

P.A. Cundall, O.D.L. Strack (1979): A discrete numerical model for granular assemblies. Geotechnique, pages 47-65. DOI 10.1680/geot.1979.29.1.47

Scholtès L. (2009): Modélisation micromécanique des milieux granulaires partiellement saturés. PhD thesis at Institut National Polytechnique de Grenoble.

Scholtès L., Chareyre B., Nicot F., Darve F. (2009): Micromechanics of granular materials with capillary effects. International Journal of Engineering Science (47), pages 64-75. DOI 10.1016/j.ijengsci.2008.07.002

Shäfer, J., Dippel, S. and Wolf, D.E. (1996): Force schemes in simulations of granular materials, Journal de Physique, Number I, Volume 6, pp. 5-20

Szabó A. és Gyenes Cs.: Térfelosztó algoritmusok és térbeli adatszerkezetek, http://goo.gl/G8Wwto (Utolsó letöltés: 2014.10.26.)

Szirmay-Kalos László, Antal György, Csonka Ferenc: Háromdimenziós grafika, animáció és játékfejlesztés ComputerBooks, 2003.

Tóth A. R. (2004): Falazott ívek vizsgálata DEM segítségével, TDK, BME, Tartószerkezetek Mechanikája Tanszék, Budapest, http://goo.gl/LfnTFa (Utolsó letöltés: 2014.10.26.)

Van Baars, S. (1996): Discrete Element Analysis of Granular Materials, Dissertation, University of Technology Delft

Varga I. (2002): Oktató program fejlesztése a „Számítógépes fizika” tantárgyhoz, Debreceni Egyetem, Elméleti Fizikai Tanszék, Debrecen, http://goo.gl/ua3bG6 (Utolsó letöltés: 2014.10.26.)

Wikipedia:

• http://hu.wikipedia.org/wiki/Reológia

• http://hu.wikipedia.org/wiki/Rezgés

• http://hu.wikipedia.org/wiki/Súrlódás

• http://hu.wikipedia.org/wiki/Körmozgás