JP3878148B2 - Electromagnetic field simulation method, electromagnetic field simulation program, and program recording medium - Google Patents

Description

（ここで、係数ｆ_ij（ｉ,ｊ＝０,１,２）は、
【数１２】

係数ｇ_ijk（ｉ＝０,１；ｊ,ｋ＝０,１,２）は、
【数１３】

係数ｈ_ij（ｉ,ｊ＝０,１）は、
【数１４】

と定義される）を満たすように誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε₁ , ε₂ を設定し、メモリに格納するステップと、(b) 離散化された時間を表す時刻ｎの初期値を設定するとともに、各格子点における電磁界の電界成分と磁界成分の初期値を設定し、メモリに格納するステップと、 (c) 各格子点に設定された誘電率と磁界成分をメモリから読み出し、時刻ｎ＋ １／２における各格子点の電界成分を計算し、メモリに格納するステップと、 (d) 電界成分をメモリから読み出し、時刻ｎ＋１における各格子点の磁界成分を計算し、メモリに格納するステップと、 (e) 時刻ｎを１増加するステップと、 (f) 時刻ｎを所定の閾値と比較し、時刻ｎが所定の閾値に達していない場合には、ステップ (c) に戻り、時刻ｎが所定の閾値に達していた場合には、電界成分と磁界成分をメモリから読み出し、出力するステップとを実行することを要旨とする。
【００１１】
第２の本発明に係る電磁界シミュレーション方法は、一つの平面内に誘電率ε₁ 、ε₂ をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求める電磁界シミュレーション方法であって、第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε₁ を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１より大きく１より小さい範囲の任意の実数値をとる定数δおよび前記格子間隔Δを用いて（１／２＋δ）Δと表される場合、前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｕ、前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｖ、前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｗ、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε₁の媒質に対応する格子点の誘電率をε₁ 、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε₂の媒質に対応する格子点の誘電率をε₂として、前記コンピュータが、(a)関係式
【数１５】

を満たすように誘電率ｕ，ｖ，ｗ，ε₁ , ε₂ を設定し、メモリに格納するステップと、(b) 離散化された時間を表す時刻ｎの初期値を設定するとともに、各格子点における電磁界の電界成分と磁界成分の初期値を設定し、メモリに格納するステップと、 (c) 各格子点に設定された誘電率と磁界成分をメモリから読み出し、時刻ｎ＋１／２における各格子点の電界成分を計算し、メモリに格納するステップと、 (d) 電界成分をメモリから読み出し、時刻ｎ＋１における各格子点の磁界成分を計算し、メモリに格納するステップと、 (e) 時刻ｎを１増加するステップと、 (f) 時刻ｎを所定の閾値と比較し、時刻ｎが所定の閾値に達していない場合には、ステップ (c) に戻り、時刻ｎが所定の閾値に達していた場合には、電界成分と磁界成分をメモリから読み出し、出力するステップとを実行することを要旨とする。
【００１２】
上記電磁界シミュレーション方法において、２次精度を有する時間差分式、および２次精度を有する空間差分式を用いた時間領域差分法によって離散化された時間の変化に応じた電磁界を順次求めることを要旨とする。
【００１３】
第３の本発明に係る電磁界シミュレーションプログラムは、一つの平面内に誘電率ε₁ 、ε₂ をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求めるために、第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε₁ を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１.５より大きく１.５より小さい範囲の任意の実数値をとる定数ｓおよび前記格子間隔Δを用いてｓΔと表される場合、前記境界を表す直線と平行であり、なおかつｙ切片が−Δである直線上に配置される第３種格子点の誘電率をａ、前記境界を表す直線と平行であり、なおかつｙ切片が−Δ／２である直線上に配置される第２種格子点の誘電率をｂ、前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｃ、前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｄ、前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｅ、前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε₁の媒質に対応する格子点の誘電率をε₁ 、前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε₂の媒質に対応する格子点の誘電率をε₂として、前記コンピュータに、(a)関係式
【数１６】

（ここで、係数ｆ_ij（ｉ,ｊ＝０,１,２）は、
【数１７】

係数ｇ_ijk（ｉ＝０,１；ｊ,ｋ＝０,１,２）は、
【数１８】

係数ｈ_ij（ｉ,ｊ＝０,１）は、
【数１９】

と定義される）を満たすように誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε₁ , ε₂ を設定し、メモリに格納するステップと、(b) 離散化された時間を表す時刻ｎの初期値を設定するとともに、各格子点における電磁界の電界成分と磁界成分の初期値を設定し、メモリに格納するステップと、 (c) 各格子点に設定された誘電率と磁界成分をメモリから読み出し、時刻ｎ＋ １／２における各格子点の電界成分を計算し、メモリに格納するステップと、 (d) 電界成分をメモリから読み出し、時刻ｎ＋１における各格子点の磁界成分を計算し、メモリに格納するステップと、 (e) 時刻ｎを１増加するステップと、 (f) 時刻ｎを所定の閾値と比較し、時刻ｎが所定の閾値に達していない場合には、ステップ (c) に戻り、時刻ｎが所定の閾値に達していた場合には、電界成分と磁界成分をメモリから読み出し、出力するステップとを実行させることを要旨とする。
【００１４】
第４の本発明に係る電磁界シミュレーションプログラムは、一つの平面内に誘電率ε₁ 、ε₂ をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求めるために、第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε₁ を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１より大きく１より小さい範囲の任意の実数値をとる定数δおよび前記格子間隔Δを用いて（１／２＋δ）Δと表される場合、前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｕ、前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｖ、前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｗ、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε₁の媒質に対応する格子点の誘電率をε₁ 、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε₂の媒質に対応する格子点の誘電率をε₂として、前記コンピュータに、(a)関係式
【数２０】

を満たすように誘電率ｕ，ｖ，ｗ，ε₁ , ε₂ を設定し、メモリに格納するステップと、(b)離散化された時間を表す時刻ｎの初期値を設定するとともに、各格子点における電磁界の電界成分と磁界成分の初期値を設定し、メモリに格納するステップと、(c)各格子点に設定された誘電率と磁界成分をメモリから読み出し、時刻ｎ＋１／２における各格子点の電界成分を計算し、メモリに格納するステップと、(d)電界成分をメモリから読み出し、時刻ｎ＋１における各格子点の磁界成分を計算し、メモリに格納するステップと、(e)時刻ｎを１増加するステップと、(f)時刻ｎを所定の閾値と比較し、時刻ｎが所定の閾値に達していない場合には、ステップ(c)に戻り、時刻ｎが所定の閾値に達していた場合には、電界成分と磁界成分をメモリから読み出し、出力するステップとを実行させることを要旨とする。
【００１５】
上記電磁界シミュレーションプログラムは、２次精度を有する時間差分式、および２次精度を有する空間差分式を用いた時間領域差分法によって離散化された時間の変化に応じた電磁界を順次求めることを要旨とする。
【００１６】
第５の本発明に係るプログラム記録媒体は、上記電磁界シミュレーションプログラムのいずれかを記録したことを要旨とする。
【００１７】
この発明におけるプログラム記録媒体とは、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＤＶＤ（Digital Versatile Disk）、光磁気ディスク、ＰＣカード等のコンピュータ読み取り可能な記録媒体を意味する。
【００１８】
【発明の実施の形態】
次に、図面を参照して本発明の実施の形態を説明する。
【００１９】
図１は、本発明の一実施形態に係る電磁界シミュレーション方法を実行する電子的な装置である電磁界シミュレーション装置の構成を示す機能ブロック図である。同図に示す電磁界シミュレーション装置１は、シミュレーション用基礎データと各種パラメタを入力するためのキーボード、マウス等の入力装置から成る入力部１１、入力された基礎データや各種パラメタを用いて電磁界シミュレーションを行うために、演算機能および制御機能を備えた中央処理装置（ＣＰＵ：Central Processing Unit）から構成される演算部１２、入力部１１からの入力情報や演算部１２からの演算結果等を記憶する記憶部１３、記憶部１３で記憶した情報をシミュレーション結果として出力するためのディスプレイ装置等から成る出力部１６を少なくとも備えている。
【００２０】
記憶部１３は、ＲＡＭ（Random Access Memory）等から構成される主記憶部１４と、ハードディスクドライブ、フレキシブルディスクドライブ、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）ドライブ、ＤＶＤ（Digital Versatile Disk）ドライブ、光磁気ディスクドライブ、ＰＣカードドライブ等から構成される補助記憶部１５とを備えており、後述する計算結果を随時格納するために必要なメモリ領域が確保されている。
【００２１】
以後に説明する本実施形態の各種処理は、一つの電子的な装置が実行する場合だけでなく、各ステップの実行を適宜分割して二つ以上の電子的な装置から構築されたシステムが全体で実行する場合も含む。そこで、本実施形態においては、一つまたは複数の電子的な装置から成る装置（コンピュータ）を「電磁界シミュレーション装置」と呼ぶ。この意味で、本実施形態には、電磁界シミュレーション処理を行うためのプログラムの指示にしたがって、例えば電磁界シミュレーション装置１上で稼動するＯＳ（Operating System）が実際の処理を行う場合等が含まれることはいうまでもない。この点は、本発明の全ての実施の形態に共通する事項である。
【００２２】
図２は、本実施形態の電磁界シミュレーションを実行する際に、誘電率が不連続的に異なる二つの領域（誘電率をε₁ 、ε₂ とする）が直線状の境界で隣接している連続媒質の平面を離散化して得られる格子状の仮想平面としての２次元Yee格子（以後、単にYee格子と記載する）の概略構成を示す説明図である。同図に示すYee格子２０は、３種類の異なるタイプの格子点２１（第１種格子点）、格子点２２（第２種格子点）、および格子点２３（第３種格子点）が交互に２次元的に配置され、各々が正方格子を単位格子とする２次元格子（同種の格子点間の格子間隔である格子定数をΔとする。ここでΔは正の定数）を構成することによって離散化された平面（仮想平面）を形成している。より具体的には、格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、格子点２２が最近接する二つの格子点２１の中点に位置し、なおかつ格子点２３が最近接する二つの格子点２１の中点のうち格子点２２が位置する中点とは異なる中点に位置するように配置した上で、格子点２１に前記平面と垂直な磁界成分を割り当てることによってその離散化された平面を構成している。したがって、このYee格子２０において、最近接する異種格子点間の距離はΔ／２である。
【００２３】
図２では、格子点２１と格子点２３とが交互に配置される列に平行な方向（図の水平方向）をｘ軸方向とし、この列のうちの一つをｘ軸としている。他方、格子点２１と格子点２２とが交互に配置される列に平行な方向（図の鉛直方向）をｙ軸とし、この列のうちの一つをｙ軸としている。このような座標軸の原点には格子点２１のいずれかが配置される。図２では、以上の設定を満たす座標系としてΣ₀を取った場合を図示しているが、この座標系の取り方が任意であることはいうまでもない。
【００２４】
図２に示す座標系Σ₀を用いると、任意の整数ｉおよびｊを用いることにより、格子点２１の座標を（ｉΔ,ｊΔ）、格子点２２の座標を（ｉΔ,（ｊ＋１／２）Δ）、格子点２３の座標を（（ｉ＋１／２）Δ,ｊΔ）とそれぞれ表すことができる。以後の説明においては、記述を簡単にするために、Δを明記しなくともその意味が十分理解できるとみなしうる場合には、各格子点のｘ、ｙ座標をそれぞれΔで割り、単なる整数ｉ、ｊの組み合わせとして、例えば上記格子点２１、２２、２３の座標をそれぞれ（ｉ,ｊ）、（ｉ,ｊ＋１／２,ｊ）、（ｉ＋１／２,ｊ）と表すことにする。
【００２５】
本実施形態においては、格子点２１には図２の平面に垂直な方向（ｚ方向）の磁界成分Ｈ_z（ｉ,ｊ）を割り当て、格子点２２には、ｘ方向の電界成分Ｅ_x（ｉ,ｊ＋１／２）を、格子点２３には、ｙ方向の電界成分Ｅ_y（ｉ＋１／２,ｊ）をそれぞれ割り当てる。これにより、Yee格子面に垂直な方向への振動を磁界のみが有する２次元ＴＭ（Transverse Magnetic）モードの電磁界をモデル化できる。ちなみに、この２次元ＴＭモードでは、ｚ方向の電界成分Ｅ_z 、ｘ方向の磁界成分Ｈ_x 、およびｙ方向の磁界成分Ｈ_yは、時間によらず常に０（一定）である。
【００２６】
二つの媒質の境界２４は、図２にも示すように、格子の配列方向であるｘ軸およびｙ軸方向に対して所定の傾斜角をなす方向を指向する。より具体的には、ｙをｘの１次関数として境界２４を表現する場合、その直線の傾きの大きさが２である。そこで、以後の説明の便宜上、この境界をなす直線の傾きを−２とし、Δおよびｓ（−１.５＜ｓ＜１.５）を用いて
ｙ＝−２ｘ＋ｓΔ （１）
とするが、本実施形態の電磁界シミュレーション方法は、傾きが＋２であってもよいし、ｙ切片が上記以外の値をとってもかまわない。これらは、座標系の取り方の任意性にも依存するものである。しかしながら、以上のような限定を加えても、本実施形態に係る電磁界シミュレーション方法の一般性が失われないことはいうまでもない。
【００２７】
式（１）で表される境界２４の近傍に配置される格子点２２および２３（ともにゼロでない電界成分を有する格子点）の誘電率は、次のように設定する：
（Ｉ-１）直線ｙ＝−２ｘ−Δ上に配置される格子点２３の誘電率をａとする。
（Ｉ-２）直線ｙ＝−２ｘ−Δ／２上に配置される格子点２２の誘電率をｂとする。
（Ｉ-３）直線ｙ＝−２ｘ上に配置される格子点２３の誘電率をｃとする。
（Ｉ-４）直線ｙ＝−２ｘ＋Δ／２上に配置される格子点２２の誘電率をｄとする。
（Ｉ-５）直線ｙ＝−２ｘ＋Δ上に配置される格子点２３の誘電率をｅとする。
【００２８】
これら以外の格子点２２および２３の誘電率は、次のように設定する：
（Ｉ-６）ｙ＜−２ｘ＋ｓΔの領域（この領域を領域Ｄ₁ とする）に配置される上記以外の格子点２２および２３の誘電率をε₁ とする。
（Ｉ-７）ｙ＞−２ｘ＋ｓΔの領域（この領域を領域Ｄ₂ とする）に配置される上記以外の格子点２２および２３の誘電率をε₂ とする。
【００２９】
なお、透磁率は、媒質に関わらず一定であるとする。
【００３０】
図３は、以上の誘電率の割り当てを一つの表３０にまとめた図である。同図に示す表３０において、例えば「ｙ＝−２ｘ−Δ」と記載されている列は、この直線上に配置される各格子点に割り当てられる誘電率を意味している。すなわち、上記（Ｉ-１）より、格子点２３には誘電率ａが割り当てられる。他方、直線ｙ＝−２ｘ−Δ上には格子点２２が存在しないため、表３０では格子点２２の誘電率を記載する欄に横線（−）を引いている。他の欄の記載についても同様である。
【００３１】
図４は、境界２４付近の各格子点に割り当てられる誘電率を概念的に示す説明図である。なお、図２と同じ記号によって記載されており、割り当てられている誘電率が明記されていない格子点２２および２３については、上記（Ｉ-６）および（Ｉ-７）にしたがって、適宜ε₁またはε₂が割り当てられているものとする。
【００３２】
領域Ｄ₁側から領域Ｄ₂側へ平面波が入射した時、磁界振幅に基づく境界での反射率Ｒおよび透過率Ｔは、格子定数Δの級数として、次のように表すことができる。
【数２１】

ここで、θ_iは入射角、θ_tは反射角をそれぞれ表している。境界２４を誘電率ε₁の領域と誘電率ε₂の領域の境界とすると、式（２）および（３）の右辺第１項はともに真値である。
【００３３】
本実施形態では、１次の誤差項の係数Ｒ₁およびＴ₁を入射角θ_iの値に関わらずゼロとする２次精度のシミュレーションを行う。そこで、式（２）および（３）で、Δの１次の係数が入射角θ_iの値に関わらずにゼロとなる条件を求めると、次の式が得られる。
【数２２】

ここで、式（４）の係数ｆ_ij（ｉ,ｊ＝０,１,２）、式（５）の係数ｇ_ijk（ｉ＝０,１；ｊ,ｋ＝０,１,２）、および式（６）の係数ｈ_ij（ｉ,ｊ＝０,１）は、上記（Ｉ-１）〜（Ｉ-５）に基づいて設定される誘電率ａ,ｂ,ｃ,ｄ,ｅを用いてそれぞれ次のように表される。
【００３４】
【数２３】

【００３５】
【数２４】

【００３６】
【数２５】

【００３７】
ちなみに、領域Ｄ₂から領域Ｄ₁へ平面波が入射した場合に２次精度を与える条件式を求めても、式（４）〜（６）と同じ式が得られる。
【００３８】
以上説明した誘電率の設定法を、「誘電率設定法Ｉ」と呼ぶ。
【００３９】
ところで、この誘電率設定法Ｉにおいて、さらに、
【数２６】

とおき、これらをｕ,ｖ,ｗについて解くことにより、次の式（７）〜（９）を得る。
【数２７】

【００４０】
これらの式（７）〜（９）を用いる場合には、境界２４を表す直線を、
【数２８】

とし、格子点２２および２３の誘電率を、次のように割り当てる：
まず、上記（Ｉ-３）〜（Ｉ-５）のｃ，ｄ，ｅをそれぞれｕ，ｖ，ｗとし、
（II-１）直線ｙ＝−２ｘ上に配置される格子点２３の誘電率をｕとする。
（II-２）直線ｙ＝−２ｘ＋Δ／２上に配置される格子点２２の誘電率をｖとする。
（II-３）直線ｙ＝−２ｘ＋Δ上に配置される格子点２３の誘電率をｗとする。
さらに、
（II-４）上記３つの条件に含まれない格子点２２および２３のうち、ｙ＜−２ｘ＋（１／２＋δ）Δの領域（この領域を領域Ｄ₁とする）に配置される格子点２２および２３の誘電率をε₁ とする。
（II-５）上記４つの条件に含まれない格子点２２および２３の誘電率をε₂ とする。この誘電率ε₂ が割り当てられる領域は、ｙ＞−２ｘ＋（１／２＋δ）Δの領域（この領域を領域Ｄ₂とする）内に含まれる。
【００４１】
条件（II-１）〜（II-５)および式（７）〜（９）に基づく誘電率設定法を「誘電率設定法II」と呼ぶことにする。この誘電率設定法IIが、誘電率設定法Ｉの特別な場合であることは、以上の説明からも明らかである。
【００４２】
図５は、上記各条件を一つの表４０にまとめて記載した図であり、表４０の見方は図３と同様である。
【００４３】
また、図６の説明図には、境界２４付近の各格子点に割り当てられる誘電率を概念的に示した。同図においても、図２と同じ記号で記載され、割り当てられている誘電率が明記されていない格子点２２および２３には、上記（II-４）および（II-５）にしたがって、適宜ε₁またはε₂が割り当てられているものとする。
【００４４】
電磁界の各成分は、以上説明したように、空間的に異なった位置に割り当てられるが、時間的な意味においても、Yee格子２０の平面に時間の次元を加えて成る３次元時空間の中の異なった位置に割り当てられる。例えば、磁界が、離散化された時刻ｔ¹,ｔ²,…,ｔ^n-1,ｔⁿ,ｔⁿ⁺¹,…（ｎは正の整数）に割り当てられるとすると、電界はｔ^1/2,ｔ^1+1/2,ｔ^2+1/2,…,ｔ^n-1/2,ｔ^n+1/2,ｔ^n+3/2,…というように、磁界成分が割り当てられた時刻の中間の時刻に割り当てられる。このように時刻を設定することによって、磁界と電界を交互に計算することが可能になる。なお、電界と磁界に割り当てられる時刻を逆に定義することも勿論可能である。
【００４５】
次に、以上の設定に基づいた電磁界シミュレーション方法について詳細に説明する。
【００４６】
本実施形態で適用する時間２次精度、空間２次精度の時間領域差分法では、電磁界シミュレーション装置１の記憶部１３内のメモリに記憶された電磁界成分から、時系列的な意味で次の時刻の電磁界が以下に示す演算によって順次求められる。
【００４７】
【数２９】

ここで、Δｔは時間ステップ間隔、Ｅ_x ^n+1/2（ｉ,ｊ＋１／２）は、格子点２２（ｉ,ｊ＋１／２）の時刻ｔ^n+1/2 における電界のｘ方向成分、Ｅ_y ^n+1/2（ｉ＋１／２.ｊ）は格子点２３の時刻ｔ^n+1/2 における電界のｙ方向成分、Ｈ_z ⁿ⁺¹（ｉ,ｊ）は格子点２１の時刻ｔⁿ⁺¹ における磁界のｚ方向成分である。また、ε（ｉ,ｊ＋１／２）は座標（ｉ,ｊ＋１／２）の格子点２２に付与される誘電率、ε（ｉ＋１／２,ｊ）は座標（ｉ＋１／２,ｊ）の格子点２３に付与される誘電率である。さらに、μは透磁率であり、前述したように媒質に関わらず一定値をとる。
【００４８】
なお、各式の右辺にある［ ］は、その［ ］内に記載された偏微分を近似する空間２次精度の差分式を意味しており、より具体的には以下の差分式で定義される：
【数３０】

【００４９】
媒質の誘電率が一様な場合のように、各電磁界成分が少なくとも２回以上偏微分可能としてよい場合は、これらの差分式が空間２次精度を与える事を各電磁界成分の空間座標（ｘ,ｙ）についてのテイラー展開を用いて示す事ができる。一方、境界２４の近傍では、各電磁界成分が境界２４と垂直な方向（ｘ方向）に対して一般に２回偏微分可能であるとは限らない。特に、境界２４に垂直な方向の電界成分は不連続であるため、これまで高精度の計算が可能であるかどうかは自明ではなかった。
【００５０】
これに対して、本実施形態においては、境界２４が、座標系Σ₀において傾き２または−２の直線で表される場合、境界２４近傍の格子点２２、２３に割り当てる誘電率を条件式（４）〜（６）または式（７）〜（９）を満たすように設定し、通常の２次精度の空間差分式（式（１４）〜（１７）を参照）を用いることにより、境界２４に入射する電磁波の入射角に関わらず境界２４での反射率と透過率を２次精度で近似可能な２次元ＴＭモードの電磁界シミュレーションを実現することができる。
【００５１】
図７は、本実施形態に係る電磁界シミュレーション方法として、時間２次精度、空間２次精度の時間領域差分法を用いる際の具体的な処理手順を示すフローチャート図である。同図に示す処理を行うにあたっては、シミュレーション領域の格子点への誘電率の設定等、シミュレーション用基礎データおよび各種パラメタが予め入力部１１から入力され、記憶部１３に記憶されており、これらの入力内容に基づいて、誘電率の具体的な設定が演算部１２によってなされているものとする。
【００５２】
図７に示す時間領域差分法においては、まず、電界と磁界の初期値を設定する（ステップＳ１０１）。ここで、初期磁界の時刻をｔ＝−１／２、初期電界の時刻をｔ＝０とする。
【００５３】
次に、初期時刻（時間ステップ）をｎ＝０と設定する（ステップＳ１０２）。いうまでもなく、図７の記号「：＝」は、「ｎを０と定義する」ことを意味している。以後の時間ステップは、上述したように、電界がｔ^1/2,ｔ^1+1/2,ｔ^2+1/2,…,ｔ^n-1/2,ｔ^n+1/2,ｔ^n+3/2,…、磁界がｔ¹,ｔ²,…,ｔ^n-1,ｔⁿ,ｔⁿ⁺¹,…と推移して時間が経過していく。
【００５４】
この後、時間の経過にしたがって電磁界成分を演算部１２において順次計算する。そこで、以後のステップは、一般のｎ（正の整数）で表される任意の時間ステップについて説明する。
【００５５】
最初に、式（１１）、（１２）を用いて時間ステップｔ＝ｎ＋１／２の電界の各成分Ｅ_x ^n+1/2（ｉ,ｊ＋１／２）、Ｅ_y ^n+1/2（ｉ＋１／２,ｊ）を計算し、その計算結果を格子点と時間ステップの組（３次元時空間の座標に相当）に対応づけて記憶部１３内のメモリに格納する（ステップＳ１０３）。
【００５６】
式（１３）を用いて時間ステップｔ＝ｎ＋１の磁界成分Ｈ_z ⁿ⁺¹（ｉ,ｊ）を計算し、その計算結果をステップＳ１０３と同様にメモリに格納する（ステップＳ１０４）。
【００５７】
電磁界成分の計算結果は、順次記憶部１３内のメモリに格納されて記憶されることになるが、前の時刻（時間ステップ）の電界または磁界を、次の時刻（時間ステップ）の電界または磁界に置換しても、以後の電磁界の計算に支障は無い。なぜならば、電界と磁界の計算は、式（１１）〜（１３）からも明らかなように、時間ステップごとに交互に実行されるからである。このようにして電磁界成分を順次更新することにより、計算時に用いるメモリ領域を削減して、負荷を軽減することが可能になる。
【００５８】
続いて、ステップＳ１０３およびＳ１０４でそれぞれ求めた電界成分Ｅ_x ^n+1/2（ｉ,ｊ＋１／２）、Ｅ_y ^n+1/2（ｉ＋１／２,ｊ）、および磁界成分Ｈ_z ⁿ⁺¹（ｉ,ｊ）から、それらに関するデータを抽出してメモリへ格納する（ステップＳ１０５）。ちなみに、この抽出された電磁界成分に関するデータは、後述するステップにおいて電磁界の各モードの強度や位相等のデータを計算する際に使用することもできる。この意味で、これらのデータを格納するメモリ領域は、電磁界成分自体を格納するメモリ領域とは異なる領域であることが好ましい。なお、実行するシミュレーションの内容によっては、このステップＳ１０５が必須ではなく、ステップＳ１０４終了後、以下に説明するステップＳ１０６に推移するようにしてもよい。
【００５９】
この後、ｎの値を１増やして時間ステップを進め（ステップＳ１０６）、ステップＳ１０７において、ｎが予め設定した閾値ｎ_maxに達していない（ｎ＜ｎ_max）と判定した場合には、ステップＳ１０３に戻って上述したステップＳ１０３〜Ｓ１０６の処理を繰り返す。なお、ステップＳ１０７の判定が演算部１２で行われることはいうまでもない。
【００６０】
他方、ステップＳ１０７でｎがｎ_maxに達したと判定した場合には、電磁界分布を出力後、電磁界の振幅や位相、反射率等を計算して計算結果を出力部１６から出力する（ステップＳ１０８）。
【００６１】
以上の処理によって得られるシミュレーション結果、およびシミュレーションのために入力される基礎データおよび各種パラメタ等は、必要に応じて記憶部１３から読み出して出力部１６から出力させることができる。
【００６２】
以上説明したように、本発明の一実施形態によれば、境界が、所定の座標系を用いてある傾きを有する直線で表される場合、境界近傍の格子点に対して条件式（４）〜（６）または式（７）〜（９）を満たすような誘電率を設定し、時間２次精度、空間２次精度の時間領域差分法を用いることにより、境界に入射する電磁波の入射角に関わらず境界での反射率と透過率を２次精度で近似可能な２次元ＴＭモードの電磁界シミュレーションを実現することができる。
【００６３】
なお、本実施形態の電磁界シミュレーション処理を実現するための電磁界シミュレーションプログラムを記録したコンピュータ読み取り可能なプログラム記録媒体をコンピュータに装着し、そのプログラム記録媒体に格納されているプログラムを読み出すことによって、コンピュータが上述した処理を実行するようにしてもよい。ここで、「コンピュータ読み取り可能な」プログラム記録媒体としては、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＤＶＤ、光磁気ディスク、ＰＣカード等を用いることができる。このようなプログラム記録媒体を提供することによって、本実施形態の電磁界シミュレーションプログラムを広く流通させることができるようになる。
【００６４】
（その他の実施形態）
ところで、本発明は、上述した一実施形態においてのみ特有の効果を奏するものと理解されるべきではない。
【００６５】
例えば、境界近傍の誘電率を設定するに際して、より広範囲の境界近傍の格子点に対して、もとの連続媒質における誘電率（ε₁ 、ε₂)とは異なる誘電率を割り当てることができる。具体的には、境界近傍において、その境界をなす直線と平行でｙ切片が異なるもののうち、境界をなす直線との距離が短い順にｍ個（ｍは２以上の任意の整数）の直線上を通過する格子点に対し、その直線ごとに誘電率の設定を行うことができる（誘電率設定法Ｉがｍ＝５の場合であり、誘電率設定法IIがｍ＝３の場合に相当）。
【００６６】
また、境界をなす直線のYee格子の軸方向に対する傾きは、上述した場合に限定されるわけではなく、軸方向に対して他の傾きを有する境界に対しても、同様に電磁界シミュレーションを実行することが可能である。
【００６７】
さらに、正方格子ではないYee格子を用いて電磁界シミュレーションを行うことも勿論可能である。
【００６８】
このように、本発明は、上記一実施形態と同様に、時間２次精度、空間２次精度の時間領域差分法を用いることによって高精度の電磁界シミュレーションを実現することのできる様々な実施の形態等を含みうるものである。
【００６９】
【実施例】
以下、本発明の一実施形態に係る電磁界シミュレーション方法を用いた電磁界シミュレーションの実施例について説明する。
【００７０】
＜シミュレーション領域＞
図８は、本実施例において適用されるシミュレーション領域の概要を示す説明図である。なお、離散化の対象となる連続媒質の２次元平面は、誘電率がε₂＝10.24の領域である半導体スラブ導波路５１と誘電率がε₁＝1.0の領域であるバックグラウンド５２が、図８に示すように配置されたものである。
【００７１】
図８に示すシミュレーション領域５０は長方形をなしており、半導体スラブ導波路５１のバックグラウンド５２との境界では、境界をなす直線が、その境界付近に設けられる適当な局所座標系において、式（１０）を満たすようになっている。すなわち、境界５３では、図中鉛直下向きをｘ軸、水平方向左向きをｙ軸、図面に垂直で図の裏側に向かう方向をｚ軸とする局所座標系Σ₁を取る一方で、境界５４では、図中鉛直上向きをｘ軸、水平方向右向きをｙ軸、図面に垂直で図の裏側に向かう方向をｚ軸とする局所座標系Σ₂を取ることにより、各境界５３および５４が式（１０）の直線として表される。
【００７２】
シミュレーション領域５０内には図２と同様のYee 格子２０が形成されており、局所座標系Σ₁もしくはΣ₂に基づいて各点の座標を定めたとき、磁界成分Ｈ_z（ｉ,ｊ）を割り当てる格子点２１（ｉ,ｊ）、電界成分Ｅ_x（ｉ,ｊ＋１／２）を割り当てる格子点２２（ｉ,ｊ＋１／２）、電界成分Ｅ_y（ｉ＋１／２,ｊ）を割り当てる格子点２３（ｉ＋１／２,ｊ）が交互に配列され、各格子点はそれぞれ同じ格子定数Δを有する正方格子をなしている。
【００７３】
＜シミュレーション用基礎データおよび各種パラメタの設定＞
シミュレーション用基礎データおよび各種パラメタについて説明する。本実施例では、各格子点に割り当てる誘電率の設定を誘電率設定法IIに基づいて設定するものとする。
【００７４】
境界５３の局所座標系Σ₁における直線の式（１０）のδの値を０とし、シミュレーション領域５０の対角線の長さを、伝搬シミュレーションを行う固有モードの整数倍、すなわち（２π／［伝搬シミュレーションを行うモードの伝搬定数］）の整数倍とした。また、シミュレーション領域５０の周囲には周期境界条件を課し、半導体スラブ導波路５１の幅を１.５μｍ（１.５×１０^-6 ｍ）とした。
【００７５】
これに対して、境界５４の局所座標系Σ₂における直線の式（１０）のδの値は、格子定数Δの値に応じて若干変化するが、その変化率がたかだか２％程度であるようにδの値を定めることができる。このδの具体的な値の例については後述する。
【００７６】
なお、以上の設定とは逆にして、境界５４の局所座標系Σ₂における直線の式（１０）のδの値を０と先に定めた後、境界５３の局所座標系Σ₁における直線の式（１０）のδの値を定めても勿論かまわない。
【００７７】
次に、伝搬シミュレーションを行う電磁波の設定について説明する。
【００７８】
まず、伝搬する電磁波の自由空間波長を１.５５μｍとした。この場合、磁界成分Ｈ_z 、および二つの方向の電界成分Ｅ_x ,Ｅ_y によって表される電磁波の伝搬モードは、６つのＴＭモード（０次〜５次）を保持することができる。各モードは、半導体スラブ導波路５１を伝搬する波の導波路境界への入射角が異なる。本実施例では、入射角が６１.４度である２次モードと入射角が２２.１度である５次モードについて伝搬シミュレーションを行った。
【００７９】
各格子点に割り当てられる電磁界成分の初期値については、解析的に得られる結果を初期値として各格子点に割り当てた。
【００８０】
図９は、半導体スラブ導波路５１とバックグラウンド５２の境界５３または５４近傍の格子点２２または２３に付与される誘電率の設定例を示す図である。同図に示す表６０においては、δ＝0.0となる境界（境界５３、５４のうち少なくとも１つの境界ではこの値をとる）において、ｕ＝0.82636, ｖ＝1.9393, ｗ＝16.061とした。
【００８１】
他方、δ≠０となる境界において、半導体スラブ導波路５１内の波長あたりのセル数（以後、導波路内波長あたりのセル数）が23.195でδ＝0.4の場合、ｕ＝0.81793, ｖ＝0.96931, ｗ＝7.1722となる。表６０では、δ＝0.4のところに、括弧書きで0.39947(=0.4-0.0053)という値が記載されているが、これは、δ≠０の場合には、格子定数Δ（導波路内波長あたりのセル数は格子の細かさの度合いを表す指標である）の値に応じて若干変化するためである。しかしながら、この変化は上述したようにたかだか２％程度なので、これが本実施例の電磁界シミュレーション結果に対して特段の影響を及ぼすことはない。
【００８２】
また、半導体スラブ導波路５１を構成する導波路内波長あたりのセル数が17.503でδ＝-0.2の場合、ｕ＝0.96970, ｖ＝2.74717, ｗ＝23.4416となる。この場合に表６０に括弧書きで記載されている値-0.19828(=-0.2+0.00172）は、導波路内波長あたりのセル数が17.503のときにシミュレーションに用いられる正確な値を示している。
【００８３】
ｕ,ｖ,ｗ以外の格子点の誘電率については、バックグラウンド５２を真空としてその誘電率ε₁を1.0とする一方、上記以外の半導体スラブ導波路５１の誘電率ε₂は10.24とした。
【００８４】
なお、本実施例のシミュレーションにおいては、Yee格子による離散化を行う際に、空間のスライス幅に対して時間のスライス幅をどの程度にするかの指標を与えるＣＦＬ数（Courant-Friedrichs-Levy number）を0.6とした。
【００８５】
＜シミュレーションの内容＞
以上説明したシミュレーション用基礎データおよび各種パラメタの設定に基づいて、２次元ＴＭモードのうちの２次モードと５次モードにつき、さまざまな格子定数Δを有するYee格子に対してシミュレーションを行った。
【００８６】
反射率Ｒおよび透過率Ｔの精度は、電磁波の位相速度の精度により評価できるので、本実施例では、この位相速度の誤差を調べた。
【００８７】
シミュレーション結果としての位相速度は、伝搬シミュレーションの電磁界分布から位相情報を抽出することによって求められる（図７のステップＳ１０８）。
【００８８】
他方、半導体スラブ導波路５１を伝搬する電磁界の伝搬モード（導波路モード）は、前述したように解析的に計算することができる。
【００８９】
そこで、これら二つの結果を用いることにより、位相速度誤差を厳密に評価した。
【００９０】
＜評価結果＞
図１０は、半導体スラブ導波路５１の導波路内波長あたりのセル数（横軸）と位相速度誤差（縦軸）の関係を示す図である。同図に示すグラフ７０では、半導体スラブ導波路５１内での導波路内波長当たりのセル数に対する位相速度誤差を百分率で表示した。
【００９１】
位相速度誤差は、２次モードが直線２１０（図で□の点を通過する実線で表示）、５次モードが直線５１０（図で○の点を通過する実線で表示）でそれぞれ示されている。これらの直線は、セル数を２倍にすると、位相速度誤差が１／４になっており、その傾きから２次精度を実現していることが分かる。
【００９２】
これに対し、曲線２２０（図で×の点を通過する破線で表示）および５２０（図で△の点を通過する破線で表示）は、従来法による電磁界シミュレーションに基づいて位相速度誤差を評価した結果を与えるものである。
【００９３】
ここで、従来法の概要を説明する。この方法は、使用する格子点の種類と位置ならびに電磁界の初期値、電磁界の時間変化を計算する差分式は本発明と同一である一方で、境界近傍の格子点への誘電率の割り当て方のみが異なっており、誘電率を割り当てるべき格子点の近傍の誘電率を算術平均し、その平均値を割り当てるものである。より具体的には、ｘ軸もしくはｙ軸に平行な長さΔの辺から成る正方形において、その正方形内部の誘電率を面積で重み付けした平均値を、正方形の対角線の交点に設けられる格子点に割り当てる。
【００９４】
この従来法は、２次元ＴＥ（Transverse Electric）モードについて適用可能であることが知られている（例えば、Allen Taflove, "Computational Electrodynamics, The Finite-Difference Time-Domain Method", p. 374, Artech House, 1995を参照）が、ここでは本実施例と比較するために２次元ＴＭモードに適用し、半導体スラブ導波路５１（を構成するYee格子上）における電磁波の伝搬シミュレーションを実行した。その際、導波路内波長当たりのセル数とδの値は、本実施例のシミュレーションと同じ値を採用した。
【００９５】
図１０からも明らかなように、従来法によって得られた曲線２２０および曲線５２０は、直線２１０や直線５１０のような直線ではないことに加え、その傾きの絶対値も直線２１０や５１０に比べて明らかに小さい。この結果から、本発明の電磁界シミュレーション方法が、従来法よりも高次の精度（２次精度）を実現していることが明らかになった。
【００９６】
以上説明した本発明の一実施例によれば、格子点２２および格子点２３のうち、境界近傍の格子点に対して式（７）〜（９）の関係式を満たす誘電率を付与し、時間差分式として式（１１）〜（１３）で定義される２次精度の差分式を用いることにより、誘電率が不連続に変化する二つの媒質の境界が、シミュレーションを行うYee格子軸に対して所定角傾斜している場合にも、２次元ＴＭモードに対する精度の高い電磁界シミュレーションを実現できることが分かる。
【００９７】
【発明の効果】
以上の説明からも明らかなように、本発明によれば、媒質の誘電率が空間的に不連続に変化する直線状の境界を有する２次元平面において、その境界が、２次元平面を格子状に離散化したときの互いに直交する格子軸のいずれに対しても平行でない場合、その離散化した格子状の平面を伝搬する電磁界成分の高精度なシミュレーションを可能にする電磁界シミュレーション方法、電磁界シミュレーションプログラム、およびプログラム記録媒体を提供することができる。
【図面の簡単な説明】
【図１】本発明の一実施形態に係る電磁界シミュレーション装置の機能ブロック図である。
【図２】本発明の一実施形態に係る電磁界シミュレーション方法で用いられるYee格子の概略構成を表す説明図である。
【図３】図２のYee格子の第２種および第３種格子点に割り当てられる誘電率の設定例（誘電率設定法Ｉ）を示す説明図である。
【図４】誘電率設定法Ｉに従って割り当てられたYee格子の境界付近の格子点の誘電率を例示する図である。
【図５】図２のYee格子の第２種および第３種格子点に割り当てられる誘電率の設定例（誘電率設定法II）を示す説明図である。
【図６】誘電率設定法Ｉに従って割り当てられたYee格子の境界付近の格子点の誘電率を例示する図である。
【図７】本発明の一実施形態に係る電磁界シミュレーション方法の概要を示すフローチャート図である。
【図８】本発明の一実施例として用いられるシミュレーション領域を示す説明図である。
【図９】図８のシミュレーション領域の半導体スラブ導波路とバックグラウンドの境界近傍の第２種および第３種格子点に付与される誘電率を示す図である。
【図１０】本発明の一実施形態に係る電磁界シミュレーション方法によって位相速度誤差を評価した結果を示す図である。
【符号の説明】
１ 電磁界シミュレーション装置
１１ 入力部
１２ 演算部
１３ 記憶部
１４ 主記憶部
１５ 補助記憶部
１６ 出力部
２０ Yee 格子
２１、２２、２３ 格子点
２４、５３、５４ 境界
５０ シミュレーション領域
５１ 半導体スラブ導波路
５２ バックグラウンド
２１０、５１０ 直線
２２０、５２０ 曲線[0001]
BACKGROUND OF THE INVENTION
The present invention relates to an electromagnetic field simulation method for analyzing an electromagnetic field using a computer, an electromagnetic field simulation program, and a program recording medium.
[0002]
[Prior art]
In recent years, “informatization of society” has made remarkable progress. Among them, electromagnetic waves as information carriers are used for information transmission in an increasingly sophisticated form. For example, optical waveguides, antennas, and the like have been required to be designed more precisely in accordance with their respective purposes. In such a situation, it is required to perform detailed simulations of changes in electromagnetic fields over time. .
[0003]
The electromagnetic field simulation method used for the above-mentioned purpose is the Finite-Difference Time-Domain method (FDTD, for example, Yee, KS, “Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media ”, IEEE Transactions on Antennas and Propagation, Vol.Ap-14, pp. 302-307, 1966.). In this simulation method, Maxwell's equation, which is a basic equation describing the time change of an electromagnetic field, is spatially differentiated using a two-dimensional lattice called a Yee lattice, and the time change of the electromagnetic field is calculated using a computer. .
[0004]
The time domain difference method approximates the partial partial differential of the electromagnetic field at the same time required to track the time change of the electromagnetic field using a discretized spatial difference formula, and based on the result of the spatial difference formula, The time change of the electromagnetic field is sequentially calculated using a time difference formula.
[0005]
This time domain difference method is classified according to the accuracy order of the difference equation used for the calculation. The accuracy order of the difference formula is such that the difference interval (time step in the case of time difference formula, cell width in the case of spatial difference formula) is τ, and when τ approaches 0, the error is τⁿIf it decreases in proportion to, it is called n-order accuracy. In the case of a spatial difference, this identification is applied in a situation where the permittivity and permeability are spatially constant, a situation where the spatial change of the permittivity and permeability is continuous, or the like.
[0006]
The most commonly used time-domain difference method is a time-domain difference method with time-secondary accuracy and space-secondary accuracy. Conventionally, various methods for achieving second-order accuracy are known when the spatial variation of the dielectric constant within the simulation region is not continuous, that is, when there is a boundary where the dielectric constant is discontinuous. For example, in Patent Document 1, when the boundary is parallel to one of the lattice axes in the discretized lattice plane, the dielectric constant for calculating the time change of the electric field component near the boundary is devised and assigned. Thus, a technique for achieving secondary accuracy is disclosed.
[0007]
[Patent Document 1]
JP 2001-243215 A
[0008]
[Problems to be solved by the invention]
However, in the time domain difference method of the above-described conventional time quadratic accuracy and spatial quadratic accuracy, there is a boundary between media having different dielectric constants in the simulation region, and these boundaries are orthogonal to each other in a discretized lattice plane. If it was not parallel to any of the lattice axes, it was not clear how to assign the dielectric constant near the boundary.
[0009]
The present invention has been made in view of such circumstances, and an object thereof is to provide a two-dimensional boundary in a two-dimensional plane having a linear boundary in which the dielectric constant of the medium changes spatially discontinuously. Electromagnetic field that enables highly accurate simulation of electromagnetic field components propagating through the discretized grid plane when the plane is discretized into a grid pattern and is not parallel to any of the orthogonal lattice axes The object is to provide a simulation method, an electromagnetic field simulation program, and a program recording medium.
[0010]
[Means for Solving the Problems]
  To achieve the above objective,Electromagnetic field simulation method according to the first aspect of the present inventionIs the dielectric constant ε in one plane₁ , Ε₂ Two mediums each having a boundary with a straight line as a boundary, and an electromagnetic field having a value of a magnetic field component parallel to the plane is zero in the plane, the discretized corresponding to the plane The virtual plane is made up of a plurality of first, second, and third types of grid points that have different electromagnetic fields assigned to each other. The second type lattice point is located at the midpoint between the two first type lattice points closest to the second type lattice point, and the second type lattice point is located among the midpoints of the two first type lattice points closest to the third type lattice point. The magnetic field component perpendicular to the plane is assigned to the first type lattice point while being arranged at a middle point different from the middle point, while the second type lattice point is assigned to the second type lattice point and the latest type lattice point. Perpendicular to the line connecting the first-type grid points in contact with each other, and In addition, an electric field component that is orthogonal to a line segment that connects the first type lattice point closest to the third type lattice point and that is parallel to the plane is assigned to the third type lattice point. The electromagnetic field simulation method for calculating the time variation of the electromagnetic field in the configured virtual plane by numerical calculation using a computer, wherein the first type lattice points and the third type lattice points are alternately arranged. One of the columns is the x axis, and one of the columns in which the first type lattice points and the second type lattice points are alternately arranged is the y axis, and the dielectric constant ε in the negative direction of the y axis₁ When the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using the coordinate system in which the medium having the position is located, the absolute value of the inclination of the straight line is 2, and the y of the straight line If the intercept is expressed as sΔ using a constant s that takes any real value in the range greater than -1.5 and less than 1.5 and the lattice spacing Δ, it is parallel to the straight line representing the boundary, and The dielectric constant of the third type lattice point arranged on the straight line whose y-intercept is −Δ is a, which is parallel to the straight line representing the boundary, and is arranged on the straight line whose y-intercept is −Δ / 2. The dielectric constant of the second type lattice point is b, parallel to the straight line representing the boundary, and the dielectric constant of the third type lattice point arranged on the straight line having a y-intercept of 0, and the straight line representing the boundary. And lattice points of the second kind arranged on a straight line with a y-intercept of Δ / 2. The electric permittivity is d, the dielectric constant of the third-type lattice point arranged on the straight line that is parallel to the straight line representing the boundary and the y intercept is Δ, and the a, b, c, d, e Dielectric constant ε in the plane corresponding to the imaginary plane among the second type or third type lattice points to which either is not assigned₁The dielectric constant of the lattice point corresponding to the medium of₁ , Among the second type or third type lattice points to which any of a, b, c, d, and e is not assigned, a dielectric constant ε in the plane corresponding to the virtual plane₂The dielectric constant of the lattice point corresponding to the medium of₂As said computer(a)Relational expression
## EQU11 ##

(Where the coefficient f_ij(I, j = 0,1,2) is
[Expression 12]

Coefficient g_ijk(I = 0, 1; j, k = 0, 1, 2) is
[Formula 13]

Coefficient h_ij(I, j = 0,1) is
[Expression 14]

Dielectric constant a, b, c, d, e, ε to satisfy₁ , ε₂ The setAnd store in memoryAnd steps to(b) Setting an initial value of time n representing discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory; (c) Read the dielectric constant and magnetic field component set at each lattice point from the memory, and time n + Calculating the electric field component of each grid point at 1/2 and storing it in memory; (d) Reading out the electric field component from the memory, calculating the magnetic field component at each lattice point at time n + 1, and storing it in the memory; (e) Incrementing time n by 1; (f) If time n is compared with a predetermined threshold and time n has not reached the predetermined threshold, (c) And when the time n has reached a predetermined threshold value, the step of reading out and outputting the electric field component and the magnetic field component from the memory;The main point is to execute.
[0011]
  Electromagnetic field simulation method according to the second aspect of the present inventionIs the dielectric constant ε in one plane₁ , Ε₂ Two mediums each having a boundary with a straight line as a boundary, and an electromagnetic field having a value of a magnetic field component parallel to the plane is zero in the plane, the discretized corresponding to the plane The virtual plane is made up of a plurality of first, second, and third types of grid points that have different electromagnetic fields assigned to each other. The second type lattice point is located at the midpoint between the two first type lattice points closest to the second type lattice point, and the second type lattice point is located among the midpoints of the two first type lattice points closest to the third type lattice point. The magnetic field component perpendicular to the plane is assigned to the first type lattice point while being arranged at a middle point different from the middle point, while the second type lattice point is assigned to the second type lattice point and the latest type lattice point. Perpendicular to the line connecting the first-type grid points in contact with each other, and In addition, an electric field component that is orthogonal to a line segment that connects the first type lattice point closest to the third type lattice point and that is parallel to the plane is assigned to the third type lattice point. The electromagnetic field simulation method for calculating the time variation of the electromagnetic field in the configured virtual plane by numerical calculation using a computer, wherein the first type lattice points and the third type lattice points are alternately arranged. One of the columns is the x axis, and one of the columns in which the first type lattice points and the second type lattice points are alternately arranged is the y axis, and the dielectric constant ε in the negative direction of the y axis₁ When the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using the coordinate system in which the medium having the position is located, the absolute value of the inclination of the straight line is 2, and the y of the straight line When the intercept is expressed as (1/2 + δ) Δ using a constant δ that takes an arbitrary real value in a range greater than −1 and less than 1, and the lattice spacing Δ, the intercept is parallel to the straight line representing the boundary; In addition, the dielectric constant of the third-type lattice points arranged on the straight line having the y-intercept of 0 is u, and is parallel to the straight line representing the boundary, and is arranged on the straight line having the y-intercept of Δ / 2. The dielectric constant of the second type lattice point is v, the dielectric constant of the third type lattice point arranged on the straight line that is parallel to the straight line representing the boundary and the y-intercept is Δ, and u, v, w Among the second-type or third-type grid points to which any of the above is not assigned, the virtual plane Permittivity within the corresponding said plane ε₁The dielectric constant of the lattice point corresponding to the medium of₁ Among the second or third type lattice points to which any of u, v, and w is not assigned, a dielectric constant ε in the plane corresponding to the virtual plane₂The dielectric constant of the lattice point corresponding to the medium of₂As said computer(a)Relational expression
[Expression 15]

To meetInvitationElectricity u, v, w, ε₁ , ε₂ The setAnd store in memoryAnd steps to(b) Setting an initial value of time n representing discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory; (c) Reading a dielectric constant and a magnetic field component set at each lattice point from the memory, calculating an electric field component at each lattice point at time n + 1/2, and storing the memory component in the memory; (d) Reading out the electric field component from the memory, calculating the magnetic field component at each lattice point at time n + 1, and storing it in the memory; (e) Incrementing time n by 1; (f) If time n is compared with a predetermined threshold and time n has not reached the predetermined threshold, (c) And when the time n has reached a predetermined threshold value, the step of reading out and outputting the electric field component and the magnetic field component from the memory;The main point is to execute.
[0012]
  In the above electromagnetic field simulation methodThe gist is to sequentially obtain an electromagnetic field corresponding to a change in time discretized by a time-domain difference method using a time difference equation having secondary accuracy and a space difference equation having secondary accuracy.
[0013]
  Electromagnetic field simulation program according to the third aspect of the present inventionIs the dielectric constant ε in one plane₁ , Ε₂ Two mediums each having a boundary with a straight line as a boundary, and an electromagnetic field having a value of a magnetic field component parallel to the plane is zero in the plane, the discretized corresponding to the plane The virtual plane is made up of a plurality of first, second, and third types of grid points that have different electromagnetic fields assigned to each other. The second type lattice point is located at the midpoint between the two first type lattice points closest to the second type lattice point, and the second type lattice point is located among the midpoints of the two first type lattice points closest to the third type lattice point. The magnetic field component perpendicular to the plane is assigned to the first type lattice point while being arranged at a middle point different from the middle point, while the second type lattice point is assigned to the second type lattice point and the latest type lattice point. Perpendicular to the line connecting the first-type grid points in contact with each other, and In addition, an electric field component that is orthogonal to a line segment that connects the first type lattice point closest to the third type lattice point and that is parallel to the plane is assigned to the third type lattice point. Any one of the columns in which the first-type grid points and the third-type grid points are alternately arranged in order to obtain the time variation of the electromagnetic field in the configured virtual plane by numerical calculation using a computer. Is the x axis, and one of the columns in which the first-type lattice points and the second-type lattice points are alternately arranged is the y-axis, and the dielectric constant ε in the negative direction of the y-axis₁ When the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using the coordinate system in which the medium having the position is located, the absolute value of the inclination of the straight line is 2, and the y of the straight line If the intercept is expressed as sΔ using a constant s that takes any real value in the range greater than -1.5 and less than 1.5 and the lattice spacing Δ, it is parallel to the straight line representing the boundary, and The dielectric constant of the third type lattice point arranged on the straight line whose y-intercept is −Δ is a, which is parallel to the straight line representing the boundary, and is arranged on the straight line whose y-intercept is −Δ / 2. The dielectric constant of the second type lattice point is b, parallel to the straight line representing the boundary, and the dielectric constant of the third type lattice point arranged on the straight line having a y-intercept of 0, and the straight line representing the boundary. And lattice points of the second kind arranged on a straight line with a y-intercept of Δ / 2. The electric permittivity is d, the dielectric constant of the third-type lattice point arranged on the straight line that is parallel to the straight line representing the boundary and the y intercept is Δ, and the a, b, c, d, e Dielectric constant ε in the plane corresponding to the imaginary plane among the second type or third type lattice points to which either is not assigned₁The dielectric constant of the lattice point corresponding to the medium of₁ , Among the second type or third type lattice points to which any of a, b, c, d, and e is not assigned, a dielectric constant ε in the plane corresponding to the virtual plane₂The dielectric constant of the lattice point corresponding to the medium of₂To the computer,(a)Relational expression
[Expression 16]

(Where the coefficient f_ij(I, j = 0,1,2) is
[Expression 17]

Coefficient g_ijk(I = 0, 1; j, k = 0, 1, 2) is
[Expression 18]

Coefficient h_ij(I, j = 0,1) is
[Equation 19]

Dielectric constant a, b, c, d, e, ε to satisfy₁ , ε₂ The setAnd store in memoryAnd steps to(b) Setting an initial value of time n representing discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory; (c) Read the dielectric constant and magnetic field component set at each lattice point from the memory, and time n + Calculating the electric field component of each grid point at 1/2 and storing it in memory; (d) Reading out the electric field component from the memory, calculating the magnetic field component at each lattice point at time n + 1, and storing it in the memory; (e) Incrementing time n by 1; (f) If time n is compared with a predetermined threshold and time n has not reached the predetermined threshold, (c) And when the time n has reached a predetermined threshold value, the step of reading out and outputting the electric field component and the magnetic field component from the memory;The gist is to execute.
[0014]
  Electromagnetic field simulation program according to the fourth aspect of the present inventionIs the dielectric constant ε in one plane₁ , Ε₂ Two mediums each having a boundary with a straight line as a boundary, and an electromagnetic field having a value of a magnetic field component parallel to the plane is zero in the plane, the discretized corresponding to the plane The virtual plane is made up of a plurality of first, second, and third types of grid points that have different electromagnetic fields assigned to each other. A midpoint where the second type grid point is located at the midpoint between the two first type grid points closest to each other, and the second type grid point is closest to the third type grid point; Are arranged at different midpoints, and a magnetic field component perpendicular to the plane is assigned to the first type lattice point, while the second type lattice point is closest to the second type lattice point. Perpendicular to the line connecting the seed grid points and flat on the plane Further, an electric field component that is orthogonal to a line segment connecting the first type lattice point closest to the third type lattice point and parallel to the plane is assigned to the third type lattice point. Then, in order to obtain the time change of the electromagnetic field in the configured virtual plane by numerical calculation using a computer, one of the columns in which the first type lattice points and the third type lattice points are alternately arranged is represented by x. On the other hand, one of the columns in which the first-type lattice points and the second-type lattice points are alternately arranged is defined as a y-axis, and a dielectric constant ε in the negative direction of the y-axis₁ When the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using the coordinate system in which the medium having the position is located, the absolute value of the inclination of the straight line is 2, and the y of the straight line When the intercept is expressed as (1/2 + δ) Δ using a constant δ that takes an arbitrary real value in a range greater than −1 and less than 1, and the lattice spacing Δ, the intercept is parallel to the straight line representing the boundary; In addition, the dielectric constant of the third-type lattice points arranged on the straight line having the y-intercept of 0 is u, and is parallel to the straight line representing the boundary, and is arranged on the straight line having the y-intercept of Δ / 2. The dielectric constant of the second type lattice point is v, the dielectric constant of the third type lattice point arranged on the straight line that is parallel to the straight line representing the boundary and the y-intercept is Δ, and u, v, w Among the second-type or third-type grid points to which any of the above is not assigned, the virtual plane Permittivity within the corresponding said plane ε₁The dielectric constant of the lattice point corresponding to the medium of₁ Among the second or third type lattice points to which any of u, v, and w is not assigned, the dielectric constant ε in the plane corresponding to the virtual plane₂The dielectric constant of the lattice point corresponding to the medium of₂To the computer,(a)Relational expression
[Expression 20]

To meetInvitationElectricity u, v, w, ε₁ , ε₂ The setAnd store in memoryAnd (b) setting an initial value of time n representing the discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory; (c) reading out the dielectric constant and magnetic field component set at each lattice point from the memory, calculating the electric field component at each lattice point at time n + 1/2, and storing it in the memory; and (d) removing the electric field component from the memory. Reading, calculating the magnetic field component of each lattice point at time n + 1 and storing it in the memory; (e) increasing the time n by 1; (f) comparing the time n with a predetermined threshold; If the predetermined threshold value has not been reached, the process returns to step (c). If the time n has reached the predetermined threshold value, the step of reading out and outputting the electric field component and magnetic field component from the memory is executed. Is the gist.
[0015]
  Electromagnetic field simulation programThe gist is to sequentially obtain an electromagnetic field corresponding to a change in time discretized by a time-domain difference method using a time difference equation having secondary accuracy and a space difference equation having secondary accuracy.
[0016]
  Program recording medium according to the fifth aspect of the present inventionIsthe aboveElectromagnetic field simulation programEitherThis is the summary.
[0017]
The program recording medium in the present invention means a computer-readable recording medium such as a hard disk, a flexible disk, a CD-ROM (Compact Disc Read Only Memory), a DVD (Digital Versatile Disk), a magneto-optical disk, and a PC card.
[0018]
DETAILED DESCRIPTION OF THE INVENTION
Next, embodiments of the present invention will be described with reference to the drawings.
[0019]
FIG. 1 is a functional block diagram showing a configuration of an electromagnetic field simulation apparatus that is an electronic apparatus that executes an electromagnetic field simulation method according to an embodiment of the present invention. An electromagnetic field simulation apparatus 1 shown in FIG. 1 is an electromagnetic field simulation using an input unit 11 including an input device such as a keyboard and a mouse for inputting basic data for simulation and various parameters, and input basic data and various parameters. In order to perform the above, a calculation unit 12 composed of a central processing unit (CPU: Central Processing Unit) having a calculation function and a control function, input information from the input unit 11, a calculation result from the calculation unit 12, and the like are stored. The storage unit 13 includes at least an output unit 16 including a display device for outputting information stored in the storage unit 13 as a simulation result.
[0020]
The storage unit 13 includes a main storage unit 14 including a RAM (Random Access Memory), a hard disk drive, a flexible disk drive, a CD-ROM (Compact Disc Read Only Memory) drive, a DVD (Digital Versatile Disk) drive, optical And an auxiliary storage unit 15 including a magnetic disk drive, a PC card drive, and the like, and a memory area necessary for storing calculation results to be described later is secured.
[0021]
Various processes of the present embodiment to be described below are not only performed by one electronic device, but the entire system constructed from two or more electronic devices by appropriately dividing the execution of each step. This includes the case of executing with Therefore, in this embodiment, a device (computer) including one or a plurality of electronic devices is referred to as an “electromagnetic field simulation device”. In this sense, this embodiment includes, for example, a case where an OS (Operating System) operating on the electromagnetic field simulation apparatus 1 performs actual processing in accordance with an instruction of a program for performing electromagnetic field simulation processing. Needless to say. This point is common to all the embodiments of the present invention.
[0022]
FIG. 2 shows two regions (dielectric constants are expressed as ε) when the electromagnetic field simulation of this embodiment is executed.₁, Ε₂Is a schematic configuration of a two-dimensional Yee lattice (hereinafter simply referred to as a Yee lattice) as a lattice-like virtual plane obtained by discretizing the planes of continuous media adjacent to each other at a linear boundary. FIG. The Yee lattice 20 shown in the figure has three different types of lattice points 21 (first type lattice points), lattice points 22 (second type lattice points), and lattice points 23 (third type lattice points) alternately. Are two-dimensionally arranged, and each constitutes a two-dimensional lattice having a square lattice as a unit lattice (a lattice constant which is a lattice interval between lattice points of the same type is Δ, where Δ is a positive constant). To form a discretized plane (virtual plane). More specifically, using the lattice points, the same kind of lattice points each form a square shape with a lattice interval Δ, and the lattice point 22 is positioned at the midpoint between the two closest lattice points 21 and the lattice point 23. Is arranged so that it is located at a midpoint different from the midpoint where the grid point 22 is located among the midpoints of the two grid points 21 closest to each other, and a magnetic field component perpendicular to the plane is assigned to the grid point 21 Constitutes the discretized plane. Accordingly, in this Yee lattice 20, the distance between the closest dissimilar lattice points is Δ / 2.
[0023]
In FIG. 2, the direction (horizontal direction in the drawing) parallel to the column in which the lattice points 21 and the lattice points 23 are alternately arranged is defined as the x-axis direction, and one of the columns is defined as the x-axis. On the other hand, the direction (vertical direction in the figure) parallel to the column in which the lattice points 21 and 22 are alternately arranged is taken as the y-axis, and one of these rows is taken as the y-axis. One of the lattice points 21 is arranged at the origin of such coordinate axes. In FIG. 2, the coordinate system satisfying the above setting is Σ₀Although the case of taking is shown in the figure, it goes without saying that the method of taking this coordinate system is arbitrary.
[0024]
Coordinate system Σ shown in FIG.₀By using arbitrary integers i and j, the coordinates of the lattice point 21 are (iΔ, jΔ), the coordinates of the lattice point 22 are (iΔ, (j + 1/2) Δ), and the coordinates of the lattice point 23 are ((I + 1/2) Δ, jΔ), respectively. In the following description, in order to simplify the description, when it can be considered that the meaning can be sufficiently understood without specifying Δ, the x and y coordinates of each lattice point are respectively divided by Δ to obtain a simple integer i. , J, for example, the coordinates of the lattice points 21, 22, 23 are represented as (i, j), (i, j + 1/2, j), and (i + 1/2, j), respectively.
[0025]
In the present embodiment, a magnetic field component H in the direction perpendicular to the plane of FIG._z(I, j) is assigned, and the grid point 22 has an electric field component E in the x direction._x(I, j + 1/2) at the grid point 23, the electric field component E in the y direction_y(I + 1/2, j) are assigned respectively. This makes it possible to model a two-dimensional TM (Transverse Magnetic) mode electromagnetic field in which only the magnetic field has vibration in a direction perpendicular to the Yee lattice plane. Incidentally, in this two-dimensional TM mode, the electric field component E in the z direction_z, Magnetic field component H in the x direction_x, And magnetic field component H in the y direction_yIs always 0 (constant) regardless of time.
[0026]
As shown in FIG. 2, the boundary 24 between the two media points in a direction that forms a predetermined inclination angle with respect to the x-axis direction and the y-axis direction, which are the arrangement directions of the lattices. More specifically, when the boundary 24 is expressed by using y as a linear function of x, the magnitude of the slope of the straight line is 2. Therefore, for the convenience of the following description, the slope of the straight line forming this boundary is set to −2, and Δ and s (−1.5 <s <1.5) are used.
y = -2x + sΔ (1)
However, in the electromagnetic field simulation method of the present embodiment, the slope may be +2, or the y-intercept may take a value other than the above. These depend on the arbitraryness of the coordinate system. However, it goes without saying that the generality of the electromagnetic field simulation method according to the present embodiment is not lost even if the above limitation is added.
[0027]
The dielectric constants of the lattice points 22 and 23 (both lattice points having non-zero electric field components) arranged in the vicinity of the boundary 24 represented by the equation (1) are set as follows:
(I-1) Let a be the dielectric constant of the lattice point 23 arranged on the straight line y = −2x−Δ.
(I-2) Let b be the dielectric constant of the lattice point 22 arranged on the straight line y = −2x−Δ / 2.
(I-3) Let c be the dielectric constant of the lattice point 23 arranged on the straight line y = −2x.
(I-4) Let d be the dielectric constant of the lattice point 22 arranged on the straight line y = −2x + Δ / 2.
(I-5) Let e be the dielectric constant of the lattice point 23 arranged on the straight line y = −2x + Δ.
[0028]
The dielectric constants of the lattice points 22 and 23 other than these are set as follows:
(I-6) region of y <−2x + sΔ (this region is defined as region D₁The dielectric constant of lattice points 22 and 23 other than the above arranged in₁And
(I-7) region of y> −2x + sΔ (this region is defined as region D₂The dielectric constant of lattice points 22 and 23 other than the above arranged in₂And
[0029]
It is assumed that the magnetic permeability is constant regardless of the medium.
[0030]
FIG. 3 is a table in which the above dielectric constant assignments are summarized in one table 30. FIG. In the table 30 shown in the drawing, for example, a column described as “y = −2x−Δ” means a dielectric constant assigned to each lattice point arranged on this straight line. That is, from the above (I-1), a dielectric constant a is assigned to the lattice point 23. On the other hand, since there is no lattice point 22 on the straight line y = −2x−Δ, a horizontal line (−) is drawn in the column in which the dielectric constant of the lattice point 22 is described in Table 30. The same applies to the descriptions in the other columns.
[0031]
FIG. 4 is an explanatory diagram conceptually showing a dielectric constant assigned to each lattice point near the boundary 24. Note that the lattice points 22 and 23 that are described by the same symbols as in FIG. 2 and for which the assigned dielectric constant is not specified are appropriately set in accordance with the above (I-6) and (I-7).₁Or ε₂Is assigned.
[0032]
Region D₁Area D from the side₂When a plane wave is incident on the side, the reflectance R and transmittance T at the boundary based on the magnetic field amplitude can be expressed as a series of the lattice constant Δ as follows.
[Expression 21]

Where θ_iIs the incident angle, θ_tRepresents the reflection angle. Boundary 24 is dielectric constant ε₁Region and dielectric constant ε₂, The first term on the right side of Equations (2) and (3) is a true value.
[0033]
In the present embodiment, the coefficient R of the first-order error term₁And T₁Is incident angle θ_iA second-order accuracy simulation is performed with zero regardless of the value of. Therefore, in equations (2) and (3), the first-order coefficient of Δ is the incident angle θ_iWhen the condition that becomes zero regardless of the value of is obtained, the following equation is obtained.
[Expression 22]

Here, the coefficient f in equation (4)_ij(I, j = 0,1,2), coefficient g in equation (5)_ijk(I = 0, 1; j, k = 0, 1, 2), and coefficient h in equation (6)_ij(I, j = 0,1) is expressed as follows using the dielectric constants a, b, c, d, e set based on (I-1) to (I-5) above. The
[0034]
[Expression 23]

[0035]
[Expression 24]

[0036]
[Expression 25]

[0037]
By the way, region D₂To region D₁Even if a conditional expression that gives secondary accuracy is obtained when a plane wave is incident on the surface, the same expressions as Expressions (4) to (6) are obtained.
[0038]
The dielectric constant setting method described above is referred to as “dielectric constant setting method I”.
[0039]
By the way, in this dielectric constant setting method I,
[Equation 26]

Then, the following equations (7) to (9) are obtained by solving these for u, v, and w.
[Expression 27]

[0040]
When these equations (7) to (9) are used, a straight line representing the boundary 24 is
[Expression 28]

And assign the dielectric constants of lattice points 22 and 23 as follows:
First, c, d, and e in (I-3) to (I-5) are u, v, and w, respectively.
(II-1) Let u be the dielectric constant of the lattice point 23 arranged on the straight line y = −2x.
(II-2) Let v be the dielectric constant of the lattice point 22 arranged on the straight line y = −2x + Δ / 2.
(II-3) The dielectric constant of the lattice point 23 arranged on the straight line y = −2x + Δ is w.
further,
(II-4) Among lattice points 22 and 23 not included in the above three conditions, a region of y <−2x + (1/2 + δ) Δ (this region is defined as region D₁The dielectric constant of the lattice points 22 and 23 arranged in₁And
(II-5) The dielectric constant of lattice points 22 and 23 not included in the above four conditions is expressed as ε₂And This dielectric constant ε₂Is assigned to an area of y> −2x + (1/2 + δ) Δ (this area is designated as area D).₂Included).
[0041]
The dielectric constant setting method based on the conditions (II-1) to (II-5) and the expressions (7) to (9) will be referred to as “dielectric constant setting method II”. It is clear from the above description that the dielectric constant setting method II is a special case of the dielectric constant setting method I.
[0042]
FIG. 5 is a diagram in which the above-described conditions are collectively described in one table 40. The way of viewing the table 40 is the same as that in FIG.
[0043]
Further, the explanatory diagram of FIG. 6 conceptually shows the dielectric constant assigned to each lattice point near the boundary 24. Also in this figure, lattice points 22 and 23 which are described with the same symbols as in FIG. 2 and whose assigned dielectric constant is not specified are appropriately set in accordance with the above (II-4) and (II-5).₁Or ε₂Is assigned.
[0044]
As described above, each component of the electromagnetic field is assigned to a spatially different position, but also in a temporal sense, in a three-dimensional space-time formed by adding a time dimension to the plane of the Yee lattice 20. Assigned to different positions. For example, the time t when the magnetic field is discretized¹, t², ..., t^n-1, tⁿ, t^{n + 1},... (N is a positive integer), the electric field is t^1/2, t^{1 + 1/2}, t^{2 + 1/2}, ..., t^{n-1 / 2}, t^{n + 1/2}, t^{n + 3/2}...,. Setting the time in this way makes it possible to calculate the magnetic field and the electric field alternately. Of course, it is possible to define the times assigned to the electric field and the magnetic field in reverse.
[0045]
Next, an electromagnetic field simulation method based on the above settings will be described in detail.
[0046]
In the time domain difference method of time secondary accuracy and spatial secondary accuracy applied in the present embodiment, the following is obtained in a time-series sense from the electromagnetic field components stored in the memory in the storage unit 13 of the electromagnetic field simulation apparatus 1. The electromagnetic fields at the times are sequentially obtained by the following calculation.
[0047]
[Expression 29]

Where Δt is the time step interval, E_x ^{n + 1/2}(I, j + 1/2) is the time t of lattice point 22 (i, j + 1/2)^{n + 1/2}X direction component of the electric field at E_y ^{n + 1/2}(I + 1 / 2.j) is the time t of the grid point 23^{n + 1/2}Y-direction component of the electric field at H,_z ^{n + 1}(I, j) is the time t of the lattice point 21^{n + 1}Is the z-direction component of the magnetic field at. Ε (i, j + 1/2) is a dielectric constant given to the lattice point 22 at coordinates (i, j + 1/2), and ε (i + 1/2, j) is a lattice point at coordinates (i + 1/2, j). The dielectric constant given to 23. Further, μ is the magnetic permeability, and takes a constant value regardless of the medium as described above.
[0048]
In addition, [] on the right side of each expression means a spatial quadratic difference expression that approximates the partial differentiation described in [], and more specifically is defined by the following difference expression. R:
[30]

[0049]
When each electromagnetic field component can be partially differentiated at least twice as in the case where the dielectric constant of the medium is uniform, the spatial coordinates of each electromagnetic field component indicate that these differential equations give the spatial secondary accuracy. It can be shown using the Taylor expansion for (x, y). On the other hand, in the vicinity of the boundary 24, each electromagnetic field component is generally not always capable of partial differentiation twice with respect to a direction perpendicular to the boundary 24 (x direction). In particular, since the electric field component in the direction perpendicular to the boundary 24 is discontinuous, it has not been obvious whether high-precision calculations have been possible so far.
[0050]
On the other hand, in this embodiment, the boundary 24 has a coordinate system Σ.₀In FIG. 2, the dielectric constant assigned to the lattice points 22 and 23 in the vicinity of the boundary 24 satisfies the conditional expressions (4) to (6) or the expressions (7) to (9). By setting and using a normal secondary accuracy spatial difference formula (see formulas (14) to (17)), the reflectance and transmittance at the boundary 24 regardless of the incident angle of the electromagnetic wave incident on the boundary 24 It is possible to realize a two-dimensional TM mode electromagnetic field simulation that can be approximated with quadratic accuracy.
[0051]
FIG. 7 is a flowchart showing a specific processing procedure when using a time domain difference method with time secondary accuracy and space secondary accuracy as the electromagnetic field simulation method according to the present embodiment. In performing the processing shown in the figure, basic data for simulation and various parameters such as setting of the dielectric constant to the lattice points of the simulation region are input from the input unit 11 in advance and stored in the storage unit 13. It is assumed that the specific setting of the dielectric constant is made by the calculation unit 12 based on the input content.
[0052]
In the time domain difference method shown in FIG. 7, first, initial values of an electric field and a magnetic field are set (step S101). Here, the time of the initial magnetic field is t = −1 / 2, and the time of the initial electric field is t = 0.
[0053]
Next, the initial time (time step) is set to n = 0 (step S102). Needless to say, the symbol “: =” in FIG. 7 means “n is defined as 0”. Subsequent time steps, as described above, the electric field is t^1/2, t^{1 + 1/2}, t^{2 + 1/2}, ..., t^{n-1 / 2}, t^{n + 1/2}, t^{n + 3/2}, ..., the magnetic field is t¹, t², ..., t^n-1, tⁿ, t^{n + 1},… And time passes.
[0054]
Thereafter, the electromagnetic field components are sequentially calculated in the calculation unit 12 as time passes. Therefore, the following steps will be described with respect to an arbitrary time step represented by general n (a positive integer).
[0055]
First, using the equations (11) and (12), each component E of the electric field at time step t = n + 1/2_x ^{n + 1/2}(I, j + 1/2), E_y ^{n + 1/2}(I + 1/2, j) is calculated, and the calculation result is stored in the memory in the storage unit 13 in association with a set of lattice points and time steps (corresponding to coordinates in a three-dimensional space-time) (step S103).
[0056]
Using Equation (13), the magnetic field component H at time step t = n + 1_z ^{n + 1}(I, j) is calculated, and the calculation result is stored in the memory in the same manner as in step S103 (step S104).
[0057]
The calculation result of the electromagnetic field component is sequentially stored and stored in the memory in the storage unit 13, but the electric field or magnetic field at the previous time (time step) is changed to the electric field or magnetic field at the next time (time step). Even if the magnetic field is replaced, there is no problem in the subsequent calculation of the electromagnetic field. This is because the calculation of the electric field and the magnetic field is alternately executed at each time step, as is apparent from the equations (11) to (13). By sequentially updating the electromagnetic field components in this way, it is possible to reduce the memory area used for the calculation and reduce the load.
[0058]
Subsequently, the electric field components E obtained in steps S103 and S104, respectively._x ^{n + 1/2}(I, j + 1/2), E_y ^{n + 1/2}(I + 1/2, j) and magnetic field component H_z ^{n + 1}From (i, j), data related to them is extracted and stored in the memory (step S105). Incidentally, the extracted data regarding the electromagnetic field component can also be used when calculating data such as the intensity and phase of each mode of the electromagnetic field in a step described later. In this sense, the memory area for storing these data is preferably a different area from the memory area for storing the electromagnetic field component itself. Depending on the contents of the simulation to be executed, this step S105 is not essential, and after step S104, the process may proceed to step S106 described below.
[0059]
Thereafter, the value of n is incremented by 1 to advance the time step (step S106). In step S107, n is a threshold value n set in advance._maxNot reached (n <n_max), The process returns to step S103 to repeat the processes of steps S103 to S106 described above. Needless to say, the determination in step S107 is performed by the calculation unit 12.
[0060]
On the other hand, n is n in step S107._maxIf the electromagnetic field distribution is determined, the electromagnetic field distribution is output, the electromagnetic field amplitude, phase, reflectance, and the like are calculated, and the calculation result is output from the output unit 16 (step S108).
[0061]
The simulation results obtained by the above processing, the basic data and various parameters input for the simulation can be read from the storage unit 13 and output from the output unit 16 as necessary.
[0062]
As described above, according to an embodiment of the present invention, when the boundary is represented by a straight line having a certain inclination using a predetermined coordinate system, the conditional expression (4) is applied to the lattice points near the boundary. The incident angle of the electromagnetic wave incident on the boundary by setting the dielectric constant satisfying the expression (6) or the expressions (7) to (9) and using the time domain difference method with the second-order temporal accuracy and the second-order spatial accuracy. Regardless of this, it is possible to realize a two-dimensional TM mode electromagnetic field simulation capable of approximating the reflectance and transmittance at the boundary with second-order accuracy.
[0063]
In addition, by mounting a computer-readable program recording medium in which an electromagnetic field simulation program for realizing the electromagnetic field simulation process of the present embodiment is recorded on a computer, and reading the program stored in the program recording medium, A computer may execute the processing described above. Here, as the “computer-readable” program recording medium, a hard disk, flexible disk, CD-ROM, DVD, magneto-optical disk, PC card, or the like can be used. By providing such a program recording medium, the electromagnetic field simulation program of this embodiment can be widely distributed.
[0064]
(Other embodiments)
By the way, the present invention should not be understood as having a specific effect only in the above-described embodiment.
[0065]
For example, when setting the dielectric constant in the vicinity of the boundary, the dielectric constant in the original continuous medium (ε₁, Ε₂) Can be assigned a different dielectric constant. Specifically, in the vicinity of the boundary, m (m is an arbitrary integer greater than or equal to 2) straight lines in the order of short distance from the straight line forming the boundary among the parallel lines forming the boundary and different y-intercepts. A dielectric constant can be set for each lattice line passing therethrough (corresponding to the case where the dielectric constant setting method I is m = 5 and the dielectric constant setting method II is m = 3).
[0066]
In addition, the inclination of the straight line that forms the boundary with respect to the axial direction of the Yee lattice is not limited to the above case, and the electromagnetic field simulation is similarly performed for a boundary having another inclination with respect to the axial direction. Is possible.
[0067]
Furthermore, it is of course possible to perform electromagnetic field simulation using a Yee lattice that is not a square lattice.
[0068]
As described above, in the present invention, as in the above-described one embodiment, various implementations capable of realizing a highly accurate electromagnetic field simulation by using the time domain difference method with time secondary accuracy and space secondary accuracy. It can include forms and the like.
[0069]
【Example】
Examples of electromagnetic field simulation using the electromagnetic field simulation method according to an embodiment of the present invention will be described below.
[0070]
<Simulation area>
FIG. 8 is an explanatory diagram showing an outline of a simulation region applied in this embodiment. Note that the two-dimensional plane of the continuous medium to be discretized has a dielectric constant of ε₂= 10.24 semiconductor slab waveguide 51 and permittivity is ε₁A background 52 which is a region of = 1.0 is arranged as shown in FIG.
[0071]
The simulation region 50 shown in FIG. 8 has a rectangular shape, and at the boundary with the background 52 of the semiconductor slab waveguide 51, a straight line forming the boundary is expressed by the equation (10) in an appropriate local coordinate system provided near the boundary. ). That is, at the boundary 53, a local coordinate system Σ with the vertical downward direction in the figure as the x-axis, the horizontal direction leftward as the y-axis, and the direction perpendicular to the drawing and going to the back side of the figure as the z-axis.₁On the other hand, at the boundary 54, a local coordinate system Σ in which the vertical upward direction in the figure is the x-axis, the horizontal direction rightward is the y-axis, and the direction perpendicular to the drawing and going to the back side of the figure is the z-axis.₂By taking, each boundary 53 and 54 is represented as a straight line in equation (10).
[0072]
A Yee lattice 20 similar to that in FIG. 2 is formed in the simulation region 50, and the local coordinate system Σ₁Or Σ₂When the coordinates of each point are determined based on the magnetic field component H_zGrid point 21 (i, j) to which (i, j) is assigned, electric field component E_xLattice point 22 (i, j + 1/2) to which (i, j + 1/2) is assigned, electric field component E_yLattice points 23 (i + 1/2, j) to which (i + 1/2, j) are assigned are alternately arranged, and each lattice point forms a square lattice having the same lattice constant Δ.
[0073]
<Basic data for simulation and various parameter settings>
The basic data for simulation and various parameters will be described. In this embodiment, it is assumed that the dielectric constant assigned to each lattice point is set based on the dielectric constant setting method II.
[0074]
Local coordinate system of boundary 53₁And the length of the diagonal line of the simulation region 50 is an integral multiple of the eigenmode in which the propagation simulation is performed, that is, (2π / [propagation constant of the mode in which the propagation simulation is performed]). An integer multiple of. Further, a periodic boundary condition is imposed around the simulation region 50, and the width of the semiconductor slab waveguide 51 is set to 1.5 μm (1.5 × 10 5^-6m).
[0075]
On the other hand, the local coordinate system Σ of the boundary 54₂The value of δ in equation (10) of the straight line varies slightly depending on the value of the lattice constant Δ, but the value of δ can be determined so that the rate of change is at most about 2%. An example of a specific value of δ will be described later.
[0076]
In contrast to the above setting, the local coordinate system Σ of the boundary 54₂First, the value of δ in the straight line equation (10) is set to 0, and then the local coordinate system Σ of the boundary 53₁Of course, the value of δ in the straight line equation (10) may be determined.
[0077]
Next, the setting of electromagnetic waves for performing propagation simulation will be described.
[0078]
First, the free space wavelength of the propagating electromagnetic wave was set to 1.55 μm. In this case, the magnetic field component H_z, And electric field component E in two directions_x, E_yThe propagation mode of the electromagnetic wave represented by can hold 6 TM modes (0th to 5th). In each mode, the incident angle of the wave propagating through the semiconductor slab waveguide 51 to the waveguide boundary is different. In this example, propagation simulation was performed for the second-order mode with an incident angle of 61.4 degrees and the fifth-order mode with an incident angle of 22.1 degrees.
[0079]
Regarding the initial value of the electromagnetic field component assigned to each grid point, the result obtained analytically was assigned to each grid point as an initial value.
[0080]
FIG. 9 is a diagram showing a setting example of the dielectric constant given to the lattice point 22 or 23 in the vicinity of the boundary 53 or 54 between the semiconductor slab waveguide 51 and the background 52. In the table 60 shown in the figure, u = 0.82636, v = 1.9393, and w = 16.061 at the boundary where δ = 0.0 (at least one of the boundaries 53 and 54 takes this value).
[0081]
On the other hand, at the boundary where δ ≠ 0, when the number of cells per wavelength in the semiconductor slab waveguide 51 (hereinafter, the number of cells per wavelength in the waveguide) is 23.195 and δ = 0.4, u = 0.817933, v = 0.96931 , w = 7.1722. In Table 60, a value of 0.39947 (= 0.4-0.0053) in parentheses is described at δ = 0.4, but this indicates that when δ ≠ 0, this is a lattice constant Δ (per waveguide wavelength). This is because the number of cells is slightly changed according to the value of (which is an index representing the degree of fineness of the lattice). However, since this change is at most about 2% as described above, this does not particularly affect the electromagnetic field simulation result of this embodiment.
[0082]
Further, when the number of cells per wavelength in the waveguide constituting the semiconductor slab waveguide 51 is 17.503 and δ = −0.2, u = 0.96970, v = 2.74717, and w = 23.4416. In this case, the value −0.19828 (= −0.2 + 0.00172) described in parentheses in Table 60 indicates an accurate value used for the simulation when the number of cells per wavelength in the waveguide is 17.503.
[0083]
Regarding the dielectric constants of lattice points other than u, v, and w, the dielectric constant ε is obtained by using the background 52 as a vacuum.₁1.0, while permittivity ε of the semiconductor slab waveguide 51 other than the above₂Was 10.24.
[0084]
In the simulation of the present embodiment, the CFL number (Courant-Friedrichs-Levy number) that gives an index of how much the slice width of time is to the slice width of space when discretization by the Yee lattice is performed. ) Was 0.6.
[0085]
<Simulation contents>
Based on the basic data for simulation described above and the settings of various parameters, simulation was performed on Yee lattices having various lattice constants Δ for the second-order mode and the fifth-order mode of the two-dimensional TM mode.
[0086]
Since the accuracy of the reflectance R and the transmittance T can be evaluated by the accuracy of the phase velocity of the electromagnetic wave, the error of this phase velocity was examined in this example.
[0087]
The phase velocity as the simulation result is obtained by extracting phase information from the electromagnetic field distribution of the propagation simulation (step S108 in FIG. 7).
[0088]
On the other hand, the propagation mode (waveguide mode) of the electromagnetic field propagating through the semiconductor slab waveguide 51 can be analytically calculated as described above.
[0089]
Therefore, the phase velocity error was strictly evaluated by using these two results.
[0090]
<Evaluation results>
FIG. 10 is a diagram showing the relationship between the number of cells per waveguide wavelength (horizontal axis) and the phase velocity error (vertical axis) of the semiconductor slab waveguide 51. In the graph 70 shown in the figure, the phase velocity error with respect to the number of cells per wavelength in the waveguide in the semiconductor slab waveguide 51 is displayed as a percentage.
[0091]
The phase velocity error is indicated by a straight line 210 (indicated by a solid line passing through the □ point in the figure) for the second order mode and a straight line 510 (indicated by a solid line passing through the point in the figure) for the fifth order mode. . These lines show that when the number of cells is doubled, the phase velocity error becomes ¼, and the secondary accuracy is realized from the inclination.
[0092]
In contrast, curves 220 (indicated by a broken line passing through a point x in the figure) and 520 (indicated by a broken line passing a point of Δ in the figure) evaluate the phase velocity error based on an electromagnetic field simulation by a conventional method. Results.
[0093]
Here, an outline of the conventional method will be described. In this method, the type and position of the grid points to be used, the initial value of the electromagnetic field, and the differential equation for calculating the time variation of the electromagnetic field are the same as in the present invention, while the dielectric constant is assigned to the grid points near the boundary. The only difference is that the dielectric constant in the vicinity of the lattice point to which the dielectric constant should be assigned is arithmetically averaged and the average value is assigned. More specifically, in a square composed of sides of length Δ parallel to the x-axis or y-axis, an average value obtained by weighting the dielectric constant inside the square by the area is used as a lattice point provided at the intersection of the diagonal lines of the square. assign.
[0094]
This conventional method is known to be applicable to a two-dimensional TE (Transverse Electric) mode (for example, Allen Taflove, “Computational Electrodynamics, The Finite-Difference Time-Domain Method”, p. 374, Artech House). , 1995) was applied to the two-dimensional TM mode for comparison with the present embodiment, and an electromagnetic wave propagation simulation in the semiconductor slab waveguide 51 (on the Yee lattice constituting it) was executed. At that time, the same number of cells per wavelength in the waveguide and the value of δ were adopted as in the simulation of this example.
[0095]
As is clear from FIG. 10, the curve 220 and the curve 520 obtained by the conventional method are not straight lines such as the straight line 210 and the straight line 510, and the absolute value of the slope is also larger than that of the straight lines 210 and 510. Obviously small. From this result, it has been clarified that the electromagnetic field simulation method of the present invention achieves higher accuracy (secondary accuracy) than the conventional method.
[0096]
According to the embodiment of the present invention described above, the dielectric constant satisfying the relational expressions (7) to (9) is given to the lattice points near the boundary among the lattice points 22 and 23, By using the differential equation of the second accuracy defined by the equations (11) to (13) as the time difference equation, the boundary between the two media whose dielectric constant changes discontinuously is compared with the Yee lattice axis for simulation. Thus, it can be seen that a highly accurate electromagnetic field simulation for the two-dimensional TM mode can be realized even when tilted at a predetermined angle.
[0097]
【The invention's effect】
As is clear from the above description, according to the present invention, in a two-dimensional plane having a linear boundary where the dielectric constant of the medium changes spatially discontinuously, the boundary is a lattice-like shape in the two-dimensional plane. Electromagnetic field simulation method that enables high-accuracy simulation of electromagnetic field components propagating in a discrete grid-like plane when not parallel to any of the mutually orthogonal lattice axes A field simulation program and a program recording medium can be provided.
[Brief description of the drawings]
FIG. 1 is a functional block diagram of an electromagnetic field simulation apparatus according to an embodiment of the present invention.
FIG. 2 is an explanatory diagram showing a schematic configuration of a Yee lattice used in an electromagnetic field simulation method according to an embodiment of the present invention.
3 is an explanatory diagram showing a setting example (dielectric constant setting method I) of dielectric constants assigned to the second type and third type lattice points of the Yee lattice in FIG. 2. FIG.
FIG. 4 is a diagram illustrating a dielectric constant of a lattice point near the boundary of a Yee lattice assigned according to a dielectric constant setting method I.
5 is an explanatory diagram showing a setting example (dielectric constant setting method II) of dielectric constants assigned to the second type and third type lattice points of the Yee lattice of FIG. 2. FIG.
6 is a diagram illustrating a dielectric constant of a lattice point near the boundary of a Yee lattice assigned according to a dielectric constant setting method I. FIG.
FIG. 7 is a flowchart showing an outline of an electromagnetic field simulation method according to an embodiment of the present invention.
FIG. 8 is an explanatory diagram showing a simulation region used as an embodiment of the present invention.
9 is a diagram showing dielectric constants given to second and third type lattice points in the vicinity of the boundary between the semiconductor slab waveguide and the background in the simulation region of FIG. 8;
FIG. 10 is a diagram showing a result of evaluating a phase velocity error by an electromagnetic field simulation method according to an embodiment of the present invention.
[Explanation of symbols]
1 Electromagnetic field simulation equipment
11 Input section
12 Calculation unit
13 Memory unit
14 Main memory
15 Auxiliary storage
16 Output section
20 Yee lattice
21, 22, 23 lattice points
24, 53, 54 border
50 Simulation area
51 Semiconductor slab waveguide
52 Background
210, 510 straight line
220, 520 curve

Claims (9)

Two media having dielectric constants ε ₁ and ε ₂ in one plane are in contact with each other with one straight line as a boundary, and an electromagnetic field whose value of the magnetic field component parallel to the plane is zero exists in the plane. sometimes,
A plurality of first-type, second-type, and third-type grid points having different electromagnetic fields assigned to each other are used for the discretized virtual plane corresponding to the plane. A square shape of Δ is located at the midpoint between the two first type grid points closest to the second type grid point, and the midpoint between the two first type grid points closest to the third type grid point. Of these, the magnetic field component perpendicular to the plane is assigned to the first type lattice point while being arranged at a midpoint different from the middle point where the second type lattice point is located. Is assigned an electric field component orthogonal to the line segment connecting the first type lattice point closest to the second type lattice point and parallel to the plane, and further, the third type lattice point and the third type lattice point It is perpendicular to the line connecting the first-type grid points that are closest to each other and is flat on the plane. Constructed by assigning a field component,
An electromagnetic field simulation method for determining a time change of an electromagnetic field in the configured virtual plane by numerical calculation using a computer,
Any of the columns in which the first-type grid points and the third-type grid points are alternately arranged is defined as the x-axis, while any of the columns in which the first-type grid points and the second-type grid points are alternately arranged. When y is the y-axis and the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using a coordinate system in which a medium having a dielectric constant ε ₁ is positioned in the negative direction of the y-axis, The absolute value of the slope of the straight line is 2, and the y-intercept of the straight line is sΔ using a constant s taking an arbitrary real value in a range larger than −1.5 and smaller than 1.5 and the lattice spacing Δ. Is expressed as
A dielectric constant of a third-type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is −Δ is a,
The dielectric constant of a second type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is −Δ / 2 is b,
The dielectric constant of the third-type lattice points arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is 0 is c,
The dielectric constant of the second type lattice point arranged on the straight line parallel to the straight line representing the boundary and having a y intercept of Δ / 2 is d,
The dielectric constant of a third type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is Δ is e,
Of the second or third type lattice points to which any one of a, b, c, d, and e is not assigned, a lattice corresponding to a medium having a dielectric constant ε _{1 in} the plane corresponding to the virtual plane The dielectric constant of the point is ε ₁ ,
Of the second or third type lattice points to which any one of a, b, c, d, and e is not assigned, a lattice corresponding to a medium having a dielectric constant ε _{2 in} the plane corresponding to the virtual plane Let the dielectric constant of the point be ε ₂ ,
The computer is
(a) Relational expression

(Where the coefficient f _ij (i, j = 0,1,2) is

The coefficient g _ijk (i = 0,1; j, k = 0,1,2) is

The coefficient h _ij (i, j = 0,1) is

A step of dielectric constant a to meet to) defined, b, c, d, e , ε 1, sets the epsilon _2, stored in the memory and,
(b) setting an initial value of time n representing the discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory;
(c) reading out the dielectric constant and the magnetic field component set at each lattice point from the memory, calculating the electric field component at each lattice point at time n + 1/2, and storing it in the memory;
(d) reading out the electric field component from the memory, calculating the magnetic field component at each lattice point at time n + 1, and storing it in the memory;
(e) increasing time n by 1;
(f) The time n is compared with a predetermined threshold value, and if the time n has not reached the predetermined threshold value , the process returns to step (c) , and if the time n has reached the predetermined threshold value, the electric field An electromagnetic field simulation method comprising: executing a component and a step of reading out and outputting the magnetic field component from a memory . Two media having dielectric constants ε ₁ and ε ₂ in one plane are in contact with each other with one straight line as a boundary, and an electromagnetic field whose value of the magnetic field component parallel to the plane is zero exists in the plane. sometimes,
A plurality of first-type, second-type, and third-type grid points having different electromagnetic fields assigned to each other are used for the discretized virtual plane corresponding to the plane. A square shape of Δ is located at the midpoint between the two first type grid points closest to the second type grid point, and the midpoint between the two first type grid points closest to the third type grid point. Of these, the magnetic field component perpendicular to the plane is assigned to the first type lattice point while being arranged at a midpoint different from the middle point where the second type lattice point is located. Is assigned an electric field component orthogonal to the line segment connecting the first type lattice point closest to the second type lattice point and parallel to the plane, and further, the third type lattice point and the third type lattice point It is perpendicular to the line connecting the first-type grid points that are closest to each other and is flat on the plane. Constructed by assigning a field component,
An electromagnetic field simulation method for determining a time change of an electromagnetic field in the configured virtual plane by numerical calculation using a computer,
Any of the columns in which the first-type grid points and the third-type grid points are alternately arranged is defined as the x-axis, while any of the columns in which the first-type grid points and the second-type grid points are alternately arranged. When y is the y-axis and the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using a coordinate system in which a medium having a dielectric constant ε ₁ is positioned in the negative direction of the y-axis, The absolute value of the slope of the straight line is 2, and the constant δ taking an arbitrary real value in the range where the y-intercept of the straight line is larger than −1 and smaller than 1, and the lattice spacing Δ (1/2 + δ) When expressed as Δ,
Let u be the dielectric constant of the third-type lattice points arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is 0.
The dielectric constant of a second type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is Δ / 2 is v,
The dielectric constant of a third-type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is Δ is w,
Among the second or third type lattice points to which any one of u, v, and w is not assigned, the dielectric constant of the lattice point corresponding to the medium having a dielectric constant ε _{1 in} the plane corresponding to the virtual plane Ε ₁ ,
Among the second or third type lattice points to which any one of u, v, and w is not assigned, the dielectric constant of a lattice point corresponding to a medium having a dielectric constant ε _{2 in} the plane corresponding to the virtual plane Ε ₂
The computer is
(a) Relational expression

The permittivity u so as to satisfy, v, w, epsilon _1, sets the epsilon _2, and storing in a memory,
(b) setting an initial value of time n representing the discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory;
(c) reading out the dielectric constant and the magnetic field component set at each lattice point from the memory, calculating the electric field component at each lattice point at time n + 1/2, and storing it in the memory;
(d) reading out the electric field component from the memory, calculating the magnetic field component at each lattice point at time n + 1, and storing it in the memory;
(e) increasing time n by 1;
(f) The time n is compared with a predetermined threshold value, and if the time n has not reached the predetermined threshold value , the process returns to step (c) , and if the time n has reached the predetermined threshold value, the electric field An electromagnetic field simulation method comprising: executing a component and a step of reading out and outputting the magnetic field component from a memory . Step (f) 3. The electromagnetic wave according to claim 1, further comprising a step of calculating at least one of an amplitude, a phase, and a reflectance of the electromagnetic field based on the electric field component and the magnetic field component output in step 1. Field simulation method. Time difference equation having the secondary accuracy, and 1 to claim, characterized in that sequentially determine the electromagnetic field in response to changes in discrete time by the FDTD method using the spatial difference equation having the secondary precision 4. The electromagnetic field simulation method according to any one of 3 . Two media having dielectric constants ε ₁ and ε ₂ in one plane are in contact with each other with one straight line as a boundary, and an electromagnetic field whose value of the magnetic field component parallel to the plane is zero exists in the plane. sometimes,
A plurality of first-type, second-type, and third-type grid points having different electromagnetic fields assigned to each other are used for the discretized virtual plane corresponding to the plane. A square shape of Δ is located at the midpoint between the two first type grid points closest to the second type grid point, and the midpoint between the two first type grid points closest to the third type grid point. Of these, the magnetic field component perpendicular to the plane is assigned to the first type lattice point while being arranged at a midpoint different from the middle point where the second type lattice point is located. Is assigned an electric field component orthogonal to the line segment connecting the first type lattice point closest to the second type lattice point and parallel to the plane, and further, the third type lattice point and the third type lattice point It is perpendicular to the line connecting the first-type grid points that are closest to each other and is flat on the plane. Constructed by assigning a field component,
In order to obtain the time change of the electromagnetic field in the configured virtual plane by numerical calculation using a computer,
Any of the columns in which the first-type grid points and the third-type grid points are alternately arranged is defined as the x-axis, while any of the columns in which the first-type grid points and the second-type grid points are alternately arranged. When y is the y-axis and the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using a coordinate system in which a medium having a dielectric constant ε ₁ is positioned in the negative direction of the y-axis, The absolute value of the slope of the straight line is 2, and the y-intercept of the straight line is sΔ using a constant s taking an arbitrary real value in a range larger than −1.5 and smaller than 1.5 and the lattice spacing Δ. Is expressed as
A dielectric constant of a third-type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is −Δ is a,
The dielectric constant of a second type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is −Δ / 2 is b,
The dielectric constant of the third-type lattice points arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is 0 is c,
The dielectric constant of the second type lattice point arranged on the straight line parallel to the straight line representing the boundary and having a y intercept of Δ / 2 is d,
The dielectric constant of a third type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is Δ is e,
Of the second or third type lattice points to which any one of a, b, c, d, and e is not assigned, a lattice corresponding to a medium having a dielectric constant ε _{1 in} the plane corresponding to the virtual plane The dielectric constant of the point is ε ₁ ,
Of the second or third type lattice points to which any one of a, b, c, d, and e is not assigned, a lattice corresponding to a medium having a dielectric constant ε _{2 in} the plane corresponding to the virtual plane Let the dielectric constant of the point be ε ₂ ,
In the computer,
(a) Relational expression

(Where the coefficient f _ij (i, j = 0,1,2) is

The coefficient g _ijk (i = 0,1; j, k = 0,1,2) is

The coefficient h _ij (i, j = 0,1) is

A step of dielectric constant a to meet to) defined, b, c, d, e , ε 1, sets the epsilon _2, stored in the memory and,
(b) setting an initial value of time n representing the discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory;
(c) reading out the dielectric constant and the magnetic field component set at each lattice point from the memory, calculating the electric field component at each lattice point at time n + 1/2, and storing it in the memory;
(d) reading out the electric field component from the memory, calculating the magnetic field component at each lattice point at time n + 1, and storing it in the memory;
(e) increasing time n by 1;
(f) The time n is compared with a predetermined threshold value, and if the time n has not reached the predetermined threshold value , the process returns to step (c) , and if the time n has reached the predetermined threshold value, the electric field An electromagnetic field simulation program for executing a component and a step of reading out and outputting the magnetic field component from a memory . Two media having dielectric constants ε ₁ and ε ₂ in one plane are in contact with each other with one straight line as a boundary, and an electromagnetic field whose value of the magnetic field component parallel to the plane is zero exists in the plane. sometimes,
A plurality of first-type, second-type, and third-type grid points having different electromagnetic fields assigned to each other are used for the discretized virtual plane corresponding to the plane. A square shape of Δ is located at the midpoint between the two first type grid points closest to the second type grid point, and the midpoint between the two first type grid points closest to the third type grid point. Of these, the magnetic field component perpendicular to the plane is assigned to the first type lattice point while being arranged at a midpoint different from the middle point where the second type lattice point is located. Is assigned an electric field component orthogonal to the line segment connecting the first type lattice point closest to the second type lattice point and parallel to the plane, and further, the third type lattice point and the third type lattice point It is perpendicular to the line connecting the first-type grid points that are closest to each other and is flat on the plane. Constructed by assigning a field component,
In order to obtain the time change of the electromagnetic field in the configured virtual plane by numerical calculation using a computer,
Any of the columns in which the first-type grid points and the third-type grid points are alternately arranged is defined as the x-axis, while any of the columns in which the first-type grid points and the second-type grid points are alternately arranged. When y is the y-axis and the y-coordinate of the point forming the boundary is expressed by a linear expression of the x-coordinate using a coordinate system in which a medium having a dielectric constant ε ₁ is positioned in the negative direction of the y-axis, The absolute value of the slope of the straight line is 2, and the constant δ taking an arbitrary real value in the range where the y-intercept of the straight line is larger than −1 and smaller than 1, and the lattice spacing Δ (1/2 + δ) When expressed as Δ,
Let u be the dielectric constant of the third-type lattice points arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is 0.
The dielectric constant of a second type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is Δ / 2 is v,
The dielectric constant of a third-type lattice point arranged on a straight line that is parallel to the straight line representing the boundary and whose y-intercept is Δ is w,
Among the second or third type lattice points to which any one of u, v, and w is not assigned, the dielectric constant of the lattice point corresponding to the medium having a dielectric constant ε _{1 in} the plane corresponding to the virtual plane Ε ₁ ,
Among the second or third type lattice points to which any one of u, v, and w is not assigned, the dielectric constant of a lattice point corresponding to a medium having a dielectric constant ε _{2 in} the plane corresponding to the virtual plane Ε ₂
In the computer,
(a) Relational expression

The permittivity u so as to satisfy, v, w, epsilon _1, sets the epsilon _2, and storing in a memory,
(b) setting an initial value of time n representing the discretized time, setting an electric field component of the electromagnetic field at each lattice point and an initial value of the magnetic field component, and storing them in a memory;
(c) reading out the dielectric constant and the magnetic field component set at each lattice point from the memory, calculating the electric field component at each lattice point at time n + 1/2, and storing it in the memory;
(d) reading out the electric field component from the memory, calculating the magnetic field component at each lattice point at time n + 1, and storing it in the memory;
(e) increasing time n by 1;
(f) The time n is compared with a predetermined threshold value, and if the time n has not reached the predetermined threshold value , the process returns to step (c) , and if the time n has reached the predetermined threshold value, the electric field An electromagnetic field simulation program for executing a component and a step of reading out and outputting the magnetic field component from a memory . Step (f) Calculating at least one of amplitude, phase and reflectance of the electromagnetic field based on the electric field component and the magnetic field component output in The electromagnetic field simulation program according to claim 5 or 6, wherein the electromagnetic field simulation program is executed. Time difference equation having the secondary accuracy, and claims 5 to wherein sequentially determining the electromagnetic field in response to changes in discrete time by the FDTD method using the spatial difference equation having the secondary precision The electromagnetic field simulation program according to any one of claims 7 to 9 . Program recording medium characterized by recording the electromagnetic field simulation program as claimed in any one of claims 5 to 8.

Images