F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Metoda konečných diferencí (Finite Differences, FD) Podobně jako jiné numerické metody, metoda konečných diferencí (MKD, někdy též metoda sítí) přibližně převádí obtížně uchopitelné diferenciální rovnice na rovnice algebraické, v tomto případě prostřednictvím jednoduché diskretizace studované oblasti a aproximace derivací diferenčními vztahy. Jsou-li výchozí diferenicální rovnice lineární, obdržíme soustavu lineárních algebraických rovnic a obtížnost teoretickou tak nejčastěji nahrazujeme obtížností numerickou (pro věrné vystižení podložních funkcí musí být diskretizace jemná a získaná soustava rovnic má zpravidla značnou dimenzi). značení: vektory a, b jsou sloupcové, řádkové jsou aT, bT složky vektoru ai matice A, B, dimenze dimA = mxn, m je počet řádků, n je počet sloupců prvky matice Aij, první index je řádkový, druhý sloupcový x[n], y(x[n]) = y[n] F8370 Moderní metody modelování ve fyzice jaro 2021 Diskretizace úlohy D. Hemzal, F. Münz hemzal@physics.muni.cz Geometrie diskretizační sítě je prakticky výhradně dána zvoleným typem souřadnic. Ve 2D tak vhodnými transformacemi mohou základními elementy kromě čtverce být také obdélníky, rovnoběžníky, nebo například výseče mezikruží v případě polárních souřadnic. Otázka řešitelnosti vzniklé algebraické soustavy je obecně dosti komplikovaná, v jednoduchém lineárním případě stačí spočítat neznámé a nezávislé rovnice. Je-li soustava nedourčená, máme volnost definovat další doplňující podmínky, kterými ji stabilizujeme (hovoří se o numerickém ukotvení úlohy) - například při hledání vlastních vektorů není z principiálních důvodů dána jejich velikost, kterou si tedy můžeme dodatečně zvolit. Je-li soustava přeurčená, nezbývá, než řešit ji přibližně. F8370 Moderní metody modelování ve fyzice jaro 2021 Řešení přeurčené lineární soustavy D. Hemzal, F. Miinz hemzal@physics.muni.cz Uvažujme přeurčenou soustavu lineárních rovnic, Ax = b (dimA = mxn, m > n). Takovou soustavu nelze obecně řešit, můžeme však požadovat její co nejlepší splnění ve smyslu nejmenších čtverců. Zavedeme-li vektor rezidua jako r = Ax — b, můžeme minimalizovat jeho normu rTr, (Ax — b)T(Ax — b) —y min. Roznásobením a derivací např. podle xT získáme podmínku ATAx — ATb = 0 , čili jsme přešli k soustavě Áx = b se čtvercovou pozitivně deflnitní maticí A = A A (dimA = nxn), kde b = A b. Jinými slovy, ve smyslu nejmenších čtverců je přeurčená soustava řešitelná vždy. F8370 Moderní metody modelování ve fyzice jaro 2021 Cauchyova úloha pro Laplasián v metodě konečných diferencí Nejčastěji užívaná centrální diference pro derivaci druhého řádu v ID má tvar u[i - 1] - 2u[i] + u[i + 1] D. Hemzal, F. Múnz hemzal@physics.muni.cz u" \i + 0(e). Výhodou kartézských souřadnic je oddělenost derivací v jednotlivých směrech: dx2 dy2 takže například ve 2D u[i-l,j]-2u[i,j]+u[i + lj] u[i, j - 1] - 2u[i, j] +u[i, j + 1] Au{i,j] =-j2-+-j2-+ O(e), což dává známé pravidlo anulování Laplasiánu při průměru sousedních hodnot u[i - 1, j] +u[i, j - 1]-4w[i, j] +u[i + +u[i,j + 1] Aíí[í,j] =-^-+ O(e). £2Au : i 3) © 0 © F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz K plnému určení dvoubodové aproximace gradientu tedy potřebujeme jeden uzel na každou stranu od vyhodnocovaného, do problému se tedy dostaneme u okrajů zkoumané oblasti. Formálně Cauchyův problém tedy vyžaduje doplnění kompatibilních okrajových podmínek, jejich diskretizovanou podobu musíme zavést i my. Budeme se zabývat dvěma typy okrajových podmínek: Dirichletova podmínka specifikuje hodnoty proměnné v některých bodech díl Jelikož je taková hodnota fixní ve smyslu neovlivnitelnosti vnitřním děním systému, lze ji například v termodynamice považovat za kontakt s termostatem (neznámá veličina je teplota systému). Neumannova podmínka specifikuje hodnoty normálové derivace proměnné v některých bodech dí} Tato hodnota vlastně specifikuje tok neznámé veličiny hranicí, takže nejběžnější podmínka un = Vu • n = 0 vlastně představuje podmínku izolace okraje (například tepelné). Dirichletova podmínka se zavádí přímo, Neumannova prostřednictvím metody zrcadlení. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Metodu zrcadlení je také třeba využít, pokud jsme provedli symetrickou redukci íl - jí vznikají uměle nové hranice, které se doplňují podmínkou Neumannova typu - na obou stranách hranice symetrie musí mít uzly (právě ze symetrie) stejné hodnoty. Na každé části dí} musí být předepsán jen jeden typ okrajové podmínky, na množině míry nula ovšem mohou být předepsány oba. Podobně, při styku dvou hraničních ploch s různými Dirichletovými podmínkami zpravdidla na jejich průsečíku zavádíme Dirichletovu podmínku o velikosti průměru výše jmenovaných. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Múnz hemzalOphysics.muni.cz Příklad - zavedení okrajových podmínek v ID Poissonově rovnici Uvažujme rovnici ^(#-l]-2#] + # + l]) jejího operátoru na ekvidistantní síti s N + 1 uzly; potom £[i] = a;[i + 1] rovnic vede na soustavu x[i] = a/N. Prvoplánové sestavení í-2 1 0 0 ... 0\ 1 -2 1 0 0 0 0 1-21 0 0 1 -2/ / 0[O] \ Í>[N-1 / /[O] \ V fW ) Tato soustava má problém v okrajových uzlech - aby diskretizace platila beze zbytku, musí být [0] = 0[iV + l] = 0. Jedná se vlastně o podmínku Dirichletova typu a vzhledem ke způsobu, kterým podmínka vznikla, se označuje za implicitní. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz Pokud bychom chtěli předepsat explicitní Dirichletovu podmínku, například [0] = T, museli bychom rovnici pro nultý uzel vynechat (hodnota v tomto uzlu je známa), a pro první uzel bychom psali což by vedlo na soustavu 1 ^(T-20[l] + 0[2]) = /[l í-2 1 0 0... 0 \ 1 -2 1 0 ... 0 0 ... 0 1 -2 1 \0 ... 0 0 1 -2/ 4>[2] j>[N-l //M - m2\ W-i v nm ) Podobně, pokud bychom navíc předepsali Neumannovu podmínku [n\ = enm - 2^, D. Hemzal, F. Múnz hemzal@physics.muni.cz a celkově soustavu 1 (-2 1 o o 1-210 o 0 ... 0 1 -2 1 \0 ... 0 0 2 -2/ [2] j>[N-l ífm-T/e\ /[2] f[N-l \f[N\-2b/ej přičemž pro 6 = 0 rozeznáváme v hlavní matici příspěvek od zrcadlení okraje. Stojí za povšimnutí, že uvedený postup automaticky monitoruje konzistenci okrajové podmínky s rovnicí. F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzal@physics.muni.cz kde /?[kg.m-3] je hustota vzorku, cř>[J.kg_1K_1] jeho tepelná kapacita a A koeficient jeho tepelné vodivosti; všechny veličiny mohou být časově i prostorově závislé. Z jednotkové analýzy levé strany vyplývá, že q je tepelný výkon ve Wattech, absorbovaný na jednotku objemu. Ze stejného důvodu musí mít tepelná vodivost jednotku Wm_1K_1; čím vyšší tepelná vodivost, tím více tepla proudí vzorkem, pokud jeho dva vzdálené konce udržujeme na odlišných teplotách. Za předpokladu ustáleného stavu se rovnice vedení tepla redukje na Na zvolené geometrii tuto rovnici doplníme fixními okrajovými podmínkami, pak můžeme úlohu rozřešit. V • (AVT) + q = 0, v homogenním vzorku, nepodrobenému dodávce tepla potom až na AT = 0. F8370 Moderní metody modelování ve fyzice jaro 2021 Ti, T2 konstantní • • • • • _ 9 1 1cu 11_ 12_ 11 3^ 14^ 15 16^ 17^ 16. 18 O o o D. Hemzal, F. Miinz hemzalOphysics.muni.cz 3 : 2T[1] + T[4] + T[8] - 4T[3] = 0 4 : T[3] + T[5] + T[9] + T[l] - 4T[4] = 0 5 : T[4] + T[6] + T[10] + T[l] - 4T[5] = 0 6 : T[5] + T[7] + T[ll] + T[l] - 4T[6] = 0 7 : 2T[6] + T[12] + T[l] - 4T[7] = 0 21 : T[l] + T[22] + 2T[18] - 4T[21] = 0 22 : T[21] + T[23] + 2T[19] - 4T[22] = 0 23 : T[22] + T[2] + 2T[20] - 4T[23] = 0 F8370 Moderní metody modelování ve fyzice jaro 2021 D. Hemzal, F. Miinz hemzalOphysics.muni.cz Solvery pro MKD MatrixMarket formát - pro zápis řídkých matic matrixFile rhsFile 22 22 62 1 3 1.0 0.0 1 1 -2.0 0.0 2 4 1.0 0.0 2 2 -2.0 0.0 3 1 1 10.0 0.0 2 1 10.0 0.0 superLU (zlinsol_print_new) matrixFile, rhsFile —> resultFile F8370 Moderní metody modelování ve fyzice jaro 2021 Shrňme si hlavní výhody a nevýhody MKD: výhody: přehlednost metody není potřeba speciálních generátorů dělení přesnost se dá regulovat volbou derivačních aproximací (/c-bodové) nevýhody: složité tvary vyžadují jemné dělení -velká náročnost na úložiště Obecně překonáno metodou konečných prvků (FEM), v konkrétních aplikacích ale použití FD s proměnným krokem může být velmi efektivní. I při použití FEM, řešení nestacionárních úloh používá v časové ose FD (časové projekce íl nemívají složité tvary) - vzniká FDTD. Každopádně je ke snížení počtu uzlů potřeba vždy maximálně využít symetrie úlohy. D. Hemzal, F. Múnz hemzal@physics.muni.cz