JP5113572B2 - シミュレーション方法及びシミュレーションプログラム - Google Patents

Description

本発明は、ＦＤＴＤ法を利用したシミュレーション方法及びシミュレーションプログラムに関するものである。

金属部材を含む構造物の一例としては、複数の金属ストリップ（金属細線）がほぼ平行に配列されて構成される金属回折格子がある。このような金属回折格子は、ワイヤグリッド偏光子、またはディスプレイ関連光学装置における偏光コントローラとしてしばしば利用される。各金属ストリップの幅または周期が、金属回折格子に入射する光の波長に比べて小さい場合、回折格子特性は、金属ストリップの幅および周期、入射ビームの偏光および入射角に依存する。そのため、所望の回折格子特性が得られるように、金属ストリップの幅及び周期等を入射ビームの特性などに応じて設計する必要がある。

金属部材を含む構造物へ電磁波が入射した場合の電磁界の解析には、解析的な手法を利用したシミュレーションが適用される。例えば、厳密結合波解析（RCWA :rigorous coupled wave analysis）などの解析的方法が無限の大きさを仮定した回折格子に対して使用されることがある。しかしながら、ＲＣＷＡに例示される解析的方法では、（１）有限の大きさを仮定すべき構造又は有限の幅を有する入射ビームの場合には使用することができない、（２）不規則な境界面を扱うことができない、（３）媒質の誘電率を複素数としているが、その実部は入射ビームの角周波数に依存していないことを仮定している、等の問題点がある。よって、例えば、上述した配列方向の長さが有限である回折格子、すなわち、実際の回折格子により近い形態を正確に解析できず、また、金属部材が、銀、アルミニウム、金などの貴金属からなる場合には、可視領域においてシミュレーションの結果で得られる電磁界分布と、実空間での電磁界分布との差が大きくなる。

一方、上記のような問題点を解決可能なシミュレーション手法として、不規則な境界面を有しており、有限の大きさの仮定が有効な構造物又は無限の仮定が可能な構造物のどちらにも適用可能なＦＤＴＤ法（Finite Differential Time Domain method）が知られている（例えば、特許文献１，２参照）。ＦＤＴＤ法は、ＦＤＴＤ領域を複数のセルに分割して、マクスウェル方程式内の空間および時間微分を有限差分によって近似することでマクスウェル方程式を直接解くものであり、様々な金属構造を解析するために使用されてきている。そして、ＦＤＴＤ法の場合、例えば、デバイ（Debye）、ドルーデ(Drude)、またはローレンツ（Lorentz）など、誘電率についての様々な解析的モデルが、周波数領域畳み込みを使用することによって、ＦＤＴＤに組み込むことができる。このような畳込みを利用した方法は、帰納的畳み込み（Recursive Convolution）法（以下、単にＲＣ法ともいう）として知られている。
特開２０００−１０５２５９号公報 特開２０００−２２７４５０号公報

しかし、従来のＦＤＴＤ法を利用したシミュレーションでは、空間的離散化の値等の条件設定によっては、シミュレーションが不安定になる傾向にあった。

そこで、本発明は、ＦＤＴＤ法によるシミュレーションをより安定して実施可能なシミュレーション方法及びシミュレーションプログラムを提供すること目的とする。

本発明のシミュレーション方法は、金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーション方法であって、構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、離散化条件決定工程で決定された空間離散化条件に基づいて、解析モデルに電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いたＦＤＴＤ法に従ってコンピュータにより算出する解析工程と、を備え、離散化条件決定工程では、帰納的畳み込み法において上記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ_１，ｃ_２であって、ＦＤＴＤ領域において所定波長を表すための格子点の数及びＦＤＴＤ領域における電磁波の解析用速度並びに金属部材の上記誘電率を含むプラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆを用いて下記式（１），（２）で表される変数ｃ_１，Ｆ，ｃ_２，Ｆ：

において、プラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆに含まれる誘電率を、所定波長の電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、変数ｃ_１，Ｆ，ｃ_２，Ｆが、ＦＤＴＤ法における計算の安定化条件としての下記式（３）及び式（４）を満たすと共に、

解析用速度が０．７より小さくなるように、格子点の数及び解析用速度をコンピュータによって決定することを特徴とする。

また、本発明に係るシミュレーションプログラムは、金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーションをコンピュータに実行させるためのシミュレーションプログラムであって、構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、離散化条件決定工程で決定された空間離散化条件に基づいて、解析モデルに電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いたＦＤＴＤ法に従ってコンピュータにより算出する解析工程と、をコンピュータに実行させ、離散化条件決定工程では、帰納的畳み込み法において上記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ_１，Ｆ，ｃ_２，Ｆであって、ＦＤＴＤ領域において所定波長を表すための格子点の数及びＦＤＴＤ領域における電磁波の解析用速度並びに金属部材の上記誘電率を含むプラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆを用いて下記式（５），（６）で表される前記変数ｃ_１，Ｆ，ｃ_２，Ｆ：

において、プラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆに含まれる上記誘電率を、所定波長の電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、変数ｃ_１，Ｆ，ｃ_２，Ｆが、ＦＤＴＤ法における計算の安定化条件としての下記式（７）及び式（８）を満たすと共に、

解析用速度が０．７より小さくなるように、格子点の数及び解析用速度をコンピュータに決定せしめることを特徴とする。

上記本発明に係るシミュレーション方法及びシミュレーションプログラムでは、所定波長λを有する電磁波に対する金属部材の誘電率の実測値及び実測値に基づいて決まる値を使用して、ＦＤＴＤ法における計算の安定化条件を満たすようにＦＤＴＤ領域の空間離散化を図っている。そのため、より正確なシミュレーションを安定して実施することが可能である。

また、本発明に係るシミュレーション方法では、実空間における電磁波の所定波長をλ、電磁波の速度をｖ及び所定波長に対応する角周波数をω_λとし、ＦＤＴＤ領域における格子点の数をλ_Ｆ、解析用速度をｖ_Ｆとし、角周波数ω_λに対する誘電率の実部及び虚部の値をε^ｒ及びε^ｉとし、角周波数ω_λに対応するＦＤＴＤ領域における解析用角周波数ω_Ｆを、

としたとき、衝突周波数ν_ｃ，Ｆ及びプラズマ角周波数ω_ｐ，Ｆは、

として与えられることが好ましい。

同様に、本発明に係るシミュレーションプログラムでは、実空間における電磁波の所定波長をλ、電磁波の速度をｖ及び所定波長に対応する角周波数をω_λとし、ＦＤＴＤ領域における格子点の数をλ_Ｆ、解析用速度をｖ_Ｆとし、角周波数ω_λに対する誘電率の実部及び虚部の値をε^ｒ及びε^ｉとし、角周波数ω_λに対応するＦＤＴＤ領域における解析用角周波数ω_Ｆを、

としたとき、衝突周波数ν_ｃ，Ｆ及びプラズマ角周波数ω_ｐ，Ｆは、

として与えられることが好ましい。

この場合、ν_ｃ、Ｆ及びω_Ｐ，Ｆは、λ_Ｆ及びｖ_Ｆの関数となる。そのため、変数ｃ_１，Ｆ及びｃ_２，Ｆが満たす条件と共に、ｖ_Ｆが０．７より小さい値をとるように、λ_Ｆ及びｖ_Ｆを決定することができる。

本発明に係るシミュレーション方法において、上記解析工程では、電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用してコンピュータにより解くことが好適である。同様に、本発明に係るシミュレーションプログラムにおいて、上記解析工程では、電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめることが好適である。

この場合、波動方程式を利用しているため計算量を低減することが可能である。

また、本発明に係るシミュレーション方法における解析工程では、電磁波におけるＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用して解く態様とすることもできる。同様に、本発明に係るシミュレーションプログラムにおける解析工程では、電磁波におけるＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめる態様とすることもできる。

更に、本発明に係るシミュレーション方法における上記解析工程では、電磁波がＴＭモード及びＴＥモードを共に含む場合には、ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用して解き、ＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用して解き、ＴＭモード及びＴＥモードに対して得られた電磁界に基づいてＦＤＴＤ領域における電磁界を算出する態様とすることもできる。同様に、本発明に係るシミュレーションプログラムにおける上記解析工程では、電磁波がＴＭモード及びＴＥモードを共に含む場合には、ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に上記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめて、ＴＭモードに対しては、マクスウェル方程式に帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、帰納的関係式を利用してコンピュータに解かさしめて、ＴＭモード及びＴＥモードに対して得られた電磁界に基づいてＦＤＴＤ領域における電磁界を算出する態様とすることができる。

本発明によれば、ＦＤＴＤ法に従ってシミュレーションを安定して実施することができる。

以下、図面を参照して本発明に係るシミュレーション方法及びシミュレーションプログラムの実施形態について説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。また、図面の寸法比率等は、説明のものと必ずしも一致していない。

本発明に係るシミュレーション方法及びシミュレーションプログラムは、金属部材を含む構造物に電磁波が入射した場合における電磁界を、ＦＤＴＤ法（Finite Difference Time Domain Method）に従って解析するためのものであり、上記構造物の設計などに好適に利用される。本発明に係る方法及びプログラムを利用した解析対象としての構造物の一例としては金属回折格子が挙げられる。本明細書では、特に断らない限り、光学素子は金属回折格子として説明する。

図１は、金属回折格子の模式図である。金属回折格子１０は、一方向（図１では、紙面に垂直な方向）に延在している複数の金属部材としての金属細線（金属ストリップ）１１が、金属細線１１の延在方向に対して略直交する方向に周期Λで配置されてなる回折格子である。説明の簡略化のために、図１に示すように、金属細線１１の配列方向をｘ軸方向とし、金属細線１１の延在方向をｚ軸方向とし、ｘ軸方向及びｚ軸方向に直交する方向をｙ軸方向とする。金属細線１１の材料としては、銀、金、アルミニウムなどが例示される。

金属回折格子１０の周期Λは、電磁波としての入射ビーム２０が有する所定波長λに対して十分小さく、例えば、５λ以下である。金属回折格子１０に入射ビーム２０が入射すると、入射ビーム２０は金属回折格子１０により回折される。周期Λが所定波長λに対して十分小さい場合、金属回折格子１０に入射ビーム２０が入射すると、ほぼ０次の回折光が生じることになる。以下の説明では、図１に示したように、金属回折格子１０に対して入射ビーム１０の入射側と反対側に回折される光を透過ビーム２１とも称し、入射ビーム２０の入射側と同じ側に回折される光を反射ビーム２２とも称す。なお、所定波長λは、電磁スペクトルにおける可視領域の波長であるとすることができる。

金属回折格子１０の特性は、周期Λ並びに金属細線１１の幅（ｘ軸方向の長さ）ｗ及び厚さ（ｙ軸方向の長さ）ｈや、回折格子法線Ｎと入射ビーム２０とのなす角である入射角θ_ｉｎや、入射ビーム２０の偏光状態などに依存する。そのため、所望の回折格子特性を実現する金属回折格子１０を設計するために本発明に係るシミュレーション方法及びプログラムが有効である。特に、金属回折格子１０のｘ軸方向の幅Ｗが３０λより小さい場合に本発明が有効である。これは、本発明に係るシミュレーション方法及びシミュレーションプログラムが、その手法として、有限な大きさを持つ物体を扱うときに有効なＦＤＴＤ法を採用していることによる。

図２は、図１に示した金属回折格子のＦＤＴＤ領域における解析モデルの模式図である。ＦＤＴＤ法では、図２に示すように、ＦＤＴＤ領域１００を複数のセルに分割することによって、マクスウェル方程式を空間的及び時間的に差分化して直接数値的に解く手法である。本実施形態では、図２に示すように、解析モデルとしての金属回折格子に対しても図１に示した金属回折格子１０と同様の符号を付すこととする。また、図２中のハッチングは、セルで分割されたＦＤＴＤ領域１００内において金属細線１１を示すための便宜的なものである。

通常、図１及び図２に示した金属回折格子１０のように一方向（図１，２中ではｚ軸方向）に延在しており、その延在方向に変化のないものを解析対象とする場合、２次元問題として扱うことができる。そこで、２次元のＦＤＴＤ法のアルゴリズムについて説明する。

ＦＤＴＤ法では、以下のマクスウェル方程式を利用する。

μは透磁率である。また、Ｅは電界、Ｈは磁界である。更に、「∇」は、ｘ_ｕ及びｙ_ｕをｘ軸方向及びｙ軸方向の単位ベクトルとしたとき、ｘ_ｕ(∂／ｘ)＋ｙ_ｕ(∂／ｙ)で定義される微分演算子である。

式（１６）中のＤは電束密度を表しており、次式で与えられる。

式（１７）中のε（ｒ，ω）は金属細線１１の誘電率である。

時間領域における電束密度Ｄは、式（１７）の両辺でフーリエ変換をとり、畳み込み定理を適用することによって次式で表される。

式（１８）において、ε（ｔ）＝Ｆ｛ε（ω）｝であり、Ｆはフーリエ変換を表す。

誘電率ε（ω）が角周波数ωのある特殊な関数で表されるとき、式（１８）の畳み込み積分は帰納的に計算可能である。これはＲＣ法として知られている。

本実施形態では、金属細線１１の誘電率ε（ω）を表す関数として、ドルーデの一次近似式（ドルーデの分散式）を採用する。この場合、誘電率ε（ω）は式（１９）のように表される。

式（１９）において、ω_ｐはプラズマ周波数、ν_ｃは衝突周波数である。ε_∞は、周波数が無限大のときの誘電率であり１に等しい。χ（ω）は電気比感受率である。

χ（ω）をフーリエ変換したもの、すなわち、Ｆ｛χ（ω）｝をχ（ｔ）と表した場合、χ（ｔ）は次式で表される。

Ｕ（ｔ）は、時間の単位ステップ関数であって、ｔ＝０の場合、Ｕ（ｔ）＝０であり、ｔ＞０の場合、Ｕ（ｔ）＝１である。

Δｔを時間ステップ間隔としたとき、時刻ｔ＝ｎΔｔ（ｎは０以上の整数）において、式（１８）の畳み込み積分を評価した場合、

が得られる。より具体的には、式（２１）は、区間［ｍΔｔ，（ｍ＋１）Δｔ］にわたってＥを定数としてとることによって、畳み込み積分から得られる。なお、ｍは０以上の整数である。

式（２１）は、以下のように書き表すことができる。

ここで、Ψは、累積変数であり下記式（２３）が成り立つ。

式（２３）は、

と表すことができる。この式（２４）において、

である。式（２５）より、

が成立する。式（２６）において、

である。また、Ψ^０＝０、Ψ^１＝Ｅ^１χ^０、Ψ^２＝Ｅ^２χ^０＋ｃ_１Ｅ^１＋ｃ_２Ψ^１とする。

式（２４）〜式（２８）より、累積変数Ψに関する以下の帰納的関係式が得られる。

式（２９）及び式（３０）では、ε_∞＝１、Δｔ＝１である。また、式（２９）及び式（３０）中、χ^０は式（２５）より、

で与えられる。更に、ｎ≦０であるすべてのｎについてΨの値は０とする。

なお、後述するように、ＦＤＴＤ法を実施する際には、式（２９）及び式（３０）に含まれておりω_ｐ及びν_ｃを用いて表されるｃ_１，ｃ_２は、ω_ｐ及びν_ｃに対応するＦＤＴＤ領域１００におけるプラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆの関数としてのｃ_１，Ｆ及びｃ_２，Ｆを採用する。

通常、２次元においては、任意の偏光状態にある入射ビーム（電磁波）２０は、２つの直交する偏光状態、すなわち、ＴＭモードとＴＥモードとに分けることが可能である。そこで、以下、ＴＭモードとＴＥモードとに分けてＦＤＴＤ法のアルゴリズムについて説明する。

[ＴＭモード]
ＴＭモードでは、式（１５）及び式（１６）を中心差分を利用して差分化して解く。

式（１６）中の時間に関する一次導関数の差分形式として式（３２）を利用する。

この場合、式（１６）に等価な差分形式の方程式は、以下のように表される。

式（３３）中のＤ_１は∇に対応する２次の中心差分形式の演算子である。また、式（３３）中のΨ^ｎ＋１−Ψ^ｎは、式（２９）により与えられる。

ここでは、式（１６）に等価な差分形式の方程式を示したが、式（１５）に関しては、従来のＦＤＴＤ法における差分形式である以下の式を用いる。

[ＴＥモード]
次にＴＥモードの場合について説明する。

ＴＥモードではマクスウェル方程式である式（１５）及び式（１６）から導かれる次の波動方程式を利用する。

本実施形態では∇・Ｅ＝０を仮定するが、高い伝導性を有する金属の場合、それらは自由電荷を長い間保持しないので、この仮定は妥当である。

式（１６）の時間に関する二次導関数の差分形式として式（３６）を利用する。

この場合、帰納的畳み込みを適用した式（３５）に等価な差分形式の方程式は、以下のように表される。

式（３７）中のＤ_１ ^２は∇^２に対応する２次の中心差分形式の演算子であり、ｖは位相速度である。また、式（３７）中のΨ^ｎ＋１＋Ψ^ｎ−１−２Ψ^ｎは式（３０）から与えられる。

本実施形態では、入射ビーム２０がＴＭモードで与えられる場合、すなわち、入射ビーム２０がＴＭ波である場合には、電磁界解析用の式（第２の電磁界解析用の式）である式（３３）及び式（３４）と、式（２９）とを利用して電界Ｅ及び磁界Ｈを算出する。入射ビーム２０がＴＥモードで与えられる場合、すなわち入射ビーム２０がＴＥ波である場合には、電磁界解析用の式（第１の電磁界解析用の式）である式（３７）と、式（３０）とを利用して電界Ｅを算出する。

次に、本実施形態におけるＦＤＴＤ領域における空間的な離散化条件の決定方法について説明する。

ＦＤＴＤ法では、図２に示したように、ＦＤＴＤ領域１００に複数のセルからなるグリッドを適用する。本実施形態のシミュレーション法で解くＦＤＴＤ法では、アルゴリズムの安定性に基づいてＦＤＴＤ領域１００の離散化条件を決定する。

式（２９）及び式（３０）より、アルゴリズムの安定性のためにｃ_１及びｃ_２は次の条件を満たすことが必要である。

ｃ_１，ｃ_２は式（２７）及び式（２８）に示すようにω_ｐ及びν_ｃを用いて表される。

一方、式（１９）より、誘電率ε（ω）の実部及び虚部をε^ｒ（ω）及びε^ｉ（ω）とすると、

が成り立つ。

よって、ω_ｐ及びν_ｃは次のように表される。

本実施形態では、入射ビーム２０が有する所定波長λに対応する角周波数をω_λとしたとき、式（４０）及び式（４１）内のε^ｒ（ω）及びε^ｉ（ω）として、角周波数ω_λに対する実測値を採用する。

所定波長λに対して誘電率の実測値が得られていない場合には、ドルーデ分散を利用して実測されている誘電率の局所的なフィッティングから上記所定波長λに対する誘電率、すなわち、ε^ｒ（ω）及びε^ｉ（ω）を得る。上記局所フィッティングは、所定波長λに最も近く誘電率が実測されている波長に対する誘電率の実測値を通るようにドルーデのε−ωプロットを実施することで得ることができる。

例えば、波長５００ｎｍ（角周波数３．７７×１０^１５／ｓｅｃ）での誘電率の実測値が得られていない一方、波長４９６ｎｍ（又は角周波数３．８×１０^１５／ｓｅｃ）と、波長５５１．１１１ｎｍ（角周波数３．４２×１０^１５／ｓｅｃ）に対してそれぞれ誘電率が得られているとする。この場合、波長５００ｎｍでの誘電率は、角周波数３．８×１０_１５／ｓｅｃに対する誘電率ε＝−９．５６４＋ｉ０．３０９をとおるドルーデのε−ωプロットから計算する（図３参照）。

より正確な誘電率を得るために、所定波長λに対して高波長側及び短波長側においてそれぞれ所定波長λに最も近く誘電率が実測されている波長に対する誘電率を基に前述したドルーデのε−ωプロットを実施し、所定波長λの誘電率としてそれぞれ得られる２つの推定値の平均値を算出して推定値とすることが好ましい。金属回折格子に入射させる光として例えばＬＥＤから出力された光等を利用される際には、多くの場合、ＬＥＤから出力される光の波長に対して誘電率が実測されていない場合があるため、上記のような算出方法は好ましい。なお、以下の説明では、実測値という場合には、上記ドルーデのε−ωプロットを実施して決まる値を含めるものとする。

また、本実施形態では、式（４２）及び式（４３）における角周波数ωをＦＤＴＤ領域における角周波数ω_Ｆに置き換えて、プラズマ角周波数ω_ｐ及び衝突周波数ν_ｃをＦＤＴＤ領域１００におけるプラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆに対応させる。この場合、プラズマ角周波数ω_ｐ，Ｆ及び衝突周波数ν_ｃ，Ｆは以下のように表される。

式（４４）及び式（４５）中のε^ｒ及びε^ｉは角周波数ω_λに対する実測値であり定数である。

そして、本実施形態では、実空間での角周波数ω_λと、式（４４）及び式（４５）中のＦＤＴＤ領域１００での角周波数（解析用角周波数）ω_Ｆとを次式によって対応づけることによって、ω_ｐ，Ｆ及びν_ｃ，Ｆを自由パラメータとしている。

式（４６）中のｖは入射ビーム２０の実空間での位相速度である。そして、λ_Ｆは、ＦＤＴＤ領域１００において波長λを表すのに必要とする格子点１１０の数であり、ｖ_ＦはＦＤＴＤ領域１００における入射ビーム２０の位相速度（解析用速度）である。

上記のように実空間におけるω_ｐ及びν_ｃに対応してＦＤＴＤ領域１００におけるω_ｐ，Ｆ及びν_ｃ，Ｆを規定した場合、式（３８），式（３９）に示す安定化条件は、

と表されることになる。

なお、ｃ_１，Ｆ及びｃ_２，Ｆは、ｃ_１，ｃ_２におけるω_ｐ及びν_ｃをω_ｐ，Ｆ及びν_ｃ，Ｆで置き換えたものであり、以下の式で表される。

また、アルゴリズムの安定性として、式（４９）及び式（５０）と共に、ＦＤＴＤ法で通常知られている以下の関係式を満たす必要がある。

前述したようにε^ｒ及びε^ｉとして角周波数ω_λに対する所定の値を利用すると共に、式（４６）を採用することによって、ω_Ｐ，Ｆ及びν_ｃ，Ｆを自由パラメータとした場合、式（４９）及び式（５０）並びに式（５１）を満たすようにλ_Ｆ及びｖ_Ｆを決定する。λ_Ｆは、所定波長λを表すための格子点１１０の数であるため、空間的離散化条件が決定されることになる。

次に、上記アルゴリズムを実施するための解析装置、シミュレーション方法及びシミュレーションプログラムについて説明する。

図４は、本発明に係るシミュレーション方法の一例を実施するための解析装置のブロック図である。

図４に示した解析装置３０は、例えばパーソナルコンピュータ装置やワークステーション等のコンピュータ装置であり、主な構成要素として、入力手段３１、メモリ３２、中央処理装置３３及び出力手段３４を備えている。

入力手段３１は、キーボードなどであり、シミュレーションを実施するためのパラメータを含む入力データを受け付ける。上記パラメータとしては、入射ビーム２０の波長λとそれに対応する角周波数ω_λを含む入射ビーム２０に関する諸条件、角周波数ω_λに対するε^ｒ及びε^ｉの実測値、金属細線１１の条件及び境界条件、及び、シミュレーションでの最大時刻ｔ_ｍａｘ（ｔ_ｍａｘ＝ｎΔｔ）等である。金属細線１１の諸条件としては、金属細線１１の幅ｗ及び厚さｈ並びに周期Λである。また、入射ビーム２０に関する条件としては上記の他に、入射ビーム２０の偏光状態（例えばＴＥモード又はＴＭモード等）、入射ビーム２０の金属回折格子１０への入射角θ_ｉｎ（金属回折格子法線Ｎと入射ビーム２０との為す角）等である。

メモリ３２は、解析装置３０を動作させるための種々のプログラムや、入力手段３１からの入力データや、後述するシミュレーション方法に沿ってＦＤＴＤ解析を実施するためのＦＤＴＤプログラム（シミュレーションプログラム）４０及びＦＤＴＤプログラム４０を実施中又は実施後の計算結果を記憶する。ＦＤＴＤプログラム４０は、式（４４）〜式（５１）に基づいて空間離散化条件、すなわち、λ_Ｆを決定せしめる離散化条件決定モジュール４１と、ＴＭモード及びＴＥモードに分けて説明したアルゴリズムに基づいてＦＤＴＤ法を実行するための解析モジュール４２とを含んで構成されている。

中央処理装置３３は、メモリ３２に記憶されているＦＤＴＤプログラム４０を読み込み、シミュレーションのための種々の演算を実行する。また、中央処理装置３３は、メモリ３２に記憶されている各種プログラムを読み込み、解析装置３０の構成要素を制御する機能も有する。

出力手段３４は、例えばディスプレイなどであり、中央処理装置３３で算出されて得られるシミュレーション結果を表示する。

次に、ＦＤＴＤ法に従ったシミュレーション方法について説明する。図５は、シミュレーション方法を示すフローチャートである。シミュレーションは、解析装置３０のメモリ３２に格納されているＦＤＴＤプログラム４０を中央処理装置３３が読み込み実行することで実施される。

図５に示すように、ステップＳ１０で操作者が入力手段３１を介してシミュレーション用のパラメータを入力する（入力工程）。入力手段３１がパラメータを受け付けると、メモリ３２がそのパラメータを記憶する。

次に、ステップＳ２０において、離散化条件を決定する（離散化条件決定工程）。この離散化条件決定工程では、中央処理装置３３が、ステップＳ１０で入力されたε^ｒ及びε^ｉの値並びに式（４４）〜式（５１）を利用して、式（４７），式（４８）及び式（５１）を満たすλ_Ｆ及びｖ_Ｆを決定する。この際、式（４７），式（４８）及び式（５１）を満たすλ_Ｆ及びν_Ｆが複数ある場合には、例えば、中央処理装置３３がそれらを比較してλ_Ｆが最小となる場合を選択するようにすることが好ましい。これにより、少ない格子点１１０の数で安定したシミュレーションが可能となるからである。

続くステップＳ３０では、ステップＳ２０で決定された空間的な離散化条件で規定されるＦＤＴＤ領域１００に、ステップＳ１０で入力された諸条件を満たす金属回折格子１０を設定し、ＦＤＴＤ領域１００における電磁界を時間ステップ間隔Δｔ毎に算出する（解析工程）。この解析工程（ステップＳ３０）においては、ステップＳ１０で入力された境界条件を必要に応じて使用する。本実施形態では、Δｔ＝１である。

ＴＥモードとＴＭモードとで使用するアルゴリズムが異なるため、ステップＳ３０では、(i)入射ビーム２０がＴＥモードとして与えられている場合、(ii)入射ビーム２０がＴＭモードとして与えられている場合、に分けて解析を実施する。

図６は、図５の解析処理のフローチャートである。解析処理では、先ず、ステップＳ３１において、入射ビーム２０の偏光状態を判定する。以下、ステップＳ１０において入射ビーム２０がＴＥモードで与えられている場合及び入射ビーム２０がＴＭモードで与えられている場合について説明する。

（ｉ）入射ビーム２０がＴＥモードで与えられている場合、ステップＳ３２Ａにおいて、セルの電磁界を初期化する。具体的には、各セルにおける電界Ｅ及び磁界Ｈの各軸方向の成分を０とする。

続いてステップＳ３３Ａに進み、式（３７）及び式（３０）を利用して電界Ｅを算出し、電界Ｅをアップデートすると共に、時刻ｔをアップデートする。なお、式（３０）におけるｃ_１，ｃ_２の値としては、ステップＳ２０で決定されたλ_Ｆ及びｖ_Ｆにより規定されるｃ_１，Ｆ及びｃ_２，Ｆの値を採用すればよい。

ステップＳ３３Ａで全てのセルにおける電界Ｅを算出すると、ステップＳ３４Ａに進み、アップデートされた時刻ｔが最大時刻ｔ_ｍａｘであるか否かを判定する。時刻ｔが最大時刻ｔ_ｍａｘより小さい場合（ステップＳ３４ＡでＹ）には、ステップＳ３３Ａに戻る。

次に、（ｉｉ）入射ビーム２０がＴＭモードで与えられている場合について説明する。この場合、ステップＳ３２Ｂにおいて、セルを初期化する。具体的には、各セルにおける電界Ｅ及び磁界Ｈの各軸方向の成分を０とする。次いで、ステップＳ３３Ｂに進み、式（３３）及び式（３４）並びに式（２９）を利用して各セルの電界Ｅ及び磁界Ｈを算出し、電界Ｅ及び磁界Ｈをアップデートすると共に、時刻ｔをアップデートする。式（２９）におけるｃ_１，ｃ_２の値に対しては、ＴＥモードの場合と同様である。

ステップＳ３３Ｂで全てのセルにおける電界Ｅ及び磁界Ｈを算出すると、ステップＳ３４Ｂに進み、アップデートされた時刻ｔが最大時刻ｔ_ｍａｘであるか否かを判定する。時刻ｔが最大時刻ｔ_ｍａｘより小さい場合（ステップＳ３４Ｂで「Ｙ」の場合）には、ステップＳ３３Ｂに戻る。なお、図６中では、ＴＥモード及びＴＭモードに分けて説明するために、時刻ｔと最大時刻ｔ_ｍａｘとの大小関係の判定ステップを２つに分けて説明しているが、ステップＳ３４Ａ及びステップＳ３４Ｂは実質的に同じステップである。

ステップＳ３４Ａ又はステップＳ３４Ｂにおいて、時刻ｔが最大時刻t_ｍａｘ以上である場合（ステップＳ３４Ａ，Ｓ３４Ｂで「Ｎ」の場合）には、図５に記載のステップＳ４０に進み、解析結果を出力手段３４により出力する（出力工程）。

また、本実施形態におけるＦＤＴＤ法の実施では、上述した各式を使用する点以外は、従来のＦＤＴＤ法のシミュレーション方法と同様であることから詳細な説明は省略している。また、金属細線１１以外の部分のセルに関しては、通常の空気中又は真空中を伝搬する電磁波に対するマクスウェル方程式又は波動方程式を扱う場合のＦＤＴＤ法と同様に、ＴＥモードに関しては下記式（５２）を利用し、ＴＭモードに関しては式（５３）及び下記式（５４）を利用すればよい。

式（５２）及び式（５３）中において、ａ＝(σ_ｄΔｔ)／（２ε_ｄ）であり、σ_ｄ及びε_ｄは隣接する金属細線１１間の媒質の電導度及び誘電率である。また、μは隣接する金属細線１１間の媒質の透磁率である。

以上説明した本実施形態のシミュレーション法及びＦＤＴＤプログラムでは、解析対象としての金属回折格子１０が有する金属細線１１の誘電率に対してドルーデの分散式を適用してアルゴリズムを定式化している。また、金属細線１１の誘電率に対して実測値を採用することでＦＤＴＤ領域１００の空間的な離散化を図っている。そして、空間離散化を図ることでアルゴリズムの安定性も保証している。その結果、金属回折格子１０に入射ビーム２０が入射した場合の電磁界を、より正確に安定してシミュレーションすることが可能である。

また、ＴＥモードに対するシミュレーションでは、マクスウェル方程式から導かれる波動方程式を利用していることから、シミュレーションのアルゴリズムを簡略化することができている。その結果、シミュレーションに要する時間の短縮化も図られている。更に、本実施形態のシミュレーション方法及びＦＤＴＤプログラムでは、ＲＣ法を採用しているため、メモリ３２に蓄えるべきデータ量を低減することが可能となっている。

次に、本実施形態で説明したシミュレーション方法及びＦＤＴＤプログラムを用いた金属回折格子１０のシミュレーション結果について説明する。ここでは、実施例１及び実施例２として２つの金属回折格子１０についてそれぞれシミュレーションを実施した。

実施例１のシミュレーションの条件は以下のとおりである。
・金属細線１１の材料：銀
・金属回折格子１０の周期Λ：１７６．１２５ｎｍ
・金属回折格子１０のフィルファクタ（ｗ／Λ）：０．５
・金属細線１１の厚さｈ：１００ｎｍ
・入射ビーム２０の波長λ：７０４．５ｎｍ（角周波数ω_λ：２．６７６×１０^１５／ｓｅｃ）
・入射ビーム２０の位相速度ｖ：３×１０^８ｍ／ｓｅｃ
・角周波数ω_λに対するε^ｒの値：−２３．４０５
・角周波数ω_λに対するε^ｉの値：０．３８７

実施例２のシミュレーションの条件は以下のとおりである。
・金属細線１１の材料：銀
・金属回折格子１０の周期Λ：４９６ｎｍ
・金属回折格子１０のフィルファクタ（ｗ／Λ）：０．５
・金属細線１１の厚さｈ：４９６ｎｍ
・入射ビーム２０の波長λ：４９６ｎｍ（角周波数ω_λ：３．８×１０^１５／ｓｅｃ）
入射ビーム２０の位相速度ｖ：３×１０^８ｍ／ｓｅｃ
・角周波数ω_λに対するε^ｒの値：−９．５６４
・角周波数ω_λに対するε^ｉの値：０．３０９

実施例１，２では、入射ビーム２０をＴＥモード及びＴＭモードとしてそれぞれ与えた場合に、入射角θ_ｉｎを変えながら図５，６に示したフローチャートに基づいてシミュレーションを実施した。なお、実施例１，２では、金属回折格子１０の大きさはｘ軸方向に無限であることを仮定している。

実施例１の条件の場合、ステップＳ２０において、λが７０４．５ｎｍであり、ε^ｒが−２３．４０５であり、ε^ｉが０．３８７であるとき、λ_Ｆは３５．２２７であり、ν_Ｆは０．５３８であった。この場合、ｃ_１，Ｆは０．２２５であり、ｃ_２，Ｆは１．００２である。また、実施例２の条件の場合、ステップＳ２０において、λが４９６ｎｍであり、ε^ｒが−９．５６４であり、ε^ｉが０．３０９であるとき、λ_Ｆは３５．２２７であり、ν_Ｆは０．５３８であった。この場合、ｃ_１，Ｆは０．０９８であり、ｃ_２，Ｆは１．００３である。

図７は、実施例１の場合のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示している。図７の横軸は入射角θ_ｉｎ（°）を示し、縦軸は透過率（％）及び反射率（％）を示している。透過率及び反射率は、入射ビーム２０の光強度に対する透過ビーム及び反射ビームの光強度の比として定義している。図８は実施例２のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示している。図８の横軸及び縦軸は図７の場合と同様である。なお、図７及び図８のシミュレーション結果は、銀から構成される無限の大きさの金属回折格子１０におけるファーフィールドでの透過率及び反射率に基づいている。ファーフィルでの結果は、ＦＤＴＤ法を実行して得られる複素電磁界を金属回折格子１０の周囲を囲む閉曲線に沿って統合することによって計算される。

図７及び図８に示すように、上記シミュレーションにより、ＴＥモード及びＴＭモードに対して透過特性及び反射特性を得ることができる。そして、上記シミュレーションにより、金属回折格子１０の構成の違いによる特性の違いが現れていることがわかる。例えば、図７に示すように、上記シミュレーションにより、実施例１の金属回折格子１０では、ＴＥモードの多くは反射しており、透過が抑制されている。また、実施例１の金属回折格子１０では、ＴＭモードの多くが透過させることができている一方、その一部を反射させることが可能となっている。このように金属回折格子１０の構成では偏光分離ができていることがわかる。また、図８に示すように、実施例２の金属回折格子１０では、ＴＭモード及びＴＥモードが共に透過可能であるが、ＴＥモードの方がより多く透過させることができており、ＴＭモードの透過率を入射角θ_ｉｎが０〜６０°で１０％〜４０％に抑制できていることがわかる。

よって、上記シミュレーション方法及びＦＤＴＤプログラムを利用すれば、金属回折格子１０に関するパラメータを変更することで、所望の透過特性及び反射特性を備えた金属回折格子１０を設計可能であることになる。

図９は、図１に示した金属回折格子が適用された導光板ユニットの構成を模式的に示す側面図である。導光板ユニット５０は、略直方体形状の導光板５１と導光板５１の表面５１ａ上に金属回折格子１０が設けられて構成されている。導光板ユニット５０は、導光板５１の側面５１ｂの側方に光源（不図示）が配置されることで面光源装置として利用することができる。

光源から出力され側面５１ｂから導光板５１内に入射された光（ビーム）は、導光板５１内を側面５１ｃ側に向けて伝搬する。この導光板５１内を伝搬する光が金属回折格子１０に入射すると、金属回折格子１０の回折により、入射された光の一部は導光板ユニット５０から出射され、他の部分は主に導光板５１内に戻される。

このように金属回折格子１０を導光板ユニット５０の一部として利用する場合には、金属回折格子１０が偏光分離素子として設計されていることが好ましい。具体的には、所定の偏光成分（例えば、Ｐ偏光成分としてのＴＭモード）が選択的に透過され、他の光が反射されて導光板５１内に戻されるように、金属回折格子１０が設計されていることが好ましい。この場合、導光板ユニット５０からは所定の偏光成分が支配的なビーム（透過ビーム２１）が出射されることになるので、導光板ユニット５０からの出射ビームを、例えば液晶表示装置におけるバックライトとして好適に使用することができる。そして、このようなバックライトを使用する液晶表示装置では、不要偏光を除去する光学素子が不要となるため、小型化や薄型化を図ることが可能である。

そして、上述したような特性を有する金属回折格子１０を設計するために、前述したシミュレーション方法及びＦＤＴＤプログラムが有効である。金属回折格子１０を、例えば液晶表示装置が有する面光源装置で使用する場合には、所定波長λは、可視領域の波長として、シミュレーションを実施すればよい。

金属回折格子１０を導光板５１上に設ける際には、図１及び図２に示したように、金属回折格子１０のみを解析対象としてシミュレーションを実施してもよいが、例えば、図９に示したように、板状の導光板５１を含む導光板ユニット５０を解析対象としても良い。

この場合、導光板５１に対応するセルでは、ＦＤＴＤ法のアルゴリズムとして以下の式を利用すればよい。すなわち、ＴＥモードの場合には式（５５）を利用し、ＴＭモードの場合には、式（５６）及び式（５７）を利用すればよい。

式（５５）及び式（５６）中においては、ａ＝(σ_ｄΔｔ)／（２ε_ｄ）であり、σ_ｄ及びε_ｄは導光板５１の電導度及び誘電率である。また、μは導光板５１の透磁率である。

また、上記実施形態では、入射ビーム２０は、ＴＥモード又はＴＭモードとして与えられるとしたが、例えば、ＴＥモード及びＴＭモードを含む非偏光状態の入射ビーム２０を与えても良い。この場合には、解析工程（ステップＳ３０）において、入射ビーム２０中のＴＥモードに対しては、ステップＳ３３Ａで説明した場合と同様にして電界Ｅを算出し、入射ビーム２０中のＴＭモードに対しては、ステップＳ３３Ｂで説明した場合と同様にして電界Ｅ及び磁界Ｈを算出し、それらを組み合わせて各セルの電界Ｅ及び磁界Ｈとすればよい。

ＴＭモード及びＴＥモードに対して算出した電界Ｅ及び磁界Ｈの組み合わせ方としては、各セルに対して算出した電界Ｅ及び磁界Ｈそれぞれの加重平均をとることができる。例えば、一つのセルに対応して、ＴＭモードに対して算出した電界をＥ_ＴＭとし、ＴＥモードに対して算出した電界をＥ_ＴＥとしたとき、そのセルにおける電界を（Ｅ_ＴＭ＋Ｅ_ＴＥ）／２とすることができる。磁界についても同様である。なお、このように加重平均をとる場合には、例えば、入射ビーム２０をＴＥモード及びＴＭモードからなる光であると仮定して算出した電界Ｅ及ぶ磁界Ｈについてそれぞれ前述のように加重平均を算出することで、非偏光状態の入射ビーム２０が入射した場合の電界Ｅ及び磁界Ｈとすることもできる。

また、金属回折格子１０は、図９を利用して説明したように導光板ユニット５０に利用することができるが、例えば、光学系内の所定の位置に配置される偏光分離素子として利用することもできる。この場合、例えば、ＴＥモード及びＴＭモードを完全分離するように金属回折格子１０を設計することが考えられる。

また、本発明に係るシミュレーションプログラムは、例えば、ＴＥモード及びＴＭモードの何れにも対応可能なプログラムとして説明したが、ＴＥモードに対するプログラム及びＴＭモードに対するプログラムすることもできる。そして、シミュレーション方法における解析工程（ステップＳ３０）では、図６に示したステップＳ３２における判定工程を省略して、図６に示したＴＥモード及びＴＭモードの各工程に沿ってシミュレーションを実施することも可能である。更に、上記実施形態では、２次元のＦＤＴＤ法について説明したが、３次元のＦＤＴＤ法に拡張することもできる。更に、上記実施形態では金属部材を含む構造物として金属回折格子１０を例示して説明したが、上記構造物は金属回折格子に限定されず、本発明に係るシミュレーション方法及びシミュレーションプログラムは、金属部材を含む任意の構造物に対して適用することができる。

本発明に係るシミュレーション方法の一実施形態適用して設計するための光学素子の模式図である。 ＦＤＴＤ領域に設定された図１に示した光学素子の模式図である。 角周波数に対する誘電率の実測値及びフィッティングカーブの一例を示す図面である。 本発明に係るシミュレーション方法の一実施形態を適用するための解析装置のブロック図である。 本発明に係るシミュレーション方法の一実施形態のフローチャートである。 図５の解析処理のフローチャートである。 実施例１のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示す図面である。 実施例２のシミュレーション結果に基づいた透過スペクトル及び反射スペクトルを示図面である。 図１に示した光学素子を含む導光板ユニットの構成を模式的に示す側面図である。

符号の説明

１０…金属回折格子（構造物）、１１…金属細線（金属部材）、２０…入射ビーム（電磁波）、１００…ＦＤＴＤ領域、１１０…格子点。

Claims (10)

1. 金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーション方法であって、
   前記構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、前記金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、
   前記離散化条件決定工程で決定された前記空間離散化条件に基づいて、前記解析モデルに前記電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いた前記ＦＤＴＤ法に従って前記コンピュータにより算出する解析工程と、
   を備え、
   前記離散化条件決定工程では、前記帰納的畳み込み法において前記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ_１，Ｆ，ｃ_２，Ｆであって、前記ＦＤＴＤ領域において前記電磁波が有する所定波長を表すための格子点の数及び前記ＦＤＴＤ領域における前記電磁波の解析用速度並びに前記金属部材の前記誘電率を含むプラズマ角周波数ω_ｐ、Ｆ及び衝突周波数ν_ｃ，Ｆを用いて下記式（１），（２）で表される前記変数ｃ_１，Ｆ，ｃ_２，Ｆ：
   

   

   
   

   

   
   において、
   前記プラズマ角周波数ω_ｐ，Ｆ及び前記衝突周波数ν_ｃ，Ｆに含まれる前記誘電率を、前記所定波長の前記電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、前記変数ｃ_１，Ｆ，ｃ_２，Ｆが、前記ＦＤＴＤ法における安定化条件に含まれる下記式（３）及び式（４）を満たすと共に、
   

   

   
   

   

   
   前記解析用速度が０．７より小さくなるように、前記格子点の数及び前記解析用速度を、前記コンピュータによって決定することを特徴とするシミュレーション方法。
2. 実空間における前記所定波長をλ、前記電磁波の速度をｖ及び前記所定波長に対応する角周波数をω_λとし、前記ＦＤＴＤ領域における前記格子点の数をλ_Ｆ、前記解析用速度をｖ_Ｆとし、前記角周波数ω_λに対する前記誘電率の実部及び虚部の値をε^ｒ及びε^ｉとし、前記角周波数ω_λに対応する前記ＦＤＴＤ領域における解析用角周波数ω_Ｆを、
   

   

   
   としたとき、前記衝突周波数ν_ｃ，Ｆ及び前記プラズマ角周波数ω_ｐ、Ｆは、
   

   

   
   

   

   
   として与えられることを特徴とする請求項１に記載のシミュレーション方法。
3. 前記解析工程では、
   前記電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータにより解くことを特徴とする請求項１又は２に記載のシミュレーション方法。
4. 前記解析工程では、
   前記電磁波におけるＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータにより解くことを特徴とする請求項１〜３の何れか一項に記載のシミュレーション方法。
5. 前記解析工程では、
   前記電磁波がＴＭモード及びＴＥモードを共に含む場合には、前記ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して解き、前記ＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して解き、前記ＴＭモード及び前記ＴＥモードに対して得られた電磁界に基づいて前記ＦＤＴＤ領域における電磁界を算出する、ことを特徴とする請求項１又は２に記載のシミュレーション方法。
6. 金属部材を含む構造物に電磁波が入射した場合の電磁界を、ＦＤＴＤ法に従って解析するためのシミュレーションをコンピュータに実行させるためのシミュレーションプログラムであって、
   前記構造物の解析モデルが配置されるＦＤＴＤ領域の空間離散化条件を、前記金属部材の誘電率を基にコンピュータにより決定する離散化条件決定工程と、
   前記離散化条件決定工程で決定された前記空間離散化条件に基づいて、前記解析モデルに前記電磁波が入射した場合の電磁界を、帰納的畳み込み法を用いた前記ＦＤＴＤ法に従って前記コンピュータにより算出する解析工程と、
   を前記コンピュータに実行させ、
   前記離散化条件決定工程では、前記帰納的畳み込み法において前記誘電率をドルーデの分散式で表すことで得られる帰納的関係式内の２つの変数ｃ_１，Ｆ，ｃ_２，Ｆであって、前記ＦＤＴＤ領域において前記電磁波の所定波長を表すための格子点の数及び前記ＦＤＴＤ領域における前記電磁波の解析用速度並びに前記金属部材の前記誘電率を含むプラズマ角周波数ω_ｐ及び衝突周波数ν_{ｃ, F}を用いて下記式（８），（９）で表される前記変数ｃ_１，Ｆ，ｃ_２，Ｆ：
   

   

   
   

   

   
   において、
   前記プラズマ角周波数ω_ｐ，Ｆ及び前記衝突周波数ν_ｃ，Ｆに含まれる前記誘電率を、前記所定波長の前記電磁波に対する実測値又は実測値に基づいて決まる値とした場合に、前記変数ｃ_１，Ｆ，ｃ_２，Ｆが、前記ＦＤＴＤ法における安定化条件に含まれる下記式（１０）及び式（１１）を満たすと共に、
   

   

   
   

   

   
   前記解析用速度が０．７より小さくなるように、前記格子点の数及び前記解析用速度を、前記コンピュータに決定せしめることを特徴とするシミュレーションプログラム。
7. 実空間における前記電磁波の前記所定波長をλ、前記電磁波の速度をｖ及び前記所定波長に対応する角周波数をω_λとし、前記ＦＤＴＤ領域における前記格子点の数をλ_Ｆ、前記解析用速度をｖ_Ｆとし、前記角周波数ω_λに対する前記誘電率の実部及び虚部の値をε^ｒ及びε^ｉとし、前記角周波数ω_λに対応する前記ＦＤＴＤ領域における解析用角周波数ω_Ｆを、
   

   

   
   としたとき、前記衝突周波数ν_ｃ，Ｆ及び前記プラズマ角周波数ω_ｐ，Ｆは、
   

   

   
   

   

   
   として与えられることを特徴とする請求項６に記載のシミュレーションプログラム。
8. 前記解析工程では、
   前記電磁波におけるＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめることを特徴とする請求項６又は７に記載のシミュレーションプログラム。
9. 前記解析工程では、
   前記電磁波におけるＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめることを特徴とする請求項６〜８の何れか一項に記載のシミュレーションプログラム。
10. 前記解析工程では、
    前記電磁波がＴＭモード及びＴＥモードを共に含む場合には、前記ＴＥモードに対しては、マクスウェル方程式から導かれる波動方程式に前記帰納的畳み込み法を適用して得られる第１の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめて、前記ＴＭモードに対しては、マクスウェル方程式に前記帰納的畳み込み法を適用して得られる第２の電磁界解析用の式を、前記帰納的関係式を利用して前記コンピュータに解かさしめて、前記ＴＭモード及び前記ＴＥモードに対して得られた電磁界に基づいて前記ＦＤＴＤ領域における電磁界を算出する、ことを特徴とする請求項６又は７に記載のシミュレーションプログラム。

Images