F8370 Moderní metody modelování ve fyzice, jaro 2015 Kateřina Zuzaňáková
Rozptyl skalární vlny na neviditelné kouli 1   Formulace problému
Řešíme problém rozptylu skalární vlny na tzv. neviditelné kouli (disku) ve 2D. Neviditelná koule o jednotkovém poloměru je popsána indexem lomu
2
Na okraji koule, tj. v r = 1, navazuje index lomu spojitě na okolní prostředí - vakuum. Volíme jednotky, ve kterých je c = 1. Šíření vlnění je pak popsáno vlnovou rovnicí
A#(r,t) - n2(r)d2^(r,t) = 0. (2)
Hledáme stacionární řešení ve tvaru ^(r, t) = ip(r) exp(—lut). Dosazením do vlnové rovnice (2) dostaneme pro funkci ip(r) Helmholtzovu rovnici
(A + A2(r)) ^(r) = 0. (3)
Simulace bude spočívat v řešení rovnice (3) v kartézských souřadnicích na obdélníku x G [—a, a],
y e [-6,6].
Na okraji x = —a volíme Dirichletovu okrajovou podmínku
4>\x=-a = 1, (4)
která odpovídá tomu, že na kouli dopadá v kladném směru osy x rovinná vlna. Na ostatních okrajích volíme Robinovu podmínku
ifjn - ikifj = 0, (5)
kde k = 2tt/\ je vlnové číslo.
2 Diskretizace
Vyjdeme z Galerkinových podmínek
j :       Y]#] / N[j] (AN[i\ +uo2n2(r)N[i}) dQ = 0. (6) Jn
Integrací per partes dostaneme pro jtý uzel rovnici
j :    V#] / N\j]VN[i]-ňdl- / (ViV[j] • V N [i] - uo2n2(r)N[j\N[i}) dQ = 0. (7)
i       Jan i Jn
Koeficienty u označíme
'3i
lan
Jj{ = j   N[j]VN[í]-ndl (8) 1
Kji = / (ViV[j] • VN[i] - uo2n2(r)N[j\N[i}) díl. Jn
S tímto označením má rovnice (7) tvar
(9)
(10)
2.1   Diagonální prvky
Nejprve spočítáme diagonální prvky matice K. S použitím aproximace indexu lomu
\r) = ^n2[k]N[k]
n
(H)
mají tvar
K» = Y, í (VnTb1-VnT[j]-a;2Vn2[fcK[fcKL7]2 Jda
T3j Jt\ k J
Gradient tvarové funkce nT[j], která je nenulová na elementu T s vrcholy j, r, s, je
1
VnT[j] = (ys - yr, xr - xs) —. Integrál ze skalárního součinu gradientů v diagonálním členu (12) je tak roven
VnT[j] ■ VnT[j] díl
[{V s - Vr)2 + (X s ~ Xrf] díl [{Vs - Vr)2 + (Xs ~ Xrf]
2 lArľ
(12)
(13)
(14)
Nyní se podívejme, jak přispěje do Kjj druhý člen v kulaté závorce (12). Do integrálu přes element T přispějí nenulovou hodnotou pouze tři členy sumy přes k:
T
^2n2[k]nT[k]nT[j]2 díl
ker
T
cj2[n2[j] / nT[jY díl + n2[r] / nT[r]nT[j]2 díl + n2[s] / nT[s]nT[j]2 díl   . (15)
t
Protože
nT[r]nT[j]2 díl = T 60
je výraz (15) roven
2 ( „2\X\ |At'
A
~6Ô~
V        20 1 J 60 1 J
Diagonální člen Kjj má tedy tvar
T3j  í 1
2|At| (n2\j\+l-{n2\r]+n2\s])
20
™2[j] + 3 [n2[r] + n2[s\)
IA
20
(16)
(17)
(18)
(19)
2
2.2   Nediagonální prvky
Nyní spočteme nediagonální prvky matice K,
K^ = Y. í ( VMj] ■ Vnr[«] - uj2Y,Ak\nT[k}nT[3\nT\i} ) díl (20)
Sčítáme přes elementy T, které obsahují uzly j a i. Pro dané j a i budou takové elementy maximálně dva. Třetí uzel elementu označíme v. První člen v kulaté závorce (20) přispěje do
K jí hodnotou
^2 í VnT[j] • VnT[i] díl = ^      - yv) (yv - yj) + (x{ - xv) (xv - Xj)]    * (21)
Příspěvek druhého členu v kulaté závorce (20) bude následující. V integrálu přes T se opět projeví pouze tři členy ze sumy přes k:
J^w2 í ^2n2[k]nT[k]nT[j]nT[i]dQ =
T3j,i T keT
^£^(»'li]+«'M + ^). (22)
Celkem dostáváme pro nediagonální prvky vyjádření
Ki't = Yl \^Vi ~ ^ (yv ~ y^ + (Xi ~ x^ (Xv ~       2IA I ~ ^2^0^ (n2^ + ^ + ~~2~^) }
T3j,i ^ \    T\ \ J )
(23)
2.3   Okrajové příspěvky
Nakonec se podíváme na okrajový integrál
Jji= /   N[j]VN[i]-ňdl, (24) Jan
který je nenulový pro uzly j na části hranice s Robinovou okrajovou podmínkou
(VN[í] -n- ikN[í]) = 0. (25)
i
S využitím této podmínky přepíšeme okrajový integrál
Jjí= í N[j\VN[i] ■ n dl = ik í N\j]N[i]dl. (26)
Jan Jan Nediagonální prvky J jí jsou nenulové, pouze pokud j a i jsou sousední hraniční uzly. J jí pak má tvar
Jji = ik^f, (27) kde L ji = a/ (xj — Xi)2 + (yj — y,i)2 je vzdálenost uzlů j a i. Diagonální prvky J j j mají tvar
Jn = ifcl (Ly + Ljr) , (28)
kde / a r jsou hraniční uzly sousedící s uzlem j.
3
2.4   Výsledná soustava rovnic
Rovnice sestavujeme pro uzly, ve kterých neznáme funkční hodnotu. Index j v rovnici (10) tedy prochází pouze uzly uvnitř simulované oblasti a uzly na části hranice s Robinovou okrajovou podmínkou, zatímco sčítací index i prochází všechny uzly. Sčítance, ve kterých je i hraniční uzel s Dirichletovou okrajovou podmínkou přesuneme na pravou stranu rovnice a dostaneme
J ji na pravé straně je nenulové pouze pro uzel i v levém horním, resp. dolním rohu simulované oblasti, tj. (—a,±6), a uzel j, který je jeho sousedním uzlem na úsečce y = ±6.
Síť bodů byla vytvořena tak, aby ve středu simulované oblasti, tj. v místě, kde se nachází čočka, byla hustější než u okrajů. Byl vytvořen program ve fortranu, který načte souřadnice uzlů a seznam elementů, naplní matici soustavy, spočítá vektor pravé strany a uloží je do textových souborů. K řešení získané soustavy rovnic pak byl použit software SuperLU.
Grafy funkce ip pro tři různé frekvence uj jsou vykresleny v obrázku 1. Vlevo je vždy umístěn graf funkce ip spočítané metodou konečných prvků. Vpravo je pro srovnání uveden odpovídající graf vypočtený semi-analytickou metodou (separací proměnných v Helmholtzově rovnici, numerickým řešením radiální rovnice uvnitř čočky a navázáním tohoto řešení na řešení vně čočky) viz [1]. Černá kružnice znázorňuje okraj čočky.
Vidíme, že funkce získané oběma způsoby si docela dobře odpovídají. Simulace se výrazně liší pouze v blízkosti horního a dolního okraje oblasti, což ale není překvapivé vzhledem k tomu, že na těchto okrajích byla při výpočtu metodou konečných prvků zvolena Robinova okrajová podmínka, která neodpovídá očekávanému šíření vlny u těchto okrajů.
Reference
[1] Tyc T., Chen H., Danner A., Xu Y., Invisible lenses with positive isotropic refractive index, Physical Review A 90, 053829 (2014).
(29)
3 Výsledky
4
a) u = 4
Obrázek 1: Funkce ip vypočtená metodou konečných prvků (vlevo) a semi-analytickou metodou (vpravo) pro různé frekvence. Kladné funkční hodnoty jsou vykresleny červeně, záporné modře a sytost barvy odpovídá velikosti funkční hodnoty.
5