Józef Droździel, Bogdan Sowiński
Politechnika Warszawska, Wydział Transportu
MODELE OSIADANIA TORU KOLEJOWEGO POD
WPŁYWEM TRWAŁYCH DEFORMACJI PODSYPKI
Rękopis dostarczono, grudzień 2017
Streszczenie: W artykule przedstawiono modele zastępcze toru z uwzględnieniem właściwości
sprężysto-plastycznych podsypki. Cechą charakterystyczną tych właściwości jest występowanie tarcia suchego między przemieszczającymi się elementami tłucznia podczas przejazdu pojazdów szynowych. Dla układów mechanicznych z tarciem suchym można zaobserwować brak powrotu do poprzedniego stanu równowagi na skutek dyssypacji energii mechanicznej w kolejnych cyklach obciążania i odciążania. W modelu uwzględniono kumulację strat energii mechanicznej w poszczególnych symulacjach przejazdu pojazdu po torze oraz zmianę parametrów sztywności i tłumienia podsypki. Zaproponowano także algorytm rozwiązywania układu równań z tarciem suchym i założeń dla procedury numerycznej realizującej taki algorytm. W końcowej części zaprezentowano wyniki symulacji przebiegu osiadania toru.
Słowa kluczowe: tor, osiadanie, modelowanie
1. WPROWADZENIE
W badaniu dynamiki drogi kolejowej można wyróżnić dwa podukłady. Pierwszy z nich to nawierzchnia kolejowa złożona z szyn, podkładów, przytwierdzeń szyn do podkładów i warstwy podsypki. Drugi zawiera szeregowo ułożone: warstwę filtracyjną (warstwę piasku lub pospółki), podtorze, i podłoże (grunt). Drugi podukład jest skrótowo dość często nazywany podtorzem. Pierwszy podukład przekazuje obciążenie na drugi. Jeżeli chodzi o pierwszy podukład to w zasadzie do analizy osiadania toru wystarczająca jest znajomość geometrii ułożenia toru, właściwości mechanicznych podsypki i początkowych nierówności geometrycznych szyn –generujących drgania w układzie.
Analiza warstw podtorza jest, zdaniem wielu autorów, podstawą do oceny jego osiadania i narastania nierówności geometrycznych toru kolejowego [7], [10], [12]. W pracy Salima [7] są zaprezentowane wyniki eksperymentu i pewne zagadnienia modelowania matematycznego w postaci równań konstytutywnych. Wee Loon Lim [12] analizował wyniki wielu doświadczeń materiałowych i pomiarów. Natomiast w pracy Simona [10] szerzej niż w [7] i [12] powiązano problem osiadania toru z ruchem pociągów. Nie pomniejszając znaczenia warstw podtorza w procesie osiadania toru, autorzy tego artykułu skupili się na przybliżonej ocenie wpływu warstwy posypki, wykazującej cechy sprężysto – plastyczne (tarcie suche) i charakteryzującej się nieliniową zależnością przenoszonych sił
od odkształceń, Z innych badań [9], [10] wynika, że warstwy podtorza większość obciążeń przenoszą w stanie sprężystym.
Proponowane modele obejmują pionowe oddziaływania dynamiczne w układzie pojazd szynowy - tor podczas ruchu po torze prostym o stałych parametrach (masy zastępcze, współczynniki sprężystości i tłumienia). Pominięte zostały specyficzne odcinki toru jak mosty, tunele, estakady, rozjazdy, przejazdy kolejowo – drogowe.
Jak wiadomo, osiadanie toru i powstawanie nierówności geometrycznych ma istotny wpływ na stan jego utrzymania. W szczególności powoduje to zmiany początkowych wymiarów, sztywności i tłumienia podsypki oraz jej zdolności do przenoszenia obciążeń na podtorze. Procesy te zależą od intensywności obciążania toru ruchem pociągów i odbywają się w różnym tempie. Eksperymentalne modele degradacji (osiadania i powstawania nierówności) toru, rozwinięte zostały przez Fröhlinga [5], Sato [8] i Shentona [9], a szeroki przegląd modeli osiadania toru zawarty jest w pracy Dahlberga [2].
W dostępnych autorom publikacjach brak jest jednak opisu zjawiska degradacji toru za pomocą odpowiednich modeli fenomenologicznych (ujmujących makroskopowo badane układy, nie wnikających szczegółowo w ich budowę). W tym artykule przedstawiono kilka modeli matematycznych długookresowej predykcji zachowania toru sprzężonego z prostym modelem pojazdu. Pokazano w nim wyniki prognozowania osiadania toru poprzez szereg symulacji. Treść artykułu jest oparta na rezultatach uzyskanych w projekcie badawczym nr 4 T12C 006 30 [11] zrealizowanym przez autorów, które po części opublikowano w pracach [1] i [3].
2. MODELE TORU Z TARCIEM SUCHYM W PODSYPCE
Przyjęto założenie, że odkształcenie warstwy podsypki po wpływem obciążenia toru może być częściowo nieodwracalne (trwałe) po ustaniu tego obciążenia. Rozmiar tego trwałego odkształcenia zależy od stanu utrzymania podsypki i przebiegu obciążania toru w czasie eksploatacji Wykonując szereg symulacji przejazdu pojazdu po badanym odcinku toru, można oszacować przebieg osiadania i zmian nierówności geometrycznych tego odcinka. Przyjęto, że parametry opisujące właściwości toru jak masy zastępcze szyn, podkładów, podsypki, parametry sztywności i tłumienia (tarcia) podsypki są stałe na badanym odcinku.Niezależnie od wariantu modelu uwzględniono w nim tarcie suche. Jak wiadomo, siła tarcia suchego jest nieciągłą funkcją prędkości względnego przesunięcia stykających się powierzchni i jest często upraszczana poprzez zastąpienie jej funkcją ciągłą, np. w postaci arcus tangens. Chcąc opisać przebieg osiadania toru nie można charakterystyki tarcia zastąpić funkcją ciągłą, bowiem przy tarciu nierozwiniętym (brak poślizgu) siła tarcia może przyjmować w pewnym zakresie różne wartości. Należy, wobec tego znaleźć efektywny algorytm rozwiązywania układu z masą zastępczą reprezentującą szyny i podkłady na sprężynie o charakterystyce liniowej i równoległym do niej tłumikiem ciernym, w którym siła tarcia będzie miała określoną wartość w obszarze jej wielowartościowości. Najprostszym modelem fizycznym, który może być brany pod uwagę jest model toru o jednym stopniu swobody stanowiący podparcie każdego koła pojazdu lub każdego zestawu kołowego (model płaski toru) i poruszający się z prędkością pojazdu (rys. 1).
Rys. 1. Elementarny model sprężysto-plastyczny równoległy (model 1.)
Bardziej złożonym wariantem jest układ modeli elementarnych sprzężonych ze sobą wzdłuż toru (rys. 2). Może on być traktowany jako model toru sprowadzony do jego linii środkowej (osi toru) lub jako model nie sprzężonych ze sobą toków szynowych. Jest to model płaski.
Rys. 2. Model toru ze sprzężonymi modelami elementarnymi wzdłuż jego linii środkowej lub toku szynowego (model 2.)
Model przestrzenny, w którym występuje sprzężenie między tokami szynowymi poprzez sztywne podkłady jest pokazany na rys. 3.
Rys. 3. Model przestrzenny toru z podkładami zamodelowanymi jako bryły sztywne (model 3.) Rozwiązanie układu równań ruchu dla modelu 1. i modelu 2. można oprzeć na rozwią-zaniu dla modelu elementarnego. Dla modelu 3. trzeba bazować na układzie o dwóch stopniach swobody, złożonego z podkładu jako sztywnej belki podpartej na końcach sprężynami i tłumikami ciernymi jak na rysunku. W modelu 3. podparcie belki stanowią elementy równoległe sprężysto-plastyczne jak w modelu 1.
Równanie ruchu i jego rozwiązanie dla modelu 1. (elementarnego) można przedstawić następująco: ) t ( f F kz z m T (1) gdzie:
z – przemieszczenie pionowe masy zastępczej,
ݖሷ – przyspieszenie pionowe masy zastępczej,
m k F
f(t) – siła wymuszająca (nacisk pionowy koła lub zestawu kołowego), FT – siła tarcia (rys. 4).
0 0 z ) z ( csign z c , c F FT T
Rys. 4. Przebieg siły tarcia w funkcji prędkości poślizgu według modelu Coulomba Podstawiając odpowiednio ) ( ) ( t f kz F z sign c F S T
równanie (1) można sprowadzić do układu równań stanu w postaci (2)
2 1 2 1 2 z z m / ) F F ( m / ) z ( sign c ) t ( f kz z S T (2) Jeśli z2 = 0 i lim FT c z20 lub lim FT c z20, to wówczas mogą zachodzić dwa przypadki: 1) FS >0 i FT < 0, przy czym FS + FT < 0,
2) FS <0 i FT > 0, przy czym FS + FT > 0
Gdy zachodzi którykolwiek z pokazanych wyżej przypadków, to do dalszych obliczeń przyjmujemy FT FS i ruch układu ustaje, gdyż z10 i z20. Przyjęcie innego założenia oznaczałoby doprowadzanie energii do układu przez siły tarcia.
Przypadek trzeci zachodzi, gdy z2 0, ale wtedy siła tarcia jest jednoznacznie określona.
Przykładowy rezultat zastosowania algorytmu pokazano na rys. 5. Wynik rozwiązania zilustrowany na rys. 5 jest testem modelu o parametrach nie związanych z właściwościami podsypki. Obrazuje on jedynie jak zachowuje się układ z tarciem suchym w warunkach wymuszenia okresowego, w którym charakterystyczne są fazy „zatrzymania” przemieszczeń.
Układ o dwóch stopniach swobody, stanowiący podstawę rozwiązania modelu 3. pokazano na rys. 6.
Składa się on ze sztywnej belki (podkładu z masami zastępczymi szyn i podsypki) o masie m i momencie bezwładności Jx względem osi wzdłużnej x, podpartego przez liniowe
FT
ݖሶ
c
elementy sprężyste i tłumiki cierne. Na belkę działają dwie równoległe siły FL(t) i FR(t), reprezentujące oddziaływania z kołami zestawu kołowego pojazdu.
Rys. 5. Przykładowy wynik zastosowania algorytmu dla układu o jednym stopniu swobody z tarciem suchym
Rys. 6. Układ o dwóch stopniach swobody jako element składowy modelu 3
Prawoskrętny układ współrzędnych kartezjańskich przyjęto, tak aby oś x była równoległa do osi toru. Zwrot tej osi jest zgodny ze zwrotem wektora prędkości pojazdu.
Równania ruchu tego układu są następujące:
) ( )] ( ) ( [ 2 )] ( ) ( [ 2 2 R L R L x R L R L F F a f b b k J F F f z k z m (3) przy czym
0 1 0 1 1 0 1 R , L R , L R , L R , L dla dla , dla ) ( (3a)
Lzb, R zb; (3b) są prędkościami poślizgu względnego powierzchni ciernych.
Siłę tarcia T i moment tarcia M, występujące w równaniach (3) oznaczamy:
( ) ( ) ( ) ( ) L R L R L R T T T f f M b f ) ( ) L) (( R L L ( L)))) ) ( ) L R L))) ((( (3c)Charakterystyki sił tarcia (TL, R) można przedstawić przy założeniu, że współczynnik tarcia statycznego jest większy niż współczynnik tarcia kinetycznego (o) lub są one sobie równe (o=). Występująca w układzie równań (3) stała f może przyjmować wartości zgodnie z formułą , 0 , 0 0 L R L R f dla f f dla , 0 L R, ,R 0 , , (3d)
W uproszczonej wersji charakterystyki sił tarcia przyjmuje się f = f0.
Mając na celu wyznaczenie sił i momentów sił w warunkach tarcia nierozwiniętego przyjęto, że L,R 0. W pracach autorów [3], [11] wyszczególniono przypadki jakie mogą zachodzić w układzie równań (3) dla sił i momentów tarcia w zależności od prędkości poślizgu względnego stykających się elementów.
Wymuszenie ruchu układu opisanego równaniami (3) stanowią siły pionowe nacisków kół na szyny przejeżdżającego pojazdu. W ruchu wymuszonym układu stan, w którym względna prędkość poślizgu L,R 0 może być stanem chwilowym lub trwać przez pewien skończony czas. W każdym z tych przypadków, siły tarcia powinny mieć takie wartości, aby spełniać w kolejnych chwilach czasu dynamiczne równania ruchu.
Z formuły L,R0 wynikają zależności:
0 0 6 b z C b z b z L L (3e) 0 0 6 b z C b z b z R R (3f)
Stałe CL i CR oznaczają przemieszczenia względne, dla których L,R 0.
Podstawiając odpowiednio (3e) i (3f) do równań (3) przy warunku L,0 lub R0,
R L R L R , L R , L R , L R L R L R , L R , L R , L F F f kz b a F b a F kC m m T F F f kz b a F b a F kC m m T ! " # $ % & ' ( ) * % & ' ( ) * ! " # $ % & ' ( ) * % & ' ( ) * 2 1 1 2 2 1 1 2 2 1 (4)
Siła T1L,R występuje dla L,0, natomiast T2L,R dla R0
Jeśli równocześnie L,0 i R0, to siły tarcia wynoszą:
% & ' ( ) * % & ' ( ) * + b a F b a F kC T L R R , L R , L 1B 2 1 2 (4a)
Występująca w równaniach (4) masa m – jest masą zastępczą szyn i podkładu w ruchu pionowym (bez obrotu), zaś masy zastępcze mL i oraz mR występują przy jednoczesnym ruchu pionowym i obrotowym i są równe:
2 b J m m x R , L +
Siły tarcia występujące w powyższych równaniach muszą spełniać warunki:
0 0,|T | f f | T | iL,R L,R dla i = 1, 2
Procedurę realizującą algorytm rozwiązania oparto na metodzie Rungego-Kutty czwartego rzędu ze zmiennym krokiem całkowania. Należy ją zaliczyć do tzw. metod „event driven” (sterowanie programem obiektowym za pomocą zdarzeń, np. zmiana znaku siły tarcia) rozwiązujących układy Filippova [4], [6]. Procedura całkuje układ równań ze stałym krokiem, aż do momentu zmiany znaku siły tarcia. Następuje wówczas zmiana kroku całkowania, aż do uzyskania takiego kroku, który dawałby z założoną dokładnością zero argumentu (
n i i iq a 1
) funkcji siły tarcia. W ogólnym przypadku zakłada się, że siła tarcia jest w postaci
0 ) ( 0 , 1 1 1 n i i i n i i i i n i i i T T q a q a csign q a c c F F
gdzie: qi - prędkości uogólnione układu.
Jest to podział przestrzeni stanu na podprzestrzenie, w których rozwiązanie jest wystarczająco gładkie do zastosowania do niego standardowej metody rozwiązywania równań różniczkowych zwyczajnych. Po wyznaczeniu zgodnie z podanym algorytmem siły tarcia procedura kontynuuje obliczenia z ponownie zwiększonym krokiem całkowania.
Praktyka wykazała, że przy znanej dokładności zera i znanym kroku dotychczasowego rozwiązania wystarczy do 3 iteracji, by „trafić w zero”.
3. SYMULACJA OSIADANIA TORU
Symulowano jazdę wagonu po torze prostym zakładając, że siły FL = FR. Dane przyjęte do obliczeń są następujące:
rozstaw szyn 2a = 1,5 m,
zastępcze punkty podparcia 2b = 1,5 m, masa zastępcza podkładu m = 400 kg,
zastępczy moment bezwładności podkładu Ix = 280 kgm2, siły wymuszenia pionowego FL i FR przyjęto w postaci
! "# $ % & ' ( ) * 2 2 1 2 , , l vt sin P F FL R
obciążenie przypadające na oś wagonu P = 4105-8105 N, rozstaw osi wózka l = 2,5 m,
prędkości wagonu v =10-50 m/s.
Symulację przeprowadzono za pomocą oprogramowania Matlab - Simulink. Schemat blokowy symulacji układu z elementami sprężysto – plastycznymi reprezentującymi podsypkę pokazano na rys. 7. Zastosowano symulację krokową, to znaczy każde uruchomienie obliczeń wiąże się z odpowiedzią na pojedynczy cykl obciążenia w postaci odkształcenia. Powstałe odkształcenie jest następnie wartością początkową dla następnego cyklu obliczeń.
W pracy [8] analizowano także symulacje układu pojazd – tor dla innych elementarnych modeli reologicznych opisujących dyssypację energii mechanicznej w podsypce, mianowicie dla: modelu sprężysto - plastycznego zmodyfikowanego (modelu sprężysto-plastycznego połączonego szeregowo ze sprężyną liniową), modelu Prandtla i zmodyfikowanego modelu Prandtla. Na podstawie oceny wyników otrzymanych z tych symulacji, przyjęto do porównania model sprężysto-plastyczny i sprężysto-plastyczny zmodyfikowany ze wskazaniem na ten ostatni, uzależniając siłę tarcia od siły nacisku osiowego.
Porównanie wyników symulacji dla ruchu płaskiego układu tor-pojazd z zastosowaniem modeli sprężysto-plastycznego i sprężysto-plastycznego zmodyfikowanego [11] przedstawiono na rys. 8 (na osi odciętych liczba cykli oznacza liczbę przejazdów zestawów kołowych). W obliczeniach przyjęto wartość zastępczego współczynnika tarcia = 0,2 oraz
jednakowe sztywności wszystkich sprężyn występujących w tych modelach k = 50 MN/m. Uzyskane z symulacji krzywe osiadania dla tych modeli można ocenić jako zbliżone do siebie, przy czym obliczenia z użyciem modelu sprężysto - plastycznego zmodyfi-kowanego są bardziej czasochłonne.
Rys. 7. Symulacja modelu z elementami sprężysto – plastycznymi w środowisku Matlab – Simulink; schemat blokowy
Prognozowanie osiadania wykonano przy różnych wartościach sztywności podparcia i różnych wartościach obciążenia statycznego na oś zestawu. Do symulacji zastosowano model, w którym pod każdym zestawem kołowym znajduje się poruszający się wraz z nim model zastępczy toru pokazany na rys. 6. Obciążenia modelowano jako siły skupione przesuwające się wzdłuż osi toru ze stałą prędkością. Symulacje przeprowadzono dla 40 000 przejazdów zestawów kołowych, co odpowiada 10 000 przejazdów 4-osiowego wagonu. Brane były pod uwagę trzy warianty statycznego obciążenia na oś: 60 kN, 100 kN i 150 kN oraz zastępczej sztywności pionowej toru 50 MN/m i 100 MN/m.
Rys. 8. Wyniki osiadania uzyskane dla modelu sprężysto-plastycznego (a) i sprężysto-plastycznego ze wzmocnieniem (b) przyspieszenie przemieszczenie predkosc Przem To Workspace1 Sila To Workspace Sign 1 s Integrator2 1 s Integrator1 1 s Integrator Fr Fl Add2 Add1 Add -K-2k/m -K-2f/m -K-1/m
Właściwości tłumienia w postaci siły tarcia zostały określone przez kąt tarcia wewnętrznego przy spoistości tłucznia równej zero i przyjęto do obliczeń = 0,2. Tor w badanych modelach składa się z szyn UIC60 podpartych na podkładach betonowych.
Symulacje przeprowadzono na odcinku toru prostego o długości 120 m co wymagało wprowadzenia do modelu 201 nieliniowych elementów podparcia podkładów (co 0,6 m), a osiadanie toru zredukowano do jego linii środkowej (płaski model toru). Porównanie wyników otrzymanych za pomocą własnego modelu z wynikami według zależności empirycznych Fröhlinga [5], odniesionych do osi toru, pokazano na rys. 9.
Rys. 9. Porównanie wyników uzyskanych na podstawie modelu Fröhlinga (linie cienkie) i modelu z podparciem sprężysto-plastycznym ze wzmocnieniem (linie pogrubione)
W obu modelach przyjęto elementy sprężysto-plastyczne zmodyfikowane ze współ-czynnikiem tarcia = 0,2 oraz sztywnością pionową podsypki k = 50 MN/m.
Przebiegi krzywych osiadania toru (rys. 9) uzyskanych z obliczeń symulacyjnych i zależności empirycznych opracowanych przez Fröhlinga [5] wskazują, że wartości osiadania, uzyskane na podstawie własnego modelu są na ogół większe od wartości uzyskanych z zależności empirycznych. Krzywe empiryczne wykazują mniejszą wrażliwość na wzrost liczby przejazdów Analizując otrzymane wyniki można zaobserwować, że wzrost obciążenia osiowego powoduje wzrost osiadania toru. Ponadto, na podstawie modelu badanego w pracy [11] pokazano, że wzrost sztywności podsypki znacząco zmniejsza wielkość osiadania, natomiast wzrost zastępczego współczynnika tarcia powoduje zwiększenie osiadania toru. Należy zaznaczyć, że opracowany model wykazuje większą wrażliwość na zmianę obciążenia osiowego niż na zmianę sztywności podsypki.
4. PODSUMOWANIE
W badanym modelu, przyczyną powodującą osiadanie toru jest tarcie suche w podsypce. Oczywiście istnieją także inne przyczyny wpływające na ten proces (np. podtorze), które nie zostały uwzględnione w proponowanych modelach.
0 10000 20000 30000 40000 liczba przejazdow zestawow
0.000 0.005 0.010 0.015 0.020 osiadanie [m] obciazenie na os 15*10^4 N obciazenie na os 10^5 N obciazenie na os 6*10^4 N obciazenie na os 15*10^4 N (Frohling) obciazenie na os 10^5 N (Frohling) obciazenie na os 6*10^4 N (Frohling)
Ze względu na brak możliwości wykonania odpowiedniego eksperymentu, ograniczono się na tym etapie do obliczeń teoretycznych i porównania ich z wynikami opartymi na równaniach empirycznych zawartych w pracy [5].
Zaproponowany model predykcji osiadania toru jest zgodnie z zamiarem autorów próbą prezentacji alternatywnego opisu tego zjawiska w stosunku do innych prac. Jest to model fenomenologiczny, który zawiera też pewne cechy modelu strukturalnego, ale nie opisuje wewnętrznego mechanizmu osiadania opartego na skomplikowanych równaniach naprężeń i odkształceń. Jest on oparty na stosunkowo prostym opisie dyssypacji energii mechanicznej spowodowanej tarciem suchym i powstawaniu trwałych odkształceń powodujących osiadanie toru.
Bibliografia
1. Chudzikiewicz A., Droździel J., Sowiński B.: Mathematical Model of Track Settlement Caused by Dry Friction in the Ballast. The Archives of Transport, vol. XXI, No 3-4, Warsaw 2009, pp. 25-38.
2. Dahlberg T.: Some Railroad Settlement Models—a Critical Review, Proc. Institution of Mechanical Engineers, Vol. 215, Part F, 2001, pp. 289-300.
3. Droździel J., Sowiński B.: Simulation of a railway track deterioration influenced by ballast stiffness with dry friction: Computers in Railways XI, Series: WIT Transactions on The Built Environment, Vol 103, Southampton, Boston, 2008, ISBN 978-1-84564-126-9, ISSN 1746-4498, pp. 693-702.
4. Filippov A.F.: Differential Equations with Discontinuous Righthand Sides. Kluwer Academic Publishers, 1988.
5. Fröhling R.D.: Low Frequency Dynamic Vehicle/Track Interaction: Modelling and Simulation, Vehicle System Dynamics Supplement Vol. 28 (1998), pp. 30-46.
6. Piroinen P.T.: Kuznetsov Y.A.: An Event-Driven Method to Simulate Filippov Systems. ACM Transactions on Mathematical Software, vol. 34, no 3 pp. 13.1-13.23, 2008.
7. Salim W.: Deformation and Degradation Aspects of Ballast and Constitutive Modelling under Cyclic Loading. The University of Wollongong, PhD thesis, 2004 (AU).
8. Sato Y.: Japanese Studies on Deterioration of Ballasted Track, Vehicle System Dynamics Vol. 24 (1995), pp. 197-208.
9. Shenton M.J.: Ballast Deformation and Track Deterioration, Track Technology. Thomas Telford Ltd. London. 1985.
10. Simon A., Spiker J.: The Mechanical Behaviour of Ballasted Railway Tracks. Delft University of Technology, The Netherlands, 2002.
11. Sowiński B. i inni: Predykcja długookresowych zmian parametrów toru kolejowego powstających w procesie eksploatacji. Raport z projektu badawczego MNiSzW nr 4 T12C 006 30, Politechnika Warszawska, Wydział Transportu, Warszawa 2008.
12. Wee Loon Lim: Mechanics of Railway Ballast Behaviour. The University of Nottingham, PhD Thesis, 2004 (UK).
MODELS OF TRACK SETTLEMENT INFLUENCED BY DURABLE BALLAST DEFORMATIONS
Summary: The article presents track equivalent models considering the elastic and plastic properties of the
ballast. A characteristic feature of these properties is the occurrence of dry friction between moving gravel elements during the passage of rail vehicles. For mechanical systems with dry friction, no return to previous equilibrium can be observed because of dissipation of mechanical energy in subsequent cycles of loading and unloading those systems. The model also considers the accumulation of mechanical energy losses in different simulated rail vehicle passage, and the change of stiffness and damping ballast parameters. An algorithm for solving a system of equations with dry friction and assumptions for a numerical procedure implementing such an algorithm was also proposed. In the final part, the results of the simulation of the course of track settlement were presented.
