I: 第12章時間領域差分法（FDTD）

FDTD（Finite-Difference Time-Domain）法は、Maxwellの回転方程式を空間と時間で差分化し、電界と磁界を交互に更新する方法です。短いパルスを一度伝搬させると、Fourier変換によって広い周波数範囲の応答を得られることが大きな特徴です。

1　Yee格子電界と磁界を半セルずらし、Stokesの法則を局所的な差分として配置します。

2　時間更新磁界を半時刻、電界を整数時刻に置き、leapfrogで交互に進めます。

3　数値分散格子解像度と伝搬方向によって、数値上の位相速度が変わることを調べます。

4　NS-FDTD差分の分母関数を基準波数へ合わせ、位相誤差を狙って減らします。

5　広帯域検証パルス幅と観測時間を分け、周波数分解能とエネルギー収支を検算します。

I-1 Yeeセル：電界と磁界を半セル・半時刻ずらす

Yee格子では、電界の接線成分をセルの辺、磁界の法線成分をセルの面中心に置きます。この配置により、一つの面を囲む電界の循環が、その面を貫く磁束の時間変化を直接更新します。

\oint_{\partial S}\vec E\cdot d\vec\ell =-\dfrac{d}{dt}\int_S\vec B\cdot\vec n\,dS, \qquad \oint_{\partial S}\vec H\cdot d\vec\ell =\dfrac{d}{dt}\int_S\vec D\cdot\vec n\,dS

時間についても $\vec E$ を整数時刻 $t^n$、$\vec H$ を半時刻 $t^{n+1/2}$ に置き、互いの既知値から交互に更新します。

\vec H^{n+1/2}=\vec H^{n-1/2} -\dfrac{\Delta t}{\mu}\nabla_h\times\vec E^n, \qquad \vec E^{n+1}=\vec E^n +\dfrac{\Delta t}{\varepsilon}\nabla_h\times\vec H^{n+1/2}

時間位置 $n\rightarrow n+1$ 視点の回転角

視点を回し、$E_x,E_y,E_z$ が辺、$H_x,H_y,H_z$ がそれぞれに垂直な面の中心に置かれることを確かめます。時間スライダーは、同じ時刻に全未知量を解くのではなく、半ステップずつ追いかける構造を示します。

I-2 2次元TM$_z$（$E_z$偏波）：Yee格子上で電磁波を進める

ここでいう「2次元」は、ベクトルまで平面内にあるという意味ではありません。$z$ 方向へ場が変化しない

$\partial/\partial z=0$ を仮定し、未知量を $E_z,H_x,H_y$ とする TM$_z$（$E_z$偏波）です。

$E_z$ は $xy$ 計算面に垂直な3次元ベクトルであり、$H_x,H_y$ は計算面内にあります。

\dfrac{\partial H_x}{\partial t} =-\dfrac{1}{\mu}\dfrac{\partial E_z}{\partial y}, \qquad \dfrac{\partial H_y}{\partial t} =\dfrac{1}{\mu}\dfrac{\partial E_z}{\partial x}

\dfrac{\partial E_z}{\partial t} =\dfrac{1}{\varepsilon} \left( \dfrac{\partial H_y}{\partial x} -\dfrac{\partial H_x}{\partial y} \right)

を使います。$E_z$ と $H_x,H_y$ を空間でも時間でも半ステップずらすことで、同じセルの周囲を回る循環と面を貫く量が対応します。

H_x^{n+1/2}=H_x^{n-1/2} -\dfrac{\Delta t}{\mu}\,\delta_yE_z^n, \qquad H_y^{n+1/2}=H_y^{n-1/2} +\dfrac{\Delta t}{\mu}\,\delta_xE_z^n

E_z^{n+1}=E_z^n +\dfrac{\Delta t}{\varepsilon} \left(\delta_xH_y^{n+1/2}-\delta_yH_x^{n+1/2}\right)

下の3D場では、$E_z$ を計算面から立ち上がる青（$-z$）・赤（$+z$）の矢印、$H_x,H_y$ を面内の緑の矢印として、Yee格子上の複数位置へ同時に表示します。薄く反復した3枚の色面は同じ $E_z(x,y)$ であり、「$z$ が変わっても同じ」ことを示します。パルスが誘電体を通ると波長と速度が変わり、境界では反射と透過が生じます。

Courant数 $S=c\Delta t/\Delta x$ 誘電体の比誘電率 $\varepsilon_r$ 外周境界時間発展動かす時刻 $n$（直接移動）

赤：$E_z>0$（$+z$）青：$E_z<0$（$-z$）緑：$H_x,H_y$（面内）正面図の ⊙ / ⊗：手前 / 奥向きの $E_z$

「正面（$+z$から）」は別の2次元図ではなく、同じ時刻・同じ数値場を $+z$ 側から見るカメラ操作です。正面では $E_z$ の矢印が視線と重なるため、手前向きを ⊙、奥向きを ⊗ でも示します。「斜視図」へ戻すと、同じ $E_z$ が面外へ伸びていたことを再確認できます。

3D場を動かして確かめること

1. 斜視図を回して $E_z$ が面外、$H_x,H_y$ が面内にあり、どれも場所ごとに値を持つベクトル場であることを確認します。

2. 一時停止して「半ステップ進む」を押し、$E_z^n$ を使って $H^{n+1/2}$、続いて $H^{n+1/2}$ を使って $E_z^{n+1}$ が更新される順序を追います。半時刻差は物理的な反応の遅れではなく、数値上の保存時刻のずれです。

3. 正面図と斜視図を往復し、時刻や場を変えずに、色面・⊙/⊗と立体矢印が同じ $E_z(x,y)$ を表すことを確かめます。

4. $\varepsilon_r$ を上げ、誘電体内で位相速度 $v=1/\sqrt{\mu\varepsilon}$ が低下することを確認します。

5. 反射境界に切り替えて波が外周で戻る様子を見た後、$S$ を安定上限 $1/\sqrt2$ より大きくし、物理波ではない格子状の増幅を確認します。

\Delta t\le \dfrac{1}{c\sqrt{\dfrac{1}{\Delta x^2}+\dfrac{1}{\Delta y^2}}}, \qquad \Delta x=\Delta y\ \Longrightarrow\ S\le\dfrac{1}{\sqrt2}

Canvasの「吸収層」は、境界へ近づくほど場を減衰させる教育用の簡易モデルです。実解析では、入射角や周波数範囲に対する反射を検証したPMLなどを使います。

I-3 数値分散：格子の向きで位相速度が変わる

格子が粗いと、数値上の位相速度が方向と波長に依存します。最短波長を $\lambda_{\min}$ とすると、まず

N_\lambda=\dfrac{\lambda_{\min}}{\Delta_{\max}}

を確認し、$N_\lambda$ を増やして位相と共振周波数が収束するか調べます。「1波長10セル」は出発点にすぎず、必要な位相精度、材料分散、曲面の階段近似によって増やします。

離散化した真空中の分散関係は、2次元では

\sin^2\!\left(\dfrac{\omega\Delta t}{2}\right) =S^2\left[ \sin^2\!\left(\dfrac{k_x\Delta x}{2}\right) +\sin^2\!\left(\dfrac{k_y\Delta y}{2}\right) \right]

となり、連続系の $\omega=c|\vec k|$ からずれます。格子細分によってこの差が減ることを、到達時刻や位相で検証します。

正方格子 $\Delta x=\Delta y=\Delta$ で伝搬角を $\theta$ とすると、$k_x=k_h\cos\theta$、$k_y=k_h\sin\theta$ です。数値的な位相速度の比は

\dfrac{v_{p,h}}{c} =\dfrac{\omega}{c k_h} =\dfrac{\omega\Delta t}{S\,k_h\Delta}

です。Canvasでは、1波長当たりのセル数 $N_\lambda$、Courant数 $S$、伝搬角を動かします。円から外れるほど方向依存性が強く、波面が長時間の伝搬で歪みます。

1波長当たり $N_\lambda$ Courant数 $S$ 伝搬角 $\theta$

I-4 NS-FDTD：基準波数で数値分散を抑える

通常の中心差分は、微分の分母にそのまま $\Delta x$ と $\Delta t$ を使います。NS-FDTD（Nonstandard FDTD）は、差分の局所配置と陽的な時間更新を保ちながら、分母を波の解析解に合わせた関数へ置き換えます。一次元の均質媒質で基準波数 $k_0$、$\omega_0=ck_0$ を選ぶと、

\delta_x^{\mathrm{NS}}f =\dfrac{f(x+\Delta x/2)-f(x-\Delta x/2)}{\phi_x}, \qquad \phi_x=\dfrac{2}{k_0}\sin\!\left(\dfrac{k_0\Delta x}{2}\right)

\delta_t^{\mathrm{NS}}f =\dfrac{f(t+\Delta t/2)-f(t-\Delta t/2)}{\phi_t}, \qquad \phi_t=\dfrac{2}{\omega_0}\sin\!\left(\dfrac{\omega_0\Delta t}{2}\right)

と置けます。平面波を代入すると、通常FDTDの係数 $S=c\Delta t/\Delta x$ の代わりに

\sin\!\left(\dfrac{\omega_h\Delta t}{2}\right) =C_{\mathrm{NS}}\sin\!\left(\dfrac{k\Delta x}{2}\right), \qquad C_{\mathrm{NS}} =\dfrac{\sin(\omega_0\Delta t/2)}{\sin(k_0\Delta x/2)}

を得ます。$k=k_0$ では $\omega_h=\omega_0$ となるため、選んだ基準波数の位相速度は一次元均質問題で厳密に一致します。

\left.\dfrac{v_{p,h}^{\mathrm{NS}}}{c}\right|_{k=k_0}=1, \qquad \dfrac{v_{p,h}^{\mathrm{NS}}}{c} =\dfrac{2}{S k\Delta x} \sin^{-1}\!\left[C_{\mathrm{NS}}\sin\!\left(\dfrac{k\Delta x}{2}\right)\right]

基準波長のセル数 $N_0=\lambda_0/\Delta x$ Courant数 $S$ 評価波数 $k/k_0$ 伝搬距離 $L/\lambda$

同じ波数をもつ連続波が $L/\lambda$ 波長だけ進む時間で比べると、数値波の累積位相差は

\Delta\phi_h =2\pi\dfrac{L}{\lambda} \left(\dfrac{v_{p,h}}{c}-1\right)

です。$k/k_0=1$ ではNS-FDTDの曲線と厳密値が重なり、右の位相ベクトルも厳密解へ戻ります。そこから評価波数を離すと誤差が再び現れ、通常FDTDとの優劣も条件によって変わります。NS-FDTDは「粗い格子が常に正解になる方法」ではなく、基準波数付近の位相誤差を差分演算子へ織り込む方法です。多次元では伝搬方向、広帯域パルスでは周波数範囲、不均質媒質では局所波数の選び方を追加検証します。

この節はColeの高精度Yee法と、矩形格子に対する分散・安定性解析を出発点にしています。原典は Cole, IEEE T-MTT 45(6), 991-996 (1997) および Kudo, Ohtani and Kashiwa (2002) です。

I-5 パルスとスペクトル：時間幅と観測時間を分ける

FDTDで広帯域特性を求めるときは、短いパルスを入射し、観測波形をFourier変換します。中心時刻 $t_0$ と幅 $\tau$ をもつGaussianパルスを

g(t)=\exp\!\left[-\left(\dfrac{t-t_0}{\tau}\right)^2\right], \qquad |G(\omega)|\propto\tau\exp\!\left[-\left(\dfrac{\omega\tau}{2}\right)^2\right]

とすると、パルスを短くするほど励起帯域は広がります。一方、周波数分解能は観測時間 $T_{\mathrm{obs}}$ で決まります。

\Delta f=\dfrac{1}{T_{\mathrm{obs}}}, \qquad f_{\mathrm{Nyquist}}=\dfrac{1}{2\Delta t}

パルス幅 $\tau/\Delta t$ 観測ステップ数 $N_t$

パルス幅を変えると連続スペクトルの広がりが変わり、観測ステップ数を変えると評価できる周波数点の間隔が変わります。「広帯域を励起する」と「共振ピークを細かく読む」は、別の数値設定です。

FDTDで必ず検証する量

CFL条件を満たし、時間刻みを小さくしても結果が変わらないこと。
格子を細かくし、共振周波数、反射係数、到達時刻が収束すること。
入射電力、反射電力、透過電力、吸収電力の収支が合うこと。
PMLを厚くしても観測量が変わらず、境界反射が十分小さいこと。
導電率・Debye・Lorentzなどの材料モデルが受動性と因果性を壊していないこと。

G: 有限要素法へ　|　H: モーメント法へ　|　教材ポータルへ戻る　|　電磁界解析ソフトウェアを見る

I: 第12章 時間領域差分法（FDTD）