JP2004004054A - Method for analyzing electromagnetic field using fdtd method, method for representing medium in analysis of electromagnetic field, simulation system, and program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a versatile technology for analyzing an electromagnetic field surely with high accuracy by representing an object being analyzed (medium) existing in an analytic space more accurately at all times. <P>SOLUTION: A conductor is represented by a line and the face of a cell having a side intersecting that conductor is divided by that conductor. The divided face including a position for determining the magnetic field strength is handled as one face as usual and the divided face not including that position is added to the adjacent face. Two electric field strengths are defined on the side intersecting the conductor. Magnetic field strength of a face BCDE is thereby determined using an integration path ABCDEFA and magnetic field strength of a face BEGH is determined using an integration path AFGHA. <P>COPYRIGHT: (C)2004,JPO

Description

【０００９】
ここで（１）式（アンペアの法則）、（２）式（ファラデーの電磁誘導の法則）、（３）式（磁界中のガウスの法則）、（４）式（電界中のガウスの法則）は独立な式ではなく、ＦＤＴＤ法において（３）、（４）式は数値誤差の評価基準として用いられる程度で、定式化には（１）、（２）式が用いられる。（１）、（２）式にＹｅｅのアルゴリズムを用いて、時空間と空間座標について中心差分法を適用し、定式化すると、３次元空間において次のようなＦＤＴＤ基本方程式を得ることができる。
【００１０】
【数２】

【００１１】
ただし、μ：透磁率、ε：誘電率、σ：抵抗率である。
（５）、（６）式中の「Ｅ」や「Ｈ」の各シンボルに下添字として付した「ｘ」や「ｙ」は、それぞれ、電界強度、磁界強度の方向を表している。ここでは、ｘ方向に割り当てた電界強度、磁界強度のみを示したが、ｙ、ｚ方向についても同様に導出できる。
【００１２】
【非特許文献１】
羽野、「ＦＤＴＤ（Ｆｉｎｉｔｅ　Ｄｉｆｆｅｒｅｎｃｅ　Ｔｉｍｅ　Ｄｏｍａｉｎ　）法」、平成１４年電気学会全国大会論文集、第５冊、ｐｐ４１１−４１４、２００２
【非特許文献２】
野田、横山、「ＦＤＴＤ法に基づく汎用サージ解析プログラムの開発」、電学論Ｂ、１２１巻５号、ｐｐ６２５−６３２、２００１
【非特許文献３】
Ｊ．Ｆａｎｇ　ａｎｄ　Ｊ．Ｒｅｎ、“Ａ　Ｌｏｃａｌｌｙ　Ｃｏｎｆｏｒｍｅｄ　Ｆｉｎｉｔｅ−Ｄｉｆｆｅｒｅｎｃｅ　Ｔｉｍｅ−Ｄｏｍａｉｎ　Ａｌｇｏｒｉｔｈｍ　ｏｆ　Ｍｏｄｅｌｉｎｇ　Ａｒｂｉｔｒａｒｙ　Ｓｈａｐｅ　Ｐｌａｎａｒ　Ｍｅｔａｌ　Ｓｔｒｉｐｓ　”、ＩＥＥＥ　ＴＲＡＮＳＡＣＴＩＯＮＳ　ＯＮ　ＭＩＣＲＯＷＡＶＥ　ＴＨＥＯＲＹ　ＡＮＤ　ＴＥＣＨＮＩＱＵＥＳ、ＶＯＬ．４１、ＮＯ．５、ＭＡＹ　１９９３
【００１３】
【発明が解決しようとする課題】
ＦＤＴＤ法では、図１に示すようなセルで解析空間を分割するのを基本とする。このため、そのセルの格子（辺）上に存在しない導体などの解析対象物体を表現する場合には、その形状を階段状に近似する階段近似法が用いられていた。しかし、階段近似法で表現される形状と実際の形状との間には大きな違いがあるのが普通である。このため、階段近似法を用いて行う従来の解析方法では、解析対象物体の共振周波数やその伝搬時間などに大きな誤差が生じて、解析を高精度に行えないという問題点があった。
【００１４】
セルサイズを小さくすることにより、階段近似法で表現される形状と実際の形状との間の違いを小さくすることができる。しかし、そのサイズを小さくすると、それに伴いセル数が増大することから、必要な計算機資源や計算時間が膨大となってしまう。このため、部分的であってもセルサイズを小さくすることは望ましくない。サージ解析の分野では、比較的に大きな解析空間を対象にすることが多いことから、言い換えればセル数が大きいのが普通であることから、特に望ましくないと言える。
【００１５】
セルの形状を変形させて解析対象物体の形状を表現する場合には、その物体に応じたモデリングが必要になり、一般性に欠ける。その物体が格子上となるように座標系を設定する方法は、複数の物体が任意の状態で解析空間内に存在していれば用いることができないのが普通である。このようなことから、セルサイズを小さくする、セルの形状を変形させるといった操作を行うことなく、解析空間内に存在する解析対象物体をより正確に表現して解析を行えるようにすることが重要であると考えられる。
【００１６】
本発明は、解析空間内に存在する解析対象物体（媒質）を常により正確に表現できるようにして、電磁界解析を確実に高精度に行えるようにするための汎用的な技術を提供することを目的とする。
【００１７】
【課題を解決するための手段】
本発明のＦＤＴＤ法を用いた電磁界解析方法は、解析空間を６面体のセルに分割し、ＦＤＴＤ（Ｆｉｎｉｔｅ　Ｄｉｆｆｅｒｅｎｃｅ　Ｔｉｍｅ　Ｄｏｍａｉｎ　）法を用いて電磁界解析を行うための方法であって、セルに分割する解析空間に存在する媒質のなかから選択した媒質を線として設定し、セルを構成する面のなかで線が辺を横切る面を該線で複数に分割して電界強度、或いは磁界強度を求める。
【００１８】
なお、上記線による面の分割は、該面を該線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、該位置を含む分割面は隣接する面に加えない形で行う、ことが望ましい。線で複数に分割される面の電界強度には、該面の周辺で計算した方向が同じ電界強度を採用する、ことが望ましい。線で複数に分割される面の磁界強度は、該線を含む積分路で求めた電界強度を用いて計算する、ことが望ましい、線で複数に分割される面と交差する方向に隣接し、且つ該線が横切る辺を有する面の磁界強度は、複数の分割面と隣接していると想定して計算する、ことが望ましい。線で複数に分割される面の電界強度は、該線を考慮せずに電界強度を全て計算した後、計算した電界強度を補正する形で求める、ことが望ましい。
【００１９】
また、分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面では、該変化を電界強度、或いは磁界強度の計算に反映させる、ことが望ましい。その反映は、変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行う、ことが望ましい。
【００２０】
本発明の電磁界解析における媒質表現方法は、解析空間を６面体のセルに分割し、ＦＤＴＤ法により行う電磁界解析で媒質を表現するために用いられる方法であって、セルに分割する解析空間に存在する媒質を選択し、該選択した媒質を、セルを構成する面の辺を横切る線として設定し表現する。
【００２１】
本発明のシミュレーション装置は、解析空間を６面体のセルに分割し、ＦＤＴＤ法を用いて電磁界解析を行うことを前提とし、解析空間に存在する媒質を線として設定する媒質設定手段と、媒質設定手段により設定した線が解析空間を分割するセルの面のなかで辺を横切る面を該線で分割して電界強度、或いは磁界強度を求める解析手段と、を具備する。
【００２２】
なお、上記線による面の分割は、該面を該線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、該位置を含む分割面は隣接する面に加えない形で行う、ことが望ましい。解析手段は、線で複数に分割される面の電界強度として、該面の周辺で計算した方向が同じ電界強度を採用する、ことが望ましい。
【００２３】
また、解析手段は、分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面では、該変化を電界強度、或いは磁界強度の計算に反映させる、ことが望ましい。その反映は、変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行う、ことが望ましい。
【００２４】
本発明のプログラムは、解析空間を６面体のセルに分割し、ＦＤＴＤ法により電磁界解析を行うシミュレーション装置として用いられるデータ処理装置に実行させることを前提とし、解析空間に存在する媒質を線として設定する機能と、設定する機能により設定した線が解析空間を分割するセルの面のなかで辺を横切る面をその線で分割して電界強度、或いは磁界強度を求める機能と、を実現させる。
【００２５】
本発明では、解析空間に存在する媒質のなかから選択した媒質（導体や構造物など）を線として設定（表現）し、セルを構成する面のなかでその線が辺を横切る面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求める。
【００２６】
その線で面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求めることにより、その線を、セルの面の辺上に配置した形で扱うことになる。その線はセルの面上に配置すれば良いから、高い自由度が得られ、汎用性は向上する。その線で表現する媒質は形状をより正確に模擬することが可能となる。これらのことから、高い汎用性を実現させつつ、電磁界解析は確実に高精度に行えるようになる。
【００２７】
上記線による面の分割を、その面をその線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、その位置を含む分割面は隣接する面に加えない形で行うようにした場合には、その分割に伴い面の数が増大することは回避される。その結果、最終的に求める強度の数、即ち解析結果として保存の対象となるデータの量の増大は回避され、必要な計算機資源の増大は抑えられる。
【００２８】
分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面において、その変化を電界強度、或いは磁界強度の計算に反映させるようにした場合には、その変化に伴って生じる誤差は回避、或いは低減させることが可能となって、その変化を計算に反映させる強度はより高精度に得られるようになる。強度の変化は、隣接する面間で伝搬していく。それにより、近傍の面で求められる強度には、媒質が同じであれば比較的に大きな違いがないのが普通である。このことから、その反映を面積の変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行うようにしたときには、計算する強度は更に高精度に得られるようになる。
【００２９】
【発明の実施の形態】
以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
＜第１の実施の形態＞
図３は、セルを構成する面とその面の辺を横切る導体の位置関係を説明する図である。始めに、図３を参照して、本実施の形態で用いられる導体（他と区別すべき媒質）の表現方法、その導体周辺における磁界と電界の近似方法について詳細に説明する。
【００３０】
本実施の形態では、解析空間は一般的な空間分割手法により６面体のセル（図１参照）に分割し、導体の形状や配置は別に設定する。それにより、セルを構成する面上の任意の位置に任意の角度で傾斜した導体の形状を線で模擬する。セルの格子上に配置できない導体を面上の線で模擬する。
【００３１】
そのように線で模擬した場合、辺を横切る導体は、それと対向する辺を横切るか、或いはそれと隣接する辺を横切ることになる。図３（ａ）は前者、即ち対向する２辺を導体が横切る場合の位置関係を示し、同図（ｂ）は後者、即ち隣接する２辺を導体が横切る場合の位置関係を示している。同図（ｂ）では、面ＡＣＤＥ以外を通る導体は破線で示してある。面上の磁界強度を計算する位置は、黒丸を内部に持つ円で示してある。
【００３２】
従来、ＦＤＴＤ法では、各辺で１つの電界だけを定義する。しかし、本実施の形態では、導体が横切る辺では、それらが交差する点を境に２つ定義する。導体に隣接する磁界の積分路は、その面の辺に沿った積分路から導体を含む積分路に変更する。そのようにして、電界強度や磁界強度を計算する面を導体により分割し、分割した面（分割面）別にそれらの強度を計算する形で格子上以外に存在する導体を考慮する。以下、図３（ａ）、図３（ｂ）に示すケース別に、導体周辺の磁界と電界の近似方法について具体的に説明する。
【００３３】
図３（ａ）に示すケースでは、導体が辺ＢＨ上の点Ａ、及び辺ＥＧ上の点Ｆで交差している。このことから、このケースでは、辺ＢＨは辺ＡＢ、ＡＨ、辺ＥＧは辺ＥＦ、ＦＧにそれぞれ２つに分けて、その中心にあるｙ方向の電界強度Ｅ_ｙを定義する。ここでは、それらの電界強度Ｅ_ｙ　を以下のように定義する。
【００３４】
【数３】

【００３５】
（１０）、（１１）式中の括弧内は、電界強度Ｅｙ　を定義した座標位置を表している。シンボル「Δｙ」はｙ方向上のセルサイズを表している。「ＡＢ」は、点Ａと点Ｂの間の長さを表している。それらは他の式でも同様である。
【００３６】
図３（ｂ）に示すケースでは、導体が辺ＡＣ上の点Ｂ、及び辺ＡＥ上の点Ｇで交差している。このことから、このケースでは、辺ＡＣは辺ＡＢ、ＢＣの２つに分け、その中心にあるｙ方向の電界強度Ｅ_ｙ　を同様に以下のように定義する。
【００３７】
【数４】

【００３８】
一方、導体の近傍に割り当てられたｘ方向の電界強度Ｅ_ｘ　は、以下のように計算する。
図３（ａ）に示すケースでは、積分路ＡＢＣＤＥＦＡと積分路ＢＣＤＥＢを用いて計算した磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）が同一でなければならないという関係から、電界強度Ｅ_ｘ　（ｉ＋１／２，ｊ，ｋ）は次式を用いて計算する。
【００３９】
【数５】

【００４０】
同様に図３（ｂ）に示すケースでは、積分路ＢＣＤＥＦＧＢと積分路ＢＣＤＥＧＢを用いて計算した磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ＋１／２，ｋ）が同一でなければならないという関係から、電界強度Ｅ_ｘ　（ｉ＋１−ＥＧ／２Δｙ，ｊ＋１，ｋ）は次式を用いて計算する。
【００４１】
【数６】

【００４２】
電界強度Ｅ_ｘ　（ｉ＋ＡＧ／２Δｘ，ｊ＋１，ｋ）についても同様の式で計算できる。
（１３）、（１４）式を導出する際に考慮したように、導体の周辺にある磁界強度は、その積分路を導体に沿った形に変更し計算する。図３（ａ）に示すケースでは、磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）を計算する際の積分路は積分路ＡＢＣＤＥＦＡとする。導体の電界強度Ｅはゼロであるから、磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）は次式を用いて計算される。
【００４３】
【数７】

【００４４】
導体が横切る面ＢＥＧＨを構成する面（分割面）ＡＢＥＦはその内部に磁界強度Ｈ_ｚ　を計算する位置を含まない。その分割面ＡＢＥＦで磁界強度Ｈ_ｚ　を計算すると、計算して残すべき磁界強度Ｈの数が面の分割によって増大する。その管理には煩雑な処理が必要となる。分割面ＡＢＥＦを面ＢＣＤＥに加えて、磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）を計算する際の積分路を積分路ＡＢＣＤＥＦＡとしたのは、そのようなことを回避するためである。他方の分割面ＡＦＧＨでは、そのようなことが生じないことから、一つの面として扱う。
【００４５】
一方、磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ＋１／２，ｋ）も同様に、積分路ＡＦＧＨＡに沿って計算される。図３（ｂ）に示すケースでも同様に、磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ＋１／２，ｋ）を計算する際の積分路は積分路ＢＣＤＥＦＧＢとなり、磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ＋３／２，ｋ）を計算する際の積分路は積分路ＡＢＧＦＨＩＡとなる。
【００４６】
導体を含む面に垂直で、その導体が交わる格子（辺）を含む面の磁界強度Ｈは、その後に計算する。その計算は、（１０）〜（１４）式により補正された電界強度Ｅを用いて行う。例えば図３（ａ）に示すケースでは、磁界強度Ｈ_ｘ　（ｉ，ｊ＋１／２，ｋ＋１／２）は電界強度Ｅ_ｙ　（ｉ，ｊ＋ＡＢ／２Δｙ，ｋ）、Ｅ_ｙ　（ｉ，ｊ＋１−ＡＨ／２Δｙ，ｋ）を用いて次式により計算する。他の磁界強度Ｈ_ｘ　（ｉ，ｊ＋１／２，ｋ＋１／２）、Ｈ_ｘ　（ｉ＋１，ｊ＋１／２，ｋ−１／２）、及びＨ_ｘ　（ｉ＋１，ｊ＋１／２，ｋ−１／２）、図３（ｂ）に示すケースでは、磁界強度Ｈ_ｘ　（ｉ，ｊ＋１／２，ｋ＋１／２）、Ｈ_ｘ　（ｉ＋１，ｊ＋１／２，ｋ−１／２）、Ｈ_ｙ　（ｉ＋１／２，ｊ＋１，ｋ−１／２）、及びＨ_ｙ　（ｉ＋１／２，ｊ＋１，ｋ−１／２）についても同様の式で計算することができる。
【００４７】
【数８】

【００４８】
格子を横切る導体と、格子上に近似された導体との接続は簡単に行うことができる。例えば図３（ａ）に示すケースで辺ＡＨは導体の内部に存在するのであれば、その辺ＡＨ上で定義する電界強度Ｅｙ　をゼロとするだけで良い。このようなことから、解析空間に存在する任意の解析対象物体（構造物）を線状の導体として表現できる。
【００４９】
上述したように、例えば図３（ａ）に示すケースでは、磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）の計算には積分路ＡＢＣＤＥＦＡを用い、辺ＡＢ上の電界強度Ｅ_ｙ　は辺ＢＣ上の電界強度Ｅ_ｙ　と定義している。そのようにすると、導体で面を分割しても、計算する強度の数が増えるようなことは回避される。つまり、必要な計算機資源が増大するようなことを回避しつつ、解析対象物体を線として模擬できるようになる。プログラム開発の面では、計算した強度を代入するための配列変数を新たに用意する必要性が回避されることから、言い換えれば、計算した強度の管理は従来と同じ方法で行えることから、その開発が容易となる。既存の解析プログラムに本発明を適用する場合には、そのためのアップデートをより容易に行えることになる。
【００５０】
辺ＢＣ上の電界強度Ｅ_ｙ　を辺ＡＢ上の電界強度Ｅ_ｙ　と定義すると、辺ＢＣ上の電界強度Ｅｙ　は辺ＡＣ上で求めた電界強度Ｅ_ｙ　として扱われることになる。しかし、ＦＤＴＤ法では中間差分が採用されていることもあって、そのように扱っても誤差は小さな範囲内に収まりやすい。磁界強度Ｈを上述したように計算し、電界強度Ｅを上述したように定義するのはそのような理由からである。計算機資源に余裕があるような場合には、導体を境に面を分割し、分割した面毎に電界強度、及び磁界強度の少なくとも一方を求めるようにしても良い。
【００５１】
導体を表現する線は、セルを構成する面上に配置しなければならない。しかし、格子上に配置する場合と比較すると、線を配置できる場所が直線上から平面上に変化したために、線を配置するうえでの制約は大幅に緩和して自由度は大きく向上することになる。直線でない線も扱えるようになる。このため、高い汎用性が得られる。たとえ階段近似を行う場合であっても、格子上に沿って線を配置する必要がないために、セルサイズを小さくすることなく、より小さな誤差で形状を表現することができる。
【００５２】
図４は、本実施の形態によるシミュレーション装置の回路構成図である。その装置は、上述したように導体を模擬して、例えばサージ解析を行うためのものである。
【００５３】
その装置は、例えばコンピュータ（データ処理装置）に本実施の形態による解析プログラムを実行させることで実現されるものである。図４に示すように、ＣＰＵ３０１、ＲＯＭ３０２、ＲＡＭ３０３、通信インターフェース３０４、例えばハードディスク装置である記憶装置３０５、可搬の記録媒体３３０の読み取りを行う記録媒体読み取り装置３０６、キーボードやポインティングデバイスなどを含めて表記した入力装置３０７、及び例えば不図示の表示装置に画像を出力する出力装置３０８がバス３０９により互いに接続されて構成されている。
【００５４】
上記解析プログラムは、記憶装置３０５、記録媒体３３０、或いは外部装置がアクセス可能な記録媒体３２０に格納されている。記録媒体３２０に解析プログラムが格納されている場合、ＣＰＵ３０１は、外部装置、ネットワーク、通信インターフェース３０４、及びバス３０９を介してその解析プログラムを取得し実行する。
【００５５】
図５は、上記解析プログラムの機能構成を示す図である。
図５中の形状ファイル５２１は、解析空間やその空間内に存在する解析対象物体などの形状や配置などの設定データを格納したファイルである。線状の導体として扱う構造物（媒質）は、例えば線の始点と終点の座標を示すデータで形状ファイル５２１中に設定されている。パラメータファイル５２２は、解析を行ううえで用いるべきパラメータの値などを格納したファイルである。電圧源や境界条件などに関する設定はそのファイル５２２を介して行われる。出力ファイル５２３は、解析結果を出力するファイルである。それらのファイル５２１〜５２３は、ユーザが指定する。
【００５６】
解析プログラム５００は、形状ファイル５２１、及びパラメータ５２２を入力ファイルとして、それに格納されているデータを用いて電磁界解析（サージ解析、など）を行い、その解析結果を出力ファイル５２３に出力するものとして実現されている。
【００５７】
解析プログラム５００の機能構成は、入力部５０１、解析部５０２，及び出力部５０２に大別される。入力部５０１は、形状ファイル５２１、及びパラメータファイル５２２に格納されたデータを取り込み、解析部５０２に提供することを主に行う。
【００５８】
解析部５０２は、境界条件処理部５０４、電界解析部５０５、及び磁界解析部５０６を備えている。境界条件処理部５０４は、設定された境界条件を電界解析部５０５、及び磁界解析部５０６に提供する。電界解析部５０５は、それが直前の時間ステップで計算した電界強度Ｅや磁界解析部５０６が計算した磁界強度Ｈを用いて現時点の時間ステップでの電界強度Ｅを計算する。他方の磁界解析部５０６は、それが直前の時間ステップで計算した磁界強度Ｈや電界解析部５０５が計算した電界強度Ｅを用いて現時点の時間ステップでの磁界強度Ｈを計算する。出力部５０３は、各解析部５０５、５０６が計算によりそれぞれ求めた強度Ｅ，Ｈを時刻に対応付けして出力ファイル５２３に出力する。
【００５９】
電界解析部５０５は、電界強度計算部５０７と電界強度補正部５０８を備える。電界強度計算部５０７は、セルの格子を導体が横切っていないことを想定して電界強度Ｅを計算する。電界強度補正部５０８は、電界強度計算部５０７が計算した電界強度Ｅを入力し、導体が格子を横切っている面の電界強度Ｅを上述したように補正する。電界強度Ｅを全て計算した後に補正を行うのは、導体が格子を横切る面では、その格子上の電界強度Ｅとして他の面の電界強度Ｅを定義していることから、他の面の更新を行う前の電界強度Ｅを定義してしまうことを確実に回避するためである。出力部５０３には、その補正を行った後の電界強度Ｅを出力する。
【００６０】
他方の磁界解析部５０６は、電界解析部５０５と同様に、磁界強度計算部５０９と磁界強度補正部５１０を備える。磁界強度計算部５０９は、セルの格子を導体が横切っていないことを想定して磁界強度Ｈを計算する。磁界強度補正部５１０は、磁界強度計算部５０９が計算した磁界強度Ｈを入力し、導体が格子を横切っている面の磁界強度Ｈを上述したように補正する。ここでも同様に、磁界強度Ｈを全て計算した後にその補正を行う。出力部５０３には、その補正を行った後の磁界強度Ｈを出力する。
【００６１】
次に、上記機能構成を有する解析プログラムをＣＰＵ３０１が実行することで実現される電磁界解析（例えばサージ解析）処理について、図６に示すそのフローチャートを参照して詳細に説明する。
【００６２】
先ず、ステップＳ１では、ユーザが指定した形状ファイル５２１、及びパラメータファイル５２２に格納されていうデータにアクセスすることにより、それらのファイル５２１、５２２で設定されている解析条件を読み込み、次の時間ステップでの計算を行ううえでの初期設定を行う。続くステップＳ２では、解析空間の持つ６つの面（それを囲む面）でそれぞれ定められた境界条件の処理を行う。その後はステップＳ３に移行する。
【００６３】
ステップＳ３では、解析空間全体で電界強度Ｅの計算を行う。その計算は、セルの格子を横切る導体が存在していないことを想定して行う。その次に移行するステップＳ４では、導体の格子の横切り方に応じて、（１０）〜（１４）式を用いて説明したように導体周辺の電界強度Ｅを補正する。その後に移行するステップＳ５では、格子上の導体に関する電界強度Ｅの補正を行い、導体内に位置する格子上の電界強度Ｅはゼロにする。ステップＳ６にはその後に移行する。
【００６４】
ステップＳ６では、解析空間内に存在する電圧源、電流源といった給電部の置かれた位置の電界強度Ｅを計算する。次のステップＳ７では、ステップＳ２〜Ｓ６を実行することで得られた電界強度Ｅや電圧値などの計算結果を出力ファイル５２３に出力する。その後はステップＳ８に移行する。
【００６５】
ステップＳ８では、解析空間全体で磁界強度Ｈを計算する。続くステップＳ９では、導体の格子の横切り方に応じて、（１５）、（１６）式を用いて説明したように導体周辺の磁界強度Ｈを補正する。その次に移行するステップＳ１０では、ステップＳ８、或いはＳ９を実行することで得た磁界強度Ｈや電流値などの計算結果を出力ファイル５２３に出力する。その後はステップＳ１１に移行して、解析を行った時間の長さである時刻ｔが最大観測時間Ｔ_ｍａｘ　より大きいか否か判定する。定められた最大観測時間Ｔ_ｍａｘ　分の解析が終了した場合、判定はＹＥＳとなり、ここで電磁界解析処理を終了する。そうでない場合には、判定はＮＯとなって上記ステップＳ１に戻り、時間ステップを１つ進めて、それ以降の処理を同様に実行する。
【００６６】
次に、本実施の形態によるシミュレーション装置を用いて行った解析結果について説明する。
図７は、導体の位置関係を示す図である。解析空間内に存在する導体やそれらの間の位置関係などを示したものである。
【００６７】
その解析空間では、大地に水平に設置された導体７０１が、その大地を覆う銅板７０２から０．６ｍの位置に置かれている。その導体７０１は、長さ４［ｍ］、半径１５×１０^−３［ｍ］のアルミパイプである。その導体７０１には、特性インピーダンス５０オームのケーブルに充電した電荷を水銀リレーで投入するタイプのパルスジェネレータ（ＰＧ）７０３が、断面積２［ｍｍ^２　］の銅線により接続されている。
【００６８】
そのような解析対象物体を持つ解析空間において、図９に示すような電圧波形をＰＧ７０３から導体７０１に印加した場合の印加端での電流・電圧波形を求める解析を検証解析として行った。その実測結果は非特許文献２から得た。また、比較・検討のために、導体７０１を格子上に配置する従来の手法（以降、「従来手法１」と記す）及び導体７０１を格子上に配置せずにそれを階段近似する従来の手法（以降、「従来手法２」と記す）による解析を行った。導体７０１は、従来手法１では図８（ａ）に示すように、ｙ、ｚの各軸方向とは直交し、ｘ軸方向とは平行な面の格子上に配置したのに対し、本実施の形態では図８（ｂ）に示すように、ｚ軸方向とは直交するｘｙ平面にｘ、ｙの各軸方向とは交差させる形で配置した。従来手法２でも基本的には図８（ｂ）に示すようにして階段状に配置した。
【００６９】
非特許文献２に記載されたＦＤＴＤ法による計算では、導体７０１の存在する位置に定義された導体７０１方向成分をゼロとし、その半径は考慮していない。解析空間を分割するセルの一辺の長さをΔｓとすると、ＦＤＴＤ法による計算半径は０．２２９８Δｓ［ｍ］として取り扱っている。
【００７０】
印加端解放時の電圧波形は、図９に示すように模擬し、銅板７０２で模擬した大地の抵抗率は１．６９×１０^−８オーム・メートルとしている。実験で用いられたカレントセンサは立ち上がり時間が２ｎｓの製品である。それにより、単位ステップ電流が流れたとしても、その応答は次式で表される１次遅れとなる。
【００７１】
【数９】

【００７２】
（１７）式を微分すると、単位インパルス応答ｈ（ｔ）は次式で表される。
【００７３】
【数１０】

【００７４】
よって、観測される波形とＦＤＴＤ法を用いた計算結果を比較する場合、実測結果に次のコンボリュージョン演算を行っている。
ｉ_０　＝ｈ（ｔ）＊ｉ（ｔ）　　　　　　　　　　　　　（１９）
その演算を行って得た実測結果と比較する形で本実施の形態による解析結果を図１０に示す。同様に、その実測結果と比較する形で従来手法１、２による解析結果を図１１、図１２にそれぞれ示す。また、本実施の形態の解析結果を従来手法１の解析結果と比較する形で図１３に示す。
【００７５】
図１０〜図１３から明らかなように、導体７０１を格子上に配置する、つまり解析空間内に存在する導体を最適な状態に配置する従来手法１は、導体７０１を階段近似する従来手法２と比較して解析結果に高い精度が得られることが分かる。本実施の形態の解析結果は、その従来手法２のそれと比較して、精度が大きく向上し、従来手法１に近い精度が得られることが分かる。その従来手法１に近い精度が得られることから、本発明を電磁界解析に適用すると、導体７０１を格子上に配置しなくとも、従来の手法では得られない高い精度を達成できることが確認できた。
＜第２の実施の形態＞
上記第１の実施の形態では、線で表現（模擬）する導体で分割された分割面のなかで磁界強度Ｈを計算する位置を含まないものは隣接する面に仮想的に加えて扱っている。そのようにして加える分割面では、隣接する面（ここでは元の面と同じ大きさである）に比べて小さいことから、理論上、それを貫く磁束はゼロと近似している。それにより、例えば図３（ａ）に示すケースでは、面ＡＢＥＦを貫く磁束はゼロとし、図３（ｂ）に示すケースでは、面ＧＥＦを貫く磁束をゼロとしている。しかし、そのように近似する場合、導体の配置によっては比較的に大きな誤差が生じてしまう可能性がある。第２の実施の形態は、その誤差をより抑えられるようにしたものである。以降、第１の実施の形態から異なる部分について詳細に説明する。
【００７６】
導体を表現する線は、図３（ａ）、（ｂ）に示すように、対向する２辺、或いは隣接する２辺を横切る。後者における説明は前者のそれと比較して容易である。このため、前者に注目して、第１の実施の形態から異なる部分を説明する。
【００７７】
第２の実施の形態では、導体を表現する線が対向する２辺を横切る場合、その線の中点を考慮する。それにより、図３（ａ）に示すようなケースでは、図１４（ａ）に示すように、導体を表現する線ＡＦの中点Ｌを考慮する。点Ｋ、Ｊは、中点Ｌとｙ方向上の位置が同じである辺ＢＨ、ＥＧ上の点である。点Ｉは、点Ｆとｙ方向上の位置が同じである辺ＢＨ上の点である。
【００７８】
線ＡＦ（以降、「導体ＡＦ」とも表記する）が導体である場合、三角形の形状をした面ＪＡＬ、ＫＦＬでは、それを貫く磁束は方向は逆で大きさは等しい。この導体ＡＦを解析空間内に表現するために、導体ＡＦを含んだ積分路ＡＢＣＤＥＦＡにファラデーの法則を用いて計算した磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）と、積分路ＡＢＣＤＥＫＪＡに，同様にファラデーの法則を用いて計算した磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）は等しくなるという条件を利用する。
【００７９】
面ＡＢＣＤＥＦ、ＡＢＣＤＥＫＪは、面ＪＡＬ、ＫＦＬの面積が等しいために、それらの面積も等しくなっている。中点Ｌをさらに考慮するのはこのためである。それにより、面ＡＢＣＤＥＦと面積、磁界強度Ｈが共に等しい別の面を利用して、その面積、つまり面ＢＣＤＥの面積に面ＡＢＥＦの面積を加えた面積での磁界強度Ｈをより高精度に算出するようにしている。
【００８０】
面ＢＣＤＥの磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）はファラデーの法則を用いて次式で計算される。
【００８１】
【数１１】

【００８２】
ここで右辺の第２項目では、
−Ｋ_３ｚ×１／（面ＢＣＤＥの面積）×（面ＢＣＤＥ上の電界強度Ｅの周回積分値）
という関係が成立している。次に、この関係を用いて、積分路ＡＢＣＤＥＦＡ、ＡＢＣＤＥＫＪＡの各各積分路についてファラデーの法則を用いて磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）を計算し、それぞれの積分路で計算された磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）は等しいという条件に従い、電界強度Ｅ_ｙ、Ｅ_ｘ　をｘ、ｙ方向で中央差分することで（２０）式が導出できる。ただし、導体ＡＦとの位置関係から、電界強度Ｅ_ｘ　の線ＫＬ上の積分値と線ＬＪ上の積分値は大きさが等しく符号が逆となり、電界強度Ｅ_ｘ　の線ＫＬ上の積分値は零（ゼロ）となる。また、電界強度Ｅ_ｙ　の線ＪＡ上の積分値と線ＫＦ上の積分値は等しくなる。
【００８３】
【数１２】

【００８４】
さらに，以上の（１９）（２０）式により計算された磁界強度Ｈは等しいので、辺ＢＥ上の電界強度Ｅ_ｘ　（ｉ＋１／２，ｊ，ｋ）は次式から計算される。
【００８５】
【数１３】

【００８６】
上述したようにして導出する磁界強度Ｈでは、考慮すべき面積を全て考慮している。それにより、磁界強度Ｈの導出に、仮想的に加えた分割面の面積分の磁束が反映されることになる。このため、（１５）式から導出されるものと比較して、より正確になる。（２１）式は、そのように磁界強度Ｈをより正確に導出することを前提にして得られる式であることから、それにより導出される電界強度Ｅもより正確なものとなる。
【００８７】
ところで、（２０）式は時間領域についてのみ中央差分を適用したファラデーの法則から磁界強度Ｈ_ｚ　（ｉ＋１／２，ｊ−１／２，ｋ）を導出するものである。電界強度Ｅ_ｘ　にはｙ方向上に１次オーダーの打ち切り誤差を含む中央差分が用いられている。表記法を変えて（２０）式の内容を示したものを（２２）式に、その中央差分を（２３）式にそれぞれ示す。
【００８８】
【数１４】

【００８９】
（２３）式が１次オーダーの打ち切り誤差を含むことは以下のように証明できる。
まず、次のテイラー展開を考える。
【００９０】
【数１５】

【００９１】
ただし、ｕ＝１＋２ＥＫ／Δｙである。（２５）式−（２４）式より
【００９２】
【数１６】

【００９３】
となり、１次オーダーの打ち切り誤差を生じることがわかる。（２６）式中の「δ」はクロネッカーのデルタであり、「ο（２）」は２次以上の項を表すシンボルである。
【００９４】
第２の実施の形態では、１次オーダーの打ち切り誤差を以下のようにしてさらに２次オーダーに減少させている。まず、次の新たなテイラー展開を考える。
【００９５】
【数１７】

【００９６】
（２７）、（２８）式より、
【００９７】
【数１８】

【００９８】
また、（２６）式＋（２９）式より、
【００９９】
【数１９】

【０１００】
となる。以上より、
【０１０１】
【数２０】

【０１０２】
が導出される。そのようにして陰的に導出される（３１）式を用いることにより、（１５）式では回避できない面積ＡＢＥＦを貫く磁束を近似的に零としたことによる誤差、さらには（２０）式に含まれる１次オーダーの打ち切り誤差が軽減される。それにより、第１の実施の形態と比較して、より正確な解析モデルを得ることができる。図１４（ｂ）に示すケースでは、磁界強度Ｈ_ｚ　^{ｎ＋１／２}（ｉ＋１／２，ｊ＋１／２，ｋ）も同様の手法により計算される。
【０１０３】
上述したような計算を図５に示すシミュレーション装置に行わせる場合、（３１）式を用いた磁界強度Ｈの計算は、図７に示す電磁界解析処理では例えばステップＳ９で行わせれば良い。（２１）式を用いた電界強度Ｅの補正は、例えばステップＳ４で行わせれば良い。このようなことから、図５に示すシミュレーション装置での第２の実施の形態の実現は、その動作の流れは基本的に第１の実施の形態におけるそれから変更させることなく行うことができる。
【０１０４】
次に、上述したようにして実現できる第２の実施の形態によるシミュレーション装置を用いて行った解析結果について説明する。その解析は、上記第１の実施の形態と同じもの、即ち図７に示すように導体が存在する解析空間解析空間において、図９に示すような電圧波形をＰＧ７０３から導体７０１に印加した場合の印加端での電流・電圧波形を求める解析を検証解析として行った。計算モデルが異なる以外、各種条件等は第１の実施の形態と全て同じである。そのようにして得られた解析結果を従来手法１による解析結果と比較する形で図１５に示す。また、従来手法１、第１の実施の形態の各モデル、更には実測結果を含め、導体７０１の始端から伝搬した電流が終端で反射され始端に戻るまでの伝搬時間（Ｐｒｏｐａｇａｔｉｏｎ　ｔｉｍｅ）［ｎｓ］、電流ピーク値（Ｐｅａｋ　ｖａｌｕｅ　ｏｆ　ｃｕｒｒｅｎｔ　）［Ａ］を図１６に示す。
【０１０５】
図１５、図１６、及び図１３から明らかなように、解析結果は従来手法１のそれと非常に良く一致している。第１の実施の形態と比較して、大きく精度が向上していることが分かる。これらのことから、第２の実施の形態で採用したモデルは、精度を向上させるうえで非常に有効であることが確認できる。特には図示していないが、（３１）式の代わりに（２０）式を用いても、第１の実施の形態と比較して精度が向上することが確認された。
【０１０６】
なお、第２の実施の形態では、導体によって分割されて隣接する面に仮想的に加える分割面の面積を正確に扱うために、言い換えればより高い精度を得るために、導体の中点（図１４（ａ）では「Ｌ」の他に「Ｋ」や「Ｊ」も対応）を考慮しているが、別の点を考慮しても良い。例えば図１４（ａ）に示すケースでは、点Ｌの代わりに点Ａ、Ｆを考慮するようにしても良い。
【０１０７】
本実施の形態（第１、及び第２の実施の形態）では、導体を表現（模擬）する線として直線を用いているが、曲線であっても良い。線で表現する対象は、導体や構造物などに限定されるわけではなく、解析空間のなかに存在する媒質のなかから必要に応じて選択すれば良い。導体が横切る辺での電界強度Ｅの近似方法などについても、本実施の形態に限定されるわけではなく、様々な変形を行うことができる。
【０１０８】
【発明の効果】
以上、説明したように本発明は、解析空間に存在する媒質のなかから選択した媒質（例えば導体や構造物など）を線として設定（表現）し、セルを構成する面のなかでその線が辺を横切る面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求める。
【０１０９】
その線で面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求めることにより、その線はセルの面の辺上に配置した形で扱うことになる。その線はセルの面上に配置すれば良いから、高い自由度が得られ、汎用性は向上する。高い自由度を実現させたことで線で表現する媒質の形状はより正確に模擬できるようになる。これらのことから、高い汎用性を実現させつつ、電磁界解析は確実に高精度に行えるようになる。
【０１１０】
上記線による面の分割を、その面をその線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、その位置を含む分割面は隣接する面に加えない形で行うようにした場合には、その分割に伴い面の数が増大することは回避される。その結果、最終的に求める強度の数、即ち解析結果として保存の対象となるデータの量の増大は回避でき、必要な計算機資源の増大は抑えることができる。
【０１１１】
分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面において、その変化を電界強度、或いは磁界強度の計算に反映させるようにした場合には、その変化に伴って生じる誤差は回避、或いは低減させることができる。このため、その変化を計算に反映させる強度はより高精度に得ることができる。強度の変化は、隣接する面間で伝搬していく。それにより、近傍の面で求められる強度には、媒質が同じであれば比較的に大きな違いがないのが普通である。このことから、その反映を面積の変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行うようにしたときには、その変化に伴って生じる誤差の回避、或いは低減はより確実に行うことができる。その結果、計算する強度は更に高精度に得ることができるようになる。
【図面の簡単な説明】
【図１】解析空間を分割するセルを説明する図である。
【図２】電磁界の時間配置を説明する図である。
【図３】セルを構成する面とその面の辺を横切る導体の位置関係を説明する図である。
【図４】本実施の形態によるシミュレーション装置の回路構成図である。
【図５】解析プログラムの機能構成を示す図である。
【図６】電磁界解析処理のフローチャートである。
【図７】導体の位置関係を示す図である。
【図８】導体の配置を説明する図である。
【図９】パルスジェネレータの電圧波形の時間変化を示す図である。
【図１０】第１の実施の形態による解析結果を示す図である。
【図１１】従来手法１による解析結果を示す図である。
【図１２】従来手法２による解析結果を示す図である。
【図１３】第１の実施の形態による解析結果を従来手法１による解析結果と比較して示す図である。
【図１４】第２の実施の形態における磁界強度、及び電界強度の計算方法を説明するための図である。
【図１５】第２の実施の形態による解析結果を示す図である。
【図１６】各モデルでの伝搬時間と電流ピーク値を示す図である。
【符号の説明】
３０１　　ＣＰＵ
３０２　　ＲＯＭ
３０３　　ＲＡＭ
３０４　　通信インターフェース
３０５　　記憶装置
３０６　　記憶媒体読み取り装置
３０７　　入力装置
３０８　　出力装置
３２０、３３０　　記録媒体
５００　　解析プログラム
５０１　　入力部
５０２　　解析部
５０３　　出力部
５０４　　境界条件処理部
５０５　　電界解析部
５０６　　磁界解析部
５０７　　電界強度計算部
５０８　　電界強度補正部
５０９　　磁界強度計算部
５１０　　磁界強度補正部
７０１　　導体
７０２　　銅板
７０３　　パルスジェネレータ[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a technique for dividing an analysis space into hexahedral cells and performing an electromagnetic field analysis using an FDTD (Finite Difference Time Domain) method.
[0002]
[Prior art]
Electromagnetic field analysis is performed in various fields including an antenna problem and an electromagnetic wave scattering problem. As the analysis method, the finite element method, the moment method, the boundary element method, and the like have been used. S. The FDTD (Finite Difference Time Domain) method, which is a finite difference method proposed by Yee, has attracted attention. This is because the algorithm has advantages such as simplicity, excellent accuracy, and wide application range.
[0003]
According to the FDTD method, the electromagnetic field is differentiated in the time domain and the analysis space, and the behavior of the electromagnetic field is numerically determined from the electrical (physical) characteristics such as the geometrical arrangement and shape of the object to be analyzed, conductivity, and magnetic permeability. This is a method that is obtained by analysis. In the time domain, discretization is performed in a short time (Δt), and the analysis space is divided into units called cells to perform discretization.
[0004]
FIG. 1 is a diagram illustrating cells that divide the analysis space. As shown in FIG. 1, the cell is set as a hexahedron (usually a cube) having a length of Δx in the x-axis direction, Δy in the y-axis direction, and Δz in the z-axis direction. The magnetic field intensity is arranged at the center of each surface, and the electric field intensity is arranged at the center of each side constituting the surface. In this way, Maxwell's equation that the rotation of an electric field creates a magnetic field and the rotation of a magnetic field creates an electric field can be directly solved by spatially shifting the magnetic field strength and the electric field strength by half a cell and making a difference. It is arranged. The cell size is generally set to 1/10 or less of the shortest wavelength of interest.
[0005]
The electric field and the magnetic field that are spatially shifted and differentiated are arranged alternately in space and time, as shown in FIG. A leapfrog algorithm that calculates the electric field strength and the magnetic field strength alternately is used. Thereby, they are calculated alternately, such as obtaining the magnetic field strength from the electric field strength and obtaining the electric field strength from the magnetic field strength. The superscripts such as "n-1 / 2" and "n" shown in FIG. 2 indicate positions arranged in the time domain. The difference in the time domain is performed so as to satisfy the Courant stability condition.
[0006]
In the FDTD method, as described above, the central difference method is used for differentiating between the time domain and the analysis space. Depending on the field in which the analysis is performed, the arrangement of the electric field and the magnetic field in the time axis and space coordinates may be reversed.
[0007]
Maxwell's equations, which are the basis of the FDTD method, are based on the law of amperial integration and the law of Faraday's electromagnetic induction. Each law is represented by a differential form and an integral form, and the Maxwell's equation generally refers to a total of four equations combined with Gauss's law regarding electric and magnetic fields. Maxwell's equation of the differential form is as follows: electric field intensity E [V / m], magnetic field intensity H [A / m], electric flux density D [C / m].²], Magnetic flux density B [T], charge density ρ [C / m³], Current density J [A / m²] Is expressed as follows.
[0008]
(Equation 1)

[0009]
Here, Equation (1) (Amper's law), Equation (2) (Faraday's law of electromagnetic induction), Equation (3) (Gauss's law in a magnetic field), Equation (4) (Gauss's law in an electric field) Is not an independent equation, but in the FDTD method, equations (3) and (4) are used as evaluation criteria for numerical errors, and equations (1) and (2) are used for formulation. When the central difference method is applied to the spatio-temporal and spatial coordinates using the Yee's algorithm in the equations (1) and (2) and formulated, the following FDTD basic equation can be obtained in a three-dimensional space.
[0010]
(Equation 2)

[0011]
Here, μ: magnetic permeability, ε: dielectric constant, σ: resistivity.
In the expressions (5) and (6), “x” and “y” added as subscripts to the symbols “E” and “H” represent the directions of the electric field strength and the magnetic field strength, respectively. Although only the electric field strength and the magnetic field strength assigned in the x direction are shown here, the same can be derived for the y and z directions.
[0012]
[Non-patent document 1]
Hano, "FDTD (Finite Difference Time Domain) Method", Proceedings of the 2002 IEEJ National Convention, 5th volume, pp411-414, 2002.
[Non-patent document 2]
Noda, Yokoyama, "Development of a general-purpose surge analysis program based on the FDTD method", Denki Kagaku B, Vol. 121, No. 5, pp625-632, 2001.
[Non-Patent Document 3]
J. Fang @ and @ J. Ren, "A Locally Conformed, Finite-Difference, Time-Domain, Algorithm, of Modeling, Arbitrary, Shape, Plana, Metal, Metal, Strip, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Satellite, Transport, Transport, Transport, Transport, Transportation, Transport, Transport, Transport, Transportation, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Transport, Vehicle Information, Technology & Information Systems. 41, NO. 5, MAY $ 1993
[0013]
[Problems to be solved by the invention]
The FDTD method is based on dividing an analysis space into cells as shown in FIG. Therefore, when expressing an object to be analyzed such as a conductor that does not exist on the grid (side) of the cell, a staircase approximation method that approximates the shape in a stepwise manner has been used. However, there is usually a great difference between the shape represented by the step approximation method and the actual shape. For this reason, the conventional analysis method using the staircase approximation method has a problem that a large error occurs in the resonance frequency of the object to be analyzed and the propagation time thereof, and the analysis cannot be performed with high accuracy.
[0014]
By reducing the cell size, the difference between the shape represented by the staircase approximation method and the actual shape can be reduced. However, when the size is reduced, the number of cells increases accordingly, and the necessary computer resources and calculation time become enormous. For this reason, it is not desirable to reduce the cell size even if it is partial. In the field of surge analysis, a relatively large analysis space is often used. In other words, the number of cells is usually large, which is particularly undesirable.
[0015]
In the case where the shape of a cell is transformed to represent the shape of an object to be analyzed, modeling according to the object is required, which lacks generality. The method of setting the coordinate system so that the object is on the grid cannot be generally used if a plurality of objects exist in the analysis space in an arbitrary state. For this reason, it is important to be able to more accurately represent and analyze the analysis target object existing in the analysis space without performing operations such as reducing the cell size or deforming the cell shape. It is considered to be.
[0016]
An object of the present invention is to provide a general-purpose technique for always accurately expressing an object to be analyzed (medium) existing in an analysis space and for reliably performing electromagnetic field analysis with high accuracy. With the goal.
[0017]
[Means for Solving the Problems]
The electromagnetic field analysis method using the FDTD method of the present invention is a method for dividing an analysis space into hexahedral cells and performing an electromagnetic field analysis using an FDTD (Finite Difference Time Domain) method. A medium selected from the media existing in the analysis space to be divided is set as a line, and among the surfaces constituting the cell, the surface where the line intersects the side is divided into a plurality of lines, and the electric field intensity or the magnetic field intensity is divided. Ask.
[0018]
The division of the plane by the line is performed by dividing a plane including a position where a magnetic field intensity is calculated among a plurality of division planes obtained by dividing the plane by the line, in addition to an adjacent plane and dividing the plane including the position. It is desirable that the surface is formed without adding to the adjacent surface. It is preferable that the electric field intensity of the surface divided into a plurality of lines adopts the same electric field intensity in the direction calculated around the surface. The magnetic field strength of the surface divided into a plurality of lines is calculated using the electric field intensity obtained by the integration path including the line, and it is desirable that the magnetic field strength is adjacent to the direction intersecting the surface divided into a plurality of lines, In addition, it is desirable that the magnetic field strength of a surface having a side crossed by the line be calculated on the assumption that the surface is adjacent to a plurality of divided surfaces. It is desirable that the electric field strength of a surface divided into a plurality of lines be obtained by calculating all the electric field strengths without taking the lines into account and then correcting the calculated electric field strengths.
[0019]
On a surface where a virtual change occurs in area by adding a dividing surface to an adjacent surface, it is desirable to reflect the change in calculation of electric field strength or magnetic field strength. The reflection is desirably performed using the electric field strength or the magnetic field strength obtained on the surface located near the surface where the change occurs.
[0020]
The method for expressing a medium in electromagnetic field analysis according to the present invention is a method used to divide an analysis space into hexahedral cells and express the medium in an electromagnetic field analysis performed by the FDTD method. Is selected, and the selected medium is set and expressed as a line crossing the side of the surface constituting the cell.
[0021]
A simulation apparatus according to the present invention is based on the premise that an analysis space is divided into hexahedral cells and electromagnetic field analysis is performed using the FDTD method, and a medium setting unit that sets a medium existing in the analysis space as a line; And analyzing means for dividing a plane which intersects a side of a cell in which a line set by the setting means divides the analysis space with the side to obtain an electric field strength or a magnetic field strength.
[0022]
The division of the plane by the line is performed by dividing a plane including a position where a magnetic field intensity is calculated among a plurality of division planes obtained by dividing the plane by the line, in addition to an adjacent plane and dividing the plane including the position. It is desirable that the surface is formed without adding to the adjacent surface. It is preferable that the analysis means adopts, as the electric field intensity of the surface divided into a plurality of lines, the electric field intensity in the same direction calculated around the surface.
[0023]
In addition, it is preferable that the analysis unit reflects the change in the calculation of the electric field strength or the magnetic field strength on a surface where a virtual change occurs in the area by adding the dividing surface to the adjacent surface. The reflection is desirably performed using the electric field strength or the magnetic field strength obtained on the surface located near the surface where the change occurs.
[0024]
The program of the present invention is based on the premise that the analysis space is divided into hexahedral cells and executed by a data processing device used as a simulation device for performing electromagnetic field analysis by the FDTD method, and the medium existing in the analysis space is represented as a line. A function of setting and a function of obtaining an electric field strength or a magnetic field strength by dividing a surface which intersects a side of a cell in which a line set by the setting function divides the analysis space with the line is realized.
[0025]
In the present invention, a medium (conductor, structure, or the like) selected from the media existing in the analysis space is set (represented) as a line, and a plurality of planes on which the line intersects the side among the surfaces constituting the cell. By dividing, at least one of the electric field strength and the magnetic field strength is obtained.
[0026]
By dividing the surface by a plurality of lines and obtaining at least one of the electric field intensity and the magnetic field intensity, the line is handled in a form arranged on the side of the cell surface. Since the line may be arranged on the surface of the cell, a high degree of freedom is obtained, and the versatility is improved. The medium represented by the line can more accurately simulate the shape. From these facts, the electromagnetic field analysis can be reliably performed with high accuracy while realizing high versatility.
[0027]
The division of the surface by the line is divided into a plurality of division surfaces obtained by dividing the surface by the line.A division surface not including the position for calculating the magnetic field strength is added to the adjacent surface, and the division surface including the position is If the processing is performed in such a manner that it is not added to adjacent planes, it is possible to avoid an increase in the number of planes due to the division. As a result, an increase in the number of finally obtained intensities, that is, an amount of data to be saved as an analysis result is avoided, and an increase in necessary computer resources is suppressed.
[0028]
When a virtual change in area is caused by adding a dividing surface to an adjacent surface, if the change is reflected in the calculation of electric field strength or magnetic field strength, the error caused by the change is This can be avoided or reduced, and the intensity of the change reflected in the calculation can be obtained with higher accuracy. The change in intensity propagates between adjacent surfaces. As a result, the strength required on the nearby surface does not usually have a relatively large difference if the medium is the same. For this reason, when the reflection is performed using the electric field strength or the magnetic field strength obtained on the surface located near the surface where the change in the area occurs, the calculated intensity can be obtained with higher accuracy. .
[0029]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings.
<First embodiment>
FIG. 3 is a diagram illustrating a positional relationship between a surface constituting a cell and a conductor crossing a side of the surface. First, a method of expressing a conductor (a medium to be distinguished from others) used in the present embodiment and a method of approximating a magnetic field and an electric field around the conductor will be described in detail with reference to FIG.
[0030]
In the present embodiment, the analysis space is divided into hexahedral cells (see FIG. 1) by a general space division method, and the shape and arrangement of conductors are set separately. Thereby, the shape of the conductor inclined at an arbitrary angle at an arbitrary position on the surface constituting the cell is simulated by a line. Conductors that cannot be placed on the cell grid are simulated by lines on the surface.
[0031]
When simulated in such a line, the conductor crossing the side will cross the side opposite to it, or cross the side adjacent thereto. FIG. 3A shows the positional relationship when the conductor traverses the former, that is, two opposing sides, and FIG. 3B shows the positional relationship when the conductor traverses the two adjacent sides. In FIG. 3B, conductors passing through other than the plane ACDE are indicated by broken lines. The position on the surface where the magnetic field strength is calculated is indicated by a circle having a black circle inside.
[0032]
Conventionally, in the FDTD method, only one electric field is defined on each side. However, in the present embodiment, two sides are defined at a point where the conductors intersect at a point where they intersect. The integration path of the magnetic field adjacent to the conductor changes from an integration path along a side of the surface to an integration path including the conductor. In this way, the surface for calculating the electric field strength and the magnetic field intensity is divided by the conductor, and the conductors other than on the grid are considered in such a manner that the strength is calculated for each divided surface (divided surface). Hereinafter, a method of approximating a magnetic field and an electric field around a conductor will be specifically described for each of the cases shown in FIGS. 3A and 3B.
[0033]
In the case shown in FIG. 3A, the conductors intersect at a point A on the side BH and a point F on the side EG. Accordingly, in this case, the side BH is divided into two sides AB and AH, and the side EG is divided into two sides EF and FG._yIs defined. Here, their electric field strength E_yIs defined as follows.
[0034]
(Equation 3)

[0035]
In the parentheses in the expressions (10) and (11), the coordinate positions defining the electric field strength Ey are shown. The symbol “Δy” indicates the cell size in the y direction. “AB” represents the length between point A and point B. They are the same in other expressions.
[0036]
In the case shown in FIG. 3B, the conductors intersect at a point B on the side AC and a point G on the side AE. Therefore, in this case, the side AC is divided into two sides AB and BC, and the electric field strength E in the y direction at the center thereof is obtained._yIs similarly defined as follows.
[0037]
(Equation 4)

[0038]
On the other hand, the electric field intensity E in the x direction assigned near the conductor_xIs calculated as follows.
In the case shown in FIG. 3A, the magnetic field intensity H calculated using the integration path ABCDEFA and the integration path BCDEB is used._zFrom the relationship that ／ (i + 1/2, j-1 / 2, k) must be the same, the electric field strength E_x(I + ／, j, k) is calculated using the following equation.
[0039]
(Equation 5)

[0040]
Similarly, in the case shown in FIG. 3B, the magnetic field intensity H calculated using the integration path BCDEFGB and the integration path BCDEGB is used._z(I + ／, j + ／, k) must be the same, so that the electric field strength E_x(I + 1−EG / 2Δy, j + 1, k) is calculated using the following equation.
[0041]
(Equation 6)

[0042]
Electric field strength E_x(I + AG / 2Δx, j + 1, k) can be calculated by the same formula.
As considered when deriving the equations (13) and (14), the magnetic field intensity around the conductor is calculated by changing its integral path to a shape along the conductor. In the case shown in FIG._zAn integration path for calculating (i + 1/2, j-1 / 2, k) is an integration path ABCDEFA. Since the electric field strength E of the conductor is zero, the magnetic field strength H_z(I + ／, j− ／, k) is calculated using the following equation.
[0043]
(Equation 7)

[0044]
The plane (divided plane) ABEF constituting the plane BEGH crossed by the conductor has a magnetic field strength H_zDoes not include the position where is calculated. The magnetic field strength H at the divided surface ABEF_zWhen is calculated, the number of magnetic field strengths H to be calculated is increased by dividing the surface. A complicated process is required for the management. The division plane ABEF is added to the plane BCDE, and the magnetic field strength H_zThe reason why the integration path ABCDEFA is used as the integration path for calculating (i + 1/2, j-1 / 2, k) is to avoid such a case. In the other divided surface AFGH, such a case does not occur, so that it is treated as one surface.
[0045]
On the other hand, the magnetic field strength H_z(I + ／, j + ／, k) is similarly calculated along the integration path AFGHA. Similarly, in the case shown in FIG._zThe integration path for calculating (i + ／, j + ／, k) is the integration path BCDEFGB, and the magnetic field strength H_zThe integration path for calculating (i + 1/2, j + 3/2, k) is the integration path ABGFHIA.
[0046]
The magnetic field strength H of the plane perpendicular to the plane containing the conductor and including the grid (side) where the conductor intersects is then calculated. The calculation is performed using the electric field intensity E corrected by the equations (10) to (14). For example, in the case shown in FIG._x(I, j + ／, k + ／) is the electric field strength E_y(I, j + AB / 2Δy, k), E_yIt is calculated by the following equation using (i, j + 1−AH / 2Δy, k). Other magnetic field strength H_x(I, j + ／, k + ／), H_x(I + 1, j + ／, k− ／), and H_x(I + 1, j + ／, k− ／), in the case shown in FIG._x(I, j + ／, k + ／), H_x(I + 1, j + ／, k− ／), H_y(I + ／, j + 1, k− ／), and H_y(I + 1/2, j + 1, k− ／) can be calculated by the same formula.
[0047]
(Equation 8)

[0048]
The connection between the conductor crossing the grid and the conductor approximated on the grid can be made easily. For example, in the case shown in FIG. 3A, if the side AH exists inside the conductor, it is only necessary to set the electric field strength Ey defined on the side AH to zero. From this, any object to be analyzed (structure) existing in the analysis space can be represented as a linear conductor.
[0049]
As described above, for example, in the case shown in FIG._z積分 (i + 1/2, j-1 / 2, k) is calculated using the integration path ABCDEFA, and the electric field intensity E on the side AB is calculated._yIs the electric field strength E on the side BC_yIs defined. In this way, even if the surface is divided by the conductor, it is possible to avoid an increase in the number of calculated intensities. That is, the object to be analyzed can be simulated as a line while avoiding an increase in necessary computer resources. In terms of program development, the need to prepare new array variables to substitute the calculated intensities is avoided, in other words, the calculated intensities can be managed in the same way as before, so the development Becomes easier. When the present invention is applied to an existing analysis program, the update for that can be performed more easily.
[0050]
Electric field strength E on side BC_yIs the electric field strength E on the side AB_y, The electric field strength Ey on the side BC is the electric field strength E obtained on the side AC._yWill be treated as. However, since the FDTD method employs an intermediate difference, the error tends to fall within a small range even when handled as such. It is for such a reason that the magnetic field strength H is calculated as described above and the electric field strength E is defined as described above. If there is room for computer resources, the surface may be divided by a conductor, and at least one of the electric field intensity and the magnetic field intensity may be obtained for each divided surface.
[0051]
The line representing the conductor must be placed on the surface that makes up the cell. However, compared to the case of placing on a grid, the place where lines can be placed has changed from a straight line to a plane, so the restrictions on placing lines are greatly relaxed and the degree of freedom is greatly improved. Become. It can handle non-straight lines. Therefore, high versatility can be obtained. Even when staircase approximation is performed, it is not necessary to arrange the lines along the grid, so that the shape can be represented with a smaller error without reducing the cell size.
[0052]
FIG. 4 is a circuit configuration diagram of the simulation device according to the present embodiment. The device simulates a conductor as described above and performs, for example, surge analysis.
[0053]
The device is realized, for example, by causing a computer (data processing device) to execute the analysis program according to the present embodiment. As shown in FIG. 4, a CPU 301, a ROM 302, a RAM 303, a communication interface 304, a storage device 305, for example, a hard disk device, a recording medium reading device 306 for reading a portable recording medium 330, a keyboard and a pointing device are included. The input device 307 described above and an output device 308 that outputs an image to, for example, a display device (not shown) are connected to each other via a bus 309.
[0054]
The analysis program is stored in the storage device 305, the recording medium 330, or the recording medium 320 accessible by an external device. When the analysis program is stored in the recording medium 320, the CPU 301 acquires and executes the analysis program via an external device, a network, the communication interface 304, and the bus 309.
[0055]
FIG. 5 is a diagram showing a functional configuration of the analysis program.
The shape file 521 in FIG. 5 is a file storing setting data such as the shape and arrangement of an analysis space and an analysis target object existing in the space. The structure (medium) to be treated as a linear conductor is set in the shape file 521 with data indicating, for example, the coordinates of the start point and the end point of the line. The parameter file 522 is a file in which parameter values and the like to be used in performing the analysis are stored. Settings relating to the voltage source, boundary conditions, and the like are made via the file 522. The output file 523 is a file for outputting an analysis result. These files 521 to 523 are specified by the user.
[0056]
The analysis program 500 performs an electromagnetic field analysis (surge analysis, etc.) using the shape file 521 and the parameters 522 as input files and data stored therein, and outputs the analysis results to an output file 523. Has been realized.
[0057]
The functional configuration of the analysis program 500 is roughly divided into an input unit 501, an analysis unit 502, and an output unit 502. The input unit 501 mainly takes in data stored in the shape file 521 and the parameter file 522 and provides the data to the analysis unit 502.
[0058]
The analysis unit 502 includes a boundary condition processing unit 504, an electric field analysis unit 505, and a magnetic field analysis unit 506. The boundary condition processing unit 504 provides the set boundary conditions to the electric field analysis unit 505 and the magnetic field analysis unit 506. The electric field analysis unit 505 calculates the electric field intensity E at the current time step using the electric field intensity E calculated at the immediately preceding time step and the magnetic field intensity H calculated by the magnetic field analysis unit 506. The other magnetic field analysis unit 506 calculates the magnetic field strength H at the current time step using the magnetic field strength H calculated at the immediately preceding time step and the electric field strength E calculated at the electric field analysis unit 505. The output unit 503 outputs the intensities E and H calculated by the respective analysis units 505 and 506 to the output file 523 in association with the time.
[0059]
The electric field analysis unit 505 includes an electric field intensity calculation unit 507 and an electric field intensity correction unit 508. The electric field strength calculator 507 calculates the electric field strength E on the assumption that the conductor does not cross the cell grid. The electric field intensity correction unit 508 receives the electric field intensity E calculated by the electric field intensity calculation unit 507, and corrects the electric field intensity E on the surface where the conductor crosses the grid as described above. The reason why the correction is performed after calculating all the electric field intensities E is that the electric field intensity E of the other surface is defined as the electric field intensity E on the lattice on the surface where the conductor crosses the lattice. This is to reliably avoid defining the electric field strength E before performing the above. The output unit 503 outputs the electric field intensity E after the correction.
[0060]
The other magnetic field analysis unit 506 includes a magnetic field intensity calculation unit 509 and a magnetic field intensity correction unit 510, like the electric field analysis unit 505. The magnetic field strength calculation unit 509 calculates the magnetic field strength H on the assumption that the conductor does not cross the cell grid. The magnetic field strength correction unit 510 receives the magnetic field strength H calculated by the magnetic field strength calculation unit 509 and corrects the magnetic field strength H on the surface where the conductor crosses the grid as described above. Here, similarly, the correction is performed after calculating all the magnetic field strengths H. The output unit 503 outputs the corrected magnetic field strength H.
[0061]
Next, an electromagnetic field analysis (for example, surge analysis) process realized by the CPU 301 executing the analysis program having the above functional configuration will be described in detail with reference to the flowchart shown in FIG.
[0062]
First, in step S1, by accessing data stored in the shape file 521 and the parameter file 522 specified by the user, the analysis conditions set in the files 521 and 522 are read, and in the next time step. Perform the initial setting for calculating. In the following step S2, processing is performed on the boundary conditions defined on the six surfaces (surfaces surrounding the analysis space) of the analysis space. Thereafter, the process proceeds to step S3.
[0063]
In step S3, the electric field strength E is calculated in the entire analysis space. The calculation is performed assuming that there are no conductors crossing the cell grid. In the next step S4, the electric field intensity E around the conductor is corrected according to the way of traversing the conductor grid as described using the equations (10) to (14). In the subsequent step S5, the electric field intensity E on the conductor on the lattice is corrected, and the electric field intensity E on the lattice located in the conductor is set to zero. Then, the process proceeds to step S6.
[0064]
In step S6, the electric field strength E at the position where the power supply unit such as the voltage source and the current source existing in the analysis space is located is calculated. In the next step S7, calculation results such as electric field strength E and voltage value obtained by executing steps S2 to S6 are output to the output file 523. After that, it moves to step S8.
[0065]
In step S8, the magnetic field strength H is calculated for the entire analysis space. In the following step S9, the magnetic field strength H around the conductor is corrected according to the way of traversing the grid of the conductor as described using the equations (15) and (16). In the next step S10, a calculation result such as the magnetic field strength H or the current value obtained by executing the step S8 or S9 is output to the output file 523. Thereafter, the process proceeds to step S11, in which the time t, which is the length of the analysis, is set to the maximum observation time T._maxIt is determined whether it is greater than. Determined maximum observation time T_maxWhen the analysis of the minute has been completed, the determination is YES, and the electromagnetic field analysis processing is ended here. If not, the determination is no and the process returns to step S1, proceeds with one time step, and performs the subsequent processing in the same manner.
[0066]
Next, the results of analysis performed using the simulation device according to the present embodiment will be described.
FIG. 7 is a diagram illustrating a positional relationship between conductors. It shows the conductors existing in the analysis space and the positional relationship between them.
[0067]
In the analysis space, a conductor 701 horizontally placed on the ground is placed at a position 0.6 m from the copper plate 702 covering the ground. The conductor 701 has a length of 4 [m] and a radius of 15 × 10^-3[M] aluminum pipe. The conductor 701 has a pulse generator (PG) 703 of a type in which electric charges charged in a cable having a characteristic impedance of 50 ohms are injected by a mercury relay, and has a cross-sectional area of 2 [mm].²].
[0068]
In the analysis space having such an object to be analyzed, an analysis for obtaining a current / voltage waveform at an application end when a voltage waveform as shown in FIG. 9 is applied from the PG 703 to the conductor 701 was performed as a verification analysis. The measurement results were obtained from Non-Patent Document 2. For comparison and study, a conventional method of arranging the conductor 701 on a lattice (hereinafter, referred to as a “conventional method 1”) and a conventional method of arranging the conductor 701 in a stepwise manner without disposing the conductor 701 on the lattice (Hereinafter, referred to as “conventional method 2”). As shown in FIG. 8A, the conductor 701 is arranged on a lattice having a plane orthogonal to the y- and z-axis directions and parallel to the x-axis direction, as shown in FIG. In the form (1), as shown in FIG. 8 (b), they are arranged on an xy plane orthogonal to the z-axis direction so as to intersect with the x- and y-axis directions. Also in the conventional method 2, they are basically arranged stepwise as shown in FIG.
[0069]
In the calculation by the FDTD method described in Non-Patent Document 2, the component in the direction of the conductor 701 defined at the position where the conductor 701 exists is set to zero, and its radius is not considered. Assuming that the length of one side of the cell dividing the analysis space is Δs, the radius calculated by the FDTD method is treated as 0.2298 Δs [m].
[0070]
The voltage waveform when the application terminal is released is simulated as shown in FIG. 9, and the resistivity of the ground simulated by the copper plate 702 is 1.69 × 10^-8Ohm meters. The current sensor used in the experiment is a product having a rise time of 2 ns. As a result, even if a unit step current flows, the response has a first-order delay represented by the following equation.
[0071]
(Equation 9)

[0072]
When the equation (17) is differentiated, the unit impulse response h (t) is expressed by the following equation.
[0073]
(Equation 10)

[0074]
Therefore, when comparing the observed waveform with the calculation result using the FDTD method, the following convolution operation is performed on the actual measurement result.
i₀= H (t) * i (t) (19)
FIG. 10 shows an analysis result according to the present embodiment in a form to be compared with an actual measurement result obtained by performing the calculation. Similarly, FIGS. 11 and 12 show the analysis results by the conventional methods 1 and 2 in comparison with the actual measurement results. FIG. 13 shows a comparison of the analysis result of the present embodiment with the analysis result of the conventional method 1.
[0075]
As is clear from FIGS. 10 to 13, the conventional method 1 in which the conductors 701 are arranged on the grid, that is, the conductors existing in the analysis space are optimally arranged, is the conventional method 2 in which the conductors 701 are approximated stepwise. It can be seen that a high accuracy is obtained in the analysis result in comparison. It can be seen that the analysis result of the present embodiment has greatly improved accuracy compared to that of the conventional method 2, and an accuracy close to that of the conventional method 1 can be obtained. Since an accuracy close to that of the conventional method 1 is obtained, it was confirmed that when the present invention is applied to electromagnetic field analysis, a high accuracy that cannot be obtained by the conventional method can be achieved without disposing the conductor 701 on the grid. .
<Second embodiment>
In the first embodiment, among the divided surfaces divided by conductors represented (simulated) by lines, those that do not include the position where the magnetic field strength H is calculated are virtually added to adjacent surfaces. . Since the divided surface added in this manner is smaller than the adjacent surface (having the same size as the original surface here), the magnetic flux penetrating therethrough is theoretically approximated to zero. Thereby, for example, in the case shown in FIG. 3A, the magnetic flux passing through the surface ABEF is set to zero, and in the case shown in FIG. 3B, the magnetic flux passing through the surface GEF is set to zero. However, in such an approximation, a relatively large error may occur depending on the arrangement of the conductors. In the second embodiment, the error can be further suppressed. Hereinafter, portions different from the first embodiment will be described in detail.
[0076]
As shown in FIGS. 3A and 3B, a line representing a conductor crosses two opposing sides or two adjacent sides. The explanation in the latter is easier than that of the former. Therefore, different points from the first embodiment will be described focusing on the former.
[0077]
In the second embodiment, when a line representing a conductor crosses two opposing sides, the midpoint of the line is considered. Thereby, in the case as shown in FIG. 3A, as shown in FIG. 14A, the midpoint L of the line AF representing the conductor is considered. Points K and J are points on the sides BH and EG at the same position in the y direction as the midpoint L. The point I is a point on the side BH where the position in the y direction is the same as the point F.
[0078]
When the line AF (hereinafter also referred to as “conductor AF”) is a conductor, the magnetic flux passing through the surfaces JAL and KFL having the triangular shape have opposite directions and the same magnitude. In order to represent this conductor AF in the analysis space, the magnetic field strength H calculated using Faraday's law on the integration path ABCDEFA including the conductor AF is used._z磁 界 (i + 1/2, j-1 / 2, k) and the magnetic field strength H similarly calculated using Faraday's law in the integration path ABCDEKJA._zThe condition that ＋１ (i + 1/2, j-1 / 2, k) is equal is used.
[0079]
The surfaces ABCDEF and ABCDEKJ have the same area because the areas of the surfaces JAL and KFL are equal. It is for this reason that the midpoint L is further taken into account. Accordingly, another surface having the same area and the same magnetic field strength H as the surface ABCDEF is used, and the magnetic field strength H in the area, that is, the area obtained by adding the area of the surface ABCEF to the area of the surface BCDE is calculated with higher accuracy. I am trying to do it.
[0080]
Field strength H of surface BCDE_z(I + ／, j− ／, k) is calculated by the following equation using Faraday's law.
[0081]
[Equation 11]

[0082]
Here, in the second item on the right side,
-K_3z× 1 / (area of surface BCDE) × (circular integral value of electric field intensity E on surface BCDE)
Is established. Next, using this relationship, the magnetic field strength H is calculated for each of the integration paths ABCDEFA and ABCDEJA using Faraday's law._z(I + 1/2, j-1 / 2, k) is calculated, and the magnetic field strength H calculated by each integration path is calculated._zAccording to the condition that (i + 1/2, j-1 / 2, k) are equal, the electric field intensity E_y, E_xExpression (20) can be derived by performing central difference of in the x and y directions. However, from the positional relationship with the conductor AF, the electric field strength E_xThe integral value on the line KL and the integral value on the line LJ of are equal in magnitude and opposite in sign, and the electric field strength E_xThe integral value on the line KL of becomes zero (zero). Also, the electric field strength E_yThe integral value on the line JA of becomes equal to the integral value on the line KF.
[0083]
(Equation 12)

[0084]
Further, since the magnetic field strengths H calculated by the above equations (19) and (20) are equal, the electric field strength E on the side BE is_x(I + ／, j, k) is calculated from the following equation.
[0085]
(Equation 13)

[0086]
In the magnetic field strength H derived as described above, all the areas to be considered are considered. As a result, the magnetic flux corresponding to the virtually added area of the divided surface is reflected in the derivation of the magnetic field strength H. For this reason, it becomes more accurate as compared with the one derived from the equation (15). Since the equation (21) is an equation obtained on the assumption that the magnetic field strength H is more accurately derived, the derived electric field strength E is also more accurate.
[0087]
By the way, from the Faraday's law that applies the central difference only in the time domain, the expression (20) shows that the magnetic field strength H_z(I + 1/2, j-1 / 2, k) is derived. Electric field strength E_xFor 差分, a central difference including a first-order truncation error in the y direction is used. Equation (20) showing the contents of equation (20) by changing the notation is shown in equation (22), and the central difference is shown in equation (23).
[0088]
[Equation 14]

[0089]
The fact that the equation (23) includes a first-order truncation error can be proved as follows.
First, consider the following Taylor expansion.
[0090]
(Equation 15)

[0091]
Here, u = 1 + 2EK / Δy. From Equation (25)-Equation (24)
[0092]
(Equation 16)

[0093]
It can be seen that a first-order truncation error occurs. In the equation (26), “δ” is Kronecker delta, and “ο (2)” is a symbol representing a second-order or higher-order term.
[0094]
In the second embodiment, the first-order truncation error is further reduced to the second-order as follows. First, consider the following new Taylor development.
[0095]
[Equation 17]

[0096]
From equations (27) and (28),
[0097]
(Equation 18)

[0098]
Also, from the expression (26) + the expression (29),
[0099]
[Equation 19]

[0100]
It becomes. From the above,
[0101]
(Equation 20)

[0102]
Is derived. By using the equation (31) which is implicitly derived in this way, the error caused by approximating the magnetic flux penetrating the area ABEF, which cannot be avoided by the equation (15), to about zero is included in the equation (20). The first-order truncation error is reduced. Thereby, a more accurate analysis model can be obtained as compared with the first embodiment. In the case shown in FIG._z ^{n + 1/2}(I + 1/2, j + 1/2, k) is calculated by the same method.
[0103]
In the case where the above-described calculation is performed by the simulation apparatus illustrated in FIG. 5, the calculation of the magnetic field strength H using Expression (31) may be performed in, for example, step S9 in the electromagnetic field analysis processing illustrated in FIG. The correction of the electric field strength E using the equation (21) may be performed, for example, in step S4. For this reason, the realization of the second embodiment in the simulation apparatus shown in FIG. 5 can be performed basically without changing the operation flow from that in the first embodiment.
[0104]
Next, an analysis result performed using the simulation apparatus according to the second embodiment, which can be realized as described above, will be described. The analysis is the same as that of the first embodiment, that is, the case where a voltage waveform as shown in FIG. 9 is applied from the PG 703 to the conductor 701 in the analysis space analysis space where the conductor is present as shown in FIG. Analysis for obtaining the current and voltage waveforms at the application end was performed as verification analysis. Various conditions and the like are all the same as in the first embodiment, except for the calculation model. FIG. 15 shows a comparison between the analysis result obtained in this way and the analysis result obtained by the conventional method 1. In addition, the propagation time (propagation @ time) [ns] from the time when the current propagated from the starting end of the conductor 701 is reflected at the end and returns to the starting end, including the conventional method 1, the models of the first embodiment, and the measurement results. FIG. 16 shows the current peak value (Peak value of current) [A].
[0105]
As is clear from FIG. 15, FIG. 16 and FIG. 13, the analysis result agrees very well with that of the conventional method 1. It can be seen that the accuracy is greatly improved as compared with the first embodiment. From these facts, it can be confirmed that the model adopted in the second embodiment is very effective in improving accuracy. Although not particularly shown, it has been confirmed that the accuracy is improved even when the expression (20) is used instead of the expression (31) as compared with the first embodiment.
[0106]
In the second embodiment, in order to accurately handle the area of a division surface virtually divided by a conductor and virtually added to an adjacent surface, in other words, in order to obtain higher accuracy, the midpoint of the conductor (see FIG. In FIG. 14A, “K” and “J” are also considered in addition to “L”), but another point may be considered. For example, in the case shown in FIG. 14A, points A and F may be considered instead of point L.
[0107]
In the present embodiment (first and second embodiments), a straight line is used as a line for expressing (simulating) a conductor, but may be a curve. The object represented by the line is not limited to a conductor or a structure, and may be selected as needed from media existing in the analysis space. The method for approximating the electric field strength E on the side crossed by the conductor is not limited to the present embodiment, and various modifications can be made.
[0108]
【The invention's effect】
As described above, according to the present invention, a medium (for example, a conductor or a structure) selected from the media existing in the analysis space is set (represented) as a line, and the line is set on the surface constituting the cell. A plane crossing the side is divided into a plurality of parts, and at least one of an electric field strength and a magnetic field strength is obtained.
[0109]
By dividing the surface by a plurality of lines and obtaining at least one of the electric field intensity and the magnetic field intensity, the line is treated as being arranged on the side of the cell surface. Since the line may be arranged on the surface of the cell, a high degree of freedom is obtained, and the versatility is improved. By realizing a high degree of freedom, the shape of the medium represented by a line can be more accurately simulated. From these facts, the electromagnetic field analysis can be reliably performed with high accuracy while realizing high versatility.
[0110]
The division of the surface by the line is divided into a plurality of division surfaces obtained by dividing the surface by the line.A division surface not including the position for calculating the magnetic field strength is added to the adjacent surface, and the division surface including the position is If the processing is performed in such a manner that it is not added to adjacent planes, it is possible to avoid an increase in the number of planes due to the division. As a result, it is possible to avoid an increase in the number of finally obtained intensities, that is, an amount of data to be saved as an analysis result, and to suppress an increase in necessary computer resources.
[0111]
When a virtual change in area is caused by adding a dividing surface to an adjacent surface, if the change is reflected in the calculation of electric field strength or magnetic field strength, the error caused by the change is It can be avoided or reduced. For this reason, the intensity that reflects the change in the calculation can be obtained with higher accuracy. The change in intensity propagates between adjacent surfaces. As a result, the strength required on the nearby surface does not usually have a relatively large difference if the medium is the same. Therefore, when the reflection is performed using the electric field strength or the magnetic field strength obtained on the surface located near the surface where the change in the area occurs, it is difficult to avoid or reduce the error caused by the change. It can be performed more reliably. As a result, the calculated intensity can be obtained with higher accuracy.
[Brief description of the drawings]
FIG. 1 is a diagram illustrating cells that divide an analysis space.
FIG. 2 is a diagram illustrating a time arrangement of an electromagnetic field.
FIG. 3 is a diagram illustrating a positional relationship between a surface constituting a cell and a conductor crossing a side of the surface.
FIG. 4 is a circuit configuration diagram of a simulation device according to the present embodiment.
FIG. 5 is a diagram showing a functional configuration of an analysis program.
FIG. 6 is a flowchart of an electromagnetic field analysis process.
FIG. 7 is a diagram showing a positional relationship between conductors.
FIG. 8 is a diagram illustrating the arrangement of conductors.
FIG. 9 is a diagram showing a time change of a voltage waveform of a pulse generator.
FIG. 10 is a diagram showing an analysis result according to the first embodiment.
FIG. 11 is a diagram showing an analysis result by the conventional method 1.
FIG. 12 is a diagram showing an analysis result by the conventional method 2.
FIG. 13 is a diagram showing an analysis result according to the first embodiment in comparison with an analysis result according to the conventional method 1.
FIG. 14 is a diagram for explaining a method of calculating a magnetic field intensity and an electric field intensity according to the second embodiment.
FIG. 15 is a diagram showing an analysis result according to the second embodiment.
FIG. 16 is a diagram showing a propagation time and a current peak value in each model.
[Explanation of symbols]
301 CPU
302 @ ROM
303 RAM
304 communication interface
305 storage device
306 Storage medium reading device
307 input device
308 output device
320, 330 recording medium
500 analysis program
501 input section
502 Analysis unit
503 output section
504 Boundary condition processing unit
505 ° electric field analysis unit
506 Magnetic field analysis unit
507 Electric field strength calculator
508 ° electric field strength correction unit
509 magnetic field strength calculator
510 ° magnetic field strength correction unit
701 conductor
702 copper plate
703 pulse generator

Claims (15)

A method for dividing an analysis space into hexahedral cells and performing an electromagnetic field analysis using an FDTD (Finite Difference Time Domain) method. Set the medium as a line,
Of the surfaces constituting the cell, the line intersecting the side is divided into a plurality of lines to determine the electric field strength, or the magnetic field strength,
An electromagnetic field analysis method using the FDTD method. The division of the plane by the line may be performed by dividing a plane including the position where the magnetic field intensity is calculated among a plurality of division planes obtained by dividing the plane by the line, in addition to an adjacent plane, and a division plane including the position. Do not add to the adjacent surface,
An electromagnetic field analysis method using the FDTD method according to claim 1. For the electric field intensity of the surface divided into a plurality of lines, the direction calculated around the surface adopts the same electric field intensity,
An electromagnetic field analysis method using the FDTD method according to claim 1. The magnetic field strength of the surface divided into a plurality of lines by the line is calculated using the electric field strength obtained by the integration path including the line,
An electromagnetic field analysis method using the FDTD method according to claim 1, 2, or 3. The magnetic field strength of a surface that is adjacent in a direction that intersects with the surface that is divided into a plurality of lines by the line, and that has a side that intersects the line, is calculated on the assumption that the surface is adjacent to the plurality of divided surfaces.
An electromagnetic field analysis method using the FDTD method according to claim 2 or 4. The electric field strength of the surface divided into a plurality of lines by the line, after calculating all the electric field strength without considering the line, is obtained in a form to correct the calculated electric field strength,
An electromagnetic field analysis method using the FDTD method according to claim 1 or 3. On a surface where a virtual change occurs in area by adding the divided surface to an adjacent surface, the change is reflected in the calculation of the electric field strength or the magnetic field strength,
An electromagnetic field analysis method using the FDTD method according to claim 2. The reflection is performed using an electric field strength or a magnetic field strength obtained on a surface located near the surface where the change occurs,
An electromagnetic field analysis method using the FDTD method according to claim 7, wherein: A method used to represent a medium in an electromagnetic field analysis performed by dividing an analysis space into hexahedral cells and performing an FDTD (Finite Difference Time Domain) method,
Selecting a medium present in the analysis space to be divided into the cells,
The selected medium is set and expressed as a line crossing a side of a surface constituting the cell,
A method for expressing a medium in electromagnetic field analysis. In a simulation apparatus that divides an analysis space into hexahedral cells and performs an electromagnetic field analysis using an FDTD (Finite Difference Time Domain) method,
Medium setting means for setting a medium existing in the analysis space as a line,
Analysis means for dividing the plane intersecting the side among the surfaces of the cells dividing the analysis space by the line set by the medium setting means, and for determining the electric field strength or the magnetic field strength,
A simulation device comprising: The division of the plane by the line may be performed by dividing a plane including the position where the magnetic field intensity is calculated among a plurality of division planes obtained by dividing the plane by the line, in addition to an adjacent plane, and a division plane including the position. Do not add to the adjacent surface,
The simulation device according to claim 10, wherein: The analysis means, as the electric field intensity of the surface divided into a plurality of lines, the direction calculated around the surface adopts the same electric field intensity,
The simulation device according to claim 10, wherein: The analysis means, on a surface where a virtual change occurs in the area by adding the division surface to an adjacent surface, reflects the change in the calculation of the electric field strength or the magnetic field strength,
The simulation device according to claim 11, wherein: The reflection is performed using an electric field strength or a magnetic field strength obtained on a surface located near the surface where the change occurs,
14. The simulation device according to claim 13, wherein: A program to divide an analysis space into hexahedral cells and to execute the data processing device used as a simulation device for performing an electromagnetic field analysis by an FDTD (Finite Difference Time Domain) method,
A function of setting a medium existing in the analysis space as a line,
A function for obtaining the electric field strength or the magnetic field strength by dividing the plane intersecting the side among the cell planes dividing the analysis space by the line set by the setting function,
The program to realize.

Images