s ss s ss = ε = = s ss s (3) と表される s の要素における s s = κ = κ, =,, (4) jωε jω s は複素比誘電率に相当する物理量であり ここで PML 媒質定数を次のように定義する すなわち κξ をPML 媒質の等価比誘電率 ξ をPML 媒質の

Transcription

1 FDTD 解析法 (Matlab 版 2 次元 PML) プログラム解説 v 概要 FDTD 解析における吸収境界である完全整合層 (Perfectl Matched Laer, PML) の定式化とプログラミングを2 次元 TE 波について解説する PMLは異方性の損失をもつ仮想的な物質であり 侵入して来る電磁波を逃さず吸収する 通常の物質と接する界面でインピーダンスが整合しており 反射波が生じないという優れた特徴をもつ PMLを利用することにより FDTD 解析において問題となる解析領域端での反射波の影響を除くことが可能となり 解析精度を大きく改善し 様々な電磁波問題を解くことができる 2.2 次元 TE 波におけるPMLの定式化ここでは 2 次元 TE 波についてのマクスウェル方程式 E = εε r J t E = r (1) E = r t に対応するPMLの定式化を行う PML 媒質内では通常の導電率 = とし その代わりに以下の 異方性テンソルで表される損失項を導入する TE 波 H E H 伝搬方向 3 次元空間におけるPMLのためのマクスウェル方程式が 異方性 PMLのテンソル s を用いて次式で与えられていることを前提とする ( この原理を理解するためには Bereger の電磁界分離による PML 定式化等を参照されたい ) すなわち E Ampere の法則 : H = εε r s J (2a) Farada の法則 : E = r s (2b) ただし ε は真空の誘電率 ε r は通常媒質の比誘電率であり テンソル s は

2 s ss s ss = ε = = s ss s (3) と表される s の要素における s s = κ = κ, =,, (4) jωε jω s は複素比誘電率に相当する物理量であり ここで PML 媒質定数を次のように定義する すなわち κξ をPML 媒質の等価比誘電率 ξ をPML 媒質の電気的な等価導電率 ξ をPML 媒質の磁気 的な等価導電率とする このような媒質定数により電磁波の伝搬方向によって損失 ( 減衰 ) が異なる 状況を実現する ( 異方性損失 ) PML では周波数分散がゼロとなる条件 ( 伝搬波形がひずまない 無 ひずみ条件 ) ε = (5) が成立していることで 周波数に依存しない吸収媒質を実現する さらに PML 媒質の誘電体とし ての構成方程式を s D = εε E s D = εε E s D = r s r s εε r E s (6) および 磁性体としての構成方程式を s B = H s B = H s B = r s r s r H s (7)

3 と置く 以上の関係から 3 次元空間における定式化を行い 最終的にそれらの方程式から 2 次元 TE 波の方程式 (1) すなわち 電磁界成分 ( E, H, H ) に相当する方程式を抽出することとする Ampere の法則 (2a) を直角座標系において書き下す まず 角周波数 ω である電磁波について ss E s J ss jωεε r E J = s J ss E s sd J = jω sd J sd J kd s D J 1 = jω kd s D J ε kd s D J と書けるが この式を逆フーリエ変換によって時間に関する微分方程式に変換する すなわち 時間 微分演算子 はフーリエ変換により複素角周波数 jω となることから (8) 式は (8) kd D J 1 kd D J = t ε kd D J (9) と変換できる これは 電磁界成分の時間変化が記述された時間に依存する方程式である 次に PML の誘電率に関する構成方程式 (6) 式からは sd = εεse sd = εεse sd = εεse r r r (1) が得られ テンソル成分を用いて書き下し

4 ( jωκ ) D = ε ε ( jωκ ) E ( jωκ ) D ε ε ( jωκ ) E ( jωκ ) D = ε ε ( jωκ ) E r ε ε = r ε ε r ε ε (11) さらに逆フーリエ変換によって時間に関する微分方程式 ( κd) D = εε r ( κe) E ε ε ( κd) D = εε r ( κe) E ε ε ( κd) D = εε r ( κe) E ε ε (12) 同様に Farada の法則 (2b) についても (5) の関係に注意し E E ss H s E E ss = jω r H s E E ss H s sb = jω sb sb kb s B 1 = jω kb sb kb s B kb = jω s B kb 1 s B ε kb s B (13) と書け これを逆フーリエ変換し

5 E E kb B E E 1 kb = B t kb B E E kb B 1 = kb B ε kb B (14) PML の透磁率に関する構成方程式 (7) 式からは sb = sh sb = sh sb = sh r r r が得られ 対応する時間に依存する微分方程式は (12) 式を参考に ( κ B) B = r ( κh) H ( κ B) B = r ( κh) H ( κ B) B = r ( κh) H となる 以上の式より 2 次元 TE 波に関する方程式を抽出すると (9) 式第 2 式より D = ( kd) J ε (15) (16) (17) (12) 式第 2 式より ( κ ) εε ( κ ) D D = r E E ε ε (18) (14) 式の第 1 式および第 3 式で電磁界成分 ( E, D, H, B, H, B ) に関する項のみ考慮し

6 E B = ( kb ) ε (19) E B = ( kb ) ε (2) (16) 式の第 1 式および第 3 式と (5) の関係を用いて ( κ B) B = r ( κh) H ε ε (21) ( κ ) ( κ ) B B = r H H ε ε (22) 3.2 次元 TE 波における PML の差分式の導出 前節 (17) 式から (22) 式を差分化し FDTD 法による時間発展方程式を導出する 差分化には2 次元 TE 波で用いたものと同様の Yee 格子を用い それぞれ E と D H と B を同じ格子点にとる 平面内の (i) TE 波 i1 E i, 1, k i H i,, k i H ik,, E ik,, H ik,, E ik,, 1 i- H i,, k o k- k k k1 (k)

7 (17) 式を差分化し H H H H D D D D 1 1 ik,, ik,, i,, k i,, k ik,, ik,, ik,, ik,, = k J ik,, D D Dt ε 2 (23) となり これを D 1 について解くと D D H H H H J 1 2εk Dt 2εDt ik,, ik,, i,, k i,, k ik,, = ik,, ik,, 2εk D t 2εk Dt D D (18) 式を差分化し 1 D D D D E E E E εε r ε ε ik,, ik,, ik,, ik,, ik,, ik,, ik,, ik,, k = k Dt 2 Dt 2 (24) (25) となり これを E 1 について解き 2εk Dt 1 E = E Dt D Dt D ( ) ( 2εk ) ( 2εk ) t 1 1 ik,, ik,, ik,, ik,, 2εk D t εε r 2εk D (19) 式の差分式は (26) E E B B B B = ik,, 1 ik,, ik,, ik,, ik,, ik,, k t 2 (27) となり これを B について解き B 2 k t 2 E E t ik,, 1 ik,, ik,, = B ik,, 2k t 2k t (28) 同様に (2) 式の差分式は E E B B B B i, 1, k ik,, i,, k i,, k i,, k i,, k = k t 2 (29) となり これを B について解き

8 B 2 k Dt 2 Dt Ei, 1, k E ik,, = t t i,, k B i,, k 2k 2k (3) その他 (21) (22) の差分化も同様に行う 導出は各自で試されたい 4.PML における材料パラメータの決定 PML の厚みを d とし その領域内で電磁波が完全に吸収されるようにするため 誘電率と導電率 に相当するパラメータ κ, ( ただし =,, ) を次のように決定する ξ ξ TE 波 方向 : 図中の PM 方向 : 解析領域 方向 : PML 領域 ( 厚さ d) L 各領域において内側から 外側への 座標を ξ とし ( ) = d m ma, ma ( m 1)l{ R()} m 1 = opt = 2ηε d 15p ε r r ただし η は空間のインピーダンス R() はξ = の境界において達成したい反射率 および m = ( κma ) κ ( ) 1 1 d, κ ma = 1 2 のように決定して電磁波を減衰させ 内側の解析領域内では =, κ = 1 とすることにより通常の空間を実現する

9 5. 解析例 5.1 正方領域における 2 次元 TE 波の伝搬 Xma=1m (4 Yee cells, id=8) Poit source PML (4 laers, pml=8) Zma=1m (4 Yee cells, kd=8) ファイル名 :fdtd_2d_fujii_matlab_pml_v211.m 解析条件およびパラメータの値 : 解析領域の大きさ ma = 1 m ma = 1 m 周波数 freq = 1 GH 分割数 (Yee 格子の数 ) 方向 4 (id = 8) 方向 4 (kd = 8) PML 層数 = 4 (pml = 8) 繰り返し計算回数 ma = 1 回波源の大きさと位置 : 点波源 中央

10 5.2 矩形導波管 または平行平板導波路における TE 波の伝搬と遮断 Xma=.2 m (2 Yee cells, id=4) Poit source Metal (Electric Boudar) PML (4 laers, pml=8) Zma=1m (5 Yee cells, kd=1) ファイル名 :fdtd_2d_fujii_matlab_pml_v211_rectagularwaveguide.m 解析条件およびパラメータの値 : 解析領域の大きさ ma =.2 m ( 遮断周波数 =.75 GH) ma = 1 m 周波数 freq =.5 GH( 遮断周波数以下 ) および 1 GH( 伝搬 ) 分割数 (Yee 格子の数 ) 方向 2 (id = 4) 方向 5 (kd = 1) PML 層数 = 4 (pml = 8) ただし = および =ma の境界のみ繰り返し計算回数 ma = 25 回波源の大きさと位置 : 点波源 左端

11 5.3 2 次元誘電体スラブ導波路中のパルス状 TE 波の伝搬と導波路分散による波形歪み Xma=.2 m (2 Yee cells, id=4) Poit source 空気 (ε r = 1) 誘電体 (ε r = 4).5m PML (4 laers (cells), pml=8) Zma=1m (15 Yee cells, kd=3) ファイル名 :fdtd_2d_fujii_matlab_pml_v211_slabwaveguide_groupvelocit.m 解析条件およびパラメータの値 : 解析領域の大きさ ma =.2 m ma = 1 m 周波数 freq = 3 GH 分割数 (Yee 格子の数 ) 方向 2 (id = 4) 方向 15 (kd = 3) PML 層数 = 4 (pml = 8) 全周囲繰り返し計算回数 ma = 5 回波源の大きさ 位置 時間波形 : 点波源 左端 波束状パルス 正弦波が Raised-cosie 型包絡線 ( エンベロープ ) で変調されており 包絡線中に正弦波が5 波長分存在するパルスが伝搬するに従い 導波路の群速度分散によってパルス波形が広がる様子を解析している さらに 正弦波の位相速度に比べ群速度が遅いため 正弦波がパルス中で相対的に前に移動していく様子が見られる