F8370 Moderní metody modelování ve fyzice jaro 2021 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 2021 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 2021 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 2021 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, [fe - Vk)x + (íCfc - ícj)^ + (xjyk - xkyj kde lAy! = 2St souvisí s plochou elementu T a platí D. Hemzal, F. Múnz hemzalOphysics.muni.cz F8370 Moderní metody modelování ve fyzice jaro 2021 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 2021 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 2021 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 2021 Příklad - Helmholtzova rovnice D. Hemzal, F. Miinz hemzal@physics.muni.cz Uvažujme stacionární formu skalární homogenní vlnové rovnice AU-^—Y = 0. c2 dť Stacionarita předpokládá časově ustálený stav, o kterém u vlnové rovnice víme, že je realizován harmonickým řešením U(x,t) = u(x) ex.p(iujt). Dosazením tohoto ansatzu vskutku odstraníme z rovnice časovou závislost a dostáváme Helmholtzovu rovnici . u2 n Au H—;ru = 0. cz Protože je vlnová rovnice lineární, lze se na Helmholtzovu rovnici dívat také jako na rovnici pro jednu frekvenční kopmonentu časově proměnného pole, nebot obecně U(x,t) = uqj(x) exp(ia;t) Obecné časově proměnné pole pak můžeme získat tak, že provedeme větší množství stacionárních simulací pro různé frekvence a ty potom v prostoru složíme s váhami, které získáme jako komponenty Fourierova rozvoje budícího pulzu. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Sestavme nyní slabou formulaci úlohy s Helmholtzovou rovnicí: J (ku+^u ) díí = 0. .2 i JI Dosazením operátoru L{u) = Au + {uj /c )u do funkcionálu MKP, přičemž u(x) = ^u[i]N[i], dostáváme j : Vu[Í\ / N[j]VN\z]dW-T4 oj VN\z]VN[j] - -^N[j]N[i dn = o. První integrál opět představuje příspěvek od okrajových podmínek a druhý od samotného operátoru Helmholtzovy rovnice. K diagonální členům Ku budou přispívat všechny elementy, obsahující uzel i. K nediagonálním členům Kíj přispějí vždy pouze dva elementy, které oba uzly i, j obsahují. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Sestavení rovnic ve 2D Pro diagonální integrál Ku platí pricemz Ku = y, / (VArw • VArw - t^IW]) d<^ TkeSUPpN[i\ rp J (^VNTk[i] ■ VNTk[i ^NTk[i]NTk[i] ) áS = i - yr)2I(2, 0) + (xr - xi)2I(0, 2) + 2(2/1 - yr)(xr - xi)I(l, 1)+ + 2(yi - yr)(xiyr - xryi)I(l, 0) + 2(xr - x{)(xiyr - xryi)I(0,1)+ + [(xiyr ~ xryif - ^[(yi - yrf + (xr - xi)2]]l(0,0) j, kde uzly /, r jsou v tomto případě zbývající dva uzly vyhodnocovaného elementu, platí Xfc(i, /, r). F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Obdobně, pro nediagonální členy platí Kíj= / + Odtud, ^■nPnP \dn = il + ir c 2 * 3 L j Rj L, = 2 f ? oj 2SV — <\(Vv- Uj)(yi - ž/u)/(2, 0) + (xj - xv)(xv - ící)7(0,2) + c + [(Uv - Uj)(xv - x i) + (yi - yv)(xj - xv)]l(l, 1)+ + [{yv ~ yj)(xiVv - xvyi) + (yi - yv)(xvyj - Xjyv)]l(l, 0) + + [(xj - xv)(xiyv - xvyi) + (xv - Xi)(xvyj - Xjyv)]l(0,1) + c + l(xvyj - Xjyv)(xiyv - xvy{)--~[(yv ~ Vj)(yi ~ Uv) + (xj - xv){xl x uj ;)]]í(0,0) kde příspěvek pro levý a pravý trojúhelník získáme jako I\ = lv=\ a Ir = Iv=r. F8370 Moderní metody modelování ve fyzice jaro 2021 Sestavovací procedura ve 2D D. Hemzal, F. Miinz fiemzal ©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í: Kik + = 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ů. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Helmholtzova rovnice - okrajové podmínky Na závěr prozkoumejme možné okrajové podmínky pro Helmholtzovu rovnici. Dirichletova podmínka předepisuje amplitudu vybrané frekvenční komponenty v daném místě hranice. Homogenní Neumannova podmínka představuje plně odrazivou hranici. Poslední podmínka, kterou je vhodné připravit, je pak hranice volně vlnění propouštějící. Tuto podmínku lze však realizovat pouze přibližně, a to například takto: uvažujme rovinnou vlnu u = uo exp(ifc • x). Má-li taková vlna volně projít hranicí, musí být její gradient Vit = u0ik exp(ifc • x) = iku rovnoběžný s normálou k této hranici. Přenásobíme-li obě strany rovnice normálou n, Vit • n = un = ik • nu vede požadavek rovnobežnosti obou vyskytujících se vektorů na podmínku k ■ n = \k\\n\ a tedy no 27ri A což je podmínka smíšeného typu. Protože tato podmínka zajištuje, že vlnění bude procházet kolmo každou příslušnou částí okraje úlohy, je velmi důležité konstruovat tvar hranice pro numerickou simulaci v souladu s předpokládaným průběhem analytického řešení úlohy. Jinou možností je k hranici, kterou má vlna simulovaný objem opouštět, přidat několik vrstev prostředí, ve kterém gradientně budeme zvyšovat útlum prostředí. Při vhodném nastavení se téměř žádná energie nevrátí zpět do simulovaného objemu a vlna efektivně prostředí volně opustí (nezávisle na směru letu vůči hranici).