Ond rej Jirouıek 13. œnora 2006mech.fd.cvut.cz/members/jirousek/download/k618y2m1/ymkp_fem.pdf ·...

Metoda konecných prvku–poznámky k prednáškám–

Ondrej Jiroušek

13. února 2006

2

Obsah

1 Úvod 71.1 Rozdelení teoretické mechaniky . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.2 Základní myšlenka MKP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.3 Velmi strucná historie MKP . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Opakování maticové algebry 112.1 Rovnost matic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2 Scítání a odecítání matic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3 Násobení skalárem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.4 Násobení matice maticí . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.5 Transpozice matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.6 Symetrické matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.7 Nulová matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.8 Jednotková matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.9 Diagonální matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.10Determinant matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.11Singulární matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.12Pozitivne definitní matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.13Inverzní matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.14Derivace a integrace matice . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Rešení soustav lineárních algebraických rovnic 153.1 Prímé metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.1.1 Gaussova eliminacní metoda . . . . . . . . . . . . . . . . . . . . . . . 15

3.1.2 Rešení pomocí LU rozkladu . . . . . . . . . . . . . . . . . . . . . . . . 16

3.1.3 GEM pomocí LDLT rozkladu . . . . . . . . . . . . . . . . . . . . . . . 18

3.2 Iteracní metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.2.1 Jacobiho iteracní metoda . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.2.2 Gaussova-Seidelova iteracní metoda . . . . . . . . . . . . . . . . . . 23

3

4 OBSAH

3.2.3 Jacobiho metoda s relaxací . . . . . . . . . . . . . . . . . . . . . . . . 24

3.2.4 Gaussova-Seidelova relaxacní metoda (SOR) . . . . . . . . . . . . . . 24

3.3 Stabilita a citlivost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.3.1 Stabilita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.3.2 Citlivost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4 Základní princip metody konecných prvku 27

5 Prímá metoda tuhosti(DIRECT STIFFNESS METHOD) 295.1 Globální soustava rovnic (MASTER STIFFNESS EQUATION) . . . . . . . . . 31

5.2 Urcení matice tuhosti prvku . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5.2.1 Poznámky k odvození matice tuhosti taženého-tlaceného prvku . . 33

5.3 Transformace souradnic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

5.4 Globalizace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

5.5 Sestavení globální matice tuhosti . . . . . . . . . . . . . . . . . . . . . . . . 35

5.6 Rešení (SOLUTION) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5.6.1 Aplikace okrajových podmínek redukcí . . . . . . . . . . . . . . . . . 37

5.6.2 Aplikace okrajových podmínek modifikací . . . . . . . . . . . . . . . 39

5.6.3 Aplikace okrajových podmínek v maticové podobe . . . . . . . . . . 39

5.7 Výpocet vnitrních sil (odvozených velicin)(POST PROCESSING) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

5.8 Implementace taženého-tlaceného prvku v Octave . . . . . . . . . . . . . . 40

5.8.1 Matice tuhosti elementu . . . . . . . . . . . . . . . . . . . . . . . . . 41

6 Variacní principy 456.1 Funkcionál celkové potenciální energie . . . . . . . . . . . . . . . . . . . . . 46

6.1.1 Vnitrní (deformacní) energie pri tahu tlaku . . . . . . . . . . . . . . . 46

6.1.2 Vnitrní (deformacní) energie pri ohybu . . . . . . . . . . . . . . . . . 47

6.1.3 Vnitrní (deformacní) energie pri kroucení . . . . . . . . . . . . . . . . 47

6.1.4 Vnejší energie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

6.1.5 Celková potenciální energie prutu . . . . . . . . . . . . . . . . . . . . 48

6.2 Lagrangeuv a Castiglianuv variacní princip . . . . . . . . . . . . . . . . . . 48

6.3 Ritzova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

6.4 Princip virtuálních prací . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

7 Tažený-tlacený prut 537.1 Variacní formulace I (x) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

7.2 Variacní formulace II (ξ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

OBSAH 5

8 Ohýbaný prut (PLANE BEAM ELEMENT) 578.1 Ohýbaný prut (Euler–Bernoulli) . . . . . . . . . . . . . . . . . . . . . . . . . 58

8.2 Ohýbaný prut s vlivem smyku (Timošenko) . . . . . . . . . . . . . . . . . . 66

8.3 Ohýbaný prut na pružném podloží (Timošenko–Winkler) . . . . . . . . . . 68

8.4 Ohýbaný prut – vliv normálové síly . . . . . . . . . . . . . . . . . . . . . . . 69

9 Izoparametrické elementy 719.1 Prirozené souradnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

9.1.1 Prirozené souradnice pro 1-D elementy . . . . . . . . . . . . . . . . . 72

9.1.2 Trojúhelníkové souradnice . . . . . . . . . . . . . . . . . . . . . . . . 73

9.1.3 Prirozené souradnice pro ctyrsten . . . . . . . . . . . . . . . . . . . . 77

9.2 Trojúhelníkový prvek s konstantní deformací (CST) . . . . . . . . . . . . . 78

9.3 Podstata izoparametrických prvku . . . . . . . . . . . . . . . . . . . . . . . 81

9.3.1 Souhrn bázových funkcí pro nejbežnejší rovinné prvky . . . . . . . 86

9.3.1.1 Constant Strain Triangle (CST, T3) . . . . . . . . . . . . . . 86

9.3.1.2 Linear Strain Element (LST, T6) . . . . . . . . . . . . . . . . 86

9.3.1.3 Linear Quadrilateral Element (Q4) . . . . . . . . . . . . . . . 86

9.3.1.4 Quadratic Quadrilateral Element (Q8) . . . . . . . . . . . . 88

10 Prvky pro rovinnou úlohu 8910.1Základní úlohy teorie pružnosti . . . . . . . . . . . . . . . . . . . . . . . . . 89

10.2Úlohy rovinné pružnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

10.2.1Rovinná deformace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

10.2.2Rovinná napjatost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

10.2.3Osová soumernost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

10.2.4Prvek pro rovinnou deformaci . . . . . . . . . . . . . . . . . . . . . . 96

10.2.5Prvek pro rovinnou napjatost . . . . . . . . . . . . . . . . . . . . . . . 98

10.2.6Statické podmínky rovnováhy (Cauchyho rovnice) . . . . . . . . . . 98

10.2.7Okrajové podmínky . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

10.2.8Odvození matice tuhosti pro CST-prvek a rovinnou úlohu . . . . . . 99

11 Deskové prvky 10111.1Predpoklady, deskové teorie . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

11.2Kirchhoffova teorie tenkých desek . . . . . . . . . . . . . . . . . . . . . . . . 101

11.2.1Kinematické rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

11.2.2Geometrické rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

11.2.3Statické rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

11.2.4Biharmonická rovnice desky . . . . . . . . . . . . . . . . . . . . . . . 104

6 OBSAH

11.2.5Funkcionál potenciální energie . . . . . . . . . . . . . . . . . . . . . . 105

11.2.6Hellinger-Reissneruv princip . . . . . . . . . . . . . . . . . . . . . . . 105

11.2.7Výpocet matice tuhosti (doplnit pozdeji) . . . . . . . . . . . . . . . . 105

11.3Mindlinova teorie tlustých desek . . . . . . . . . . . . . . . . . . . . . . . . 107

12 Úvod do dynamiky v MKP 10912.1Variacní princip v dynamice - Hamiltonuv princip . . . . . . . . . . . . . . 109

12.2Matice hmotnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

12.2.1Globalizace matice hmotnosti . . . . . . . . . . . . . . . . . . . . . . 113

12.3Vlastní kmitání s více stupni volnosti . . . . . . . . . . . . . . . . . . . . . . 114

12.3.1Metoda konstant tuhosti . . . . . . . . . . . . . . . . . . . . . . . . . 114

12.3.2Príklad (metoda konstant tuhosti) . . . . . . . . . . . . . . . . . . . . 117

12.4Rešení kmitání pomocí MKP . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

12.5Vlastnosti vlastních tvaru . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

12.5.1Ortogonalita vlastních tvaru . . . . . . . . . . . . . . . . . . . . . . . 119

12.6Tlumené kmitání . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

12.7Rešení vlastního kmitání . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

12.7.1Stodolova iteracní metoda . . . . . . . . . . . . . . . . . . . . . . . . 122

12.7.2Inverzní iterace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

12.7.3Metoda iterace podprostoru . . . . . . . . . . . . . . . . . . . . . . . . 123

13 Rešení úloh nelineární mechaniky 12513.1Prístupy k rešení nelineárních problému . . . . . . . . . . . . . . . . . . . . 125

13.2Inkrementání technika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

13.3Iteracní (Newtonovy) techniky . . . . . . . . . . . . . . . . . . . . . . . . . . 128

13.4Smíšené techniky (STEP-ITERATIVE OR MIXED) . . . . . . . . . . . . . . . . . 128

14 Numerická integrace 13114.1Rovnomerné delení (Newton-Cotesovo schema) . . . . . . . . . . . . . . . . 131

14.2Nerovnomerné delení (Gaussova integracní formule) . . . . . . . . . . . . . 132

Kapitola 1

Úvod

Tento materiál vznikl z prednášek pro posluchace kurzu “Metoda konecných prvku" naFakulte dopravní, CVUT v Praze. Cílem tohoto kurzu je seznámit studenty se základ-ními myšlenkami metody konecných prvku, predevším v oblasti teoretické pružnosti adynamiky konstrukcí. Protože je kurz urcen pro posluchace tretího až pátého rocníkubakalárského i magisterského studia, bylo nutno kurz doplnit i o nekolik kapitol, patrí-cích do jiných teoretických predmetu. Postupne tak do tohoto materiálu pribyly strucnékapitoly se základy maticové algebry, kapitola pojednávající o metodách rešení soustavlineárních algebraických rovnic i kapitola o numerické integraci. Tyto kapitoly, patrícído kurzu numerické analýzy, nejsou bežne prednášeny na FD CVUT.

Skripta postupne vznikala tak, aby na ctenáre kladla pouze minimální znalosti z teoriepružnosti, ci matematické analýzy. V prubehu semestru byly doplnovány partie, kterénebyly známé všem posluchacum kurzu. Postupne tak vznikly úvodní partie, které opa-kují nejduležitejší pojmy z teorie pružnosti, nebo dynamiky konstrukcí.

Metoda konecných prvku je ve popsána nejprve pomocí prímé metody tuhosti, kteráje uvedena pro nejjednodušší typ namáhání a nejjednodušší prvek, tedy pro prípadtahu (tlaku) a prutový prvek. Na tomto prípadu jsou uvedeny postupy používané vevšech ostatních prípadech namáhání, jako je transformace sourednic, globalizace avýpocet odvozených (sekundárních) velicin. Je zde také uvedena implementace taženého(tlaceného) prvku v Octave (Matlabu).

V následující kapitole jsou uvedeny základní myšlenky variacního poctu, uveden Lagran-geuv a Castiglianuv variacní princip a jednoduše popsána Ritzova metoda. Její použitíje zdokumentováno na príkladu nalezení rovnice ohybové cáry jednoduché staticky ur-cité konstrukce. Je zde rovnež uvedeno odvození matice tuhosti taženého (tlaceného)prutu pomocí principu virtuální práce (PVP).

Protože matice tuhosti taženého (tlaceného prvku) byla odvozena pomocí PVP již v ka-pitole variacního poctu, je další kapitola rovnou venována ohýbanému prvku. Je zdeprovedeno odvození bázových funkcí, matice tuhosti prvku je zde již odvozena pomocíenergetického principu.

Další kapitola se zabývá izoparametrickými elementy, jejich podstatou. Je zde prove-deno odvození trojúhelníkových souradnic, dále prirozených souradnic napríklad proctyrsten. Podrobneji jsou zde vysvetleny izoparametrické prvky. Aplikace bázových funkcíizoparametrických elementu pro odvození matice tuhosti trojúhelníkového prvku s kon-stantní deformací je uvedeno v záveru této kapitoly.

V následující kapitole jsou probrány základní úlohy teorie pružnosti a na príkladu rovin-ných prípadu pružnosti,tedy rovinné deformace a rovinné napjatosti odvozeny matice

7

8 KAPITOLA 1. ÚVOD

tuhosti elementu pro tyto dve základní úlohy matematické pružnosti. V záveru kapitolyje provedeno strucné shrnutí další rovinný elementu, elementy vyššího rádu.Kapitola 12 se venuje MKP v dynamice konstrukcí. Úvodem jsou zopakovány základnípojmy volného i vynuceného, netlumeného i tlumeného kmitání. Úvod do kmitání sou-stav s více stupni volnosti je vysvelen pomocí metody konstant tuhosti, je zde vysvetlenai Stodolova iteracní metoda. Je zde provedeno odvození matice hmotnosti elementu. Zá-kladní vztahy pro MKP jsou dále odvozeny pomocí Hamiltonova variacního principu.Metody pro rešení vlastního kmitání, jako jsou Rayleig-Ritzova metoda, metoda statickékondenzace, nebo iteracní metody jsou provedeny pouze strucne. Rešení vynucenéhokmitání je uvedeno pouze pomocí rozvoje do vlastních tvaru.Kapitola 13 pojednává strucne o metodách rešení nelineárních problému. Inkremen-tální i iteracní techniky jsou nastíneny pouze v hlavních myšlenkách, smíšené metodydokonce jen informativne.Poslední kapitola je venována numerické integraci. Je zde provedeno pomerne detailneodvození Newton-Cotesových vzorcu a Gaussovy integracní fomrule.

1.1 Rozdelení teoretické mechanikyMechanika jako soucást fyziky muže být rozdelena do trí cástí: teoretická mechanika,aplikovaná a výpoctová. My se samozrejme budeme zabývat tou poslední. Výpoctovoumechaniku lze ovšem i dále delit, napr na nanomechaniku, mikromechaniku, mezo-mechaniku, mechaniku pevných teles, kapalin... My se zameríme pouze na mechanikupevných teles.

Mechanika

teoretickaaplikovana

vypoctova

nano amicro

kontinua

pevnych teleskapalinmulti−

systemu

Mechanika je v základním kurzu delena na dve cásti, na statiku a dynamiku. Dyna-mika na rozdíl od statiky zkoumá vliv setrvacných sil. Vlastnosti zkoumaného systémumužou být i ve statice casove závislé, setrvacné síly se však zanedbávají. Statické pro-blémy mohou být skritne statické, nebo kvazistatické. Pri kvazistatické analýze je trebapocítat s casem, napr. v problémech jako je únavová analýza (cyklování materiálu), se-dání základu, plasticita závislá na rychlosti deformace. Cas v techto príklad však jetreba uvažovat pouze jako relativní parametr. Z oblasti dynamiky se budeme zabývatpouze vlastním i vynuceným kmitáním konstrukcí, odezvou konstrukce na harmonickébuzení a stanovováním vlastních frekvecí a módu kmitání.Dalším merítkem, jak mechaniku delit je na lineární a nelineární. Zdroju nelinearit jevíce, nelinearita muže být materiálová, nebo geometrická. Budeme se zabývat problémynelineární mechaniky, odlišnostmi rešení, nelineárními materiálovými modely i mate-matickými metodami rešení soustav rovnic.

1.2 Základní myšlenka MKPVe strucnosti, základní myšlenkou metody konecných prvku je reprezentace konstrukcepomocí jednoduchých cástí, tzv. elementu. Tyto elementy jsou navzájem spojené v tzv.

1.3. VELMI STRUCNÁ HISTORIE MKP 9

uzlových bodech. Jednoduché funkce jsou voleny k aproximaci hledané veliciny, v tomtotextu nejcasteji posunutí, po jednotlivých prvcích. Rešením konstrukce pomocí MKPjsou v tomto prípade hodnoty uzlových posunutí, tedy posunutí v konecném poctu uz-lových bodu. Jednoduché funkce jsou voleny tak, aby vhodným zpusobem aproximovalyposunutí po prvku, nejcasteji to bývají polynomy n-tého stupne, mohou to však být inapríklad trigonometrické funkce.

Variacní principy mechaniky, napr. princip minima potenciální energie, jsou používányk rešení podmínek rovnováhy na každém prvku. Napríklad pro prípad elastického ma-teriálu, potenciální energie je dána souctem vnitrní deformacní energie (práce vnitrníchsil na deformaci) a potenciální energie vnejších sil, tedy zatížení. Jestliže má být kon-strukce v rovnováze, tato potenciální energie musí nabývat minimální hodnoty.

Z tohoto strucného popisu plynou základní požadavky kladené na ctenáre. Ackolivvšechna témata jsou v tomto textu pokryta, vetšinou se jedná o pomerne strucný vý-klad a proto se doporucuje, aby ctenár byl seznámen alespon v minimálním rozsahu snásledujícími tématy:

• základy matematické teorie pružnosti

• variacní principy mechaniky

• maticová algebra

• rešení soustav lineárních algebraických rovnic

• numerická integrace

Algoritmizace nekterých základních postupu metody konecných prvku je v textu na-stínena pomocí nekterého z vyšších programovacích jazyku. Ackoliv stále ješte mnohovedeckého kódu je psáno ve Fortranu, soucasný trend je používat vyšší programovacíjazyky. V projektech s otevreným zdrojovým kódem je nejcasteji používaným jazykem(2005) jazyk C, následuje perl, C+, java, PHP. K výukovým úcelum je dokonce vhodnejšípoužívat jazyky jako Matlab, Maple, ci Mathematicu. V textu je užito open-sourcovéekvivalenty Matlabu, systému octave (http://www.octave.org); pro nastínení mož-nosti použití systému pro symbolickou algebru systém macsyma (http://macsyma.something.org).

1.3 Velmi strucná historie MKP

Obvyklá otázka zní: kdo vymyslel metodu konecných prvku? Na to není úplne jedno-duchá odpoved’. Casto se uvádí jako první R. Courrant, který poprvé užil myšlenku “pocástech spojitých polí” . Radeji se vyhneme diskuzi o “objeviteli” MKP a shrnme alesponnekteré nejvýznamejší milníky v ére MKP:

• 1943 Courrant (zavedení “po cástech spojitých polí”)

• 1956 Turner, Clough, Martin a Topp

• 1960 Clough (poprvé použil pojmu “konecný element”, rovinné problémy)

• 70-léta aplikace na sálových pocítacích

• 80-léta mikropocítace, PC, pre- a post- procesory, význam komercních MKP sys-tému roste

10 KAPITOLA 1. ÚVOD

• 90-léta analýzy velmi velkých strukturálních problému, aplikace v letectví, kosmo-nautice

M. J. Turner formuloval prímou metodu tuhosti (DSM, Direct Stiffness Method) v obdobímezi roky 1950 a 1962. Díky jemu zacala jako první v praxi používat DSM firma Boeing,zatímco ostatní letecké spolecnosti stále ješte používaly pro problémy aeroelasticity si-lovou metodu. Zavedení izoparametrických modelu, tvarových funkcí a tzv. patch testuje pricítáno B. M. Ironsovi. Dalším prukopníkem na poli MKP byl R. J. Melosh, kterýpopsal spojitost s Ritzovou metodou a systematizoval variacní zpusob odvození maticetuhosti elementu. Konecne E. L. Wilson vyvinul první “open source” konecnoprvkovýsystém.

Jako skutecný zacátek éry MKP je považováno publikování clánku “Stiffness and de-flection analysis of complex atructures, J. Aero Sci., 1956” v roce 1956 (Turner, Clough,Martin a Topp) i když pojmu “konecný element” užil poprvé až R. W. Clough v roce 1960v clánku, kde použil MKP pro rešení rovinného problému.

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

Kapitola 2

Opakování maticové algebry

Soustavu lineárních algebraických rovnic:

a11x1 + a12x2 + . . .+ a1nxn = b1

a21x1 + a22x2 + . . .+ a2nxn = b2 (2.1)...

an1x1 + an2x2 + . . .+ annxn = bn

s neznámými xi budeme zapisovat ve forme:

Ax = B (2.2)

kde

A = [aij ] =

a11 a12 . . . a1n

a21 a22 . . . a2n

. . . . . . . . . . . .an1 an2 . . . ann

(2.3)

x = xi =

x1

x2

...xn

(2.4)

b = bi =

b1b2...bn

(2.5)

A se nazývá maticí soustavy (ctvercová velikosti nxn), x a b jsou (sloupcové) vektory, bje tzv. vektor pravých stran a x obsahuje neznámé.

11

12 KAPITOLA 2. OPAKOVÁNÍ MATICOVÉ ALGEBRY

2.1 Rovnost matic

Matice A a B stejných rozmeru mxn jsou si sobe rovny, jestliže všechny jejich cleny jsousi rovny, t.j. aij = bij pro všechna i,j.

2.2 Scítání a odecítání matic

Scítání a odecítání je definováno pouze pro matice stejných rozmeru a to následovne:

C = A+B . . . cij = aij + bij (2.6)D = A−B . . . dij = aij − bij

2.3 Násobení skalárem

λA = [λaij ] (2.7)

2.4 Násobení matice maticí

Pro matice A (mxn) a B (nxp) je definován soucin AB následovne:

C = AB . . . cij =

m∑

k=1

aikbkj (2.8)

kde i = 1, 2, ...,m a j = 1, 2, ..., p. Nezapomínejme, že obecne AB 6= BA ale (AB)C = A(BC),tj. násobení matic je asociativní.

2.5 Transpozice matice

Transpozice matice A = [aij ] je dána následovne:

AT = [aji] (2.9)

pripomenme, že platí: (AB)T = BTAT .

2.6 Symetrické matice

Matice je nazývá symetrická, jestliže platí: AT = A nebo aij = aji pro všechna i,j.

2.7. NULOVÁ MATICE 13

2.7 Nulová matice

Nulovou matici budeme znacit 0. Muže být jakéhokoliv rozmeru, napr. nulová maticerozmeru 2x3

0 =

[0 0 00 0 0

]

(2.10)

2.8 Jednotková matice

I =

1 0 . . . 00 1 . . . 0. . . . . . . . . . . .0 0 . . . 1

(2.11)

Všimneme si, že platí: AI = A, Ix = x.

2.9 Diagonální matice

Je taková matice, ketrá má nenulové prvky pouze na hlavní diagonále. Napr.:

D =

11 0 0 00 −3 0 00 0 0 00 0 0 17

Takovouto matici budeme zapisovat zkráceným zápisem D = diag[11,−3, 0, 17].

2.10 Determinant matice

Determinant ctvercové matice A je skalární císlo, oznacované detA nebo |A|. Výpocetdeterminantu matice 2x2, nebo 3x3 je triviální, použijeme tzv. Sarussovo pravidlo:

det

[a11 a12

a21 a22

]

= a11a22 − a12a21 (2.12)

det

a11 a12 a13

a21 a22 a23

a31 a32 a33

= a11a22a33 + a12a23a31 + a21a32a13 − a13a22a31 − a12a21a33 − a23a32a11

(2.13)

2.11 Singulární matice

Ctvercová matice A je singulární, jestliže její determinant je roven nule detA = 0, cožznamená, že soustava rovnic nemá jedno jediné rešení. (Soustava rovnic má práve jednorešení práve tehdy, když detA 6= 0.

14 KAPITOLA 2. OPAKOVÁNÍ MATICOVÉ ALGEBRY

2.12 Pozitivne definitní matice

ctvercová matice A (nxn) se nazývá pozitivne definitní jestliže pro každý nenulový vektorx rozmeru n platí:

xTAx > 0

Pozitivne definitní matice není singulární.

2.13 Inverzní matice

Pro ctvercovou a regulární matici A (detA 6= 0 ) je inverzní matice A−1vytvorena tak, žeplatí:

AA−1 = A−1A = I (2.14)

Matice adjungovaná k matici A je definována jako AdjA = CTij, kde C:

Cij = (−1)i+jMij (2.15)

kde Mij je determinant redukované matice získané eliminací i-tého rádku a j-téhosloupce matice A. Pomocí adjungované matice sematice inverzní dá vypocítat násle-dovne:

A−1 =1

detAAdjA (2.16)

Dá se snadno ukázat, že (AB)−1 = B−1A−1.

Jestliže matice soustavy lineárních algebraických rovnic 2.1 A je regulární, rešení tétosoustavy se dá vyjádrit jako:

x = A−1b (2.17)

tudíž hlavním úkolem pri hledání rešení soustavy je nalezení inverzní matice soustavyA−1 .

2.14 Derivace a integrace matice

Mejme A(t) = [aij(t)]. Potom derivací matice A podle casové promenné t rozumíme:

d

dtA(t) =

[daij(t)

dt

]

(2.18)

a její integrací:

∫

A(t)dt =

[∫

aij(t)dt

]

(2.19)

Kapitola 3

Rešení soustav lineárníchalgebraických rovnic

V tomto odstavci se omezíme pouze na soustavy lineárních algebraických rovnic, vekterých bude matice A ctvercová a regulární, rádu n:

Ax = b (3.1)

Matice A je regulární práve tehdy, když detA 6= 0 a rovnice 3.1 má pak práve jednorešení. Metody pro rešení rovnice 3.1 mužeme rozdelit do dvou skupin:

• Prímé metody: tyto metody poskytují presné rešení rovnice 3.1 a to v konecnémpoctu kroku. Patrí sem Gaussova eliminacní metoda a metody využívající LU roz-kladu matice A.

• Iteracní metody: rešení je dosaženo pomocí série iterací, tedy výpocetních kroku,kterými se prinližejeme k presnému rešení. Výsledný vektor x je tedy pouze pri-bližný. Zmíníme se strucne o Jacobiho metode a Gauss-Seidelove metode.

Obecne muže být težké ríct, která z metod bude nejvýhodnejší pro rešení dané soustavyrovnic, itecní metody budou vhodnejší v prípade, že matice A je rídká, nebo dokonce vý-razne pásová. Matice tuhosti celého problému sestavená pri rešení úloh pomocí metodykonecných prvku je obvykle výrazne pásová, proto mají iteracní metody v MKP velkývýznam. Príme metody nelze rovnež nekdy použít z duvodu rozsáhlosti dané úlohy, kdyse matice nevejde do operacní pameti pocítace.

3.1 Prímé metody

3.1.1 Gaussova eliminacní metoda

Tato metoda je obecne známá a navíc není problém si ji odvodit. Rovnici 3.1 rozepišmepomocí jednotlivých clenu:

a11x1 + a12x2 + . . .+ a1nxn = b1

a21x1 + a22x2 + . . .+ a2nxn = b2 (3.2)...

an1x1 + an2x2 + . . .+ annxn = bn

15

16 KAPITOLA 3. REŠENÍ SOUSTAV LINEÁRNÍCH ALGEBRAICKÝCH ROVNIC

Nyní eliminujme prvky pod hlavní diagonálou tak, aby nám na posledním rádku vyšlajedna rovnice o jedné neznámé, napr. ve tvaru annxn = ξbnn. Prvním krokem bude eli-minace první neznámé x1 ze všech ostatních rovnic. První rovnici vynásobenou clenem−a21

a11pricteme k druhému rádku, první rovnici vynásobenou clenem −a31

a11pricteme k

tretímu rádku a tak dále postupne k poslední rovnici (n− 1 operací):

a11x1 + a12x2 + . . .+ a1nxn = b1

a(1)22 x2 + . . .+ a

(1)2n xn = b

(1)2 (3.3)

...a(1)n2 x2 + . . .+ a(1)

nnxn = b(1)n

Postupne eliminujme další rádky. Cleny, kterými pokaždé násobíme jednotlivé rovnice,budeme nazývat multiplikátory a obecne v k-tém kroku budeme k-tý rádek soustavynásobit multiplikátorem

mk−1ik = −a

k−1ik

ak−1kk

(3.4)

a tím nulovat prvky v k-tém sloupci. Po n− 1 krocích eliminujeme puvodní soustavu nasoustavu s horní trojúhelníkovou maticí:

a11x1 + a12x2 + . . .+ a1nxn = b1

a(1)22 x2 + . . .+ a

(1)2n xn = b

(1)2 (3.5)

...a(n−1)

nn xn = b(n−1)n

Z posledního rádku je zrejmé

xn =b(n−1)n

a(n−1)nn

(3.6)

Ostatní neznámé xi získáme zpetnou substitucí v poradí xn−1, xn−2, ..., x1.Algoritmus Gaussovy eliminacní metody je možno snadno zapsat napr. takto (bez testumožnosti delení, tedy ak−1

kk 6= 0):

pro k = 1, 2, ..., n− 1pro i = k + 1, k + 2, ..., n

mik=−ak−1ik

ak−1kk

pro j = k + 1, k + 2, ..., n+ 1ak

ij = ak−1ij +mk−1

ik ak−1kj

konec jkonec i

konec k

3.1.2 Rešení pomocí LU rozkladu

Pri rešení soustavy rovnic pomocí GEM se vlastne matice A prevedla na tzv. horní trojú-helníkovou matici, tedy matici, její cleny pod hlavní diagonálou jsou nulové. Podívejmese tedy podrobneji na rozklad matice A pomocí trojúhelníkové matice. Predpokládejme,že tento rozklad existuje a že tedy matici A mužeme zapsat ve tvaru:

A = L · U (3.7)

3.1. PRÍMÉ METODY 17

kde L je dolní trojúhelníková matice a U je horní trojúhelníková matice (Upper a Lower).To znamená (pro jednoduchost predpokládejme nyní matice 4x4:

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44

=

l11l21 l22l31 l32 l33l41 l42 l43 l44

·

u11 u12 u13 u14

u22 u23 u24

u33 u34

u44

(3.8)

Požijme tohoto rozkladu pro rešení rovnice 3.1:

Ax = (L · U) · x = L · (U · x) = b (3.9)

substitucíL · y = b (3.10)

a následným rešenímU · x = y (3.11)

Výhodou tohoto postupu samozrejme je, že rešení soustavy rovnic s horní trojúhelníko-vou maticí je triviální (pouze zpetná substituce). Pro rešení rovnice 3.11 použijeme tedyzpetnou substituci:

xn =bnunn

(3.12)

xi =1

uii

yi −i−1∑

j=1

uijxj

i = n− 1, n− 2, ..., 1 (3.13)

pro rešení rovnice 3.10 doprednou substituci:

y1 =b1l11

(3.14)

yi =1

lii

bi −i−1∑

j=1

lijyj

i = 2, 3, ..., n (3.15)

Metod na urcení LU rozkladu matice A je více a navíc LU rozklad matice není urcen

jednoznacne. Napríklad pro matici A =

(4 108 2

)

lze snadno najít ekvivalentní LU

rozklady ve tvaru:

A =

(4 108 2

)

=

(2 04 2

)(2 50 −9

)

=

(2 04 −1

)(2 50 18

)

Nejednoznacnost LU rozkladu je zrejmá již z toho, že soustava

LU = A

obsahuje n2 rovnic, ale je potreba urcit n(n+ 1) neznámých lij a uij. Lze tedy volne volitn neznámých. Obvykle se jako volí diagonální prvky matice L a to jako jednotkové, tedy

lij = 1, i = 1, 2, ..., n


Jedním z nejcasteji používaných algoritmu pro nalezení LU rozkladu matice A je Cho-leského algoritmus pro symetrické pozitivne definitní matice. Je-li matice A pozitivnedefinitní, dá se ukázat, že existuje horní trojúhelníková matice U , která má na diago-nále pouze kladné prvky tak, že platí

A = UTU (3.16)

Choleského algoritmus se pak dá zapsat takto:

pro j = 1, 2, ..., n

ujj =√

ajj −∑j−1

r=1 u2jr

pro i = j + 1, j + 2, ..., n

uij = 1ujj

(

aij −∑j−1

s=1 uisujs

)

konec ikonec j

Záverem ješte poznamenejme, že uvedený LU rozklad je ekvivalentní s Gaussovou ele-minacní metodou a proto je existuje LU rozklad matice A práve tehdy, když je možnéaplikovat Gaussovu eliminacní metodu. Rešení soustavy Ax = b pomocí LU rozkladu jevýhodné zejména tehdy, když vektor b je vícenásobný, tedy máme rešit soustavu rovnic,kde je více pravých stran. V MKP tato situace nastává pomerne casto, vektor b repre-zentuje vektor zatížení, tedy je možné snadno rešit príklady s nekolika zatežovacímistavy.

3.1.3 GEM pomocí LDLT rozkladu

Jak již bylo receno, metoda konecných prvku vede na rešení soustavy lineárních al-gebraických rovnic. Matice A bude obvykle symetrická, pozitivne definitní matice. Prosymetrické matice se GEM obvykle modifikuje a rešení soustavy Ax = b se pak na-zývá rešení Gaussovou eliminacní metodou s použitím LDLT rozkladu (horní trojúhel-níková matice je v prípade symetrie transponovanou dolní trojúhelníkovou maticí, tedyU = LT ). Podívejme se na rešení soustavy rovnic GEM pomocí LDLT rozkladu trochupodrobneji.

Predpokládejme tedy, že máme soustavu linárních algebraických rovnic se symetrickoumaticí A ve tvaru:

Ax = b (3.17)

v metode konecných prvku bude tato soustava reprezentovat známou rovnici Kr = f ,matice A tedy bude mít význam matice tuhosti konstrukce K a vektor pravé strany bbude reprezentovat vektor transformovaného zatížení f .

Rešení soustavy rovnic 3.17 budeme hledat pomocí LDLT rozkladu matice A, tedy roz-kladu na dolní trojúhelníkovou matici L, diagonální matici D a horní trojúhelníkovoumatici U = LT . V matici LT budou pod diagonálou pouze nulové prvky, výsledné rešeníbude jednoduché provést “odspodu”, nebot’ poslední rádek je vlastne jedna rovnice projednu neznámou. Ostatní neznámé získáme zpetným dosazením do zbývajících rovnic,nebot’ postupne v každé rovnici se bude vyskytovat jen jediná neznámá.

Hledejme tedy rozklad matice A pomocí posloupnosti následujících transformací:

L−1n−1L

−1n−2...L

−12 L−1

1 A = U (3.18)

3.1. PRÍMÉ METODY 19

kde U je horní trojúhelníková matice ve tvaru:

U =

1..

1−li+1,i 1−li+2,i 1

.−ln,i 1

(3.19)

v tomto zápisu se prvky matice U nazývají Gaussovy násobitele (potrebné pro nulovánívšech sloupcu pod hlavní diagonálou) a z úvahy nulové hodnote prvku pod diagonálouve výsledné matici lze psát:

li+j,j =ai+j,i

aii(3.20)

Máme-li horní trojúhelníkovou matici U , pro výraz 3.18 potrebujeme její inverzi. Uvedomíme-li si, že matice L je dolní trojúhelníková, její inverzi získáme pouhou zmenou znaménka,tedy vynásobením všech poddiagonálních prvku −1. Výraz 3.18 lze tedy snadno prepsata pro rozklad matice tuhosti získáme jednoduchý výraz:

A = L1L2...Ln−1U = LU (3.21)

Horní trojúhelníkovou matici U mužeme snadno rozložit na soucin diagonální matice Da horní trojúhelníkové, která bude mít na diagonále samé jednicky, oznacme ji U :

U = DU (3.22)

Pro jistotu matice ješte rozepíšeme:

D =

U11

U22

..Unn

(3.23)

U =

1 L1,2 L1,3 . L1,n

1 L2,3 . L2,n

. Li,j Li,n

. Ln−1,n

1

(3.24)

Dosadíme-li rozklad U = DU do 3.21, získáme:

A = LU = LDU (3.25)

Nyní využijeme symetrie matice tuhosti, tedy platnosti AT = A a provedeme transpozicirozkladu matice A (traspozice diagonální matice je ta samozrejme samá matice):

AT =(LDU

)T= U

TDTLT = U

TDLT (3.26)

Nezapomenme, že vzhledem ke zmínené symetrie matice tuhosti musí platit A = AT :

LDU = UTDLT (3.27)


Z tohoto zápisu vyplývá rovnost mezi maticemi U a LT . Využijme tedy této rovnosti adosad’me U = LT do 3.25. Obdržíme tedy výsledný LDLT rozklad symetrické matice A:

A = LDLT (3.28)

Máme-li matici A v kýženém LDLT rozkladu, lze tento rozklad pomerne jednoduše po-užít pro rešení soustavy rovnic Ax = b. Najdeme nejprve pomocný vektor b, pro kterýbude platit:

Lb = b (3.29)

Tento pomocný vektor b získáme snadno pomocí modifikace vektoru pravých stran stej-nými koeficienty, které již máme k dispozici z rozkladu matice A (jsou to stejné Gaus-sovy koeficienty):

b = L−1b = L−1n−1L

−1n−2...L

−12 L−1

1 b (3.30)

a dosad’me tento výraz do rovnice Ax = b (a zároven rozložme A pomocí LDLT rozkladu):

LDLTx = Lb (3.31)

Nyní tuto rovnici mužeme zleva vynásobit postupne inverzní maticí k L:

DLTx = b (3.32)

a posléze inverzní maticí D−1:LTx = D−1b (3.33)

Invezní matici k diagonální najdeme snadno a matice LT je horní trojúhelníková, to zna-mená, že poslední rovnice z této soustavy je rovnice o jedné neznámé. Ostatní neznámézískáme zpetnou substitucí.

3.2 Iteracní metody

Výhodou iteracních metod je zejména to, že nekladou tak velké nároky na pamet’ pocí-tace jako metody prímé. Matice A se nemusí celá ukládat do operacní pameti a protovelké úlohy, které není možno pocítat pomocí GEM, lze rešit iterativne.

Základní myšlenkou všech iteracních metod je prevedení rovnice Ax = b na tvar vhodnýk iterování, tedy tvar, ve kterém bude na levé strane pouze neznámý vektor x a pravástrana bude obsahovat modifikovanou matici A, napr. matici M :

x = Mx+ n (3.34)

nyní lze snadno napsat iteracní formuli ve tvaru:

x(k+1) = Mx(k) + n (3.35)

To znamená, že jsme zapsali (k + 1) priblížení hledaného vektoru x v závislosti na jehohodnote z predešlého iteracního kroku (k). Obecne lze popsat libovolnou iteracní me-todu v nekolika málo krocích:

1. volba pocátecního vektoru x(0) =

x(0)1 , x

(0)2 , ..., x

(0)n ,

3.2. ITERACNÍ METODY 21

2. výpocet (k + 1) (v prvním kroku tedy první) aproximace vektoru rešení v závislostina predchozí iteraci

x(k+1) = Mx(k) + n

3. test zastavovací podmínky (rozhodnutí o ukoncení ci pokracování iteracního pro-cesu):

x(k+1) − x(k) < δ

kde δ je zvolená presnost rešení.

4. odhadnutí chyby rešení; obecne lze zapsat∣∣∣x(k+1) − x

∣∣∣ < ε

slovne se dá vyjádrit, že (k + 1) iterace aproximuje presné rešení x s chybou ε.

Iteracní metody lze ješte rozdelit na stacionární a nestacionární. Pri nestacionárníchiteracních metodách se v prubehu výpoctu mení obecne iteracní matice M i iteracnívektor n, popr. pouze jeden z nich. V tomto textu se budeme zabývat pouze stacionár-ními metodami.

3.2.1 Jacobiho iteracní metoda

Zkusme si odvodit algoritmus JIM. Jde vlastne o to, prpsat rešenou soustavu rovnic

a11x1 + a12x2 + . . .+ a1nxn = b1

a21x1 + a22x2 + . . .+ a2nxn = b2 (3.36)...

an1x1 + an2x2 + . . .+ annxn = bn

na tvar vhodný k iterování, tedy potrebujeme na obou stranách rovnice vektor x. Tedypotrebujeme v i-té rovnici vyjádrit i-tou neznámou xi, tedy vydelme každou rovnici dia-gonálním prvkem aii:

x1 =1

a11(b1 − a12x2 − a13x3 − ...− a1nxn)

x2 =1

a22(b2 − a21x1 − a23x3 − ...− a2nxn)

...xn =

1

ann(bn − an1x1 − an2x2 − ...− ann−1xn−1)

Nyní vidíme, že jsme skutecne dostali tvar vhodný k iterování. Snadno lze totiž oznacit:

M =

0 −a12

a11· · · −a1n

a11

−a21

a220 · · · −a2n

a22

......

. . ....

− an1

ann− an2

ann· · · 0

a n =

b1a11b2a22

...bn

ann

(3.37)

a dospet tak k Jacobiho iteracní formuli ve tvaru:

x(k+1) = Mx(k) + n (3.38)


Všimneme-li si lépe rovnice 3.37, zjistíme, že ke stejnému výsledku lze dospet pomocíLDU rozkladu matice A, tedy rozkladu na dolní trojúhelníkovou matici L, diagonálnímatici D = [dii] a horní trojúhelníkovou matici U , pro jednoduchost v prípade matice(4, 4):

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44

=

l21l31 l32l41 l42 l43

+

d11

d22

d33

d44

+

u12 u13 u14

u23 u24

u34

(3.39)Soustavu rovnic lze tedy pomocí LDU rozkladu zapsat jednoduše takto:

(L+D + U)x = b (3.40)

Toto lze snadno upravit na tvar vhodný k iterování:

Lx+Dx+ Ux = b

Dx = −(L+ U)x+ b

x = −D−1(L+ U)x+D−1b

Matice M a vektor n v obecném zápisu iteracní metody budou mít tedy tvar

M = −D−1(L+ U)

n = D−1b

Jacobiho iteracní metodu lze tedy psát ve tvaru:

x(k+1) = −D−1(L+ U)x(k) +D−1b (3.41)

Algoritmus JIR mužeme pro matici A(n, n) a vektor b(n) zapsat ve dvou cyklech (m jepocet iterací, zde tedy jednoduše, bez zastavovací podmínky):

pro k = 1, 2, ...,mpro i = 1, 2, ..., n

x(k+1)i = 1

aii

(

bi −∑i−1

j=1 aijx(k)j −∑n

j=i+1 aijx(k)j

)

konec ikonec k

Duležitou otázkou u všech iteracních metod je otázka jejich konvergence. Ke konver-genci si uvedeme dve základní vety na záver této kapitoly.

Poznámka V metode konecných prvku budeme rešit obycejne rovnici Ax = b, kde ma-tice A bude symetrická podle hlavní diagonály. V tom prípade rozklad matice A nahorní trojúhelníkovou U , diagonální D a dolní trojúhelníkovou L bude snazší o to,že horní trojúhelníková matice je transpozicí dolní, tedy U = LT . Zmínený rozkladtedy budeme nazývat shodne s literaturou LDLT rozkladem matice A. Jednodušelze zapsat A = LDLT , napr. v prípade (4, 4) matice:

a11 a12 a13 a14

a21 a22 a23 a24

a31 a32 a33 a34

a41 a42 a43 a44

=

l21l31 l32l41 l42 l43

+

d11

d22

d33

d44

+

lT12 lT13 lT14lT23 lT24

lT34

(3.42)

3.2. ITERACNÍ METODY 23

3.2.2 Gaussova-Seidelova iteracní metoda

Jacobiova metoda se v praxi témer nepoužívá a to zejména proto, že Gauss-Seidelovametoda konverguje skoro vždy, když konverguje Jacobiho metoda, navíc muže konver-govat i v prípade, že Jacobiho metoda nekonverguje a rychlost konvergence je u GSMvetšinou vyšší, než konvergence JIM.

GSM narozdíl od JIM používá již vypocítané složky vektoru x(k+1) v dalším výpoctu.To je i duvodem, proc se GSM rovnež nazývá metoda postupných oprav. V (k + 1)-nímiteracním kroku to bude tedy vypadat následovne:

x(k+1)1 =

1

a11

(

b1 − a12x(k)2 − a13x

(k)3 − ...− a1nx

(k)n

)

x(k+1)2 =

1

a22

(

b2 − a21x(k+1)1 − a23x

(k)3 − ...− a2nx

(k)n

)

...x(k+1)

n =1

ann

(

bn − an1x(k+1)1 − an2x

(k+1)2 − ...− ann−1x

(k+1)n−1

)

Použitím LDU rozkladu matice A zapišme opet rešenou soustavu rovnic ve tvaru vhod-ném k iterování:

Lx+Dx+ Ux = b

(L+D)x+ Ux = b

(L+D)x = −Ux+ b

z toho plyne

x = (L+D)−1Ux+ (L+D)

−1b (3.43)

tuto rovnici mužeme rozrešit vzhledem ke (k + 1)-ní iteraci x(k+1)

x(k+1) = − (L+D)−1Ux(k) + (L+D)

−1b (3.44)

GSM má velkou výhodu a tou je možnost dokázat její konvergenci v prípade, že maticeA je symetrická a pozitivne definitní. Další výhodou GSM je menší pamet’ová nároc-nost oproti JIM. Pro GSM jsou duležité nekteré její další vlastnosti, zejména podmínkykonvergence. Uvedeme si je zde pouze bez dukazu, podrobnosti viz napr. [rektorys Varmetody, spíš nekde jinde].

Veta Je-li matice A symetrická a pozitivne definitní, Gaussova-Seidelova metoda prorešení soustavy rovnic Ax = b konverguje nezávisle na volbe pocátecního vektorux(0).

Veta (Postacující podmínka konvergence) Posloupnost iterací

x(k+1) = Mx(k) + n

konverguje pro libovolný pocátecní vektor x(0) a libovolný vektor n, když normaiteracní matice M splnuje nerovnost

‖M‖ ≤ q < 1 (3.45)


3.2.3 Jacobiho metoda s relaxací

Prepíšeme-li JIM (3.41) na tvar

x(k+1) = x(k) −D−1(L+D + U)x(k) +D−1b

to znamená, že (k + 1)-ní iteraci získáme jako soucet k-té iterace a tzv. rezidua r(k+1):

x(k+1) = x(k) + r(k+1),kde r(k+1) = D−1(

b−Ax(k))

vynásobíme-li zmínené reziduum r(k+1) vhodne zvoleným císlem ω, mužeme získat no-vou iteracní metodu, která muže konvergovat rychleji.Jacobiho iteracní metoda s relaxací se pak dá zapsat

x(k+1) = x(k) + ωD−1(

b−Ax(k))

(3.46)

kde k = 1, 2, ... a ω se nazývá relaxacní parametr. Ve složkách dostáváme vztah

x(k+1)i = x

(k)i +

ω

aii

bi −i−1∑

j=1

aijx(k)j −

n∑

j=1

aijx(k)j

(3.47)

3.2.4 Gaussova-Seidelova relaxacní metoda (SOR)

Tzv. superrelaxacní metoda se používá k urychlení konvergence. Gaussovu iteracní me-todu prepišme opet pomocí rezidua do tvaru

x(k+1)i = x

(k)i + r

(k+1)i

kde

r(k+1)i =

1

aii

bi −i−1∑

j=1

aijx(k+1)j −

n∑

j=1

aijx(k)j

i-té reziduum r(k+1)i vynásobme vhodným parametrem ω. Obdržíme

x(k+1)i = x

(k)i + ωr

(k+1)i = (1− ω)x

(k)i

ω

aii

bi −i−1∑

j=1

aijx(k+1)j −

n∑

j=1

aijx(k)j

(3.48)

Gaussovu-Seidelovu relaxacní metodu (Successive OverRelaxation, SOR) lze tedy zapsat

x(k+1) = x(k) + ωD−1(

b−Mx(k+1) − (D +N)x(k))

(3.49)

Dá se ukázat, že aby metoda SOR konvergovala, je nutné, aby parametr ω splnovaljednoduchou konvergencní podmínku

0 < ω < 2 (3.50)

Pro úplnost zapišme ješte algoritmus SOR:

pro k = 1, 2, ...,mpro i = 1, 2, ..., n

x(k+1)i = (1− ω)x

(k)i + ω

aii

(

bi −∑i−1

j=1 aijx(k+1)j −∑n

j=i+1 aijx(k)j

)

konec ikonec k

3.3. STABILITA A CITLIVOST 25

3.3 Stabilita a citlivost

3.3.1 Stabilita

Duležitým aspektem pri numerickém rešení problému algebry jsou stabilita a citlivostdané metody.

Rozsáhlé úlohy MKP vedou na velké množství rovnic ve výsledné soustave lin. alg. rov-nic. Nekteré úlohy mohou mít špatné numerické vlastnosti, napr. úloha muže být citlivána malé zmeny vstupních dat. To znamená, že malá zmena ve vstupních datech mužemít velký vliv na sledované veliciny ve výsledku.

Proto bývá nutné (i když je to casto pri numerických výpoctech opomíjeno), provést tzv.citlivostní analýzu daného problému.

Vše na príkladech

Teorie - odkaz na ...

Príklad Gaussovou eliminací vyrešte soustavu rovnic Ax = b, kde

A =

[0.0001 1

1 −1

]

a b =

[10

]

Uvažujte aritmetiku s konecnosu presností, kde strojová presnost je εM = 10−4.

3.3.2 Citlivost

Citlivost na perturbaci vstupních dat, tedy jak se projeví malá zmena vstupních dat(zatížení, materiálové charakteristiky, rozmery konstrukce...) na sledovaných výsledcích(velicinách jako napetí, deformace, maximální posun...).


Kapitola 4

Základní princip metodykonecných prvku

Principem metody je reprezentace složité geometrie zkoumané oblasti, at’ už rovinnénebo prostorové, pomocí velmi jednoduchých cástí (konecných elementu) spojených vtzv. uzlových bodech. Pod vnejším zatížením se zkoumané teleso deformuje a pro repre-zentaci skutecných posunutí jednotlivých uzlových bodu volí jednoduché funkce, kterétato posunutí aproximují. Pro stanovení soustavy podmínek rovnováhy se používá vari-acních principu mechaniky (princip minima potenciální energie).Potenciální energie zatíženého pružného telesa je dána souctem vnitrní energie (vý-sledek deformace, tzv. deformacní energie) a potenciální energie vnejších sil (zatížení).Jestliže je teleso v rovnováze, tato potenciální energie nabývá minimální hodnoty.Tento princip si ukážeme na jednoduchém príklade pružiny zatížené vnejší silou.

F [N]

ku

Obrázek 4.1: Pružina zatížená vnejší silou

Vyjádreme nyní vztahy pro vnitrní a vnejší energii tohoto systému:

• vnitrní, deformacní energie v deformované pružine

U =1

2(sila.v.pruzine)× (posunuti) =

1

2(ku)u =

1

2ku2

• potenciální energie vnejší síly F

Wp = (sila)× (posunuti) = −Fu

27

28 KAPITOLA 4. ZÁKLADNÍ PRINCIP METODY KONECNÝCH PRVKU

• celková potenciální energie systému

Π = U +Wp =1

2ku2 − Fu

• minimum potenciální energie

∂Π

∂u= ku− F = 0⇒ ku = F

nám dává podmínku rovnováhy pro zatíženou pružinu.

Malinko predbehneme a popíšeme si postup metody konecných prvku.

1. Diskretizace kontinua

2. Volba bázových funkcí posunutí (DISPLACEMENT MODEL)

3. Odvození matice tuhosti jednotlivého elementu pomocí variacního principutento krok vede k podmínce rovnováhy pro jednotlivý element, kterou budemeuvádet ve tvaru:

[k]ere = Re (4.1)

[k]e ...matice tuhosti elementu (ELEMENT STIFFNESS MATRIX)re...vektor uzlových posunutí(NODAL DISPLACEMENT VECTOR)Re...vektor uzlových sil(NODAL FORCE VECTOR)Prvky matice tuhosti se nekdy nazývají též príspevkové koeficienty (INFLUENCECOEFFICIENTS). Prvek matice tuhosti kij vyjadruje totiž sílu v míste a smeru ne-známé i od jednotkového zatížení v míste a smeru odpovídajícímu neznámé j.

4. Globalizace. Sestavení celkové matice tuhosti (matice tuhosti soustavy, globálnímatice tuhosti, GLOBAL STIFFNESS MATRIX) pro celé kontinuum. Sestavení glo-bálního vektoru zatížení (GLOBAL FORCE VECTOR)

[K]r = R (4.2)

5. Výpocet neznámých posunutí (globální vektor posunutír). Toto nelze samozrejmeprovést bez geometrických okrajových podmínek (uložení, podeprení kontrukce).

6. Výpocet odvozených velicin, napr. deformací a napetí ze známých uzlových posu-nutí.

Kapitola 5

Prímá metoda tuhosti(DIRECTSTIFFNESS METHOD)

Prutové konstrukce jsou patrne nejjednodušším príkladem sloužícím k demonstracizákladních principu MKP. Jsou to soustavy kloubove spojených prutu schopných pre-nášet pouze osové zatížení (tah, tlak). Jednotlivé pruty pritom mohou být myšlene na-hrazeny pružinami o tuhosti k = EA

L .

y

x

y

x

Obrázek 5.1: Príhradové konstrukce = hmotné body spojené kyvnými pruty

Patrne nejvetší výhodou bude fakt, že pri malém poctu prutu veškeré výpocty budemožno provést rucne bez pomoci pocítace a dále pak relativne snadná implementace(ukázáno bude pomocí MATLABu).Ukázková príhradovina sestává ze trí prutu a ze trí stycníku. V terminologii MPK sepruty nazývají elementy (ELEMENTS) a stycníky uzlové body (NODES, NODAL POINTS).Charakteristiky prutu budou oznaceny horním indexem v závorkách, napr. délky prutubudou L(i), moduly pružnosti E(i) a prurezové plochy A(i) .Prurezové plochy i modulypružnosti se predpokládají konstantní po celé délce prutu. Pruty (resp. elementy) budouobecne oznaceny horním indexem e, napr. prurezová plocha elementu je A(e).Geometrie konstrukce je vztažena ke globální soustave sourednic x, y (GLOBAL COORDINATESYSTEM, OVERALL COORDINATE SYSTEM). Externí síly (zatížení) mohou pusobit pouze

29

30 KAPITOLA 5. PRÍMÁ METODA TUHOSTI(DIRECT STIFFNESS METHOD)

y

x

5

3

Fx3 = 20 Fy3 = 10

4

Obrázek 5.2: Zadání príkladu

y

x

E1, A1, L1

E3, A3, L3

E2, A2, L2

fy3, uy3

fx3, ux3

fx2, ux2

fy2, uy2fy1, uy1

fx1, ux2

Obrázek 5.3: Koncové síly a koncová posunutí konstrukce

5.1. GLOBÁLNÍ SOUSTAVA ROVNIC (MASTER STIFFNESS EQUATION) 31

ve stycnících konstrukce. Globální vektor zatížení predstavuje sloupcový vektor o šestirádcích:

f =

fx1

fy1

fx2

fy2

fx3

fy3

(5.1)

Dalším duležitým prvkem bude vektor posunutí. Z mechaniky víme, že celkový stavpríhradové konstrukce je úplne dán posuny jednotlivých stycníku. Budeme-li tedy znátglobální vektor posunutí, budeme schopni vypocítat i napetí v jednotlivých prutech.Globální vektor posunutí je tedy:

r =

ux1

uy1

ux2

uy2

ux3

uy3

(5.2)

Tyto neznámé se také nazývají stupne volnosti (DEGREES OF FREEDOM, STATE VARIABLES).Dalším úkolem bude zajisté specifikovat okrajové podmínky (BOUNDARY CONDITIONS).Tyto hodnoty vlastne známe dopredu a mužeme je rovnou specifikovat v globálním vek-toru posunutí, ovšem z výpocetního hlediska toto muže pockat až jako poslední krokpred výpoctem soustavy rovnic.

5.1 Globální soustava rovnic (MASTER STIFFNESS EQUATION)

Tato soustava rovnic vyjadruje vztah mezi uzlovými silami f a stycníkovými posuvyr celé konstrukce pred urcením okrajových podmínek. Výsledná soustava lineárníchrovnic bude mít tvar:

fx1

fy1

fx2

fy2

fx3

fy3

=

K11 K12 K13 K14 K15 K16

K21 K22 K23 K24 K25 K26

K31 K32 K33 K34 K35 K36

K41 K42 K43 K44 K45 K46

K51 K52 K53 K54 K55 K56

K61 K62 K63 K64 K65 K66

ux1

uy1

ux2

uy2

ux3

uy3

(5.3)

v maticovém zápisu pakf = Kr (5.4)

kde K je globální matice tuhosti, která je symetrická,

5.2 Urcení matice tuhosti prvku

Prímá metoda spocívá ve trech krocích: (1) oddelení jednotlivých prvku, (2) lokalizace,(3) odvození matice tuhosti prvku.


V prvním kroku vedle globální soustavy souradnic zavedeme pro každý prvek ješte lo-kální souradnicový systém, v nemž odvodíme vztahy mezi vektorem posunutí a vek-torem uzlových sil prostrednictvím lokální matice tuhosti. Lokální soustavy souradnicbudeme znacit x(e), y(e) a úhel mezi globální osou x a lokální osou x(e) budeme znacitφ(e). Pocátek lokálního souradného systému je umísten do stredu každého prutu.

l

ui uj

vi vj

x

y

Obrázek 5.4: Tažený-tlacený prvek v lokálním souradném systému

Uzlové síly a uzlová posunutí jsou mezi sebou svázány pomocí:

f = Kr (5.5)

což pro prehlednost ješte rozepíšeme:

fxi

fyi

fxj

fyj

=

Kxixi Kxiyi Kxixj Kxiyj

Kyixi Kyiyi Kyixj Kyiyj

Kjixi Kxjyi Kxjxj Kxjyj

Kyjxi Kyjyi Kyjxj Kyjyj

uxi

uyi

uxj

uyj

(5.6)

f se nazývá vektor uzlových sil elementu a r vektor uzlových posunutí elementu a K jematice tuhosti elementu. Nejjednosušším zpusobem, jak stanovit prvky matice tuhostielementu je využít znalostí z mechaniky materiálu. Prut lze nahradit lineární pružinouo tuhosti

ks =EA

L(5.7)

tudíž vztah síla-posunutí nabývá tvaru:

F = ksu =EA

Lu (5.8)

kde F je vnitrní normálová síla v prutu a u je prodloužení (zkrácení) prutu. Normálovousílu i prodloužení lze zapsat pomocí uzlových sil a uzlových posunutí:

F = fxj − fxi d = uxj − uxi (5.9)

tyto rovnice vyjadrují podmínku rovnováhy a kinematickou podmínku kompatibility.Spojením 5.8 a 5.9 získáváme maticový zápis

fxi

fyi

fxj

fyj

=EA

L

1 0 −1 00 0 0 0−1 0 1 00 0 0 0

uxi

uyi

uxj

uyj

= Kr (5.10)

5.2. URCENÍ MATICE TUHOSTI PRVKU 33

tudíž

K =EA

L

1 0 −1 00 0 0 0−1 0 1 00 0 0 0

(5.11)

K je matice tuhosti taženého-tlaceného prvku v lokálních souradnicích.

5.2.1 Poznámky k odvození matice tuhosti taženého-tlaceného prvku

Odvození matice tuhosti prvku a význam jednotlivých clenu matice tuhosti lze osvetlitvelmi názorne na následujícím príkladu. Vztah mezi uzlovými silami

f

a vektoremuzlových posunutí r pomocí matice tuhosti prvku je dán:

fxi

fyi

fxj

fyj

=

K11 K12 K13 K14

K21 K22 K23 K24

K31 K32 K33 K34

K41 K42 K43 K44

uxi

uyi

uxj

uyj

(5.12)

Prepišme si co vlastne znamená první rádek tohoto zápisu:

fx1 = K11ux1 +K12uy1 +K13ux2 +K14uy2 (5.13)

Nyní položme všechna koncová posunutí prutu krome ux1 rovna nule. Rovnice se námzredukuje na fx1 = K11ux1. Nyní ihned vidíme význam jednotlivých prvku matice tu-hosti, nebot’ platí: K11 = fx1

ux1.

Mužeme tedy slovne psát: Clen Kij matice tuhosti prvku je síla (moment) v míste asmeru neznémé i od jednotkového posunutí (pootocení) v míste a smeru neznáméj. Všechna ostatní posunutí (a pootocení) jsou rovna nule (držena).

Toto je znázorneno na obrázku 5.5. Posunutí v libovolném míste prutu lze spocítat jakou(x) = Fx

EA , celkové protažení (zkrácení) prutu je u = FlEA . Nyní není nic jednoduššího než

dosadit do výrazu pro K11:

K11 =fx1

ux1=

fx1

fx1lEA

=EA

l(5.14)

Na tomto míste snad ješte doplnme poznámku k symetrii matice tuhosti. Ta je videtrovnež z významu jejích jednotlivých prvku. Rozepišme si napríklad ješte její tretí rádek

fy1 = K31ux1 +K32uy1 +K33ux2 +K34uy2 (5.15)

a porovnejme s 5.13. Zatímco K13 = fx1

ux2se dá slovy vyjádrit jako: “jak velkou silou

musíme pusobit v míste 1 aby v míste 3 vzniklo jednotkové posunutí, K31 =fx2

ux1je: “síla

vzniklá v míste 3 od jednotkového posunutí v míste 1. Toto ilustruje názorne obrázek5.5.


l

ux1

fx1

EA

Obrázek 5.5: Význam prvku matice tuhosti prutu

5.3 Transformace souradnic

Pro posunutí vyjádrená v globální, resp. lokální souradnicové soustave platí následujícíjednoduchý transformacní vztah:

uxi = uxi cosφ+ uyi sinφ

uyi = −uxi sinφ+ uyi cosφ (5.16)uxj = uxj cosφ+ uyj sinφ

uyj = −uxj sinφ+ uyj cosφ

zapsáno v maticové forme:

uxi

uyi

uxj

uyj

=

c s 0 0−s c 0 00 0 c s0 0 −s c

uxi

uyi

uxj

uyj

(5.17)

kde c = cosφ a s = sinφ, matice se nazývá transformacní matice pro pole posunutí a znacíse T . Dle obrázku 5.6 je zrejmé, že se uzlové síly transformují dle fxi = fxi−fyi sinφ, cožvede k alternativnímu zápisu:

fxi

fyi

fxj

fyj

=

c −s 0 0s c 0 00 0 c −s0 0 s c

fxi

fyi

fxj

fyj

(5.18)

ve kterém se transformacní matice nazývá transformacní matice pro uzlové síly, znacitji budeme Tf a platí T = T T

f .

5.4 Globalizace

Od této chvíle budeme používat již dríve zmínený index (e), který znací velicinu týkajícíse elementu. Podmínka rovnováhy elementu v globálních souradnicích je dána

f (e) = K(e)r(e) (5.19)

Transformace z globálního souradného systému je dána pomocí

5.5. SESTAVENÍ GLOBÁLNÍ MATICE TUHOSTI 35

y

ux2 xux1

α

α

uy1

uy2

fy1

fx1

fx2

fy2

uy2

uy1

ux1

ux2

y

x

Obrázek 5.6: Transformace souradnic

r(e) = T (e)r(e) f (e) = (T (e))T f(e) (5.20)

Zavedením techto výrazu do 5.5 a porovnáním s 5.19 získáme duležitý transformacnívztah mezi maticí tuhosti elementu v globálním souradnicovém systému K(e) a maticítuhosti elementu v lokálním souradnicovém systému K

(e):

K(e) = (T (e))TK(e)T (e) (5.21)

Násobíme-li matici tuhosti elementu zleva transformovanou a zprava puvodní transfor-macní maticí T (e), získáváme po provedení techto operací

K(e) =E(e)A(e)

L(e)

c2 sc −c2 −scsc s2 −sc −s2−c2 −sc c2 sc−sc −s2 sc s2

(5.22)

kde K(e) je matice tuhosti elementu v globálních souradnicích.

ÚKOL: Matici tuhosti v globálních souradnicích lze snadno vycíslit pro náš uve-dený príklad. Stejne snadné je i zapsat podmínky rovnováhy pro všechny tri ele-menty v globálních souradnicích.

5.5 Sestavení globální matice tuhosti

Úkolem pri sestavování globální matice tuhosti je umístení príspevku matic tuhostijednotlivých elementu. Tento proces se obecne nazývá (MERGING). Tento proces se dá


fyzikálne interpretovat jako opetné sestavení jednotlivých elementu (nyní prutu) do celékonstrukce. V globální matici tuhosti scítáme príspevky jednotlivých elementu do tu-hosti celé konstrukce. Toto musí splnovat matematicky následující dve pravidla:

• Kompatibilita deformací: stycníkové posuny všech prutu spojených v jednom styc-níku musí být stejné.

• Rovnováha ve stycníku: v každém stycníku je splnena silová podmínka rovnováhy,t.j. vnitrní síly všech prutu jdoucích do stycníku vyrovnávají vnejší zatížení tohotostycníku.

Dle obrázku 5.7 musí naše prutová konstrukce splnovat následující:

(3)

(2)

f(3)3

f(3)3

f(2)3

f(2)3

f3

Obrázek 5.7: Rovnováha ve stycníku - f3 je vnejší zatížení stycníku

u(2)x3 = u

(3)x3 u

(2)y3 = u

(3)y3 (5.23)

fx3 = f(2)x3 + f

(3)x3 fy3 = f

(2)y3 + f

(3)y3

Postup sestavení globální matice tuhosti (matice tuhosti konstrukce) se dá symbolickyzapsat jako:

f (1) = K(1)r f (2) = K(2)r f (3) = K(3)r (5.24)

f = f (1) + f (2) + f (3) = (K(1) +K(2) +K(3))r = Kr (5.25)

5.6. REŠENÍ (SOLUTION) 37

Jestliže toto vše provedeme, naše výsledné rovnice (MASTER STIFFNESS EQUATIONS)budou mít tvar:

fx1

fy1

fx2

fy2

fx3

fy3

=

20 10 −10 0 −10 −1010 10 0 0 −10 −10−10 0 10 0 0 00 0 0 5 0 −5−10 −10 0 0 10 10−10 −10 0 −5 10 15

ux1

uy1

ux2

uy2

ux3

uy3

(5.26)

Tento postup sestavování globální matice tuhosti se však v praxi nepoužívá a to z násle-dujícího duvodu. Aby se soucet matic tuhosti jednotlivých elementu dal provést, muselyby tyto matice mít hodnost rovnou hodnosti globální matice tuhosti, v našem prípadehodnost 6. Bylo by tedy nutné pro každou matici tuhosti elementu alokovat pamet’ pronxn císel.

Z tohoto duvodu se provádí sestavování globální matice tuhosti pomocí tzv. kódovýchcísel, kdy se každému prvku nejprve priradí kódová císla, nekdy se hovorí o tzv. ta-bulce stupnu volnosti elementu. Matice tuhosti jednotlivých prvku nejsou zvetšeny navelikost globální matice tuhosti (navíc se vetšinou rovnou prihlíží k symetrii matic, projednoduchost v textu vynecháme). Pro náš prípad budou kódová císla u matice tu-hosti jednotlivých elementu vypadat následovne (císlo vždy odpovídá pozici príslušnéhostupne volnosti elementu v globálním vektoru uzlových posunutí):

C(1) = 1, 2, 3, 4 C(1) = 3, 4, 5, 6 C(1) = 1, 2, 5, 6

Pomocí techto kódových císel se potom príspevky tuhostí jednotlivých elementu dotuhosti celé konstrukce promítnou “umístením” prvku jednotlivých matic tuhostí ele-mentu takto:

K(1)11 +K

(3)11 K

(1)12 +K

(3)12 K

(1)13 K

(1)14 K

(3)13 K

(3)14

K(1)21 +K

(3)21 K

(1)22 +K

(3)22 K

(1)23 K

(1)24 K

(3)23 K

(3)24

K(1)31 K

(1)32 K

(1)33 +K

(2)11 K

(1)34 +K

(2)12 K

(2)13 K

(2)14

K(1)41 K

(1)42 K

(1)43 +K

(2)21 K

(1)44 +K

(2)22 K

(2)23 K

(2)24

K(3)31 K

(3)32 K

(2)31 K

(2)32 K

(2)33 +K

(3)33 K

(2)34 +K

(3)34

K(3)41 K

(3)42 K

(2)41 K

(2)42 K

(2)43 +K

(3)43 K

(2)44 +K

(3)44

(5.27)

Tento postup bude snad patrnejší z grafického znázornení. Schematicky pomocí barev-ného obrázku (5.8) je znázorneno, na kterých místech v globální matice tuhosti se vnašem prípade scítají príspevky matic tuhosti jednotlivých elementu.

5.6 Rešení (SOLUTION)

Nyní nám zbývá poslední krok: specifikovat okrajové podmínky a vyrešit soustavu rov-nic. Tím získáme globální vektor posunutí r. Temto neznámým se ríká primární,ostatní veliciny (deformace, napetí) se budou odvozovat na základe tohoto vektoru.

5.6.1 Aplikace okrajových podmínek redukcí

Jedním ze zpusobu, jak do soustavy rovnic aplikovat okrajové podmínky, je použít tzv.redukci, t.j. na místa známých posunu (v našem prípade to budou tri nuly - dve za


Obrázek 5.8: Grafické znázornení sestavení globální matice tuhosti

pevný kloub a jedna za posuvnou podporu). Bez specifikace okrajových podmínek nemásoustava rovnic rešení, nebot’ matice tuhosti je singulární (rádky a sloupce matice Kjsou lineární kombinací). Fyzikálne se toto dá interpretovat pohybem telesa jako tuhéhocelku (RIGID BODY MOTION), konstrukce “plave” v rovine x, y.Aplikujeme tedy okrajové podmínky (uložení konstrukce) následovne:

ux1 = uy1 = uy2 = 0 (5.28)

Dále mužeme specifikovat vektor vnejších sil f:fx2 = 0 fx3 = 20 fy3 = 10 (5.29)

Okrajové podmínky nám eliminují tri rádky a tri sloupce, což vede na následující sou-stavu rovnic:

10 0 00 10 100 10 15

ux2

ux3

uy3

=

02010

(5.30)

Tato rovnice se nazývá (REDUCED MASTER STIFFNESS SYSTEM). Tuto soustavu rovnicsnadno rucne vypocteme (napr. Gaussovou eliminací) a výsledný vektor posunutí je:

ux2

ux3

uy3

=

00.4−0.2

(5.31)

Toto je pouze cástecný vektor posunutí, celkový vektor by byl:

r =

0000

0.4−0.2

(5.32)

5.6. REŠENÍ (SOLUTION) 39

5.6.2 Aplikace okrajových podmínek modifikací

Aplikace okrajových podmínek redukcí má zjevnou nevýhodu, která se projeví pri pro-gramování (preusporádání matice, která je již v pameti). Tento zpusob aplikace okrajo-vých podmínek se tedy nepoužívá a okrajové podmínky se aplikují tzv. modifikací, prikteré do rádku a sloupcu matice tuhosti, které odpovídají nulovým neznámým posu-nutím, zapíšeme nulu, pouze na místo diagonálních clenu zapíšeme jednicku; v našemprípade je ux1 = uy1 = uy2 = 0, tedy soustava bude vypadat následovne:

fx1

fy1

fx2

fy2

fx3

fy3

=

1 0 0 0 0 00 1 0 0 0 00 0 10 0 0 00 0 0 1 0 00 0 0 0 10 100 0 0 0 10 15

ux1

uy1

ux2

uy2

ux3

uy3

(5.33)

Rešením 5.33 nyní získáme úplný vektor globálních posunutí, vcetne nulových hodnotposunutí v místech, kde byly predepsány okrajové podmínky. Tento postup má zjevnouvýhodu pri aplikaci nenulových okrajových podmínek, v našem prípade pri aplikacipredepsaných posunutí (pokles podpory). Motivací je opet vyhnutí se preusporádánímatice tuhosti celé konstrukce. Aplikace nenulových okrajových podmínek se pak pro-vede ve dvou krocích. Nejprve modifikujme globální matici tuhosti K tak, aby rádky spredepsanými posuny byly triviální rovnice vyjadrující vlastne pouze hodnotu prede-psaného posunu v daném míste. Pro jednoduchost uvažujme nenulová posunutí danátakto ux1 = 0.1 uy1 = −0.2 uy2 = −0.3. První, druhá a ctvrtá rovnice tedy pouze vyjadrujíhodnotu predepsaného posunutí v techto místech (napr. 1× ux1 = 0.1):

1 0 0 0 0 00 1 0 0 0 00 0 10 0 0 00 0 0 1 0 00 0 0 0 10 100 0 0 0 10 15

ux1

uy1

ux2

uy2

ux3

uy3

=

0.1−0.2

0−0.32010

(5.34)

Rešením této soustavy získáme modifikovaný vektor globálních posunutí r(M). V dalšímkroku mužeme tuto rovnici modifikovat tak, že všechny známé cleny z levé strany pre-suneme na pravou stranu a získáme tak modifikovaný vektor uzlových sil f (M). Rešenímtéto modifikované soustavy rovnic získáme úplný vektor uzlových posunutí.

5.6.3 Aplikace okrajových podmínek v maticové podobe

Aplikace nenulových okrajových podmínek se dá snadneji vyjádrit v maticové podobe.Nejprve rozdelme soustavu globálních rovnic Kr = f na dve soustavy rovnic:

[K11 K12

K21 K22

] [r1

r2

]

=

[f1

f2

]

(5.35)

Rozepíšeme-li si první rádek této soustavy:

K11r1 + K12r2 = f1 (5.36)


bude zrejmé, že pokud vektor r2 bude obsahovat pouze známá, predepsaná posunutía zároven vektor f1 pouze známé koncové síly, mužeme neznámé posuny vektoru r1

snadno urcit:

K11r1 = f1 −K12r2 (5.37)

Tato rovnice vyjadruje vlastne redukovaný systém rovnic. Jestliže budou okrajové pod-mínky nulové (homogenní), bude vektor r2, který obsahuje predepsané posuny nulovýa rovnice se zjednodušší na:

K11r1 = f1 (5.38)

Tento postup aplikace okrajových podmínek modifikací se nejcasteji zapisuje v maticovépodobe následovne:

[K11 0

0 I

] [r1

r2

]

=

[f1 −K12r2

r2

]

(5.39)

Toto je pouze redukovaný systém 5.37 rozšírený o triviální rovnici Ir2 = r2. Výhodoutohoto postupu je, že opet získáme úplný vektor posunutí vcetne predepsaných posunutír2.

5.7 Výpocet vnitrních sil (odvozených velicin)(POST PROCESSING)

Posledním krokem je výpocet odvozených velicin (DERIVERED or RECOVERED QUANTITIES).Odvozenými se tyto veliciny nazývají proto, protože jsou vypocteny na základe primár-ních neznámých, kterými jsou uzlová posunutí.

Jestliže vektor výsledných uzlových posunutí vynásobíme zleva globální maticí tuhostiK, získáme vektor uzlových sil, ve kterém jsou obsaženy nejen známé vnejší síly (zatí-žení), ale i reakce v podporách: fx1 = −20, fy1 = −20, fy2 = 10. Je snadné zkontrolovatpodmínky rovnováhy konstrukce jako celku.

Další odvozenou velicinou jsou normálové síly v prutech a jim odpovídající normálovénapetí. Normálové síly v jednotlivých prutech N (e) a normálová napetí σ(e) lze snadnospocítat na základe znalosti protažení jednotlivých prutu ∆(e):

∆(e) = u(e)xj − u

(e)xi N (e) =

E(e)A(e)

L(e)∆(e) σ(e) =

N (e)

A(e)(5.40)

Poznámka: Tato kapitola sloužila jako úvod do koncepce metody konecných prvku.Postupu MKP bylo vysvetleno na velmi jednoduchém typu prvku (tažený-tlacenýprut). V dalším výkladu bude znacení ponekud zjednodušeno, bude vynechán in-dex (e) znacící element, ale napríklad matice tuhosti elementu bude znacena ma-lým [k],k a globální matice tuhosti velkým [K],K.

5.8 Implementace taženého-tlaceného prvku v Octave

Výše uvedený príklad taženého-tlaceného prvku se pokusíme co nejjednodušeji imple-mentovat pomocí nekterého z jazyku symbolické algebry. Nejpoužívanejšími jazyky provedecko-technické výpocty jsou Matlab, jeho open-source implementace Octave, dále je

5.8. IMPLEMENTACE TAŽENÉHO-TLACENÉHO PRVKU V OCTAVE 41

možno použít napr. systém Macsyma, nebo další komercní systémy jako je Mathema-tica, ci Maple.

My se omezíme na implementaci za použití systému Octave, který je zdarma ke staženína adrese www.octave.org. Výhodou tohoto jazyka je zejména jeho témer stoprocentníkompatibilita s Matlabem, i když tech sto procent je diskutabilních. Programovací stylpoužitý v príkladech je volen tak, aby byl co nejnázornejší, nikoli nejefektivnejší, cinejpohodlnejší. Pro názornost to urcite postací.

Vezmeme jednoduchý príklad z mechniky prutových konstrukcí. Príhradovina složenáze trí prutu je namáhána jedinou šiknou silou dle obrázku . Pomocí prusecné metody, cistycníkové metody spocteme normálové síly v jednotlivých prutech a z podmínek rovno-váhy reakce ve vnejších vazbách. Výsledné síly jsou následující: reakce Ax = −35.255kN,Ay = −48.296kN a By = 12.941kN a osové síly v prutech S1 = +7.471kN, S2 = −14.953kNa S3 = +55.768kN. Výsledky zkontrolujeme pomocí MKP, kde z posunutí jednotlivýchuzlových bodu spocítáme tyto odvozené veliciny.

! "$#α = 45

o

Obrázek 5.9: Príhradová konstrukce zatížená šikmou silou

5.8.1 Matice tuhosti elementu

Jako první bude potreba vypocítat matici tuhosti taženého-tlaceného prutu v globálníchsouradnicích. Za tímto úcelem byla navržena funkce TrussStiffnesMatrix.m, která vracíjako výsledek matici tuhosti elementu rádu 4x4. Implementace je velmi názorná:


function k = TrussStiffnesMatrix (E, A, x1, y1, x2, y2)dx=abs(x2-x1);dy=abs(y2-y1);L=sqrt(dx^2+dy^2);c=dx/L;s=dy/L;k = (E*A/L)*[[ c^2, c*s, -c^2, -c*s],

[c*s, s^2, -s*c, -s^2],[-c^2, -s*c, c^2, s*c],[-s*c, -s^2, s*c, s^2]];

endfunction

Funkce tedy pro zadané souradnice prvku x vrací matici tuhosti elementu již v glo-bálních souradnicích. Matice tuhosti elemetu c. 1, 2, 3 vypocítáme pomocí volání tétofunkce ve tvaru:

k1 = TrussStiffnesMatrix (2.1e11, 1/12*80*150^3, 0, 0, 10, 0);k2 = TrussStiffnesMatrix (2.1e11, 1/12*60*120^3, 10, 0, 5, 10*sin(60*pi/180));k3 = TrussStiffnesMatrix (2.1e11, 1/12*60*120^3, 5, 10*sin(60*pi/180), 0, 0);

Tvar matice tuhosti elementu c. 3 (šikmý prut) zkontrolujme pomocí (matice musí býtsymetrická dle hl. diagonály):

octave:4> k3k3 =

4.9497e+05 4.9497e+05 -4.9497e+05 -4.9497e+054.9497e+05 4.9497e+05 -4.9497e+05 -4.9497e+05

-4.9497e+05 -4.9497e+05 4.9497e+05 4.9497e+05-4.9497e+05 -4.9497e+05 4.9497e+05 4.9497e+05

Pomocí funkce MergeKeIntoK.m budeme nyní sestavovat matici tuhosti celého kon-strukce (globální matici tuhosti). Toto provedeme tak,

function K_glob = MergeKeIntoK(i,j,k,l,K_loc,K_glob)ii=0;list=[i,j,k,l];for i=list

jj=0;ii+=1;for j=list

jj+=1;K_glob(i,j)+=K_loc(ii,jj);

endforendfor

endfunction

Pred vlastní sestavováním globální matice tuhosti je treba tuto inicializovat. Toto pro-vedeme pomocí nulové matice:

k=zeros(6,6);

Tento príkaz vytvorí nulovou matici k velikosti (6, 6), do které budeme postupne nacítatmatice tuhosti jednotlivý elemetu:

5.8. IMPLEMENTACE TAŽENÉHO-TLACENÉHO PRVKU V OCTAVE 43

k=MergeKeIntoK (1,2,3,4,k1,k);k=MergeKeIntoK (3,4,5,6,k2,k);k=MergeKeIntoK (1,2,5,6,k3,k);

Nyní nadefinujme funkci ApplyBCs, která bude mít za úkol aplikování okrajových pod-mínek (zatím pouze nulová posunutí). Toto provedeme modifikací globální matice tu-hosti k. Nyní ovšem nebudeme moci použít postupu popsaného v odstavci 5.6, ale po-užijeme postupu následujícího: “vyškrtneme” opet odpovídající rádky a sloupce maticeodpovídající nulovým posunutím. Tímto by ovšem matice tuhosti celé konstrukce bylasingulární, pro zachování regulární matice je potreba místo diagonálních prvku kii na-psat jednicky; k tomu nám poslouží následující procedura:

function K_mod = ApplyBCs (which,K_glob)for i=which

for j=1:rows(K_glob)K_glob(i,j)=K_glob(j,i)=0.0;if (i==j)

K_glob(i,j)=1.0;endif

endforendforK_mod=K_glob;

endfunction

Proceduru zavoláme pomocí

k_mod=ApplyBCs([1,2,4],k);

Nyní matice k_mod obsahuje globální matici tuhosti celé konstrukce v globálním sou-radném systému s elminovanými rádky a sloupci, které odpovídají nulovým posunutímv globálním souradném systému. Globální vektor posunutí r získáme rešením lineárnísoustavy algebraických rovnic Kr = f , kde f je globální vektor zatížení. Globální vektorzatížení bude mít tvar:

f=[0,0,0,0,50*sin(45*pi/180),50*cos(45*pi/180)]

Rešení soustavy lineárních algebraických rovniclze v systému octave získat jednodušepomocí soucinu inverzní matice k−1

moda transponovaného vektoru zatížení fT , tedy po-mocí následujícího kódu:

r=inverse(k_mod)*f’

Vektor r nyní obsahuje výsledná posunutí uzlových bodu. Pro náš príklad dostaneme:

rT = 0.00000, 0.00000, 0.00000, 0.00000,−1.04706, 0.60452

Posledním krokem je tzv. postprocessing, t.j. výpocet odvozených velicin, v našem prí-pade se bude jednat o deformace, napetí a normálové síly v jednotlivých prutech kon-strukce. Za tímto úcelem lze použít následující procedury:


function n_force = GetIntForces (E, A, x1, y1, x2, y2, u)dx=abs(x2-x1);dy=abs(y2-y1);L=sqrt(dx^2+dy^2);c=dx/L;s=dy/L;# ind=reverse_index(???) -> mapovani ix, iy, jx, jy na glob. vektor;# mozna volat funkci rovnou s u_elem, pocitanym predem v lokal. s.s.u_elem = [ c*u[ix]+s*u[iy], -s*u[ix]+c*u[iy], c*u[jx]+s*u[jy], -s*u[jx]+c*u[jy] ]; u[jy] ];e_elem = (u_elem[3]-u_elem[1])/L;n_force = E*A*e_elem;

endfunction

Vektor deformace (t.j. normálové deformace (protažení) jednotlivých prutu) lze vypocítatpomocí rozdílu uzlových posunutí delených délkou prutu ε(x) =

ui−uj

l a císelne námvyjde:

ε =

2.9647e− 06−9.8829e− 063.6884e− 05

Z vektoru deformace lze snadno spocítat vektor napetí σ, v našem prípade taženého-tlaceného prutu se jedná o normálová napetí v jednotlivých prutech, která lze spocítatpomocí Hookova zákona σ(x) = Eε(x):

σ =

0.62258−2.07547.7456

Normálové síly v jednotlivých prutech snadno získáme pomocí vztahuN(x) =∫

(A) σ(x)dA,v našem prípade prizmatického prutu vynásobením normálového napetí prurezovouplochou jednotlivých prutu; císelne:

N =

7.471−14.94355.768

Kapitola 6

Variacní principy

Približné metody používané v klasické mechanice (ale nejen v mechanice):

• klasické variacní metody: zejména Ritzova metoda. Tyto metody pracují s pozloup-nostmi tzv. bázových funkcí, které jsou ruzné od nuly témer všude na celé vyšet-rované oblasti Ω.

• metoda sítí: dnes již také klasická metoda, svého casu používaná na rešení desek.Jedná se v podstate o ...

• metoda konecných prvku: jedná se vlastne o zobecnenou variantu Ritzovy metody,kdy bázové funkce jsou ruzné od nuly pouze v nekterých prvcích. Výhodou metodyoproti klasické Ritzove je její nezávislost na tvaru zkoumané oblasti Ω a snadná de-finice okrajových podmínek. Rovnež matice výsledné soustavy má výrazne pásovýcharakter, což usnadnuje rešení (resp. umožnuje použití spec. metod).

Princip variacního poctu - mezi funkcemi urcitých vlastností najít takovou funkci, prokterou daný integrál (který je závislý na dané funkci a jejích derivacích) nabývá extrémníhodnoty. Hledané funkci se ríká extremála.

I =

∫ b

a

F (x, y, y′

) dx = max min

(6.1)

prípustné (konkurencní) funkceLagrange - diferenciáltotální diferenciál

F =∂f

∂xdx+

∂f

∂ydy +

∂f

∂zdz =

n∑

i=1

∂f

∂xidxi = ∇T f.h (6.2)

kde ∇f (operátor nabla) - gradient funkce f(x)variace funkce - velmi malý prírustek funkcevariace v mechanice u(x), ε(x), p(x), σ(x) - tyto potom nazýváme virtuální posuny, defor-mace, zatížení, napetí.Lze ukázat odvození matice tuhosti prvku pomocí variacních principu mechaniky (prin-cip minima potenciální energie, (PV posunutí, PV sil), Hellinger-Reissner, modifikovanévariacní principy...)Pro úplnost pripomenme nezbytné pojmy z klasické mechaniky (ci teorie pružnosti), sekterými budeme pracovat. Jedná se zejmána o pojmy:

45

46 KAPITOLA 6. VARIACNÍ PRINCIPY

6.1 Funkcionál celkové potenciální energie

V klasické Ritzove metode se hledaná funkce, v našem prípade funkce pruhybu w(x)aproximuje pomocí lineární kombinace bázových funkcí wi. Úkolem je najít funkci pru-hybu w(x), která vyhovuje diferenciální rovnici - v našem prípade diferenciální rovniciohybové cáry prutu:

−EIw′′

= M(x) (6.3)

Hledanou funkci w(x) aproximujeme tedy lineární kombinací bázových funkcí ψi:

w(x) =n∑

i=1

aiψi (6.4)

Bázové funkce ψi volíme tak, aby vyhovovaly okrajovým podmínkám a neznámé koefici-enty ai urcíme pomocí minimalizace funkcionálu. V našem prípade budeme minimali-zovat potenciální energii Π , tedy budeme hledat její parciální derivace podle neznámýchkoeficientu ai:

∂Π

∂a1= 0

∂Π

∂a2= 0 (6.5)

... (6.6)∂Π

∂an= 0

Minimalizace potenciální energie tedy vede k rešení soustavy n rovnic o n neznámých6.5. Více bude zrejmé na príkladu, zopakujme však nejdríve nekteré nutné pojmy zklasické teorie pružnosti. Budeme pracovat s pojmy jako deformacní energie, prácevnejších sil a diferenciální rovnice ohybové cáry.

6.1.1 Vnitrní (deformacní) energie pri tahu tlaku

Predpoklad: lineární izotropní materiál (samozrejme homogenní). Hustota deformacníenergie λ (plocha pod krivkou v diagramu σ − ε) je u lineárního materiálu zrejme:

λ =1

2σε (6.7)

navíc pro elastický materiál platí Hookuv zákon:

σ = Eε (6.8)

Pouze pripomenme, že normálové napetí pri tahu (tlaku) je rovno σ = N(x)A(x) .

Potom celková deformacní energie naakumulovaná v deformovaném objemu V se dávyjádrit dvema rovnocennými integrálními vztahy:

U =

∫

(V )

λ dV =

∫

(V )

1

2σε dV =

∫ l

0

∫

(A)

σ2

2EdAdx =

1

2

∫ l

0

σ2

EAdx (6.9)

6.1. FUNKCIONÁL CELKOVÉ POTENCIÁLNÍ ENERGIE 47

U =

∫

(V )

λ dV =

∫

(V )

1

2σε dV =

∫ l

0

∫

(A)

Eε2

2dAdx =

1

2

∫ l

0

EAε2 dx (6.10)

Pri použití prvního z nich a dosazení výrazu pro normálové napetí pri tahu/tlaku:

σ =N(x)

A(x)(6.11)

snadno dospejeme k výrazu:

U =

∫ l

0

∫

(A)

N2

2EA2dAdx =

1

2

∫ l

0

N2(x)

EAdx

6.1.2 Vnitrní (deformacní) energie pri ohybu

Opet predpokládejme lineární izotropní materiál. Pri odvozování výrazu pro deformacníenergii se zmení pouze výraz pro normálové napetí pri ohybu. To je pri tomto namáhánídáno pomocí:

σ =M

Iy (6.12)

Potom mužeme psát:

U =

∫

(V )

λ dV =

∫

(V )

1

2σε dV =

∫ l

0

∫

(A)

σ2

2EdAdx =

1

2

∫ l

0

∫

(A)

M2

2EI2y2 dAdx (6.13)

Výraz∫

(A) y2 dA je zrejme výraz pro moment setrvacnosti, tudíž lze výraz 6.13 zjednodu-

šit na výsledný

U =1

2

∫ l

0

M2 (x)

EIdx (6.14)

6.1.3 Vnitrní (deformacní) energie pri kroucení

Predpoklady zustávají stejné, pouze se mení zatížení. Prizmatický prut je namáhánkroutícím momentem Mk. Pri kroucení vzniká smykové napetí, které je dáno pomocí:

τ =Mk

Ipρ (6.15)

kde Ip polární moment setrvacnosti prurezu (napr. pro plný kruh je Ip = 12πr

4) a ρ jevzdálenost vyšetrovaných vláken od stredu (polární souradnice, polomer). Deformacníenergii urcíme opet integrací hustoty deformacní energie pre celý objem:

U =

∫

(V )

λ dV =

∫

(V )

1

2τγ dV =

∫

(V )

τ2

2GdV =

∫ l

0

∫

(A)

M2k

2GI2p

ρ2 dAdx (6.16)

Výraz∫

(A)ρ2 dA je polární moment setrvacnosti Ip, opet lze zkrátit na výsledné:

U =1

2

∫ l

0

M2k (x)

GIpdx


q0

F1 F2 F3

Obrázek 6.1: Ohýbaný nosník

6.1.4 Vnejší energie

Mejme nosník namáhaný pouze ohybem. Zatížení sestává z osamelých bremen F1, F2, ..., Fna spojitého zatížení q(x).

• spojité zatížení q(x)

• silové zatížení (osamelá bremena) Fi

práce vnejších sil (osamelých bremen a spojitého zatížení) je potom:

W =

n∑

i=1

Fiwi +

∫ l

0

q(x)w(x) dx (6.17)

6.1.5 Celková potenciální energie prutu

Celková potenciální energie je dána rozdílem deformacní energie vnitrních sil (napetí) aenergie vnejších sil (práce zatížení).

Π = U −W (6.18)

Pro náš problém ohybaného prutu dostáváme pro potenciální energii výraz:

Π = U −W =1

2

∫ l

0

M2(x)

EIdx−

n∑

i=1

Fiwi −∫ l

0

q(x)w(x) dx (6.19)

Pomocí variacního principu hledáme funkci pruhybu w(x) pomocí minimalizace poten-ciální energie Π. Mezi klasické variacní principy mechaniky patrí Lagrangeruv variacníprincip (vede na princip virtuálních posunutí) a Castiglianuv variacní princip (vede naprincip virtuálních sil).

6.2 Lagrangeuv a Castiglianuv variacní princip

Lagrangeuv variacní princip Pri variaci složek posunutí, splnujících geometrické pod-mínky uvnitr oblasti Ω i na její hranici Γ, je variace potenciální energie systémunulová.

6.3. RITZOVA METODA 49

Dusledkem Lagrangeova variacního principu je veta o minimu potenciální energie sys-tému, která ríká:

Veta o minimu potenciální energie Ze všech kinematicky prípustných stavu (t.j. stavukteré vyhovují okrajovým podmínkám, napr. podeprení prutu) u(x) nastane ten,pro který je nabývá potenciální energie minimální hodnoty. To znamená, že funk-cionál je stacionární.

δΠ = δU − δW = 0 (6.20)

Duálním principem k Lagrangeove variacnímu principu je Castiglianuv variacní princip.Tento vede k vete o minimu komplementární (doplnkové) energie a slovy se dá vyjádrittakto:

Castiglianuv variacní princip Pri variaci složek sil a napetí, splnujících statické pod-mínky rovnováhy uvnitr telesa Ω i na jeho hranici Γ, je variace doplnkové energiesystému nulová.

Dusledkem Castiglianova variacního principu je veta o minimu komplementární ener-gie:

Veta o minimu doplnkové energie Ze všech možných stavu napetí pružného telesa,které vyhovují podmínkám rovnováhy uvnitr i na hranici telesa, nastane práveten, pro který doplnková energie systému nabývá minimální hodnoty.

6.3 Ritzova metoda

Metoda, která reší tvar ohybové cáry práve pomocí Lagrangeova variacního principu senazývá Ritzova metoda. Tato metoda reší úlohy pružnosti energeticky (nalezení pruhy-bové cáry nosníku). Výsledné rešení je približné, ale s dosažením požadované presnosti.Presnost záleží na volbe aproximace ohybové cáry, tedy na volbe bázových funkcí wi(x).

Postup: Neznámou funkci pruhybu w(x) aproximujeme pomocí bázové funkce ψ, castejipomocí posloupnosti funkcí ψi:

w(x) =n∑

i=1

aiψi (6.21)

Bázové funkce ψi volíme tak, aby splnovaly okrajové podmínky na hranici Γ a neznámékoeficienty ai urcíme z podmínek minima potenciální energie:

∂Π

∂ai=∂ (Eext +Eint)

∂ai=∂ (U −W )

∂ai= 0 (6.22)

Nejjednodušší bude popsat algoritmus Ritzovy metody na jednoduchém príkladu hle-dání pruhybové cáry nosníku.

Príklad Pomocí Ritzovy metody stanovte maximální pruhyb konzoly délky l, ohybovétuhosti EI zatížené na volném konci osamelým bremenem F .


Rešení nalezneme postupne ve dvou krocích. Nejprve volme bázovou funkci co nejjedno-dušší, rekneme jako kvadratickou funkci (budeme potrebova spojitou druhou derivaci):

w(x) =n∑

i=1

aiψi = a1x2

Sestrojme funkcionál potenciální energie jako rozdíl práce vnejších sil a deformacníenergie vnitrních sil:

Π = Ee +Ei = W − U =EI

2

∫ l

o

(

w′′

(x))2

dx − Fw

=EI

2

∫ l

0

4a21 dx− Fa1l

2 =EI

24a2

1l− Fa1l2

Nyní hledejme extrém Π vzhledem k hledanému parametru a1:∂Π

∂a1= 4EIla1 − F l2 = 0

a1 =F l2

4EIl=

F l

4EI

Hledaná funkce pruhybu je tedy:

w(x) =F l

4EIx2

Tento výsledek dává oproti presné hodnote maximálního pruhybu wmax = Fl3

3EI 25%chybu, pokusme se tedy výsledek zpresnit použitím lépe zvolené bázové funkce. Hle-dejme tedy tentokrát pruhyb jako lineární kombinaci za použitím pouhých dvou clenu,napr. za použití snadno derivovatelných mocninných funkcí:

w(x) = a1x2 + a2x

3

Opet sestrojme funkcionál Π (druhá derivace w′′

(x) = 2a1 + 6a2x):

Π =EI

2

∫ l

0

4a21 + 24a1a2x+ 36a2

2x2dx− Fa1l

2 − Fa2l3

=EI

2

(4a2

1l + 12a1a2l2 + 12a2

2l3)− Fa1l

2 − Fa2l3

Parciální derivace Π podle všech hledaných parametru mají být nulové:∂Π

∂a1=

EI

2

(8la1 + 12a2l

2)− F l2 = 0

∂Π

∂a2=

EI

2

(12a1l

2 + 24la32

)− F l3 = 0

Rešení tšchto dvou rovnic dává hodnoty pro parametry a1 = Fl2EI a a2 = − F

6EI . Dosazeníma1, a2 do rovnice pruhybu, získáme hledané rešení ohybové cáry:

w(x) = a1x2 + a2x

3 =F l

2EIx2 − F

6EIx3

Toto rešení nám již dává pro hodnotu maximálního pruhybu zcela presnou hodnotu:

wmax = w(l) =F l3

2EI− F l3

6EI=

F l3

3EI

6.4. PRINCIP VIRTUÁLNÍCH PRACÍ 51

6.4 Princip virtuálních prací

Z Lagrangeova variacního principu minima potenciální energie vyplývá deformacní me-toda, obecne známá metoda pro výpocet staticky neurcitých konstrukcí. Užití Castig-lianova variacního principu minima doplnkové energie vede zase k silové metode. Obaprincipy vycházejí z obecne platného principu virtuální práce (PVP).

V mechanice jsme byli zvyklí rozdelovat PVP na na princip virtuálních posunutí (PVp)a princip virtuálních sil (PVs), dle toho, kterou z velicin považujeme za virtuální, tedymyšlenou, nekonecne malé hodnoty.

Strucne se tedy pro PVs dá zapsat:

δW =n∑

i=1

δFiri +m∑

j=1

δMjφj (6.23)

nebo, v prípade PVp:

δW =

n∑

i=1

Fiδri +

m∑

j=1

Mjδφj (6.24)

Princip virtuálních prací zjednodušene praví: “Je-li teleso v rovnováze, virtuální prácevnitrních sil se rovná virtuální práci vnejších sil. Vnitrní síly budou v prípade statickéanalýzy napetí, vnejší síly zatížení. Pro sestavení rovnice vyjadrující virtuální práci jedobré si uvedomit, že napetí “pracují” na deformacích, síly “pracují” na posunutícha momenty “pracují” na pootoceních. Podmínka rovnosti virtuálních prací vnejších avnitrních sil se pak jednoduše zapíše takto:

∫

(V )

δεTσ dV =

∫ l

0

δuTqx dx+ δrT R∫

(V )

([B]T δr)TE[B]r dV =

∫ l

0

([N ]δr)T qx dx+ δrT R (6.25)

δrT∫

(V )

[B]TE[B] dV

︸︷︷︸

[K]

r = δr (

∫ l

0

[N ]T qx dx+ R)︸︷︷︸

f

Kr = f Kr = f

Tato rovnice vyjadruje vztah mezi vektorem koncových sil na prutu f a vektorem kon-cových posunutí r. Tyto jsou mezi sebou svázány pomocí matice K, tedy matice tuhostielementu. V našem prípade taženého-tlaceného prvku platí, že jeho matice tuhosti je:

K =

∫

(V )

[B]TE [B] dV =

∫ l

0

∫

(A)

[B]TE [B] dAdx =

∫ l

0

[B]TEA [B] dx (6.26)

Vzpomenme, že EA je normálová tuhost prutu. Malinko predbehneme a dodejme, žepozdeji pri odvozování matice tuhosti pro jiný prvek dojdeme k obdobnému zápisu,napríklad pro prípad ohýbaného prutu bude matice tuhosti:

K =

∫ l

0

[B]TEI [B] dx (6.27)


kde soucin EI je ohybová tuhost prutu. Postupne bude možno zpusob sestavení maticetuhosti elementu zobecnit a dojdeme k obecnému zpusobu vyjádrení matice tuhostielementu:

K =

∫

(Ω)

[B]T

[D] [B] dΩ (6.28)

kde Ω bude integracní oblast, tedy oblast celého elementu (délka, plocha, objem) amatice [D] bude matice materiálové tuhosti. O tom ale více v následujících kapitolách.

Kapitola 7

Tažený-tlacený prut

7.1 Variacní formulace I (x)

Pokusme se tedy nyní o odvození matice tuhosti elementu pro tažený-tlacený prvekpomocí principu virtuálních prací. Tažený-tlacený prut je v rovine dán pomocí sourad-nic dvou uzlových bodu. V každém uzlovém bodu i a j definujme pouze jedno koncovéposunutí ui a uj . Pro tažený-tlacený prvek platí je tedy vektor koncových posunutí:

r =

ui

uj

(7.1)

Pomocí techto koncových posunutí a matice bázových funcí je aproximováno libovolnéposunutí uvnitr prvku, tedy posunutí u(x):

u(x) =[N1 N2

]

ui

uj

(7.2)

ve zkráceném zápisu

u(x) = [N ]r u = Nr (7.3)

Hledáme tedy bázové funkce Ni, pomocí kterých budeme aproximovat libovolné posu-nutí u(x) ze známých koncových posunutí

ui uj

T .

x = lx = 0

N1(x) = 1 −

x

l

x

x = 0 x = l

x

N2(x) = x

l

Obrázek 7.1: Bázové funkce taženého-tlaceného prutu na intervalu < 0, l >

Z obrázku 7.1 je zrejmé, že posunutí v míste x prutu se dá vyjádrit pomocí posunutíuzlových bodu jako:

u(x) = (1− x

l)ui +

x

luj (7.4)

53

54 KAPITOLA 7. TAŽENÝ-TLACENÝ PRUT

Matice bázových funkcí pro tažený-tlacený prut bude mít tedy tvar:

[N ] =[

1− xl

xl

](7.5)

Abychom mohli vyjádrit deformacní energii prutu (U =∫

(V )Eε2

2 dV ), budeme potrebovatvyjádrit vektor deformace pomocí koncových posunutí:

ε(x) =du(x)

dx=d[N ]

dxr =

[− 1

l1l

]ui

uj

(7.6)

ε(x) = [B]r ε= Br (7.7)

Pro prípad lineárního izotropního materiálu bude vektor napetí:

σ(x) = Eε(x) = E(εx − εx0) = E([B]r − εx0) (7.8)

kde εx0 je vektor pocátecních deformací, napr. od teplotních zmen εx0 = ∆l/l = α.∆t.l/l =α∆t

Použijme pro náš tažený-tlacený prvek Lagrangeova principu minima potenciální ener-gie. K sestavení funkcionálu Π budeme potrebovat výrazy pro vnitrní energii a pro po-tenciální energii vnejších sil. Vnitrní (deformacní) energii prutu mužeme zapsat pomocíhustoty deformacní energie (energie naakumulované v jednotce objemu dokonale pruž-ného materiálu stlaceného pusobením vnejších sil) λ:

U (e) =

∫

(V )

λ dV =

∫

(V )

1

2σε dV =

∫

(V )

Eε2

2dV =

=1

2

∫ l

0

εTEAε dx =1

2

∫ l

0

rTBTEABr dx =

=1

2

[

u(e)i u

(e)j

] − 1l

1l

EA[− 1

l1l

]

u(e)i

u(e)j

= (7.9)

=1

2

[

u(e)i u

(e)j

] EA

l2

[1 −1−1 1

]

dx

︸︷︷︸

[K](e)

u(e)i

u(e)j

=

=1

2

r(e)(T )

K(e)r(e)

Tažený-tlacený prut je zatížen pouze spojitým zatížením q(x) pusobícím shodne s podél-nou osou prutu. Práce vykonaná temito vnejšími silami je:

W (e) =

∫ l

0

q(x)u(x) dx =

∫ l

0

qNT r(e) dx =(

r(e))T∫ l

0

q

[1− x

lxl

]

dx =

=(

r(e))T

f (e) (7.10)

Potenciální energie je rozdíl vnitrní (deformacní) energie a práce vnejších sil (zatížení):

Π(e) = U (e) −W (e) =1

2

(

r(e))T

Kr(e) −(

r(e))T

f (e) (7.11)

7.2. VARIACNÍ FORMULACE II (ξ) 55

Princip minima potenciální energie nám ríká, že pri zatížení vnejšími silami zaujme pruttakovou polohu, pri které nabyde jeho potenciální energie minimální hodnoty:

(δr)T ∂Π

∂r= (δr)

T[K(e)r(e) − f (e)] = 0 =⇒ K(e)r(e) = f (e) (7.12)

Matice tuhosti taženého-tlaceného prvku je tedy:

[K](e) =

∫ l

0

EA

l2

[1 −1−1 1

]

dx (7.13)

[K](e)

=EA

l

[1 −1−1 1

]

(7.14)

a vektor transformovaného (zobecneného zatížení):

f (e) =

∫ l

0

qNT dx =

∫ l

0

q

[1− x

lxl

]

dx = q

[

x− x2

2lx2

2l

]x=l

x=0

=

[ql2ql2

]

(7.15)

Rovnice 7.11 a 7.12 vyjadrují duležité vztahy, které použijeme i dále v tomto textu,zapišme si je radeji ješte jednou prehledneji, tedy bez indexu (e) oznacujících, že sejedná o zápis platící pro prvek:

Π =1

2rTKr − rT f

Kr = f (7.16)

7.2 Variacní formulace II (ξ)

Pro úplnost ješte odvodíme matici tuhosti taženého-tlaceného prutu pomocí bázovýchfunkcí volených na intervalu 〈0, 1〉 namísto intervalu 〈0, l〉. Postup se nám bude hoditpozdeji u bázových funkcí složitejších elementu. Výhodné je totiž (zejména z hlediskanumerické integrace) volit práve interval jednotkový. Postup si mužeme ukázat naprí-klad pomocí transformace intervalu 〈0, l〉 na interval 〈0, 1〉. Bázové funkce N1 a N2 pre-jdou ve tvar:

N1 = 1− ξN2 = ξ

Pri integraci nesmíme zapomenout na Jacobián transformace J. Souradnice vyjádrenév jednotlivých souradných systémech jsou svázány pomocí vztahu:

ξ =x

l

Jacobián transformace bude tedy mít tvar:

J =dξ

dx=

1

l

56 KAPITOLA 7. TAŽENÝ-TLACENÝ PRUT

Nyní je už velmi snadné vyjádrit matici tuhosti taženého-tlaceného prutu:

K =

∫

BTDB dΩ =

∫ 1

0

[−11

]

EA[−1 1

] 1

ldξ =

=EA

l

∫ 1

0

[1 −1−1 1

]

dξ =EA

l

[ξ −ξ−ξ ξ

]1

0

=EA

l

[1 −1−1 1

]

Pro úplnost ješte doplnme odvození matice tuhosti taženého-tlaceného prutu pomocíprirozených souradnic. Více se prirozenými souradnicemi budeme venovat ve zvláštníkapitole, nyní si jenom strucne ukážeme postup. Prirozené souradnice se volí tak, žejejich hodnota v uzlových bodech elementu nabývá hodnot +1 a -1. Pro prípad taženého-tlaceného prutu musíme tedy provést transformaci souradnic z intervalu < 0; l > nainterval < −1; +1 >. Hledejme tedy lineární funkci ξ ve tvaru: ξ = ax + b. Tato funkcemusí splnovat následující podmínky: ξ(x = 0) = −1 a ξ(x = l) = 1. Z techto podmínekvyjde: b = −1 a a = 2

l . Hledanou transformaci lze tedy zapsat ve tvaru:

ξ =2

lx− 1

Nuní odvodíme lineární bázové funkce N1 a N2, pomocí kterých budeme aproximovatposunutí po prutu (nyní místo na intervalu < 0; l >na intervalu < −1; +1 >. Pomocíobrázku 7.2 snadno nahlédneme, že hledané funkce budou mít tento tvar:

N1(ξ) =1

2− 1

2ξ =

1

2(1− ξ)

N2(ξ) =1

2+

1

2ξ =

1

2(1 + ξ)

ξ ξ = +1ξ = −1

N1 = 1

2(1 − ξ)

ξξ = −1 ξ = +1

N2 = 1

2(1 + ξ)

Obrázek 7.2: Bázové funkce taženého-tlaceného prutu na intervalu < −1; 1 >

Nyní sestavíme metici tuhosti prvku pomocí výrazu K =∫

(Ω)BTDB dΩ. K tomuto bu-

deme potrebovat matici B, kterou získáme derivací bázových funkcí podle prirozenésouradnice ξ:

B =dN

dξ=

[

−1

2,1

2

]

Jacobián transformace bude:J =

dξ

dx=

2

l

Potom matice tuhosti nabyde již známého tvaru:

K =

∫

Ω

BTDB dΩ =

∫ 1

−1

[− 1

212

]

EA[− 1

212

] 2

ldξ =

=2EA

l

∫ 1

−1

[14 − 1

4− 1

414

]

dξ =2EA

l

[24 − 2

4− 2

424

]

=EA

l

[1 −1−1 1

]

Kapitola 8

Ohýbaný prut (PLANE BEAMELEMENT)

V první cásti mechaniky (statika) se zabýváme statickou analýzou konstrukcí a její velkácást je venována prubehum vnitrních sil prutových konstrukcí. Nejcastejším úkolem jetedy stanovit prubeh vnitrních sil (normálové N(x) posouvající T (x) a ohybového M(x),popr. kroutícího momentu Mx(x).

Pro odvození základních vztahu mezi geometrickými velicinami budeme vycházet zedvou teorií. První teorii budeme nazývat Euler-Bernoulliho která se rovnež nazývá kla-sická, ci inženýrská. Tato teorie predpokládá, že prurez prutu zustává i po deformacirovinný a je kolmý k deformované strednici prutu. Tato teorie je platná pro štíhlé pruty(štíhlost prutu je pomer délky a prurezové plochy).Druhá teorie bývá nejcasteji nazývána Timošenkova a stejne jako Euler-Bernoullihoteorie predpokládá zachování rovinnosti prícného rezu, avšak tento již není po zatíženíprutu kolmý k deformované strednici. Timošenkova teorie nezanedbává vliv smyku nadeformaci prícného rezu a je duležitá pro ohyb masivních nosníku.V následujícím odstavci odvodíme tedy matici tuhosti ohýbaného prutu založeného naEuler-Bernoulliho teorii, t.j. s temito predpoklady:

• rovinná symetrie: existuje podélná rovina symetrie, podle níž je nejen prícný rezprutu, ale i pusobící zatížení symetrické. Z této symetrie vyplývá i symetrie vnitr-ních sil podél stejné roviny.

• promenlivost prícného rezu: prícný rez prutu je po délce prutu nemenný, nebo jedán hladkou funkcí A(x). My budeme v následujících úvahách považovat prut zaprizmatický, t.j. bude platit: dA(x)/dx = 0.

• kolmost prícných rezu: rezy zustávají po deformaci rovinné a kolmé k deformo-vané strednici prutu (tenké pruty, u masivních prutu je nutno brát v úvahu vlivposouvající síly a Timošenkovu teorii).

• deformacní energie: na deformacní energii prutu má vliv pouze ohybový moment,vliv posouvající síly je v energetické bilanci zanedbán.

• linearizace: pocítáme pouze s teorií prvního rádu, t.j. deformace jsou dostatecnemalé. Viz teorie infitezimálních deformací, teorie prvního rádu.

• elastické (pružné) chování materiálu: materiál prutu považujeme za lineární izot-ropní a homogenní.

57

58 KAPITOLA 8. OHÝBANÝ PRUT (PLANE BEAM ELEMENT)

Uvažujme prizmatický prut délky l který je namáhán pouze ohybem. V krajních bodechprutu definujme dve uzlové neznámé – posunutí ui (resp. uj ) a pootocení prícného rezuφi (resp. φj ) dle obr. 8.1.

l

vi vj

ξ ∈< −1; 1 >

y

φi φj

x ∈< 0, l >

Obrázek 8.1: Uzlové neznámé pro ohýbaný prut

Dle Euler-Bernoulliho hypotézy o zachováníx kolmosti prícných rezu k deformovanéstrednici prutu, lze posunutí ve smeru osy x vyjádrit pouze jako funkci vzdálenosti odnormálové osy y a úhlu pootocení prícného rezu φ. To je duvod, proc na prvku neuvažu-jeme posunutí ve smeru osy x jako neznámou a mužeme psát následující geometrickévztahy:vektor posunutí:

u(x, y)v(x, y)

=

−y ∂v(x)∂x

v(x)

=

−yv′y

=

−yφv(x)

(8.1)

Pri ohybu bude vznikat pouze deformace εx. Vektor deformace bude mít tedy pouzejednu složku εx, index x tedy vynecháme a zapíšeme:

ε =∂u

∂x= −y ∂

2v

∂x2= −y d

2v

dx2= −yκ (8.2)

kde κ = d2vdx2 je krivost prutu.

vektor napetí:

σ = Eε = −Ey d2v

dx2= −Eyκ (8.3)

ohybový moment:

M =

∫

(A)

−yσ dx = Ed2v

dx2

∫

(A)

y2 dA = EIκ (8.4)

8.1 Ohýbaný prut (Euler–Bernoulli)Nyní odvodíme matici tuhosti ohýbaného prutu. Prut má dva uzlové body, v každém znich dve defomacní neznámé - vi,j posunutí ve smeru osy y, φi,j - pootocení koncovýchprurezu (obr. 8.2).Posunutí libovolného místa prutu urcíme pomocí matice bázových (tvarových, SHAPEFUNCTIONS) funkcí [N ] a vektoru koncových posunutí prutu r:

v =[N1 N2 N3 N4

]

vi

φi

vj

φj

= [N ]r(e) (8.5)

8.1. OHÝBANÝ PRUT (EULER–BERNOULLI) 59

l

x

z

φ1

φ2

v1v2

Obrázek 8.2: Euler–Bernoulliho ohýbaný prut - geometrické vztahy

bázové funkce Ni napíšeme opet pomocí prirozené souradnice ξ = 2xl − 1. Pripomenme,

že prirozená souradnice má hodnotu -1 pro x = 0 a hodnotu 1 pro x = l. Tvar bázovýchfunkcí pro ohýbaný prut je na obrázku 8.3 a odvodíme si je pomocí následující úvahy.Dle jejich prubehu je zrejmé, že je lze zapsat jako kubické polynomy (viz. infexní body)a napr. první bázovou funkci budeme hledat ve tvaru:

vi = 0, φi = 1

vi = 1, φi = 0

N1(ξ)

N2(ξ)

vj = 0, φj = 1vi = 0, φi = 0

vi = 0, φi = 0

N3(ξ)

N4(ξ)

vj = 0, φj = 0

vj = 0, φj = 0

vj = 1, φj = 0

Obrázek 8.3: Bázové (tvarové) funkce ohýbaného prutu

N1(ξ) = aξ3 + bξ2 + cξ + d

Koeficienty a, b, c a d nalezneme pomocí okrajových podmínek:

N1(−1) = 1

N′

1(−1) = 0

N1(+1) = 0

N′

1(+1) = 0

Derivace N1 dle ξ je rovna N ′

1(ξ) = 3aξ2 +2bξ+c, rešíme tedy následující soustavu rovnic:

−a+ b− c+ d = 1

3a− 2b+ c = 0

a+ b+ c+ d = 0

3a+ 2b+ c = 0


Koeficienty vyjdou: a = 14 , b = 0, c = − 3

4 , d = 12 . Hledaná bázová funkce je tedy N1(ξ) =

14ξ

3 − 34ξ + 1

2 = 14 (1− ξ)2 (2 + ξ).

Nyní odvodíme bázovou funkci N(ξ), t.j. bázovou funkci pro prípad jednotkového po-otocení v uzlu i. Jednotkové pootocení v uzlu i pro interval 〈0, l〉 znamená pootocenívelikosti φ1 = l

2 pro interval 〈−1,+1〉. Hledáme tedy funkci N2(ξ) = aξ3 + bξ2 + cξ + d snásledujícími okrajovými podmínkami:

N2(−1) = 0

N′

2(−1) =l

2N2(+1) = 0

N′

2(+1) = 0

Rešíme tedy následující soustavu rovnic (nyní již s výhodou maticove, nebot’ matice levéstrany zustává nemenná):

−1 1 −1 13 −2 1 01 1 1 13 2 1 0

abcd

=

0l200

Koeficienty vyjdou (použijme napr. opet Octave): a = l8 ,b = − l

8 ,c = − l8 a d = l

8 . FunkceN2(ξ) pro jednotkové pootocení uzlu i má tedy tvar: N2(ξ) = 1

8 l(ξ3 − ξ2 − ξ + 1) = 1

8 l(1 −ξ)2(1 + ξ). Zbývající funkce N3(ξ) a N4(ξ) ponecháme na ctenári, uvedeme zde pouzevýsledek:

N1(ξ) =1

4(1− ξ)2(2 + ξ)

N2(ξ) =1

8l(1− ξ)2(1 + ξ) (8.6)

N1(ξ) =1

4(1 + ξ)2(2− ξ)

N2(ξ) = −1

8l(1 + ξ)2(1− ξ)

vektor deformace:

ε(x) = −yκ = −y d2v

dx2= −yN ′′

r = −y 4

l2d2v(ξ)

dξ2= −y 4

l2d2N(ξ)

dξ2r = −y[B]r(e) (8.7)

κ = [B] r (8.8)

Pri odvození tohoto vztahu jsme použili pravidlo o derivaci složené funkce. Osvežme sizákladní kurz matematiky pomocí následujícího odstavce. Máme funkci f(x) a funkciξ = 2x

l − 1. Pri derivaci funkce f(x) musíme použít zmínenou vetu o derivaci složenéfunkce:

df(x)

dx=df

dξ

dξ

dx=df

dξ

2

l

Druhá derivace funkce f(x) dle promenné x tedy bude:

d2f

dx2=d(df(x)

dx )

dξ

dξ

dx=d2f

dξ22

l

2

l=

4

l2d2f

dξ2


Matice B tedy bude obsahovat 4l2 násobky druhých derivací bázových funkcí Ni(ξ) dle ξ:

[B] =4

l2d2N

dξ2=

4

l2

[d2N1

dξ2d2N2

dξ2d2N3

dξ2d2N4

dξ2

]

(8.9)

Derivací 8.6 snadno odvodíme:

[B] =1

l

[

6 ξl , 3ξ − 1, −6 ξ

l , 3ξ + 1]

(8.10)

Matici tuhosti ohýbaného prutu odvodíme podobným zpusobem jako u prípadu taženého-tlaceného prutu z principa minima potenciální energie. Deformacní energii ohýbanéhoprutu lze odvodit pomocí hustoty deformacní energie λ, v prípade lineárního elastickéhomateriálu bude λ = 1

2σε:

U =

∫

(V )

λ dV =1

2

∫

(v)

σε dV =

∫

(V )

σ2

2EdV =

∫ l

0

∫

(A)

M2(x)

2EI2y2 dAdx =

∫ l

0

M2(x)

2EIdx

=

∫ l

0

E2I2(

v′′

)2

2EIdx =

1

2

∫ l

0

EI

(d2v

dx2

)2

dx =1

2

∫ l

0

EIκ2 dx

Nyní proved’me pár úvah ohledne κ = Br: za prvé, κ je skalár, to znamená, že κ = κT :

Br = rTBT

S využitím tohoto vycíslíme hodnotu κ2:

κ2 = (Br)2 = BrBr = rTBT rTBT = rTBTBr

S využitím techto poznatku nyní upravne výraz pro deformacní energii ohýbanéhoprutu:

U =1

2

∫ l

0

EIrTBTBr dx

Celková potenciální enegie Π bude dána rozdílem deformacní energie prutu a prácevnejších sil:

Π =1

2

∫ l

0

EIrTBTBr dx−∫ l

0

qv dx

=1

2rT

∫ l

0

BTEIB dxr −∫ l

0

qNr dx

=1

2rTKr − rT f

V tomto nám již známém výrazu je f vektor transformovaného zatížení, který jsmezískali temito úpravami (práce vnejších sil je opet skalár, musí tedy platit W = W T :

W = W T =

(∫ l

0

qN dxr

)T

= rT

(∫ l

0

qN dx

)T

= rT

∫ l

0

NT qT dx = rT

∫ l

0

NT q dx

Celková energie prutu je tedy dána pomocí:

Π = U −W =1

2rTKr − rT f (8.11)


kde

K =

∫ l

0

EIBTB dx =

∫ 1

−1

EIBTB1

2l dξ (8.12)

je matice tuhosti ohýbaného prutu a lze ji numericky snadno vycíslit:

K(e) =1

2EIl

∫ 1

−1

1

l2

6 ξl

3ξ − 1−6 ξ

l3ξ + 1

[

6 ξl 3ξ − 1 −6 ξ

l 3ξ + 1]dξ =

=EI

2l3

∫ 1

−1

36ξ2 6ξ(3ξ − 1)l −36ξ2 6ξ(3ξ + 1)l(3ξ − 1)2l2 −6ξ(3ξ − 1)l (9ξ2 − 1)l2

36ξ2 −6ξ(3ξ + 1)l(3ξ + 1)2l2

dξ (8.13)

=EI

l3

12 6l −12 6l4l2 −6l 2l2

12 −6l4l2

a f je vektor transformovaného zatížení:

f =

∫ l

0

NT q dx =1

2ql

∫ 1

−1

NT dξ =1

2ql

∫ 1

−1

14 (1− ξ)2(2 + ξ)18 l(1− ξ)2(1 + ξ)14 (1 + ξ)2(2− ξ)− 1

8 l(1 + ξ)2(1− ξ)

dξ = ql

12112 l12− 1

12 l

(8.14)

Úkol: overte výpoctem hodnoty koncových sil ohýbaného prutu z obr. 8.4, t.j. vypocteteprvky matice tuhosti ohýbaného prutu prímou metodou. K výpoctu použijte silovounebo deformacní metodu, popr. mužete využít energetického principu.

6EIl2

4EIl

6EIl2

12EIl2

12EIl2

6EIl2

2EIl

6EIl2

vi = 0, φi = 1

vi = 1, φi = 0

vj = 0, φj = 0

vj = 0, φj = 0

Obrázek 8.4: Koncové síly ohýbaného prutu


Pro jistotu si toto mužeme provést metodou fiktivního nosníku. Pruhyb v libovolnémmíste nosníku je roven fiktivnímu ohybovému momentu delenému ohybovou tuhostínosníku:

v(x) =Mf (x)

EI

obdobne pootocení je rovno fiktivní posouvající síle delené ohybovou tuhostí:

φ(x) =Tf (x)

EI

Nosník uvolníme dle obrázku 8.5 kde jsou patrné i okrajové podmínky: vi = 1 a φi =0. Tyto použijeme jako okrajové podmínky pro rešení staticky urcité konzoly zatíženéstaticky neurcitými reakcemi z uvolneného vetknutí vlevo.

AB

vi = 1, φi = 0

Al

Al2

2

Bl

B

l

2

2

3l

Obrázek 8.5: Fiktivní nosník - stanovení koncových sil

Fiktivní nosník zatížíme plochou ohybového momentu, stanovíme náhradní bremena(viz. obr. 8.5) a mužeme psát:

v(0) =Mf

EI=

1

EI

(

−Bl l2

+Al2

2

2

3l

)

= 1

φ(0) =Tf

EI=

1

EI

(

Bl − Al2

2

)

= 0

Rešením techto rovnic snadno získáme hodnoty neznámých reakcí A,B:

A =12EI

l3, B =

6EI

l2


Zbývající reakce urcíme ze svislé a momentové podmínky rovnováhy na konzole zatíženésilami A,B.

Poznámka

Pro úplnost doplnme ve strucné podobe odvození matice tuhosti ohýbaného prutu po-mocí bázových funkcí Ni(x) definovaných na intervalu 〈0; l〉. Bázové funkce budeme opethledat jako kubické polynomy:

Ni(x) = Ax3 +Bx2 + Cx+D

Postupne budeme hledat konstanty A,B,C,D pro prípady jednotkového posunutí (po-otocení) v uzlech i a j. S využitím nekterého ze systému pro rešení úloh numerickéalgebry (Maple, Mathematica, Maxima) hledáme konstanty A,B,C,D tak, aby funkceNi(x) vyhovovala okrajovým podmínkám prutu. Pro první bázovou funkce N1(x) mámeokrajové podmínky: N1(x) = 1, N ′

1(x) = 0, N1(l) = 0 a N ′

1(l) = 0. Každá okrajová podmínkapredstavuje jeden rádek matice A. Rešíme tedy soustavu 4 rovnic ve tvaru Ax = b, kdematice A bude mít tvar:

A =

0 0 0 10 0 1 0l3 l2 l 13l2 2l l 0

a vektor b bude postupne

b1 =

1000

b2 =

0100

b3 =

0010

b4 =

0001

Ukažme si použití systému Maxima (aplikace okrajových podmínek: N1(0) = 1,N ′

1(0) = 0,N1(l) = 0, N ′

1(l) = 0) pro stanovení koeficientu A,B,C,D první bázové funkce N1(x):

(1) A: matrix([0,0,0,1],[0,0,1,0],[l^3,l^2,l,1],[3*l^2,2*l,1,0]);

(% o2)

0 0 0 10 0 1 0l3 l2 l 13l2 2l 1 0

(2) b: matrix([1],[0],[0],[0]);

(% o2)

1000

(3) (A ^^ -1) . b;

(% o3)

2l3

− 3l2

01

(4)

První bázová funkce N1(x) pro posunutí vi = 1 je tedy:

N1(x) =2x3

l3− 3x2

l2+ 1


Zbývající bázové funkce stanovíme postupnou zámenou vektoru b1 za b2, b3, b4. Zde jevýsledek:

N2(x) =x3

l2− 2x2

l+ x

N3(x) =3x2

l2− 2x3

l3

N4(x) =x3

l2− x2

l

Matici [B] vyjadruje vztah mezi vektorem deformace ε a vektorem koncových posunutí,získáme ji tedy jako druhé derivace bázových funkcí:

B =[

d2N1

dx2d2N2

dx2d2N3

dx2d2N4

dx2

]

=[

12xl3 − 6

l26xl2 − 4

l6l2 − 12x

l36xl2 − 2

l

]

matici tuhosti ohýbaného prutu získáme opet integracíK =∫

ΩBTDBdΩ; nezapomenme,

že nyní je však oblast Ω interval 〈0; l〉:

K =

∫ l

0

BTEIB dx = EI

∫ l

0

12xl3 − 6

l26xl2 − 4

l6l2 − 12x

l36xl2 − 2

l

[12xl3 − 6

l26xl2 − 4

l6l2 − 12x

l36xl2 − 2

l

]=

= EI

12l3

6l2 − 12

l36l2

6l2

4l − 6

l22l

− 12l3 − 6

l212l3 − 6

l26l2

2l − 6

l24l

Došli jsme tedy ke stejnému výsledku jako v prípade integrace na intervalu 〈−1; 1〉.Uved’me si ješte celý postup odvození matice tuhosti ohýbaného prutu pomocí bázovýchfunkcí volených na intervalu 〈0; l〉 pomocí systému Maxima:

A: matrix([0,0,0,1],[0,0,1,0],[l^3,l^2,l,1],[3*l^2,2*l,1,0]);b: matrix([1],[0],[0],[0]);n_1: (A ^^ -1) . b;b: matrix([0],[1],[0],[0]);n_2: (A ^^ -1) . b;b: matrix([0],[0],[1],[0]);n_3: (A ^^ -1) . b;b: matrix([0],[0],[0],[1]);n_4: (A ^^ -1) . b;bazove funkce tedy budou:N1: n_1[1]*x^3+n_1[2]*x^2+n_1[3]*x+n_1[4];B1: diff(N1,x,2);N2: n_2[1]*x^3+n_2[2]*x^2+n_2[3]*x+n_2[4];B2: diff(N2,x,2);N3: n_3[1]*x^3+n_3[2]*x^2+n_3[3]*x+n_3[4];B3: diff(N3,x,2);N4: n_4[1]*x^3+n_4[2]*x^2+n_4[3]*x+n_4[4];B4: diff(N4,x,2);integraci B^T.B ziskame matici tuhosti:B: matrix(B1,B2,B3,B4);


K: integrate(B . transpose(B), x, 0, l);spravnost prubehu bazovych funkci muzeme overit jejich vykreslenim:l:1;plot2d(N1,[x,0,l]);plot2d(N2,[x,0,l]);plot2d(N3,[x,0,l]);plot2d(N4,[x,0,l]);

8.2 Ohýbaný prut s vlivem smyku (Timošenko)

Timošenkova teorie narozdíl od Euler-Bernoulliho teorie uvažuje vliv posouvající síly nadeformacní energii prutu. Základní kinematické predpoklady se zmení – prícné rezy sicezustávají i nadále rovinnými, nejsou však již kolmé k pretvorené strednici (obr. 8.6).

γ1

v′

1 v′

2

γ2

v′

1=

∂v

∂x

l

v1

φ1

φ2

v2

Obrázek 8.6: Ohýbaný prut s vlivem smyku - geometrické vztahy

Úhel odklonu prícného rezu prutu po deformovaci od teoretického kolmého prícnéhorezu oznacíme γ. Celkové zkosení prícného rezu θ se bude rovnat souctu úhlu φ (Euler-Bernoulliho hypotéza, prícné rezy kolmé k deformované strednici) a úhlu γ vznikléhood pusobení posouvající síly T (x):

θ = φ+ γ =∂v

∂x+ γ = u

′

+ γ (8.15)

Ostatní geometrické vztahy zustávají shodné s ohýbaným prutem rešeným dle Euler-Bernoulliho hypotézy, tedy:

u (x, y) = −yθ (8.16)v (x, y) = v(x)

Pro odvození matice tuhosti budeme potrebovat deformacní energii prutu. Prut je na-máhán ohybovým momentem M(x) a posouvající silou T (x):

U =1

2

∫ l

0

M2(x)

EI+T 2(x)

GAsdx (8.17)

8.2. OHÝBANÝ PRUT S VLIVEM SMYKU (TIMOŠENKO) 67

Dosad’me vztahy mezi vnitrními silami a rovnicí ohybové cáry, tedy: −EIv′′

= M(x) a−EIv′′′

= T (x) do 8.17:

U =1

2

∫ l

0

EI(

v′′

)2

+E2I2

GAs

(

v′′′

)2

dx (8.18)

Pripomenme, že As je redukovaná, efektivní plocha prícného rezu, která odolává smyku.Její velikost se spocítá z rovnosti deformacní energie pri smyku 1

2τγ = 12

τ2

GAs– tato

musí být rovna energii zpusobenou rozdelením smykového napetí po prurezu. Napr.pro obdélníkový prurez vychází efektivní prurezová plocha A

s = 56bh.

Nyní zaved’me smykovou štíhlost Φ jako dvanáctinásobek pomeru ohybové tuhosti EI asmykové tuhosti GAs delený kvadrátem délky:

Φ =12EI

GAsl2

Potreba násobení 12 a delení kvadrátem délky se projeví až pri scítání ohybové a smy-kové cásti matice tuhosti, pro smykovou štíhlost Φ = 0 se totiž Timošenkuv prut staneEuler-Bernoulliho prutem.

Pomocí výrazu pro smykovou štíhlost upravme výraz pro deformacní energii prutu do-sazením do 8.18:

U =1

2

∫ l

o

EI(

v′′

)2

+Φl2

12EI(

v′′′

)2

dx (8.19)

Do tohoto výrazu mužeme dosadit v′′

= Br. Nezapomenme, že díky úvahám o skalárnípovaze deformacní energie mužeme kvadrát

(

v′′

)2

psát jako rTBTBr:

U =1

2

∫ l

0

rTBTEIBr +Φl2

12rT(

B′

)T

EIB′

r dx =

=1

2

(

rT

∫ l

0

BTEIB dx r + rT

∫ l

0

Φl2

12

(

B′

)T

EIB′

dx r

)

= (8.20)

=1

2rT(Kohyb +Ksmyk

)r

Matice tuhosti ohýbaného prutu s vlivem posouvající síly je tedy:

K = Kohyb +Ksmyk

Kohyb =

∫ l

0

BTEIB dx (8.21)

Ksmyk =

∫ l

0

Φl2

12

(

B′

)T

EIB′

dx

V techto výrazech B′ je vektor obsahující tretí derivace bázových funkcí dle x:

B′

=

12l3

6l2 − 12

l36l2

(8.22)

Na záver pomocí výše odvozených vztahu vycísleme matici tuhosti Timošenkova prutu.První cást matice tuhosti (pouze vliv ohybu) známe, vycísleme tedy pouze vliv posouva-


jící síly Ksmyk :

Ksmyk =Φl2

12

∫ l

0

12l36l2

− 12l36l2

EI

12l3

6l2 − 12

l36l2

dx =

= ΦEI

12l3

6l2 − 12

l36l2

6l2

3l − 6

l23l

− 12l3 − 6

l212l3 − 6

l26l2

3l − 6

l23l

Po sectení Ksmyk a Kohyb dospejeme k matici tuhosti Timošenkova prutu, vyjádrenépomocí smykové štíhlosti Φ:

K = EI (1 + Φ)

12l3

6l2 − 12

l36l2

6l2

1l (4 + Φ) − 6

l21l (2− Φ)

− 12l3 − 6

l212l3 − 6

l26l2

1l (2− Φ) − 6

l21l (4 + Φ)

= (8.23)

=EI

l(1 + Φ)

12l2

6l − 12

l26l

6l 4 + Φ − 6

l 2− Φ− 12

l3 − 6l

12l3 − 6

l6l2 2− Φ − 6

l 4 + Φ

(8.24)

8.3 Ohýbaný prut na pružném podloží (Timošenko–Winkler)

Dalším duležitým jednodimenzionálním prvkem je ohýbaný prut na pružném, Winkle-rovském podloží. Pružné podloží lze modelovat jako soustavu pružin podepírajících prutpo celé jeho délce. Winklerovské podloží je zjednodušeným prípadem, protože zanedbávájak vliv vícedimenzionální pružnosti i trení.

Reakce pusobící na infinitesimální délce prutu dx je rovna dfw = −kwv(x) dx kde kw

je tuhost Winklerova podloží na jednotku délky. Její jednotky jsou[

Nm/m

]tedy

[Nm2

].

Hustota deformacní energie stlaceného podloží délky dx je λ = 12kwv(x)

2dx. Deformacníenergie prutu na pružném podloží je potom:

U = U + Uw

kdeUw =

1

2

∫ l

0

kwv2(x) dx (8.25)

Matici tuhosti prutu na Winklerovském podloží získáme prictením tuhosti podloží kmatici tuhosti prutu. Celková deformacní energie ohýbaného prutu na pružném podloží(Winkler–Timošenko) bude:

U =1

2

∫ l

o

EI(

v′′

)2

+Φl2

12EI(

v′′′

)2

dx+ kw

∫ l

0

v2 dx (8.26)

Proved’me toto pro prípad Timošenkova prutu. Matici tuhosti podloží získáme snadno:

Kw = kw

∫ l

0

v2(x) dx = kw

∫ l

0

BTB dx (8.27)

8.4. OHÝBANÝ PRUT – VLIV NORMÁLOVÉ SÍLY 69

Po vycíslení dostaneme:

Kw = kw

12l3

6l2 − 12

l36l2

6l2

4l − 6

l22l

− 12l3 − 6

l212l3 − 6

l26l2

2l − 6

l24l

(8.28)

Celková matice tuhosti ohýbaného prutu na Winklerovském podloží je souctem maticetuhosti prutu a matice tuhosti W. podloží:

K =EI

l(1 + Φ)

12l2

6l − 12

l26l

6l 4 + Φ − 6

l 2− Φ− 12

l3 − 6l

12l3 − 6

l6l2 2− Φ − 6

l 4 + Φ

+ kw

12l3

6l2 − 12

l36l2

6l2

4l − 6

l22l

− 12l3 − 6

l212l3 − 6

l26l2

2l − 6

l24l

(8.29)

8.4 Ohýbaný prut – vliv normálové sílyUvažujme prut namáhaný kombinací ohybového momentu M , posouvající síly T a nor-málové síly N . Nosník bude nyní mít v každém uzlu tri stupne volnosti – posunutí vesmeru osy x, posunutí ve smeru osy z a pootocení prícného rezu φ.Vektor uzlových posunutí r bude mít tedy celkem šest složek:

r =

u1

v1φ1

u2

v2φ2

Deformacní energie prutu bude:

U =1

2

∫ l

0

M2(x)

EI+T 2(x)

GAs+N2(x)

EAdx (8.30)

Nyní neexistuje žádná geometrická vazba mezi posuny a pootoceními. Matici tuhostiprutu získáme snadno sectením matice tuhosti Timošenkova ohýbaného prutu a taženého–tlaceného prutu (rozšírené pro zbývající stupne volnosti):

K =EI

l(1 + Φ)

0 0 0 0 0 00 12

l26l 0 − 12

l26l

0 6l 4 + Φ 0 − 6

l 2− Φ0 − 12

l3 − 6l 0 12

l3 − 6l

0 6l2 2− Φ 0 − 6

l 4 + Φ0 0 0 0 0 0

+EA

l

1 0 0 −1 0 00 0 0 0 0 00 0 0 0 0 0−1 0 0 1 0 00 0 0 0 0 00 0 0 0 0 0

(8.31)Na prut muže pusobit spojité zatížení libovolného smeru f0(x). Toto zatížení lze rozložitna zatížení pusobící v ose prutu n0(x) a zatížení kolmé ke strednici q0(x). Pro jednodu-chost uvažujme spojité zatížení konstantní hodnoty f0 = q0 + n0. Vektor transformova-ného zatížení bude:

f =

−n0l2

q0l2

q0l2

12n0l2

q0l2

− q0l2

12

(8.32)


Kapitola 9

Izoparametrické elementy

Nejcastejším typem používaných prvku se staly tzv. izoparametrické prvky, t.j. prvky unichž se k interpolaci jak geometrie, tak posunutí po prvku používá stejných interpo-lacních funkcí. Toto lze názorne zapsat:

x = [N ]T xn a u = [N ]T r (9.1)

Interpolacní (bázové) funkce [N ] se zapisují pomocí prirozené souradnice ξ, jejíž konceptbude vysvetlen v následujícím odstavci.

9.1 Prirozené souradnice

Jak už bylo receno v predchozích kapitolách, používáme dva typy souradných systému- lokální a globální. Lokální souradný systém je definován pro každý element a po-mocí transformacních zákonu jsou hodnoty prepocítávány do globálního souradnéhosystému. V tomto systému je obvykle definována geometrie celého zkoumaného telesa,zadáváno zatížení apod.

Prirozený souradný systém je lokální systém, který dovoluje popsat polohu libovolnéhobodu elementu pomocí nekolika bezrozmerných císel (v rovine dvou, v prostoru trí),jejichž velikost je vždy menší nebo rovna jedné. Tento systém je navíc vymyšlen tak,aby hodnota souradnice v uzlovém bodu elementu byla jednotková. Zavedení priroze-ných souradnic zjednodušuje formulaci konecných elementu a podstatne zjednodušujenumerickou integraci, nebot’ prevádí integraci na obecném intevalu na integraci na in-tervalu < −1; +1 >.

ll2 l1

(1, 0) P (ξ1, ξ2) (0, 1)x1 P (x) x2

Obrázek 9.1: Ke konceptu prirozených souradnic (1-D)

71

72 KAPITOLA 9. IZOPARAMETRICKÉ ELEMENTY

9.1.1 Prirozené souradnice pro 1-D elementy

Prirozené souradnice jsou takové souradnice, které dávají v uzlovém bodu elementujednotkovou hodnotu souradnice. Bude platit, že soucet všech prirozených souradnicje roven jedné, t.j.

n∑

i=1

ξi = 1

Pocet prirozených souradnic bude vždy o jednu vyšší, než pocet souradnic v odpovídají-cím kartézském souradném systému, t.j. pro jednodimenzionální prvky budeme potre-bovat dve prirozené souradnice ξ1 a ξ2, pro rovinné prvky tri souradnice ξ1, ξ2, ξ3 a vprostoru dokonce ctyri ξ1, ξ2, ξ3, ξ4. To je dáno už tím, že máme navíc jednu podmínkupro prirozené souradnice a to, že jejich soucet je roven jedné.

Podívejme se nyní na odvození prirozených souradnic pro tažený-tlacený prut z obrázku9.1. Transformacní vztah mezi kartézskou souradnicí libovolného bodu P x a priroze-nými souradnicemi ξ1, ξ2 lze zapsat:

1x

=

[1 1x1 x2

]ξ1ξ2

(9.2)

První rovnice pouze vyjadruje, že soucet prirozených souradnic je roven 1. Inverzívztahu 9.2 dostáváme:

ξ1ξ2

=

[x2/ (x2 − x1) −1/ (x2 − x1)−x1/ (x2 − x1) 1/ (x2 − x1)

]1x

=

[x2/l −1/l−x1/l 1/l

]1x

=(9.3)

=1

l

[x2 −1−x1 1

]1x

(9.4)

Názorne si ukážeme použití prirozeného souradného systému pro tažený-tlacený prut.Bázové funkce v prirozeném souradném systému mají tvar:

ξ1 =x2

l− x

l=x2 − xl

ξ2 =x

l− x1

l=x− x1

l

Všimneme si, že platí

ξ1 =l1l

a ξ2 =l2l

(9.5)

tedy, že prirozenou souradnici mužeme rovnež získat jako pomer délky "zbývajícího"intervalu k celkové délce. Využijme nyní prirozených souradnic pro 1-D prvek k odvo-zení matice tuhosti taženého-tlaceného prutu. Vzpomenme, že matice tuhosti prutu jeK =

∫

(Ω)BTDB dΩ. Matice B vyjadruje vztah mezi vektorem deformace ε a vektorem

koncových posunu r. Posunutí po prvku u(x) je aproximováno pomocí prirozenýchbázových funkcí N1 = ξ1 a N2 = ξ2 dle vztahu 9.2. Protože je vektor deformace ε rovenderivaci posunutí u(x) je matice B rovna:

B =dN

dξ=[

dξ1

dxdξ2

dx

]=[− 1

l1l

](9.6)

9.1. PRIROZENÉ SOURADNICE 73

Nyní mužeme již snadno odvodit matici tuhosti elementu integrací na oblasti elementuΩ:

K =

∫

(Ω)

BTDB dΩ =

∫ l

0

[− 1

l1l

]

EA[− 1

l1l

]dx =

= EA

∫ l

0

[1l2 − 1

l2

− 1l2

1l2

]

dx =EA

l

[1 −1−1 1

]

9.1.2 Trojúhelníkové souradnice

Podívejme se nyní na prirozený souradnicový systém pro rovinné elementy. V rovine setyto prirozené souradnice casto nazývají plošnými souradnicemi. Poloha bodu bude protrojúhelníkový prvek definována tremi parametrickými souradnicemi ξ1, ξ2, ξ3 (z nichžpouze dve budou nezávislé). Pro ctyrúhelník budeme definovat dve prirozené souradniceξ a ζ (viz. obr. 9.2).

y

x

A3

A1

A2

P (L1, L2, L3)

L1 = 0

L3 = 0

L2 = 0

(x1, y1)

(1, 0, 0)

(x2, y2)

(0, 1, 0)

(x3, y3)

(0, 0, 1)y

x

(1,−1)

(x3, y3)

(x4, y4) ζ

ξ

(−1, 1)

(x2, y2)

(1, 1)

(1,−1)

(x1, y1)

Obrázek 9.2: Prirozené souradnice v rovine

Plošné souradnice ξ1, ξ2, ξ3 lze definovat rovnež pomocí pomeru ploch trojúhelníku, nakteré rozdeluje libovolný bod P (x, y) puvodní trojúhelník, t.j.

ξ1 =A1

Aξ2 =

A2

Aξ3 =

A3

AA = A1 +A2 +A3 (9.7)

Plochu trojúhelníku s vrcholy o souradnicích (x1, y1), (x2, y2), (x3, y3) lze spocítat pomocídeterminantu matice obsahující souradnice jeho vrcholových bodu:

A =1

2

∣∣∣∣∣∣

1 1 1x1 x2 x3

y1 y2 y3

∣∣∣∣∣∣

(9.8)


Prirozené souradnice potom budou:

ξ1 =A1

A=

∣∣∣∣∣∣

1 1 1x x2 x3

y y2 y3

∣∣∣∣∣∣

∣∣∣∣∣∣

1 1 1x1 x2 x3

y1 y2 y3

∣∣∣∣∣∣

ξ2 =A2

A=

∣∣∣∣∣∣

1 1 1x1 x x3

y1 y y3

∣∣∣∣∣∣

∣∣∣∣∣∣

1 1 1x1 x2 x3

y1 y2 y3

∣∣∣∣∣∣

ξ3 =A3

A=

∣∣∣∣∣∣

1 1 1x1 x2 xy1 y2 y

∣∣∣∣∣∣

∣∣∣∣∣∣

1 1 1x1 x2 x3

y1 y2 y3

∣∣∣∣∣∣

(9.9)

ξ1 = 0, ξ2 = 0, ξ3 = 0(x2, y2)

ξ1 = 1, ξ2 = 0, ξ3 = 0(x1, y1)

ξ = 1

ξ = 0

ξ1 = 0, ξ2 = 0, ξ3 = 0

(x3, y3)

ξ = 3

4

ξ = 1

4

ξ = 1

2

P (ξ1, ξ2, ξ3)

Obrázek 9.3: Trojúhelníkové souradnice

Soucet takto definovaných plošných souradnic je opet jednotkový, jinými slovy soucetplošných souradnic kteréhokoliv bodu trojúhelníku je roven jedné. Je videt, že totovskutku platí:

ξ1 + ξ2 + ξ3 =A1

A+A2

A+A3

A= 1 (9.10)

Vztah mezi kartézskými souradnicemi x, y a plošnými (v tomto prípade je budeme nazý-vat trojúhelníkovými) souradnicemi ξ1, ξ2, ξ3 je dán pomocí:

1xy

=

1 1 1x1 x2 x3

y1 y2 y3

ξ1ξ2ξ3

(9.11)

Prirozené souradnice ξi mužeme rovnež získat prímo inverzí vztahu 9.11. Inverzi pro-vedeme pomocí systému pro symbolickou algebru Maxima. Nadefinejeme symbolickoumatici transformace a budeme hledat její inverzi. Použijeme príkazy matrix a invert.

maxima] T: matrix([1,1,1],[x1,x2,x3],[y1,y2,y3]);


(D1)

1 1 1x1 x2 x3

y1 y2 y3

(C1) T1: invert(T);

(D2)

x2y3−x3y2

x2y3−x1y3−x3y2+x1y2+x3y1−x2y1

y2−y3


x3−x2

x2y3−x1y3−x3y2+x1y2+x3y1−x2y1x3y1−x1y3


y3−y1


x1−x3

x2y3−x1y3−x3y2+x1y2+x3y1−x2y1x1y2−x2y1


y1−y2


x2−x1


V inverzní matici vidíme, že všechny cleny mají stejný jmenovatel, mužeme tedy maticizjednodušit pomocí ev (vytkneme determinant matice):

(C2) T_det:determinant(T);

(C3) N: T_det*T1;

(D3)

x2y3 − x3y2 y2 − y3 x3 − x2

x3y1 − x1y3 y3 − y1 x1 − x3

x1y2 − x2y1 y1 − y2 x2 − x1

Tento determinant je zároven roven dvojnásobku plochy trojúhelníka. Hledaný inverznívztah mezi trojúhelníkovými a kartézskými souradnicemi lze tedy zapsat ve tvaru:

ξ1ξ2ξ3

=

1

2A

x2y3 − x3y2 y2 − y3 x3 − x2

x3y1 − x1y3 y3 − y1 x1 − x3

x1y2 − x2y1 y1 − y2 x2 − x1

1xy

(9.12)

Stejného výsledku lze dosaáhnout i pomocí plošných souradnic. Prirozená souradniceξ1 je rovna pomeru ploch A1a A. Ukažme si pouze odvození první prirozené souradniceξ1. Budeme potrebat plochy trojúhelníku A1 a A:

A1 =1

2

∣∣∣∣∣∣

1 1 1x x2 x3

y y2 y3

∣∣∣∣∣∣

A =1

2

∣∣∣∣∣∣

1 1 1x1 x2 x3

y1 y2 y3

∣∣∣∣∣∣

Opet pomocí maximy:

(C4) A_1:matrix([1,1,1],[x,x2,x3],[y,y2,y3]);

(D4)

1 1 1x x2 x3

y y2 y3

Prirozenou souradnici ξ1 získáme delením ploch A1a A (resp. jejich dvojnásobku):

(C5) ev(A_1.invert(A),detout);

(D5)

1 1 1x x2 x3

y y2 y3

·

0

B

B

@

x2y3 − x3y2 y2 − y3 x3 − x2

x3y1 − x1y3 y3 − y1 x1 − x3

x1y2 − x2y1 y1 − y2 x2 − x1

1

C

C

A


Tedy stejný výsledek odpovídající prvnímu rádku získaného pomocí inverze matice trans-formace. Nyní použijme ješte Maximu k získání matice B. Budeme potrebovat derivacefunkcí ξi dle x a y. Sestavme nejprve funkce ξi:

(C6) xi_1: N[1,1]+N[1,2]*x+N[1,3]*y;

(C7) xi_2: N[2,1]+N[2,2]*x+N[2,3]*y;

(C8) xi_3: N[3,1]+N[3,2]*x+N[3,3]*y;


N1 = ξ1 = x2 y3 + x (y2 − y3)− x3 y2 + (x3 − x2) y

N2 = ξ2 = x (y3 − y1)− x1 y3 + x3 y1 + (x1 − x3) y

N3 = ξ3 = x1 y2 + x (y1 − y2)− x2 y1 + (x2 − x1) y

Pro lineární trojúhelníkový prvek (s uzlovými body pouze ve vrcholech) bude platit, žei-tá prirozená souradnice je prímo i-tou bázovou funkcí Ni(x, y). Abychom získali maticiB = ∂N budeme potrebovat derivace bázových funkcí Ni(x, y) podle x a y:

B =

dN1

dx 0 dN2

dx 0 dN3

dx 0

0 dN1

dy 0 dN2

dy 0 dN3

dydN1

dxdN1

dydN2

dxdN2

dydN3

dxdN3

dy

(9.13)

(C9) dxi_1_dx: diff(xi_1,x);

(C10) dxi_1_dy: diff(xi_1,y);





Matici B (pro jednoduchost vynecháme clen 12A ) pak získáme snadno:

(C15) B: matrix([dxi_1_dx, 0, dxi_2_dx, 0, dxi_3_dx, 0],

(C16) [0, dxi_1_dy, 0, dxi_2_dy, 0, dxi_3_dy],

(C17) [dxi_1_dy, dxi_1_dx, dxi_2_dy, dxi_2_dx, dxi_3_dy, dxi_3_dx]);

B =

y2 − y3 0 y3 − y1 0 y1 − y2 00 x3 − x2 0 x1 − x3 0 x2 − x1

x3 − x2 y2 − y3 x1 − x3 y3 − y1 x2 − x1 y1 − y2

Nyní by bylo již snadné odvodit matici tuhosti trojúhelníkového prvku. Potrebovalibychom pouze znát matici materiálové tuhosti D (v dalším textu bude uvedena ma-tice materiálové tuhosti pro prípad rovinné napjatosti a rovinné defomace). Obecne mámatice materiálové tuhosti pro rovinný prvek tento tvar:

D =

D11 D12 D13

D21 D22 D23

D31 D32 D33

Matici tuhosti trojúhelníkového prvku bychom získali integrací:

K =

∫

(Ω)

BTDB dΩ =

y2 − y3 0 x3 − x2

0 x3 − x2 y2 − y3

y3 − y1 0 x1 − x3

0 x1 − x3 y3 − y1

y1 − y2 0 x2 − x1

0 x2 − x1 y1 − y2

D11 D12 D13

D21 D22 D23

D31 D32 D33

(9.14)

y2 − y3 0 y3 − y1 0 y1 − y2 00 x3 − x2 0 x1 − x3 0 x2 − x1

x3 − x2 y2 − y3 x1 − x3 y3 − y1 x2 − x1 y1 − y2

∫

(Ω)

dΩ


9.1.3 Prirozené souradnice pro ctyrsten

Analogicky k plošným souradnicím lze snadno definovat objemové souradnice. Pro ctyr-sten z obrázku 9.4 dostáváme transformacní vztah mezi kartézskými a prirozenýmisouradnicemi:

yz

x

(x2, y2, z2)(0, 1, 0, 0)

(x3, y3, z3)(0, 0, 1, 0)

P (x, y, z)P (L1, L2, L3, L4)

(0, 0, 0, 1)(x4, y4, z4)

(x1, y1, z1)(1, 0, 0, 0)

Obrázek 9.4: Prirozenýné souradnice pro ctyrsten (3-D)

1xyz

=

1 1 1 1x1 x2 x3 x4

y1 y2 y3 y4z1 z2 z3 z4

ξ1ξ2ξ3ξ4

(9.15)

Objemové souradnice ξi mohou být definovány analogicky k plošným souradnicím jakoobjemové souradnice, tedy pomerem objemu protilehlého ctyrstenu k objemu celkovéhoctyrstenu:

ξi =Vi

Vpro i = 1, 2, 3, 4 (9.16)

Objem tetraedru s vrcholy o souradnicích (x1, y1), (x2, y2), (x3, y3), (x4, y4) je:

V =

∣∣∣∣∣∣∣∣

1 1 1 1x x2 x3 x4

y y2 y3 y4z z2 z3 z4

∣∣∣∣∣∣∣∣

(9.17)


prirozená souradnice ξ1 bude tedy dána pomereme determinantu:

ξ1 =V1

V=

∣∣∣∣∣∣∣∣

1 1 1 1x x2 x3 x4

y y2 y3 y4z z2 z3 z4

∣∣∣∣∣∣∣∣

∣∣∣∣∣∣∣∣

1 1 1 1x1 x2 x3 x4

y1 y2 y3 y4z1 z2 z3 z4

∣∣∣∣∣∣∣∣

(9.18)

Zbývající prirozené souradnice ξi lze odvodit obdobne. Prímocarejší je však odvodit pri-rozené souradnice pro ctyrsten pomocí inverze vztahu 9.15. Použijme Maximu a získámesnadno následující transformacní vztah:

8

>

>

<

>

>

:

ξ1ξ2ξ3ξ4

9

>

>

=

>

>

;

=1

6V

2

6

6

4

x2 (y3z4 − y4z3) − x3 (y2z4 − y4z2) + x4 (y2z3 − y3z2) −y3z4 + y2z4 + y4z3 − y2z3 − y4z2 + y3z2−x1 (y3z4 − y4z3) + x3 (y1z4 − y4z1) − x4 (y1z3 − y3z1) y3z4 − y1z4 − y4z3 + y1z3 + y4z1 − y3z1x1 (y2z4 − y4z2) − x2 (y1z4 − y4z1) + x4 (y1z2 − y2z1) −y2z4 + y1z4 + y4z2 − y1z2 − y4z1 + y2z1

−x1 (y2z3 − y3z2) + x2 (y1z3 − y3z1) − x3 (y1z2 − y2z1) y2z3 − y1z3 − y3z2 + y1z2 + y3z1 − y2z1

x3z4 − x2z4 − x4z3 + x2z3 + x4z2 − x3z2 −x3y4 + x2y4 + x4y3 − x2y3 − x4y2 + x3y2−x3z4 + x1z4 + x4z3 − x1z3 − x4z1 + x3z1 x3y4 − x1y4 − x4y3 + x1y3 + x4y1 − x3y1x2z4 − x1z4 − x4z2 + x1z2 + x4z1 − x2z1 −x2y4 + x1y4 + x4y2 − x1y2 − x4y1 + x2y1

−x2z3 + x1z3 + x3z2 − x1z2 − x3z1 + x2z1 x2y3 − x1y3 − x3y2 + x1y2 + x3y1 − x2y1

3

7

7

5

8

>

>

<

>

>

:

1x

y

z

9

>

>

=

>

>

;

(9.19)

Pro lineární ctyrsten (uzly pouze ve vrcholech), bude pro bázové funkce platit:

Ni = ξi

Derivací bázových funkcí Ni = ξi postupne dle x, y a z lze sestavit matici B = ∂N :

B =

dN1

dx 0 0 dN2

dx 0 0 dN3

dx 0 0 dN4

dx 0 0

0 dN1

dy 0 0 dN2

dy 0 0 dN3

dy 0 0 dN4

dy 0

0 0 dN1

dz 0 0 dN2

dz 0 0 dN3

dz 0 0 dN4

dz

0 ∂N1

∂z∂N1

∂y 0 ∂N2

∂z∂N2

∂y 0 ∂N3

∂z∂N3

∂y 0 ∂N4

∂z∂N4

∂y

∂N1

∂z 0 ∂N1

∂x∂N2

∂z 0 ∂N2

∂x∂N3

∂z 0 ∂N3

∂x∂N4

∂z 0 ∂N4

∂x

∂N1

∂y∂N1

∂x 0 ∂N2

∂y∂N2

∂x 0 ∂N3

∂y∂N3

∂x 0 ∂N4

∂y∂N4

∂x 0

(9.20)

a integrací K =∫

(Ω)BTDB dΩ získat matici tuhosti prvku.

9.2 Trojúhelníkový prvek s konstantní deformací (CST)

Využití konceptu prirozených souradnic si ukážeme na príklade odvození matice tuhostilineárního trojúhelníkového prvku. Je to nejjednodušší prvek pro rešení rovinných úlohpružnosti (rovinná deformace a rovinná napjatost). Díky své definici pomocí trí uzlovýchbodu je posunutí po prvku aproximováno lineárne a díky tomu je deformace prvkukonstantí. Casto tedy bývá nazýván prvek CST, tedy Constant Strain Triangle.

Základní rovnice získáme pomocí principu minima potenciální energie. Deformacníenergie bude:

U =

∫

(V )

1

2σT ε dV =

1

2

∫

(V )

(Dε)Tε dV =

1

2

∫

(Ω)

hεTDεdΩ (9.21)

9.2. TROJÚHELNÍKOVÝ PRVEK S KONSTANTNÍ DEFORMACÍ (CST) 79

kde h je tloušt’ka prvku. Práce vnejších sil bude rovna práci plošných sil b =bx by

T

na posunutí u =ux uy

T a práci sil, pusobících na hranici Γ (tyto síly bývají ozna-covány TRACTION FORCES, t:

W =

∫

(Ω)

huT b dΩ +

∫

Γ

rT t dΓ (9.22)

Potenciál deformacní energie:

Π = U −W =1

2

∫

(Ω)

hεTDεdΩ−(∫

(Ω)

huT b dΩ +

∫

Γ

rT t dΓ

)

=

=1

2

∫

(Ω)

hrTBTDBr dΩ−(∫

(Ω)

hrTNT b dΩ +

∫

Γ

rT t dΓ

)

=

=1

2rT

∫

(Ω)

hBTDB dΩr − rT

(∫

(Ω)

hNT b dΩ +

∫

Γ

t dΓ

)

=

=1

2rTKr − rT f

kde K =∫

(Ω) hBTDB dΩ je matice tuhosti prvku pro rovinnou deformaci/napjatost a

f =∫

(Ω) hNT b dΩ +

∫

Γ t dΓ je vektor transformovaného zatížení.

Posunutí libovolného bodu uvnitr prvku budeme aproximovat pomocí posunutí uzlo-vých bodu a trojúhelníkových souradnic ξi ve tvaru:

u(x, y) = u1ξ1 + u2ξ2 + u3ξ3

v(x, y) = v1ξ1 + v2ξ2 + v3ξ3 (9.23)

Toto zapíšeme pomocí maticového zápisu:

uv

=

[ξ1 0 ξ2 0 ξ3 00 ξ1 0 ξ2 0 ξ3

]

u1

v1u2

v2u3

v3

(u = Nr) (9.24)

Matici B, tedy matici vztahu mezi vektorem deformace a uzlových posunutí získámederivací matice bázových funkcí N . K tomuto využijeme opet geometrické rovnice:

ε = ∂u = ∂Nr = Br (9.25)

kde matice ∂ je opet tzv. operátorová matice, tedy matice symbolických derivací. V na-šem prípade rovinného problému má rozmer (3,2) a lze jí zapsat následujícícm zpuso-bem:

∂ =

∂∂x 00 ∂

∂x∂∂y

∂∂x

(9.26)

Pomocí operátorové matice ∂ vyjádreme geometrické rovnice pomocí vektoru koncovýchposunutí:


ε = ∂u = ∂Nr =

∂∂x 00 ∂

∂x∂∂y

∂∂x

[ξ1 0 ξ2 0 ξ3 00 ξ1 0 ξ2 0 ξ3

]

u1

v1u2

v2u3

v3

=

=

∂ξ1

∂x 0 ∂ξ2

∂x 0 ∂ξ3

∂x 0

0 ∂ξ1

∂y 0 ∂ξ2

∂y 0 ∂ξ3

∂y∂ξ1

∂y∂ξ1

∂x∂ξ2

∂y∂ξ2

∂x∂ξ3

∂y∂ξ3

∂x

u1

v1u2

v2u3

v3

=

=1

2A

y23 0 y31 0 y12 00 x32 0 x13 0 x21

x32 y23 x31 y31 x21 y12

u1

v1u2

v2u3

v3

(ε = Br)

Kde matice B je symbolická operátorová matice. Vztah mezi napetím a deformací je dánpomocí materiálové matice tuhosti D. V našem prípade rovinného problému vztah mezinapetím a pretvorením zapíšeme pomocí:

σ =

σx

σy

τxy

=

D11 D12 D13

D21 D22 D23

D31 D32 D33

εx

εy

γxy

(σ = Dε) (9.27)

Nyní mužeme matici tuhosti CST prvku odvodit dle následujícího vztahuK =∫

(Ω)hBTDB dΩ.

Je zrejmé, že pokud bude tloušt’ka prvku konstantní, lze integraci snadno provést (aprovést násobení výše uvedených trí matic). Matice CST prvku bude mít rozmer (6,6) anabyde následujícího tvaru:

K =

∫

(Ω)

hBTDB dΩ = h

∫

(Ω)

dΩBTDB =hA

4A2BTDB =

=h

4A2

y23 0 x32

0 x32 y23y31 0 x13

0 x13 y31y12 0 x21

0 x21 y12

D11 D12 D13

D21 D22 D23

D31 D32 D33

y23 0 y31 0 y12 00 x32 0 x13 0 x21

x32 y23 x13 y31 x21 y12

Nyní dokoncíme zbývajících pár kroku nezbytných k dodefinování vztahu pro CST pr-vek. Pro vyjádrení funkcionálu potenciální energie Π = U −W nám zbývá vyjádrit výrazpro práci vnejších sil. Uvažujme na prvku pouze objemové síly, definované pomocí slo-žek bx a by (jednotky N/m3):

b =

bxby

(9.28)

9.3. PODSTATA IZOPARAMETRICKÝCH PRVKU 81

vektor koncových sil získáme obdobným zpousobem, jako v prípade taženého-tlacenéhoprutu:

f =

fx

fy

=

∫

Ω

hNT b dΩ =

∫

Ω

h

ξ1 00 ξ1ξ2 00 ξ2ξ3 00 ξ3

bxby

dΩ (9.29)

Pro jednoduchost uvažujme tloušt’ku prvku h konstantní, stejne tak i složky zatížení bxa by. Potom s prihlédnutím k definici plošných souradnic lze psát:

∫

Ω

ξ1 dΩ =

∫

Ω

ξ2 dΩ =

∫

Ω

ξ3 dΩ =1

3A (9.30)

Po úprave:

f =

fx

fy

=Ah

3

bxbybxbybxby

(9.31)

9.3 Podstata izoparametrických prvku

Izoparametrické prvky používají stejné interpolacní funkce k aproximaci jejich geome-trie i posunu uvnitr prvku. V prípade CST prvku lze toto zapsat pomocí následujícíhomaticového zápisu:

1xyuv

=

1 1 1x1 x2 x3

y1 y2 y3u1 u2 u3

v1 v2 v3

ξ1ξ2ξ3

(9.32)

u = N1u1 +N2u2 +N3u3 = ξ1u1 + ξ2u2 + ξ3u3 (9.33)v = N1v1 +N2v2 +N3v3 = ξ1v1 + ξ2v2 + ξ3v3

V prípade trojúhelníkového prvku zjistíme, že pri prechodu z lineárního prvku na prvkyvyšších rádu (bázové funkce budou polynomy vyšších stupnu) jsou zpresnena pouzeposunutí, geometrie zustává definována stejne. Tyto elementy se nazývají superparame-trické. V prípade kvadratického trojúhelníkového prvku (LST - linear strain triangle) sejedná o nálsedující:

1xyuv

=

1 1 1 1 1 1x1 x2 x3 x4 x5 x6

y1 y2 y3 y4 y5 y6u1 u2 u3 u4 u5 u6

v1 v2 v3 v4 v5 v6

ξ1ξ2ξ3ξ4ξ5ξ6

(9.34)


Obecný zápis

1xyzuvwT

=

1 1 . . . 1x1 x2 . . . xn

y1 y2 . . . yn

z1 z2 . . . zn

u1 u2 . . . un

v1 v2 . . . vn

w1 w2 . . . wn

T1 T2 . . . Tn

N1

N2

...

Nn

(9.35)

Poslední rádek predstavuje teplotní pole po prvku, které je po prvku aproximováno zezadaných uzlových hodnot pomocí stejných aproximacních funkcí. Stejným zpusobemmuže být aproximována napr. tloušt’ka deskového, ci skorepinového prvku.

Konstrukci bázových funkcí pro mezilehlé uzlové body si ukážeme nejprve na prípadeLST prvku, pozdeji ukážeme to samé s použitím kvadratického ctyrúhelníkového prvkuQ8. U CST prvku platilo, že bázové funkce byly vlastne prirozené souradnice:

N1 = ξ1 N2 = ξ2 N3 = ξ3 (9.36)

(x2, y2)

(x3, y3)

(x1, y1)

1

ξ1

(x1, y1)

(x2, y2)

1

N1(ξ)

(x3, y3)

Obrázek 9.5: Bázové funkce: prechod ke kvadratickým báz. funkcím

Podíváme-li se na obrázek 9.5 zjistíme, že bázové funkce pro rohové uzlové body LSTprvku lze zapsat pomocí prirozené souradnice ξ1 ve tvaru:

N1 = Aξ21 +Bξ1 + C

pomocí okrajových podmínek:

N1(0) = 0 N1(1

2) = 0 N1(1) = 1

získáme vztah pro první bázovou funkci LST prvku: N1(ξ) = 2ξ21 − ξ1. Zbývající bázovéfunkce snadno zjistíme odmenou okrajových podmínek:

N1 = ξ1 (2ξ1 − 1)

N2 = ξ2 (2ξ2 − 1) (9.37)N3 = ξ3 (2ξ3 − 1)


(x1, y1)

(x2, y2)

(x3, y3)

1

N1(ξ)ξ1

241

ξ2

Obrázek 9.6: Kvadratická báz. funkce N4(ξ) pro LST prvek

Bázové funkce pro mezilehlé uzlové body lze rovnež s pomocí obrázku 9.6 odhadnout,napr. funkci N4(ξ) budeme hledat ve tvaru:

N4(ξ) = Aξ1ξ2

kostantu A získáme aplikací podmínky pro jednotkové posunutí v míste uzlu, tedy:

N4(1

2) = 1⇒ Aξ1(

1

2)ξ2(

1

2) = 1

A1

2

1

2= 1⇒ A = 4

zbývající bázové funkce N5a N6 obdobným zpusobem, výsledek je:

N4 = 4ξ1ξ2

N5 = 4ξ2ξ3 (9.38)N6 = 4ξ3ξ1

Matice bázových funkcí pro LST prvek má tedy tvar:

N =ξ1 (2ξ1 − 1) , ξ2 (2ξ2 − 1) , ξ3 (2ξ3 − 1) , 4ξ1ξ2, 4ξ2ξ3, 4ξ3ξ1

(9.39)

Nyní se v rychlosti podívejme na ctyrúhelníkový izoparametrický prvek Q8, definovanýpomocí osmi uzlových bodu. Pomerne jednoduše lze odhadnout bázové funkce v prípadejeho lineární varianty Q4 s uzlovými body pouze ve vrcholech. Bázové funkce jsou velmijednoduché, lze je konstruovat jako soucin dvou funkcí zvlášt’ pro souradnici ξ a η.Lineární funkce ve smeru souradnice ξ je dána pomocí 1

2 (1− ξ) , stejne tak ve smeru ηbude ona lineární funkce 1

2 (1− η). Jejich soucinem získáme první bázovou funkci prolineární ctyrúhelníkový prvek NQ4

1 (ξ, η).

Obdobne získáme zbývající bázové funkce pro lineární prvek Q4, pouze meníme kvad-ranty a tím i príslušné znaménko. Pro prvek Q4 mužeme tedy zapsat následující lineární


funkce:

NQ4

1 (ξ, η) =1

4(1− ξ)(1− η)

NQ4

2 (ξ, η) =1

4(1 + ξ)(1− η) (9.40)

NQ4

3 (ξ, η) =1

4(1 + ξ)(1 + η)

NQ4

4 (ξ, η) =1

4(1− ξ)(1 + η)

Pomerne jednoduchou geometrickou úvahou lze z lineárních bázových funkcí konstru-ovat kvadratické bázové funkce pro prvek Q8. Postup je zrejmý z obrázku 9.7.

η

(x1, y1) (x2, y2)

(x3, y3)(x4, y4)

N5(ξ)

1 ξ

1

η

(x1, y1) (x2, y2)

(x3, y3)(x4, y4)

N8(ξ)

ξ

Obrázek 9.7: Kvadratická báz. funkce N5(ξ) a N5(ξ) pro prvek Q8

Bázovou funkci NQ85 (ξ, η) získáme jako soucin dvou funkcí: ve smeru souradnice ξ kva-

dratické funkce 1− ξ2 a ve smeru souradnice η lineární funkce 12 (1− η):

NQ8

5 (ξ, η) =1

2(1− ξ2)(1− η)

Obdobne pro funkci N8(ξ, η), tentokrát bude však kvadratická funkce ve smeru sourad-nice η:

NQ8

8 (ξ, η) =1

2(1− ξ)(1− η2)

Pro rohové uzly a bázové funkce N1(ξ, η)...N4(ξ, η) lze postupovat od lineární bázovéfunkce NQ4

1 (ξ, η) již definované pro lineární prvek Q4.

NQ8

step1 =1

4(1− ξ) (1− η)

Ve smeru souradnice η nejprve odecteme polovinu funkce NQ8

8 (polovinu musíme odecístproto, abychom pro souradnici η = 0 získali vždy 0):

NQ8

step2 = Nstep1 −1

2NQ8

8 =1

4

[(1− ξ) (1− η)− (1− ξ)

(1− η2

)](9.41)

stejný postup lze uplatnit i pro druhý smer, tedy pro souradnici ξ:

NQ8

1 = Nstep2−1

2N5 = N1−

1

2NQ4

5 −1

2NQ4

8 =1

4

[(1− ξ)(1− η)− (1− ξ)(1− η2)− (1− ξ2)(1− η)

]

(9.42)

Všechny geometrické úvahy jsou zrejmé z obr. 9.8.


η

(x1, y1) (x2, y2)

(x3, y3)(x4, y4)

1

Nstep1

1(ξ)

ξ

Nstep11 = 1

4 (1− ξ) (1− η)

η

(x1, y1) (x2, y2)

(x3, y3)(x4, y4)

1

Nstep2

1(ξ)

ξ

Nstep21 = Nstep1

1 − 12N8

η

(x1, y1) (x2, y2)

(x3, y3)(x4, y4)

1

Nstep3

1(ξ)

ξ

N1 = Nstep21 − 1

2N5

Obrázek 9.8: Postup tvorby kvadratické báz. funkce N1(ξ) pro prvek Q8


9.3.1 Souhrn bázových funkcí pro nejbežnejší rovinné prvky

Záverem si shrnme, jak budou vypadat bázové funkce pro nejbežnejší typy rovinnýchelementu. Prirozené souradnice nyní i pro trojúhelník oznacme pomocí trí písmen reckéabecedy, konkrétne ξ, η a ζ:

9.3.1.1 Constant Strain Triangle (CST, T3)

Tento element je znázornen na obr. 10.5.

v

u

v2

u2

v3

u3

v1

u1

(x2, y2)

(x3, y3)

(x, y)

(x1, y1)

Obrázek 9.9: Trojúhelníkový prvek s konstantní deformací (Constant Strain Triangle)

N1 = ξ

N2 = η (9.43)N3 = ζ

9.3.1.2 Linear Strain Element (LST, T6)


N1 = ξ(2ξ − 1)

N2 = η(2η − 1)

N3 = ζ(2ζ − 1) (9.44)N4 = 4ξη

N5 = 4ηζ

N6 = 4ηξ

9.3.1.3 Linear Quadrilateral Element (Q4)



v1

u1

v3

u3

v6

u6 v

u

v4

u4

v5

u5

v2

u2(x1, y1)

(x2, y2)

(x3, y3)

(x, y)

Obrázek 9.10: Trojúhelníkový prvek s lineární deformací (Linear Strain Triangle)

η

ξ

(x3, y3)

(x4, y4)

u2

v3

u3v4

v1

u4

u1

ξ = −1ξ = +1

v2

η = +1

η = −1(x1, y1)

(x2, y2)

Obrázek 9.11: Ctyrúhelníkový prvek s konstantní deformací (Linear Quadrilateral Ele-ment)


N1(ξ, η) =1

4(1− ξ)(1− η) N2(ξ, η) =

1

4(1 + ξ)(1− η) (9.45)

N3(ξ, η) =1

4(1 + ξ)(1 + η) N4(ξ, η) =

1

4(1− ξ)(1 + η)

9.3.1.4 Quadratic Quadrilateral Element (Q8)


η

ξ

(x3, y3)

(x4, y4)

u2

v3

u3

v1

u1

ξ = −1ξ = +1

v2

η = +1

η = −1(x1, y1)

(x2, y2)

u8

u5

u6

u7

v4

u4

v8

v5

v6

v7

Obrázek 9.12: Ctyrúhelníkový prvek s lineární deformací (Quadratic Quadrilateral Ele-ment)

N1(ξ, η) =1

4(1− ξ)(η − 1)(ξ + η + 1) N2(ξ, η) =

1

4(1 + ξ)(η − 1)(η − ξ + 1) (9.46)

N3(ξ, η) =1

4(1 + ξ)(1 + η)(ξ + η − 1) N4(ξ, η) =

1

4(1− ξ)(1 + η)(ξ − η + 1)

Kapitola 10

Prvky pro rovinnou úlohu

10.1 Základní úlohy teorie pružnosti

Pripomenme si základní rovnice teoretické pružnosti (geometrické, statické a fyzikálnírovnice). V prostoru máme urcit vektor posunutí

u =

uvw

vektor deformace:

ε =

εx

εy

εz

γyz

γxz

γxy

a vektor napetí:

σ =

σx

σy

σz

τyz

τxz

τxy

Prvními rovnice, které si pripomeneme, jsou geometrické rovnice, které vyjadrují vztahmezi vektorem posunutí u a vektorem deformace ε:

εx =∂u

∂xγyz =

∂w

∂y+∂v

∂z

εy =∂v

∂yγxz =

∂w

∂x+∂u

∂z(10.1)

εz =∂w

∂zγxy =

∂v

∂x+∂u

∂y

89

90 KAPITOLA 10. PRVKY PRO ROVINNOU ÚLOHU

Vztah mezi napetím a pretvorením je u lineárního elastického materiálu dán rozšírenýmHookeovým zákonem, jenž je z teorie pružnosti znám ve tvaru:

εx =σx

Ex− νxy

σy

Ey− νxz

σz

Ezεx =

1

E[σx − ν(σy + σz)]

εy =σy

Ey− νyx

σx

Ex− νyz

σz

Ezεy =

1

E[σy − ν(σx + σz)]

εz =σz

Ez− νzx

σx

Ex− νzy

σy

Eyεx =

1

E[σz − ν(σy + σx)] (10.2)

γyz =1

Gyzτyz

γxz =1

Gxzτxz

γyx =1

Gyxτyx

tyto rovnice se dají jednoduše zapsat pomocí maticového zápisu ve tvaru:

ε = [C]σ ε = Cσ (10.3)

kde [C] je matice poddajnosti pro obecne anizotropní materiál. Inverzí vztahu 10.3 zís-káváme vztah pro výpocet vektoru napetí, ve kterém [D] je matice tuhosti materiálu:

σ = [D] ε (10.4)

Platí, že [D]−1

= [C]. Ješte pro úplnost “vycíslíme” matici poddajnosti (izotropní mate-riál):

[C] = [D]−1 =

1Ex

−νxy

Ey

−νxz

Ez 0 0 0−νyx

Ex

1Ey

−νyz

Ez 0 0 0−νzx

Ex

−νzy

Ey

1Ez

0 0 0

0 0 0 1Gyz

0 0

0 0 0 0 1Gxz

0

0 0 0 0 0 1Gxy

(10.5)

Další rovnice jsou Cauchyovy rovnice, plynoucí z podmínky rovnováhy na infinitesimál-ním objemu dV = dx× dy × dz.

∂σx

∂x+∂τxy

∂y+∂τxz

∂z+ bx = 0 (10.6)

∂τxy

∂x+∂σy

∂x+∂τyz

∂z+ by = 0 (10.7)

∂τxz

∂x+∂τyz

∂y+∂σz

∂z+ bz = 0 (10.8)

zkrácene:[∂]σ+ b = 0 (10.9)

10.1. ZÁKLADNÍ ÚLOHY TEORIE PRUŽNOSTI 91

z

x

σy

τyz

τyx

σz

τzy

τxz +∂τyz

∂xdx

σx +∂σx

∂xdx

σy +∂σy

∂ydy

σx

τxy

τzx +∂τzx

∂zdz

τzy +∂τzy

∂zdz

σz +∂σz

∂zdz

y

τyx +∂τyx

∂ydy

τyz +∂τyz

∂ydy

τxy +∂τxy

∂xdx

τxz

τzx

dV = dxdydz

Obrázek 10.1: Podmínky rovnováhy na elementu


kde b je vektor objemových sil (jednotky[

kgm3

]

) a [∂] je tzv. operátorová matice:

[∂] =

∂∂x 0 0 0 ∂

∂z∂∂y

0 ∂∂y 0 ∂

∂z 0 ∂∂x

0 0 ∂∂z

∂∂y

∂∂x 0

(10.10)

Provedením momentových výminek rovnováhy kolem os rovnobežných s osami sourad-nic a procházejících težištem kvádru obdržíme rovnice vyjadrující tzv. vetu o vzájem-nosti smykových napetí:

γxy = γyx γyz = γzy γxz = γzx (10.11)

V prostoru hlavních napetí (Haigh-Westergarduv prostor) bude vektor napetí mít tvar:

σT =σ1 σ2 σ3 0 0 0

(10.12)

Doplnme ješte pro úplnost výrazy pro invarianty napetí (budeme potrebovat pri teoriíchporušení):

J1 = σx + σy + σz = σ1 + σ2 + σ3

J2 = σyσz + σxσz + σxσy − τ2yz − τ2

xz − τ2xy = σ2σ3 + σ1σ3 + σ1σ2 (10.13)

J3 = σxσyσz + 2τyzτxzτxy − σxτ2yz − σyτ

2xz − σzτ

2xy = σ1σ2σ3

a invarianty deformace:

I1 = εx + εy + εz = ε1 + ε2 + ε3

I2 = εyεz + εxεz + εxεy −1

4(γ2

yz + γ2xz + γ2

xy) = ε2ε3 + ε1ε3 + ε1ε2 (10.14)

I3 = εxεyεz +1

4(γyzγxzγxy − εxγ

2yz − εyγ

2xz − εzγ

2xy) = ε1ε2ε3

Pomocí invariantu napetí/deformace lze jednoduše rozložit napjatost v bode na na-pjatost objemovou (volumetrickou) a deviatorickou. Také jsou tyto napjatosti nazýványnapjatost zmeny objemu a napjatost zmeny tvaru. Symbolicky zapíšeme pro vektor na-petí:

σT = σvT + σDTσvT =

J1/3 J1/3 J1/3 0 0 0

(10.15)

σDT =σx − J1/3 σy − J1/3 σz − J1/3 τyz τxz τxy

a vektor deformace:

εT = εvT + σDTεvT =

I1/3 I1/3 I1/3 0 0 0

(10.16)

εDT =εx − I1/3 εy − I1/3 εz − I1/3 τyz τxz τxy

10.2 Úlohy rovinné pružnosti

Velmi používaným zjednodušením problému pružnosti je snížení dimenzionality danéhoproblému. Velmi casto zjednodušení prostorového problému na dvoudimenzionální máza následek dramatické snížení casové nárocnosti na výpocet.

10.2. ÚLOHY ROVINNÉ PRUŽNOSTI 93

10.2.1 Rovinná deformace

Nejsrozumitelnejším príkladem rovinné deformace je dlouhá stena, napr. operná zed’z obrázku 10.2. Geometrie steny ani její zatížení se nemení podél její délky a snadnozjistíme, že všechny námi zkoumané veliciny jsou závislé pouze na souradnicích x a y.

Obrázek 10.2: Rovinná deformace - operná zed’

Podíváme-li se na libovolný prícný rez, zjistíme, že posunutí ve smeru osy z jsou nulováa tudíž i deformace εz, γyz, γzx budou nulové. Vektor deformace se tímto zredukuje abude mít pouze tri nenulové složky:

εx =∂u

∂xεy =

∂v

∂yγxy =

∂u

∂y+∂v

∂x(10.17)

Po dosazení εz = 0 do rozšíreného Hookova zákona:

εz =1

E[σx − ν (σy + σz)] = 0 (10.18)

dostáváme:σz = ν(σx + σy) (10.19)

To znamená, že napetí ve smeru osy z je na rozdíl od deformace v tomto smeru nenulové,dá se však spocítat z hodnot normálových napetí ve dvou smerech roviny. Ostatní složkynapetí, se dají vyjádrit z rozšíreného Hookova zákona dosazením σz = ν(σx + σy) dovýrazu pro εx a εy:

εx =1

Eσx − ν [σy + ν(σx + σy)] =

1

E

[σx − νσy − ν2σx − ν2σy

]=

=1

E

[σx

(1− ν2

)− νσy (1 + ν)

]=

1− ν2

E

(

σx −ν

1− ν σy

)

εy = · · · =1− ν2

E

(

− ν

1− ν σx + σy

)

Pro rovinnou deformaci dostáváme tedy matici materiálové poddajnosti:

C =

1−ν2

E − ν(ν+1)E 0

− ν(ν+1)E

1−ν2

E 00 0 1

G

(10.20)


Inverzí matice materiálové poddajnosti [C] získáme matici materiálové tuhosti [D]. Tremipríkazy pomocí Maximy:

(C1) C: matrix([(1-nu^2)/E,-nu*(1+nu)/E,0],[-nu*(1+nu)/E,(1-nu^2)/E,0],[0,0,1/G]);

(C2) D: factor(invert(C));

(C3) D_mod: ((nu+1)*(2*nu-1))*D;

získáme

[D] =

(ν−1)E(ν+1)(2ν−1)

νE(ν+1)(2ν−1) 0

νE(ν+1)(2ν−1)

(ν−1)E(ν+1)(2ν−1) 0

0 0 G

=

1

(ν + 1) (2ν − 1)

(ν − 1)E νE 0νE (ν − 1)E 00 0 (ν + 1) (2ν − 1)G

Do této rovnice ješte dosad’me vztah mezi normálovým modulem pružnosti E a modu-lem pružnosti ve smyku, který je dán G = E/2(1+ ν). Pro rovinnou deformaci (a lineárníizotropní materiál) tedy platí tyto fyzikální rovnice:

σx

σy

τxy

=

E

(1 + ν)(1− 2ν)

1− ν ν 0ν 1− ν 00 0 1−2ν

2

εx

εy

γxy

(10.21)

10.2.2 Rovinná napjatost

Na rozdíl od rovinné deformace, kde jsme meli konstrukci, její jeden rozmer byl v po-rovnání s ostatními rozmery velký, rovinná napjatost zahrnuje konstrukce, které sevyznacují jedním rozmerem zanedbatelným v porovnání se zbylými rozmery. Príklademtakové konstrukce jsou desky. Tenká deska zatížená pouze ve své strednicové rovine zobr. 10.3 je nejcastejší ukázka rovinné napjatosti.

σz

= 0

Obrázek 10.3: Rovinná napjatost - tenká deska zatížená ve své strednicové rovine

Uvažujme prípad, kdy deska není zatížena zatížením pusobícím kolmo k její strednicovérovine. Složky vektory napetí τyz a τzx budou nulové na povrchu desky a složka σz budenulová po celé tloušt’ce desky. Potom zbývající složky vektoru napetí σx, σy, τxy mohou


být brány po výšce jako prumerné hodnoty zcela nezávislé na souradnici z. Složky γyz,γzx jsou nulové, zatímco εz se dá vyjádrit pomocí zbývajících nenulových složek vektorudeformace:

εz = − ν

1− ν (εx + εy) (10.22)

vektor deformace bude shodný s vektorem deformace z prípadu rovinné deformace(pouze si musíme uvedomit, že deformace ve smeru osy z není nulová):

εx =∂u

∂xεy =

∂v

∂yγxy =

∂u

∂y+∂v

∂x(10.23)

Do rozšíreného Hookova zákona dosad’me σz = 0. Matice poddajnosti bude mít tvar:

C =

1E − ν

E 0− ν

E1E 0

0 0 1G

(10.24)

Pomocí Maximy nalezneme opet matici materiálové tuhosti [D] inverzí matice materiá-lové poddajnosti [C]:(C1) C: matrix([1/E,-nu/E,0],[-nu/E,1/E,0],[0,0,1/G]);

D =

− E(ν−1)(ν+1) − νE

(ν−1)(ν+1) 0

− νE(ν−1)(ν+1) − E

(ν−1)(ν+1) 0

0 0 G

(C2) D: factor(invert(C));

vytkneme výraz − 1(ν−1)(ν+1) pred matici:

(C3) D_mod: -1*(nu-1)(nu+1)*D;

D = − 1

(ν − 1) (ν + 1)

E νE 0νE E 00 0 − (ν − 1) (ν + 1)G

(10.25)

Konstitutivní rovnice pro rovinnou deformaci jsou:

σx

σy

τxy

=

E

1− ν2

1 ν 0ν 1 00 0 1−ν

2

εx

εy

γxy

(10.26)

10.2.3 Osová soumernost

Znacná cást problému teoretické mechaniky zahrnuje analýzu osove soumerných prvku.Jestliže navíc zatížení takového prvku je rovnež symetrické podél stejné osy, lze s výho-dou využít této symetrie. Príkladem takovéto situace muže být napríklad válcová pod-pera mostu, na níž pusobí pouze hydrostatický tlak okolní vody.Je samozrejme velmi výhodné tento problém uvažovat ve válcovém souradném systému,zobrazeném na obrázku. Z duvodu symetrie, složky vektoru napetí budou nezávislé nasouradnici φ a navíc všechny derivace vzhledem k φ vymizí. Z tohoto plyne nulovostsložek v, γrφ, γφz, τrφ, τφz. Geometrické rovnice pro nenulové složky vektoru deformacebudou mít tvar:

εr =∂u

∂rεφ =

u

rεz =

∂w

∂zγrz =

∂u

∂z+∂w

∂r(10.27)


r, u

z, w

Θr, u

Obrázek 10.4: Osová soumernost

Konstitutivní rovnice:

σr

σz

σφ

τrz

=

1− ν ν ν 01− ν ν 0

1− ν 0sym 1−ν

2

εr

εz

εφ

γrz

(10.28)

10.2.4 Prvek pro rovinnou deformaci

Posunutí (u, v) jsou na rovinném elementu interpolována ze známých uzlových posunutí(ui, vi) pomocí bázových funkcí Ni

uv

=

[N1 0 N2 0 . . .0 N1 0 N2 . . .

]

u1

v1u2

v2...

u = Nr (10.29)

kde matice N je matice bázových funkcí, u vektor posunutí a r vektor uzlových posunutí.Vektor posunutí (po prvku) je pouze funkcí uzlových posunutí, ríkáme, že je pomocíbázových funkcí aproximován z uzlových posunutí.

Pripomenme, že pri rovinné deformaci platí εz = γyz = γxz = 0 σz 6= 0 = ν(σx + σy).Geometrické rovnice nabydou tvaru

εx

εy

γxy

=

∂/∂x 00 ∂/∂y

∂/∂y ∂/∂x

uv

ε = ∂Nr = Bu (10.30)


kde B = ∂N je geometrická matice (STRAIN-DISPLACEMENT MATRIX) a ∂ je operátorovámatice. Nyní odvodíme obecný vzorec pro matici tuhosti elementu z energetického prin-cipu. Deformacní energie v elementu se dá vyjádrit jako:

U =1

2

∫

(V )

σT ε dV =1

2

∫

(V )

(σxεx + σyεy + τxyγxy) dV

=1

2

∫

(V )

(Dε)T ε dV =1

2

∫

(V )

εTDε dV =

1

2

∫

(V )

uTB

TDBu dV (10.31)

=1

2u

T

∫

(V )

BTDB dV u =

1

2u

Tku

Obecný výraz pro matici tuhosti elementu má tedy tvar:

k =

∫

(V )

BTDB dV (10.32)

Matice D (pouze v prípade 1-D napjatosti je to skalár) se nazývá matice materiálovétuhosti a v prípade rovinné defomace, bude mít tvar 10.33. Pripomenme si ovšem, žetento vztah platí pouze pro prípad lineárního homogenního izotropního materiálu.

D =E

(1 + ν)(1− 2ν)

1− ν ν 0ν 1− ν 00 0 1−2ν

2

(10.33)

Tento vztah nejlépe odvodíme inverzí rozšíreného Hookova zákona s uvážením pocátec-ních deformací (napr. od teploty)

εx

εx

γxy

=

1/E −ν/E 0−ν/E 1/E 0

0 0 1/G

σx

σy

τxy

+

εx0

εy0

γxy0

ε = Cσ (10.34)

kde

C =

1/E −ν/E 0−ν/E 1/E 0

0 0 1/G

se nazývá matice materiálové poddajnosti a samozrejme platí C = D−1.

Fyzikální rovnice s uvážením vlivu pocátecních deformací (napr. opet vliv teploty) budoumít tvar

σ =E

(1 + ν)(1− 2ν)

1− ν ν 0ν 1− ν 00 0 1−2ν

2

εx

εy

γxy

−

εx0

εy0

γxy0

σ = D(ε− ε0)

(10.35)

vektor pocátecních deformací vlivem zmeny teploty bude

εx0

εy0

γxy0

=

α∆tα∆t0

(10.36)


10.2.5 Prvek pro rovinnou napjatost

Tento prvek bude v zásade stejný jako prvek pro rovinnou deformaci, jen matice ma-teriálové tuhosti bude mít jiný tvar. Fyzikální rovnice dostaneme snadno nahrazenímvztahu pro Younguv modul pružnosti v tahu/tlaku a ve smyku dle

E′

= E(1− ν2)

ν′

= ν(1− ν) (10.37)G

′

= G

Matice materiálové tuhosti pro prípad rovinné napjatosti bude tedy

D =E

1− ν2

1 ν 0ν 1 00 0 1−ν

2

(10.38)

Pro prehlednost ješte doplnme fyzikální rovnice

σx

σy

τxy

=

E

1− ν2

1 ν 0ν 1 00 0 1−ν

2

εx

εy

γxy

−

εx0

εy0

γxy0

σ = D(ε− ε0) (10.39)

10.2.6 Statické podmínky rovnováhy (Cauchyho rovnice)

Na zacátku této kapitoly jsme pripomneli sestavení podmínek rovnováhy ve 3-D naelementárním objemu dV = dxdydz. Uvážíme-li objemové síly v rovine (síly na jednotkuplochy, napr. tíhové síly) budou mít tyto tzv. Cauchyho rovnice tvar (jsme v rovine):

∂σx

∂x+∂τxy

∂y+ bx = 0 (10.40)

∂τxy

∂x+∂σy

∂y+ by = 0

10.2.7 Okrajové podmínky

Okrajové podmínky jsou specifikovány jako vnejší síly a predepsané posunutí. Posunutína hranici mohou být jak nulová (homogenní okrajové podmínky), tak mít stanovenounenulovou hodnotu. V metode konecných prvku jsou všechny typy zatížení (spojité si-lové zatížení, osamelé síly, spojité momentové zatížení ci osamelý moment) prevádenana jednotlivé síly pusobící v uzlových bodech.

Okrajové podmínky se v prípade rovinného problému (u = u, v)dají symbolicky zapsatjako:

u = u v = v (10.41)tx = tx ty = ty

kde okrajové silové podmínky znacíme tx, ty z duvodu jejich názvu – temto silám se ríkáTRACTION FORCES a jsou to vlastne hodnoty napetí na hranici, t.j. napr. tx = σxnx, kdenx znací normálu.


10.2.8 Odvození matice tuhosti pro CST-prvek a rovinnou úlohu

Tento element je znázornen na obr. 10.5. Prvek je definován vektorem uzlových posu-nutí:

rT =u1 v1 u2 v2 u3 v3

(10.42)

v

u

v2

u2

v3

u3

v1

u1

(x2, y2)

(x3, y3)

(x, y)

(x1, y1)

Obrázek 10.5: trojúhelníkový prvek konstatní deformace (Constant Strain Triangle)

Posuny uvnitr prvku jsou aproximovány pomocí posunutí uzlových bodu a bázovýchfunkcí ve tvaru:

N1

N2

N3

=

1

2A

x2y3 − x3y2 y2 − y3 x3 − x2

x3y1 − x1y3 y3 − y1 x1 − x3

x1y2 − x2y1 y1 − y2 x2 − x1

1xy

(10.43)

N1 +N2 +N3 = 1 Ni =

1 v uzlu i0 v ostatnich uzlech

Geometrické rovnice lze zapsat pomocí složek vektoru deformace:

εx =∂u

∂x= ∂Nr =

3∑

i=1

∂Ni

∂xui

εy =∂v

∂y= ∂Nr =

3∑

i=1

∂Ni

∂yvi

γxy =∂v

∂x+∂u

∂x=

3∑

i=1

(∂Ni

∂xvi +

∂Ni

∂yui

)

Což lze symbolicky zapsat pomocí následujícího maticového zápisu (tyto rovnice již byly


pro CST prvek odvozeny v predešlé kapitole, není tedy treba se k nim detailneji vracet):

εx

εy

γxy

=

1

2A

y23 0 y31 0 y12 00 x32 0 x13 0 x21

x32 y23 x31 y31 x21 y12

u1

v1u2

v2u3

v3

(10.44)

Z tohoto zápisu je zrejmé, že všechny tri složky vektoru deformace jsou po prvku kon-stantní (odtud název prvku). Jakmile máme vztah mezi vektorem deformace ε a kon-cových posunutí r, máme matici B a mužeme vyjádrit matici tuhosti CST prvku inte-grací:

K =

∫

Ω

BTDB dΩ

Kapitola 11

Deskové prvky

11.1 Predpoklady, deskové teorie

Deska je konstrukce, jejíž dva rozmery jsou výrazne vetší než zbývající tretí rozmer. Dlepomeru tloušt’ky desky k jejím zbývajícím rozmerum t

a , tb rozdelujeme desky na:

• tlusté desky: 15 − 1

10

– s vlivem smyku– Reissner-Mindlinova teorie– Timoshenkova teorie ohýb. prutu

• tenké desky: 110 − 1

50

– bez vlivu smyku– Kirchhoffova teorie– Euler-Bernoulli teorie ohýb. prutu

• velmi tenké desky: < 150

– geometrická nelinearita– von Karmánova teorie– teorie druhého rádu (podmínky rovnováhy na pretvoreném prutu)

11.2 Kirchhoffova teorie tenkých desek

Predpoklady:

1. geometrická linearita: malé pruhyby a malé deformace

2. lineární materiál: platnost Hookova zákona σ = Eε

3. tenká deska: ta,b ∈

⟨15 ,

150

⟩

Jsou-li splneny uvedené predpoklady, lze výpocty založit na tzv. Bernoulliho hypotéze:

101

102 KAPITOLA 11. DESKOVÉ PRVKY

• normály se nedeformují (neprodlužují)

εz = 0

• zanedbání Poissonova efektuσz = 0

11.2.1 Kinematické rovnice

strednicová rovina:

• pruhybyw = w(x, y)

• pootocení materiálového bodu

φx =∂w

∂y

φy =∂w

∂x

• posuny bodu P (x, y, z):

u = −z ∂w∂x

= −zφy

v = −z ∂w∂y

= −zφx

w

11.2.2 Geometrické rovnice

Vztahy mezi vektorem posunutí a vektorem deformace.

εx =∂u

∂x= −z ∂

2w

∂x2= −zκx

εy =∂v

∂y= −z ∂

2w

∂y2= −zκy

εz =∂w

∂z= 0

γyz =∂w

∂y+∂v

∂z= 0− z ∂

2w

∂y∂z= 0

γxz =∂w

∂x+∂u

∂z= 0− z ∂

2w

∂x∂z= 0

γxy =∂v

∂x+∂u

∂y= −z ∂

2w

∂y∂x− z ∂

2w

∂x∂y= −2zκxy

vektor deformace bude mít tedy pouze tri nenulové složky:

ε =

εx

εy

γxy

=

−zκx

−zκy

−2zκxy

(11.1)

11.2. KIRCHHOFFOVA TEORIE TENKÝCH DESEK 103

11.2.3 Statické rovnice

Zde ve tvaru moment~krivost.Predpokládejme opet lineární elastický materiál, platnost Hookova zákona. Fyzikálnírovnice mužeme psát ve tvaru:

σx

σy

τxy

=

D11 D12 D13

D22 D23

sym D33

εx

εy

γxy

= −z [D]

κx

κy

2κxy

Dosad’me nyní tyto dyzikální rovnice do podmínek rovnováhy pro vnitrní síly. Ohybovýmoment na jednotku délky získáme snadno jako výslednice príslušného napetí:

mxdy =

∫ h/2

−h/2

σx(−z)dydz ⇒ mx = −∫ h/2

−h/2

zσxdz

Obdobne i ostatní složky ohybových momentu:

mx = −∫ h/2

−h/2

zσxdz

my = −∫ h/2

−h/2

zσydz

mx = −∫ h/2

−h/2

zτxydz

Tyto rovnice prepišme do maticové formy a zároven dosad’me za napetí σ = −zDκ:

mx

my

mxy

= −

h/2∫

−h/2

z(−z)

D11 D12 D13

D22 D23

sym D33

κx

κy

2κxy

dz =

=

h/2∫

−h/2

z2dz

D11 D12 D13

D22 D23

sym D33

κx

κy

2κxy

=

[z3

3

]h2

−h2

[D] κ =

=Eh3

12 (1− ν2)

1 ν 0ν 1 00 0 1−ν

2

κx

κy

2κxy

výraz Eh3

12(1−ν2) se nazývá desková tuhost. Celá matice Eh3

12(1−ν2)

1 ν 0ν 1 00 0 1−ν

2

je pak ma-

tice tuhosti izotropního materiálu pri rovinné napjatosti. Známe-li momenty mx, my amxy snadno dopocteme napetí z podmínek rovnováhy:

mx =1

2σxh

2

2

3h =

1

6σxh

2 ⇒ σx =6mx

h2

obdobne i pro zbývající napetí:

σx = ±6mx

h2

σy = ± 6my

h2

τxy = ±6mxy

h2


Nyní proved’me na infinitesimálním elementu podmínky rovnováhy. Nejprve ve svislémsmeru:

−txdy − tydx+

(

tx +∂tx∂x

dx

)

dy +

(

ty +∂ty∂y

dy

)

dx+ q0dxdy = 0

po odectení nekterých clenu a vydelení rovnice dA = dxdy vychází:

∂tx∂x

+∂ty∂y

+ q0 = 0 (11.2)

Obdobným zpusobem z momentových podmínek rovnováhy získáme:

∂mx

∂x+∂mxy

∂x+ tx = 0

∂mxy

∂x+∂my

∂y+ ty = 0 (11.3)

Z rovnice 11.3 vyjádreme posouvající síly tx a ty:

tx = −∂mx

∂x− ∂mxy

∂y

ty = −∂mxy

∂x− ∂my

∂y

a dosad’me do 11.2:

−∂2mx

∂x2− ∂2mxy

∂x∂y− ∂2my

∂y2− ∂2mxy

∂x∂y+ q0 = 0

Úpravou získáme duležitý vztah:

∂2mx

∂x2+ 2

∂2mxy

∂x∂y+∂2my

∂y2= q0 (11.4)

11.2.4 Biharmonická rovnice desky

Zapišme symbolicky geometrické rovnice:

κ = Pw

kinematické rovnice:M = Dκ

a konecne statické rovnice:P TM = q

V techto rovnicích je P symbolický operátorový vektor obsahující druhé derivace:

P T =

∂2

∂x2,∂2

∂y2,∂2

∂y2

vektor κ je vektor deskových krivostí:

κT = κx, κy, κxy

11.2. KIRCHHOFFOVA TEORIE TENKÝCH DESEK 105

a vektor momentu M :MT = Mx, My, Mxy

Upravme nyní pomocí techto výrazu statické rovnice:

q = P TM = P TDκ = P TDPw = ∇2D∇2w = D∇4w

kde ∇4 je biharmonický operátor definovaný jako soucet parciálních derivací:

∇4 =∂4

∂x4+ 2

∂4

∂x2∂y2+

∂4

∂y4

Získali jsme tedy velmi duležitý vztah pro desku, který se nazývá desková rovnice:

D∇4w = q (11.5)

Je vhodné si ihned všimnout podobnosti s úplnou diferenciální rovnicí ohybové cáryu Euler-Bernoulli teorie ohýbaného prutu: EIwIV = q a srovnat výraz pro ohybovoutuhost prutu EI s deskovou tuhostí D.

11.2.5 Funkcionál potenciální energie

Π(w) = U(w) −W (w)

U(w) =1

2

∫

(Ω)

MTκ dΩ =1

2

∫

(Ω)

κTDκdΩ =1

2

∫

(Ω)

wTP TDPw dΩ

W (w) =

∫

(Ω)

qw dΩ +

∫

(Γ)

tw − mφ dΓ

11.2.6 Hellinger-Reissneruv princip

funkcionál

Π(w,M) = U(w,M)−W (w,M)

deformacní energie:

U(w,M) =

∫

(Ω)

(

MTκ− 1

2MTD−1M

)

dΩ

11.2.7 Výpocet matice tuhosti (doplnit pozdeji)

Trojúhelníkový prvek DKT (Discrete Kirchhoff Theory)

Ve funkcionálu potenciální energie se objevují nejvýše první deriva hledaných funkcí.Proto je v tomto prípade na hledané aproximace funkcí w, φx a φy kladena pouze pod-mínka C0 spojitosti (spojitost funkcních hodnot). cást energie, která prísluší posouva-jícím silám se pri konstrukci funkcionálu zanedbá a v diskrétních bodech se zavede


požadavek zachování normál po deformaci. Tento predpoklad je správný, nebot’ i ex-perimenty prokázaly, že pro tenké desky je energie príslušná smyku zanedbatelná vesrovnání s ohybovou energií.vekor uzlových posunutí:

r =w1 φx1 φy1 w2 φx2 φy2 w3 φx3 φy3

T

z podmínek nulové smykové deformace platí následující geometrické vztahy:

φx =∂w

∂y

φy =∂w

∂x

pootocení φx a φy jsou aproximována pomocí bázových funkcí a uzlových pootocení:

φx = N1φx1 +N2φx2 +N3φx3

φy = N1φy1 +N2φy2 +N3φy3

protože se jedná o lineární trojúhelníkový prvek, budou aproximacní funkce Ni rovnyprímo trojúhelníkovým souradnicím odvozeným v 9.1.2, tedy:

Ni = ξi =Ai

A

Aby byla splnena Kirchhoffova hypotéza, pro aproximaci funkcí ∂w∂s a φn se volí kvadra-

tické polynomy, které jsou ztotožneny ve trech bodech každé strany. Je zrejmé, že nahranicích mezi prvky je zachována spojitost funkcí w, ∂w

∂s , φs i φn.

φx

φy

=

[cosα − sinαsinα cosα

]φn

φs

(11.6)

∂w∂s∂w∂n

=

[cosα − sinαsinα cosα

]φx

φy

(11.7)

Aproximace pootocení φx a φy se pak zapíše pomocí aproximacních funkcí Nx1, Ny1,...,Nx6, Ny6:

Nx1 = −1.5 (d6N6 − d5N5)

Nx2 = N1 − e5N5 − e6N6

Nx3 = b5N5 + b6N6

......

...

Matice tuhosti se pak vypocte dle:

K =

∫

(A)

BTDoB dA

kde Do matice tuhosti izotropního materiálu pri rovinné napjatosti daná pomocí:

Do =Eh3

12 (1− ν2)

1 ν 0ν 1 00 0 1−ν

2

11.3. MINDLINOVA TEORIE TLUSTÝCH DESEK 107

11.3 Mindlinova teorie tlustých desek


Kapitola 12

Úvod do dynamiky v MKP

Až dosud jsme se zabývali analýzou konstrukcí zatížených staticky, t.j. zatížení se me-nilo velmi pomalu, nebylo tedy treba brát v úvahu vliv setrvacných sil, které jsou du-sledkem zrychlení konstrukce, nebo její cásti. Dynamika se zabývá vlivem setrvacnýchsil, kdy zatížení pusobící na konstrukci zpusobují její rozkmitání. Temito zatíženími bý-vají nejcasteji

• pohyblivé zatížení

• rotující stroje

• náraz hmotných predmetu

• úcinky vetru

• zemetresení

• úcinky exploze

Úlohy dynamiky lze v zásade rozdelit do dvou hlavních kategorií: volné kmitání a vynu-cené kmitání. Ze základního kurzu dynamiky (MEC2) si pamatujeme, že oproti statickéanalýze nám do výpoctu vstupuje další duležitá velicina - hmotnost (a s ní spojené se-trvacné síly, které se objeví v pseudostatických podmínkách rovnováhy). V MKP totoznamená, že v maticovém zápisu podmínek rovnováhy se nám objeví další matice - ma-tice hmotnosti. Tu budeme uvažovat ve dvou základních variantách, jako konzistentnímatici hmotnosti (CONSISTENT MASS MATRIX) a nekonzistentí matici hmotnosti (LUMPEDMASS MATRIX). První z nich lze odvodit pomocí variacního principu, v dynamice použí-váme Hamiltonuv princip, což je variacní princip zohlednující historii posunutí.

12.1 Variacní princip v dynamice - Hamiltonuv princip

Jak už bylo receno, v dynamice používáme jako variacní princip Hamiltonuv, kterývyužívá funkcionál L, tzv. Lagrangián, definovaný pomocí:

L = K − U −Wp (12.1)

kde K je kinetická energie telesa, U je jeho deformacní energie a Wp je práce vnejšíchsil. Výraz pro kinetickou energii zavedeme pomocí její hustoty dK:

dK =1

2ρuT u dV (12.2)

109

110 KAPITOLA 12. ÚVOD DO DYNAMIKY V MKP

kde ρ je objemová hustota a tecka nad u znací derivaci dle casu, v tomto prípade tedyrychlost. Práci vnejších sil zapíšemem pomocí jejího prírustku:

dWp = −uTX dV − uTT dS (12.3)

kde X jsou objemové síly, t.j. síly na jednotku objemu (príkladem muže být tíhová síla), aT jsou síly pusobící na hranici, tzv. TRACTION FORCES. Deformacní energie se dá vyjádritpomocí materiálové matice tuhosti [D], nebo materiálové matice poddajnosti [C]:

dU =1

2σT ε dV =

1

2σT [C]σ dV (12.4)

=1

2εT [D]ε dV

Funkcionál z rovnice 12.1 pro lineární elastický materiál bude mít podobu:

L =1

2

∫

(V )

(ρuTu − εT [D]ε+ 2uTX

)dV +

∫

(S)

uTT dS (12.5)

Hamiltonuv princip lze slovne vyjádrit takto:Ze všech možných historií posunutí mezi dvemi casovými okamžiky t1 a t2, které vyhovujípodmínkám kompatibility a kinematickým okrajovým podmínkám, nastane ta z nich, prokterou nabývá Lagrangián minimální hodnoty, t.j. integrál rozdílu kinetické a potenciálníenergie je stacionární.

δ

∫ t2

t1

Ldt = 0 (12.6)

Dosadíme-li do 12.1 výrazy pro vektor posunutí u = [N ]r a vektor deformace ε =[∂]u = [∂][N ]r = [B]r, získáme:

L =1

2

∫

(V )

(rT [B]T [D][B]r − ρrT [N ]T [N ]r − 2rTX

)dV −

∫

(S)

rT [N ]T T dS(12.7)

Použijeme-li variacního principu z 12.6, obdržíme:∫ t2

t1

(

δrT∫

(V )

[B]T [D][B] dV r − δrT∫

(V )

ρ[N ]T [N ] dV r− (12.8)

−δrT∫

(V )

[N ]T X dV − δrT∫

(S)

[N ]T T dS)

dt = 0

nyní integrujme druhý clen podle casu pomocí per partes:∫ t2

t1

δrT∫

(V )

ρ[N ]T [N ] dV r dt =

[

δrT∫

(V )

ρ[N ]T [N ] dV r]t2

t1

− (12.9)

−∫ t2

t1

δrT∫

(V )

ρ[N ]T [N ] dV r dt

Podle Hamiltonova principu musí být platit: δr(t1) = δr(t2) = 0, to znamená, žeprvní clen za rovnítkem ve výrazu 12.9 musí být roven nule. Dosadíme-li zbývající do12.8 a vytkneme-li δrT , obdržíme:

∫ t2

t1

δrT(∫

(V )

ρ[N ]T [N ] dV r+

∫

(V )

[B]T [D][B] dV r− (12.10)

−∫

(V )

[N ]T X dV −∫

(S)

[N ]T T dS)

dt = 0

12.2. MATICE HMOTNOSTI 111

Protože posunutí δr jsou virtuální (nenulová), musí celý výraz ve velké závorce býtnula. Tímto dostáváme pohybové rovnice pro element:

[m]r] + [k]r] = f (12.11)

kde [m] je konzistentní matice hmotnosti elementu definovaná pomocí:

[m] =

∫

(V )

ρ[N ]T [N ] dV (12.12)

[k] je matice tuhosti elementu a f je vektor uzlového zatížení.

12.2 Matice hmotnosti

V predchozím odstavci jsme pomocí Hamiltonova principu odvodili pohybové rovnice proelement, pricemž jsme odvodili také výraz pro konzistentní matici hmotnosti elementu.V MKP se však casto používá nekonzistentní matice hmotnosti a to zejména z duvoduúspory pameti. Nekonzistentní matice hmotnosti je totiž diagonální a krome zminova-ných úspor pameti také výrazne snižuje nárocnost výpoctu (násobení diagonální maticíje snadné).

Uved’me nejprve príklad výpoctu konzistentní matice hmostnosti pro jednoduchý prípadtaženého-tlaceného prvku. Pripomenme, že bázové funkce pro tento prvek mají tvar:N1 = 1− ξ a N2 = ξ. Konzistentní matice hmotnosti bude tedy:

m =

∫ 1

−1

ρ

1− ξξ

1− ξ ξAl dξ = ρAL

∫ 1

−1

[

(1− ξ)2 ξ − ξ2ξ − ξ2 ξ2

]

dξ

= ρAL

[13

16

16

13

]

=

[m3

m6

m6

m3

]

Pro stanovení nekonzistentní matice hmotnosti použijeme následujícího postupu. Hmot-nost prvku soustredíme do jeho uzlových bodu. Toto lze provést pouze u jednoduchýchprvku - napr. u prutového prvku. Pro náš prípad taženého-tlaceného prvku je to velmisnadné - polovinu celkové hmotnosti prvku pridelíme do každého uzlu:

[m] =

[ρAl2 0

0 ρAl2

]

=

[m2 00 m

2

]

U složitejších prvku toto nelze jednoduše provést a proto se používá tzv. metody kon-denzace - vypocítáme konzistentní matici hmotnosti a ponechají se pouze její diagonálnícleny, ostatní se nahradí nulami. Tímto postupem ovšem není zachována celková hmot-nost prvku a proto je treba diagonální prvky prenásobit takovým koeficientem, kterýnám celkovou hmotnost zarucí. Tento koeficient bude:

ξ = ∆Mel

∑ni=1mii

kde ∆ je dimenzionalita problému (∆=1 pro jednorozmerný prvek, ∆=2 pro dvourozmernýa ∆=3 pro prostorový prvek, Mel je hmotnost elementu a mii jsou diagonální prvky ma-tice hmotnosti príslušející pouze posuvným stupnum volnosti.

Uved’me jednoduchý príklad prutového elementu se ctyrmi stupni volnosti - dva rotacnía dva posuvné v každém uzlovém bode. Do každého uzlu máme soustredit polovicní


hmotnost prvku, t.j. ρAl/2. Nekonzistentní matice hmotnosti prvku tedy bude (ve smerurotacních stupnu volnosti je 0):

[m] =ρAl

2

1 0 0 00 0 0 00 0 1 00 0 0 0

Doplnme výraz pro konzistentní matici hmotnosti pro tentýž prvek (lze snadno odvoditintegrací m =

∫

V ρNTNdV známých bázových funkcí Ni pro i = 1, 2, 3, 4). Nezapomenme

na Jacobián transformace J = dxdξ = 1

2 l:

[m] = ρA

∫ 1

−1

NTNl

2dξ =

ρAl

420

156 22l 54 −13l22l 4l2 13l −3l2

54 13l 156 −22l−13l −3l2 −22l 4l2

Tato matice hmotnosti zachovává translacní i rotacní setrvacnost (stupne volnosti naprvku jsou dány vektorem koncových neznámých: r =

vi φi vj φj

T . Sestavení di-agonální matice hmotnosti je zrejmé pro clenym11 a m33, tedy cleny u posuvných stupnuvolnosti. V každém uzlu musí být ve smeru translacního stupne volnosti polovicní hmot-nost: 1

2ρAl (toto vychází samozrejme i souctem “translacních” prvku konzistentní maticehmotnosti - m11a m13 - dostáváme tak ρAl 156+54

420 = 12ρAl). Pro prvky ve smeru rotacních

stupnu volnosti se v prípade tohoto prvku obvykle volí nezáporná hodnota ve tvaruαml2, kde α∈

⟨0; 1

50

⟩. (Moment setrvacnosti prutu ke koncovému, resp. strednímu bodu

je 13ml

2, resp. 112ml

2). Nekonzistentní matice hmotnosti ohýbaného prutu tedy bude:

[m] = ρAl

12 0 0 00 αl2 0 00 0 1

2 00 0 0 αl2

Ukažme si ješte rozdíl mezi nekonzistentní a konzistentní maticí hmotnosti na prípadutrojúhelníkového prvku konstatní deformace (CST). Nekonzistentní matice hmotnostiprvku CST:

[m] =ρhA

3

1 0 0 0 0 00 1 0 0 0 00 0 1 0 0 00 0 0 1 0 00 0 0 0 1 00 0 0 0 0 1

a jeho konzistentní matice hmotnosti:

[m] = ρh

∫

(A)

[N ]T [N ] dA =ρhA

12

2 1 1 0 0 01 2 1 0 0 01 1 2 0 0 00 0 0 2 1 10 0 0 1 2 10 0 0 1 1 2

Pro prvky vyššího rádu je postup sestavení nekonzistentní matice hmotnosti obdobný,uved’me si príklad prutového elementu s mezilehlým uzlem. V každém uzlu je defino-váno pouze jedno uzlové posunutí: r =

u1 u2 u3

T . Stupne volnosti jsou sestaveny

12.2. MATICE HMOTNOSTI 113

tak, že nejprve jsou ocíslovány krajní uzly 1, 2 a nakonec prostrední uzel 3. Pomocíbázových funkcí N1 (ξ), N2 (ξ), N3 (ξ) snadno odvodíme konzistentní matici hmotnosti:

[m] = ρA

∫ 1

−1

NTNl

2dξ =

ρAl

30

4 −1 2−1 4 22 2 16

Chceme-li odvodit nekonzistentní matici hmotnosti pro tento prvek, logicky priradímemenší cást hmotnosti prvku do koncových bodu a vetší cást do prostredního uzlu. Vpredmetu “kinematika a dynamika" se podobným zpusobem prevádel prostý nosníkna kmitání s jedním zpusobem volnosti: do krajních bodu se umístila šestina celkovéhmotnosti prutu 1

6ρAl a strední cásti se priradily zbývající dve tretiny hmotnosti prutu23ρAl. Stejného výsledku dosáhneme aplikací postupu prevedení konzistentní maticehmotnosti na nekonzistentní. Secteme postupne všechny cleny v jednotlivých rádcích avýsledek umísteme na diagonálu:

[m] = ρAl

4−1+230 0 00 −1+4+2

30 00 0 2+2+16

30

=

ρAl6 0 0

0 ρAl6 0

0 0 2ρAl3

=

m6 0 00 m

6 00 0 2

3m

12.2.1 Globalizace matice hmotnosti

Podobne jako tomu bylo u matice tuhosti, matice hmotnosti jsou nejprve sestaveny vlokálním souradném systému definovaným pro element a teprve pri tzv. globalizaci jeúcinek matic hmotnosti jednotlivých elementu promítnut do globální matice hmotnosti(matice hmotnosti konstrukce). Globalizace se dá provést opet pomocí transformacnímatice:

M = T TMT (12.13)kde M je matice hmotnosti v lokálním souradném systému, T je transformacní maticeodvozená v a M je matice hmotnosti elementu v globálním souradném systému.Transformacní matice T je v principu transformacní matice odvozená v kapitole 5.6, snekolika vyjímkami. Za prvé tranasformacní matice T pro transformaci posunutí mužebýt obdélníková. Napríklad pri odvozování globální matice tuhosti taženého–tlacenéhoprvku mohl mít vektor posunutí pouze dva cleny: r =

u1 u2

. Potom pri globalizaci

do 2D ci 3D je transformacní matice rozmeru 2×2, ci 3×3. Transformace pro maticituhosti funguje správne, nebot’ lokální matice tuhosti má nulové cleny v odpovídají-cích smerech. Tato transformace by však nefungovala správne pro transformaci ma-tice hmotnosti. Záver tedy je, že lokální matice hmotnosti musí mít všechny translacnícleny. Pro tažený–tlacený prvek musíme tedy definovat vektor uzlových posunutí vcetneposunu kolmých ke strednici prutu (vi):

u =u1 v1 u2 v2

T (12.14)

Konzistentní a diagonální matice hmotnosti v lokálních souradnicích budou:

Mcons =1

6ρAl

2 0 1 00 2 0 11 0 2 00 1 0 2

Mdiag =1

2ρAl

1 0 0 00 1 0 00 0 1 00 0 0 1

(12.15)

Globalizaci proved’me pomocí 12.13. Nejprve pro diagonální matici hmotnosti:

Mdiag = T TMdiagT =1

2ρAlT T I4T =

1

2ρAlT TT =

1

2ρAlI4 (12.16)


Je zrejmé, že globální matice hmotnosti je shodná s lokální maticí hmotnosti. Pro kon-zistentní matici hmotnosti bude transformace:

Mcons = T TMconsT =

c −s 0 0s c 0 00 0 c −s0 0 s c

1

6ρAl

2 0 1 00 2 0 11 0 2 00 1 0 2

c s 0 0−s c 0 00 0 c s0 0 −s c

(12.17)

Kde c = cosα a s = sinα. Po provedení maticového násobení zjistíme, že dojdeme ke stej-nému záveru jako u diagonální matice hmotnosti (T TMconsT = M cons) a to díky clenum(c2 + s2

)a (cs− sc) a blokové forme matice hmotnosti. Toto lze dokázat následovne. Na

chvíli serad’me vektor koncových posunutí takto: u =u1 u2 v1 v2

T – tedy nejprveposunutí ve smeru osy x a pak posunutí ve smeru osy y. Dále zapišme matici hmot-nosti pomocí dvou submatic M . Matici transformace mužeme díky preskupení vektoruposunutí prepsat rovnež pomocí submatic s kosíny a síny:

M cons =

[M 0

0 M

]

T =

[cI2 sI2−sI2 cI2

]

(12.18)

kdeM =

1

6ρAl

[2 11 2

]

(12.19)

Proved’me nyní transformaci:

Mcons = T TMconsT =

[cI2 −sI2sI2 cI2

] [M 0

0 M

][cI2 sI2−sI2 cI2

]

=

=

[ (c2 + s2

)M (cs− sc) M

(sc− cs) M(c2 + s2

)M

]

=

[M 0

0 M

]

(12.20)

Z uvedeného plyne, že transformace tzv. diagonální matice s opakujícími se bloky jerovna té samé matici. Duležitý záver pro tranformaci matice hmotnosti je následující:Bloková matice hmotnosti (lokální) se globalizuje ve stejnou matici tehdy, když všechnystupne volnosti jsou translacní a když všechny patrí do stejného globálního souradnéhosystému. Této vlastnosti je samozrejme vhodné využít a pro matice hmotnosti preskocittransformaci z lokálního do globálního souradného systému.

Uvedeného postupu nelze použít v následujících prípadech:

• element má jiné než translacní stupne volnosti

• bloky v matici hmotnosti nemají stejnou velikost (napr. element s krivými hranami,skorepinový element)

V techto prípadech je nutno použít obvyklého transformacního vztahu 12.13.

12.3 Vlastní kmitání s více stupni volnosti

12.3.1 Metoda konstant tuhosti

Uvažujme vlastní kmitání pružného proste podepreného nosníku z obrázku 12.1. Kmi-tání budeme rešit jako kmitání soustavy se dvema stupni volnosti. Do tretin rozpetí

12.3. VLASTNÍ KMITÁNÍ S VÍCE STUPNI VOLNOSTI 115

k21v1(t)

v1(t) v2(t)

k11v1(t)

k12v2(t) k22v2(t)

m1 m2

Obrázek 12.1: kmitání nosníku se dvema stupni volnosti


nosníku umístíme hmotné body a do nich sostredíme odpovídající cást hmotnosti nos-níku.

Zapišme pohybové rovnice (podmínky rovnováhy sil pusobících na jednostlivé hmotnébody ve smeru kmitání):

m1v1 + k11v1 + k12v2 = 0 (12.21)m2v2 + k21v1 + k22v2 = 0

Rešení této soustavy homogenní diferenciálních rovnic druhého rádu hledejme ve tvaru:

v1(t) = v01sinωt

v2(t) = v02sinωt

Dvakrát derivujme dle casu a dosad’me zpet do 12.21:

−m1ω2v0

1sinωt+ k11v01sinωt+ k12v

02sinωt = 0

−m2ω2v0

2sinωt+ k21v01sinωt+ k22v

02sinωt = 0

Obe rovnice mužeme delit sinωt:

[k11 −m1ω

2]v01 + k12v

02 = 0 (12.22)

k21v02 +

[k22 −m2ω

2]v02 = 0

Podmínkou existence netriviálního rešení soustavy lineárních algebraických rovnic jerovnost determinantu matice levé strany nule:

(k11 −m1ω

2) (k22 −m2ω

2)− k12k21 = 0

Toto je polynom ctvrtého stupne pro ω. Snadno se lze presvedcit, že této rovnici odpoví-dají dve ruzná rešení vetší než nula.

0 < ω1 < ω2

Pro kmitání s více stupni volnosti platí, že lze stanovit maximálne tolik vlastních frek-vencí (a jim odpovídajících vlastních tvaru), kolik má daná soustava stupnu volnosti.Dosadíme-li vypoctené hodnoty vlastních frekvencí ωi zpet do 12.22, zjistíme, že tytorovnice jsou lineárne závislé, t.j není možné spocítat hodnoty výchylek odpovídající jed-notlivým vlastním frekvencím, pouze jejich vzájemný pomer. Toto si ukážeme názornejina príkladu v následující kapitole.

12.3. VLASTNÍ KMITÁNÍ S VÍCE STUPNI VOLNOSTI 117

%&%%&%'&''&'

(&((&()&))&)

k1

k2

y2(t)

m2

m1

y1(t) *+**+*,+,,+,

-+--+-.+..+.

k1(y1 − y2)

k2y2

k1(y2 − y1)

m2y2

m1y1

Obrázek 12.2: kmitání se dvema stupni volnosti

12.3.2 Príklad (metoda konstant tuhosti)

Soustava dvou hmotných bodu m1a m2 kmitá na pružinách o tuhostech k1resp. k2.Urcete vlastní kruhové frekvence volného kmitání a sestavte pohybové rovnice.Rešení:V case promenné vertikální výchylky obou hmotných bodu oznacme y1(t) a y2(t). Nahmotné body pusobí tíhové síly, vratné síly zpusobené stlacením, ci protažením pružina setrvacné síly. Posune-li se hmota m1o výchylku y1a zároven hmota m2 o y2, vzniknev pružine 1vratná síla o velikosti k1(y1 − y2). Protože se jedná o vlastní kmitání, tíhovésíly pusobící na hmotné body by se objevily na pravé strane pohybových rovnic a navlastní kruhové frekvence kmitání nemají vliv (v zápisu vynecháváme). Pohybové rovnice(pseudostatické podmínky rovnováhy ve smeru kmitání):

m1y1 + k1(y1 − y2) = 0

m2y2 + k1(y2 − y1) + k2y2 = 0

prepišme maticove:[m1 00 m2

]y1y2

+

[k1 −k1

−k1 k1 + k2

]y1y2

=

00

rešení této soustavy homogenních diferenciálních rovnic druhého rádu budeme hledatve tvaru:

y1(t) = y01sinωt

y2(t) = y02sinωt

pro dosazení budeme potrebovat druhé derivace dle casu:

y1(t) = −ω2y01sinωt

y2(t) = −ω2y02sinωt

po dosazení a jednoduché uprave:[k1 −m1ω

2 −k1

−k1 (k1 + k2)−m2ω2

]y1y2

=

00

(12.23)


podmínka rešitelnosti soustavy lineárních algebraických rovnic:

det

[k1 −m1ω

2 −k1

−k1 (k1 + k2)−m2ω2

]

= 0

(k11 −m1ω

2) (k1 + k2 −m2ω

2)− k2

1 = 0

po dosazení dostaneme:(10− 5ω2

) (10 + 20− 10ω2

)− 100 = 0

1⇒ ω1 =√

1 = 1s−1

ω2 =

4⇒ ω2 =√

4 = 2s−1

Hodnoty vlastních kruhových frekvencí soustavy ω1, ω2 dosad’me do 12.23. Obdržímesoustavu lineárních algebraických rovnic pro neznámé výchylky (resp. amplitudy vý-chylek) y0

1 a y02:

[10− 5 −10−10 30− 10

]y01

y02

=

00

Snadno zjistíme, že se jedná o soustavu závislých rovnic (druhá rovnice je lineárníkombinací prvé), ze které nelze spocítat hodnoty výchylek y0

1 a y02, pouze jejich pomer:

5y1 − 10y2 = 0

−10y1 + 20y2 = 0

y01

y02

=2

1

Výsledné pohybové rovnice vlastního kmitání sostavy hmotných bodu budou mít tedytvar:

y1(t) = 2sin1t

y2(t) = 1sin2t

12.4 Rešení kmitání pomocí MKP

Aplikací Hamiltonova principu snadno dojdeme k výchozím rovnicím pro rešení jednot-livých druhu kmitání. V základním kurzu dynamiky jsme kmitání rozdelovali na ctyriskupiny, na kmitání vlastní netlumené, vlastní tlumené, vynucené netlumené a vynucenétlumené. Jednoduchou analogií k pseudostatickým podmínkám rovnováhy dospejeme kzákladním rovnicím pro jednotlivé typy kmitání.

Kmitání vlastní netlumené:[M ] r+ [K] r = 0 (12.24)

vlastní tlumené:[M ] r+ [C] r+ [K] r = 0 (12.25)

vynucené netlumené:[M ] r+ [K] r = F1 sinω1t (12.26)

12.5. VLASTNOSTI VLASTNÍCH TVARU 119

vynucené tlumené:[M ] r+ [C] r+ [K] r = F1 sinω1t (12.27)

Význam jednotlivých clenu v techto rovnicích je zrejmý, pro úplnost strucne zopaku-jeme:

• [M ] matice hmotnosti

• r vektor uzlových zrychlení

• [C] matice tlumení

• r vektor uzlových rychlostí

• F1 vektor amplitud budících sil

• ω1 kruhová frekvence budící síly

12.5 Vlastnosti vlastních tvaru

Pro vlastní netlumené kmitání platí maticový zápis podmínek rovnováhy:

[M ] r+ [K] r = 0

kde M je matice hmotnosti celé soustavy a K je matice tuhosti soustavy. Rešení hle-dáme ve známém tvaru, který zapíšeme do sloupcového vektoru r(t):

r(t) =r0sinωt

r(t) = −ω2r0sinωt

Po dvojím derivování tohoto výrazu dle casu a dosazení do podmínek rovnováhy zís-káme:

− [M ]ω2r0sinωt+ [K]

r0sinωt = 0

z cehož plyne duležitá rovnice pro vlastní kmitání:

([K]− ω2 [M ]

) r0

= 0 (12.28)

12.5.1 Ortogonalita vlastních tvaru

Z rovnice 12.28 platí pro každou vlastní frekvenci ωi a ωj:

−ω2iMri +Kri = 0

−ω2jMrj +Krj = 0


Vynásobme první rovnici zleva rTj a druhou rovnici (rovnež zleva) rT

i :

−ω2i r

Tj Mri + rT

j Kri = 0 (12.29)−ω2

j rTi Mrj + rT

i Krj = 0

Transpozicí vztahu rTj Mri a obdobne vztahu rT

j Kri dostaneme:

(rTj Mri

)T= rT

i MT rj

(rTj Kri

)T= rT

i KT rj

Vzhledem k symetrii matic M a K tudíž platí:

rTj Mri = rT

i Mrj (12.30)rTj Kri = rT

i Krj

Po uvážení techto vlastností a sectení rovnic 12.29 dostaneme:(ω2

i − ω2j

)rTi Mrj = 0

Vlastní kruhová frekvence ωi je ruzná od ωj. Potom mužeme psát:

rTi Mrj = 0

Toto se dá slovy vyjádrit takto: Vlastní vektory jsou ortogonální vzhledem k matici hmot-nosti. Strucne se ríká, že vlastní tavry jsou M-ortogonální.

Nyní vynásobme první rovnici ve vztahu 12.29 1ω2

i

a obdobne druhou 1ω2

j

:

−rTj Mri +

1

ω2i

rTj Kri = 0

−rTi Mrj +

1

ω2j

rTi Krj = 0

rovnice secteme a uvažme 12.30:(ω2

i − ω2j

)rTi Krj = 0

Toto se dá slovy vyjádrit takto: Vlastní vektory jsou ortogonální i vzhledem k matici tu-hosti. Strucne: vlastní vektory jsou K-ortogonální.

Jak už vyplynulo z predchozí kapitole, velikosti vlastních tvaru nemají žádný fyzikálnívýznam, proto je zvykem vlastní tvary normovat vhledem k matici hmotnosti. Toto seprojeví jako velmi výhodné pri rešení kmitání rozvojem do vlastních tvaru. Chceme-livlatní tvary normovat vzhledem k matici hmotnosti, budeme chtít, aby platilo:

rTi Mrj = I (12.31)

kde I je jednotková matice. Pro normovaný vlastní tvar bude tedy platit:

ri,norm =ri

√

rTi Mri

(12.32)

12.6. TLUMENÉ KMITÁNÍ 121

dosad’me 12.31 do první rovnice v 12.29:

−ω2i I + rT

j Kri = 0

ihned je videt, že:rTi Krj = ω2

i

Nejcasteji se s vlastními tvary pracuje v tomto tvaru, t.j. vlastní tvary jsou ortonormálnívzhledem k matici hmotnosti, strucne M-ortonormální.

Zapamatujme si duležitý poznatek o volném kmitání soustav s více stupni volnosti:

• Velikosti posunutí a z nich odvozené veliciny (napr. napetí) nemají u vlastníhokmitání smysl, ci fyzikální význam.

12.6 Tlumené kmitání

Uvažujme vynucené tlumené kmitání. Toto kmitání je urcené pohybovou rovnicí:

[M ]r+ [C]r+ [K]r = F (t) (12.33)

Kde [C] je matice viskózního tlumení. Obecne by se konzistentní matice tlumení elementudala získat obdobne jako konzistentní matice hmotnosti pomocí vztahu:

C =

∫

V

NTµNdV

V praxi se však tento vztah nevyužívá, nebot’ koeficienty lineárního tlumení µ je velmisložité získat, ale protože vliv tlumení je obecne menší než vliv setrvacnosti a tuhosti, jemožné matici tlumení urcit ve tvaru lineární kombinace matice hmotnosti a tuhosti:

[C] = α[M ] + β[K] (12.34)

kde α a β jsou proporciální koeficienty zohlednující vliv rychlosti uzlových bodu a rych-losti deformace. Tyto proporcionální koeficienty získáme pomocí experimentálního me-rení. Tomuto typu tlumení se ríká Rayleighovo tlumení.

12.7 Rešení vlastního kmitání

Úlohou vlastního kmitání je nalézt vlastní kruhové frekvence a jim príslušné vlastnítvary kmitání. Pro prípad vlastního kmitání nabyde pohybová rovnice tvaru:

[M ]r+ [K]r = 0 (12.35)

Rešení této diferenciální rovnice hledejme ve tvaru:

r(t) = r0 . sinωt (12.36)

dvojím derivováním dle casu a zpetným dosazením:

− [M ] rω2sinωt+ [K] r sinωt = 0(−ω2 [M ] + [K]

)r = 0 (12.37)


prevedeme soustavu diferenciálních rovnic druhého rádu na problém vlastních císel promatice K a M . Netriviální rešení této rovnice nalezneme splnením podmínky rešitelnostisoustavy rovnic ve tvaru:

det(−ω2 [M ] + [K]) = 0 (12.38)

Nejcasteji je tato úloha rešena nekterou z iteracních metod, pro ilustraci uvedeme Sto-dolovu iteracní metodu, metodu inverzní iterace a metodu iterace podprostoru.

12.7.1 Stodolova iteracní metoda

vynásobme rovnici 12.37 zleva maticí poddajnosti [δ]:(

[δ] [K]− [δ]ω2(j) [m]

)r(j)

= 0(ω2 [δ] [m]− [E]

) r(j)

= 0 (12.39)

nebot’ platí: [δ] [K] = [E], tedy soucin matice poddajnosti a matice tuhosti je jednotkovámatice. Tímto jednoduchým postupem jsme získali základní vztah metody konstantpoddajnosti, která je analogickou úlohou k metode konstant tuhosti. Metodou konstantpoddajnosti se dále nebudeme zabývat, vynásobme pouze tuto rovnici nyní clenem 1/ω2:

(

[δ] [m]− 1

ω2(j)

[E]

)

r(j)

= 0

tuto rovnici je nyní již snadné upravit na tvar vhodný k iterování:

1

ω2(j)

r(j)

= [δ] [m]r(j)

aproximaci vlastního tvaru pomocí tvaru predchozího mužeme tedy zapsat pomocí:r(j)i+1

= ω2(j) [δ] [m]

r(j)i

Toto je základní výraz pro Stodolovu iteracní metodu, jejíž postup si strucne uvedeme vnásledujících krocích:

1. volba pocátecního tvaru kmitánír(j)0- tento tvar lze zvolit libovolný (pozor však

na to, že jednotlivé cleny nesmí být navzájem lineárne závislé, napr. tvar 1, 0, 0, ...., 0Tje vhodný. Casteji se však používá tvar vypocítaný jako deformace od vlastní tíhykonstrukce.

2. výpocet normovaného tvaru:

r(j)0

=

r(j)0

(r(j)0)T

[m]r(j)0

3. Gramm-Schmidtova ortogonalizace: všechny vlastní tvary kmitání musí být vzá-jemne ortogonální, proto je potreba zajistit, aby byl napr. tretí vlastní tvar kolmý k

12.7. REŠENÍ VLASTNÍHO KMITÁNÍ 123

prvnímu a druhému vlastnímu tvaru. Tento krok je samozrejme vynechán pri vý-poctu prvního vlastního tvaru (není zatím k cemu ortogonalizovat). Tak tedy napr.pri výpoctu tretího vlastního tvaru musí být splneny tyto podmínky ortogonality:

(r(1)0)T

[m]r(3)0

= 0(r(2)0)T

[m]r(3)0

= 0

pri výpoctu tretího vlastního tvaru je tedy potreba vyloucit první a druhý tvarpomocí podmínek ortogonality. Toto lze zapsat pomocí:

[m]r(j)0

= [m]r(j)0 −

j−1∑

k=1

C(k) [m]r(k)

0

kdeC(k) =

(r(j)0)T

[m]r(k)

0

4. výpocet vlastní kruhové frekvence príslušné vlastnímu tvaru:

ω2i+1 =

(r(j)i+1

)T

[m]r(j)i

(r(j)i+1

)T

[m]r(j)i+1

5. zastavovací podmínka: ∣∣∣∣

(

ωi+1(j)

)2

−(

ωi(j)

)2∣∣∣∣

ω2i+1(j)

5 10−2s

kde s je pocet desetinných míst potrebných pro danou presnost.

12.7.2 Inverzní iterace

Metoda inverzní iterace je založená na Stodolove iteracní metode. Opet volíme pocátecníaproximaci vlastního tvaru, kterou oznacíme r0. Vyjdeme ze vztahu:

[K] r = ω2 [M ] r (12.40)

Tento vztah platí pro libovolný vlastní tvar ri a jemu príslušnou vlastní frekvenci ωi.Nyní vypocteme amplitudu setrvacných sil Fsetrv

12.7.3 Metoda iterace podprostoru


Kapitola 13

Rešení úloh nelineárnímechaniky

Dosud jsme predpokládali lineární elastický izotropní materiál v nemž pri zatížení vzni-kají pouze malé deformace. Tato predstava o materiálu je však znacne vzdálena tomu,jak se reálný materiál skutecne pri zatežování chová. Platnost Hookova zákona je uvetšiny reálných materiálu omezena pouze na malé deformace.

...

Nelineární problémy je možno rozdelit do trí skupin a to:

• materiálová (fyzikální) nelinearita

• geometrická nelinearita

• nelinearita zpusobená zmenou stavu (kontaktní úloha)

Materiálová nelinearita spocívá v nelineárních konstitutivních vztazích. Napetí neníprímo úmerné deformaci. Pri geometrické nelinearite nemužeme deformace považovatza nekonecne malé. Kategorie geometrických nelinearit zahrnuje velká posunutí a velkédeformace.

13.1 Prístupy k rešení nelineárních problému

Budeme se pouze zabývat temi nejbežnejšími zpusoby rešení nelineárních problému,tudíž si budeme moci dovolit zvolit jeden ze trí výše uvedených typu nelinearit a vyložitproblematiku pouze pomocí tohoto jednoho prípadu, napr. pomocí materiálové neli-nearity. Materiálovou nelinearitu je nejspíš nejjednodušší si predstavit graficky. Pouzenaznacíme jakým zpusobem se popsané metody dají rozšírit na další dve kategorie,zejména pak na geometrickou nelinearitu.

Rešení nelineárních úloh pomocí metody konecných prvku je vetšinou provádeno po-mocí trí základních technik:

• inkrementální (postupné) - INCREMENTAL OR STEPWISE METHODS

• iteracní (Newtonovy) - ITERATIVE OR NEWTON METHODS

125

126 KAPITOLA 13. REŠENÍ ÚLOH NELINEÁRNÍ MECHANIKY

• smíšené - STEP-ITERATIVE OR MIXED PROCEDURES

V prípade materiálové nelinearity budou mít podmínky rovnováhy pro element tvar:

[k][r] = f (13.1)

V teto rovnici však matice tuhosti elementu [k] je však funkcí nelineárních materiálo-vých vlastností [D(σ)]. Toto mužeme zapsat symbolicky ve tvaru:

[k] = [k(r, f)] (13.2)

Toto je schematicky znázorneno na obr. 13.1.

kr = f

r

f

ε

σ = D(σ)ε

σ

Obrázek 13.1: Materiálová nelinearita (a) zatížení vs. posunutí (b) napetí vs. deformace

13.2 Inkrementání technika

Základem inkrementálního postupu je rozdelení zatížení na nekolik malých kroku.Casto jsou tyto kroky voleny o stejné velikosti, ale obecne nemusí být. Prírustek zatíženíje aplikován v každém kroku a výpocet v každém kroku spocívá v rešení soustavy line-árních rovnic. Jinými slovy, matice tuhosti [K] je stejná v každém zatežovacím kroku,ale liší se mezi jednotlivými kroky. Výpocet každého kroku je založen na prírustu po-sunu r, tento inkrementální postup je opakován až do okamžiku dosažení celkovéhozatížení. Dá se také ríci, že inkrementální postup prevádí nelineární úlohu na úlohu pocástech lineární.

Chceme-li popsat inkrementální techniku, uvažujme celkové zatížení dané vektoremF. Pocátecní (referencní) stav je dán vektorem pocátecních sil F0 a vektorem pocá-tecních posunutí r0. Nejcasteji tyto vektory budou oba nulové, v prípade, že zacínámez nezatíženého, nedeformovaného stavu. Rozdelíme celkové zatížení na N prírustu, t.j.celkové zatížení aplikované na konstrukci bude:

F = F0+

N∑

j=1

∆Fj (13.3)

13.2. INKREMENTÁNÍ TECHNIKA 127

Po i-tém prírustku je zatížení dáno vektorem:

Fi = F0+

i∑

j=1

∆Fj (13.4)

a vektor posunutí:

ri = r0+

i∑

j=1

∆rj (13.5)

Abychom mohli spocítat prírustek posunutí (posunutím je rízen celý výpocet), použijemevýpoctu matice tuhosti na základe predešlého prírustku:

[ki−1]∆ri = ∆Fi pro i = 1, 2, 3...,M (13.6)

kde[ki−1] = [ki−1(ri−1, Fi−1)] (13.7)

Pocátecní hodnota matice tuhosti [k0] se vypocte z materiálových konstant na základepracovního diagramu na zacátku zatížení. Schema inkrementálního postupu je znázor-neno na obr. 13.2.

r

∆ri = k−1

i−1∆fi

1

k0

∆fi

/1032547698:6;0=<?>A@10CBED769F7690CB

G D76;H703IJD76;F7690CBfi

fi−1

ki−1

f

ri−1 ri

Obrázek 13.2: Inkrementální prístup

Obycejne se v inkrementálním prístupu využívá tecného modulu pružnosti k výpoctu[D(σ)] a k výpoctu matice tuhosti [K]. Tato matice tuhosti se pak nazývá tecná maticetuhosti.

Inkrementální metoda je analogická s numerickými metodami používanými k integracisoustav lineárních i nelineárních diferenciálních rovnic, jako je Eulerova metoda, nebometoda Runge-Kutta. Ke zpresnení výše uvedené metody se nejcasteji používá zjemnenídelení (rozdelení zatížení na více prírustku). Musíme si ovšem uvedomit, že tecná maticetuhosti musí být prepocítána v každém kroku, to znamená velké nároky na výpocetnícas.


13.3 Iteracní (Newtonovy) techniky

Pri iteracních metodách, na rozdíl od inkrementální techniky, je v každém kroku apli-kováno celkové zatížení. Protože v každém kroku používáme približnou hodnotu maticetuhosti v každém kroku, není dosaženo rovnováhy a po každém kroku musí dojít kvyrovnání. Toho je dosaženo tím, že po každé iteraci je vypocítán tzv. vektor nevyrovna-ného zatížení a v dalším kroku je využit k výpoctu prírustku posunutí. Tento postup jeopakován do té doby, než je dosaženo rovnováhy s požadovanou presností.

Oznacme opet F0 vektor pocátecního zatížení a r0 vektor pocátecních posunutí. Vobecném prípade tyto nebudou vždy nulové. Po i-tém cyklu iteracního procesu budezatížení:

Fi = F − Fe,i−1 (13.8)

kde F je celkové zatížení a Fe,i−1 je zatížení vyrovnané po predchozím kroku. Prí-rustek posunutí je behem i-tého kroku spocítán pomocí:

[k(i)]∆ri = Fi (13.9)

Celkové posunutí po i-tém kroku je vypocítano z:

ri = r0+i∑

j=1

∆rj (13.10)

Na konci i-tého iteracního kroku je vypocítán vektor nevyrovnaných sil Fe,i jako za-tížení nezbytné k zachování posunutí ri. Tato procedura je opakována do té doby,než prírustky posunutí nebo nevyrovnané síly jsou nulové, respektive menší než danálimitní hodnota.

V každém iteracním kroku je treba vypocítat matici tuhosti. Nejbežnejším zpusobemje pocítat tecnou matici tuhosti na koci predchozího iteracního kroku, t.j. jako tecnu kzatežovací krivce F − r v bode ri−1, Fi−1.Tento postup znamená vycíslení tecné matice tuhosti v každém iteracním kroku, cožmuže být znacne casove nárocné a proto se casto používá modifikovaná iteracní tech-nika, která neprepocítává v každém kroku matici tuhosti, nýbrž po celou dobu iteracepocítá s pocátecní tecnou maticí tuhosti. Rozdíl mezi klasickou a modifikovanou iteracnítechnikou je patrný z obr. 13.3.

13.4 Smíšené techniky (STEP-ITERATIVE OR MIXED)

Tyto techniky využívají kombinace inkrementálního a iteracního postutpu. Jeden z mož-ných prístupu je patrný z obr. 13.5 kde zpresnení výpoctu je dosaženo v každém prí-rustku provedením nekolika modifikovaných iteracních kroku, t.j. bez prepocítávánítecné matice tuhosti.

13.4. SMÍŠENÉ TECHNIKY (STEP-ITERATIVE OR MIXED) 129

F3

rr4r3r2r1

∆r1 ∆r2

Fe2

F4

k0

1

Fe1

F

F2

ki−1

Obrázek 13.3: Iteracní metoda

r3r2r1r

F

k0

F

Obrázek 13.4: Modifikovaná iteracní metoda

r

∆Fi+1

Fi

∆Fi

Fi−1

ri−1 ri

F

Obrázek 13.5: Smíšená (mixed) technika


Kapitola 14

Numerická integrace

Pro výpocet matice tuhosti jakéhokoliv prvku je treba provést numerickou integracisoucinu matic. Integrál matice je roven matici integrálu jednotlivých prvku dané ma-tice. Numerickou integraci budeme cílit na interval < −1; +1 >, pouze zrídka výsledekzobecníme pro interval < a; b >. Pripomenme, že matice tuhosti prvku se dá zapsat po-mocí

∫

ΩBTDB dΩ, nám však vzhledem k výše uvedenému postací se omezit na integrál

obecné funkce f(ξ), f(ξ, η), popr. f(ξ, η, ζ):∫ 1

−1

f(ξ) dξ

∫ 1

−1

∫ 1

−1

f(ξ, η) dξ dη

∫ 1

−1

∫ 1

−1

∫ 1

−1

f(ξ, η, ζ) dξ dη dζ (14.1)

Funkci f(ξ) nahradíme na intervalu < −1,+1 > polynomem, který oznacíme F (ξ). Apro-ximaci funkce f(ξ) stanovíme tak, že ve vybraných bodech ξi z intervalu < −1; +1 >bude platit f(ξi) = F (ξi), t.j. polynomická funkce bude temito body procházet a hledanýintegrál

∫ 1

−1 f(ξ) dξ urcíme jako∫ 1

−1 F (ξ) dξ. Presnost takto vypocteného integrálu se budezvyšovat s poctem zvolených bodu na daném intervalu.K aproximaci funkce f(ξ) se casto používají Lagrangeovy polynomy. Zvolme na inter-valu 〈−1, 1〉 (n+ 1)bodu, kterými mužeme proložit polynom n-tého stupne. Pro libovolnýstupen bude mít Lagrangeuv polynom následující tvar:

hi(ξ) =(ξ − ξ0)(ξ − ξ1)...(ξ − ξi−1)(ξ − ξi+1)...(ξ − ξn)

(ξi − ξ0)(ξi − ξ1)...(ξi − ξi−1)(ξi − ξi+1)...(ξi − ξn)(14.2)

Aproximacní funkci F (ξ) lze vyjádrit pomocí:

F (ξ) = f0h0(ξ) + f1h1(ξ) + ...+ fnhn(ξ) (14.3)

Dle zpusobu delení intervalu si rozdelíme a detailneji popíšeme dva možné zpusobynumerické integrace. První prístup rozdeluje daný interval na n stejne velkých podin-tervalu, zatímco druhý zpusob se snaží o racionálnejší umístení integracních bodu, tedyjedná se o nerovnomerné delení intervalu.

14.1 Rovnomerné delení (Newton-Cotesovo schema)

Rozdelme tedy nejprve daný interval rovnomerne, t.j. pro jednotlivé body ξi na intervalu〈−1, 1〉 bude platit:

ξ0 = −1, ξ1 = −1 +2

n, ξ2 = −1 +

4

n, ..., ξn−1 = −1 +

2(n− 1)

n, ξn = 1

131

132 KAPITOLA 14. NUMERICKÁ INTEGRACE

Presná hodnota integrálu∫ 1

−1f(ξ)dξ potom bude:

∫ 1

−1

f(ξ) dξ =

n∑

i=0

fi

∫ 1

−1

hi(ξ) dξ +Rn (14.4)

Výraz∫ 1

−1 hi(ξ) dξ lze vycíslit pro jednotlivé stupne Lagrangeových mnohoclenu a tytovýrazy se potom nazývají Newton-Cotesovými konstantami Cn

i . Vzorec 14.4 lze potomprepsat ve tvaru:

∫ 1

−1

f(ξ) dξ =

n∑

i=0

fiCni +Rn (14.5)

Výhodou Newton-Cotesova schematu je fakt, že jej lze snadno zobecnit pro libovolnýinterval 〈a; b〉. Výraz pro približnou (tedy beze zbytku) hodnotu integrálu

∫ b

a f(x) dx budemít tvar:

∫ b

a

f(x) dx = (b− a)n∑

i=0

fiCni (14.6)

14.2 Nerovnomerné delení (Gaussova integracní formule)

Nyní se pokusme o následující úvahu. Daný interval 〈−1, 1〉 nebudeme delit rovnomerne,ale integracní body (body aproximace dané funkce) se budeme snažit volit tak, abyvýsledná približná hodnota integrálu byla co nejbližší jeho presné hodnote. Nyní tedynebudou neznámými pouze váhové koeficienty v aproximacním predpisu, ale i body nadaném intervalu. V tomto prípade však pro vyjádrení daného integrálu potrebujemedvojnásobný pocet parametru. Daný integrál zapíšeme tedy jako následující soucet:

∫ 1

−1

f(ξ) dξ ≈ α1f(ξ1) + α2f(ξ2) + ...+ αnf(ξn) (14.7)

kde αi jsou tedy váhové koeficienty a ξi jsou body na intervalu 〈−1, 1〉. Duležité je siuvedomit, že obojí jsou nyní neznámé. V predchozím odstavci jsme aproximovali hleda-nou funkci polynomem (n − 1) stupne. Nyní máme však k dispozici dvojnásobný pocetparametru, aproximacní funkce F (ξ) muže tedy být až polynom (2n − 1) stupne. Hle-dáme tedy další aproximacní funkci, která zvýší stupen aproximacního polynomu, alezároven zustane nulovou ve všech n bodech, kterými procházejí naše Lagrangeovy poly-nomy. Nejsnadnejší je rozšírit funkci (ξ−ξ1)(ξ−ξ2)...(ξ−ξn) obecným polynomem n-téhostupne:

F (ξ) = (ξ − ξ1)(ξ − ξ2)...(ξ − ξn)(A0 +A1ξ +A2ξ2 +A3ξ

3 + ...+Anξn) (14.8)

Hledáme tedy aproximaci funkce f(ξ) ve tvaru:

f(ξ) ≈ F (ξ) = (ξ − ξ1)(ξ − ξ2)...(ξ − ξn)︸︷︷︸

p(ξ)

n∑

i=0

Aiξi (14.9)

Integrujme:∫ 1

−1

f(ξ) dξ ≈∫ 1

−1

F (ξ) dξ =

n∑

i=0

Ai

∫ 1

−1

p(ξ)ξi−1 dξ (14.10)

14.2. NEROVNOMERNÉ DELENÍ (GAUSSOVA INTEGRACNÍ FORMULE) 133

Budeme se nyní snažit vybrat souradnice ξi tak, aby se všech n integrálu∫ 1

−1p(ξ)ξi−1 dξ

rovnalo 0. Tyto podmínky vedou na soustavu nelineárních rovnic:∫ 1

−1

p(ξ)ξi−1 dξ = 0 (i = 1, 2, ..., n) (14.11)

Tuto vlastnost mají tzv. Legendrovy polynomy a ξi jsou tedy koreny Legendrových po-lynomu. Budeme-li tedy znát tyto koreny ξi, lze váhové koeficienty αi spocítat snadnopomocí:

αi =

∫ 1

−1

hi(ξ) dξ (14.12)

Ukažme si výpocet korenu i váhových koeficientu pro pár prvních delení intervalu〈−1, 1〉. Pro prípad jediného integracního bodu uprostred intervalu (pro sudý pocet dílkubude vždy jeden integracní bod uprostred intervalu, tedy v nule). Spocítejme tedy h1(ξ)a α1:

h1(ξ) = 1

α1 =

∫ 1

−1

1 dξ = 2

Nyní hledejme Gaussovy integracní doby pro prípad delení na dva intervaly. V tomtoprípade využijme symetrie ξ1 = −ξ2 a rešme soustavu nelineárních rovnic ve tvaru:

∫ 1

−1

(ξ − ξ1)(ξ − ξ2) dξ = 0

∫ 1

−1

(ξ − ξ1)(ξ − ξ2)ξ dξ = 0

Druhá podmínka je díky symetrii intervalu identicky splnena. Rešme tedy první pod-mínku a dosad’me podmínku symetrie:

∫ 1

−1

(ξ − ξ1)(ξ − ξ2) dξ =

∫ 1

−1

(ξ + ξ2)(ξ − ξ2) dξ =

∫ 1

−1

(ξ2 − ξ22) dξ =

=

[ξ3

3

]1

−1

− ξ22 [ξ]1−1 =

2

3− ξ222 = 0

Z této podmínky plyne, že:

ξ2 =1√3

ξ1 = − 1√3

Váhové koeficienty α1 = α2 získáme již snadno integrací:

α1 = α2 =

∫ 1

−1

ξ − ξ2ξ1 − ξ2

dξ = −√

3

2

∫ 1

−1

ξ − 1√3dξ = 1

134 KAPITOLA 14. NUMERICKÁ INTEGRACE

Pro prípad delení intervalu na tri díly si uved’me již jen podmínky, ze kterých urcímepotrebné váhové koeficienty i body na intervalu:

∫ 1

−1

(ξ − ξ1)(ξ − ξ2) dξ = 0

∫ 1

−1

(ξ − ξ1)(ξ − ξ2)ξ dξ = 0

∫ 1

−1

(ξ − ξ1)(ξ − ξ2)ξ2 dξ = 0

První i tretí podmínka jsou opet splneny identicky a integrací druhé získáme následujícípodmínku pro ξ3:

∫ 1

−1

(ξ2 − ξ23

)ξ2 dξ =

[ξ5

5− ξ23

ξ3

3

]1

−1

=2

5− 2

3ξ23 = 0

Z tohoto plyne, že ξ3 =√

35 .

Váhové koeficientyopet získámé integrací následujícího:

α1 =

∫ 1

−1

(ξ − ξ2)(ξ − ξ3)(ξ1 − ξ2)(ξ1 − ξ3)

dξ =5

9

α2 =

∫ 1

−1

(ξ − ξ1)(ξ − ξ3)(ξ2 − ξ1)(ξ2 − ξ3)

dξ =8

9

Date post:	26-Nov-2020
Category:	Documents
Upload:	others
View:	1 times
Download:	0 times

Ond rej Jirouıek 13. œnora 2006mech.fd.cvut.cz/members/jirousek/download/k618y2m1/ymkp_fem.pdf ·...

Documents