F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal, F. Miinz hemzal@physics.muni.cz Slabá formulace diferenciální úlohy Předpokládejme diferenciální rovnici v silné formě L(u) = 0. Ve slabé formě ji lze psát / L(u)dn = 0. Platnost rovnice v silné formě implikuje platnost formy slabé, opačně tomu rozhodně není. Aby slabá forma implikovala silnou, muselo by platit / vL(u)dn = 0, pro všechny (vhodně integrovatelné) funkce v. Matematický smysl slabé formulace spočívá v převedení hledané funkce u mimo diferenciální operátor L(u), takže podmínky na diferencovatelnost hledaného řešení jsou oslabeny a převedeny na podmínky diferencovatelnosti pomocných funkcí v. Formálně se převodu dá dosáhnout (opakovanou) aplikací per partes až k / uĽ^(v)dQ, = 0 (+ okrajové členy), Jo, kde Z/t je oprátor sdružený s L. V praxi můžeme vyvádění u zastavit, kde potřebujeme: například pro operátory Laplaceova typu, které jsou druhého řádu, zpravidla ponecháme požadavek na spojitost hledané funkce i funkcí aproximačních a per partes aplikujeme pouze jednou: / v/\dVt= I vVu-dS- / Vu-Vvdtt. Jn Jan Jn F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal, F. Miinz hemzal@physics.muni.cz Vratme se zpět k zavedení slabé formy. Jedná se o soubor nekonečně mnoha podmínek, ale můžeme se omezit na vybranou třídu funkcí v, v rámci které chceme slabé řešení optimalizovat. Pro potřeby MKP samozřejmě volíme v = ^2v[i]N\ i čili snažíme se slabé řešení optimalizovat na množině funkcí přípustných jako řešení problému MKP; dostáváme Y^v[i] J N[i]L(ú)án = 0. i Mají-li v[i] být libovolná, dostáváme tadiční Galerkinovy podmínky i : f N[i\L(u)dn = 0. Prověřme ještě vliv jednoduchých okrajových podmínek: okrajový integrál můžeme napsat jako / vVu ■ ndS, takže pro triviální Neumannovu podmínku na hranici, iin|ao = Vit^o • n = 0, nedostáváme žádný příspěvek. Dirichletova podmínka také nepřináší žádný okrajový příspěvek, ale z jiného důvodu. Po dosazení ansatzu MKP do obou vyskytujících se funkcí máme i : v / N[i\VN[j] -ndS. ^ Jan j U Dirichletovské hranice sice není derivace hledané funkce nijak omezena, ale sestavujeme pouze pro uzly i ležící uvnitř simulovaného objemu a hodnota jejich tvarových funkcí na hranici je nulová. F8370 Moderní metody modelování ve fyzice jaro 2023 Metoda konečných prvků ve 2D Věnujme se nejprve tvaru aproximačních funkcí N [i]. D. Hemzal, F. Múnz hemzal@physics.muni.cz Definice zůstává stejná: N[i] se skládá z po částech lineárních oblastí (v našem případě částí rovin) tak, že je jednotková v uzlu i a nulová na všech stěnách elementů obsahujícího uzel i, které samy uzel i neobsahují. Každá ze stěn N [i], označovaná jako NT[i] je tedy tvořena částí roviny ax + by + cz + d = 0, jejíž koeficienty jsou v rámci T{i,j, k] definovány z podmínek N\i N[í\\xm=0 N[i\\x[k]=0. (vzhledem k linearitě N [i] to zaručuje i nulovost na spojnici j-k) V praxi obvykle integrujeme v rovině xy, takže tvarové funkce volíme přímo jako NT[i] = ax + by + c. F8370 Moderní metody modelování ve fyzice jaro 2023 Soustava rovnic omezující NT[i] má tvar řešení Kramerovým pravidlem dává NT{i ax[i] + bx[i] + c = 1 ax[j] + bx[j] + c = 0 ax[k] + bx[k] + c = 0, [(ž/j - Vk)x + (íCfc - íCj)^ + (xjyk ~ xkyj kde IAtI = ISt souvisí s plochou elementu T a platí D. Hemzal, F. Múnz hemzalOphysics.muni.cz F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal, F. Miinz hemzal@physics.muni.cz Momentové integrály ve 2D Je vidět, že pro operátory Laplaceova typu (lineární, druhého řádu) každé dosazení do funkcionálu MKP povede na integrand ne vyšší než druhé mocniny v proměnných x, y, přičemž parametry integrálu budou známé hodnoty - souřadnice vrcholů elementů. To ovšem také znamená, že celé sestavování rovnic se rozpadne na výpočet mnoha integrálů typu kde a + b < 2. Tyto integrály nazýváme momentovými. Zjevně (0,0) = ST. Integrály přes prvek mohou zjevně být rozepsány jako Nalezněme nyní momentové integrály (trojúhelníkovým) segmentem T. (k, j) nad obecně položeným Xa Pl Xb P2 F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal, F. Miinz hemzal@physics.muni.cz Obecná přímka p zadaná v rovině xy body [x\, y{\ a [#2,2/2] může být zapsána ve tvaru p : 2/1-2/2 , %m ~ x2yi y =-x H-- X\ — Xi X\ — Xi což pro momentové integrály přináší V a ~Vc g. j xaVc~xcVa xa xc xa xc V a ~Vc g. j xaVc~xc.Va xa xc xa xc iT(kj) xky^ dydx + x kyj dyd x Xa Va-Vb x 1 ^aVb-xbya Xb Vb-Vc x 1 xbyc-xcyb xa x\) xa x\) x\) xc x\) xc Po provedení integrace můžeme všechny momenty až do druhého řádu psát 7T(0,0) = ST 7T(l,0) = ^(xa + xb + xc)ST ^(0,1) = ^{ya + yb + yc)ST T1 / (2, 0) = ~\Xa -\- Xu -\- Xc ~\- XaX\) ~\- X\)XC -\- XaXc)oT o 7T(0, 2) = ^(2/a + 2/6 + 2/c + VaVb + 2/62/c + yayc)ST JT(1, 1) = [2(xaya + xbyb + xcyc) + Xa{yb + 2/c) + ^6(2/a + 2/c) + xc(ya + 2/fe)J<$r Všimněme si, že momentové integrály byly upraveny až na tvar, ve kterém jsou zcela symetrické - tímto způsobem je překonána speciálnost volby polohy bodů při odvození momentových integrálů (například singularity, které by se projevily při svislé hraně v posledním obrázku). F8370 Moderní metody modelování ve fyzice jaro 2023 D. Hemzal, F. Miinz hemzal@physics.muni.cz Okrajové příspěvky ve 2D FEM Příspěvek okrajových podmínek Dirichletova a homogenního Neumannova typu v matici soustavy nevystupuje, v obecnějším případě jejich lineární kombinace jej lze nalézt explicitně. Předpokládejme smíšenou okrajovou podmínku ve tvaru au + un = 0, čili (aN{j]+VN{j]-n) = 0. Potom v rámci okrajového integrálu můžeme psát Tu[j] í NtWVNtÍj] ■ ndl = -aTuÍj] í NT^NT[j]dl. ■ JdT ■ JdT Vyhodnotme poslední integrál nad hranicí i-j elementu T(i,j, k). Na zvolené hraně / elementu platí 'X 'X _ _3 atT\ l 'X 1 'X l 'X 'X a •X i -I 7 % Oj J Oj J Parametrizaci hrany pak volíme jako x = Xj + {x i — xj)t, y = yj + {yi — yj)r, takže celkem dostáváme / NTm^j]dl = - ľ {X-X>KX-Xt)Ll3dT = -Lij ľ r(r - l)dr = ^ JdT JO \xi~xj) JO D kde Lij = \J (x i — x j)2 + {jji — y j)2 představuje délku hranice v uvažovaném elementu. Obdobně dostaneme ÍNT\i]NTm = LtJ ŕr2dr = ^f. F8370 Moderní metody modelování ve fyzice jaro 2023 Sestavovací procedura ve 2D D. Hemzal, F. Miinz hemzal©physics.muni.cz pro každý element T : pro každý vnitřní uzel i elementu T: 1) Ku + = IT 2a) je-li uzel j vnitřní: Ki« + = Iv=k je-li uzel j vnější: bi + = u[j]Iv=k 2b) je-li uzel k vnitřní: + = I V=J je-li uzel k vnější: bi + = u[k]L, další vnitřní uzel elementu T další element. Pořadí indexů v matici K je třeba nemíchat (index aktuálního uzlu je při přičítání neustále první). V kroku 2) je uzel v vždy ten, který nevystupuje v indexech matice K. Indexy i, j, k (tam kde se jedná o vnitřní uzly) nejsou pořadová čísla odpovídajících uzlů, ale pořadová čísla těchto uzlů v seznamu vnitřních uzlů.