(57)【特許請求の範囲】
【請求項１】
解析空間を６面体のセルに分割し、ＦＤＴＤ（Finite Difference Time Domain ）法を用いて電磁界解析を行うための方法であって、
前記セルに分割する前記解析空間に存在する媒質のなかから選択した媒質を線として設定し、
前記セルを構成する面のなかで前記線が辺を横切る面を該線で複数に分割して電界強度、或いは磁界強度を求める、
ことを特徴とするＦＤＴＤ法を用いた電磁界解析方法。
【請求項２】
前記線による面の分割は、該面を該線で分けて得られる複数の分割面のなかで前記磁界強度を計算する位置を含まない分割面は隣接する面に加え、該位置を含む分割面は隣接する面に加えない形で行う、
ことを特徴とする請求項１記載のＦＤＴＤ法を用いた電磁界解析方法。
【請求項３】
前記線で複数に分割される面の電界強度には、該面の周辺で計算した方向が同じ電界強度を採用する、
ことを特徴とする請求項１記載のＦＤＴＤ法を用いた電磁界解析方法。
【請求項４】
前記線で複数に分割される面の磁界強度は、該線を含む積分路で求めた電界強度を用いて計算する、
ことを特徴とする請求項１、２、または３記載のＦＤＴＤ法を用いた電磁界解析方法。
【請求項５】
前記線で複数に分割される面と交差する方向に隣接し、且つ該線が横切る辺を有する面の磁界強度は、複数の分割面と隣接していると想定して計算する、
ことを特徴とする請求項２、または４記載のＦＤＴＤ法を用いた電磁界解析方法。
【請求項６】
前記線で複数に分割される面の電界強度は、該線を考慮せずに電界強度を全て計算した後、計算した電界強度を補正する形で求める、
ことを特徴とする請求項１、または３記載のＦＤＴＤ法を用いた電磁界解析方法。
【請求項７】
前記分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面では、該変化を前記電界強度、或いは磁界強度の計算に反映させる、
ことを特徴とする請求項２〜５の何れか１項に記載のＦＤＴＤ法を用いた電磁界解析方法。
【請求項８】
前記反映は、前記変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行う、
ことを特徴とする請求項７記載のＦＤＴＤ法を用いた電磁界解析方法。
【請求項９】
解析空間を６面体のセルに分割し、ＦＤＴＤ（Finite Difference Time Domain ）法により行う電磁界解析で媒質を表現するために用いられる方法であって、
前記セルに分割する前記解析空間に存在する媒質を選択し、
該選択した媒質を、前記セルを構成する面の辺を横切る線として設定し表現する、
ことを特徴とする電磁界解析における媒質表現方法。
【請求項１０】
解析空間を６面体のセルに分割し、ＦＤＴＤ（Finite Difference Time Domain ）法を用いて電磁界解析を行うシミュレーション装置において、
前記解析空間に存在する媒質を線として設定する媒質設定手段と、
前記媒質設定手段により設定した前記線が前記解析空間を分割する前記セルの面のなかで辺を横切る面を該線で分割して電界強度、或いは磁界強度を求める解析手段と、
を具備することを特徴とするシミュレーション装置。
【請求項１１】
前記線による面の分割は、該面を該線で分けて得られる複数の分割面のなかで前記磁界強度を計算する位置を含まない分割面は隣接する面に加え、該位置を含む分割面は隣接する面に加えない形で行う、
ことを特徴とする請求項１０記載のシミュレーション装置。
【請求項１２】
前記解析手段は、前記線で複数に分割される面の電界強度として、該面の周辺で計算した方向が同じ電界強度を採用する、
ことを特徴とする請求項１０記載のシミュレーション装置。
【請求項１３】
前記解析手段は、前記分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面では、該変化を前記電界強度、或いは磁界強度の計算に反映させる、
ことを特徴とする請求項１１記載のシミュレーション装置。
【請求項１４】
前記反映は、前記変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行う、
ことを特徴とする請求項１３記載のシミュレーション装置。
【請求項１５】
解析空間を６面体のセルに分割し、ＦＤＴＤ（Finite Difference Time Domain ）法により電磁界解析を行うシミュレーション装置として用いられるデータ処理装置に実行させるプログラムであって、
前記解析空間に存在する媒質を線として設定する機能と、
前記設定する機能により設定した前記線が前記解析空間を分割する前記セルの面のなかで辺を横切る面を該線で分割して電界強度、或いは磁界強度を求める機能と、
を実現させるためのプログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、解析空間を６面体のセルに分割し、ＦＤＴＤ（Finite Difference Time Domain ）法を用いて電磁界解析を行うための技術に関する。
【０００２】
【従来の技術】
電磁界解析は、アンテナ問題や電磁波散乱問題を始めとする各種分野で行われている。その解析手法として、有限要素法やモーメント法、或いは境界要素法などが用いられていたが、近年、Ｋ．Ｓ．Ｙｅｅにより提案された有限差分法であるＦＤＴＤ（Finite Difference Time Domain ）法が注目されている。これは、アルゴリズムが簡単であること、優れた精度を持つこと、適用範囲が広いこと、などの長所を有しているためである。
【０００３】
そのＦＤＴＤ法は、電磁界を時間領域と解析空間において差分化し、解析対象物体の幾何学的な配置や形状、導電率、透磁率等の電気的（物理的）特性から電磁界の挙動を数値解析により求める手法である。時間領域では微小時間（Δｔ）で離散化を行い、解析空間はセルと呼ばれる単位に分割して離散化を行う。
【０００４】
図１は、解析空間を分割するセルを説明する図である。図１に示すように、セルは、ｘ軸方向はΔｘ、ｙ軸方向はΔｙ、ｚ軸方向はΔｚの長さを持つ６面体（通常は立方体）として設定される。磁界強度は、各面の中心に配置され、電界強度は、面を構成する各辺の中心に配置されている。そのようにして、磁界強度と電界強度を空間的に半セル分ずらして差分化することにより、電界の回転が磁界をつくり、磁界の回転が電界をつくる、というマクスウェルの方程式を直接的に解ける配置となっている。セルサイズは、一般的には問題とする最短波長の１／１０以下に設定される。
【０００５】
空間的にずらして差分化される電界と磁界は、図２に示すように、時空間では交互に配置される。電界強度と磁界強度を交互に計算するリープフロッグアルゴリズムが用いられている。それにより、電界強度から磁界強度を求め、磁界強度から電界強度を求めるというようにそれらは交互に計算される。図２中に表記された「ｎ−１／２」や「ｎ」などの上添字は、時間領域上で配置された位置を表している。時間領域の差分は、Courant 安定条件を満たすように行われる。
【０００６】
ＦＤＴＤ法では、上述したように、時間領域と解析空間の差分化に中心差分法が使われている。解析を行う分野によって、電界と磁界の時間軸、空間座標中の配置が逆となる場合もある。
【０００７】
ＦＤＴＤ法の基本となるマクスウェルの方程式は、アンペアの周回積分の法則とファラデーの電磁誘導の法則を基本とする。それぞれの法則は微分形と積分形とで表され、上記マクスウェルの方程式は、通常、電界、磁界に関するガウスの法則と組み合わせた計４つの方程式のことを言う。微分形のマクスウェルの方程式は、電界強度Ｅ［Ｖ／ｍ］、磁界強度Ｈ［Ａ／ｍ］、電束密度Ｄ［Ｃ／ｍ2 ］、磁束密度Ｂ［Ｔ］、電荷密度ρ［Ｃ／ｍ3 ］、電流密度Ｊ［Ａ／ｍ2 ］を用いて次のように表される。
【０００８】
【数１】
0003703812.tif 000002
【０００９】
ここで（１）式（アンペアの法則）、（２）式（ファラデーの電磁誘導の法則）、（３）式（磁界中のガウスの法則）、（４）式（電界中のガウスの法則）は独立な式ではなく、ＦＤＴＤ法において（３）、（４）式は数値誤差の評価基準として用いられる程度で、定式化には（１）、（２）式が用いられる。（１）、（２）式にＹｅｅのアルゴリズムを用いて、時空間と空間座標について中心差分法を適用し、定式化すると、３次元空間において次のようなＦＤＴＤ基本方程式を得ることができる。
【００１０】
【数２】
0003703812.tif 000003
【００１１】
ただし、μ：透磁率、ε：誘電率、σ：抵抗率である。
（５）、（６）式中の「Ｅ」や「Ｈ」の各シンボルに下添字として付した「ｘ」や「ｙ」は、それぞれ、電界強度、磁界強度の方向を表している。ここでは、ｘ方向に割り当てた電界強度、磁界強度のみを示したが、ｙ、ｚ方向についても同様に導出できる。
【００１２】
【非特許文献１】
羽野、「ＦＤＴＤ（Finite Difference Time Domain ）法」、平成１４年電気学会全国大会論文集、第５冊、ｐｐ４１１−４１４、２００２
【非特許文献２】
野田、横山、「ＦＤＴＤ法に基づく汎用サージ解析プログラムの開発」、電学論Ｂ、１２１巻５号、ｐｐ６２５−６３２、２００１
【非特許文献３】
J.Fang and J.Ren、“A Locally Conformed Finite-Difference Time-Domain Algorithm of Modeling Arbitrary Shape Planar Metal Strips ”、IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES、VOL.41、NO.5、MAY 1993
【００１３】
【発明が解決しようとする課題】
ＦＤＴＤ法では、図１に示すようなセルで解析空間を分割するのを基本とする。このため、そのセルの格子（辺）上に存在しない導体などの解析対象物体を表現する場合には、その形状を階段状に近似する階段近似法が用いられていた。しかし、階段近似法で表現される形状と実際の形状との間には大きな違いがあるのが普通である。このため、階段近似法を用いて行う従来の解析方法では、解析対象物体の共振周波数やその伝搬時間などに大きな誤差が生じて、解析を高精度に行えないという問題点があった。
【００１４】
セルサイズを小さくすることにより、階段近似法で表現される形状と実際の形状との間の違いを小さくすることができる。しかし、そのサイズを小さくすると、それに伴いセル数が増大することから、必要な計算機資源や計算時間が膨大となってしまう。このため、部分的であってもセルサイズを小さくすることは望ましくない。サージ解析の分野では、比較的に大きな解析空間を対象にすることが多いことから、言い換えればセル数が大きいのが普通であることから、特に望ましくないと言える。
【００１５】
セルの形状を変形させて解析対象物体の形状を表現する場合には、その物体に応じたモデリングが必要になり、一般性に欠ける。その物体が格子上となるように座標系を設定する方法は、複数の物体が任意の状態で解析空間内に存在していれば用いることができないのが普通である。このようなことから、セルサイズを小さくする、セルの形状を変形させるといった操作を行うことなく、解析空間内に存在する解析対象物体をより正確に表現して解析を行えるようにすることが重要であると考えられる。
【００１６】
本発明は、解析空間内に存在する解析対象物体（媒質）を常により正確に表現できるようにして、電磁界解析を確実に高精度に行えるようにするための汎用的な技術を提供することを目的とする。
【００１７】
【課題を解決するための手段】
本発明のＦＤＴＤ法を用いた電磁界解析方法は、解析空間を６面体のセルに分割し、ＦＤＴＤ（Finite Difference Time Domain ）法を用いて電磁界解析を行うための方法であって、セルに分割する解析空間に存在する媒質のなかから選択した媒質を線として設定し、セルを構成する面のなかで線が辺を横切る面を該線で複数に分割して電界強度、或いは磁界強度を求める。
【００１８】
なお、上記線による面の分割は、該面を該線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、該位置を含む分割面は隣接する面に加えない形で行う、ことが望ましい。線で複数に分割される面の電界強度には、該面の周辺で計算した方向が同じ電界強度を採用する、ことが望ましい。線で複数に分割される面の磁界強度は、該線を含む積分路で求めた電界強度を用いて計算する、ことが望ましい、線で複数に分割される面と交差する方向に隣接し、且つ該線が横切る辺を有する面の磁界強度は、複数の分割面と隣接していると想定して計算する、ことが望ましい。線で複数に分割される面の電界強度は、該線を考慮せずに電界強度を全て計算した後、計算した電界強度を補正する形で求める、ことが望ましい。
【００１９】
また、分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面では、該変化を電界強度、或いは磁界強度の計算に反映させる、ことが望ましい。その反映は、変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行う、ことが望ましい。
【００２０】
本発明の電磁界解析における媒質表現方法は、解析空間を６面体のセルに分割し、ＦＤＴＤ法により行う電磁界解析で媒質を表現するために用いられる方法であって、セルに分割する解析空間に存在する媒質を選択し、該選択した媒質を、セルを構成する面の辺を横切る線として設定し表現する。
【００２１】
本発明のシミュレーション装置は、解析空間を６面体のセルに分割し、ＦＤＴＤ法を用いて電磁界解析を行うことを前提とし、解析空間に存在する媒質を線として設定する媒質設定手段と、媒質設定手段により設定した線が解析空間を分割するセルの面のなかで辺を横切る面を該線で分割して電界強度、或いは磁界強度を求める解析手段と、を具備する。
【００２２】
なお、上記線による面の分割は、該面を該線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、該位置を含む分割面は隣接する面に加えない形で行う、ことが望ましい。解析手段は、線で複数に分割される面の電界強度として、該面の周辺で計算した方向が同じ電界強度を採用する、ことが望ましい。
【００２３】
また、解析手段は、分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面では、該変化を電界強度、或いは磁界強度の計算に反映させる、ことが望ましい。その反映は、変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行う、ことが望ましい。
【００２４】
本発明のプログラムは、解析空間を６面体のセルに分割し、ＦＤＴＤ法により電磁界解析を行うシミュレーション装置として用いられるデータ処理装置に実行させることを前提とし、解析空間に存在する媒質を線として設定する機能と、設定する機能により設定した線が解析空間を分割するセルの面のなかで辺を横切る面をその線で分割して電界強度、或いは磁界強度を求める機能と、を実現させる。
【００２５】
本発明では、解析空間に存在する媒質のなかから選択した媒質（導体や構造物など）を線として設定（表現）し、セルを構成する面のなかでその線が辺を横切る面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求める。
【００２６】
その線で面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求めることにより、その線を、セルの面の辺上に配置した形で扱うことになる。その線はセルの面上に配置すれば良いから、高い自由度が得られ、汎用性は向上する。その線で表現する媒質は形状をより正確に模擬することが可能となる。これらのことから、高い汎用性を実現させつつ、電磁界解析は確実に高精度に行えるようになる。
【００２７】
上記線による面の分割を、その面をその線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、その位置を含む分割面は隣接する面に加えない形で行うようにした場合には、その分割に伴い面の数が増大することは回避される。その結果、最終的に求める強度の数、即ち解析結果として保存の対象となるデータの量の増大は回避され、必要な計算機資源の増大は抑えられる。
【００２８】
分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面において、その変化を電界強度、或いは磁界強度の計算に反映させるようにした場合には、その変化に伴って生じる誤差は回避、或いは低減させることが可能となって、その変化を計算に反映させる強度はより高精度に得られるようになる。強度の変化は、隣接する面間で伝搬していく。それにより、近傍の面で求められる強度には、媒質が同じであれば比較的に大きな違いがないのが普通である。このことから、その反映を面積の変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行うようにしたときには、計算する強度は更に高精度に得られるようになる。
【００２９】
【発明の実施の形態】
以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
＜第１の実施の形態＞
図３は、セルを構成する面とその面の辺を横切る導体の位置関係を説明する図である。始めに、図３を参照して、本実施の形態で用いられる導体（他と区別すべき媒質）の表現方法、その導体周辺における磁界と電界の近似方法について詳細に説明する。
【００３０】
本実施の形態では、解析空間は一般的な空間分割手法により６面体のセル（図１参照）に分割し、導体の形状や配置は別に設定する。それにより、セルを構成する面上の任意の位置に任意の角度で傾斜した導体の形状を線で模擬する。セルの格子上に配置できない導体を面上の線で模擬する。
【００３１】
そのように線で模擬した場合、辺を横切る導体は、それと対向する辺を横切るか、或いはそれと隣接する辺を横切ることになる。図３（ａ）は前者、即ち対向する２辺を導体が横切る場合の位置関係を示し、同図（ｂ）は後者、即ち隣接する２辺を導体が横切る場合の位置関係を示している。同図（ｂ）では、面ＡＣＤＥ以外を通る導体は破線で示してある。面上の磁界強度を計算する位置は、黒丸を内部に持つ円で示してある。
【００３２】
従来、ＦＤＴＤ法では、各辺で１つの電界だけを定義する。しかし、本実施の形態では、導体が横切る辺では、それらが交差する点を境に２つ定義する。導体に隣接する磁界の積分路は、その面の辺に沿った積分路から導体を含む積分路に変更する。そのようにして、電界強度や磁界強度を計算する面を導体により分割し、分割した面（分割面）別にそれらの強度を計算する形で格子上以外に存在する導体を考慮する。以下、図３（ａ）、図３（ｂ）に示すケース別に、導体周辺の磁界と電界の近似方法について具体的に説明する。
【００３３】
図３（ａ）に示すケースでは、導体が辺ＢＨ上の点Ａ、及び辺ＥＧ上の点Ｆで交差している。このことから、このケースでは、辺ＢＨは辺ＡＢ、ＡＨ、辺ＥＧは辺ＥＦ、ＦＧにそれぞれ２つに分けて、その中心にあるｙ方向の電界強度Ｅyを定義する。ここでは、それらの電界強度Ｅy を以下のように定義する。
【００３４】
【数３】
0003703812.tif 000004
【００３５】
（１０）、（１１）式中の括弧内は、電界強度Ｅy を定義した座標位置を表している。シンボル「Δｙ」はｙ方向上のセルサイズを表している。「ＡＢ」は、点Ａと点Ｂの間の長さを表している。それらは他の式でも同様である。
【００３６】
図３（ｂ）に示すケースでは、導体が辺ＡＣ上の点Ｂ、及び辺ＡＥ上の点Ｇで交差している。このことから、このケースでは、辺ＡＣは辺ＡＢ、ＢＣの２つに分け、その中心にあるｙ方向の電界強度Ｅy を同様に以下のように定義する。
【００３７】
【数４】
0003703812.tif 000005
【００３８】
一方、導体の近傍に割り当てられたｘ方向の電界強度Ｅx は、以下のように計算する。
図３（ａ）に示すケースでは、積分路ＡＢＣＤＥＦＡと積分路ＢＣＤＥＢを用いて計算した磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）が同一でなければならないという関係から、電界強度Ｅx （ｉ＋１／２，ｊ，ｋ）は次式を用いて計算する。
【００３９】
【数５】
0003703812.tif 000006
【００４０】
同様に図３（ｂ）に示すケースでは、積分路ＢＣＤＥＦＧＢと積分路ＢＣＤＥＧＢを用いて計算した磁界強度Ｈz （ｉ＋１／２，ｊ＋１／２，ｋ）が同一でなければならないという関係から、電界強度Ｅx （ｉ＋１−ＥＧ／２Δｙ，ｊ＋１，ｋ）は次式を用いて計算する。
【００４１】
【数６】
0003703812.tif 000007
【００４２】
電界強度Ｅx （ｉ＋ＡＧ／２Δｘ，ｊ＋１，ｋ）についても同様の式で計算できる。
（１３）、（１４）式を導出する際に考慮したように、導体の周辺にある磁界強度は、その積分路を導体に沿った形に変更し計算する。図３（ａ）に示すケースでは、磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）を計算する際の積分路は積分路ＡＢＣＤＥＦＡとする。導体の電界強度Ｅはゼロであるから、磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）は次式を用いて計算される。
【００４３】
【数７】
0003703812.tif 000008
【００４４】
導体が横切る面ＢＥＧＨを構成する面（分割面）ＡＢＥＦはその内部に磁界強度Ｈz を計算する位置を含まない。その分割面ＡＢＥＦで磁界強度Ｈz を計算すると、計算して残すべき磁界強度Ｈの数が面の分割によって増大する。その管理には煩雑な処理が必要となる。分割面ＡＢＥＦを面ＢＣＤＥに加えて、磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）を計算する際の積分路を積分路ＡＢＣＤＥＦＡとしたのは、そのようなことを回避するためである。他方の分割面ＡＦＧＨでは、そのようなことが生じないことから、一つの面として扱う。
【００４５】
一方、磁界強度Ｈz （ｉ＋１／２，ｊ＋１／２，ｋ）も同様に、積分路ＡＦＧＨＡに沿って計算される。図３（ｂ）に示すケースでも同様に、磁界強度Ｈz （ｉ＋１／２，ｊ＋１／２，ｋ）を計算する際の積分路は積分路ＢＣＤＥＦＧＢとなり、磁界強度Ｈz （ｉ＋１／２，ｊ＋３／２，ｋ）を計算する際の積分路は積分路ＡＢＧＦＨＩＡとなる。
【００４６】
導体を含む面に垂直で、その導体が交わる格子（辺）を含む面の磁界強度Ｈは、その後に計算する。その計算は、（１０）〜（１４）式により補正された電界強度Ｅを用いて行う。例えば図３（ａ）に示すケースでは、磁界強度Ｈx （ｉ，ｊ＋１／２，ｋ＋１／２）は電界強度Ｅy （ｉ，ｊ＋ＡＢ／２Δｙ，ｋ）、Ｅy （ｉ，ｊ＋１−ＡＨ／２Δｙ，ｋ）を用いて次式により計算する。他の磁界強度Ｈx （ｉ，ｊ＋１／２，ｋ＋１／２）、Ｈx （ｉ＋１，ｊ＋１／２，ｋ−１／２）、及びＨx （ｉ＋１，ｊ＋１／２，ｋ−１／２）、図３（ｂ）に示すケースでは、磁界強度Ｈx （ｉ，ｊ＋１／２，ｋ＋１／２）、Ｈx （ｉ＋１，ｊ＋１／２，ｋ−１／２）、Ｈy （ｉ＋１／２，ｊ＋１，ｋ−１／２）、及びＨy （ｉ＋１／２，ｊ＋１，ｋ−１／２）についても同様の式で計算することができる。
【００４７】
【数８】
0003703812.tif 000009
【００４８】
格子を横切る導体と、格子上に近似された導体との接続は簡単に行うことができる。例えば図３（ａ）に示すケースで辺ＡＨは導体の内部に存在するのであれば、その辺ＡＨ上で定義する電界強度Ｅy をゼロとするだけで良い。このようなことから、解析空間に存在する任意の解析対象物体（構造物）を線状の導体として表現できる。
【００４９】
上述したように、例えば図３（ａ）に示すケースでは、磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）の計算には積分路ＡＢＣＤＥＦＡを用い、辺ＡＢ上の電界強度Ｅy は辺ＢＣ上の電界強度Ｅy と定義している。そのようにすると、導体で面を分割しても、計算する強度の数が増えるようなことは回避される。つまり、必要な計算機資源が増大するようなことを回避しつつ、解析対象物体を線として模擬できるようになる。プログラム開発の面では、計算した強度を代入するための配列変数を新たに用意する必要性が回避されることから、言い換えれば、計算した強度の管理は従来と同じ方法で行えることから、その開発が容易となる。既存の解析プログラムに本発明を適用する場合には、そのためのアップデートをより容易に行えることになる。
【００５０】
辺ＢＣ上の電界強度Ｅy を辺ＡＢ上の電界強度Ｅy と定義すると、辺ＢＣ上の電界強度Ｅy は辺ＡＣ上で求めた電界強度Ｅy として扱われることになる。しかし、ＦＤＴＤ法では中間差分が採用されていることもあって、そのように扱っても誤差は小さな範囲内に収まりやすい。磁界強度Ｈを上述したように計算し、電界強度Ｅを上述したように定義するのはそのような理由からである。計算機資源に余裕があるような場合には、導体を境に面を分割し、分割した面毎に電界強度、及び磁界強度の少なくとも一方を求めるようにしても良い。
【００５１】
導体を表現する線は、セルを構成する面上に配置しなければならない。しかし、格子上に配置する場合と比較すると、線を配置できる場所が直線上から平面上に変化したために、線を配置するうえでの制約は大幅に緩和して自由度は大きく向上することになる。直線でない線も扱えるようになる。このため、高い汎用性が得られる。たとえ階段近似を行う場合であっても、格子上に沿って線を配置する必要がないために、セルサイズを小さくすることなく、より小さな誤差で形状を表現することができる。
【００５２】
図４は、本実施の形態によるシミュレーション装置の回路構成図である。その装置は、上述したように導体を模擬して、例えばサージ解析を行うためのものである。
【００５３】
その装置は、例えばコンピュータ（データ処理装置）に本実施の形態による解析プログラムを実行させることで実現されるものである。図４に示すように、ＣＰＵ３０１、ＲＯＭ３０２、ＲＡＭ３０３、通信インターフェース３０４、例えばハードディスク装置である記憶装置３０５、可搬の記録媒体３３０の読み取りを行う記録媒体読み取り装置３０６、キーボードやポインティングデバイスなどを含めて表記した入力装置３０７、及び例えば不図示の表示装置に画像を出力する出力装置３０８がバス３０９により互いに接続されて構成されている。
【００５４】
上記解析プログラムは、記憶装置３０５、記録媒体３３０、或いは外部装置がアクセス可能な記録媒体３２０に格納されている。記録媒体３２０に解析プログラムが格納されている場合、ＣＰＵ３０１は、外部装置、ネットワーク、通信インターフェース３０４、及びバス３０９を介してその解析プログラムを取得し実行する。
【００５５】
図５は、上記解析プログラムの機能構成を示す図である。
図５中の形状ファイル５２１は、解析空間やその空間内に存在する解析対象物体などの形状や配置などの設定データを格納したファイルである。線状の導体として扱う構造物（媒質）は、例えば線の始点と終点の座標を示すデータで形状ファイル５２１中に設定されている。パラメータファイル５２２は、解析を行ううえで用いるべきパラメータの値などを格納したファイルである。電圧源や境界条件などに関する設定はそのファイル５２２を介して行われる。出力ファイル５２３は、解析結果を出力するファイルである。それらのファイル５２１〜５２３は、ユーザが指定する。
【００５６】
解析プログラム５００は、形状ファイル５２１、及びパラメータ５２２を入力ファイルとして、それに格納されているデータを用いて電磁界解析（サージ解析、など）を行い、その解析結果を出力ファイル５２３に出力するものとして実現されている。
【００５７】
解析プログラム５００の機能構成は、入力部５０１、解析部５０２，及び出力部５０２に大別される。入力部５０１は、形状ファイル５２１、及びパラメータファイル５２２に格納されたデータを取り込み、解析部５０２に提供することを主に行う。
【００５８】
解析部５０２は、境界条件処理部５０４、電界解析部５０５、及び磁界解析部５０６を備えている。境界条件処理部５０４は、設定された境界条件を電界解析部５０５、及び磁界解析部５０６に提供する。電界解析部５０５は、それが直前の時間ステップで計算した電界強度Ｅや磁界解析部５０６が計算した磁界強度Ｈを用いて現時点の時間ステップでの電界強度Ｅを計算する。他方の磁界解析部５０６は、それが直前の時間ステップで計算した磁界強度Ｈや電界解析部５０５が計算した電界強度Ｅを用いて現時点の時間ステップでの磁界強度Ｈを計算する。出力部５０３は、各解析部５０５、５０６が計算によりそれぞれ求めた強度Ｅ，Ｈを時刻に対応付けして出力ファイル５２３に出力する。
【００５９】
電界解析部５０５は、電界強度計算部５０７と電界強度補正部５０８を備える。電界強度計算部５０７は、セルの格子を導体が横切っていないことを想定して電界強度Ｅを計算する。電界強度補正部５０８は、電界強度計算部５０７が計算した電界強度Ｅを入力し、導体が格子を横切っている面の電界強度Ｅを上述したように補正する。電界強度Ｅを全て計算した後に補正を行うのは、導体が格子を横切る面では、その格子上の電界強度Ｅとして他の面の電界強度Ｅを定義していることから、他の面の更新を行う前の電界強度Ｅを定義してしまうことを確実に回避するためである。出力部５０３には、その補正を行った後の電界強度Ｅを出力する。
【００６０】
他方の磁界解析部５０６は、電界解析部５０５と同様に、磁界強度計算部５０９と磁界強度補正部５１０を備える。磁界強度計算部５０９は、セルの格子を導体が横切っていないことを想定して磁界強度Ｈを計算する。磁界強度補正部５１０は、磁界強度計算部５０９が計算した磁界強度Ｈを入力し、導体が格子を横切っている面の磁界強度Ｈを上述したように補正する。ここでも同様に、磁界強度Ｈを全て計算した後にその補正を行う。出力部５０３には、その補正を行った後の磁界強度Ｈを出力する。
【００６１】
次に、上記機能構成を有する解析プログラムをＣＰＵ３０１が実行することで実現される電磁界解析（例えばサージ解析）処理について、図６に示すそのフローチャートを参照して詳細に説明する。
【００６２】
先ず、ステップＳ１では、ユーザが指定した形状ファイル５２１、及びパラメータファイル５２２に格納されていうデータにアクセスすることにより、それらのファイル５２１、５２２で設定されている解析条件を読み込み、次の時間ステップでの計算を行ううえでの初期設定を行う。続くステップＳ２では、解析空間の持つ６つの面（それを囲む面）でそれぞれ定められた境界条件の処理を行う。その後はステップＳ３に移行する。
【００６３】
ステップＳ３では、解析空間全体で電界強度Ｅの計算を行う。その計算は、セルの格子を横切る導体が存在していないことを想定して行う。その次に移行するステップＳ４では、導体の格子の横切り方に応じて、（１０）〜（１４）式を用いて説明したように導体周辺の電界強度Ｅを補正する。その後に移行するステップＳ５では、格子上の導体に関する電界強度Ｅの補正を行い、導体内に位置する格子上の電界強度Ｅはゼロにする。ステップＳ６にはその後に移行する。
【００６４】
ステップＳ６では、解析空間内に存在する電圧源、電流源といった給電部の置かれた位置の電界強度Ｅを計算する。次のステップＳ７では、ステップＳ２〜Ｓ６を実行することで得られた電界強度Ｅや電圧値などの計算結果を出力ファイル５２３に出力する。その後はステップＳ８に移行する。
【００６５】
ステップＳ８では、解析空間全体で磁界強度Ｈを計算する。続くステップＳ９では、導体の格子の横切り方に応じて、（１５）、（１６）式を用いて説明したように導体周辺の磁界強度Ｈを補正する。その次に移行するステップＳ１０では、ステップＳ８、或いはＳ９を実行することで得た磁界強度Ｈや電流値などの計算結果を出力ファイル５２３に出力する。その後はステップＳ１１に移行して、解析を行った時間の長さである時刻ｔが最大観測時間Ｔmax より大きいか否か判定する。定められた最大観測時間Ｔmax 分の解析が終了した場合、判定はＹＥＳとなり、ここで電磁界解析処理を終了する。そうでない場合には、判定はＮＯとなって上記ステップＳ１に戻り、時間ステップを１つ進めて、それ以降の処理を同様に実行する。
【００６６】
次に、本実施の形態によるシミュレーション装置を用いて行った解析結果について説明する。
図７は、導体の位置関係を示す図である。解析空間内に存在する導体やそれらの間の位置関係などを示したものである。
【００６７】
その解析空間では、大地に水平に設置された導体７０１が、その大地を覆う銅板７０２から０．６ｍの位置に置かれている。その導体７０１は、長さ４［ｍ］、半径１５×１０-3［ｍ］のアルミパイプである。その導体７０１には、特性インピーダンス５０オームのケーブルに充電した電荷を水銀リレーで投入するタイプのパルスジェネレータ（ＰＧ）７０３が、断面積２［ｍｍ2 ］の銅線により接続されている。
【００６８】
そのような解析対象物体を持つ解析空間において、図９に示すような電圧波形をＰＧ７０３から導体７０１に印加した場合の印加端での電流・電圧波形を求める解析を検証解析として行った。その実測結果は非特許文献２から得た。また、比較・検討のために、導体７０１を格子上に配置する従来の手法（以降、「従来手法１」と記す）及び導体７０１を格子上に配置せずにそれを階段近似する従来の手法（以降、「従来手法２」と記す）による解析を行った。導体７０１は、従来手法１では図８（ａ）に示すように、ｙ、ｚの各軸方向とは直交し、ｘ軸方向とは平行な面の格子上に配置したのに対し、本実施の形態では図８（ｂ）に示すように、ｚ軸方向とは直交するｘｙ平面にｘ、ｙの各軸方向とは交差させる形で配置した。従来手法２でも基本的には図８（ｂ）に示すようにして階段状に配置した。
【００６９】
非特許文献２に記載されたＦＤＴＤ法による計算では、導体７０１の存在する位置に定義された導体７０１方向成分をゼロとし、その半径は考慮していない。解析空間を分割するセルの一辺の長さをΔｓとすると、ＦＤＴＤ法による計算半径は０．２２９８Δｓ［ｍ］として取り扱っている。
【００７０】
印加端解放時の電圧波形は、図９に示すように模擬し、銅板７０２で模擬した大地の抵抗率は１．６９×１０-8オーム・メートルとしている。実験で用いられたカレントセンサは立ち上がり時間が２ｎｓの製品である。それにより、単位ステップ電流が流れたとしても、その応答は次式で表される１次遅れとなる。
【００７１】
【数９】
0003703812.tif 000010
【００７２】
（１７）式を微分すると、単位インパルス応答ｈ（ｔ）は次式で表される。
【００７３】
【数１０】
0003703812.tif 000011
【００７４】
よって、観測される波形とＦＤＴＤ法を用いた計算結果を比較する場合、実測結果に次のコンボリュージョン演算を行っている。
ｉ0 ＝ｈ（ｔ）＊ｉ（ｔ）                          （１９）
その演算を行って得た実測結果と比較する形で本実施の形態による解析結果を図１０に示す。同様に、その実測結果と比較する形で従来手法１、２による解析結果を図１１、図１２にそれぞれ示す。また、本実施の形態の解析結果を従来手法１の解析結果と比較する形で図１３に示す。
【００７５】
図１０〜図１３から明らかなように、導体７０１を格子上に配置する、つまり解析空間内に存在する導体を最適な状態に配置する従来手法１は、導体７０１を階段近似する従来手法２と比較して解析結果に高い精度が得られることが分かる。本実施の形態の解析結果は、その従来手法２のそれと比較して、精度が大きく向上し、従来手法１に近い精度が得られることが分かる。その従来手法１に近い精度が得られることから、本発明を電磁界解析に適用すると、導体７０１を格子上に配置しなくとも、従来の手法では得られない高い精度を達成できることが確認できた。
＜第２の実施の形態＞
上記第１の実施の形態では、線で表現（模擬）する導体で分割された分割面のなかで磁界強度Ｈを計算する位置を含まないものは隣接する面に仮想的に加えて扱っている。そのようにして加える分割面では、隣接する面（ここでは元の面と同じ大きさである）に比べて小さいことから、理論上、それを貫く磁束はゼロと近似している。それにより、例えば図３（ａ）に示すケースでは、面ＡＢＥＦを貫く磁束はゼロとし、図３（ｂ）に示すケースでは、面ＧＥＦを貫く磁束をゼロとしている。しかし、そのように近似する場合、導体の配置によっては比較的に大きな誤差が生じてしまう可能性がある。第２の実施の形態は、その誤差をより抑えられるようにしたものである。以降、第１の実施の形態から異なる部分について詳細に説明する。
【００７６】
導体を表現する線は、図３（ａ）、（ｂ）に示すように、対向する２辺、或いは隣接する２辺を横切る。後者における説明は前者のそれと比較して容易である。このため、前者に注目して、第１の実施の形態から異なる部分を説明する。
【００７７】
第２の実施の形態では、導体を表現する線が対向する２辺を横切る場合、その線の中点を考慮する。それにより、図３（ａ）に示すようなケースでは、図１４（ａ）に示すように、導体を表現する線ＡＦの中点Ｌを考慮する。点Ｋ、Ｊは、中点Ｌとｙ方向上の位置が同じである辺ＢＨ、ＥＧ上の点である。点Ｉは、点Ｆとｙ方向上の位置が同じである辺ＢＨ上の点である。
【００７８】
線ＡＦ（以降、「導体ＡＦ」とも表記する）が導体である場合、三角形の形状をした面ＪＡＬ、ＫＦＬでは、それを貫く磁束は方向は逆で大きさは等しい。この導体ＡＦを解析空間内に表現するために、導体ＡＦを含んだ積分路ＡＢＣＤＥＦＡにファラデーの法則を用いて計算した磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）と、積分路ＡＢＣＤＥＫＪＡに，同様にファラデーの法則を用いて計算した磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）は等しくなるという条件を利用する。
【００７９】
面ＡＢＣＤＥＦ、ＡＢＣＤＥＫＪは、面ＪＡＬ、ＫＦＬの面積が等しいために、それらの面積も等しくなっている。中点Ｌをさらに考慮するのはこのためである。それにより、面ＡＢＣＤＥＦと面積、磁界強度Ｈが共に等しい別の面を利用して、その面積、つまり面ＢＣＤＥの面積に面ＡＢＥＦの面積を加えた面積での磁界強度Ｈをより高精度に算出するようにしている。
【００８０】
面ＢＣＤＥの磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）はファラデーの法則を用いて次式で計算される。
【００８１】
【数１１】
0003703812.tif 000012
【００８２】
ここで右辺の第２項目では、
−Ｋ3z×１／（面BCDEの面積）×（面BCDE上の電界強度Ｅの周回積分値）
という関係が成立している。次に、この関係を用いて、積分路ＡＢＣＤＥＦＡ、ＡＢＣＤＥＫＪＡの各各積分路についてファラデーの法則を用いて磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）を計算し、それぞれの積分路で計算された磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）は等しいという条件に従い、電界強度Ｅy、Ｅx をｘ、ｙ方向で中央差分することで（２０）式が導出できる。ただし、導体ＡＦとの位置関係から、電界強度Ｅx の線ＫＬ上の積分値と線ＬＪ上の積分値は大きさが等しく符号が逆となり、電界強度Ｅx の線ＫＬ上の積分値は零（ゼロ）となる。また、電界強度Ｅy の線ＪＡ上の積分値と線ＫＦ上の積分値は等しくなる。
【００８３】
【数１２】
0003703812.tif 000013
【００８４】
さらに，以上の（１９）（２０）式により計算された磁界強度Ｈは等しいので、辺ＢＥ上の電界強度Ｅx （ｉ＋１／２，ｊ，ｋ）は次式から計算される。
【００８５】
【数１３】
0003703812.tif 000014
【００８６】
上述したようにして導出する磁界強度Ｈでは、考慮すべき面積を全て考慮している。それにより、磁界強度Ｈの導出に、仮想的に加えた分割面の面積分の磁束が反映されることになる。このため、（１５）式から導出されるものと比較して、より正確になる。（２１）式は、そのように磁界強度Ｈをより正確に導出することを前提にして得られる式であることから、それにより導出される電界強度Ｅもより正確なものとなる。
【００８７】
ところで、（２０）式は時間領域についてのみ中央差分を適用したファラデーの法則から磁界強度Ｈz （ｉ＋１／２，ｊ−１／２，ｋ）を導出するものである。電界強度Ｅx にはｙ方向上に１次オーダーの打ち切り誤差を含む中央差分が用いられている。表記法を変えて（２０）式の内容を示したものを（２２）式に、その中央差分を（２３）式にそれぞれ示す。
【００８８】
【数１４】
0003703812.tif 000015
【００８９】
（２３）式が１次オーダーの打ち切り誤差を含むことは以下のように証明できる。
まず、次のテイラー展開を考える。
【００９０】
【数１５】
0003703812.tif 000016
【００９１】
ただし、ｕ＝１＋２ＥＫ／Δｙである。（２５）式−（２４）式より
【００９２】
【数１６】
0003703812.tif 000017
【００９３】
となり、１次オーダーの打ち切り誤差を生じることがわかる。（２６）式中の「δ」はクロネッカーのデルタであり、「ο（２）」は２次以上の項を表すシンボルである。
【００９４】
第２の実施の形態では、１次オーダーの打ち切り誤差を以下のようにしてさらに２次オーダーに減少させている。まず、次の新たなテイラー展開を考える。
【００９５】
【数１７】
0003703812.tif 000018
【００９６】
（２７）、（２８）式より、
【００９７】
【数１８】
0003703812.tif 000019
【００９８】
また、（２６）式＋（２９）式より、
【００９９】
【数１９】
0003703812.tif 000020
【０１００】
となる。以上より、
【０１０１】
【数２０】
0003703812.tif 000021
【０１０２】
が導出される。そのようにして陰的に導出される（３１）式を用いることにより、（１５）式では回避できない面積ＡＢＥＦを貫く磁束を近似的に零としたことによる誤差、さらには（２０）式に含まれる１次オーダーの打ち切り誤差が軽減される。それにより、第１の実施の形態と比較して、より正確な解析モデルを得ることができる。図１４（ｂ）に示すケースでは、磁界強度Ｈz n+1/2（ｉ＋１／２，ｊ＋１／２，ｋ）も同様の手法により計算される。
【０１０３】
上述したような計算を図５に示すシミュレーション装置に行わせる場合、（３１）式を用いた磁界強度Ｈの計算は、図７に示す電磁界解析処理では例えばステップＳ９で行わせれば良い。（２１）式を用いた電界強度Ｅの補正は、例えばステップＳ４で行わせれば良い。このようなことから、図５に示すシミュレーション装置での第２の実施の形態の実現は、その動作の流れは基本的に第１の実施の形態におけるそれから変更させることなく行うことができる。
【０１０４】
次に、上述したようにして実現できる第２の実施の形態によるシミュレーション装置を用いて行った解析結果について説明する。その解析は、上記第１の実施の形態と同じもの、即ち図７に示すように導体が存在する解析空間解析空間において、図９に示すような電圧波形をＰＧ７０３から導体７０１に印加した場合の印加端での電流・電圧波形を求める解析を検証解析として行った。計算モデルが異なる以外、各種条件等は第１の実施の形態と全て同じである。そのようにして得られた解析結果を従来手法１による解析結果と比較する形で図１５に示す。また、従来手法１、第１の実施の形態の各モデル、更には実測結果を含め、導体７０１の始端から伝搬した電流が終端で反射され始端に戻るまでの伝搬時間（Propagation time）［ｎｓ］、電流ピーク値（Peak value of current ）［Ａ］を図１６に示す。
【０１０５】
図１５、図１６、及び図１３から明らかなように、解析結果は従来手法１のそれと非常に良く一致している。第１の実施の形態と比較して、大きく精度が向上していることが分かる。これらのことから、第２の実施の形態で採用したモデルは、精度を向上させるうえで非常に有効であることが確認できる。特には図示していないが、（３１）式の代わりに（２０）式を用いても、第１の実施の形態と比較して精度が向上することが確認された。
【０１０６】
なお、第２の実施の形態では、導体によって分割されて隣接する面に仮想的に加える分割面の面積を正確に扱うために、言い換えればより高い精度を得るために、導体の中点（図１４（ａ）では「Ｌ」の他に「Ｋ」や「Ｊ」も対応）を考慮しているが、別の点を考慮しても良い。例えば図１４（ａ）に示すケースでは、点Ｌの代わりに点Ａ、Ｆを考慮するようにしても良い。
【０１０７】
本実施の形態（第１、及び第２の実施の形態）では、導体を表現（模擬）する線として直線を用いているが、曲線であっても良い。線で表現する対象は、導体や構造物などに限定されるわけではなく、解析空間のなかに存在する媒質のなかから必要に応じて選択すれば良い。導体が横切る辺での電界強度Ｅの近似方法などについても、本実施の形態に限定されるわけではなく、様々な変形を行うことができる。
【０１０８】
【発明の効果】
以上、説明したように本発明は、解析空間に存在する媒質のなかから選択した媒質（例えば導体や構造物など）を線として設定（表現）し、セルを構成する面のなかでその線が辺を横切る面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求める。
【０１０９】
その線で面を複数に分割して電界強度、及び磁界強度の少なくとも一方を求めることにより、その線はセルの面の辺上に配置した形で扱うことになる。その線はセルの面上に配置すれば良いから、高い自由度が得られ、汎用性は向上する。高い自由度を実現させたことで線で表現する媒質の形状はより正確に模擬できるようになる。これらのことから、高い汎用性を実現させつつ、電磁界解析は確実に高精度に行えるようになる。
【０１１０】
上記線による面の分割を、その面をその線で分けて得られる複数の分割面のなかで磁界強度を計算する位置を含まない分割面は隣接する面に加え、その位置を含む分割面は隣接する面に加えない形で行うようにした場合には、その分割に伴い面の数が増大することは回避される。その結果、最終的に求める強度の数、即ち解析結果として保存の対象となるデータの量の増大は回避でき、必要な計算機資源の増大は抑えることができる。
【０１１１】
分割面を隣接する面に加えることにより面積に仮想的な変化が生じる面において、その変化を電界強度、或いは磁界強度の計算に反映させるようにした場合には、その変化に伴って生じる誤差は回避、或いは低減させることができる。このため、その変化を計算に反映させる強度はより高精度に得ることができる。強度の変化は、隣接する面間で伝搬していく。それにより、近傍の面で求められる強度には、媒質が同じであれば比較的に大きな違いがないのが普通である。このことから、その反映を面積の変化が生じる面の近傍に位置する面で求められる電界強度、或いは磁界強度を用いて行うようにしたときには、その変化に伴って生じる誤差の回避、或いは低減はより確実に行うことができる。その結果、計算する強度は更に高精度に得ることができるようになる。
【図面の簡単な説明】
【図１】解析空間を分割するセルを説明する図である。
【図２】電磁界の時間配置を説明する図である。
【図３】セルを構成する面とその面の辺を横切る導体の位置関係を説明する図である。
【図４】本実施の形態によるシミュレーション装置の回路構成図である。
【図５】解析プログラムの機能構成を示す図である。
【図６】電磁界解析処理のフローチャートである。
【図７】導体の位置関係を示す図である。
【図８】導体の配置を説明する図である。
【図９】パルスジェネレータの電圧波形の時間変化を示す図である。
【図１０】第１の実施の形態による解析結果を示す図である。
【図１１】従来手法１による解析結果を示す図である。
【図１２】従来手法２による解析結果を示す図である。
【図１３】第１の実施の形態による解析結果を従来手法１による解析結果と比較して示す図である。
【図１４】第２の実施の形態における磁界強度、及び電界強度の計算方法を説明するための図である。
【図１５】第２の実施の形態による解析結果を示す図である。
【図１６】各モデルでの伝搬時間と電流ピーク値を示す図である。
【符号の説明】
３０１    ＣＰＵ
３０２    ＲＯＭ
３０３    ＲＡＭ
３０４    通信インターフェース
３０５    記憶装置
３０６    記憶媒体読み取り装置
３０７    入力装置
３０８    出力装置
３２０、３３０    記録媒体
５００    解析プログラム
５０１    入力部
５０２    解析部
５０３    出力部
５０４    境界条件処理部
５０５    電界解析部
５０６    磁界解析部
５０７    電界強度計算部
５０８    電界強度補正部
５０９    磁界強度計算部
５１０    磁界強度補正部
７０１    導体
７０２    銅板
７０３    パルスジェネレータ
【図１】
0003703812.tif 000022
【図２】
0003703812.tif 000023
【図３】
0003703812.tif 000024
【図４】
0003703812.tif 000025
【図５】
0003703812.tif 000026
【図６】
0003703812.tif 000027
【図７】
0003703812.tif 000028
【図８】
0003703812.tif 000029
【図９】
0003703812.tif 000030
【図１０】
0003703812.tif 000031
【図１１】
0003703812.tif 000032
【図１２】
0003703812.tif 000033
【図１３】
0003703812.tif 000034
【図１４】
0003703812.tif 000035
【図１５】
0003703812.tif 000036
【図１６】
0003703812.tif 000037