電磁界における移動境界問題の数値解析法
黒田
道子
†a)Numerical Technique for the Analysis of Electromagnetic Field with Moving
Boundaries
Michiko KURODA
†a)あらまし 移動体通信や可変機能をもつデバイスの開発など，近年電磁界数値解析の分野において，移動する 物体や移動する波源などの任意運動，任意形状に対する数値解析の必要性が強く求められている．ここでは，複 数の格子を用いる重合格子法をローレンツ変換も考慮してFDTD 法に導入し，電磁界移動境界問題の数値解析 法を提案している．この手法の精度を検証し，物体が移動する場合，波源が移動する場合についての数値計算結 果を理論値と比較し，良好な精度が得られたので報告する． キーワード 移動境界，重合格子法，ローレンツ変換，FDTD 法，移動体通信
1.
ま え が き
近年，移動体通信や可変機能をもつデバイスの開発 が活発に行われているが，これらの設計，解析のため には電磁界における移動境界問題の数値解析法が強く 求められている．電磁界解析分野における移動境界問 題は，これまでは理論解析法[1]∼[3]やFDTD法に 移動位置だけを導入する方法[4]が提案されているが， 単純な運動や一定方向への運動に限られていた．この ため，任意形状，任意運動に対しては限界があった． ここでは，FDTD法を用いた移動する物体や境界を もつ場合，また波源が移動する場合の数値解析法につ いて述べる． 移動境界問題を電磁界数値解析に導入するには困難 な問題がある．それは電磁波の速度が光速であるとい う点である．このため，低速運動問題を取り扱う場合 は，時間離散間隔を光速に合わせているため必要な移 動距離に到達するまでに計算時間がかかること，高速 運動問題を考える場合はローレンツ変換を考慮しなけ ればならないことである．筆者らは，これらの問題点を 解決し，FDTD法を用いて移動境界問題を直接数値解 析できる方法についてこれまで提案してきた[5]∼[9]． †東京工科大学，八王子市Tokyo University of Technology, 1404–1 Katakura, Hachioji-shi, 192–0982 Japan a) E-mail: [email protected] ここでは，重合格子法[10]を高速運動の場合にも適 用できるようローレンツ変換をFDTD法に導入した 数値解析法について論述する．ここで用いる重合格子 法は主格子上に副格子を重ね，副格子を移動させ，補 間法を用いて各格子上の物理量のやり取りを行うこと で移動する境界の数値解析を可能とする方法である． 本論文では，この手法について述べた後，電磁界数 値解析において高速移動体の場合に問題となるローレ ンツ変換を高速の場合に適用できるように新たに考案 した重合格子法，ローレンツ変換とFDTD法を組み 合わせた数値計算アルゴリズムについて説明する．こ の手法を用いた数値解析法の精度を検証し，数値結果 と理論値との比較を行う．また，主な例として，物体 が移動する場合と波源が移動する場合についての数値 例を示す．
2.
重合格子法とローレンツ変換
重合格子法とは，複数の格子上でお互いに物理量を やり取りしながら解析を行う手法である[10]．ローレ ンツ変換を考慮した電磁界における移動境界問題では， 静止系と運動系がある．これに重合格子法を適用す るには，主格子(Main mesh)を静止系に副格子(Sub mesh)を運動系として，運動系である副格子を移動さ せ格子間での物理量のやり取りを行うことにより移動 境界問題を電磁界において直接数値解析できるように なる．図 1 相対運動する二つの慣性系 Fig. 1 Two inertial frames in relative motion.
2. 1 電磁界のローレンツ変換 ローレンツ変換とは二つの慣性系間の座標（時間座 標と空間座標）を結びつけるマクスウェルの方程式を 満足する座標変換方式である[1]．図1で示すように運 動系K が静止系Kに対して速度ベクトルvで任意 方向へ移動しているとき，空間因子と時間因子は，以 下のベクトルの表示式で表すことができる． dr= γ(dr − vdt) + (γ − 1)v × v × dr v2 (1) dt= γ dt −v • dr c2 (2) ただし，cは光速である． γ = 1 1 − v2/c2 = 1 1 − β2 ここで，dr，dtはK系，dr，dt はK系における 空間因子と時間因子である．二つの慣性系が上記の関 係にあるとき，その電磁界変換式は式(3)∼(6)で与え られる． E= E + γ(E⊥+ v × B) (3) H= H + γ(H⊥− v × D) (4) B= B + γ B⊥− v × E c2 (5) D= D + γ D⊥+v × H c2 (6) ここで，E，H，D，Bは運動方向に対して平行 な成分であり，E⊥，H⊥，D⊥，B⊥は運動方向に対 して垂直な成分である． 2. 2 重合格子法 重合格子法とは，複数の格子間での物理量をお互い にやり取りしながら解析を行う手法である．移動境界 問題に対しては，静止系を主格子(Main mesh)，運動 図 2 線形補間モデル Fig. 2 Interpolation model.
系を副格子(Sub mesh)とし，副格子を移動させ，格 子間で物理量のやりとりを行うことにより移動境界問 題を解決できる．副格子は複数導入することも可能で ある． 各格子間での物理量のやりとりは補間法を用いて行 う．補間法には様々な手法があるが，ここではFDTD 法と重合格子法を組み合わせて使用し，空間離散間隔 や時間離散間隔が微小であるため，線形補間を用いる． 図2は，重合格子法の2次元モデルを示す．破線が主 格子，実線が副格子である．主格子上の物理量φP は， 副格子上の4点の物理量φ1，φ2，φ3，φ4 より線形補 間を用いて， φP =(φ1x2+φ2x1)y2+(φ3x2+φ4x1)y1 (x1+x2)(y1+y2) (7) で求められる． 2. 3 FDTD法へのローレンツ変換の導入 つぎに，FDTD法へローレン変換を導入する手法 について述べる．FDTD法は空間格子，時間格子とも に格子間隔は一定である．ところが，ローレンツ変換 を用いることで，式(1)，(2)のように空間離散間隔も 時間離散間隔も変化してしまう．このため，空間離散 間隔，時間離散間隔を一定に保ちながらローレンツ変 換を導入することが求められる．この問題を解決する 計算手法を考案したので以下に示す． ロ ー レ ン ツ 変 換 を 導 入 し ，重 合 格 子 法 を 用 い た FDTD法の計算アルゴリズムを図3に示す．(a)，(d) を静止系Kの主格子，(b)，(c)を運動系Kの副格 子としている．図3においてローレンツ変換を用い て線形補間を行うとき，式(1)，(2)により時間離散 間隔が変化してしまうので，FDTD法のアルゴリズ ムを保つために，時間離散間隔の整合を取ることが 重要になってくる．運動系K では，空間のみでなく
図 3 時間と空間のアルゴリズム
Fig. 3 Time and space algorithm for FDTD and Lorentz transformation by overset grid gen-eration method. 時間の進み方が変化するので，静止系Kでの時間離 散間隔をΔtとしたとき K 系での時間離散間隔は Δt=1 − β2Δtとなる．ここでは，x方向へ速度 vで移動する場合について，図3のアルゴリズムで説 明する． 図3 (a)での電界Eをローレンツ変換により(b)の K系に変換し，(b)から(c)へFDTD法で計算を行っ た後，(c)から(d)へローレンツ逆変換を行う．K系 での電界と磁界の空間離散間隔をΔx/2，Δy/2とす る．時間，空間アルゴリズムに従って電界Eをローレ ンツ変換すると(d)でΔt進んだEの値を得ることが できる．同様の変換をHに対しても行うが，ローレ ンツ変換を用いているため，(d)でΔt進んだHの値 を得ることができない．(d)でΔt進んだHの値を得 るため，K系上に座標変換する磁界Hはあらかじめ 空間位置Xと時間T となるよう空間も時間も補間を 行った後，(a)から(b)へローレンツ変換を行う．この 補間を行うことにより，(d)において電界E，磁界H ともにΔt進んだ物理量を得ることができ，FDTD法 にローレンツ変換を導入することができるようになる．
3.
重合格子法を用いた
FDTD
法の精度
提案した手法の精度について検証することは，この 手法の有用性を示す上で重要である．ここでは，主格 子と副格子の空間離散間隔の比と計算精度について検 証する． 精度を検証するためのモデルを図4に示す．正方 図 4 2次元モデル Fig. 4 2D numerical model.形の主格子の上に正方形の副格子を重ね，両方の中 心はそろえている．主格子長はLm× Lm，副格子長 はLs× Ls，Im(m, n)，Is(i, j)は各々主格子の格子 数，副格子の格子数を表す．ここで，主格子の空間離 散間隔と副格子の空間離散間隔の比を以下の式R (sub rate)で表す． R = Ls Is Lm Im (8) 精度を検証するため，主格子数は400 × 400で固定 し，副格子数を36 × 36，42 × 42，50 × 50，62 × 62， 84 × 84，125 × 125 と変化させて，主格子と副格子 の空間離散間隔比Rに対して精度を検証していく． 主格子の空間離散間隔は，Δxm = Δym= λ/60（λ は波長）で一定である．ここでは，副格子長 Ls は 変化せず，副格子数が変化するため，副格子の空間 離散間隔 Δxs，Δys は変化する．Lm = 0.02(m)， Ls= 2.5 × 10−3(m)，Δt = 5.0 × 10−14(sec)として 計算を行う．副格子数を36×36から250×250まで変 化させるため，それに対応して副格子の空間離散間隔 が変化する．FDTD法ではCourant Conditionを満 足する必要があるため，副格子数が最大の250 × 250， つまり副格子の空間離散間隔が最小の値を取るときに， この条件を満足するようにΔtを選んでいる． 主格子上のm = 25から100GHzの平面波を入力 し，主格子上に副格子を重ねた場合と主格子のみの場 合について，以下の式を用いて誤差を検証する． Relative error =Ey− E0 E0 (9) ここで，Eyは主格子の上に副格子を重ねた場合の電 界，E0 は主格子のみの場合の電界を示す． 副格子の空間離散間隔と主格子の空間離散間隔比R
(a)
(b)
(c)
図 5 (a)R = 0.6 の場合の電界分布．(b) R = 1 の場合 の電界．(c)R = 1.25 の場合の電界分布．
Fig. 5 (a) Electric field,R = 0.6. (b) Electric field, R = 1.0. (c) Electric field, R = 1.25.
図 6 副格子と主格子空間離散間隔比R と誤差の関係
Fig. 6 Absolute error of the amplitude ofEy forR.
(subrate)と計算精度についての結果を，図5に示す． 図5 (a)は R = 0.6 の場合，(b)はR = 1 の場合， (c)はR = 1.25の場合の電界分布である．これらの 図より，(a)，(b)の場合には電界の乱れは見られない が，R = 1.25になると誤差が出てきているのがわか る．また，図6にRと誤差の関係を示す．この結果よ り，Rが1を超えると急激に誤差が大きくなることが わかる．これは，線形補間において，副格子の空間離 散間隔が大きくなることにより誤差が生じたものと思 われる．R ≤ 1の場合は，副格子の空間離散間隔が小 さい値を取るため，誤差が少なく良好な精度で計算が できることが分かった． 以上の結果より，副格子の空間離散間隔は小さい方 が精度がよく，副格子と主格子の空間離散間隔の比R は1以下がよい精度が得られることが分かった．
4.
重合格子法を用いた移動する物体
ここでは，移動する物体に対する電磁波について， 数値計算を行う．移動する誘電体に対する2次元モデ ルを図7に示す．方形誘電体が入射波に対して，遠ざ かる場合と近づく場合について，数値計算を行う．主 格子の上に重ねた副格子を移動する誘電体物体に設定 する．主格子の格子数(m, n) を400 × 400，副格子 数(i, j) を250 × 250とし，主格子の空間離散間隔 をΔxm= Δym = λ/60，副格子の空間離散間隔を Δxs = Δys = λ/300 とした．ここで，λ は波長で ある．時間離散間隔は，Δt = λ × 10 −4 2c とし，主格 子上のm = 25から誘電体に向けて100 GHzの平面 図 7 入射波に対して移動する誘電体Fig. 7 Numerical model of the dielectric body moving for the incident wave.
(a) (b) (c) 図 8 (a)物体が静止している場合の電界分布v = 0．(b) 入 射波に対して近づく場合の電界分布 (v = 0.1c)．(c) 入 射波に対して遠ざかる場合の電界分布 (v = 0.1c)．
Fig. 8 (a) EM field scattering from static dielectric body withv = 0. (b) EM field scattering from dielectric
body approaching to the incident wave (φ = π) with v = 0.1c. (c) EM field scattering from di-electric body moving away from the incident wave, (φ = 2π) with v = 0.1c. 波を入射した．副格子は，移動する誘電体物体とし， εr= 3.5の比誘電率の場合の電界分布を示す．入射波 へ向かって進む場合(φ = π)と遠ざかる場合(φ = 2π) について，移動速度を v = 0.1c(m/s)として電界分 布を求めた結果を図8に示す．(a)は静止している場 合，(b)は入射波に対して近づいていく場合，(c)は入 射波から遠ざかる場合についての電界分布を示す．こ れらの結果より，反射波，透過波に違いがみられ，運 動の影響が出ていることがよくわかる． この結果を理論値と比較してみる．誘電体からの反 射波や透過波の振幅比は理論式(10)，(11)で表され 表 1 入射波に対する透過波の振幅比 Table 1 Ratio of the amplitude for transmitted wave
to the incident wave.
表 2 入射波に対する反射波の振幅比 Table 2 Ratio of the amplitude for reflected wave to
the incident wave.
る．[1] Ampt Ampi = 2(1 + v/c√εr) (1 +√εr)(1 + v/c) (10) Ampr Ampi = 1 − v/c 1 + v/c· 1 −√εr 1 +√εr (11) ここで，Ampt は透過波の振幅，Ampi は入射波の振 幅，Ampr は反射波の振幅を表す． 表1，2に理論値と数値解析結果の振幅比の比較を 示す．この結果より，本手法による数値結果と式(10)， (11)の結果は，よく一致していることがわかり，物体 が移動する場合の電磁波の数値解析結果の精度を確認 することができた．表1，2の結果において，近づく 場合は，主格子と副格子の位置的距離が近くなるため， 線形補間の精度が上がり，遠ざかる場合は，主格子と 副格子の位置的距離が遠くなるため，近づく場合に比 べて誤差が大きくなったと考えられる．
5.
重合格子法を用いた移動する波源
これまでは，移動する物体への電磁波の解析法につ て述べてきたが，ここでは波源が移動する場合につい ての数値計算を行う[11]．図9にモデル図を示す．副図 9 計算精度検証モデル
Fig. 9 Numerical model to verify the accuracy.
図 10 R (sub rate) の違いによる相対誤差と標準偏差 Fig. 10 Sub rateR versus relative error and standard
deviation. 格子上に波源を置き，副格子を移動させることで，波 源の移動ができるようにしている．波源が移動する場 合の精度を検証するため，この図を用いて，主格子と 副格子の空間離散間隔の比を表す式(8)を用いてR を変化させ，以下の式(12)を用いて，相対誤差及び 標準偏差を求める．相対誤差は式(13)に示す．主格 子をΔxm= Δym= λ/40，副格子間隔をλ/40から λ/340まで変化させ，Δt = λ 500c，入射波750MHz， 移動速度v = 0.01cで計算を行った． f =1 + v/c 1 − v/cf0 (12) Relative error
=theoretical results−numerical results theoretical results (13) 図10にその結果を示す．横軸に主格子と副格子の 比R (sub rate)，縦軸に相対誤差と標準偏差を示す． 式(12)に示す理論値と計算値の相対誤差と標準偏差 を用いて比較する．波源が移動しているため，副格子 の空間離散間隔Δxを小さくすると，格子数が増し， 図 11 十字路形モデル
Fig. 11 Numerical model having a shape of cross section.
(a)
(b)
図 12 (a)静止波源の場合の電界分布．(b) 移動波源の場 合の電界分布．
Fig. 12 (a) Electric Field with stationary input wave. (b) Electric Field with moving input wave.
線形補間の回数が増え，その結果，相対誤差が大きく なると考えられる．標準偏差と誤差を検討した結果， 標準偏差が一番小さく，相対誤差は丸めの範囲でほぼ 同等のR = 1の場合を最適とした． 次に，波源が移動する場合について，図11の交差 点における車両の移動を想定したモデル図に基づき数 値解析を行った．十字形に交差した自由空間の周りを 比誘電率εr= 7.0，導電率σ = 0.0389[S /m]の誘電
体で囲み，点AからBまで波源を移動させる．点C， Pは観測点とする．静止波源は，点Bにおいて静止 状態で電磁波を放射するものとしている．解析条件と して，全長50λ，自由空間幅10λ，主格子の空間離散 間隔と副格子の空間離散間隔の比をR = 1 とするた め，空間離散間隔Δx = Δy = λ/40，時間離散間隔 Δt = λ/80とし，入射波750MHz，波源の移動速度は v = 0.01cとした．主格子の格子数は，2000 × 2000， 副格子の格子数は，50 × 50とした．点Bで静止の場 合の電界分布と点Aから点Bへ移動する場合の電界 分布を図12に示す．これは，移動波源が点Aから出 発し点Bに到着した時間で表示している．これらの図 を比較すると，十字路中央の上下に，波源の移動の影 響が出ているのがよくわかる．
6.
む す び
移動体通信の高速化や可変機能をもつ光・マイクロ 波デバイスの発展に伴い，電磁界数値解析の分野にお いて移動境界問題の数値解析法が強く求められてきて いる．本論文では，重合格子法とFDTD法を組み合わ せた電磁界数値解析法について論述した．特に，ロー レンツ変換を電磁界数値解析に導入する手法について 詳しい説明をした． 電磁界数値解析分野への数値解析法の導入において， 困難な点は電磁界が光速であることである．FDTD法 においては，Courant条件を満足するため，時間離散 間隔Δtは小さな値を取り，また計算機容量が限られ ているため，ドップラー効果を明確に現すためには， ここで提案した数値例においても，入射周波数を高く とらなければならなかったこと，物体や波源の移動速 度を高速に設定しなければならなかったことなどまだ まだ，解決しなければならない問題点がある． しかし，これまで厳密解法が主であった電磁界数値 解析分野に移動境界数値解析法を提案することができ たことは，任意形状，任意運動，更には3次元問題を 解析することが可能になり，大きな進歩であると考え る．今後は，並列処理や高速大容量の計算機などを用 いた計算処理法によりこの手法はより一般的になり利 用範囲も拡大していくことを期待する． 文 献[1] V. Bladel, Relativity and Engineering, Springerverlag, Berlin, 1984.
[2] T. Shiozawa, “Electromagnetic scattering by a mov-ing small particle,” J. Appl. Phys., vol.39, no.7,
pp.2933–2998, June 1968.
[3] K. Tanaka, “Scattering of electromagnetic waves by a rotating perfectly conducting cylinder with arbi-trary cross section: Point-matching method,” IEEE Trans. Antennas. Propag., vol.28, no.6, pp.796–803, Nov. 1980.
[4] F. Harfoush, A. Taflove, and G.A. Kriegamann, “A numerical technique for analyzing electromagnetic wave scattering from moving surfaces in one and two dimensions,” IEEE Trans. Antennas. Propag., vol.37, no.1, pp.55–63, Jan. 1989.
[5] M. Kuroda and S. Kuroda, “The body fitted grid generation with moving boundary and its application for optical phase modulation,” IEICE Trans Electron, vol.E76-C, no.3, pp.480–485, March 1993.
[6] H. Iwamatsu and M. Kuroda, “Over set grid gen-eration method coupled with FDTD method while considering the Doppler effect,” IEEJ Trans. FM, vol.129, no.10, pp.699–703, Oct. 2009.
[7] S. Sahrani and M. Kuroda, “A novel approach for the analysis of electromagnetic field with rotating body,” Aces Journal, vol.26, no.8, pp.651–659, Aug. 2011. [8] S. Sahrani and M. Kuroda, “FDTD analysis with
overset grid generation method for rotating body and evaluation of its accuracy,” IEICE Trans. Electron, vol.E96-C, no.1, pp.35–41, Jan. 2013.
[9] S. Sahrani and M. Kuroda, “Numerical technique for the analysis of electromagnetic field by moving dielec-tric body,” IEEJ Trans. FM, vol.133, no.5, pp.255– 259, May 2013.
[10] J.F. Tjompson, Z.U.A. Warsi, and C.W. Mastin, Nu-merical Grid Generation, North-Holland, New York 1985.
[11] T. Akata, S. Sahrani, and M. Kuroda, “Numerical analysis of the EM field from a moving source and the application for a moving vehicle,” EuMC2014, Oct. 2014. （平成 28 年 9 月 27 日受付，11 月 29 日再受付， 29年 4 月 12 日公開） 黒田 道子 （正員：フェロー） 1973年静岡大学工学部電気工学科卒， 1975年早稲田大学大学院理工学研究科修 士課程了，1978 年早稲田大学大学院理工 学研究科博士課程了．工学博士．1986 年∼ 1987年オハイオ州立大学客員研究員．1990 年東京工科大学工学部情報工学科助教授， 1998年同大学教授，2004 年同大コンピュータサイエンス学部 教授，2011 年∼2013 年コンピュータサイエンス学部長，2014 年東京工科大学名誉教授．電磁界移動境界問題の数値解析に取 り組む．電子情報通信学会フェロー．IEEE Senior Member．
