株式会社 科学技術研究所
ホーム > FDTD法とは > FDTD法とは 第1章

FDTD法とは

 

総論

FDTD法とは電界E[V/m]、磁界H[A/m]、電束密度D[C/m2]、磁束密度B[T]、電流密度I[A/m2]とした場合に 以下の(1),(2)式で表されるMaxwell方程式を時間的(Time)、空間的(Domain)に離散化し陽的な時間進行法を用いて経時的な電磁界強度を求める手法です。ここでは、FDTD法について具体的に分かりやすくご説明します。尚、電流密度は J で示すケースも多いですが、本ページでは虚数単位 j と混同しないように I で示します。*FDTD法とは、Finite-Difference Time-Domain methodの略号です。
Keyword:FDTD,FDTD法,電界,磁界,電束密度,磁束密度,Maxwell方程式,FDTD method,Finite-Difference,Time-Domain method,electromagnetic field,Maxwell equation.

 

Overview

FDTD method for an electromagnetic field is a numerical scheme to obtain the time domain distribution of electromagnetic field by solving the explicitly discretized Maxwell's equations which are represented by the following equations (1) and (2). Here, we define an electric field as E[V/m], a magnetic field as H[A/m], an electric flux density D as [C/m2], a magnetic flux density B as [T], and a current density I as [A/m2].
( A current density is often indicated by J, but in this site, it is indicated by I to avoid confusion with the imaginary unit j.)
(FDTD method is an abbreviation for Finite-Difference Time-Domain method. )

格子は多くの場合、図1のような電界及び磁界ベクトルの配置が半格子ずれたスタガード格子(Staggered Grid)を用います。 FDTD法では電磁波の波形を再現するために格子は解析対象とする波長よりも十分に小さい必要があります。 かつ時間的な陽的解法ですので、時間進行法に用いるタイムステップは格子幅Δh及び光速cを用いてCFL条件(Courant-Friedrichs-Lewy Condition; クーラン条件)から求められるΔh/cが上限値となります。 従って、電磁波の伝搬の様子を空間・時間的に可視化できるという大きなメリットがあり、また解析領域の大きさと解析対象波長が近い場合には非常に効率的に解が得られます。 一方で解析領域が波長の100倍程度になると格子点数が 膨大なものとなり解析が困難になるというデメリットがあります。

スタガード格子

図1 スタガード格子

可視光(波長400nm~800nm程度)を扱う場合にはサブマイクロ~マイクロメートルの解析領域を対象とするような問題で非常に効率的な解析が可能です。 特にナノスケールの構造をもった物体の電磁波の反射・散乱・透過特性の解析や プラズモン解析などでは多くの実績があります。 また電子レンジのようなマイクロ波加熱問題も波長と解析領域の大きさが近い場合が多く効率よく解析ができます。 さらに波長の長い電波は建物の壁面を考慮することはできませんが、建物を一つの物体と見立てれば大空間における解析も可能です。

可視化の様子(赤:電界、青:磁界)

図2 電磁界の可視化様子(赤:電界、青:磁界)
上図は、フリーソフト「KeyFDTDちゃん」の計算サンプルです。

お問い合わせ・資料ダウンロードのお申し込みはこちら
 

お問い合わせ・資料ダウンロード

 

FDTD・FDTD法の詳解目次
第1章 FDTD電磁波解析の基礎用語
第1節 支配方程式・微分・差分
第2節 メッシュ・離散化・計算時間
第3節 励振源・境界条件・収束判定・解析の実効値
第4節 解析事例(マルチモードキャビティによる加熱)
 
第2章 FDTD電磁波解析を用いた研究計画
第1節 シミュレーションプラン
第2節 結果までの難所
第3節 プラン作成実践
第4節 解析事例(マルチモードキャビティによる培地加熱)
 
第3章 FDTD電磁波解析と仮定
第1節 マクスウェル方程式が仮定していること
第2節 分極の振る舞いに関する仮定
第3節 FDTD法に関する仮定
第4節 解析事例(2つの位相制御ポートを持つ円筒キャビティ)
まとめ FDTD法で仮定されていることと注意点・対処
 
第4章 FDTD電磁波解析結果の報告・発表・論文化
第1節 妥当性の主張文例
第2節 シミュレーション結果の提示
第3節 報告書・プレゼンテーションについて
第4節 解析事例(NaCl水溶液のマイクロ波加熱)
まとめ シミュレーションの妥当性に関する注意点・対処
 
お問い合わせ・資料ダウンロードのお申し込みはこちら
 

お問い合わせ・資料ダウンロード

 

第1章 FDTD電磁波解析の基礎用語

 

第1章 第1節 支配方程式・微分・差分

 

Maxwell方程式

 
Maxwell方程式を微分型で書いたのが上の(eq.1.a)~(eq.1.d)です。時間tと位置xにおける物理量を表すため、従属変数はすべて時間と位置ベクトルの関数として表されます。式中の文字変数はそれぞれ電界E[V/m]、磁界H[A/m]、電束密度D[C/m2]、磁束密度B[T]、電流密度I[A/m2]、電荷密度ρ[C/m3]です。
(eq.1.a)は磁力線が閉曲線であることを示しており、また右辺が0なので磁気単極子が存在しないことを表しています。(eq.1.b)は磁場の時間変化が電界の回転を発生(電磁誘導)させることを示しています。(eq.1.c)はある体積を見たとき、その表面から出る電束密度ベクトルが、内部に含まれる電荷密度に等しいことを表しています。(eq.1.a)と対比させると右辺の値が存在しますが、これは電荷単極子が存在することを意味しています。(eq.1.d)は電束密度の時間変化と電流密度が磁場の回転を発生することを意味しています。
この方程式系には四則演算と時間と空間に関する1回微分計算のみが含まれているため後述の離散化が容易です。またこれらの方程式から電磁波の波動方程式も導出できます。

 

コントロールボリュームと微分演算子

 
コントロールボリュームは物理量を観測する空間のことです。日本語では検査体積ですが、コントロールボリューム(CV)の方が一般的です。
Maxwell方程式はCVの電磁界特性を微分演算子を用いて記述したものです。微分演算子は数値解析やシミュレーションを行う場合、理解が必要です。厳密な意味での数学的定義の理解は後回しにして、その演算子が物理的に意味することを理解することが重要です。上記説明で∇・はCVの収支と紹介していますが発散の意味もあります。

 

微分と差分

 
微分の定義を数学で最初に学ぶのは高校生の頃です。その導入は時間-移動距離のグラフで平均変化率の極限が微分だという説明だった記憶があります。
コンピュータで数値計算するにあたり微分演算子はどのように扱うでしょうか。数値計算では、数学の学習とは反対に微分を差分(平均変化率)で近似します。この差分近似により、コンピュータが得意な四則演算のみで微分を計算できます。
FDTD法は微分型のMaxwell方程式を離散化(差分近似)する点に特徴があります。ベクトルの回転に空間微分が隠れている点も踏まえ、Maxwell方程式には空間微分と時間微分が含まれています。従ってFDTD法の電磁波解析ソフトウェアには空間の離散化と時間の離散化の手続きが必要です。
お問い合わせ・資料ダウンロードのお申し込みはこちら
 

お問い合わせ・資料ダウンロード

 

第1章 第2節 メッシュ・離散化・計算時間

メッシュ生成ー物性値と形状再現

 
空間の離散化に先だってメッシュ生成の手続きが必要です。メッシュ生成は、1.メッシュ分割、2.物性値割当て、という2ステップの手続きです。始めに、解析領域全体を或る有限サイズのメッシュに分割します<メッシュ分割>。次に、各メッシュの交点に物性値を割当てます<物性値割当て>。
励振源及び解析対象もメッシュ分割します。上図は等間隔且つ直交なメッシュで分割し、誘電率を割当てた例です。
1次元の場合メッシュの形状は線分です。解析領域を複数の線分に分割します(i=1~i=2の線分、i=2~i=3の線分、...、i=5~i=6の線分)。そして線分同士の交点に物性値を割当てます(i=1に真空の誘電率、i=2に真空の誘電率、...、i=6に誘電体の誘電率)。コントロールボリューム(CV)の観点で言えば、i=1を中心とした真空のCVがあり、i=2を中心とした真空のCVがあり、...と理解してください。i=3とi=4の中間には真空と誘電体の界面があります。
2次元の場合メッシュの形状は平面です。解析領域を複数の平面に分割し、平面同士の交点に物性値を割当てます(〇の位置に真空の誘電率、×の位置に誘電体の誘電率)。円筒誘電体もメッシュ分割され、点線で示すCVの群で表現されます。
凸形状で近似されることが重要です。解析によってはこの解像度では不足する場合もあり得ますし、十分かもしれません。この見極めは解析の目的を念頭に置いて検討する必要があります。

空間の離散化

 
空間の離散化では、電磁波を記述する上で欠かせない電界及び磁界をメッシュに沿って配置していきます。 空間微分を含む式を良く見ること、登場する物理量の関係に注意することが重要です。
Z方向の1次元空間を伝搬する電磁波を扱いましょう。Z方向に伝搬する電磁波には「電界成分Ey、磁界成分Hxをもつ電磁波」と「電界成分Ex、磁界成分Hyをもつ電磁波」の2通りが考えられますが、ここでは「電界成分Ey、磁界成分Hxをもつ電磁波」に着目します。
Maxwell方程式の(eq.1-b, eq.1-d)を1次元に適用し、電界Eyと磁界Hxに関する式を求めると次の2式が得られます。(eq.7)から磁界Hxの単位時間の変化はEyの空間微分で求められます。また(eq.8)から電界Eyの単位時間の変化はHxの空間微分で求められます。

差分近似の発想に倣って(eq.7)、(eq.8)を式変形した結果が(eq.9)、(eq.10)です。(eq.9)は、2つ並んだ電界から磁界が計算できることを表します。(eq.10)は、2つ並んだ磁界及び自身の過去の値から電界が計算できることを表します。
電磁波は伝搬するもの、という着想から得られるのが上記の図です。Ey1、Hx1+1/2、Ey2、Hx2+1/2、...と互い違いに配置します。Ey1で励振した電磁波が(eq.9)によってHx1+1/2に、更に(eq.10)によってEy2に、そのまた(eq.9)によってHx2+1/2に、と伝搬するわけです。
FDTD法では電界と磁界を空間・時間的にずらして配置(スタガード配置)して計算を効率的に行います。

時間の離散化

 
(eq.9), (eq.10)の左辺は電界EyとHxの時間変化を表していますが、これも離散化が必要です。また前述の通り、電界Eyと磁界Hxは時間的にもずらして配置するので電界が時刻n=1, 2, …に存在するとき磁界は時刻n=1+1/2, 2+1/2, …に配置します。
時刻を離散化して時刻n=1と時刻n=2の間隔がΔtとすると左辺を差分化して(eq.11), (eq.12)に書き換えられます。
(eq.12)でEyin+1/2の項が表れていますがこの値をEyin+1/2=0.5×(Eyin+1+Eyin)として、(eq.11)についてHxi+1/2n+1/2、(eq.12)についてEyin+1を求める式に書き換えれば時間についても離散化された式が得られます。
時間を離散化した結果、時間間隔Δtが現れたことに注意して下さい。Δtは一般的にタイムステップと呼ばれシミュレーションでは重要です。

計算時間の考え方

 
FDTD法ではシミュレーションにかかる時間に関して概ね上の3式で表されることを考慮します。(eq.13)はCFL条件と呼ばれ、タイムステップΔtの間に電磁波がメッシュの間を飛び越えないという条件です。なおこの式の左辺を右辺で割った値はCFL数と呼ばれ、(eq.13)は「CFL数が1以下でなければならない」とも言い換えられます。メッシュ間隔dx, dy, dzが小さくなるとΔtが小さくなることを忘れてはいけません。
FDTD法は物理量を求める際に、過去の情報のみを用いる陽解法なのでCFL数は1以下でなければなりません。物理量を求める際に過去と現在の連立方程式を解く陰解法を用いるFDTD法の場合はCFL数が1以上の場合もありますが、電磁波解析ではこの手法は現段階で一般的ではありません。
(eq.14)は1タイムステップ進めるために必要な解析時間tstepと解析領域の各方向メッシュ数ni, nj, nkの関係を示し、tstepは全メッシュ数(ni×nj×nk)に比例することが分かります。
(eq.15)は解析開始から幅W、奥行D、高さHの解析領域内部に電磁波が行き渡って、定常状態になるまでの解析ステップ数nstepとW、D、H、Δtの関係を示しています。
解析領域を一定として、各方向のメッシュ幅をdx=dy=dz=dhとすると、各方向のメッシュ数はni=W/dh, nj=D/dh, nk=H/dhとなります。またタイムステップΔtはdh/cです。定常解を求めるために必要な全解析時間ttotalは1ステップにかかる計算時間(eq.14)とステップ数(eq.15)の積で表され(eq.16)の関係が成り立ちます。
つまり解析領域が同じ大きさでメッシュ幅を各方向均一のメッシュを使用した場合、定常状態を得るための解析時間は「メッシュ幅の4乗に反比例」します。端的に言うとメッシュの密度を2倍にすると定常状態を得るまでの計算時間は16倍になります。
この問題を避けるためには以下の3点について検討が必要です。

計算時間短縮のための3つの検討

  • メッシュは3方向で同じ幅である必要はなくX, Y, Z方向でメッシュ数を調整します。但しメッシュ幅が極端に異なると誤差の原因となるためそのアスペクト比は5倍以内に収まるように注意します。
  • 少しでもメッシュ幅を大きくします。均一メッシュの場合はメッシュ幅を1.25倍にすると解析時間は40%程度まで短縮できます。
  • 不均一メッシュを使用するケースもありますが、誤差の挙動が複雑なので一概には勧めていません。商用ソフトでは解析領域内にメッシュの細かい領域を作るマルチグリッド機能を備えているものがあり、これを使用するほうが有利です。
 
お問い合わせ・資料ダウンロードのお申し込みはこちら
 

お問い合わせ・資料ダウンロード

 

第1章 第3節 励振源,境界条件,収束判定,解析の実効値

励振源-周波数とモード

 
励振源
4種類の励振源
FDTD法のシミュレーションでは周波数の制限はありません。このため電波~マイクロ波帯~テラヘルツ帯~可視光帯~紫外線の幅広い範囲で活用できます。
励振源の設定で重要なのは励振モードの設定です。例えばマイクロ波帯では多くの場合、上記の4種類の励振源が使用されます。ソフトウェアによってはより多くのモードが準備されていることもありますが、この4種が基本です。
アンテナを励振源とする場合は、モノポールアンテナやダイポールアンテナの場合は点源と誘電体を組み合わせてアンテナそのものをモデル化します。構造の細かなアンテナはそのままモデル化することが困難な場合が多く、如何に等価なアンテナをモデル化するかがポイントになります。

境界条件

 
境界条件
3種類の境界条件
電磁波シミュレーションで用いられる境界条件は主に上に示した3種類です。MUR1及びPML(Perfect Matched Layer)は電磁波の反射を0にする条件として定式化されているため吸収境界条件と呼ばれますが実際には解放をモデル化する場合に使用されます。
完全反射条件PEC(Perfect Electric Conductor)は入射した電磁波を完全反射する条件です。導波管やキャビティ壁面はこの境界条件でしばしばモデル化されます。理想的な(導電率∞の)導体を仮定しているために、壁面での損失をシミュレートすることはできないので注意が必要です。
また周期的な構造を持つ電波吸収体のようなものをシミュレートする場合には周期境界条件を使用します。周期境界条件を用いると1周期の解析領域を計算するだけで、無限に同じ構造が続いていることを模擬できます。構造が距離を置いて存在する場合には、構造と解析領域境界の距離は構造同士の距離の半分になるように配置することに注意が必要です。

収束判定

 
収束判定
収束判定の仕方
一般的にシミュレーションは電磁界が存在しない解析領域に時刻0において電磁波の励振を開始します。またFDTD法は時々刻々の物理量を得る手法(時間進行法)であるため解析領域内に電磁波が十分に行き渡った状態になるまでは過渡的な状態の結果が得られます。
過渡的な結果が得られるのはFDTD法の大きなメリットである一方、次項の実効値を求める場合には過渡的な状態では物理的な意味を持たない結果が得られます。このため実効値を求めるために積算を行う場合は、過渡的な状態が定常的な状態にまで進んでいる(=収束している)かの確認が必要です。
この確認には、解析領域の代表的位置における時系列の電磁界強度を使用します。上図はマイクロ波加熱シミュレーションの例です。キャビティ中心を代表的位置に採用し、時間平均した電界絶対値|E|をチャート表示しました。時間平均とは、時刻tから1周期前の時刻までの物理量を平均することです。
チャートからは次のことが読み取れます。1.解析初期において電磁波強度0から始まり強度が上昇する。2.解析終期において電磁波強度は或る値に漸近する。
電磁波強度の漸近は「電磁波が解析領域に十分行き渡り、それ以上タイムステップを進めても時間平均した電磁波強度分布が変化しない」状態を意味します。これを定常状態と呼称します。
FDTD法においては、定常状態とみなせるまでタイムステップを進めた後に物理量の実効値を求めます。

実効値の導出

 
実効値の導出
時々刻々の値を積分(積算)して実効値を求める
FDTD法は時間進行法であることを前述しましたが、これは得られる物理量はすべて時々刻々の値であることを意味しています。電磁波を連続照射する実験や実機では電磁波の電界成分の時々刻々の値ではなく実効値を用います。
FDTD法では直接実効値を求めることは出来ないので、時々刻々の値を積分(積算)して実効値を求めるのが一般的です。有限要素法では実効値を反復計算によって求めるため、有限要素法ではこのプロセスはありません。もし有限要素法の経験のある方は、FDTDでは実効値を積分で求めることを忘れないよう注意が必要です。
またこのプロセスは前述した収束状態で行わなければ、物理的に意味がありません。積分(積算)を行う場合には必ず収束性の確認が必要です。

解析結果から求められる値

 
上の表でĒi,j,kなどの変数の上に‾がついているものは積算を行って得られる値です。FDTD法ではメッシュ各点の物理量を時系列で扱っているために、非常に多くの情報を得ることが出来ます。とはいえ、すべての情報を保存するのはディスク容量の無駄でもあり現実的ではありません。
電界によるマイクロ波加熱問題では、電界ベクトルの実効値及びSARのデータが必須で必要に応じて実効Poyntingベクトルのデータを保存します。
一方、アンテナ解析の場合には、電磁界強度の実効値に加え位相のデータが必要です。更にアンテナそのもののインピーダンスを検討する場合には複素電流ベクトルのデータが必要であり、保存するべきデータが多くなります。
電波伝送問題では、DUR(Desired Undesired Ratio)が必要ですが、電界ベクトルの実効値と電界成分の位相データがあれば複素電界を求めることが出来るため、複素数の形でDURを計算できます。
お問い合わせ・資料ダウンロードのお申し込みはこちら
 

お問い合わせ・資料ダウンロード

 

第1章 第4節 解析事例(マルチモードキャビティによる加熱)

実験とシミュレーションの比較

 
実験では一般的に具体的な実験系を構成し実験を行い具体的な結果を得ます。一方シミュレーションでは具体的な系を一旦数値モデルに落とし込み(抽象化)、そのモデルから数値計算によって具体的な結果を得ます。
同じ実験系について実験とシミュレーションを行った場合、結果が異なることはさほど珍しいことではありません。シミュレーションでは理想的な系を想定しますが、実験系では機器の錆び、汚れ、歪みなどが頻繁にあり、励振源も理想的な正弦波を発振する機器は存在しないからです。
適切に行われたシミュレーションで実験と異なる結果が得られる場合は、シミュレーションで抽象化出来ていない物理が存在すると考えられます。前述の実験系の非理想的な条件もその一つであり、Maxwell方程式で抽象化出来ない何らかの現象の可能性もあります。実験とシミュレーションが一致しないところには新しい発見の可能性が埋まっているのです。尤もシミュレーションが適切に行われていることが前提ですが。

解析事例(マルチモードキャビティによる加熱)

 
以降で図化処理等について説明するために例題を紹介します。幅300×奥行200×高さ200[mm]のマルチモードキャビティに2.45GHzのTE10導波管が接続されて、キャビティ中心のワークを加熱するものとします。キャビティ内部は空気で満たされており、ワークの物性はε’-jε” = 4-j2で半径20[mm]、高さ20[mm]の円筒形としました。
境界条件は導波管のキャビティの反対側を吸収境界条件MUR1としてキャビティ面は完全導体壁PECとしました。
解析はマルチモードキャビティ内の電磁波伝搬と定在波形成を定性的に把握するものとします。この目的を考慮してメッシュ数60×50×60を使用して、CFL数=0.5の条件でタイムステップを求め20000ステップの解析を行いました。

図化処理(コンター)

 
シミュレーション結果を文書化する際に頻繁に使用されるのが電磁界の強度やSARを色で表現したコンターです。コンターは平面上の分布を視覚的にかつ定量的に把握できる点で非常に優れています。
コンター使用時の注意点はカラーで出力した画像をモノクロまたは2色、3色刷りの印刷物で使用しないことです。カラーのコンターは赤色に大きな値、青紫色に小さな値を割り当てることが多いです。これをモノトーンにすると大きな値と小さな値共に濃い色で印刷されるため、非常に紛らわしい図表になってしまうのです。(下の図)論文や印刷物がモノトーン印刷の場合それに合わせて最初から図化処理する必要があります。(上の図)

図化処理(ラインコンター)

 
ラインコンターは等値線のことです。コンターと同様、シミュレーション結果の図化処理では頻繁に使用されます。コンターに比べると値の把握が直感的には出来にくいが、モノトーン画像でもコンターで発生したような問題は発生しにくいのが特長です。またコンターよりも値の傾きを等値線の間隔から把握しやすいという大きなメリットもあります。
値の把握には上の図のように線の濃さを変更するケースや、ソフトウェアによっては等高線に小さな文字で値を追加できるものもあります。上記の図は線の濃さを変更した図です。
等高線に値を追加した図を使用する場合、印刷時の大きさを考慮して文字の大きさを調整する必要があります。また.pngや.jpeg等のデータを用いた場合には文字も画像として扱われるために印刷時に判読できない可能性が高くなる点に注意が必要です。

図化処理(ポインティングベクトル、流線)

 
ポインティングベクトル・流線(マルチモードキャビティによる加熱)
電磁波解析ではベクトルの描画は電磁界の成分について描画することが多いですがPoyntingベクトルを描画するとエネルギーの伝搬様式を視覚的に捉えられます。また流線は電磁波解析で描画されることはほとんどありませんが、Poyntingベクトルの流線は加熱問題でエネルギー移動把握の有力なツールです。

 

お問い合わせ・資料ダウンロードのお申し込みはこちら

お問い合わせ・資料ダウンロード