VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ BRNO UNIVERSITY OF TECHNOLOGY

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ BRNO UNIVERSITY OF TECHNOLOGY

FAKULTA STAVEBNÍ ÚSTAV STAVEBNÍ MECHANIKY FACULTY OF CIVIL ENGINEERING INSTITUTE OF STRUCTURAL MECHANICS

MODELOVÁNÍ VOLNÉHO PRUTU ZATÍŽENÉHO SLEDUJÍCÍM ZATÍŽENÍM MODELLING OF FREE BEAM LOADED BY FOLLOWER LOAD

BAKALÁŘSKÁ PRÁCE BACHELOR'S THESIS

AUTOR PRÁCE

JAN MAŠEK

AUTHOR

VEDOUCÍ PRÁCE SUPERVISOR

BRNO 2014

Ing. PETR FRANTÍK, Ph.D.

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ FAKULTA STAVEBNÍ Studijní program Typ studijního programu Studijní obor Pracoviště

B3607 Stavební inženýrství Bakalářský studijní program s prezenční formou studia 3647R013 Konstrukce a dopravní stavby Ústav stavební mechaniky

ZADÁNÍ BAKALÁŘSKÉ PRÁCE Student

Jan Mašek

Název

Modelování volného prutu zatíženého sledujícím zatížením

Vedoucí bakalářské práce

Ing. Petr Frantík, Ph.D.

Datum zadání bakalářské práce Datum odevzdání bakalářské práce V Brně dne 30. 11. 2013

30. 11. 2013 30. 5. 2014

............................................. prof. Ing. Drahomír Novák, DrSc. Vedoucí ústavu

................................................... prof. Ing. Rostislav Drochytka, CSc., MBA Děkan Fakulty stavební VUT

Podklady a literatura Literatura dle vývoje a pokynů vedoucího práce. Brepta, R., Půst, L., Turek, F.: Mechanické kmitání, Technický průvodce 71, nakladatelství Sobotáles, Praha, 1994. Zásady pro vypracování Nastudování potřebných znalostí dle pokynů vedoucího práce. Zorientování se v problematice. Vytvoření numerických modelů a jejich aplikace. Předepsané přílohy

............................................. Ing. Petr Frantík, Ph.D. Vedoucí bakalářské práce

Abstrakt C´ılem pˇredkl´ adané pr´ ace je vytvoˇren´ı nelineárn´ıho dynamického modelu volného prutu zat´ıˇzeného sleduj´ıc´ı silou. Re´ aln´ ym pˇredobrazem modelu je ˇst´ıhlá pruˇzná raketa zat´ıˇzen´ a tahem reaktivn´ıho motoru. Vzhledem k povaze u ´lohy mus´ı model umoˇzn ˇovat velké deformace. Dalˇs´ımi poˇzadavky na aplikovan´ y model je vystiˇzen´ı materiálového tlumen´ı, tlumen´ı vlivem interakce modelu s okoln´ım prostˇred´ım a zaveden´ı vlivu gravitaˇcn´ıho pole. Model bude v budoucnu pouˇzit pro studium pokritického chován´ı uvaˇzované konstrukce, proto je jednou z priorit n´ızk´ a ˇcasov´ a n´ aroˇcnost simulace. Vlastn´ı vytvoˇren´ a formulace numerického modelu bude implementována v programovac´ım jazyku Java. Pro zobrazen´ı pr˚ ubˇehu simulace a sledován´ı stavov´ ych promˇenn´ ych bude vytvoˇreno odpov´ıdaj´ıc´ı grafické prostˇred´ı. Správnost odvozeného modelu bude ovˇeˇrena porovnán´ım vybran´ ych hodnot s analytick´ ym ˇreˇsen´ım.

Kl´ıˇ cov´ a slova Fyzik´ aln´ı diskretizace kontinua, metoda tuh´ ych d´ılc˚ u, Lagrangeovská mechanika, pohybové rovnice, nekonzervativn´ı zat´ıˇzen´ı, velké deformace.

4

Abstract The aim of the presented thesis is to create a non-linear dynamical model of a free rod loaded by a follower force. The model is inspired by a slender flexible missile loaded by an end thrust. Because of the nature of the problem, the model has to be capable of large deflections. Another requirements on the model are to implement material the damping and the damping due to the interaction of the model with a surrounding medium and the influence of gravitational field. In the future, the model will be used for examination of the post-critical behavior of such construction. Therefore low computational demands of the simulation are required. The derived formulation of the numerical model will be implemented using the Java programming language. For observation of the simulation process and for monitoring of the state variables, an appropriate graphic interface will be created. The accuracy of the derived model will be verified by the comparison of selected values to the analytical solution.

Keywords Physical discretization of continuum, rigid finite element method, Lagrangian mechanics, equations of motion, nonconservative load, large deflections.

5

Bibliografick´ a citace pr´ ace ˇ MASEK, Jan. Modelov´ an´ı volného prutu zat´ıˇzeného sleduj´ıc´ım zat´ıˇzen´ım. Brno, 2014. 64 s., ´ 15 s. pˇr´ıl. Bakal´ aˇrsk´ a pr´ ace. Vysoké uˇcen´ı technické v Brnˇe, Fakulta stavebn´ı, Ustav stavebn´ı mechaniky. Vedouc´ı pr´ ace Ing. Petr Frant´ık, Ph.D.

6

Prohl´ aˇ sen´ı Prohlaˇsuji, ˇze jsem svou bakal´ aˇrskou práci zpracoval samostatnˇe, pod odborn´ ym veden´ım vedouc´ıho pr´ ace Ing. Petra Frant´ıka, Ph.D., a ˇze jsem uvedl vˇsechny pouˇzité informaˇcn´ı zdroje.

V Brnˇe dne 30.5.2014

.......................................... podpis autora Jan Maˇsek

7

Podˇ ekov´ an´ı: Rád bych podˇekoval vedouc´ımu své bakaláˇrské práce panu Ing. Petru Frant´ıkovi, Ph.D. za otevˇren´ y a konstruktivn´ı odborn´ y a pˇrátelsk´ y osobn´ı pˇr´ıstup v pr˚ ubˇehu naˇs´ı spolupráce. Pˇredevˇs´ım d´ıky jeho neziˇstné pomoci jsem si osvojil schopnosti a zp˚ usob myˇslen´ı nezbytné k vypracov´ an´ı této pr´ ace. ´ Dále bych r´ ad podˇekoval vyuˇcuj´ıc´ım a vˇedeck´ ym pracovn´ık˚ um Ustavu stavebn´ı mechaniky, kteˇr´ı ve mnˇe sv´ ym nadˇsen´ım a peˇclivou prac´ı probudili podobné zal´ıben´ı v mechanice. V neposledn´ı ˇradˇe chci podˇekovat sv´ ym rodiˇc˚ um, kteˇr´ı mˇe vˇzdy ne´ unavnˇe podporovali a poskytnuli mi tak skvˇelé z´ azem´ı po celou dobu studia. Podˇekován´ı patˇr´ı i mému bratrovi, bez jehoˇz pomoci a podpory by pro mne studium bylo mnohem obt´ıˇznˇejˇs´ı.

V Brnˇe dne 30.5.2014

.......................................... podpis autora Jan Maˇsek

8

Obsah 1 Teoretick´ yu ´ vod 1.1 Klasick´ a mechanika . . . . . . . . . . . . . . . . . . . . . . ´ 1.1.1 Uvod a historie . . . . . . . . . . . . . . . . . . . . 1.1.2 Pravo´ uhl´ a souˇradná soustava . . . . . . . . . . . . 1.1.3 Pol´ arn´ı souˇradn´ a soustava . . . . . . . . . . . . . . 1.1.4 Newtonovy pohybové zákony . . . . . . . . . . . . 1.1.5 Z´ akon zachov´ an´ı hybnosti . . . . . . . . . . . . . . 1.1.6 Raketov´ y pohon . . . . . . . . . . . . . . . . . . . 1.2 Energie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Potenci´ aln´ı energie . . . . . . . . . . . . . . . . . . 1.2.2 Kinetick´ a energie . . . . . . . . . . . . . . . . . . . 1.2.3 Z´ akon zachov´ an´ı energie . . . . . . . . . . . . . . . 1.3 Lagrangeovsk´ a mechanika . . . . . . . . . . . . . . . . . . 1.3.1 Zobecnˇené souˇradnice systému . . . . . . . . . . . 1.3.2 Konfiguraˇcn´ı prostor, zobecnˇené rychlosti systému 1.3.3 Lagrangeova funkce . . . . . . . . . . . . . . . . . 1.3.4 Virtu´ aln´ı posunut´ı, virtuáln´ı práce . . . . . . . . . 1.3.5 D’Alembert˚ uv princip . . . . . . . . . . . . . . . . 1.3.6 Hamilton˚ uv pricip nejmenˇs´ı akce . . . . . . . . . . 1.3.7 Lagrangeovy rovnice . . . . . . . . . . . . . . . . . 1.4 Pohybové rovnice a jejich ˇreˇsen´ı . . . . . . . . . . . . . . . 1.4.1 Pohybové rovnice . . . . . . . . . . . . . . . . . . . 1.4.2 Eulerova metoda . . . . . . . . . . . . . . . . . . . 1.4.3 Semi-implicitn´ı Eulerova metoda . . . . . . . . . . 1.4.4 Runge-Kuttova metoda . . . . . . . . . . . . . . . 1.5 Zdroje nelinearit . . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Geometrick´ a nelinearita . . . . . . . . . . . . . . . 1.5.2 Materi´ alov´ a nelinearita . . . . . . . . . . . . . . . 1.5.3 Nelinearita zp˚ usobená okrajov´ ymi podm´ınkami . . 1.6 Konzervativn´ı a nekonzervativn´ı s´ıly . . . . . . . . . . . . 1.6.1 Konzervativn´ı s´ıly . . . . . . . . . . . . . . . . . . 1.6.2 Nekonzervativn´ı s´ıly . . . . . . . . . . . . . . . . . 1.7 Metoda tuh´ ych d´ılc˚ u . . . . . . . . . . . . . . . . . . . . .

9

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13 13 13 14 16 17 18 19 21 21 22 22 24 24 24 24 25 26 26 27 28 28 28 29 30 31 31 32 32 33 33 33 34

OBSAH

2 Odvozen´ı numerick´ eho modelu ´ 2.1 Uvod . . . . . . . . . . . . . . . . . . . . . 2.1.1 Princip odvozen´ı . . . . . . . . . . 2.2 Energie . . . . . . . . . . . . . . . . . . . 2.2.1 Kinetick´ a energie modelu . . . . . 2.2.2 Potenci´ aln´ı energie modelu . . . . 2.3 Pohybové rovnice . . . . . . . . . . . . . . 2.3.1 Maticov´ y tvar . . . . . . . . . . . . 2.4 Aplikovan´ y model volného prutu . . . . . 2.4.1 Sleduj´ıc´ı s´ıla . . . . . . . . . . . . 2.4.2 Materi´ alové tlumen´ı . . . . . . . . 2.4.3 Tlumen´ı vlivem okoln´ıho prostˇred´ı 2.4.4 Gravitaˇcn´ı zrychlen´ı . . . . . . . . 2.4.5 Aplikovan´ y model . . . . . . . . .

. . . . . . . . . . . . .

36 36 37 37 37 40 41 44 47 47 48 49 52 52

. . . . .

53 53 55 56 56 58

. . . . . .

59 59 60 60 61 65 67

5 Z´ avˇ er 5.1 V´ yhledy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70 71

Literatura

72

Seznam symbol˚ u

74

Seznam pˇ r´ıloh

76

Pˇ r´ılohy A: Implementace modelu v jazyce Java . . . . . . . . . . . . . . . . . . . . . . . . . B: Kompaktn´ı disk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77 77 90

3 V´ ysledky numerick´ ych simulac´ı 3.1 Vlastn´ı frekvence prutu . . . . 3.2 Konzola zat´ıˇzen´ a vlastn´ı t´ıhou 3.3 Kritick´ a s´ıla . . . . . . . . . . . 3.3.1 Voln´ y prut . . . . . . . 3.3.2 Konzola . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . .

. . . . .

4 Implementace modelu 4.1 Implementace metod ˇreˇsen´ı dynamického systému . 4.1.1 Eulerova metoda . . . . . . . . . . . . . . . 4.1.2 Semi-implicitn´ı Eulerova metoda . . . . . . 4.1.3 Runge-Kuttova metoda . . . . . . . . . . . 4.2 Srovn´ an´ı metod ˇreˇsen´ı dynamického systému . . . 4.3 Grafické prostˇred´ı modelu . . . . . . . . . . . . . .

10

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

. . . . . . . . . . . . .

. . . . .

. . . . . .

Seznam obr´ azk˚ u 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 1.10 1.11 1.12 1.13 1.14 1.15

Zatˇr´ıdˇen´ı klasické mechaniky. . . . . . . . . . . . . . . . . . . . . . Poloha bodu P v souˇradné soustavˇe. . . . . . . . . . . . . . . . . . Poloha bodu P vzhledem k pohybujic´ı se souˇradné soustavˇe x0 y 0 z 0 . Poloha bodu P v pol´ arn´ı souˇradné soustavˇe. . . . . . . . . . . . . . Raketa zat´ıˇzen´ a tahem motoru. . . . . . . . . . . . . . . . . . . . . Translaˇcn´ı a rotaˇcn´ı pruˇzina. . . . . . . . . . . . . . . . . . . . . . Tuhé kyvadlo s jedn´ım stupnˇem volnosti. . . . . . . . . . . . . . . Trajektorie v´ yvoje dynamického systému. . . . . . . . . . . . . . . Konzola zat´ıˇzen´ a osamˇelou silou. . . . . . . . . . . . . . . . . . . . Neline´ arn´ı pracovn´ı diagram materiálu. . . . . . . . . . . . . . . . Pˇr´ıklady nˇekter´ ych neline´ arn´ıch vazeb. . . . . . . . . . . . . . . . . Zn´ azornˇen´ı konzervativn´ıch vlastnost´ı gravitaˇcn´ı s´ıly. . . . . . . . . Pruˇzn´ a ˇst´ıhl´ a raketa zat´ıˇzená tahem motoru. . . . . . . . . . . . . Idealizace konstrukce metodou tuh´ ych d´ılc˚ u. . . . . . . . . . . . . . Varianty kontaktn´ıch kumulátor˚ u pˇretvoˇren´ı. . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

13 14 16 16 19 22 22 27 31 32 33 34 34 34 35

2.1 2.2 2.3 2.4 2.5

Model volného prutu zat´ıˇzeného sleduj´ıc´ı silou. . . . . Model volného prutu. . . . . . . . . . . . . . . . . . . Prvn´ı d´ılec modelu. . . . . . . . . . . . . . . . . . . . . P˚ usoben´ı sleduj´ıc´ı s´ıly. . . . . . . . . . . . . . . . . . . Odvozen´ı z´ avislosti u ´tlumu na smˇeru pohybu modelu.

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

36 37 38 47 50

3.1 3.2 3.3 3.4 3.5 3.6

Konvergence hodnoty prvn´ı vlastn´ı frekvence volného prutu. Z´ avislost v´ ypoˇcetn´ı doby na poˇctu d´ılc˚ u. . . . . . . . . . . . Konzola zat´ıˇzen´ a vlastn´ı t´ıhou. . . . . . . . . . . . . . . . . Konvergence hodnoty kritické s´ıly volného prutu. . . . . . . N´ ar˚ ust v´ ypoˇcetn´ı doby s poˇctem stupˇ n˚ u volnosti. . . . . . . Konvergence hodnoty kritické s´ıly konzoly. . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

54 54 55 57 57 58

4.1 4.2 4.3 4.4 4.5 4.6

Maxim´ aln´ı ˇcasov´ y krok v závislosti na poˇctu d´ılc˚ u a metodˇe. V´ ypoˇcetn´ı doba v z´ avislosti na poˇctu d´ılc˚ u a metodˇe. . . . . . Z´ akladn´ı grafické prostˇred´ı modelu. . . . . . . . . . . . . . . . Zobrazen´ı stavové promˇenné v závislosti na ˇcase. . . . . . . . Zobrazen´ı stavov´ ych promˇenn´ ych ve vzájemné závislosti. . . . Zobrazen´ı promˇenn´ ych stejného typu. . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

65 66 68 69 69 69

11

. . . . .

. . . . .

Seznam tabulek 3.1 3.2 3.3 3.4 3.5 3.6

Srovn´ an´ı vypoˇcten´ ych v´ ysledk˚ u s analytick´ ym ˇreˇsen´ım vlastn´ıch frekvenc´ı. . . Konvergence hodnoty prvn´ı vlastn´ı frekvence a odpov´ıdaj´ıc´ı v´ ypoˇcetn´ı doba. Porovn´ an´ı v´ ysledk˚ u pr˚ uhybu a pootoˇcen´ı voln´ıho konce konzoly. . . . . . . . . V´ ysledky simulac´ı speci´ aln´ıho modelu volného prutu. . . . . . . . . . . . . . . V´ ysledky simulac´ı proveden´ ych v aplikaci FyDiK. . . . . . . . . . . . . . . . . V´ ysledky simulac´ı speci´ aln´ıho modelu konzoly. . . . . . . . . . . . . . . . . .

53 54 55 56 56 58

4.1 4.2

Maxim´ aln´ı ˇcasov´ y krok v závislosti na poˇctu d´ılc˚ u a metodˇe. . . . . . . . . . V´ ypoˇcetn´ı doba v z´ avislosti na poˇctu d´ılc˚ u a metodˇe. . . . . . . . . . . . . . .

65 66

12

Kapitola 1

Teoretick´ yu ´ vod Vu ´vodn´ı ˇc´ asti pr´ ace je zapotˇreb´ı pohovoˇrit o teori´ıch a pojmech, jenˇz jsou v dalˇs´ıch ˇcástech práce pouˇzity. Pˇrestoˇze lze u ˇcten´ aˇre pˇredpokládat znalosti fyziky, matematiky a stavebn´ı mechaniky, je na m´ıstˇe mu nab´ıdnout kompaktn´ı teoretick´ y u ´vod, kter´ y nav´ıc pˇredejde pˇr´ıpadnému nedorozumˇen´ı.

1.1 1.1.1

Klasick´ a mechanika ´ Uvod a historie

Mechanika je vˇedn´ı obor zab´ yvaj´ıc´ı se studiem pohybu tˇeles v prostoru a ˇcase. Jej´ı souˇcást, klasická mechanika, zkoum´ a pohyb tˇeles a objekt˚ u v makroskopickém mˇeˇr´ıtku 1 pˇri rychlosti v´ yraznˇe niˇzˇs´ı neˇz rychlost svˇetla, viz. obrázek 1.1. Rychlost

10 -9m

3×10 8 ms-1

Klasická mechanika

Relativistická mechanika

10 -9m

Rozměry

3×10 8 ms-1

Kvantová mechanika

Kvantová teorie pole

Obr´ azek 1.1: Zatˇr´ıdˇen´ı klasické mechaniky.

S mechanick´ ymi jevy se ˇclovˇek pot´ yká od sam´ ych poˇcátk˚ u vlastn´ı existence. Napˇr´ıklad statická rovnov´ aha, princip p´ aky, kladky, naklonˇená rovina a dalˇs´ı jednoduché stroje, byly ˇ pˇredmˇetem b´ ad´ an´ı vˇedc˚ u a vyn´ alezc˚ u jiˇz ve starovˇekém Recku. Jedn´ım z nejv´ yznamnˇejˇs´ıch ´des, kter´ byl zcela jistˇe Archime y nav´ıc definoval pojmy jako tˇeˇziˇstˇe a statick´ y moment. 1

Rozmˇery tˇelesa jsou dosti velké pro pozorov´ an´ı pouh´ ym okem.

13

´ u ´ vod Teoreticky

Klasick´ a mechanika v dneˇsn´ı podobˇe se zaˇcala utváˇret v 16. stolet´ı s nástupem tzv. modern´ı empirické vˇedy, jenˇz se op´ırá pˇredevˇs´ım o experiment a mˇeˇren´ı. Zde je vhodné poznamenat, ˇze patrnˇe mnohem vˇetˇs´ı nesnáz´ı neˇz nevyvinutá technologie byl pro tehdejˇs´ı vˇedce mnohdy zatrvrzel´ y odpor spoleˇcnosti. Mezi nejvˇetˇs´ı osobnosti té doby patˇr´ı napˇr´ıklad Galileo Galilei a Blaise Pascal. Prav´ ym otcem klasické mechaniky je potom Isaac Newton, jenˇz v roce 1687 ve své práci Philosophiae Naturalis Principia Mathematica publikoval z´ akon vˇseobecné gravitace a tˇri pohybové zákony, viz [18].

1.1.2

Pravo´ uhl´ a souˇ radn´ a soustava

Poloha bodu v pravo´ uhl´ e souˇ radn´ e soustavˇ e Kaˇzdému bodu P v trojrozmˇerném prostoru m˚ uˇze b´ yt pˇriˇrazen vektor r, kter´ y urˇcuje jeho polohu vzhledem ke zvolenému poˇcátku O souˇradné soustavy xyz, viz obrázek 1.2. Vektor r lze formulovat jako: b r = x bi + y bj + z k,

(1.1)

b jsou potom tzv. jednotkové vektory, které reprezentuj´ı osy kde x, y, z jsou re´ aln´ a ˇc´ısla, bi, bj a k pravo´ uhlého souˇradného systému. Pro tyto vektory plat´ı:       0 0 1 b = 0 . bi = 0 , bj = 1 , k (1.2) 1 0 0

osa z P

r z O

osa x y x

osa y Obr´ azek 1.2: Poloha bodu P v souˇradné soustavˇe.

14


Rychlost bodu v pravo´ uhl´ e souˇ radn´ e soustavˇ e Rychlost pohybu bodu P v prostoru lze chápat jako velikost okamˇzité zmˇeny vektoru r popisuj´ıc´ıho jeho polohu. Vektor rychlosti v lze potom zapsat jako: dr ∆r = lim , dt t→0 ∆t kde ∆r je zmˇena vektoru r za ˇcasov´ y okamˇzik ∆t, tedy: v=

(1.3)

∆r = r(t + ∆t) − r(t).

(1.4)

Pˇred dosazen´ım do v´ yrazu (1.3) pˇripomeˇ nme, ˇze obecnˇe plat´ı: d da db (a + b) = + . dt dt dt Potom je jiˇz moˇzné vyj´ adˇrit vektor rychlosti v: dr dx b dy b dz b = i+ j+ k, dt dt dt dt ˇcemuˇz odpov´ıd´ a:

(1.5)

v=

(1.6)

b v = vx bi + vz bj + vy k.

(1.7)

Pro dalˇs´ı potˇreby jeˇstˇe poznamenejme, ˇze tedy plat´ı: vx =

dy dz dx , vy = , vz = . dt dt dt

(1.8)

Pohybliv´ a souˇ radn´ a soustava Pˇri pozorov´ an´ı chov´ an´ı pohybuj´ıc´ıho se tˇelesa m˚ uˇze b´ yt v´ yhodné jeho polohu nebo polohu jeho ˇcásti vztahovat k vhodnˇe um´ıstˇené pohyblivé souˇradné soustavˇe. Uvaˇzujme pohybuj´ıc´ı se pravo´ uhlou souˇradnou soustavu x0 y 0 z 0 a soustavu xyz, která je v klidu. Polohu bodu P vzhledem k poˇcátku O soustavy xyz urˇcuje vektor r, polohu poˇcátku O0 soustavy x0 y 0 z 0 potom vektor s, viz obrázek 1.3. Polohu bodu P v˚ uˇci pohybuj´ıc´ı se soustavˇe x0 y 0 z 0 lze popsat vektorem r0 : r0 = r − R(t).

(1.9)

Vektor rychlosti pohybu v0 bodu P , tedy zmˇenu jeho polohy v ˇcase, lze vyjádˇrit snadno: dr0 dr dR(t) = − . dt dt dt Zkrácenˇe tedy: v0 =

(1.10)

v0 = v − V,

(1.11)

kde V je vektor rychlosti pohybu poˇcátku O0 soustavy x0 y 0 z 0 vzhledem k soustavˇe xyz. Pˇri ˇreˇsen´ı u ´lohy nejen v mechanice je vhodné zvolit poˇcátek souˇradného systému tak, aby bylo pozorov´ an´ı chov´ an´ı modelu co moˇzná nejsnazˇs´ı. Jak bude ukázáno pozdˇeji, v pˇredkládaném modelu je pouˇzita stacionárn´ı souˇradná soustava pro sledován´ı polohy modelu v prostoru a z´ aroveˇ n pohybliv´ a souˇradná soustava, jej´ıˇz v´ yhodné pouˇzit´ı bude vysvˇetleno pozdˇeji.

15


z

z' P

r' r s

O'

x'

y'

O

x

y Obr´ azek 1.3: Poloha bodu P vzhledem k pohybujic´ı se souˇradné soustavˇe x0 y 0 z 0 .

1.1.3

Pol´ arn´ı souˇ radn´ a soustava

Jak bylo uk´ az´ ano, v pravo´ uhlé souˇradné soustavˇe lze snadno popsat polohu i rychlost bodu. U nˇekter´ ych u ´loh je ale v´ yhodnˇejˇs´ı zvolit polárn´ı souˇradnou soustavu, jenˇz m˚ uˇze b´ yt pˇrehlednˇejˇs´ı a ˇreˇsen´ı m˚ uˇze znaˇcnˇe zjednoduˇsit. Nam´ısto pravo´ uhl´ ych souˇradnic je poloha bodu P v polárn´ı souˇradné soustavˇe urˇcena u ´hlem ϕ a pr˚ uvodiˇcem r, viz obr´ azek 1.4. y

r

ω

yP

P

r φ

0

xP

x

Obr´ azek 1.4: Poloha bodu P v polárn´ı souˇradné soustavˇe.

Vektor ω je vektorem u ´hlové rychlosti kolmo na pr˚ uvodiˇc r. Vektor r potom znázorˇ nuje zmˇenu délky pr˚ uvodiˇce r v ˇcase pˇri konstantn´ım u ´hlu ϕ. Lze tedy psát: dϕ , dt d |r| |r| = . dt |ω| =

(1.12)

Rychlost pohybu bodu P v rovinˇe lze popsat vektorem v, kter´ y vznikne sloˇzen´ım vektor˚ u ω a r, tedy: v = ω + r.

(1.13)

16


Pro pˇrevod mezi pravo´ uhl´ ymi a polárn´ımi souˇradnicemi bodu P lze snadno odvodit vztahy: xP = |r| cos ϕ,

(1.14)

yP = |r| sin ϕ, pˇriˇcemˇz: q |r| = x2P + yP2 , xP . ϕ = arctan yP

1.1.4

(1.15)

Newtonovy pohybov´ e z´ akony

Jak jiˇz bylo zm´ınˇeno, Isaac Newton v roce 1687 vyslovil tˇri pohybové zákony, jenˇz jsou základem klasické mechaniky a dynamiky. Zákony popisuj´ı vztahy mezi pohybem tˇelesa a silami na tˇeleso p˚ usob´ıc´ımi. Prvn´ı Newton˚ uv z´ akon Prvn´ı Newton˚ uv z´ akon, naz´ yvan´ y také jako Zákon setrvaˇcnosti, je moˇzné v pˇrekladu zapsat: ”Jestliˇze na tˇeleso nep˚ usob´ı ˇz´ adné vnˇejˇs´ı s´ıly nebo je výslednice vnˇejˇs´ıch sil nulov´ a, potom tˇeleso setrv´ av´ a v klidu nebo v rovnomˇerném pˇr´ımoˇcarém pohybu.” Lze si vˇsimnout, ˇze z´ akon zmiˇ nuje pouze vnˇejˇs´ı s´ıly, které p˚ usob´ı na tˇeleso. Je tedy zˇrejmé, ˇze vnitˇrn´ı s´ıly mezi ˇc´ astmi tˇelesa nemaj´ı na pohyb tˇelesa jako celku, pˇresnˇeji na pohyb jeho tˇeˇziˇstˇe, vliv. Matematicky lze Z´ akon setrvaˇcnosti formulovat napˇr´ıklad takto: X dv Fext = 0 ⇔ = 0. (1.16) dt Druh´ y Newton˚ uv z´ akon Druh´ y Newton˚ uv z´ akon, zn´ am´ y také jako Zákon s´ıly, ˇr´ıká: ”P˚ usob´ı-li na tˇeleso s´ıla, pak se tˇeleso pohybuje se zrychlen´ım, které je pˇr´ımo u ´mˇerné p˚ usob´ıc´ı s´ıle a nepˇr´ımo u ´mˇerné hmotnosti tˇelesa.” Tvrzen´ı z´ akona lze matematicky zapsat jako: Fext , (1.17) m kde a je vektor zrychlen´ı tˇelesa a m je hmotnost tˇelesa. Je zˇrejmé, ˇze v´ yraz (1.17) m˚ uˇzeme pˇrepsat do zn´ amé podoby: a=

Fext = m a.

(1.18)

Rovnici (1.18) lze jeˇstˇe d´ ale upravit: dv d(mv) dp = = , dt dt dt tedy m˚ uˇzeme ˇr´ıci, ˇze s´ıla F je rovna ˇcasové zmˇenˇe hybnosti p. Fext = m

17

(1.19)


Tˇ ret´ı Newton˚ uv z´ akon Prvé dva Newtonovy z´ akony se zab´ yvaly odezvou tˇelesa na p˚ usoben´ı vnˇejˇs´ıch sil. Tˇret´ı Newton˚ uv z´ akon, naz´ yvan´ y také jako Zákon akce a reakce, hovoˇr´ı o vzájemném p˚ usoben´ı dvou tˇeles: ”Proti kaˇzdé akci vˇzdy p˚ usob´ı stejn´ a reakce; jinak: vz´ ajemn´ a p˚ usoben´ı dvou tˇeles jsou vˇzdy stejnˇe velk´ a a m´ıˇr´ı na opaˇcné strany.” Tˇret´ı pohybov´ y z´ akon vyj´ adˇren´ y matematicky potom vypadá následovnˇe: FA,B = −FB,A ,

(1.20)

kde FA,B je s´ıla vyv´ıjen´ a tˇelesem A na tˇeleso B a FB,A je s´ıla, kterou tˇeleso B zpˇetnˇe p˚ usob´ı na tˇeleso A. Tedy: X FA,B + FB,A = 0 ; Fext = 0. (1.21)

1.1.5

Z´ akon zachov´ an´ı hybnosti

Uvaˇzujme uzavˇrenou soustavu2 obsahuj´ıc´ı dvˇe na sebe navzájem silovˇe p˚ usob´ıc´ı hmotn´ a tˇelesa A a B. Podle Tˇret´ıho Newtonova zákona plat´ı, ˇze p˚ usob´ıc´ı s´ıly maj´ı stejnou velikost a opaˇcn´ y smˇer. Rovnici (1.20) lze dále upravit: mA aA = −mB aB , mA

dvA dvB = −mB , dt dt dpA dp = − B, dt dt

(1.22)

d (p + pB ) = 0, dt A z ˇcehoˇz je zˇrejmé, ˇze je-li silov´ a výslednice soustavy nulov´ a, celkov´ a hybnost soustavy z˚ ust´ av´ a konstantn´ı, viz [13]. Pro soustavu obsahuj´ıc´ı n hmotn´ ych bod˚ u analogicky plat´ı: P=

˙ = P

n X i=1

i=1

n X

n X

i=1 n X

pi =

n X

p˙i =

i=1

mi vi ,

(1.23)

n

X dvi X mi = mi ai = Fi , dt

(1.24)

i=1

˙ = 0 ⇔ P = konst, Fi = 0 ⇔ P

(1.25)

i=1

kde pi jsou d´ılˇc´ı hybosti hmotn´ ych bod˚ u soustavy a P je celková hybnost soustavy. 2 Uzavˇren´ a soustava tˇeles nebo hmotn´ ych bod˚ u je takov´ a soustava, kde na sebe obsaˇzené entity mohou p˚ usobit silovˇe, popˇr. si vymˇen ˇovat energii, ale nemohou s okol´ım vymˇen ˇovat hmotu.

18


1.1.6

Raketov´ y pohon

Pˇredlohou pˇredkl´ adaného dynamického modelu je ˇst´ıhlá raketa zat´ıˇzená tahem raketového motoru, proto je na m´ıstˇe zm´ınit i princip fungován´ı raketového pohonu. Raketov´ y pohon se ˇrad´ı mezi z´ akladn´ı typy reaktivn´ıho pohonu. Dalˇs´ımi typy reaktivn´ıho pohonu jsou: • proudov´ y pohon, • náporov´ y pohon, • pulzaˇcn´ı pohon, • iontov´ y pohon. Princip reaktivn´ıho pohonu vycház´ı z Tˇret´ıho Newtonova zákona a Zákona zachov´ an´ı hybnosti, viz [19]. Reaktivn´ı motory vyv´ıj´ı tah d´ıky tryskán´ı média z hnac´ı trysky motoru. V pˇr´ıpadˇe raketového motoru se jedná o spaliny vznikaj´ıc´ı hoˇren´ım paliva. Uvaˇzujme raketu o hmotnosti m pohybuj´ıc´ı se rychlost´ı v, viz obrázek 1.5. Rychlost unikaj´ıc´ıch spalin oznaˇcme jako vex .

vex

v m

Obr´ azek 1.5: Raketa zat´ıˇzená tahem motoru.

Pˇri hoˇren´ı paliva a trysk´ an´ı spalin raketa neustále ztrác´ı hmotnost. Hybnost rakety p v ˇcase t lze zapsat jako: p(t) = m(t)v(t).

(1.26)

V ˇcase (t + dt) se hmotnost rakety sn´ıˇz´ı na (m(t) − dm). Hybnost rakety se zmˇen´ı na: p(t + dt) = (m(t) − dm) (v(t) + dv) .

(1.27)

Ztracené spaliny maj´ı hmotnost dm a rychlost (v(t) − vex ). Celková hybnost rakety a ztracen´ ych spalin P v ˇcase (t + dt) je potom: P(t + dt) = (m(t) − dm) (v(t) + dv) + dm(v(t) − vex ).

(1.28)

Zmˇena celkové hybnosti dP m´ a tedy tvar: dP = P(t + dt) − P(t) = m dv + dm vex .

(1.29)

Pokud na soustavu nep˚ usob´ı ˇzádné vnˇejˇs´ı s´ıly, jej´ı hybnost z˚ ustává konstantn´ı a jej´ı zmˇena v ˇcase dP je nulov´ a. Lze tedy tvrdit, ˇze: m dv = −dm vex .

(1.30)

19


Pokud vydˇel´ıme obˇe strany rovnice (1.30) ˇclenem dt, z´ıskáme v´ yraz: m

dv dm =− vex . dt dt

(1.31)

Tedy: m a = −m ˙ vex ,

(1.32)

kde (−m) ˙ je rychlost, s jakou raketa ztrác´ı hmotnost hoˇr´ıc´ıho paliva. Podle Druhého Newtonova z´ akona, viz rovnice (1.18), lze levou stranu v´ yrazu (1.32) oznaˇcit za hnac´ı s´ılu, tedy tah raketového motoru Fthrust : Fthrust = −m ˙ vex .

(1.33)

Separace promˇenn´ ych v rovnici (1.30) dále vede na v´ yraz: dv = −

dm vex . m

(1.34)

Pokud je rychlost trysk´ an´ı spalin vex konstantn´ı, integrac´ı obdrˇz´ıme v´ ysledek: m 0 , v − v0 = vex ln m

(1.35)

kde v0 je poˇc´ ateˇcn´ı rychlost rakety a m0 je poˇcáteˇcn´ı celková hmotnost rakety. Bliˇzˇs´ı pohled na v´ yraz (1.35) ukazuje rychlostn´ı omezen´ı raketového pohonu. Pomˇer m0 /m je maxim´ aln´ı ve chv´ıli, kdy je veˇskeré palivo spáleno a m obsahuje pouze hmotnost rakety a n´ akladu. I v krajn´ım pˇr´ıpadˇe, kdy by poˇcáteˇcn´ı hmotnost rakety byla tvoˇrena z 90% hmotnost´ı paliva, je nejvyˇsˇs´ı pomˇer m0 /m roven 10, tedy maximáln´ı z´ıskaná rychlost (v − v0 ) bude pˇribliˇznˇe 2.3 vex . Z uvedeného vypl´ yv´ a, ˇze je vhodné maximalizovat rychlost tryskán´ı spalin vex a d´ ale s v´ yhodou navrhnout raketu o v´ıce stupn´ıch, jeˇz mohou b´ yt v pr˚ ubˇehu letu odhozeny, ˇc´ımˇz se sn´ıˇz´ı celkov´ a hmotnost rakety ve vztahu ke stupˇ n˚ um následuj´ıc´ım.

20


1.2

Energie

Energie je jednou ze z´ akladn´ıch fyzikáln´ıch veliˇcin. Energie je skalárn´ı veliˇcinou, je tedy urˇcena pouze velikost´ı. V pˇr´ırodˇe lze energii popsat v mnoha podobách, viz [19]. Uved’me ty hlavn´ı: • kinetick´ a energie, • potenci´ aln´ı energie, – gravitaˇcn´ı potenci´ aln´ı energie, – potenci´ aln´ı energie vzniklá pruˇznou deformac´ı, – tlakov´ a potenci´ aln´ı energie, • elektrick´ a energie, • magnetick´ a energie, • vnitˇrn´ı energie, – tepeln´ a energie, – jadern´ a energie, – chemick´ a energie.

1.2.1

Potenci´ aln´ı energie

Potenciáln´ı energie Ep (v zahraniˇcn´ı literatuˇre oznaˇcovaná jako V ) je energie, jenˇz se m˚ uˇze kumulovat v systému, viz [19] a [20]. Pˇr´ıkladem m˚ uˇze b´ yt gravitaˇcn´ı potenciáln´ı energie, kdy je potenci´ aln´ı energie u ´mˇern´ a hmotnosti tˇelesa m, gravitaˇcn´ımu zrychlen´ı g a v´ yˇsce tˇeˇziˇstˇe tˇelesa nad srovn´ avac´ı rovinou h, tedy: Ep = m g h.

(1.36)

Dalˇs´ı moˇznost´ı kumulace energie je pruˇzná deformace tˇelesa. Pˇri protaˇzen´ı nebo stlaˇcen´ı translaˇcn´ı pruˇziny tuhosti k a poˇc´ ateˇcn´ı délky l o délku ∆l, viz obrázek 1.6, p˚ usob´ı na pruˇzinu s´ıla Fts : Fts = k ∆l.

(1.37)

Pro potenci´ aln´ı energii translaˇcn´ı pruˇziny Ep,ts plat´ı: Z ∆l 1 Ep,ts = Fts dl = k ∆l2 . 2 0

(1.38)

Podobnˇe v pˇr´ıpadˇe rotaˇcn´ı pruˇziny, viz obrázek 1.6, spojuj´ıc´ı dva tuhé d´ılce pˇri deformaci ou ´hel ϕ p˚ usob´ı moment Mrs : Mrs = k ϕ.

(1.39)

V´ yraz pro potenci´ aln´ı energii Ep,rs nashromáˇzdˇenou v rotaˇcn´ı pruˇzinˇe má potom tvar: Z ϕ 1 Ep,rs = Mrs dϕ = k ϕ2 . (1.40) 2 0

21


Mrs k

l+Δl

k

φ

Fts

Obr´ azek 1.6: Translaˇcn´ı a rotaˇcn´ı pruˇzina.

1.2.2

Kinetick´ a energie

Kinetick´ a energie je moˇzn´ a nejbˇeˇznˇejˇs´ı forma energie. Pro hmotn´ y bod o hmotnosti m a rychlosti v lze v´ yraz pro jeho kinetickou energii Ek (v zahraniˇcn´ı literatuˇre oznaˇcovanou také T ) formulovat jako: 1 Ek = m v 2 . 2

1.2.3

(1.41)

Z´ akon zachov´ an´ı energie

Zákon zachov´ an´ı energie ˇr´ık´ a, ˇze celková energie izolované soustavy3 z˚ ustává konstantn´ı. Energie tedy nem˚ uˇze b´ yt nijak vyrobena ani zniˇcena, pouze se m˚ uˇze pˇremˇen ˇovat z jedné formy do jiné. Uvaˇzujme kyvadlo o jednom stupni volnosti, viz obrázek 1.7, s nehmotn´ ym tuh´ ym závˇesem délky l, na jehoˇz konci je soustˇredˇená hmota m. Na hmotu m p˚ usob´ı gravitaˇcn´ı zrychlen´ı o velikosti g. x

φ0 l m

h

y

Obr´ azek 1.7: Tuhé kyvadlo s jedn´ım stupnˇem volnosti.

3 Izolovan´ a soustava tˇeles nebo hmotn´ ych bod˚ u je takov´ a soustava, kde na sebe obsaˇzené entity mohou p˚ usobit silovˇe, ale nemohou s okol´ım vymˇen ˇovat energii nebo hmotu.

22


Kyvadlo je v ˇcase t0 = 0 vych´ yleno od svislého rovnováˇzného stavu o u ´hel ϕ0 . Jak bylo uvedeno, jeho potenci´ aln´ı energie Ep (t0 ) má velikost: Ep (t0 ) = m g h = m g l sin ϕ(t0 )

(1.42)

Rychlost z´ avˇesu kyvadla v lze vyjádˇrit pomoc´ı jeho u ´hlové rychlosti jako: v = l ω,

(1.43)

kde ω je u ´hlov´ a rychlost kyvadla, tedy: ω=

dϕ . dt

(1.44)

Pro kinetickou energii kyvadla Ek plat´ı vztah (1.41). Po u ´pravˇe z´ıskáváme v´ yraz: 1 Ek = m l2 ω 2 . 2

(1.45)

V ˇcase t0 m´ a kyvadlo nulovou u ´hlovou rychlost, jeho kinetická energie Ek je proto nulov´ a. ´ Uhlová v´ ychylka kyvadla ϕ (a tedy v´ yˇska h) je maximáln´ı, coˇz znamená i maximáln´ı hodnotu potenciáln´ı energie. Celkov´ a energie kyvadla E je prost´ ym souˇctem jeho potenciáln´ı a kinetické energie: E = Ep + Ek = konst.

(1.46)

Pro energii kyvadla v ˇcase t0 tedy plat´ı: E = Ep (t0 ) + Ek (t0 ) = m g l sin ϕ(t0 ) + 0.

(1.47)

Po uvolnˇen´ı kyvadla zaˇcne v´ ychylka ϕ i v´ yˇska kyvadla h klesat, naopak rychlost kyvadla bude r˚ ust. V libovolném okamˇziku t1 lze vystihnout energii kyvadla v´ yrazem (1.46), plat´ı tedy: 1 E = Ep (t1 ) + Ek (t1 ) = m g l sin ϕ(t1 ) + m l2 ω(t1 )2 . 2

(1.48)

V okamˇziku t2 , kdy je kyvadlo v nejniˇzˇs´ı poloze, je v´ yˇska kyvadla h nad srovnávac´ı rovinou nulov´ aau ´hlov´ a rychlost kyvadla ω je maximáln´ı. Veˇskerá energie kyvadla je tedy ve formˇe kinetické energie. Lze proto psát: 1 E = Ep (t2 ) + Ek (t2 ) = 0 + m l2 ω(t2 )2 . 2

(1.49)

Protoˇze uvaˇzovan´ y mechanick´ y systém je konzervativn´ı a nedocház´ı tedy ke ztrát´ am energie, popsan´ y cyklus pohybu kyvadla a pˇremˇen ˇován´ı energie z potenciáln´ı na kinetickou se neust´ ale opakuje.

23


1.3

Lagrangeovsk´ a mechanika

V pr˚ ubˇehu let 1772 aˇz 1788 Joseph-Louis Lagrange publikoval svoji formulaci klasické Newtonovské mechaniky. Lagrangeova formulace je ekvivalentn´ım popisem stavu systému v porovnán´ı s Newtonovsk´ ym pˇr´ıstupem. Na rozd´ıl od Newtona, vyuˇz´ıvaj´ıc´ıho k popisu pohybu vektory, Lagrange pouˇz´ıv´ a skalárn´ı veliˇciny - kinetickou a potenciáln´ı energii systému. Lagrange˚ uv pˇr´ıstup rovnˇeˇz pouˇz´ıv´ a efektivnˇejˇs´ıho popisu stavu systému, totiˇz konfiguraˇcn´ı prostor zobecnˇen´ ych souˇradnic. Lagrangeova formulace je, jak bude pozdˇeji ukázáno, velmi elegantn´ım pˇr´ıstupem popisu dynamiky systému, viz [7] Pomoc´ı jednoduché poˇcáteˇcn´ı rovnice, obsahuj´ıc´ı pouze skalárn´ı veliˇciny, lze relativnˇe snadno sestrojit pohybové rovnice popisuj´ıc´ı chován´ı systému.

1.3.1

Zobecnˇ en´ e souˇ radnice syst´ emu

Jak bylo zm´ınˇeno, Lagrangeovsk´ y popis stavu systému nepouˇz´ıvá k popisu stavu systému nutnˇe kartézskou souˇradnou soustavu. K popsán´ı stavu systému lze zvolit libovolné parametry, které v kaˇzdém okamˇziku jednoznaˇcnˇe popisuj´ı stav sledovaného systému. Pouˇzit´ı vhodnˇe zvolen´ ych zobecnˇen´ ych souˇradnic b´ yvá v´ yhodnˇejˇs´ı také z d˚ uvodu znaˇcné rozmanitosti u ´loh v mechanice. Pokud m´ a sledovan´ y systém n stupˇ nu volnosti4 , pro jednoznaˇcn´ y popis stavu systému je zˇrejmˇe zapotˇreb´ı alespoˇ n n zobecnˇen´ ych souˇradnic. Poznamenejme, ˇze pouˇzit´ı v´ıce zobecnˇen´ ych souˇradnic jiˇz d´ ale nepˇrináˇs´ı ˇzádné v´ yhody.

1.3.2

Konfiguraˇ cn´ı prostor, zobecnˇ en´ e rychlosti syst´ emu

Zobecnˇené souˇradnice systému tvoˇr´ı tzv. konfiguraˇcn´ı prostor vˇsech pˇr´ıpustn´ ych poloh neboli konfigurac´ı systému. Takto popsán, konfiguraˇcn´ı prostor neposkytuje u ´plnou informaci o fyzik´ aln´ım stavu systému, viz [7]. V dynamické u ´loze je zapotˇreb´ı znát nav´ıc rychlosti hmotn´ ych bod˚ u, tedy zobecnˇené rychlosti systému. Oznaˇc´ıme-li zobecnˇené souˇradnice systému jako q0 , q1 , . . . , qn−1 , pro jejich zobecnˇené rychlosti q˙0 , q˙1 , . . . , q˙n−1 plat´ı: q˙i =

1.3.3

dqi . dt

(1.50)

Lagrangeova funkce

Uvaˇzujme hmotn´ y bod M o hmotnosti m pohybuj´ıc´ı se prostorem. Kinetická energie hmotného bodu Ek m´ a zˇrejmˇe tvar, viz 1.41. Kinetická energie je tedy funkc´ı rychlosti hmotného bodu: 1 1 1 Ek = m v2 = m r˙ 2 = m x˙ 2 + y˙ 2 + z˙ 2 = Ek (x, ˙ y, ˙ z) ˙ . 2 2 2

(1.51)

Potenci´ aln´ı energie hmotného bodu Ep nab´ yvá tvaru, viz 1.36. Závis´ı tedy na poloze hmotného bodu v prostoru: Ep = m g h(r) = Ep (x, y, z).

(1.52)

4

Stupeˇ n volnosti je smˇer posunu nebo pootoˇcen´ı, jimiˇz se hmotn´ y bod nebo tˇeleso m˚ uˇze pohybovat. V prostoru obecnˇe n = 3N − v, kde N je poˇcet hmotn´ ych bod˚ u a v je poˇcet vazeb.

24


Lagrangeova funkce (nebo zkr´ acenˇe Lagrangián) L je prost´ ym rozd´ılem dvou skalár˚ ukinetické a potenci´ aln´ı energie systému. Tedy: L = Ek − Ep ,

(1.53)

v zahraniˇcn´ı literatuˇre pak: L = T − V.

(1.54)

Lagrangeova funkce L je funkc´ı rychlosti a polohy hmotného bodu, tedy vˇsech promˇenn´ ych popisuj´ıc´ıch dynamick´ y stav systému: L = L (t, x, y, z, x, ˙ y, ˙ z) ˙ .

(1.55)

Pro obecn´ y systém popsan´ y pomoc´ı n zobecnˇen´ ych souˇradnic a rychlost´ı: ˙ . L = L (t, q0 , q˙0 , q1 , q˙1 , . . . , qn−1 , q˙n−1 ) = L (t, q, q)

1.3.4

(1.56)

Virtu´ aln´ı posunut´ı, virtu´ aln´ı pr´ ace

Virtuáln´ı posunut´ı hmotného bodu δr je libovolné myˇslené posunut´ı, které vyhovuje okrajov´ ym podm´ınk´ am u ´lohy, viz [11] a [12]. M˚ uˇze j´ım b´ yt virtuáln´ı posun δs nebo virtuáln´ı pootoˇcen´ı δϕ. Virtu´ aln´ı pr´ ace δW je pr´ ace skuteˇcn´ ych sil nebo moment˚ u na virtuáln´ım pˇrem´ıstˇen´ı, tedy napˇr´ıklad: δW = F δs,

(1.57)

δW = M δϕ.

Virtuáln´ı pr´ ace je skal´ arn´ı veliˇcinou, nebot’ vzniká skalárn´ım souˇcinem vektor˚ u. Uvaˇzujme hmotn´ y bod, na kter´ y p˚ usob´ı v´ yslednice vnˇejˇs´ıch sil F. Vektor s´ıly F lze rozloˇzit ve smˇeru souˇradnicov´ ych os x, y, z na d´ılˇc´ı vektory Fx , Fy , Fz . Podobnˇe virtuáln´ı posun δs se sklád´ a ze sloˇzek δsx , δsy , δsz . Virtuáln´ı práce δW v´ yslednice vnˇejˇs´ıch sil F na virtuáln´ım posunut´ı δs je potom d´ ana vztahem: δW = Fx δsx + Fy δsy + Fz δsz .

(1.58)

Pro statickou rovnov´ aˇznou soustavu sil dále plat´ı Lagrange˚ uv princip virtuáln´ıch prac´ı: ”Virtu´ aln´ı pr´ ace δW rovnov´ aˇzné soustavy sil Fi (i = 0, 1, . . . , n) p˚ usob´ıc´ı na hmotný bod nebo tuhé tˇeleso je pˇri libovolném virtu´ aln´ım posunut´ı δr rovna nule.” Matematicky zaps´ ano: δW =

n X

Fi δri = 0.

(1.59)

i=0

25


1.3.5

D’Alembert˚ uv princip

Sestavme pohybovou rovnici pro hmotn´ y bod M o hmotnosti m, na kter´ y p˚ usob´ı v´ yslednice vnˇejˇs´ıch sil F a v´ yslednice reakˇcn´ıch sil R: F + R = m a = m ¨r,

(1.60)

kde r je polohov´ y vektor hmotného bodu M . Pravou stranu rovnice (1.60) lze oznaˇcit jako setrvaˇcnou s´ılu Z, pro kterou plat´ı: ˙ Z = −m a = −p.

(1.61)

Po dosazen´ı tedy lze ps´ at: F + R + Z = 0,

(1.62)

coˇz je matematick´ ym z´ apisem D’Alembertova principu, jenˇz má p˚ uvod v Newtonov´ ych pohybov´ ych z´ akonech: ”Dynamický systém je v rovnov´ aze pod vlivem sil vnˇejˇs´ıch, reakˇcn´ıch a setrvaˇcných.” Dále uˇzijme principu virtu´ aln´ıch prac´ı. Pokud budeme uvaˇzovat ideáln´ı vazby, pro které plat´ı: δWR = R δs = 0,

(1.63)

z´ıskáváme v´ yraz: δW = F δs + Z δs = 0,

(1.64)

˙ = 0, δW = F δs − (m ¨r) δs = F δs − pδs

(1.65)

tedy:

coˇz je formulace D’Alembertova prinicpu ve formˇe virtuáln´ıch prac´ı: ”Rozd´ıl virtu´ aln´ıch prac´ı, konaných vnˇejˇs´ımi silami p˚ usob´ıc´ımi na hmotný bod na virtu´ aln´ım posunut´ı δs, a virtu´ aln´ıch prac´ı, konaných silami setrvaˇcnými na témˇze posunut´ı, je roven nule. ”

1.3.6

Hamilton˚ uv pricip nejmenˇ s´ı akce

Definujme nejprve akci systému S jako integrál Lagrangeovy funkce L, viz (1.56), mezi dvˇema ˇcasov´ ymi okamˇziky ta a tb : Z tb ˙ dt S(ta , tb ) = L (t, q, q) (1.66) ta

Princip nejmenˇs´ı akce se zakl´ adá na tvrzen´ı, ˇze ze vˇsech pˇr´ıpustn´ ych trajektori´ı systému je uskuteˇcnˇena trajektorie nejv´ yhodnˇejˇs´ı, viz [7]. Napˇr´ıklad z hlediska energetické nároˇcnosti, nutného ˇcasu5 a podobnˇe. U mechanického systému pˇredpokládáme, ˇze ze vˇsech moˇzn´ ych trajektorii´ı probˇehne ta, pro kterou integrál (1.66) nab´ yvá minimáln´ı hodnoty. ´ Uveden´ y integr´ al se naz´ yv´ a funkcionál6 . Ukolem je tedy naj´ıt minimáln´ı hodnotu funkcionálu, totiˇz extrém ve stacion´ arn´ım bodˇe. 5

Pˇr´ıkladem je tzv. Fermat˚ uv princip v optice: Svˇetlo se v prostoru ˇs´ıˇr´ı z jednoho bodu do druhého po takové dr´ aze, aby doba potˇrebn´ a k pˇrekon´ an´ı této dr´ ahy byla minim´ aln´ı. 6 Funkcion´ al je zobrazen´ı, které pˇriˇrazuje funkci ˇc´ıslo - energii, ˇcas a podobnˇe.

26


1.3.7

Lagrangeovy rovnice

Hamilton˚ uv princip nejmenˇs´ı akce dále pouˇzijeme k odvozenn´ı Lagrangeov´ ych rovnic. Jak bylo uvedeno, v´ yvoj stavu systému bude prob´ıhat po trajektorii s minimáln´ı akc´ı. Zaved’me virtuáln´ı posunut´ı δq(t) jako infinitezimáln´ı rozd´ıl myˇslené a uskuteˇcnˇené trajektorie systému v libovolném ˇcase t, viz obr´ azek 1.8: δq(t) = qvirt (t) − qreal (t).

(1.67)

tb

uskutečněná trajektorie qreal t

δq t neuskutečněná (virtuální) trajektorie qvirt

ta

Obr´ azek 1.8: Trajektorie v´ yvoje dynamického systému.

Dále uved’me, ˇze virtu´ aln´ı i uskuteˇcnˇená trajektorie zaˇc´ıná a konˇc´ı v totoˇzném bodˇe konfiguraˇcn´ıho prostoru: δq(ta ) = δq(tb ) = 0,

(1.68)

a nakonec z´ amˇennost operac´ı derivace podle ˇcasu a variace δ: d ˙ δq = δ q. dt

(1.69)

Dále budeme hledat minimum integrálu akce S. Stanovme podm´ınky extremálnosti: Z tb ˙ dt = 0, δ L (t, q, q) ta tb

Z

ta Z tb ta

˙ dt = 0, δL (t, q, q) ∂L ∂L δqk + δ q˙k dt = 0. ∂qk ∂ q˙k

(1.70)

Nyn´ı proved’me integraci per partes druhého ˇclenu v´ ysledku (1.70): Z

tb

ta

∂L d δqk − ∂qk dt

∂L ∂ q˙k

δqk

∂L dt + δqk ∂ q˙k

27

tb = 0. ta

(1.71)


S pˇrihlédnut´ım k tvrzen´ı (1.68) lze posledn´ı ˇclen v´ yrazu (1.71) oznaˇcit za nulov´ y. M˚ uˇzeme potom ps´ at: Z tb ∂L d ∂L δqk dt = 0. (1.72) − ∂qk dt ∂ q˙k ta Rovnost (1.72) mus´ı platit pro kaˇzdé ˇcasy ta a tb a pro kaˇzdé virtuáln´ı posunut´ı δqk . Zároveˇ n virtu´ aln´ı posunut´ı δqk jsou na sobˇe nezávislá, protoˇze poˇcet zobecnˇen´ ych souˇradnic je roven poˇctu stupˇ n˚ u volnosti. Pro kaˇzdé k potom mus´ı platit: ∂L d ∂L = 0, k = 0, 1, . . . , n − 1 , (1.73) − ∂qk dt ∂ q˙k kde n je poˇcet zobecnˇen´ ych souˇradnic systému. Rovnice (1.73) se naz´ yvaj´ı Lagrangeovy rovnice, viz [20]. Pˇredstavuj´ı nutné podm´ınky extremálnosti (minimality) funkcionálu akce S. Matematicky se jedná o obyˇcejné diferenciáln´ı rovnice druhého ˇr´ adu pro uskuteˇcnˇenou trajektorii qk , viz [14] a [15].

1.4

Pohybov´ e rovnice a jejich ˇ reˇ sen´ı

1.4.1

Pohybov´ e rovnice

Pohybov´ a rovnice systému je matematick´ y pˇredpis popisuj´ıc´ı moˇzné trajektorie v´ yvoje systému v závislosti na ˇcase, viz [14] a [2]. Pro ˇreˇsen´ı pohybov´ ych rovnic je potˇreba znát poˇcáteˇcn´ı podm´ınky, zde poˇc´ ateˇcn´ı hodnoty stavov´ ych promˇenn´ ych systému. 7 Pro ˇreˇsen´ı disipativn´ıch systém˚ u bˇeˇznˇe postaˇcuj´ı tzv. explicitn´ı numerické metody, viz [16]. Konzervativn´ı8 systémy zpravidla vyˇzaduj´ı robustnˇejˇs´ı numerické metody kv˚ uli nároˇcnosti na numerickou stabilitu ˇreˇsen´ı.

1.4.2

Eulerova metoda

Eulerova metoda je nejjednoduˇsˇs´ı explicitn´ı metodou pro ˇreˇsen´ı obyˇcejn´ ych diferenciáln´ıch rovnic prvn´ıho ˇr´ adu. Eulerova metoda je jednokroková9 a jednobodová10 numerická metoda. Uvaˇzujme diferenci´ aln´ı rovnici ve tvaru: dx = F (t, x(t)). dt Pˇredpis Eulerovy metody je potom:

(1.74)

x(t + h) = x(t) + h F (t, x(t)),

(1.75)

kde h je krok numerické metody a F (t, x(t)) je smˇernice aproximaˇcn´ı pˇr´ımky v bodˇe x(t). Smˇernici F (t, x(t)) lze obdrˇzet jako numerickou prvn´ı derivaci. Pˇribliˇznˇe potom plat´ı: f 0 (x) ≈

f (x + h) − f (x) . h

(1.76)

7

U disipativn´ıho mechanického systému se celkov´ a energie systému nezachov´ av´ a. Disipace energie m˚ uˇze b´ yt zp˚ usobena napˇr´ıklad tlumen´ım systému. 8 V konzervativn´ıch mechanick´ ych systémech z˚ ust´ av´ a celkov´ a energie systému konstantn´ı. 9 Pro ˇreˇsen´ı k-krokové numerické metody je zapotˇreb´ı zn´ at k pˇredchoz´ıch stav˚ u - krok˚ u. 10 Pˇri ˇreˇsen´ı pomoc´ı l-bodové numerické metody je potˇreba spoˇc´ıtat l hodnot funkce F (t, x(t)).

28


Po u ´pravˇe lze tedy ps´ at: f (x + h) ≈ f (x) + h f 0 (x).

(1.77)

Chyba Eulerovy metody je ˇr´ adu O(h2 ). Pro ˇreˇsen´ı dynamického systému s n stavov´ ymi promˇenn´ ymi nab´ yv´ a pˇredpis Eulerovy metody tvaru: xi (t + h) = xi (t) + h Fi (t, x1 (t), x2 (t), . . . , xn (t)).

1.4.3

(1.78)

Semi-implicitn´ı Eulerova metoda

Mnohem stabilnˇejˇs´ıho ˇreˇsen´ı pˇri pouˇzit´ı Eulerovy metody lze dosáhnout jednodouchou u ´pravou vycházej´ıc´ı z vlastnost´ı dynamick´ ych systém˚ u. Dynamické systémy popisuj´ıc´ı chován´ı mechanick´ ych konstrukc´ı maj´ı p˚ uvod v diferenciáln´ıch rovnic´ıch druhého ˇrádu a pomoc´ı jednoduch´ ych substituc´ı se upravuj´ı na diferenciáln´ı rovnice prvn´ıho ˇrádu. Napˇr´ıklad diferenci´ aln´ı rovnici ve tvaru: m

d2 x dx +c + k x = 0, dt2 dt

(1.79)

lze substituc´ı v=

dx , dt

(1.80)

upravit na soustavu dv c k = − v − x, t m m dx = v. dt

(1.81)

Pˇri ˇreˇsen´ı Eulerovou metodou, viz (1.75), dostáváme diferenˇcn´ı vztahy: c k v(t) + x(t)), m m x(t + h) = x(t) + h v(t). v(t + h) = v(t) − h (

(1.82)

Pokud nejdˇr´ıve vyˇreˇs´ıme hodnotu stavové promˇenné v v ˇcase t + h, lze pˇredpis (1.82) upravit na: c k v(t) + x(t)), m m x(t + h) = x(t) + h v(t + h). v(t + h) = v(t) − h (

(1.83)

Popsan´ a u ´prava Eulerovy metody tedy spoˇc´ıvá v pouˇzit´ı novˇe vypoˇctené rychlosti v pro v´ ypoˇcet nového posunut´ı x. Stejné stability v´ ypoˇctu dosáhneme opaˇcnou u ´pravou, totiˇz uˇzit´ım novˇe vypoˇcteného posunut´ı pro v´ ypoˇcet nové rychlosti.

29


1.4.4

Runge-Kuttova metoda

Runge-Kuttova metoda je jednokroková a ˇctyˇrbodová numerická metoda. Pro v´ ypoˇcet je tedy zapotˇreb´ı zn´ at jeden pˇredchoz´ı stav a hodnoty funkce F (t, x(t)) ve ˇctyˇrech bodech. Pr˚ ubˇeh funkce je aproximov´ an polynomem ˇctvrtého stupnˇe. Pˇredpis Runge-Kuttovy metody je následuj´ıc´ı: k1 = F t, x(t) , h h , k2 = F t + , x(t) + k1 2 2 h h (1.84) k3 = F t + , x(t) + k2 , 2 2 k4 = F t + h, x(t) + k3 h , h x(t + h) = x(t) + k1 + 2k2 + 2k3 + k4 . 6 Chyba Runge-Kuttovy metody je ˇrádu O(h5 ). Pro ˇreˇsen´ı soustav rovnic dynamického systému s n stavov´ ymi promˇenn´ ymi lze pˇredpis zobecnit: k1i = Fi t, x1 (t), x2 (t), . . . , xn (t) h h h h k2i = Fi t + , x1 (t) + k11 , x2 (t) + k12 , . . . , xn (t) + k1n 2 2 2 2 h h h h (1.85) k2i = Fi t + , x1 (t) + k21 , x2 (t) + k22 , . . . , xn (t) + k2n 2 2 2 2 k4i = Fi t + h, x1 (t) + k31 h, x2 (t) + k32 h, . . . , xn (t) + k3n h h xi (t + h) = x(t) + k1i + 2k2i + 2k3i + k4i , i = 1, 2, . . . , n. 6

30


1.5

Zdroje nelinearit

Fyzikáln´ı nelinearitou rozum´ıme existenci nelineárn´ı závislosti mezi charakteristikami vnˇejˇs´ıho zat´ıˇzen´ı, napjatosti a deformace konstrukce. U reáln´ ych konstrukc´ı se vˇzdy vyskytuje nelineárn´ı chov´ an´ı, nicménˇe u mnoh´ ych u ´loh je za urˇcit´ ych podm´ınek linearizace ˇreˇsen´ı pˇr´ıpustn´ a a vede ke snazˇs´ımu vyˇreˇsen´ı problému se zanedbatelnou odchylkou oproti skuteˇcnému ˇreˇsen´ı.

1.5.1

Geometrick´ a nelinearita

Geometrick´ a nelinearita ˇreˇsen´ı je zp˚ usobena závislost´ı tvaru (geometrie) konstrukce na zat´ıˇzen´ı nebo okrajov´ ych podm´ınk´ ach. Vznik geometrické nelinearity ukaˇzme na následuj´ıc´ım pˇr´ıkladu. Uvaˇzujme prut délky l o jednom stupni volnosti, kterému v pootoˇcen´ı brán´ı pruˇzina tuhosti k, viz obr´ azek 1.9. Na volném konci prutu p˚ usob´ı osamˇelá s´ıla F . k

F l

φ

F

l cos φ Obr´ azek 1.9: Konzola zat´ıˇzená osamˇelou silou.

Momentov´ a podm´ınka k bodu vetknut´ı má tvar: F l cos ϕ = k ϕ.

(1.86)

Pˇri malém pootoˇcen´ı nosn´ıku lze uvaˇzovat: cos ϕ ≈ 1,

(1.87)

coˇz vede k line´ arn´ı z´ avislosti mezi p˚ usob´ıc´ı silou F a pootoˇcen´ım ϕ. Pro malé deformace lze tento pˇredpoklad pouˇz´ıt, aniˇz bychom se dopustili závaˇzné chyby. S rostouc´ı deformac´ı se ale lineárn´ı ˇreˇsen´ı od pˇresného v´ yznamnˇe odchyluje a nen´ı jiˇz pˇr´ıpustné. Pro pˇresnˇejˇs´ı aproximaci lze nelineárn´ı funkci rozloˇzit do ˇrady, z n´ıˇz se dále pouˇzij´ı prvn´ı dva nebo v´ıce ˇclen˚ u. Pro zm´ınˇen´ y pˇr´ıpad plat´ı: cos ϕ = 1 −

ϕ2 ϕ4 ϕ6 + − + ... 2! 4! 6!

(1.88)

Po dosazen´ı do momentové podm´ınky obdrˇz´ıme: ϕ2 F l 1− = k ϕ. 2!

31

(1.89)


Tato aproximace se jev´ı jako vyhovuj´ıc´ı pro vˇetˇs´ı rozsah u ´hlu ϕ. Pˇri dalˇs´ım nárustu nelinearity se pochopitelnˇe i toto ˇreˇsen´ı zaˇc´ıná liˇsit od pˇresn´ ych hodnot. Uu ´loh, kde vznikaj´ı natolik v´ yznamná pˇretvoˇren´ı, ˇze analytické ˇreˇsen´ı je sloˇzité nebo nemoˇzné, nezb´ yv´ a neˇz ˇreˇsit u ´lohu numericky. Tedy rozdˇelit ˇreˇsen´ı na podkroky a v kaˇzdém uvaˇzovat odpov´ıdaj´ıc´ı geometrii konstrukce.

1.5.2

Materi´ alov´ a nelinearita

Materiálov´ a nelinearita vznik´ a v pˇr´ıpadˇe, pokud materiál modelu nen´ı lineárnˇe pruˇzn´ y s platnost´ı Hookeova z´ akona: σ = E ε.

(1.90)

Závislost mezi napjatost´ı a deformac´ı tedy nen´ı lineárn´ı, viz obrázek 1.10, a vlastnosti materiálu jsou funkc´ı jeho deformace: σ = E(ε) ε.

(1.91)

F, σ

Δl,ε

Obr´ azek 1.10: Nelineárn´ı pracovn´ı diagram materiálu.

Podobnˇe jako v pˇredchoz´ım pˇr´ıpadˇe je tˇreba u ´lohu ˇreˇsit numericky a v kaˇzdém kroku pouˇz´ıt materi´ alové charakteristiky odpov´ıdaj´ıc´ı deformaci modelu.

1.5.3

Nelinearita zp˚ usoben´ a okrajov´ ymi podm´ınkami

Pˇri ˇreˇsen´ı mechanick´ ych u ´loh jsou skuteˇcná uloˇzen´ı a zat´ıˇzen´ı ˇcasto nahrazována idealizovan´ ymi vazbami a silami. Tato zjednoduˇsen´ı ovˇsem nesm´ı b´ yt na u ´kor v´ ystiˇznosti modelu. Ve skuteˇcnosti se uloˇzen´ı a zat´ıˇzen´ı dˇeje jako vzájemn´ y kontakt mezi tˇelesy. Tato skuteˇcnost má za n´ asledek promˇenlivé charakteristiky okrajov´ ych podm´ınek.

32


Jako pˇr´ıklad uved’me zmˇenu polohy teoretick´ ych podpor a smˇeru p˚ usoben´ı jejich reakc´ı, viz obrázek 1.11a, promˇenlivost styˇcn´ ych ploch na kontaktu tˇeles (1.11b) nebo promˇennou tuhost podpory (1.11c). Bˇeˇzné je také napˇr´ıklad nadzdviˇzen´ı konstrukce z podpory nebo naopak kontakt s dalˇs´ım objektem v d˚ usledku jej´ı deformace. a)

b)

c)

Obr´ azek 1.11: Pˇr´ıklady nˇekter´ ych nelineárn´ıch vazeb.

1.6 1.6.1

Konzervativn´ı a nekonzervativn´ı s´ıly Konzervativn´ı s´ıly

Konzervativn´ı silou je kaˇzd´ a s´ıla, která je závislá pouze na své poloze. Je tedy moˇzné v kaˇzdém okamˇziku urˇcit jej´ı potenciál. Jako pˇr´ıklad konzervativn´ı s´ıly lze uvést gravitaˇcn´ı s´ılu, viz sekce 1.2. Pro konzervativn´ı s´ıly plat´ı: • Pˇri pohybu tˇelesa po uzavˇrené kˇrivce, viz obrázek 1.12a, je práce konaná konzervativn´ı silou rovna nule: I W = F dr = 0. (1.92) C

• Pˇri pohybu tˇelesa z jednoho bodu do druhého nen´ı celková práce vykonaná konzervativn´ı silou z´ avisl´ a na trajektorii, viz obrázek 1.12b. • Konzervativn´ı s´ıla je opaˇcn´ a ke gradientu svého potenciálu: F = −∇ Φ.

1.6.2

(1.93)

Nekonzervativn´ı s´ıly

Pro nekonzervativn´ı s´ıly neplat´ı uvedené vlastnosti sil konzervativn´ıch. Pro nekonzervativn´ı s´ıly tedy nelze urˇcit potenci´ al, nebot’ r˚ uzné trajektorie vedou k jiné vykonané práci. Stejnˇe tak jimi vykonan´ a pr´ ace na uzavˇrené kˇrivce nen´ı rovna nule. Zde jako pˇr´ıklad lze uvést tˇrec´ı nebo tlum´ıc´ı s´ıly. Nekonzervativn´ı s´ıly nejsou z´ avislé pouze na poloze svého p˚ usobiˇstˇe, ale také na okamˇzité konfiguraci systému. Nekonzervativn´ı s´ıly, které mˇen´ı své p˚ usoben´ı v závislosti na geometrii konstrukce, se naz´ yvaj´ı sleduj´ıc´ı s´ıly, viz [10]. Reáln´ ym pˇr´ıkladem sleduj´ıc´ıch sil je tryskán´ı média z konzolovˇe uloˇzené trubky nebo právˇe pˇredmˇet této pr´ ace - tah reaktivn´ıho motoru p˚ usob´ıc´ı teˇcnˇe v˚ uˇci konci rakety, viz obrázek 1.13.

33


b) h

a) h

B

A

x

A

x

Obr´ azek 1.12: Zn´ azornˇen´ı konzervativn´ıch vlastnost´ı gravitaˇcn´ı s´ıly.

F Obr´ azek 1.13: Pruˇzná ˇst´ıhlá raketa zat´ıˇzená tahem motoru.

1.7

Metoda tuh´ ych d´ılc˚ u

Metoda tuh´ ych d´ılc˚ u, naz´ yvan´ a také metoda simplexn´ıch prvk˚ u, viz [8] a [2], uˇz´ıvá speciáln´ıho 11 zp˚ usobu diskretizace kontinua . Kontinuum je rozdˇeleno na tuhé d´ılce vzájemnˇe spojené prvky schopné deformace, tzv. kumulátory pˇretvoˇren´ı, viz obrázek 1.14b. Narozd´ıl od reálné konstrukce a metody koneˇcn´ ych prvk˚ u nez˚ ustává model hladce spojit´ y, viz obrázek 1.14a. V´ yhodou metody tuh´ ych d´ılc˚ u je bezproblémové vystihnut´ı geometrické nelinearity. V pˇr´ıpadˇe dynamické formulace nav´ıc obdrˇz´ıme pˇr´ımo soustavu pohybov´ ych rovnic v kanonickém tvaru. a)

b)

Obr´ azek 1.14: Idealizace konstrukce metodou tuh´ ych d´ılc˚ u. 11

Jako fyzik´ aln´ı diskretizace je myˇsleno vytvoˇren´ı modelu konstrukce, kter´ y je svou podstatou diskrétn´ı, lze jej ale ch´ apat jako re´ alnou konstrukci, viz [4].

34


Kaˇzd´ y tuh´ y d´ılec m´ a vlastn´ı pomˇernou hmotnost a vlastn´ı tˇeˇziˇstˇe, je tedy nositelem kinetické energie modelu. Deformace - potenciáln´ı energie modelu - se akumuluje na kontaktech tuh´ ych d´ılc˚ u v kumul´ atorech pˇretvoˇren´ı. Tyto jsou bˇeˇznˇe uvaˇzovány jako dokonale pruˇzné. Pruˇzné spojen´ı tuh´ ych d´ılc˚ u by mˇelo b´ yt voleno s ohledem na povahu u ´lohy tak, aby dostateˇcnˇe vystihovalo podstatu modelu, ale nenavyˇsovalo zbyteˇcnˇe sloˇzitost odvozen´ı a nebo v´ ypoˇcetn´ı n´ aroˇcnost ˇreˇsen´ı. Vybrané moˇzné varianty jsou ukázány na obrázku 1.15: • spojen´ı tuh´ ych d´ılc˚ u rotaˇcn´ı pruˇzinou (1.15a), • spojen´ı translaˇcn´ı pruˇzinou (1.15b), • spojen´ı kombinac´ı rotaˇcn´ıch a translaˇcn´ıch pruˇzin (1.15c), • spojen´ı obecnou pruˇzinou umoˇzn ˇuj´ıc´ı nav´ıc zahrnut´ı vlivu smyku (1.15d). a)

b)

c)

d)

Obr´ azek 1.15: Varianty kontaktn´ıch kumulátor˚ u pˇretvoˇren´ı.

35

Kapitola 2

Odvozen´ı numerick´ eho modelu 2.1

´ Uvod

Zámˇerem pr´ ace je vytvoˇren´ı nelineárn´ıho dynamického modelu volného prutu schopného velk´ ych deformac´ı. Pro tento u ´kol pouˇzijeme zvláˇstn´ı pˇr´ıpad metody tuh´ ych d´ılc˚ u, viz [8]. V´ yhodou tohoto modelu je, jak pozdˇeji uvid´ıme, jeho jednoduchost a nenároˇcnost na v´ ypoˇcetn´ı dobu. Ta je u numerick´ ych model˚ u pˇr´ımo u ´mˇerná nejvyˇsˇs´ı frekvenci, kterou model kmit´ a. U ˇst´ıhlého ocelového prutu jsou ˇrádovˇe nejvyˇsˇs´ı frekvence dosahovány pˇri normálovém kmit´ an´ı (kmit´ an´ı pruˇzn´ ych kumulátor˚ u normálového pˇretvoˇren´ı - translaˇcn´ıch pruˇzin). Na chov´ an´ı takové konstrukce má vˇsak normálové kmitán´ı na rozd´ıl od kmit´ an´ı pˇr´ıˇcného zanedbateln´ y vliv. Vol´ıme tedy model normálovˇe dokonale tuh´ y, ˇc´ımˇz se vyhneme vzniku vyˇsˇs´ıch frekvenc´ı. Jedin´ ymi kumulátory pˇretvoˇren´ı v modelu tak z˚ ustávaj´ı rotaˇcn´ı pruˇziny spojuj´ıc´ı tuhé d´ılce. Voln´ y prut uvaˇzujeme konstantn´ıho pr˚ uˇrezu s hmotnost´ı M ,

y

F

0

x

Obr´ azek 2.1: Model volného prutu zat´ıˇzeného sleduj´ıc´ı silou.

délkou L a ohybovou tuhost´ı EI. V´ ypoˇctov´ y model rozdˇeluje prut na n tuh´ ych d´ılc˚ u vzájemnˇe spojen´ ych klouby. V m´ıstech kloub˚ u nahrazuj´ı ohybovou tuhost lineárn´ı rotaˇcn´ı pruˇziny, nebot’ d´ıky dostateˇcné ˇst´ıhlosti prutu uvaˇzujeme materiál jako lineárnˇe pruˇzn´ y. Lze ukázat, ˇze pro prvky modelu plat´ı: m=

M L EI , l= , k= , n n l

(2.1)

kde m je hmotnost tuhého d´ılce, l je délka tuhého d´ılce a k je tuhost rotaˇcn´ı pruˇziny.

36

´ho modelu Odvozen´ı numericke

2.1.1

Princip odvozen´ı

Pˇri odvozov´ an´ı se budeme ˇr´ıdit z´ asadami klasické mechaniky, viz [19]. Nejdˇr´ıve sestav´ıme Lagrangeovu funkci L, viz sekce 1.3.7. L = Ek − Ep ,

(2.2)

kde Ek je kinetick´ a energie modelu a Ep je potenciáln´ı energie modelu. Pohybové rovnice modelu potom obdrˇz´ıme uˇzit´ım Eulerovy-Lagrangeovy rovnice, viz (1.73): ∂L d ∂L = , (2.3) dt ∂ q˙i ∂qi kde t je ˇcas, qi je zobecnˇen´ a souˇradnice systému a q˙i je zobecnˇená rychlost systému.

2.2 2.2.1

Energie Kinetick´ a energie modelu

Uvaˇzujme rovinn´ y model volného prutu sloˇzen´ y ze tˇr´ı tuh´ ych d´ılc˚ u, viz obrázek 2.2. D´ılce jsou stejnˇe dlouhé a rotaˇcn´ı pruˇziny maj´ı stejnou tuhost. Vztahy, které zde z´ıskáme, n´ am dále poslouˇz´ı pro zobecnˇen´ı celého problému. y

vy,3 vy,2 vy,1 0

vy,0

1

2

vx,3

φ2 ω2

φ1 ω1 vx,2

φ0 ω0 v x,1

vx,0 0

x Obr´ azek 2.2: Model volného prutu.

37


Jako prvn´ı formulujme v´ yraz pro kinetickou energii prvn´ıho d´ılce Ek1 . Pro vytnut´ y diferenciáln´ı element o délce dr plat´ı: 1 1 dEk0 = dm v12 = ρ dr v12 , (2.4) 2 2 kde dm je hmotnost diferenci´ aln´ıho elementu, ρ je hustota elementu a v1 je absolutn´ı rychlost elementu. y

l r

vy,0

φ0

dm

vx,0 0

x Obr´ azek 2.3: Prvn´ı d´ılec modelu.

Absolutn´ı rychlost elementu z´ıskáme sloˇzen´ım dvou nezávisl´ ych pohyb˚ u - translaˇcn´ıho pohybu prvn´ıho bodu o rychlosti v0 a rotaˇcn´ıho pohybu prvn´ıho d´ılce o u ´hlové rychlosti ω0 . D´ılˇc´ı sloˇzky translaˇcn´ı rychlosti vx1 a vy1 m˚ uˇzeme tedy zapsat ve tvaru: vx1 = vx0 − ω0 r sin(ϕ0 ),

(2.5)

vy1 = vy0 + ω0 r cos(ϕ0 ). Pro sloˇzky translaˇcn´ı rychlosti vx0 a vy0 plat´ı: d d x0 , vy0 = y0 , dt dt a pro u ´hlové rychlosti ωi obecnˇe:

(2.6)

vx0 =

d ϕi , i = 0, 1 . . . n − 1. dt S uˇzit´ım v´ yrazu: q 2 + v2 , v1 = vx1 y1 ωi =

(2.7)

(2.8)

m˚ uˇzeme po u ´pravˇe a dosazen´ı do rovnice (2.4) psát: 1 1 dEk0 = ρ dr v12 = ρ dr (vx0 − ω0 r sin(ϕ0 ))2 + (vy0 + ω0 r cos(ϕ0 ))2 . 2 2 Celkovou kinetickou energii prvn´ıho d´ılce Ek1 z´ıskáme integrac´ı po celé délce: Z l 1 Ek0 = ρ vx0 − ω0 r sin(ϕ0 ))2 + (vy0 + ω0 r cos(ϕ0 ))2 dr 2 0 1 1 2 2 2 2 = m vx0 + vy0 + ω0 l − vx0 ω0 l sin(ϕ0 ) + vy0 ω0 l cos(ϕ0 ) . 2 3

38

(2.9)

(2.10)


V pˇr´ıpadˇe druhého d´ılce postupujeme obdobnˇe. Celkovou rychlost elementu na druhém d´ılci obdrˇz´ıme sloˇzen´ım translaˇcn´ı rychlosti poˇcáteˇcn´ıho bodu v0 , obvodové rychlosti koncového bodu prvn´ıho d´ılce a rotaˇcn´ıho pohybu diferenciáln´ıho elementu na druhém d´ılci ou ´hlové rychlosti ω1 . Pro sloˇzky vx2 a vy2 z´ıskáváme vztahy: vx2 = vx0 − ω0 l sin(ϕ0 ) − ω1 r sin(ϕ1 ),

(2.11)

vy2 = vy0 + ω0 l cos(ϕ0 ) + ω1 r cos(ϕ1 ).

V´ yraz pro kinetickou energii diferenciáln´ıho elementu na druhém d´ılci má tedy následuj´ıc´ı podobu: 1 1 dEk1 = ρ dr v22 = ρ dr (vx0 − ω0 l sin(ϕ0 ) − ω1 r sin(ϕ1 ))2 + 2 2 +(vy0 + ω0 l cos(ϕ0 ) + ω1 r cos(ϕ1 ))2 ,

(2.12)

Pro z´ısk´ an´ı kinetické energie druhého d´ılce Ek1 opˇet zintegrujme v´ yraz (2.12) pˇres celou délku d´ılce: Z l 1 1 1 2 2 2 v2 dr = m vx0 + vy0 + ω02 l2 + ω1 2l2 − 2 vx0 ω0 l sin(ϕ0 ) + Ek1 = ρ 2 0 2 3 + 2 vy0 ω0 l cos(ϕ0 ) − vx0 ω1 l sin(ϕ1 ) + vy0 ω1 l cos(ϕ1 ) + (2.13) i + ω0 ω1 l2 cos(ϕ1 − ϕ0 ) . Oproti v´ yrazu (2.10) budou nadále pˇrib´ yvat souˇciny u ´hlov´ ych rychlost´ı, které v dalˇs´ım odvozen´ı pohybov´ ych rovnic zp˚ usob´ı nár˚ ust sloˇzitosti v´ ypoˇctu. S uˇzit´ım stejn´ ych pˇredpoklad˚ u jako v pˇredchoz´ım pˇr´ıpadˇe m˚ uˇzeme dále odvodit v´ yraz pro kinetickou energii tˇret´ıho d´ılce Ek2 . D´ılˇc´ı sloˇzky celkové rychlosti diferenciáln´ıho elementu na tˇret´ım d´ılci vzniknou obdobnˇe: vx3 = vx0 − ω0 l sin(ϕ0 ) − ω1 l sin(ϕ1 ) − ω2 r sin(ϕ2 ), vy3 = vy0 + ω0 l cos(ϕ0 ) + ω1 l cos(ϕ1 ) + ω2 r cos(ϕ2 ). Po integraci z´ısk´ av´ ame kinetickou energii tˇret´ıho d´ılce Ek2 ve tvaru: Z l 1 1 1 2 2 2 Ek2 = ρ v3 dr = m vx0 + vy0 + ω02 l2 + ω12 l2 + ω22 l2 − 2 vx0 ω0 l sin(ϕ0 ) + 2 0 2 3 + 2 vy0 ω0 l cos(ϕ0 ) − 2 vx0 ω1 l sin(ϕ1 ) + 2 vy0 ω1 l cos(ϕ1 ) − − vx0 ω2 l sin(ϕ2 ) + vy0 ω2 l cos(ϕ2 ) + 2 ω0 ω1 l2 cos(ϕ1 − ϕ0 ) + i + ω0 ω2 l2 cos(ϕ2 − ϕ0 ) + ω1 ω2 l2 cos(ϕ2 − ϕ1 ) .

39

(2.14)

(2.15)


Nyn´ı formulujme obecn´ y tvar kinetické energie s-tého tuhého d´ılce Ek,s . Vˇsimneme-li si rozvoje v´ yraz˚ u (2.10), (2.13) a (2.15), m˚ uˇzeme psát: ( s−1 s−1 X X 1 1 22 2 2 2 2 Ek,s = m vx0 + vy0 + l ωi + ωs l − 2 vx0 l ωi sin(ϕi ) + 2 3 i=0

+ 2 vy0 l + 2 l2

s−1 X

ωi cos(ϕi ) i=0 "k−1 s−2 X X

+l

+ ωs l [−vx0 sin(ϕs ) + vy0 cos(ϕs )] +

ωk ωi cos(ϕi − ϕk ) +

k=0 2

i=0

s−1 X

i=1

s−1 X

(2.16)

# ωk ωi cos(ϕi − ϕk ) +

i=k+1

) ωk ωs cos(ϕs − ϕk ) .

k=0

V´ yraz pro celkovou kinetickou energii modelu Ek má tedy podobu: ( n−1 X 1 1 2 2 2 2 ωi + (n − i − 1) − Ek = m n vx0 + vy0 + l 2 3 i=0

− l vx0

+ l vy0

n−1 X

[ωi sin(ϕi ) (1 + 2(n − i − 1))] +

i=0 n−1 X

[ωi cos(ϕi ) (1 + 2(n − i − 1))] +

(2.17)

i=0

+

n−1 X

" 2 l2

j=0

+l

2

n−2 X

j−2 X k−1 X k=0

ωk ωi cos(ϕi − ϕk ) +

i=1

j−1 X

! ωk ωi cos(ϕi − ϕk ) +

i=k+1

#) ωk ωn−1 cos(ϕn−1 − ϕk )

.

k=0

2.2.2

Potenci´ aln´ı energie modelu

Z´ıskán´ı v´ yrazu pro potenci´ aln´ı energii Ep bude snazˇs´ı, nebot’ v uvaˇzovaném modelu se potenciáln´ı energie akumuluje pouze v rotaˇcn´ıch pruˇzinách mezi tuh´ ymi d´ılci. Pruˇzinˇe spojuj´ıc´ı s-t´ y a s+1 d´ılec pˇridˇel´ıme index s, s+1. Napjatost pruˇziny m˚ uˇzeme zapsat jako: Ms,s+1 = k ∆ϕ = k (ϕs+1 − ϕs ) ,

(2.18)

pro jej´ı potenci´ aln´ı energii Ep,s tedy plat´ı: 1 1 Ep,s = Ms,s+1 ∆ϕ = k (ϕs+1 − ϕs )2 . 2 2

(2.19)

Celkovou potenci´ aln´ı energii Ep z´ıskáme sumac´ı vˇsech d´ılˇc´ıch v´ yraz˚ u: Ep =

X

Ep,i

n−2 1 X = k (ϕs+1 − ϕs )2 . 2

(2.20)

s=0

40


2.3

Pohybov´ e rovnice

Jak jiˇz bylo zm´ınˇeno, pro odvozen´ı pohybov´ ych rovnic pouˇzijeme základn´ı vztahy (2.2) a (2.3). Z´ aroveˇ n budeme m´ıt na pamˇeti substituce (2.6) a (2.7). Odvozen´ı tedy provedeme podle rovnic: d ∂L ∂L = , (2.21) dt ∂vx0 ∂x0 d dt

d dt

∂L ∂vy0 ∂L ∂ωs

=

=

∂L , ∂y0

(2.22)

∂L , s = 0, 1, 2, . . . , n − 1. ∂ϕs

(2.23)

Nejprve provedeme derivace kinetické energie Ek podle stavov´ ych promˇenn´ ych vx0 , vy0 a ωs . Podle v´ yrazu (2.17) m˚ uˇzeme ps´ at: ( ) n−1 X ∂Ek 1 = m 2 n vx0 − l [ωi sin(ϕi ) (1 + 2(n − i − 1))] , (2.24) ∂vx0 2 i=0

( ) n−1 X ∂Ek 1 = m 2 n vy0 + l [ωi cos(ϕi ) (1 + 2(n − i − 1))] , ∂vy0 2

(2.25)

i=0

1 ∂Ek 1 2 + (n − s − 1) − l vx0 sin(ϕs ) (1 + 2(n − s − 1)) + = m 2 l ωs ∂ωs 2 3 + l vy0 cos(ϕs ) (1 + 2(n − s − 1)) + 2 l

2

s−1 X

ωi cos(ϕs − ϕi )+

(2.26)

i=0 2

+ l ωn−1 cos(ϕs − ϕn−1 )

o

,

V´ ysledky, které jsme obdrˇzeli, d´ ale zderivujeme podle ˇcasu: ( n−1 X d ∂Ek 1 [ω˙i sin(ϕi ) (1 + 2(n − i − 1))] − = m 2 n v˙ x0 − l dt ∂vx0 2 i=0 ) n−1 X 2 −l ωi cos(ϕi ) (1 + 2(n − i − 1)) ,

(2.27)

i=0

d dt

∂Ek ∂vy0

( n−1 X 1 = m 2 n v˙ y0 + l [ω˙i cos(ϕi ) (1 + 2(n − i − 1))] − 2 i=0 ) n−1 X 2 −l ωi sin(ϕi ) (1 + 2(n − i − 1)) , i=0

41

(2.28)


d dt

∂Ek ∂ωs

1 1 2 = m 2 l ω˙ s + (n − s − 1) − l v˙ x0 sin(ϕs ) (1 + 2(n − s − 1)) − 2 3 − l vx0 ωs cos(ϕs ) (1 + 2(n − s − 1)) + l v˙ y0 cos(ϕs ) (1 + 2(n − s − 1)) − − l vy0 ωs sin(ϕs ) (1 + 2(n − s − 1)) + 2 l

2

s−1 X

ω˙ i cos(ϕs − ϕi ) −

i=0

− 2 l2

s−1 X

ωi ωs sin(ϕs − ϕi ) + ωi2 sin(ϕs − ϕi ) + l2 ω˙ n−1 cos(ϕs − ϕi )−

i=0 2

2 − l ωs ωn−1 sin(ϕn−1 − ϕs ) + l2 ωn−1 sin(ϕn−1 − ϕs )

o

. (2.29)

Nyn´ı si pˇriprav´ıme derivace kinetické energie Ek podle promˇenn´ ych x0 , y0 a ϕs . Z podm´ınek rovnováhy lze odvodit, ˇze: ∂Ek ∂Ek = 0, = 0. ∂x0 ∂y0

(2.30)

Derivace podle u ´hl˚ u pootoˇcen´ı d´ılc˚ u ϕs má potom tvar: ∂Ek 1 n = m − l vx0 ωs cos(ϕs ) (1 + 2(n − s − 1)) − l vy0 ωs sin(ϕs ) (1 + 2(n − s − 1)) − ∂ϕs 2 s−1 n−2 X X − 2 l2 ωi ωs sin(ϕs − ϕi ) + 2 l2 ωs ωi sin(ϕi − ϕs ) + i=0

i=s+1

+ l2 ωs ωn−1 sin(ϕn−1 − ϕs ) , (2.31)

∂Ek 1 n = m − l vx0 ωn−1 cos(ϕn−1 ) − l vy0 ωn−1 sin(ϕn−1 ) − ∂ϕn−1 2 ) n−2 X − l2 ωi ωn−1 sin(ϕn−1 − ϕi ) . i=0

42

(2.32)


Posledn´ım krokem je z´ıskat derivace potenciáln´ı energie Ek podle u ´hl˚ u pootoˇcen´ı d´ılc˚ u ϕs . Vyjdeme z v´ yrazu (2.20) a z´ısk´ ame v´ ysledek: ∂Ep = k (ϕ0 − ϕ1 ), ∂ϕ0

(2.33)

∂Ep = k (−ϕs−1 + 2ϕs − ϕs+1 ), ∂ϕs

(2.34)

∂Ep = k (ϕn−1 − ϕn−2 ), ∂ϕn−1

(2.35)

∂Ep = 0, ∂x0

(2.36)

∂Ep = 0. ∂x0

(2.37)

V tuto chv´ıli je jiˇz vˇse potˇrebné pˇripraveno a podle v´ yraz˚ u (2.21), (2.22) a (2.23) m˚ uˇzeme vytvoˇrit soustavu rovnic: d ∂L ∂Ek ∂Ep − , (2.38) = dt ∂vx0 ∂x0 ∂x0 d dt

d dt

∂L ∂vy0 ∂L ∂ωs

=

=

∂Ek ∂Ep − , ∂y0 ∂y0

(2.39)

∂Ek ∂Ep − , s = 0, 1, 2, . . . , n − 1. ∂ϕs ∂ϕs

43

(2.40)


2.3.1

Maticov´ y tvar

Pro pˇrehlednost z´ apisu pˇrevedeme soustavu rovnic do maticového tvaru: M

d (v + ω) = Qω 2 − Kϕ, dt d ϕ = ω, dt d s = v. dt

(2.41)

kde v, ω a ϕ jsou vektory stavov´ ych promˇenn´ ych systému - translaˇcn´ıch rychlost´ı, u ´hlov´ ych rychlost´ı a u ´hl˚ u jednotliv´ ych d´ılc˚ u. Matice M je symetrick´ a a reprezentuje hmotné momenty setrvaˇcnosti d´ılc˚ u. M˚ uˇzeme si vˇsimnout, ˇze pˇri mal´ ych hodnotách vzájemného pootoˇcen´ı d´ılc˚ u (mal´ ych deformac´ıch prutu) z˚ ustane pˇribliˇznˇe konstantn´ı, k ˇcemuˇz by doˇslo v pˇr´ıpadˇe linearizace ˇreˇsen´ı (momenty setrvaˇcnosti d´ılc˚ u nejsou z´ avislé na jejich pootoˇcen´ı), viz [2]. Rovnˇeˇz je zde zˇrejmé rozdˇelen´ı matice, kdy se prvn´ı dva ˇrádky vztahuj´ı k translaˇcn´ım rychlostem a zbyl´ y obsah matice se t´ yk´ a rychlost´ı u ´hlov´ ych. Matice Q obsahuje transformaˇcn´ı momenty setrvaˇcnosti. Nen´ı jiˇz ˇctvercová, avˇsak m˚ uˇzeme vidˇet, ˇze spodn´ı ˇctvercov´ a ˇc´ ast je antisymetrická s nulami na diagonále. Na rozd´ıl od matice M se zde hodnoty pˇri mal´ ych deformac´ıch prutu bl´ıˇz´ı nule. V´ yznam matice Q tedy nar˚ ustá s deformac´ı prutu a uplatn´ı se pro nelineárn´ı ˇreˇsen´ı. Maticov´ y tvar pohybov´ ych rovnic uzav´ırá pásová matice ohybov´ ych tuhost´ı rotaˇcn´ıch pruˇzin K. Pro dodrˇzen´ı form´ alnˇe správného tvaru matice nesm´ıme zapomenout doplnit dva prázdné ˇr´ adky (derivace podle vx0 a vy0 ). Matice tuhosti má potom tvar:   0 0 0 ... 0 0      0 0 0 . . . 0 0        1 −1           −1 2 −1   .  (2.42) K = k    −1 2 −1       .. .. ..   . . .           −1 2 −1     −1 1

44

 −l sin ϕ0 −l sin ϕ1 2n 0 (1 + 2(n − 1)) (1 + 2(n − 2))     l cos ϕ0 l cos ϕ1  2n  (1 + 2(n − 1)) (1 + 2(n − 2))     1  2l2 ( + (n − 1)) 2l2 cos(ϕ1 − ϕ0 )  3   1  1  M = m 2l2 ( + (n − 2)) 2  3          symetrie      

45 .. .

...

...

...

...

−l sin ϕn−1



     l cos ϕi  l cos ϕn−1  (1 + 2(n − i − 1))     2 2 2l cos(ϕi − ϕ0 ) l cos(ϕn−1 − ϕ0 )      2l2 cos(ϕi − ϕ1 ) l2 cos(ϕn−1 − ϕ1 ) ,     .. ..  . .     1 2 2 2l ( + (n − i − 1)) l cos(ϕn−1 − ϕi )   3      2 2 l 3

−l sin ϕi (1 + 2(n − i − 1))

(2.43)


l cos ϕ0 l cos ϕ1 l cos ϕ2  (1 + 2(n − 1)) (1 + 2(n − 2)) (1 + 2(n − 3))     l sin ϕ0 l sin ϕ1 l sin ϕ2   (1 + 2(n − 1)) (1 + 2(n − 2)) (1 + 2(n − 3))     l2 sin(ϕ1 − ϕ0 ) l2 sin(ϕ2 − ϕ0 )  0  (1 + 2(n − 2)) (1 + 2(n − 3))     2 l2 sin(ϕ2 − ϕ1 )  l sin(ϕ0 − ϕ1 ) 0 1  (1 + 2(n − 2)) (1 + 2(n − 3)) Q = m 2    2  l sin(ϕ0 − ϕ2 ) l2 sin(ϕ1 − ϕ2 )  0  (1 + 2(n − 3)) (1 + 2(n − 3))    .. .. ..   . . .    2  l sin(ϕ0 − ϕi ) l2 sin(ϕ1 − ϕi ) l2 sin(ϕ2 − ϕi )  (1 + 2(n − i − 1)) (1 + 2(n − i − 1)) (1 + 2(n − i − 1))   l2 sin(ϕ0 − ϕn−1 ) l2 sin(ϕ1 − ϕn−1 ) l2 sin(ϕ2 − ϕn−1 )



l2 sin(ϕi − ϕ2 ) (1 + 2(n − i − 1)) .. . 0 l2 sin(ϕi − ϕn−1 )

.

..

... ...

...

...

l2 sin(ϕi − ϕ0 ) (1 + 2(n − i − 1))

...

l2 sin(ϕi − ϕ1 ) (1 + 2(n − i − 1))

l sin ϕi (1 + 2(n − i − 1))

...

...

l cos ϕi (1 + 2(n − i − 1))

      l sin ϕn−1      2 l sin(ϕn−1 − ϕ0 )       2 l sin(ϕn−1 − ϕ1 ) .     2 l sin(ϕn−1 − ϕ2 )     ..   .     2 l sin(ϕn−1 − ϕi )     0

l cos ϕn−1



(2.44)


46


2.4

Aplikovan´ y model voln´ eho prutu

V pˇredchoz´ı ˇc´ asti jsme se zab´ yvali odvozen´ım konzervativn´ıho modelu volného prutu. Pro naˇse potˇreby ale model jeˇstˇe postr´ ad´ a mnoˇzstv´ı potˇrebn´ ych souˇcást´ı. V prvé ˇradˇe je tˇreba doplnit do modelu sleduj´ıc´ı zat´ıˇzen´ı, jenˇz je nezbytnou sloˇzkou u ´lohy. Dále je zapotˇreb´ı postihnout r˚ uzné varianty disipace energie, nebot’ varianta konzervativn´ıho modelu nen´ı pro simulaci reálného chov´ an´ı v´ ystiˇzn´ a. Posledn´ı nezbytnou souˇcást´ı je zahrnut´ı vlivu gravitaˇcn´ıho pole.

2.4.1

Sleduj´ıc´ı s´ıla

Sleduj´ıc´ı s´ıla je hlavn´ı souˇc´ ast´ı modelu, nebot’ nahrazuje tah reaktivn´ıho motoru v reálné u ´loze. Sleduj´ıc´ı s´ıla mˇen´ı sv˚ uj smˇer v závislosti na deformaci konstrukce, viz [10]. V naˇsem pˇr´ıpadˇe tak, ˇze p˚ usob´ı pod konstantn´ım u ´hlem v˚ uˇci konci volného prutu. Pˇridán´ı u ´ˇcink˚ u sleduj´ıc´ı s´ıly tedy d´ av´ a pohybov´ ym rovnic´ım (2.41) tvar: M

d (v + ω) = Qω 2 − Kϕ + F, dt

(2.45)

kde F je vektor u ´ˇcink˚ u sleduj´ıc´ı s´ıly. Sleduj´ıc´ı s´ılu oznaˇcme jako Ff ollower a u ´hel sevˇren´ y s posledn´ım d´ılcem modelu potom jako υ, viz obr´ azek 2.4. D´ılˇc´ı sloˇzky sleduj´ıc´ı s´ıly Ff ollower,x a Ff ollower,y potom m˚ uˇzeme vyjádˇrit jako: Ff ollower,x = Ff ollower cos(ϕn−1 + υ),

(2.46)

Ff ollower,y = Ff ollower sin(ϕn−1 + υ).

y

Ffollower

υ

φ

n-1

0

x Obr´ azek 2.4: P˚ usoben´ı sleduj´ıc´ı s´ıly.

47


Vektor vlivu sleduj´ıc´ı s´ıly F vznikne sloˇzen´ım translaˇcn´ıho p˚ usoben´ı sleduj´ıc´ı s´ıly a moment˚ u sleduj´ıc´ı s´ıly v˚ uˇci jednotliv´ ym kloub˚ um:       Fx 0 0             F      0 0  y                       0  sin(ϕ0 )   cos(ϕ0 )                    0 sin(ϕ ) cos(ϕ )      . 1 1 F =   +Ff ollower,x l  − F l (2.47) f ollower,y           .      .. ..  .      . .  .                        0 sin(ϕn−2 ) cos(ϕn−2 )             0 sin(ϕn−1 ) cos(ϕn−1 )

2.4.2

Materi´ alov´ e tlumen´ı

Je zˇrejmé, ˇze v´ ystiˇzn´ a formulace materiálového tlumen´ı vyˇzaduje experimentáln´ı ovˇeˇren´ı pˇri vhodn´ ych podm´ınk´ ach. Pˇredloˇzen´ a definice bude proto pozdˇeji ovˇeˇrena experimentem. Pˇri formulaci materi´ alového tlumen´ı vycház´ıme z nejjednoduˇsˇs´ı definice lineárn´ıho viskózn´ıho tlumen´ı: M = −c m ω,

(2.48)

kde c je koeficient u ´tlumu, m je hmotnost d´ılce a ω je vzájemná relativn´ı rychlost sousedn´ıch d´ılc˚ u. Pro popisovan´ y model tedy plat´ı: M = −cint m

d ∆ϕ = −cint m ωr , dt

(2.49)

kde cint je koeficient materi´ alového u ´tlumu a ωr je rychlost rozev´ırán´ı nebo sv´ırán´ı dvou sousedn´ıch d´ılc˚ u. Pro tlum´ıc´ı moment mezi s-t´ ym a s+1 d´ılcem tedy zˇrejmˇe plat´ı: M = −cint m (ωs+1 − ωs ).

(2.50)

Matici materi´ alového u ´tlumu oznaˇcme Di . Pohybové rovnice (2.41) tedy z´ıskaj´ı následuj´ıc´ı podobu: M

d (v + ω) = Qω 2 − Kϕ − Di ω, dt

(2.51)

kde ω je vektor u ´hlov´ ych rychlost´ı d´ılc˚ u.

48


Matice materi´ alového tlumen´ı Di tedy nab´ yvá tvaru:   0 0 0 ... 0 0     0  0 0 . . . 0 0        1 −1          −1 2 −1      . Di = cint m     −1 2 −1       .. .. ..   . . .           −1 2 −1     −1 1

2.4.3

(2.52)

Tlumen´ı vlivem okoln´ıho prostˇ red´ı

Do chov´ an´ı modelu chceme zahrnout i disipaci energie vlivem okoln´ıho prostˇred´ı. Vektor tlumen´ı vlivem okoln´ıho prostˇred´ı oznaˇc´ıme De . Upravme tedy pohybové rovnice (2.41) do tvaru: M

d (v + ω) = Qω 2 − Kϕ − De . dt

(2.53)

Nejprve definujme tlum´ıc´ı s´ılu obecnˇe jako: D = c m v,

(2.54)

kde c je koeficient u ´tlumu, m je hmotnost d´ılce a v je rychlost d´ılce. Abychom postihli závislost velikosti u ´tlumu na smˇeru pohybu, vztah (2.54) zap´ıˇseme jako: D = cext ai m v,

(2.55)

kde cext je koeficient u ´tlumu vlivem okoln´ıho prostˇred´ı a ai je koeficient zohledˇ nuj´ıc´ı vliv rozd´ılného natoˇcen´ı d´ılc˚ u v˚ uˇci smˇeru jejich pohybu, pro kter´ y plat´ı: ϕ + ϕ i−1 i , ai = sin ϕv,i − 2

(2.56)

kde ϕv,i je u ´hel vektoru absolutn´ı rychlosti d´ılce, kter´ y z´ıskáme podle vztahu: vy,i ϕv,i = arctan . vx,i

49

(2.57)


Pˇrehlednˇe je odvozen´ı zn´ azornˇeno na obrázku 2.5.:

y

φ φ i-1 i

vi

2

φ

v,i

φ

i

φ

i-1

0 Obr´ azek 2.5: Odvozen´ı závislosti u ´tlumu na smˇeru pohybu modelu.

Pˇri pohledu na v´ yraz (2.56) si vˇsimnˇeme, ˇze v pˇr´ıpadˇe malého rozd´ılu pr˚ umˇerného u ´hlu natoˇcen´ı sousedn´ıch d´ılc˚ u a smˇeru jejich pohybu bude hodnota koeficientu ai velmi mal´ a. Naopak s rostouc´ım rozd´ılem tˇechto u ´hl˚ u koeficient ai nab´ yvá maximáln´ıch hodnot a tlumen´ı je tak nejv´ yraznˇejˇs´ı. D´ıky této formulaci se intenzita tlumen´ı vlivem okoln´ıho prostˇred´ı bude v kaˇzdém m´ıstˇe modelu adekv´ atnˇe liˇsit. Do v´ yrazu pro tlum´ıc´ı s´ılu nav´ıc vloˇz´ıme ˇclen si , kter´ ym zavedeme vliv stabilizace rakety pomoc´ı kˇridélek. V m´ıstech modelu, kde bude poˇzadováno toto pˇr´ıdavné tlumen´ı, je hodnota parametru si vˇetˇs´ı neˇz jedna. Tˇren´ı o pláˇst’ rakety bude zanedbáno. Obecnˇe tedy plat´ı: D = cext si ai m v,

(2.58)

Pro tlum´ıc´ı s´ıly p˚ usob´ıc´ı v kloubech modelu nakonec m˚ uˇzeme psát: m vx,0 , 2 m vy,0 , Dy,0 = cext s0 a0 2 Dx,i = cext si ai m vx,i , i = 1, 2, . . . , n − 1,

Dx,0 = cext s0 a0

Dy,i = cext si ai m vy,i , i = 1, 2, . . . , n − 1, m Dx,n = cext sn an vx,n , 2 m Dy,n = cext sn an vy,n . 2

50

(2.59)


Vektor tlumen´ı vlivem okoln´ıho prostˇred´ı De se skládá z translaˇcn´ıch u ´ˇcink˚ u v´ yslednice tlum´ıc´ıch sil a moment˚ u tlum´ıc´ıch sil k jednotliv´ ym kloub˚ um modelu. V´ ysledn´ y tvar je tedy následuj´ıc´ı: Pn

i=0 Dx,i







0



0



      P      n       0 0  i=0 Dy,i                  Pn Pn       0 D sin(ϕ ) D cos(ϕ )     0 0 i=1 x,i  i=1 y,i              P P       n n 0 sin(ϕ ) D cos(ϕ ) D      1 1 De =  i=2 x,i  + l  i=2 y,i  −l .             .. .. ..       . . .                   Pn Pn         sin(ϕk−1 ) i=k Dx,i  cos(ϕk−1 ) i=k Dy,i  0             0 sin(ϕn−1 )Dx,n cos(ϕn−1 )Dy,n

51

(2.60)


2.4.4

Gravitaˇ cn´ı zrychlen´ı

Posledn´ı zb´ yvaj´ıc´ı ˇc´ ast aplikovaného modelu je vliv p˚ usoben´ı gravitaˇcn´ıho pole. Pˇridejme tedy do pohybov´ ych rovnic (2.41) vektor gravitaˇcn´ıho zrychlen´ı A: M

d (v + ω) = Qω 2 − Kϕ + A. dt

(2.61)

Velikost t´ıhové s´ıly Fg p˚ usob´ıc´ı v tˇeˇziˇsti tuhého d´ılce odvod´ıme snadno z pˇredpokladu: Fg = m g,

(2.62)

kde m je hmotnost d´ılce a g je intenzita gravitaˇcn´ıho pole. Podobnˇe pro setrvaˇcnou s´ılu Fa p˚ usob´ıc´ı v tˇeˇziˇsti tuhého d´ılce plat´ı: Fa = m ag ,

(2.63)

kde m je setrvaˇcn´ a hmotnost d´ılce a ag je zrychlen´ı p˚ usob´ıc´ı na d´ılec. Vektor gravitaˇcn´ıho zrychlen´ı A je tedy sloˇzen z translaˇcn´ıch u ´ˇcink˚ u setrvaˇcn´ ych sil a z moment˚ u setrvaˇcn´ ych sil ke kloub˚ um modelu:     0 0          M a   0 g                1  0   (n − 2 ) cos(ϕ0 )              1 (n − − 1) cos(ϕ ) 0 .    1 (2.64) A= − a m l 2 g          .   ..   .   .   .               1  0   (n − 2 − i) cos(ϕi )          1 0 2 cos(ϕn−1 )

2.4.5

Aplikovan´ y model

V´ ysledn´ a podoba aplikovaného modelu vznikne sloˇzen´ım vˇsech zm´ınˇen´ ych d´ılˇc´ıch sloˇzek, viz (2.45), (2.51), (2.53) a (2.61). Koneˇcn´ y tvar aplikovaného modelu je tedy následuj´ıc´ı: M

d (v + ω) = Qω 2 − Kϕ − Di − De ω + F + A, dt d ϕ = ω, dt d s = v. dt

52

(2.65)

Kapitola 3

V´ ysledky numerick´ ych simulac´ı 3.1

Vlastn´ı frekvence prutu

Pro zjiˇstˇen´ı vlastn´ıch frekvenc´ı konstrukce je zapotˇreb´ı vhodnˇe zvolit bud´ıc´ı impuls. V pˇr´ıpadˇe volného prutu je obt´ıˇzné prut rozkmitat se souˇcasn´ ym zachován´ım rovnováhy sil. Mohou tedy vzniknout neˇz´ adouc´ı posuny nebo pootoˇcen´ı, které komplikuj´ı sledován´ı vibrac´ı. Voln´ y prut m˚ uˇzeme vybudit symetricky, v takovém pˇr´ıpadˇe ovˇsem obdrˇz´ıme pouze vlastn´ı frekvence se symetrick´ ymi amplitudami. Z toho d˚ uvodu byl model upraven také do podoby, kdy je poˇc´ atek prutu pevnˇe vetknut a konstrukce tedy p˚ usob´ı jako konzola. Simulace byly provedeny pro ocelov´ y prut o délce L = 2 m, pr˚ uˇrezové ploˇse A = 0.001 m2 , momentu setrvaˇcnosti I = 8.33×10−9 m4 , modulu pruˇznosti E = 210×109 Pa a hmotnosti m = 15 700 g. Prut byl rozdˇelen na 40 tuh´ ych d´ılc˚ u. Postupy analytického ˇreˇsen´ı, jakoˇz i koˇreny frekvenˇcn´ıch rovnic, byly pˇrevzaty, viz [17]. Krok v´ ypoˇctu byl nastaven na 8×10−4 s, simulovan´ y ˇcas byl 50 s. V´ ysledky frekvenˇcn´ı anal´ yzy m˚ uˇzeme vidˇet v tabulce 3.1. Frekvence f1 , voln´ y prut f3 , voln´ y prut f5 , voln´ y prut f1 , konzola f2 , konzola f3 , konzola f4 , konzola f5 , konzola

Analytické ˇreˇsen´ı [Hz] 13.2876 71.8056 177.3140 2.1966 13.7662 38.5457 75.5343 124.8640

Simulace [Hz] 13.3896 72.7841 176.1263 2.1352 13.7785 38.4750 75.6227 125.2371

Odchylka [%] +0.772 +1.363 −0.670 −2.685 +0.211 −0.059 +0.241 +0.299

Tabulka 3.1: Srovn´ an´ı vypoˇcten´ ych v´ ysledk˚ u s analytick´ ym ˇreˇsen´ım vlastn´ıch frekvenc´ı.

53

´ sledky numericky ´ ch simulac´ı Vy

Tabulka 3.2 potom ukazuje konvergenci hodnot k prvn´ı vlastn´ı frekvenci volného prutu a ˇcasovou n´ aroˇcnost simulace. V´ ysledky jsou znázornˇeny na obrázku 3.1. Poˇcet d´ılc˚ u Frekvence [Hz] 10 12.6271 15 13.1993 20 13.3815 25 13.3893 30 13.3895 35 13.3896 40 13.3896

V´ ypoˇcetn´ı doba [s] 3.624 6.325 9.928 16.904 25.137 36.510 50.581

Tabulka 3.2: Konvergence hodnoty prvn´ı vlastn´ı frekvence a odpov´ıdaj´ıc´ı v´ ypoˇcetn´ı doba.

13.5

128,82

Frekvence [Hz]

13.3 13.1 12.9 12.7 12.5 10

15

20 25 30 Počet dílců

35

40

Obr´ azek 3.1: Konvergence hodnoty prvn´ı vlastn´ı frekvence volného prutu.

Výpočetní doba [s]

50 40 30 20 10 0 10

15


35

40

Obr´ azek 3.2: Závislost v´ ypoˇcetn´ı doby na poˇctu d´ılc˚ u.

54


3.2

Konzola zat´ıˇ zen´ a vlastn´ı t´ıhou

Dalˇs´ım parametrem, kter´ y byl ovˇeˇren, je pr˚ uhyb a pootoˇcen´ı volného konce konzoly zat´ıˇzené vlastn´ı t´ıhou q, viz obr´ azek 3.3. Toto ovˇeˇren´ı je zároveˇ n ovˇeˇren´ım správnosti u ´ˇcink˚ u gravitaˇcn´ıho pole. V´ ypoˇcet byl proveden pro ocelov´ y prut o délce L = 5 m, ohybové tuhosti EI = 100N m2 a hmotnosti m = 50 g. Konzolového uloˇzen´ı prutu je dosaˇzeno fixac´ı polohy prvn´ıho kloubu modelu a pootoˇcen´ı prvn´ıho d´ılce. Z toho d˚ uvodu je skuteˇcná délka konzoly L = (5 − 5/n), kde n je poˇcet d´ılc˚ u. Simulace prob´ıhala s modelem o 50 d´ılc´ıch. Vzorce pro analytické ˇreˇsen´ı byly pˇrevzaty, viz [11]. Pr˚ uhyb volného konce konzoly δ lze analyticky vyjádˇrit: q l4 δ= = 8EI

0,05×9.81 5−5/50

× (5 − 5/50)4 = 0.0721 m.

8 × 100

(3.1)

Pro pootoˇcen´ı volného konce konzoly Θ plat´ı vztah: q l3 = Θ= 6EI

0,05×9.81 5−5/50

× (5 − 5/50)3 = 0.0196 rad

6 × 100

(3.2)

q

δ Θ

l

Obr´ azek 3.3: Konzola zat´ıˇzená vlastn´ı t´ıhou.

Srovn´ an´ı v´ ysledk˚ u z numerického v´ ypoˇctu s analytick´ ym ˇreˇsen´ım nab´ız´ı tabulka 3.3.

Pr˚ uhyb δ [m] Pootoˇcen´ı Θ [rad]

Analytické ˇreˇsen´ı Simulace 0.0721 0.0735 0.0196 0.0198

Odchylka [%] +1.904 +1.020

Tabulka 3.3: Porovn´ an´ı v´ ysledk˚ u pr˚ uhybu a pootoˇcen´ı voln´ıho konce konzoly.

55


3.3

Kritick´ a s´ıla

Ztráta stability konstrukce a hodnota kritické s´ıly jsou velmi vhodné nástroje pro ovˇeˇren´ı správnosti chov´ an´ı modelu. Jev ztr´ aty stability u numerického modelu potvrzuje kvalitativnˇe stejné neline´ arn´ı chov´ an´ı jako u re´ alné konstrukce. Protoˇze ke ztrátˇe stability docház´ı náhle, je samotn´ a hodnota kritické s´ıly v´ yborn´ ym kvantitativn´ım ovˇeˇren´ım chován´ı modelu, viz [3].

3.3.1

Voln´ y prut

Simulace pro z´ısk´ an´ı hodnoty kritické s´ıly volného prutu byly provádˇeny pro ocelov´ y prut o délce L = 5 m a ohybové tuhosti EI = 100 Nm2 . Simulovan´ y ˇcas byl 200 sekund. Z´ıskané hodnoty kritické s´ıly byly d´ ale porovnány s v´ ysledky ze simulac´ı proveden´ ych pro stejn´ y model v aplikaci FyDiK, viz [6]. V´ ysledky z numerick´ ych simulac´ı byly dále porovnány s analytick´ ym ˇreˇsen´ım, viz [9]: Fcr ≈ 109.54

EI 100 = 109.54 2 = 438.16 N. L2 5

(3.3)

V´ ysledky numerick´ ych simulac´ı speciáln´ıho modelu - obdrˇzené hodnoty kritické s´ıly, jejich odchylky oproti analytickému ˇreˇsen´ı (3.3) a odpov´ıdaj´ıc´ı v´ ypoˇcetn´ı doba - jsou pˇrehlednˇe zobrazeny v tabulce 3.4. Poˇcet d´ılc˚ u Kritick´ a s´ıla [N] 5 326.12 10 411.83 15 426.81 20 432.05 25 434.49 30 435.80 35 436.59 40 437.14 45 437.48 50 437.91

Odchylka [%] −25.571 −6.009 −2.590 −1.394 −0.838 −0.539 −0.358 −0.233 −0.155 −0.057

V´ ypoˇcetn´ı doba [s] 1.255 2.397 4.326 7.344 12.042 18.783 27.486 38.478 51.829 71.481

Tabulka 3.4: V´ ysledky simulac´ı speciáln´ıho modelu volného prutu.

V´ ysledky obdrˇzené pˇri simulac´ıch v aplikaci FyDiK ukazuje podobnˇe strukturovaná tabulka 3.5. Poˇcet d´ılc˚ u Kritick´ a s´ıla [N] 10 400.73 20 433.27 30 439.52 40 441.73 50 442.75

Odchylka [%] −8.543 −1.116 +0.310 +0.815 +1.048

V´ ypoˇcetn´ı doba [s] 27.666 84.761 216.790 444.110 825.040

Tabulka 3.5: V´ ysledky simulac´ı proveden´ ych v aplikaci FyDiK.

56


Grafické zn´ azornˇen´ı konvergence kritické s´ıly nab´ız´ı obrázek 3.4, nár˚ ust v´ ypoˇcetn´ı doby v závislosti na poˇctu d´ılc˚ u modelu potom obrázek 3.5.

Kritická síla [N]

440.0 430.0 420.0 410.0 400.0 390.0 10 15 20 25 30 35 40 45 50 Počet dílců Obr´ azek 3.4: Konvergence hodnoty kritické s´ıly volného prutu.

Výpočetní doba [s]

80.0 60.0 40.0 20.0 0.0 10 15 20 25 30 35 40 45 50 Počet dílců Obr´ azek 3.5: N´ ar˚ ust v´ ypoˇcetn´ı doby s poˇctem stupˇ n˚ u volnosti.

57


3.3.2

Konzola

Hodnota kritické s´ıly byla ovˇeˇrena i na konzolovém prutu. Model z˚ ustává totoˇzn´ y, tedy oce2 lov´ y prut o délce L = 5 m a ohybové tuhosti EI = 100 Nm . V´ ysledky numerické simulace budou porovn´ any s analytick´ ym ˇreˇsen´ım, viz [1] a [3]: Fcr ≈ 20.05

EI . L2

(3.4)

Jak bylo zm´ınˇeno, délka vlastn´ı konzoly je vˇzdy L = (5−5/n). Kritická s´ıla tedy nav´ıc závis´ı na poˇctu d´ılc˚ u a liˇs´ı se tak pro kaˇzd´ y pˇr´ıpad. V´ ysledky numerick´ ych simulac´ı ve srovnán´ı s analytick´ ym ˇreˇsen´ım m˚ uˇzeme vidˇet v tabulce 3.6. Kritick´ a s´ıla, simulace [N] 93.95 88.23 85.85 84.51 83.75 83.21 82.98 82.83 82.70 82.65

Poˇcet d´ılc˚ u 5 10 15 20 25 30 35 40 45 50

Kritická s´ıla, analytické ˇreˇsen´ı [N] 125.31 99.01 92.06 88.87 87.02 85.83 84.99 84.37 83.89 83.51

Odchylka [%] −33.382 −12.221 −7.241 −5.152 −3.908 −3.144 −2.419 −1.845 −1.435 −1.037

Tabulka 3.6: V´ ysledky simulac´ı speciáln´ıho modelu konzoly.

Grafické zn´ azornˇen´ı konvergence v´ ysledk˚ u nab´ız´ı obrázek 3.6 94.0

Kritická síla [N]

92.0 90.0 88.0 86.0 84.0 82.0 5

10

15

20


40

45

Obr´ azek 3.6: Konvergence hodnoty kritické s´ıly konzoly.

58

50

Kapitola 4

Implementace modelu Implementace numerického modelu byla provedena v programovac´ım jazyku Java s kompilátorem Java Development Kit, viz [21]. Implementace prob´ıhala v prostˇred´ı Eclipse IDE for Java Developers, viz [22]. Pro frekvenˇcn´ı anal´ yzu kmitán´ı prutu je pouˇzita Rychlá Fourierova transformace (FFT), viz [4].

4.1

Implementace metod ˇ reˇ sen´ı dynamick´ eho syst´ emu

Dále jsou uvedeny zdrojové k´ ody ukazuj´ıc´ı pr˚ ubˇeh v´ ypoˇctu jednoho kroku jednotliv´ ymi numerick´ ymi metodami. Pro zpˇr´ıstupnˇen´ı kódu ˇctenáˇri nejdˇr´ıve následuje struˇcn´ y popis jednotliv´ ych metod: • assembleMatrixM() - sestav´ı matici hmotn´ ych moment˚ u setrvaˇcnosti M. • assembleMatrixQ() - sestav´ı matici transformaˇcn´ıch moment˚ u setrvaˇcnosti Q. • assembleMatrixK() - sestav´ı matici tuhost´ı rotaˇcn´ıch pruˇzin K. • assembleMatrixD() - sestav´ı matici materiálového u ´tlumu Di . • createEquatuion() - sestav´ı zb´ yvaj´ıc´ı vektory v soustavˇe rovnic: De , F, A. • getSolution() - vyˇreˇs´ı soustavu rovnic Gaussovou eliminaˇcn´ı metodou, vr´ at´ı vektor ˇreˇsen´ı. • updateAngles() - aktualizuje u ´hly pootoˇcen´ı tuh´ ych d´ılc˚ u. • updateCoords() - aktualizuje souˇradnice kloub˚ u spojuj´ıc´ıch tuhé d´ılce. • updateOmegas() - aktualizuje u ´hlové rychlosti tuh´ ych d´ılc˚ u. • updateVelocities() - aktualizuje translaˇcn´ı rychlosti kloub˚ u spojuj´ıc´ıch tuhé d´ılce. • updateCG() - vypoˇc´ıt´ a polohu tˇeˇziˇstˇe modelu z aktuáln´ıch souˇradnic. • updateKineticEnergy() - vypoˇc´ıtá aktuáln´ı kinetickou energii modelu. • updatePotentialEnergy() - vypoˇc´ıtá aktuáln´ı potenciáln´ı energii rotaˇcn´ıch pruˇzin.

59

Implementace modelu

• updateOverallEnergy() - vrát´ı celkovou energii modelu. • updateTime() - pˇriˇcte ˇcasov´ y krok k pˇredchoz´ımu ˇcasu simulace.

4.1.1

Eulerova metoda

Z implementace Eulerovy metody je zˇrejmá jednoduchost numerického ˇreˇsen´ı nelineárn´ı dynamické u ´lohy. Jedin´ a n´ aroˇcn´ a ˇc´ ast je ˇreˇsen´ı soustavy rovnic Gaussovou eliminaˇcn´ı metodou. 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730

/* Euler method */ public void stepEuler () { assembleMatrixM () ; assembleMatrixQ () ; assembleMatrixK () ; assembleMatrixD () ; // createEquation () ; // getSolution ( leftSide , rightSide ) ; // updateAngles () ; updateCoords () ; // updateOmegas () ; updateVelocities () ; // updateCG () ; // updateKineticEnergy () ; up da te Pot en ti nal En er gy () ; updateOverallEnergy () ; // updateTime () ; }

4.1.2

Semi-implicitn´ı Eulerova metoda

Jedin´ ym rozd´ılem oproti klasické Eulerovˇe metodˇe je zámˇena poˇrad´ı metod updateAngles() a updateCoords() s metodami updateOmegas() a updateVelocities(), viz sekce 1.4.3. 747 748 749 750 751 752 753 754 755 756 757 758 759

/* Semi - implicit Euler method */ public void stepImprovedEuler () { assembleMatrixM () ; assembleMatrixQ () ; assembleMatrixK () ; assembleMatrixD () ; // createEquation () ; // getSolution ( leftSide , rightSide ) ; // updateOmegas () ;

60

Implementace modelu

760 761 762 763 764 765 766 767 768 769 770 771 772

updateVelocities () ; // updateAngles () ; updateCoords () ; // updateCG () ; // updateKineticEnergy () ; up da te Pot en ti nal En er gy () ; updateOverallEnergy () ; // updateTime () ; }

4.1.3

Runge-Kuttova metoda

Runge-Kuttova metoda je ˇctyˇrbodová numerická metoda. Proto je v rámci jednoho kroku potˇreba ˇctyˇrikr´ at sestavit a vyˇreˇsit soustavu rovnic. Je zˇrejmé, ˇze tato skuteˇcnost bude vést ke znaˇcné ˇcasové n´ aroˇcnosti ˇreˇsen´ı. 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834

/** Runge - Kutta method */ public void stepRungeKutta () { Vector RKsolutionK1 = new Vector () ; Vector RKsolutionK2 = new Vector () ; Vector RKsolutionK3 = new Vector () ; Vector RKsolutionK4 = new Vector () ; // /* Saving values xi ( t ) */ for ( int i =0; i < angles . length ; i ++) { anglesT [ i ]= angles [ i ]; omegasT [ i ]= omegas [ i ]; } for ( int i =0; i < coordinates . length ; i ++) { coordinatesT [ i ]= coordinates [ i ]; velocitiesT [ i ]= velocities [ i ]; } // /* k1 x ( t ) * / / / / / / / / / / / / / / / / / / / / / / / / / / / assembleMatrixM () ; assembleMatrixQ () ; assembleMatrixK () ; assembleMatrixC () ; // createEquation () ; // getSolution ( leftSide , rightSide ) ; // updateOmegas () ; updateVelocities () ; // updateAngles () ; updateCoords () ;

61

Implementace modelu

835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890

// RKsolutionK1 . allocateFor ( solution ) ; // for ( int row =0; row < solution . getRowCount () ; row ++) { RKsolutionK1 . addValue ( row , solution . getValue ( row ) ) ; } /* k2 ( xi ( t ) + k1i * h /2) * / / / / / / / / / / / / / / / / / / / / / / / / / / / /* input for k2 */ for ( int row =0; row < solution . getRowCount () ; row ++) { solution . setValue ( row , RKsolutionK1 . getValue ( row ) /2 d ) ; } // for ( int i =0; i < angles . length ; i ++) { angles [ i ]= anglesT [ i ]; omegas [ i ]= omegasT [ i ]; } for ( int i =0; i < coordinates . length ; i ++) { coordinates [ i ]= coordinatesT [ i ]; velocities [ i ]= velocitiesT [ i ]; } // /* solving for k2 */ updateOmegas () ; updateVelocities () ; updateAngles () ; updateCoords () ; // // assembleMatrixM () ; assembleMatrixQ () ; assembleMatrixK () ; assembleMatrixC () ; // createEquation () ; // getSolution ( leftSide , rightSide ) ; // RKsolutionK2 . allocateFor ( solution ) ; // for ( int row =0; row < solution . getRowCount () ; row ++) { RKsolutionK2 . addValue ( row , solution . getValue ( row ) ) ; } // // /* k3 ( xi ( t ) + k2i * h /2 * / / / / / / / / / / / / / / / / / / / / / / / / / / / /* input for k3 */ for ( int row =0; row < solution . getRowCount () ; row ++) { solution . setValue ( row , RKsolutionK2 . getValue ( row ) /2 d ) ; }

62

Implementace modelu

891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946

// for ( int i =0; i < angles . length ; i ++) { angles [ i ]= anglesT [ i ]; omegas [ i ]= omegasT [ i ]; } for ( int i =0; i < coordinates . length ; i ++) { coordinates [ i ]= coordinatesT [ i ]; velocities [ i ]= velocitiesT [ i ]; } // /* solving for k3 */ // updateOmegas () ; updateVelocities () ; updateAngles () ; updateCoords () ; // // assembleMatrixM () ; assembleMatrixQ () ; assembleMatrixK () ; assembleMatrixC () ; // createEquation () ; // getSolution ( leftSide , rightSide ) ; // RKsolutionK3 . allocateFor ( solution ) ; // for ( int row =0; row < solution . getRowCount () ; row ++) { RKsolutionK3 . addValue ( row , solution . getValue ( row ) ) ; } // // /* k4 ( xi ( t ) + k3i * h ) * / / / / / / / / / / / / / / / / / / / / / / / / / / / for ( int row =0; row < solution . getRowCount () ; row ++) { solution . setValue ( row , RKsolutionK3 . getValue ( row ) ) ; } // for ( int i =0; i < angles . length ; i ++) { angles [ i ]= anglesT [ i ]; omegas [ i ]= omegasT [ i ]; } for ( int i =0; i < coordinates . length ; i ++) { coordinates [ i ]= coordinatesT [ i ]; velocities [ i ]= velocitiesT [ i ]; } // /* solving for k4 */ updateOmegas () ;

63

Implementace modelu

947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002

updateVelocities () ; updateAngles () ; updateCoords () ; // // assembleMatrixM () ; assembleMatrixQ () ; assembleMatrixK () ; assembleMatrixC () ; // createEquation () ; // getSolution ( leftSide , rightSide ) ; // RKsolutionK4 . allocateFor ( solution ) ; // for ( int row =0; row < solution . getRowCount () ; row ++) { RKsolutionK4 . addValue ( row , solution . getValue ( row ) ) ; } // // /* final solution xi ( t + h ) = xi ( t ) + h /6( k1i +2 k2i +2 k3i + k4i ) */ for ( int row =0; row < solution . getRowCount () ; row ++) { solution . setValue ( row , ( RKsolutionK1 . getValue ( row ) + 2 d * RKsolutionK2 . getValue ( row ) + 2 d * RKsolutionK3 . getValue ( row ) + RKsolutionK4 . getValue ( row ) ) /6 d ) ; } // for ( int i =0; i < angles . length ; i ++) { angles [ i ]= anglesT [ i ]; omegas [ i ]= omegasT [ i ]; } for ( int i =0; i < coordinates . length ; i ++) { coordinates [ i ]= coordinatesT [ i ]; velocities [ i ]= velocitiesT [ i ]; } // updateOmegas () ; updateVelocities () ; // updateAngles () ; updateCoords () ; // updateCG () ; // updateKineticEnergy () ; up da te Pot en ti nal En er gy () ; updateOverallEnergy () ; // updateTime () ; }

64

Implementace modelu

4.2

Srovn´ an´ı metod ˇ reˇ sen´ı dynamick´ eho syst´ emu

Pro u ´plnost je na m´ıstˇe nab´ıdnout vzájemné porovnán´ı pouˇzit´ ych numerick´ ych metod pro ˇreˇsen´ı dynamického systému. Maxim´ aln´ı ˇ casov´ y krok Maximáln´ım ˇcasov´ ym krokem rozum´ıme takov´ y krok hmax , pˇri kterém je numerické ˇreˇsen´ı jeˇstˇe stabiln´ı. Porovn´ an´ı v´ ysledk˚ u nab´ız´ı tabulka 4.1 a obrázek 4.1. Poˇcet d´ılc˚ u 10 20 30 40 50

hmax , Euler [s] 0.000 2 0.000 06 0.000 02 0.000 01 0.000 002

hmax , S-I Euler [s] 0.02 0.005 0.002 0.001 0.000 6

hmax , Runge-Kutta [s] 0.025 0.005 0.002 0.001 0.000 7

Tabulka 4.1: Maxim´ aln´ı ˇcasov´ y krok v závislosti na poˇctu d´ılc˚ u a metodˇe.

1

10

20

40

50

Euler Runge-Kutta S-IgEuler

0.1 Maximálnígkrokg[s]

Početgdílců 30

0.01 0.001 0.0001 0.00001

0.000001 Obr´ azek 4.1: Maxim´ aln´ı ˇcasov´ y krok v závislosti na poˇctu d´ılc˚ u a metodˇe.

Z v´ ysledk˚ u je patrné, ˇze klasická Eulerova metoda je o dva aˇz tˇri ˇrády ménˇe stabiln´ı neˇz semi-implicitn´ı Eulerova metoda a metoda Runge-Kutta, pro které jsou v´ ysledky velmi podobné a z´ avis´ı na konkrétn´ıch podm´ınkách simulace.

65

Implementace modelu

ˇ Casov´ a n´ aroˇ cnost v´ ypoˇ ctu V´ ypoˇcetn´ı doba pˇri maxim´ aln´ım kroku hmax závis´ı také na v´ ypoˇcetn´ı nároˇcnosti jednotliv´ ych metod. Jak bude uk´ az´ ano, v´ ypoˇcetn´ı nároˇcnost znev´ yhodˇ nuje ˇctyˇrbodovou metodu RungeKutta. V´ ysledky nab´ız´ı tabulka 4.2 a dále ilustruje obrázek 4.2. Poˇcet d´ılc˚ u hmax , Euler [s] 10 18.506 20 181.768 30 1 385.120 40 16 253.156 50 151 288.201

hmax , S-I Euler [s] 0.694 2.739 14.253 63.996 208.206

hmax , Runge-Kutta [s] 1.328 8.126 31.485 111.889 263.767

Tabulka 4.2: V´ ypoˇcetn´ı doba v závislosti na poˇctu d´ılc˚ u a metodˇe.

100000

VýpočetníKdobaK[s]

10000 1000 100 10 Euler Runge-Kutta S-IKEuler

1 0.1 10

20

30 PočetKdílců

40

50

Obr´ azek 4.2: V´ ypoˇcetn´ı doba v závislosti na poˇctu d´ılc˚ u a metodˇe.

V´ ysledn´ a ˇcasov´ a n´ aroˇcnost v´ ypoˇctu m´ırnˇe favorizuje semi-implicitn´ı Eulerovu metodu. Nicménˇe je na m´ıstˇe pˇripomenout, ˇze oproti metodˇe Runge-Kutta má semi-implicitn´ı Eulerova metoda pouze empirick´ y p˚ uvod a nemus´ı tak b´ yt vˇzdy spolehlivá.

66

Implementace modelu

4.3

Grafick´ e prostˇ red´ı modelu

Pro zobrazov´ an´ı pr˚ ubˇehu simulace bylo dále vytvoˇreno grafické prostˇred´ı, které umoˇzn ˇuje sledovat chov´ an´ı modelu v pr˚ ubˇehu simulace. Zároveˇ n nab´ız´ı mnoˇzstv´ı dalˇs´ıch moˇznost´ı zobrazován´ı, jak bude uk´ az´ ano d´ ale. Z´ akladn´ı prostˇ red´ı Základn´ı grafické prostˇred´ı modelu je znázornˇeno na obrázku 4.3. Graficky zobrazuje: • vlastn´ı model - tuhé d´ılce spojené rotaˇcn´ımi pruˇzinami, • tˇeˇziˇstˇe modelu, • lok´ aln´ı souˇradn´ y systém - hlavn´ı osy setrvaˇcnosti, • p˚ usob´ıc´ı sleduj´ıc´ı s´ılu a zn´ azornˇen´ı gravitaˇcn´ıho pole. Textovˇe jsou zobrazov´ any: • fyzik´ aln´ı vlastnosti modelu - délka, hmotnost, ohybová tuhost a odpov´ıdaj´ıc´ı parametry tuh´ ych d´ılc˚ u a rotaˇcn´ıch pruˇzin, • parametry tlumen´ı a p˚ usob´ıc´ıho zat´ıˇzen´ı • souˇradnice tˇeˇziˇstˇe a smˇery hlavn´ıch os setrvaˇcnosti, • momenty setrvaˇcnosti a deviaˇcn´ı moment, • nastaven´ı simulace - metoda ˇreˇsen´ı, ˇcasov´ y krok, ˇcas simulace, v´ ypoˇcetn´ı doba a poˇcet krok˚ u.

Dalˇ s´ı moˇ znosti zobrazen´ı V pr˚ ubˇehu v´ ypoˇctu lze v re´ alném ˇcase sledovat libovoln´ y poˇcet libovoln´ ych stavov´ ych promˇenn´ ych tˇemito zp˚ usoby: • v z´ avislosti na ˇ case, viz obrázek 4.4, • ve vz´ ajemn´ e z´ avislosti na dalˇs´ı zvolené promˇenné, viz obrázek 4.5, • zobrazen´ı promˇ enn´ ych stejn´ eho typu, napˇr´ıklad u ´hlov´ ych rychlost´ı ω, viz obr´ azek 4.6.

67

Implementace modelu

Obr´ azek 4.3: Základn´ı grafické prostˇred´ı modelu.

68

Implementace modelu

Obr´ azek 4.4: Zobrazen´ı stavové promˇenné v závislosti na ˇcase.

Obr´ azek 4.5: Zobrazen´ı stavov´ ych promˇenn´ ych ve vzájemné závislosti.

Obr´ azek 4.6: Zobrazen´ı promˇenn´ ych stejného typu.

69

Kapitola 5

Z´ avˇ er V teoretické ˇc´ asti pr´ ace je ˇcten´ aˇri nab´ıdnut kompaktn´ı teoretick´ yu ´vod. Jsou popsány základn´ı zákony klasické mechaniky a dynamiky a je osvˇetlen energetick´ y - Lagrangeovsk´ y - pˇr´ıstup k popisu dynamického systému. V návaznosti Lagrangeovské odvozen´ı pohybov´ ych rovnic systému jsou pops´ any numerické metody jejich ˇreˇsen´ı. V u ´vodu práce byly dále zm´ınˇeny zdroje nelinearit v mechanickém systému a vlastnosti konzervativn´ıch a nekonzervativn´ıch sil. Jako posledn´ı je pops´ an princip metody tuh´ ych d´ılc˚ u, která je pouˇzita pro vytvoˇren´ı numerického modelu. V praktické ˇc´ asti pr´ ace je prezentován postup odvozen´ı nelineárn´ıho modelu volného prutu za pouˇzit´ı metody tuh´ ych d´ılc˚ u. Model je dále doplnˇen o nezbytné souˇcásti jako je zat´ıˇzen´ı sleduj´ıc´ı silou, vliv gravitaˇcn´ıho pole a materiálové tlumen´ı. Zvláˇstn´ı pozornost byla vˇenov´ ana podrobnému vystiˇzen´ı tlumen´ı vlivem okoln´ıho prostˇred´ı. S ohledem na re´ alnou pˇredlohu u ´lohy byly zavedeny vlivy rozd´ılného smˇeru natoˇcen´ı a pohybu d´ılc˚ u, a nakonec i moˇznost stabilizace rakety pomoc´ı kˇridélek. Tlumen´ı vlivem okoln´ıho prostˇred´ı m´ a v´ yznamn´ y vliv na chován´ı modelu na hranici ztráty stability a následné pokritické chov´ an´ı. Vzhledem k v´ yhodné formulaci modelu volného prutu byla provedena u ´prava pohybov´ ych rovnic a byl vytvoˇren samostatn´ y model konzoly. V dalˇs´ı ˇc´ asti pr´ ace jsou prezentovány v´ ysledky numerick´ ych simulac´ı model˚ u: • anal´ yza vlastn´ıch frekvenc´ı volného prutu a konzoly, • pr˚ uhyb a pootoˇcen´ı konzoly zat´ıˇzené vlastn´ı t´ıhou, • hodnota kritické sleduj´ıc´ı s´ıly volného prutu a konzoly. Pˇri simulac´ıch se potvrdil pˇredpoklad n´ızké ˇcasové nároˇcnosti v´ ypoˇctu. I simulace s velmi podrobn´ ymi modely lze prov´ adˇet mnohem rychleji neˇz v reálném ˇcase. Byla provedena vlastn´ı implementace numerického modelu v programovac´ım jazyce Java a zároveˇ n bylo vytvoˇreno grafické prostˇred´ı umoˇzn ˇuj´ıc´ı podrobnˇe sledovat pr˚ ubˇeh simulace.

70

´ ve ˇr Za

5.1

V´ yhledy

S pˇrihlédnut´ım k rozs´ ahl´ ym moˇznostem, které nab´ız´ı modely a jejich implementace, je plánov´ ano budouc´ı pouˇzit´ı vytvoˇren´ ych model˚ u pro: • kompletn´ı ovˇeˇren´ı vlastn´ıch frekvenc´ı volného prutu d´ıky lokáln´ı souˇradné soustavˇe tvoˇrené hlavn´ımi osami setrvaˇcnosti, • ovˇeˇren´ı vlastn´ıch tvar˚ u volného prutu a konzoly, • podrobné studium pokritického chován´ı konstrukc´ı, • studium vlivu tlumen´ı vlivem vnˇejˇs´ıho prostˇred´ı na pokritické chován´ı konstrukc´ı. Vzhledem k trvaj´ıc´ımu z´ ajmu autora o programován´ı se dále nab´ız´ı moˇznost paralelizace u ´lohy pouˇzit´ım programovac´ıho jazyka Java nebo nVidia CUDA.

71

Literatura [1] BECK, M. Die Knicklast des einseitig eingespannten, tangential gedr˝ uckten Stabes, Zeitschrift Angewandte Mathematik und Physik, Volume 3, Issue 3, 1952, p. 225-228. ´ [2] FRANTÍK, P. Neline´ arn´ı projevy mechanick´ ych konstrukc´ı, dizertaˇcn´ı práce, Ustav stavebn´ı mechaniky FAST VUT v Brnˇe, Brno 2004, 117 s. [3] FRANTÍK, P. Post-critical behaviour of beam loaded by follower force. CD proceedings ˇ a republika, 2009, 9 s. of national conference Engineering Mechanics 2009, Svratka, Cesk´ ISBN 80-86246-35-2 [4] FRANTÍK, P., MACUR, J. Kritická s´ıla imperfektovan´ ych systém˚ u. Modelov´ an´ı ˇ a republika, a mˇeˇren´ı neline´ arn´ıch jev˚ u v mechanice, sborn´ık pˇr´ıspˇevk˚ u, Neˇctiny, Cesk´ 2006, str. 67-70, ISBN 80-02-01827-3 [5] FRANTÍK, P. Implementace Rychlé Fourierovy transformace, Java bal´ık. Dostupné z: www.kitnarf.cz/java/. [6] FRANTÍK, P. Aplikace FyDik. Dostupné z: http://fydik.kitnarf.cz/ [7] GOLDSTEIN, H. Classical mechanics, 3 edition. Addison-Wesley, 2001, 680 s. ISBN 978-020-1657-029. ´ [8] HENRYCH, J. Upln´ a soustava finitn´ıch metod mechaniky a moˇznosti dalˇs´ıho rozvoje, ˇ studie 6.85 CSAV, Praha 1985, 165 s. [9] GURAN, A., OSSIA, K. On the stability of a flexible missile under an end thrust. Math. Comput. Modelling, Vol. 14, 1990, 5 s. [10] HERMANN, G. Dynamics and stability of mechanical systems with follower forces. Stanford University, Stanford, California, 1971, 234 s. ˇ AK, ´ J., KYTYR ´ J. Statika stavebn´ıch konstrukc´ı: základy stavebn´ı mechaniky, [11] KADLC staticky urˇcité prutové konstrukce. 3. vyd. Brno: VUTIUM, 2010, 349 s. ISBN 978-80214-3419-6. ˇ AK, ´ ´ J. Statika stavebn´ıch konstrukc´ı. Tˇret´ı dostisk druhého vyd. [12] KADLC J., KYTYR V Brnˇe: VUTIUM, 2009, 431 s. ISBN 978-80-214-3428-8. [13] LEECH, J. W. Klasick´ a mechanika (orig. Butler & Tanner Ltd., Frome, 1958), SNTL nakladatelstv´ı tech. lit., Praha 1970, 133 s.

72

Literatura

´ [14] MACUR, J. Uvod do teorie dynamick´ ych systém˚ u a jejich simulace. 1. vyd. Brno: VUT, 1995, 87 s. ISBN 80-214-0698-4. [15] MACUR, J. Dynamické systémy. [online]. [cit. http://www.fce.vutbr.cz/studium/materialy/Dynsys/

2014-05-03].

Dostupné

z:

[16] MACUR, J. Simulace dynamick´ ych systém˚ u v jazyce Java. [online]. [cit. 2014-05-03]. Dostupné z: http://www.fce.vutbr.cz/studium/materialy/simul/ ˇ [17] MELCER, J. Kuch´ arov´ a, D.: Dynamika stavebných konˇstrukci´ı, Zilinsk´ a univerzita, ˇ Zilina 2004, 240 s. [18] NEWTON, Isaac. The principia: mathematical principles of natural philosophy. S.l.: Snowball Pub, New York 2010. 991 s.ISBN 16-079-6240-3. [19] TAYLOR, John R. Classical mechanics. Sausalito, Calif.: University Science Books, c2005, xiv, 786 s. ISBN 18-913-8922-X. [20] TORBY, B. Energy Methods. Advanced Dynamics for Engineers, CBS College Publishing, New York 1984, 426 s. [21] http://www.oracle.com/technetwork/java/ [22] https://www.eclipse.org/downloads/

73

Literatura

Seznam symbol˚ u x, y, z r î, ˆj, k ˆ v t vx , vy , vz O ϕ ω xP , yP Fext a m p P vex Ep , T Ek , V k Fts Mrs Ep,ts Ep,rs g l qi q˙i L δr δW δs

souˇradné osy pravo´ uhlé souˇradné soustavy polohov´ y vektor bodu jednotkové vektory vektor rychlosti bodu ˇcas sloˇzky vektoru rychlosti poˇc´ atek souˇradné soustavy u ´hel od horizont´ aly u ´hlov´ a rychlost souˇradnice bodu P v pravo´ uhlé souˇradné soustavˇe vnˇejˇs´ı s´ıla p˚ usob´ıc´ı na hmotn´ y bod zrychlen´ı hmotnost hmotného bodu hybnost hmotného bodu celkov´ a hybnost soustavy rychlost trysk´ an´ı spalin potenci´ aln´ı energie kinetick´ a energie tuhost rotaˇcn´ı nebo translaˇcn´ı pruˇziny s´ıla p˚ usob´ıc´ı na translaˇcn´ı pruˇzinu moment p˚ usob´ıc´ı na rotaˇcn´ı pruˇzinu potenci´ aln´ı energie akumulovaná v translaˇcn´ı pruˇzinˇe potenci´ aln´ı energie akumulovaná v rotaˇcn´ı pruˇzinˇe gravitaˇcn´ı zrychlen´ı délka d´ılce, kyvadla zobecnˇen´ a souˇradnice systému zobecnˇen´ a rychlost systému Lagrangeova funkce virtu´ aln´ı posunut´ı virtu´ aln´ı pr´ ace virtu´ aln´ı posun

74

Literatura

δϕ R Z S qreal qvirt h k1 , k2 , k3 , k4 E ε σ W ∇ Φ L M I n Ms,s+1 Ff ollower c cint cext a D Fg ag f δ Θ Fcr

virtu´ aln´ı pootoˇcen´ı v´ yslednice reakˇcn´ıch sil setrvaˇcn´ a s´ıla akce systému uskuteˇcnˇen´ a trajektorie neuskuteˇcnˇen´ a trajektorie ˇcasov´ y krok koeficienty R-K metody modul pruˇznosti pomˇerné pˇrevoˇren´ı napˇet´ı pr´ ace gradient potenci´ al s´ıly délka prutu hmotnost prutu moment setrvaˇcnosti poˇcet d´ılc˚ u napjatost pruˇziny s, s + 1 sleduj´ıc´ı s´ıla koeficient tlumen´ı koeficient materiálového tlumen´ı koeficient tlumen´ı vlivem vnˇejˇs´ıho prostˇred´ı koeficient zohledˇ nuj´ıc´ı rozd´ılné natoˇcen´ı d´ılc˚ u a smˇeru jejich pohybu tlum´ıc´ı s´ıla gravitaˇcn´ı s´ıla zrychlen´ı p˚ usob´ıc´ı na d´ılec frekvence pr˚ uhyb volného konce konzoly pootoˇcen´ı volného konce konzoly kritick´ a s´ıla

Matice a vektory M Q K Di F v ω De A

matice hmotn´ ych moment˚ u setrvaˇcnosti matice transformaˇcn´ıch moment˚ u setrvaˇcnosti matice tuhost´ı rotaˇcn´ıch pruˇzin matice materi´ alového u ´tlumu vektor zat´ıˇzen´ı sleduj´ıc´ı silou vektor translaˇcn´ıch rychlost´ı vektor u ´hlov´ ych rychlost´ı vektor tlumen´ı vlivem vnˇejˇs´ıho prostˇred´ı vektor vlivu gravitaˇcn´ıho pole

75

Literatura

Seznam pˇ r´ıloh Pˇr´ıloha A: Implementace modelu v jazyce Java Pˇr´ıloha B: Kompaktn´ı disk

76

Literatura

Pˇ r´ılohy Pˇ r´ıloha A: Implementace modelu v jazyce Java 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

package cz . dl . freerodmodel ; import java . io . Fi leNo tFoun dExc eptio n ; import java . io . ObjectInputStream . GetField ; import java . io . PrintWriter ; import import import import

cz . kitnarf . fft . AmplitudeAnalyser ; cz . kitnarf . fft . FFT ; cz . kitnarf . fft . Peak ; cz . kitnarf . matrix .*;

public class FreeRodModel { /* Variables */ private int elementCount ; // private double mass ; private double length ; private double stiffness ; private double innerDampingCoef ; private double exterDampingCoef ; private double spikeDampingCoef ; private double tailDampingCoef ; private double mediumDampingCoef ; // private double elementMass = mass / elementCount ; private double elementLength = length / elementCount ; private double elementStiffness = stiffness / elementLength ; // private double step ; private double time ; // private double kineticEnergy ; private double potentialEnergy ; private double overallEnergy ; // private double followerForce ; private double followerAngle ; // private double gravityAcceleration ; /* Solution */ private Vector solution ;

77

Literatura

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

// RK private private private private

Vector Vector Vector Vector

RKsolutionK1 ; RKsolutionK2 ; RKsolutionK3 ; RKsolutionK4 ;

/* Arrays */ double [] angles double [] anglesSin double [] anglesCos double [] double [] double [] double [] double [] double [] double [] double []

= new double [ elementCount ]; = new double [ elementCount ]; = new double [ elementCount ];

omegas = new double [ elementCount ]; coordinates = new double [2* elementCount +2]; velocities = new double [2* elementCount +2]; centerOfGravity = new double [2]; dampingForcesX = new double [ elementCount +1]; dampingForcesY = new double [ elementCount +1]; aeroDampingCoef = new double [ elementCount +1]; velocityAngles = new double [ elementCount +1];

/* RK arrays */ double [] omegasT = new double [ elementCount ]; double [] anglesT = new double [ elementCount ]; double [] coordinatesT = new double [2* elementCount +2]; double [] velocitiesT = new double [2* elementCount +2]; /* Matrices */ int rows = elementCount +2; int columns = rows ; double [][] matrixM = new double [][] matrixQ = new double [][] matrixK = new double [][] matrixC = new

double double double double

[ rows ][ columns ]; [ rows ][ elementCount ]; [ rows ][ elementCount ]; [ rows ][ elementCount ];

/**/ Matrix leftSide = new Matrix () ; Vector rightSide = new Vector () ;

/** Initial state */ public void setInitialState ( double vx0Init , double vy0Init , double [] omegasInit , double [] anglesInit , double massInit , double lengthInit , double stiffnessInit , double innerDampingCoefInit , double followerForceInit , double gravityAccelerationInit , double followerAngleInit , double exterDampingCoefInit , double spikeDampingCoefInit , double tailDampingCoefInit , double me dium Dampi ngCo efIni t ) { velocities [0]= vx0Init ; velocities [1]= vy0Init ; // for ( int i =0; i < omegas . length ; i ++) omegas [ i ]= omegasInit [ i ]; // for ( int i =0; i < angles . length ; i ++) angles [ i ]= anglesInit [ i ];

78

Literatura

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

for ( int i =0; i < anglesSin . length ; i ++) anglesSin [ i ]= Math . sin ( angles [ i ]) ; for ( int i =0; i < anglesCos . length ; i ++) anglesCos [ i ]= Math . cos ( angles [ i ]) ; // mass = massInit ; length = lengthInit ; stiffness = stiffnessInit ; innerDampingCoef = innerD ampingCoefInit ; followerForce = followerForceInit ; followerAngle = followerAngleInit ; gravityAcceleration = gr a vi t yA cc e le ra t io n In it ; exterDampingCoef = exterD ampingCoefInit ; spikeDampingCoef = spikeD ampingCoefInit ; tailDampingCoef = tailDampingCoefInit ; mediumDampingCoef = me diumD ampi ngCoe fIni t ; } /** Constructor */ public FreeRodModel ( int elementCountInit , double stepInit ) { step = stepInit ; elementCount = elementCountInit ; } /** Calculates current model time */ private void updateTime () { time += step ; } /** Calculates current angular velocities */ private void updateOmegas () { for ( int i =0; i < omegas . length ; i ++) { omegas [ i ]+= step * solution . getValue ( i +2) ; } } /** Calculates current angles */ private void updateAngles () { for ( int i =0; i < angles . length ; i ++) { angles [ i ]+= step * omegas [ i ]; } for ( int i =0; i < anglesSin . length ; i ++) anglesSin [ i ]= Math . sin ( angles [ i ]) ; for ( int i =0; i < anglesCos . length ; i ++) anglesCos [ i ]= Math . cos ( angles [ i ]) ; } /** Calculates current coordinates of mass points */ private void updateCoords () { for ( int i =0; i < coordinates . length ; i ++) {

79

Literatura

153 154

if (i <2) coordinates [ i ]+= step * velocities [ i ]; /* x */ else if ( i %2==0) coordinates [ i ]= coordinates [i -2]+ elementLength * anglesCos [( i - i /2 -1) ]; /* y */ else if ( i %2!=0) coordinates [ i ]= coordinates [i -2]+ elementLength * anglesSin [( i - i /2 -2) ];

155 156 157 158 159 160 161 162 163 164 165 166 167

168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199

} } /** Calculates current velocities */ private void updateVelocities () { /* velocities [0] , velocities [1]*/ for ( int i =0; i <2; i ++) velocities [ i ]+= step * solution . getValue ( i ) ; /* other velocities */ for ( int i =2; i < velocities . length ; i ++) { /* x */ if ( i %2==0) velocities [ i ]= velocities [i -2] - omegas [( i - i /2 -1) ]* elementLength * anglesSin [( i - i /2 -1) ]; /* a % b == 0 A divisible by B if zero . */ /* y */ else if ( i %2!=0) velocities [ i ]= velocities [i -2]+ omegas [( i - i /2 -2) ]* elementLength * anglesCos [( i - i /2 -2) ]; } } /** Calculates current center of gravity coordinates */ private void updateCG () { /* Center of gravity coordinates */ /* x */ double topX = coordinates [0]+ coordinates [ coordinates . length -2]; for ( int i =2; i < coordinates . length -2; i +=2) topX +=2* coordinates [ i ]; double bottomX =2* elementCount ; centerOfGravity [0]= topX / bottomX ; /* y */ double topY = coordinates [1]+ coordinates [ coordinates . length -1]; for ( int i =3; i < coordinates . length -1; i +=2) topY +=2* coordinates [ i ]; double bottomY =2* elementCount ; centerOfGravity [1]= topY / bottomY ; } /** Calculates current kinetic energy */ private void updateKineticEnergy () { kineticEnergy =0; kineticEnergy += elementCount *( velocities [0]* velocities [0]+ velocities [1]* velocities [1]) ; for ( int i =0; i < elementCount ; i ++) kineticEnergy += omegas [ i ]* omegas [ i ]* elementLength * elementLength *(1/3 d +( elementCount -i -1) ) ; for ( int i =0; i < elementCount ; i ++) kineticEnergy -= velocities [0]* omegas [ i ]* elementLength * anglesSin [ i ]*(1+2*( elementCount -i -1) ) ; for ( int i =0; i < elementCount ; i ++) kineticEnergy += velocities [1]* omegas [ i ]* elementLength * anglesCos [ i ]*(1+2*( elementCount -i -1) ) ; for ( int i =0; i < elementCount ; i ++) { for ( int l =0; l <= i ; l ++) {

80

Literatura

200 201 202 203 204

for ( int k =0; k <= l ; k ++) { if (k < l ) { if ( l == i ) kineticEnergy += elementLength * elementLength * omegas [ k ]* omegas [ l ]* Math . cos ( angles [ l ] - angles [ k ]) ; else kineticEnergy +=2* elementLength * elementLength * omegas [ k ]* omegas [ l ]* Math . cos ( angles [ l ] - angles [ k ]) ; } }

205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252

} } kineticEnergy *=0.5* elementMass ; } /** Calculates current potential energy */ private void upd at eP ote nt in alE ne rg y () { potentialEnergy =0; for ( int i =0; i < elementCount -1; i ++) { potentialEnergy +=( angles [ i +1] - angles [ i ]) *( angles [ i +1] - angles [ i ]) ; } potentialEnergy *=0.5* elementStiffness ; } /** Calculates current overall energy */ private void updateOverallEnergy () { overallEnergy =0; overallEnergy = kineticEnergy + potentialEnergy ; }

/** Fills the M matrix */ private void assembleMatrixM () {

for ( int r =0; r < rows ; r ++) { for ( int c =0; c < columns ; c ++) { /* filling the first symmetric half */ /* velocities [0] row */ if ( r ==0) { if ( c ==0) matrixM [0][0] =2 d * elementCount ; if ( c ==1) matrixM [0][1] =0 d ; if (c >=2) matrixM [0][ c ] = - elementLength * anglesSin [c -2]*(1+2*( elementCount -( c -1) ) ) ; if ( c == columns -1) matrixM [0][ c ] = - elementLength * anglesSin [c -2]; }

81

Literatura

253 254 255 256 257 258 259

/* velocities [1] row */ if ( r ==1) { if ( c ==0) matrixM [1][0] =0 d ; if ( c ==1) matrixM [1][1] =2 d * elementCount ; if (c >=2) matrixM [1][ c ] = elementLength * anglesCos [c -2]*(1+2*( elementCount -( c -1) ) ) ; if ( c == columns -1) matrixM [1][ c ] = elementLength * anglesCos [c -2]; } /* generic omega rows */ if (r >=2 && r < rows -1) { if ( r == c ) matrixM [ r ][ r ] =2 d * elementLength * elementLength *(1/3 d +( elementCount +( r -(2* r -1) ) ) ) ; if (c > r && c < columns -1) matrixM [ r ][ c ] = elementLength * elementLength * Math . cos ( angles [c -2] - angles [r -2]) *(1+2*( elementCount -( c -1) ) ) ; if (c > r && c == columns -1) matrixM [ r ][ c ] = elementLength * elementLength * Math . cos ( angles [c -2] - angles [r -2]) ; } /* last omega row */ if ( r == rows -1 && c == r ) matrixM [ r ][ c ] =2/3 d * elementLength * elementLength ;

260 261 262 263 264 265 266

267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299

} } ////////////////////////////////////////// /* generating the other symmetric half */ for ( int r =0; r < rows ; r ++) { for ( int c =0; c < columns ; c ++) { if ( r != c ) matrixM [ c ][ r ]= matrixM [ r ][ c ]; } } } /** Fills the Q matrix */ private void assembleMatrixQ () { for ( int r =0; r < rows ; r ++) { for ( int c =0; c < rows -2; c ++) { if ( r ==0) matrixQ [ r ][ c ]= elementLength * anglesCos [ c ]*(1+2*( elementCount -( c +1) ) ) ; if ( r ==1) matrixQ [ r ][ c ]= elementLength * anglesSin [ c ]*(1+2*( elementCount -( c +1) ) ) ; if (r >=2 && r < rows -1 ) {

82

Literatura

300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354

if (c >= r -2) matrixQ [ r ][ c ]= elementLength * elementLength * Math . sin ( angles [ c ] - angles [r -2]) *(1+2*( elementCount -( c +1) ) ) ; } } } /* generating the other antisymmetric half */ for ( int r =2; r < rows ; r ++) { for ( int c =0; c < rows -2; c ++) { if (r -2!= c ) matrixQ [ c +2][ r -2]= - matrixQ [ r ][ c ]; } } } /** Fills the K matrix */ private void assembleMatrixK () { for ( int r =2; r < rows ; r ++) { for ( int c =0; c < elementCount ; c ++) { if ( r ==2) { if ( c ==0) matrixK [ r ][ c ]=1; if ( c ==1) matrixK [ r ][ c ]= -1; } if (r >2 && r < rows -1) { if ( c == r -3) matrixK [ r ][ c ]= -1; if ( c == r -2) matrixK [ r ][ c ]=2; if ( c == r -1) matrixK [ r ][ c ]= -1; } if ( r == rows -1) { if ( c == r -3) matrixK [ r ][ c ]= -1; if ( c == r -2) matrixK [ r ][ c ]=1; } } } } /** Fills the C matrix */ private void assembleMatrixC () { for ( int r =2; r < rows ; r ++) { for ( int c =0; c < elementCount ; c ++) { if ( r ==2) { if ( c ==0) matrixC [ r ][ c ]=1; if ( c ==1) matrixC [ r ][ c ]= -1; }

83

Literatura

355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409

if (r >2 && r < rows -1) { if ( c == r -3) matrixC [ r ][ c ]= -1; if ( c == r -2) matrixC [ r ][ c ]=2; if ( c == r -1) matrixC [ r ][ c ]= -1; } if ( r == rows -1) { if ( c == r -3) matrixC [ r ][ c ]= -1; if ( c == r -2) matrixC [ r ][ c ]=1; } } } }

/** Solves the equation , gets the solution vector */ private void createEquation () { /* Creating the matrices and vectors for solving */ /* Matrix M */ Matrix kitMatrixM = new Matrix () ; kitMatrixM . allocateFor ( matrixM ) ; kitMatrixM . setValues ( matrixM ) ; kitMatrixM . multiply (0.5* elementMass ) ; // leftSide = kitMatrixM ; // System . out . println ( kitMatrixM ) ; /* Matrix Q */ Matrix kitMatrixQ = new Matrix () ; kitMatrixQ . allocateFor ( matrixQ ) ; kitMatrixQ . setValues ( matrixQ ) ; kitMatrixQ . multiply (0.5* elementMass ) ; // System . out . println ( kitMatrixQ ) ; /* Matrix K */ Matrix kitMatrixK = new Matrix () ; kitMatrixK . allocateFor ( matrixK ) ; kitMatrixK . setValues ( matrixK ) ; kitMatrixK . multiply ( elementStiffness ) ; // System . out . println ( kitMatrixK ) ; /* Matrix C */ Matrix kitMatrixC = new Matrix () ; kitMatrixC . allocateFor ( matrixC ) ; kitMatrixC . setValues ( matrixC ) ; kitMatrixC . multiply ( innerDampingCoef * elementMass ) ; // System . out . println ( kitMatrixC ) ; /* Vector omega ^2 */ Vector kitOmegaSq = new Vector () ; kitOmegaSq . allocate ( elementCount , 1) ; for ( int r =0; r < elementCount ; r ++) kitOmegaSq . addValue (r , omegas [ r ]* omegas [ r ]) ;

84

Literatura

410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461

// System . out . println ( kitOmegaSq ) ; /* Vector of angles */ Vector kitFi = new Vector () ; kitFi . allocate ( elementCount , 1) ; for ( int r =0; r < elementCount ; r ++) // System . out . println ( kitFi ) ;

kitFi . addValue (r , angles [ r ]) ;

/* Vector of angular velocities */ Vector kitOmega = new Vector () ; kitOmega . allocate ( elementCount , 1) ; for ( int r =0; r < elementCount ; r ++) kitOmega . addValue (r , omegas [ r ]) ; // System . out . println ( kitOmega ) ; double followerForceX = followerForce * Math . cos ( angles [ elementCount -1]+ followerAngle ) ; double followerForceY = followerForce * Math . sin ( angles [ elementCount -1]+ followerAngle ) ;

/* Vector Fx * l * sin ( fi ) */ Vector kitFollowerX = new Vector () ; kitFollowerX . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if ( r ==0) kitFollowerX . addValue (r , - followerForceX ) ; if (r >1) kitFollowerX . addValue (r , followerForceX * elementLength * anglesSin [r -2]) ; }

/* Vector - Fy * l * cos ( fi ) */ Vector kitFollowerY = new Vector () ; kitFollowerY . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if ( r ==1) kitFollowerY . addValue (r , - followerForceY ) ; if (r >1) kitFollowerY . addValue (r , - followerForceY * elementLength * anglesCos [r -2]) ; } /* Vector of moments of gravity forces X */ Vector kitGravityXmoment = new Vector () ; kitGravityXmoment . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if (r >1) { double value =0; kitGravityXmoment . addValue (r , value ) ; } }

/* Vector of moments of gravity forces Y */

85

Literatura

462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488

489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505

Vector kitGravityYmoment = new Vector () ; kitGravityYmoment . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if (r >1) { kitGravityYmoment . addValue (r , - elementMass * gravityAcceleration * elementLength * anglesCos [r -2]*( elementCount -0.5 -( r -2) ) ) ; } } // System . out . println ( kitGravityYmoment ) ;

/* Vector of translational gravity forces */ Vector kitGravityTrans = new Vector () ; kitGravityTrans . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if ( r ==0) kitGravityTrans . addValue (r , 0) ; if ( r ==1) kitGravityTrans . addValue (r , - mass * gravityAcceleration ) ; } /* Vector of overall gravity effect */ Vector kitGravityOverall = new Vector () ; kitGravityOverall . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) kitGravityOverall . addValue (r , kitGravityXmoment . getValue ( r ) + kitGravityYmoment . getValue ( r ) + kitGravityTrans . getValue ( r ) ) ; // System . out . println ( kitGravityOverall ) ; /* Aero damping coefficients */ for ( int i =0; i <= elementCount ; i ++) { if ( velocities [2* i ] >0 && velocities [2* i +1] >0 ) velocityAngles [ i ]= Math . atan ( velocities [2* i +1]/ velocities [2* i ]) ; if ( velocities [2* i ] <0 && velocities [2* i +1] >0 ) velocityAngles [ i ]= Math . PI - Math . atan ( - velocities [2* i +1]/ velocities [2* i ]) ; if ( velocities [2* i ] <0 && velocities [2* i +1] <0 ) velocityAngles [ i ]= Math . PI + Math . atan ( - velocities [2* i +1]/ - velocities [2* i ]) ; if ( velocities [2* i ] >0 && velocities [2* i +1] <0 ) velocityAngles [ i ]= Math . atan ( velocities [2* i +1]/ - velocities [2* i ]) ; if ( velocities [2* i ]==0 && velocities [2* i +1] >0 ) velocityAngles [ i ]=0.5* Math . PI ; if ( velocities [2* i ]==0 && velocities [2* i +1] <0 ) velocityAngles [ i ]= -0.5* Math . PI ; if ( velocities [2* i ] >0 && velocities [2* i +1]==0 ) velocityAngles [ i ]=0; if ( velocities [2* i ] <0 && velocities [2* i +1]==0 ) velocityAngles [ i ]= Math . PI ; if ( i ==0) aeroDampingCoef [ i ]= Math . abs ( Math . sin ( Math . abs ( velocityAngles [ i ] angles [ i ]) ) ) ;

506

86

Literatura

507 508

if ( i == elementCount ) aeroDampingCoef [ i ]= Math . abs ( Math . sin ( Math . abs ( velocityAngles [ i ] angles [i -1]) ) ) ;

509 510 511 512 513 514 515 516 517

if ( i !=0 && i != elementCount ) aeroDampingCoef [ i ]= Math . abs ( Math . sin ( Math . abs ( velocityAngles [ i ] -( angles [i -1]+ angles [ i ]) /2 d ) ) ) ; /* The damping forces */ for ( int s =0; s <= elementCount ; s ++) { if ( s ==0) dampingForcesX [ s ]= - exterDampingCoef *( spikeDampingCoef * mediumDampingCoef * aeroDampingCoef [ s ]) *0.5* elementMass * velocities [ s *2]; if (s >=1 && s < elementCount ) dampingForcesX [ s ]= - exterDampingCoef *( mediumDampingCoef * aeroDampingCoef [ s ]) * elementMass * velocities [ s *2]; if ( s == elementCount ) dampingForcesX [ s ]= - exterDampingCoef *( tailDampingCoef * mediumDampingCoef * aeroDampingCoef [ s ]) *0.5* elementMass * velocities [ s *2];

518

519

520 521 522 523 524 525

526

527

528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548

} for ( int s =0; s <= elementCount ; s ++) { if ( s ==0) dampingForcesY [ s ]= - exterDampingCoef *( spikeDampingCoef * mediumDampingCoef * aeroDampingCoef [ s ]) *0.5* elementMass * velocities [ s *2+1]; if (s >=1 && s < elementCount ) dampingForcesY [ s ]= - exterDampingCoef *( mediumDampingCoef * aeroDampingCoef [ s ]) * elementMass * velocities [ s *2+1]; if ( s == elementCount ) dampingForcesY [ s ]= - exterDampingCoef *( tailDampingCoef * mediumDampingCoef * aeroDampingCoef [ s ]) *0.5* elementMass * velocities [ s *2+1]; } /* Vector of external translational damping forces */ Vector kitExtDampingTrans = new Vector () ; kitExtDampingTrans . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if ( r ==0) { double value =0; for ( int i =0; i <= elementCount ; i ++) value += dampingForcesX [ i ]; kitExtDampingTrans . addValue (r , value ) ; } if ( r ==1) { double value =0; for ( int i =0; i <= elementCount ; i ++) value += dampingForcesY [ i ]; kitExtDampingTrans . addValue (r , value ) ; }

87

Literatura

549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583

584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600

}

/* Vector of moments of external damping forces X */ Vector kitExtDampingXmoment = new Vector () ; kitExtDampingXmoment . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if (r >1) { double value =0; for ( int i =r -1; i <= elementCount ; i ++) value += dampingForcesX [ i ]; kitExtDampingXmoment . addValue (r , elementLength * anglesSin [r -2]* value ); } }

/* Vector of moments of external damping forces Y */ Vector kitExtDampingYmoment = new Vector () ; kitExtDampingYmoment . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) { if (r >1) { double value =0; for ( int i =r -1; i <= elementCount ; i ++) value += dampingForcesY [ i ]; kitExtDampingYmoment . addValue (r , elementLength * anglesCos [r -2]* value ); } }

/* Vector of overall external damping effect */ Vector kitExtDampingOverall = new Vector () ; kitExtDampingOverall . allocate ( elementCount +2 , 1) ; for ( int r =0; r < elementCount +2; r ++) kitExtDampingOverall . addValue (r , - kitExtDampingXmoment . getValue ( r ) + kitExtDampingYmoment . getValue ( r ) + kitExtDampingTrans . getValue ( r ) ) ;

//************************// /* Solving the equation */ /* Multiplying matrices and vectors */ Vector QtimesOmegaSq = Vector . multiply ( kitMatrixQ , kitOmegaSq ) ; Vector KtimesFi = Vector . multiply ( kitMatrixK , kitFi ) ; Vector CtimesOmega = Vector . multiply ( kitMatrixC , kitOmega ) ; /* Summation of right side */ Vector rightSide1 = new Vector () ; rightSide1 . allocate ( elementCount +2 , 1) ; for ( int r =0; r < rightSide1 . getRowCount () ; r ++) { // if (r <3) rightSide . setValue (r , 0) ;

88

Literatura

601

602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653

rightSide1 . setValue (r , QtimesOmegaSq . getValue ( r ) - KtimesFi . getValue ( r ) - CtimesOmega . getValue ( r ) + kitFollowerX . getValue ( r ) + kitFollowerY . getValue ( r ) + kitGravityOverall . getValue ( r ) + kitExtDampingOverall . getValue ( r ) ) ; } rightSide = rightSide1 ; } private void getSolution ( Matrix leftSide , Vector rightSide ) { /* Solution vector */ GaussElimination solver = new GaussElimination () ; solution = solver . solve ( leftSide , rightSide ) ; // System . out . println ( solution ) ; }

public double getCoordinate ( int index ) { return coordinates [ index ]; } public double getVelocity ( int index ) { return velocities [ index ]; } public double [] getCoordinates () { return coordinates ; } public double [] getAngles () { return angles ; }

public double [] g e t C e n t e r O f G r a v i t y C o o r d i n a t e s () { return centerOfGravity ; } public double [] getOmegas () { return omegas ; } public double getPotentialEnergy () { return potentialEnergy ; } public double [] getAeroCoefficients () { return aeroDampingCoef ; } public double [] getVelocityAngles () { return velocityAngles ; }

89

Literatura

654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681

public double [] getExternalDampingX () { return dampingForcesX ; } public double [] getExternalDampingY () { return dampingForcesY ; } public double [] getVelocitiesX () { double [] velocitiesX = new double [ elementCount +1]; for ( int i =0; i < velocitiesX . length ; i ++) { velocitiesX [ i ]= velocities [2* i ]; } return velocitiesX ; } public double [] getVelocitiesY () { double [] velocitiesY = new double [ elementCount +1]; for ( int i =0; i < velocitiesY . length ; i ++) { velocitiesY [ i ]= velocities [2* i +1]; } return velocitiesY ; }

Pˇ r´ıloha B: Kompaktn´ı disk Na pˇriloˇzeném kompaktn´ım disku se nacház´ı: • elektronick´ a verze bakal´ aˇrské práce, • kompletn´ı zdrojové k´ ody numerického modelu a grafického rozhran´ı.

90

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ BRNO UNIVERSITY OF TECHNOLOGY

Recommend Documents