Index of /rozprawy2/11331

Pełen tekst

(1)AKADEMIA GÓRNICZÓ-HUTNICZA im. Stanisława Staszica w Krakowie Wydział Inzynierii Mechanicznej i Robotyki Katedra Mechaniki i Wibroakustyki. PRACA DÓKTÓRSKA ANALIZA WŁASCIWÓSCI DYNAMICZNYCH PRZENÓSNIKÓW WIBRACYJNYCH DZIAŁAJĄCYCH W ÓPARCIU Ó ZJAWISKÓ DYNAMICZNEJ ELIMINACJI DRGAN mgr inż. Marek Gajowy Promotor: prof. dr hab. inż. Jerzy Michalczyk Kraków 2017 1.

(2) Spis treści: 1. WPROWADZENIE – CEL PRACY ..................................................................................... 5 2. Zasada eliminacji dynamicznej drgań (Frahm 1909) [1] ....................................................... 7 2.1. Model dynamicznego eliminatora drgań [2] ................................................................... 7 2.2. Szerokość kotliny antyrezonansowej ............................................................................. 10 3. MODEL MATEMATYCZNY PRZENOŚNIKA ANTYREZONANSOWEGO................ 12 3.1. Równania ruchu ............................................................................................................. 12 3.2. Siły kontaktowe przy zderzeniach – model nadawy (Michalczyk J. 2008) .................. 15 3.3. Model silnika asynchronicznego ................................................................................... 18 3.4. Dobór parametrów – strojenie maszyny ........................................................................ 19 4. PROGRAM SYMULACYJNY ........................................................................................... 21 4.1. Krok całkowania ............................................................................................................ 21 4.2. Procedura przeprowadzania symulacji .......................................................................... 21 5. WYZNACZENIE PODSTAWOWYCH PARAMETROW DYNAMICZNYCH .............. 22 5.1. Symulacja przenośnika nieobciążonego nadawą ........................................................... 22 5.2. Symulacja przenośnika obciążonego nadawą ................................................................ 25 5.2.1. Prędkość nadawy w funkcji ilości nadawy na przenośniku .................................... 29 5.2.2. Wydajność w funkcji ilości nadawy ....................................................................... 30 5.2.3. Moc w funkcji wydajności ...................................................................................... 31 5.2.4. Siły przekazywane na podłoże w funkcji wydajności ............................................ 34 6. ANALIZA MODALNA MODELOWANEGO PRZENOŚNIKA ...................................... 35 6.1 Linearyzacja układu nieliniowego .................................................................................. 35 6.2 Częstości drgań własnych przenośnika .......................................................................... 36 6.3 Formy drgań własnych ................................................................................................... 37 6.4 Wnioski ........................................................................................................................... 39 7. WYZNACZENIE AMPLITUD DRGAŃ MAKSYMALNYCH PODCZAS WYBIEGU . 40 7.1. Klasyczne metody wyznaczania amplitud w rezonansie przejściowym ....................... 40 7.1.1. Przyczyny błędów metod klasycznych. .................................................................. 45 7.1.2. Metoda nomograficzna (Cieplok 2009) [25, 26] .................................................... 46 7.2. Metoda bilansu energetycznego I [13,29] ..................................................................... 49 7.3. Zmodyfikowana metoda bilansu energetycznego ......................................................... 54 7.4. Wyniki i wnioski ........................................................................................................... 56 2.

(3) 8. SAMOSYNCHRONIZACJA I ROZSYNCHRONIZOWANIE PRZENOŚNIKA ............ 57 8.1. Warunek Blechmana na wystąpienie samosynchronizacji ............................................ 57 8.2. Ocena stabilności ruchu synchronicznego wibratorów z kryterium Blechmana ........... 58 8.3. Badania symulacyjne ..................................................................................................... 64 8.3.1. Stabilność ruchu w przypadku zróżnicowanych oporów na łożyskach silników .. 64 8.3.2. Stabilność synchronizacji przy zmianie kąta rozstawu wibratorów. ...................... 66 8.3.3. Wrażliwość synchronizacji od zmiany odległości wibratorów od środka ciężkości korpusu .............................................................................................................................. 67 8.3.4. Wpływ ilości nadawy na synchronizację. ............................................................... 67 8.4. Wnioski .......................................................................................................................... 69 9. WARIANTY KONSTRUKCYJNE – ZMIENNY STOSUNEK MASY RYNNY DO MASY KORPUSU ................................................................................................................... 70 9.1. Wstęp ............................................................................................................................. 70 9.2. Analiza porównawcza właściwości przenośnika przy stałej masie korpusu i zmiennej masie rynny........................................................................................................................... 70 9.2.1. Prędkość transportowania nadawy .......................................................................... 71 9.2.2. Wydajność transportowania w funkcji ilości nadawy na przenośniku ................... 72 9.2.3. Stosunek zapotrzebowania mocy do wydajności w funkcji wydajności ................ 73 9.2.4. Sumaryczna siła przekazywana na fundament........................................................ 74 9.2.5. Porównanie szerokości kotliny antyrezonansowej przy stałym stosunku masy rynny 𝑀𝑟 do współczynnika sprężystości względnej 𝑘𝑓 ................................................. 75 9.2.6. Porównanie amplitud w rezonansie przejściowym przy stałym stosunku masy rynny 𝑀𝑟 do współczynnika sprężystości względnej 𝑘𝑓 ................................................ 76 9.3. Analiza porównawcza właściwości przenośnika przy stałej masie rynny i zmiennej masie korpusu. ...................................................................................................................... 77 9.3.1. Prędkość transportowania nadawy .......................................................................... 78 9.3.2. Wydajność transportowania w funkcji ilości nadawy na przenośniku ................... 79 9.3.3. Stosunek zapotrzebowania mocy do wydajności w funkcji wydajności ................ 80 9.3.4. Sumaryczna siła przekazywana na fundament........................................................ 80 10. PODSUMOWANIE ........................................................................................................... 82 11. LITERATURA ................................................................................................................... 84. 3.

(4) Wykaz oznaczeń: ky. – współczynnik sprężystości wzdłużna pakietu sprężyn śrubowych,. kx. – współczynnik sprężystości poprzeczna pakietu sprężyn śrubowych,. kf. – współczynnik sprężystości poprzeczna pakietu resorów,. ke. – współczynnik sprężystości,.  , wym ,. – częstość wymuszająca,. n , i. – częstość drgań własnych. e. – mimośród,. m1,2,w. – masa niewyważona,. Mr,. – masa rynny,. M. – macierz mas,. K. – macierz sztywności, współczynnik sprężystosci,. Mk,M. – masa korpusu, masa chroniona,. Jk. – główny centralny moment bezwładności korpusu,. Jr. – główny centralny moment bezwładności rynny,. J1 ,J2. – główny centralny moment bezwładności wibratora,. A0 ,A,B. – amplituda drgań,. F0. – amplituda siły wymuszającej,. 𝑥, 𝑦, 𝑓. – przemieszcznie liniowe,. 𝛼, 𝜑1 , 𝜑2. – przemieszcznie kątowe,. 𝑘𝑝. – współczynnik podrzutu,. 𝑃. – moc silnika asynchronicznego,. 𝑀𝑢𝑡. – moment utyku silnika napędowego,. 𝑀𝑒𝑙. – moment elektryczny silnika napędowego,. 𝑞. – wektor przemieszczeń,. R. – współczynnik restytucji,. p. – wykładnik Hertza Sztajermana,. 𝑠. – poślizg silnika asynchroniczmego,. 𝐴𝑟. – amplituda drgań rynny. Δ𝜑. – różnica faz wibratorów. 𝑄. – wydajność. 𝑅𝑚𝑎𝑥. – siła sumaryczna na podłoże. 𝐸, 𝑉. – energia kinetyczna, energia potencjalna 4.

(5) 1. WPROWADZENIE – CEL PRACY Przenośniki wibracyjne posiadają szereg cech przydatnych przy wykorzystaniu ich dla transportu bliskiego materiałów sypkich w różnych gałęziach przemysłu. Można tu wymienić: odporność na wysoką, rzędu 1000℃, temperaturę, możliwość odbioru ciepła materiału transportowanego, możliwość izolowanego transportu materiałów wydzielających trujące gazy i pyły , możliwość schładzania lub nagrzewu nadawy, prowadzenia reakcji chemicznych itp. W chwili obecnej dużym zainteresowaniem na rynku przenośników wibracyjnych cieszą się nowe konstrukcje, których schemat dynamiczny oparty jest na zjawisku eliminacji dynamicznej. Przenośniki te, według danych producentów posiadają szereg zalet, m. innymi wyższe amplitudy drgań przy tym samym napędzie i wyjątkowo niskie wartości sił dynamicznych przekazywanych na podłoże. Poza patentami konstrukcji przenośników opartych na zasadzie eliminatora Frahma, np.[3,4] istnieje stosunkowo niewiele prac poświęconych dynamice tych maszyn, np.[5,6,23]. W szczególności brak jest prac badających zachowanie tego typu maszyn podczas zmian ich obciążenia nadawą, lub prace te zawierają błędy [35] , jak np. [8]. Ponadto prezentowane konstrukcje należą głównie do przenośników o małych gabarytach i wydajności. Celem pracy jest analiza właściwości eksploatacyjnych przenośników wibracyjnych, których zasada działania oparta jest na zjawisku eliminacji dynamicznej Frahma, w szczególności ocena wpływu obciążenia przenośnika nadawą na przebieg i wydajność procesu roboczego, reakcje dynamiczne na podłoże oraz amplitudy drgań w stanach przejściowych. W pracy zbudowano model symulacyjny takiego przenośnika obciążonego nadawą i przeprowadzono badania wydajności transportowania oraz sił przekazywanych na podłoże w stanach pracy ustalonej i amplitud maksymalnych drgań w okresach rozruchu i wybiegu. W pracy zajęto się również zjawiskiem samosynchronizacji wibratorów [11, 24] występującym w działaniu tych urządzeń. Skupiono się w tym przypadku na badaniu wrażliwości samosynchronizacji na zmiany parametrów geometrycznych zamocowania wibratorów oraz nierównomierność obciążenia silników. Ważnym elementem rozprawy stało się dopracowanie metody szacowania amplitudy w rezonansie przejściowym podczas wybiegu maszyny o wielu stopniach swobody na bazie bilansu energetycznego [13]. Jest to metoda posiadająca znaczną przewagę nad klasycznymi sposobami wyznaczania amplitud podczas wybiegu, bazującymi na konieczności wyznaczenia przyspieszenia wirujących wibratorów, co stanowi poważny problem i źródło dużych błędów we wspomnianych metodach [29]. Dzięki zmodyfikowanej metodzie bilansu energetycznego 5.

(6) można uzyskać wartości amplitud z błędem względnym rzędu od 2 do 36%, co stanowi zdecydowaną przewagę nad metodami klasycznymi, gdzie błąd oszacowania przekracza często kilkaset procent.. 6.

(7) 2. Zasada eliminacji dynamicznej drgań (Frahm 1909) [1] 2.1. Model dynamicznego eliminatora drgań [2] Rozpatrzmy dynamikę układu chronionego {𝑀, 𝐾, 𝐶}, którego przyczyną ruchu jest siła harmoniczna 𝐹𝑜 e𝑖𝜈𝑡 , z eliminatorem dynamicznym {𝑚𝑒 , 𝑘𝑒 , 𝑐𝑒 }, przedstawionego schematycznie na rysunku 1.. Rys. 1. Eliminator dynamiczny. Równania ruchu układu przedstawiają zależności (1): 𝑀𝑥̈ + 𝐶𝑥̇ + 𝐾𝑥 + 𝑐𝑒 (𝑥̇ − 𝑦̇ ) + 𝑘𝑒 (𝑥 − 𝑦) = 𝐹𝑜 e𝑖𝜈𝑡 𝑚𝑒 𝑦̈ − 𝑐𝑒 (𝑥̇ − 𝑦̇ ) − 𝑘𝑒 (𝑥 − 𝑦) = 0. (1) Rozwiązanie układu równań możemy zapisać w postaci: 𝑥 = 𝐴e𝑖𝜈𝑡 , 𝑥̇ = 𝑖𝜈𝐴e𝑖𝜈𝑡 , 𝑥̈ = −𝜈 2 𝐴e𝑖𝜈𝑡 𝑦 = 𝐵e𝑖𝜈𝑡 , 𝑦̇ = 𝑖𝜈𝐵e𝑖𝜈𝑡 , 𝑦̈ = −𝜈 2 𝐵e𝑖𝜈𝑡. (2) gdzie A i B – amplitudy zależne od parametrów dynamicznych układu i parametrów wymuszenia. Po podstawieniu poszukiwanych rozwiązań (2) do (1) otrzymano: 𝑀(−𝜈 2 𝐴e𝑖𝜈𝑡 ) + 𝐶𝑖𝜈𝐴e𝑖𝜈𝑡 + 𝐾𝐴e𝑖𝜈𝑡 + 𝑐𝑒 (𝑖𝜈𝐴e𝑖𝜈𝑡 − 𝑖𝜈𝐵e𝑖𝜈𝑡 ) + 𝑘𝑒 (𝐴e𝑖𝜈𝑡 − 𝐵e𝑖𝜈𝑡 ) = 𝐹𝑜 e𝑖𝜈𝑡 𝑚𝑒 (−𝜈 2 𝐵e𝑖𝜈𝑡 ) − 𝑐𝑒 (𝑖𝜈𝐴e𝑖𝜈𝑡 − 𝑖𝜈𝐵e𝑖𝜈𝑡 ) − 𝑘𝑒 (𝐴e𝑖𝜈𝑡 − 𝐵e𝑖𝜈𝑡 ) = 0. (3) 7.

(8) Równanie powyższe można przedstawić w postaci macierzowej. [. 𝑖𝜈𝑡 𝑖𝜈𝑡 −𝑀𝜈 2 + (𝐶 + 𝑐𝑒 )𝑖𝜈 + (𝐾 + 𝑘𝑒 ) −𝑐𝑒 𝑖𝜈 − 𝑘𝑒 ] ∙ [ 𝐴e𝑖𝜈𝑡 ] = [ 𝐹𝑜 e ] 2 −𝑐𝑒 𝑖𝜈 − 𝑘𝑒 −𝑚𝑒 𝜈 + 𝑐𝑒 𝑖𝜈 + 𝑘𝑒 𝐵e 0. (4) Jako rozwiązanie otrzymano amplitudy zespolone A i B:. 𝐴=−. 𝐾𝑚𝑒. 𝑣2. − 𝐾𝑘𝑒 + 𝑀𝑘𝑒. 𝑣2. − 𝑀𝑚𝑒. 𝑣4. 𝐹𝑜 𝑐𝑒 𝑖𝑣 − 𝐹𝑜 𝑚𝑒 𝑣 2 + 𝐹𝑜 𝑘𝑒 + 𝑘𝑒 𝑚𝑒 𝑣 2 + 𝑐𝑒 𝑖𝑚𝑒 𝑣 3 − 𝐶𝑐𝑒 𝑖 2 𝑣 2 − 𝐶𝑖𝑘𝑒 𝑣 − 𝐾𝑐𝑒 𝑖𝑣 + 𝐶𝑖𝑚𝑒 𝑣 3 + 𝑀𝑐𝑒 𝑖𝑣 3. (5). 𝐵=−. 𝐹𝑜 𝑐𝑒 𝑖𝑣 + 𝐹𝑘𝑒 𝐾𝑚𝑒 𝑣 2 − 𝐾𝑘𝑒 + 𝑀𝑘𝑒 𝑣 2 − 𝑀𝑚𝑒 𝑣 4 + 𝑘𝑒 𝑚𝑒 𝑣 2 + 𝑐𝑒 𝑖𝑚𝑒 𝑣 3 − 𝐶𝑐𝑒 𝑖 2 𝑣 2 − 𝐶𝑖𝑘𝑒 𝑣 − 𝐾𝑐𝑒 𝑖𝑣 + 𝐶𝑖𝑚𝑒 𝑣 3 + 𝑀𝑐𝑒 𝑖𝑣 3. (6) Przykładowo przebiegi amplitud drgań: układu chronionego bez eliminatora A0, układu chronionego z eliminatorem A oraz eliminatora drgań B dla wybranych parametrów układu i wymuszenia pokazano ma rysunku 2.. Rys.2. Zależność bezwymiarowej amplitudy drgań mas układu w funkcji bezwymiarowej częstości wymuszającej [2]. 𝐹1 =. 𝛾=. 𝐹𝑜 = 1, 𝑀𝑔. 𝑐𝑒 = 0.05, 𝐶. 𝜉=. 𝐶 2√𝐾𝑀. 𝑍1 =. 𝑍 , 𝑀𝑔 𝐾. = 0.1,. 𝜇=. 𝑚𝑒 = 0.1, 𝑀. 𝑍 = {|𝐴0 |, | 𝐴|, |𝐵|},. 𝜀=. 𝛿=. 𝑘𝑒 = 0.12, 𝐾. 𝑣 𝐾 , 𝜔𝑜 = √ 𝜔𝑜 𝑀. (7). Dla małych tłumień w układzie 𝑐𝑒 = 𝐶 ≅ 0 amplitudy drgań A0, A i B można zapisać 8.

(9) 𝐴0 = 𝐹0. 𝐴 = 𝐹0 |. 𝐵 = 𝐹0 |. 1 , |(𝐾 − 𝑀𝑣 2 )|. (8). 𝑘𝑒 − 𝑚𝑒 𝑣 2. |,. (9). |,. (10). (𝐾 + 𝑘𝑒 − 𝑀𝑣 2 ) ∙ (𝑘𝑒 − 𝑚𝑒 𝑣 2 ) − 𝑘𝑒 2 𝑘𝑒. (𝐾 + 𝑘𝑒 − 𝑀𝑣 2 ) ∙ (𝑘𝑒 − 𝑚𝑒 𝑣 2 ) − 𝑘𝑒 2. Charakterystyki amplitudowo-częstotliwościowe przedstawiono na rysunku. Rys.3. Charakterystyki amplitudowo-częstotliwościowe dla układu nietłumionego [2]. Z postaci rozwiązań zilustrowanych na rysunku 3., wynika, że w przypadku gdy: 𝑘𝑒 𝑣=√ 𝑚𝑒. (11). amplituda drgań A masy chronionej M będzie miała wartość zerową (krzywa czerwona). Warto przy tym zauważyć, że amplituda drgań masy 𝑚𝑒 ma wtedy wartość:. 𝐵=. 𝐹0 𝑘𝑒. (12) 𝑘𝑒. Warunek 𝑣 = √. 𝑚𝑒. jest warunkiem dynamicznej eliminacji drgań za pomocą dołączonego 𝐹. dodatkowego układu mechanicznego – eliminatora dynamicznego. Zależność 𝐵 = 𝑘0. 𝑒. 9.

(10) pokazuje, że w warunkach eliminacji drgań układ dołączony nie znajduje się w stanie drgań rezonansowych. Porównując przebiegi amplitud drgań masy M z eliminatorem przedstawione na rysunkach 2 i 3 łatwo zauważyć, że dla tłumień różnych od zera amplituda A osiąga 𝑘. minimum dla 𝑣 ≈ √𝑚𝑒 . Stąd przy założeniu, że masa eliminatora drgań 𝑚𝑒 jest stała i jest 𝑒. zwykle znacznie mniejsza od masy układu głównego. 𝑚𝑒 𝑀. ≤ 0.1 możemy wyznaczyć optymalną. sztywność eliminatora ke.. 2.2. Szerokość kotliny antyrezonansowej Istotne znaczenie dla bezpieczeństwa pracy układu posiada szerokość przedziału pomiędzy częstościami rezonansowymi, zwana kotliną antyrezonansową. Aby określić jej szerokość z rozwiązania układu liniowego nietłumionego (1) wyznaczono amplitudę dla dwóch współrzędnych x i y w postaci zależności (9) i (10). Częstości drgań własnych można otrzymać poprzez przyrównanie mianownika w powyższych równaniach do zera. (𝐾 + 𝑘𝑒 − 𝑀𝑣 2 ) ∙ (𝑘𝑒 − 𝑚𝑒 𝑣 2 ) − 𝑘𝑒 2 = 0. (13). dokonując wcześniej przekształceń i podstawień a więc Dokonując dzielenia równania powyższego przez 𝑘𝑒 otrzymujemy następującą postać (𝐾 + 𝑘𝑒 − 𝑀𝑣 2 ) ∙ (𝑘𝑒 − 𝑚𝑒 𝑣 2 ) − 𝑘𝑒 2 = 0 | ∶ 𝑘𝑒 (𝐾 + 𝑘𝑒 − 𝑀𝑣 2 ) ∙ (. (14). 𝑘𝑒 𝑣2 𝑘𝑒 2 − )− =0 𝑘𝑒 𝑘𝑒 𝑘𝑒 𝑚𝑒. (15). po uproszczeniu otrzymano postać 𝑘𝑒 𝑣2 (𝐾 + 𝑘𝑒 − 𝑀𝑣 2 ) ∙ ( − ) − 𝑘𝑒 = 0 | ∶ 𝐾 𝑘𝑒 𝑘𝑒 𝑚𝑒. (16). Dokonując kolejnego podzielenia równania tym razem przez K otrzymano 𝐾. (𝐾 +. 𝑘𝑒 𝐾. −. 𝑀𝑣 2 𝐾. 𝑘. ) ∙ (𝑘𝑒 − 𝑒. 𝑣2 𝑘𝑒 𝑚𝑒. )−. 𝑘𝑒 𝐾. =0. (17). co po uproszczeniu i przekształceniu dało. 10.

(11) (1 +. 𝑘𝑒 𝑣2 𝑣2 𝑘𝑒 − ) ∙ (1 − )− =0 𝐾 𝑘𝑒 𝐾 𝐾 𝑀 𝑚𝑒. (18). Jak udowodniono w literaturze [1], sprężystość K nie wpływa na eliminację w przypadku takiego układu mechanicznego, co pozwala nam przyjąć ten sam stosunek 𝐾. = 𝜔𝑛2 lub też przekształcając 𝑀 Definiując stosunek mas (1 + 𝜇 −. 𝑚𝑒 𝑀. 𝑘𝑒 𝐾. =. 𝑘𝑒 𝑚𝑒. =. 𝑚𝑒 𝑀. = 𝜇 otrzymujemy ostatecznie. 𝑣2 𝑣2 ) ∙ (1 − )−𝜇 =0 𝜔𝑛2 𝜔𝑛2. (19). Pierwiastkami tego równania są wartości. 𝑣 2 𝜇 𝜇2 √ ( ) = (1 + ) ∓ 𝜇 + 𝜔𝑛 2 4. (20). Uzyskano zależność szerokości kotliny antyrezonansowej jedynie od stosunku mas, pokazaną na wykresie na rys. 4. 2 1.8 1.6 1.4. 𝑣 / 𝜔n𝑛. 1.2 1 0.8 0.6 0.4 0.2 0. 0. 0.05. 0.1. 0.15. 0.2. 0.25. 0.3. 0.35. 0.4. 0.45. 0.5. 𝜇. Rys. 4. Wykres przebiegu dwóch częstości rezonansowych. 𝑣 𝜔𝑛. analizowanego układu od stosunku mas. 𝑚𝑒 𝑀. =. 𝜇. Jak widać na rysunku 4. stosunek masy eliminatora do masy głównej nie może być zbyt niski, gdyż powoduje to zmniejszenie kotliny antyrezonansowej a więc układ może być narażony na zwiększone amplitudy podczas pracy z nierównomiernym obciążeniem lub innymi zakłóceniami mogącymi przesunąć punkt pracy w częstość bliską częstości rezonansowej.. 11.

(12) 3. MODEL MATEMATYCZNY PRZENOŚNIKA ANTYREZONANSOWEGO 3.1. Równania ruchu W pracy poddano analizie antyrezonansowy przenośnik wibracyjny o konstrukcji pokazanej na rysunku 5a. Na jego podstawie zbudowano model przenośnika wibracyjnego z 6 stopniami swobody (𝑥, 𝑦, 𝛼, 𝑓, 𝜑1 , 𝜑2 ) oraz 3-kolumnowej i 4-warstwowej nadawy (12 elementów). Przyjęty schemat dynamiczny przenośnika i nadawy pokazano na rys. 5b. Przenośnik składa się z dwóch mas, tzw. masy korpusu 𝑀𝑘 , masy rynny 𝑀𝑟 , która w tym modelu pełni rolę eliminatora dynamicznego oraz dwóch wibratorów inercyjnych z indywidualnym napędem indukcyjnym. Korpus został podparty sprężyście na sprężynach śrubowych przymocowanych do podłoża, natomiast rynna przenośnika na listwach resorujących zamocowanych do korpusu. Wibratory bezwładnościowe biegnące przeciwbieżnie są źródłem wypadkowej siły wymuszającej działającej na korpus na kierunku drgań.. Rys. 5a. Schemat konstrukcyjny analizowanego przenośnika antyrezonansowego.. 12.

(13) Rys.5b. Model obliczeniowy przenośnika wibracyjnego. 13.

(14) Wibratory zostały zamocowane w taki sposób, aby symetralna odcinka łączącego środki przegubów wibratorów przecinała środek masy korpusu oraz rynny. W badaniach symulacyjnych przyjęto następujące wartości stałych geometrycznych i masowych: Wymiary geometryczne przenośnika: 𝑙 = 8 𝑚 – długość czynna rynny (odległość transportowania), 𝐿 = 2 𝑚, 𝐻 = 0,48 𝑚, 𝐿1 = 0,97 𝑚, ℎ1 = 0,33 𝑚, 𝐿2 = 0,77𝑚, ℎ2 = 0,67𝑚, ℎ𝑟 = 1,1𝑚, 𝐿𝑟 = 1,92 𝑚, 𝑒 = 0,06𝑚 – długość mimośrodu (przyjęta z warunku, aby współczynnik podrzutu wynosił 𝑘𝑝 = 3,3. Masy i momenty bezwładności: 𝑀𝑟 = 1000 𝑘𝑔, 𝑀𝑘 = 2500 𝑘𝑔 – odpowiednio masa rynny i masa korpusu; Σ𝑚𝑛 = 12 ÷ 3000 𝑘𝑔 – masa (zmienna) nadawy; 𝑚1,2 = 50 𝑘𝑔 – masy niewyważone wibratorów; 𝐽𝑘 = 12200 𝑘𝑔𝑚2 , 𝐽𝑟 = 5000 𝑘𝑔𝑚2 , 𝐽1 = 𝐽2 = 0,18 𝑘𝑔𝑚2 –. odpowiednio. centralne. momenty bezwładności korpusu, rynny i masy niewyważonej wibratora, Parametry elementów sprężystych: 𝑘𝑓 = 10962000 𝑁/𝑚, 𝑘𝑦 = 2328000 𝑁/𝑚, 𝑘𝑥 = 0.5𝑘𝑦 , 𝑏𝑥 = 230,14. 𝑁𝑠 𝑁𝑠 𝑁𝑠 , 𝑏𝑦 = 460,28 , 𝑏𝑓 = 666,87 𝑚 𝑚 𝑚. Parametry nadawy: 𝑅 = 0,01 – współczynnik restytucji zderzeń pomiędzy kolejnymi warstwami nadawy, 𝜇𝑝𝑛 = 0,4 – współczynnik tarcia między pierwszą warstwą nadawy a rynną, 𝜇𝑛𝑛 = 0,7 – współczynnik tarcia wewnętrznego dla nadawy, 1. 𝑑 = 3 𝑙 – długość pojedynczej warstwy nadawy Parametry dynamiczne silnika asynchronicznego: 𝑃 = 4 𝑘𝑊 – moc silnika asynchronicznego, 𝑀𝑢𝑡 = 97,4 𝑁𝑚 – moment utyku silnika napędowego, 𝜔𝑠 = 104,7 𝑟𝑎𝑑/𝑠,. 𝜔𝑢𝑡 = 79,6 𝑟𝑎𝑑/𝑠. –. synchronicznego i utyku. 14. odpowiednio. prędkość. kątowa. biegu.

(15) Do analizy ruchu układu przyjęto w stanie równowagi statycznej maszyny bez nadawy układ bezwzględny centralny względem korpusu maszyny o osiach 𝑥, 𝑦, kąt obrotu korpusu 𝛼, przemieszczenie względne rynny względem korpusu 𝑓 oraz bezwględne kąty obrotu wibratorów 𝜑1 i 𝜑2 . Zespół równań opisujący ruch maszyny, wyprowadzony metodą równań Lagrange’a II rodzaju, zapisać można w postaci macierzowej: (21). [𝑀] ∙ [𝑞̈ ] = [𝑄] gdzie. 𝑞̈ =. 𝑑2 𝑑𝑡 2. [𝑥, 𝑦, 𝛼, 𝑓, 𝜑1 , 𝜑2 ]𝑇. (22). Macierz mas M przedstawia formuła (23) natomiast wektor wyrazów wolnych Q zaprezentowano w postaci (24).. 3.2. Siły kontaktowe przy zderzeniach – model nadawy (Michalczyk J. 2008) W trakcie pracy przenośnika wibracyjnego dochodzi do licznych zderzeń pomiędzy nadawą oraz rynną przenośnika. Podczas kolizji powstają bardzo duże siły kontaktowe, które są opisane według zmodyfikowanego modelu Hertza, zaproponowanego w pracy [9]. Proces przebiegu siły kontaktowej w funkcji odkształcenia zderzanych brył δ oraz w zależności od fazy zderzenia został zobrazowany na rys. 6. Pole powierzchni ograniczone krzywymi {1} i {2}, jest miarą strat energii związanych m.in. z odkształceniem plastycznym zderzanych brył, wydzieleniem się ciepła i kruszeniem lokalnych zabrudzeń. Straty te zwykle opisuje się parametrem nazywanym współczynnikiem restytucji R.. Rys.6. Siły kontaktowe w funkcji odkształcenia w dwóch fazach zderzenia. W przyjętym modelu nadawę podzielono na kolumny i warstwy zderzające się podczas ruchu przenośnika jedynie z warstwami bezpośrednio stykającymi się w pionie. 15.

(16) 𝑴=. 𝑚1 + 𝑚2 + 𝑀𝑘 + 𝑀𝑟. 0. ℎ1 𝑚1 + ℎ2 𝑚2 − ℎ𝑟 𝑀𝑟. 𝑀𝑟 𝑐𝑜𝑠𝛽. 𝑚1 cos(𝛽 + 𝜑1 ) 𝑒. 𝑚2 cos(𝛽 − 𝜑2 ) 𝑒. -. 𝑚1 + 𝑚2 + 𝑀𝑘 + 𝑀𝑟. −𝐿1 𝑚1 − 𝐿2 𝑚2 + 𝐿𝑟 𝑀𝑟. 𝑀𝑟 𝑠𝑖𝑛𝛽. 𝑚1 sin(𝛽 + 𝜑1 ) 𝑒. 𝑚2 sin(𝛽 − 𝜑2 ) 𝑒. -. -. 𝐽𝑘 + 𝐽𝑟 + 𝐿21 𝑚1 + 𝐿22 𝑚2 + 𝐿2𝑟 𝑀𝑟 + + ℎ12 𝑚1 + ℎ22 𝑚2 + ℎ𝑟2 𝑀𝑟. 𝐿𝑟 𝑀𝑟 𝑠𝑖𝑛𝛽 − ℎ𝑟 𝑀𝑟 𝑐𝑜𝑠𝛽. −𝐿1 𝑚1 𝑠𝑖𝑛(𝛽 + 𝜑1 )𝑒 + +ℎ1 𝑚1 cos(𝛽 + 𝜑1 ) 𝑒. −𝐿2 𝑚2 𝑠𝑖𝑛(𝛽 − 𝜑2 )𝑒 + +ℎ2 𝑚2 cos(𝛽 − 𝜑2 ) 𝑒. -. -. -. 𝑀𝑟. 0. 0. -. -. -. -. 𝐽1 + 𝑚1 𝑒 2. 0. -. -. -. -. -. 𝐽2 + 𝑚2 𝑒 2. (23). 𝑚1 𝑠𝑖𝑛(𝛽 + 𝜑1 )𝜑̇ 12 𝑒 − 𝑚2 sin(𝛽 − 𝜑2 ) 𝜑̇ 22 𝑒 − 𝑘𝑥 (𝑥 + 𝐻𝛼) − 𝑏𝑥 (𝑥̇ + 𝐻𝛼̇ ) − 𝑇101 − 𝑇102 − 𝑇103 𝑚2 𝑐𝑜𝑠(𝛽 − 𝜑2 )𝜑̇ 22 𝑒 − 𝑚1 𝑐𝑜𝑠(𝛽 + 𝜑1 ) 𝜑̇ 12 𝑒 − 𝑘𝑦 𝑦 −. 𝑸=. 𝑏𝑦 (𝑦̇ − 𝐿𝛼̇ ) 𝑏𝑦 (𝑦̇ + 𝐿𝛼̇ ) − − 𝐹101 − 𝐹102 − 𝐹103 2 2. 𝐿1 𝑚1 𝑐𝑜𝑠(𝛽 + 𝜑1 ) 𝜑̇ 12 𝑒 − ℎ2 𝑚2 sin(𝛽 − 𝜑2 ) 𝜑̇ 22 𝑒−𝐿2 𝑚2 𝑐𝑜𝑠(𝛽 − 𝜑2 )𝜑̇ 22 𝑒 + ℎ1 𝑚1 𝑠𝑖𝑛(𝛽 + 𝜑1 )𝜑̇ 12 𝑒+ 𝑏𝑦 (𝑦̇ + 𝐿𝛼̇ )𝐿 𝑏𝑦 (𝑦̇ − 𝐿𝛼̇ )𝐿 𝑘𝑦 (𝑦 + 𝐿𝛼)𝐿 𝑘𝑦 (𝑦 − 𝐿𝛼)𝐿 −𝑏𝑥 (𝑥̇ + 𝐻𝛼̇ )𝐻 − + − 𝑘𝑥 (𝑥 + 𝐻𝛼)𝐻 − + + 𝑇101 ℎ𝑟 + 𝑇102 ℎ𝑟 + 𝑇103 ℎ𝑟 − 𝐹102 𝑑 − 𝐹103 2𝑑 2 2 2 2 −𝑏𝑠 𝑓̇−𝑘𝑠 𝑓 −𝑚1 𝑔 sin(𝛽 + 𝜑1 ) 𝑒 + 𝑀𝑒𝑙1 + 𝑀01 −𝑚2 𝑔 sin(𝛽 − 𝜑2 ) 𝑒 + 𝑀𝑒𝑙2 + 𝑀02. (24). 16.

(17) Analityczny opis siły kontaktowej przedstawionej graficznie na rys. 6 przedstawia wzór (25).  1  R2  F j ,( j 1,k )  ( j 1,k   j ,k ) p k 1  1  sgn( j 1,k   j ,k ) sgn( j 1,k   j ,k )  2  . . . (25). k – stała sprężystości kontaktowej brył, ( j 1,k  j ,k ) – łączna deformacja brył, ( j 1,k  j ,k ) – prędkość względna brył, R – współczynnik restytucji, p – wykładnik Hertza Sztajermana zależny od kształtu powierzchni brył w punkcie zderzenia. Siła tarcia działająca na poszczególne warstwy nadawy na kierunku poziomym opisuje równanie (26). T j ,( j 1,k )  Fj ,( j 1,k ) sgn( j ,k   j 1,k ). (26). Fj,(j-1,k) – składowa normalna nacisku warstwy j-tej na j-1 w kolumnie k-tej, Tj(j-1,k) – składowa styczna nacisku warstwy j-tej na j-1 w kolumnie k-tej będąca siłą tarcia. Jeśli kolejne warstwy nadawy (w danej kolumnie) j i j-1 nie stykają się, wówczas siła kontaktowa na kierunku normalnym Fj,(j-1,k) i stycznym Tj,(j-1,k) pomiędzy tymi warstwami jest równa zero. Równania ruchu na kierunkach x i y poszczególnych warstw nadawy, z uwzględnieniem oddziaływania przenośnika na dolne warstwy nadawy, mają postać:. mnj ,k   T j ,( j 1,k )  T j 1,( j ,k ). (27). mnj ,k  mnj ,k g  Fj ,( j 1,k )  Fj 1,( j ,k ) Cechą charakterystyczną przyjętego opisu sił kontaktowych jest możliwość opisania zjawiska zderzenia jako procesu dynamicznego przy zachowaniu wartości Newtonowskiego współczynnika restytucji R określonego w modelu "bezczasowym".. 17.

(18) 3.3. Model silnika asynchronicznego W budowie maszyn wibracyjnych najpowszechniej stosowane są silniki indukcyjne. Jednym z najprostszych modeli matematycznych indukcyjnych silników elektrycznych jest model Klossa a więc charakterystyka statyczna silnika (28), opisująca moment elektromechaniczny na wale silnika asynchronicznego. Funkcję tę charakteryzuje zależność tylko od jednej zmiennej – prędkości kątowej wirnika.. 𝑀𝑒𝑙 (𝜔) =. 2𝑀𝑢𝑡 (𝜔𝑠 − 𝜔𝑢𝑡 )(𝜔𝑠 − 𝜔) (𝜔𝑠 − 𝜔𝑢𝑡 )2 + (𝜔𝑠 − 𝜔)2. (28). Moc nominalną silnika oraz pozostałe stałe wyznaczono z poniższych wzorów:. 𝜔𝑠 =. 2𝜋 ∙ 50 𝑙𝑝𝑏. 𝑀𝑛𝑜𝑚 =. 𝑃𝑛𝑜𝑚 𝜔𝑛𝑜𝑚. 𝑀𝑢𝑡 = 𝑀𝑛𝑜𝑚 ∙ 𝑝. 𝑠𝑛𝑜𝑚 =. 𝜔𝑠 − 𝜔𝑛𝑜𝑚 𝜔𝑠. 𝑠𝑢𝑡 = 𝑠𝑛𝑜𝑚 (𝑝 + √𝑝2 − 1). 𝑠𝑢𝑡 =. 𝜔𝑠 − 𝜔𝑢𝑡 𝜔𝑠. gdzie:. 𝑀𝑛𝑜𝑚. – moment znamionowy silnika elektrycznego. ω. – prędkość kątowa wirnika. p. – przeciążalność. 𝜔𝑠 , 𝜔𝑢𝑡 , 𝜔𝑛𝑜𝑚 – prędkość wirnika synchroniczna, utyku, nominalna 𝑠𝑢𝑡 , 𝑠𝑛𝑜𝑚. – poślizg utyku, nominalny. lpb. – liczba par biegunów. Charakterystyka momentowa silnika klatkowego różni się znacząco od charakterystyki opisanej wzorem Klossa tylko w stanach szybkozmiennych np. podczas rozruchu, jednak jak wykazano w pracy [30] zastosowanie charakterystyki statycznej, zamiast dynamicznej wiążącej 18.

(19) się z koniecznością wprowadzenia kolejnych stopni swobody i wydłużenia czasu symulacji, dało zadowalające wyniki również w stanie rozruchu co skłoniło autora do pozostania przy charakterystyce statycznej związanej ze wzorem Klossa. Funkcję Klossa można zmodyfikować uzupełniając ją o czynnik poprawiający charakterystykę momentu w fazie rozruchu. 𝑀𝑒𝑙 (𝜔)∗ = 𝑀𝑒𝑙 (𝜔) + 𝑀𝑝 (𝜔). (29). 𝑀𝑒𝑙 (𝜔) – podstawowy wzór Klossa Człon poprawkowy 𝑀𝑝 (𝜔) przyjąć można np.: 𝜔. 𝑀𝑝 (𝜔) = 𝑀𝑝0 ∙ 𝑒 −Ω. (30). Stałe 𝑀𝑝0 , Ω funkcji uzupełniającej mogą zostać tak dobrane, aby charakterystyka momentowa poprawiona (29) odpowiadała wartościom katalogowym silników klatkowych.. 3.4. Dobór parametrów – strojenie maszyny W badanym modelu przenośnika wibracyjnego eliminator dynamiczny stanowi dołączona do masy korpusu masa rynny Mr. Warunkiem zaistnienia eliminacji drgań korpusu, 𝑘. a więc jej pracy w antyrezonansie jest spełnienie zależności 𝜔𝑤𝑦𝑚 = √𝑀𝑓 , gdzie 𝜔𝑤𝑦𝑚 = 𝑟. 104,7. 𝑟𝑎𝑑 𝑠. – częstość siły wymuszającej generowanej przez wibrator bezwładnościowy ( prędkość. kątowa rotoru silnika) Wyprowadzono konieczny współczynnik sprężystości przedstawiony w (31): 2 𝑘𝑓 = 𝑀𝑟 ∙ 𝜔𝑤𝑦𝑚 = 10962000 𝑁/𝑚. (31). Podparcie korpusu zapewniono przez sprężyny śrubowe o wartościach stałych sprężystości: 𝑘𝑦 = 2328000 𝑁/𝑚, 𝑘𝑥 = 0.5𝑘𝑦 .. 19.

(20) Wykorzystano również współczynnik rozproszenia energii materiału podparcia sprężystego odpowiednio maszyny i rynny 𝜓1 = 0.13, 𝜓2 = 0.04 a dzięki nim obliczono ekwiwalentny współczynnik tłumienia wiskotycznego (32) 𝜓 𝑘𝑥. 𝑏𝑥 = 2𝜋𝜔1. 𝑤𝑦𝑚. 𝜓 𝑘𝑦. , 𝑏𝑦 = 2𝜋𝜔1. 𝑤𝑦𝑚. 𝜓 𝑘𝑓. , 𝑏𝑓 = 2𝜋𝜔2. (32). 𝑤𝑦𝑚. 20.

(21) 4. PROGRAM SYMULACYJNY 4.1. Krok całkowania Do przeprowadzania symulacji i zarządzania nimi wykorzystano dwa środowiska programistyczne: Matlab i Pascal. Język programowania pakietu Matlab jest językiem programowania wysokiego poziomu, zawiera wiele gotowych bibliotek, które ułatwiają użytkownikowi rozwiązywanie problemów z wielu dziedzin nauki. Pascal jest natomiast językiem niższego poziomu w stosunku do Matlaba. Programy napisane w językach niskopoziomowych są zwykle lepiej zoptymalizowane, działają szybciej, ale napisanie ich wymaga od programisty pewnych umiejętności i czasu. W przypadku programu symulacyjnego badającego ruch przenośnika nieobciążonego nadawą posłużono się programem Matlab a w nim dostarczoną funkcją ode45 korzystającą z algorytmu Runghego-Kutty IV rzędu ze zmiennym krokiem całkowania układu równań różniczkowych. Metoda ta stanowi dużą przewagę nad stałokrokowym odpowiednikiem napisanym w programie Pascal z racji na szybsze zwracanie wyniku. Program napisany w Pascalu ze stałokrokowym algorytmem R-K IV-rzędu stanowi przewagę w przypadku symulacji układu obciążonego nadawą. Obliczanie przemieszczeń i prędkości w przypadku następujących po sobie zderzeń nadawy i przenośnika (pojawianie się krótkotrwałych dużych wartości sił kontaktowych) generuje w Matlabie wykorzystującym funkcję ode45 problem w postaci wydłużonego czasu ustalania odpowiedniego kroku całkowania przez co program staje się niestabilny wydłużając czas oczekiwania na rozwiązanie.. 4.2. Procedura przeprowadzania symulacji Przebieg procesu przeprowadzania symulacji i analizy wyników został podzielony na trzy etapy. Można je porównać z etapami analizy wytrzymałościowej w programach typu MES. Są to: • preprocessing – przygotowywanie modelu matematycznego do obliczeń numerycznych • solving – obliczenia numeryczne • postprocessing – wizualizacja i analiza wyników obliczeń Preprocessing, solving oraz postprocessing są realizowane w przypadku przenośnika nieobciążonego nadawą w środowisku Matlab, natomiast do obliczeń numerycznych dla przenośnika obciążonego nadawą użyto programu napisanego w języku PASCAL. Do rysowania wykresów użyto programu GnuPlot. 21.

(22) 5. WYZNACZENIE PODSTAWOWYCH PARAMETROW DYNAMICZNYCH 5.1. Symulacja przenośnika nieobciążonego nadawą W pierwszej kolejności przeprowadzono analizę działania przenośnika nieobciążonego nadawą. Wzrost amplitudy drgań na początku rozruchu został wygenerowany przez przechodzenie układu przez częstości rezonansowe. Amplituda przemieszczeń współrzędnych 𝑥, 𝑦 oraz 𝛼 wraz z upływem czasu zmierzała do wartości bliskich zeru co wskazuje na eliminacje dynamiczną drgań, gdyż wibratory wirują będąc zamocowane na korpusie. Masa rynny i odpowiednio nastrojone resory zapewniają stale działającą siłę równoważącą siłę generowaną przez niewyważone masy wirujące. Rynna uzyskuje w stanie ustalonym amplitudę drgań rzędu 𝐴𝑟(𝑠𝑦𝑚𝑢𝑙𝑎𝑐𝑗𝑎) = |𝑓(t) + 𝑥(t)𝑐𝑜𝑠𝛽 + 𝑦(t)sin𝛽|𝑚𝑎𝑥 = 5,87 ∙ 10−3 𝑚. Wartość ta jest porównywalna z wynikiem oszacowania amplitudy drgań rynny z zasady zachowania środka masy, którą można stosować w przypadku braku lub niewielkich sił działających z zewnątrz na układ punktów materialnych, tzn. przy założeniu dostatecznie „miękkiego” układu podarcia korpusu. Masa korpusu praktycznie unieruchomiona nie jest składową wzoru wynikłego z zasady zachowania środka masy. 𝐴𝑟(𝑠𝑧𝑎𝑐𝑜𝑤𝑎𝑛𝑎) =. 2 𝑚𝑤 𝑒 𝑀𝑟. =. 2∙50∙0,06 1000. = 6 ∙ 10−3 𝑚. (33). Oszacowana 𝐴𝑓(𝑠𝑧𝑎𝑐𝑜𝑤𝑎𝑛𝑎) różni się nieznacznie od wyników symulacji, gdyż nie uwzględnia mikrodrgań korpusu wynikających z obecności tłumienia uwzględnionego w układzie symulowanym. Bliskość rezultatów wskazuje na poprawność działania symulacji przenośnika. Na rysunkach 7, 8, 9, 10. przedstawiono przemieszczenia i prędkości na. poszczególnych współrzędnych uogólnionych.. 22.

(23) 0.2. v y [m/s]. 𝑣𝑦 [𝑚/𝑠]. 0.4. 0. x. 𝑣𝑥v [𝑚/𝑠] [m/s]. 0.1 0.05. -0.05 -0.1. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 0 -0.2 -0.4. 100. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. -3. 2. x 10. 0.01. y [𝑚]. y [m]. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. -0.01. 100. Rys.7. Przemieszczenie i prędkość pozioma i pionowa korpusu przenośnika w funkcji czasu. 0.1.  [rad/s]. 𝑟𝑎𝑑 ] 𝑠. 0.05. 𝜔𝛼 [. 0 -0.05 -0.1. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. -3. 2. x 10.  [rad]. 𝛼 [𝑟𝑎𝑑]. 1 0 -1 -2. 0. Rys.8. Kąt wahań i predkość kątowa korpusu przenośnika w funkcji czasu. 1. Vf [m/s]. 0.5. 𝑣𝑓 [𝑚/𝑠]. 0. 0 -0.5 -1. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 0.01 0.005. 𝑓 [𝑚]. -2. 0 -0.005. -1. f [m]. x [m]. x [𝑚]. 1. 0.005. 0 -0.005 -0.01. Rys.9. Przemieszczenie i prędkość względna rynny przenośnika w funkcji czasu.. 23.

(24) 100. 50. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 50. 0. 100. 15000. 15000. 10000. 10000.  2 [rad]. 𝜑2 [𝑟𝑎𝑑]. 𝜑1 [𝑟𝑎𝑑]. 0.  2 [rad]. 150. 𝑟𝑎𝑑 𝜔𝜑2 [ ] 𝑠  2 [rad/s].  1 [rad/s]. 𝑟𝑎𝑑 𝜔𝜑1 [ ] 𝑠. 150. 5000 0. 0 -5000. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. 5000. -5000. 0. 0. Rys.10. Prędkość kątowa i przemieszczenie kątowe wibratorów 1 i 2 w funkcji czasu.. Sytuacja przedstawiona na powyższych rysunkach jest efektem samosynchronizacji dwóch biegnących przeciwbieżnie wibratorów. Silniki włączone w tej samej chwili z różnicą fazową Δ𝜑 = 2,094 𝑟𝑎𝑑 (rys. 11) po niemal 1 sekundzie doznają synchronizacji dzięki momentom wibracyjnym działającym ze strony zamocowań sprowadzając różnicę fazową do wartości oscylującej nieznacznie wokół 0 a więc tak, aby wibratory tworzyły dwie siły, których wypadkowa pokrywa się z kierunkiem symetralnej łączącej środki mas korpusu i rynny. Dzieje się tak dzięki zjawisku samosynchronizacji, której dokładne omówienie przedstawione jest w kolejnych rozdziałach.. 2.5. 2. Δ𝜑 d[𝑟𝑎𝑑] [rad]. 1.5. 1. 0.5. 0. -0.5. 0. 10. 20. 30. 40. 50 Czas [s]. 60. 70. 80. 90. 100. Rys.11. Różnica przemieszczeń kątowych ∆𝜑 wirników 1 i 2 w funkcji czasu.. 24.

(25) 5.2. Symulacja przenośnika obciążonego nadawą W następnej kolejności przeprowadzono analizę działania przenośnika obciążonego nadawą. Przyjęto następujące parametry do opisu sił kontaktowych podczas zderzeń brył: Współczynnik tarcia 𝜇𝑛𝑛 wynikajacy z naturalnego kąta usypu, który przyjęto 𝜃 = 35° 𝜇𝑛𝑛 = 𝑡𝑔(𝜃) = 0,7 Wpółczynik tarcia między pierwszą warstwą nadawy a rynną 𝜇𝑝𝑛 = 0,4 Wartość współczynnika restytucji przyjęto na poziomie: 𝑅 = 0,01 – współczynnik restytucji zderzeń pomiędzy kolejnymi warstwami nadawy Współczynnik sprężystości strefy kontaktowej przyjęto na poziomie 𝑘𝑠𝑡𝑟𝑒𝑓𝑦 𝑘𝑜𝑛𝑡𝑎𝑘𝑡𝑜𝑤𝑒𝑗 = 1,5 ∙ 108 𝑁/𝑚 Podobnie jak w przypadku przenośnika nieobciążonego nadawą w fazie rozruchu wystąpił wzrost amplitudy drgań generowany przez rozpędzający się silnik, którego częstość przechodzi przez kolejne częstości rezonansowe. Na rysunkach 12, 13, 14, 15 przedstawiono przemieszczenia i prędkości po współrzędnych uogólnionych przenośnika obciążonego nadawą w ilości 𝑚𝑛 = 840 𝑘𝑔. Rynna w warunkach pracy ustalonej uzyskuje bezwzględną amplitudę drgań rzędu: 𝐴𝑟(𝑠𝑦𝑚𝑢𝑙𝑎𝑐𝑗𝑎) = |𝑓(t) + 𝑥(t)𝑐𝑜𝑠𝛽 + 𝑦(t)sin𝛽|𝑚𝑎𝑥 = 5,91 ∙ 10−3 𝑚.. Wartość ta jest wciąż bardzo zbliżona do amplitudy drgań rynny oszacowanej z zasady zachowania środka masy. Zwiększenie nieznaczne amplitudy w przypadku przenośnika obciążonego wiąże się z oddziaływaniem uderzającej nadawy w rynnę przenośnika. Amplituda przemieszczenia bezwzględnego rynny obliczona na drodze symulacji cyfrowej w przypadku przenośnika zarówno nieobciążonego jak i obciążonego nadawą jest mniejsza od szacowanej ze względu na tłumienie podparcia maszyny, które powoduje, że masa korpusu doznaje niezerowych przemieszczeń w obszarze antyrezonansowym. 25.

(26) 0.4. 0.1. 0.2. 𝑣𝑦 v[𝑚/𝑠] [m/s]. 0. 0. y. v x [m/s]. 𝑣𝑥 [𝑚/𝑠]. 0.2. -0.1 -0.2. -0.2 -0.4. 0. 5. 10. 15 t [s]. 20. 25. 30. 0. 5. 10. 15 t [s]. 20. 25. 30. 5. 10. 15 t [s]. 20. 25. 30. -3. 0.01. 5. x 10. 0 -0.005 -0.01. 0. y [m]. y [𝑚]. x [m]. x [𝑚]. 0.005. 0. 5. 10. 15 t [s]. 20. 25. -5. -10. 30. 0. Rys. 12. Przemieszczenie i prędkość pozioma i pionowa korpusu przenośnika obciążonego nadawą (𝑚𝑛 = 840 𝑘𝑔) w funkcji czasu. 0.1. 𝜔𝛼 [𝑟𝑎𝑑/𝑠]  [rad/s]. 0.05 0 -0.05 -0.1. 0. 5. 10. 15 t [s]. 20. 25. 30. 5. 10. 15 t [s]. 20. 25. 30. -3.  [rad]. α [𝑟𝑎𝑑]. 2. x 10. 0. -2. -4. 0. Rys. 13. Przemieszczenie kątowe i prędkość kątowa korpusu przenośnika obciążonego nadawą (𝑚𝑛 = 840 𝑘𝑔) w funkcji czasu. 1. v f [m/s]. 𝑣𝑓 [𝑚/𝑠]. 0.5 0 -0.5 -1. 0. 5. 10. 15 t [s]. 20. 25. 30. 0. 5. 10. 15 t [s]. 20. 25. 30. 0.01. 𝑓f [𝑚] [m]. 0.005 0 -0.005 -0.01. Rys. 14. Przemieszczenie i prędkość względna rynny przenośnika obciążonego nadawą (𝑚𝑛 = 840 𝑘𝑔) w funkcji czasu.. 26.

(27) 150. 50. 100. 50. 0 0. 5. 10. 15 t [s]. 20. 25. 30. 4000. 4000. 3000. 3000.  2 [rad]. 𝜑2 [𝑟𝑎𝑑].  1 [rad]. 2 [rad/s]. 100. 0. 𝜑1 [𝑟𝑎𝑑]. 𝜔𝜑2 [𝑟𝑎𝑑/𝑠]. 1 [rad/s]. 𝜔𝜑1 [𝑟𝑎𝑑/𝑠]. 150. 2000 1000 0. 0. 5. 10. 15 t [s]. 20. 25. 0. 5. 10. 15 t [s]. 20. 25. 30. 0. 5. 10. 15 t [s]. 20. 25. 30. 2000 1000 0. 30. Rys. 15. Prędkość i przemieszczenie kątowe wibratora 1 i 2 przenośnika obciążonego nadawą (𝑚𝑛 = 840 𝑘𝑔) w funkcji czasu.. Silniki włączone w tej samej chwili z różnicą fazową ∆𝜑 = 2,094 𝑟𝑎𝑑 w krótkim czasie doznają synchronizacji jednak w odróżnieniu od przenośnika nieobciążonego nadawą różnica kątów fazowych ustaliła się na poziomie ∆𝜑 ≈ −12,55 𝑟𝑎𝑑 a więc w przybliżeniu jest to podwojona wartość. miary jednego obrotu wyrażonego w radianach (rys. 16). Wartość ta. odpowiadała ponownie sytuacji w której wirujące wibratory tworzyły dwie siły, których wypadkowa pokrywa się z kierunkiem symetralnej łączącej środki mas korpusu i rynny. Zatem wibratory osiągają bieg synchroniczny podczas pracy przenośnika z obciążeniem 𝑚𝑛 = 840 𝑘𝑔, jednak nie idą podczas rozruchu cały czas współbieżnie.. 0. -2. [𝑟𝑎𝑑]  [rad] dΔ𝜑. -4. -6. -8. -10. -12. -14 0. 5. 10. 15 Czas [s]. 20. 25. 30. Rys. 16. Różnica przemieszczeń kątowych ∆𝜑 wirników 1 i 2 przenośnika obciążonego nadawą (𝑚𝑛 = 840 𝑘𝑔) w funkcji czasu.. 27.

(28) Ruch nadawy pokazano na rysunku 17.b. Parametry przenośnika przyjęto w taki sposób, aby zapewnić jednotaktowy ruch nadawy, co obrazuje rysunek. Na jeden okres ruchu rynny przenośnika przypada jeden okres ruchu nadawy. Warunkiem wystąpienia jednotaktowości ruchu jest spełnienie warunku aby współczynnik podrzutu materiału transportowanego oscylował wokół wartości 𝑘𝑝 = 3,3. 𝐴𝑟(𝑠𝑦𝑚𝑢𝑙𝑎𝑐𝑗𝑎) ∙ 𝜔𝑤𝑦𝑚 2 ∙ 𝑠𝑖𝑛𝛽 𝑘𝑝 = 𝑔 ∙ 𝑐𝑜𝑠𝛿. (34). gdzie 𝐴𝑓(𝑠𝑦𝑚𝑢𝑙𝑎𝑐𝑗𝑎) = 0.00591 m – amplituda przemieszczenia rynny, 𝜔𝑤𝑦𝑚 = 104,7 𝑟𝑎𝑑/𝑠 – częstość wymuszająca w stanie ustalonym, 𝛽 = 30° – kąt pochylenia kierunku siły wymuszającej wibratorów do kierunku stycznej do powierzchni rynny przenośnika, 𝛿 = 0° – kąt pochylenia rynny. Podstawiając powyższe dane otrzymano wartość współczynnika 𝑘𝑝 bliską wartości pożądanej.. 𝑘𝑝 =. 0.00591∙104,72 ∙𝑠𝑖𝑛30° 9,81∙cos 0°. ≈ 3.3. Jednotaktowość przemieszczania nadawy (jak wykazano w pracy [31]) sprawia, że moc rozpraszana w procesie transportowania jest najmniejsza.. 28.

(29) -3. 4. x 10. przemieszczenie rynny warstwa1. 2. warstwa2 0. warstwa3 warstwa4. 𝑦𝑖 [𝑚]. y i [rad]. -2 -4 -6 -8 -10 -12. 0. 5. 10. 15 Czas [s]. 20. 25. 30. -3. x 10 -2. przemieszczenie rynny warstwa1. -3. warstwa2 warstwa3. y i [rad]. -4. warstwa4. -5. -6. -7. -8 21.1. 21.12 21.14 21.16 21.18 Czas [s]. 21.2. 21.22 21.24 21.26. a). b). Rys. 17. Przemieszczenie rynny przenośnika oraz kolejnych 4 warstw nadawy z kolumny środkowej.. 5.2.1. Prędkość nadawy w funkcji ilości nadawy na przenośniku. Jednym z głównych parametrów opisujących właściwości każdego przenośnika wibracyjnego jest szybkość transportowania, gdyż to ona determinuje projektantów tych urządzeń w pierwszej kolejności. W przypadku analizowanej maszyny prędkość ta jest uśrednioną wartością prędkości transportowania wszystkich 12 warstw. Na rys. 18 przedstawiono wykres prędkości transportowania dla różnej ilości nadawy transportowanej na rynnie przenośnika.. 29.

(30) 0.5. 𝑣𝑡𝑟 [m/s] vtr[𝑚/𝑠]. 0.4. 0.3. 0.2. 0.1. 0. 0. 500. 1000. 1500 mn [kg]. 2000. 2500. 3000. 𝑚𝑛 [𝑘𝑔]. Rys. 18. Zależność średniej szybkości transportowania 𝑣𝑡𝑟 (uśrednienie po warstwach) od obciążenia przenośnika nadawą 𝑚𝑛 .. Wyniki te potwierdzają rezultaty uzyskane w pracach: [10,6], gdzie wykazano, że pomimo obciążenia przenośnika znaczną ilością nadawy przy odpowiednio dobranym współczynniku podrzutu o wartości ok. 3, przenośnik zachowuje w dalszym ciągu wysokie prędkości transportowania.. 5.2.2. Wydajność w funkcji ilości nadawy. Funkcja pozwalająca określić wydajność przenośnika w zależności od masy nadawy znajdującej się na przenośniku i średniej prędkości transportowania ma postać:. 𝑄=. 𝑚𝑛 𝑘𝑔 ∙ 𝑣ś𝑟 [ ] 𝑙 𝑠. (35). gdzie 𝑚𝑛. – masa nadawy na rynnie przenośnika,. 𝑄. – wydajność. 𝑙 = 3𝑑 – długość rynny 𝑣ś𝑟. – prędkość średnia transportowania. 30.

(31) 140. 120. 𝑄 [𝑘𝑔/𝑠] Q [kg/s]. 100. 80. 60. 40. 20. 0. 0. 500. 1000. 1500 mn[𝑘𝑔] [kg] 𝑚 𝑛. 2000. 2500. 3000. Rys.19. Wykres wydajności przenośnika 𝑄 w funkcji masy znajdującej się na nim nadawy 𝑚𝑛 .. Jak wynika z wykresu w badanym zakresie wraz ze wzrostem ilości nadawy na przenośniku maszyna uzyskuje coraz większe wartości wydajności. Oczywiście jest to okupione zużyciem większej ilości energii co zostanie pokazane na kolejnych wykresach.. 5.2.3. Moc w funkcji wydajności. Energia konsumowana podczas procesu transportowania wiąże się z podtrzymaniem ruchu nadawy i przenośnika oraz nadawaniu prędkości transportowania kolejnym porcjom nadawy lądującym na przenośniku. W czasie całej symulacji wartość obciążenia a więc ilość nadawy na całej długości przenośnika jest taka sama. Zakłada się więc dostarczanie w sposób ciągły tych samych porcji nadawy, które opuszczają rynnę przenośnika. Moc zużywana podczas procesu liczona jest na podstawie mocy chwilowej silników, które są źródłem energii dostarczanej maszynie. W obliczeniach ujęto tylko moc w stanie ustalonym. Moc chwilową silników 𝑃𝑐ℎ wyraża wzór: 𝑃𝑐ℎ = 𝑀𝑒𝑙1 ∙ 𝜔1. (36). gdzie 𝑀𝑒𝑙 𝑖 – moment i-tego silnika elektrycznego, 𝜔𝑖 – prędkość i-tego silnika. 31.

(32) Na wykresie (rys. 20) przedstawiono moc chwilową silnika 1 przy obciążeniu przenośnika nadawą w ilości 𝑚𝑛 = 3000 𝑘𝑔. 9000 8000 7000. P1[W]. 𝑃1 [𝑊]. 6000 5000 4000 3000 2000 1000 0. 0. 5. 10. 15 t [s]. 20. 25. 30. t [𝑠]. Rys. 20. Moc chwilowa silnika 1 obciążonego nadawą w ilości mn = 3000 kg.. Zależność pozwalającą wyliczyć moc konsumowaną przez cały proces 𝑃𝑝𝑟𝑜𝑐 w czasie Δ𝑡 przedstawia wzór wiążący pracę dwóch silników: Δ𝑡. 𝑃𝑝𝑟𝑜𝑐. Δ𝑡. Δ𝑡. Δ𝑡. ∫ 𝑃𝑐ℎ1 ∙ 𝑑𝑡 + ∫0 𝑃𝑐ℎ2 ∙ 𝑑𝑡 ∫0 𝑀𝑒𝑙1 ∙ 𝜔1 ∙ 𝑑𝑡 + ∫0 𝑀𝑒𝑙2 ∙ 𝜔2 ∙ 𝑑𝑡 [𝑊] = 0 = Δ𝑡 Δ𝑡. (37). Z dokonanych obliczeń uzyskano zależność mocy 𝑃 zużywanej przez proces transportu w funkcji ilości nadawy na przenośniku (rys. 21) . Wykres ten wskazuje na wzrost zapotrzebowania mocy przy rosnącym obciążeniu nadawą rynny przenośnika. Bardziej istotny jednak z punktu widzenia eksploatatora jest punkt pracy maszyny, a więc przy jakiej wydajności transportowanie jest bardziej opłacalne.. 32.

(33) 7. 6. 5. P[kW]. P [𝑘𝑊]. 4. 3. 2. 1. 0 0. 500. 1000. 1500 mn [kg]. 2000. 2500. 3000. Rys. 21. Wykres zapotrzebowania mocy P w funkcji obciążenia rynny nadawą mn .. W tym celu dokonano operacji dzielenia mocy 𝑃 zużywanej na proces przez wydajność przenośnika 𝑄 w tym procesie a wynik przedstawiono w funkcji wydajności (rys. 22). Wynika z niego, że wraz ze wzrostem wydajności maleje zużycie energii na transport jednostkowej ilości nadawy, przy czym dla bardzo małych wydajności zjawisko to jest bardzo silne, zaś przy wzroście wydajności zmierza asymptotycznie do wartości stałej. 1.4. 1.2. 0.6. P/Q [kW*s/kg]. 𝑘𝑊 ∙ 𝑠 ] 𝑘𝑔. 0.8. 𝑃/𝑄 [. 1. 0.4. 0.2. 0. 0. 20. 40. 60 80 Q [kg/s]. 100. 120. 140. 𝑘𝑔 𝑄 [ ] 𝑠. Rys.22. Stosunek zapotrzebowania mocy do wydajności w funkcji wydajności.. Transport staje się najbardziej opłacalny, jeśli wartość wydajności przekracza 𝑄 = 25 𝑘𝑔/𝑠. Nie określono w pracy maksymalnej wartości wydajności. Ma ona jednak ograniczenia ze względu na obciążenie dynamiczne podłoża i parametry wytrzymałościowe rynny i całego przenośnika i wiąże się z odpowiednim doborem mocy napędu. 33.

(34) 5.2.4. Siły przekazywane na podłoże w funkcji wydajności. Przenośnik wibracyjny wykorzystujący eliminacje dynamiczną posiada tę zaletę, że nie generuje znaczących sił na podłoże. W przypadku nieobciążonego przenośnika dokonuje się redukcja sił działających na korpus, a więc siły wymuszającej przenoszonej z dwóch wibratorów oraz sił generowanych przez resory, na których zawieszona jest rynna na korpusie. Drgania korpusu wówczas praktycznie ustają. Pozostają jednak mikrodrgania wynikające z tłumienia występującego w układzie. Podczas procesu transportowania nadawy, a więc przy dodatkowym obciążeniu maszyny dochodzi do częstych zderzeń nadawy z rynną przenośnika co prowadzi do pojawienia się dodatkowych sił działających na rynnę. Generowana zmienna, inna niż poprzednio siła w resorach nie jest w stanie zrównoważyć sił działających od zamocowań wibratorów co prowadzi do wzrostu drgań korpusu i sił przekazywanych na podłoże. Ponadto, zwiększanie ilości nadawy na rynnie powoduje pogorszenie synchronizacji wibratorów – patrz rozdz. 8.3.4., co wpływa na zwiększenie drgań korpusu i sił przenoszonych na podłoże. Na rysunku 23. przedstawiono wykres siły przekazywanej na podłoże od obciążenia rynny nadawą. 4. 4. x 10. 3. 2.5. 2.5. 2. 2. Rmax [N]. 𝑅𝑚𝑎𝑥 [𝑁]. Rmax [N]. 𝑅𝑚𝑎𝑥 [𝑁]. 3. 1.5. 1.5. 1. 1. 0.5. 0.5. 0 0. 500. 1000. 1500 mn [kg]. 2000. 2500. x 10. 0 0. 3000. 20. 40. 60 80 Q [kg/s]. 𝑚𝑛 [𝑘𝑔]. 𝑄 [𝑘𝑔/𝑠]. a). b). 100. 120. 140. Rys. 23. Wykres siły przekazywanej na podłoże w zależności: a) od ilości nadawy na rynnie b) od wydajności. Poniżej zaprezentowano sposób obliczenia średniej wartości obciążenia podłoża dla stanu pracy ustalonej: 𝑅𝑚𝑎𝑥 = √(Σ𝑅𝑥 )2 + (Σ𝑅𝑦 ). 2. [𝑁]. (38). gdzie Σ𝑅𝑥 – sumaryczna siła dynamiczna działająca na sprężyny podporowe na kierunku x, Σ𝑅𝑦 – sumaryczna siła dynamiczna działająca na sprężyny podporowe na kierunku y.. 34.

(35) 6. ANALIZA MODALNA MODELOWANEGO PRZENOŚNIKA Z punktu widzenia zastosowań technicznych ważne jest również zachowanie przenośnika podczas rezonansu przejściowego przy rozruchu i wybiegu. Narastające wówczas wielokrotnie amplitudy drgań mogą prowadzić do uszkodzenia maszyny lub nadmiernego obciążenia konstrukcji podporowej. Ponieważ rozruch odbywa się stosunkowo szybko, większe zagrożenie związane jest na ogół z wybiegiem swobodnym wibratora po wyłączeniu silnika. Aby prześledzić to zagadnienie dla przypadku przenośnika antyrezonansowego należy wyznaczyć częstości i formy jego drgań własnych.. 6.1 Linearyzacja układu nieliniowego Z uwagi na niewielkie rozproszenie energii w sprężynach , pochodzące głównie od tłumienia materiałowego i konstrukcyjnego, możliwe jest analizowanie drgań własnych układu jako nietłumionych. Pozostałe nieliniowości modelu maszyny bez nadawy związane są z wpływem drgań korpusu na ruch wibratorów oraz z przyspieszeniem Coriolisa w ruchu złożonym rynny. Stosunki wymiarowo-masowe w typowych konstrukcjach maszyn wibracyjnych pozwalają na ogół na pominięcie w analizie drgań własnych ruchu wibratora i skupienie jego masy w osi obrotu [11]. Podobnie analizując wartość przyspieszenia Coriolisa, można dla typowych stosunków wymiarowych maszyny, pominąć je w porównaniu do pozostałych składników przyspieszenia. Postępując w ten sposób można równania dynamiczne (21) do (24) sprowadzić do postaci: [𝑀] ∙ {𝑞̈ } + [𝐾] ∙ {𝑞} = {0}. (39). gdzie: {0} – wektor zerowy, [𝑀] – macierz mas, [𝐾] – macierz sprężystości Rozwiązanie ogólne równania jednorodnego przyjmuje postać:. {𝑞(𝑡)} = {𝑞0 } ∙ sin(𝜔𝑡 + 𝛾). (40). gdzie: {𝑞0 } – wektor amplitud współrzędnych : x, y, α, f, 𝜔 – częstość własna, 35.

(36) 𝛾 – kąt fazowy drgań. Stąd otrzymujemy równanie macierzowe: ([𝐾] − 𝜔2 ∙ [𝑀]) ∙ {𝑞(𝑡)} = {0}. (41). Istnienie rozwiązania niezerowego tego równania jest możliwe jeśli macierz ([𝐾] − 𝜔2 ∙ [𝑀]) jest osobliwa, tzn. 𝑑𝑒𝑡([𝐾] − 𝜔2 ∙ [𝑀]) = 0. (42). Powyższa zależność stanowi równanie 4 stopnia na ω2 i prowadzi do wyznaczenia 4 częstości własnych (niekoniecznie różnych). Częstości –ω odrzucamy jako pozbawione sensu fizycznego. Po uwzględnieniu postaci macierzy mas i sprężystości lewa strona tego równania przyjmuje postać:. k x  ( M k  M r ) 2  0 K    2  M    k x H  hr M r  2  2   M r cos( ). . . k x H  hr M r  2. 0 k y  ( M k  M r )  Lr M r  2.  Lr M r . 2. 2. k x H 2  k y L2  [ J k  J r  M r L2r  M r hr2 ] 2.  M r sin( ). 2. [ M r Lr sin( )  M r hr cos( )] 2.  M r cos( ) 2.    [ M r Lr sin( )  M r hr cos( )] 2    M r 2  k s   M r sin( ) 2. ( 43) 6.2 Częstości drgań własnych przenośnika Po podstawieniu przyjętych wcześniej wartości liczbowych otrzymujemy rozwiązanie w postaci 4 częstości własnych: ω1 = 16,42 𝑟𝑎𝑑/𝑠 ω2 = 21,82 𝑟𝑎𝑑/𝑠 ω3 = 27,67 𝑟𝑎𝑑/𝑠 ω4 = 124,57 𝑟𝑎𝑑/𝑠. 36.

(37) 6.3 Formy drgań własnych Ponieważ zerowanie się wyznacznika głównego macierzy oznacza, że równania są zależne liniowo, dlatego nie jest możliwe otrzymanie konkretnych wartości amplitud. Podstawienie do macierzy danej częstości własnej umożliwi wyznaczenie odpowiadającej jej formy drgań a więc stosunku amplitud poszczególnych współrzędnych. Dla konkretnych wartości częstości drgań otrzymano formy drgań w postaci: 𝐴 {𝑞 0 } = ( 𝐵 ) 𝐶 𝐷 Obliczenie. (44). pierwszej. formy. drgań dla pierwszej ω1 = 16,42 𝑟𝑎𝑑/𝑠.  2,197 105  0 2 K     M   q0    8,578 105  5  2,336 10. . . częstości. 0. 8,578 105. 1,384 106.  5,18 105.  5,18 105. 3,614 106.  1,349 105. 0. drgań. własnych.  2,336 105   A   0        1,349 105   B   0      C  0 0      1,069 107   D   0  (45). Korzystając z rozkładu macierzy LU otrzymano formę drgań odniesioną do amplitudy względnego przemieszczenia f : Dla ω1 = 16,42 𝑟𝑎𝑑/𝑠, 48,32 𝐴 −4,43 𝐵 (𝐶) = ( )∙𝐷 −12,10 𝐷 1. ( 45a). Postępując podobnie dla pozostałych częstości otrzymano: dla ω2 = 21,82 𝑟𝑎𝑑/𝑠, 18,22 𝐴 12,51 𝐵 (𝐶) = ( )∙𝐷 8,80 𝐷 1. (45b). 37.

(38) dla ω3 = 27,67 𝑟𝑎𝑑/𝑠, −10,75 𝐴 45,26 𝐵 (𝐶) = ( )∙𝐷 −11,10 𝐷 1. (45c). dla ω4 = 124,57 𝑟𝑎𝑑/𝑠, −0,25 𝐴 −0,15 𝐵 (𝐶) = ( )∙𝐷 −0,00014 𝐷 1. (45d). Dokonano sprawdzenia prawidłowości form drgań. Po podstawieniu warunków początkowych zgodnych z poszczególnymi formami drgań dla współrzędnych w układzie pozbawionym wymuszenia inercyjnego oraz tłumienia otrzymano drgania o częstości zgodnej z częstościami własnymi. 4 3 2. y [m]. 1 0 -1 -2 -3 -4 0. 0.5. 1 t [s]. 1.5. 2. Rys. 24. Przemieszczenie pionowe korpusu przy pierwszej formie drgań wykazało częstość ω1 = 16,42 𝑟𝑎𝑑/𝑠. 10. y [m]. 5. 0. -5. -10. 0. 0.5. 1. 1.5. 2. t [s]. Rys. 25. Przemieszczenie pionowe korpusu przy drugiej formie drgań wykazało częstość ω2 = 21,82 𝑟𝑎𝑑/𝑠. 38.

(39) 40 30 20. y [m]. 10 0 -10 -20 -30 -40 0.5. 1. 1.5 t [s]. 2. 2.5. Rys. 26. Przemieszczenie pionowe korpusu przy trzeciej formie drgań wykazało częstość ω3 = 27,67 𝑟𝑎𝑑/𝑠. 0.15. 0.1. y [m]. 0.05. 0. -0.05. -0.1. -0.15 0. 0.1. 0.2. 0.3. 0.4. 0.5 t [s]. 0.6. 0.7. 0.8. 0.9. 1. Rys. 27. Przemieszczenie pionowe korpusu przy czwartej formie drgań wykazało częstość ω4 = 124,57 𝑟𝑎𝑑/𝑠. 6.4 Wnioski Można zauważyć, że trzy pierwsze spośród czterech częstości drgań własnych położonych jest blisko siebie. Kotlina antyrezonansowa stanowi obszar miedzy częstościami ω3 = 27,67 𝑟𝑎𝑑/𝑠 oraz ω4 = 124,57 𝑟𝑎𝑑/𝑠 a wartość antyrezonansowa częstości wymuszającej ω𝑤𝑦𝑚 = 104,7 𝑟𝑎𝑑/𝑠.. 39.

(40) 7. WYZNACZENIE AMPLITUD DRGAŃ MAKSYMALNYCH PODCZAS WYBIEGU. 7.1. Klasyczne metody wyznaczania amplitud w rezonansie przejściowym Zjawisko rezonansu przejściowego występuje w miękko podpartych układach mechanicznych poddanych quasi-harmonicznej sile wymuszającej z częstotliwością zmieniającą się w sposób ciągły i przechodzącą przez częstości drgań własnych. Urządzenia, których dotyczy ten problem to m.in.. systemy wibroizolacji, niewyważone maszyny. wirnikowe strojone nadrezonansowo tj. przenośniki wibracyjne, przesiewacze, stoły wibracyjne, kraty wstrząsowe, młyny. Rezonans przejściowy pojawia się głównie podczas rozruchu i wybiegu maszyny. Ze względu na praktyczne znaczenie problemu, rezonans przejściowy był przedmiotem wielu badań teoretycznych. Rezonans ten jest bardzo niepożądany i niebezpieczny, gdyż prowadzi do dużych obciążeń dynamicznych elementów maszyn i konstrukcji [34]. W celu przedstawienia problemu narastających drgań w układach mechanicznych posłużono się przykładem maszyny wibracyjnej lub układu wibroizolacji niewyważonej maszyny wirnikowej, której schemat przedstawiony jest na rysunku 28.. Rys. 28. Model maszyny wibracyjnej: m – masa korpusu maszyny, m1 – masa niewyważona, e – mimośród masy niewyważonej, J – moment bezwładności zredukowany do punktu obrotu, Mel – moment napędowy (moment elektryczny), y – przemieszczenie korpusu, φ – kąt obrotu masy niewyważonej. 40.

(41) Załóżmy, że układ w stanie ustalonym pracuje w zakresie nadrezonansowym, tzn.. ω𝑢𝑠𝑡 ≫ ω0 = √. 𝑘 𝑚 + 𝑚1. gdzie ω0 – oznacza częstość własną układu nietłumionego. Układ taki podczas rozruchu i wybiegu wirnika przechodzi przez częstość rezonansową. Generuje to wzmożone drgania. Rezonans przejściowy jest zjawiskiem niepożądanym i niebezpiecznym, gdyż prowadzi do dużych obciążeń dynamicznych punktów podparcia, więzów kinematycznych i elementów konstrukcyjnych maszyn. Pierwsze analizy opisujące jakościowo rezonans przejściowy opierały się na przyjęciu a priori przebiegu siły wymuszającej jako zadanej funkcji czasu. Przyjmowano, że jej częstość jest funkcją liniowo zależną jedynie od czasu, zakładając jednocześnie stałą amplitudę tej siły lub uzależniając ją jedynie od prędkości i przyspieszenia kątowego wibratora bezwładnościowego. Równanie opisujące ruch korpusu maszyny z rys. 28. przedstawia się w tym drugim przypadku w następujący sposób: (𝑚1 + 𝑚)𝑦̈ + 𝑏𝑦̇ + 𝑘𝑦 = −𝑚𝑒𝜑̈ cos(𝜑) + 𝑚𝑒𝜑̇ 2 sin(𝜑) = 𝑓(𝑡). (46). gdzie 𝜑 = 𝜑(𝑡) − zakładana postać, zwykle liniowa, zmienności prędkości kątowej wirnika. Pierwsza praca poświęcona rezonansowi przejściowemu została opublikowana w 1932 roku przez F. Lewisa. Równanie powyższe dla stałej amplitudy siły wymuszającej zostało przekształcone do postaci. 𝑦̈ + 2𝛽𝜔𝑜 ∙ 𝑦̇ + 𝜔𝑜 2 ∙ 𝑦 =. gdzie. 𝜔𝑜 = √𝑚. 𝑘. 1 +𝑚. 𝑓(𝑡) 𝑚1 + 𝑚. (47). ,. 41.

(42) 𝛾 = 2(𝑚. 𝑏. 1 +𝑚)𝜔𝑜. Rozwiązanie powyższego równania przy zerowych warunkach początkowych zostało przedstawione w postaci 𝑦=. 𝑡 𝑡 1 [𝑒 𝑎𝑡 ∫ 𝑓(𝜏)𝑒 −𝑎𝑡 𝑑𝜏 − 𝑒 𝑏𝑡 ∫ 𝑓(𝜏)𝑒 −𝑏𝑡 𝑑𝜏] 2𝑖𝑣(𝑚1 + 𝑚) 0 0. (48). gdzie 𝑎 = −𝛾𝜔0 + 𝑖𝑣 𝑏 = −𝛾𝜔0 − 𝑖𝑣 𝑣 = 𝜔0 √1 − 𝛾 2 𝑖 = √−1 Kilkanaście lat później A. Kac [20] wyznaczył rozwiązanie równania dla siły wymuszającej 𝑓(𝑡) w postaci 𝜀𝑡 2 𝑓(𝑡) = 𝑚1 𝑒(𝜀𝑡) cos ( − 𝜑0 ) 2 2. (49). gdzie 𝜑0 − faza poczatkowa siły. Kac uzyskał w ten sposób równanie, którego rozwiązanie jest utrudnione ze względu na pracochłonne obliczanie całek składowych. W praktyce zadowalamy się wyznaczeniem przebiegu obwiedni drgań, której dokonał w swojej pracy Goliński [18]. Na podstawie pracy Lewisa [21], w której siła wymuszająca rozpatrywana jest jako siła sinusoidalna o stałym module ale liniowo narastającej prędkości, A. Kac wyznaczył zależność na wartość częstości rezonansowej (rys. 29). 𝜔𝑟𝑒𝑧 = 𝜔𝑜 (1 + Δ𝑛 ) =. 𝑡𝑔𝛿 ∙ √1,5 ∙ 𝜋 ∙ 𝑡𝑔𝛿 (√𝑡𝑔𝛿 + 0,28𝛾). (50). 2. 42.

(43) gdzie 𝑡𝑔𝛿 =. 𝜀 𝜔𝑜 2. – wskaźnik przyspieszenia. Rys. 29. Porównanie obwiedni amplitudy drgań wymuszonych w warunkach quasi-stacjonarnych (1) oraz w warunkach rezonansu przejściowego (2) odniesionych do wartości 𝐴𝑢 =. 𝑚1 𝑒 𝑚1 +𝑚. , tj. do wartości asymptoty. amplitudy drgań przy 𝜔 → ∞, gdzie A – amplituda drgań korpusu, 𝑚1 – masa niewyważenia [18]. Na bazie prac Lewisa [21] i Katza [20] opracowano nomogram (rys. 31) przedstawiający współczynnik zwielokrotnienia amplitudy 𝛼 dla masy m zależny od 𝑡𝑔𝛿 [19]:. 𝑡𝑔𝛿 = 𝜀 𝛾=. 𝑚 𝑘. 𝑏 2𝑚𝜔𝑜. 𝛼=𝐴. 𝜀. – współczynnik przyspieszenia. = 𝜔2. 𝑜. =. 𝐴𝑚𝑎𝑥. 𝑏 𝑘 𝑚. 2𝑚√. =. 𝑏 2√𝑘𝑚. – współczynnik tłumienia – współczynnik zwielokrotnienia amplitudy. 𝜔=∞,𝜉=0. 43.

(44) Nomogram jest rysunkiem przedstawiającym szereg wykresów będących funkcją kilku wybranych parametrów zestawionych w jednym układzie współrzędnych.. Rys. 31. Współczynnik zwielokrotnienia amplitudy drgań rezonansowych 𝛼 w funkcji wskaźnika przyspieszenia dla różnych wartości liczb tłumienia 𝛾. [19]. Metody przedstawione dotychczas cechuje ta wada, że uwzględniają one jedynie składową normalną rzeczywistej siły wymuszającej, która zawiera składową styczną. Ponadto przyspieszenie kątowe wibratora, które w tych metodach przyjmowane jest jako wartość stała a w rzeczywistości podlega znaczącym zmianom [29]. Późniejsze opracowania opierały się na podobnych założeniach co prace Lewisa i Kaca z wyjątkiem analizy Markerta i Seidlera (22).. Równanie pozwalające obliczyć. amplitudy w rezonansie przejściowym zawierało tym razem czynniki wiążące składową styczną siły wymuszającej. Rozwiązaniem równania przedstawionego w tej pracy jest formuła pozwalająca obliczyć maksymalną amplitudę w rezonansie, której główny czynnik przedstawiono poniżej:. 𝐴𝑚𝑎𝑥 ≈. 𝑐𝜔𝑛. (51). √𝜀. gdzie 𝜀 – przyspieszenie kątowe wibratora, 𝜔𝑛 – częstość rezonansowa 44.

(45) Głównym błędem popełnianym we wszystkich dotychczasowych pracach było przyjmowanie stałej wartości przyspieszenia kątowego, którego dokładne oszacowanie w obszarze okołorezonansowym jest bardzo trudne. Podobne prace prowadzili również Banaszewski i Turkiewicz [32, 33]. Efektem ich prac było uzyskanie formuły na amplitudę maksymalną w rezonansie przejściowym zależną od wartości spodziewanej przyspieszenia kątowego  niewyważonego wirnika stanowiącego źródło sił wymuszających. Dla układów słabo tłumionych zależność ta odpowiada na ogół w przybliżeniu formule:. 𝐴𝑚𝑎𝑥 =. 𝑎. ( 52 ). √𝜀. gdzie: a – wielkość stała.. 7.1.1. Przyczyny błędów metod klasycznych.. Rzeczywistym modelem odwzorowującym pracę maszyny z rysunku 28 jest układ o dwóch stopniach swobody a równania dynamiczne opisujące ten układ zaprezentowano poniżej:. (𝑚1 + 𝑚)𝑦̈ + 𝑏𝑦̇ + 𝑘𝑦 = −𝑚𝑒𝜑̈ cos(𝜑) + 𝑚𝑒𝜑̇ 2 sin(𝜑) 𝐽𝑧𝑟 𝜑̈ = 𝑀𝑧 − 𝑀𝑤. (53). gdzie 𝑀𝑧 → momenty: napędowy (podczas rozruchu) i moment oporu, 𝑀𝑤 = −𝑚𝑒𝑦̈ cos(𝜑) → moment wibracyjny Składnik 𝑀𝑤 nazwany momentem wibracyjnym [15] swoje źródło ma w oddziaływaniu sił reakcji w przegubie łączącym wibrator z korpusem. Podczas rezonansu jego wartość znacząco rośnie co ma decydujący wpływ na wartość przyspieszenia kątowego wibratora. Rzeczywista wartość przyspieszenia przy wybiegu maszyny jest zdecydowanie większa od wartości przyspieszenia wynikającego jedynie z oporów ruchu a dzieje się tak poprzez generowany hamujący moment wibracyjny. Oddziaływanie momentu wibracyjnego na prędkość kątowa wibratora przedstawiono na poniższym rysunku [16]: 45.

(46) a b. czas 𝑡 [𝑠][s] Rys. 30. Wpływ momentu wibracyjnego na prędkość kątową wibratora podczas wolnego wybiegu, 𝑀𝑧 < 0 a) wybieg bez momentu wibracyjnego, b) wybieg z uwzględnionym momentem wibracyjnym. Wartości amplitud uzyskiwanych z metod przedstawionych dotychczas w tym rozdziale cechowały się dużymi błędami prowadząc do kilkukrotnego przeszacowania amplitud podczas wybiegu [29] co stało się przyczyną prowadzenia dalszych analiz i tworzenia nowych sposobów szacowania amplitud w rezonansie przejściowym.. 7.1.2. Metoda nomograficzna (Cieplok 2009) [25, 26]. Przegląd prac uwzględniających opisane wyżej zjawisko zaczniemy od pozycji [26], w której rozpatrzono przebieg rezonansu przejściowego dla układu symetrycznego pokazanego na rys.31. We wspomnianej pracy wykorzystano transformację równań ruchu modelu maszyny symetrycznej do wirującego z prędkością masy niewyważonej układu współrzędnych. Stało się możliwe zdefiniowanie jednostek względnych i dwukrotne zmniejszenie liczby parametrów opisujących ruch układu.. 46.

(47) Rys. 31. Model maszyny wibracyjnej: mk – masa korpusu maszyny, m– masa niewyważona, e – mimośród masy niewyważonej, Jzr – moment bezwładności zredukowany do punktu obrotu, Mel – moment napędowy (moment elektryczny), x,y – przemieszczenie korpusu, φ – kąt obrotu masy niewyważonej. Zbiór sześciu parametrów fizycznych mk, me, Jzr, Mel, k, b potrzebnych do zapisania równań dynamicznych we współrzędnych naturalnych został zredukowany do trzech parametrów , q,  w zapisie przy zastosowaniu parametrów względnych:. m2e2 ( mk  m ) J zr.  q.  . (54). M el 1 J zr 02. (55). b 2 ( mk  m )k. (56). k mk  m. (57). gdzie:. o . Uzyskano przekształcone równania, które pozwoliły na czytelną graficzną reprezentację dużej liczby rozwiązań równań ruchu układu uzyskanych przy wykorzystaniu programu do obliczeń numerycznych. Na rysunku 32 przestawiono uzyskany w ten sposób nomogram dla fazy wybiegu swobodnego maszyny symetrycznej. Wykres ten może być również. 47.

(48) wykorzystany dla maszyn o trajektorii liniowej, co wymaga odczytu dla dwukrotnie pomniejszonej wartości  w stosunku do danej zależnością (54).. Rys. 32. Współczynnik zwielokrotnienia amplitudy drgań korpusu maszyny przy wybiegu swobodnym wibratora.. Dzięki metodzie zaproponowanej przez jej autora (Cieplok, 2009) dla maszyn o konstrukcji symetrycznej (rys. 31) oraz tych, które wykonują ruch prostoliniowy (rys. 28) otrzymujemy stosunkowo dokładne oszacowanie wzmocnień amplitudowych zarówno dla rozbiegu jak i wybiegu maszyny. Pozwala również na badanie wybiegów z uwzględnieniem oporów ruchu. Istnieją ograniczenia powyższego sposobu szacowania amplitud w rezonansie przejściowym. Pierwsza uwaga dotyczy symetryczności układu przedstawionego na rys. 31. Obiekty rzeczywiste o ruchu płaskim nie odpowiadają na ogół takiemu schematowi. Wynika to ze zróżnicowania współczynników sprężystości na kierunku x i y. Wykazano, że układy mechaniczne cechują się dużą wrażliwością na odmienne wartości tych stałych. Trudno zatem zastosować w tych przypadkach model maszyny symetrycznej. Drugim problemem są trudności zastosowania metody do układów, gdzie ruch maszyny odbywa się po trajektorii odcinkowej gdy dla osiągnięcia takiego ruchu wykorzystuje się dwa przeciwbieżne wibratory inercyjne a ich współfazowość uzyskuje się dzięki synchronizacji swobodnej [11]. Jak wykazano w pracy [27], bieg współfazowy przestaje być stabilny przy 48.

(49) wybiegu po przekroczeniu częstości drgań wahadłowych korpusu. Wibratory w obszarze tej częstości charakteryzuje tendencja do rozfazowania. Zatem schemat układu zmienia się zasadniczo w stosunku do przyjętego w omawianej metodzie .. 7.2. Metoda bilansu energetycznego I [13,29] Ponieważ omówione wcześniej metody ograniczone są do analizy układów o jednym stopniu swobody, z uwzględnieniem lub bez uwzględnienia stopnia swobody związanego z ruchem obrotowym wibratora oraz ze względu na występującą w analizowanym układzie utratę współfazowego biegu wibratorów w obszarze okołorezonansowym – rys.33 ,. Rys.33. Różnica przemieszczeń kątowych silników nr 1 i 2 przy wybiegu.. nie znajdują zastosowania inne metody oceny amplitud maksymalnych w rezonansie przejściowym [12], gdyż zakładają one współfazowy bieg wibratorów i nie pozwalają na wyznaczenie maksymalnych amplitud drgań kątowych maszyny. Dlatego do analizy amplitud maksymalnych zastosowano metodę energetyczną [ 13, 29]. Metoda ta oparta jest na złożeniu, że energia zespołu wibratorów , z jaką wchodzą one w rezonans przejściowy, przekształcana jest głównie na wzrost energii drgań układu. Zjawisko to. powoduje. gwałtowne. załamanie. prędkości. okołorezonansowej , pokazane na rys. 34.. 49. obrotowej. wibratorów. w. strefie.

(50) 100. 80. 1 [rad/s]. 60. 40. 20. 0. -20 90. 100. 110. 120. 130. 140. 150. 160. t [s]. Rys.34. Zjawisko utraty prędkości kątowej wibratora w zakresie okołorezonansowym. Warto nadmienić, że metoda oparta na bilansie energetycznym została pierwszy raz wykorzystana w pracy [Agranovskaja & Blechman, 1969] z tym, że energia kinetyczna zespołu wibratorów nie została przyrównana do energii kinetycznej pozostałej części maszyny (korpus i rynna) lecz do maksymalnej energii potencjalnej układu zawieszenia korpusu. Trudnością tej metody stała się konieczność określenia prędkości kątowej przy której następuje przemiana energii a prędkość ta nie odpowiada bezpośrednio częstościom drgań własnych układu (Lewis, 1932). Metoda ta również nie wskazuje sposobu jej określenia. Bilans energetyczny analizowanej maszyny a więc równość energii kinetycznej ruchu obrotowego wibratorów i energii kinetycznej układu drgającego przedstawia się zatem w następujący sposób [13]:. n. 1 1 T  J zr  02i  q max i Mq max i 2 2. (58). gdzie n – liczba wibratorów, Jzr – moment bezwładności układu napędowego zredukowany do wału wirnika, ω0i – prędkość kątowa, przy której następuje wymiana energii (w przybliżeniu i-ta częstość własna układu), qmax i – wektor współrzędnych opisujących drgania układu z i-tą formą drgań , M – macierz mas zlinearyzowanego układu drgającego, wyznaczona na podstawie (59).. 50.