Numerické metody mechaniky II

Numerické metody mechaniky II Martin Mánek Ostrava 2024 2 Obsah 1. Metoda konečných prvků v mechanice kontinua ........................................................................... 3 2. Nelineární úlohy mechaniky ............................................................................................................ 5 2.1. Lineární mechanika ................................................................................................................ 5 2.2. Nelineární mechanika ............................................................................................................ 5 2.3. Geometrické nelinearity ........................................................................................................ 5 2.3.1. Velké posuvy ................................................................................................................. 5 2.3.2. Velké deformace ........................................................................................................... 5 2.3.3. Tenzory deformace ............................................................................................................ 6 2.3.4. Tenzory napětí ................................................................................................................... 6 2.4. Materiálové nelinearity .......................................................................................................... 6 2.4.1. Lineárně elastické modely ............................................................................................. 6 2.4.2. Hyperelastické materiálové modely .............................................................................. 7 2.4.3. Viskoelastické materiály ................................................................................................ 8 2.4.4. Plastické chování materiálů ........................................................................................... 9 2.5. Kontakt ................................................................................................................................. 10 2.5.1. Metoda Lagrangeových multiplikátorů ....................................................................... 11 2.5.2. Penaltový přístup ........................................................................................................ 11 2.6. Metody řešení základní rovnice ........................................................................................... 11 3. Modální analýza ............................................................................................................................ 12 4. Tlumení .......................................................................................................................................... 13 4.1. Materiálové tlumení ............................................................................................................ 13 4.2. Viskózní tlumení ................................................................................................................... 14 5. Nestacionární úlohy mechaniky .................................................................................................... 15 5.1. Implicitní algoritmus ............................................................................................................ 15 5.2. Explicitní algoritmus ............................................................................................................. 16 5.3. Porovnání algoritmů ............................................................................................................ 17 6. Ztráta stability ............................................................................................................................... 18 7. Teplotní úlohy ................................................................................................................................ 19 Použitá literatura ................................................................................................................................... 20 3 1. Metoda konečných prvků v mechanice kontinua Základní myšlenkou metody konečných prvků je rozložení těles na konečný počet menších častí – elementů, prvků. Oproti analýze chování celé součásti nebo soustavy součástí je popis chování těchto prvků poměrně jednoduchý. V deformační variantě dochází k náhradě posuvů náhradními funkcemi ve tvaru polynomů. Postup metody lze popsat v několika krocích: • Diskretizace neboli rozdělení geometrie na podoblasti – konečné prvky, elementy. • Formulace chování jednotlivých elementů. • Opětovné složení do výsledné soustavy rovnic popisující celý systém pomocí rovnic pro jednotlivé elementy. • Zadání okrajových a počátečních podmínek. • Řešení soustavy získaných rovnic a získání primárních neznámých – pro deformační variantu se jedná o posuvy. • Řešení dalších veličin odvozených s primárních neznámých, např. napětí a přetvoření atd. Formulace chování elementu Pro matematický popis chování bývá použit maticový zápis, proto je ho užito i zde. Spojitě rozložené posuvy prvku {𝑢} jsou vyjádřeny pomocí matice tvarových funkcí [𝐴] a zobecnělých posuvů {𝑞}. {𝑢}=[𝐴]{𝑞}. (1.1) Diskretizace přetvoření pro spojitě rozložená přetvoření {𝜀} jsou vyjádřena pomocí matice [𝐵] a zobecnělých posuvů {𝑞} v jednotlivých uzlech. Matice [𝐵] vychází ze vztahů mezi posuvy a přetvořením daných teorií pružnosti. {𝜀}=[𝐵]{𝑞}. (1.2) Pro výpočet napětí lze pro lineární chování izotropního homogenního materiálu uvažovat platnost Hookeova zákona. S využitím matice [𝐶] lze napětí nad elementem vyjádřit jako: {𝜎}=[𝐶]{𝜀}. (1.3) S využitím rovnice 1.2 dostaneme vztah: {𝜎}=[𝐶][𝐵]{𝑞}. (1.4) Pro odvození základní rovnice MKP využijeme princip virtuálních prací. Virtuální posuvy a přetvoření lze pro danou diskretizaci vyjádřit ve tvarech? {𝛿𝑢}=[𝐴]{𝛿𝑞}. (1.5) {𝛿𝜀}=[𝐵]{𝛿𝑞}. (1.6) V dalším postupu jsou zanedbány objemové síly, ale jsou uvažovány setrvačné síly na základě d´Alambertova principu s využitím hustoty 𝜌. {𝑋}=−𝜌{𝑢󰇘}. (1.7) Matice tvarových funkcí není funkcí času, proto lze vyjádřit zrychlení na základě zobecnělých zrychlení v uzlech. {𝑢󰇘}=[𝐴]{𝑞󰇘}. (1.8) Aplikací předchozích vztahů dostaneme rovnici platnou pro každý prvek v následujícím tvaru: {𝛿𝑞}𝑇∫[𝐵]𝑇{𝜎}𝑑𝑉 𝑉={𝛿𝑞}𝑇∫𝜌[𝐴]𝑇[𝐴]𝑑𝑉{𝑞󰇘} 𝑉+{𝛿𝑞}𝑇∫[𝐴]𝑇{𝑝}𝑑𝑆 𝑉+{𝛿𝑞}𝑇{𝐹𝑂𝑆}. (1.9) 4 Rovnice může být také napsána jako rovnost sil v uzlech, které odpovídají deformační energii prvku tedy vnitřním silám {𝑓𝑖} a vnějších sil {𝑓𝑒}, které jsou tvořeny silami objemovými {𝑓𝑉}, plošnými {𝑓𝑆} a osamělými {𝑓𝑂𝑆}. Platí pak: {𝑓𝑖}={𝑓𝑒}, (1.10) kde: {𝑓𝑖}=∫[𝐵]𝑇{𝜎}𝑑𝑉, 𝑉 (1.11) {𝑓𝑒}={𝑓𝑉}+{𝑓𝑆}+{𝑓𝑂𝑆}, (1.12) {𝑓𝑉}=∫𝜌[𝐴]𝑇[𝐴]𝑑𝑉 {𝑞󰇘}, 𝑉 (1.13) {𝑓𝑆}=∫[𝐴]𝑇{𝑝}𝑑𝑆. 𝑆 (1.14) Rovnováha vnitřních a vnějších sil musí platit pro lineární i nelineární příklady pružnosti. Při platnosti Hookeova zákona je možné napsat rovnici vycházející z této rovnováhy: [𝑘]{𝑞}={𝑓𝑒}, (1.15) kde [𝑘] představuje matici tuhosti daného prvku. [𝑘]=∫[𝐵]𝑇[𝐶][𝐵]𝑑𝑉. 𝑉 (1.16) V případě nestacionárních úloh je potřeba vyjádřit rovnici pro objemové síly v závislosti na matici hmotnosti prvku [𝑚]. {𝑓𝑉}=[𝑚]{𝑞󰇘}, (1.17) [𝑚]=∫𝜌[𝐴]𝑇[𝐴]𝑑𝑉. 𝑉 (1.18) Pohybová rovnice prvku má pak tvar: [𝑚]{𝑞󰇘}+[𝑘]{𝑞}={𝑓}, (1.19) kde: {𝑓}={𝑓𝑆}+{𝑓𝑂𝑆}. (1.20) Dalším krokem je opětovné složení jednotlivých elementů do celkové geometrie. S využitím znalosti polohy jednotlivých elementů lze z matic tuhosti a hmotnosti pro jednotlivé prvky sestavit globální matice tuhosti a hmotnosti celé konstrukce. V oblasti lineární pružnosti se předpokládají malé posuvy a přetvoření, lineární chování materiálu a neměnnost okrajových podmínek. Výsledný vztah pro metodu konečných prvků má pro nestacionární úlohy tvar: [𝑀]{𝑄󰇘}+[𝐾]{𝑄}={𝐹}. (1.21) Pro stacionární úlohy zůstane pouze tvar: [𝐾]{𝑄}={𝐹}. (1.22) 5 2. Nelineární úlohy mechaniky V mechanice je možné se často setkat s úlohami, kde již není možné dodržet veškeré předpoklady lineární pružnosti. V nelineární mechanice pružného tělesa existují příklady, které obsahují nelineární deformační charakteristiku, velké posuvy nebo deformace, které významně mění tvar původní geometrie. Rovněž spoje mezi tělesy často nelze popsat prostou lineární charakteristikou. 2.1. Lineární mechanika Pro lineární mechaniku je typická vysoká míra idealizace řešeného problému. Je však zaručena řešitelnost těchto lineárních modelů. Konečný stav modelu je závislý pouze na zadaných hodnotách zatížení, nikoliv na způsobu jejich dosažení a pořadí. Lze využít princip superpozice, kdy je výsledný stav složen z výsledků dílčích stavů. Základní předpoklady jsou: • Rovnováha je vyjádřena v nezdeformovaném stavu tělesa. • Mezi posuvy a deformacemi platí lineární vztahy. • Materiál je lineárně elastický. • Všechny vazby jsou lineární. 2.2. Nelineární mechanika V úlohách nelineárních na rozdíl od lineárních neexistuje věta o existenci řešení, konvergenci a stabilitě. Můžeme se tak setkat s úlohami, které v zadaných podmínkách nelze vyřešit. Rovněž u nich záleží na posloupnosti zátěžných stavů, neplatí tak princip superpozice. Prakticky můžeme rozlišit několik druhů nelinearit. • Geometrické nelinearity. • Materiálové nelinearity. • Strukturální nelinearity (kontakty). 2.3. Geometrické nelinearity Geometrické nelinearity jsou způsobeny velkými posuvy a natočeními, ty mohou způsobit i velké přetvoření. Z tohoto pohledu můžeme rozlišovat dva případy geometrických nelinearit, velké posuvy (a natočení) a velká přetvoření. Chování materiálu může být v obou případech lineární i nelineární. 2.3.1. Velké posuvy Konstrukce může vlivem zatížení dosáhnout velkého posunutí nebo natočení např. rybářský prut, avšak deformace zůstávají malé. Pro popis lze použít Cauchyho geometrické rovnice jako v lineární mechanice, ale rovnováha konstrukce musí být vyjádřena pro zdeformovaný tvar. 2.3.2. Velké deformace V tomto případě se k velkým posuvům a deformacím přidávají navíc ještě velké deformace (přetvoření). V lineární pružnosti při předpokladu malých deformací jsou zanedbány nelineární členy geometrických vztahů a používá se inženýrský (Cauchyho) tenzor přetvoření. Hranice použitelnosti se uvádí do hodnoty deformace 1 %. To je dostatečné pro většinu konstrukcí, nicméně v simulacích tváření nebo při použití poddajných materiálů jako pryže, plasty, pěny mohou deformace dosahovat desítek i stovek procent. Je pak nutno uvažovat i nelineární členy geometrických vztahů. 6 2.3.3. Tenzory deformace Tenzory deformace používané v mechanice kontinua vycházejí ze dvou základních metod popisu kontinua. První je Lagrangeův přístup, který používá jako referenční hodnoty počáteční polohu a čas. Druhým přístupem je Eulerův, který za referenční hodnoty považuje okamžité souřadnice. První bývá převážně používán v mechanice tuhých těles a druhý v mechanice tekutin. Na základě Langrangeova přístupu je možné odvodit Greenův tenzor deformace (také Green-Lagrangeův, Green-St. Venantův): 𝜀𝑗𝑘 𝐺=1 2(𝜕𝑢𝑗 𝜕𝑎𝑘+𝜕𝑢𝑘 𝜕𝑎𝑗+𝜕𝑢𝑖 𝜕𝑎𝑗𝜕𝑢𝑖 𝜕𝑎𝑘). (2.1) Z Eulerova přístupu vychází Almasiho tenzor deformace: 𝜀𝑗𝑘 𝐺=1 2(𝜕𝑢𝑗 𝜕𝑥𝑘+𝜕𝑢𝑘 𝜕𝑥𝑗+𝜕𝑢𝑖 𝜕𝑥𝑗𝜕𝑢𝑖 𝜕𝑥𝑘). (2.2) Dojde-li pouze k malým hodnotám posunutí a deformací, součiny parciálních derivací jsou řádově menší než první 2 členy v závorkách. V takovém případě je možné je zanedbat a vztahy přejdou v Cauchyho tenzor deformace. Dalších možností, jak vyjádřit deformaci je celá řada, např. inženýrská deformace, logaritmická deformace atd. 2.3.4. Tenzory napětí V lineárních případech se napětí počítá ze síly působící na plošku nedeformovaného tělesa. Při zahrnutí nelinearity v důsledku velkých deformací se musí zohlednit i změna geometrie tělesa. Prvním používaným tenzorem je Cauchyho tenzor napětí, který se vztahuje k okamžitému zatížení a geometrii tělesa. Tím se přidružuje k Almasiho tenzoru deformace. Greenovu tenzoru deformace je přidružen tzv. Langrangeův neboli 1. Piollův tenzor napětí. Má však tu nevýhodu, že není symetrický, proto se v praxi využívá 2. Piolla – Kirchhoffův tenzor napětí. Tento tenzor nemá fyzikální význam a je založen na působení fiktivní síly na těleso před deformací. Tenzor je invariantní vůči posuvům a rotacím tělesa jako tuhého celku. Pro posouzení napjatosti ve zkoumané konstrukci je třeba pracovat s Cauchyho tenzorem. Mezi oběma tenzory existují vztahy, které umožňují jednotlivé složky mezi sebou přepočítat. 2.4. Materiálové nelinearity V technické praxi se můžeme setkat s materiály s různým chováním. Lineární pružnost předpokládá platnost Hookeova zákona a elastické chování materiálu. Některé materiály, přestože jsou zatěžovány v elastické oblasti vykazují významně nelineární chování např. pryže nebo plasty. Chování se liší i v závislosti na charakteru zatížení. Viskoelastické materiály i při malých posuvech a deformacích mohou vykazovat jistou míru disipace energie závislou například na rychlosti zatěžování nebo frekvenci daného zatížení. Další běžnou nelinearitou může být překročení meze kluzu, kde dochází ke vzniku plastických deformací a materiál již nelze popsat prostým lineárním vztahem. Pro popis materiálových nelinearit slouží velké množství materiálových modelů, z nichž některé jsou popsány v tomto textu. 2.4.1. Lineárně elastické modely Lineárně elastický izotropní model materiálu je nejjednodušší materiálový model v mechanice pružných těles. Řídí se Hookovým zákonem ve všech směrech. Definován je pomocí 2 konstant – modulu pružnosti v tahu E a Poissonovým číslem μ. Je použitelný zejména pro kovové materiály. Pro vyhodnocení simulace je nutné ještě stanovit dovolenou hodnotu napětí. Lineárně elastický ortotropní model se obdobně jako předchozí model řídí Hookeovým zákonem, avšak hodnoty modulu pružnosti a Poissonova čísla nejsou ve všech směrech stejné. Pro výpočet je třeba zadat chování ve 3 osách, moduly pružnosti v tahu Ex, Ey, Ez, Poissonova čísla μxy, μxz, μyz 7 a moduly pružnosti ve smyku Gxy, Gxz a Gyz. Je také nutné znát dovolenou hodnotu napětí ve 3 směrech pro posouzení pevnosti. Tento model se může být použit pro dřevo, kosti nebo kompozitní materiály. Lineárně elastický anizotropní model na rozdíl od předchozích modelů nemá žádnou rovinu symetrie Pro úplnou definici je potřeba 36 materiálových konstant. Matice může být symetrická, v takovém případě stačí 21 nezávislých parametrů. Tento model je vhodný pro modelování vysoce anizotropních materiálů. 2.4.2. Hyperelastické materiálové modely Hyperelastické modely jsou vhodné pro materiály, které vykazují elastické (vratné) chování ve velkém rozsahu deformací. Deformace mohou dosahovat desítek až stovek procent a závislost napětí a deformace bývá silně nelineární. Typickými materiály jsou pryže, pěny, měkké biologické tkáně, polymery v kaučukovitém stavu. Materiály mohou mít rovněž významnou závislost chování na teplotě. U pryží se často předpokládá nestlačitelnost, kdy objemový modul je řádově vyšší než modul pružnosti při jednoosém zatížení. Se stlačitelností je potřeba uvažovat až v případě, kdy jsou omezeny tvarové změny např. stlačování těsnění do příliš malé drážky, kdy se těsnění opírá o všechny stěny drážky. Většina modelů využívá funkce hustoty deformační energie, ze které jsou následně počítána napětí. Při definici jsou většinou využity deformační invarianty založené na třech protaženích λ. 𝐼1=𝜆1 2+𝜆2 2+𝜆3 2, (2.3) 𝐼2=𝜆1 2𝜆2 2+𝜆2 2𝜆3 2+𝜆1 2𝜆3 2, (2.4) 𝐼3=𝜆1 2𝜆2 2𝜆3 2, (2.5) 𝜆=𝐿 𝐿0, (2.6) kde L je deformovaná délka vzorku, L0 je počáteční délka vzorku. Pro nestlačitelný materiál vychází třetí invariant nulový. Neo-Hooke model Tento materiálový model je nejjednodušším z hyperelastických modelů. Poskytuje dobré výsledky pro protažení do λ = 1,3 u gumových nebo biologických materiálů. Kvůli nedostatku parametrů neposkytuje přesný popis experimentálního chování a nedokáže popsat nárůst tuhosti při velkých deformacích, nehrozí však nesmyslné výsledky při simulaci jiných stavů, než na které byly nastaveny jeho parametry. Model je dán touto funkcí: 𝑊=𝐶10(𝐼1−3)+𝐾 2(𝐽−1)2, (2.7) kde C10 je konstanta definovaná jako polovina počátečního smykového modulu, K je počáteční objemový modul a J je poměr mezi deformovaným a původním objemem. Pro nestlačitelný materiál je J = 1 a rovnice přejde do tvaru: 𝑊=𝐶10(𝐼1−3). (2.8) Mooney-Rivlin model Mooney-Rivlinův model je rozšířením Neo-Hooke modelu o další parametry a invarianty. Tyto parametry umožňují přesnější popis experimentálních dat a přesnou predikci chování pro simulaci stejného charakteru zatížení jako v experimentu. S přibývajícími parametry však model ztrácí fyzikální základ a dochází k prostému proložení naměřených dat křivkou. Model je proto vhodné používat pouze pro simulace, kde zatížení odpovídá provedenému experimentu, jinému charakteru zatížení model nemusí odpovídat a v některých případech může dojít fyzikálně nesmyslným výsledkům. Pro 3 parametry má model tvar: 8 𝑊=𝐶10(𝐼1−3)+𝐶01(𝐼2−3)+𝐶11(𝐼1−3)(𝐼2−3)+𝐾 2(𝐽−1)2, (2.9) kde C10, C01 a C11 jsou konstanty modelu, které lze získat proložením naměřené závislosti mezi napětím a deformací křivkou. Yeoh model Yeohův model je založen na předchozím modelu, ale obsahuje pouze členy obsahující první invariant. Oproti Mooney-Rivlinovu modelu neposkytuje tak přesnou shodu s experimentálními daty zejména pro nízké hodnoty deformace. Pokud jsou parametry modelu kladné, model má monotónně zpevňující průběh, je tak možné jej použít i pro jiné než naměřené zátěžné stavy. Počet parametrů není přímo omezen a jeho rovnici je možné napsat jako: 𝑊=∑𝐶𝑖0(𝐼1−3)𝑖+∑1 𝑑𝑘(𝐽−1)2𝑘, 𝑁 𝑘=1 𝑛 𝑖=1 (2.10) kde Ci0 a dk jsou materiálové konstanty modelu. Arruda-Boyce model Tento model je založen na struktuře gum (polymerní řetězce), na kterých lze dosáhnout velmi dobrých výsledků. Pro nestlačitelné chování je definován pouze 2 parametry. Kvůli tomu neposkytuje tak přesný popis jako předchozí modely, ale je použitelný i pro jiné zatížení, než na které byl nastaven. Model je definován rovnicí: 𝑊=𝐶[1 2(𝐼1−3)+1 20𝜆𝐿 2(𝐼12−9)+11 1050𝜆𝐿 4(𝐼13−27)+19 7000𝜆𝐿 6(𝐼14−81) +519 673750𝜆𝐿 8(𝐼15−243)]+1 𝑑(𝐽2−1 2−ln𝐽), (2.11) kde C je počáteční smykový modul, λL je limitní protažení řetězců a d je parametr stlačitelnosti. Další modely: Výpočetní softwary nabízí celou řadu dalších materiálových modelů, patří sem například Extended Tube model, Gentův model, Henckyho model, Ogdenův model atd. Některé z modelů byly odvozeny přímo pro pěnové materiály např. Blatz-Ko foam model, Ogden hyperfoam model atd. Při použití hyperelastických materiálových modelů je potřeba vždy zvážit použitelnost modelu pro dané použití. Některé modely mohou výborně popisovat zatížení, na kterém byly nastaveny jejich parametry, pro jiný typ zatížení mohou však dávat nepřesné nebo přímo fyzikálně nesmyslné výsledky. 2.4.3. Viskoelastické materiály Viskoelastické materiály kombinují chování elastické a viskózní, které je časově závislé. Po aplikaci zatížení je deformace vlivem elasticity konstantní a viskózní deformace roste v čase. Modely mohou být použity buď pro simulace relaxace materiálu, kde po zatížení dojde k nárůstu napětí, které se vlivem viskózní složky v čase snižuje (relaxace) nebo pro simulace cyklického zatížení, kde vlivem viskozity dochází k disipaci energie. Oproti elastickému materiálu tak vzniká hysterezní smyčka ve vztahu mezi napětím a deformací. Základními modely jsou Maxwellův a Kelvin-Voigtův model, další modely pak vznikají jejich kombinací. Oba modely je možné si zjednodušeně představit jako kombinaci pružiny a tlumiče. Maxwellův model Tento model je možné si přestavit jako sériové zapojení pružiny s tlumičem. Po zatížení dochází ke skokové změně deformace vlivem elastické složky, při dalším působení zatížení narůstá viskózní složka deformace. Po odstranění zatížení mizí elastická složka deformace, avšak zůstává viskózní složka, která 9 je trvalá. Při simulaci relaxace dojde při skokové změně deformace ke skokovému nárůstu napětí (Obrázek 1), které v čase kompletně zrelaxuje. Model je popsán rovnicí: 𝑑𝜀𝑡 𝑑𝑡=𝑑𝜀𝑑 𝑑𝑡+𝑑𝜀𝑠 𝑑𝑡=𝜎 𝜂+1 𝐸𝑑𝜎 𝑑𝑡, (2.12) kde εS je elastická deformace (pružina), εd je viskózní složka deformace (tlumič), E je modul pružnosti materiálu a η je viskozita materiálu. Kelvin-Voigtův model Oproti předchozímu modelu si lze tento model představit jako paralelní zapojení pružiny a tlumiče. Na skokovou změnu zatížení nenastává skoková změna deformace, ale deformace je brzděna viskózní složkou. Po odstranění zatížení nezůstává žádná zbytková deformace (Obrázek 1). Model popisuje rovnice: 𝜎(𝑡)=𝐸𝜀(𝑡)+𝜂𝑑𝜀(𝑡) 𝑑𝑡 . (2.13) Obrázek 1: Odezva viskoelastických modelů, a) Maxwellův model, b) Kelvin-Voigtův model 2.4.4. Plastické chování materiálů Ve výpočtech se můžeme setkat s případy, kdy hodnota napětí přesáhne mez kluzu materiálu a dojde ke vzniku trvalé (plastické) deformace, které po odlehčení nezmizí. Vztah mezi napětím a deformací se stává nelineární. Rovněž je nutné znát podmínky zatěžování jako historii nebo rychlost. V některých případech není vliv rychlosti zatěžování významný a lze použít model, který je nezávislý na rychlosti. Naproti tomu například v analýzách tečení materiálu, kde plastická deformace roste s časem nebo naopak při analýzách nárazů bývá nezbytné využití časově závislého modelu. Hranici mezi pružným a nepružným chováním je vymezena plochou plasticity, která je popsána skalární funkcí – podmínkou plasticity. Při použití ideálně pružně plastického materiálu je plocha plasticity neměnná, kovy však v plastické oblasti vykazují zpevnění. Toto zpevnění můžeme rozdělit na 3 druhy: • izotropní, • kinematické, • kombinované. Při izotropním zpevnění se mění velikost plochy plasticity. Velikost plochy plasticity je dána poloměrem Y. Toto zpevnění je vhodné pro modelování monotónního zatěžování např. kování, tažení, ale nelze jím vhodně popsat cyklické namáhání. Toto zpevnění lze popsat jako: 𝑓=𝑓(𝜎)−𝑌=0. (2.14) 16 V praktických aplikacích se častěji než metoda dopředných diferencí používají jiné metody. V současnosti se z implicitních schémat používají prakticky pouze schémata navržená Hilberem, Hughesem a Taylorem. Vlastnosti implicitního algoritmu jsou: • Posuvy, rychlosti a zrychlení v čase tn+1 jsou vyjádřeny pomocí hodnot v časových krocích tn a tn+1. Z toho plyne název algoritmu – implicitní. • Schéma je nepodmíněně stabilní, stabilita není závislá na délce časového kroku. Stabilitou se rozumí udržet malé odchylky řešení, pokud se počáteční podmínky liší jen o malou hodnotu. Při nevhodné délce časového kroku může řešení dosahovat fyzikálně nesmyslných výsledků, výpočet však dále konverguje. Typické je zhroucení výpočtu během několika časových kroků. • Schéma lze použít i pro řešení stacionární úlohy, kdy se vypustí matice hmotnosti. • V každém časovém kroku je třeba řešit soustavu 5.9 a včetně dynamické matice tuhosti. Tomu se lze vyhnout v případě lineární úlohy s konstantním časovým krokem. Pro tento případ stačí dynamickou matici tuhosti sestavit v prvním kroku. Při použití implicitního algoritmu je snahou používat co největší časové kroky, což vyžaduje použití tenzorů velkých deformací při popisu kinematiky. To vede k nutnosti provádět v každém časovém kroku iterace pro splnění pohybové rovnice. Pro řešení se používá Newton-Raphsonova metoda. 5.2. Explicitní algoritmus Řešena je stejná pohybová rovnice jako v případě odvození implicitního algoritmu. V časovém kroku tn má rovnice tvar: [𝑀]{𝑄󰇘}𝑛+[𝐾]{𝑄}𝑛={𝐹}𝑛. (5.10) Pomocí metody centrálních diferencí je zrychlení možné vyjádřit jako funkci posuvů: {𝑄󰇘}𝑛={𝑄}𝑛+1−2{𝑄}𝑛+{𝑄}𝑛−1 ∆𝑡2. (5.11) Po dosazení a úpravě pohybové rovnice platí: ([𝑀] ∆𝑡2){𝑄}𝑛+1={𝐹}𝑛−[𝐾]{𝑄}𝑛+[𝑀]2{𝑄}𝑛−{𝑄}𝑛−1 ∆𝑡2. (5.12) Vlastnosti explicitního algoritmu jsou: • Posuvy v čase tn+1 jsou získány z pohybové rovnice pro předchozí časový krok. Z toho plyne název algoritmu – explicitní. • Algoritmus je podmíněně stabilní, podmínkou stability je dostatečně malý časový krok, jehož velikost musí být menší, než je délka kritického časového kroku tc. Kritický časový krok odpovídá délce průchodu napěťové vlny nejmenším elementem. Jeho velikost je možné určit z charakteristického rozměru elementu (různý pro různé typy elementů) h. • Nelze zanedbat matici hmotnosti, algoritmus tak není přímo použitelný pro řešení stacionárních úloh. Toto omezení lze obejít umělým zvýšením hustoty materiálu. Vznikají tak sice nereálné setrvačné síly a kinetická energie, ty jsou ale řádově nižší vůči přetvárným silám a deformační energii. • Výhoda explicitního algoritmu nastává při použití diagonální matice hmotnosti. Soustava rovnic se rozpadne na samostatné nezávislé rovnice. Z každé rovnice lze přímo vyjádřit neznámou a ∆𝑡<∆𝑡𝑐, (5.13) ∆𝑡𝑐=ℎ √𝐸 𝜌 . (5.14) 17 není nutné obtížně sestavovat dynamickou matici hmotnosti jako u implicitního algoritmu, to umožňuje výrazně rychlejší řešení rovnic. • Vlivem integrace lineárních prvků pouze v jednom bodě může při diskretizaci čtyřúhelníkovými nebo šestistěnnými prvky vznikat tzv. hourglass energie. Prvky se deformují do tvaru přesýpacích hodin bez vnějšího zatížení, které by toto chování způsobovalo. To může vést k nereálným výsledkům nebo zhroucení výpočtu vlivem přidané energie. Předejít tomu lze zmenšením velikosti elementu, některou z metod tlumení hourglass energie, použitím trojúhelníkových nebo čtyřstěnných elementů nebo plnou integrací elementů, tu však nenabízí všechny řešiče. 5.3. Porovnání algoritmů Z popisu algoritmů vyplývá, že každý je ze svého principu vhodný pro jiné použití. Volbou vhodného algoritmu je možné značně zkrátit časovou náročnost analýzy. Implicitní algoritmus je z důvodu možnosti využití velmi dlouhých časových kroků vhodný pro pomalejší a dlouhotrvající děje s menší mírou nelinearity, kde mezi jednotlivými časovými kroky nedochází k velkým změnám. To mohou být např. simulace relaxace materiálu, pohyb mechanismů nebo předepínání konstrukcí. Explicitní algoritmus je kvůli malým časovým krokům vhodný zejména pro krátkodobější děje s velkou mírou nelinearit. Typickým použitím jsou simulace nárazů, explozí, průstřelů. V těchto simulacích bývají zahrnuty materiálové nelinearity v podobě velkých plastických deformací nebo s vlivem rychlosti deformace na chování materiálu. Mezi tělesy může být značné množství kontaktů. S přibývajícími nelinearitami významně neroste čas výpočtu. Algoritmus je možné použít i pro simulace tváření materiálů, které bývají silně nelineární z důvodu velké plastické deformace a užití kontaktů. 18 6. Ztráta stability Vzpěr neboli ztráta stability je geometrická nestabilita, kdy dochází k rychlému vybočení konstrukce vlivem zatížení. Při vybočení dojde k náhlé ztrátě tuhosti konstrukce. Ve většině případů k němu dochází vlivem zatížení v tlaku, např. při sešlápnutí plechovky nebo tlakovém zatěžování štíhlých nosníků. Na vzpěr by tam měly být kontrolovány konstrukce se štíhlými částmi zatíženými tlakem. Ztráta stability může mít globální charakter, kdy ztratí stabilitu celá konstrukce nebo lokální charakter, kdy dojde ke ztrátě stability pouze v určité části konstrukce a zbytek může dále přenášet zatížení. V případě globální ztráty stability konstrukce není dále schopna zatížení přenášet a dochází zpravidla k destrukci. Vzpěr nebývá většinou ovlivněn setrvačnými účinky, ke změně deformační energie na kinetickou dochází až po ztrátě stability. Pro řešení může tak být použita pouze část základní rovnice MKP pro stacionární úlohy. [𝐾]{𝑄}={𝐹}. (6.1) Ve většině případů je úkolem zjistit bezpečnost vůči ztrátě stability, případně limitní (kritickou) sílu, kterou může být konstrukce zatížena. Řešení je možné provést obdobně jako u modální analýzy převedením na problém vlastních čísel: ([𝐾]+𝜆𝑖[𝐾𝐺]){𝜓𝑖}={0}, (6.2) kde [K] je matice tuhosti, λi zátěžný faktor udávající poměr mezi kritickým zatížením, které způsobí ztrátu stability a aplikovaným zatížením, [KG] napěťová matice tuhosti závislá na zatížení konstrukce a {ψi} tvar módu ztráty stability. V případě symetrického profilu dochází v momentě ztráty stability k nejednoznačnému řešení, kdy není možné určit, na kterou stranu dojde k vybočení (Obrázek 5). Řešení zde má tzv. bifurkační bod, ve kterém dojde k rozdvojení zátěžné cesty. To při řešení způsobuje 2 stejné výsledky. Obrázek 5: Kritická síla při vzpěru K řešení ztráty stability může být využita lineární analýza, nelineární analýza nebo transientní analýza. V lineární analýze dojde k vybočení a ztrátě tuhosti při dosažení kritické síly, v případě nelineární analýzy dochází k vybočení ještě před dosažením kritické síly. Výsledky získané lineární analýzou jsou velmi nekonzervativní, zátěžný faktor neboli bezpečnost vůči ztrátě stability by tak měl být dostatečně vysoký podle typu a použití konstrukce. Z analýzy lze vyhodnotit pouze bezpečnost vůči ztrátě stability a její tvar, nelze vyhodnotit posunutí, napětí nebo deformace. Obrázek 6: Analýza ztráty stability tvaru 19 7. Teplotní úlohy Teplotní úlohy jsou velmi rozšířenou oblastí výpočtů využívajících MKP. Cílem analýz bývá určení rozložení teplotního pole, určení teplotních gradientů, tepelných toků uvnitř a na hranicích těles a celkové teplené bilance soustavy. Z pohledu času je možné tepelné úlohy rozdělit na stacionární a nestacionární. Nestacionárními analýzami mohou být náběhy tepelných zařízení, rozběhy motorů, ohřev nebo chlazení materiálu. Ve stacionárních úlohách jsou analyzovány ustálené stavy konstrukcí, např. ustálený chod motoru nebo tepelné ztráty budov. Tepelné analýzy mohou být kombinovány s napěťově-deformačními analýzami. Teplotní pole může být použito pro výpočet napětí v materiálu, který vlivem teploty dilatuje. Opačně lze zkoumat např. tváření, kde se velká část deformační energie mění v teplo. Teplo se může šířit 3 způsoby: • Vedením (kondukcí) – teplo si předávají částice, které spolu sousedí. • Prouděním (konvekcí) – teplo se přenáší pohybem tekutin (plynu nebo kapalin), teplá proudící tekutina teplo předává a chladnější odebírá. Proudění tepla často nastává v kombinaci s vedením tepla. Proudění se dělí na volné, které vzniká vlivem různé hustoty teplé a studené tekutiny. Tekutina s nižší hustotou pak stoupá a s vyšší hustotou klesá. Druhou možností je nucené proudění, kdy je tekutina poháněna ventilátorem nebo čerpadlem. • Zářením (sáláním, radiací) – teplo se přenáší elektromagnetickým zářením. Tento přenos probíhá i ve vakuu. Časově závislý problém vedení tepla ve 3D lze popsat diferenciální rovnicí: 𝑘(𝜕2𝑇 𝜕𝑥2+𝜕2𝑇 𝜕𝑦2+𝜕2𝑇 𝜕𝑧2)+𝑄=𝜌𝑐𝜕𝑇 𝜕𝑡, (7.1) kde k je součinitel tepelné vodivosti, T je teplota, Q tepelný výkon, ρ hustota materiálu, c měrná tepelná kapacita a t čas. V případě řešení stacionárního problému je pravá strana nulová. Zkráceně lze rovnici zapsat jako: 𝑘∇2𝑇+𝑄=𝜌𝑐𝜕𝑇 𝜕𝑡. (7.2) Množství tepla, které se prostředím šíří nebo předává je možné definovat měrným tepelným tokem. Pro vedení tepla je závislý na součiniteli tepelné vodivosti a gradientu teploty. 𝑞=−𝑘∇𝑇. (7.3) Teplo, které se předává mezi tekutinou a stěnou je závislé na rozdílu teplot stěny a tekutiny a součiniteli přestupu tepla α. Jeho velikost se může významně změnit v závislosti na teplotě, typu povrchu nebo proudění tekutiny. Vztah pro tepelný tok mezi tekutinou a stěnou je: 𝑞=𝛼∆𝑇. (7.4) V aplikacích, kde se dosahuje velmi vysokých teplot (tváření za tepla) může být významné i sdílení tepla zářením. Tepelný tok je pak závislý na emisivitě tělesa ε, Stefan-Boltzmanově konstantě σ a rozdílu čtvrtých mocnin teploty zářiče T a okolního prostředí T0. 𝑞=𝜀𝜎(𝑇4−𝑇04). (7.5) Zahrnutím radiace se problém stává nelineární a vzhledem ke čtvrté mocnině u teplot je možný vnik významné chyby v důsledku nepřesné teploty zářiče. Rovnici vedení tepla v diskrétní podobě je možné definovat podobně jako základní rovnici MKP pro řešení deformace a napjatosti. [𝐶𝑇]{𝑈𝑇󰇗}+[𝐾𝑇]{𝑈𝑇}={𝐹𝑇}, (7.6) kde CT je globální matice tepelné vodivosti, KT globální matice tepelné kapacity a FT vektor tepelného zatížení a UT jsou neznámé teploty v uzlech. 20 Použitá literatura [1] Fusek, M., R. Halama. 2011. MKP a MHP. Ostrava: Vysoká škola báňská – Technická univerzita Ostrava, Západočeská univerzita v Plzni. [2] Ondrouch, J., J. Podešva. 2012. Aplikovaný mechanik jako součást týmů konstruktérů a vývojářů: teorie a příklady k předmětu Technické kmitání. Ostrava: Vysoká škola báňská – Technická univerzita Ostrava. Dostupné také z: http://projekty.fs.vsb.cz/147/ucebniopory/978-80-248-27629.pdf. [3] Visser, G. Modal Analysis: What It Is and Is Not [online]. [cit. 2024-11-25]. Dostupné z: https://simteq.co.za/blog/modal-analysis/. SIMTEQ Engineering (Pty) Ltd. [4] ANSYS, INC. 2021. Buckling: Local Buckling vs Buckling, Methods for Solving [online]. [cit. 202411-25]. Dostupné z: https://innovationspace.ansys.com/courses/wpcontent/uploads/2021/01/Lesson3-Buckling.pdf. [5] Schauer, P. 2006. Aplikovaná fyzika: Přenos tepla. Brno: Vysoké učení technické v Brně. [6] Petruška, J. 2007. Inovace předmětů studijních programů strojního inženýrství v oblasti teplotního namáhání: Řešení úloh vedení tepla a teplotní napjatosti pomocí MKP. Brno: Vysoké učení technické v Brně. [7] ANSYS, INC. 2020. Fundamentals of Damping [online]. [cit. 2024-11-25]. Dostupné z: https://innovationspace.ansys.com/courses/wp-content/uploads/sites/5/2020/12/2.6.2Fundamentals-of-Damping-New-Template.pdf. Transformace formy a obsahu vysokoškolského vzdělávání na VŠB-TUO NPO_VŠB-TUO_MSMT16605/2022 17. LISTOPADU 2172/15 708 00 OSTRAVA-PORUBA univ[email protected]z www.vsb.cz Toto dílo je licencováno pod CC BY 4.0