Givoli-Neta-Patlashenkoの高次吸収境界条件を用いた時間領域FEMの平行平板誘電体スラブ導波路の散乱係数計算
Givoli,Neta,Patlashenkoが提出した高次吸収境界条件は、導波管のスカラ波の吸収には効果がありました。
今回は境界上の媒質が変化する場合について試してみました。文献では媒質は一様として導出されているためそのまま使用できません。先ずは媒質が変化する場合を定式化してみました。
元の文献は、こちら。
Dan Givoli , Beny Neta & Igor Patlashenko
"Finite element analysis of time-dependent semi-infinite wave-guides with high-order boundary treatment"
http://www.math.nps.navy.mil/~bneta/higdon_fe_pap.pdf
February 27, 2003
International Journal for Numerical Methods in Engineering
Volume 58, Issue 13, p.1955- 1983,
7 December 2003
境界を線要素で分割したとき、要素内では媒質定数は一定であるとすると、各々の要素内では、次のスカラヘルムホルツ方程式が成立します。
q d2 u / dt2 - c0^2 p ∇^2 u = 0 (2)
ここで、TEモード(TMzモード)のとき
u = Ez, p = 1 / μr , q = εr
です。
Higdon ABCは文献と同じく
(d / dx + (1 / c1) d / dt) u = Φ1
(d / dx + (1 / c2) d / dt) Φ1 = Φ2
(d / dx + (1 / cJ) d / dt) ΦJ-1 = 0
ここで
Φ0 = u
ΦJ = 0
とすると、
(d / dx + (1 / cj) d / dt) Φj - 1 = Φj j = [1, J] (17)
さて、Φjは(2)を満たすものとすると
p d2 Φj / dx2 + p d2 Φj / dy2 - (1 / c0^2) q d2 Φj / dt2 = 0 (19)
ここで、
p d2 Φj / dx2 = p (d / dx - (1 / cj + 1) d / dt) (d / dx + (1 / cj + 1) d / dt)
+ (p / cj + 1-2) d2 Φj / dt2 (20)
より
p (d / dx - (1 / cj) d / dt) (d / dx + (1 / cj) d / dt) Φj - 1
+ (p / cj + 1^2 - q / c0^2) d2 Φj - 1 / dt2
+ p d2 Φj - 1 / dy2 = 0 (21)
(17)より
p (d / dx - (1 / cj) d / dt) Φj + (p / cj^2 - q / c0^2) d2 Φj - 1 / dt2
+ p d2 Φj - 1 / dy2 = 0 (22)
また、(17)のj + 1に関する式より
p (d / dx + (1 / cj+1) d / dt) Φj = p Φj + 1 (23)
(23) - (22)より
p ( 1 / cj+1 + 1 / cj) d Φj / dt =
p Φj + 1 + (p / cj^2 - q / c0^2) d2 Φj -1 / dt2 + p d2 Φj - 1 / dy2 (24)
従って、境界上が不均一な媒質からなる場合のGivoli-Neta-PatlashenkoのABCは。
β0 d2 u / dt2 + p d u / dx = p Φ1 (25)
βj d Φj / dt - αj d2 Φj - 1 / dt2 - p d2 Φj - 1 / dy2 = p Φj + 1 (26)
j = [ 1, (J - 1) ]
αj = p / cj^2 - q / c0^2
β0 = p / c1
βj = p (1 / cj + 1 / cj + 1)
となります。
吸収境界条件の速度cj は境界上でy方向には一定とした場合、cjは導波路のモードの位相速度を設定することができます。
cjがy方向にも変化する設定も可能だと思われますが、まだ試していないので何とも言えないです。
平行平板誘電体スラブ導波路グレーティングの散乱係数の計算
比誘電率ε core = 2.4025、導波路幅aの誘電体導波路()が導波路幅W = 30 a平行平板導波路に装荷されているモデルを考えます。
誘電体導波路に８個の比誘電率ε grating = 2.1025、長さaの誘電体が間隔aで埋め込まれているとします。
これは、周波数領域FEM（２次三角形要素)で計算した結果です。これから時間領域FEMで計算するモデルでは分割数は２倍とり、1次三角形要素で分割しました。また、左右の一様誘電体導波路部分を長くとりました。
比誘電率分布
入射波としては、変調ガウシアンパルスを用いました。
Ez(t) = E0(y) cos (ωc ( t - T0)) exp ( - (t - T0)^2 / (2 Tp^2))
T0 = 0.48 (1.0 / fc) ( 5 / 2)
Tp = T0 / (2 √(2 ln(2) )
fc：搬送波の周波数
ωc = 2Πfc
E0(y) 搬送波周波数のときの誘電体装荷平行平板導波路の基本モード分布
ここで、fcに対応する規格化周波数は2W / λ = 12に選びました。T0の0.48や搬送波の周波数は不連続のない誘電体導波路の電界時間変化を計算して反射波が少なくなるように選びました。
5次ABCを用いたときの計算結果は次の通りです。
ABCの次数を変化させたときの結果は次の通りです。
1次ABC
2次ABC
3次ABC
4次ABC
5次ABC
6次ABC
良好な結果が得られていると思います。
【追記 2013-06-22】
DelFEM4Netで上記サンプルを実装しました。よろしければご覧ください。
DelFEM4Net_SampleAppのwg2d_td_higher_order_abcに含まれています。スペースキーで問題が切り替わります。計算終了後gキーでグラフのみの画面を表示できます。
DelFEM4Net https://code.google.com/p/delfem4net/
DelFEM4Net_SampleApp(2013-06-22) https://code.google.com/p/delfem4net/downloads/detail?name=DelFEM4Net_SampleApp_20130622.zip
依存DLL
※ILNumerics(32Bit)のdllです。正式には本家よりインストールしてください。
ILNumerics http://ilnumerics.net/
本家を見たら以前見たときとかなり変わってました。なんかCommunity Edition(無料)とProduct Editionが明記されていますね。Community Editionの場合、GPL v.3に準拠する必要があるようです。
