JP3818874B2 - 電磁波解析装置および電磁波解析プログラム - Google Patents
電磁波解析装置および電磁波解析プログラム Download PDFInfo
- Publication number
- JP3818874B2 JP3818874B2 JP2001193224A JP2001193224A JP3818874B2 JP 3818874 B2 JP3818874 B2 JP 3818874B2 JP 2001193224 A JP2001193224 A JP 2001193224A JP 2001193224 A JP2001193224 A JP 2001193224A JP 3818874 B2 JP3818874 B2 JP 3818874B2
- Authority
- JP
- Japan
- Prior art keywords
- calculation
- calculated
- simulation time
- electromagnetic field
- field
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Lifetime
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明はマックスウェル方程式を時間と空間について差分法で解く電磁波解析装置および電磁波解析プログラムに関し、特に空間を不均一なサイズのセルに分割してマックスウェル方程式の解を求める電磁波解析装置および電磁波解析プログラムに関する。
【０００２】
【従来の技術】
電磁波の過渡的な挙動を、電子計算機を使った数値計算によって解析する方法の一つに時間領域差分法（Finite-Difference Time-Domain Method、以下ＦＤＴＤ法という）がある。この手法は、マックスウェル方程式を時間と空間について差分法で解く方法であり、適用範囲が広いこと、計算精度や計算効率が良好であることなどから、今日幅広く利用されている。ＦＤＴＤ法については、「K.S.Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media,"IEEE Trans. AP-14, pp.302-307, 1966.」や「A.Taflove, Computational Electrodynamics, MA, Artech House, 1995.」に詳しい。
【０００３】
ＦＤＴＤ法は、マックスウェル方程式の２つの回転の式を、時間ｔと空間（ｘ−ｙ−ｚ空間）について差分法で解く方法である。
以下に、マックスウェル方程式について説明する。なお、以下の式において、Ｅｘは電界のｘ軸方向の成分、Ｅｙは電界のｙ軸方向の成分、Ｅｚは電界のｚ軸方向の成分、Ｈｘは磁界のｘ軸方向の成分、Ｈｙは磁界のｙ軸方向の成分、Ｈｚは磁界のｚ軸方向の成分、μは透磁率、εは誘電率、σは電気伝導率である。
【０００４】
１次元（Ｘ軸方向に伝搬する電磁波）のマックスウェル方程式は、次の式で表される。
【０００５】
【数１】
【０００６】
【数２】
【０００７】
２次元（Ｘ軸、Ｙ軸方向に伝搬する電磁波）のマックスウェル方程式は、ＴＭ波(Transverse Magnetic wave)については次の式で表される。
【０００８】
【数３】
【０００９】
【数４】
【００１０】
【数５】
【００１１】
２次元（Ｘ軸、Ｙ軸方向に伝搬する電磁波）のマックスウェル方程式は、ＴＥ波(Transverse Electric wave)については次の式で表される。
【００１２】
【数６】
【００１３】
【数７】
【００１４】
【数８】
【００１５】
３次元のマックスウェル方程式は、次の式で表される。
【００１６】
【数９】
【００１７】
【数１０】
【００１８】
【数１１】
【００１９】
【数１２】
【００２０】
【数１３】
【００２１】
【数１４】
【００２２】
ＦＤＴＤ法では、通常、解析領域を直方体のセル（差分格子）に分割することにより、空間的な離散化を行う。このとき、セルの大きさが小さいほど、電磁波の過渡的な挙動が急激に変わる領域においても、少ない誤差で計算することができる。
【００２３】
ところが、空間を均一なセルで一様に分割した場合、しばしば多大なセル数を必要とし、計算時間や計算機メモリの増加を引き起こす。そこで、大きさの異なるセルを組み合わせて用いるいわゆる不均一なセルで解析領域を分割することにより、セル数を節約することがたびたび行われる。
【００２４】
【発明が解決しようとする課題】
しかし、解析領域をサイズが不均一なセルで分割した場合、セルサイズが相異なる領域の境界（以下、セルサイズ境界という）において、誤差が増加することが知られている。
【００２５】
すなわち、セルサイズ境界では、隣接するセルに基づいた差分近似の中心位置とセルサイズ境界の位置とがずれてしまう。この位置のずれに基づく差分近似の打ち切り誤差が、電磁波解析における誤差の増加を招いている。
【００２６】
セルサイズが同一の領域では、差分近似の打ち切り誤差は２次のオーダであるが、セルサイズ境界における差分近似の打ち切り誤差は１次のオーダである。これは、解析対象よっては、計算精度を著しく低下させる要因となる。
【００２７】
本発明はこのような点に鑑みてなされたものであり、不均一セルを用いた電磁波解析の誤差を低減させた電磁波解析装置を提供することを目的とする。
【００２８】
【課題を解決するための手段】
本発明では上記課題を解決するために、図１に示すような電磁波解析装置が提供される。本発明に係る電磁波解析装置は、解析領域１を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するものである。
【００２９】
本発明の電磁波解析装置は、セルサイズ判断手段２、第１の計算手段３、第２の計算手段４および出力手段５を有している。
セルサイズ判断手段２は、解析領域１に設定された各計算位置Ｐ１，Ｐ２に関し、周囲のセルのサイズの同一性を判断する。第１の計算手段３は、セルサイズ判断手段２において周囲のセルのサイズが同じであると判断された計算位置Ｐ１における電磁界を、第１の計算手法で計算する。第２の計算手段４は、セルサイズ判断手段２において周囲のセルのサイズが異なると判断された計算位置Ｐ２における電磁界を、第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算する。出力手段５は、第１の計算手段と第２の計算手段とにより計算された各計算位置における電磁界の値を出力する。
【００３０】
このような電磁波解析装置によれば、サイズの同じセルサイズが均一な領域内の計算位置においては、セルサイズ判断手段２によって、周囲のセルのサイズが同一であると判断される。すると、その計算位置の電磁界は、第１の計算手段３により、第１の計算手法で計算される。一方、異なるサイズのセル同士の隣接する計算位置では、セルサイズ判断手段２によって、周囲のセルのサイズが異なると判断される。すると、その計算位置の電磁界は、第２の計算手段４により、第１の計算手法よりも誤差の少ない第２の計算手法で計算される。第１の計算手段３および第２の計算手段４で計算された電磁界の値は、出力手段５によって出力される。
【００３１】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して説明する。
図１は、本発明の原理構成図である。本発明に係る電磁波解析装置は、電磁波の挙動を数値計算する場合に、マックスウェル方程式を時間と空間について差分法で解くものである。差分法では、解析領域１を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するものである。
【００３２】
本発明の電磁波解析装置は、セルサイズ判断手段２、第１の計算手段３、第２の計算手段４および出力手段５を有している。
セルサイズ判断手段２は、解析領域１に設定された各計算位置Ｐ１，Ｐ２に関し、周囲のセルのサイズの同一性を判断する。
【００３３】
たとえば、セルサイズ判断手段２は、解析領域１に設定された座標軸毎に、隣接するセルの各軸方向の幅の同一性を判断する。すなわち、セルサイズ判断手段２は、ある軸方向に隣接するセル同士のその軸方向へのセルの幅が同じであれば、その軸方向への電磁波解析においては、周囲のセルのサイズが同一であるものと判断する。一方、セルサイズ判断手段２は、ある軸方向に隣接するセル同士のその軸方向へのセルの幅が異なれば、その軸方向への電磁波解析においては、周囲のセルのサイズが異なるものと判断する。
【００３４】
第１の計算手段３は、セルサイズ判断手段２において周囲のセルのサイズが同じであると判断された計算位置Ｐ１における電磁界を、第１の計算手法で計算する。たとえば、第１の計算手段３は、電界または磁界の偏微分係数を算出する際に、陰的スキーム（陰解法に基づいた計算手法）を用いず、既に算出されている磁界の値から算出する。
【００３５】
第２の計算手段４は、セルサイズ判断手段２において周囲のセルのサイズが異なると判断された計算位置Ｐ２における電磁界を、第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算する。
【００３６】
たとえば、第２の計算手段４は、周囲のセルのサイズが異なる計算位置の電磁界の偏微分係数を算出する際に、高次の陰的スキームを用いて算出する。このとき、第２の計算手段４は、周囲のセルのサイズが異なる計算位置の周囲の電磁界の微分係数を前記第１の計算手段３から取得し、第１の計算手段３で計算された電磁界の微分係数を用いて、陰的スキームの解を求める。
【００３７】
出力手段５は、第１の計算手段３と第２の計算手段４とにより計算された各計算位置における電磁界の値を出力する。
このような電磁波解析装置によれば、セルサイズが均一な領域内の計算位置においては、セルサイズ判断手段２によって、周囲のセルのサイズが同一であると判断される。すると、その計算位置の電磁界は、第１の計算手段３により、第１の計算手法で計算される。一方、異なるサイズのセル同士の隣接する計算位置では、セルサイズ判断手段２によって、周囲のセルのサイズが異なると判断される。すると、その計算位置の電磁界は、第２の計算手段４により、第１の計算手法よりも誤差の少ない第２の計算手法で計算される。第１の計算手段３および第２の計算手段４で計算された電磁界の値は、出力手段５によって出力される。
【００３８】
このように、本発明では、セルの大きさが相異なる領域が存在する場合、セルの大きさが均一の領域における電磁界の計算式と、異なる大きさのセルの境界における電磁界の計算式とが異なる。セルの中心を磁場の計算位置とし、セル同士の境界を電界の計算位置とした場合、たとえば、以下の計算式で電磁界を求め出すことができる。
【００３９】
ここでは、簡単のため１次元問題について考える。また、電気伝導率σは０とする。
電磁波解析シミュレーションにおける時刻Δｔ（時間刻み幅）×ｎ（ｎは０以上の整数）までの電磁界の過渡的変化が計算されているときに、時刻Δｔ×（ｎ＋１）までの電磁界を求めるときの計算式を以下に示す。
【００４０】
なお、時刻Δｔ×ｎを単にｎと表し、時刻Δｔ×（ｎ＋１）をｎ＋１と表し、時刻Δｔ×（ｎ＋１／２）をｎ＋１／２と表すものとする。各物理量の右上に、その物理量の計算時刻（シミュレーション上の時刻）を示す。
【００４１】
また、ｉは０以上の整数であり、ｉ番目のセルの始点位置を示している。ｉ＋Δｘ（ｘ軸方向の空間刻み幅：セルサイズ）は、ｉ＋１番目のセルの始点位置である。ｉ＋１／２は、ｉ番目のセルの中心位置である。各物理量の右下に、その物理量の計算位置を示す。
【００４２】
まず、時刻ｎ＋１／２における磁界は、
【００４３】
【数１５】
【００４４】
で算出される。式（１５）における電界の偏微分係数は、
【００４５】
【数１６】
【００４６】
で算出される。
また、時刻ｎ＋１の電界は、
【００４７】
【数１７】
【００４８】
で算出される。式（１７）における磁界の偏微分係数は、セルサイズが均一な領域においては、
【００４９】
【数１８】
【００５０】
で算出される。式（１８）の差分近似式の打ち切り誤差は２次のオーダである。打ち切り誤差の詳細は後述する。
また、セルサイズ境界における偏微分係数は、
【００５１】
【数１９】
【００５２】
で算出される。なお、Ａ、Ｂ、Ｃは、以下の通りである。
【００５３】
【数２０】
【００５４】
【数２１】
【００５５】
【数２２】
【００５６】
式（１９）の差分近似式の打ち切り誤差は３次のオーダである。打ち切り誤差の詳細は後述する。
なお、セルサイズ境界における偏微分係数を求める場合、始めにセルサイズ境界以外の場所における空間についての偏微分係数を式（１６）および式（１８）を用いて差分近似する。次にセルサイズ境界上における空間についての偏微分係数を式（１９）を用いて差分近似する。
【００５７】
なお、式（１５）から式（２２）に示した計算は、ある特定の軸方向（ここではｘ方向）に着目し、Ｅ(電界)がセルの端、Ｈ(磁界)がセルの中心に配置した場合に成立する。Ｅ(電界)がセルの中心、Ｈ(磁界)がセルの端に配置してある場合は、式（１５）から式（２２）の各式のＥとＨならびにεとμを入れ替えればよい。
【００５８】
このように、本発明によれば、セルの大きさが相異なる領域の境界上における空間についての偏微分係数は、式（１９）を用いて差分近似しているため、差分近似の打ち切り誤差は２次のオーダとなる。なお、それ以外の場所における空間についての偏微分係数を差分近似する場合には、従来と同様の式（１６）および式（１８）を用いるが、これは中心差分であることが保証されるので、差分近似の打ち切り誤差は常に２次のオーダとなる。
【００５９】
これにより、打ち切り誤差の増加を防ぐことができる。すなわち、不均一セルを用いた電磁波解析全体としての計算精度が向上する。
ここで、上記計算式の導出方法について説明する。
【００６０】
ＦＤＴＤ法による差分近似を行う場合、関数ｆ（ｘ）について、次のテーラー級数を考える。
【００６１】
【数２３】
【００６２】
【数２４】
【００６３】
式（２３）から式（２４）を減算し、ｆ’について解くと、
【００６４】
【数２５】
【００６５】
となる。
ここで、δは、クロネッカーのデルタである。すなわち、（Δｘi−Δｘi-1）＝０のときδ（Δｘi−Δｘi-1）＝０であり、（Δｘi−Δｘi-1）≠０のときδ（Δｘi−Δｘi-1）＝１である。Ｏ（Δｘm）は、Δｘに関してｍ次のオーダの項を表している。Δｘは、十分小さな値（１未満）であるため、ある項の次数が高いほど、その項の値は小さな値となる。
【００６６】
なお、式（２３）〜式（２５）は、空間についての差分近似を考察しているが、時間についても同様である。すなわち、ｘ座標を時間ｔに置き換え、Δｘを電磁波解析における時間刻み幅Δｔに置き換え、計算位置ｉを時刻ｎに置き換えればよい。ここで、時間刻み幅Δｔは常に一定であるものとする。すると、（Δｔi−Δｔi-1）＝０、すなわちδ＝０である。
【００６７】
以上のテーラー級数の式に基づいて、まず、セルサイズが均一な場合の計算式を求める。
図２は、均一セルにおける電磁界計算位置の配置図である。図２では、電界計算位置２１を黒棒で示し、磁界計算位置２２を白棒で示している。図２に示すように、電界計算位置２１と磁界計算位置２２とが、Ｘ軸上に交互に配置されている。この例では、隣接する２つの電界計算位置２１に挟まれる領域がセルである。従って、電界計算位置２１の間隔が空間刻み幅（セルサイズ）Δｘである。
【００６８】
均一セルの場合には、複数の電界計算位置２１は等間隔に配置されている。隣接する電界計算位置２１の中間に、磁界計算位置２２が配置されている。なお、ｉ番目のセルの中央に配置された磁界計算位置２２は、「ｉ＋１／２」で表すものとする。
【００６９】
空間のセルのサイズが均一な場合、（Δｘi−Δｘi-1）＝０となる。すなわち、δ＝０である。従って、均一セルの場合には、式（１）、式（２）に示した１次元のマックスウェル方程式に対して、時間と空間とに関して式（２５）に基づく差分法を適用すると、
【００７０】
【数２６】
【００７１】
【数２７】
【００７２】
【数２８】
【００７３】
【数２９】
【００７４】
となる。
このように、解析対象となる空間を均一なサイズのセルで一様に分割した場合、電磁界成分の時間と空間についての偏微分係数は、式（２７）および式（２９）に示す中心差分で近似され、打ち切り誤差は２次のオーダである。
【００７５】
一方、大きさの異なるセルが存在する場合、セルの大きさが相異なる領域の境界上での偏微分係数は、空間についての差分近似が中心差分からずれる。
図３は、不均一セルにおける電磁界計算位置の配置図である。図３では、電界計算位置２３を黒棒で示し、磁界計算位置２４を白棒で示している。図３に示すように、電界計算位置２３と磁界計算位置２４とが、Ｘ軸上に交互に配置されている。この例では、隣接する２つの電界計算位置２３に挟まれる領域がセルである。
【００７６】
図３では、電界計算位置２３の間隔が不定であり、不均一なセルが形成されている。図３の例では、ｉ−１番目のセルの電界計算位置（ｉ−１）とｉ番目のセルの電界計算位置（ｉ）との距離よりも、ｉ番目のセルの電界計算位置（ｉ）とｉ＋１番目のセルの電界計算位置（ｉ＋１）との距離の方が長い。すなわち、ｉ−１番目のセルのセルサイズΔｘi-1とｉ番目のセルのセルサイズΔｘiとは異なる大きさである。すなわち、電界計算位置（ｉ）は、セルサイズ境界である。
【００７７】
セルサイズ境界である電界計算位置（ｉ）における磁界の偏微分係数は、ｉ−１番目のセルの中心に設定された磁界計算位置（ｉ−１／２）とｉ番目のセルの中心に設定された磁界計算位置（ｉ＋１／２）とに基づいて計算される。ところが、電界計算位置（ｉ）は、磁界計算位置（ｉ−１／２）と磁界計算位置（ｉ＋１／２）との中点とずれている。
【００７８】
そのため、ｉ番目のセルの電界計算位置（ｉ）での磁界の偏微分係数を求める式（２９）は、以下の式（３０）に置き換えられる。
【００７９】
【数３０】
【００８０】
式（３０）に示すように、セルサイズが不均一な場合、打ち切り誤差が１次のオーダになってしまう。この打ち切り誤差の増加は、解析対象によっては著しい計算誤差の原因となる。
【００８１】
そこで、本発明では、打ち切り誤差を低減する方法として、高次の差分近似を用いる方法を部分的に適用する。高次の差分近似については、「W.J.Goedheer and J.H.H.M.Potters,"A compact finite difference scheme on a non-equidistance mesh,"Journal of Computational Physics, Vol. 61, pp.269-279, 1985.」や「A.Taflove, Advances in Computational Electro- dynamics, MA, Artech House, pp.71, 1998.」に説明されている。
【００８２】
本発明では、高次の陰的スキームによる差分近似を、以下のようにして導出した。
まず、次のテーラー級数を考える。
【００８３】
【数３１】
【００８４】
【数３２】
【００８５】
【数３３】
【００８６】
【数３４】
【００８７】
｛式（３１）｝−｛式（３２）｝から、
【００８８】
【数３５】
【００８９】
を得る。同様に、｛式（３３）｝−｛式（３４）｝から、
【００９０】
【数３６】
【００９１】
を得る。｛式（３５）｝−｛式（３６）｝×｛（Δｘi−Δｘi-1）／８｝から、
【００９２】
【数３７】
【００９３】
を得る。｛式（３３）｝×Δｘi-1＋｛式（３４）｝×Δｘiから、
【００９４】
【数３８】
【００９５】
を得る。｛式（３７）｝×Δｘi-1Δｘi−｛式（３８）｝×｛（−２Δｘi 2＋５ΔｘiΔｘi-1−２Δｘi-1 2）／２４｝から、
【００９６】
【数３９】
【００９７】
を得る。式（３９）の両辺を｛（Δｘi-1＋Δｘi-1）Δｘi-1 Δｘi-1｝／２で割って整理すると、
【００９８】
【数４０】
【００９９】
となる。なお、Ａ，Ｂ，Ｃは、式（２０）〜式（２２）に示した通りである。
ここで、均一セルであればδ＝０となり、電界および磁界のそれぞれに関して、
【０１００】
【数４１】
【０１０１】
【数４２】
【０１０２】
が得られる。なお、
【０１０３】
【数４３】
【０１０４】
である。
式（４０）、式（４１）および式（４２）による電磁波解析を行えば、打ち切り誤差は、式（４１）と式（４２）とを用いた均一セルでは４次のオーダとなる。不均一セルでは、セルサイズ境界において式（４０）が適用され、部分的に３次のオーダとなる。すなわち、従来の方法に比べて打ち切り誤差を低減させることが可能となる。
【０１０５】
なお、式（４０）、式（４１）および式（４２）は陰解法である。すなわち、１つの式の中の複数の項に、同じ時刻の異なる計算位置の偏微分係数が含まれている。そのため、任意の偏微分係数を算出するには、以下の式（４４）で表現される連立一次方程式を解く必要が生じる。
【０１０６】
【数４４】
【０１０７】
このような連立方程式の解を、解析領域のすべてに関して求めると、計算時間の増加を招いてしまう。そこで、式（４０）を局所的に適用することを考える。
はじめに、セルサイズ境界以外の偏微分係数を式（２９）により、２次のオーダの打ち切り誤差で算出する。式（２９）は、すでに求められている磁界の値から、直接算出することができる。
【０１０８】
次に、セルサイズ境界上の偏微分係数を式（４０）により算出する。このとき、式（４０）の左辺の第１項および第３項は、式（２９）によりすでに計算されているので、式（４０）の左辺の第２項を直ちに計算することができる。
【０１０９】
なお、式（４０）の左辺の第１項、第３項は、打ち切り誤差が２次のオーダで計算されているので、式（４０）の左辺の第２項の打ち切り誤差も２次のオーダである。
【０１１０】
式（４０）の左辺の第２項を算出する式は、以下のようになる。
【０１１１】
【数４５】
【０１１２】
式（４５）から３次のオーダの項の計算を省くことで、式（１９）が得られる。
このように、本発明では、セルサイズ境界に限定して高次元の差分近似法を適用する。これにより、計算時間の増加を最小限に抑えながら、打ち切り誤差を低減させることが可能となる。
【０１１３】
なお、上記の説明では簡単のために、１次元（ｘ軸方向に伝搬する電磁波）の場合について説明したが、この手法は、２次元や３次元にも容易に拡張することができる。
【０１１４】
また、セルサイズが連続して変化する場合についても適用可能である。ただし、セルサイズが連続して変化する場合には、式（１９）を計算するために連立一次方程式を解く必要がある。
【０１１５】
以下に、本発明を実現するためのコンピュータについて具体的に説明する。
図４は、本発明の実施の形態に用いるコンピュータのハードウェア構成例を示す図である。コンピュータ１００は、ＣＰＵ(Central Processing Unit)１０１によって装置全体が制御されている。ＣＰＵ１０１には、バス１０７を介してＲＡＭ(Random Access Memory)１０２、ハードディスクドライブ（ＨＤＤ:Hard Disk Drive）１０３、グラフィック処理装置１０４、入力インタフェース１０５、および通信インタフェース１０６が接続されている。
【０１１６】
ＲＡＭ１０２は、ＣＰＵ１０１に実行させるＯＳ(Operating System)のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、ＲＡＭ１０２には、ＣＰＵ１０１による処理に必要な各種データが格納される。ＨＤＤ１０３は、ＯＳやアプリケーションプログラムが格納される。
【０１１７】
グラフィック処理装置１０４には、モニタ１１が接続されている。グラフィック処理装置１０４は、ＣＰＵ１０１からの命令に従って、画像をモニタ１１の画面に表示させる。入力インタフェース１０５には、キーボード１２とマウス１３とが接続されている。入力インタフェース１０５は、キーボード１２やマウス１３から送られてくる信号を、バス１０７を介してＣＰＵ１０１に送信する。
【０１１８】
通信インタフェース１０６は、ネットワーク１０に接続されている。ネットワーク１０は、たとえばインターネットのような広域ネットワークである。通信インタフェース１０６は、ネットワーク１０を介して、他のコンピュータとの間でデータの送受信を行う。
【０１１９】
以上のようなハードウェア構成によって、本実施の形態の処理機能を実現することができる。たとえば、電磁波解析プログラムをコンピュータ１００に実行させることにより、コンピュータ１００を電磁波解析装置として機能させることができる。
【０１２０】
図５は、電磁波解析装置として動作するコンピュータの機能ブロック図である。コンピュータ１００には、入力装置１２ａと出力装置１１ａとが接続されている。
【０１２１】
入力装置１２ａは、キーボード１２やマウス１３などのデータを入力するための機器である。ユーザは、入力装置１２ａを介して、電磁波解析の初期値を入力することができる。
【０１２２】
出力装置１１ａは、モニタ１１などのデータ出力機器である。コンピュータ１００で計算された電磁界の過渡的な挙動を示す情報は、出力装置１１ａを介して出力される。
【０１２３】
コンピュータ１００は、初期値設定部１１１，電磁界データ記憶部１１２、セルサイズ判定部１１３、および電磁界計算部１１４で構成される。
初期値設定部１１１は、入力装置１２ａから入力された電磁波解析の初期値を、電磁界データ記憶部１１２に設定する。初期値としては、解析領域の定義（セルサイズ、誘電率、透磁率、電気伝導率など）や、電磁界の初期状態などが設定される。
【０１２４】
電磁界データ記憶部１１２は、たとえば、コンピュータ１００のＲＡＭ１０２内の記憶領域である。電磁界データ記憶部１１２は、初期値設定部１１１によって設定された初期値や、電磁界計算部１１４によって計算された電磁界の過渡的な挙動を示す情報を保持する。
【０１２５】
セルサイズ判定部１１３は、電磁界データ記憶部１１２に設定されたセルサイズの情報を取得し、セルサイズ境界を判定する。セルサイズ判定部１１３は、セルサイズ境界の場所を、電磁界計算部１１４に通知する。
【０１２６】
電磁界計算部１１４は、電磁界データ記憶部１１２に設定された電磁界データに基づいて、電磁界の過渡的な挙動を計算する。電磁界計算部１１４は、計算結果を、逐一電磁界データ記憶部１１２に格納する。また、所定の時間分の計算が終了すると、電磁界計算部１１４は、電磁界の過渡的な挙動を示す情報を、解析結果として出力装置１１ａに対して出力する。解析結果は、時間刻み幅毎の各時刻における所定の計算位置での磁界や電界の値である。
【０１２７】
図６は、本実施の形態による電磁波解析処理のフローチャートである。以下に、図６に示す処理をステップ番号に沿って説明する。
［ステップＳ１１］初期値設定部１１１は、入力装置１２ａからの入力に応じて初期設定を行う。初期設定の内容は、セル数、セルサイズ、各セルの材料定数（誘電率ε、透磁率μ、電気伝導率σなど）、境界条件（吸収境界条件、電気壁境界条件など）、励振条件（波源の位置、波形など）、解析終了時間、時間刻み幅などの設定ならびに電磁界計算のための係数の算出である。
【０１２８】
なお、このとき、電磁波解析シミュレーションの時刻Ｔを０に設定する。
［ステップＳ１２］セルサイズ判定部１１３は、初期設定で設定された各セルに関して、セルサイズを判定する。セルサイズの判定により、両側のセルサイズが等しい電界配置地点（電界の計算位置）と、両側のセルサイズが異なる電界配置地点とが区別される。
【０１２９】
［ステップＳ１３］電磁界計算部１１４は、両側のセルサイズが等しい電界配置地点における磁界の偏微分係数を、式(１８)によって計算する。
［ステップＳ１４］電磁界計算部１１４は、両側のセルサイズが異なる電界配置地点における磁界の偏微分係数を、式(１９)によって計算する。
【０１３０】
［ステップＳ１５］電磁界計算部１１４は、ステップＳ１３，Ｓ１４で計算した磁界の偏微分係数を用いて、式(１７)によって電界の計算を行う。
［ステップＳ１６］電磁界計算部１１４は、境界条件に依存した電界を、別途その条件式に基づいて計算を行う。その後、電磁界計算部１１４は、時刻ＴをΔＴ／２だけ進める。
【０１３１】
［ステップＳ１７］電磁界計算部１１４は、式（１５）および（１６）による磁界の計算を行う。
［ステップＳ１８］電磁界計算部１１４は、境界条件に依存した磁界を、別途その条件式に基づいて計算を行う。その後、電磁界計算部１１４は、時刻ＴをΔＴ／２だけ進める。
【０１３２】
［ステップＳ１９］電磁界計算部１１４は、時刻Ｔが、あらかじめ設定した解析終了時刻Ｔmaxを超えたか否か判断する。時刻Ｔが解析終了時刻Ｔmaxを超えた場合には、処理がステップＳ２０に進められる。時刻Ｔが解析終了時刻Ｔmaxを超えていない場合には、処理がステップＳ１３に進められる。
【０１３３】
［ステップＳ２０］電磁界計算部１１４は、計算結果を出力装置１１ａへ出力する。
以上のような処理により、様々な環境での電磁波の挙動を計算することができる。以下に、電磁波解析の計算事例について説明する。
【０１３４】
［平面波の伝播］
まず、平面波の伝搬についての計算事例を説明する。
図７は、平面波の伝搬を計算するモデルを示す図である。図７では、右方向がｘ軸の正、上方向がｙ軸の正、紙面手前方向がｚ軸の正である。ここでは、２次元ＴＥ波を扱うものとする。図７において格子状に表されたセルの辺上に電界Ｅｘ、Ｅｙが配置され、セルの中心に磁界Ｈｚが配置されている。
【０１３５】
解析領域３０は、セルサイズの大きな領域３１とセルサイズの小さな領域３２とで構成される。領域３１のセルサイズはΔｘC×Δｙ、領域３２のセルサイズはΔｘF×Δｙである。セルサイズの異なる２つの領域３１，３２は、セルサイズ境界３３で接している。
【０１３６】
Ｙ方向の終端には電気壁境界条件を設定し、Ｘ方向の終端はセルサイズ境界から十分遠方になるようにしてある。これは、Ｘ方向の終端からの反射の影響を除去するためである。
【０１３７】
セルサイズ境界は、電磁波解析の計算のために解析領域に便宜的に設けられたものであり、解析対象自身の特性ではない。そのため、解析領域に電磁波を伝搬させた場合、計算上、セルサイズ境界からの反射波が検出されないことが理想である。
【０１３８】
そこで、セルサイズ境界から−Ｘ方向に十分に離れた地点において、電界Ｅｙをガウス波で励振し、パルスを＋Ｘ方向へ伝搬させ、セルサイズ境界における反射量を計算した。Δｙ＝Δｘc＝1.0mm、ΔｘF＝0.1 mmの場合（Δｘc／ΔｘF＝１０）の反射特性は、以下の通りである。
【０１３９】
図８は、セルサイズ境界における反射特性を示す図である。図８では、横軸に周波数（ＧＨｚ）を示し、縦軸に反射量（ｄＢ）を示している。反射量の数値は、−１２０ｄＢに近いほどセルサイズ境界からの反射が少ないことを意味している。従来手法の反射特性は実線４１で表し、本発明の実施の形態に基づく手法の反射特性は、点線４２で表している。
【０１４０】
図８からも分かるように、本実施の形態によれば、従来手法に比べ反射量が２０〜４０ｄＢ低減できていることがわかる。これは、セルサイズ境界上の電界Ｅｙを計算する際に、打ち切り誤差が低減され領域の不連続性が緩和されたからであると考えられる。
【０１４１】
なお、この特性は、Δｙ＝Δｘc＝1.0mm、ΔｘF＝0.025mmの場合（Δｘc／ΔｘF＝４０）でもほとんど変わらない。
［パッチアンテナの反射特性］
次に、パッチアンテナの反射特性の解析結果について説明する。
【０１４２】
図９は、パッチアンテナの平面図である。図９の例では、解析領域５０内にパッチアンテナ５１が配置されている。
解析領域５０は、縦横118.0mmの正方形の領域である。パッチアンテナ５１の先端部は、縦横共に11.8mmの大きさであり、解析領域５０の４方の縁から50.0mm離れている。
【０１４３】
パッチアンテナ５１は、幅1.8mmの信号線５２で外部に接続される。パッチアンテナ５１は、信号線５２との接続部において、信号線５２の両側に深さ4.0mm、幅1.0mmの溝が形成されている。
【０１４４】
なお、パッチアンテナ５１が形成されている基板の厚さは0.8mm、比誘電率は3.274であるものとする。
図９に示したパッチアンテナ５１は、社団法人電子情報通信学会マイクロ波シミュレータ研究会において規範問題とされているマイクロストリップ給電パッチアンテナである。たとえば、「並木武文, 坂口拓史, 伊藤公一, “FDTD法を用いたパッチアンテナ解析における計算精度についての一考察,”信学技報, AP99-12, pp.17-22, May 1999.」や「田口光雄, “アンテナ設計から見た電磁界シミュレータの使い方と評価,” 信学誌, vol.83, no.11, pp.878-883, Nov. 2000.」において、パッチアンテナ５１を用いたシミュレーション例が詳しく説明されている。
【０１４５】
本実施の形態では、このパッチアンテナ５１を計算するために、不均一セルを使って領域分割し、計算モデルを作成する。
図１０は、パッチアンテナの計算モデルの一例を示す図である。図中、右方向がｘ軸の正、上方向がｙ軸の正、紙面手前方向がｚ軸の正である。
【０１４６】
ここで、放射パッチであるパッチアンテナ５１の端部近傍は、電磁界の特異性を考慮して非常に微細なセルに分割してある。セルの分割手法に関しては、「E.M.Daniel and C.J.Railton, "Fast finite difference time domain analysis of microstrip patch antennas,"IEEE AP-S Digest, pp.414-417, 1991.」や「大西輝夫, 柏達也, 内藤行雄, 細矢良雄, “同軸給電パッチアンテナのFD-TD解析における効率化の一検討,”信学総大, B-1-128, 1997.」で説明されている。
【０１４７】
Ｚ方向についても同様に、パッチアンテナ５１の近傍は非常に微細なセルに分割している。図１０の計算モデルにおける全セル数は76380個である。
図１１は、図１０の計算モデルにおける最小ならびに最大セルサイズを示す図である。ｘ軸方向のセルサイズΔＸの最小値(Minimum cell)は、0.200mmである。ｘ軸方向のセルサイズΔＸの最大値(Maximum cell)は、2.000mmである。ｙ軸方向のセルサイズΔＹの最小値(Minimum cell)は、0.050mmである。ｙ軸方向のセルサイズΔＹの最大値(Maximum cell)は、2.000mmである。ｚ軸方向のセルサイズΔＺの最小値(Minimum cell)は、0.069mmである。ｚ軸方向のセルサイズΔＺの最大値(Maximum cell)は、1.200mmである。
【０１４８】
また、解析領域の終端にはMurの１次の吸収境界条件を設定した。Murの１次の吸収境界条件に関しては、「G.Mur,"Absorbing Boundary Conditions for the finite-difference approximation of the time-domain electromagnetic field equations,"IEEE Trans.EMC-23, pp.377-382, 1981.」で説明されている通りである。
【０１４９】
このようなパッチアンテナ５１のサンプルを用意し、そのサンプルを用いて反射特性を測定すると共に、本実施の形態における解析手法と従来の解析手法とにより反射特性を計算した。
【０１５０】
図１２は、図１０に示す計算モデルの反射特性の計算結果と測定結果とを示す図である。図１２では、横軸に周波数（ＧＨｚ）を示し、縦軸に反射量（ｄＢ）を示している。サンプルによる測定結果を破線６１で表し、本発明の実施の形態に基づく手法による計算結果を点線６２で表し、従来手法による計算結果を実線６３で表している。
【０１５１】
図１２からも分かるように、本発明の実施の形態に基づく手法で電磁波解析を行えば、従来に比べて、測定結果に近い計算結果を得ることができる。なお、測定値の信頼性については、十分に検討してある。測定値の信頼性の検討内容は、「並木武文, 坂口拓史, 伊藤公一, “FDTD法を用いたパッチアンテナ解析における計算精度についての一考察,”信学技報, AP99-12, pp.17-22, May 1999.」（以下、この文献を「計算精度考察例」と呼ぶ）に説明されているものと同様である。
【０１５２】
ここで、図１２に示した反射特性から共振周波数を抽出し、測定値に対する計算値の相対誤差を算出した。
図１３は、共振周波数と相対誤差を示す図である。測定値に基づく共振周波数(Resonant frequency)は、7.0100ＧＨｚである。従来の手法の計算値に基づく共振周波数は6.8624ＧＨｚであり、相対誤差(Relative error)は2.10％である。本発明の実施の形態に係る手法の計算値に基づく共振周波数は6.9756ＧＨｚであり、相対誤差は0.49%である。
【０１５３】
図１２、図１３から明らかなように、図１０の計算モデルを用いた場合、本手法によって大幅な精度改善が図れることがわかる。
今回の計算で必要としたＣＰＵ時間とメモリサイズを以下に示す。
【０１５４】
図１４は、計算に使用されたＣＰＵ時間とメモリサイズとを示す図である。従来の手法で使用したＣＰＵ時間は493秒であり、使用したメモリサイズは5.5ＭＢである。本発明の実施の形態の手法で使用したＣＰＵ時間は666秒であり、使用したメモリサイズは6.4ＭＢである。
【０１５５】
なお、計算には、ワークステーション(CPU：UltraSPARC II 360MHz)を使用した。本発明による手法の場合、通常のＦＤＴＤ法と比較して、ＣＰＵ時間は約1.35倍、メモリは1.16倍増加している。しかし、その精度改善効果からすれば問題となる程度ではないと考えられる。
【０１５６】
前記の計算精度考察例における通常のＦＤＴＤ法を用いた同様の検討では、不均一セルを使用した場合は、共振周波数を相対誤差1.0%以下で算出することが出来なかった。そして、共振周波数を相対誤差1.0%以下で算出するためには、莫大な数の微細な均一セルを使用しなければならなかった。この場合、ほぼ同性能のワークステーションを使用して、ＣＰＵ時間が600分以上、メモリが250ＭＢ必要であった。これらの結果を鑑みれば、本発明による手法の有効性は明らかである。
【０１５７】
なお、上記の処理機能は、コンピュータによって実現することができる。その場合、電磁波解析装置が有すべき機能の処理内容を記述したプログラムが提供される。そのプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリなどがある。磁気記録装置には、ハードディスク装置（ＨＤＤ）、フレキシブルディスク（ＦＤ）、磁気テープなどがある。光ディスクには、ＤＶＤ(Digital Versatile Disc)、ＤＶＤ−ＲＡＭ(Random Access Memory)、ＣＤ−ＲＯＭ(Compact Disc Read Only Memory)、ＣＤ−Ｒ(Recordable)／ＲＷ(ReWritable)などがある。光磁気記録媒体には、ＭＯ(Magneto-Optical disc)などがある。
【０１５８】
プログラムを流通させる場合には、たとえば、そのプログラムが記録されたＤＶＤ、ＣＤ−ＲＯＭなどの可搬型記録媒体が販売される。また、プログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することもできる。
【０１５９】
プログラムを実行するコンピュータは、たとえば、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、自己の記憶装置に格納する。そして、コンピュータは、自己の記憶装置からプログラムを読み取り、プログラムに従った処理を実行する。なお、コンピュータは、可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することもできる。また、コンピュータは、サーバコンピュータからプログラムが転送される毎に、逐次、受け取ったプログラムに従った処理を実行することもできる。
【０１６０】
（付記１） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算する電磁波解析装置において、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断するセルサイズ判断手段と、
前記セルサイズ判断手段において周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算する第１の計算手段と、
前記セルサイズ判断手段において周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算する第２の計算手段と、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する出力手段と、
を有することを特徴とする電磁波解析装置。
【０１６１】
（付記２） 前記第２の計算手段は、前記第２の計算手法として、陰的スキームを用いることを特徴とする付記１記載の電磁波解析装置。
（付記３） 前記第２の計算手段は、前記陰的スキームにより、周囲のセルのサイズが異なる計算位置の電磁界の偏微分係数を算出することを特徴とする付記２記載の電磁波解析装置。
【０１６２】
（付記４） 前記第１の計算手段は、周囲のセルのサイズが同じである計算位置の電磁界の偏微分係数を計算し、
前記第２の計算手段は、周囲のセルのサイズが異なる計算位置の周囲の電磁界の微分係数を前記第１の計算手段から取得し、前記第１の計算手段で計算された電磁界の微分係数を用いて、前記陰的スキームの解を求めることを特徴とする付記３記載の電磁波解析装置。
【０１６３】
（付記５） 前記セルサイズ判断手段は、解析領域に設定された座標軸毎に、隣接するセルの各軸方向の幅の同一性を判断することを特徴とする付記１記載の電磁波解析装置。
【０１６４】
（付記６） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析プログラムにおいて、
コンピュータに、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算し、
周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算し、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する、
処理を実行させることを特徴とする電磁波解析プログラム。
【０１６５】
（付記７） 前記第２の計算手法として、陰的スキームを用いることを特徴とする付記６記載の電磁波解析プログラム。
（付記８） 第２の計算手法による計算の際には、前記陰的スキームにより、周囲のセルのサイズが異なる計算位置の電磁界の偏微分係数を算出することを特徴とする付記７記載の電磁波解析プログラム。
【０１６６】
（付記９） 前記第１の計算手法による計算の際には、周囲のセルのサイズが同じである計算位置の電磁界の偏微分係数を計算し、
前記第２の計算手段は、周囲のセルのサイズが異なる計算位置の周囲の電磁界の微分係数を前記第１の計算手段から取得し、前記第１の計算手段で計算された電磁界の微分係数を用いて、前記陰的スキームの解を求めることを特徴とする付記８記載の電磁波解析プログラム。
【０１６７】
（付記１０） 前記セルサイズ判断手段は、解析領域に設定された座標軸毎に、隣接するセルの各軸方向の幅の同一性を判断することを特徴とする付記６記載の電磁波解析プログラム。
【０１６８】
（付記１１） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体において、
前記コンピュータに、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算し、
周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算し、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する、
処理を実行させることを特徴とする記録媒体。
【０１６９】
（付記１２） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析方法において、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算し、
周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算し、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する、
ことを特徴とする電磁波解析方法。
【０１７０】
【発明の効果】
以上説明したように本発明では、周囲のセルのサイズが異なる計算位置における電磁界を、他の計算位置よりも少ない誤差の計算手法で計算するようにしたため、セルのサイズが切り替わるセルサイズ境界における計算誤差を減少させることができ、むやみに計算量を増大させることなく、誤差の少ない電磁波解析を行うことができる。
【図面の簡単な説明】
【図１】本発明の原理構成図である。
【図２】均一セルにおける電磁界計算位置の配置図である。
【図３】不均一セルにおける電磁界計算位置の配置図である
【図４】本発明の実施の形態に用いるコンピュータのハードウェア構成例を示す図である。
【図５】電磁波解析装置として動作するコンピュータの機能ブロック図である。
【図６】本実施の形態による電磁波解析処理のフローチャートである。
【図７】平面波の伝搬を計算するモデルを示す図である。
【図８】セルサイズ境界における反射特性を示す図である。
【図９】パッチアンテナの平面図である。
【図１０】パッチアンテナの計算モデルの一例を示す図である。
【図１１】図１０の計算モデルにおける最小ならびに最大セルサイズを示す図である。
【図１２】図１０に示す計算モデルの反射特性の計算結果と測定結果とを示す図である。
【図１３】共振周波数と相対誤差を示す図である。
【図１４】計算に使用されたＣＰＵ時間とメモリサイズとを示す図である。
【符号の説明】
１ 解析領域
２ セルサイズ判断手段
３ 第１の計算手段
４ 第２の計算手段
５ 出力手段
１０ ネットワーク
１１ モニタ
１２ キーボード
１３ マウス
１００ コンピュータ
１０１ ＣＰＵ
１０２ ＲＡＭ
１０３ ハードディスクドライブ
１０４ グラフィック処理装置
１０５ 入力インタフェース
１０６ 通信インタフェース
１０７ バス
【発明の属する技術分野】
本発明はマックスウェル方程式を時間と空間について差分法で解く電磁波解析装置および電磁波解析プログラムに関し、特に空間を不均一なサイズのセルに分割してマックスウェル方程式の解を求める電磁波解析装置および電磁波解析プログラムに関する。
【０００２】
【従来の技術】
電磁波の過渡的な挙動を、電子計算機を使った数値計算によって解析する方法の一つに時間領域差分法（Finite-Difference Time-Domain Method、以下ＦＤＴＤ法という）がある。この手法は、マックスウェル方程式を時間と空間について差分法で解く方法であり、適用範囲が広いこと、計算精度や計算効率が良好であることなどから、今日幅広く利用されている。ＦＤＴＤ法については、「K.S.Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media,"IEEE Trans. AP-14, pp.302-307, 1966.」や「A.Taflove, Computational Electrodynamics, MA, Artech House, 1995.」に詳しい。
【０００３】
ＦＤＴＤ法は、マックスウェル方程式の２つの回転の式を、時間ｔと空間（ｘ−ｙ−ｚ空間）について差分法で解く方法である。
以下に、マックスウェル方程式について説明する。なお、以下の式において、Ｅｘは電界のｘ軸方向の成分、Ｅｙは電界のｙ軸方向の成分、Ｅｚは電界のｚ軸方向の成分、Ｈｘは磁界のｘ軸方向の成分、Ｈｙは磁界のｙ軸方向の成分、Ｈｚは磁界のｚ軸方向の成分、μは透磁率、εは誘電率、σは電気伝導率である。
【０００４】
１次元（Ｘ軸方向に伝搬する電磁波）のマックスウェル方程式は、次の式で表される。
【０００５】
【数１】
【数２】
２次元（Ｘ軸、Ｙ軸方向に伝搬する電磁波）のマックスウェル方程式は、ＴＭ波(Transverse Magnetic wave)については次の式で表される。
【０００８】
【数３】
【数４】
【数５】
２次元（Ｘ軸、Ｙ軸方向に伝搬する電磁波）のマックスウェル方程式は、ＴＥ波(Transverse Electric wave)については次の式で表される。
【００１２】
【数６】
【数７】
【数８】
３次元のマックスウェル方程式は、次の式で表される。
【００１６】
【数９】
【数１０】
【数１１】
【数１２】
【数１３】
【数１４】
ＦＤＴＤ法では、通常、解析領域を直方体のセル（差分格子）に分割することにより、空間的な離散化を行う。このとき、セルの大きさが小さいほど、電磁波の過渡的な挙動が急激に変わる領域においても、少ない誤差で計算することができる。
【００２３】
ところが、空間を均一なセルで一様に分割した場合、しばしば多大なセル数を必要とし、計算時間や計算機メモリの増加を引き起こす。そこで、大きさの異なるセルを組み合わせて用いるいわゆる不均一なセルで解析領域を分割することにより、セル数を節約することがたびたび行われる。
【００２４】
【発明が解決しようとする課題】
しかし、解析領域をサイズが不均一なセルで分割した場合、セルサイズが相異なる領域の境界（以下、セルサイズ境界という）において、誤差が増加することが知られている。
【００２５】
すなわち、セルサイズ境界では、隣接するセルに基づいた差分近似の中心位置とセルサイズ境界の位置とがずれてしまう。この位置のずれに基づく差分近似の打ち切り誤差が、電磁波解析における誤差の増加を招いている。
【００２６】
セルサイズが同一の領域では、差分近似の打ち切り誤差は２次のオーダであるが、セルサイズ境界における差分近似の打ち切り誤差は１次のオーダである。これは、解析対象よっては、計算精度を著しく低下させる要因となる。
【００２７】
本発明はこのような点に鑑みてなされたものであり、不均一セルを用いた電磁波解析の誤差を低減させた電磁波解析装置を提供することを目的とする。
【００２８】
【課題を解決するための手段】
本発明では上記課題を解決するために、図１に示すような電磁波解析装置が提供される。本発明に係る電磁波解析装置は、解析領域１を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するものである。
【００２９】
本発明の電磁波解析装置は、セルサイズ判断手段２、第１の計算手段３、第２の計算手段４および出力手段５を有している。
セルサイズ判断手段２は、解析領域１に設定された各計算位置Ｐ１，Ｐ２に関し、周囲のセルのサイズの同一性を判断する。第１の計算手段３は、セルサイズ判断手段２において周囲のセルのサイズが同じであると判断された計算位置Ｐ１における電磁界を、第１の計算手法で計算する。第２の計算手段４は、セルサイズ判断手段２において周囲のセルのサイズが異なると判断された計算位置Ｐ２における電磁界を、第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算する。出力手段５は、第１の計算手段と第２の計算手段とにより計算された各計算位置における電磁界の値を出力する。
【００３０】
このような電磁波解析装置によれば、サイズの同じセルサイズが均一な領域内の計算位置においては、セルサイズ判断手段２によって、周囲のセルのサイズが同一であると判断される。すると、その計算位置の電磁界は、第１の計算手段３により、第１の計算手法で計算される。一方、異なるサイズのセル同士の隣接する計算位置では、セルサイズ判断手段２によって、周囲のセルのサイズが異なると判断される。すると、その計算位置の電磁界は、第２の計算手段４により、第１の計算手法よりも誤差の少ない第２の計算手法で計算される。第１の計算手段３および第２の計算手段４で計算された電磁界の値は、出力手段５によって出力される。
【００３１】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して説明する。
図１は、本発明の原理構成図である。本発明に係る電磁波解析装置は、電磁波の挙動を数値計算する場合に、マックスウェル方程式を時間と空間について差分法で解くものである。差分法では、解析領域１を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するものである。
【００３２】
本発明の電磁波解析装置は、セルサイズ判断手段２、第１の計算手段３、第２の計算手段４および出力手段５を有している。
セルサイズ判断手段２は、解析領域１に設定された各計算位置Ｐ１，Ｐ２に関し、周囲のセルのサイズの同一性を判断する。
【００３３】
たとえば、セルサイズ判断手段２は、解析領域１に設定された座標軸毎に、隣接するセルの各軸方向の幅の同一性を判断する。すなわち、セルサイズ判断手段２は、ある軸方向に隣接するセル同士のその軸方向へのセルの幅が同じであれば、その軸方向への電磁波解析においては、周囲のセルのサイズが同一であるものと判断する。一方、セルサイズ判断手段２は、ある軸方向に隣接するセル同士のその軸方向へのセルの幅が異なれば、その軸方向への電磁波解析においては、周囲のセルのサイズが異なるものと判断する。
【００３４】
第１の計算手段３は、セルサイズ判断手段２において周囲のセルのサイズが同じであると判断された計算位置Ｐ１における電磁界を、第１の計算手法で計算する。たとえば、第１の計算手段３は、電界または磁界の偏微分係数を算出する際に、陰的スキーム（陰解法に基づいた計算手法）を用いず、既に算出されている磁界の値から算出する。
【００３５】
第２の計算手段４は、セルサイズ判断手段２において周囲のセルのサイズが異なると判断された計算位置Ｐ２における電磁界を、第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算する。
【００３６】
たとえば、第２の計算手段４は、周囲のセルのサイズが異なる計算位置の電磁界の偏微分係数を算出する際に、高次の陰的スキームを用いて算出する。このとき、第２の計算手段４は、周囲のセルのサイズが異なる計算位置の周囲の電磁界の微分係数を前記第１の計算手段３から取得し、第１の計算手段３で計算された電磁界の微分係数を用いて、陰的スキームの解を求める。
【００３７】
出力手段５は、第１の計算手段３と第２の計算手段４とにより計算された各計算位置における電磁界の値を出力する。
このような電磁波解析装置によれば、セルサイズが均一な領域内の計算位置においては、セルサイズ判断手段２によって、周囲のセルのサイズが同一であると判断される。すると、その計算位置の電磁界は、第１の計算手段３により、第１の計算手法で計算される。一方、異なるサイズのセル同士の隣接する計算位置では、セルサイズ判断手段２によって、周囲のセルのサイズが異なると判断される。すると、その計算位置の電磁界は、第２の計算手段４により、第１の計算手法よりも誤差の少ない第２の計算手法で計算される。第１の計算手段３および第２の計算手段４で計算された電磁界の値は、出力手段５によって出力される。
【００３８】
このように、本発明では、セルの大きさが相異なる領域が存在する場合、セルの大きさが均一の領域における電磁界の計算式と、異なる大きさのセルの境界における電磁界の計算式とが異なる。セルの中心を磁場の計算位置とし、セル同士の境界を電界の計算位置とした場合、たとえば、以下の計算式で電磁界を求め出すことができる。
【００３９】
ここでは、簡単のため１次元問題について考える。また、電気伝導率σは０とする。
電磁波解析シミュレーションにおける時刻Δｔ（時間刻み幅）×ｎ（ｎは０以上の整数）までの電磁界の過渡的変化が計算されているときに、時刻Δｔ×（ｎ＋１）までの電磁界を求めるときの計算式を以下に示す。
【００４０】
なお、時刻Δｔ×ｎを単にｎと表し、時刻Δｔ×（ｎ＋１）をｎ＋１と表し、時刻Δｔ×（ｎ＋１／２）をｎ＋１／２と表すものとする。各物理量の右上に、その物理量の計算時刻（シミュレーション上の時刻）を示す。
【００４１】
また、ｉは０以上の整数であり、ｉ番目のセルの始点位置を示している。ｉ＋Δｘ（ｘ軸方向の空間刻み幅：セルサイズ）は、ｉ＋１番目のセルの始点位置である。ｉ＋１／２は、ｉ番目のセルの中心位置である。各物理量の右下に、その物理量の計算位置を示す。
【００４２】
まず、時刻ｎ＋１／２における磁界は、
【００４３】
【数１５】
で算出される。式（１５）における電界の偏微分係数は、
【００４５】
【数１６】
で算出される。
また、時刻ｎ＋１の電界は、
【００４７】
【数１７】
で算出される。式（１７）における磁界の偏微分係数は、セルサイズが均一な領域においては、
【００４９】
【数１８】
で算出される。式（１８）の差分近似式の打ち切り誤差は２次のオーダである。打ち切り誤差の詳細は後述する。
また、セルサイズ境界における偏微分係数は、
【００５１】
【数１９】
で算出される。なお、Ａ、Ｂ、Ｃは、以下の通りである。
【００５３】
【数２０】
【数２１】
【数２２】
式（１９）の差分近似式の打ち切り誤差は３次のオーダである。打ち切り誤差の詳細は後述する。
なお、セルサイズ境界における偏微分係数を求める場合、始めにセルサイズ境界以外の場所における空間についての偏微分係数を式（１６）および式（１８）を用いて差分近似する。次にセルサイズ境界上における空間についての偏微分係数を式（１９）を用いて差分近似する。
【００５７】
なお、式（１５）から式（２２）に示した計算は、ある特定の軸方向（ここではｘ方向）に着目し、Ｅ(電界)がセルの端、Ｈ(磁界)がセルの中心に配置した場合に成立する。Ｅ(電界)がセルの中心、Ｈ(磁界)がセルの端に配置してある場合は、式（１５）から式（２２）の各式のＥとＨならびにεとμを入れ替えればよい。
【００５８】
このように、本発明によれば、セルの大きさが相異なる領域の境界上における空間についての偏微分係数は、式（１９）を用いて差分近似しているため、差分近似の打ち切り誤差は２次のオーダとなる。なお、それ以外の場所における空間についての偏微分係数を差分近似する場合には、従来と同様の式（１６）および式（１８）を用いるが、これは中心差分であることが保証されるので、差分近似の打ち切り誤差は常に２次のオーダとなる。
【００５９】
これにより、打ち切り誤差の増加を防ぐことができる。すなわち、不均一セルを用いた電磁波解析全体としての計算精度が向上する。
ここで、上記計算式の導出方法について説明する。
【００６０】
ＦＤＴＤ法による差分近似を行う場合、関数ｆ（ｘ）について、次のテーラー級数を考える。
【００６１】
【数２３】
【数２４】
式（２３）から式（２４）を減算し、ｆ’について解くと、
【００６４】
【数２５】
となる。
ここで、δは、クロネッカーのデルタである。すなわち、（Δｘi−Δｘi-1）＝０のときδ（Δｘi−Δｘi-1）＝０であり、（Δｘi−Δｘi-1）≠０のときδ（Δｘi−Δｘi-1）＝１である。Ｏ（Δｘm）は、Δｘに関してｍ次のオーダの項を表している。Δｘは、十分小さな値（１未満）であるため、ある項の次数が高いほど、その項の値は小さな値となる。
【００６６】
なお、式（２３）〜式（２５）は、空間についての差分近似を考察しているが、時間についても同様である。すなわち、ｘ座標を時間ｔに置き換え、Δｘを電磁波解析における時間刻み幅Δｔに置き換え、計算位置ｉを時刻ｎに置き換えればよい。ここで、時間刻み幅Δｔは常に一定であるものとする。すると、（Δｔi−Δｔi-1）＝０、すなわちδ＝０である。
【００６７】
以上のテーラー級数の式に基づいて、まず、セルサイズが均一な場合の計算式を求める。
図２は、均一セルにおける電磁界計算位置の配置図である。図２では、電界計算位置２１を黒棒で示し、磁界計算位置２２を白棒で示している。図２に示すように、電界計算位置２１と磁界計算位置２２とが、Ｘ軸上に交互に配置されている。この例では、隣接する２つの電界計算位置２１に挟まれる領域がセルである。従って、電界計算位置２１の間隔が空間刻み幅（セルサイズ）Δｘである。
【００６８】
均一セルの場合には、複数の電界計算位置２１は等間隔に配置されている。隣接する電界計算位置２１の中間に、磁界計算位置２２が配置されている。なお、ｉ番目のセルの中央に配置された磁界計算位置２２は、「ｉ＋１／２」で表すものとする。
【００６９】
空間のセルのサイズが均一な場合、（Δｘi−Δｘi-1）＝０となる。すなわち、δ＝０である。従って、均一セルの場合には、式（１）、式（２）に示した１次元のマックスウェル方程式に対して、時間と空間とに関して式（２５）に基づく差分法を適用すると、
【００７０】
【数２６】
【数２７】
【数２８】
【数２９】
となる。
このように、解析対象となる空間を均一なサイズのセルで一様に分割した場合、電磁界成分の時間と空間についての偏微分係数は、式（２７）および式（２９）に示す中心差分で近似され、打ち切り誤差は２次のオーダである。
【００７５】
一方、大きさの異なるセルが存在する場合、セルの大きさが相異なる領域の境界上での偏微分係数は、空間についての差分近似が中心差分からずれる。
図３は、不均一セルにおける電磁界計算位置の配置図である。図３では、電界計算位置２３を黒棒で示し、磁界計算位置２４を白棒で示している。図３に示すように、電界計算位置２３と磁界計算位置２４とが、Ｘ軸上に交互に配置されている。この例では、隣接する２つの電界計算位置２３に挟まれる領域がセルである。
【００７６】
図３では、電界計算位置２３の間隔が不定であり、不均一なセルが形成されている。図３の例では、ｉ−１番目のセルの電界計算位置（ｉ−１）とｉ番目のセルの電界計算位置（ｉ）との距離よりも、ｉ番目のセルの電界計算位置（ｉ）とｉ＋１番目のセルの電界計算位置（ｉ＋１）との距離の方が長い。すなわち、ｉ−１番目のセルのセルサイズΔｘi-1とｉ番目のセルのセルサイズΔｘiとは異なる大きさである。すなわち、電界計算位置（ｉ）は、セルサイズ境界である。
【００７７】
セルサイズ境界である電界計算位置（ｉ）における磁界の偏微分係数は、ｉ−１番目のセルの中心に設定された磁界計算位置（ｉ−１／２）とｉ番目のセルの中心に設定された磁界計算位置（ｉ＋１／２）とに基づいて計算される。ところが、電界計算位置（ｉ）は、磁界計算位置（ｉ−１／２）と磁界計算位置（ｉ＋１／２）との中点とずれている。
【００７８】
そのため、ｉ番目のセルの電界計算位置（ｉ）での磁界の偏微分係数を求める式（２９）は、以下の式（３０）に置き換えられる。
【００７９】
【数３０】
式（３０）に示すように、セルサイズが不均一な場合、打ち切り誤差が１次のオーダになってしまう。この打ち切り誤差の増加は、解析対象によっては著しい計算誤差の原因となる。
【００８１】
そこで、本発明では、打ち切り誤差を低減する方法として、高次の差分近似を用いる方法を部分的に適用する。高次の差分近似については、「W.J.Goedheer and J.H.H.M.Potters,"A compact finite difference scheme on a non-equidistance mesh,"Journal of Computational Physics, Vol. 61, pp.269-279, 1985.」や「A.Taflove, Advances in Computational Electro- dynamics, MA, Artech House, pp.71, 1998.」に説明されている。
【００８２】
本発明では、高次の陰的スキームによる差分近似を、以下のようにして導出した。
まず、次のテーラー級数を考える。
【００８３】
【数３１】
【数３２】
【数３３】
【数３４】
｛式（３１）｝−｛式（３２）｝から、
【００８８】
【数３５】
を得る。同様に、｛式（３３）｝−｛式（３４）｝から、
【００９０】
【数３６】
を得る。｛式（３５）｝−｛式（３６）｝×｛（Δｘi−Δｘi-1）／８｝から、
【００９２】
【数３７】
を得る。｛式（３３）｝×Δｘi-1＋｛式（３４）｝×Δｘiから、
【００９４】
【数３８】
を得る。｛式（３７）｝×Δｘi-1Δｘi−｛式（３８）｝×｛（−２Δｘi 2＋５ΔｘiΔｘi-1−２Δｘi-1 2）／２４｝から、
【００９６】
【数３９】
を得る。式（３９）の両辺を｛（Δｘi-1＋Δｘi-1）Δｘi-1 Δｘi-1｝／２で割って整理すると、
【００９８】
【数４０】
となる。なお、Ａ，Ｂ，Ｃは、式（２０）〜式（２２）に示した通りである。
ここで、均一セルであればδ＝０となり、電界および磁界のそれぞれに関して、
【０１００】
【数４１】
【数４２】
が得られる。なお、
【０１０３】
【数４３】
である。
式（４０）、式（４１）および式（４２）による電磁波解析を行えば、打ち切り誤差は、式（４１）と式（４２）とを用いた均一セルでは４次のオーダとなる。不均一セルでは、セルサイズ境界において式（４０）が適用され、部分的に３次のオーダとなる。すなわち、従来の方法に比べて打ち切り誤差を低減させることが可能となる。
【０１０５】
なお、式（４０）、式（４１）および式（４２）は陰解法である。すなわち、１つの式の中の複数の項に、同じ時刻の異なる計算位置の偏微分係数が含まれている。そのため、任意の偏微分係数を算出するには、以下の式（４４）で表現される連立一次方程式を解く必要が生じる。
【０１０６】
【数４４】
このような連立方程式の解を、解析領域のすべてに関して求めると、計算時間の増加を招いてしまう。そこで、式（４０）を局所的に適用することを考える。
はじめに、セルサイズ境界以外の偏微分係数を式（２９）により、２次のオーダの打ち切り誤差で算出する。式（２９）は、すでに求められている磁界の値から、直接算出することができる。
【０１０８】
次に、セルサイズ境界上の偏微分係数を式（４０）により算出する。このとき、式（４０）の左辺の第１項および第３項は、式（２９）によりすでに計算されているので、式（４０）の左辺の第２項を直ちに計算することができる。
【０１０９】
なお、式（４０）の左辺の第１項、第３項は、打ち切り誤差が２次のオーダで計算されているので、式（４０）の左辺の第２項の打ち切り誤差も２次のオーダである。
【０１１０】
式（４０）の左辺の第２項を算出する式は、以下のようになる。
【０１１１】
【数４５】
式（４５）から３次のオーダの項の計算を省くことで、式（１９）が得られる。
このように、本発明では、セルサイズ境界に限定して高次元の差分近似法を適用する。これにより、計算時間の増加を最小限に抑えながら、打ち切り誤差を低減させることが可能となる。
【０１１３】
なお、上記の説明では簡単のために、１次元（ｘ軸方向に伝搬する電磁波）の場合について説明したが、この手法は、２次元や３次元にも容易に拡張することができる。
【０１１４】
また、セルサイズが連続して変化する場合についても適用可能である。ただし、セルサイズが連続して変化する場合には、式（１９）を計算するために連立一次方程式を解く必要がある。
【０１１５】
以下に、本発明を実現するためのコンピュータについて具体的に説明する。
図４は、本発明の実施の形態に用いるコンピュータのハードウェア構成例を示す図である。コンピュータ１００は、ＣＰＵ(Central Processing Unit)１０１によって装置全体が制御されている。ＣＰＵ１０１には、バス１０７を介してＲＡＭ(Random Access Memory)１０２、ハードディスクドライブ（ＨＤＤ:Hard Disk Drive）１０３、グラフィック処理装置１０４、入力インタフェース１０５、および通信インタフェース１０６が接続されている。
【０１１６】
ＲＡＭ１０２は、ＣＰＵ１０１に実行させるＯＳ(Operating System)のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、ＲＡＭ１０２には、ＣＰＵ１０１による処理に必要な各種データが格納される。ＨＤＤ１０３は、ＯＳやアプリケーションプログラムが格納される。
【０１１７】
グラフィック処理装置１０４には、モニタ１１が接続されている。グラフィック処理装置１０４は、ＣＰＵ１０１からの命令に従って、画像をモニタ１１の画面に表示させる。入力インタフェース１０５には、キーボード１２とマウス１３とが接続されている。入力インタフェース１０５は、キーボード１２やマウス１３から送られてくる信号を、バス１０７を介してＣＰＵ１０１に送信する。
【０１１８】
通信インタフェース１０６は、ネットワーク１０に接続されている。ネットワーク１０は、たとえばインターネットのような広域ネットワークである。通信インタフェース１０６は、ネットワーク１０を介して、他のコンピュータとの間でデータの送受信を行う。
【０１１９】
以上のようなハードウェア構成によって、本実施の形態の処理機能を実現することができる。たとえば、電磁波解析プログラムをコンピュータ１００に実行させることにより、コンピュータ１００を電磁波解析装置として機能させることができる。
【０１２０】
図５は、電磁波解析装置として動作するコンピュータの機能ブロック図である。コンピュータ１００には、入力装置１２ａと出力装置１１ａとが接続されている。
【０１２１】
入力装置１２ａは、キーボード１２やマウス１３などのデータを入力するための機器である。ユーザは、入力装置１２ａを介して、電磁波解析の初期値を入力することができる。
【０１２２】
出力装置１１ａは、モニタ１１などのデータ出力機器である。コンピュータ１００で計算された電磁界の過渡的な挙動を示す情報は、出力装置１１ａを介して出力される。
【０１２３】
コンピュータ１００は、初期値設定部１１１，電磁界データ記憶部１１２、セルサイズ判定部１１３、および電磁界計算部１１４で構成される。
初期値設定部１１１は、入力装置１２ａから入力された電磁波解析の初期値を、電磁界データ記憶部１１２に設定する。初期値としては、解析領域の定義（セルサイズ、誘電率、透磁率、電気伝導率など）や、電磁界の初期状態などが設定される。
【０１２４】
電磁界データ記憶部１１２は、たとえば、コンピュータ１００のＲＡＭ１０２内の記憶領域である。電磁界データ記憶部１１２は、初期値設定部１１１によって設定された初期値や、電磁界計算部１１４によって計算された電磁界の過渡的な挙動を示す情報を保持する。
【０１２５】
セルサイズ判定部１１３は、電磁界データ記憶部１１２に設定されたセルサイズの情報を取得し、セルサイズ境界を判定する。セルサイズ判定部１１３は、セルサイズ境界の場所を、電磁界計算部１１４に通知する。
【０１２６】
電磁界計算部１１４は、電磁界データ記憶部１１２に設定された電磁界データに基づいて、電磁界の過渡的な挙動を計算する。電磁界計算部１１４は、計算結果を、逐一電磁界データ記憶部１１２に格納する。また、所定の時間分の計算が終了すると、電磁界計算部１１４は、電磁界の過渡的な挙動を示す情報を、解析結果として出力装置１１ａに対して出力する。解析結果は、時間刻み幅毎の各時刻における所定の計算位置での磁界や電界の値である。
【０１２７】
図６は、本実施の形態による電磁波解析処理のフローチャートである。以下に、図６に示す処理をステップ番号に沿って説明する。
［ステップＳ１１］初期値設定部１１１は、入力装置１２ａからの入力に応じて初期設定を行う。初期設定の内容は、セル数、セルサイズ、各セルの材料定数（誘電率ε、透磁率μ、電気伝導率σなど）、境界条件（吸収境界条件、電気壁境界条件など）、励振条件（波源の位置、波形など）、解析終了時間、時間刻み幅などの設定ならびに電磁界計算のための係数の算出である。
【０１２８】
なお、このとき、電磁波解析シミュレーションの時刻Ｔを０に設定する。
［ステップＳ１２］セルサイズ判定部１１３は、初期設定で設定された各セルに関して、セルサイズを判定する。セルサイズの判定により、両側のセルサイズが等しい電界配置地点（電界の計算位置）と、両側のセルサイズが異なる電界配置地点とが区別される。
【０１２９】
［ステップＳ１３］電磁界計算部１１４は、両側のセルサイズが等しい電界配置地点における磁界の偏微分係数を、式(１８)によって計算する。
［ステップＳ１４］電磁界計算部１１４は、両側のセルサイズが異なる電界配置地点における磁界の偏微分係数を、式(１９)によって計算する。
【０１３０】
［ステップＳ１５］電磁界計算部１１４は、ステップＳ１３，Ｓ１４で計算した磁界の偏微分係数を用いて、式(１７)によって電界の計算を行う。
［ステップＳ１６］電磁界計算部１１４は、境界条件に依存した電界を、別途その条件式に基づいて計算を行う。その後、電磁界計算部１１４は、時刻ＴをΔＴ／２だけ進める。
【０１３１】
［ステップＳ１７］電磁界計算部１１４は、式（１５）および（１６）による磁界の計算を行う。
［ステップＳ１８］電磁界計算部１１４は、境界条件に依存した磁界を、別途その条件式に基づいて計算を行う。その後、電磁界計算部１１４は、時刻ＴをΔＴ／２だけ進める。
【０１３２】
［ステップＳ１９］電磁界計算部１１４は、時刻Ｔが、あらかじめ設定した解析終了時刻Ｔmaxを超えたか否か判断する。時刻Ｔが解析終了時刻Ｔmaxを超えた場合には、処理がステップＳ２０に進められる。時刻Ｔが解析終了時刻Ｔmaxを超えていない場合には、処理がステップＳ１３に進められる。
【０１３３】
［ステップＳ２０］電磁界計算部１１４は、計算結果を出力装置１１ａへ出力する。
以上のような処理により、様々な環境での電磁波の挙動を計算することができる。以下に、電磁波解析の計算事例について説明する。
【０１３４】
［平面波の伝播］
まず、平面波の伝搬についての計算事例を説明する。
図７は、平面波の伝搬を計算するモデルを示す図である。図７では、右方向がｘ軸の正、上方向がｙ軸の正、紙面手前方向がｚ軸の正である。ここでは、２次元ＴＥ波を扱うものとする。図７において格子状に表されたセルの辺上に電界Ｅｘ、Ｅｙが配置され、セルの中心に磁界Ｈｚが配置されている。
【０１３５】
解析領域３０は、セルサイズの大きな領域３１とセルサイズの小さな領域３２とで構成される。領域３１のセルサイズはΔｘC×Δｙ、領域３２のセルサイズはΔｘF×Δｙである。セルサイズの異なる２つの領域３１，３２は、セルサイズ境界３３で接している。
【０１３６】
Ｙ方向の終端には電気壁境界条件を設定し、Ｘ方向の終端はセルサイズ境界から十分遠方になるようにしてある。これは、Ｘ方向の終端からの反射の影響を除去するためである。
【０１３７】
セルサイズ境界は、電磁波解析の計算のために解析領域に便宜的に設けられたものであり、解析対象自身の特性ではない。そのため、解析領域に電磁波を伝搬させた場合、計算上、セルサイズ境界からの反射波が検出されないことが理想である。
【０１３８】
そこで、セルサイズ境界から−Ｘ方向に十分に離れた地点において、電界Ｅｙをガウス波で励振し、パルスを＋Ｘ方向へ伝搬させ、セルサイズ境界における反射量を計算した。Δｙ＝Δｘc＝1.0mm、ΔｘF＝0.1 mmの場合（Δｘc／ΔｘF＝１０）の反射特性は、以下の通りである。
【０１３９】
図８は、セルサイズ境界における反射特性を示す図である。図８では、横軸に周波数（ＧＨｚ）を示し、縦軸に反射量（ｄＢ）を示している。反射量の数値は、−１２０ｄＢに近いほどセルサイズ境界からの反射が少ないことを意味している。従来手法の反射特性は実線４１で表し、本発明の実施の形態に基づく手法の反射特性は、点線４２で表している。
【０１４０】
図８からも分かるように、本実施の形態によれば、従来手法に比べ反射量が２０〜４０ｄＢ低減できていることがわかる。これは、セルサイズ境界上の電界Ｅｙを計算する際に、打ち切り誤差が低減され領域の不連続性が緩和されたからであると考えられる。
【０１４１】
なお、この特性は、Δｙ＝Δｘc＝1.0mm、ΔｘF＝0.025mmの場合（Δｘc／ΔｘF＝４０）でもほとんど変わらない。
［パッチアンテナの反射特性］
次に、パッチアンテナの反射特性の解析結果について説明する。
【０１４２】
図９は、パッチアンテナの平面図である。図９の例では、解析領域５０内にパッチアンテナ５１が配置されている。
解析領域５０は、縦横118.0mmの正方形の領域である。パッチアンテナ５１の先端部は、縦横共に11.8mmの大きさであり、解析領域５０の４方の縁から50.0mm離れている。
【０１４３】
パッチアンテナ５１は、幅1.8mmの信号線５２で外部に接続される。パッチアンテナ５１は、信号線５２との接続部において、信号線５２の両側に深さ4.0mm、幅1.0mmの溝が形成されている。
【０１４４】
なお、パッチアンテナ５１が形成されている基板の厚さは0.8mm、比誘電率は3.274であるものとする。
図９に示したパッチアンテナ５１は、社団法人電子情報通信学会マイクロ波シミュレータ研究会において規範問題とされているマイクロストリップ給電パッチアンテナである。たとえば、「並木武文, 坂口拓史, 伊藤公一, “FDTD法を用いたパッチアンテナ解析における計算精度についての一考察,”信学技報, AP99-12, pp.17-22, May 1999.」や「田口光雄, “アンテナ設計から見た電磁界シミュレータの使い方と評価,” 信学誌, vol.83, no.11, pp.878-883, Nov. 2000.」において、パッチアンテナ５１を用いたシミュレーション例が詳しく説明されている。
【０１４５】
本実施の形態では、このパッチアンテナ５１を計算するために、不均一セルを使って領域分割し、計算モデルを作成する。
図１０は、パッチアンテナの計算モデルの一例を示す図である。図中、右方向がｘ軸の正、上方向がｙ軸の正、紙面手前方向がｚ軸の正である。
【０１４６】
ここで、放射パッチであるパッチアンテナ５１の端部近傍は、電磁界の特異性を考慮して非常に微細なセルに分割してある。セルの分割手法に関しては、「E.M.Daniel and C.J.Railton, "Fast finite difference time domain analysis of microstrip patch antennas,"IEEE AP-S Digest, pp.414-417, 1991.」や「大西輝夫, 柏達也, 内藤行雄, 細矢良雄, “同軸給電パッチアンテナのFD-TD解析における効率化の一検討,”信学総大, B-1-128, 1997.」で説明されている。
【０１４７】
Ｚ方向についても同様に、パッチアンテナ５１の近傍は非常に微細なセルに分割している。図１０の計算モデルにおける全セル数は76380個である。
図１１は、図１０の計算モデルにおける最小ならびに最大セルサイズを示す図である。ｘ軸方向のセルサイズΔＸの最小値(Minimum cell)は、0.200mmである。ｘ軸方向のセルサイズΔＸの最大値(Maximum cell)は、2.000mmである。ｙ軸方向のセルサイズΔＹの最小値(Minimum cell)は、0.050mmである。ｙ軸方向のセルサイズΔＹの最大値(Maximum cell)は、2.000mmである。ｚ軸方向のセルサイズΔＺの最小値(Minimum cell)は、0.069mmである。ｚ軸方向のセルサイズΔＺの最大値(Maximum cell)は、1.200mmである。
【０１４８】
また、解析領域の終端にはMurの１次の吸収境界条件を設定した。Murの１次の吸収境界条件に関しては、「G.Mur,"Absorbing Boundary Conditions for the finite-difference approximation of the time-domain electromagnetic field equations,"IEEE Trans.EMC-23, pp.377-382, 1981.」で説明されている通りである。
【０１４９】
このようなパッチアンテナ５１のサンプルを用意し、そのサンプルを用いて反射特性を測定すると共に、本実施の形態における解析手法と従来の解析手法とにより反射特性を計算した。
【０１５０】
図１２は、図１０に示す計算モデルの反射特性の計算結果と測定結果とを示す図である。図１２では、横軸に周波数（ＧＨｚ）を示し、縦軸に反射量（ｄＢ）を示している。サンプルによる測定結果を破線６１で表し、本発明の実施の形態に基づく手法による計算結果を点線６２で表し、従来手法による計算結果を実線６３で表している。
【０１５１】
図１２からも分かるように、本発明の実施の形態に基づく手法で電磁波解析を行えば、従来に比べて、測定結果に近い計算結果を得ることができる。なお、測定値の信頼性については、十分に検討してある。測定値の信頼性の検討内容は、「並木武文, 坂口拓史, 伊藤公一, “FDTD法を用いたパッチアンテナ解析における計算精度についての一考察,”信学技報, AP99-12, pp.17-22, May 1999.」（以下、この文献を「計算精度考察例」と呼ぶ）に説明されているものと同様である。
【０１５２】
ここで、図１２に示した反射特性から共振周波数を抽出し、測定値に対する計算値の相対誤差を算出した。
図１３は、共振周波数と相対誤差を示す図である。測定値に基づく共振周波数(Resonant frequency)は、7.0100ＧＨｚである。従来の手法の計算値に基づく共振周波数は6.8624ＧＨｚであり、相対誤差(Relative error)は2.10％である。本発明の実施の形態に係る手法の計算値に基づく共振周波数は6.9756ＧＨｚであり、相対誤差は0.49%である。
【０１５３】
図１２、図１３から明らかなように、図１０の計算モデルを用いた場合、本手法によって大幅な精度改善が図れることがわかる。
今回の計算で必要としたＣＰＵ時間とメモリサイズを以下に示す。
【０１５４】
図１４は、計算に使用されたＣＰＵ時間とメモリサイズとを示す図である。従来の手法で使用したＣＰＵ時間は493秒であり、使用したメモリサイズは5.5ＭＢである。本発明の実施の形態の手法で使用したＣＰＵ時間は666秒であり、使用したメモリサイズは6.4ＭＢである。
【０１５５】
なお、計算には、ワークステーション(CPU：UltraSPARC II 360MHz)を使用した。本発明による手法の場合、通常のＦＤＴＤ法と比較して、ＣＰＵ時間は約1.35倍、メモリは1.16倍増加している。しかし、その精度改善効果からすれば問題となる程度ではないと考えられる。
【０１５６】
前記の計算精度考察例における通常のＦＤＴＤ法を用いた同様の検討では、不均一セルを使用した場合は、共振周波数を相対誤差1.0%以下で算出することが出来なかった。そして、共振周波数を相対誤差1.0%以下で算出するためには、莫大な数の微細な均一セルを使用しなければならなかった。この場合、ほぼ同性能のワークステーションを使用して、ＣＰＵ時間が600分以上、メモリが250ＭＢ必要であった。これらの結果を鑑みれば、本発明による手法の有効性は明らかである。
【０１５７】
なお、上記の処理機能は、コンピュータによって実現することができる。その場合、電磁波解析装置が有すべき機能の処理内容を記述したプログラムが提供される。そのプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリなどがある。磁気記録装置には、ハードディスク装置（ＨＤＤ）、フレキシブルディスク（ＦＤ）、磁気テープなどがある。光ディスクには、ＤＶＤ(Digital Versatile Disc)、ＤＶＤ−ＲＡＭ(Random Access Memory)、ＣＤ−ＲＯＭ(Compact Disc Read Only Memory)、ＣＤ−Ｒ(Recordable)／ＲＷ(ReWritable)などがある。光磁気記録媒体には、ＭＯ(Magneto-Optical disc)などがある。
【０１５８】
プログラムを流通させる場合には、たとえば、そのプログラムが記録されたＤＶＤ、ＣＤ−ＲＯＭなどの可搬型記録媒体が販売される。また、プログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することもできる。
【０１５９】
プログラムを実行するコンピュータは、たとえば、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、自己の記憶装置に格納する。そして、コンピュータは、自己の記憶装置からプログラムを読み取り、プログラムに従った処理を実行する。なお、コンピュータは、可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することもできる。また、コンピュータは、サーバコンピュータからプログラムが転送される毎に、逐次、受け取ったプログラムに従った処理を実行することもできる。
【０１６０】
（付記１） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算する電磁波解析装置において、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断するセルサイズ判断手段と、
前記セルサイズ判断手段において周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算する第１の計算手段と、
前記セルサイズ判断手段において周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算する第２の計算手段と、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する出力手段と、
を有することを特徴とする電磁波解析装置。
【０１６１】
（付記２） 前記第２の計算手段は、前記第２の計算手法として、陰的スキームを用いることを特徴とする付記１記載の電磁波解析装置。
（付記３） 前記第２の計算手段は、前記陰的スキームにより、周囲のセルのサイズが異なる計算位置の電磁界の偏微分係数を算出することを特徴とする付記２記載の電磁波解析装置。
【０１６２】
（付記４） 前記第１の計算手段は、周囲のセルのサイズが同じである計算位置の電磁界の偏微分係数を計算し、
前記第２の計算手段は、周囲のセルのサイズが異なる計算位置の周囲の電磁界の微分係数を前記第１の計算手段から取得し、前記第１の計算手段で計算された電磁界の微分係数を用いて、前記陰的スキームの解を求めることを特徴とする付記３記載の電磁波解析装置。
【０１６３】
（付記５） 前記セルサイズ判断手段は、解析領域に設定された座標軸毎に、隣接するセルの各軸方向の幅の同一性を判断することを特徴とする付記１記載の電磁波解析装置。
【０１６４】
（付記６） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析プログラムにおいて、
コンピュータに、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算し、
周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算し、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する、
処理を実行させることを特徴とする電磁波解析プログラム。
【０１６５】
（付記７） 前記第２の計算手法として、陰的スキームを用いることを特徴とする付記６記載の電磁波解析プログラム。
（付記８） 第２の計算手法による計算の際には、前記陰的スキームにより、周囲のセルのサイズが異なる計算位置の電磁界の偏微分係数を算出することを特徴とする付記７記載の電磁波解析プログラム。
【０１６６】
（付記９） 前記第１の計算手法による計算の際には、周囲のセルのサイズが同じである計算位置の電磁界の偏微分係数を計算し、
前記第２の計算手段は、周囲のセルのサイズが異なる計算位置の周囲の電磁界の微分係数を前記第１の計算手段から取得し、前記第１の計算手段で計算された電磁界の微分係数を用いて、前記陰的スキームの解を求めることを特徴とする付記８記載の電磁波解析プログラム。
【０１６７】
（付記１０） 前記セルサイズ判断手段は、解析領域に設定された座標軸毎に、隣接するセルの各軸方向の幅の同一性を判断することを特徴とする付記６記載の電磁波解析プログラム。
【０１６８】
（付記１１） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体において、
前記コンピュータに、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算し、
周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算し、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する、
処理を実行させることを特徴とする記録媒体。
【０１６９】
（付記１２） 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析方法において、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
周囲のセルのサイズが同じであると判断された計算位置における電磁界を、第１の計算手法で計算し、
周囲のセルのサイズが異なると判断された計算位置における電磁界を、前記第１の計算手法よりも少ない誤差で計算可能な第２の計算手法で計算し、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する、
ことを特徴とする電磁波解析方法。
【０１７０】
【発明の効果】
以上説明したように本発明では、周囲のセルのサイズが異なる計算位置における電磁界を、他の計算位置よりも少ない誤差の計算手法で計算するようにしたため、セルのサイズが切り替わるセルサイズ境界における計算誤差を減少させることができ、むやみに計算量を増大させることなく、誤差の少ない電磁波解析を行うことができる。
【図面の簡単な説明】
【図１】本発明の原理構成図である。
【図２】均一セルにおける電磁界計算位置の配置図である。
【図３】不均一セルにおける電磁界計算位置の配置図である
【図４】本発明の実施の形態に用いるコンピュータのハードウェア構成例を示す図である。
【図５】電磁波解析装置として動作するコンピュータの機能ブロック図である。
【図６】本実施の形態による電磁波解析処理のフローチャートである。
【図７】平面波の伝搬を計算するモデルを示す図である。
【図８】セルサイズ境界における反射特性を示す図である。
【図９】パッチアンテナの平面図である。
【図１０】パッチアンテナの計算モデルの一例を示す図である。
【図１１】図１０の計算モデルにおける最小ならびに最大セルサイズを示す図である。
【図１２】図１０に示す計算モデルの反射特性の計算結果と測定結果とを示す図である。
【図１３】共振周波数と相対誤差を示す図である。
【図１４】計算に使用されたＣＰＵ時間とメモリサイズとを示す図である。
【符号の説明】
１ 解析領域
２ セルサイズ判断手段
３ 第１の計算手段
４ 第２の計算手段
５ 出力手段
１０ ネットワーク
１１ モニタ
１２ キーボード
１３ マウス
１００ コンピュータ
１０１ ＣＰＵ
１０２ ＲＡＭ
１０３ ハードディスクドライブ
１０４ グラフィック処理装置
１０５ 入力インタフェース
１０６ 通信インタフェース
１０７ バス
Claims (5)
- 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算する電磁波解析装置において、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断するセルサイズ判断手段と、
シミュレーション時刻が過去の電磁界に基づいて計算対象時刻の電磁界を計算する第１の計算手法を用い、前記セルサイズ判断手段において周囲のセルのサイズが同じであると判断された計算位置における前記計算対象時刻の電磁界を計算する第１の計算手段と、
陰解法に基づいた計算手法である陰的スキームを利用して電磁界を計算する第２の計算手法を用い、周囲のセルの前記計算対象時刻の電磁界については前記第１の計算手段の計算結果を適用することにより、前記セルサイズ判断手段において周囲のセルのサイズが異なると判断された計算位置における前記計算対象時刻の電磁界を計算する第２の計算手段と、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する出力手段と、
を有することを特徴とする電磁波解析装置。 - 前記第２の計算手段は、周囲のセルの前記計算対象時刻の磁界の計算結果を前記第１の計算手段から取得し、取得した磁界の計算結果に基づいて前記計算対象時刻の前記計算位置の磁界の偏微分係数を計算し、シミュレーション時刻が過去の電磁界と前記計算位置の磁界の偏微分係数とから前記計算位置における前記計算対象時刻の電磁界を計算することを特徴とする請求項１記載の電磁波解析装置。
- 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算する電磁波解析装置において、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断するセルサイズ判断手段と、
時間刻み幅がΔｔであり、シミュレーションのステップ数ｎまでの電磁界の過渡的変化が既に求められている状況下でステップ数（ｎ＋１）における電磁界を計算するときに、前記セルサイズ判断手段において周囲のセルのサイズが同じであると判断された計算位置について、シミュレーション時刻（ｎ−１／２）×Δｔにおける磁界とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界を計算し、算出したシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界からシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の偏微分係数を計算し、算出したシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の偏微分係数とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１）×Δｔにおける電界を計算する第１の計算手段と、
前記セルサイズ判断手段において周囲のセルのサイズが異なると判断された計算位置について、シミュレーション時刻（ｎ−１／２）×Δｔにおける磁界とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界を計算し、算出したシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界と前記第１の計算手法で算出された周囲のセルのシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の微分係数とから計算位置のシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の微分係数を計算し、算出した計算位置の磁界の微分係数とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１）×Δｔにおける電界を計算する第２の計算手段と、
前記第１の計算手段と前記第２の計算手段とにより計算された前記各計算位置における電磁界の値を出力する出力手段と、
を有することを特徴とする電磁波解析装置。 - 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析プログラムにおいて、
コンピュータに、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
シミュレーション時刻が過去の電磁界に基づいて計算対象時刻の電磁界を計算する第１の計算手法を用い、周囲のセルのサイズが同じであると判断された計算位置における前記計算対象時刻の電磁界を計算し、
陰解法に基づいた計算手法である陰的スキームを利用して電磁界を計算する第２の計算手法を用い、周囲のセルの前記計算対象時刻の電磁界については前記第１の計算手法による計算結果を適用することにより、周囲のセルのサイズが異なると判断された計算位置における前記計算対象時刻の電磁界を計算し、
計算された前記各計算位置における電磁界の値を出力する、
処理を実行させることを特徴とする電磁波解析プログラム。 - 解析領域を複数のセルに分割し、セル毎の電磁界の過渡的な挙動を数値計算するための電磁波解析プログラムにおいて、
コンピュータに、
前記解析領域に設定された各計算位置に関し、周囲のセルのサイズの同一性を判断し、
時間刻み幅がΔｔであり、シミュレーションのステップ数ｎまでの電磁界の過渡的変化が既に求められている状況下でステップ数（ｎ＋１）における電磁界を計算するときに、周囲のセルのサイズが同じであると判断された計算位置について、シミュレーション時刻（ｎ−１／２）×Δｔにおける磁界とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界を計算し、算出したシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界からシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の偏微分係数を計算し、算出したシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の偏微分係数とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１）×Δｔにおける電界を計算し、
周囲のセルのサイズが異なると判断された計算位置について、シミュレーション時刻（ｎ−１／２）×Δｔにおける磁界とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界を計算し、算出したシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界と電界計算時に算出された周囲のセルのシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の微分係数とから前記計算位置のシミュレーション時刻（ｎ＋１／２）×Δｔにおける磁界の微分係数を計算し、算出した前記計算位置の磁界の微分係数とシミュレーション時刻ｎ×Δｔにおける電界とからシミュレーション時刻（ｎ＋１）×Δｔにおける電界を計算し、
計算された前記各計算位置における電磁界の値を出力する、
処理を実行させることを特徴とする電磁波解析プログラム。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2001193224A JP3818874B2 (ja) | 2001-06-26 | 2001-06-26 | 電磁波解析装置および電磁波解析プログラム |
| US10/022,539 US6662125B2 (en) | 2001-06-26 | 2001-12-20 | Electromagnetic wave analyzer and program for same |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2001193224A JP3818874B2 (ja) | 2001-06-26 | 2001-06-26 | 電磁波解析装置および電磁波解析プログラム |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2003006181A JP2003006181A (ja) | 2003-01-10 |
| JP3818874B2 true JP3818874B2 (ja) | 2006-09-06 |
Family
ID=19031554
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2001193224A Expired - Lifetime JP3818874B2 (ja) | 2001-06-26 | 2001-06-26 | 電磁波解析装置および電磁波解析プログラム |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US6662125B2 (ja) |
| JP (1) | JP3818874B2 (ja) |
Families Citing this family (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4720964B2 (ja) * | 2001-05-31 | 2011-07-13 | 日本電気株式会社 | Ｆｅｍ解析方法、プログラム、およびシステム |
| JP2004320885A (ja) * | 2003-04-16 | 2004-11-11 | Hitachi Ltd | 回転機の磁界解析方法及びプログラム |
| WO2005057434A1 (ja) * | 2003-12-10 | 2005-06-23 | Murata Manufacturing Co., Ltd. | 電磁界解析装置、電磁界解析プログラムおよびそのプログラムを記録した記録媒体 |
| EP1617309B1 (en) * | 2004-07-15 | 2011-01-12 | Fujitsu Limited | Simulation technique with local grid refinement |
| JP4520822B2 (ja) * | 2004-10-29 | 2010-08-11 | 富士通株式会社 | 電磁波解析装置、電磁波解析方法および電磁波解析プログラム |
| JP2006209590A (ja) | 2005-01-31 | 2006-08-10 | Ricoh Co Ltd | 電磁界解析装置および解析方法、ならびに解析プログラム |
| JP4262762B2 (ja) * | 2006-12-11 | 2009-05-13 | シャープ株式会社 | 電磁界解析プログラム |
| GB0723224D0 (en) * | 2007-11-27 | 2008-01-09 | Fujitsu Ltd | A multi-level gris embedding process with uniform and non-uniform grid refinement for multigurid-fdtd electromagnetic solver |
| GB0723222D0 (en) * | 2007-11-27 | 2008-01-09 | Fujitsu Ltd | A very stable multigrid fdtd solver |
| JP2009294386A (ja) * | 2008-06-04 | 2009-12-17 | Toshiba Corp | リソグラフィシミュレーション方法 |
| JP5149321B2 (ja) * | 2010-03-24 | 2013-02-20 | 株式会社東芝 | 電磁場シミュレーション方法、電磁場シミュレーション装置、半導体装置の製造方法 |
| US9416642B2 (en) | 2013-02-01 | 2016-08-16 | Halliburton Energy Services, Inc. | Modeling subterranean rock blocks in an injection treatment simulation |
| US9798042B2 (en) | 2013-02-01 | 2017-10-24 | Halliburton Energy Services, Inc. | Simulating an injection treatment of a subterranean zone |
| US9239407B2 (en) | 2013-08-27 | 2016-01-19 | Halliburton Energy Services, Inc. | Injection treatment simulation using condensation |
| CN107918696B (zh) * | 2017-10-20 | 2021-05-11 | 西安电子科技大学 | 相控阵天线的多场耦合分析方法及计算机程序 |
| JP7121111B2 (ja) * | 2018-03-07 | 2022-08-17 | 住友重機械工業株式会社 | 磁場解析装置、解析方法、及びプログラム |
| JP2021110997A (ja) * | 2020-01-07 | 2021-08-02 | 住友重機械工業株式会社 | シミュレーション方法、シミュレーション装置、及びプログラム |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO1994003413A1 (en) * | 1992-08-10 | 1994-02-17 | New Mexico State University Technology Transfer Corporation | A digitally-configurable analog vlsi chip and method for real-time solution of partial differential equations |
| US5880838A (en) * | 1996-06-05 | 1999-03-09 | California Institute Of California | System and method for optically measuring a structure |
| US6226599B1 (en) | 1997-03-05 | 2001-05-01 | Fujitsu Limted | Electromagnetic wave analyzer apparatus |
| JP3633765B2 (ja) | 1997-11-19 | 2005-03-30 | 富士通株式会社 | シミュレーション装置及びシミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体 |
| JP3639116B2 (ja) | 1998-06-12 | 2005-04-20 | 富士通株式会社 | 電磁波解析装置及び電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体 |
| US6665849B2 (en) * | 1999-06-09 | 2003-12-16 | Interuniversitair Microelektronica Centrum Vzw | Method and apparatus for simulating physical fields |
-
2001
- 2001-06-26 JP JP2001193224A patent/JP3818874B2/ja not_active Expired - Lifetime
- 2001-12-20 US US10/022,539 patent/US6662125B2/en not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| US6662125B2 (en) | 2003-12-09 |
| JP2003006181A (ja) | 2003-01-10 |
| US20030050760A1 (en) | 2003-03-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3818874B2 (ja) | 電磁波解析装置および電磁波解析プログラム | |
| Lee et al. | A 3-D spectral-element method using mixed-order curl conforming vector basis functions for electromagnetic fields | |
| Bracken et al. | S-domain methods for simultaneous time and frequency characterization of electromagnetic devices | |
| US20110270593A1 (en) | Method and apparatus for simulating physical fields | |
| JP5332731B2 (ja) | 電磁界シミュレータ及び電磁界シミュレーション装置 | |
| Kumar et al. | Analytical modeling and performance benchmarking of on-chip interconnects with rough surfaces | |
| JP3639116B2 (ja) | 電磁波解析装置及び電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体 | |
| Dadash et al. | Analytical adjoint sensitivity formula for the scattering parameters of metallic structures | |
| JP2019164102A (ja) | 電磁波測定点算出装置及び放射妨害波測定装置 | |
| Wu et al. | A new method of modeling three-dimensional MIC/MMIC circuits: The space-spectral domain approach | |
| JPH11295365A (ja) | 電磁界解析方法 | |
| Jithesh et al. | A review on computational EMI modelling techniques | |
| Demeester et al. | Construction of the Dirichlet to Neumann boundary operator for triangles and applications in the analysis of polygonal conductors | |
| Christopoulos et al. | The application of transmission-line modeling (TLM) to electromagnetic compatibility problems | |
| JP3703812B2 (ja) | Ｆｄｔｄ法を用いた電磁界解析方法、電磁界解析における媒質表現方法、シミュレーション装置、及びプログラム | |
| Chen et al. | A model for the prediction of the shielding effectiveness of cylindrical enclosure | |
| Jedlicka et al. | Coupling through tortuous path narrow slot apertures into complex cavities | |
| US20230243878A1 (en) | Computer-readable recording medium storing program, information processing apparatus, and obtainment method of conductor loss and dielectric loss | |
| Zhou et al. | Application of subspace projection approaches for reduced-order modeling of electromagnetic systems | |
| Nilavalan et al. | Quantifying numerical dispersion in non-orthogonal FDTD meshes | |
| Becker et al. | Finite difference time domain modeling of high speed electronic packages | |
| Ndagijimana et al. | The Use of a Full-Wave 3D Technique in a CAD Oriented Analysis: The TLM Method | |
| Chebolu et al. | The analysis of microwave antennas using the FDTD method | |
| US11501049B1 (en) | Systems and methods for modeling interactions of power and signals in a multi-layered electronic structure | |
| Grayaa et al. | Full-wave analysis of shielded planar circuits using different models of sources |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20050524 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20050715 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060124 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060324 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20060613 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20060613 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 Ref document number: 3818874 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100623 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110623 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120623 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120623 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130623 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140623 Year of fee payment: 8 |
|
| EXPY | Cancellation because of completion of term |