F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz FEM ve 3D a v obecné dimenzi Začněme s konstrukcí tvarových funkcí NT[a] = ax + (3y + -{z + 8 ve 3D. Výchozími elementy T ve 3D jsou jehlany, definované čtyřmi vrcholy, takže podobně jako v případě nižších dimenzí dostáváme soustavu Xb xc \Xd Va Vb Vc Vd Za Zc Zd 1\ 1 1 1/ M 7 w /1\ o o w Tento tvar není přímo vhodný pro zápis v obecné dimenzi v tom smyslu, že ji s přibývající dimenzí nelze přímočaře rozšiřovat o nové sloupce a řádky. Protože nás však bude zajímat především determinant této matice, můžeme si dovolit manipulovat s jejími řádky, například do tvaru A = Xb Xb Xc Xd V a ~ Vb Vb - y c y c - y d yd Za Zb Zc Zb Zc Zd Zd 0\ 0 0 1/ Nyní je rozvojem podle posledního sloupce patrné, že / A Xa Xb \XC Xb Xc Xd ya yb yc yb yc yd Za Zb Zc a tato matice se dá s přibývajícími dimenzemi rozšiřovat velmi pohodlným způsobem. Jako bonus stačí počítat determinant z matice řádu o jedna menšího, než ve výchozím tvaru. Navíc existuje obecný vztah mezi objemem d-rozměrného simplexu a maticí A, V(T) = i|A|, jak se dá snadno ověřit na ID a 2D případu. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Múnz hemzal@physics.muni.cz K vyjádření tvarových funkcí NT[a] = ax + (3y + + ô je ještě třeba pomocí Kramerova pravidla spočítat determinanty Aa až Ag, nejpohodlněji asi v původní soustavě; například po záměně prvního sloupce dostáváme Aa = 1 Va Za 1 0 Vb Zh 1 0 yc zc 1 0 y d zd 1 Zb 1 Vb Zb - zc 0 Ve Zc 1 — Vc - y d zc ~ zd 0 Vá Zd 1 Vd Zd 1 yb -yc zb- zc y c -vá zc- zd ve kterém poznáváme plochu projekce stěny protilehlé vrcholu a do roviny y z. Obdobně dopandou A^ a A7, pro As dostaneme Aa = Xa xc Xd Va Vb Vc Vd Za Zb Zc Zd Xb Xc Xd Vb Vc Vd Zb Zc Zd Můžeme tedy tvarovou funkci přehledně zapsat jako 1 NT[a A(T) [A^x(^a)x + A^y(^a)y + A^z(^a)z + A$(^a F8370 Moderní metody modelování ve fyzice jaro 2021 FEM, pokročilé vztahy Vratme se k (pro jednoduchost dvojrozměrnému) integrálu D. Hemzal, F. Múnz hemzal@physics.muni.cz J J N[a]dxdy, pro tvarovou funkci N [a] = ax -\- f3y-\- 7 nad elementem T. Až dosud jsme tento integrál rozdělovali na momentové příspěvky Jj N[a]áxáy = a/T(l,0) + /3/T(0,1) + 7/T(0, 0). Pokusme se nyní o explicitní dosazení všech komponent; především, pro element s vrcholy a, 6, c máme A = xayb + xbyc + xcya - xayc - xbya - xcyb, dále Aa = (yb — yc), = [xc — xb), A7 = (xbyc — xcyb) a pro zúčastněné momentové integrály platí /T(o,o) = IAI 7T(1, 0) = ^-(xa + xb + xc) IAI 7T(0,1) = ^{ya + Vb + Vc) IAI Po dosazení tedy celkem dostáváme ärľ 1, , (Vb ~ yc)(xa + xb + xc) + (xc - xb)(ya + yb + yc) + 3(xbyc - xcyb) |A| 1 T7V[a]da;d^ =-3Ä-T = 6|A| a k výpočtu tohoto typu integrálů tak lze vlastně náročné sestavení momentových integrálů zcela obejít. F8370 Moderní metody modelování ve fyzice jaro 2021 Obecněji se dá ukázat, že nad libovolným elementem T(a,b,c) platí (Ar[o])ť(JV[6]y(JV[c])*dxd» = (i + ^!+2). D. Hemzal, F. Miinz hemzal@physics.muni.cz AI tento vztah lze poměrně jendoduše zobecnit do obecné dimenze (povšimněme si, že nad 2D elementem - trojúhelníkem - se vyskytují nejvýše tři různé N\ IAI Další typ integrálu, který se ve fyzikálních úlohách vyskytuje, je JJtVN\i] ■ VN\j]áxáy = (a[í\a\j] +/3[í\P\J V integrandu může také vystupovat libovolná funkce; vyřešíme ve 2D případu integrál 't f(x,y)(N[a]y(N[b]y(N[c])kdxdy. Chceme-li postupovat v souladu s myšlenkami MKP, vyjádříme vystupující funkci pomocí jejího FEM ansatzu f(x, y) = J^i /[fl-Wlfl? takže celý integrál má tvar £/M / N[í\(N[a])\N[li[y(N[c])kdxdy l Jt a každý maticový příspěvek by tedy měl vznikat sumováním / přes všechny elementy. Z vlastností tvarových funkcí ovšem samozřejmě vypývá, že do celé sumy přispějí jen ty N[l], jejichž / náleží aktuálnímu elementu T. Celkem tedy If = [(i + l)\j\k\f[a] + + l)\k\f[b] + i\j\{k + 1)!/[C]] -——. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Shrnutí V rámci zobecňování zaveďme nová označení pro vrcholy (a,b,c, ...) elementu T a jejich souřadnice v prostoru s dimenzí d. 1) Oindexujme nejprve vrcholy římskými indexy: ai, i = l..d + 1, souřadnice oindexujme řeckými indexy a G (0, d): x[0] = 1, x[l] = x, x[2] = y, atd. Pro souřadnice vrcholů tak celkem dostáváme výrazy typu ai[a] a je například di[0] = 1, 02[3] = bz atd. 2) Dále zaveďme nové označení I[a, (3] momentových integrálů až do druhého řádu včetně: indexy a, (3 G (0, d) nyní udávají, které souřadnice vstupují do výpočtu momentu: 7[a>/3] = JJ x[a]x\p]&n, Porovnáním s původním označením máme například 7[1,1] = 7(2,0), /[l, 3] = 7(1,0,1), 7[0, 2] = 7(0,1), apod. Výhodou nového zápisu je, že je schopen jediným kompaktním zápisem 'd+l d+1 d+1 lT[a,P] = (d+l)(d + 2) ^ a>% [«] • ai 1^1 + ai Ialai 1^1 IA .i=l i=l d\ postihnout všechny momenty do druhého řádu pro libovolnou dimenzi úlohy. Přitom oi[l] -a2[l] ai[2] -a2[2 o2[l]-o3[l] a2[2]-a3[2 K vyjádření tvarových funkcí je ještě třeba pomocí Kramerova pravidla spočítat ještě částečné determinanty, d NT[a] = A t A0(-a) + ^A^a) kde Ar b[l] b[i c[l] c[2 Pro integrál z tvarových funkce nad elementem T pak platí 7! T(iv[i])-(iv[2]r2...d^=(|/| + d)! A kde pro multiindex 7 = (il, i2,...) platí 7! =Uj ij\ a |7| = ■ i,-. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Na závěr zbývá vyjádřit momentové příspěvky okrajových integrálů N[j]VN[i\dW. Pro jednoduchost demonstrujme postup na 3D příspěvku QT(i,j, k) = J J xYzkáS- na takový integrál vede nehomogenní Neumannova podmínka i podmínka smíšeného typu. Je-li hraniční stěna parametrizována jako z = ax + f3y + 7, můžeme shora uvedený integrál (prvního druhu) napsat jako Qt(íJ, k) x y Jzk\ll+ (^f) + ) dxdy, ox J \oy kde koeficienty parametrizace stěny se určí opět Kramerovým pravidlem a = aa/a, atd. z f Za Zb Ixa ya lN Xb Vb \xc Ve 1, \7, Potom 1 + (I) + (i) = yrwTTP=^(^r + (a^)2 + i^f- Integrujeme tedy vlastně přes průmět stěnového elementu (zde trojúhelníka) do roviny x-y, jeho skutečný sklon je pak zohledněn vystupující odmocninou. To s sebou přináší nepříjemnost v tom smyslu, že i když stěna má zaručeně nenulovou plochu, její průmět již může být degenerovaný do úsečky (pokud bude svislá) a potom A-,z = 0 a shora uvedený vztah bude (zatím) divergovat. Problém se vyřeší úplným vyhodnocením vystupujících členů. Vyhodnotme nejprve momenty typu Qt(í,J,0); pro ně zjevně platí Qt(í,3,0) = j^—J(A^y + (A^y + (A^ í íxVdxdy = ^(A^z)2 + (A^ + (^IzgíŮ a vztah bude vždy nedegenerováný. Dosazením z = ax + (3y + 7 pak pro momenty až do druhého řádu podobně dostáváme Qt(0, 0,1) = ^(A^)2 + (A^)2 + (A^,)2 ±(za + zb + ZC) Qt(1, 0,1) = \J(A^z)2 + (A^)2 + (A^y)2 ^ [2(xaza + xbzb + xczc) + xa(zb + zc) + xb(za + zc) + xc(za + zb) Qt(0, 1,1) = ^/(A^z)2 + (A^)2 + (A^y)2 ^ [2(2/a2;a + 2/6z6 + yczc) + 2/a(z6 + zc) + 2/6(za + zc) + yc(za + 2:5) QT(0, 0, 2) = ^/(A^)2 + (A^)2 + (A^y)2 -^(^ + zl + zl + zazb + *6*c + zazc). 12 K divergencím již dojít nemůže a navíc jsou 3D okrajové momenty plně vyjádřitelné pomocí jejich 2D objemových variant: r2D| kde ovšem za a, (3 může libovolně figurovat také z. Uvedený vztah je pak snadno rozšiřitelný do obecné dimenze.