超広帯域電磁界解析のための周波数依存性FDTD法

招待論文

超広帯域電磁界解析のための周波数依存性

FDTD

法

チャカロタイ

ジェドヴィスノプ

†a)

_{和氣加奈子}

†

_渡辺

_聡一

†

陳

強

††

澤谷

邦男

†††

Frequency-Dependent FDTD Method for Ultra-Wideband Electromagnetic

Analyses

Jerdvisanop CHAKAROTHAI

†a)

, Kanako WAKE

†

, Soichi WATANABE

†

,

Qiang CHEN

††

, and Kunio SAWAYA

†††

あらまし 本論文では，著者らが提案した高速逆ラプラス変換 (Fast Inverse Laplace Transform; FILT) 及び Prony 法を利用した周波数依存型時間領域有限差分 (Frequency-Dependent Finite-Difference Time-Domain; (FD)2_{TD) 法について一般定式化を行い，広帯域電磁界による散乱問題及びアンテナ問題に対して提案手法を適} 用し，その有効性を示す．まず Jonscher-Raicu 分散性を有する均一誘電体球に対して散乱解析を行い，Mie の 解析解との比較により，提案手法の妥当性を確認する．次に，均一誘電直方体近傍における広帯域アンテナの反 射係数を求め，生体組織近傍における広帯域アンテナの設計に対して提案手法の有効性及び有用性を示す．最後 に，人体近傍において設計した広帯域アンテナの性能評価について有用性及び計算効率性を明らかにした． キーワード 時間領域有限差分法，広帯域，周波数依存性，散乱断面積，広帯域アンテナ，人体モデル

1.

ま え が き

時間領域における電磁界解析は1966年に発表され たYeeの先駆的な研究により著しく進歩してきた[1]． Yeeにより提案されたアルゴリズムは等方性媒質に 対して，時間領域でマクスウェル方程式を解くもの であり，時間領域有限差分(Finite-difference time-domain; FDTD)法と呼ばれる[2]．周波数領域解法 であるモーメント法[3]や有限要素法[4]に対して， FDTD法の最も有利な点の一つは時間領域で，一度 の解析で，広帯域の解が得られることである．FDTD 法のロバスト性，多用途性及び実装の容易さ等により， これまで高速電子回路解析[5]，非線形現象解析[6]， †_{国立研究開発法人情報通信研究機構，小金井市}

National Institute of Information and Communications Technology, 4–2–1 Nukui-Kitamachi, Koganei-shi, 184–8795 Japan

††_{東北大学大学院工学研究科通信工学科専攻，仙台市}

Tohoku University, 6–6–05 Aramaki Aza Aoba, Aoba-ku, Sendai-shi, 980–8579 Japan

†††_{東北大学産学連携機構イノベーション戦略推進センター，仙台市}

Head office of Enterprise Partnerships, Tohoku University, 468–1 Aoba, Aramaki-Aza, Aoba-ku, Sendai-shi, 980–0845 Japan a) E-mail: jerd@nict.go.jp フォトニック結晶[7]，電磁界ばく露評価[8]∼[10]等 の様々な電磁界問題に対して適用されてきた．しかし ながら，Yeeによるアルゴリズムは広帯域にわたる周 波数に依存する電磁特性を有する分散性媒質に対して モデル化できない欠点がある．また特定の周波数にお ける一定な複素誘電率または複素透磁率を用いた解析 では，比較的狭帯域でしか解が得られず，応用範囲が 限られていた． FDTD法を分散性媒質の解析へ拡張するための手法 として，回帰的畳み込み(Recursive convolution; RC) 法[11], [12]，補助微分方程式(Auxiliary differential equation; ADE)法[13]∼[15]及びZ変換法[16]∼[18] が挙げられる．これらはDebye分散またはLorentz 分散モデルによって表される周波数依存性媒質の解 析に対して有効であるが，いずれの手法も生体組織 のように通常複素誘電率がCole-Cole分散モデルに よって表される媒質への適用が困難である．これは Cole-Cole分散式では，時間領域において分数次数微 分の定式化が必要になるためである．この問題を克服 するために，これまで様々な手法が提案された．例え ば，Cole-Cole分散を複数のDebye分散項の和によっ て表す手法[19], [20]や双次Z 変換[21]等が用いられ

た．複数のDebye分散項における緩和時間等のパラ メータ決定には，非線形最適化手法[19], [20]やPade 近似[22]等が用いられた．またRiemann-Liouville理 論を用い，より一般化されたHavriliak-Negami分散 式やJonscher-Raicu分散式のFDTD定式化もある が，定式化が複雑である[23]∼[25]．一方，これまで 著者らが一般的な周波数分散性媒質に対する解析の ために，高速逆ラプラス変換(Fast inverse Laplace transform; FILT) [26]∼[28]及びProny法を組み合

わせた手法を提案した[29]．FILTによって周波数領

域における複素誘電率を時間領域インパルス応答へ

変換した後，Prony法によってZ 領域での表現を求

める．これによりFDTD法における定式化が簡単で

かつ容易に実装が可能となる．著者らはこれまで複素 比誘電率がDebye，Davidson-Cole，Cole-Cole及び

Havriliak-Negamiの分散モデルによって表される周 波数分散性媒質に対して提案手法が適用可能であるこ とを示したとともに，広帯域パルス電磁界に対するば く露評価問題に対して提案手法を適用し，有効性を示 した[29], [30]． 本論文では，これまで著者らが提案した周波数依存 性FDTD法をより一般化されたJonscher-Raicu分 散式[31]への適用可能性及び妥当性について検討する 他，磁気分散性を有する磁性体フェライトへ適用し， 人体近傍における広帯域アンテナの性能評価に対して 提案手法の有用性及び計算効率を示す．提案手法の利 点としては，分数次数の微分方程式を時間領域で定式 化する必要がなく，FILT及びProny法を用いること で，直接FDTD計算内の更新式における係数を導出 することができ，様々な分散モデルへの適用が可能で ある点と文献[24], [25]に比べてFDTD実装が容易な 点が挙げられる．また提案手法は広帯域電磁界解析に おいて，他の手法と比べてもロバスト性，多用途性， 効率性の観点からも有利である．本論文の構成は以下 のようになる．まず2.では，提案手法における定式 化及び更新式中の係数の決定法について述べる．次に 3.及び4.では，それぞれ広帯域パルス電磁界による 散乱解析及び広帯域アンテナ問題について様々な計算 結果を示しながら，本提案手法の有用性及び効率性を 示す．

2.

提 案 手 法

本章では，著者らが提案した周波数依存性FDTD 法について一般定式化を行い，Jonscher-Raicu分散 式に対する係数の決定について述べる． 2. 1 マクスウェル方程式のFDTD定式化 媒質内におけるマクスウェル方程式は時間領域にお いて，一般的に次式によって表すことができる． ∇ × H = ∂D ∂t, ∇ × E = − ∂B ∂t. (1) ここで，E，H，D及びBはそれぞれ電界[V/m]，磁 界[A/m]，電束密度[C/m2_] 及び磁束密度[T]を表し ており，全て空間及び時間の関数である．電界Eと電 束密度D及び磁界H と磁束密度Bは構成方程式に よって関連付けられており，周波数領域において一般 的に， D(ω) = ε0ε˜r(ω)E(ω) = ε0 ⎡ ⎣εr∞+ σ jωε0 + NE  g=1 χEg(ω) ⎤ ⎦E(ω), (2) B(ω) = μ0μ˜r(ω)H(ω) = μ0 ⎡ ⎣μr∞+ σ ∗ jωμ0 + NH  h=1 χHh(ω) ⎤ ⎦H(ω). (3) と表される．なお，ここでは，媒質は異方性をもたな いとする．上式の中で，ε0，μ0，εr∞，μr∞，σ及び σ∗ はそれぞれ真空中の誘電率[F/m]，真空中の透磁 率[H/m]，無限周波数における比誘電率及び比透磁 率，導電率[S/m]及び磁気導電率[Wb/(A·m)]であ る．NE 及びNH はそれぞれ電気感受率χE g と磁気感 受率χHh の項数である．χEg 及びχHh は媒質によって 異なるモデルが用いられる．例えば，生体組織の場合， 電気感受率をCole-Coleモデルで表されるのが実験値 とよく一致していると報告されており[32], [33]，光領 域における金属膜の電気感受率はドュルーデ・ローレ ンツ分散式によって表せる[34]．フェライト等の磁性 材料に対しては，磁気感受率がローレンツ分散モデル によって表される[35]．式(2)及び(3)より複素比誘 電率ε˜r 及び複素比透磁率μ˜r が周波数に対する関数 となっているため，時間領域においてそれぞれ電束密 度D及び磁束密度Bを求めるためには，非常に時間 のかかる畳み込み積分を計算しなければならない．一 方，時間領域における畳み込み積分はz 領域では，単 なる掛け算となるため，これを利用して複素比誘電率 ˜ εr 及び複素比透磁率μ˜rをz 領域で表現できれば，畳

み込み積分に必要な計算コストを大幅に削減できる． 本研究では，これまで著者らが提案した高速逆ラプラ ス変換及びProny法を用いた手法によって複素比誘 電率ε˜r 及び複素透磁率μ˜r を容易にz領域に変換す ることができることを示した[29]．詳細な導出手順は 文献[29]を参照されたい．なお，定式化を簡単にする ために，NE= NH= 1と仮定し，z変換を適用する と，ε˜r 及びμ˜r は次式のように表せる． ˜ εr(z) = εr∞+1 + z −1 1− z−1 σΔt 2ε0 +XE1(z), (4) ˜ μr(z) = μr∞+1 + z −1 1− z−1 σ∗Δt 2μ0 +XH1 (z). (5) ここで，Δtは時間ステップ間隔であり，XE1 及びXH1 は次式のように表される． X∗ 1(z) = N_l∗  l=1 A∗l 1− p∗_lz−1 + N_k∗  k=1 Bk∗− Ck∗z−1 1− r∗_kz−1_{+ q}∗ kz−2 (6) Nl∗及びNk∗（∗ = E or H）は実極の数及び複素極の 対の数を表す．また係数A∗l，p∗l，Bk∗，Ck∗，r∗k及び qk∗は文献[29]に示す手順で計算でき，全て実数とし て求まるため，FDTD法における更新式の中に直接 代入して計算できる． FDTD法における電界及び磁界の更新式を導出する ために，式(4)及び(5)をz 領域で表現された式(2) 及び(3)に代入し，最終的には，電界及び磁界につい て以下の更新式が得られる． En = 1 GE 1 ⎡ ⎣GE 2E n−1₊Δt ε0∇ × H n−1 2 + N_lE  l=1  1− pEl  Pn−1l + N_kE  k=1  1− rEk  Qn−1 k + q E kQ n−2 k ⎤ ⎦, (7) Hn+1₂ = 1 GH 1 ⎡ ⎣GH 2 H n−1 2−Δt μ0∇ × E n + N_lH  l=1  1− pHl  Rn−1 2 l + N_kH  k=1  1− rHk  Sn−1₂ k + q H kS n−3₂ k ⎤ ⎦. (8) ただし， GE1 = εr∞+σΔt 2ε0 + NE_l  l=1 AEl + N_kE  k=1 BkE, (9) GH1 = μr∞+σ ∗_Δt 2μ0 + NH_l  l=1 AHl + N_kH  k=1 BkH, (10) GE2 = εr∞− σΔt 2ε0 + N_kE  k=1 CkE, (11) GH2 = μr∞− σ ∗_Δt 2μ0 + N_kH  k=1 CkH. (12) また補助変数Pnl，Qnk，R n+1₂ l 及びS n+1₂ k は以下の 更新式によって計算される． Pn l = p E lP n−1 l + A E l E n , (13) Qn k= r E kQ n−1 k − q E kQ n−2 k + B E kE n_{− C}E kE n−1 , (14) Rn+1₂ l = p H l R n−1₂ l + A H l H n+1₂ , (15) Sn+1₂ k = r H kS n−1 2 k − q H kS n−3 2 k + B H kH n+1₂ − CH k Hn− 1 2_, (16) ここで，文献[29], [30]における定式化と異なり，電束 密度Dn及び導電率に関する補助変数In を含まない 形で定式化を行っている点に注意されたい．これによ り補助変数Dn 及びIn のメモリ使用量を削減でき， かつ計算効率も向上させることができる．式(7)–(16) の計算フローチャートを図1に示す．まずFILT及び Prony法によって求めた係数A∗l，p∗l，Bk∗，Ck∗，r∗k， qk∗の値を読み込んだ後，式(9)–(12)により更新式の 中に用いられる各係数GE1，GH1 ，GE2，GH2 を計算し， FDTD法の計算ループに入る．計算ループの中では， まず式(7)によって電界Enを計算した後，式(13)及 び(14)により補助変数Pnl 及びQnk をそれぞれ更新 する．それから電界について完全整合層(PML)の計 算を行う．磁界についても同様に式(8)磁界Hn+12 を 計算した後，式(15)及び(16)により補助変数Sn+12 及びRn+1_{2 をそれぞれ更新し，最後に磁界について} PMLを計算する．そして，計算ステップ数NT に達

図 1 提案手法の計算フローチャート するまで，計算ループを繰り返す． 本 提 案 手 法 の 使 用 メ モ リ 量 に つ い て 以 下 に 述 べ る．三次元モデルの場合，解析空間の大きさ NA = Nx× Ny× Nz とすると，En 及びHnに必要なメ モリ量はそれぞれ x，y，及び z 成分の変数に要す る 6NA である．補助変数の Pnl，Qnk，R n+1/2 l 及 びSn+1/2k に必要なメモリ量は，N B を媒質が占め るFDTDセル数とすると，3× NB× (NlE+ 2NkE+ NlH+ 2N H k )である．したがって，全使用メモリ量は 6NA+ 3× NB× (NlE+ 2N E k + N H l + 2N H k )であ る．ここでそれぞれの媒質には分散性を考慮する項数 NE l ，NkE，NlH 及びNkH が違う値をもっていること にも注意されたい． 本手法の最大の利点として，FILT及びProny法を 利用することで，Cole-Cole分散式で示されるような 時間領域における分数次数微分の定式化が不要であり， かつ様々な分散性媒質に対して適用可能である．また 従来のFDTD法への拡張が容易であることが挙げら れる．一方，提案手法の欠点として，複数の数値手法 を用いているため，解の不安定性を引き起こす可能性 があるが，文献[29]に示す安定性解析を行った結果， 本研究で検討している全ての解析モデルで選択した時 間ステップ間隔においては，安定的な解を与えること が分かった． 2. 2 更新式における係数の決定 更新式における係数A∗l，p∗l，B∗k，Ck∗，r∗k，q∗kは FDTD計算を始める前に，あらかじめ求めておく必要 がある．本研究で解析に用いている媒質は主に Cole-Cole分散モデルによって表現できる生体組織である が，より一般的な複素誘電率の表現として， Jonscher-Raicuモデルがある[31]．著者らはこれまでDebye，

Cole-Cole，Davidson-Cole，Havriliak-Negamiの分

散式について提案手法の妥当性を検討してきた[29]．こ こで，Jonscher-Raicu分散モデルへの適用可能性につ いて初めて検討を行うこととなる．また透磁率に分散 性を有するフェライトのモデル化の妥当性検討も行う． まずJonscher-Raicu分散式は周波数領域において， 次式のように表すことができる． ˜ εr(ω) = εr∞+ σ jωε0 + Δχ [(jωτa)1−α+ (jωτb)1−β]γ. (17) ここで，Δχはそれぞれ比誘電率の変化分，τa とτb は緩和時間，α，β，γは0∼1の値をもつ実定数であ る．なお，文献[31]においては，τa及びτbは同じ値 をもっているが，本研究では，それぞれ別々の値であ ることに注意されたい．ここでは，誘電率に対しての み分散性を有する媒質であると仮定するが，透磁率の 場合にも同様に適用することができる．式(17)におけ る第一項及び第二項はFDTD法の中で簡単に定式化 できるため，第三項に対してのみFILT及びProny法 を適用し，次のようにz 領域における表現を求めた． まずFILTを適用する前に，式(17)を複素周波数で 表現し，緩和時間τaとτbを時間ステップ間隔Δtで 規格化する．これは，FILTを用いて時間領域インパル ス応答を求める際，便宜上，規格化時間t/Δt = n = 1, 2, 3, . . .において時間領域の解を求められるようにす るためである．したがって，次式で表される電気感受 率χ(s)に対してFILTを適用する． χ(s) = ⎡ Δχ ⎣_sτa Δt 1−α +  sτb Δt 1−β γ (18) 計算に用いたJonscher-Raicu分散式の中の各パラメー タを表1に示す． 式(18)に 対 し てFILT及 びProny 法 を 適 用 し ，

Δt = 3.85 psのときの各係数 AE l，pEl，BEk，CkE， rE k 及びqEk を求めた．その結果，NlE= 9，NkE= 0 となり，AE l 及びpEl のみ値をもち，それらの値を表2 にまとめた． 定式化の妥当性を検討するために，文献[30]に示す 一次元モデルにおいて，媒質の反射係数を求めた．解 析空間のセル数は5000で，解像度は2 mmである． 入射電界をT0 = 0.385 ns，a0 = 0.146 nsのガウシ アンパルスを用いる．本論文における計算は全てワー クステーション(Intel Xeon E5-2697v3 @ 2.6 GHz，

256-GB memory，22 parallel threads)上で行われた．

図2にJonscher-Raicu分散モデルを有する媒質の反 射係数を10 MHz∼4 GHzまでの周波数範囲内で示 す．反射係数は10 MHzにおいて0.93であるが，周 波数が高くなるにつれて低下し，4 GHzにおいて0.34 表 1 Jonscher-Raicu分散式における各パラメータ 表 2 Jonscher-Raicu分散式に対する更新式中の係数 図 2 Jonscher-Raicu媒質の反射係数 である．図2より本提案手法によって求めた反射係数 はよく理論値と一致していることから，一次元におけ る定式化の妥当性を確認できた． 一方，生体組織等に用いられるCole-Cole分散モデ ルは式(17)の一般的なJonscher-Raicu分散式の中に α = 1及びγ = 1とした式に等しいため，同様に提 案手法を適用できる．文献[30]にはCole-Cole分散モ デルに対する提案手法の妥当性について議論している ため，そちらも参照されたい．

3.

散 乱 解 析

本章では誘電体球による散乱解析を行い， Jonscher-Raicu分散モデルに対する提案手法の妥当性を示す． 3. 1 解析モデル まず散乱解析に用いた誘電体球を図3に示す．誘電 体球の半径は rd = 10 cmで，媒質は前章で述べた Jonscher-Raicu分散モデルで表されるものとした．誘 電体球モデルは立方体セルで離散化していることから， 球の表面は階段近似としている．FDTD計算に必要な 項数は9個で，更新式の係数の値を表2に示す．図3に 示されるように，平面波はy軸と平行な偏波を有する ものであり，+z方向から誘電体球に入射している．入 射電界は1 V/mとした．また解析モデルの解像度及び 時間ステップ幅をそれぞれΔx = Δy = Δz = 2 mm， Δt = 3.85 psとした．解析領域の大きさは141 x 141 x 141セルである．領域外に向かって伝搬する電磁波を 吸収するために，8層の完全整合層を適用している． 総計算ステップ数は10000ステップである．誘電体球 の場合，10000ステップにかかった計算時間は8分40 秒であった．なお，10000ステップの計算後の解析空 間内における電磁エネルギー量はその最大値と比較し て，−63.2 dBと非常に小さいため，解が収束してい ることを確認した． 3. 2 解 析 結 果 図4には誘電体球の中心点からz軸に沿って+4 cm 図 3 散乱解析に用いる Jonscher-Raicu 分散性誘電体球

及び−4 cm離れた点において記録された電界の時間 波形を示す．図4に示されている時間領域における Mieの解析解は周波数領域におけるMieの解析解の逆 フーリエ変換によって得られたものである[36]．図4 から提案手法によって求めた誘電体球内部の電界波形 はMieによる解析解と非常によく一致しているため， 提案手法の妥当性を確認した．また図4に示すように 電磁界が誘電体中で伝搬していくにつれ，遅延応答に おいて，Mieの解析解との差が表れていることが示さ れている．これは，FDTD法において離散化された空 間内の電磁波の伝搬速度が実際の伝搬速度と異なるた めである[2]．この現象は特に誘電体内の伝搬に対して 顕著である[37]． 図5に10 MHz∼2 GHzの範囲内の各周波数にお ける散乱断面積及びMieの解析解を示す．ここで，全 ての数値結果は提案手法によって一度の時間領域で のFDTD解析によって得られたことに注意されたい． 図5より提案手法の解析結果はMieの解析と非常に 図 4 誘電体球の中心からz 軸上 −4 cm 及び +4cm 離 れた点における電界波形の計算結果及び Mie 解析解 図 5 一度の FDTD 計算によって得られた各周波数にお ける散乱断面積及びその Mie 解析解 よく一致していることから，ここでも散乱問題に対し て，提案手法の有効性を示すことができた．

4.

人体近傍での広帯域アンテナの解析

本章では，提案手法を用いて人体近傍の広帯域アン テナの設計に対する有用性及び計算効率性を示す． 4. 1 広帯域アンテナの構造及び解析モデル まず解析に用いた広帯域アンテナの構造を図6に 示す．アンテナはカプセル内視鏡や体内インプラン ト等のような小型医療デバイスとの通信を目的とし たものである．設計目標は人体と接触している状態 で，300 MHz∼3 GHzの周波数帯で電圧定在波比 (VSWR)が2以下である．このアンテナの形状はフェ ライトを除けば，自己補対アンテナの一種である[38]． アンテナは，大きさa × a = 0.03 × 0.03 m2 の2個 の正方形パッチ及び磁性体フェライトから構成されて おり，給電点におけるギャップ幅を wg = 2 mmと し，内部抵抗を50 Ωとした．フェライトは大きさ b × b = 0.1 × 0.1 m2_，厚さ hf = 6 mmである．フェ ライトは特に低周波帯におけるアンテナ特性の改善及 び高周波における放射電磁界の抑制のために用いる． フェライトと正方形パッチとの間の間隔はhd= 2 mm とした．ここで用いたTDK社製フェライトの複素比 透磁率の実測値を参考にし，次式で示されているロー レンツ分散式の中の各パラメータを推定した．なお， 複素透磁率の測定周波数の上限は1 GHzである． μr(ω) = μr+ jμr = μr∞+ Δχ 1 + 2jωδτL− ω2τL2 , (19) ここで，μr∞，Δχ，δ及びτLはそれぞれ無限周波数 における比透磁率(= 1.0)，ローレンツ分散における比 透磁率の変化分(= 1600)，減衰係数(= 14)及びロー レンツ分散の緩和時間(= 0.83 ns)である．またフェ ライトの比誘電率及び導電率は13及び0.002 S/mと した．図7には，上述のパラメータで求めた複素比透 図 6 数値解析に用いる広帯域アンテナの計算モデル

図 7 フェライトの複素比透磁率の実部と虚部（μr∞ = 1.0，Δχ = 1600，δ = 17，τL= 0.83 ns） 図 8 広帯域アンテナを含めた解析モデル 磁率及び実際に測定した複素比透磁率の実部及び虚部 を示しており，両者はよく一致していることから推定 されたパラメータの妥当性を確認できた．なお，本研 究で設計したアンテナは数値解析における提案手法の 有用性を示すために用いるものであり，アンテナは最 適な形状ではないことに注意されたい． 図8 (a)及び(b)に本研究で検討している解析モ デルを示す．図8 (a)の場合，広帯域アンテナ近傍に は，筋肉組織の電気定数を有する誘電直方体が配置 されている．アンテナ設計のために，生体組織の電 気定数の周波数依存性は考慮されなければならない 重要な要素であり[39]，本研究では電気定数が4項の Cole-Coleモデルによって表されるとする[32]．フェ ライトの周波数分散性を考慮するために，時間ステッ プ幅を Δt = 3.85 psとし，式(19)に対して同様に FILT及びProny法を適用し，FDTD計算における フェライト媒質のための更新式中の係数を求め，表3 に示す．一方，筋肉の誘電率を有する誘電直方体の大 表 3 フェライトに対する更新式中の係数 きさは人体胴体を模擬するために0.2 × 0.2 × 0.1 m3 (d = 0.2 m，dz= 0.1 m)とし，アンテナと誘電直方 体との間の間隔をla = 2 mmまたは0 mm（接触し ている状態）とした．筋肉組織に対して，FILT及び Prony法によって求めたFDTD更新式に必要な項数 は21である．その他の人体を構成する全て51種類 の生体組織もFDTD解析を行う前に，更新係数AlE 及びplE を求めた後，文献[30]に示す一次元モデル において反射係数を算出した．その結果，いずれの生 体組織においても10 MHz∼4 GHzの範囲内で反射 係数が理論値と比較して2%以内で一致しているため， これらの係数を実際のFDTD計算に用いる．なお， Cole-Cole分散モデルに対する提案手法の妥当性につ いて論じられている文献[30]も参照されたい． 図8 (a)の解析モデルの大きさは，141×141×101セ ルで，解像度は2 mmである．入射電圧はT0= 385 ps， a0= 146 ns及び1 Vの振幅を有するガウシアンパル スとして与える[30]．前節と同様に8層のPMLを用 いる．計算ステップ数は4000である．周波数分散性 を有するフェライト及び筋肉誘電直方体が存在する場 合においても，4000ステップにかかった計算時間は前 節と同様な計算環境上で，約77秒である．一方，図 8 (b)に示す解析モデルには，解剖学的な数値人体モデ ル近傍に広帯域アンテナを配置した．アンテナとの距 離はlb= 0 としているが，人体表面は平らではない ため，アンテナが完全に接触している状態ではないこ とに注意されたい．配置場所はそれぞれカプセル内視 鏡及び植込み医療機器の位置を想定し，人体腹部及び 胸部である．人体モデルは情報通信研究機構が開発し た「TARO」の上半身だけのモデルを用いた[40]．解 析モデルの大きさは360× 200 × 549セルである．解 像度及び給電に用いる入射電圧波形は図8 (a)のモデ ルと同じである．10000ステップにかかった計算時間 はそれぞれCase A及びCase Bの場合，3912秒及び 4382秒であった． 4. 2 解 析 結 果 図8 (a)の解析モデルを用い，広帯域アンテナの反 射係数を求めた結果を図9に示す．図9に示すよう

図 9 筋肉の誘電直方体近傍の広帯域アンテナの反射係数 図 10 各解析周波数における誘電直方体内のx 成分の電 界分布 に，誘電直方体がない場合の共振周波数は1 GHz及び 2 GHz付近である．帯域幅は1 GHzにおいて28%で， 比較的狭帯域である．一方，筋肉の誘電直方体から la= 2 mm離れてアンテナを配置した場合の帯域幅は 121%である．誘電体と接触している場合（la = 0）， アンテナの帯域幅は更に164%に増大し，VSWR = 2 以下の周波数範囲は297 MHz∼3.04 GHzであり，設 計目標(300 MHz∼3 GHz)を満足している．また フェライトがない場合，図9に示すように，帯域幅が 360 MHz∼3.25 GHzの160%と小さくなり，設計目 標を達成できなかった．したがって，低周波帯におけ るアンテナの特性を改善するためには，フェライトの 設置が必要であることがわかる．以上の結果より，広 範囲にわたる周波数分散特性を有する媒質を含むアン テナの特性解析にも本提案手法が有効であることが示 された． 図10にアンテナが誘電直方体と接触している場合 の300 MHz，500 MHz，1 GHz及び2 GHzにおけ 図 11 誘電直方体の中心点におけるx 成分の電界波形 る誘電体内の x成分の電界分布を示している．横軸 にはアンテナからの距離を示す．図10には，提案手 法で求めた結果及び筋肉の電気定数を用いてそれぞれ の周波数においてYeeアルゴリズムを用いて解析した 結果を示している．図10から分かるように，一度で 提案手法により求めた電界分布は，それぞれの解析周 波数で求めた場合とよく一致しており，本提案手法に よる広帯域の解を求めることが可能であることを示し た．次に，複素誘電率の周波数依存性を正確に考慮す る必要性を示すために，筋肉の誘電直方体の複素誘電 率を一定にした場合(図中の凡例：Yee @ 300 MHz) 及び提案手法で計算した場合の誘電体内部の中心点に おけるx成分の電界の時間波形を図11に示す．図11 よりCole-Coleモデルを用いた場合，電界のパルス波 形は滑らかに変化するのに比べ，300 MHzにおける 一定な電気定数を用いた場合，電界波形はCole-Cole モデルのものとかなり異なり，パルスの振れ幅も大き くなる．電界のピーク値についてもCole-Coleモデル と比べ，約42%大きい．ここでは，電界のx成分の みを示しているが，他の成分についても同様な傾向が 見られる．したがって，生体組織内において広帯域パ ルス電磁界を利用した無線通信の性能評価等を行うた めには，誘電率の周波数依存性を正確に考慮する必要 があることが示唆されている． 次に，図8 (b)の解析モデルを用い，人体近傍におい て広帯域アンテナの特性を評価した．図12に設計した 広帯域アンテナを人体の胸部(Case A)及び腹部(Case B)に配置した場合の反射係数を示す．合わせて前述 で，筋肉の誘電直方体近傍に置いた場合(la= 0 mm) の反射係数も示してある．結果より設計した広帯域ア ンテナの反射係数は人体近傍の配置場所によってかな

図 12 人体の胸部 (Case A) 及び腹部 (Case B) 近傍に おける広帯域アンテナの反射係数 図 13 設計した広帯域アンテナを人体近傍に配置した場 合の心臓 (Case A) 及び小腸 (Case B) 内のx 成 分の電界の時間波形 り変化し，設計された目標から外れることが分かった． 特に胸部に配置した場合，VSWRが2以下の下限周 波数が594 MHzにシフトする結果となった．したがっ て，この広帯域アンテナを安定的に動作させるために は，人体とアンテナとの間の整合層が必要であると考 えられる． 図13には，提案手法によって求めた心臓(Case A) 及び小腸(Case B)内のx成分の電界波形を示す．ア ンテナからの距離はCase Aの場合，5.8 cmで心臓壁 内に観測点を設けたのに対して，Case Bの場合の距 離は6.0 cmである．図13より心臓内のパルス波形の 幅は小腸内のそれよりも小さいことが見られる．これ は図12に示すようにアンテナを人体胸部に置いた場 合，低周波帯における反射係数が大きくなり，アンテ ナ特性が悪化していることから，アンテナから低周波 成分の電磁界が抑制されたためである．またパルスの 継続時間は約4 nsで，両者ともおおむね同じである 表 4 各解析モデルに対する計算時間 ことが分かった．このように，提案手法によって解剖 学的人体モデル内における広帯域パルス電磁界の波形 を正確に求めることが可能であることを示した．今後， 広帯域パルス電磁界を用いた人体内外通信やカプセル 内視鏡や体内小型医療デバイスの位置検出等の様々な 応用に本提案手法を用いることが考えられる． 4. 3 提案手法の計算効率 最後に，提案手法の計算効率について考える．表4 には本研究で用いた解析モデルに対する計算時間をま とめた．計算ステップ数は全ての解析モデルにおいて， 10000ステップとした．まず散乱問題については，提 案手法によるJonscher-Raicu型の分散性を有する誘 電体球の解析にかかる計算時間は375秒であった．同 様なモデルで時間領域の波形を離散的にフーリエ変換 して14点の周波数における散乱断面積の結果を求め る場合，計算時間は約433秒であった．一方，一定な 電気定数を代入してYeeアルゴリズムにより一つの周 波数での解を求めるのに必要な計算時間は324秒であ るので，単純に14点の周波数における解を求めるため には，324× 14 = 4536秒必要であるが，周波数分散 性を考慮した提案手法では，433秒で計算を終えるこ とができる．更に，解析周波数の点数が多くなればな るほど，本提案手法の計算効率がより高くなる．表4 より放射問題についても同様なことがいえる．誘電直 方体近傍の広帯域アンテナの解析では，単純に時間領 域での解（時間波形）を求めるだけであれば，10000 ステップの計算でかかった時間は279秒である．提案 手法による解析は一定な電気定数を用いた解析の計算 時間(183秒)よりも長くかかるが，アンテナの反射特 性を広帯域で求めることができる．また20点の周波 数における電界分布を求めるために追加的に必要な計 算時間はわずか4秒であり，本提案手法の計算効率の

高さが示せた．最も有利な点としてこれまで生体内に おける広帯域パルス電磁界の正確な評価ができなかっ たが，本研究で示したように，人体近傍の広帯域アン テナからの広帯域パルス電磁界の伝搬や体内における 電磁界の時間波形を提案手法によって求めることがで きる．今後，提案手法の広帯域無線通信技術への様々 な応用が期待できる．

5.

む す び

本研究では，著者らがこれまで提案してきた高速 逆ラプラス変換及びProny法を用いた周波数依存性 FDTD法について，一般定式化を行い，電磁界の散 乱問題及び放射問題に対して提案手法の有効性及び有 用性を示した．まずJonscher-Raicu分散式への本手 法の適用可能性について検討し，誘電体球の散乱解析 においてMieの解析解とよく一致しているため，手 法の妥当性を確認した．次に，放射問題については， 広帯域アンテナの設計に提案手法を用い，人体近傍に おけるアンテナの特性を様々な条件で評価し，手法の 有用性を示した．その結果，広帯域アンテナは近傍に 置いた誘電体または人体に強く影響を受けることが分 かった．また提案手法を用いて不均一人体モデル内に おける広帯域パルス電磁界の時間波形を求めることが できることを示した．最後に，本提案手法の計算効率 について論じ，複数点の周波数での解を求める場合に おいても，本提案手法が非常に有効であることが確認 された． 今後の研究課題として，提案手法を用いて人体内外 の無線通信設計や体内小型医療デバイスの位置推定等 の検討を行っていくことが考えられる． 謝辞 本研究は平成30年度日本学術振興会から科 学研究費補助金若手研究（課題番号18K18376）の助 成を受けて行われた． 文 献

[1] K.S. Yee, “Numerical solution of initial bound-ary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag., vol.14, no.3, pp.302–307, 1966.

[2] A. Taflove and S.C. Hagness, Computational electro-dynamics: the finite-difference time-domain method, Artech House, 2005.

[3] R.F. Harrington, Field Computation by Moment Methods, Wiley-IEEE Press, 1993.

[4] J.M. Jin, The Finite Element Method in Electromag-netics, John Wiley & Sons, 1993.

[5] V.A. Thomas, M.E. Jones, M. Piket-May, A. Taflove,

and E. Harrigan, “The use of SPICE lumped circuits as sub-grid models for FDTD analysis,” IEEE Mi-crow. Guided Wave Lett., vol.4, no.5, pp.141–143, 1994.

[6] P.M. Goorjian and A. Taflove, “Direct time integra-tion of Maxwell’s equaintegra-tions in nonlinear dispersive media for propagation and scattering of femtosecond electromagnetic solitons,” Opt. Lett., vol.17, no.3, pp.180–182, 1992.

[7] 岩本 敏，荒川泰彦，“FDTD 法のフォトニック結晶研究

への応用，”レーザー研究，vol.36, no.10, pp.614–620, 2008.

[8] J. Wang, O. Fujiwara, S. Kodera, and S. Watanabe, “FDTD calculation of whole-body average SAR in adult and child models for frequencies from 30 MHz to 3 GHz,” Phys. Med. Biol., vol.51, no.17, pp.4119– 4127, 2006.

[9] P.J. Dimbylow, A. Hirata, and T. Nagaoka, “In-tercomparison of whole-body averaged SAR in Eu-ropean and Japanese voxel phantoms,” Phys. Med. Biol., vol.53, no.20, pp.5883–5897, 2008.

[10] J. Chakarothai, K. Wake, and S. Watanabe, “Conver-gence of a single-frequency FDTD solution in numer-ical dosimetry,” IEEE Trans. Microw. Theory Tech., vol.64, no.3, pp.707–714, 2016.

[11] R. Luebbers, F.P. Hunsberger, K.S. Kunz, R.B. Standler, and M. Schneider, “A frequency-dependent finite-difference time-domain formulation for disper-sive materials,” IEEE Trans. Electromagn. Compat., vol.32, no.3, pp.222–227, 1990.

[12] R.J. Luebbers and F. Hunsberger, “FDTD for Nth-order dispersive media,” IEEE Trans. Antennas Propag., vol.40, no.11, pp.1297–1301, 1992. [13] T. Kashiwa, Y. Ohtomo, and I. Fukai, “A

finite-difference time-domain formulation for transient propagation in dispersive media associated with Cole-Cole’s circular ARC law,” Microw. Opt. Technol. Lett., vol.3, pp.416–419, 1990.

[14] R.M. Joseph, S.C. Hagness, and A. Taflove, “Di-rect time integration of Maxwell’s equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses,” Opt. Lett., vol.16, no.18, pp.1412–1414, 1991. [15] O.P. Gandhi, B.-Q. Gao, and J.-Y. Chen, “A

frequency-dependent finite-difference time-domain formulation for induced current calculations in hu-man beings,” Bioelectromagnetics, vol.13, no.6, pp.543–555, 1992.

[16] D.M. Sullivan, “Frequency-dependent FDTD meth-ods using Z transforms,” IEEE Trans. Antennas Propag., vol.40, no.10, pp.1223–1230, 1992. [17] D.M. Sullivan, “Z-transform theory and the FDTD

method,” IEEE Trans. Antennas Propag., vol.44, no.10, pp.28–34, 1996.

invari-ance method in Z-transform theory for frequency-dependent FDTD methods,” IEEE Trans. Antennas Propag., vol.57, no.9, pp.2683–2690, 2009.

[19] M. Mrozowski and M.A. Stuchly, “Parameterization of media dispersive properties for FDTD,” IEEE Trans. Antennas Propag., vol.45, no.9, pp.1438–1439, 1997.

[20] D.F. Kelley, T.J. Destan, and R.J. Luebbers, “Debye function expansions of complex permittivity using a hybrid particle swarm-least squares optimization approach,” IEEE Trans. Antennas Propag., vol.55, no.7, pp.1999–2005, 2007.

[21] S. Su, W. Dai, D.T. Haynie, and N. Simicevic, “Use of the z-transform to investigate nanopulse penetration of biological matter,” Bioelectromagnetics, vol.26, no.5, pp.389–97, 2005.

[22] I.T. Rekanos and T.G. Papadopoulos, “FDTD mod-eling of wave propagation in Cole-Cole media with multiple relaxation times,” IEEE Antennas Wirel. Propag. Lett., vol.9, pp.67–69, 2010.

[23] L. Mescia, P. Bia, and D. Caratelli, “Fractional derivative based FDTD modeling of transient wave propagation in Havriliak-Negami media,” IEEE Trans. Microw. Theory Tech., vol.62, no.9, pp.1920– 1929, 2014.

[24] D. Caratelli, L. Mescia, P. Bia, and O.V. Stukach, “Fractional-calculus-based FDTD algorithm for ul-trawideband electromagnetic characterization of ar-bitrary dispersive dielectric materials,” IEEE Trans. Antennas Propag., vol.64, pp.3533–3544, 2016. [25] P. Bia, D. Caratelli, L. Mescia, R. Cicchetti, G.

Maione, and F. Prudenzano, “A novel FDTD formu-lation based on fractional derivatives for dispersive Havriliak-Negami media,” Signal Process., vol.107, pp.312–318, 2015.

[26] T. Hosono, “Numerical inversion of Laplace trans-form and some applications to wave optics,” Radio Sci., vol.16, no.6, pp.1015–1019, 1981.

[27] S. Kishimoto, S. Ohnuki, Y. Ashizawa, K. Nakagawa, and W.C. Chew, “Time domain analysis of nanoscale electromagnetic problems by a boundary integral equation method with fast inverse Laplace transform,” J. Electromagnet. Waves, vol.26, pp.997– 1006, 2012.

[28] S. Ohnuki, Y. Kitaoka, and T. Takeuchi, “Time-domain solver for 3D electromagnetic problems using the method of moments and the fast inverse Laplace transform,” IEICE Trans. Electron., vol.E99-C, no.7, pp.797–800, July 2016.

[29] J. Chakarothai, “Novel FDTD scheme for anal-ysis of frequency-dependent medium using fast inverse Laplace transform and Prony’s method,” IEEE Trans. Antennas Propag., vol.67, Jan. 2019. (Early Access) DOI: 10.1109/TAP.2018.2878077 [30] J. Chakarothai, S. Watanabe, and K. Wake,

“Numer-ical dosimetry of electromagnetic pulse exposures us-ing FDTD method,” IEEE Trans. Antennas Propag., vol.66, no.10, pp.5397–5408, 2018.

[31] V. Raicu, “Dielectric dispersion of biological mat-ter: Model combining Debye-type and “universal” re-sponses,” Phys. Rev. E, vol.60, no.4, pp.4677–4680, 1999.

[32] S. Gabriel, R.W. Lau, and C. Gabriel, “The dielec-tric properties of biological tissues: III. Paramedielec-tric models for the dielectric spectrum of tissues,” Phys. Med. Biol., vol.41, no.11, pp.2271–2293, 1996. [33] S. Kensuke, W. Kanako, and W. Soichi,

“Develop-ment of best fit Cole-Cole parameters for measure-ment data from biological tissues and organs between 1 MHz and 20 GHz,” Radio Sci., vol.49, no.7, pp.459– 472, 2014.

[34] W.H.P. Pernice, F.P. Payne, and D.F.G. Gallagher, “A general framework for the finite-difference time-domain simulation of real metals,” IEEE Trans. An-tennas Propag., vol.55, no.3, pp.916–923, 2007. [35] M. Tanaka, J. Chakarothai, and H. Inoue, “An

FD-FDTD modeling of ferrite,” IEICE Trans. Electron. (Japanese Edition), vol.J86-C, no.9, pp.1030–1031, Sept. 2003.

[36] J.A. Stratton, Electromagnetic Theory, New Jersey, John Wiley & Sons, 2007.

[37] J. Chakarothai, Q. Chen, and K. Sawaya, “Three-dimensional electromagnetic scattering analysis of lossy/lossless dielectrics using constrained interpo-lation profile method,” AP-RASC’10, no.KB2-5, Toyama, Japan, 2010.

[38] Y. Mushiake, “Self-complementary antennas,” IEEE Antennas Propag. Mag., vol.34, no.6, pp.23–29, 1992. [39] P. Kosmas, C.M. Rappaport, and E. Bishop, “Model-ing with the FDTD method for microwave breast can-cer detection,” IEEE Trans. Microw. Theory Tech., vol.52, no.8, pp.1890–1897, 2004.

[40] T. Nagaoka, S. Watanabe, K. Sakurai, E. Kunieda, S. Watanabe, M. Taki, and Y. Yamanaka, “Devel-opment of realistic high-resolution whole-body voxel models of Japanese adult males and females of av-erage height and weight, and application of models to radio-frequency electromagnetic-field dosimetry,” Phys. Med. Biol., vol.49, no.1, pp.1–15, 2004.

（2018 年 9 月 4 日受付，12 月 23 日再受付， 2019年 4 月 15 日公開）

チャカ ロ タ イ ジェド ヴィス ノ プ （正員） 平 15 秋田大卒．平 22 東北大大学院博 士課程了．現在，情報通信研究機構電磁波 研究所電磁環境研究室に所属．広帯域電 磁界解析手法や無線電力伝送システムに 対するばく露評価などの生体電磁環境研 究に従事．平 26 URSI 若手科学賞，平 27 電気学会基礎・材 料・共通部門優秀論文発表賞，平 30 本会エレクトロニクスシ ミュレーション研究会優秀論文発表賞受賞．電気学会，信学会， Bioelectromagnetics Society，IEEE，ACES 各会員． 和氣加奈子 （正員） 平 12 都立大大学院工学研究科電気工学 専攻博士課程了．同年郵政省通信総合研究 所（現，国立研究開発法人情報通信研究機 構）入所．現在，電磁波研究所電磁環境研 究室主任研究員．生体電磁環境研究に従事． 平 11 URSI 若手科学賞，平 21 通ソ論文賞， 平 27 第 59 回前島密賞，平 29 第 27 回電気学術振興賞（進歩 賞）等受賞．電気学会，信学会，IEEE，Bioelectromagnetics Society等の会員． 渡辺 聡一 （正員） 平 8 都立大（現，首都大学東京）大学院 工学研究科電気工学専攻博士課程了．同年 郵政省通信総合研究所（現，国立研究開発 法人情報通信研究機構）入所．工博．現在， 電磁波研究所電磁環境研究室研究マネー ジャ．平 8 URSI 若手科学賞，平 9 信学会 論文賞，平 10 信学会学術奨励賞，平 11 信学会功労感謝状，平 17 Physics in Medcine and Biology誌 The Roberts Prize， 平 21 信学会通信ソサイエティBest Letter Award，平 24 文 部科学大臣表彰科学技術賞，平 25 日本 ITU 協会賞，平 26 前島密賞各受賞，平 17 より国際非電離放射線防護委員会主 委員会委員，平 25 より信学会シニア会員．電気学会，IEEE， Bioelectromagnetics Society各会員． 陳 強 （正員：フェロー） 昭 63 年西安電子科技大学卒．平 6 東北 大大学院博士課程了．同大学助手，助教授， 准教授を経て，平 25 より同大学院工学研 究科通信工学専攻電磁波工学分野教授．ア ンテナ，マイクロ波・ミリ波，電磁界の測 定法及び数値解析法の研究等に従事．平 5 本学会学術奨励賞受賞．平 8，18，22 と 24 通ソ活動功労賞． 平 26 通ソ功労顕彰状．平 20 本学会論文賞，第 2 回喜安善市 賞．信学会環境電磁工学研究専門委員会幹事，アンテナ・伝播 研究専門委員会幹事，光応用電磁界計測時限研究専門委員会初 代委員長，無線電力伝送研究専門委員長を歴任．平 29 本学会 フェロー． 澤谷 邦男 （正員：フェロー） 昭 46 東北大学・工・通信卒．昭 51 同 大学大学院博士課程了．同大学助手，助教 授，教授を経て現在同大学産学連携機構特 任教授．プラズマ中のアンテナ，電磁波の 散乱・回折理論，移動通信用アンテナの研 究に従事．工博．IEEE フェロー，映像情 報メディア学会会員．昭 56 本会学術奨励賞，昭 63 同論文賞， 平 18 同通ソ論文賞，平 21 同論文賞・喜安善市賞，平 26 同通 ソ論文賞，平 29 同功績賞受賞．