JP4984464B2 - 電磁界シミュレータおよび電磁界シミュレートプログラム - Google Patents

Description

本発明は、電磁界をシミュレートする電磁界シミュレータ、およびコンピュータに電磁界をシミュレートさせる電磁界シミュレートプログラムに関する。

従来より、通信アンテナの性能などを解析するために、通信電波などを対象とした電磁界をシミュレートする電磁界シミュレータや電磁界シミュレートプログラムが知られている。また、このような電磁界シミュレータや電磁界シミュレートプログラムで電磁界をシミュレートする代表的な手法として、有限差分時間領域法（ＦＤＴＤ法：Ｆｉｎｉｔｅ Ｄｉｆｆｅｒｅｎｃｅ Ｔｉｍｅ Ｄｏｍａｉｎ法）が知られている（例えば、非特許文献１および非特許文献２参照。）。この手法は、電磁界の時間変化を記述する基本方程式であるＭａｘｗｅｌｌの方程式を、空間的・時間的に差分化し、電磁界の時間変化を追跡するものである。この手法によれば、空間及び時間の離散化に用いるグリッドの間隔（ステップ）が十分に小さく設定されることにより、電磁界の時間変化が詳細にシミュレートされる。このようなＦＤＴＤ法の利点としては、計算原理が単純で演算の高速化が容易である点や、原理的に時間波形が計算されるため過渡的な電磁波特性が評価できる点や、３次元計算が容易である点などがあげられる。

電磁界成分の時間変化をＦＤＴＤ法で算出する式を以下に示す。（Ａｌｌｅｎ Ｔａｆｌｏｖｅ， Ｓｕｓａｎ Ｃ． Ｈａｇｎｅｓｓ ： “Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｅｌｅｃｔｒｏｄｙｎａｍｉｃｓ”より）

Ｅｘ，Ｅｙ，Ｅｚ：電界各成分
Ｈｘ，Ｈｙ，Ｈｚ：磁界各成分
Ｊsourceｘ：電流源
Ｍsourceｘ：磁流源
ε：誘電率
σ：導電率
σ＊：磁気伝導率
Δｘ，Δｙ，Δｚ：空間グリッド幅
Δｔ：時間分割幅
電磁界成分の右下に付された添字ｉ，ｊ，ｋは空間グリッドの番号であって、シミュレート空間の座標を表している。また、電磁界成分に付された添字ｎ，ｎ＋１は時間ステップを示す番号で、時間ステップが１進むと時間がΔｔだけ進む。
Yee, K.S., "Numerical solution of initial boundary value problems involving Maxwell'sequations in isotropicmedia", IEEE Transa. Antennas Propagat., Vol.Ap-14,p.302-307,1996. 宇野亨，「ＦＤＴＤ法による電磁界およびアンテナ解析」，１９９８，コロナ社

差分式の計算方法であるＦＤＴＤ法で十分な計算精度が得られるためには、シミュレーションにおける時間ステップおよび空間ステップが、それぞれ、光の周期および波長の１／１０〜１／２０程度に小さく設定されていることが望ましい。また、実際のシミュレーション時には、解析するモデルの各グリッドについて電界と磁界のＸ，Ｙ，Ｚ成分それぞれの瞬間値がＲＡＭなどからなる主記憶装置に保存されることが必要である。計算機の記憶容量には制限があるために、一般に、モデルの大きさは電磁波の波長の数倍に制限される。

光は波長が１ミクロン以下の電磁波であるため、ＦＤＴＤ法で計算が可能な領域が数ミクロンに限られ、これまではＦＤＴＤ法を光のシミュレーションに用いる事例は少なかったが、近年、光で情報記憶媒体に対する情報アクセスを行う技術が発達し、光の波長より小さい構造を持つ光学素子を用いる近接場光学などが注目されてくるに連れて、光の波長の数倍程度というミクロな領域における光の振る舞いを解析することが増えてきた。このような近接場光学では回折の効果が大きく、光を光線として扱うことができないため、ＦＤＴＤ法などを用いた電磁界シミュレータや電磁界シミュレートプログラムを利用して、光を電磁波として解析することが有効であり、ＦＤＴＤ法を光学モデルに適用した近接場光学の研究や光学素子の開発設計が盛んに行われている。

ところで、通信電波などのシミュレーションでも光のシミュレーションでも、電磁波の源を設定することが必要となるが、通信電波などの場合には、アンテナや回路に流れる電流を源として設定することができる。即ち、上述した式（１）、式(２)の電流源Ｊsourceｘに回路などの電流を代入し、磁流源Ｍsourceｘを「０」とすることで所望のシミュレーションが実現する。

一方、ＦＤＴＤ法で光学問題を計算する場合には、解析領域(シミュレート空間)の外側から光が入射する場合が多いため、入射する光を再現するためには、シミュレート空間の一部分で電界を入射波の周波数で強制振動させる方式が用いられる。電界を強制的に振動させた部分を波源と呼ぶ。この波源は現実の解析対象には存在せず、解析を行うために仮想的に設定したものである。

最も簡単な実現方式は、波源の部分で式(1)を適用せずに、以下に示す式(a)に従って電界Exを強制振動させる。

また、波源の部分以外では式(1)、式(2)の演算を、電流源Ｊsourceｘも磁流源Ｍsourceｘも「０」として実行することでシミュレート空間内の電界、磁界に波源部分の電界Ex振動が伝わり、光の伝搬がシミュレートされる。しかし、このように電界Exを強制振動させると、波源が配置されたことにより光が反射される現象がおきる。

図１は、電界Exを強制振動させる方式でのシミュレーション例を示す図である。

この図には、上方に直線状の波源１が配置され、波源１と相対する位置には反射面２が配置されたモデルが示されている。またこのモデルでは、波源１には、反射面２に向かって収束する球面波に相当する電流源が与えられている。

図１に示すシミュレーション結果には、波源１と反射面２との間に、反射面２と平行な干渉縞が現れている。これは、反射面２で反射されて波源１に戻った光が、波源１でさらに反射されて多重反射を引き起こしたことによるもので、仮想的な波源１の存在が計算に影響を及ぼした、間違った結果である。

波源１によるこのような反射を避ける手法は既に知られており、上記の式(a)に替えて

という式（３）を用いることが提案されている。ここでＥ０は波源における入射波の電界振幅、ωは入射波の角周波数、φは波源における入射波の初期位相であり、入射波の電界振幅Ｅ０は座標の関数となっている。

この式(３)では、波源に設定される電界として、入射波として振動する電界成分に、波源外からの光の伝搬を計算した場合の電界成分を足し合わせた値が用いられている。このため、波源中でも光の伝搬が計算されるので反射が発生しない。

図２は、波源中での光の伝搬が考慮されたシミュレーション例を示す図である。

ここでも、図１で用いられたモデルと同じモデルが用いられている。

この図２に示すシミュレーション結果では、波源１による反射がなくなったため、波源１と反射面２の間には、波源１から反射面２に進む光と、反射面２から波源１に向かう反射光との１回だけの干渉によって円弧状の縞模様が現れている。

上述した式(a)のように電流源を強制振動させて反射が起きる波源はハードソースと称され、式(３)のように入射波の振動に伝搬項を加えた波源はソフトソースと称される。光学解析ではソフトソースを用いるのが一般的である。

ところが、式(３)のソフトソースを用いて設定した波源の場合、シミュレーションの前提として期待された入射波の振幅がシミュレーション結果には正しく反映されず、入射光強度を正しく設定することができないという問題がある。

図３は、入射光強度を確認するためのシミュレーションの一例を示す図であり、図４は、入射光強度を確認するためのシミュレーションの別の例を示す図である。

図３では、波源における入射波の電界を振幅１Ｖ／ｍの正弦波で振動させた場合に屈折率１．５の媒質中を伝搬する光の電界振動がシミュレートされている。

計算条件として、光の伝搬方向のグリッドの大きさを１０ｎｍ、時間ステップ１３．３４３ａｓに設定した第１の条件と、グリッドの大きさ２０ｎｍ、時間ステップ１８．７９２ａｓに設定した第２の条件とが用いられており、図３には、各計算条件におけるシミュレーション結果が示されている。ここに示されたシミュレーション結果では、波源の振幅や物理的な解析モデルは不変であるにも関わらず、計算条件の違いにより振幅が変化している。しかも、どちらのシミュレーション結果で得られた電界振幅も、本来得られるべき振幅１Ｖ／ｍとは一致していない。

図４では、入射波の電界としてピーク電界１Ｖ／ｍ、パルス幅１６．７ｎｓのガウス型パルスが設定され、波源での屈折率が異なる４種類の場合それぞれについて、伝搬する光の電界振動がシミュレートされている。なお、それら４種類の場合について、いずれもグリッドサイズは１０ｎｍ、空間ステップは１０ａｓで共通である。

この図４に示すシミュレーション結果でも、図３に示す正弦波の場合と同様に、波源における電界の設定が一定にも係わらず、屈折率が異なると、伝搬した電界振動の振幅が異なっている。

本発明は上記事情に鑑み、電磁波の振幅を正確にシミュレートすることができる電磁界シミュレータおよび電磁界シミュレートプログラムを提供することを目的とする。

上記目的を達成する本発明の電磁界シミュレータは、
ある時刻における電磁界の空間分布を用いて、次の時刻における電磁界の空間分布を算出することを繰り返す電磁界シミュレータにおいて、
電磁界のシミュレーションにおける時空間的な粗さを設定する粗密度設定部と、
シミュレート空間上の各位置における電磁的な物理係数を設定する係数設定部と、
電磁界の空間分布に対する境界条件として電磁波の波源位置と強度とを設定する波源設定部と、
上記波源設定部で設定された波源位置と強度とを有する電磁波がシミュレーションの結果上で再現されるように、上記粗密度設定部で設定された時空間的な粗さおよび上記係数設定部で設定された電磁的な物理係数に基づいて、シミュレート空間上で波源の位置に付与される電磁界成分の強度を算出する強度算出部とを備えたことを特徴とする。

本発明の電磁界シミュレータによれば、波源の位置に付与される電磁界成分の強度が時空間的な粗さおよび電磁的な物理係数に基づいて算出されるので、所望の電磁波振幅がシミュレーション上で正確に再現される。

本発明の電磁界シミュレータは、電磁界のシミュレート中に、シミュレート空間上の波源位置に、上記強度算出部で得られた強度の電磁界成分と波源位置外から伝搬してくる電磁界成分との和を付与し続ける波源形成部を備えることが好ましい。

このような波源形成部を備えることにより、波源による電磁波の反射が回避される。

また、本発明の電磁界シミュレータは、上記強度算出部が、上記係数設定部で波源の位置に設定された誘電率εおよび透磁率μと、上記粗密度設定部で設定された時間粗さΔｔおよび電磁波進行方向の空間粗さδｚとに基づいて、波源の位置に付与される電界成分の強度Ｅ０を
Ｅ０＝２Δｔ／（δｚ√εμ）
と算出するものであることが好ましい。

このような式を用いることで、所望の強度の光を正確に再現する電界成分を容易に得ることができ、シミュレーション上で所望の光を再現することができる。

上記目的を達成する本発明の電磁界シミュレートプログラムは、
コンピュータに組み込まれて、そのコンピュータに、
ある時刻における電磁界の空間分布を用いて、次の時刻における電磁界の空間分布を算出することを繰り返させる電磁界シミュレートプログラムにおいて、
上記コンピュータ上に、
電磁界のシミュレーションにおける時空間的な粗さを設定する粗密度設定部と、
シミュレート空間上の各位置における電磁的な物理係数を設定する係数設定部と、
電磁界の空間分布に対する境界条件として電磁波の波源位置と強度とを設定する波源設定部と、
上記波源設定部で設定された波源位置と強度とを有する電磁波がシミュレーションの結果上で再現されるように、上記粗密度設定部で設定された時空間的な粗さおよび上記係数設定部で設定された電磁的な物理係数に基づいて、シミュレート空間上で波源の位置に付与される電磁界成分の強度を算出する強度算出部とを構成させることを特徴とする。

本発明の電磁界シミュレートプログラムによれば、本発明の電磁界シミュレータを構成する要素をコンピュータによって容易に実現することができる。

なお、本発明の電磁界シミュレートプログラムがコンピュータ上に構成する粗密度設定部などといった要素は、１つの要素が１つのプログラム部品によって構築されるものであってもよく、１つの要素が複数のプログラム部品によって構築されるものであってもよく、複数の要素が１つのプログラム部品によって構築されるものであってもよい。また、これらの要素は、そのような作用を自分自身で実行するものとして構築されてもよく、あるいは、コンピュータに組み込まれている他のプログラムやプログラム部品に指示を与えて実行するものとして構築されても良い。

以上説明したように、本発明によれば、電磁波の振幅を正確にシミュレートすることができる。

以下図面を参照して本発明の実施の形態を説明する。

ここでは、本発明の電磁界シミュレートプログラムの一実施形態がコンピュータに組み込まれて実行されることによってそのコンピュータ上に本発明の電磁界シミュレータの一実施形態が構成される例について説明する。

図５は、本発明の一実施形態が適用されるコンピュータを示す外観斜視図である。

このコンピュータ１００は、ＣＰＵ、ＲＡＭメモリ、ハードディスク等を内蔵した本体部１０１、本体部１０１からの指示により蛍光面１０２ａに画面表示を行うＣＲＴディスプレイ１０２、このコンピュータ内にユーザの指示や文字情報を入力するためのキーボード１０３、蛍光面１０２ａ上の任意の位置を指定することによりその位置に応じた指示を入力するマウス１０４を備えている。

本体部１０１は、さらに、外観上、フレキシブルディスク２２（図５には図示せず；図６参照）やＣＤ−ＲＯＭ２０が装填されるフレキシブルディスク装填口１０１ａおよびＣＤ−ＲＯＭ装填口１０１ｂを有しており、その内部には、装填されたフレキシブルディスクやＣＤ−ＲＯＭ２０をドライブする、フレキシブルディスクドライブ１１４、ＣＤ−ＲＯＭドライブ１１５（図６参照）も搭載されている。

本実施形態では、ＣＤ−ＲＯＭ２００に本発明の電磁界シミュレートプログラムの一実施形態が記憶されており、このＣＤ−ＲＯＭ２００がＣＤ−ＲＯＭ装填口１０１ｂから本体部１０１内に装填され、ＣＤ−ＲＯＭドライブ１１５によりそのＣＤ−ＲＯＭ２００に記憶された電磁界シミュレートプログラムがこのコンピュータ１００のハードディスク内にインストールされる。このコンピュータ１００のハードディスク内にインストールされた電磁界シミュレートプログラムが起動されると、このコンピュータ１００上には、本発明の電磁界シミュレータの一実施形態が構築される。

図６は、図５に示すコンピュータのハードウェア構成図である。

ここに示すようにコンピュータ１００には、中央演算処理装置（ＣＰＵ）１１１、ＲＡＭ１１２、ハードディスクコントローラ１１３、フレキシブルディスクドライブ１１４、ＣＤ−ＲＯＭドライブ１１５、マウスコントローラ１１６、キーボードコントローラ１１７、およびディスプレイコントローラ１１８が備えられており、それらはバス１１０で相互に接続されている。

フレキシブルディスクドライブ１１４、ＣＤ−ＲＯＭドライブ１１５は、図５を参照して説明したように、フレキシブルディスク２１０、ＣＤ−ＲＯＭ２００が装填され、装填されたフレキシブルディスク２１０、ＣＤ−ＲＯＭ２００をアクセスするものである。

また、ここには、ハードディスクコントローラ１１３によりアクセスされるハードディスク２２０、マウスコントローラ１１６により制御されるマウス１０４、キーボードコントローラ１１７により制御されるキーボード１０３、およびディスプレイコントローラ１１８により制御されるＣＲＴディスプレイ１０２も示されている。

前述したように、ＣＤ−ＲＯＭ２００には電磁界シミュレートプログラムが記憶されており、ＣＤ−ＲＯＭドライブ１１５により、そのＣＤ−ＲＯＭ２００から電磁界シミュレートプログラムが読み込まれ、バス１１０を経由し、ハードディスクコントローラ１１３によりハードディスク２２０内に格納される。実際の実行にあたっては、そのハードディスク２２０内の電磁界シミュレートプログラムはＲＡＭ１１２上にロードされ、ＣＰＵ１１１により実行される。

図７は、本発明の電磁界シミュレートプログラムの一実施形態を示す図である。ここでは、この電磁界シミュレートプログラム３００は、ＣＤ−ＲＯＭ２００に記憶されている。

この電磁界シミュレートプログラム３００は、図５に示すコンピュータ１００内で実行され、そのコンピュータ１００を、電磁界をシミュレートする電磁界シミュレータとして動作させるものであり、計算モデル入力部３１０と計算条件入力部３２０と波源入力部３３０と振幅補正部３４０とＦＤＴＤ法計算部３５０と計算結果解析部３６０と計算結果出力部３７０とを有する。

この電磁界シミュレートプログラム３００の各要素の詳細については後述する。

図８は、本発明の電磁界シミュレータの一実施形態の機能ブロック図である。

この電磁界シミュレータ４００は、図７の電磁界シミュレートプログラム３００が、図５に示すパーソナルコンピュータ１００にインストールされて実行されることにより構成されるものである。

この電磁界シミュレータ４００は計算モデル入力部４１０と計算条件入力部４２０と波源入力部４３０と振幅補正部４４０とＦＤＴＤ法計算部４５０と計算結果解析部４６０と計算結果出力部４７０とから構成されており、計算モデル入力部４１０、計算条件入力部４２０、波源入力部４３０、振幅補正部４４０、ＦＤＴＤ法計算部４５０、計算結果解析部４６０、および計算結果出力部４７０は、図７に示す電磁界シミュレートプログラム３００を構成する、計算モデル入力部３１０、計算条件入力部３２０、波源入力部３３０、振幅補正部３４０、ＦＤＴＤ法計算部３５０、計算結果解析部３６０、および計算結果出力部３７０それぞれによってパーソナルコンピュータ１００上に構成されている。

このように、図８に示す電磁界シミュレータ４００の各要素は、図７に示す電磁界シミュレートプログラム３００の各要素に対応するが、図８の各要素は、図５に示すパーソナルコンピュータ１００のハードウェアとそのパーソナルコンピュータで実行されるＯＳやアプリケーションプログラムとの組合せで構成されているのに対し、図７に示す各要素はそれらのうちのアプリケーションプログラムのみにより構成されている点が異なる。

以下、図８に示す電磁界シミュレータ４００の各要素の詳細を説明することによって、図７に示す電磁界シミュレートプログラム３００の各要素も合わせて説明する。

計算モデル入力部４１０は、本発明にいう係数設定部の一例に相当し、シミュレーション空間に対して屈折率などといった物性値を設定することによってシミュレーションの対象となる物理的なモデルを入力するものである。

計算条件入力部４２０は、本発明にいう粗密度設定部の一例に相当し、シミュレーションにおける時空間的な計算条件、即ち、時間ステップ間隔および空間グリッド幅を入力するものである。

波源入力部４３０は、本発明にいう波源設定部の一例に相当し、上述した物理的なモデルに対する入射光を、電界振幅、角周波数、および初期位相の入力によって設定するものである。

これら計算モデル入力部４１０、計算条件入力部４２０、および波源入力部４３０は、ハードウェアとしては、図６に示すキーボード１０３、マウス１０４、およびＣＰＵ１１１が担っている。

振幅補正部４４０は、本発明にいう強度算出部の一例に相当し、物理的なモデルと計算条件とに基づいて、シミュレーションの結果に入射光が再現されるように電界振幅を補正する。

ＦＤＴＤ法計算部４５０は、電磁界を前述のＦＤＴＤ法でシミュレートするものであり、本発明にいう波源形成部の一例としての役割も担うものである。

ここで、ＦＤＴＤ法による電磁界のシミュレートと電界振幅の補正とについて詳細に説明する。このＦＤＴＤ法では、シミュレート空間が多数のグリッドに分割され、各グリッド上に電磁界の成分が配置される。また、グリッド上に配置された成分を記憶するための３次元配列変数が図６に示すＲＡＭ１１２に用意される。

図９は、グリッド上に配置された電磁界の成分を表す図である。

ここには、空間が多数のグリッドに分割されたうちの１つのグリッド５００が示されており、このグリッド５００のサイズは、空間が分割された空間ステップΔｘ，Δｙ，Δｚに等しい。グリッド５００の頂点５１０には、３次元配列変数の配列番号が対応付けられており、この図５の左下の頂点５１０が配列番号（ｉ、ｊ、ｋ）に対応している。

ＦＤＴＤ法では、グリッド５００の各辺の中点５２０に、辺に平行な電界成分が配置され、グリッド５００の各面の中心点５３０に、面に垂直な磁界成分が配置される。このため、これらの電磁界成分を記憶するための３次元配列変数の配列番号には値「１／２」が用いられて、辺の中点５２０は（ｉ＋１／２、ｊ、ｋ）や（ｉ、ｊ＋１／２、ｋ）などという配列番号で示され、面の中心点５３０は（ｉ＋１／２、ｊ＋１／２、ｋ）や（ｉ、ｊ＋１／２、ｋ＋１／２）などという配列番号で示される。このような配列番号は空間座標に相当する。

このような電磁界成分の時間変化を算出する式を以下に示す。

Ｅｘ，Ｅｙ，Ｅｚ：電界各成分
Ｈｘ，Ｈｙ，Ｈｚ：磁界各成分
ε：誘電率
σ：導電率
σ＊：磁気伝導率
Δｘ，Δｙ，Δｚ：空間グリッド幅
Δｔ：時間分割幅
Ｅ０：波源における入射波の電界振幅
Ａ：電界振幅の補正係数
ω：入射波の角周波数
φ：波源における入射波の初期位相
電磁界成分の右下に付された添字ｉ，ｊ，ｋはグリッドの番号に相当する配列番号であって、シミュレート空間の座標を表している。また、電磁界成分に付された添字ｎ，ｎ＋１は時間ステップを示す番号で、時間ステップが１進むと時間がΔｔだけ進む。

上記式（３）’は、「発明が解決しようとする課題」の欄で示したソフトソースの式(３)を改良した式であり、電界振幅Ｅ０に対して補正係数Ａが掛かっている。

上記式（３）’が示すように、時間ステップ「ｎ＋１／２」における電界のｘ成分は、１ステップ過去の時間ステップ「ｎ−１／２」における電界のｘ成分と、１／２ステップ過去の時間ステップ「ｎ」における磁界のｙ成分およびｚ成分と、時間ステップ「ｎ」における電流源のｘ成分とから算出される。このとき、磁界のｙ成分およびｚ成分としては、図9に示すグリッド５００上で電界のｘ成分Ｅｘを取り囲んでいる磁界のｙ成分およびｚ成分が用いられる。上記式（２）が示すように、磁界のｘ成分についても同様に、過去の電磁界成分などから算出される。なお、電界成分と磁界成分とでは、時間ステップが互いに１／２ステップずれており、上記式（３）’式（２）は交互に実行されることとなる。

これら式（３）’式（２）は電磁界のｘ成分を算出する式であるが、電磁界のｙ成分やｚ成分についても同様の式によって算出される。これらの式に従った演算が図８のＦＤＴＤ法計算部４５０で実行されることによって、図６に示すＲＡＭ１１２に用意された３次元配列変数には、電磁界ベクトルの場を表したデータが記憶されることとなる。

次に、式（３）’における補正係数Ａの算出について説明する。

ソフトソースでは波源内を伝搬する光も同時に計算するため、時間ステップ間隔やセルサイズなどの計算条件、および波源が配置された場所の屈折率などの物性値により結果が変化する。従来のソフトソースの計算式（３）を見ると、右辺の電界と磁界にかかる係数にはε，σ，σ＊といった物性値やΔｘ，Δｔといった、計算条件で定義される変数が含まれている。このことは、次の時間ステップの電界を計算するときに、参照する電界や磁界の寄与が物理モデルと計算条件との双方の影響を受けることを意味する。従って波源に対して強制的に与えられる電界の振幅についても物理モデルおよび計算条件に応じた補正を施す必要がある。

この補正係数を求めるための基礎の式として、ここでは、波源の位置に仮想的な電流を発生させる式(１)を用いて考察を進める。式（１）と式（２）（磁流源は「０」）より、電流から磁界が発生し、さらにこの磁界から電界が発生する物理現象がシミュレートされる。実際には波源は真空中やガラスなど絶縁体の中に配置する場合が多いので、この電流はあくまで仮想的なものである。

また、計算の便宜上、入射光として、Ｚ方向に進行する直線偏光で電界のＸ成分がＥ０、Ｙ成分が「０」である光を考える。また、波源の位置には吸収がない、すなわち導電率σが「０」であるとする。電流が流れるとアンペールの法則に従って、電流を中心に右ねじ回りの磁界が発生する。磁界と電界の関係を式で表すと、

となる。次に、このアンペールの法則を、図９で説明したグリッド上に当てはめる。

図１０は、アンペールの法則を、図９で説明したグリッド上に当てはめた様子を表す図である。

この図１０にはＦＤＴＤ法における解析空間のグリッド上で、電流Ｊｘ’と、その電流Ｊｘ’を取り巻く磁界Ｈｙ，Ｈｚが示されている。これらの電流Ｊｘ’および磁界Ｈｙ，Ｈｚにアンペールの法則の式（４）を当てはめると次の式（５）が得られる。

ここで磁界の向きは軸の＋方向が正の値である。

一方、Ｚ方向に進み、電界振幅がＥｘ＝Ｅ０，Ｅｙ＝０となる平面波の磁界成分は、媒質の誘電率と透磁率を用いて次の関係式（６）を満たす。

上式（５）から電流と磁界の関係が得られ、式（６）から磁界と電界の関係が得られるので、これら式（５），式（６）を併せると、電界振幅Ｅ０を発生する電流の大きさは、

と求められる。この式（７）を、仮想的な電流源を用いる上述した式（１）に当てはめて導電率σを「０」とすると電流の項は、

となる。このように求められた電流の項と、ソフトソースの上記改良式（３）’の入射波の項とを比較することにより補正係数Ａが求められ、ＦＤＴＤ法の電界ｘ成分の式は

となる。この式（９）のうち補正係数Ａに相当する部分については図８に示す振幅補正部４４０でシミュレーションの実行前に計算が行われる。そして、ＦＤＴＤ法計算部４５０では、その計算された補正係数Ａの値が上記改良式（３）’に代入されて用いられることにより、実質的に式（９）によるシミュレーションが行われることとなる。なお、ここでの説明は簡単のため平面波で説明したが、どのような形態の波も局所的に見れば平面波で近似できるため、他の形態の入射波についても補正係数Ａの値は式（９）の第３項から得られる。また、波源における吸収はなく、導電率σが「０」であるという前提も用いたが、これも、波源に設定される部分で本来生じる吸収については別の項で改めて計算されることとなるので、波源のみで、かつ入射光のみが吸収なしであるとしてもシミュレーション全体に対しては特に影響はない。

このように補正係数Ａが加味されたシミュレーションが行われることにより、シミュレーションの結果には、図８に示す波源入力部４３０で入力されたとおりの入射光が、計算条件に関わらず正確に再現されることとなる。

図８に示す計算結果解析部４６０は、ＦＤＴＤ法計算部４５０によってシミュレートされた電磁界のデータを用いて電磁界の解析を行うものである。

計算結果出力部４７０は、計算結果解析部４６０による解析の結果をハードディスク２２０に保存する。またグラフに加工して、図５に示すＣＲＴディスプレイ１０２に表示させるものである。

ここで、電磁界シミュレータ４００の全体的な動作手順についてフローチャートを用いて改めて整理する。

図１１は、電磁界シミュレータ４００の全体的な動作手順を示すフローチャートである。

この図１１に示すフローチャートの説明に際しては、図８に示す要素を、特に図番を断らずに参照する場合がある。

電磁界シミュレータ４００の動作が開始されると、上記式（３）’式（２）で用いられる誘電率εや導電率σ等といった各種の物理係数（材料定数）がシミュレート対象のモデルに基づいて計算モデル入力部４１０で入力される（ステップＳ０１）。

次に、図９に示すグリッド５００における空間グリッド幅Δｘ，Δｙ，Δｚ、およびシミュレートされる電磁界の時間ステップ間隔Δｔが計算条件として計算条件入力部４２０で入力される（ステップＳ０２）。

さらに、シミュレート対象に対して入射する光の振幅Ｅ０、角周波数ω、初期位相φが波源入力部４３０で入力される（ステップＳ０３）。このとき振幅Ｅ０については、座標の関数として入力されることで平面波や球面波などが設定され、さらに時間の関数として入力されることでパルス波などが設定される。

その後、振幅補正部４４０によって、振幅Ｅ０に対する上述した補正係数Ａが算出され（ステップＳ０４）、その算出された補正係数Ａの値も含めて、上述した改良式（３）’の各係数値がＦＤＴＤ法計算部４５０に設定される（ステップＳ０５）。

次に、モデルにおける時間ｔが時間分割幅Δｔだけ進み（ステップＳ０６）、上記式（３）’式（２）に従ったＦＤＴＤ法の演算がＦＤＴＤ法計算部４５０で実行される（ステップＳ０７）。このとき、電界成分を算出する式（３）’の方が１／２ステップだけ先なので、演算は式（３）’、式（２）の順で実行される。また、入射波の振幅が時間の関数である場合にはその関数に従った時間変化の計算も実行される。なお、演算によって得られた電界成分のうち、後の解析で必要となるものについては、図６に示すハードディスク２２０などに適宜に保存される。

このような式（３）’、式（２）の演算が終了すると、次に、モデルにおける時間ｔが、シミュレート計算を終了する終了時刻ｔｍａｘを超えたか否かが判定され（ステップＳ０８）、終了時刻ｔｍａｘ以前であると判定された場合にはステップＳ０６に戻って上記の手順が繰り返される。

一方、ステップＳ０８で、モデルにおける時間ｔが終了時刻ｔｍａｘを超えたと判定された場合には、ＦＤＴＤ法計算部４５０はシミュレーション計算を終了し、シミュレーション結果の解析が計算結果解析部４６０によって行われる（ステップＳ０９）。ここでの解析内容については、従来から周知の各種解析が適宜採用されるものとして詳細な説明は省略する。

最後に、計算結果出力部４７０によって解析結果がグラフ化されて図５に示すＣＲＴディスプレイ１０２に表示され、またその結果をハードディスク２２０に保存し電磁界シミュレータ４００の動作が終了する。

以下、本実施形態の電磁界シミュレータ４００に具体的なモデルが適用されたシミュレーション例と解析例を用いて、電磁界シミュレータ４００によるシミュレーションの効果を説明する。

図１２は、図３に示すシミュレーション例を本実施形態に適用した場合の結果を示す図である。

上述したように、このシミュレーションでは、波源における入射波の電界を振幅１Ｖ／ｍの正弦波で振動させた場合に屈折率１．５の媒質中を伝搬する光の電界振動がシミュレートされている。また、計算条件として、光の伝搬方向のグリッドの大きさを１０ｎｍ、時間ステップ１３．３４３ａｓに設定した第１の条件と、グリッドの大きさ２０ｎｍ、時間ステップ１８．７９２ａｓに設定した第２の条件とが用いられている。各条件における補正係数の値は、真空の光速をｃ＝２．９９７９ｘ１０８［ｍ／ｓ］とすると、第１の条件では０．５３３５１となり、第２の条件では０．３７５５８となる。これらの補正係数の値が用いられてシミュレートされた結果では、図１２に示すように、どちらの計算条件でも同様の振幅の波形が得られており、しかも、波源で設定した１Ｖ／ｍの振幅と同じ振幅を有する伝搬光が正確に再現されている。

図１３は、図４に示すシミュレーション例を本実施形態に適用した場合の結果を示す図である。

上述したように、このシミュレーションでは、入射波の電界としてピーク電界１Ｖ／ｍ、パルス幅１６．７ｎｓのガウス型パルスが設定され、波源での屈折率が異なる４種類の場合それぞれについて、伝搬する光の電界振動がシミュレートされている。また、それら４種類の場合について、いずれもグリッドサイズは１０ｎｍ、空間ステップは１０ａｓで共通である。

この場合の補正係数は、屈折率に応じて異なる値となり、各屈折率に対応した補正係数がそれぞれ適用されることにより、図１３に示すように、４種類の場合すべてについて振幅が１Ｖ／ｍで同一となる解析結果が得られる。

このように、上述した補正係数の算出式は波源の特性によらず適用可能であり、正弦波だけでなく任意の波形に適用可能であることがわかる。

図１４および図１５は、厚さ４００ｎｍの金の薄膜に円形の微小開口があるモデルに対するシミュレーション例を示す図であり、図１４では開口径が直径１００ｎｍ、図１５では開口径が直径２００ｎｍである。

図１４，図１５の双方について、金薄膜３の上方の波源４から波長４０５ｎｍの平面波が入射した場合のシミュレーションが行われている。また、双方とも波源４の位置の屈折率は同じで、グリッドサイズも光伝播方向は同じ２０ｎｍとなっている。ただし薄膜３の面内方向のグリッドサイズは、開口の大きさに合わせて、開口径１００ｎｍのモデルでは２ｎｍ、開口径２００ｎｍのモデルでは５ｎｍに設定されている。

また、時間ステップについては、空間グリッドの最小間隔と時間ステップとの間で、ＣＦＬ（Ｃｏｕｒａｎｔ， Ｆｒｉｅｄｒｉｃｈ ａｎｄ Ｌａｖｙ）条件あるいはＣｏｕｒａｎｔ条件と呼ばれる関係を満たす必要があり、時間ステップが空間グリッドの最小間隔に対して十分に小さなステップとなっていないとＦＤＴＤ法による計算値が発散してしまい、正常な計算を行うことができない。このため、図１４に示す開口径１００ｎｍのモデルでは、面内方向のグリッドサイズが２ｎｍなので、時間ステップは４．７１ａｓ以下であることが必要であり、開口径２００ｎｍのモデルでは、面内方向のグリッドサイズが５ｎｍなので、時間ステップは１１．６ａｓ以下であることが必要がある。また、シミュレーションの便宜のため、光の周期が時間ステップの整数倍となるような設定が望ましく、計算時間の短縮などのために時間ステップは長い方が望ましい。このため、光の周期１．３５１ｆｓに対して、開口径が１００ｎｍのモデルでは時間ステップが４．６１ａｓ（１／２９３周期）に設定され、開口径が２００ｎｍのモデルでは時間ステップが１１．３５ａｓ（１／１１９周期）に設定されている。

これら２つのモデルに対するシミュレーションの解析結果を比較する。

図１６は、従来のシミュレーション技術で上記２つのモデルをシミュレーションした時の解析結果を比較した図である。

ここでは、解析結果として、開口中心に沿う各点における全電界の２乗値を、光の伝搬方向にプロットしたグラフが示されている。このグラフの横軸は光の伝搬方向の座標を示しており、波源はｚ＝＋４００ｎｍに位置し、薄膜はｚ＝０〜−４００ｎｍの間に位置している。光は＋Ｚ方向から−Ｚ方向に伝搬している。

電界の２乗値は光強度に比例するため、本来ならばこのグラフからは、開口を通過してくる光の強度が開口径に依存して変化する様子を評価できるはずである。しかし、図１６に示す例では、波源から薄膜へと伝搬する光（ｚ＝０〜＋４００ｎｍの間の光）の強度が２つのモデルで異なっているため、開口を通過してくる光の強度を比較することは意味がない。

正確な比較のためには、２つのモデルで時間ステップを同一にして、同一の計算条件で解析しなければならず、開口径２００ｎｍのモデルではＣＦＬ条件で制限される時間ステップのおよそ１／２で計算することとなり、計算時間が約倍になってしまう。また、先に開口径が２００ｎｍのモデルにおける解析結果を得ていて、その後、開口径が１００ｎｍのモデルに対する計算を行う必要ができたような場合には、開口径が１００ｎｍのモデルに対する計算を、開口径が２００ｎｍのモデルと同じ時間ステップで実行する事ができない（ＣＦＬ条件を満たさない）ため、開口径が２００ｎｍのモデルを再計算する事になる。

図１７は、本実施形態の電磁界シミュレータで上記２つのモデルをシミュレーションした時の解析結果を比較した図である。

上述した補正係数によって電界振幅が補正されているため、この図１７では、薄膜に入射する前の＋Ｚ領域における光強度が、２つのモデルでほぼ同じになっていることが分かる。従って、本実施形態の電磁界シミュレータによるシミュレーションの解析結果では、薄膜の円形開口に光が入射した後の減衰を、計算条件が異なる２つのモデル間で正確に比較することができる。また、これら２つのモデル以外のモデルを後から追加して計算する場合であっても正確な比較が可能である。

このように計算条件が異なっていても正確な比較ができるので、各モデルにおける計算に際しては、各モデルに適した計算条件を用いることができ、計算時間などの効率がよい。

なお、本明細書で例示したような簡単なモデルの場合には、シミュレーションの結果に対して補正を行うことも可能であるが、誘電率などに非線形性があるような場合には、シミュレーション前の振幅補正が不可欠となる。

また、上記説明では、本発明の電磁界シミュレートプログラムの実施形態によって、本発明の電磁界シミュレータの実施形態を構成する要素がコンピュータ上に構成される例が示されているが、本発明の電磁界シミュレータは、ハードウェアで構成された要素を有するものであってもよい。

また、上記説明では、電界成分について振幅が補正される例が示されているが、本発明にいう強度算出部は、磁界成分について補正後の強度を算出するものであってもよい。

電流源Ｊsourceｘを強制振動させる方式でのシミュレーション例を示す図である。 波源中での光の伝搬が考慮されたシミュレーション例を示す図である。 入射光強度を確認するためのシミュレーションの一例を示す図である。 入射光強度を確認するためのシミュレーションの別の例を示す図である。 本発明の一実施形態が適用されるコンピュータを示す外観斜視図である。 図５に示すコンピュータのハードウェア構成図である。 本発明の電磁界シミュレートプログラムの一実施形態を示す図である。 本発明の電磁界シミュレータの一実施形態の機能ブロック図である。 グリッド上に配置された電磁界の成分を表す図である。 アンペールの法則を、図９で説明したグリッド上に当てはめた様子を表す図である。 電磁界シミュレータの全体的な動作手順を示すフローチャートである。 図３に示すシミュレーション例を本実施形態に適用した場合の結果を示す図である。 図４に示すシミュレーション例を本実施形態に適用した場合の結果を示す図である。 厚さ４００ｎｍの金の薄膜に開口径が直径１００ｎｍの円形微小開口があるモデルに対するシミュレーション例を示す図である。 厚さ４００ｎｍの金の薄膜に開口径が直径２００ｎｍの円形微小開口があるモデルに対するシミュレーション例を示す図である。 従来のシミュレーション技術で上記２つのモデルをシミュレーションした時の解析結果を比較した図である。 本実施形態の電磁界シミュレータで上記２つのモデルをシミュレーションした時の解析結果を比較した図である。

符号の説明

１ 波源
２ 反射面
３ 薄膜
４ 波源
１００ コンピュータ
１０１ 本体部
１０２ａ 蛍光面
１０２ ＣＲＴディスプレイ
１０３ キーボード
１０４ マウス
１０１ａ フレキシブルディスク装填口
１０１ｂ ＣＤ−ＲＯＭ装填口
１１０ バス
１１１ 中央演算処理装置（ＣＰＵ）
１１２ ＲＡＭ
１１３ ハードディスクコントローラ
１１４ フレキシブルディスクドライブ
１１５ ＣＤ−ＲＯＭドライブ
１１６ マウスコントローラ
１１７ キーボードコントローラ
１１８ ディスプレイコントローラ
２００ ＣＤ−ＲＯＭ
２１０ フレキシブルディスク
２２０ ハードディスク
３００ 電磁界シミュレートプログラム
３１０ 計算モデル入力部
３２０ 計算条件入力部
３３０ 波源入力部
３４０ 振幅補正部
３５０ ＦＤＴＤ法計算部
３６０ 計算結果解析部
３７０ 計算結果出力部
４００ 電磁界シミュレータ
４１０ 計算モデル入力部
４２０ 計算条件入力部
４３０ 波源入力部
４４０ 振幅補正部
４５０ ＦＤＴＤ法計算部
４６０ 計算結果解析部
４７０ 計算結果出力部
５００ グリッド

Claims (4)

1. ある時刻における電磁界の空間分布を用いて、次の時刻における電磁界の空間分布を算出することを繰り返す電磁界シミュレータにおいて、
   電磁界のシミュレーションにおける時空間的な粗さを設定する粗密度設定部と、
   シミュレート空間上の各位置における電磁的な物理係数を設定する係数設定部と、
   電磁界の空間分布に対する境界条件として電磁波の波源位置と強度とを設定する波源設定部と、
   前記波源設定部で設定された波源位置と強度とを有する電磁波がシミュレーションの結果上で再現されるように、前記粗密度設定部で設定された時空間的な粗さおよび前記係数設定部で設定された電磁的な物理係数に基づいて、シミュレート空間上で波源の位置に付与される電磁界成分の強度を算出する強度算出部とを備えたことを特徴とする電磁界シミュレータ。
2. 電磁界のシミュレート中に、シミュレート空間上の波源位置に、前記強度算出部で得られた強度の電磁界成分と波源位置外から伝搬してくる電磁界成分との和を付与し続ける波源形成部を備えたことを特徴とする請求項１記載の電磁界シミュレータ。
3. 前記強度算出部が、前記係数設定部で波源の位置に設定された誘電率εおよび透磁率μと、前記粗密度設定部で設定された時間粗さΔｔおよび電磁波進行方向の空間粗さδｚとに基づいて、前記波源設定部で設定された電界成分の強度Ｅ０に対して、波源の位置に付与される電界成分の強度Eを
   Ｅ＝２Δｔ／（δｚ√εμ）Ｅ０
   と算出するものであることを特徴とする請求項１記載の電磁界シミュレータ。
4. コンピュータに組み込まれて、該コンピュータに、
   ある時刻における電磁界の空間分布を用いて、次の時刻における電磁界の空間分布を算出することを繰り返させる電磁界シミュレートプログラムにおいて、
   前記コンピュータ上に、
   電磁界のシミュレーションにおける時空間的な粗さを設定する粗密度設定部と、
   シミュレート空間上の各位置における電磁的な物理係数を設定する係数設定部と、
   電磁界の空間分布に対する境界条件として電磁波の波源位置と強度とを設定する波源設定部と、
   前記波源設定部で設定された波源位置と強度とを有する電磁波がシミュレーションの結果上で再現されるように、前記粗密度設定部で設定された時空間的な粗さおよび前記係数設定部で設定された電磁的な物理係数に基づいて、シミュレート空間上で波源の位置に付与される電磁界成分の強度を算出する強度算出部とを構成させることを特徴とする電磁界シミュレートプログラム。

Images