F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Münz hemzalOphysics.muni.cz Metoda konečných prvků (Finite Elements Method) metoda konečných diferencí: - průběh funkce je přesný - vyjádření derivace je přibližné metoda konečných prvků: - průběh funkce je přibližný - vyjádření derivace je přesné výhody MKP: dobře popisuje složité oblasti nevýhody MKP: těžko se získávají testovací 'analytická' řešení Zkoumaná oblast je rozdělena na elementy TJ; v rámci každého z elementů se fyzikální parametry úlohy zpravidla předpokládají konstantní (na (společné) hranici dvou elementů tedy často mají skok). Hledaná veličina ^(x), se diskretizuje jako ip[i] v uzlech v místech vrcholů elementů TJ. Uzly, které jsou prvky hranice díl nazýváme hraniční, ostatní uzly označujeme jako vnitřní X Xe Xö X4 X3 Xs Xi- 26 27 21 t 16 \ 11 i g i 7 22 17 12 7 2 28 _29 .30 25 23 18 13 8 24 19 14 20 15 10 Í>i Í>2 f 3 f 4 f 5 F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Aproximace se provádí s pomocí tvarových funkcí N [i] (lineárních, kvadratických ... obecně řádu k) vztahujících uzel i s přilehlými elementy T všeobecně platí, že se snažíme používat nejnižší řád aproximace, který ještě umožní najít netriviální řešení zvolené diferenciální rovnice; například pro rovnice druhého řádu stačí k = 1. lepší aproximace dosáhneme překročením nej nižšího řádu, v praxi často znamená možnost hrubšího dělení íž; jednotlivé elementy vš; získávají větší počet parametrů a tak výpočetní úspora nemusí b velká. řádům k > 3 se vyhýbáme kvůli jejich tendenci oscilovat. a následně vytvořených aproximačních funkcí N [i] jako sjednocení %--všech NT[i] zkonstruovaných z tvarových funkcí, které zasahují do elementů T obsahujících uzel i. u zjevně tedy N[i] a N[j] mají nenulovou společnou podporu pouze pokud uzly i a j náleží stejnému elementu. F8370 Moderní metody modelování ve fyzice D. Hemzal, F. Miinz jaro 2021 hemzal@physics.muni.cz Metoda konečných prvků v ID V případě ID úlohy vznikají elementy T prostým (byt nehomogenním) dělením intervalu G R na subintervaly. I wi wj ní i U všech typů elementů je třeba dbát, aby byly nedegenerované, tj. aby JT díl ^ 0, v ID se podmínka redukuje na vyloučení subinter-valů nulové délky. o Lineární tvarové funkce N~[i], resp. N+\i] spojují uzel i s ne-jbližším sousedem vlevo a vpravo a mají formu úseček nabývajících hodnoty 1 v uzlu i a hodnoty 0 v uzlu i — 1, resp. i + 1: N~[i\ :y= X~X][l~ \ xe (x[i-l],x[i\) x[i\ — x [i — 1J :y= X^x][\+\ xe{x[í[,x[i + l]) x[i\ — X[l + 1J Aproximační funkce N [i] je potom sjednocením N [i] = N~[i] |J N+[i]. U vyšších řádů tvarových funkcí se pro každý element přidávají vnitřní parametry, ovlivňující chování složitějších křivek na každém z intervalů (počet přidaných parametrů je dán balancováním stupňů volnosti). F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Jádro metody vzniká spojením dvou fundamentálních kroků: 1) uvědoměním, že hledanou funkci ty můžeme (ve spojité oblasti) vyjádřit jako ty(x) = ^2^[i]N[i] pouze s využitím hodnot ty [i] v diskretizovaných uzlech. zároveň, díky volbě aproximačních funkcí nezasahujících za nej bližších k sousedů, se na příspěvku v daném intevalu (x[i], x[i + 1]) podílí jen k tvarových funkcí (resp. k jejich aproximačních větví) přitom ty [i] jsou z hlediska diferenciálního operátoru konstantami. Skutečně, pro k = 1 mezi dvěma body platí Tímto způsobem je do formulace vnesena nespojitost (v prvních derivacích pro k = 1, atd.), ale metoda jako celek je na to připravena druhým krokem, 2) reformulovaním tzv. slabé variační úlohy pro řešenou soustavu rovnic ty(x) ty [i] - ty[i + l] x[i + 1] — x[i] X + ty[i]x[i + 1] - ty [i + l]x[i] x[i + 1] — x[i] F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzalOphysics.muni.cz Slabá formulace variační úlohy Budeme se zabývat řešením Cauchyova problému pro lineární diferenciální operátor L (obsahující derivace) fc-tého řádu. Uvažujme tedy oblast s hranicí díl a na ní rovnici Lu = 0, včetně příslušných okrajových podmínek /(w, Vit,... Aplikovat metodu konečných prvků jako Lu = 0 prostřednictvím u = w[ž]7V[i] formálně nelze kvůli nespojitosti tvarových funkcí. Mohli bychom ovšem uvažovat slabší podmínku [ L(u)dn = 0. Jelikož se jedná o integraci v Riemannově smyslu, jsou povoleny nespoj itosti na množině míry nula. Řešení původní rovnice splní i slabší intergrální podmínku, opačně to ale neplatí. Z důvodu aproximace MKP slabší řešení navíc původní formulaci vyhovět ani nemůže a musíme obecně očekávat R = L(u) Ý 0. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Máme pochopitelně zájem o co nejlepší splnění výchozí rovnice, a v tomto ohledu exituje několik postupů, které je možné aplikovat: metoda nejmenších čtverců. Požadujeme mm, 'o a aplikací standardního aparátu metody nejmenších čtverců dostaneme dR J ■ 2 JnRdu[j díl = 0. metoda vážených reziduí. Reziduum se minimalizuje vzhledem ke vhodným váhovým funkcím podmínkami J ■ 'n WjRdtt = 0 metoda nejmenších čtverců: Wj = ^_r ou[j metoda kolokační: Wj = (5(|x — x[j]|) - fixujeme přesné splnění rovnice v uzlových bodech metoda Galerkinova: Wj = N\j] - nejlepší aproximace na (pod)prostoru span(7V[j F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz MKP — příklad sestavení v ID, Poissonova rovnice Uvažujme Poissonovu rovnici L(u) = Au - 4np = 0. Je-li u = ^2iU[i]N[i], dostáváme pro Galérkinovy podmínky J- í N[j]L (Tu\z]N[i\) = Tu[i\ í N[j]V(VN[i\)dn-A7T í N[j]pdn = o, kde ovšem můžeme psát N\j]V(VN[i]) = V(N[j] V N [i]) — VAT [7] VAT [i] a na první člen použít Gaussovu větu. Dostaneme j : ľ«[i] í N[j]VN[i\dW ~y2u[i\ í VA^[j]VA^[ž]d^ = 4tt í N[j]pdn, i JdťL i Jn Jn což je algebraická soustava Ku = b, kde příspěvek objemových členů je Kíj= í VN[j]VN[i\dn bj = 47T í N[j]pdn. Omezíme-li se Dirichletovými a homogenními Neumannovými podmínkami, nepřináší okrajový integrál žádný vklad. Přesněji, triviální Neumannova podmínka nevyžaduje žádnou modifikaci rovnic, podmínka Dirichletova svými hodnotami vstoupí do integrálu u okrajových elementů, ale jelikož je hodnota na hranici zadána, všechny takové členy se převádí jako konstanta na pravou stranu (v samotných Dirichletovsky zadaných uzlech rovnice samozřejmě nesestavujeme vůbec). F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz V uvedeném zápisu by sestavení rovnic procházelo ve velké smyčce všechny uzly, v každém cyklu by se dohledávaly všechny elementy, které uzel obsahují. To není praktické, zejména ve vyšších dimenzích výstup z generátoru sítě obsahuje primárně seznam elementů a jeho neustálé prohledávání pro aktuální uzel by bylo extrémně zdlouhavé. Problém má jednoduché řešení tkvící v samotném zavedení MKP, kdy (pro k = 1) tvarové funkce nezasahují mimo jeden element. Například diagonální prvek Kjj = / \7N[j]\7N[j]díl má dva příspěvky, ŕ X [j ] ŕ X [j -\-1 ] Kjj= / (VN-{j])2dx+ / (VN+[j])2dx, Jxlj-l] Jxlj] -x[j] ľx[j+l] (ViV-[j])2da;+ / 'x[j-l] Jx[j] které lze pohodlně vygenerovat při procházení jednotlivých elementů, Ku = -—;-K-.-tt + 1 LJJ x[j]-x[j-l] ' x[j + l}-x[j\' Podobně (pro i ^ j) bude Kíj = / VN[i]\7N[j]díl nenulové pouze pro i = j ± 1, /»_|_ 2_ j -j /* ^ [«7 ] 1 K™=L VAr+[j]w'[j+l]d^-^+i]-^ľ k^=U^-^^=-^w. Základní kontrolou našeho sestavení je symetrie matice K, v obecnějším případě bývá matice soustavy před zavedením okrajových podmínek hermiteovská. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Komponenty vektoru pravé strany bj = 4-7t / N[j]pdíl se rovněž rozpadají na jednoduché příspěvky rx\j] 'x[j- px\j+l] N~[j]pdx + 47t / N+[j]pdx, 11 Jx\ '^^x[j]-x[j -1} x\j + l] -x[j P[a\-ô--h PiP\-ô- přičemž způsobů, jak konkrétně zapracovat hodnoty p je více. Můžeme například považovat p za konstantní na jednotlivých elementech, potom ŕ X [j ] ŕ X [j -\-1 ] bj = 47ľp[a] / N~[j]dx + 4irp[/3] / N+[j]dx = 4tt Jx[j-1] Jx[j] Nebo můžeme, více v souladu s myšlenkou MKP, uvažovat p = p[i]N[i], odkud ŕ X [j ] ŕ X [j -\-1 ] 6j=47rVp[i] / N-\j]N[i]dx + 4Tr Vp[i] / 7V+[j]iV[ž]da;. Nyní opěr s výhodou využijeme nepřekrývající se aproximační funkce, např. / iV-[j]iV[i]da; = pL/] / (iV-[j])2da; + pL/-l] / N~[j]N+[j - l]dx Jx[j-1] Jx[j-l] Jx[j-l] a / iV+[j]iV[i]da; = py] / (iV+[j])2da; + p[j + 1] / 7V+[j]7V-[j + l]da; i Jx[j] Jx[j] takže celkem b\j] p[j]xti +11 ~ ~ ^ + p[j - i]x^ ~ ~ ^ + ^[j + r +1' ~ x^ 4tt 3 6 6 Vypišme si všechny získané koeficienty ve zjednodušením případě ekvidistatního kroku e: 6 S e z pravé strany rozeznáváme v sestavené matici diskretizovaný Laplasián, na pravé straně zbývá váhovaný průměr hodnot v okolí aktuálního uzlu, přičemž většina váhy spočívá na hodnotě v aktuálním uzlu. Důvod této změny (ve FD by vystupoval pouze aktuální uzel) tkví v Galerkinově optimalizaci - rovnici se snažíme řešit nejen v uzlových bodech, ale i v oblasti mezi nimi.