修 士 論 文 の 和 文 要 旨
研究科・専攻 大学院 情報理工学研究科 情報・ネットワーク工学専攻 博士前期課程
氏 名 星野 真伸 学籍番号 1931131
論 文 題 目 CIP-BS 法の低数値分散性を利用した大規模な散乱問題解析法の提案
要 旨
本研究ではConstrained Interpolation Profile-Basis Set (CIP-BS)法の粗い離散化であって も時間ステップを減少させることで位相誤差を低減可能という性質を利用し、波長に対し広大 な解析領域を考慮する問題への適用とその評価を行う。 まずCIP-BS 法の性質がどの程度の粗い離散化まで効果を発揮するのかを明らかにするため 分散関係式を導出し数値分散解析を行う。その結果、波長あたりのセル数であるサンプリング 密度一定の下で時間ステップを減少させたとき、数値位相速度は時間ステップの2 乗に比例す る傾きで厳密な速度に接近し、ある値より時間ステップが小さくなるとサンプリング密度の −5.5乗に比例する値に収束することが明らかになった。つまり、時間ステップやサンプリング 密度を適切な組み合わせで設定しなければ精度は上がらず計算時間だけが増大してしまう。 次に解析領域内に散乱体が分布する状況に対応するため、不均一な離散化をCIP-BS 法に適 用し数値解析を行いその誤差評価から散乱体モデリングに関するパラメータについて考察す る。解析状況として無限長円筒誘電体に対する散乱問題を考える。散乱場の振幅誤差は誘電体 直径あたりのセル数の−1.4乗に比例し、入射波長が変化しても傾きのオーダーは一定であった が、位相誤差は入射波長の差異に対し異なる分布を見せた。この結果から散乱体モデリングの ためのセル数による影響は振幅に生じやすいと推測する。 最後に波長に対し大規模な解析領域を有する散乱問題として3 × 3アレー状に配置した円筒
誘電体に対する散乱問題の数値解析を行い Finite Difference Time Domain (FDTD)法と比較
することで計算効率を評価する。同程度の位相精度で比較したときCIP-BS 法の更新式を解く
ために共役勾配法を用いた場合、高精度な領域ではFDTD 法に比べ省メモリであった。一方、
計算時間を比較したとき、CIP-BS 法に LU 分解を用いても FDTD 法が高速であったが高精度
化に伴う計算時間の増加割合は CIP-BS 法の方が小さいため、CIP-BS 法はより広大な解析領
令和
2
年度 修士論文
CIP-BS
法の低数値分散性を利用した
大規模な散乱問題解析法の提案
学籍番号
1931131
星野真伸
電気通信大学大学院
情報理工学研究科 情報・ネットワーク工学専攻
電子情報学プログラム
主任指導教員 安藤 芳晃 准教授
指導教員
野村 英之 教授
提出日
令和
3
年
1
月
25
日
概要
本研究では CIP-BS 法の粗い離散化でも時間ステップを減少することで位相誤差を低減可能であると いう性質に着目し、波長に対し広大な解析領域を考慮する問題への適用とその評価を行う。 まず CIP-BS 法の特性がどの程度の粗い離散化でも効果を発揮するのかを調査するため数値分散解析 を行う。その結果、サンプリング密度 νx一定のもとで時間ステップ ∆t を減少させていくと、νxが大き いときは νxに関係無く数値位相速度は O(∆t2)で厳密値に接近していくため、パラメータ決定に際し νx 大きいときに所望の精度を得ることができる ∆t を選択し、その中で νxを最小になる様に選択すべきで あると提案する。 CIP-BS法に不均一離散化を適用し、そのモデリング精度の評価のため T Mzモードにおける無限長円 筒散乱体の散乱問題を解析し散乱場の誤差を調査する。結果、位相誤差も詳細な離散化により誤差の低 減が見られるが直径あたりのセル数 Nmが 12 より大きいとき Nmの増加と共に誤差の低減割合は減少 する。即ち振幅誤差にモデリングの影響が強く現れるため、モデリングのためのセル数は振幅の精度が 所望の値となる様にすべきと考える。 最後に波長に対して広大な解析領域を持つ散乱問題として 3 × 3 のアレー状に整列した誘電体に対す る散乱問題の数値解析を行う。CIP-BS 法と FDTD 法を用いて解析を行い、同じ位相精度で比較した結 果、CG 法を用いた CIP-BS 法が省メモリであり、LU 分解を用いた CIP-BS 法が最もメモリを要した。 計算時間においては FDTD 法が高速であったが、高精度化に伴う計算時間の増加割合は CIP-BS 法が少 なかったためより広大な解析領域に適すると考えられる。目 次
第 1 章 序論 3 1.1 背景 . . . 3 1.2 目的 . . . 3 第 2 章 CIP-BS 法 4 2.1 1次元 CIP-BS 法の基底関数 . . . 4 2.2 基底関数の多次元化 . . . 6 2.3 CIP-BS法の計算スキーム . . . 6 第 3 章 数値分散解析 9 3.1 数値分散関係式の導出 . . . 9 3.2 解析結果 . . . 11 3.2.1 サンプリング密度への依存性 . . . 11 3.2.2 高次化の影響 . . . 13 3.2.3 時間ステップと数値位相速度の関係 . . . 15 第 4 章 散乱体のモデリングに関するパラメータの検討 17 4.1 CIP-BS法基底関数の拘束条件を満たす不均一離散化 . . . 17 4.2 数値解析 . . . 18 4.2.1 問題設定及び数値解析に用いるパラメータ . . . 18 4.2.2 計算結果 . . . 18 第 5 章 大規模散乱シミュレーション 23 5.1 数値解析 . . . 23 5.1.1 問題設定及び数値解析に用いるパラメータ . . . 23 5.1.2 計算結果 . . . 24 第 6 章 結論 29 付 録 A 基底関数同士の内積 32第
1
章 序論
1.1
背景
Constrained Interpolation Profile(CIP)法は流体力学の分野において移流方程式の数値解法として開 発された [1]。CIP 法は解だけでなくその導関数まで計算し高精度化を図る手法であり、双曲線型の偏微 分方程式に適用可能であるという汎用性から電磁気分野にも応用され成果を挙げている [2]。しかし電磁 界解析において主に扱われる多次元問題に適用する際に方向分離を必要とし、本来存在しないはずの高 次微分項を計算しなければならないため厳密に方程式を扱うことが出来ない。そこで CIP 法の補間関数 の拘束条件を満たす関数を基底関数とし、多次元問題に直接適用可能な CIP - Basis Set(CIP-BS) 法が 開発された [3] [4]。CIP-BS 法は比較的新しい手法ではあるものの精力的な研究の成果により電磁界解 析へ適用され、更に個別技術の開発も進み実用に足る手法となった [5][6]。CIP-BS 法は一つの離散点に 複数の未知係数を割当て、またそれらを係数とする直交していない基底関数を用いて Galerkin 法を適用 するためバンド幅の広い疎行列を解く必要があり、計算量は用いる連立方程式の数値解法に強く依存す る。特に反復計算において高速で正確な計算が可能な LU 分解はバンド幅による影響を受けやすく [7]、 CIP-BS法を適用する問題は吟味しなければならない。そこで CIP-BS 法の特性である他の手法に比べ 粗い離散化でも時間ステップを小さくすることで位相誤差を低減できるという点に着目する [6]。一般に 位相誤差低減には波長あたりのセル数であるサンプリング密度を増加させる必要があり、解析領域が波 長に比べ広大である問題ではサンプリング密度増加に伴う計算資源の増大が障害となる。CIP-BS 法の 特性を鑑みれば、そのような問題を粗い離散化でも高精度に計算を行うことが可能であり、効率的な解 析が期待できる。また解析領域中にグリッドに沿わない形状の散乱体が存在するときモデリングのため に詳細な離散化を必要とするため、CIP-BS 法の不均一離散化についても考察する。
1.2
目的
本研究では波長に比べ解析領域が広大である状況に対し効率的な数値解析手法を提案する。まず CIP-BS法の粗い離散化でも時間ステップを減少させることで位相誤差を低減できるという特性がどの程度の 粗さに耐えうるのかを調査するため数値分散解析を行う。すにで一般的に用いられている電磁界解析手 法である Finite Difference Time Domain(FDTD) 法では更新式を単一周波数の平面波について解くこと で位相速度などの諸性質が明らかとなっている [8]。そこで CIP-BS 法においても同様に数値分散関係式 を導出して解析を行い、その結果から所望の精度が得られ、効率の良いパラメータの決定法を提案する。 次に効率的な計算資源利用のため不均一離散化を検討する。CIP-BS 法基底関数は不均一離散化境界上 においても定義可能であることを利用し CIP-BS 法による不均一離散化を提案する。散乱問題のシミュ レーションを行いその数値解析結果からグリッドに沿わない形状の散乱体を十分なモデリングに必要な セル数を調査する。 最後に波長に対し広大な解析領域を考える必要のある問題に対し CIP-BS 法を適用し計算効率を調査 する。CIP-BS 法はバンド幅が広い疎行列を解く必要があるため、連立一次方程式の数値解法に計算効 率が比較的強く依存する。そこで直接法である LU 分解と反復法である CG 法のそれぞれを用いたとき の経過時間と使用メモリを調査する。また FDTD 法と比較し CIP-BS 法の計算効率を評価する。第
2
章
CIP-BS
法
本章では CIP-BS 法の基底関数が満たすべき性質と導出過程を述べ、導出した基底関数を用いた更新 式の定式化について述べる。以下では電磁場の 0 次導関数が登場するがこれは電磁場の値そのものを意 味する。2.1
1
次元
CIP-BS
法の基底関数
x
( )
F
x
iF
i
(0)i
F
(1)i+1
F
(0)i+1
F
(1)x
i+1x
図 2.1: K=1 の 1 次元 CIP 法補間関数の拘束条件の模式図 CIP-BS法では電磁場とその 1 次から K 次導関数までを計算するため、電磁場を展開する際それら導 関数も用いる。1 次元の解析領域を任意の間隔で離散化し、i 番目のノードの座標を xiとする。CIP-BS 法における電磁場 F (x) の展開式は、x = xiの F (x) の p 次導関数の値である Fi(p)を係数とする基底関 数 φ(p) i (x)を用いて次式で表される。 F(x) = Nx X i=0 K X p=0 Fi(p)φ(p)i (x) (2.1.1) ここで φ(p)i は基底関数である。この展開式と CIP 法補間関数が一致する様に φ (p) i は決定される。CIP 法補間関数は xi≤ x ≤ xi+1で図 2.1 に示す様な拘束条件を満たし、それは次式で表される。 ∂dF ∂xd(x) = Fi(d) (x = xi) Fi+1(d) (x = xi+1) (2.1.2)ここで φ(p)i の台を隣接する 1 セルに制限すると、xi≤ x ≤ xi+1で非零の値を持つのは φ(p)i と φ(p)i+1であ る。式 (2.1.2) に式 (2.1.1) を代入することで φ(p) i と φ (p) i+1の拘束条件が得られるが φ (p) i+1の条件は i = i −1 とすることで xi−1 ≤ x ≤ xiでの φ(p)i の条件とみなすことができる。したがって φ (p) i の満たすべき条件 は次式で表される。 ∂dφ(p) i ∂xd (xj) = 1 (p = d, j = i) 0 otherwise (2.1.3)
xi≤ x ≤ xi+1において φ(p)i が満たすべきは xj= xi, xi+1の式 (2.1.3) であり、これは 2(K + 1) 元連立 方程式である。つまり φ(p)+ が 2K + 1 次多項式ならば一意に決定できるので φ (p) i = P2K+1 n=0 an(x− xi)n とする。xj= xiの式 (2.1.3) より K 次以下の項は p 次の項だけであるため φ(p)i = 1 p!(x− xi) p+ 2K+1 X n=K+1 an(x− xi)n (2.1.4) となる。残った xj= xi+1の式 (2.1.3) より ∆x = xi+1− xiと置くと i ≥ K + 1 の aiに関する連立方程 式は次式となる。 1 1 · · · 1 (2k + 1) 2k · · · (K + 1) .. . ... 2K+1PK 2KPK · · · K+1PK ∆x2K+1−pa2K+1 ∆x2K−pa2K .. . ∆xK+1−paK+1 = −p!1 −p p! .. . −1 0 .. . 0 (2.1.5)
これを解くことで xi≤ x ≤ xi+1での φ(p)i が求まる。同様に xj= xi, xi−1の式 (2.1.3) より xi−1≤ x ≤ xi での φ(p)
i が求めるが、xj = xiの式 (2.1.3) より φ (p)
i が式 (2.1.4) で表されることに変わりは無く、xj= xi−1 の式 (2.1.3) は xj = xi+1の式 (2.1.3) が xi+1が xi−1 に置き換わっただけであるため、∆x が負となっ てしまうが ∆x = xi−1 − xi とすれば xi−1 ≤ x ≤ xiでの φ (p) i の係数 anを決定する連立方程式は式 (2.1.5)と全く同じ式に帰着する。以下に本研究で主に用いる K = 1 の基底関数 φ(0) i 、φ (1) i を示す。こ こで ∆x+= xi+1− xi、∆x− = xi− xi−1である。図 2.2、2.3 はそれぞれ φ
(0) i 、φ (1) i の概形である。各 ノード点で式 (2.1.3) の拘束条件を満たしていることが確認できる。 φ(0)i (x) = 2x−xi ∆x+ 3 − 3x−xi ∆x+ 2 + 1 (x≥ xi) −2x−xi ∆x− 3 − 3x−xi ∆x− 2 + 1 (x < xi) (2.1.6) φ(1)i (x) = ∆x+ x−xi ∆x+ 3 − 2x−xi ∆x+ 2 +x−xi ∆x+ (x≥ xi) ∆x− x−xi ∆x− 3 + 2x−xi ∆x− 2 +x−xi ∆x− (x < xi) (2.1.7) 0 1 xi - 1 xi xi + 1 図 2.2: 基底関数 φ(0)i の概形 0 xi - 1 xi xi + 1 図 2.3: 基底関数 φ(1)i の概形
2.2
基底関数の多次元化
一般的に電磁界解析では多次元問題を取り扱うが、CIP-BS 法を多次元問題に適用するためには基底関 数を多次元化しなくてはならない。そこで 1 次元と同様に 2 次元の電磁場 F (x, y) の展開式から基底関 数が満すべき要件を考え、基底関数を導出する。xy 平面上の解析領域を x 方向に Nxセル、y 方向に Ny セルに任意の間隔で離散化し x 方向に ix番目、y 方向に iy番目のノードの座標を (xix, yiy)とする。電 磁場 F (x, y) はノード (xix, yiy)上の x について px次、y について py次の F (x, y) の導関数の値 F (px,py) ix,iy を係数に持つ二次元の基底関数 Φ(px,py) ix,iy (x, y)を用いて次式で展開される。 F(x, y) = Nx X ix=0 Ny X iy=0 K X px=0 K X py=0 F(px,py) ix,iy Φ (px,py) ix,iy (x, y) (2.2.1) 展開式を x について px階、y について py階微分した導関数は (x, y) = (xix, xiy)で ∂px+pyF ∂xpx∂ypy(xix, yiy) = F (px,py) ix,iy (2.2.2) とならなければならないため基底関数 Φ(px,py) ix,iy は次の条件を満たす。 ∂dx+dyΦ(px,py) ix,iy ∂xdx∂ydy (xjx, yjy) = 1 ( dx= px, dy= py, jx= ix, jy = iy ) 0 otherwise (2.2.3) ここで Φ(px,py) ix,iy (x, y) = f (x)g(y)と因数分解可能であるすると式 (2.2.3) は f(x)、g(y) のそれぞれの条件 として以下の様に記述できる。 ∂dxf ∂xdx(xjx) = 1 ( dx= px, jx= ix) 0 otherwise (2.2.4) ∂dyg ∂ydy(yjy) = 1 ( dy= py, jy= iy ) 0 otherwise (2.2.5) 式 (2.2.4)、(2.2.5) は 1 次元の基底関数が満たすべき条件である式 (2.1.3) と一致するため、これらの条 件から導出される関数は 1 次元 CIP-BS 法の基底関数である。即ち 2 次元の基底関数 Φ(px,py) ix,iy は 1 次元 の基底関数 φ(p)i を用いて次式で定義される。 Φ(px,py) ix,iy (x, y) = φ (px) ix (x)φ (py) iy (y) (2.2.6)2.3
CIP-BS
法の計算スキーム
本節では電磁気学の支配方程式である Maxwell の方程式を CIP-BS 法によって解析する方法を述べる。 Maxwellの方程式の内、Faraday の式と Ampere の式は以下の 2 式である。∇ × E = −µ∂H∂t (2.3.1) ∇ × H = ε∂E∂t (2.3.2) ここで E は電場、H は磁場、ε は誘電率、µ は透磁率である。本研究では主に 2 次元 TM 波を取り扱う ため該当成分だけを取り出すと ∂Ez ∂y = −µ ∂Hx ∂t (2.3.3) ∂Ez ∂x = µ ∂Hy ∂t (2.3.4) ∂Hy ∂x − ∂Hx ∂y = ε ∂Ez ∂t (2.3.5)
となる。時間についての離散化幅である時間ステップを ∆t とし式 (2.3.1) は時刻 t = n∆t を中心に、式 (2.3.2)は時刻 t = (n +1 2)∆tを中心に、それぞれ時間差分化を行うと次式となる。 µHn+ 1 2 x = µHn− 1 2 x − ∆t ∂En z ∂y (2.3.6) µHn+ 1 2 y = µHn− 1 2 y + ∆t ∂En z ∂x (2.3.7) εEzn+1 = εEn z + ∆t ∂Hn+12 y ∂x − ∂Hn+12 x ∂y ! (2.3.8) ここでの上付き文字は時刻を表す。空間については x 方向に ∆x 間隔で Nxセル、y 方向に ∆y 間隔で Nyセルに離散化する。ここで第 2.2 節で導出した CIP − BS 法の基底関数 Φ を用いると式 (2.3.8) の電 磁場の各成分は次の様に展開される。 εX i,p Ez|n+1,(p)i Φ(p)i (x, y) =ε X i,p Ez|n,(p)i Φ(p)i (x, y) + ∆t X i,p Hy| n+1 2,(p) i ∂Φ(p)i ∂x (x, y)− X i,p Hx|n+ 1 2,(p) i ∂Φ(p)i ∂y (x, y) (2.3.9) ここで下付きの i は i = {(ix, iy) ∈ N2|0 ≤ ix ≤ Nx, 0 ≤ iy ≤ Ny}、上付きの p は p = {(px, py)∈ N2 |0 ≤ px≤ K, 0 ≤ py≤ K} であり展開係数が座標 (ix∆x, iy∆y)の x ついて px次、y について py次 の電磁場の導関数であることを示す。式 (2.3.10) において未知係数は未来の時刻の電界 En+1,p z,i である。 基底関数 Φ(p) i は係数と同数存在しているため、全ての基底関数を試験関数とし式 (2.3.10) の両辺と内積 をとり連立方程式を得る。例としてある基底関数 Φ(q)j (x, y)と式 (2.3.10) の両辺と内積をとると X i,p Ez|n+1,(p)i DΦ(q)j εΦ (p) i E =X i,p Ez|n,(p)i DΦ(q)j εΦ (p) i E + ∆t X i,p Hy|n+12,(p) i * Φ(q)j ∂Φ(p)i ∂x + −X i,p Hx|n+12,(p) i * Φ(q)j ∂Φ(p)i ∂y + (2.3.10) となり、基底関数の台が中心ノードに隣接する 1 セルであることを考慮すると X a,p Ez|n+1,(p)j+a DΦ(q)j εΦ (p) j+a E =X a,p Ez|n,(p)j+a DΦ(q)j εΦ (p) j+a E + ∆t X a,p Hy|n+ 1 2,(p) j+a * Φ(q)j ∂Φ(p)j+a ∂x + −X a,p Hx|n+ 1 2,(p) j+a * Φ(q)j ∂Φ(p)j+a ∂y +! (2.3.11) となる。ここで a = {(ax, ay)∈ N2| − 1 ≤ ax≤ 1, −1 ≤ ay≤ 1} である。式 (2.3.6)、(2.3.7) について も同様に展開し内積をとると次式となる。 X a,p Hx|n+1,(p)j+a DΦ(q)j µΦ (p) j+a E =X a,p Hx|n,(p)j+a DΦ(q)j µΦ (p) j+a E − ∆tX a,p Ez|n+ 1 2,(p) j+a * Φ(q)j ∂Φ(p)j+a ∂y + (2.3.12) X a,p Hy|n+1,(p)j+a D Φ(q)j µΦ (p) j+a E =X a,p Hy|n,(p)j+a D Φ(q)j µΦ (p) j+a E + ∆tX a,p Ez| n+1 2,(p) j+a * Φ(q)j ∂Φ(p)j+a ∂x + (2.3.13)
式 (2.3.11)、(2.3.12)、(2.3.13) と同様に全ての基底関数を試験関数としそれぞれと内積をとると次の連 立方程式が得られる。 MµH n+1 2 x = MµHn− 1 2 x − ∆tMyEzn (2.3.14) MµH n+1 2 y = MµHn− 1 2 y + ∆tMxEzn (2.3.15) MεEzn+1= MεEzn+ ∆t MxH n+1 2 y − MyH n+1 2 x (2.3.16) ここで Hx、Hn n y、Eznは未知数ベクトルである。K = 1 の時、これらのベクトルを次式で定義する。 Hxn= thHx |n,(0,0)0,0 , Hx| n,(0,1) 0,0 , Hx| n,(1,0) 0,0 , Hx| n,(1,1) 0,0 , Hx| n,(0,0) 0,1 , · · · i (2.3.17) Hyn= thHy |n,(0,0)0,0 , Hy| n,(0,1) 0,0 , Hy| n,(1,0) 0,0 , Hy| n,(1,1) 0,0 , Hy| n,(0,0) 0,1 , · · · i (2.3.18) Ezn= t h Ez| n,(0,0) 0,0 , Ez| n,(0,1) 0,0 , Ez| n,(1,0) 0,0 , Ez| n,(1,1) 0,0 , Ez| n,(0,0) 0,1 , · · · i (2.3.19) また Mµ、Mε、Mx、My は係数行列である。式 (2.3.17)、(2.3.18)、(2.3.19) より要素 Fi(p)は (K + 1)2Nyix+(K +1)2iy+(K +1)px+py列の要素であるため I(p)
i = (K +1)2Nyix+(K +1)2iy+(K +1)px+py と置き、試験関数を Φ(q) j とすると Mµ、Mε、Mx、Myは次式で表せる。 [Mµ]I(q) j ,I (p) i = D Φ(q)j µΦ (p) i E (2.3.20) [Mε]I(q) j ,I (p) i =DΦ(q)j εΦ (p) i E (2.3.21) [Mx]I(q) j ,I (p) i = * Φ(q)j ∂Φ(p)i ∂x + (2.3.22) [My]I(q) j ,I (p) i = * Φ(q)j ∂Φ(p)i ∂y + (2.3.23) 各行列の試験関数と基底関数の台が重なる要素、即ち ix− 1 ≤ jx≤ ix+ 1、iy− 1 ≤ jy ≤ iy+ 1である 内積以外は 0 であるため疎行列であることが分かる。 以上の様に更新式は定式化され式 (2.3.14)、(2.3.15)、(2.3.16) を繰返し解くが、時間に関して蛙跳び 差分を行ったため、式 (2.3.14)、(2.3.15) と式 (2.3.16) は交互に解く必要がある。
第
3
章 数値分散解析
数値解析手法にはその手法特有の誤差要因である数値分散が存在する。これにより位相速度に誤差が 生じ、存在しないはずの減衰が現れるなどの影響を及ぼす。この影響の大小は数値解析に用いるパラメー タに依存する。本節では CIP-BS 法の数値分散関係式を導出し数値分散の影響を明らかにすることで適 切なパラメータの決定法を提案する。3.1
数値分散関係式の導出
以下では解析領域全体が真空で光速 c0、誘電率 ε0、透磁率 µ0であり、x 方向を ∆x 間隔、y 方向を ∆y間隔で離散化する状況を考える。CIP-BS 法の解として x 軸に対し θ 方向に伝搬する周波数 ω の平 面波 F (x, y) を考える。この平面波の x について px次、y について py次の導関数の座標 (ix∆x, iy∆y)、 時刻 t = n∆t での値を次の様に置く。∂px+pyF
∂xpx∂ypy(ix∆x, iy∆y) =
1 ∆xpx∆ypy
˜
F(px,py)ejωn∆t−j2π˜λk(cos θ ix∆x+sin θ iy∆y) (3.1.1)
ここで ˜k は規格化数値波数であり数値波数ベクトルは2π˜k λ (cos θ, sin θ)と表す。また ˜F (px,py)は ∆x、∆y によって無次元化された導関数の振幅である。基底関数を無次元化すると次式となる。 ˜ Φ(p)i (x, y) = 1 ∆xpx∆ypyΦ (p) i (x, y) (3.1.2) ここで i = (ix, iy)、p = (px, py)である。式 (3.1.1) と式 (3.1.2) を用いて時刻 t = n∆t の平面波 F(n∆t, x, y)を展開すると次式となる。 F(n∆t, x, y) =X i,p ˜
F(p)ejωn∆t−j2π˜λk(cos θ ix∆x+sin θ iy∆y)Φ˜(p)
i (x, y) (3.1.3)
x方向サンプリング密度 νx= λ
∆x、y 方向サンプリング密度 νy =∆yλ 、クーラン数 ξ = c0∆t q 1 ∆x2+ 1 ∆y2、 を定義すると以下の式が成り立つ。 ω∆t = 2π λc0∆t = 2π λ ξ q 1 ∆x2+∆y12 =q2πξ ν2 x+ νy2 (3.1.4) パラメータ νx、νy、ξ を用いて η と ζ(i) を次式で定義する。 η= ejω∆t= e j√2πξ ν2x+ν2 y (3.1.5)
ζ(i) = e−j2π˜λk(cos θ ix∆x+sin θ iy∆y)= e−j2π
˜
kix cos θνx +iy sin θνy
(3.1.6) 式 (3.1.5)、(3.1.6) を式 (3.1.3) に代入すると次式となる。 F(n, ∆x, y) = ηnX p ˜ F(p)X i ζ(i) ˜Φ(p)i (x, y) (3.1.7)
となり、式 (3.1.7) を用いて時間差分式の式 (2.3.8) の電磁場の各方向成分を展開すると次式となる。 ηn+1X p ˜ Ez(p) X i ζ(i) ˜Φ(p)i (˜x,y) =η˜ nX p ˜ Ez(p) X i ζ(i) ˜Φ(p)i (˜x,y)˜ + ηn+12ξsin β X p Z ˜Hy(p) X i ζ(i) ∂ ˜Φ (p) i ∂x˜ (˜x,y)˜ − ηn+1 2ξcos β X p Z ˜Hx(p) X i ζ(i) ∂ ˜Φ (p) i ∂y˜ (˜x,y)˜ (3.1.8) ここで Z =qµ0 ε0 は真空中の波動インピーダンス、˜H (p) x 、˜Hy(p)、˜Ez(p)は式 (3.1.7) の ˜F(p)に対応する無次元 化された電磁場の導関数の振幅、β はセルの対角線と x 軸とのなす角度で sin β = ∆y √ ∆x2+∆y2 = νx √ν2 x+ν2y 、 cos β = √ ∆x ∆x2+∆y2 = νy √ ν2 x+νy2 であり、また ˜x = x ∆x、˜y = y ∆yで ∂ ˜Φ(p)i ∂x˜ = ∆x ∂ ˜Φ(p)i ∂x 、 ∂ ˜Φ(p)i ∂y˜ = ∆y ∂ ˜Φ(p)i ∂y で ある。式 (3.1.8) の両辺と試験関数 Φq j の内積をとり両辺を ηnζ(j)で割ると次の様になる。 ηX p ˜ E(p)z X a ζ(a) D ˜Φ(q)j Φ˜ (p) j+a E =X p ˜ Ez(p) X a ζ(a) D ˜Φ(q)j Φ˜ (p) j+a E + η12ξsin β X p Z ˜Hy(p) X a ζ(a) * ˜ Φ(q)j ∂ ˜Φ(p)j+a ∂x˜ + − η12ξcos β X p Z ˜Hx(p) X a ζ(a) * ˜ Φ(q)j ∂ ˜Φ(p)j+a ∂y˜ + (3.1.9) ここで a = {(ax, ay)∈ N2 | − 1 ≤ ax≤ 1, −1 ≤ ay ≤ 1} である。無次元化された基底関数同士の内積は 付録に示す内積の ∆x+と ∆x−を 1 と置き換えたものと一致する。式 (2.3.8) と同様に式 (2.3.6)、(2.3.7) についても式 (3.1.8) から式 (3.1.9) の操作を行うと以下の連立方程式を得る。 η ˜MZHx= ˜MZHx− η12ξcos β ˜MyEz
η ˜MZHy= ˜MZHy+ η12ξsin β ˜MxEz
η ˜M Ez= ˜M Ez+ η12ξsin β ˜MxZHy− η 1 2ξcos β ˜MyZHx (3.1.10) ここで Hx、Hy、Ezは未知数ベクトルであり (K + 1)px+ py番目の要素がそれぞれ ˜Hx(p)、 ˜Hy(p)、 ˜Ez(p) である。また ˜M、 ˜Mx、 ˜Myは ((K + 1)qx+ qy,(K + 1)px+ py)≡ (m, n) が次式となる。 [ ˜M]m,n=X a ζ(a) D ˜Φ(q)j Φ˜ (p) j+a E (3.1.11) [ ˜Mx]m,n=X a ζ(a) * ˜ Φ(q)j ∂ ˜Φ(p)j+a ∂x˜ + (3.1.12) [ ˜My]m,n=X a ζ(a) * ˜ Φ(q)j ∂ ˜Φ(p)j+a ∂y˜ + (3.1.13) 式 (3.1.10) を Ezについて整理すると次の様な方程式となる。 (η− 1)2M˜ − ξ2sin2 β η ˜MxM˜−1M˜x− ξ2cos2β η ˜MyM˜−1M˜y Ez= 0 (3.1.14) 式 (3.1.14) が非自明解を持つためには解は一意に決定出来てはならない。即ち次式の様に係数行列の行 列式は 0 とならなくてはならない。 (η− 1) 2M˜ − ξ2sin2β η ˜MxM˜−1M˜x− ξ2cos2β η ˜MyM˜−1M˜y = 0 (3.1.15) この式 (3.1.15) が CIP-BS 法の数値分散関係式である。式 (3.1.11)、(3.1.12)、(3.1.13) より係数行列は あるノード番号 j に関するものである事が分かるが、均一に離散化した状況を考えているため j を変化 させても内積の値は変化せず一般性は保たれる。
3.2
解析結果
本節では数値分散関係式である式 (3.1.15) より数値分散におけるパラメータの影響を明らかにし適切 なパラメータ決定法を提案する。 以下では νx= νyの正方形のセルで領域が離散化された状況を考え、式 (3.1.15) を解き数値分散の性質 を明らかにする。比較のため式 (3.2.1) で表される FDTD 法の数値分散関係式より [8]、得られた FDTD 法の数値分散の影響を示す。 νx2+ νy2 ξ2 sin 2 πξ q ν2 x+ νy2 = νx2sin2 π˜k νx cos θ ! + νy2sin2 π˜k νy sin θ ! (3.2.1)3.2.1
サンプリング密度への依存性
FDTD法はクーラン条件によって与えられる時間ステップの上限に近いほど、即ちクーラン数 ξ が 1 に近いほど位相速度は厳密な値に接近する [8]。そこで K = 1、ξ 一定とした時の CIP-BS 法と ξ = 0.99 とした時の FDTD 法の数値位相速度とサンプリング密度の関係を図 3.1 に示す。図 3.1 では縦軸が数値 位相速度 v を真空中の光速 c0で割った規格化位相速度 v c0 であり、1 に近づくほど v が厳密な値に近い事 を表す。図 3.1 より高橋らの報告と同様に [6]、CIP-BS 法は ξ が小さくなるほど位相速度は厳密な値に 近づくことが確認できる。またいずれの手法でも νxが大きくなるにつれ位相速度が厳密な値に近づく事 が分かるが、その傾きは ξ によって異なる。まず ξ = 0.1 では FDTD 法と同様に O(ν−2 x )の傾きを持ち、 ξ= 0.001では O(ν−5.5 x )となり、早く収束している。そして ξ = 0.01 は θ = 0 度では νx= 10、θ = 45 度では νx= 5をそれぞれ境に O(ν−2 x )と O(νx−5.5)の傾きが共存しており、O(ν−5.5x )となっている領域 では ξ = 0.001 の位相速度分布と一致している。即ち、ξ を減少させ続けても O(ν−5.5 x )の傾きを持つ分 布に収束するため、適切な νxと ξ の組み合わせでパラメータを設定する必要がある。また O(ν−2 x )の傾 きを持つ領域では伝搬方向の違いによる変化は無いこと。よって O(ν−2 x )の傾きを持つ分布上の位相速 度を取る様にパラメータを設定することで誤差の大きい伝搬方向の成分に全体の精度が依存する様な状 況を回避できる。 数値波数が虚部を持つとき、伝搬距離に応じた数理減衰が生じる。図 3.2 に K = 1 の CIP-BS 法にお ける数値減衰とサンプリング数の関係を示す。伝搬方向 θ = 45 度の場合、νx ≥ 1 範囲では数値減衰が 確認されなかった。数値減衰は νx<1.4でしか確認出来ず、実際のシミュレーションでは使用しない範 囲であるため減衰は考慮する必要はないと考えられる。1+10
-8
1+10
-7
1+10
-6
1+10
-5
1+10
-4
1+10
-3
1.01
1.1
2
4
6
8 10 12
20
1-10
-8
1-10
-7
1-10
-6
1-10
-5
1-10
-4
1-10
-3
0.99
0.9
Normalized phase velocity
v/c
0
for CIP-BS
Normalized phase velocity
v/c
0
for FDTD
Sampling density
ν
x
=
ν
y
FDTD ξ = 0.99
CIP-BS ξ = 0.1
CIP-BS ξ = 0.01
CIP-BS ξ = 0.001
(a)伝搬方向 θ =0 度の場合1+10
-8
1+10
-7
1+10
-6
1+10
-5
1+10
-4
1+10
-3
1.01
1.1
2
4
6
8 10 12
20
1-10
-8
1-10
-7
1-10
-6
1-10
-5
1-10
-4
1-10
-3
0.99
0.9
Normalized phase velocity
v/c
0
for CIP-BS
Normalized phase velocity
v/c
0
for FDTD
Sampling density
ν
x
=
ν
y
FDTD ξ = 0.99
CIP-BS ξ = 0.1
CIP-BS ξ = 0.01
CIP-BS ξ = 0.001
(b)伝搬方向 θ =45 度の場合 図 3.1: K = 1 の CIP-BS 法と数値位相速度 v c0 とサンプリング密度 νxの関係0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1
1.2
1.4
1.6
1.8
2
Attenuation [Np/
λ
]
Sampling density
ν
x
=
ν
y
CIP-BS
ξ = 0.1
CIP-BS
ξ = 0.01
CIP-BS
ξ = 0.001
図 3.2: K = 1 の CIP-BS 法における伝搬方向 θ = 0 度の数値減衰とサンプリング密度 νxの関係3.2.2
高次化の影響
CIP-BS法において高次化の試みは今までに成されてきており、K = 2 の CIP-BS 法は時間依存の Schrodinger方程式への適用が行われ、高精度化に成功している [2][3]。そこで本節では Maxwell の方程 式では高次化はいかなる意味を持つのかを明らかにするため K = 1, 2 の時の数値分散の影響を比較す る。図 3.3 に K = 1, 2 の CIP-BS 法におけるサンプリング密度 νx= νyと数値位相速度の関係を示す。 図 3.3 において実線は数値位相速度が厳密な値 c0より速い領域 (左縦軸を参照)、破線は c0より遅い領 域 (右縦軸を参照) を意味している。K = 1 の位相速度は図 3.1 で見た通り、ξ が大きいとき O(ν−2 x )の 傾きを持ち、ξ が減少すると O(ν−5.5 x )の分布に収束していく。一方 K = 2 の場合は、ξ が大きい場合は K= 1の時と同様 O(ν−2 x )の傾きをもつが、K = 1 より高精度まで ξ の減少と共に厳密な位相速度に接 近し、ξ = 10−7の分布に収束する。その収束の仕方は ξ が減少していくにつれ c0より位相速度が遅く なる区間がより大きい νxまで徐々に広がり、ξ = 10−7まで減少すると観測した νxの区間全体が c0よ り遅くなっている。観測した区間内で最も厳密値に近いのは ξ = 10−7、νx∈ [7, 10] の時であり、この区 間のパラメータを用いることが効率的であると考えられるが、そのときの位相速度の相対誤差は 10−13 とかなり小さなものとなっている。しかし単精度の仮数部は 7 桁、倍精度でも 15 桁程であり丸め誤差 を考えると過剰な高精度である。また ξ が大きい場合 K = 1 の O(ν−2 x )の傾きをもつ分布と一致してお り、高次化の意味がない。そして K = 1 のときの未知数はセル数の 4 倍に比例していたが、K = 2 の倍 はセル数の 9 倍に比例する。LU 分解の計算時間が未知数の 2 乗に比例することを考えれば、K = 2 の ときの計算時間は K = 1 のときの81 16 倍程度であり、4 倍以上の計算時間を要する。したがって高次化 で高精度になるとは限らず、なったとしても効率が良いとは考えにくい。もし用いるのであれば 4 倍精 度を必要とする様な高精度なシミュレーションを行う際に ξ = 10−7、νx∈ [7, 10] 程度のパラメータを選 ぶと効率が良いと考えられる。1+10
-14
1+10
-12
1+10
-10
1+10
-8
1+10
-6
1+10
-4
1+10
-2
2
4
6 8 10 12
20
1-10
-14
1-10
-12
1-10
-10
1-10
-8
1-10
-6
1-10
-4
1-10
-2
Normalized phase velocity
v/c
0
Normalized phase velocity
v/c
0
Sampling density
ν
x
=
ν
y
(a)伝搬方向 θ =0 度の場合1+10
-14
1+10
-12
1+10
-10
1+10
-8
1+10
-6
1+10
-4
1+10
-2
2
4
6 8 10 12
20
1-10
-14
1-10
-12
1-10
-10
1-10
-8
1-10
-6
1-10
-4
1-10
-2
Normalized phase velocity
v/c
0
Normalized phase velocity
v/c
0
Sampling density
ν
x
=
ν
y
(b)伝搬方向 θ =45 度の場合 K = 2, ξ = 10-1 K = 2, ξ = 10-2 K = 2, ξ = 10-3 K = 2, ξ = 10-4 K = 2, ξ = 10-6 K = 2, ξ = 10-7 K = 1, ξ = 10-1 K = 1, ξ = 10-2 K = 1, ξ = 10-3 図 3.3: K = 1, 2 の CIP-BS 法と数値位相速度 v c0 とサンプリング密度 νxの関係 (実線は c0より速い領 域であり左縦軸を参照、破線は c0より遅い領域であり右縦軸を参照)3.2.3
時間ステップと数値位相速度の関係
第節では報告と同様に位相速度が時間ステップ ∆t の大きさを決定する因子の一つであるクーラン数 ξが小さくなるほど精度が向上し、またサンプリング密度増加に伴う精度向上割合も改善を見せた。そ こで本節では ∆t と数値位相速度を調査し適切なパラメータの決定法を提案する。 K= 1の CIP-BS 法と FDTD 法における νx一定の下での数値位相速度 v と ∆t の関係を図 3.4、3.5 に示す。図 3.4 より時間ステップ ∆t が小さくなるにつれ数値位相速度が厳密な値に近づく事が分かる が、FDTD 法では図 3.5 に示す通り ∆t がクーラン条件により決定される ∆t の上限から離れるほど位相 速度は厳密値から離れるためこの性質は CIP-BS 法特有の性質であると考えられる [8]。図 3.4 の黒い破 線は O(∆t2)の傾きをもつ近似曲線であり ∆t が大きい時、v は λ あたりのセル数 νxに依らずに O(∆t2) で厳密値 c0に近づく。一方 ∆t が減少していくと黒い破線から外れオーダーが小さくなっていくが、サ ンプリング密度 νxが多いほどが高い精度まで O(∆t2)の傾きで厳密値に近付いていく。そして ∆t を十 分に小さくし位相速度が収束する値図 3.1 で確認した通り ν−5.5 x に比例する。黒い破線上では νxに依ら ず ∆t のみに不均一離散化境界での反射波が抑制できると考える。また FDTD 法では大きな伝搬方向依 存性が存在するが、CIP-BS 法では黒い破線は伝搬方向に依らずほぼ一致している。そのためより大きな 時間ステップで黒い破線から位相速度分布が離れてしまう伝搬方向 0 度の波の数値位相速度を黒い破線 と一致するようパラメータを設定することで、精度の悪い伝搬方向の波に誤差が依存することを避ける ことが可能である。それゆえパラメータを決定する際は黒い破線を参照し所望の精度が得られる ∆t を 決定し、その ∆t で黒い破線と一致する νxの内最も小さいものを選択することが望ましい。1+10
-71+10
-61+10
-51+10
-41+10
-30.001
λ/c
00.01
λ/c
0Normalized phase velocity
v/c
0Time step
∆t
Ο(∆t2) νy = νx = 3 4 5 6 7 8 9 10 (a)伝搬方向 θ =0 度の場合1+10
-71+10
-61+10
-51+10
-41+10
-30.001
λ/c
00.01
λ/c
0Normalized phase velocity
v/c
0Time step
∆t
Ο(∆t2) νy = νx = 3 4 5 6 7 8 9 10 (b)伝搬方向 θ =45 度の場合 図 3.4: CIP-BS 法における数値位相速度 v c0 と時間ステップ ∆t の関係0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
0.02
λ/c
00.04
λ/c
00.06
λ/c
0Normalized phase velocity
v/c
0Time step
∆t
(a)伝搬方向 θ =0 度の場合0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
0.02
λ/c
00.04
λ/c
00.06
λ/c
0Normalized phase velocity
v/c
0Time step
∆t
ν
x= ν
y=10
9
8
7
6
5
4
(b)伝搬方向 θ =45 度の場合 図 3.5: FDTD 法における数値位相速度 v c0 と時間ステップ ∆t の関係第
4
章 散乱体のモデリングに関するパラメータ
の検討
散乱問題を数値解析する際、セルに沿わない散乱体の形状のモデリングのため必要以上に離散化を 行う必要があり計算資源の浪費が生じる。これは広大な解析領域を取り扱う場合顕著な影響を受ける。 CIP-BS法では一つのノードに対し 2(K + 1) 個の未知数が対応するためセル数の増大に伴う計算量の増 加が他の手法に比べ大きい。また第 3.2 節で明らかにした通り離散化を詳細に行ったとしても精度の向上 につながるとは限らない。そこで CIP-BS 法に不均一離散化を適用し効率的な計算資源の利用を試みる。4.1
CIP-BS
法基底関数の拘束条件を満たす不均一離散化
不均一に離散化する場合、セルの大きさが異なる境界が生じるが、CIP-BS 法の基底関数は式 (2.1.6)、 (2.1.7)で表され i 番目のノードに隣接するセルの大きさ ∆x+及び ∆x−が異なるとしても定義可能で ある。また CIP 法補間関数の拘束条件である式 (2.1.2) よりセルの両端の電磁場とその導関数をもって 補間関数は拘束されなくてはならない。FDTD 法などの不均一離散化手法であるサブグリッド法では図 4.1の様な離散化がなされるが、この方法では不均一離散化境界においてセルの両端にノードが配置され ていないため関数が定義出来ない。そのため図 4.2 の様に長方形のセルを用いセルのノード同士が重な る様に離散化する。Fine sampling
図 4.1: サブグリッド法で用いられる離散化Fine sampling
図 4.2: CIP-BS 法基底関数の拘束条件を満たす離散 化2d
2d
rε =2
Incident
wave
E
zH
tA
B
C
D
d
k
図 4.3: 無限長円筒散乱体に対する散乱問題の状況設 定の模式図 m 8 cells 8 cells 8 cells 8 cells5d
5d
10 cells/ d
10 cells/ d
cells/
N
A B C D TFSF NPMLd
図 4.4: 数値解析に用いるパラメータの模式図4.2
数値解析
本節では円筒散乱体の散乱問題について不均一離散化を適用した CIP-BS 法の数値解の誤差を調査し 散乱体のモデリングに必要なセル数を調査する。4.2.1
問題設定及び数値解析に用いるパラメータ
図 4.3 に示す通り、直径 d、比誘電率 εr= 2の無限長円筒散乱体に対する散乱問題を考える。この散 乱体に図の矢印の向きに波長 λ = d, 2d の平面波を入射しそれぞれの散乱場を確認する。観測地点は円筒 散乱体の中心を重心とした 1 辺 2d の正方系 ABCD である。ここで辺 DA で後方散乱、辺 BC で前方 散乱を観測できる様に配置する。 CIP-BS法で数値計算を行う際、図 4.4 に示す通り解析領域を 1 辺 5d の正方形としその中心に円筒誘 電体を配置する。入射波は電界が次式となる中心周波数 ω =2π λc0のパルス波とする。 Ez(x, y, t) = E0sin(ωt− kx) exp −(x− c0t) 2 2σ2 p (4.2.1) ここで E0は振幅、k = 2π λ は波数、σp= 0.32dである。d あたりのセル数は自由空間では 10 セル散乱体内 部では Nmセル、となる様に離散化し、時間ステップ ∆t は自由空間中の自由空間中のセルサイズ ∆x = d 10 を用いて ∆t = 0.01∆x c0 √ 2 と設定する。不均一離散化境界の 2 セル外側を囲む様に Total Field/Scattered Field(TFSF)境界を配置し入射波と散乱場を分離し、解析領域の外側に 8 層の Nearly Perfectly Matched Layer(NPML)を配置し反射波を抑制する [6][9]。4.2.2
計算結果
図 4.5、4.6 に散乱場の厳密解と Ns= 10の CIP-BS 法の数値解を示す。図 4.5 は入射波の波長が λ = d、 図 4.6 は λ = 2d の時の散乱場である。
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
A
B
C
D
A
Normalized amplitude
Observation points
Exact solution
CIP-BS N
m
= 10
(a)入射波に対する規格化振幅-π
-π/2
0
π
/2
π
A
B
C
D
A
Phase [rad]
Observation points
(b)位相 図 4.5: λ = d の入射波に対する散乱場の (a) 振幅及び (b) 位相0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
A
B
C
D
A
Normalized amplitude
Observation points
Exact solution
CIP-BS N
m
= 10
(a)入射波に対する規格化振幅-π
-π/2
0
π
/2
π
A
B
C
D
A
Phase [rad]
Observation points
(b)位相 図 4.6: λ = 2d の入射波に対する散乱場の (a) 振幅及び (b) 位相この散乱場に現れた誤差を定量的に調査するため、振幅を式 (4.2.2) で表される平均二乗相対誤差(Mean Square Relative Error:MSRE)、位相を式 (4.2.3) で表される平均二乗誤差 (Mean Square Error:MSE) で 評価する。 s 1 L ˆ L FCIP-BS− Fexact Fexact 2 dl (4.2.2) s 1 L ˆ L (FCIP-BS− Fexact) 2 dl (4.2.3) ここで L は観測地点の全体の長さ、FCIP-BSが数値解、Fexactが厳密解である。この振幅と位相の誤差と Nmの関係を図 4.7 に示す。ここでの誤差要因は数値分散による誤差とモデリングの粗さによる誤差が 考えられる。振幅の誤差は入射波長に関わらず O(N−1.4 m )の傾きをもつため、同一の誤差要因が支配的 であると考える。第 3.2 節より 1 波長あたり 10 セル以上で離散化を行っているため数値減衰は存在しな い。また位相速度の変化が振幅に影響を及ぼす場合、本来観測するべき周波数とは別の成分を観測した ということになるが、式 (5.1.1) のスペクトルは多項式オーダではないため Nmの変化によって数値分散 に影響が現れたとは考え難い。そのため振幅誤差はモデリングの粗さによる要因が支配的であると考え られる。 一方、位相誤差は傾きが一定ではなく、Nm<12のとき入射波長に関わらず O(N−1 m )の傾きを持ち、 Nm≥ 15 のとき λ = d では O(N−0.5 m )、λ = 2d では O(Nm−0.8)の傾きを持つ。いずれの入射波長の場合 も傾きは小さくなっているがこれは数値分散による影響であると考える。図 3.4 より ∆t 一定の下でサン プリング密度が増加した場合、位相速度は図中の黒い破線に収束している。図 4.7(b) でも数値分散の影響 が支配的となれば誤差の収束が生じると考えられ Nmの増大によって傾きが減少したのはその影響であ ると考える。また入射波長によって位相誤差の傾きの変化に差が生じたのは λ = d の数値解より λ = 2d の解の方が 1 波長あたりのセル数が多いため数値分散は少ないためであると考える。事実、数値分散解 析より自由空間中をグリッドに沿った方向に伝搬する波の数値位相速度は λ = d のとき ∆t =0.0007λ c0 な ので v = (1 + 10−6)c0、λ = 2d のとき ∆t = 0.00035λ c0 なので v = (1 + 0.2 × 10 −6)c0であり λ = 2d の方 が厳密値に近い。即ち位相誤差は傾きが O(N−1 m )である領域では数値分散よりモデリングの誤差が支配 的であり、Nmを大きくしていくと傾きが小さくなり数値分散が支配的となったと考える。 以上の結果から位相誤差に比べ振幅誤差はモデリングの影響が強く現れるため、モデリングに要する セル数を決定する場合は振幅の精度を基準に決定すべきであると考えられる。
0.1
0.2
0.4
0.6
0.8
1.0
2.0
10
15
20
25
30
MSRE of amplitude [
×10
-2
]
Num. of cells per diameter of the scatterer N
m
[cells]
Incident wave length
λ
λ = d
λ = 2d
(a)振幅の平均二乗相対誤差0.4
0.6
0.8
1.0
2.0
10
15
20
25
30
MSE of phase [
×10
-2
rad]
Num. of cells per diameter of the scatterer N
m
[cells]
Incident wave length
λ
λ = d
λ = 2d
(b)位相の平均二乗誤差
第
5
章 大規模散乱シミュレーション
第 3.2 節にて数値位相速度は O(∆t2)で厳密解に近づくことが明らかとなった。そしてそれは ∆t が大 きければセルサイズに依らずに ∆t だけで精度が決定されるため粗い離散化でも高精度な計算が実現で きる。そこで本節では波長に対し広大な解析領域を考慮する散乱問題に CIP-BS 法を適用し現在一般的 に用いられている手法である FDTD 法と比較することで利点を調査する。Incident
wave
=2
H
tE
z15λ
k
6λ
6λ
6λ
6λ
r
ε
λ
15λ
D
C
B
A
図 5.1: 周期的に並んだ無限長円筒散乱体に対する散 乱問題の模式図 λ 17 λ 17 D A B C NPML or PML 8 cells 8 cells 8 cells 8 cells TFSF 図 5.2: 数値解析に用いるパラメータの模式図 (FDTD 法は均一に離散化)5.1
数値解析
5.1.1
問題設定及び数値解析に用いるパラメータ
CIP-BS法と FDTD 法の比較のため、図 5.1 に示す様な周期的に並んだ無限長円筒散乱体の散乱問題 を解析する。波長を λ としたとき無限長円筒散乱体の直径は λ、比誘電率は εr= 2であり、これが 3 × 3 のアレー状に中心間距離 6λ の間隔で並んでいる。ここに図 5.1 の矢印の方向に平面波を入射する。観測 地点はアレー中央の散乱体中心を重心とする 1 辺 15λ の正方形 ABCD 辺上である。 以上の散乱問題の数値解析を行う際は図 5.2 の様にパラメータを設定する。解析領域を 1 辺 17λ の正 方形としこの重心とアレーの中央の散乱体中心が一致する様に散乱体を配置する。入射波電界は次式で 表される中心周波数 ω =2π λc0のパルス波とする。Ez(x, y, t) = E0sin(ωt− kx) exp −(x− c0t) 2 2σ2 p (5.1.1)
ここで E0は振幅、k =2π λ は波数、σp= 0.32λである。 CIP-BS法ではサンプリング密度が自由空間中で Nsセル散乱体内部で Nd = 20セルとなる様に離散 化する。各散乱体の 1 セル外側を囲むように TFSF 境界を配置し入射波と散乱場を分離する。解析領域 の外側に 8 層の NPML を配置し反射波を抑制する。更新式のための連立方程式解法には LU 分解か共役 勾配 (Conjugate Gradient:CG) 法を用いる。CG 法において前処理行列は係数行列の対角成分を取り出 した行列を用いる。 一方、FDTD 法では λ あたりのセル数を Ns= Ndセルとし均一に離散化する。時間ステップは安定 条件を満たし、良好な精度で計算できる様に自由空間中のセルサイズ ∆x = λ Ns を用いて、∆t = 0.99∆x c0 √ 2 と設定する [8]。各散乱体の 4 セル外側を囲むように TFSF 境界を配置し入射波と散乱場を分離し、解析 領域の外側に 8 層の PML を配置し反射波を抑制する [8]。
5.1.2
計算結果
前節で設定した散乱問題の数値解を定量的に評価し、FDTD 法と CIP-BS 法が同精度であるようなパラ メータを明らかにする。そしてそのパラメータでのメモリ使用量及び計算時間を比較することで CIP-BS 法の評価を行う。 図 5.3 に評価対象である正方形 ABCD 辺上の散乱場の振幅と位相を示す。CIP-BS 法の結果は LU 分 解を用いたものである。Ns=1000セルとした時の FDTD 法の散乱場を参照解として振幅と位相の誤差 を評価する。CIP-BS 法と FDTD 法による数値解の振幅の平均二乗相対誤差 (MSRE) をそれぞれ図 5.4、 5.5に示す。CIP-BS 法の Ns= 4の振幅誤差は他のサンプリング密度に比べ 0.5×10−2程大きい。これは セルが大きく TF 領域が広くなる事で生じる TFSF 境界での入射波の分離精度の低下や不均一離散化境界 での反射などの誤差要因が強く現れていると考えられる。また Ns= 6, 8の振幅の誤差は ∆t > 8 × 10−3 で指数的な誤差の低下が見られ、Ns = 6では ∆t = 4 × 10−3 λ c0、Ns = 8では ∆t = 6 × 10 −2 λ c0 程で 振幅誤差は 0.6 × 10−2付近に収束している。誤差の指数的な減少は位相速度が ∆t が減少したことで数 値位相速度が厳密値に近くなり、観測する周波数成分が対象成分に接近する際に入射波のスペクトルが 現れたものであると考えられる。また図 5.4 では散乱体内部のサンプリング密度を Nd =20としていた が、Nd =30に変更し、Ns = 6、∆t = 0.071λ c0 と設定して CIP-BS 法を用いて解析した結果、振幅誤 差は 0.51 × 10−2となった。モデリングの粗さに対する振幅誤差の依存性は図 4.7(a) に示しており、そ の結果では円筒散乱体の直径あたりのセル数の −1.4 乗に比例していたが、今回の結果では −0.15 乗に 比例する程度であるため、モデリングの粗さは支配的でなく Nd = 20で十分であると考えられる。図 5.5は FDTD 法の振幅の誤差を示していおり、振動を繰り返す複雑な分布が確認できる。CIP-BS 法で も Ns= 6の誤差より Ns= 8の誤差の方が低い領域が見て取れ、Nsの変化によって振動することが予 想できる。しかし CIP-BS 法の ∆t の変動には振幅誤差は敏感な反応は見せないため、高精度化に際し て他の手法より Nsの変化が少ない CIP-BS 法は解の安定性は高いと言える。0
0.5
1
1.5
2
2.5
A
B
C
D
A
Normalized amplitude
Observation points
FDTD, N
s=1000
CIP-BS,N
s=6,
∆t = 0.013λ/c
0 (a)入射波に対する規格化振幅-π
-π/2
0
π
/2
π
A
B
C
D
A
Phase[rad]
Observation points
(b)位相 図 5.3: 周期的に並んだ散乱体を囲む正方形 ABCD 辺上での散乱場の (a) 振幅及び (b) 位相0.5
1.0
1.5
2.0
2.5
2
4
6
8
10
12
MSRE of amplitude [
×10
-2
]
Time step
∆t [×10
-3
λ/c
0
]
N
s
= 4
N
s
= 6
N
s
= 8
図 5.4: CIP-BS 法での振幅の平均二乗相対誤差0.1
0.5
1.0
5.0
100
200
300
400
MSRE of amplitude [
×10
-2
]
Sampling density N
s
of FDTD
図 5.5: FDTD 法での振幅の平均二乗相対誤差CIP-BS法と FDTD 法による数値解の位相の平均二乗誤差 (MSE) を図 5.6 に示す。振幅と同様に CIP-BS法では LU 分解と CG 法を用いたが図 5.6 では FDTD 法と Ns= 6の CIP-BS 法の結果の位相 誤差の近似曲線が一致する様にプロットしている。図 5.6 より CIP-BS 法の位相誤差は O(∆t0.9)の傾き を持ち、FDTD 法の位相誤差は O(N−1.5 s )の傾きを持つことが確認できる。また CIP-BS 法の Ns= 4 の誤差だけが ∆t の減少に伴い減少割合が小さくなっていることが確認できる。これは図 3.4 で見られた 数値位相速度の収束が生じたものであると考えられる。即ちこの位相誤差は数値分散による影響が支配 的であり、この区間では CIP-BS 法の Nsを 6 以上大きくしても精度は向上しない。そこでこのスケール で FDTD 法と Ns= 6の CIP-BS 法の使用メモリと計算時間を比較するとそれぞれ図 5.7、5.8 となる。 ここで使用メモリ量は仮想メモリ量であるが、どの手法においてもスワップは発生していないため計算 時間にメモリの書き出し速度による差は考慮しない。CIP-BS 法は高精度化に当たってセル数を増加さ せ無いため使用メモリ量は増加していない事が図 5.7 から確認できる。CIP-BS 法の更新式を LU 分解 で解くと FDTD 法に比べ多くのメモリを利用しなくてはならないが、CG 法で解けば高精度な領域では FDTD法より少ないメモリで計算可能である。計算時間を見ると CG 法に比べ LU 分解が高速であるこ とわかるが、それでも FDTD 法と比べると時間を要する。しかし計算時間の増加は FDTD 法は O(N3 s) 傾きを持ち、CIP-BS 法は O(∆t−1)の傾きを持つが図 5.8 中で比較すると FDTD の法の方が大きな傾き であることが分かる。それ故、解析領域に対し散乱体の分布密度がより小さく、またより広大な解析領 域を有する状況で CIP-BS 法は適していると考えられる。
10
-2
10
-1
2
4
6
8
10
100
200
300
Mean square error of phase [rad]
Time step
∆t of CIP-BS [×10
-3
λ/c
0
]
Sampling density N
s
of FDTD
FDTD
CIP-BS N
s
=
4
CIP-BS N
s
=
6
CIP-BS N
s
=
8
図 5.6: CIP-BS 法及び FDTD 法数値解の位相の平均二乗誤差0.2
0.4
0.6
0.8
1
2
4
6
2
4
6
8
10
100
200
300
Used virtual memory [GB]
Time step
∆t of CIP-BS [×10
-3
λ/c
0
]
Sampling density N
s
of FDTD
FDTD
CIP-BS, N
s
= 6, CG
CIP-BS, N
s
= 6, LU
図 5.7: 使用メモリ量10
100
1000
2
4
6
8
10
100
200
300
Calculation time [min]
Time step
∆t of CIP-BS [×10
-3
λ/c
0
]
Sampling density N
s
of FDTD
FDTD
CIP-BS, N
s= 6, CG
CIP-BS, N
s= 6, LU
第
6
章 結論
本研究では CIP-BS 法の粗い離散化でも時間ステップを減少することで位相誤差を低減可能であると いう性質に着目し、波長に対し広大な解析領域を考慮する問題への適用とその評価を行った。 まず CIP-BS 法の特性がどの程度の粗い離散化でも効果を発揮するのかを調査するため数値分散解析 を行った。K = 1 の時クーラン数 ξ 一定の下でサンプリング数 νxと数値位相速度の関係を調査した。νx を減少させると、ξ が大きいときは O(ν−2 x )の傾きで厳密な速度に接近し、ξ が小さくなると O(νx−5.5) の傾き持つ分布に収束した。O(ν−2 x )の分布は伝搬方向依存性が見られなかった。数値減衰は νx <1.4 でしか確認出来ず実際のシミュレーションでは用いない範囲であるため減衰は無視できると考えた。ま たサンプリング密度 νx一定のもとで時間ステップ ∆t と数値位相速度の関係を調査した結果、∆t を減 少させていくと、νxが大きいときは νxに関係無く数値位相速度は O(∆t2)で厳密値に接近していくが、 νxが小さくなると傾きは小さくなり数値位相速度は収束し、その収束した値の νxの-5.5 乗に比例する ものである。よってパラメータ決定に際しパラメータ決定に際し νx大きいときに所望の精度を得ること ができる ∆t を選択し、その中で νxを最小になる様に決定すべきであると提案した。 そして高次化した CIP-BS 法についても数値分散解析を行った。K = 2 としたとき、ξ 一定のもとで νxと数値位相速度の関係を確認した。K = 2 の時も K = 1 のときと同様、ξ が大きければ O(νx−2)の傾 きで厳密値に接近し、K = 1 の同じ ξ でも O(ν−2 x )の傾きを有するとき、K = 1 と K = 2 の数値位相速 度の分布は一致し、ξ が大きいときは K = 1 で十分であると考えた。また K = 2 のとき ξ を減少させて いくと収束をみせ ξ = 10−7、νx∈ [7, 10] が程度を選択すると効率が良いことが明らかとなったが、丸め 誤差を考えると実用的でないと判断した。 次に解析領域内にグリッドに沿わない散乱体が存在するとき、モデリングのためにその領域だけは詳 細に離散化する必要があるため CIP-BS 法に不均一離散化を適用した。そのモデリング精度の評価のた め T Mzモードにおける無限長円筒散乱体の散乱問題を解析し散乱場の誤差を調査した。結果、振幅誤 差は散乱体を詳細に離散化することでその直径上のセル数 Nmの −1.4 乗に比例する低減を見せた。こ れは入射波の波長を変化させても変わらないためモデリングによる影響であると考えられる。一方、位 相誤差も詳細な離散化により誤差の低減が見られるが、Nm>12のとき Nmの増加と共に誤差の低減割 合は減少した。これは Nmの増加と共に数値分散の影響が支配的になるために生じたものであると考え た。以上の結果より振幅誤差にモデリングの影響が強く現れるため、モデリングのためのセル数は振幅 の精度が所望の値となる様にすべきと考える。 最後に波長に対して広大な解析領域を持つ散乱問題として 3 × 3 のアレー状に整列した誘電体に関 する散乱問題の数値解析を行う。CIP-BS 法と FDTD 法を用いて解析を行い散乱場の精度を比較した。 CIP-BS法の場合、∆t の減少と共に振幅誤差は指数的に減少した後、収束した。より詳細にモデリング を行っても収束後の誤差は変化することはなかったため、直径あたり 20 セル程度の離散化で十分である と考えられる。また同じ位相精度で比較したとき、高精度な領域ならば CG 法を用いた CIP-BS 法が最 も省メモリであった。計算時間においては FDTD 法が高速であったが、高精度化に伴う計算時間の増加 割合は CIP-BS 法が少なかったためより広大な解析領域に適すると考えられる。謝辞
本研究を進めるにあたりご指導くださった安藤芳晃准教授を初め、多くの助言を頂いた研究室の皆様 に感謝申し上げます。
参考文献
[1] H. Takawaki, A. Nishiguchi and T. Yabe, “Cubic interpolated pseudo-particle method (CIP) for solving hyperbolic-type equations ”, Jour. Comp Phys. vol 61,pp.261-268,1985.
[2] 矢部孝, 尾形陽一, 内海隆行, CIP 法原子から宇宙までを解くマルチスケール解法, 森北出版, 2007 [3] T. Utsumi,T. Yabe, J. Koga, T. Aoki and A. Sekine, “Accurate Basis Set by the CIP Method for
the Solutionsof the Schrodinger Equation ”, Comp Phys. Commun.,Vol.157, pp.121-138, 2004 [4] T. Utsumi, T. Yabe, T. Aoki, J. Koga, and M. Yamagiwa, “Solution of hyperbolic equations with
the CIP-BS method ”JSME Int. J. B,vol.47, pp.768776, 2004
[5] Y. Ando, A. Noba and S. Murakoshi “CIP-BS method for solving Maxwell’s equations ”, ETMS 2010 Digest, pp.111-113, 2010.
[6] Y.Ando and Y.Takahashi, “CIP Basis Set method for electromagnetic simulation ”, IEICE Trans. Electron., vol. E97-C, no.1, pp.26-32, 2014.
[7] G.H. Golub, C.F.V. Loan, Matrix Computations, Third Ed., The John Hopkins University Press.,London, 1996.
[8] A. Taflove and S. C. Hagness, Computational electrodynamics : the finite-difference time-domain method,Second ed. , Artech House, Boston, 2000
[9] 高橋 祐介, 安藤 芳晃, “CIP-BS 法における各種個別技術の開発, 電気通信大学”, 電気通信大学, 修 士論文, 2011 年
