UČEBNÍ TEXTY VYSOKÝCH ŠKOL Vysoké učení technické v Brně Fakulta strojního inženýrství Doc. RNDr. Libor Čermák, CSc. Algoritmy metody konečných prvků Předmluva k revidovanému elektronickému vydání Tato skripta jsou určena studentům matematického inženýrství FSI VUT v Brně pro studium předmětu Numerické metody III. Jedná se o revidovaný text stejně pojmenovaných (již vyprodaných) skript, která vyšla v PC-DIR Real, s.r.o., Brno, v roce 2000. Brno, září 2005 Libor Čermák 3 Předmluva Metoda konečných prvků (MKP) je univerzálním nástrojem pro efektivní řešení rozmanitých inženýrských problémů, které vyžadují provádět výpočty v oborech „teoretických" inženýrských disciplín jako je pružnost a pevnost, termomechanika, hydromechanika atp. K masivnímu rozšíření MKP došlo začátkem 70-tých let. Od té doby zůstává MKP v popředí zájmu inženýrů, matematiků i fyziků. Problémům spojeným s MKP je věnováno ohromné množství publikací, vznikají rozsáhlé programové systémy, konají se četné konference. Krása MKP spočívá v tom, že její podstatu lze vysvětlit na jednoduchém matematicky formulovaném modelovém problému. To je také hlavním cílem výuky MKP v mezioborovém studiu matematického inženýrství na FSI VUT a MU v Brně. Zatímco učební text [30] Prof. RNDr. Alexandra Ženíška, DrSc. objasňuje matematické základy MKP, tento „doplňující" text se věnuje především „technologii" MKP, tj. poměrně detailnímu popisu „realizačních" algoritmů. Učební text je rozčleněn do tří kapitol zabývajících se algoritmizací MKP postupně v jedné, dvou a třech prostorových proměnných. Jedna prostorová proměnná zdaleka neumožňuje MKP ukázat svou „plnou sílu". Kapitola první má proto význam především pedagogický: prostřednictvím okrajového diferenciálního problému druhého řádu (tah-tlak osově namáhaného prutu) a čtvrtého řádu (ohyb nosníku podle Kirchhoffovy teorie) jsou na jednoduchém definičním oboru (úsečka) objasněny všechny podstatné kroky MKP-algoritmizace: přechod od formulace klasické (tj. diferenciám' rovnice a okrajové podmínky) k formulaci slabé, konečněprvková aproximace, elementární matice a vektory, sestavení soustavy rovnic a zpracování jejího řešení. Těžištěm skripta je kapitola druhá, tedy úlohy rovinné. Modelový okrajový problém druhého řádu (typu stacionárního vedení tepla) je diskretizován nejdříve užitím lineárního trojúhelníkového prvku. Podrobně jsou vysvětleny datové struktury popisující triangulaci, značná pozornost je věnována sestavovacím algoritmům. Uvedena je také algoritmizace dvou nestacionárních úloh (vedení tepla, kmitání membrány), stručně je zmíněna rovněž úloha vlastních čísel. Dis-kretizace dvourozměrné úlohy pružnosti MKP demonstruje schopnost této metody vypořádat se s problémem popsaným soustavou dvou parciálních diferenciálních rovnic (Laméovy rovnice). Následuje popis nejpoužívanějších izoparametrických trojúhelníkových a čtyřúhelníkových konečných prvků včetně jejich aplikace na řešení úlohy typu stacionárního vedení tepla. Problémy reálného světa jsou nelineární, což tento text zohledňuje zařazením popisu metod pro řešení nelineárních úloh diskretizovaných MKP (například pro úlohy typu nelineárního vedení tepla, stacionárního i nestacionárního). Efektivní řešení konvekčně- difúzni ch úloh s dominantní kon-vekcí (KDUsDK) patří k nejobtížnějším úlohám numerické matematiky a MKP zde teprve v poslední době začíná konkurovat zatím více používané metodě konečných objemů (MKO). V tomto textu je popsána: pro stacionární KDUsDK upwind technika se současným použitím jak MKO tak MKP a pro nestacionární KDUsDK metoda charakteristik v kombinaci s MKP. Poslední kapitola uvádí třírozměrné izoparametrické konečné prvky (čtyřstěny, pětistěny a šestistěny) a jejich použití demonstruje opět na problému typu stacionárního vedení tepla. Kromě již zmíněného učebního textu [30] lze z česky psaných materiálů doporučit k hlubšímu studiu knihy [10], [27], [22], [4], [28], [11], [23], [21] a skripta [14], [15]. Z programových systémů MKP jmenujme alespoň ANSYS [1] a NEXIS [18] českého původu. Pro výukové účely je velmi vhodný PDE toolbox MATLABu [12]. Za chyby a přepisy, které se bohužel ve skriptu jistě vyskytnou, se dopředu omlouvám. Budu vděčný všem čtenářům, kteří mě na ně upozorní. Brno, září 2000 Libor Čermák 4 Obsah 1 Jednorozměrné úlohy 6 Základní pojmy a označení ................................. 6 1.1 Okrajový problém pro 0DR2............................. 6 a) Klasická formulace ................................. 6 b) Slabá formulace................................... 8 c) Metoda konečných prvků.............................. 11 1.2 Okrajový problém pro ODR4............................. 17 a) Klasická formulace ................................. 17 b) Slabá formulace................................... 18 c) Metoda konečných prvků.............................. 20 2 Rovinné úlohy 25 2.1 Základní pojmy a označení ........... ................... 25 2.2 Klasická formulace................................... 26 2.3 Greenova formule.................................... 28 2.4 Slabá formulace..................................... 28 2.5 Triangulace, po částech lineární funkce........................ 29 2.6 Diskrétní slabá formulace............................... 32 2.7 Elementární matice a vektory............................. 34 a) Elementární matice a vektor na elementu e.................... 34 b) Elementární matice a vektor na straně S..................... 37 c) Sestavení globální matice a vektoru........................ 38 2.8 Několik poznámek ................................... 44 2.9 Minimalizační formulace................................ 47 2.10 Nestacionární úloha vedení tepla........................... 48 2.11 Dynamika........................................ 51 2.12 Rovinná napjatost a rovinná deformace ....................... 54 a) Klasická formulace ................................. 54 b) Slabá formulace................................... 57 c) Diskrétní slabá formulace.............................. 58 d) Elementární matice a vektory........................... 59 e) Sestavení globální matice a vektoru........................ 62 f) Závěrečné poznámky ................................ 63 2.13 Izoparametrické prvky................................. 63 2.14 Nelineární úlohy .................................... 76 a) Stacionární úloha.................................. 76 b) Nestacionární úloha................................. 82 2.15 Konvekčně-difúzní úlohy s dominantní konvekcí................... 83 a) Stacionární úloha, upwind metoda......................... 83 b) Nestacionární úloha, metoda charakteristik.................... 90 3 Prostorové úlohy 95 Literatura 106 5 1. Jednorozměrné úlohy Základní pojmy a označení Prostor funkcí spojitých v intervalu (0,£) označíme C(0,£). Podobně C(0,£), C(0,£) a C(0,£) postupně označují prostory spojitých funkcí v intervalech (0,£), (0,£) a (0,£). Dále Ck(0,£), Ck(0,£), Ck(0,£), Ck(0,£) označují prostory funkcí, které jsou v uvažovaných intervalech spojité a mají v nich spojité derivace až do řádu k včetně. Prostor funkcí po částech spojitých v intervalu (0,£) označíme PC(0,£). Přitom / G PC(0,£) znamená, že existuje dělení 0 = to < t\ < ■ ■ ■ < tn = £ takové, že f E C(íj_i,íj) a že existují konečné jednostranné limity \imx^t._1+f (x) a \imx^t._ f (x), i = l,...,n. Symbolem PCk(0,£) značíme prostor funkcí, které jsou v intervalu (0, £) spojité spolu se svými derivacemi až do řádu k — l včetně, a jejichž k-tá derivace je v intervalu (0, £) po částech spojitá. Lebesgueův prostor funkcí integrovatelných s kvadrátem v intervalu (0,£) označíme L/2(0,£). Symbolem Hk(0,£) označíme Sobolevův prostor funkcí / takových, že /, /',..., f^ G L/2(0,£). Je přirozené položit H°(0,£) = L2(0,£). Zřejmě PCk(0,£) C Hk(0,£). Je známo, že Hk{0,£) C Ck~ľ(0,£), a proto se nedopustíme žádné velké nepřesnosti, když si pod funkcí z prostoru Hk(0,£) budeme představovat funkci z prostoru PCk(0,£). 1.1. Okrajový problém pro ODR2 a) Klasická formulace Naším cílem je nalézt funkci vyhovující diferenciální rovnici a splňující v bodě x = 0 buďto okrajovou podmínku u(0) = g0 nebo okrajovou podmínku p(0K(0) = a0u(0) - f30 (1.3b) a v bodě x = £ buďto okrajovou podmínku u(£) = g e (1.4a) nebo okrajovou podmínku ■p(£)u'(£) = aeu(£) - pe. (1.4b) 6 Okrajové podmínky (1.3a) a (1.4a) se nazývají Díríchletovy, okrajová podmínka (1.3b) resp. (1.4b) se pro a0 = 0 resp. cti = 0 nazývá Neumannova a pro a0 ^ 0 resp. a£ ^ 0 Newtonova. Úloha (1.1)-(1.4) může popisovat například problém tahu-tlaku prutu, tedy prutu namáhaného pouze tahem popřípadě tlakem. V tom případě je u(x) posunutí střednicové čáry prutu, p(x) = E{x)A{x), kde E(x) je Youngův modul pružnosti a A(x) je plocha průřezu prutu, q(x) je měrný odpor podloží, na němž prut spočívá, f(x) je intenzita zatížení, g0, gi jsou předepsaná posunutí konců prutu, a0, ae Jsou tuhosti pružin v koncových bodech prutu a /3q, fit jsou zadané síly působících na koncích prutu. Tatáž úloha popisuje také průhyb tenkého prutu, někdy se také hovoří o průhybu struny. Jinou aplikací popsanou formálně stejnými rovnicemi a okrajovými podmínkami je například stacionárni úloha vedeni tepla v tyčí. Pak u(x) je teplota, p(x) je koeficient tepelné vodivosti, q(x) je koeficient přestupu tepla „povrchem" tyče do obklopujícího prostředí, f(x) = fz + que je součtem tepelného zdroje fz(x) a tepelného toku que {ue(x) je teplota obklopující „povrch" tyče), go, gt jsou předepsané teploty okrajů tyče, «o, ae jsou koeficienty přestupu tepla okraji tyče a /3q, fit jsou zadané tepelné toky na okrajích tyče. Rovnice (1.2) pak má tvar -ipu')' + q(u - ue) = fz. Newtonovy okrajové podmínky se v úloze vedení tepla obvykle píší ve tvaru p(0K(0) = aoMO) - Uq), -p(í)u'(í) = at{u{í) - ue), kde Uq, Ui jsou teploty okolí konců tyče. Úloha (1.1)-(1.4) může být modelem i pro jiné problémy, například z oblasti potenciálního proudění tekutin, elektrostatického potenciálu atd. Funkci u G C2(0,£), vyhovující rovnici (1.2) a splňující jednu z okrajových podmínek (1.3) a jednu z okrajových podmínek (1.4), nazveme klasickým řešením úlohy (1.1)- (1.4) . Pro existenci klasického řešení je třeba předpokládat dostatečnou hladkost dat, tedy splnění „podmínek hladkosti" p e C^OJ), q, f e C(0,£). (1.5a) Dále budeme v souladu s fyzikálním významem dat předpokládat splnění „fyzikálních podmínek" p(x) >p0>0, q(x) > 0 v (0,£), a0 > 0, ctl > 0. (1.5b) Pro jednoznačnost řešení je třeba ještě připojit „podmínky uložení", a sice že je splněna alespoň jedna z následujících tří podmínek : a) platí buďto okrajová podmínka (1.3a) nebo (1.4a); b) q(x) > go > 0 na části intervalu (0,£); c) «o > 0 nebo ae > 0. Jsou-li tedy splněny podmínky (1.5), má úloha (1.1)-(1.4) jediné řešení. 7 Poznámka 1. Není-li splněna žádná z podmínek (1.5c), je úloha (1.1)-(1.4) tvaru -(pu')' = f, p(0K(0) = -a), -p(l)u'(l) = -pt (1.6) a nazývá se Neumannova úloha. Její řešení je ľ z(x) r u(x) = C + / dx, kde z(x) = p(x)u'(x) = — /3q — / f(s)ds (1-7) J p{x) Jo a kde C je libovolná konstanta. Jestliže —z{£) = —p{£)u'{£) = —j3e, tedy je-li splněna podmínka rovnováhy + & + / /(s)ds = 0, (1.8) ./o má Neumannova úloha (1.7) nekonečně mnoho řešení navzájem se lišících o konstantu C. Pokud však podmínka rovnováhy (1.8) splněna není, Neumannova úloha (1.6) řešení vůbec nemá. □ b) Slabá formulace Funkcí v G Cľ(0,£) nazveme testovací, je-li rovna nule v tom krajním bodě intervalu (0,£), v němž je předepsána Dirichletova okrajová podmínka. Pro konkrétnost se omezíme na okrajové podmínky (1.3a) a (1.4b), takže testovací funkce v splňuje v(0) = 0. Násobme rovnici (1.2) testovací funkcí v a integrujme přes (0,£). Integrací per-partes členu Jg [—(pu')']v dx a následným užitím okrajové podmínky (1.4b) a vztahu v(0) = 0 obdržíme pí pí fvdx= / \—{pu')' + qu] v dx = — pu'v\x_Q + / \pu'v' + quv] dx = o Jo x~ Jo ľ = [ct£u(í) — fii\v{€) + / [pu'v' + quv] dx. Jo Odvodili jsme tedy, že řešení u úlohy (1.2), (1.3a), (1.4b) musí splňovat kromě Dirichletovy okrajové podmínky u(0) = Qq také rovnost [pu'v' + quv] dx + aeu(£)v(£) = í fvdx + (3ev(£) (1.9) Jo pro každou funkci v G Cľ(0,£), v(0) = 0. Okrajová podmínka (1.4b) Newtonova typu, která se stala součástí integrální rovnice (1.9) a je tak automaticky splněna, se nazývá přirozenou okrajovou podmínkou. Dirichletovu okrajovou podmínku (1.3a), která součástí rovnice (1.9) není a jejíž explicitní splnění proto musíme vyžadovat, nazýváme podstatnou nebo také hlavní okrajovou podmínkou. Rovnice (1.9) je dobře definována i v případě, kdy funkce u a v jsou z prostoru X = iJ1(0,£). Testovací funkce pak volíme z prostoru V = {v G X\ v(0) = 0} testovacích funkcí a řešení u z množiny W = {v G X\ v(0) = Qq} přípustných řešení. Dále označíme a(u,v) = I [pu'v' + quv] dx + a£u(£)v(£), 0 (1.10) ŕ Jo Pak úlohu najít u G W splňující a(u,v) = L(v) \/v G V (1.11) nazýváme slabou (nebo také variační) formulací problému (1.2), (1.3a), (1.4b). Řešení úlohy (1.11) nazveme slabým řešením. Zeslabené podmínky hladkosti p,q,f ePC(0,i) (1.5a') spolu s předpoklady (1.5b) a (1.5c) zaručují jednoznačnou existenci slabého řešeni. Slabá formulace má v úloze tahu-tlaku prutu význam principu virtuálních posunutí a samotné testovací funkce v E V mají význam virtuálních posunutí Su přípustných řešení u G W. Formulace problému založená na principu virtuálních posunutí je obecnější než formulace popsaná diferenciální rovnicí. Je tomu tak proto, že „diferenciální formulaci" (1.2), (1.3a), (1.4b) lze při platnosti podmínek (1.5a) a pro u G C2(0,£) odvodit z formulace slabé postupem opačným k tomu, jímž jsme získali rovnost (1.11). Poznámka 2. Uveďme si tvar V, W, a(u,v) a L (v) pro všechny možné kombinace okrajových podmínek. (aa) Okrajové podmínky (1.3a), (1.4a) (1.12aa) V = {veX\ v(0) = v(£) = 0}, W = {v G X | v(0) = g0, v(£) pí pí a(u,v)= / [pu'v' + quv] dx, L{v) = i fvdx. Jo Jo (ab) Okrajové podmínky (1.3a), (1.4b) V = {v G X | v(0) = 0}, W = {v G X | v(0) = g0}, pí pí a(u, v) = [pu'v' + quv] dx + a£u(£)v(£), L{v) = j fvdx + fítv(í). (ba) Okrajové podmínky (1.3b), (1.4a) V = {v G X | v(£) = 0}, W = {v G X | v(£) = ge}, a(u, v) = [pu'v' + quv] dx + a0u(0)v(0), L (v) = I fv dx + (30v(0) . (1.12ab) o (1.12ba) o (bb) Okrajové podmínky (1.3b), (1.4b) V = X, W = X, a(u, v) = I [pu'v' + quv] dx + a0u(0)v(0) + a^u{£)v{£), L{v)= / f v dx + /30v(0) + pev(£). Jo (1.12bb) Ve všech těchto čtyřech případech je zaručena jednoznačná existence slabého řešení úlohy (1.11). □ 9 Poznámka 3. Jak jsme již uvedli, v úloze tahu-tlaku prutu má fy resp. fy význam síly a «o resp. on má význam tuhosti pružiny působící v levém resp. pravém konci prutu. Ve slabé formulaci je účinek sil reprezentován jejich virtuální prací fyv(0) resp. fyv(í) jako součást virtuální práce L (v) vnějších sil a vliv pružin je zohledněn započtením virtuálních prací aou(0)v(0) resp. aeu(£)v(£) do virtuální práce a(u, v) vnitřních sil. Je proto přirozené očekávat, že síla fy působící ve vnitřním bodě c přispěje k virtuální práci vnějších sil prací fyv(c) a pružina tuhosti ac umístěná ve vnitřním bodě c přispěje k virtuální práci vnitřních sil prací acu(c)v(c). Ukažme si, že tomu tak skutečně je. Bodové zatížení fy lze považovat za idealizované vyjádření spojitého zatížení intenzity fe(x) = fy/2e působícího na úseku (c — e,c + s), kde e je velmi malé číslo. Pak totiž jQ f£ dx = fy a současně užitím věty o střední hodnotě pro v G C(0,£) dostaneme lime^0 Jq feV(^x = fyv(c). Zcela analogicky lze bodový účinek reprezentovaný pružinou o tuhosti ac považovat za idealizované spojité pružné uložení na úseku (c — e, c + e), na kterém intenzita pružného odporu qe(x) = ac/2e. Pak jQ qe dx = ac a současně lime^0 f0 qeuvdx = olcu(c)v(c). Jsou-li fy síly a on tuhosti pružin působících v bodech q, můžeme a(u, v) a L (v) zapsat ve tvaru a(u, v) = [pu'v' + quv] dx + ceiu(ci)v(ci), Jo i (1-13) L(v)= I fvdx + S^fyv(ci), Jo i přičemž v závislosti na okrajových podmínkách, v souladu s (1.12), jsou do ^ q;íií(q)í;(q) zahrnuty také členy aQu(0)v(0), a^u{l)v{£) a do X^A'Kq) také členy fyv(0), fyv(l). Příslušné slabé řešení úlohy (1.11) existuje a je určeno jednoznačně. Úloha (1.11) pro a(u, v), L(v) podle (1.13) aV,W podle (1.12) má význam i pro jiné aplikace než je úloha tahu-tlaku prutu. Tak například v úloze vedení tepla reprezentuje fy — qíj-u(cj) intenzitu bodového tepelného zdroje působícího v bodě q. Výše zmíněnou úlohu lze však uvažovat také zcela abstraktně, bez vztahu ke konkrétní technické aplikaci. Přitom z předpokladu q > 0, viz (1.5b), dostáváme pro zaručení existence slabého řešení požadavek a,i > 0 ve všech bodech q, pro jednoznačnost slabého řešení stačí v případě okrajových podmínek (1.3b), (1.4b) a pro q = 0 žádat > 0 alespoň v jednom bodě q. Pro úlohu (1.13) tedy dostáváme fyzikální podmínky p(x) > po > 0, q(x) > 0 v (0,£), on > 0 pro všechna i, (1.5b') a podmínky uložení vyžadující, že je splněna alespoň jedna z následujících tří podmínek : a) platí buďto okrajová podmínka (1.3a) nebo (1.4a); (1.5c') b) q{x) > qo > 0 na části intervalu (0,£); c) olí > 0 pro nějaké i. □ 10 Poznámka 4. Dirichletovu okrajovou podmínku u(0) = Qq resp. u(l) = ge lze přibližně realizovat pomocí Newtonovy okrajové podmínky p(0K(0) = a0[u(0) - g0] resp. - p{i)u'{i) = ae[u(i) - ge] pro velké a0 resp. cti .Tuto skutečnost si ilustrujme na dvou příkladech. a) V úloze tahu-tlaku prutu se nejčastěji setkáváme s homogenní Dirichletovou okrajovou podmínkou u(0) = 0 resp. u(l) = 0. Je zřejmé, že nulové posunutí lze zajistit připojením příslušného konce prutu k velmi tuhé pružině, což odpovídá velké hodnotě «0 resp. . b) V úloze vedení tepla má a0 resp. cti význam koeficientu přestupu tepla. Čím větší je hodnota tohoto koeficientu, tím více se teplota u(0) resp. u(l) přiblíží k teplotě okolí, kterou reprezentuje právě qq resp. ge. Ukazuje se tedy, že při řešení praktických úloh plně vystačíme s okrajovými podmínkami (1.3b) a (1.4b), takže se můžeme zabývat úlohou (1.11) pouze pro V = W = X (1.14) a pro a(u,v) a L (v) určené rovnicemi (1.13). Jsou-li splněny podmínky (1.5a'), (1.5b') a (1.5c':b-c), má úloha (1.11), (1.13), (1.14) právě jedno slabé řešení. □ c) Metoda konečných prvků Na intervalu (0,£) zvolíme dělení 0 = xq < x\ ■ ■ ■ < = t a na každé úsečce (xí-i,Xí) délky hi = Xí — Xí-i hledáme přibližné řešení U(x) ve tvaru lineárního polynomu procházejícího body (xí-i, Ui-i) a (xí,Uí), takže />< . _ />< />< _ />> . U(x) = Ui-iWi-^x) + UíWí(x), kde Wi-i(x) Wi(x) hi hi Funkce U (x) je tedy na celém intervalu (0, l) po částech lineárni funkcí^ určenou předpisem n U{x) = Y^U%Wl{x), (1.15) i=0 kde Wi(x) se jsou tzv. bázové funkce, lineární na každé úsečce (xk-i,Xk) a takové, že 1 pro i = j, 0 pro i ý j- Wi(Xj) 0 = Xq Xi+1 xn-i xn = £ Obr. 1.1. Lineární Lagrangeovy bázové funkce 11 Úsečku (xk-i,Xk), na které je definována lineární funkce určená svými hodnotami v uzlech Xk-i a Xk, nazýváme lineárním Lagrangeovým konečným prvkem nebo také elementem. Nechť h = max! /\ 7, ^ " ^-1 1 ^ " ^-1 / (pU Wi) = fc^—j—- - = PH—j— , r+1(pí/y) = (----) ( - — ) = p,l+i- hi+i ) V hi+i) 2 hi+i ľ(qUwi) = ^hilqi-^+Ui-! ■ 0 + q^U ■ 1] = ^h^Ui, ľ+1(qUwi) = \hi+1[qi+Ui ■ 1 + qi+i-Ui+1 ■ 0] = \hi+1qi+Ui, ŕ(fwi) = lMU-i,+ ■ o + /ť_ • i] = ifci/i-, r+1(M) = • i + /i+i,- • o] = \hi+1fi+. Dosadíme-li odtud do (1.18), obdržíme Uí — Ui-i Uí — Ui+1 i Pi-i-7-+Pí+k-Z-+ žihiqi-+ hi+1qi+\Ui + oiiUi = 2 ht +2 hi+1 (1.21-i) = \ [h"ifi- + h>i+ifi+] + /3j pro i = 1,..., N — 1. Podobně odvodíme piU°~Ul + iMo^o + a0í/0 = ^i/o+ + Ä (1-21-0) pro i = 0 a Pn-±un , UN 1 + \hNqN-UN + ajv#jv = \hNfN_ + /3_/v (1.21-N) pro i = N. Soustava rovnic (1.21) má symetrickou pozitivně definitní a tedy regulární matici. Řešení Uq, ..., U^ soustavy rovnic (1.21) určuje MKP-slabé řešení U (x). Pro chybu u — U a její derivaci platí za předpokladu u G C2(0,£) odhad dj max max I-—r(u — U)\ = 0(h2~3) pro j = 0,1. (1.22) li(x), Q\ = ©j_i, 0*2 = &i, A[ = Aj_i, A2 = Ať. Pak na (x^Xi) je r(Pu'v') = ŕ([e[w[ + e>*]'p [A\w\ + Aiwi]') = [ 1 2) \r([wi]'P[w\]') r([wi]'P[wi]'] AI Al kde ©' ©i ©2 kde Kt2 = \hi Podobně odvodíme ľ(qUv) = [@J]TKí2A\ položíme K* = K*1 + K*2, a dále odvodíme P(fv) = [&]TF\ kde P = \hi ífy'+ a A1 Qi-i,+ 0 0 qt- Al 15 Z rovnice 0 =ah(U,v) - Lh(v) N N N N í=l y^ľ(pU'v' + qUv) + y^jaiU{xi)v{xi _i=l i=0 N N =Yy&Yvěx - fí+y 0*na* - ä] i=l i=0 a rovnice (1.23) dostaneme rovnost N N 0T [KA - F] = ^[©^[KW - F] + y, 0i[«iAť - a i=0 (1.25) i=0 z níž plyne postup, jak pomocí elementárních matic K* = {k%rs}^ s=l, elementárních vektorů F* = (F{,Fl)T a čísel [3i sestavit globálni matici K a globálni vektor F: stačí srovnat členy se stejnými indexy u parametrů 0 a A (pro určení prvků matice K) nebo jen u parametru 0 (pro určení prvků vektoru F) na levé a na pravé straně rovnice (1.25). Soustava rovnic (1.24) se symetrickou třídiagonální maticí je tvaru / a0 b0 0 0 ... b0 ai bi 0 ... 0 bi a2 b2 ... \ / A0 \ Ai A2 V bN-2 &n-1 b^-i 0 &Ar_i aN J &n-1 \ An J ( Fo \ Fl F2 Fn-i \FN J a pro její nenulové koeficienty odvodíme ao = &íi + "o, 22 aN = k22 + aN, k i+l 11 ctí, b b0 = kl2, Fo = Fl + (3Q, pi+i Fn = F2 + (3n- U'í+1 rp _ rpi ^12 ) ri — r2 Pi, i 1,...,N Snadno ověříme, že jsme opravdu dostali rovnice (1.21). Jestliže chceme pro q = 0 získat U{x,i) = u{x,i) přesné, viz poznámka 6, musíme při výpočtu F* integrovat přesně. Pro f'l{x) = fl na elementu konstantní dostaneme a pro lineární fl(x) = f'l(xi-i)w\(x) + f% (xí)w2(x) obdržíme ť_ ! Í2f\xl_l) + f\xi)\ ~^l\ř{xl_1) + 2f\xi) I 16 1.2. Okrajový problém pro ODR4 a) Klasická formulace Hledáme funkci u{x) G C4(0, £), (1.26) která vyhovuje diferenciální rovnici (pu")" - (ru')' + qu = f v(0,£). (1.27) Pokud jde o okrajové podmínky, nabízí se nám více možností než pro rovnici druhého řádu. V každém z krajních bodů c = 0 a c = l vybereme dvě z následujících čtyř podmínek u(c) = uc, (1.28.1) u'(c) = po > 0, r(x) > 0, q(x) > 0 v (0,£), a0 > 0, /30 > 0, ae > 0, /3e > 0. (1.30b) Pro jednoznačnost řešení je třeba ještě připojit „podmínky uložení". Existuje celá řada rovnocenných podmínek uložení, uveďme si některé z nich: a) v jednom z koncových bodů jsou předepsány buďto podmínky (1.29.1) nebo podmínky (1.29.2) s 7C > 0 nebo podmínky (1.29.3) s ac > 0 nebo podmínky (1.29.4) s ac > 0, 7C > 0; b) v obou koncových bodech jsou předepsány buďto podmínky (1.29.2) nebo podmínky (1.29.4) s ac > 0; c) v jednom koncovém bodě jsou předepsány buďto podmínky (1.30c) (1.29.2) nebo podmínky (1.29.4) s ac > 0 a ve druhém koncovém bodě jsou předepsány buďto podmínky (1.29.3) nebo podmínky (1.29.4) s 7c > 0; d) na části intervalu (0,£) platí současně r(x) > r0 > 0, q{x) > qo > 0. Jsou-li tedy splněny podmínky (1.30), má úloha (1.26) - (1.29) jediné řešení. b) Slabá formulace Funkci v G C2(0, £) nazveme testovací, jestliže v(c) = 0 v tom krajním bodě c, v němž je předepsána podmínka (1.28.1), a jestliže v'(d) = 0 v tom krajním bodě d, v němž je předepsána podmínka podmínka (1.28.2). Abychom byli konkrétní, zvolíme si okrajové podmínky u(0) = uQ, = 3, jak plyne z definice a,h(U,v) a {wí(x)} obrázky 2 a 3), a při dostatečně přesné numerické integraci je rovněž pozitivně definitní a tedy regulární. Vyřešením soustavy rovnic KA (1.42) získáme hledané parametry a2j = U(xí) a a2j+1 = U'(xí), i = 0,..., N. Matici K a vektor F sestavíme pomocí elementárních matic K* a elementárních vektorů F* příslušných elementům (xí_i,Xí). MKP-slabé řešení U a testovací funkce v je na elementu (xí-i,Xí) tvaru 4 4 U{x) = v{x) = YJ&jw)(x), 21 kde a\ = u(Xi_i) = a2i-2 , ©í = V{xi_l) = @2i-2 , ai = í/'(^i-l) = a2l-i, ©2 = ViXi-i.) = 02í_! , a\ = u (xi) = a2t, ei = v(xi) = e2,t, a\ = u'(xí) = a2,1+1 , e\ = v'(xí) = e2,t+1 jsou parametry a w{(x) = ň1(j~*i-1'SJ pro ^(0 = 1-3^ + 2^, w\(x) = hlň2(^~*i-1^ pro ň2(o = š-2e+e, wl{x) = Ň3 (^ř1) pro Mo = ^e-2e, wi(X) = kň, (^^j pro ň4(o = -e+e jsou „elementární" bázové funkce. Vyjádříme 4 4 ľ(pU"v") = r(£e>*]'v £>[«;*]'') = [&]Va\ 3=1 1=1 kde K*1 = {fcj}}Jí=1 pro = r([w}]"p Kf) a kde & = (01, 02, 0*, &\)t, A* = (ai, a2, a*, a\)t. Podobně vyjádříme 4 4 r(rW) = r(^0,k]v^ azkť) = [eť]rKť2Ať, j=i i=i kde K*2 = }JÍ=1 pro A;}2 = r(k]VK]'), 4 4 r(g^) = i'£ejffl;?5] aiw\) = [e*]TKí3A\ i=i z=i kde K*3 = {fcj?}},i=1 pro k* = Ť(w) qwf), položíme K* = K*1 + K*2 + K*3, a dále vyjádříme 4 r(^) = r(^0>}/) = [eí]TFíl, kde F^F^F^F^Fff pro Ff = r(w)f), i2 ľ{mv') = r(^e}[w}]'m) = [@*]TF kde Fí2 = (Ff,Ff,Ff,Ff)T pro F]2 = r([w^]'m) 22 a položíme F1 = F*1 + Fl2. Jsou-li p(x) = p% r(x) = r% = g% /(x) = m{x) = m2 i) konstantní, snadným výpočtem obdržíme K p hl K i2 30hi f 12 Qhi -12 Ah2 -Qhi 2h2 -12 -6/íí 12 -Qhi v 6hi 2fc? -Qhi Ah2 J f 36 3hi 36 3ht \ 3/ij Ah2 -3hi -h2 i 36 -3hi 36 -3hi V 3hi -h2 -3hi Ah2 ) K i3 420 156 22fcť 22/1,- 4/i? V 54 13^ -3h2 54 156 -22fc,- -13/1, \ 4ä? / / 6 \ f% 12 yi2 m í "I \ 0 1 V o/ Jsou-li na elementu (x,í_i,Xí) funkce /(x) = a m(x) ( 21 f (x\) + 9f(x2) \ hi 6 FJ1 = ^ 60 ^[3f (4) + 2/ť (4)] 9/i(4) + 21/i(4) V -^[2f (4) + 3f (4)] / F" = — 12 m1 (x) lineární, pak / —6rrí(x\) — 6ml(x2) \ rrŕ(x\) — rrŕ(x2) 6rrí (x\) + 6ml(x2) \ — rrŕ(x\) + ml(x'2) J Obecně však stačí integrály počítat jen přibližně formulí, která je přesná pro polynomy třetího stupně. Lze proto použít Simpsonovu formuli P{g) = ^[/(^-i) + Af(xt_i) + f{xi)} nebo Gaussovu formuli i\g) = hhilfixii - ^hi) + f(Xii + fhi)], kde x- i \hi. Z rovnic (1.41) a (1.39) užitím výše uvedených vyjádření členů P(pU"v" + rU'v' + qUv) a P(fv + mv') dostaneme rovnost 0T[KA - F] = n n = ^[0i]T[KiAi - Fť] + Y,{Q2í[(XíA2í - A] + e2í+1[7íA2í+1 Si]}, (1.43) i=0 23 z níž plyne postup, jak pomocí elementárních matic K* = {kzrs}fs=1, elementárních vektorů F* = {F*}f=l a čísel fy, 7«, Si sestavit globální matici K a globální vektor F: stačí srovnat členy se stejnými indexy u parametrů 0 a A (pro určení prvků matice K) nebo jen u parametru 0 (pro určení prvků vektoru F) na levé a pravé straně rovnice (1.43). Následuje algoritmus sestavení matice K a vektoru F. 1) Matici K a vektor F vynulujeme. 2) Pro každý element (xí_i,Xí), i = 1,...,N, definujeme čtyři čísla Q(l) = 2« — 2, Q(2) = 2i - 1, Q(3) = 2i, Q(4) = 2i + 1 a provedeme: a) klrs pro r, s = 1, 2, 3,4 přičteme k prvku matice K s indexy [Q(r), Q(s)]; b) fl pro r = 1, 2, 3,4 přičteme k prvku vektoru F s indexem Q (r). 3) Pro každý uzel xÍ7 i = 0,..., N, provedeme: a) on přičteme k prvku matice K s indexy [2i, 2i]; b) 7j přičteme k prvku matice K s indexy [2i + 1, 2i + 1]; c) fy přičteme k prvku vektoru F s indexem 2i; d) Si přičteme k prvku vektoru F s indexem 2i + 1. Protože matice tuhosti K je symetrická a má nenulové koeficienty jen v hlavní diagonále a ve třech sousedních horních a dolních subdiagonálách, stačí sestavovat prvky kij pro j = i,...,max(i + 3,2iV+l). Pro chybu u — U a její derivace platí za předpokladu u G C4(0,£) odhad dj max max I-—t {u — U)\ = 0{lr~^) pro j = 0,1,2,3. l 1 celé; b) oblasti fži,..., ÍŽjv takové, že Q = {jf=1 Í2« H Q j = 0 pro i ý Ji c) funkce fu ..., fN, f,t e C(í)j), / = ft v í)ť, i = 1,..., N. Prostor všech po částech spojitých funkcí budeme značit PC(Q). Symbolem PCk(Q) označíme prostor všech funkcí / takových, že f e Ck~ľ(Q) a Dk f e PC (Q), kde Dk f označuje k-té derivace funkce /. Lebesgueův prostor funkcí integrovatelných s kvadrátem v Q označíme L2(f2). Sobolevův prostor všech funkcí / takových, že /, Df,..., Dkf e L2(f2), označíme Hk(íl). Je-li S část hranice, pak prostor funkcí spojitých na S označíme C(S) a prostor funkcí po částech spojitých na S označíme PC(S). Připomeňme, že PC(íl) C L2(íl) a PCk(íl) C Hk(íl). Proto funkce z prostoru PC(íl) může posloužit jako příklad funkce z prostoru L2(Í2) a podobně funkce z prostoru PC1(Í2) jako příklad funkce z prostoru Hľ(íl). 2.2. Klasická formulace Naším cílem je určit funkci u(x,y) e c2 (Q) n c1 (Q u r2) n c (H) (2.1) vyhovující diferenciální rovnici -V • [pVíí] + qu = f víl (2.2) a splňující okrajové podmínky u = g na Ti, (2-3) —p— = au — 13 na T2. (2.4) on Tečka ve vztahu (2.2) značí skalární součin a protože _ . í d d \ T . ^r^n í d d \ T (du du\T V = grad= I—,— , je - V • [pVu\ = - —, — = \ox oyJ \ox oyJ \ox oyJ ô ( du \ 0 f du \ Tx{PTx)-Ty{PTy)^-{pU^-{pU^ O hranici V předpokládejme r = T\ U Y2, T\ H T2 = 0. Dále n = (nx,ny)T = (ni,ri2)T je jeánotkový vektor vnější normály hranice a du du du — = n • \7u = nx— + ny — dn dx dy 26 je derivace ve směru vnější normály. Okrajová podmínka (2.3) se nazývá Díríchletova, okrajová podmínka (2.4) se pro a = 0 nazývá Neumannova a pro q/0 Newtonova. Úloha může popisovat například problém dvourozměrného stacionárního vedení tepla. V tom případě je u(x,y) teplota, p(x,y) koeficient tepelné vodivosti, q(x,y) je koeficient přestupu tepla „plochou" Q, f(x,y) = Q + que je součtem tepelného zdroje Q(x,y) a tepelného toku que „plochou" Q (ue(x,y) je teplota okolí „plochy" Q), g(x,y) předepsaná teplota na hranici Ti, a(x,y) koeficient přestupu tepla a j3(x,y) zadaný tepelný tok na hranici T2. Rovnice (2.2) pak má tvar -V • [p Vu] + q(u -ue) = Q v a (2.2') Newtonova okrajová podmínka se v úloze vedení tepla obvykle píše ve tvaru du -p— = a(u - u0) na T2, (2.4 ) on kde ua(x,y) je teplota okolí (hranice r2). Pokud tutéž úlohu interpretujeme jako úlohu průhybu membrány, je u(x,y) výchylka, p(x,y) reprezentuje tuhost membrány, q(x,y) odpor prostředí, f(x,y) zatížení, g(x,y) předepsaný pokles podepřené části Ti hranice, a(x, y) tuhost pružinového uložení a /3(x, y) zatížení na části T2 hranice. Úloha (2.1)-(2.4) může být modelem i pro další problémy, například pro potenciální proudění tekutin, elektrostatický potenciál, kroucení tyče, difúzni šíření příměsi atd. Funkce u(x,y) s vlastností (2.1) splňující (2.2)-(2.4) se nazývá klasické řešení. Pro jeho existenci je třeba předpokládat dostatečnou hladkost dat, a sice pe^fň), f,qeC(ň), geCif,), a,/3eC(Ť2). (2.5a) Dále budeme v souladu s obvyklým fyzikálním významem dat předpokládat p(x,y) > Po > 0, q(x,y) > 0 v Q, a(x,y) > 0 na T2. (2.5b) Pro jednoznačnost řešení je třeba dále připojit podmínku mesiTi > 0 nebo Qo > 0 na Qq C Q, mes2f2o > 0, nebo (2.5c) o:(x,y) > a0 > 0 na r20 C T2, nieslo > 0, kde mesiTi resp. mesir2o je délka části Ti resp. T2o hranice a mes2f2o je plocha podoblasti Qq- Pro existenci řešení musíme předpokládat ještě něco navíc, například že Q je regulárni oblast s hladkou hranicí a T = Ti, nebo že T = To, nebo že ^J 6 ' ' (2.5d) íl je konvexní polygon a T = Ti. Soubor podmínek (2.5a) - (2.5d) nám zaručuje jednoznačnou existenci klasického řešení problému (2.1)-(2.4). 27 Poznámka 12. Nesplnění podmínky (2.5c) znamená řešit úlohu Ou —V • [pVíí] = / v fž, — p— = —/3 na T. on Pomocí Greenovy formule (2.6) lze ukázat, že tato úloha má řešení jen tehdy, platí-li f dxdy + / /3ds = 0. n Jt V tom případě je řešení nekonečně mnoho a navzájem se liší o konstantu. □ 2.3. Greenova formule Nechť f(x1,x2), g(x1,x2) G H1 (Q). Potom platí / ^-gdx= í fgmás- í f^-dx. (2.6) Jn clxi Jv Jn clxi Zde dx = dxi dx2, ds je diferenciál oblouku T a rti je i-tá složka jednotkového vektoru vnější normály hranice. 2.4. Slabá formulace Nechť v G C1 (Q). Vynásobením rovnice (2.2) funkcí v, integrací přes Q a použitím Greenovy formule dostaneme V • [p Vu]v dxdy = p[uxnx + uyny]v ds + / p[uxvx + uyvy] dxdy = (2-7) Jn = — í p—— vds+ í [pVu ■ Vv]dxdy = í (f — qu)vdxdy. Jt vn Jn Jn Funkci v(x,y) G C1(Í2) nazveme testovací, splňuje-li v = 0 na Ti. Dosadíme-li do (2.7) za v testovací funkci a užijeme (2.4), obdržíme rovnici / [pVw • + quv] dxdy + / auv ds = / f v dxdy + (3v ds. Jn Jt2 Jn Jt2 Označme si a(u, v) = [p Vm • + quv] dxdy + / auv ds, Jn Jt2 (2.8) (2.9) L{v)= / j"vdxdy + / /3vds, (2.10) 'n Jt2 položme X = H1 (Q) a dále označme V = {veX\v = 0 na TJ, W = {veX\v = g na rj. (2.11) a(u,v) je symetrická bilineární forma a L (v) je lineární funkcionál. V nazveme prostorem testovacích funkcí a W nazveme množinou přípustných řešení. Zdůrazněme, že W není lineárním prostorem funkcí, neboť pro g ^ 0 a v±,V2 G W neplatí v\ + v2 G W. Slabou (variační) formulací úlohy (2.1)-(2.4) rozumíme úlohu najít u G W splňující a(u,v) = L(v) \/v G V. (2.12) Funkci u vyhovující (2.12) nazýváme slabým řešením (úlohy (2.1)-(2.4)). Jednoznačnou existenci slabého řešení lze zaručit za výrazně slabších předpokladů než byly předpoklady (2.5a) - (2.5d), které jsme potřebovali pro jednoznačnou existenci řešení klasického. Pokud jde o hladkost dat, stačí místo (2.5a) předpokládat p,q,fePC(Q), geCiľJ, a,f3ePC(T2). (2.5a') Tyto požadavky jsou realistické a odpovídají tomu, s čím se setkáváme při řešení praktických úloh. Předpoklady (2.5b) a (2.5c) ponecháme v platnosti. To není omezující, neboť tyto předpoklady jsou v praktických úlohách splněny. Konečně předpoklad (2.5d) nahradíme přirozeným požadavkem Q je regulární oblast s hranicí po částech hladkou, (2.5ď) který opět plně vyhovuje praktickým potřebám. Shrňme tedy, že jednoznačná existence slabého řešení je zaručena, jsou-li splněny předpoklady (2.5a'), (2.5b), (2.5c) a (2.5ď). 2.5. Triangulace, po částech lineární funkce Předpokládejme, že Q je polygon. Q pokryjeme triangulaci7 skládající se z trojúhelníkových elementů e takových, že uzávěry každých dvou různých trojúhelníků jsou buďto disjunktní nebo mají společný vrchol nebo stranu a že n = y e. Trojúhelníky triangulace budeme nazývat také prvky. Vrcholy trojúhelníků budeme nazývat uzly. Tringulaci zvolíme tak, aby body průniku r\ H T2 byly uzly. Počet všech prvků triangulace označíme PP, počet všech uzlů PU, počet uzlů ležících na Ti označíme PB a počet zbývajících uzlů (ležících uvnitř Q a na T2) označíme PN. Strany elementů budeme značit S, množinu všech stran ležících na hranici T2 označíme S a počet všech stran S G S označíme PS. Nejdelší stranu trojúhelníků triangulace označíme h. Symbol h budeme v dalším používat jako měřítko jemnosti triangulace. Kromě toho, použijeme-li písmeno h jako index u nějaké veličiny vyznačíme tím, že tato veličina závisí na triangulaci T. Triangulaci popíšeme pomocí následujících souborů dat: 1. PRVKY[i,j], i = 1,..., PP, j = 1, 2, 3, jsou čísla uzlů elementu q G T; 29 2. STRANY[i,j], i = 1,..., PS, j = 1, 2, jsou čísla uzlů strany Si G S; 3. BODY[i], i = 1,..., PB, jsou čísla uzlů ležících na 1\; 4. X[i], Y[i], i = 1,..., PU, jsou rr-ové a y-ové souřadnice uzlů. 12 14 13 18 Obr. 6. TViangulovaná polygonálni oblast Z čistě formálních důvodů, umožňujících snadnější matematické vyjadřování, budeme v dalším textu předpokládat, že na hranici I\ leží uzly s čísly PN+1,..., PU. Pokud tomu tak není, jako třeba v příkladu podle obrázku 6, snadno si poradíme : uzly přečíslujeme. Můžeme k tomu použít kódovací tabulku KC, která ke každému uzlu I přiřadí tak zvané kódové číslo K C [I]: {Vytvoření pole kódových čísel KC[1..PU]} for I:=l to PU do KC[I]=0; for I:=l to PB do KC[BODY[I]]=-I; J:=0; for I:=l to PU do if KC[I] = 0 then begin J:=J+1; KC[I]:=J end; Vidíme, že uzly neležící na I\ mají přiřazena kódová čísla 1,..., P N a uzly ležící na I\ čísla —1,..., —PB. Je-li KC[I] > 0, je novým číslem uzlu I číslo KC[I], pro K C [I] < 0 je novým číslem uzlu I číslo P N — K C [I]. Funkci, která je na každém trojúhelníku e G T lineární a globálně, tj. na Q, spojitá, nazveme funkcí po částech lineární. Každá taková funkce v(x,y) je jednoznačně určena svými hodnotami v{xi,yi) ve vrcholech Pj(xj,í/j) triangulace T. Prostor všech po částech lineárních funkcí označíme X^. Funkce z X^ jsou v Q spojité a mají v Q po částech spojité derivace, tedy Xh C PC1 (Q) C X. 30 PRVKY STRANY BODY 1 4 5 1 5 2 2 5 6 2 6 3 3 6 7 4 8 9 4 9 5 5 9 10 5 10 6 6 10 11 6 11 7 7 11 12 13 9 8 13 14 9 14 10 9 14 15 10 15 11 10 15 16 11 16 12 11 18 14 13 18 19 14 19 15 14 19 20 15 15 20 21 15 21 16 21 22 16 22 17 16 22 23 17 1 2 3 4 5 6 7 8 9 10 11 12 1 2 2 3 3 7 7 12 12 16 16 17 17 23 23 22 22 21 21 20 20 19 19 18 13 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 X Y 0 0 0 0,5 0 1 0,5 0 0,5 0,5 0,5 1 0,5 1,5 1 0 1 0,5 1 1 1 1,5 1 2 1,5 0 1,5 0,5 1,5 1 1,5 1,5 1,5 1,75 2 0 2 0,5 2 1 1,75 1,25 1,75 1,5 1,75 1,75 KC 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: Tab. 1: Okrajové podmínky 31 Speciálním případem funkcí z Xh jsou funkce Wi(x,y), které jsou v Pí rovny jedné a v ostatních vrcholech jsou rovny nule. Těmto funkcím se říká bázové. Každá funkce v G Xh může být pomocí svých hodnot v uzlech a pomocí bázových funkcí vyjádřena ve tvaru PU v(x,y) = ^2v(xi,yi)wi(x,y). (2.13) i=i Xft je prostor dimenze PU. Podstatnou a charakteristickou vlastností metody konečných prvků je ta skutečnost, že bázové funkce v ní používané jsou nenulové jen na velmi malé poddoblasti oblasti Q. Říkáme, že bázové funkce Wi mají malý nosič. (Nosičem funkce ip je uzávěr množiny všech bodů, v nichž je funkce ip nenulová). Dále definujme prostor testovacích funkcí Vh = {veXh\ v(P3) = 0 VPj G Tl} (2.14) a množinu přípustných funkcí Wh = {veXh\ v{Pj) = g{P]) VP, G ľx}. (2.15) Pro v G Vh platí PN v(x,y) = ^2v(xi,yi)wi(x,y) (2.16) i=i a pro v G Wh je PN PU v(x,y) = ^2v(xi,y,i)wi(x,y) + ^ 9(xí,yi)wÁx,y)- (2-17) i=l í=PN+l Vh je podprostor prostoru V a dimenze Vh je rovna PN. Množina přípustných řešení Wh pro g ý 0 opět zřejmě není lineárním prostorem a Wh je podmnožinou W pouze tehdy, když pro každou funkci v G Wh je v = g na Ti. 2.6. Diskrétní slabá formulace Označme Ie(g) numericky nebo přesně spočtený je g(x, y) áxáy a Is(g) numericky nebo přesně spočtený js g(x, y) ás. {Ie{g) a Is{g) jsou zřejmě lineární funkcionály). Pak diskrétní slabá (variační) formulace úlohy (2.1)-(2.4) zní najít U G Wh splňující a,h(U,v) = Lh(v) Wv G Vh, (2.18) 32 kde ah(U,v) = ^P(pVU ■ Vv + qUv) + ^Is(aUv), Lh(v) = j2nfv) + J2lS^- ses (2.19) (2.20) 2eT Ses Hodnotu MKP-slabého řešení U v uzlu Pj označíme A j = U (x j, y j) a hodnotu testovací funkce v v uzlu Pi označíme 0« = v(xi,yi). Podle (2.16) a (2.17) je PN PU U(x,y) = ^AjWj(x,y) + ^ g{xj,yj)wj{x,y), j=l j=PN+l PN (2.21) i=l Všimněte si, že '^2^=1AjWj(x,y) G Vh a z = J2^=PN+i9(xjjyj)wj(x^y) e Wh- Vzhledem k bilineárnosti formy a,h(U,v) a lineárnosti funkcionálu Lh(v) je ^p^7 0 =ah(U,v) - Lh(v) PN PN i=l kde jsme si označili i=i {j=i pn r pn pi/ Lh(wi)- ^ ah{wj,wi)g(xj,yj) j=PN+l Q1 ÍKA - Fl (2.22) pi/ hj = a>h(wj,Wi), Fi = Lh(wi)- ^ ah(wj,wi)g(xj,yj), j=PN+l K = {%}gl1, F = (F1,...,FPW)T, A = (A1;..., AP7V)T, 0 = (©!,..., 0p^)t. Protože @ je libovolný vektor, musí platit KA = F. (2.23) (2.24) Z této soustavy lineárních rovnic spočítáme neznámé Ai,..., Ap^. K se nazývá matice tuhosti a F vektor zatížení. Uveďme některé vlastnosti matice tuhosti. 1. K je symetrická: ah(u>j,Wi) = ah(u>i,Wj); 2. K je řídká: ah(u>j,Wi) = 0, pokud uzly Pj a Pí nejsou vrcholy téhož trojúhelníka; 33 3. K je pozitivně definitní: © K© = a,h(v,v); integrujeme-li přesně, je a,h(v,v) = / [pVíí • + qv2] dxdy + / av2 ds > 0 pro í)/0 Jn Jt2 a dále a,h(v,v) = 0 právě když v = 0 (stačí použít(2.5b) a (2.5c)); matice K však bude pozitivně definitní pro dostatečně malé h i tehdy, když integrujeme numericky užitím formulí, které uvedeme v následujícím odstavci 2.7; 4. Při vhodném očíslování proměnných je K pásová matice. Soustavu rovnic (2.24) lze bez problémů vyřešit: pásovost a pozitivní definitnost matice K umožňuje efektivní nasazení některé z modifikací Gaussovy eliminační metody bez výběru hlavních prvků, řídkost a pozitivní definitnost matice K zase vybízí k řešení soustavy rovnic výkonnou iterační metodou. 2.7. Elementární matice a vektory Matici tuhosti K a vektor zatížení F sestavíme pomocí elementárních matic Ke a elementárních vektorů Fe pro e G T a elementárních matic K.s a elementárních vektorů Fs pro S G S. Matici K budeme nazývat také globální maticí tuhosti a vektor F globálním vektorem zatížení. a) Elementární matice a vektor na elementu e Uvažme jeden konkrétní element e triangulace T s vrcholy P^(x^,y^), i = 1,2,3. Pro uzel Pf je i lokálním číslem uzlu na elementu e. Je-li c = c(e) číslo elementu e, pak globálním číslem uzlu Pf je číslo I = PRVKY[c, i]. Pf a Pj jsou tedy jen různá označení téhož uzlu. Nechť wf(x,y) je restrikce bázové funkce wj(x,y) na element e. Označme si w1(x, y) = a\x + b\y + 4, i = 1, 2, 3, Ne = Ne(x,y) = (wl(x,y),w%(x,y),wl(x,y)) Ae = (A^, a2, Ag)T, kde A,e = U(xlyt), i % = 1,2,3 (2.25) (2.26) (2.27) @e = (Ql,Qe2,Qe3)T, kde 0? = 1 = 1,2,3. Pak slabé řešení U a testovací funkci v lze na elementu e vyjádřit ve tvaru U = NeAe v = Ne@e. (2.29) Dále označme / dw\ dw2 dw% \ (2.30) 34 Pak V£/ = BeAe a X7v = Be@e. Nyní již vyjádříme ľ(pVU ■ Vv + qUv) = ľ([Vv]TpVU + vqU) = = Je([0e]T[Be]TpBeAe + [0f [NfgNeAe) = [0e]TKeAe kde Ke = Kel + Ke2 a Kel = /e([BfpBe), Ke2 = F([NfgNe). Při výpočtu prvků matice Kel se obvykle používá formule ľ(g) =pl(e)g(xeT,yeT), (2.31) (2.32) (2.33) (2.34) kde pl(e) je plocha elementu e a x^ = \{x\ +^2 + 2:3), y^ = \{Vi + Jsou souřadnice těžiště elementu e. Formule (2.34) integruje přesně polynomy prvního stupně. Pomocí této formule snadno určíme K el pl(e)p(xeT,yeT) \ 3al + ala3 + b\b\ \ &2&1 CÍ2a2 + ^2^2 a2fl3 + ^2^3 OgO^ + ^3^2 a3a3 + ^3^3 J (2.35) Zbývá spočítat koeficienty af, 6^ vlastností bázových funkcí 1,2,3, a vyjádřit pl(e). K tomu účelu stačí využít 0 pro i ý 3i 1 pro i = j, hJ 1,2,3. (2.36) Tak dostaneme / x\ yf 1 \ V V2 vi ( a\ ď2 a% \ °l °2 \ "3 C3 / / 1 0 0 \ 0 1 o 0 0 1 a pomocí Cramerova pravidla zjistíme = (vl - ye3)/de, b\ = (x% - xe2)/de, r.e _ (x2ye3 - ytxt)/de, ď2 = (%e - yí)/de, b% = (x\ - xe3)/de, c2 — (x\y\ - yixl)/de, a% = foí - yl)/de, b% = (%2 - xl)/de, ), (2.55) kde iy = \{xf + x'2), y^ = \{yf + ) jsou souřadnice středu strany S. Pak dostaneme Ks = idsa(4,ž/r) 1 1 , Fs = ids/3(4,^)( n ) (2.56) ) ) ), ¥S = \dsp{xsT,ysT)(^li a matice K.s je plná. Někdy se při výpočtu matice K.s používá speciální postup, Ks = /S([NS]TqNs) = a(x%,y%) f [N5]TN5ds, J s který vede opět na plnou matici Ks = \dsa{xsT^)(y2i (2.57) c) Sestavení globální matice a vektoru Zkombinujeme-li (2.22), (2.19), (2.20), (2.32), (2.48), (2.51) a (2.52) vidíme, že pro slabé řešení U a libovolnou testovací funkci v platí 0 = ah(U, v) - Lh{v) = 0T [KA - F] = = ]!T[ee]T [KeAe - fe] + Y}®sľ íksas - FT • (2'58) eeT ses Z této rovnosti lze odvodit pravidla pro sestavení globální matice K a vektoru F z lokálních matic Ke, K.s a lokálních vektorů Fe, Fs. Uvedeme si dva postupy. Algoritmus 1 (eliminační) Víme už, že uzly neležící na Ti mají přiřazena kódová čísla 1,..., PN a uzly ležící na r\ čísla —1,..., —PB. Předpokládejme, že předepsané hodnoty řešení v uzlech BODY[I] jsou uloženy v pozici G[I] pole G[1..PB]. Následuje popis sestavovacího algoritmu. {Nulování globální matice K ~ K[1..PN, 1..PN] a vektoru F ~ F[1..PN]} for I:=l to PN do begin for J:=l to PN do K[I,J]=0; F[I]:=0 end; {Zpracování příspěvků z elementů e G T} for E: = l to PP do {prvky} begin {Výpočet elementární matice Ke ~ KE[1..3,1..3] a vektoru Fe ~ FPfl.^] } for I:=l to 3 do {řádky} 38 begin II:=KC[PRVKY[E,I]]; if II > 0 then {rovnice se sestavuje} begin for J:=l to 3 do {sloupce} begin JJ:=KC[PRVKY[E,J]]; if J J > 0 then K[II,JJ]:=K[II,JJ]+KE[I,J] {příspěvek do matice K} else F[II]:=F[II]-KE[I,J]*G[-JJ] {příspěvek do vektoru F} end; {sloupce} F[II]:=F[II]+FE[I] {příspěvek do vektoru F} end {sestavování rovnice} end {řádky} end; {prvky} { Zpracování příspěvků ze stran S G S} for S:=l to PS do {strany} begin {Výpočet elementární matice Ks ~ KS[1..2,1..2] a vektoru F5 ~ FS[1..2]} for I:=l to 2 do {řádky} begin II:=KC[STRANY[S,I]]; if II > 0 then {rovnice se sestavuje} begin for J:=l to 2 do {sloupce} begin JJ:=KC[STRANY[S,J]]; if JJ > 0 then K[II,JJ]:=K[II,JJ]+KS[I,J] {příspěvek do matice K} else F[II]:=F[II]-KS[I,J]*G[-JJ] {příspěvek do vektoru F} end; {sloupce} F[II]:=F[II]+FS[I] {příspěvek do vektoru F} end {řádky} end {sestavování rovnice} end; {strany} Předpokládejme, že řešení soustavy rovnic s maticí soustavy K[1..PN, 1..PN] a vektorem pravé strany F[1..PN] dostaneme v poli Z[1..PN]. Pak pole D[1..PU] hodnot slabého řešení v uzlech dostaneme zpětným překódováním. for I:=l to PU do begin J:=KC[I]; if J > 0 then D[I]:=Z[J] else D[I]:=G[-J] end; Algoritmus 1 nazýváme eliminační proto, že rovnice příslušné uzlům hranice i\ vůbec nesestavujeme a že parametry příslušné uzlům hranice i\ vyeliminujeme tím, že za ně 39 dosadíme předepsané hodnoty řešení a příslušné příspěvky z elementárních matic převedeme na pravé strany rovnic. Algoritmus 2 (pružinový) Při praktickém výpočtu diskrétního slabého řešení se někdy používá postup, který si nyní popíšeme. Je známo, že Dirichletova okrajová podmínka u = g může být přibližně splněna prostřednictvím okrajové podmínky Newtonova typu du kde u je velké kladné číslo. Lze ukázat, že pro o —> +00 platí u —> g na r\. Zahrňme tedy podmínku (2.59) do diskrétní slabé formulace (2.18). To se projeví tím, že bude Xh=Vh=Wh, že počet neznámých parametrů PN bude roven počtu uzlů PU a dále že do a,h(U,v) přibudou členy Ps(aUv) a do L h (v) členy Is(agv) pro strany S ležící na r\. Integraci nových hraničních členů proveďme lichoběžníkovou formulí (2.53). Sestavovací algoritmus 2 odvodíme z algoritmu 1. Nejdříve zpracujeme elementy e G T a strany S E S a sestavíme tak matici a vektor, které označíme K a F. Matici K a vektor F dále upravíme tak, že pro každý uzel Pi ležící na hranici Ti • přičteme adi k prvku ku matice K a • přičteme odig(x,i,yi) k prvku Fi vektoru F. Přitom di se rovná jedné polovině ze součtu délek těch (jedné nebo dvou) stran hranice Ti, které obsahují uzel P«. Protože členy adi mají být hodně velké, lze je nahradit členy Qi = ttku, kde k je velké číslo, například k = 1020. Protože Qí je výrazně větší než ku, a protože lze předpokládat, že Qí je současně výrazně větší než Fi, stačí pro každý uzel Pí ležící na hranici Ti • nahradit prvek ku matice K číslem nku a • nahradit prvek Fi vektoru F číslem Kkag(xi,yi). Následuje popis algoritmu 2. {Nulování globální matice K ~ K[1..PU, 1..PU] a vektoru F ~ F[1..PU]} for I:=l to PU do begin for J:=l to PU do K[I,J]=0; F[I]:=0 end; {Zpracování příspěvků z elementů e G T} for E: = l to PP do {prvky} begin {Výpočet elementární matice Ke ~ KE[1..3,1..3] a vektoru Fe ~ Fi^l.^] } for I:=l to 3 do {řádky} begin II:=PRVKY[E,I]; 40 for J:=l to 3 do {sloupce} begin JJ:=PRVKY[E,J]; K[II,JJ]:=K[II,JJ]+KE[I,J] {příspěvek do matice K} end; {sloupce} F[II]:=F[II]+FE[I] {příspěvek do vektoru F} end {řádky} end; {prvky} { Zpracování příspěvků ze stran S G S} for S:=l to PS do {strany} begin {Výpočet elementární matice Ks ~ KS[1..2,1..2] a vektoru Fs ~ FS[1..2]} for I:=l to 2 do {řádky} begin II:=STRANY[S,I]; for J:=l to 2 do {sloupce} begin JJ:=STRANY[S,J]; K[II,JJ]:=K[II,JJ]+KS[I,J] {příspěvek do matice K} end; {sloupce} F[II]:=F[II]+FS[I] {příspěvek do vektoru F} end {řádky} end; {strany} { Zpracování Dirichletovy okrajové podmínky} KAPA:=le20; for B:=l to PB do {body hranice rx} begin II:=BODY[I]; Q:=K[II,II]*KAPA; K[II,II]:=Q; {modifikace diagonálního prvku} F[II]:=G[I]*Q {modifikace pravé strany} end; {body hranice rx} Vyřešením soustavy rovnic s maticí soustavy K[1..PU, 1..PU] a vektorem pravé strany F[l..PU] dostaneme přímo pole D[1..PU] hodnot slabého řešení v uzlech. Algoritmus 2 nazýváme pružinový proto, že v úlohách pružnosti má a vystupující ve vztahu (2.59) význam tuhosti pružinového uložení. Příklad Uvažme úlohu, jejíž geometrie a triangulace je popsána pomocí obrázku 6 a tabulky 1 a o jejíž datech předpokládáme P = P, 9 = 0, f =Jy, 9 = g(l-x)2, a = a, /3 = ]3ny, kde p, f,g,ä a (3 jsou zadané konstanty a ny je y—vá složka jednotkového vektoru vnější normály hranice. Naším cílem je zápis několika rovnic soustavy (2.24). 41 Začněme tím, že uvedeme elementární matici prvku 1. Z (2.35) dostaneme p 2 / 1 -1 \ 0 -1 2 -1 0 \ -1 1 / element typu 1 Snadno ověříme, že elementární matice prvků 1, 3, 5, 6, 8, 10 a 12 (označme je jako elementy typu 1) jsou stejné: pro p(x,y) = p konstantní má element e a jeho obraz v zobrazení x = ax + b, ý = ay + c, (a, 6, c jsou konstanty) tutéž matici Kel, jak lze zjistit analýzou vztahu (2.35). Prvky 2, 4, 7, 9 a 11 označme jako elementy typu 2, jejich elementární matice je 2 V 0 -1 0 1 -1 -1 \ -1 2/ element typu 2 Prvky 13, 15, 17, 19, 20, 22 a 27 typu 3 mají elementární matici / 1 0 -1 \ 2 V 0 -1 element typu 3 a prvky 14, 16, 18, 21, 23, 26 a 28 typu 4 mají elementární matici / 1 -1 0 \ -1 1 2 -1 V 0 -1 2 -1 element typu 4 / Výpočtem ověříme, že elementární matice prvku 24 je stejná jako u prvků typu 2 a elementární matice prvku 25 je stejná jako u prvků typu 4. Elementární vektor Fe počítaný užitím formule (2.49) má na všech prvcích tvar kfpKe) ( ví \ \yij 0, 03125 pro prvky 26, 27, 28, kde pl(e) = l 0, 0625 pro prvky 24, 25, 0,125 pro ostatní prvky je plocha prvku e a y f je y-ová souřadnice jeho i-tého vrcholu. Elementární matice K.s a vektory Fs počítané podle (2.54) jsou tvaru K' 1 0 0 1 kde ď ■s ( 0,5 pro strany 1, 2, 11, 12, 0,5\/2 pro strany 3, 4, 5, 0, 25 pro stranu 6, 7, 8, 9, , 0,25\/2 pro stranu 10 42 je délka strany S, í 0,5 pro strany 3, 4, 5, Fs = ^ď^í J, kde dsn = < 0,25 pro stranu 7, 10, I 0 pro zbývající strany je součin dsny délky strany S a y-ové složky jednotkového vektoru vnější normály na této straně. Sestavme nejdříve rovnici příslušnou neznámé v uzlu 6. Z kódovací tabulky vyčteme KC[6]=4, takže půjde o čtvrtou rovnici. Sestavíme ji pomocí příspěvků z prvků 3, 4, 5, 9, 10 a 11 obsahujících uzel 6. Začněme například prvkem 3. To je prvek typu 1. Z pole PRVKY zjistíme, že uzel 6 je třetím uzlem prvku 3 (neboť PRVKY[3,3]=6) a proto použijeme třetí řádek elementární matice a elementárního vektoru. Protože KC[PRVKY[3,1]] = 1, KC[PRVKY[3, 2]] = 3, KC[PRVKY[3, 3]] = 4, příspěvek do levé strany rovnice bude ip[0-Ai-l-A3 + l-A4] a příspěvek do pravé strany rovnice i/0,125-1, neboť obsah prvku 6 je 0,125 a y-ová souřadnice uzlu 6 je 1. Stejně zpracujeme i zbývající prvky a výsledek zaznamenáme do tabulky. prvek typ řádek levá strana pravá strana 3 i 3 0,5p[-A3 + A4] 0,125/3/ 4 2 2 0,5p[A4- A2] 0,125/3/ 5 1 2 0,- 5p [-A2 + 2A4 - A5] 0,125/3/ 9 2 3 0,- 5p [-A3 - A7 + 2A4] 0,125/3/ 10 1 1 0,5p[A4- A7] 0,125/3/ 11 2 1 0,5p [A4 - A5] 0,125/3/ Sloučením jednotlivých příspěvků dostaneme rovnici p [-A2 - A3 + 4A4 - A5 - A7] = 0,25/. Stejný výsledek obdržíme standardní diskretizací rovnice p Au = / metodou sítí (stačí si uvědomit, že velikost oka sítě h = 0,5, takže h2 = 0, 25, a že v uzlu 6 je / = /). Další rovnicí, kterou si odvodíme, bude rovnice příslušná neznámé v uzlu 5. Protože KC[5]=3, jde o třetí rovnici. Sestavíme ji pomocí příspěvků z prvků 1, 2, 3, 7, 8 a 9. Narozdíl od odvození předchozí rovnice nyní uzel 4, společný vrchol prvků 1 a 7, leží na hranici Ti a je v něm tedy předepsána hodnota řešení u(P±) = g(P&) = g(l—0,5)2 = 0,257?. 43 Na hranici r\ leží také uzel 1, pro který u(Pi) jednotlivých prvků zaznamenáme do tabulky g(P1) = g(l - 1) 0. Příspěvky prvek typ řádek levá strana pravá strana 1 l 3 0, 5pA3 0,5pg(P4)- f o, 125 • 0, 5/3/ 2 2 2 0,5p [A3 - Ai] o, 125 • 0, 5/3/ 3 1 2 0, 5p[-A1 + 2A3- A4] o, 125 • 0, 5/3/ 7 2 3 0,5p [-A6 + 2A3] 0,5pg(P4)- f o, 125 • 0, 5/3/ 8 1 1 0,5p[A3- A6] o, 125 • 0, 5/3/ 9 2 1 0,5p[A3-A4] o, 125 • 0, 5/3/ a zjistíme, že člen g{P\) se mezi nimi neobjevil. Sestavená rovnice má tvar p [-Ax + 4A3 - A4 - A6] = 0, 25 [pg + 0, 5/]. Také tuto rovnici lze dostat standardní diskretizací používanou v metodě sítí. Nakonec sestavíme rovnici příslušnou neznámé v uzlu 16. Protože KC[16] = 12, jde o rovnici číslo 12. Sestavíme ji pomocí příspěvků z prvků 18, 19, 25, 26 a 27 užitím tabulky prvek typ řádek levá strana pravá strana 18 4 2 0,5p[-An+2A12-A8] 0,125 19 3 1 0,5p[Ai2-A8] 0,125 25 4 3 0,5p[-Ai6 + Ai2] 0,0625 26 4 3 0,5p[-A17 + a12] 0,03125 27 3 3 0,5p[-Ai7 - A13 + 2A12] 0,03125 1,5/3/ 1,5/3/ 1,5/3/ 1,5/3/ 1,5/3/ Protože strana 5 spojující uzly 12 a 17 je částí hranice r2, přičteme k levé straně člen A;f 2A12 = 0, 25v/2čřA12 a k pravé straně člen F£ = 0, 25/3. Na hranici T2 leží také strana 6 spojující uzly 16 a 17. Na této straně je /3 = 0, takže k pravé straně nepřičítáme nic, k levé straně však přičteme příspěvek kf xAi2 = 0,125ľrAi2. Celkem dostaneme rovnici -p A8 - 0, 5p An + [3, 5p + 0,125(1 + 2y/2) ä]A12 - 0, 5p A13-0, 5p A16 -pA17 = 0,125 [1, 5/ + 2^]. Získat analog této rovnice metodou sítí je nesnadné. 2.8. Několik poznámek Poznámka 13. V technických aplikacích nás často zajímá hodnota gradientu řešení. Na elementu e je gradient konstantní, VU = BeAe (2.60) 44 a odpovídá nejlépe hodnotě Vw v těžišti elementu e. Za aproximaci Vw v uzlu P se obvykle bere aritmetický průměr hodnot Vř7 ze všech elementů obsahujících uzel P jako svůj vrchol. □ Poznámka 14. Za předpokladu dostatečné hladkosti slabého řešení lze ukázat, že pro chybu u — U slabého řešení a jeho konečněprvkové po částech lineární aproximace platí max|w- U\ = 0{h2) a max ||V(u - U)\\ = 0{h) Ve G T, kde || • || označuje délku vektoru. □ Poznámka 15. Tam, kde se slabé řešení u prudce mění, je aproximace U pomocí po částech lineárních funkcí málo přesná, pokud triangulaci nezvolíme dostatečně přesně. V praxi často známe tato kritická místa. Pak stačí zahustit síť v kritické části oblasti Q, zatímco v nekritické části oblasti Q jemnou síť nevolíme, abychom zbytečně nezvětšovali počet neznámých. Takovýto způsob volby sítě je jednou ze značných výhod metody konečných prvků. □ Poznámka 16. Triangulaci volíme tak, aby žádný vnitřní úhel trojúhelníků triangulace nebyl příliš velký. Je-li to možné, používáme trojúhelníky s úhly rovnými nejvýše 7r/2 a pokud možno se vyhýbáme také úhlům příliš malým. Jsou-li vnitřní úhly všech trojúhelníků triangulace netupé (tj. menší nebo rovné 7r/2), pak uvážíme-li výklad za vztahem (2.42) a použijeme-li elementární matice Ke2 tvaru (2.44) a elementární matice K.s tvaru (2.54), zjistíme, že matice tuhosti K má na hlavní diagonále kladné prvky a mimo hlavní diagonálu prvky nekladné. Navíc lze ukázat, že K je íreducíbílně diagonálně dominantní matice (stručně IDDM), viz [7], což znamená, že: a) K je ireducibilní, tj. pro libovolné dva indexy i0,j, 1 < i0,j < PN, existují nenulové koeficienty kioil,kilÍ2,... ,kie_ljie, kde ie=j; j3) K je diagonálně dominantní, tj. pn IM>X^|%| pro i = l,..., PA^; 7) existuje index i0, pro který pn l^*o*ol ^ ^ ] l^*oiľ IDDM s kladnými diagonálními a nekladnými mimodiagonálními prvky je tzv. regulární M-matice, viz [13], [8]. Je známo, že symetrická regulární M-matice je pozitivně definitní. Je-li tedy matice tuhosti K symetrická regulární M-matice, je pozitivně definitní nezávisle na velikosti h. □ Poznámka 17. Triangulace jsou obvykle generovány automaticky. Lokální zhušťování se často provádí bisekcí: spojením středů stran dostaneme z originálního trojúhelníka typu A čtyři shodné menší trojúhelníky typu B. Napojení na okolní hrubou triangulaci 45 zajistí dva přechodové trojúhelníky typu C vzniklé rozdělením trojúhelníka typu A těžnicí. Lokální zjemnění triangulace je ilustrováno na obrázku 7. □ Poznámka 18. Doposud jsme předpokládali, že Q je polygon. To nám umožnilo bez problémů vykrýt oblast Q triangulací T. Jak ale postupovat, když Q je regulární oblast s hranicí po částech hladkou a hladké hraniční oblouky jsou křivé? I v takovém případě pokryjeme oblast Q triangulací T. Přitom T rozumíme stejně jako dříve množinu trojúhelníkových elementů e« takových, že ěj děj pro i ^ j; je buďto množina prázdná nebo je to společná strana nebo společný vrchol. Označíme ňh = \Jě Hranici oblasti Q h značme dílh nebo také IV Naším cílem samozřejmě je, aby náhrada oblasti Q pomocí byla pokud možno co nej dokonalejší. Toho lze docílit, přijmeme-li následující předpoklady: 1) uzly triangulace a těžiště jejich elementů leží v fž; 2) uzly ležící na hranici leží také na hranici T; 3) společné body průniku Tif^ jsou uzly triangulace. Z podmínky 2 mimo jiné vyplývá, že při vykrývaní oblasti Q elementy e nesmíme Obr. 7. Zjemnění triangulace bisekcí uvnitř náhradní oblasti vytvořit omylem nějaké „díry". Je zřejmé, že nemusí platit Qh C Q ani Q C ÍV Avšak obě oblasti Q a se mohou lišit jen v okolí svých hranic T aľ/,a kdyby se křivé oblouky hranice T vypřímily, pak by se oblasti Q h a Q ztotožnily. Podmínka 3 nám zase umožňuje přirozeným způsobem přiřadit k části hranice Ti její aproximaci IV a k části T2 její aproximaci T2h tak, aby = IV U T2h a přitom TVnr2ft =T1nT2. Algoritmus popisující vytvoření globální matice K a globálního vektoru F se změní jen nepatrně. Malý problém totiž představuje ta skutečnost, že Q h nemusí být částí Q a že ľ / Td, takže funkce p, q a / nemusejí být definovány na celé oblasti ílh, funkce g nemusí být definována na a funkce a a /3 nemusejí být definovány na r2^. S funkcí g žádné problémy nevzniknou, neboť hodnoty této funkce potřebujeme znát jen v uzlech hranice Tíh a tyto uzly jsou podle předpokladu 2 současně uzly hranice IV U zbývajících funkcí stačí znát jejich hodnoty v uzlech formulí numerické integrace. Vzhledem k podmínce 1 a 2 žádné problémy nenastanou, jsou-li uzly integrace současně uzly triangulace a nebo jsou-li těžišti elementů e. Takže jedinou potíž představují vztahy (2.56) a (2.57) obsahující hodnoty funkcí a a /3 ve středu (rrf,, yfl) strany S. Ale zde je pomoc snadná: hodnotu ve středu nahradíme aritmetickým průměrem hodnot v koncových bodech strany. □ 46 2.9. Minimalizační formulace Definujme funkcionál U(w) = \a(w,w) - L(w) (2.61) a zabývejme se minimalizační úlohou najít u G W splňující U(u) < U(w) Vu> G W. (2.62) Ukážeme, že „slabá" úloha (2.12) a minimalizační úloha (2.62) mají stejná řešení. Předpokládejme nejdříve, že u je řešením (2.12). Nechť w G W a položme v = w — u, takže v EV a, w = u + v. Máme Yi{w) = U(u + v) = \ a{u + v, u + v) — L{u + v) = = | a(u, u) + | a(it, i?) + | a(-u, it) + | a(-u, v) — L(u) — L(v) = = | a{u, u) — L{u) + a{u, v) — L (v) + \ a(v, v) = = Ľ.(u) + ±a(v,v)>Ľ.(u), poněvadž a(u, v) = a(v, u), a(u, v) — L(v) = 0 a a(v, v) > 0, takže u je řešením (2.62). Obráceně, nechť u je řešením (2.62). Potom pro každé v E V a každé reálné číslo t je n(it) < Il(u + tv). Tedy diferencovatelná funkce tp{t) = U(u + tv) = \ a{u + tv,u + tv) — L{u + tv) = = | a{u, u) + | a(it, to) + | a(to, it) + | a(to, to) — L(it) — L(tv) = = | a{u, u) — L{u) + t[a{u, v) — L (v)] + | t2a{v, v) má minimum pro ŕ = 0 a tudíž 0 už matice soustavy rovnic (2.69) diagonální není a proto vždy užijeme bezpodmínečně stabilní ^-metodu. Nejpřesnější výsledky obdržíme pro Crankovo-Nicolsonovo schéma, kterému odpovídá 9 = \ (metoda je 2. řádu přesnosti), nej stabilnější je však implicitní Eulerovo schéma, kterému přísluší 9 = 1. Za předpokladu existence dostatečně hladkého slabého řešení lze ukázat, že platí , f 0(h2 + At) pro 9 Ý h , ^T)\--U\ = {0{h2 + Aŕ) pro, = i (2.71) kde At = maxj AU. Proto se nejčastěji používá Crankovo-Nicolsonovo schéma. Jen v případě nesouladu počátečních a okrajových podmínek nebo při prudké změně okrajových podmínek je účelné pro několik časových kroků použít implicitní Eulerovo schéma, které je sice méně přesné, zato však stabilnější. 50 2.11. Dynamika Budeme se zabývat úlohou určit funkci u{x,y,t), (x,y) G íl, t G (0,T), vyhovující diferenciální rovnici ^-V-bVíí] + g« = / pro (x,y)eQ,te(0,T) (2.72) a splňující okrajové podmínky u = g pro (x,y) G ľi, í G (0,T), (2.73) -p— = cm - /3 pro (x, y) G r2, í G (0,T) (2.74) on a počáteční podmínky du(x,y,0) — u(x,y,0) = ip{x,y), -—-= ijj(x,y) pro (x,y) G íl. (2.75) Úloha (2.72) - (2.75) je dobrým modelem kmitající membrány. Přitom ŕ je čas, u{x,y,t) výchylka, (du/dt rychlost a d2u/dt2 zrychlení), g(x, y, t) hustota, p{x,y,t) charakterizuje tuhost membrány, q{x,y,t) odpor okolního prostředí, g{x,y,t) předepsaný pokles podepřené části r\ hranice, a(x, y, t) tuhost pružinového uložení a j3(x, y, t) zatížení na části T2 hranice. Při odvozování slabé formulace postupujeme podobně jako v úloze vedení tepla a dostaneme (gutt, v) + a{u, v) = L (v) Wv G V a pro ŕG(0,T), (2.76a) u = g pro (x, y) G Ti, t G (0,T), (2.76b) du(x,y,0) — u(x,y, 0) = oo a jim odpovídající nenulové vlastní funkce ui, u2,... s vlastnostmi a(ui,Uj) = (gUi,Uj) = 0 pro i ^ j. Příslušná diskrétní formulace je tvaru najít AGtaí/Gl4, Í//0, splňující ah(U,v) = A(gU,v)h \/v G Vh. (2.89) Lze ukázat, že existuje P N kladných vlastních čísel 0 < Ax < A2 < • • • < Ap^ a jim odpovídajících nenulových vlastních funkcí Ui, U2, ■ ■ ■ Up^ takových, že pro i ^ j platí ah(Ui,Uj) = (gUi,Uj)h = 0. Maticový zápis (2.89) je najít A G M a ó G RPN, ô ^0, splňující Kô = AMô. (2.90) Vlastní vektory ôi příslušné vlastní číslům Aj lze určit tak, že ôjK.ôj = ôjMôj = 0 pro i / j a á. Máj = 1. Přitom souvislost mezi vlastním vektorem ôi = (5j,i, • • •, 5í,pn)t 53 a vlastní funkcí U i je dána vztahem Ui(x,y) = Y2j=iSijWj(x,y). Zobecněný problém vlastních čísel (2.90) řešíme pomocí vhodné numerické metody. V praxi nás přitom obvykle zajímá jen několik nej menších vlastních čísel a jim odpovídajících vlastních funkcí. Mezi nejčastěji používané metody patří metoda iteraci v podprostorech, viz [2],[4], Lanczosova metoda, viz [4] nebo Arnoldiho metoda, viz [12],[24]. Za předpokladu existence dostatečně hladkého slabého řešení platí maximi - Ui\ = 0(h2), \\ - Ať| = 0(h2) pro i = 1,PN. n Ukažme si, jak můžeme pomocí vlastních čísel A j a vlastních vektorů ô í určit řešení úlohy (2.77) v případě, že matice M a K nazávisejí na t. Řešení hledejme ve tvaru A(ŕ) = Yq(ŕ), kde Y = (<5i,..., Ôpn) a q(r) = (qi(t),..., qpN(t))T je třeba určit. Dosazením do MA + KA = F dostaneme MY q + KYq = F. Vynásobíme-li tuto rovnici zleva maticí YT, dostaneme q + Aq = YTF, neboť YTMY = I je jednotková matice a YTKY = A, kde A = diag{A1;..., Ap^} je diagonální matice s vlastními čísly Ai,..., Apn na hlavní diagonále. Dostaneme tedy PN na sobě nezávislých obyčejných diferenciálních rovnic q\(t) + AjOj(r) = hi(t), kde hi(t) je i-tá složka vektoru YTF(r). Vynásobením rovnice A(0) = Yq(0) = tp resp. A(0) = Yq(0) = V zleva maticí YTM dostaneme q(0) = YTM<^ resp. q(0) = YTM?/>, takže funkce g«(r) je určena počátečními podmínkami Oj(0) = rÍ7 (ji(0) = Sj, kde resp. Sj je i-tá složka vektoru YTM<^ resp. YTMV'. V případě úlohy (2.77') a proporcionálního útlumu, tedy pro C = %M + a^K, opět klademe A (ŕ) = Yq(ŕ) a i-tou složku g«(r) vektoru q(r) získáme snadno řešením počáteční úlohy q\(t) + (aM + aKAi) qt(t) + Aiqi(t) = h,t(t), qi(0) = r,t, qi(0) = s,t, (2.91) kde hi(t), Ti a Sj jsou stejně jako dříve i-té složky vektorů YTF(ŕ), YTM<^ a YTM?/>. V některých případech lze získat uspokojující aproximaci řešení pomocí jen několika nejmenších vlastních čísel a vektorů. Položíme Y = ... ,6m) pro vhodné m < PN, řešením úloh (2.91) pro i = 1,... ,m určíme g«(ŕ) a dostaneme A (ŕ) = Yq(ŕ), kde q(ŕ) = (qi(t),...,qm(t))T. 2.12. Rovinná napjatost a rovinná deformace Doposud jsme se zabývali skalárními úlohami, jejichž řešením byla jediná funkce. V tomto odstavci si uvedeme významnou aplikaci, jejíž řešení tvoří dvojice funkcí. a) Klasická formulace V rovinné úloze pružnosti nebo deformace vektor napětí cr = (o-x,o-y,TXy)T (2.92) splňuje Cauchyovy rovnice statické rovnováhy ^ + ^ + /i = 0, ^ + ^ + /2 = 0 vro(x,y)en, (2.93) ox oy ox oy 54 kde ax(x,y) resp. ay(x,y) je normálové napětí v řezu x = konst. resp. y = konst., rxy je tečné napětí v řezu x = konst. nebo y = konst. a f\(x,y) resp. f2(x,y) je složka vektoru f=(/i,/2)r (2-94) objemového zatížení působícího ve směru osy x resp. y. Hookův zákon (fyzikální rovnice) zapíšeme ve tvaru 99) I qAT(1 + z/) pro rovinnou deformaci, Přitom a je koeficient tepelné roztažnosti, AT je změna teploty a v je Poissonův součinitel příčného zúžení (0 < z/ < 0,5). Pro izotropní materiál jsou z prvků matice tuhosti D nenulové pouze členy dn=d22, d12=d21 a d33 a platí pro ně v případě rovinné napjatosti _ _ E _ _ Ev _ E d\\ = d22 =---, d\2 = d2\ =---, (Í33 = ———■—-, (2.100) 1 — uz 1 — uz 2(1 + v) a v případě rovinné deformace , _ , _ E(l - v) Ev E du — «22 — 7——rf.-^r-T, «12 — «21 — 7——rf.-^r-T, «33 — TTf-i I—v (2.101) (1 + z/)(l - 2u) (1 + z/)(l - Iv) 2(1 + z/) kde i? je Youngův modul pružnosti v tahu resp. tlaku. Z geometrických rovnic dui du2 dui du2 £x=c^^ £y=<3y-> 7- = ^+9^' (2-102) dostaneme vztah mezi vektorem deformace a vektorem posunutí u=(Ul,u2)T, (2.103) 55 v němž ui(x,y) resp. U2(x,y) je posun bodu (x,y) ve směru osy x resp. y. Dosazením ze vztahů (2.95), (2.102) do (2.93) dostaneme soustavu dvou parciálních diferenciálních rovnic d_ dx dy a dui , a "11 --r "12 d 31 ■ dx dui dx d- 32 du2 f dis ( dux du2 dy v dy dx du2 m33 ( dui ^_ du2 dy dx d d h ~ 7^ [duel + cři2£° + cři37°J - [^3i4 + d32S d^lly] d ' dui «31 n " dx , a 9u2 f «32 n " dy Vd33 / <9mi du2 dx V dx d «21 n - a 9u2 Vd22 n -dy M23 ^ (9-ui ^_ du2 dy v dy dx (2.104) Í2 d d [d31e°x + d32e°y + g^sT^J - ^ [feS + «'22^° + ^237^] • dx Pro izotropní materiál se rovnice (2.104) nazývají Laméovy statické rovnice. Zavedeme-li Laméovy konstanty Ev A v Ev pro rovinnou napjatost, pro rovinnou deformaci, E (1 + v){l - 2v) pak dosazením do (2.100) a (2.101) snadno ověříme, že dn = d22 = ^ + 2/x, di2 = A, d33 = //, d13 = d23 = 0, 2(1 + v) (2.105) (2.106) a vyjádříme-li počáteční poměrné prodloužení e° pomocí (2.99), dostaneme Laméovy rovnice ve tvaru -//Au — (A + //) grad div u = f — 2(A + //) grad et kde Au Aw2 A-u,; <92w,- d2m dx2 dy2 div u dui du2 —1 H---■ dx dy Okrajové podmínky uvažujeme dvojího druhu. Na části Ti hranice předepišme geometrické okrajové podmínky ui = gi, u2 = g2 na Ti a na části r2 hranice statické okrajové podmínky axnx + rxyny = pi, rxynx + ayny = p2 na T2. (2.107) (2.10Í 56 Zde g±(x,y) a g2(x,y) jsou předepsané posuny na části I\ hranice a p±(x, y) a p2(x,y) jsou rr-ová a y-ová složka povrchového zatížení. Označme g = (g±,g2)T a p = (pi,P2)T- Budeme předpokládat, že geometrické okrajové podmínky jsou předepsány alespoň na části hranice, tedy že r\ ^ 0. Rovnice (2.104) spolu s okrajovými podmínkami (2.107) a (2.108) představují klasickou (diferenciální) formulaci úlohy rovinné napjatosti či deformace. b) Slabá formulace Definujme prostor X dvojic funkcí v(x,y) = (vi(x,y),v2(x,y))T předpisem X = {v K G H1 (Q), v2 G H1 (Q)}, (2.109) prostor testovacích funkcí V = {vGX|v = 0 na Ti} (2.110) a množinu přípustných řešení W = {vGX|v = g naľi}. (2.111) Vynásobme první z Cauchyových rovnic (2.93) funkcí v\ a druhou funkcí v2, kde v = (vi,v2) G V, obě rovnice sečtěme a integrujeme přes Q. Užitím Greenovy formule (2.6), statických okrajových podmínek (2.108), Hookova zákona (2.95), geometrických rovnic (2.102) a označení e(v)=(^,??,^ + P-)T (2.112) dostaneme dor dr, xy dx dy h Vl drxy doy dx dy Í2 v2 > dxdy dn n L {[axnx + txyny] vx + [rxynx + ayny] v2} ds- dvi dvi dv2 ' dx dvo dxdy + / [/i^i + f2v2] dxdy dx dy dx dy = / [Pivi + P2V2] ds — e (v) • cr dxdy + / v • f dxdy = = v • p ds — / e(v) • D [e(u) — e°] dxdy + / v • f dxdy. Jt2 Jn Proto slabou (variační) formulaci úlohy rovinné napjatosti či deformace lze psát ve tvaru najít u G W splňující a(u, v) = L(y) Vv G V, (2.113) kde a(u,v)= / e(v) • Be(u)dxdy, (2.114) Jn L(v) = / v-ídxdy+ / e(v) • Ds° dxdy + / v-pds. (2.115) Jn Jn Jt2 57 Jestliže funkce G PC(í!) pro i, j = 1,2,3, funkce /1? /2, 7^ G PC(fž), funkce Pi, p2 G PC(r2) a funkce gi, g2 G C(ľi), pak existuje jediné slabé řešení úlohy (2.113). c) Diskrétní slabá formulace Oblast Q triangulujeme stejně jako v odstavci 2.5 a definujeme konečněprvkový prostor funkcí Xh = {v\v1eXh,v2eXll}. (2.116) Dimenze prostoru je rovna 2 PU. V prostoru X^ volíme bázové funkce wn(x,y) = (wj(x, ř/),0)r, víj2{x,y) = (0,Wj{x,y))T, j = l,...,PU. (2.117) Pak každou funkci v(x,y) G XZi můžeme zapsat ve tvaru PU v(x,y) = ^2 \vl(xj,yj)w:ji(x,y) + v2(xj,yj)wj2(x,y)} . (2.118) Konečněprvkový prostor testovacích funkcí definujeme předpisem Yh = {v G X, I v(P;) = 0 VP, G r\} (2.119) a množinu konečněprvkových přípustných funkcí určíme jako množinu wh = {v G X, I v(P,) = g(P,) VP, G ľx}. (2.120) Pro v G Vft platí PN v(x,y) = [vÁxj>yj)™ji(x>y) + v2(xj,Vj)wj2(x,y)] (2.121) a pro v G je PN v(^ ž/) = -'A^'Wjm-''- V) + '2Í-0-//;:iw;-'(-ľ-ž/)] + 3=1 (2.122) pi/ + j] [9i(xj,yj)^ji(x,y) + 92(xj,yj)wj2(x,y)}. j=PN+l Pak diskrétní slabou formulaci úlohy rovinné napjatosti či deformace lze zapsat ve tvaru najít U G W,, splňující ah(U,\) = Lh(v) Vv G Vft, (2.123) kde (U,v) = X)/B(e(v)-De(U)), (2.124) °& eeJ Ses 58 Hodnotu složky Ui řešení U v uzlu Pj označíme Aji = Ui(xj,Vj) a hodnotu složky vi testovací funkce v v uzlu Pi označíme Qn = vi(xi,yi). Podle (2.119) a (2.120) je PN ufo y) = [ajiwji(x' y) + Ai2wi2(x, y)] + PU + Yl Í9i(xj,yj)wji(x,y) +92(xj,yj)wj2(x,y)}, (2.126) j=PN+l PN v(x, y)=Y [®íiwn(xi y) + ©í2wí2(x, y)]. Vzhledem k bilineárnosti formy (^(U, v) a lineárnosti funkcionálu L^(v) lze dosazením z (2.126) do (2.124) a (2.125) odvodit PN (PN \ 0 = ah(u, v) - Lh(v) = ^[e,]T J ^KíjAj ~fA= 0t[KA - F], (2.127) i=i l i=i J kde jsme užili označení v - ív \PN k - f ah(™ji,™íi) a/l(wj-2,wii) V ^(WjijW^) aft(wi2,wí2) F = (Fľ,...,F?„f a F,= ( ^f]e,..., [Vw; }e) = [Je]- Le = [det Je}- Ge~Le kde (vivr,..., vivpee j dŇl 9Ň;c \ dŇ: Pe (2.169) (2.170) Označme ještě Pak z (2.33), (2.157), (2.169) a (2.48) plyne Kel = f([BfpBe) = Je([Be]TpeBe|detJe|) = = Ie([Ĺe]T[Ge]Tpe GeLe|det J6)"1), Ke2 = /e([Ne]rř/Ne) = Je([Ňe]TgeŇe|det Je|), (2.171) (2.172) (2.173) 72 Fe = Je([Ne]J /) = Je([Ne]J f |det Je|). (2.174) Závěrem uveďme, jaké kvadraturní formule Je(-) je vhodné použít pro výpočet elementárních matic Kel, Ke2 a elementárního vektoru Fe. Pro prvek T3 jsme takové formule již uvedli v části a) odstavce 2.7. Zopakujme tedy, že při integraci prvků matice Kel užijeme formuli ie(g) = \g(\,\) (2.175) a při integraci prvků matice Ke2 a vektoru Fe formuli Ie{g) = i[£(0, 0) + g(l, 0) + g(0,1)]. (2.176) Obě formule (2.175) a (2.176) integrují přesně polynomy ^rf pro 0 < i + j < 1. Pro zbývající typy prvků doporučíme vždy jen jednu kvadraturní formuli pro výpočet obou matic Kel, Ke2 i vektoru Fe. Pro prvek T6 je vhodné použít formuli Ie(g) = \\g(\, 0) + g{\, \) + g(0, ±)], (2.177) která integruje přesně polynomy ^rf pro 0 < i+j < 2. Pro prvek Q4 použijeme Gaussovu součinovou formuli 2x2 Ie(g) = g(—a, — a) + g(a, — a) + g(a, a) + g(—a, a) pro a = (2.178) která integruje přesně polynomy ^rf pro 0 < i, j < 3. Konečně pro prvky Q8 a Q9 je vhodné použít Gaussovu součinovou formuli 3x3 ÍeÍ9) = |^(0, 0) + f [g(-a, 0) + g(0, -a) + g(a, 0) + g(0, a)] + (2.179) 25 ' + j^[g(—a,—a) + g(a,—a) + g(a, a) + g(—a, a)] pro a=y0, 6, která integruje přesně polynomy ^rf pro 0 < i, j < 5. Přejděme nyní k odvození elementární matice k.s a vektoru Fs na straně S. K tomu účelu si označíme N5 = Ns(x,y) = (wS(x,y),...,w%s(x,y)), As = (A?,..., A^)T, kde Af = U(xf, yf), i = l,...,ps, @s = (0f,..., e^)T, kde ef = í;(rrf,yf), i = 1,... ,ps. Pak slabé řešení U a testovací funkci v lze na straně S* vyjádřit ve tvaru U = N5 A5, v = Ns@s. Označme ještě Ňs = (iVr,...,JV£). 73 Pak ze vztahu (2.51), v němž místo a píšeme a1, vztahu (2.52), v němž místo /3 píšeme /37 a ze vztahu (2.158) plyne Ks = /SqN5]tc/N5) = J5([Ň5]Ta/5Ň5J5), (2.180) Fs = /S^nY^) = IS([NS]TPISJS). (2.181) Zbývá dodat, jaké kvadraturní formule Is(-) použijeme pro výpočet elementární matice k.s a elementárního vektoru fs. Pro strany E2 (příslušné prvkům T3 a Q4) užijeme lichoběžníkovou formuli ÍS(g) = lM0) + m) (2-182) integrující přesně polynomy prvního stupně a pro strany E3 (příslušné prvkům T6, Q8 a Q9) užijeme buďto Simpsonovu formuli J5(^) = ^(0)+4^(i) + ^(l)] (2.183) nebo Gaussovu formuli se dvěma uzly ÍS(9) = m3-r)+9(3-±fi)}. (2.184) Obě formule (2.183) i (2.184) integrují přesně polynomy stupně třetího. Globální matici soustavy a globální vektor zatížení sestavíme modifikací algoritmů uvedených v části c) odstavce 2.7. Po vyřešení soustavy rovnic získáme hodnoty U(xi,yi), i = 1,...,PU. Na každém prvku proto známe Af=U(xf,yf), i = l,...,pe, a můžeme spočítat hodnotu konečněprvkového řešení U a jeho gradientu \7U v libovolném bodu P*(x*,y*) prvku e. Protože bázové funkce w^(x,y) obecně neumíme explicitně vyjádřit, musíme přejít na referenční element é. Najdeme vzor P*(^*,i]*) bodu P* v zobrazení (2.153), tedy rozřešíme obecně nelineární soustavu rovnic ** = *7*), y* = ye(C,v*)) (2-185) pro neznámé £*, rf, a pak určíme f/(^,^) = f>(r,^) = Ňe(r,^)Ae, VU(x*,y*) = WB(£W) = Be(£W)Ae = {[det Je]_1GeLe}(^*,?]*)Ae. Pokud se spokojíme s hodnotou gradientu \7U v některém z uzlových bodů Pf(x^,yf), pak pro určení jeho vzoru žádnou soustavu rovnic řešit nemusíme, neboť souřadnice uzlu PíÍ^íiVí) známe. Je-li soustava (2.185) lineární, rozřešíme ji snadno, a je-li nelineární, pak k jejímu vyřešení obvykle postačí provést jen několik málo kroků Newtonovy metody. Za předpokladu dostatečné hladkosti slabého řešení lze ukázat, že pro chybu u — U slabého řešení a jeho konečněprvkové aproximace platí max>- U\ = 0(hr) a max ||V(u - U)\\ = 0{hr-v) Ve G T, (2.186) Í2hní2 ěníí 74 kde r=2 pro prvky T3, Q4 a r=3 pro prvky T6, Q8, Q9. Je-li q = 0 a jestliže používáme čtyřúhelníkové prvky, které se málo liší od rovnoběžníků, pak ve speciálních bodech P těchto prvků dostáváme přesnější hodnoty Vř7: [V(it — U)](P ) = 0(hr). Říkáme, že v těchto bodech nastává superkonvergence gradientů. Pro prvek Q4 super konvergence nastává v bodu, který je obrazem uzlu Gaussovy kvadraturní formule 1 x 1 na referenčním elementu é, tedy v bodu P (xe(0, 0), ye(0, 0)), a pro prvky Q8 a Q9 nastává superkonvergence v bodech, které jsou obrazy uzlů Gaussovy kvadraturní formule 2x2 na referenčním elementu é, tedy v bodech Pe(xe(±a, ±a),ye(±a, ±a)), kde a = Zajímá-li nás gradient ve vrcholu P triangulace T, pak ho spočteme jako aritmetický průměr z hodnot gradientů v „přiléhajících" Gaussových bodech prvků obsahujících tento vrchol P. Superkonvergence gradientů v Gaussových bodech čtyřúhelníkových elementů se běžně využívá v úloze rovinné napjanosti či deformace. (V bilineární formě a(u, v) rovinné úlohy napjatosti či deformace, viz (2.114), se nevyskytuje analog členu quv obsaženého v bilineární formě a(u,v) úlohy (2.1)-(2.4), viz (2.9). Proto odpadá analogie výše uvedeného požadavku q = 0 a k dosažení superkonvergence v Gaussových bodech stačí, abychom používali čtyřúhelníkové prvky tvarově blízké rovnoběžníkům.) Poznámka 22. Prvek Q9 je přibližně stejně přesný jako prvek Q8, ve srovnání s ním však má parametr Ag = Ue(xg,yg) navíc. To ale není zase tak velký nedostatek, neboť rovnici pro parametr 0g můžeme získat hned při zpracovávání elementu e, má totiž tvar Yýj=i k9jAj + k99A9 = Fš'> odtud lze parametr Ag vyjádřit, Aeg = [F9e - £)®=1 k^A^/k^, a dosadit ho do zbývajících rovnic. Také toto dosazení lze provést „lokálně": zvolíme ®9 = [— Si=i kí9®í]/k99 a d° rovnice 0 = dh(U, v) — Lh(v) vložíme za element e příspěvek [0e]t[keAe - fe] = [0e]t[keAe - fe], kde ée = (e?,...,e§)T, äe = (a?,...,adt, k = {%}fJ=1, fe = (Ff,...,F8ef (šikovně zvolený „vázaný" parametr ©g zaručuje symetrii matice ke). Redukovanou matici ke a redukovaný vektor fe lze získat snadno užitím Gaussovy eliminační metody: v soustavě rovnic keAe = fe odeliminujeme prvky A;|g a i = 1,..., 8, posledního sloupce a řádku a v prvních osmi řádcích a sloupcích upravené matice soustavy najdeme matici ke a v prvních osmi řádcích upraveného vektoru pravé strany najdeme vektor fe. Elementární matici ke a vektor fe rozešleme do globální matice k a globálního vektoru f obvyklým způsobem. Po vyřešení globální soustavy rovnic na elementu e dopočítáme Ag = [Fg — Y^8j=i k9j^j]/k99- Právě popsaná technika eliminace „vnitřních" parametrů elementu (zde jednoho) bývá označována jako metoda kondenzace vnitřních parametrů. Zobecnění této metody, spočívající v „předeliminaci" parametrů vnitřních vzhledem k celé skupině elementů, bývá označováno jako statická kondenzace. Princip je následující: konstrukci rozčleníme na subkonstrukce, vnitřní parametry subkonstrukcí odeliminujeme (částečná eliminace ne příliš rozsáhlých soustav s řídkými maticemi), ze zbytků matic a vektorů subkonstrukcí příslušných „hraničním" parametrům subkonstrukcí sestavíme globální soustavu rovnic (bohužel již s plnou maticí), vyřešíme ji a získáme hraniční parametry a nakonec na subkonstrukcích dopočítáme jejich vnitřní parametry. Předpokladem dosažení očekávaného efektu, tj. snížení celkového objemu výpočtů ve srovnání s případem, kdy řešíme konstrukci jako jeden celek, je kromě promyšlené strukturalizace konstrukce zejména kvalitní implementace zmíněné techniky statické kondenzace. □ 75 2.14. Nelineární úlohy Mnohý problém technické praxe vystihneme podstatně realističtěji, když k jeho popisu použijeme nelineární parciální diferenciální rovnici a případně také nelineární okrajové podmínky. Tak například úloha (2.1) - (2.4) je dobrým modelem stacionární úlohy vedení tepla při relativně nízkých teplotách. Při vyšších teplotách je však už třeba předpokládat, že funkce p, q, f, a a /3 závisejí také na teplotě u. V nestacionární úloze vedení tepla (2.63) - (2.66) je třeba při vyšších teplotách uvážit také závislost funkce c na teplotě u. Poznamenejme, že některé praktické problémy nelze vůbec jako lineární formulovat, a to ani za zjednodušujících, avšak stále ještě realistických předpokladů. V dalším budeme předpokládat, že funkce p, q i / mohou záviset také na prvních parciálních derivacích ux a uy hledaného řešení u. Diskretizace metodou konečných prvků probíhá u nelineárních úloh formálně stejně jako u úloh lineárních: ve slabé formulaci nahradíme slabé řešení u MKP-slabou aproximací U. Závislost „dat" p, q, f, a, j3, c na u (nebo i na ux, uy) způsobí, že soustava algebraických rovnic (2.24) případně soustava obyčejných diferenciálních rovnic (2.68) nebo (2.77') je nelineární (vystupují v ní hodnoty funkcí p, q, f, a, j3, c závislé na U (nebo i na Ux, Uy) a tedy závislé na A popřípadě A(t)). Postupy, kterými se vypořádáváme s nelinearitou v úlohách stacionárních se obvykle liší od postupů používaných ke zvládnutí nelinearity v nestacionárních úlohách. a) Stacionární úloha Stejně jako v lineární úloze sestavíme elementární matice Ke a elementární vektory Fe pro e G T a elementární matice K.s a elementární vektory Fs pro S G S. Nyní však Ke a Fe budou záviset na Ae a Ks a Fs na As. Po sestavení globální matice soustavy K(A) a globálního vektoru pravé strany F(A) dostaneme nelineární soustavu rovnic R(A) = K(A)A — F(A) = 0 (2.187) pro neznámý vektor parametrů A. Soustavu (2.187) řešíme některou z metod pro řešení nelineárních soustav rovnic. Vyjdeme z počáteční aproximace A*-0-1 a zvolenou metodou počítáme další aproximace A*-^, k = 1,... tak dlouho, až pro k=K lze považovat za dostatečně přesnou aproximaci řešení A. K určíme tak, aby buďto reziduum R(A^) bylo dostatečně malé, nebo aby byl dostatečně malý rozdíl dvou po sobě jdoucích aproximací A(-x*)-A(-x~1*), případně aby byly dostatečně malé oba dva tyto vektory. Nejjednodušší metodou řešení soustavy (2.187) je metoda prosté iterace: K(A(fc))A(fc+1) = F(A(fc)). (2.188) Metoda prosté iterace je velmi jednoduchá, konverguje však spíše výjimečně, a proto její použití je omezené. Nejčastěji se soustava (2.187) řeší Newtonovou metodou: H(A(fc))Y^ = -R(A(fc)), (2.189) A(k+i) = A(k) + YW) (2.190) kde H(A) = {dRi(A)/dAj}™=1 (2.191) je Jacobiova matice zobrazení R(A) = (/^(A),..., Rp^(A))T. Globální matici H a globální vektor R sestavujeme z elementárních matic He a elementárních vektorů Re pro 76 e G T a elementárních matic H5 a elementárních vektorů R5 pro S E S. Předpokládejme, že pracujeme s izoparametriekými prvky zavedenými v odstavci 2.13. Označme Re(Ae) = Ke(Ae)Ae - Fe(Ae) (2.192) lokální reziduum na prvku e a He(Ae) = {0i2?(AB)/3A!}£.=1 (2.193) Jacobiovu matici zobrazení Re(Ae) (i?f (Ae) je i-tý člen vektoru Re(Ae)). Užitím (2.33) a (2.172)-(2.174) dostaneme Re(Ae) = [Je(peKpe) + Ie(qeKqe)]Ae - Je(/eF/e), (2.194) kde Kpe = [Ĺe]T[Ge]TGeLe|det Je|_1, (2.195) Kge = [Ňe]TŇe|det Je|, (2.196) F/e = [Ňe]T|det Je|. (2.197) Pro Jacobiovu matici He zobrazení Re(Ae) platí He = We + We - H/e, (2.198) kde RPe = Kel + re(ŔPeAe[^l) (2199) 1 Pe Hs.=1 (2.206) (2.207) Jacobiovu matici zobrazení RS(A'S) (Rf(As) je i-tý člen vektoru RS(A'S)). Užitím (2.180) a (2.181) dostaneme R5(A5) = Ís(áISKaS)As - Is0ISF^s) kde KaS = [ŇS]TŇsjs, F13'3 = [NS]Tjs. Pro Jacobiovu matici Hs zobrazení RS(A'S) platí H5 = HaS - H13'3, kde H dá is dáIS <9Af " '' d AI "Ps d/3IS d/3 is <9Af '"'9A|s Přitom z (2.159) a (2.156) plyne ])• (2.208) (2.209) (2.210) (2.211) (2.212) (2.213) PS áIS(^,As) = Y,<*tä,vf, Af)Äf (0 i=l Ps Pis(t As) = Y,Ptä,Vi, *i)Nf (0- (2.214) (2.215) Speciálně pro prvek T3 a elementární matice a vektory určené v lineárním případě vztahy (2.40), (2.44), (2.49) a (2.54) jsou matice Hpe, Hqe, řťe, Ha5 a H'35 tvaru Rpe = Kel + KpeAe[ <9p(Pf) <9p(P£) <9p(P£) <9A^ ' dAe2 1 dA% (2.199') kde K pe 2|de -r — s tyč se re se \ _re _ fe fe ť -se - ť a kde p(Pf) je hodnota funkce p v těžišti Pf prvku e, \dq{F \ dA We = Ke2 + l |de| Ae í M^ľ) j ? (2.200') 6 '"'j ' i,j=l ^Íe = l\de\\d-W-X , (2-201') kde q(Pf) je hodnota funkce g ve vrcholu Pf prvku e, kde f(Pf) je hodnota funkce / ve vrcholu Pf prvku e, 2 ^A«| . (2.212-) ^ 3 > i, 1 = 1 Ha5 = K5 i 0a? i yS \ in Vi/-Vi^lv-i i-v-j- o din 1 mr\ r\< t tc\ n rvi n "OS kde a(Pi) je hodnota funkce a ve uzlu P* strany S*, ^ i j i,j=l kde j3(Pf) je hodnota funkce /3 ve uzlu Pf strany S. V úloze vedení tepla funkce p, q, f, a a /3 závisejí na prvku e obvykle pouze na teplotě u. V tom případě = |P'(|[A? + AI + AI]), Mgf) = q\At)5„ 3 3 df(Pt) fVAeU da(Pf) s d/3(Pf) s QAe = f (A,)%' gA5 = a (A, )% gA5 = p (a, 3 kde 5jj = 1 a ^ = 0 pro i ^ j. Jacobiovu matici HS(AS) zobrazení RS(AS), s = e,S, lze spočítat také jen přibližně. Užitím jednoduché formule numerického derivování spočteme i-tf sloupec matice HS(AS) například ze vzorce OlViA") II íA;. ...,Af + e,..., Asps) - IliA;...., A?,..., Asps) dAst e kde e je vhodně zvolené malé kladné číslo. Globální matici H=H(A^) a globální vektor R=R(A(-fc*)) sestavíme modifikací eli-minačního algoritmu uvedeného v části c) odstavce 2.7. Jediný podstatný rozdíl spočívá v tom, že prvky elementárních matic He resp. Hs přispívají pouze do globální matice. Konkrétně, je-li fof • resp. hfj prvek elementární matice He resp. Hs a jestliže lokálním indexům i a j příslušejí kódová čísla I a J, pak je-li J<0 nebo J<0, prvek fof- resp. neovlivní ani globální matici H ani globální vektor R. Jinými slovy, zpracovávají se jen takové prvky elementárních matic a vektorů, jejimž indexům příslušejí kladná kódová čísla: pro 7>0aJ>0se resp. hf, se přičte k hu a pro I > 0 se R\ resp. Rf se přičte 79 k Ri (h h jsou prvky matice Haf?/ jsou prvky vektoru R). Poznamenejme, že globálni matice H je obecně nesymetrická (jak je zřejmé například z tvaru (2.199')). Je známo, že Newtonova metoda s výjimkou speciálních případů konverguje jen tehdy, je-li počáteční aproximace A*-0-1 „dostatečně blízká" řešení A. Abychom překonali tento nedostatek Newtonovy metody, používá se její modifikace založená na strategii tlumení. V tzv. tlumené Newtonově metodě se místo soustavy (2.189) řeší soustava H(A(fc))Y(fc) = -wR(Aw), (2.216) kde uj < 1 je vhodně zvolený tlumící faktor zajišťující pokles absolutní hodnoty rezidua: ||R(A(fc+1))|| < ||R(A(fc))||, || • || označuje nějakou vektorovou normu, například délku vektoru. Jednu z možností výběru tlumícího faktoru navrhli Armijo a Goldstein (viz [19], také [12]): za uj doporučují vzít nej větší z čísel 1, |, ^,..., pro které platí ||R(A(fc+1))|| < (1 - |)||R(A(fc))||. (2.217) uj je tedy zvoleno tak, aby reziduum pokleslo alespoň 1 — f krát. Významnou vlastností této strategie je ta skutečnost, že pokud A^ konverguje k řešení A, pak uj —> 1, a tedy v blízkosti řešení je dosaženo kvadratického řádu konvergence jako u standardní Newtonovy metody. Jestliže však podmínku (2.217) nelze splnit pro uj > ujq, kde ujq je zvolená „dostatečně malá" konstanta, musíme bohužel konstatovat, že ani tlumená Newtonova metoda nekonverguje - počáteční aproximace A*-0-1 je prostě příliš špatná. Tento problém lze obvykle překonat metodou parametrizace, kterou popíšeme později. Lineární soustava rovnic (2.216) se často řeší iterační metodou, neboť zejména při menší nelinearitě je řešení této soustavy blízké 0, takže známe dobrou počáteční aproximaci pro iterační výpočet. Zvláště výhodné je použití zobecněné Ríchardsonovy metody (zkráceně GR-metody) eventuálně urychlené například pomocí zobecněné metody sdružených gradientů, viz [13]. V tomto textu se omezíme na použití „čisté" GR-metody. Pro lineární soustavu rovnic Ay = b lze GR-metodu popsat takto: 1) Je počáteční aproximace, r*-0*1 = b — Ay*-0-1; 2) pro j = l,... ,J počítej: Md(i) = r(j-i)j y(j) =y0'-i) + d0'); r(i) = j.Cj-1) _ Ad^, kde M je vhodná iterační matice splňující tyto 2 požadavky: 1) M je dobrá aproximace matice A, 2) řešení soustavy rovnic s maticí M je „laciné". 80 Všimněte si, že pro M = A je y^ řešením soustavy Ay = b. Ukažme si nyní, jak lze vhodně aplikovat GR-metodu na řešení rovnice (2.216). Nechť je nějaká aproximace matice H(A^) a C^Vl^ = je její LU-rozklad na dolní trojúhelníkovou matici £^ a horní trojúhelníkovou matici VSk\ Obvykle se volí Jí(fc) = H(A(ř)) pro vhodné / = l(k) < k. (2.218) Doporučený způsob výběru parametru l(k) uvedeme později. Y^ počítáme GR-metodou takto: y(°) = 0, = -WR(A«), y0') =y0'-i)+d0'); L j = i,...,j, (2.219) Ttí) = r0--D -H(A«)d« J Y(k) = y(J). Počet kroků J volíme malý. V dalším předpokládejme, že matice IK^*1 je určena vztahem (2.218). Pak pro J=l dostáváme známou modifikaci Newtonovy metody, v níž je určeno vztahem H(A(ř))Y(fc) = -wR(Aw). Věnujme se nyní volbě l(k). Pro l(k) = 0 je zřejmě Jí(fc) = H(A(0)) pro všechna k. V tom případě stačí v průběhu celého řešení rovnice (2.187) provést jen jediný Lf/-rozklad, a to právě matice H(A(-°-)). Pro k > 0 pak už provádíme jen „laciné" řešení soustav rovnic se stále stejnou dolní a horní trojúhelníkovou maticí. U této modifikace Newtonovy metody však lze očekávat jen pomalou nebo dokonce vůbec žádnou konvergenci. Proto se obvykle postupuje tak, že Lf/-rozklad provádíme častěji, například pro k = 0,m, 2m,..., kde m je zvolené přirozené číslo. Pak l(k) = m[k 4- m], kde symbol 4- značí celočíselné dělení. Tak například pro m = 1 je l(k) = k, pro m = 2 je 1(0) = 0, 1(1) = 0, 1(2) = 2, /(3) = 2, 1(4) = 4, /(5) = 4 atd. Jak metodu prosté iterace tak i tlumenou Newtonovu metodu je vhodné, a v zájmu dosažení konvergence často dokonce nutné, kombinovat s metodou parametrizace, viz [?]. K úloze R(A) = 0 přiřadíme třídu parametrických úloh íR(A, s) = 0 spojitě závislých na parametru s G (0,1) tak, aby íR(A, 1) = R(A) a aby úloha íR(A, 0) = 0 byla jen slabě nelineární nebo dokonce lineární. Zvolíme přiměřeně jemné dělení 0 = s0 < sľ < ■ ■ ■ < sq = 1 a řešíme postupně úlohy íR(A, sq) = 0 pro q = 0,..., Q. Řešení těchto úloh označme A*-9*1. Řešení A*-0-1 slabě nelineární úlohy získáme snadno. Pro q > 0 jsou již sice příslušné úlohy silněji nelineární, při jejich řešení však můžeme vyjít z velmi dobrých počátečních aproximací A*-9'0-1 = A*-9-1-1. Speciálním případem metody parametrizace je tzv. přírůstková metoda, běžně používaná například při řešení úloh nelineární pružnosti a plasticity. Ve standardní přírůstkové metodě se předpokládá, že „zatížení" F(A) = F na řešení A nezávisí, a že „míra nelinea-rity" matice K(A) závisí na „velikosti" ||F|| zatížení F. Parametrickou úlohu íR(A, sq) = 0 proto volíme ve tvaru K(A)A = F(9), kde F(9) = F(0) + s9(F - F(0)), (2.220) 81 a stejně jako dříve předpokládáme, že úlohu %(A, 0) = 0 příslušnou „zatížení" F^0^ umíme snadno vyřešit. Metodu označujeme jako přírůstkovou proto, že řešení A*-9*1 dostaneme z předchozího řešení A*-9-1-1 „přítížením" o přírůstek zatížení AF^ = — F^-1). b) Nestacionární úloha Omezíme se na nelineární nestacionární úlohu vedení tepla. Pak matice C, K a vektor F v úloze (2.68) jsou závislé na řešení A (ŕ). Užitím ^-metody opět dostaneme soustavu rovnic (2.69). Ta je však nyní nelineární, neboť matice Cí+e, Kí+Él a vektor Fí+e závisejí na Aí+e = A* + d(Aí+1 - A%): Ct+d = C(ti+e, Aí+e), W+e = K(ti+e, Aí+e), F*+d = F(ti+e, Aí+e). (2.69f') Aí+1 je tedy řešením nelineární soustavy rovnic íť+1 = R(Aí+e, Ai+1) = 0, (2.221) kde (2.222) R(w, A) = [C(ti+e, w) + AtmiU+e, w)] A- - [C(ti+e, w) - Att{l - 8)K(ti+e, w)]Ať - AtiF(ti+e,w). Nelineární rovnici (2.221) lze línearízovat: z lineární soustavy rovnic R(w, Y) = 0, kde í A* pro 9 ^ 1 nebo i = 0, w = < o ■ n -, 2 n 2.223 [ §A* - iA1-1 pro 9 = \ a i > 0, V ; spočteme Y a položíme Aí+1 = Y. Linearizované řešení lze případně vylepšit provedením několika kroků metody prosté iterace: zvolíme Y^0^ = Y, z lineárních soustav rovnic R(A* + 0(Y(j) - A*), Y(j+1)) = 0 (2.224) počítáme postupně Y^+1\ j = 0,..., J — 1, a nakonec položíme Aí+1 = Y^J^. Uvedená stategie bude úspěšná, je-li nelinearita poměrně slabá a časový krok At,i dosti krátký. Při silné nelinearitě, která je typická například u úloh se změnou fáze, nezbude než řešit rovnici (2.221) Newtonovou metodou, případně některou z jejích modifikací uvedených v části a) tohoto odstavce. Globální rezidum Rí+1 a jeho Jacobiovu matici Hí+1 opět sestavujeme z lokálních reziduí Rí+1'e a jejich Jacobiových matic Hí+1'e pro e G T a z lokálních reziduí a jejich Jacobiových matic FF+1,S pro S G S. Jedním ze sčítanců vytvářejících matici Hí+1'e je také matice Hí+1'ce příslušná funkci c, ol\s ) rs=1 kde Ci+e>e = Ce(ti+e, Ai+e'e), přičemž Ai+e'e = A*'e + 0(Ai+1'e - A*'e) a Aj'e jsou aproximace teplot Ae(tj) v uzlech prvku e v čase tj, j = i, i + 1. V lineárním případě pro prvek T3 je matice Ce uvedena v odstavci 2.10 a pro izoparametrické prvky snadno odvodíme Ce = f([NfcNe) = Je([Ňe]TceŇe|det Je|). 82 Proto kde matice Cce = K 96 je určena vztahem (2.196) a čl+e,e je hodnota funkce c(rr, y, t, u) pro x = xe(^,r]), y = ye(Í,r]), t = ti+e a u = J2l=i ^l+e'eŇ^(£, rj), přičemž A\+e,e jsou složky vektoru Aí+Él'e. Pro prvek T3 je ď[4+e'eA;+1'e] 6 H «A*+l,e Hí+1'ce = i \ď kde c*+ÉI'e = c(xer,yer,ti+Q, A%r+e,e). Matice Hí+1'ce je diagonální a jestliže funkce c závisí na prvku e pouze na teplotě u, pak její diagonální prvky jsou = \ \ď\ [c(A;+e'e) + 0Al+l^d{Al+e^)l r = 1,2, 3. 2.15. Konvekčně-difúzní úlohy s dominantní konvekcí Konvekčně-difúzni úlohy s dominantní konvekcí (stručně KDUsDK) samy o sobě, nebo jako podstatná součást úloh složitějších, jsou předmětem trvalého zájmu inženýrů i matematiků. Byla vymyšlena celá řada nejdříve jednoduchých a postupně značně rafinovaných metod a nové metody stále vznikají. Jde evidentně o živou problematiku, která je stále ve vývoji. Řešení KDUsDK je věnováno nesmírné množství odborných prací, z nichž mnohé lze najít například v knize [16]. My se v tomto textu omezíme na dvě základní techniky používané při řešení KDUsDK: pro stacionární úlohu uvedeme upwind metodu a pro nestacionární úlohu metodu charakteristik. Při řešení KDUsDK se jako základní diskre-tizační technika stále značně používá diferenční metoda nebo metoda konečných objemů na pravidelných sítích. Metoda konečných prvků se prosazovala poměrně obtížně, neboť specifické techniky typu upwind nebo charakteristiky se „roubují" na metodu konečných prvků obtížněji a nezřídka také méně efektivně. My však v duchu tohoto učebního textu zůstaneme věrni metodě konečných prvků, jen při upwind technice si vypomůžeme metodou konečných objemů. Obě techniky, tj. upwind a charakteristiky, vyložíme v nejjed-nodušší formě tak, aby vynikly jejich podstatné rysy. Základní myšlenky obou metod lze řadou důmyslných postupů výrazně zefektivnit. Jejich rozmanitá vylepšení zde však úmyslně popisovat nebudeme, neboť chceme udržet přiměřený rozsah tohoto odstavce a také nechceme čtenáře mást těžko sledovatelnými technickými detaily. a) Stacionární úloha, upwind metoda Hledejme funkci u(x,y) splňující podmínky (2.1), vyhovující diferenciální rovnici V • [tu] - V • [pVíí] + qu = f v Q (2.225) a splňující okrajové podmínky (2.3) a (2.4). Úloha popisuje například transport látky (chemické příměsi) v proudovém poli (tekutině). V tom případě je u(x,y) hmotnost příměsi v jednotkovém objemu tekutiny, r = q\ je součin hustoty g(x,y) a rychlosti v = (^1,^2) tekutiny s rr-ovou složkou v\(x,y) a y-avou 83 složkou v2(x, y), takže r\ = gv± je x-ová a r2 = gv2 y-ová složka vektoru r měrné objemové hybnosti, p(x,y) je koeficient difúze a f(x,y) — q(x,y)u(x,y) je objemový zdroj příměsi. Tatáž úloha popisuje také transport entalpíe v tekutině. Pak u(x,y) je entalpie (entalpie u a teplota T jsou svázány vztahem du = cdT, kde c je tepelná kapacita), r je opět roven součinu hustoty g a rychlosti v, p(x, y) = X/c je podíl tepelné vodivosti A a tepelné kapacity c a f(x,y) — q(x,y)u(x,y) je objemový výkon zdrojů energie. Obecně lze úlohu (2.1), (2.225), (2.3) a (2.4) považovat za matematický popis zachováni jisté fyzikální veličiny, která je ovlivňována difúzi, reprezentovanou difúzním členem —V • [p Vit], konvekci, reprezentovanou konvekčním členem V - [rit] a generaci případně absorpci, reprezentovanou zdrojovým členem / — qu. Dále budeme předpokládat, že vektor r = gv splňuje tzv. rovnicí kontinuity V-r^ + ^ = %^ + ^= 0. (2.226) Ox oy dx oy Pak V • [ru] = r ■ Vw = r\ux + r2uy a rovnici (2.225) lze zapsat v ekvivalentním tvaru Ou Ou -V • [pVu] + n— + r2— + qu = f v a (2.225') ox oy Abychom zaručili jednoznačnou existenci klasického řešení konvekčně-difúzní úlohy, doplňme soubor podmínek (2.5a) - (2.5d) o podmínku T_ = {(x, y) G díl | n(x, y) ■ r(x, y) < 0} C ľi (2.5e) zajišťující, aby na části T_ hranice, kterou do oblasti Q tekutina vtéká, byla zkoumaná veličina u známa. Navíc předpokládejme r1}r2 G C1(Í2). Míru vlivu konvekce a difúze na transport zkoumané veličiny vyjadřuje tzv. globálni Reynoldsovo číslo R = \\r\\£/p, kde ||r|| je délka vektoru měrné objemové hybnosti r, l je charakteristický rozměr oblasti Q a p je koeficient difúze. Pro „velká" Reynoldova čísla R (ve skutečnosti jde o hodnoty funkce R(x,y)) hovoříme o konvekčně-difúzních úlohách s dominantní konvekci: transport je rozhodujícím způsobem ovlivňován konvekci, vliv difúze je okrajový. Slabou formulaci odvodíme obvyklým způsobem a obdržíme úlohu najít u G W splňující a(u, v) + b(u, v) = L (v) \/v G V, (2.227) kde a(u,v) je určeno vztahem (2.9), L (v) vztahem (2.10) a b(u,v) je bilineární forma b(u,v)= I V-[ru]váxáy. (2.228) Jn Diskretizaci proveďme užitím prvku T3, takže a(u,v) aproximujeme pomocí dh(u,v), viz (2.19), a L (v) pomocí Lh(v), viz (2.20). Aproximaci konvekčního členu musíme provést obezřetně. Označme Re = ^ (2.229) P , 84 tzv. lokálni Reynoldsovo číslo. V (2.229) je h=he charakteristický rozměr prvku e. Pro jednoduchost předpokládejme, že he je nejdelší strana trojúhelníka e. Způsob diskretizace b(u,v) závisí na velikosti Re. Abychom si ozřejmili úskalí, která mohou nevhodnou diskretizací členu b(u,v) vzniknout, zabývejme se jednorozměrnou modelovou konvečně-difúzní úlohou -pu" + rv! = 0 pro x G (0,1), u(0) = a, u(l) = /3, (2.230) kde p > 0, t^O, aa/3 jsou konstanty. Přesné řešení je 1 — expf-rrl u(x) = a + (P - a)--f—. (2.231) Interval (0,1) rozdělíme na N stejných dílků délky h=l/N a úlohu (2.230) diskretizu-jeme. Metoda konečných prvků aproximující klasické řešení u pomocí funkce U po částech lineární vede na soustavu rovnic —Ui+i + 2Uí - Ui-i Uí+i - Ui-i P-^2--hr-2h-= U0 = a,UN = /3. (2.232) Tutéž soustavu rovnic obdržíme, když úlohu (2.230) diskretizujeme diferenční metodou, v níž difúzni člen — pu" aproximujeme standardně vztahem -pu \x=Xi~p--2- a konvekční člen ru' aproximujeme rovněž standardně užitím centrální diference ruf\x=Xi^rUi+1-hUi-\ (2.233) Jak se snadno přesvědčíme, řešení úlohy (2.232) pro Re = rh/p ^ 2 je , 1 D 1 'l Ui = A + B 2--n.e 1 — \Re (2.234) kde konstanty A a B jsou určeny okrajovými podmínkami. Pro \Re\ > 2 zřejmě Ui osciluje, což je v rozporu s monotónností přesného řešení u, viz (2.231), přičemž „velikost" oscilací je přímo úměrná \/3 — a\ a \Re\. Nežádoucí oscilace nenastanou, když konvekční člen aproximujeme pomocí jednostranné diference: Ui - Ui-x . . Ui - í/i_i r-ul ~ r--- = r--- pro r > 0, lx_ii h h (2.235) , Ui+i — Ui . .Ui — Ui+i rul w r- = r - pro r < 0. ix-xí ^ ' ' h Pro diskrétní řešení pak platí U = A++ B+(1 + Rey pror>0, (2.236) Ui = A~+ B- (1 - R^-' pro r < 0, 85 přičemž konstanty A+, B+, A~ a B~ určíme z okrajových podmínek. Aproximace pomocí jednostranné diference jez fyzikálního hlediska přirozená: informaci o řešení v uzlu Xi čerpáme ze znalosti řešení proti „proudu", proti „větru". Proto takovou jednostrannou aproximaci konvekčního člene nazýváme upwind aproximací a hovoříme o upwind schématu (upwind metodě). Řešení (2.236) je monotónní funkcí parametru i a je tedy „fyzikálně" realistické rovněž pro velká Re. Vraťme se nyní k diskretizaci konvekčního členu b(u, v). Je-li Re na elementech e malé, aproximujeme b(u,v) pomocí bh(U,v) standardním způsobem: bh(U,v) = Y,n[nUx + r2Uy]v) = ^[ee]TKe3Ae kde Ke3 = i \ď\ {aínixlvt) + 6|r2(x?,tf)}* (2.237) (index i je řádkový, j sloupcový, a|, b'j viz (2.37)). Poznamenejme, že prvky kff ele mentární matice K jsme získali numerickou integrací členů dwe, dwe. dx dy užitím formule (2.43). Matice Ke3 je nesymetrická. Nesymetrická je proto také lokální matice Ke=Kel+Ke2+Ke3 a globální matice K. Matice K je pro dostatečně malé h regulární. Pro velká lokální Reynoldsova čísla Re však řešením soustavy rovnic (2.24) většinou obdržíme neuspokojivé výsledky vyznačující se fyzikálně neopodstatněným a nesprávným oscilatorickým průběhem konečněprvkového řešení U. Diskretizaci konvekčního členu je proto třeba provést jinak. Postupovat budeme podobně jako v práci [6], tj. užijeme kombinaci metody konečných prvků, metody konečných objemů a upwind techniky. Ke každému uzlu Pj přiřadíme tzv. duální element Di (označovaný také jako konečný objem nebo stručněji box). Popišme si jeho konstrukci. Nechť e je jeden z elementů, které mají uzel Pj za vrchol. Označme zbývající vrcholy elementu e jako Pj, Pfc a dále označme střed strany PjP/ jako Pij, střed strany pP^ jako P& a těžiště elementu e jako Pijk- Pak průnik duálního elementu Di se zvoleným elementem e je roven čtyřúhelníku Dijk s vrcholy Pj, P^, P^k, Pik, viz obrázek 17. Duální element Di je tedy sjednocením čtyřúhelníku D^k příslušných všem trojúhelníkům e=PiPjPk obsahujícím vrchol Pj. Nechť S (i) označuje množinu všech uzlů sousedících s uzlem Pí (dva uzly jsou sousední, pokud jsou vrcholy téhož elementu). Ke každé straně PíPj, j G S (i), přísluší dva přímé úseky T Obr. 17. Dijk=D,i n e «=1,2, hranice dDi duálního elementu Pj: leží-li alespoň jeden z uzlů p nebo Pj uvnitř Q, pak jeden z těchto přímých úseků je úsečkou spojující body P^ a P^ a druhý je úsečkou spojující body P^ a P^, kde P^ je těžiště elementu s vrcholy p, Pj a p, viz obrázek 18; pokud oba uzly Pj a Pj leží na hranici díl, je jeden z těchto přímých úseků úsečkou spojující body P^ a P^k a druhý je úsečkou 86 spojující body P^ a Pí, viz obrázek 19. Označme alf, délku úsečky Tf, a nf- jednotkový Obr. 18. Di a D j pro vnitřní stranu PíPj Přistupme nyní k vlastní aproximaci členu b(u, v). Užitím Greenovy formule (2.6) a rovnice kontinuity (2.226) postupně dostaneme pu „ pu b{U,v)=YÍ V-[rU^dxdy ^Yv(pi) Í V • [r í/] dardy = i=i i=i Pí/ » = ľ^) V-[r{i7-f/(Pi)}]d^dž/ = i=i Pí/ « = VVpi) / {v-v)[U-U{Pi)]ás= (2.23Í i=i Pí/ 2 „ i=l i es (i) a=l jTij jeS(i) a=l "íij pu 2 í=i jes(i) ol=i 87 kde n je jednotkový vektor vnější normály díl a dfj'Bfj je vhodně zvolená aproximace jTa (r • n1j)[U — U (Pí)] ás. Tak například „přirozená" aproximace m = \r(Plj).^j]\U(PlJ)-U(Pl)] (2.239) pro U (Pij) = \U(Pi) + U(Pj)]/2 není pro velká Re vhodná a dává podobně neuspokojivé výsledky jako aproximace (2.237). Proto postupujeme jinak. Rozlišíme dva případy. 1) Pokud úsečka Tfj neleží na hranici, užijeme „upwind" aproximaci členu "Bfj, která ve vztahu (2.239) místo U (Pij) bere U (Pí) pro r(P^) -n?. > 0 nebo U(Pj) pro r(P^) -n?. < 0, takže = min{0, r(Pý-) • n£}[£/(P,) - í/(Pť)] P^ T" £ 0a 2) Leží-li úsečka Lf* na hranici, klademe £?• = 0 pro T" G 9a (2.240) (2.241) Zdůvodněni. Podíváme-li se na obrázek 19 vidíme, že k hraničnímu uzlu p přísluší dvě hraniční úsečky T2j a L2ř a jim odpovídající členy !Bf- a !B2ř. Z (2.238) plyne j2 ^2 d2rP>2 ~ Pii (r -n)[U- U(Pí)] ds. P 1 ?,i Integrál je však přirozené aproximovat nulou, což je formálně totéž jako položit 23 ®2 =0. □ P 23 Pí Obr. 19. ľ)i pro PP,- G <9$1 Obr. 20. Lokální značení 88 Člen bh(U,v) definovaný vztahy (2.238) a (2.240) vyjádříme pomocí elementárních matic ke3 pro e G T a K.S3 pro S G S: bh(U,v) = j][0e]tke3Ae. (2.242) Věnujme se dále vyjádření elementární matice ke3. Na trojúhelníku e užijeme lokální značení, viz obrázek 20. Vrcholy označíme Pf, P|, P3e, střed strany P1eP2e označíme P^ nebo P2i a podobně střed strany P2eP3e označíme P|3 nebo P3e2 a střed strany P3eP1e označíme P^ nebo Pfg. Těžiště trojúhelníka označíme Pf23- Nechť nf- je jednotkový vektor vnější normály trojúhelníka P^P^P^ na straně P^-Pf^- Označíme-li délku úsečky P^P^s jako df-, platí "i,- = -nj-i = ±(ž/ľ23 - Vtj, 4- - ^23)/^ (2-243) kde znaménko zvolíme tak, aby vektory nf • a P/PJ svíraly úhel menší nebo roven |, tedy aby nf • • PfPj > 0. Dále označíme m?. = min{0,d?.r(P^)-n?.}. (2.244) Pak z (2.238), (2.240) - (2.244) odvodíme, že elementární matice ke3 je tvaru k e3 / — mf2 — mf3 mf2 mf3 \ m21 ~m21 - m23 m23 V mll mÍ2 ~mll ~ mÍ2 ) (2.245) Všimněte si, že součet prvků v každém řádku matice ke3 je roven nule a dále, že na hlavní diagonále matice ke3 jsou prvky nezáporné a mimo ni prvky nekladné. Jestliže označíme k3 globální matici sestavenou z elementárních matic ke3, pak i tato matice má obě výše zmíněné vlastnosti: součet prvků v každém jejím řádku je roven nule, prvky na hlavní diagonále jsou nezáporné a zbývající prvky jsou nekladné. Proto platí ke3(A + c)=ke3A, kde c je vektor stejných prvků c, což je diskrétní analog vztahu V • [r(u + c)] = V • [rit]. Elementární matice ke je nyní součtem tří matic, ke=kel+ke2+ke3. Protože matice ke3 je nesymetrická, matice ke a globální matice k je rovněž nesymetrická. Lze ukázat, že pro dostatečně malé h]e matice k regulární. Použijeme-li trojúhelníkové prvky s netupými vnitřními úhly, matice ke2 volíme ve tvaru (2.44) a matice příslušné Newtonově okrajové podmínce volíme ve tvaru (2.54), pak podobně jako v poznámce 16 (odstavec 2.8) můžeme konstatovat, že matice k má kladné prvky na hlavní diagonále a nekladné mimo ni, že k je IDDM, tedy že k je M-matice, tedy že k je regulární nezávisle na velikosti h. Za předpokladu dostatečné hladkosti slabého řešení lze ukázat, že pro chybu u — U slabého řešení a jeho konečněprvkové po částech lineární aproximace, získané upwind aproximací konvekčního členu podle (2.238), (2.240) a (2.241), platí max \u — U\ = 0(h). n Poznamenejme, že pro malá Re je použitelná standardní aproximace podle (2.237) a v tom případě pro chybu u — U platí stejný odhad, jako kdybychom řešili úlohu bez konvekčního členu, viz poznámka 14 (odstavec 2.8). 89 b) Nestacionární úloha, metoda charakteristik Hledejme funkci u{x,y,t), (x,y) G íl, t G (0, T), vyhovující diferenciální rovnici + V-[tu]-V ■\pVu] + qu = f pro (x, y) G íl, t G (0, T), (2.246) splňující okrajové podmínky (2.64), (2.65) a počáteční podmínku (2.66). Nestacionární rovnice (2.246) opět charakterizuje zachování fyzikální veličiny u, při němž se kromě difúzního toku, konvekčního toku a objemové generace případně absorpce uvažuje také časová změna ut. Rovnice (2.246) je dobrým modelem nestacionárního šíření chemické příměsi nebo entalpie v proudovém poli. Stejně jako ve stacionárním případě budeme předpokládat platnost podmínky (2.5e). Rovnice kontinuity nabývá pro nestacionární úlohu tvaru dg „ dg dr\ dr2 dg d(gvi) d(gv2) — + V • r = — H--- H--- = — + + v = 0. (2.226') dt dt dx dy dt dx dy Pokud data g, r, p, q, f, g, a a /3 nezávisejí na čase t, pak ustálené řešení nestacionární úlohy, pro které ut = 0, vyhovuje rovnici (2.225). To umožňuje řešit stacionární úlohy technikami vyvinutými pro řešení úloh nestacionárních, tedy například metodou charakteristik, kterou nyní popíšeme. Postupovat budeme podobně jako v práci [3] a [20]. Ke každému bodu x = (x,y) G íl a času t G (0,T) přiřadíme tzv. charakteristiku X(r) = X(x, t;r) určenou jako řešení dvou obyčejných diferenciálních rovnic f r(X(r),r) pro X(r) G íl, , , . X(t) = < s počáteční podmínkou Xm = x. (2.247) L 0 pro X(t) G díl, Rovnici (2.247) řešíme „pozpátku", tj. pro r < t. Pro r1}r2 G PCľ(íl) je charakteristika určena jednoznačně. Všimněte si, že díky podmínce X(r) = 0 charakteristika nikdy neopustí íl. Charakteristika X(x, t; r) je tedy rovnicí křivky (trajektorií), po níž se pohybuje pomyslný hmotný bod, který je unášen spolu s tekutinou „rychlostí" r (ve skutečnosti je r součinem rychlosti v a hustoty g), a který v čase t dorazí do bodu x. V technické mluvě lze charakteristikou X(x, t; r) označit proudnici, která v čase t prochází bodem x. Pro časovou derivaci složené funkce u(X(t),t) platí Ľu du(X(t),t) _ du(X(r),r) mK ' dt dr Odtud a z (2.226') plyne, že rovnici (2.246) lze přepsat do tvaru ^ + r(x, t) ■ Vm(x,t). (2.24Č Ľu g^-V-[PVu] + qu = f, xGfi, ÍG(0,T). (2.249) Tuto rovnici diskretizujeme ^-metodou, tj. užijeme formuli typu --h Q Dt V At ti+i 90 pro 9 G (0,1). (2.248) aproximujeme takto: Dit — (x,ŕi+i) « [w(X(x,ŕm;ŕm),ŕm) - it(X(x, ŕm; ŕ;), íj)]/Z\íí, a protože z (2.247) plyne X(x, ti+1; ti+1) = x, dostáváme — (x,tť+1) « K+1(x) - uť(Xť(x))]/2\tť, (2.251) kde X*(x) = X(x, ŕj) a kde «•'(£) = u(£,tj) pro £ = x, X*(x), j = i, i + 1. Proto z (2.249) užitím (2.250) a (2.251) dostaneme -—-+ (1 - 9)Qt + 0Qt+1 = 0, (2.252) kde ui+1 « íi(x,ri+i), íľ « ^(X^x),^), Q-? = [^']_1{-V • [p>W] + gV - p} pro j = i,i + 1, a kde ^, p>, gJ a J-5 označují hodnoty těchto funkcí v (x,tj), u% »2 it(x, ŕj). Zdůrazněme, že zatímco ve vztahu (2.251) «•'(£) = it(£,i/), v rovnici (2.252) jsme tímtéž symbolem označili už jen aproximaci u(£,tj). Další aproximací (1 — 9)Ql + 9Qt+1 k, Qt+e, kde £Í+ÉI, pí+ÉI, gí+ÉI a jsou hodnoty těchto funkcí v (x, ti+g), ti+g = ti + #Aŕj, a ul+e = (1 — 9)ul + 9uí+1, dostaneme - V • [pt+9Vut+d] + qt+9ut+d = f+9 pro x G íl (2.253) Rovnici (2.253) doplníme okrajovými podmínkami v? = gj, pro x G Ti, j = i, i + 1, (2.254) l+ed^ = ai+eui+e _ p+e x G r (2>255) on kde aí+ÉI a /3Í+ÉI jsou hodnoty funkcí a a /3 v (x,ŕj+e) a kde i = 0,...,Q — 1 (Q zde označuje počet časových kroků). Úlohu (2.253) - (2.255) diskretizujeme v proměnné x metodou konečných prvků užitím elementu T3. Vypustíme slabou formulaci a zapíšeme přímo diskrétní slabou formulaci. Aproximujeme ú1 ~ U% = Yn=i ^ljwj i kde A* = Í7j je hledaná aproximace u(xj,yj,ti) (pro i = 0 klademe ř/j1 = Lp(xj,yj)), a dostáváme úlohu najít U*+1 splňující (gž+e^~^, ^ + ah(Ut+d, v) = Lh{v) Vv G \4 (2.256) kde PU ^(x) = ^f/>,(XÚx)), (2.257) přičemž X^(x) je přesné nebo jen přibližné, numericky spočtené, řešení X*(x)=X(x, ti+i; ti) úlohy (2.247), a kde dále 91 (w,v)h = YJF(wv), (2.258) Ui+e = (1 - 9)Ul + #£/m, (2.259) ah(w, v) = ľ(pl+eVw ■ Vv + qi+ewv) + ^ J^+V;), (2.260) eeT ses = £ P(f+0v) + £ J5(/3í+%). (2.261) seT ses Dosadíme-li do rovnice (2.256) za f/í+ÉI z (2.259) a převedeme členy obsahující hledanou funkci Ul+1 na levou stranu, dostaneme (gi+eUi+1,v)h + Ati6ah(Ui+1,v) = (gl+eÚ\v)h - Att(l - 9)ah(U\v) + AttLh(v). S metodou charakteristik souvisí pouze první člen na pravé straně této rovnice, zbývající členy jsou (až na jeden formální rozdíl) stejné jako v nestacionární úloze vedení tepla studované v odstavci 2.10 (funkci g zde odpovídá funkce c v odstavci 2.10). Je proto přirozené použít stejné formule numerické integrace jako v odstavci 2.10 a při výpočtu Ie(gl+eU'lv) užít formuli (2.43). Tedy Ie(pl+e^7w ■ Vi>) počítáme formulí (2.34), všechny zbývající integrály na elementu e počítáme formulí (2.43) a integrály na straně S počítáme formulí (2.53). Snadno odvodíme ľ(gt+dUlv) = [&e]TCe't+d, kde Ce't+d - — 6 / e,í+dTTe,í \ e'i+0 t jeA QÍ Uh,2 \ e,i+8fTe,i i Přitom gf+e = g(xej,yej,ti+e) a ř/JJ = U'l(X.eQ, kde X^. = X^(xp a x^e = {x^yf) jsou souřadnice bodu PJ, j = 1,2,3. X^'1- lze spočítat pomocí n > 1 kroků vhodné Rungovy-Kuttovy formule: pro 9 ý \ stačí metoda řádu 1 a pro 9 = | je vhodné užít metodu řádu 2. Určení U^%- je algoritmicky poměrně složitá operace vyžadující obecně nemalý objem výpočtů. Proto je účelné pro každý uzel P^ o souřadnicích = (xi,yi) spočítat X^ = X^(x^) a pak U%h^ = ř7l(X^) a na elementu e pro uzel P!f = Pg položit U^j = U%h^. Protože rovnice sestavujeme jen pro uzly neležící na Ti, Ulhí stačí určit jen pro P^ ^ Ti. Uveďme si nyní podrobně postup výpočtu Ulhí. Označme Ar = Ati/n a položme ts = sAr, s = 0,..., n. Dále nechť Y° = x^ pro zvolený uzel P^ ^ r\. Pak postupně pro s = 1,..., n počítejme: a) pro 9 ý 2 explicitní Eulerovou metodou (stručně EEM) řádu 1 Y| = Yf-1 - Arkf, kde kf = r(Y|-\ ri+1 - r^); 92 b) pro 0 = \ Heunovou metodou (stručně HM) řádu 2 Y| = Y,-1 - i^rCk?1 + kf), kde kf = r(Yf-\ŕi+1 - rs_!), kf = r(Yrx - Ark?\ŕi+1 - rs) Nakonec položíme = Y™ a vypočteme Úlkí = U%(X.he). Funkce r bývá často známa jen v uzlech P j a v časech t i. To nám umožňuje předpokládat, že r je pro (x, y) G ě a t G tvaru 3 í- — í t — t- r(x,y,t) = ^rej{ť)wej{x,y), kde r^ŕ) = r(^,y?,ŕj) + ^^r(^e,ŕm). i=i * * Při výpočtu Y^ je třeba vyčíslovat funkci r jednou nebo dvakrát podle toho, zda používáme EEM nebo HM. Abychom mohli kf případně k|2 spočítat, potřebujeme vědět, v jakém elementu e se bod Y|0 = Y^T1 případně Y^ = Y^T1 — Z\rkf nachází, a pak spočítat r(Yse5,ts5) = Ej=irKís5M(Y«) pro ó = 0,1 a pro ŕs0 = - *ai = - Ts-Podobně, chceme-li vyčíslit Ulhí = U%(X.he), potřebujeme znát element e, v němž bod leží, a pak spočítat U%(X.he) = Xľ/=i U%(xj,yj)Wj(X.he). Efektivní stanovení bodů Y|5 a výpočet hodnoty U%(X.he) je základem úspěšné implementace metody charakteristik. Určení Ulhí lze výrazně urychlit, když volíme Ati tak, aby byla splněna tzv. Courantova-Fríedríchsova-Lewyova podmínka (stručně CFL podmínka) _max ||r||Z\rj < min/i, ííx{0,T) ^ a když ~KZM určíme v jednom kroku EEM. Pak ~KZM = — Z\ŕjr(x^, í«+i) a tento bod leží v jednom z trojúhelníků s vrcholem P^. Výslednou soustavu rovnic pro výpočet neznámých parametrů Aí+1 lze zapsat ve tvaru [Ci+e + Ati0Kt+9]Ai+1 = Či+e - AUil - ^)Kí+eAí + AtiFi+e, (2.262) přičemž globální matici K.í+e a globální vektor Fí+e sestavíme pomocí elementárních matic Ke,í+e = Kei,i+e + Ke2,í+e^ z nichž Kei,i+e přísluší členu Ie(pt+9Vw ■ Vv) a Ke2'í+e členu Je(gí+É,-uw), elementárních matic K.s^+e příslušných členům Is{crJrewv), elementárních vektorů Ye,l+e příslušných členům Ie{fl+ev) a elementárních vektorů Fs,z+e příslušných členům Ps(j3í+ev), a dále globální matici Cí+e sestavíme z elementárních matic Qe^+e příslušných členům Ie{ol+ewv) a globální vektor Cl+e sestavíme z elementárních vektorů Qe,i+e příslušných členům Ie{nl+eU'lv). Matice soustavy (2.262) je symetrická, pro dostatečně malé h pozitivně definitní, a jsou-li vnitřní úhly triangulace netupé, je pozitivně definitní pro každé h. Pozitivní defi-nitnost matice soustavy (2.262) je velkou předností metody charakteristik, srovnáváme-li ji s metodami založenými na upwind aproximacích konvekčního členu, neboť ty vždy vedou k maticím nesymetrickým. Za předpokladu dostatečně hladkého slabého řešení u lze ukázat, že platí max \u{xj,yj,U) -U)\ = 0{h + Ať + h2/At), 93 kde q = 1 pro 9 ^ a q = 2 pro 9 = ^. Pro 9 = 0 je matice soustavy (2.262) diagonální a tedy výpočet Aí+1 je „laciný". I když je EEM jen podmíněně stabilní, pro „malý" difúzni člen není podmínka omezující délku časového kroku nijak zvlášť přísná (pro p C 1 jde prakticky o CFL podmínku), a proto se pro řešení KDUsDK někdy EEM používá. Diskretizujeme-li rovnici (2.246) standardně pomocí EEM, dostaneme rovnici / . - U'1 \ \Ql-2^->VJ +ah{U\v) + bh{U\v) = Lh(v), v níž lze konvekční člen bh(U\v) aproximovat upwind technikou, viz (2.242), a dostat tak zhruba „stejně kvalitní" výpočetní algoritmus jako je metoda charakteristik spojená s EEM. Proto je také tento postup někdy při řešení KDUsDK používán. 94 3. Prostorové úlohy Použití metody konečných prvků při řešení úloh ve trojrozměrném prostoru (stručně ve 3D) si budeme demonstrovat na problému popsaném stejně jako v kapitole 2 diferenciální rovnicí (2.2) a okrajovými podmínkami (2.3), (2.4). Nyní však je Q oblast ve 3D a u, V-, 0.1 í i Qi a a P Jsou funkce tří prostorových proměnných x, y a z. Také operátor V a vektor vnější normály n je třísložkový, V = (d/dx,d/dy,d/dz)T a n = (nx,ny,nz)T. Pro odvození slabé formulace potřebujeme Greenovu formuli. Ta má i ve 3D tvar (2.6), jen f a, g jsou funkce tří proměnných x±, X2, x%, dx = dr^id^d^ a ds je diferenciál plochy T. Vynásobíme-li diferenciální rovnici (2.2) testovací funkcí v G C1 (Q), v = 0 na Ti, integrujeme přes Q, použijeme Greenovu formuli a přirozenou okrajovou podmínku (2.4), dostaneme rovnici a(u,v) = L (v), kde bilineární forma a(u,v) a lineární funkcionál L (v) je určen vzorci (2.9) a (2.10), v nichž místo dxdy píšeme dxdydz. Slabá formulace je tvaru (2.12), jednoznačná existence slabého řešení je zaručena podmínkami (2.5a') a 3D analogií podmínek (2.5b), (2.5c) (gaa jsou funkcemi 3 nezávislých proměnných) a (2.5ď) (Q je mnohostěn s rovinnými stěnami, tedy polyedr, nebo obecněji „nedegenerovaný" mnohostěn se zakřivenými hladkými stěnami (tj. plochami se spojitě se měnící tečnou rovinou)). V MKP se při řešení 3D úloh obvykle používají izoparametrické prvky. Tvarově jde nejčastěji o šestistěny (speciálně o čtyřboké hranoly nebo dokonce o kvádry). Někdy je však k vykrytí oblasti nutné použít také pětistěny (speciálně trojboké hranoly) nebo dokonce čtyřstěny (speciálně trojboké jehlany). Protože šestistěn lze rozdělit na dva pětistěny a pětistěn zase na tři čtyřstěny, stačí pracovat jen s nej jednoduššími čtyřstěnnými prvky. Přesto se pětistěny a zejména šestistěny hojně používají. Na obrázku 21 je čtyřstěnný prvek T4 se čtyřmi uzly a s přímými hranami, na obrázku 22 je čtyřstěnný prvek T10 s deseti uzly a s křivými hranami: každá hrana je určena třemi uzly. (Označení čtyřstěnných prvků písmenem T je odvozeno od prvního písmene anglického slova „tetrahedron".) Obr. 21. Prvek T4 Obr. 22. Prvek T10 Na obrázku 23 je pětistěnný prvek P6 se šesti uzly, na obrázku 24 je pětistěnný prvek P15 s patnácti uzly, (písmeno P od prvního písmene anglického slova „pentahedron") na obrázku 25 je šestistěnný prvek H8 s osmi uzly a na obrázku 26 je šestistěnný prvek H20 s dvaceti uzly (písmeno H od prvního písmene anglického slova „hexahedron"). Referenční 95 prvek é je pro prvky T4 a T10 uveden a pro prvky H8 a H20 na obrázku 29. na obrázku 27, pro prvky P6 a P15 na obrázku 28 Obr. 23. Prvek P6 Obr. 24. Prvek P15 Obr. 25. Prvek H8 Obr. 26. Prvek H20 Na referenčním prvku é definujeme bázové funkce iV,f(£, r/, s vlastnostmi 1 pro i = j 0 pro i ý j a pro i,j = l,...,pe kde pe = 4 pro prvek T4, pe = 10 pro prvek T10, pe = 6 pro prvek P6, pe = 15 pro prvek P15, pe = 8 pro prvek H8 a pe = 20 pro prvek H20 a kde Q) jsou souřadnice 96 uzlu Pj. Pro prvek t4 je Ň? = l-Z-ri-(, pi (0,0,0) (0,1,0) p«e Ae = (Ío,o) „ P7e = (0,i,0) Pf pe ^2' 0, p\ (1,0,0) 4 - (0,0,1) (3.3.0) (0,0, i) (0,^3) 10 P|(1,0,0) Obr. 27. Čtyřstěnný referenční prvek pro prvek T10 je = 2(1 - ■Qíh-Z-v-O, Ňi = 2^ - ҧ = 277(77 Ňi = 2C(C - Ňl = 4£(1 -£-v -o, = 4^77, Ň<7 = 477(1 -£-v -o, Ňi = 4c(1 - Ňi = iv- ^10 = 4^c, »7-0, ne(o,i,i) Pf = (0,0, -1) pe _ (1,0, č? = (0,1, -1) pe _ (0,0, pe = (1,0, 1) pe _ (0,1, pe = (i,0, -1) pe _ f1 1 V2' 2' pe = (o,i -1) pe _ m.0 — (i o, pe -Ml V2- 2' 1) pe _ ^12 — (o,i pe ^13 = (0,0, 0) pe _ (1,0, pe ^15 = (0,1, 0) P|(l,0,-1) Obr. 28. Pětistěnný referenční prvek 97 pro prvek P6 je #f = ±//(i-C), #f = |0i + 0, pro prvek P15 je T'/(l + C), = (l-e-r/)(-e- Ňl = ^-K-1)(1 -0, Ni = ^ + K-i)(i + 0, = 2£(l-£-T7)(l -C), Ňl = 2r/(l-£-r/)(l -0, Ňl, = 2^(1 + 0, #13 = (l-£-r/)(l- #15 = >/(i-C2), 16, 17 -+- 13 20 I ÝP l l i pe pe M -Ml ------o- - ^(1,1,1) P 15 P\ 14 ^ pe ' r12 pe P 18 19 Pí P 10 Pí Pf(-1,-1,-1) P ^> 0 pro (£, ^, C) e é, je násada í>e vyjádřena formálně stejně jako souřadnice xe, ye a ze popisující geometrii, viz (3.1), tedy máme izoparametrický prvek. Oblast Q vykryjeme konečnými prvky (tj. triangulujeme ji). Použijeme buďto prvky T4, P6 a H8 (s uzly jen ve vrcholech) nebo prvky T10, P15 a H20 (s uzly také na hranách). Množinu T všech prvků nazveme triangulací, sjednocení všech prvků triangulace označíme ílh- Hranici oblasti Q h značíme ľV Částem r\ a T2 hranice T vhodným způsobem přiřadíme jejich aproximace a T2h-V moderních programech MKP je triangulace generována automaticky: uživatel specifikuje oblast, materiálové vlastnosti, okrajové podmínky a typ prvků, zbytek „zařídí program". Konečněprvkový prostor definujeme jako prostor funkcí v(x,y, z) definovaných na Q h takových, že jejich restrikce na element e je tvaru (3.5). Stěny prvků budeme značit S, množinu všech stěn ležících na hranici T2h označíme S. Prvky T4, P6, H8 s přímými hranami mají dva typy stěn, trojúhelníkové F3 a čtyřúhel-níkové F4, viz obrázky 30 a 31, A z r>S '°pf -i>r2 S 1 I 1 -^ 'ps y Obr. 30. Stěna F3 100 prvky T10, P15 a H20 s křivými hranami mají rovněž dva typy stěn, trojúhelníkové F6 a čtyřúhelníkové F8, viz obrázky 32 a 33. (Stěny značíme F podle prvního písmene anglického překladu „face" českého slova stěna.) Obr. 32. Stěna F6 Obr. 33. Stěna F8 Referenční prvek Š příslušný trojúhelníkovým stěnám F3, F6 a čtyřúhelníkovým stěnám F4, F8 je zakreslen na obrázku 34 a 35. Obr. 34. Referenční prvek stěn F3 a F6 Obr. 35. Referenční prvek stěn F4 a F8 Na referenčních prvcích definujeme bázové funkce: pro stěnu F3 pro stěnu F6 = 2(l-£-í7)(i-É-*7), Ňi = 4£(l-£ - v) Ňi = 2£(£-±), Ňi = 4£t7, Ňi = 277(77-i), Ňi = 477(1 - e - v) 101 pro stěnu F4 Ňf = \(1-0(1-V), N* = 1(1+0(1 + ??), Ňi = \(l+t)(l-v), Ňi = \(l-t)(l + v), a pro stěnu Fí iVf = í(l -v)(-€-v -1), Ng = 1(1 -a(i- - v) Ňi = í(l + 0(1 - v)(€ -v- 1), ҧ = 1(1 + 0(i- v2) Ňi = i(l + 0(1 1), Ň7S = 1(1 -a(i- iVf = Ki -0(1 -1), Ňi = 1(1 -0(i- v2) Geometrii stěny S určuje zobrazení Ps PS y = ys^r1) = Y.yiŇi^i wo(š,v)eš, (3.6) PS Q Z = z i=l kde (xf ,yf,zf) jsou souřadnice uzlu Pf, i = 1,... ,p$, ps = 3 pro stěnu F3, p,s = A pro stěnu FA, p,s = 6 pro stěnu F6, p,s = 8 pro stěnu P8. Na stěně S definujeme bázové funkce wf(x,y,z) = Ňts(^,r]), i = l,...,ps, (3.7) kde (£,rj) G Š je vzorem bodu [x,y,z) E S v zobrazení (3.6). Bázová funkce wf(x,y,z) je zřejmě restrikcí bázové funkce we-(x,y,z) pro P? = Pf na stěnu S* . Restrikci konečně-prvkové funkce v z prostoru na stranu S lze vyjádřit ve tvaru Ps v = vs(x,y,z) = Yvíwí(x,y,z), kde ví = v(xf,yf,zf), i = l,...,ps, (3.8) i=i pro (rr, y, z) G S* a ve tvaru PS i=i pro (£, ?y) G iS. Pro funkci g(x,y,z) definovanou na stěně S označme gs = gs(Š, v) = g(xS(Š, v),ys(Š, v), zS(Š,v)) Pro (£,v) e š. 102 K funkcím a a /3 přiřadíme interpolanty a1 a /37, jejichž restrikce aIS a j3IS na stěnách S G S jsou tvaru PS = ^a{x^,y?,z?)w?{x,y,z), PS í_1 (3.10) PS v ' í=l Protože 5-(x,í/,2;)dxdí/d2; = Jj^,r],()\det3e\ d^di]d(, přibližnou hodnotu integrálu na elementu e počítáme užitím kvadraturní formule Ie(-) na referenčním elementu é, Je(í?) = /e(r|detJe|). (3.11) Podobně postupujeme při přibližném výpočtu integrálu na stěně S. Využijeme vztahu / g(x,y,z)ds= / g(£,rj) J5(£,rj) d£dr/, J s J š kde Js = VEG - F2, dxs \ 2 ( dys \ 2 f dzs \ 2 E dxs dxs dys dys dzs dzs (3.12) <9£ di] <9£ di] <9£ drj ' dxs \ 2 f dys \2 í dzs di] J \ di] J \ di] a počítáme Is(g) = Is(gsJs), (3.13) kde Is(-) je kvadraturní formule na referenčním elementu Š. Diskrétni slabá formulace má tvar (2.18), množina přípustných řešení Wh a prostor testovacích funkcí Vh viz (2.15) a (2.14), bilineární forma a,h(U,v) a lineárni funkcionál L h (v) viz (2.160) a (2.161). Elementárni matice a vektory se získají zcela analogicky jako ve 2D případě, Kel = /e([Le]T[Je]-y [Je]~TLe|det Je|), (3.14) Ke2 = Je([Ňe]Tge Ňe|det Je|), (3.15) 103 Fe = Je([Ňe]T/e |det Je|), (3.16) Ks = jS^s^t&is jqsjs^ (3>17) Fs = p([Ň5]T/3/5J5), (3.18) kde Ňe = (Ň?,...,Ň°e), Le = (VŇ?,...,VŇ°e) přičemž V = (d/d£, d/&q, d/d()T, [Je]_T={[Je]T}_1 je matice inverzní k matici [Je]T a Ň5 = (Áf,..., Ň^). Zbývá uvést vhodné formule numerické integrace. Pro prvek T4 lze užít formuli Ie(g) = \g(\M) (3.19) nebo formuli Ho) = átô(0, 0, 0) + g(l, 0, 0) + g(0,1,0) + g(0, 0,1)]. (3.20) Obě formule (3.19) i (3.20) integrují přesně polynomy f^rfC^ pro 0