(57)【特許請求の範囲】
【請求項１】
コンピュータを、交番磁界発生源から送信される電磁界の波長が孔を含む金属構造体の大きさより大きい場合の、前記金属構造体周辺の電磁界分布を計算する手段として機能させるためのプログラムであって、コンピュータに、前記金属構造体表面を複数の微小面積からなる金属要素に分割し、各金属要素上に流れるループ電流値を変数として設定する第１のステップ、前記交番磁界発生源から送信された前記金属要素上での磁束密度を外部磁束密度として計算する第２のステップ、前記外部磁束密度と前記ループ電流値とを用いて、各金属要素における磁場の法線方向成分が０となる方程式を、前記各金属要素に対してそれぞれ導出し、前記孔の外周を流れる電流値を変数として設定して前記孔に対しファラデーの法則を用いた方程式を導出し、導出された前記孔に対する前記方程式と、前記各金属要素に対して導出された方程式とを、前記孔の外周を流れる電流と前記各金属要素上を流れるループ電流との相互作用を考慮して連成させて、前記各金属要素と前記孔とに対する連立方程式を構築する第３のステップ、前記連立方程式を解いて、前記各金属要素上に流れるループ電流値と前記孔の外周を流れる電流値とを求める第４のステップ、及び得られた前記ループ電流値と前記孔の外周を流れる電流値とを用いて、前記金属構造体周辺の磁場分布を計算する第５のステップを実行させるための電磁界分布計算プログラム。
【請求項２】
コンピュータを、孔に要素を生成し、生成した孔要素に仮想的なループ電流値を変数として設定し、前記孔要素上での外部磁束密度と前記仮想的なループ電流値とを用いて、前記孔要素における磁場の法線方向成分が０となる方程式を導出すると共に、前記仮想的なループ電流値が前記孔の外周に流れる電流値と同じ値であるとして、前記孔に対する方程式を導出する手段として機能させるための請求項１記載の電磁界分布計算プログラム。
【請求項３】
コンピュータを、金属構造体が閉導体を含む場合、前記閉導体を構成する特定の金属要素のループ電流値を所定の値に固定して連立方程式を解く手段として機能させるための請求項１または２記載の電磁界分布計算プログラム。
【請求項４】
コンピュータに、第４のステップにおいて連立方程式を解く際に、前記連立方程式を複数のブロックに分割するステップＡ、各ブロック毎に他のブロックからの寄与を考慮して各ブロックの電流分布を求めるステップＢ、及び得られた各ブロックの電流分布が収束するまで、または得られた各ブロックの電流分布を用いて第５のステップで計算する電磁界分布が収束するまで、前記ステップＢの計算を反復するステップＣを実行させるための請求項１〜３のいずれかに記載の電磁界分布計算プログラム。
【請求項５】
コンピュータを、金属構造体表面を複数の微小面積からなる金属要素に分割すると共に、前記金属構造体表面を複数のブロックに分割し、各ブロック毎に、それぞれ第２のステップ及び第３のステップを施して各ブロックに対する連立方程式を構築し、第４のステップでは、各ブロック毎に構築した前記連立方程式を解いて各ブロックの電流分布を求める手段として機能させるための請求項１〜３のいずれかに記載の電磁界分布計算プログラム。
【請求項６】
コンピュータに、第４のステップにおいて各ブロックの電流分布を求める際に、各ブロック毎に他のブロックからの寄与を考慮して各ブロックの電流分布を求めるステップＢ、及び得られた各ブロックの電流分布が収束するまで、または得られた各ブロックの電流分布を用いて第５のステップで計算する電磁界分布が収束するまで、前記ステップＢの計算を反復するステップＣを実行させるための請求項５記載の電磁界分布計算プログラム。
【請求項７】
コンピュータを、各ブロックの電流分布の初期値を求める際に、ブロック間の相互作用を無視して上記初期値を求める手段として機能させるための請求項４または６記載の電磁界分布計算プログラム。
【請求項８】
コンピュータを、各ブロックの電流分布の初期値を求める際に、外部磁場の発信源に最も近いブロックの電流分布は他のブロックからの寄与を無視し、上記ブロック以外の電流分布は他のブロックからの寄与を考慮して解く手段として機能させるための請求項４または６記載の電磁界分布計算プログラム。
【請求項９】
コンピュータを、各金属要素における外部磁束密度の法線方向成分の大きい順に並べた連立方程式を、前記外部磁束密度の法線方向成分の大きさに応じて複数のブロックに分割し、各ブロックの電流分布の初期値を求める際に、前記外部磁束密度の法線方向成分の最も大きいブロックの電流分布は他のブロックからの寄与を無視し、上記ブロック以外の電流分布は他のブロックからの寄与を考慮して解く手段として機能させるための請求項４記載の電磁界分布計算プログラム。
【請求項１０】
コンピュータに、孔を含む金属構造体表面を複数の微小面積からなる金属要素に分割し、各金属要素上に流れるループ電流値を変数として設定するステップ、前記孔に要素を生成し、生成した孔要素に仮想的なループ電流値を変数として設定するステップ、前記金属構造体表面を複数のブロックに分割するステップ、各ブロックに含まれる前記金属要素及び前記孔要素に対し、それぞれ交番磁界発生源から送信された前記各要素上での磁束密度を外部磁束密度として計算するステップ、各ブロックに含まれる各金属要素に対しては、前記外部磁束密度と前記ループ電流値とを用い、各金属要素における磁場の法線方向成分が０となる方程式を導出し、各ブロックに含まれる孔要素に対しては、前記孔要素上での外部磁束密度と前記仮想的なループ電流値とを用いて、前記孔要素における磁場の法線方向成分が０となる方程式を導出すると共に、前記仮想的なループ電流値が前記孔の外周に流れる電流値と同じ値であるとして、前記孔に対する方程式を導出し、導出された前記孔に対する前記方程式と、前記各金属要素に対して導出された前記方程式とを、前記孔の外周を流れる電流と前記各金属要素上を流れるループ電流との相互作用を考慮して連成させて、前記各金属要素と前記孔とに対する連立方程式を、各ブロック毎に構築するステップ、各ブロックに対して閉導体が存在するか否かを判定し、閉導体が存在する場合は、前記閉導体を構成する特定の金属要素のループ電流値を所定の値に固定するステップ、各ブロックの連立方程式を独立に解き、各ブロックの電流分布を計算するステップ、得られた各ブロックの電流分布を用い、各ブロック間の相互作用を考慮して各ブロックの電流分布の補正を行うステップ、及び補正された電流分布を用いて、前記金属構造体周辺の電磁界分布を計算するステップを実行させるための請求項１記載の電磁界分布計算プログラム。
【請求項１１】
交番磁界発生源から送信される電磁界の波長が孔を含む金属構造体の大きさより大きい場合の、前記金属構造体周辺の電磁界分布を計算する装置であって、複数の微小面積からなる金属要素に分割された金属構造体表面の金属要素データを入力する入力手段、各金属要素上に流れるループ電流値を変数として設定するループ電流設定手段、交番磁界発生源から送信された前記金属要素上での磁束密度を外部磁束密度として計算する外部磁束密度計算手段、前記外部磁束密度と前記ループ電流値とを用いて、各金属要素における磁場の法線方向成分が０となる方程式を、前記各金属要素に対してそれぞれ導出し、前記孔の外周を流れる電流値を変数として設定して前記孔に対しファラデーの法則を用いた方程式を導出し、導出された前記孔に対する前記方程式と、前記各金属要素に対して導出された方程式とを、前記孔の外周を流れる電流と前記各金属要素上を流れるループ電流との相互作用を考慮して連成させて、前記各金属要素と前記孔とに対する連立方程式を構築する構築手段、前記連立方程式を解いて、前記各金属要素上に流れるループ電流値と前記孔の外周を流れる電流値とを求める求解手段、及び得られた前記ループ電流値と前記孔の外周を流れる電流値とを用いて、前記金属構造体周辺の電磁界分布を計算する計算手段を備えたことを特徴とする電磁界分布計算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、交番磁界発生源から送信された電磁波が金属構造体で散乱される場合の電磁界分布を計算する計算方法、及びこの計算方法を用いた装置に関するものである。また、電磁界分布計算プログラムに関するものである。
【背景技術】
【０００２】
たとえば、自動車内に設置された送信アンテナから送信された電磁波が車体構造で散乱または反射される場合の電磁界分布、あるいは送信アンテナから送信された電磁波が車体構造にループ電流を流した結果生成される電磁界分布を計算する場合、従来の電磁界分布計算方法では、たとえば、下記の式のような磁界積分方程式を用いて解析を行っていた（例えば、非特許文献１参照。）。【０００３】
【数１】
0004665969.tif 000002
【０００４】
ここで、Ｈi（ｒ）は金属構造体表面の観測座標ｒにおける入射磁界ベクトルであり、Ｋ（ｒ）は観測座標ｒにおいて金属構造体表面に生成される表面電流密度である。前記磁界積分方程式は、金属構造体に外部磁界Ｈi（ｒ）が印加された場合に、金属表面に生成される表面電流密度Ｋ（ｒ）を解く方程式である。ｎに＾記号がついた変数は金属構造体表面の法線方向の単位ベクトルである。したがって、前記方程式の左辺は入射磁界ベクトルの接線方向成分を与えている。Ωは観測座標から見込んだｒにおける金属構造体が張る立体角、ωは高周波磁界の角周波数、μは媒質の透磁率、∇はベクトル解析における微分演算子、ｄＳ’は面積の微小要素、ＧA（ｒ｜ｒ’）は座標ｒ’に位置する微小電流ベクトルが座標ｒに生成するベクトルポテンシャルを与えるグリーン関数である。なお、右辺第２項は、前記表面電流密度による金属構造体表面の観測座標ｒにおける磁界の接線成分を与えている。以上から、前記の積分方程式は磁界の接線成分に関するものである。【０００５】
【非特許文献１】Ａ．Ａ．Ｋ．Ｍｏｈｓｅｎ著、Ｍａｇｎｅｔｉｃ  Ｆｉｅｌｄ  Ｉｎｔｅｇｒａｌ  Ｅｑｕａｔｉｏｎ  ｆｏｒ  Ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ  Ｓｃａｔｔｅｒｉｎｇ  ｂｙ  Ｃｏｎｄｕｃｔｉｎｇ  Ｂｏｄｉｅｓ  ｏｆ  Ｒｅｖｏｌｕｔｉｏｎ  ｉｎ  Ｌａｙｅｒｅｄ  Ｍｅｄｉａ、  Ｐｒｏｇｒｅｓｓ  Ｉｎ  Ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃｓ  Ｒｅｓｅａｒｃｈ、  ＰＩＥＲ  ２４巻、１９９９年（第２２頁、式３）
【発明の開示】
【発明が解決しようとする課題】
【０００６】
従来の電磁界分布計算方法で用いられていた前記磁界積分方程式は、非特許文献１の第２２頁、１〜６行目にも書かれているように、閉じた金属構造体に対してのみ適用できることが知られている。例えば、１枚の金属シートは開いた構造体であるが、そのシートのエッジ位置で、前記式の右辺の第１項を評価してみると、まず立体角Ωは、４πをわずかに下回る値となることがわかる。一方、エッジ位置の電流密度は、一般に非常に増大することが知られている。この結果、第１項の（１−Ω／４π）はゼロに近い値であり、Ｋ（ｒ）は非常に大きな値になるはずである。このようなゼロに近い値と非常に大きい値の積が右辺の第１項であることは、前記積分方程式を解いて、金属シート上の電流分布を計算した場合に、少なくともエッジ位置で計算精度が大きく損なわれることが容易に推定できる。
以上の考察から、従来の電磁界分布計算方法で用いられていた磁界積分方程式では、金属シートのような開いた構造体を高精度に解析することが困難であった。【０００７】
本発明は、かかる問題点を解決するためになされたものであり、車など、表面が金属シートからなる構造体のように開いた金属構造を含む場合でも、高精度、かつ高速に電磁界分布が解析計算できる電磁界分布計算方法を得ることを目的としている。
また、前記電磁界電磁界分布計算が可能な装置及び電磁界分布計算プログラムを提供することを目的としている。
【課題を解決するための手段】
【０００８】
この発明に係る電磁界分布計算プログラムは、コンピュータを、交番磁界発生源から送信される電磁界の波長が孔を含む金属構造体の大きさより大きい場合の、前記金属構造体周辺の電磁界分布を計算する手段として機能させるためのプログラムであって、コンピュータに、前記金属構造体表面を複数の微小面積からなる金属要素に分割し、各金属要素上に流れるループ電流値を変数として設定する第１のステップ、前記交番磁界発生源から送信された前記金属要素上での磁束密度を外部磁束密度として計算する第２のステップ、前記外部磁束密度と前記ループ電流値とを用いて、各金属要素における磁場の法線方向成分が０となる方程式を、前記各金属要素に対してそれぞれ導出し、前記孔の外周を流れる電流値を変数として設定して前記孔に対しファラデーの法則を用いた方程式を導出し、導出された前記孔に対する方程式と、前記各金属要素に対して導出された方程式とを、前記孔の外周を流れる電流と前記各金属要素上を流れるループ電流との相互作用を考慮して連成させて、前記各金属要素と前記孔とに対する連立方程式を構築する第３のステップ、前記連立方程式を解いて、前記各金属要素上に流れるループ電流値と前記孔の外周を流れる電流値とを求める第４のステップ、及び得られた前記ループ電流値と前記孔の外周を流れる電流値とを用いて、前記金属構造体周辺の電磁界分布を計算する第５のステップを実行させるためのものである。
【０００９】
また、この発明に係る電磁界分布計算装置は、交番磁界発生源から送信される電磁界の波長が孔を含む金属構造体の大きさより大きい場合の、前記金属構造体周辺の電磁界分布を計算する装置であって、複数の微小面積からなる金属要素に分割された金属構造体表面の金属要素データを入力する入力手段、各金属要素上に流れるループ電流値を変数として設定するループ電流設定手段、交番磁界発生源から送信された前記金属要素上での磁束密度を外部磁束密度として計算する外部磁束密度計算手段、前記外部磁束密度と前記ループ電流値とを用いて、各金属要素における磁場の法線方向成分が０となる方程式を、前記各金属要素に対してそれぞれ導出し、前記孔の外周を流れる電流値を変数として設定して前記孔に対しファラデーの法則を用いた方程式を導出し、導出された前記孔に対する前記方程式と、前記各金属要素に対して導出された方程式とを、前記孔の外周を流れる電流と前記各金属要素上を流れるループ電流との相互作用を考慮して連成させて、前記各金属要素と前記孔とに対する連立方程式を構築する構築手段、前記連立方程式を解いて、前記各金属要素上に流れるループ電流値と前記孔の外周を流れる電流値とを求める求解手段、及び得られた前記ループ電流値と前記孔の外周を流れる電流値とを用いて、前記金属構造体周辺の電磁界分布を計算する計算手段を備えたものである。
【発明の効果】
【００１１】
この発明は、エッジのある開いた構造体に対しても高精度に電流解を計算できる積分方程式を構築しているため、高精度で、かつ高速に電磁界分布計算が可能となる効果がある。
【図面の簡単な説明】
【００１２】
【図１】本発明の実施の形態１による電磁界分布計算により得られた磁場分布を示す図である。
【図２】本発明の実施の形態１に係わるループループ電流モデルを示す図である。
【図３】本発明の実施の形態１に係わる磁場分布を示す図である。
【図４】本発明の実施の形態１による電磁界分布計算方法を示すフローチャート図である。
【図５】本発明の実施の形態１による電磁界分布計算装置を示すブロック図である。
【図６】本発明の実施の形態２に係わるループループ電流モデルを示す図である。
【図７】本発明の実施の形態４による電磁界分布計算方法を示すフローチャート図である。
【図８】本発明の実施の形態４に係わるブロック分割を説明する図である。
【図９】本発明の実施の形態５による電磁界分布計算方法を示すフローチャート図である。
【図１０】本発明の実施の形態５に係わるブロック分割を説明する図である。
【図１１】本発明の実施の形態６による電磁界分布計算方法を示すフローチャート図である。
【図１２】本発明の実施の形態７による電磁界分布計算方法を示すフローチャート図である。
【図１３】本発明の実施の形態８による電磁界分布計算方法を示すフローチャート図である。
【図１４】本発明の実施の形態９に係わるアンテナ解析モデルを示す斜視図である。
【図１５】本発明の実施の形態９に係わるアンテナ解析モデルを示す断面構成図である。
【図１６】本発明の実施の形態９による解析方法を説明する図である。
【図１７】本発明の実施の形態９による解析方法により得られたアンテナの電流分布を示す図である。
【図１８】本発明の実施の形態１０に係わるループ電流モデルを示す斜視図である。
【図１９】本発明の実施の形態１０に係わる辺電流モデルを示す斜視図である。
【図２０】本発明の実施の形態１０に係わる磁場分布を示す図である。
【図２１】本発明の実施の形態１０に係わる磁力分布を示す図である。
【図２２】本発明の実施の形態１０による電磁界分布計算方法を示すフローチャート図である。
【図２３】本発明の実施の形態１１による電磁界分布計算方法を示すフローチャート図である。
【図２４】本発明の実施の形態１１による他の電磁界分布計算方法を示すフローチャート図である。
【図２５】本発明の実施の形態１２による電磁界分布計算方法を説明する図である。
【図２６】本発明の実施の形態１２による電磁界分布計算方法を説明する図である。
【図２７】本発明の実施の形態１２による電磁界分布計算方法を説明する図である。
【図２８】本発明の実施の形態１２による電磁界分布計算方法を示すフローチャート図である。
【図２９】本発明の実施の形態１２による電磁界分布計算方法を説明する図である。
【図３０】本発明の実施の形態１２による他の電磁界分布計算方法を示すフローチャート図である。
【図３１】本発明の実施の形態１２による他の電磁界分布計算方法を説明する図である。
【図３２】本発明の実施の形態１２によるさらに他の電磁界分布計算方法を示すフローチャート図である。
【図３３】本発明の実施の形態１４による電磁界分布計算方法を示すフローチャート図である。
【図３４】本発明の実施の形態１４による他の電磁界分布計算方法を示すフローチャート図である。
【図３５】本発明の実施の形態１４による他の電磁界分布計算方法を示すフローチャート図である。
【図３６】本発明の実施の形態１４による他の電磁界分布計算方法を示すフローチャート図である。
【図３７】本発明の実施の形態１４による他の電磁界分布計算方法を示すフローチャート図である。
【図３８】本発明の実施の形態１４による他の電磁界分布計算方法を示すフローチャート図である。
【図３９】本発明の実施の形態１５による電磁界分布計算方法を示すフローチャート図である。
【図４０】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４１】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４２】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４３】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４４】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４５】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４６】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４７】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４８】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図４９】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図５０】本発明の実施の形態１５による電磁界分布計算方法を説明する図である。
【図５１】本発明の実施の形態１５による他の電磁界分布計算方法を示すフローチャート図である。
【図５２】本発明の実施の形態１５による他の電磁界分布計算方法を示すフローチャート図である。
【図５３】本発明の実施の形態１５による他の電磁界分布計算方法を示すフローチャート図である。
【図５４】本発明の実施の形態１５による他の電磁界分布計算方法を示すフローチャート図である。
【符号の説明】
【００１３】
１  自動車金属要素構成データ、２  要素データ入力手段、３  ループ電流設定手段、４  外部磁束密度計算手段、５  連立方程式の構築手段、６  連立方程式の求解手段、７  磁束密度分布の計算手段、１００  半波長アンテナ、２００  第１の面上のシールドロッド群、３００  第２の面上のシールドロッド群。
【発明を実施するための最良の形態】
【００１４】
実施の形態１．
スマートエントリシステムは従来のキーレスエントリー機能とイモビライザ機能を併せ持つシステムであり、住宅や車等に搭載されるシステムである。例えば車の場合、キー操作をすることなくドアロックを解除し、ワンタッチでエンジンをスタートすることができる。このシステムでは携帯型専用リモコンと車載送受信機との間で通信を行うが、通常、車から携帯機への送信はＬＦ帯域（３０ｋＨｚ〜３００ｋＨｚ）の磁場が用いられる。本実施の形態１では、このスマートエントリシステムを搭載した自動車に対して本発明による電磁界分布計算を行うものを示す。【００１５】
図１は本発明の実施の形態１による電磁界分布計算により得られた磁場分布を示す図であり、たとえば自動車１１のドア部に取り付けられた送信アンテナ１２と携帯型専用リモコンが通信できる通信可能範囲１３を模式的に示した図である。得られた磁場分布は、たとえばスマートエントリシステムの通信可能範囲を設計するために用いることができる。また、所定の領域内の磁場分布における最大磁場、最小磁場を設計する際にも用いることが可能である。【００１６】
以下に、本発明の実施の形態１による電磁界分布計算方法を示す。交番磁界発生源（ここでは送信アンテナ）から発生した交番磁界は、自動車の金属表面に生成されるループ電流によって散乱される。図２は、たとえば自動車内に装着された送信アンテナから送信された交番磁界１４が自動車などの金属構造体表面１５に生成するループ電流分布をモデル化したものである。送信アンテナから発生する交番磁界は１００ｋＨｚ程度の周波数をもっており、この場合、前記金属構造体に流れるループ電流の表皮厚さは５０ないし１００ミクロン程度である。前記表皮厚さは、金属表面と車内外の磁場評価位置までの距離に比べて充分小さいので、表面だけにループ電流が流れていると近似することができる。さらに、自動車の大きさは前記周波数の波長に比べて充分小さいため、車の金属体上では電流の位相変化は無視でき、定在波は生じない。したがって、金属表面上で電荷の蓄積はなく、空間への変位電流は無視できることになる。結果的に、電流の発散はなく、ループ電流を用いてループ電流分布を高精度にモデリングできる。【００１７】
図２において、各ループ電流をＩ1、Ｉ2、・・・とし、磁束密度を評価する評価点の座標を各ループ内に１つずつ設定した。図２では、簡単のために円形ループ電流を表示した。実際の自動車の場合はＣＡＤデータとして三角形要素と四角形要素が使われているので、これらの要素データを使えばよい。
なお、四角形以上の多角形要素を用いることも可能である。【００１８】
図３は、金属表面１５において交番磁界は必ず面に平行に走行するため、法線方向成分Ｂnが常にゼロ（Ｂn＝０）になるという電磁気学の境界条件を示している。【００１９】
図４は、本実施の形態１の電磁界分布計算方法を流れ図として示したものである。ステップＳ１で、Ｎ個の微小面積に分割された各金属要素上に流れる未知のループ電流値をＮ個の変数Ｉk（ｋ＝１〜Ｎ）としてそれぞれ定義する。ステップＳ２で、Ｎ個の金属要素上に各々、評価点Ｐk（座標ｒk、ｋ＝１〜Ｎ）を設定し、各評価点での外部磁束密度の法線方向成分と、前記各評価点Ｐkにおいて、全てのループ電流Ｉk（ｋ＝１〜Ｎ）が生成する磁束密度の法線方向成分の総和との和が０なる方程式を、各評価点毎に導出してＮ個の方程式からなる連立方程式を構築する。ここで、外部磁束密度とは、たとえば、車内に設置された送信アンテナで送信された磁束密度分布を言い、例えば、下記式で表されるビオサバールの法則によって計算を行うことができる。ここで、Ｒは金属要素上の位置座標、Ｉaiは送信アンテナにおける電流値、ｒaiは送信アンテナの電流Ｉaiの位置座標を示す。【００２０】
【数２】
0004665969.tif 000003
【００２１】
前記ビオサバールの法則によって求められる磁束密度Ｂiを、下記式に示すように、すべての送信アンテナについて足し合わせる。【００２２】
【数３】
0004665969.tif 000004
【００２３】
ステップＳ３で、連立方程式を解いて各金属要素のループ電流値Ｉk（ｋ＝１〜Ｎ）を求める。ステップＳ４で、求めた電流値Ｉk（ｋ＝１〜Ｎ）を用いて、各評価点における外部磁束密度と、各評価点における、全てのループ電流が生成する磁束密度の総和との和を計算し、金属構造体周辺の空間の磁場分布を計算する。【００２４】
図５は本実施の形態１の電磁界分布計算装置のブロック図を示している。自動車金属構造要素データ１が装置に入力されると、図４の各ステップに対応する処理が計算機内のプログラムで実行される。したがって、図５の各手段２〜７で表現された機能ブロックは計算機内のプログラム内に存在する。【００２５】
前記計算に用いる連立方程式について詳細に説明する。送信アンテナで生成された磁界と、周囲の金属表面で生成されたループ電流による磁界を合成した磁界が実際に空間で観測される磁界である。図２に示す金属表面上のループ電流モデルにおいて、アンテナで生成された交番磁場は車体の金属部にループ電流を流し、金属部表面では磁力線は金属面に平行になる。定式化すれば、下記の式（１）になる。ここで、Ｂ（ｒ）は金属表面の磁場評価点Ｐ（座標ｒ）における磁束密度、Ｂ０（ｒ）は磁場評価点Ｐ（座標ｒ）における送信アンテナによる外部磁束密度、Ｂ（ｒ’、ｒ）は金属表面上の座標ｒ’位置での電流ループによる磁場評価点Ｐ（座標ｒ）における磁束密度、ｎは金属表面の法線方向の単位ベクトルである。
Ｂ（ｒ）・ｎ＝Ｂ０（ｒ）・ｎ＋ΣＢ（ｒ’、ｒ）・ｎ＝０            （１）【００２６】
繰り返しになるが、式（１）は送信アンテナによる磁界と、周囲の金属表面で生成されたループ電流による磁界を合成した磁界が実際に空間で観測されることと、これらのトータル磁場は、金属表面では法線方向成分が常にゼロになることを意味する。積分方程式法による定式では、このような境界条件に注目して未知変数の数と評価点の数を同一にして、連立方程式を作る。図２に示すように、金属シート上に多数の円電流を分布させ、それぞれの円電流の値を変化させることにより、任意のループ電流分布を表現することができる。金属構造体が車両の場合、この考えは、実車の金属部の要素に拡張することができる。すなわち、要素が三角形と四角形の混合要素の場合、これらの各要素を循環するループ電流を考えれば、このループ電流値を変数に取ればよい。与えられたＬＦアンテナに電流を流して磁場を生成し、次に金属表面上の各要素の中心で、磁場の法線方向成分がゼロになるように各要素のループ電流値を解くことができ、下記の式（２）、（３）で定式化できる。ここで、ｇ（ｒ’、ｒ）は座標ｒ’の位置に対応するループに単位電流を流した場合の、座標ｒの評価点における磁束密度、Ｉ（ｒ’）は座標ｒ’に対応するループ電流値である。
Ｂ（ｒ’、ｒ）＝ｇ（ｒ’、ｒ）Ｉ（ｒ’）                  （２）
Σ［ｇ（ｒ’、ｒ）・ｎ］Ｉ（ｒ’）＝−Ｂ０（ｒ）・ｎ      （３）【００２７】
計算においては、金属シートを除去して、空間にループ電流を配置すればよいため、ビオサバールの法則で磁場を計算できる。
なお、この方法は誘導磁場解析用であり、変位電流を考慮する問題には適用できないことに留意する。同じ積分方程式法でも、電界積分方程式法（狭義のモーメント法）では、電流変数は金属面内の直交する２成分であり、さらに複素数であるため、今回のループ電流を変数としてビオサバールの法則で計算する場合の４倍の未知数が発生する。即ち、本発明の電磁界分布計算方法におけるアルゴリズムでは、変位電流を考慮する必要がないため、求める変数は虚数部分が無く、さらに磁場の法線方向を境界条件として用いているため、１方向の成分のみの変数であるので、電界積分方程式法に比べて未知数が１／４にできるため、メモリ容量が少なく、計算時間も短くなることが理解できる。【００２８】
なお、本発明の電磁界分布計算方法は、電磁波の散乱体である金属構造物の大きさが、電磁波の波長の１／１０以下の大きさで、かつ金属の導電率と電磁波の周波数より決まる表皮厚さに比べ、金属の厚さが厚い場合に適応できる。【００２９】
表皮厚さについて、さらに詳細に述べる。
表皮厚さをδ、透磁率をμ、周波数をｆ、金属の導電率をσとすると、δ＝√（１／πｆμσ）となる。例えば、１３４ｋＨｚにおいて、比透磁率μが３００で、導電率σが１０7［Ｓ／ｍ］の金属を考えると、表皮厚さδは０．０２５ｍｍである。表皮厚さとは、電磁波がｅ-1に減衰する厚さである。仮に、金属の厚さが表皮厚さの５倍だとすると、透過する電磁波の割合はｅ-5＝０．６７％となる。通常の測定器を用いると測定精度は１％程度である。したがって、１％程度の解析精度が求められる場合、表皮厚さの５倍以上、すなわち前記の例では、０．１２５ｍｍ以上の厚さであれば適用範囲である。さらに好ましくは、表皮厚さの１０倍以上で、この場合透過する電磁波の割合はｅ-10＝０．００４５％となる。【００３０】
さらに、波長λに対する金属構造物の大きさについても、以下に詳細に述べる。
電磁波は、進行方向の座標をｚとすると、ｅｘｐ（ｉ２πｚ／λ）で伝わることが一般に知られている。本発明では、この電磁波の位相変化を無視できると仮定しているため、適応範囲が限定される。例えば、１３４ｋＨｚにおいて、波長は、２２３９ｍである。仮に金属構造体の大きさが、波長の１／５０である４５ｍの場合、電磁波の位相は、ｅｘｐ（ｉ２π／５０）＝０．９９２＋０．１２５ｉとなる。本発明では、この位相を１とみなすため、誤差は０．８％程度である。金属構造体の大きさが、λ／１０である場合の誤差は２０％程度である。一般に、通常の測定器の測定精度は１％程度であるので、１％程度の解析精度が求められる場合は、金属構造体の大きさが波長の１／５０以下であれば本発明を適応できる。さらに好ましくは、波長の１／２００となる。この場合、０．０５％の誤差となる。【００３１】
実施の形態２．
図６は本発明の実施の形態２に係わるループ電流モデルを示す図である。
本実施の形態２は実施の形態１と同様、自動車内に装着された送信アンテナから送信された交番磁界が自動車などの金属構造体表面に生成するループ電流によって散乱される場合の電磁界分布を計算する方法であるが、本実施の形態における計算では、図６に示すように、磁束密度を評価する座標を各ループ内に３つずつ設定した。即ち、各金属要素にそれぞれ３点、合計３Ｎ個の評価点Ｐki（座標ｒki、ｋ＝１〜Ｎ、ｉ＝１〜３）を設けた。
さらに、本実施の形態２の電磁界分布計算方法では、実施の形態１で述べたと同様、金属表面では磁界の法線方向成分が常にゼロになるという点に着目し、各金属要素内の３個の評価点における磁場の法線方向成分の和が０となるというアルゴリズムに基いて電磁界分布を計算している。即ち、図４におけるステップ２で、各金属要素内の３個の評価点における磁場の法線方向成分の和が０なる方程式を、各金属要素毎に導出して、Ｎ個の方程式からなる連立方程式を構築する。ここで、金属要素内の各評価点における磁場の法線方向成分は、各評価点での外部磁束密度の法線方向成分と、各評価点Ｐkiにおいて、全てのループ電流Ｉk（ｋ＝１〜Ｎ）が生成する磁束密度の法線方向成分の総和との和である。ステップＳ３で、構築した連立方程式を解いて各金属要素のループ電流値Ｉk（ｋ＝１〜Ｎ）を求め、ステップＳ４で、求めた電流値Ｉk（ｋ＝１〜Ｎ）を用いて、各評価点における外部磁束密度と、各評価点において全てのループ電流が生成する磁束密度の総和との和を計算し、金属構造体周辺の空間の磁場分布を計算する。【００３２】
以下に、本実施の形態２における前記計算に用いる連立方程式について詳細に説明する。図６に示す金属表面上のループ電流モデルにおいて、アンテナで生成された交番磁場は車体の金属部にループ電流を流し、金属部表面では磁力線は金属面に平行になる。Ｎ個に分割された金属要素内の３個の評価点における磁場の法線方向成分の和が０なることを定式化すれば、下記の式（４）になる。
ΣiＢ（ｒki）・ｎ＝ΣiＢ０（ｒki）・ｎ＋ΣiΣjＢ（ｒ’j、ｒki）・ｎ
＝０                                                  （４）
ここで、Ｂ（ｒki）は磁場評価点Ｐki（座標ｒki）における磁束密度であり、Ｂ０（ｒki）は磁場評価点Ｐki（座標ｒki）における送信アンテナによる外部磁束密度、Ｂ（ｒ’j、ｒki）は金属表面上の座標ｒ’j位置での電流ループによる磁場評価点Ｐki（座標ｒki）における磁束密度、ｎは金属表面の法線方向の単位ベクトルである。【００３３】
実施の形態１と同様、座標ｒ’jの位置に対応する電流ループに単位電流を流した場合の、座標ｒkiにおける磁束密度をｇ（ｒ’、ｒki）、ｒ’j位置に対応するループ電流値をＩ（ｒ’j）とすると、
Ｂ（ｒ’j、ｒki）＝ｇ（ｒ’j、ｒki）Ｉ（ｒ’j）                  （５）
となる。したがって、本実施の形態２においては前記アルゴリズムが
ΣiΣj［ｇ（ｒ’j、ｒki）・ｎ］Ｉ（ｒ’j）＝−ΣiＢ０（ｒki）・ｎ      （６）
なる方程式により定式化できる。【００３４】
なお、本実施の形態では、各金属要素内の３個の評価点における磁場の法線方向成分の和を例に説明したが、複数点であれば何点でも可能である。
また、各金属要素内の３個の評価点における磁場の法線方向成分の和以外にも、金属要素内の３個の評価点における磁場の法線方向成分の二乗和＝０で連立方程式を構築してもよい。【００３５】
実施の形態３．
実施の形態１、２では、ステップ２において、各金属要素内に１つ又は３つの評価点を設けて、連立方程式を構築したが、本実施の形態では、各金属要素に対して面積分を行い、各金属要素全体を評価点として連立方程式を構築する。即ち、図４におけるステップ２で、各金属要素における外部磁束密度の法線方向成分を前記各金属要素上で面積分した値と、前記各金属要素において、全てのループ電流Ｉk（ｋ＝１〜Ｎ）が生成する磁束密度の法線方向成分の総和を、前記各金属要素上で面積分した値との和が０となる連立方程式を各金属要素毎に導出して、Ｎ個の方程式からなる連立方程式を構築する。ステップＳ３では、構築した連立方程式を解いて各金属要素のループ電流値Ｉk（ｋ＝１〜Ｎ）を求め、ステップＳ４で、求めた電流値Ｉk（ｋ＝１〜Ｎ）を用いて、各金属要素における外部磁束密度と、各金属要素において全てのループ電流が生成する磁束密度の総和との和を計算し、金属構造体周辺の空間の磁場分布を計算する。【００３６】
以下に、本実施の形態３における前記計算に用いる連立方程式について詳細に説明する。Ｎ個に分割された各金属要素において、磁場の法線方向成分が０になることを定式化すれば、下記の式（７）になる。∫skはｋ番目の金属要素の面積分、Ｓkはｋ番目の金属要素の面積範囲を示す。
∫SkＢ（ｒ）・ｎ ｄＳ＝∫SkＢ０（ｒ）・ｎ ｄＳ
＋∫SkΣjＢ（ｒ’j、ｒ）・ｎ ｄＳ
＝０                                              （７）
実施の形態１と同様、座標ｒ’jの位置に対応する電流ループに単位電流を流した場合の、座標ｒにおける磁束密度をｇ（ｒ’、ｒ）、ｒ’j位置に対応するループ電流値をＩ（ｒ’j）とすると、
Ｂ（ｒ’j、ｒ）＝ｇ（ｒ’j、ｒ）Ｉ（ｒ’j）                            （８）
となる。したがって、本実施の形態３においては前記アルゴリズムが
∫SkΣj［ｇ（ｒ’j、ｒ）・ｎ］Ｉ（ｒ’j） ｄＳ＝−∫SkＢ０（ｒki）・ｎ ｄＳ
（９）
なる方程式により定式化できる。【００３７】
実施の形態４．
図７は本発明の実施の形態４による電磁界分布計算方法を示すフローチャート図である。本実施の形態は、実施の形態１で示した図４のステップＳ３における連立方程式の求解方法を工夫したものである。
図７において、ステップＳ３１では、ステップＳ２で立てた連立方程式を複数のブロックに分割して、ブロック連立方程式を立てる。例えば、自動車の磁界分布の解析において、金属要素を図８に示すような３つのブロックＡ、Ｂ、Ｃに分割する。ブロックの分割の方法は、例えば３ブロックに分割する場合、ｚ座標に２つの閾値をもうけ、その閾値よりもｚ座標が大きい金属要素をブロックＡとし、また異なる閾値よりもｚ座標が小さい要素をブロックＢとする。残りの金属要素をブロックＣとする。分割数は特に３分割に限らず、６分割や１３分割など計算資源に応じて変更する。分割に用いる閾値は、ｚ座標以外にｘ座標やｙ座標、さらに複合的にｘｙｚ座標を組み合わせたものが考えられる。【００３８】
ステップＳ３２では、ブロック連立方程式における各ブロックを、式（３）の右辺項（外部磁場の寄与）が大きいブロックの順番に並び替える。ただし、ブロック内の要素の順番については特に操作を行わない。例えば自動車内に送信アンテナが取り付けられている場合には、送信アンテナからの距離が近いブロックの順に並び替える。図８の例では、ブロックＢが送信アンテナに一番近いのでブロック１とする。次に近いブロックＡをブロック２とし、ブロックＣをブロック３とする。【００３９】
前記ステップＳ３１、Ｓ３２の動作を具体的に示す。実施の形態１で示した式（３）を行列形式で書き下すと、【００４０】
【数４】
0004665969.tif 000005
【００４１】
と書き表すことができる。上式は例えば２００００元の連立方程式である。ただし、２００００はこの場合の金属要素の数Ｎである。
通常の計算機の資源では、２００００元の連立方程式を解くことは困難であるため、この連立方程式をブロックに分割する。分割の方法は、例えば下記のように行及び列をそれぞれ点線部分で３分割する。【００４２】
【数５】
0004665969.tif 000006
【００４３】
各ループ電流による磁束密度の法線方向成分をブロック分割した行列をＤと呼び、各ループ電流をブロック分割した行列をＪと呼び、各外部磁場ベクトル（外部磁束密度）の法線方向成分をブロック分割した行列をＥと呼ぶと、式（１１）は式（１２）のようなブロック行列方程式に書き換えることができる。ただし、式（１２）は３分割された行のブロックの順番が、送信アンテナからの距離が近いブロックの順に並び替えられたものであり、Ｄ11、Ｄ12、Ｄ13、Ｊ1、Ｅ1は最も送信アンテナに近いブロック１（ブロックＢ）の評価点に関する分割行列、Ｄ21、Ｄ22、Ｄ23、Ｊ2、Ｅ2は次に送信アンテナに近いブロック２（ブロックＡ）の評価点に関する分割行列、Ｄ31、Ｄ32、Ｄ33、Ｊ3、Ｅ3はブロック３（ブロックＣ）の評価点に関する分割行列である。【００４４】
【数６】
0004665969.tif 000007
【００４５】
次に、ステップＳ３３では、式（１２）のブロック行列方程式における分割行列Ｊ1、Ｊ2、Ｊ3（各ブロックのループ電流分布（分布）に対応）の初期値Ｊ1O、Ｊ2O、Ｊ3Oを次式より求める。なお、次式において、行列Ｄの逆行列をＤ-1と記載する。【００４６】
【数７】
0004665969.tif 000008
【００４７】
前記式（１３−１）は、外部磁場の発信源から近いブロックのループ電流は他のブロックからの寄与を無視して解くことを意味しており、式（１３−２）は、次に外部磁場の発信源から近いブロックのループ電流は最初のブロックのみ考慮して解くことを意味しており、式（１３−３）は、次に外部磁場の発信源から近いブロックを最初の２ブロックのみ考慮して解くことを意味している。このような初期値の求め方は次のような理由による。外部磁場が大きいということは、その金属要素に流れるループ電流も大きい。従って、各ブロックを、外部磁場の寄与の大きい順に並べた時、最初のブロックのループ電流の初期値を求める際には、他のブロックからの寄与を無視する。次のブロックのループ電流の初期値を求める際には、外部磁場と最初のブロックに初期値として与えたループ電流を与えそこから発生する磁場をあたかも外部磁場として考慮し、ループ電流を計算する。次のブロックのループ電流を求める際には、外部磁場と、最初の２つのブロックから発生する初期値として与えたループ電流を与えそこから発生する磁場をあたかも外部磁場として考慮し、ループ電流を計算する。このように逐次的に初期値を求めることによって良い初期値を求めることができる。良い初期値が与えられている場合、一般に次の工程（ステップＳ３４）で行う反復計算において、少ない反復回数で収束する。【００４８】
ステップＳ３４における反復計算について説明する。ステップＳ３３で得られた初期値を用い、以下の式（１４）に従って解の精度の更新を行う。反復計算にあたっては、各ブロックを外部磁場の寄与の大きい順に並べ、初期値を求めたときと同様の工夫で、寄与の大きいものを優先して計算する。なお、ｋは反復回数を示す。これにより、高速に収束する。【００４９】
【数８】
0004665969.tif 000009
【００５０】
従来のブロック反復法では、２００００元の連立方程式に対して、数十回の反復を要したが、反復の際の初期値を求め方、及び反復の順番を工夫することによってより少ない回数の反復で高精度に収束させることができる。
また、ブロック反復法は、従来よく知られている方法であるが、電磁気の散乱問題は一般には直接法で解かれており、従来、電磁気の散乱問題に反復法を利用することはなされていなかった。
本発明では上述のように工夫したブロック反復法を用い、図８の車モデルを解析したところ、６回の反復計算で３％以下の誤差に収束した。従来の直接解法で解いた場合に比べ１０倍高速で、使用メモリが１／１０で解析可能である。【００５１】
ステップＳ４では、ステップＳ３４で得られたループ電流値Ｉk（ｋ＝１〜Ｎ）を用いて、各評価点における外部磁束密度と、各評価点における、全てのループ電流が生成する磁束密度の総和との和を計算し、金属構造体周辺の空間の磁場分布を計算する。
ステップＳ５では、前回の反復計算で得られた磁場分布と比較して収束判定を行い、収束していない場合は、ステップＳ３４に戻り、さらに反復計算を繰り返す。収束している場合は、ステップＳ６で、収束した磁場分布を解とする。【００５２】
なお、図７のフローチャートでは、収束の判定に磁場分布を用いたが、ループ電流分布を収束の判定に用いることもできる。反復計算において、ある反復におけるループ電流の値と次の反復におけるループ電流の値が十分近ければ収束したとみなすことが出来る。【００５３】
実施の形態５．
図９は本発明の実施の形態５による電磁界分布計算方法を示すフローチャート図である。本実施の形態は、実施の形態４と同様、図４のステップＳ３における連立方程式の求解方法を工夫したものであるが、実施の形態４との違いは金属要素を複数のブロックに分割する方法が異なっている。本実施の形態では、金属要素を図１０に示すように外部磁束密度の大きさに応じて複数のブロックＡ１、Ａ２、Ａ３に分割する。例えば、３分割の場合、全ての金属要素を外部磁場の寄与が大きい順番に並べたあとに、上から順にＡ１（ブロック１）、Ａ２（ブロック２）、Ａ３（ブロック３）の３ブロックに分割する。
図９において、ステップＳ３０では、各金属要素を外部磁束密度の寄与が大きい金属要素の順番に並べ替える。
ステップ２では、並び替えられた順に、金属要素の評価点での外部磁束密度の法線方向成分と、前記評価点において、全てのループ電流が生成する磁束密度の法線方向成分の総和との和が０なる方程式を導出して、Ｎ個の評価点に対する連立方程式を構築する。
ステップＳ３１ではステップＳ２で立てた連立方程式を複数のブロックに分割して、ブロック連立方程式を立てる。
ステップＳ３３〜Ｓ３４では実施の形態４と同様の操作を行う。即ち、各ブロックのループ電流の初期値を求める際に、最初のブロック（ここでは外部磁束密度の最も大きいブロック１（Ａ１））のループ電流は他のブロック２（Ａ２）、およびブロック３（Ａ３）からの寄与を無視して解き、得られたループ電流値を用いて各ブロックのループ電流を反復計算する。【００５４】
本実施の形態では、外部磁場の寄与が大きい順番に金属要素を並べたあとに上から順に３ブロックに分割したため、図１０のブロック分割の結果と図８のブロック分割の結果とは異なっている。実施の形態４におけるブロック分割（図８）は、人間がまとまった構造として認識できるブロック分割であるに対して、本実施の形態５のブロック分割（図１０）は計算機による計算を伴わないとブロック分割を行うことができないが、実施の形態５のブロック分割は、例えば球体など人間がまとまった構造として認識できない場合にも機械的にブロック分割を行うことができる。【００５５】
なお、本実施の形態では、外部磁場の寄与が大きい順番に金属要素を並べたあとに、並び替えられた順に各金属要素の評価点での方程式を導出して連立方程式を構築し、上から順に３ブロックに分割したが、各金属要素に対するＮ個の連立方程式を構築し、外部磁束密度の法線方向成分の大きい順に前記Ｎ個の連立方程式を並び替え、並び替えた連立方程式を、前記外部磁束密度の法線方向成分の大きさに応じて複数のブロックに分割しても同様の結果が得られる。【００５６】
前記実施の形態では３分割で行ったが、分割数は特に３分割に限らず、６分割や１３分割など計算資源に応じて変更する。例えば、３００００要素の金属要素からなる車両モデルを用い、６分割すると、１ブロックあたり５０００要素になる。１ブロックあたりの使用メモリは、５０００×５０００×８＝２００Ｍｂｙｔｅとなる。一方、従来の方法では、３００００×３００００×８＝７．２ＧＢ必要になる。現在普及しているパソコンの多くは３２ｂｉｔであるため、使用最大メモリは約２〜４ＧＢであり、メモリも高価なので、２００ＭＢｙｔｅ程度のメモリ空間で計算できることは極めて有用な効果である。【００５７】
なお、前記実施の形態４、５においては実施の形態１で示した連立方程式の解法について説明したが、実施の形態２、３で述べた連立方程式を解く際にも、実施の形態４、５で示した分割ブロックの考え方を適用して、ループ電流値を求め、磁場分布を計算することが可能である。【００５８】
実施の形態６．
図１１は本発明の実施の形態６による電磁界分布計算方法を示すフローチャート図である。本実施の形態は、実施の形態４、５と同様、連立方程式（または行列方程式）を作成後、連立方程式を複数のブロックに分割して各ブロックごとのループ電流分布の初期値を計算し、さらに各ブロック間の相互作用を考慮して反復計算によって前記ループ電流分布を補正することによって、正しい電磁界分布に近づけていくものであるが、実施の形態４、５との違いは、各ブロックごとのループ電流分布の初期値の求め方が異なっている。本実施の形態では、前記初期値を求める際に、式（１３−１）〜（１３−３）の代わりに、下記の式（１３０）を用いる。【００５９】
【数９】
0004665969.tif 000010
【００６０】
図において、ステップＳ１００では、金属構造体をＮ個の金属要素に分割し、Ｎ個の微小面積に分割された各金属要素上に流れる未知のループ電流値をＮ個の変数Ｉk（ｋ＝１〜Ｎ）としてそれぞれ定義する。ステップＳ２００では、Ｎ個の金属要素上に各々設定された評価点Ｐk（座標ｒk、ｋ＝１〜Ｎ）での外部磁束密度の法線方向成分と、前記各評価点Ｐkにおいて、全てのループ電流Ｉk（ｋ＝１〜Ｎ）が生成する磁束密度の法線方向成分の総和との和が０なる方程式を、各評価点毎に導出してＮ個の方程式からなる連立方程式（行列方程式）を構築する。ステップＳ３１０では、ステップＳ２００で立てた行列方程式を複数のブロックに分割し、ブロック行列方程式を立てる。ステップＳ３３０では、前記ブロック行列方程式における各ブロックのループ電流分布の初期値を、前記式（１３０）を用い、ブロック間の相互作用を無視して求める。ステップＳ３４０では、ステップＳ３３０で得られた初期値を用い、ループ電流分布の補正を行う。なお、前記ブロック行列方程式は、実施の形態４（図７のステップＳ３２）、または実施の形態５（図９のステップＳ３０）で説明したように、外部磁界の大きさを考慮して配列されたものであり、補正にあたっては、実施の形態４、５と同様、式（１４）を用い、各ブロック間の相互作用を考慮して各ブロックのループ電流分布を反復計算する。ステップＳ４００では、得られたループ電流分布を用いて、金属構造体周辺の空間の磁場分布を計算する。ステップＳ５００では、得られた磁場分布と前回の反復計算で得られた磁場分布とを比較して収束判定を行い、収束していない場合は、ステップＳ３４０に戻り、さらに反復計算を繰り返す。収束している場合は、ステップＳ６００で、収束した磁場分布を出力する。【００６１】
本実施の形態では、ブロック間の相互作用を無視してループ電流分布の初期値を求めており、初期値の精度は実施の形態４、５に比べ悪くなるが、実際にプログラミングを行うさいに、ソースコードが単純になるという効果がある。【００６２】
なお、図１１のフローチャートでは、収束の判定に磁場分布を用いたが、ループ電流分布を収束の判定に用いることもできる。反復計算において、ある反復におけるループ電流の値と次の反復におけるループ電流の値が十分近ければ収束したとみなすことが出来る。【００６３】
実施の形態７．
図１２は本発明の実施の形態７による電磁界分布計算方法を示すフローチャート図である。本実施の形態は、実施の形態４、５、６と同様、行列方程式を作成後、行列方程式を複数のブロックに分割して、各ブロックごとのループ電流分布の初期値を計算し、さらに、各ブロック間の相互作用を考慮して反復計算によって補正することによって、正しい電磁界分布に近づけていくものであるが、実施の形態４、５、６との違いは、各ブロック間の相互作用を考慮した補正の方法が異なっている。
本実施の形態では、補正に際して、式（１４）の代わりに、下記の式（１４０）を用いる。なお、αは緩和係数（通常１≦α＜２）である。【００６４】
【数１０】
0004665969.tif 000011
【００６５】
図において、ステップＳ１００、Ｓ２００、Ｓ３１０、Ｓ４００、Ｓ５００、Ｓ６００は、実施の形態４〜６と同様の操作を行う。ステップＳ３３１では、式（１３−１）〜（１３−３）、または式（１３０）を用いて各ブロックごとのループ電流分布の初期値を計算する。ステップＳ３４１では、式（１４０）を用い、各ブロック間の相互作用を考慮して各ブロックごとのループ電流分布を反復計算する。【００６６】
本実施の形態では、補正に際して、緩和係数α（通常１≦α＜２）を導入することにより、ブロック間の相互作用を高速に収束させることができ、計算時間の短縮が可能となる効果がある。【００６７】
なお、図１２のフローチャートでは、収束の判定に磁場分布を用いたが、ループ電流分布を収束の判定に用いることもできる。反復計算において、ある反復におけるループ電流の値と次の反復におけるループ電流の値が十分近ければ収束したとみなすことが出来る。【００６８】
実施の形態８．
図１３は本発明の実施の形態８による電磁界分布計算方法を示すフローチャート図である。本実施の形態は、実施の形態４と同様、金属構造体の構造に基いて分割されたブロック行列方程式を用いて、電磁界分布を計算するものであるが、実施の形態４との違いは、先に金属要素を構造に基いてブロック分割してから、各ブロックの行列方程式を作成し、作成した各ブロックの行列方程式をハードディスクに保存する。その後、それぞれのブロック毎にループ電流分布の初期値を計算し、さらに、各ブロック間の相互作用を考慮して反復計算によって補正することによって、正しい電磁界分布に近づけていく。【００６９】
図において、ステップＳ１００、Ｓ４００、Ｓ５００、Ｓ６００は、実施の形態４〜７と同様の操作を行う。ステップＳ３１１では金属要素を複数のブロックＡ、Ｂ、Ｃに分割する。分割の際には、例えば、車において、ドア、ボンネット、天井、床等に分割し、それぞれの金属要素が５０００〜７０００要素以下となるようにする。ステップＳ２０１では、各ブロックを送信アンテナからの距離が近いブロックの順に並び替え、３分割の場合、ブロック１〜３とする。さらに、行列方程式（１２）における各ブロック１〜３の行列Ｄ11、Ｄ22、Ｄ33、Ｅ1、Ｅ2、Ｅ3、Ｊ1、Ｊ2、Ｊ3を作成し、各ブロックの行列方程式を構築する。ステップＳ３３０では、作成した行列Ｄ11、Ｄ22、Ｄ33、Ｅ1、Ｅ2、Ｅ3に対し、式（１３０）を用いて各ブロックごとのループ電流分布の初期値Ｊ1O、Ｊ2O、Ｊ3Oを計算する。ステップＳ２０５では、各ブロック間の相互作用の行列を作成する。すなわち、３分割の場合、行列方程式（１２）における行列Ｄ12、Ｄ13、Ｄ21、Ｄ23、Ｄ31、Ｄ32を作成する。ステップＳ３４０では、前回得られたループ電流値と、ステップＳ２０１、Ｓ２０５で作成した行列Ｄ11、Ｄ22、Ｄ33、Ｅ1、Ｅ2、Ｅ3、Ｄ12、Ｄ13、Ｄ21、Ｄ23、Ｄ31、Ｄ32とを用い、式（１４）により、各ブロックのループ電流分布値を反復計算する。【００７０】
解くべき問題の規模が大規模（具体例：金属要素が２００００要素）であるとき、行列方程式の行列（２００００×２００００）は、３．２ＧＢｙｔｅのメモリ容量を必要とし、一般に普及している３２ｂｉｔの計算機ではメモリに収まらない。６４ｂｉｔの計算機でも３ＧＢｙｔｅ以上のメモリを実装したマシンは高価で一般には普及していない。通常、メモリに収まらない場合には、ハードディスクに行列を記憶するが、先にブロック分割を行っておけば、ブロックごとに行列を作成してハードディスクに記憶することが可能となり、実際にプログラミングを行う際のソースコードが単純になるという効果がある。【００７１】
なお、図１３のフローチャートでは、反復計算に際し、式（１４）を用いたが、式（１４０）を用いても良い。この場合、行列方程式（１２）における行列Ｄ12、Ｄ13、Ｄ21、Ｄ23、Ｄ31、Ｄ32は使用しないのでステップＳ２０５は必要ない。
また、図１３のフローチャートでは、初期値の計算に際し、式（１３０）を用いたが、式（１３−１）〜式（１３−３）を用いても良い。この場合、行列方程式（１２）における行列Ｄ21、Ｄ31、Ｄ32を使用するので、ステップＳ２０５を初期値の計算の前に施す必要がある。【００７２】
また、図１３のフローチャートでは、収束の判定に磁場分布を用いたが、ループ電流分布を収束の判定に用いることもできる。反復計算において、ある反復におけるループ電流の値と次の反復におけるループ電流の値が十分近ければ収束したとみなすことが出来る。【００７３】
実施の形態９．
実施の形態４〜８で示したブロック分割の方法、及び反復法による計算は、アンテナの近傍に波長以上の大きさを有する散乱体が存在する場合の電磁界解析に適用することも可能であり、メモリの節約と高速化に効果がある。
図１４に半波長ダイポールアンテナとＬ型シールド体で構成された解析モデルを示す。図１５は図１４の断面構成図である。Ｌ型シールド体は長さ２λのロッドを１５本用いて、λ／１０間隔でＬ型配置した構成である。ダイポールアンテナとシールドロッドの距離を３λ／１０とし、ロッド直径ｄをλ／１００πとする。波長ダイポールアンテナ上の電流分布及び各導体に流れる電流分布は、電界積分方程式法で求めることができる。【００７４】
以下では、長さの単位を波長で規格化する。半波長送信アンテナを９分割し、中央要素で給電条件を満足させる。１５本の２λ長のシールドロッドもそれぞれ３６分割すれば、分割されたすべての要素がλ／１８で同一になる。座標ｚ’における電流をＩ（ｚ’）とすれば、座標ｚにおける電界の接線成分Ｅｚ（ｚ）は次式で与えられる。【００７５】
【数１１】
0004665969.tif 000012
【００７６】
ここで、εは誘電率、βは２π／λ（λは波長）、Ψ（ｒ、ｒ’）は次式で与えられる。
Ψ（ｒ、ｒ’）＝ｅｘｐ（−Ｉβ｜ｒ−ｒ’｜）／（４π｜ｒ−ｒ’｜）      （１６）【００７７】
給電点以外の金属表面では、電界の接線成分はゼロである。式（１５）を用いることにより、要素数Ｎと同じ数の未知電流Ｉ（ｚ）を含む式をＮ個用意することができる。
実施の形態４で示したブロック分割による計算方法を適用する。本実施の形態では、半波長アンテナ１００、第１の面上のシールドロッド群２００、第２の面上のシールドロッド群３００からなる３つのグループを考えると、式（１７）が成立する。【００７８】
【数１２】
0004665969.tif 000013
【００７９】
式（１７）において、Ｉｌｉｓｔ１は半波長アンテナ１００上の９分割した要素の電流リスト、Ｉｌｉｓｔ２はシールドロッド群２００上の要素の電流リスト、Ｉｌｉｓｔ３はシールドロッド群３００上の要素の電流リストである。ｅＶｅｃは半波長アンテナ１００上の電界の接線成分のリストで、給電点以外ではゼロである。式（１７）のブロック行列ｆ１１、ｆ１２、ｆ１３、ｆ２１、ｆ２２、ｆ２３、ｆ３１、ｆ３２、ｆ３３は式（１５）、（１６）から数値積分で求める。
ブロック分割法では、給電点を含むグループから計算を始める。まず、半波長ダイポールアンテナ上の電流分布Ｉｌｉｓｔ１の初期値を、シールドロッド群２００、３００を無視して求める。したがって、
Ｉｌｉｓｔ１＝ｆ１１-1・ｅＶｅｃ
で求める。ここで、ｆ１１-1は行列ｆ１１の逆行列である。
次に、アンテナ上の電流を前記のように計算したＩｌｉｓｔ１として、アンテナに近い方の面上のシールドロッド上の電流を求める。図１４の場合は、第１の面と第２の面のそれぞれとアンテナの距離は同一であるため、どちらを最初に計算しても構わない。以下では、第１の面上の電流Ｉｌｉｓｔ２を次式で計算している。
Ｉｌｉｓｔ２＝ｆ２２-1（０−ｆ２１・Ｉｌｉｓｔ１）
最後に、前記で求めたアンテナ上の電流Ｉｌｉｓｔ１及び第１の面上の電流Ｉｌｉｓｔ２を用いて、第２の面上の電流Ｉｌｉｓｔ３を次式で計算している。
Ｉｌｉｓｔ３＝ｆ３３-1（０−ｆ３１・Ｉｌｉｓｔ１−ｆ３２・Ｉｌｉｓｔ２）【００８０】
この計算の結果、アンテナ上の電流Ｉｌｉｓｔ１、第１の面上のシールドロッド電流Ｉｌｉｓｔ２、第２の面上のシールドロッド電流Ｉｌｉｓｔ３の近似解を求めることができる。図１６にアンテナ上の電流Ｉｌｉｓｔ１の実数部を要素番号に対してドットとしてプロットした。実線は従来の電界積分方程式法、すなわちブロック分割しないで、全体の電流分布を一度に求めた場合の電流分布の実数部を示す。ブロック分割法の最初の計算でほぼ正しい解が得られているが、誤差があることが認められる。
続いて、以下の３つの式をこの順に計算し、さらにこの３つの式の計算を複数回、反復する。なお、これらの式は、式（１７）から導かれた。【００８１】
Ｉｌｉｓｔ１＝ｆ１１-1（ｅＶｅｃ−ｆ１２・Ｉｌｉｓｔ２−ｆ１３・Ｉｌｉｓｔ３）
Ｉｌｉｓｔ２＝ｆ２２-1（０−ｆ２１・Ｉｌｉｓｔ１−ｆ２３・Ｉｌｉｓｔ３）
Ｉｌｉｓｔ３＝ｆ３３-1（０−ｆ３１・Ｉｌｉｓｔ１−ｆ３２・Ｉｌｉｓｔ２）【００８２】
３回反復後のアンテナ上の電流Ｉｌｉｓｔ１の実数部を要素番号に対してドットとしてプロットした結果を図１７に示す。全体の電流分布を一度に求めた場合（実線）とブロック分割法による解（ドット）が十分一致していることがわかる。
シールドロッド上の電流も同様に十分一致した。【００８３】
従来の電界積分方程式法は、全体を一度に行列方程式として解くため、大規模な計算となり、必要なメモリ容量が１０ＧＢを超え、計算が不可能になることがあったが、ブロック分割法を適用することにより、小規模な問題に分割して逐次計算を進めることが可能になる。その結果、計算機のメモリ容量も１〜２ＧＢで充分であり、汎用のパソコンなどで大規模な計算を実行できる。さらに、計算時間も大幅に低減する効果がある。【００８４】
実施の形態１０．
前記各実施の形態では金属構造体が開いた金属構造体（開導体）である場合を示したが、本実施の形態では閉導体を含む場合について説明する。
金属構造体が閉導体の場合、その閉導体を構成する金属要素のある１つのループ電流を所定の値に固定することで計算が不安定にならずに解くことができる。以下に閉導体における解析について詳細に説明する。【００８５】
本実施の形態では、例えば、図１８に示すような単純な四面体を考える。このような閉じた構造体に対し、四面体の各面を金属要素とし、実施の形態４と同様に行列方程式を作成すると次式のようになる。【００８６】
【数１３】
0004665969.tif 000014
【００８７】
ただし、Ｊ1〜Ｊ4は四面体の各面におけるループ電流、Ｅ1〜Ｅ4は外部磁場が各面に作る磁場の法線方向成分である。また、Ｄ11〜Ｄ44は四面体の各面における、各ループ電流による磁束密度の法線方向成分である。
図１８のように閉じた構造体の場合は、磁場ガウスの法則により、外部磁場の面積分が０、即ち
∫Ｅ・ｄＳ＝０
なる関係が成立するので、式（１８）の連立方程式は独立でなくなり、（Ｎ−１）次独立となる。即ち、四面体の場合、連立方程式の数は４本であるが、独立な方程式の数は３本となる。従来の方法でこの問題を解くと非常に不安定な数値解を与える。【００８８】
ところが、図１９に示すように、辺電流Ｌ1〜Ｌ6を定義すると、四面体の辺電流Ｌ1〜Ｌ6とループ電流Ｊ1〜Ｊ4は式（１９）によって関係つけられる。
Ｌ1＝Ｊ1−Ｊ2
Ｌ2＝−Ｊ1−Ｊ3
Ｌ3＝−Ｊ2＋Ｊ4
Ｌ4＝Ｊ3＋Ｊ4                （１９）
Ｌ5＝Ｊ2＋Ｊ3
Ｌ6＝Ｊ1−Ｊ4【００８９】
ここで、αをある定数とし、
Ｊ1’＝Ｊ1＋α
Ｊ2’＝Ｊ2＋α
Ｊ3’＝Ｊ3−α                （２０）
Ｊ4’＝Ｊ4＋α
とすると、式（１９）は
Ｌ1＝Ｊ1’−Ｊ2’
Ｌ2＝−Ｊ1’−Ｊ3’
Ｌ3＝−Ｊ2’＋Ｊ4’
Ｌ4＝Ｊ3’＋Ｊ4’              （２１）
Ｌ5＝Ｊ2’＋Ｊ3’
Ｌ6＝Ｊ1’−Ｊ4’
と変換される。【００９０】
このような変換を行っても辺電流としてはまったく同一のものを得る。すなわち、ループ電流は式（２０）のようなバイアスを加える変換に対して同一の辺電流を与える。これは、回路の電圧がバイアの任意性があることに類似している。したがって、閉じた導体を含む場合、あるひとつのループ電流を固定することによって解くことができることがわかる。具体的に４面体の場合について考えると、例えば、Ｊ4を０に固定した場合、式（２２）と変形できる。【００９１】
【数１４】
0004665969.tif 000015
【００９２】
式（２２）は４本の独立した連立方程式であり、これを解くことによって閉じた導体における各要素のループ電流値を求めることができ、その結果、磁場分布を解析することが可能となる。【００９３】
なお、本実施の形態においては、金属構造体が閉じた導体の場合について説明したが、金属構造体が閉じた導体を含む場合（具体例：車のガソリンタンク）についても同様の議論は適用可能である。閉じた導体のうち１つの金属要素のループ電流値をある値に固定（例えば０に固定）することによって行列方程式が不定とならず正確な電磁場分布が計算できる。【００９４】
図２２は本発明の実施の形態１０による電磁界分布計算方法を示すフローチャート図である。図において、ステップＳ１００では、金属構造体を要素に分割し、分割された各金属要素上に流れる未知のループ電流値をそれぞれ変数として定義する。ステップＳ２００では、金属要素上に各々設定された評価点での外部磁束密度の法線方向成分と、前記各評価点において、全てのループ電流が生成する磁束密度の法線方向成分の総和との和が０なる方程式を、各評価点毎に導出して連立方程式（行列方程式）を構築する。ステップＳ７００では、金属構造体が閉じた導体を含んでいるか否かを判定する（この例では判定をコンピュータを用いて自動的に行った。共通辺を持たない辺を持つ三角形を除去していくことを、形状が変化しなくなるまで繰り返す。最終的に多面体構造が無くなれば閉じていない。多面体構造が残れば、その多面体構造が閉じているという判定となる。）。閉じた導体を含まない場合は、ステップＳ３００にて、実施の形態１と同様にしてループ電流分布の計算を行う。閉じた導体を含む場合は、ステップＳ７１０で、閉導体内の１要素を一定値に固定してから、ステップＳ３００にてループ電流分布の計算を行う。ステップＳ４００では、得られたループ電流値を用いて金属構造体周辺の空間の磁場分布を計算する。ステップＳ６００では、得られた磁場分布を出力する。【００９５】
閉じた導体に磁場をかけた場合の磁場分布を、前記方法によって解析した結果を示す。
非常に遠方にあるヘルムホルツコイルを磁場発生源とし、その中心に金属球があるモデルについて考える。非常に遠方に磁場発生源があるため、金属球付近では外部磁場は、ｙ軸方向に成分をもち、かつ一様となっている。このような単純なモデルの場合、本発明のような数値計算を行わなくとも電磁気学の定理によって磁場分布を求めることができる。
金属球の半径をＲ、金属中心からみた対頂角をθ、Ｂｘ及びＢｙを磁場のＸ成分及びＹ成分とし、電磁気学の定理より磁場分布を求めると、
Ｂｘ＝２ｓｉｎθｃｏｓθ−ｓｉｎθｃｏｓθ／２／Ｒ3
Ｂｙ＝ｃｏｓθ2−ｓｉｎθ2−ｃｏｓθ2／Ｒ3−ｓｉｎθ／２／Ｒ3
となる。この数式から得られる磁力線図を図２０に示す。
同様のモデルを本実施の形態の手法で解析した場合の磁力線図を図２１に示す。
図２０の磁力線図と図２１の磁力線図とが同様の磁力線を与えることから、本実施の形態による手法が正しい結果をもたらすことがわかる。【００９６】
実施の形態１１．
図２３は本発明の実施の形態１１による電磁界分布計算方法を示すフローチャート図である。本実施の形態は、実施の形態１０と同様、金属構造体が閉導体を含む場合も考慮した行列方程式の計算例であるが、本実施の形態では、閉導体も含めた行列方程式の計算に際し、行列方程式を構造に基いてブロック分割し、ブロック毎のループ電流分布の初期値を求めると共に、各ブロックの相互作用を考慮して前記ループ電流分布を補正することにより、ループ電流分布の計算を行うものである。
図において、ステップＳ１００では、金属構造体をＮ個の金属要素に分割し、ステップＳ２００では、各金属要素に対するＮ個の方程式からなる連立方程式（行列方程式）を構築する。次に、ステップＳ７００では、金属構造体が閉じた導体を含んでいるか否かを判定する。閉じた導体を含まない場合は、ステップＳ３１０にて、ステップＳ２００でたてた行列方程式を複数のブロックに分割し、ブロック行列方程式を立てる。閉じた導体を含む場合は、ステップＳ７１０で、閉導体内の１要素を一定値に固定してから、ステップＳ３１０にてブロック行列方程式をたてる。ステップＳ３３１では、前記ブロック行列方程式における各ブロックのループ電流分布の初期値を、前記式（１３）または式（１３０）を用いて求める。ステップＳ３４２では、ステップＳ３３１で得られた初期値を用い、前記式（１４）、式（１４０）により、各ブロック間の相互作用を考慮したループ電流分布の補正を行う。ステップＳ４００では、得られたループ電流分布を用いて、金属構造体周辺の空間の磁場分布を計算する。ステップＳ５００では、得られた磁場分布と前回の反復計算で得られた磁場分布とを比較して収束判定を行い、収束していない場合は、ステップＳ３４２に戻り、さらに反復計算を繰り返す。収束している場合は、ステップＳ６００で、収束した磁場分布を出力する。【００９７】
このようにすることにより、閉導体を含むものにおいても、計算時間の短縮が可能となる効果がある。【００９８】
なお、本実施の形態では、電磁界分布を計算する際に、行列方程式を作成→閉導体検出→ブロック分割→ループ電流分布の計算の順で計算を行ったが、「行列方程式を作成」、「閉導体検出」、及び「ブロック分割」の順番を入れ替えても良い。
図２４は前記順番を入れ替えた変形例であり、閉導体検出→行列方程式を作成→ブロック分割→ループ電流分布の計算の順で計算を行ったものである。【００９９】
なお、図２３〜図２４に示す例以外にも、「行列方程式を作成」、「閉導体検出」、及び「ブロック分割」の順番を任意に入れ替えて計算しても良く、同様の効果がある。
順番を入れ替えた変形例（図を省略）を以下に示す。
・行列方程式を作成→ブロック分割→閉導体検出→ループ電流分布の計算
・閉導体検出→ブロック分割→行列方程式を作成→ループ電流分布の計算
・ブロック分割→行列方程式を作成→閉導体検出→ループ電流分布の計算
・ブロック分割→閉導体検出→行列方程式を作成→ループ電流分布の計算【０１００】
実施の形態１２．
本実施の形態では、金属構造体が孔を有する場合について説明する。金属構造体に孔が存在する場合（具体例：自動車の窓枠、フロントウィンドウ）、内側環状電流を模擬できないため、正しく解析できない問題がある。
以下に、孔の開いた円板導体からの電磁波の散乱問題について説明する。例えば、図２５（ａ）に示すように、円板導体（以下、円板リングと記す）を、実施の形態１で述べたように複数の金属要素に分割する。円板リングでは、本来、ループ電流は淵部に集中し、図２５（ｂ）に示すように、ループ電流の大部分は内側の淵と外側の淵に同じ向きに流れる。ところが、実施の形態１で述べたように、分割した各金属要素上にループ電流が流れていると仮定するとこの状態を模擬できない。すなわち、図２５（ｃ）に示すように、内側の淵と外側の淵には逆向きに電流が流れてしまう。
実際に、図２６（ａ）に示すように、ソースコイル１６（ｒ1＝半径０．１ｍ）を円板リング１７（内径２ｒ2＝０．１ｍ、外径２ｒ3＝０．３ｍ）に近接配置（間隔Ｌ＝０．１ｍ）して、円板リング１７の散乱問題を実施の形態１で述べたような方法で解析した。その結果、図２６（ｂ）の実線に示すように、ほとんど円板リングのシールド効果はないという物理的に誤った電磁界解析結果を得る。なお、図２６（ｂ）の点線は有限要素法で同様のモデルの電磁界解析を行った解析結果であり、計算時間がかかるが、実際の電磁界分布とほぼ一致する結果が得られている。点線と実線とでは異なった結果となっていることがわかる。
そこで、本実施の形態では、金属構造体に孔が存在する場合、新たな自由度として、図２７（ａ）に示すように、円板リングの内側の淵（孔の外周）に流れる電流を追加する。そして、孔の外周に流れる孔電流に対して、ファラデーの法則を適用し、外側の金属導体部分と共に連成解析を行う。
その結果、図２７（ｂ）の実線に示すように、正しい電磁界分布がえられた。有限要素法を用いて確認を行った（図２７（ｂ）の点線）ところ、両者は良く一致した結果となる。【０１０１】
金属構造体が孔を含む場合の具体的な解析法について、さらに詳細に説明する。
本実施の形態１２は、孔に対しファラデーの法則を用いた方程式を導出し、導出された前記孔に対する方程式と、各金属要素に対して導出された方程式とを連成させて、前記各金属要素と前記孔とに対する連立方程式を構築し、この連立方程式を解いて、各金属要素上に流れるループ電流値と孔の外周を流れる電流値とを計算して電流分布を求め、金属構造体周辺の電磁界分布を解析するものである。このようにすることにより、金属構造体に孔がある場合の電磁界分布が高速に計算可能になる。【０１０２】
図２８は本発明の実施の形態１２による電磁界分布計算方法を示すフローチャート図である。本実施の形態では、たとえば、図２９に示すような円板金属構造の電磁波散乱問題を考える。まず、ステップＳ１００で、円板金属構造体をＮ個の三角要素に分割する。次に、ステップＳ２００で、実施の形態１の式（１０）と同様にして、Ｎ個の金属要素に対する行列方程式を作成する。なお、式（１０）の符号を、ここでは式（２３）に示す符号に置き換える。【０１０３】
【数１５】
0004665969.tif 000016
【０１０４】
次に、ステップＳ８００で、孔を検出する。図２９のような単純な構造の場合には、目視によっても検出でき、外部より孔に関する情報を入力することも可能であるが、より複雑な構造の場合には、金属要素構成データより自動的に孔を検出することが望ましい。【０１０５】
以下に孔を自動的に検出するアルゴリズムを示す。
まず、各分割された金属要素（三角要素）に対してエッジ検出を行う。エッジとは他の金属と接していない要素の辺のことである。次にエッジを連結し、エッジの連結から構成される多角形を作成する。図２９の例の場合、二つの円（正確には３６角形）が検出される。この例では、二つの円（正確には３６角形）の内側のみに孔が存在する。検出された多角形内部の任意の点（例えば多角形を複数の分割し、分割された全ての要素の重心をとる）が金属要素上にあるか否かによって孔が存在するかどうかの判別が可能である。【０１０６】
孔が存在する場合、ステップＳ２０２にて、孔に対し、
［Φ］＝−［Φ0］      （２４）
なるファラデーの法則に基いて、行列方程式を作成する。ここで、Φは、孔に鎖交する磁束、Φ0は、孔に鎖交する外部磁束である。ストークスの定理を適用することによって、式（２４）は式（２５）に変形される。【０１０７】
【数１６】
0004665969.tif 000017
【０１０８】
ここで、Ａは、孔の外周を流れる電流が作成する単位電流あたりのベクトルポテンシャルであり、Ａ0は、外部コイルが作成するベクトルポテンシャルである。Ｉholeは孔の外周を流れる孔電流、∫ｄｌは孔の外周に亘る積分を表す。【０１０９】
次に、ステップＳ２０３にて、孔の外周を流れる電流と各金属要素上を流れるループ電流との相互作用を考慮して、金属要素と孔を連成させた連成行列方程式を作成する。孔がひとつのみ存在する例を示すが、孔が複数存在しても同様の考え方が適用できる。【０１１０】
【数１７】
0004665969.tif 000018
【０１１１】
ここで、Ａ1〜ＡNは、Ｎ個の金属要素がそれぞれ、孔に対して作成する単位電流あたりのベクトルポテンシャルである。この単位電流あたりのベクトルポテンシャルを周回積分することにより、単位電流あたりの孔に鎖交する磁束が求まる。Ｇhole1,N+1〜ＧholeN,N+1は、孔電流ＩholeがＮ個の各金属要素上に作成する磁場の法線成分である。
また、Ｇmetal1,1〜ＧmetalN,N、Ｉmetal1〜ＩmetalN、Ｂ０metal1〜Ｂ０metalNは、実施の形態１のｇ（ｒ’、ｒ）・ｎ、Ｉ（ｒ’）、Ｂ０（ｒ）・ｎに相当する。【０１１２】
ステップＳ３００では、前記式（２６）を解くことにより、各金属要素上に流れるループ電流値と孔の外周を流れる電流値とが算出され、電流分布が求められる。ステップＳ４００、Ｓ６００では、前記各実施の形態と同様、ビオサバールの法則によって電磁場分布を算出し、算出した電磁場分布を出力する。【０１１３】
なお、孔が複数（例えば、孔Ａ1、孔Ａ2）ある場合には、式（２６）は式（２７）のように拡張される。【０１１４】
【数１８】
0004665969.tif 000019
【０１１５】
式（２７）において、Ａ10、Ａ20はそれぞれ外部コイルが孔Ａ1、Ａ2に作るベクトルポテンシャル、Ａ11〜Ａ1N、Ａ21〜Ａ2Nは、それぞれ、Ｎ個の金属要素が孔Ａ1、Ａ2に作成する単位電流あたりのベクトルポテンシャルである。この単位電流あたりのベクトルポテンシャルを周回積分することにより、単位電流あたりの孔Ａ1、Ａ2に鎖交する磁束が求まる。また、Ｍ12とＭ21は、孔間の相互インダクタンスである（Ｍ12＝Ｍ21）。【０１１６】
本実施の形態１２の電磁界分布計算方法を実現する装置は、図５に示したブロック図と同様の構成である。孔の検出、孔電流の設定、および孔を含む連立方程式の構築は、図５の構築手段５においてなされる。金属構造要素データ１が装置に入力されると、図２８の各ステップに対応する処理が計算機内のプログラムで実行される。【０１１７】
図３０は本発明の実施の形態１２による他の電磁界分布計算方法を示すフローチャート図である。
まず、ステップＳ１００で、円板金属構造体をＮ個の三角要素に分割し、ステップＳ２００で、Ｎ個の金属要素に対する行列方程式（式（２３））を作成する。
次に、ステップＳ８００で、孔を検出し、孔が存在する場合、ステップＳ８１０にて、孔要素の作成を行う。図３１が作成例である。自動的に三角要素を作成するアルゴリズムは、デローニ法やペーブ法が知られている。図３１は、デローニ法によって自動的に作成した孔要素の例である。
次に、ステップＳ２０２で、Ｍ個の分割された孔要素に対しても金属要素のときと同様の方程式を作成する。すなわち、下記式（２８）を作成する。【０１１８】
【数１９】
0004665969.tif 000020
【０１１９】
ただし、Ｂ０hole1〜Ｂ０holeMは各孔要素における外部磁束密度の法線成分、Ｉhole1〜ＩholeMは各孔要素におけるループ電流である。ここで、１つの孔に属するすべての孔要素に一様のループ電流Ｉholeが流れていると仮定する。すなわち、孔要素内のループ電流の値をすべて同じＩholeと置く。
Ｉhole1＝Ｉhole2＝・・・＝ＩholeM＝Ｉhole    （２９）【０１２０】
すべて同じ値と置くと、アンペールの法則によって孔の外周にＩholeという電流が流れていることと等価になる。この方法によって、孔の周に流れる電流が再現できるため、孔がある場合にも正しい解析が可能となる。
式（２９）を用いて、式（２８）を変形すると【０１２１】
【数２０】
0004665969.tif 000021
【０１２２】
となる。ただし、ｎ＝１〜Ｍ、ｍ＝１〜Ｍ、Ｂ０hole1＝Ｂ０hole2＝・・・＝Ｂ０holeM＝Ｂ０hole
とする。
孔要素が複数存在する場合には、複数の式（３０）式が存在する。
式（３０）式をみると、左辺は、孔の外周が作成する総磁束であり、右辺は外部コイルが作成する総磁束となっている。これはすなわち、ファラデーの法則である。
このように、孔に要素を生成し、生成された孔要素に対し、金属要素と同様の行列方程式を作成し、さらに、１つの孔に属するすべての孔要素に一様のループ電流が流れていると仮定すれば、孔に対して、ファラデーの法則を用いた方程式を導出することができる。
孔に対する方程式を構築する際に、図３０に示す方法を用いた場合、図２８に示す方法に比べると、孔に要素を作成するステップが必要ではあるが、ベクトルポテンシャルを計算する計算時間がいらず、計算が容易となる。また、孔要素を扱う際に、金属要素と同様の取扱いができるため、ソースコードが単純になるというメリットがある。【０１２３】
次に、ステップＳ２０３では、金属要素と孔要素を連成させた行列方程式を作成する。孔がひとつのみ存在する例を示すが、孔が複数存在しても同様の考え方が適用できる。【０１２４】
【数２１】
0004665969.tif 000022
【０１２５】
ステップＳ３００では、前記式（３１）を解くことにより、各電流値が算出され、電流分布が求められる。ステップＳ４００、Ｓ６００では、前記各実施の形態と同様、ビオサバールの法則によって電磁場分布を算出し、算出した電磁場分布を出力する。【０１２６】
なお、孔が複数ある場合、例えば、２つの孔ｈｏｌｅ１（Ｍ１個の孔要素に分割）、ｈｏｌｅ２（Ｍ２個の孔要素に分割）の場合には、式（３１）は式（３２）のように拡張される【０１２７】
【数２２】
0004665969.tif 000023
【０１２８】
ただし、ｎ１＝１〜Ｍ１、ｍ１＝１〜Ｍ１、ｎ２＝１〜Ｍ２、ｍ２＝１〜Ｍ２、Ｂ０hole1、Ｉhole1はｈｏｌｅ１における外部磁束密度及びループ電流、Ｂ０hole2、Ｉhole2はｈｏｌｅ２における外部磁束密度及びループ電流である。【０１２９】
図３２は本発明の実施の形態１２による、さらに他の電磁界分布計算方法を示すフローチャート図である。
図２８、図３０においては、金属要素に対する行列方程式をたてた後に、孔を検出し、孔に対する方程式を立てると共に、金属要素に対する方程式と孔に対する方程式を連成させて、金属構造体全体の連立方程式をたてた。図３２においては、先にステップＳ８００で孔を検出し、ステップＳ８１０で孔要素を生成後、ステップＳ２０３で、孔要素と金属要素とに対し、一度に連成行列方程式の作成を行う。ステップＳ２０３で作成する行列方程式は、図３０のステップＳ２０３で作成したものと同じ方程式となる。すなわち、金属要素に対しても孔要素に対しても、磁界の法線成分が０という境界条件を課した方程式となる。【０１３０】
なお、図３２のステップＳ２０３において、図２８のステップＳ２０３と同様の連成行列方程式を構築しても良い。この場合は、ステップＳ８１０は不要であり、孔要素を作成する必要はない。【０１３１】
実施の形態１３．
前記実施の形態１２では金属構造体が孔を含む場合の例を示したが、前記実施の形態１２のステップＳ２０３で作成する行列方程式の計算に際しては、前記実施の形態４〜７に示した方法と同様の方法を用いて、前記行列方程式をブロック分割し、ブロック毎の電流分布の初期値を求めると共に、各ブロックの相互作用を考慮して前記電流分布を補正することにより、電流分布の計算を行うようにしてもよい。
このようにすることにより、孔を含むものにおいても、計算時間の短縮が可能となる効果がある。【０１３２】
なお、本実施の形態では、実施の形態１２に対し、ブロック分割を行って計算する例、すなわち、
・孔検出→行列方程式を作成（孔と金属要素とに対する行列方程式の作成）→ブロック分割→電流分布の計算
あるいは、
・行列方程式を作成（金属要素に対する行列方程式の作成）→孔検出（孔を含めた連成行列方程式の作成を行う）→ブロック分割→電流分布の計算
の順で計算を行ったが、実施の形態８と同様の方法を用いて、先に金属構造体を複数のブロックに分割し、各ブロック毎に連立方程式を構築しても良い。すなわち、「行列方程式を作成」、「孔検出」、及び「ブロック分割」の順番を入れ替えても計算しても良く、同様の効果がある。
順番を入れ替えた変形例を以下に示す。
・行列方程式を作成（金属要素に対する行列方程式の作成）→ブロック分割→孔検出（孔を含むブロックに対して連成行列方程式の作成を行う）→電流分布の計算
・孔検出→ブロック分割→行列方程式を作成（ブロック毎に孔と金属要素とに対する行列方程式の作成）→電流分布の計算
・ブロック分割→行列方程式を作成（金属要素に対する行列方程式の作成）→孔検出（孔を含むブロックに対して連成行列方程式の作成を行う）→電流分布の計算
・ブロック分割→孔検出→行列方程式を作成（ブロック毎に孔と金属要素とに対する行列方程式の作成）→電流分布の計算【０１３３】
また、本実施の形態１３においては、前記各実施の形態と同様、磁場分布を用いて収束の判定を行っても良いし、電流分布を用いて収束の判定を行っても良い。【０１３４】
実施の形態１４．
図３３は本発明の実施の形態１４による電磁界分布計算方法を示すフローチャート図である。本実施の形態では、孔と閉導体と含む場合にも対応できる計算方法を示すものである。
図において、ステップＳ１００で、金属構造体を要素に分割する。次に、ステップＳ８００で、実施の形態１２と同様にして、孔を検出し、孔が存在しない場合は、ステップＳ７００に進み、孔が存在する場合は、ステップＳ８１０にて、孔内部に要素を作成し、ステップＳ２０３で孔要素と金属要素とに対する行列方程式の作成を行う。
次に、ステップＳ７００では、実施の形態１０と同様にして、閉じた導体を検出し、閉導体が存在しない場合は、ステップＳ３００に進み、閉導体が存在する場合は、ステップＳ７１０にて、閉導体内の１要素を一定値に固定する。【０１３５】
次に、ステップＳ３００にて、最終的に得られた行列方程式を解き、電流分布の計算を行う。本実施の形態において得られる最終的な行列方程式は、金属構造体に孔または閉導体が含まれる場合にも対応した方程式であり、種々の形状に対応した電流分布が計算できる。ステップＳ４００では、得られた電流分布を用いて金属構造体周辺の空間の磁場分布を計算する。ステップＳ６００では、得られた磁場分布を出力する。【０１３６】
なお、前記実施の形態では、電磁界分布を計算する際に、孔検出→行列方程式を作成→閉導体検出→電流分布の計算の順で計算を行ったが、「孔検出」、「行列方程式を作成」、及び「閉導体検出」の順番を入れ替えても良い。
図３４〜図３８は前記順番を入れ替えた変形例である。各図において、計算順所は以下のとおりである。
・図３４：孔検出→閉導体検出→行列方程式を作成（孔と金属要素とに対する行列方程式の作成）→電流分布の計算
・図３５：閉導体検出→孔検出→行列方程式を作成（孔と金属要素とに対する行列方程式の作成）→電流分布の計算
・図３６：閉導体検出→行列方程式を作成（金属要素と対する行列方程式の作成）→孔検出（孔を含めた連成行列方程式の作成を行う）→電流分布の計算
・図３７：行列方程式を作成（金属要素に対する行列方程式の作成）→閉導体検出→孔検出（孔を含めた連成行列方程式の作成を行う）→電流分布の計算
・図３８：行列方程式を作成（金属要素と対する行列方程式の作成）→孔検出（孔を含めた連成行列方程式の作成を行う）→閉導体検出→電流分布の計算【０１３７】
なお、前記図３６〜図３８では、行列方程式をステップＳ２００、ステップＳ２０２、及びステップＳ２０３で作成するが、前記ステップＳ２００では、金属要素に対する行列方程式を作成し、前記ステップＳ２０２では孔要素に対する行列方程式を作成し、前記ステップＳ２０３では金属要素と孔要素に対する行列方程式を作成している。
また、本実施の形態では、孔に孔要素を作成し、ステップＳ２０２、及びステップＳ２０３では式（３０）、式（３１）に従った行列方程式を作成したが、孔要素を作成せず、式（２５）、式（２６）（または式（２７））に従った行列方程式を作成しても良い。【０１３８】
本実施の形態による電磁界分布計算方法を用いることにより、開いた導体のみならず、閉じた導体や孔を含む導体等、種々の形状の金属構造体の電磁界分布解析が可能となる効果がある。【０１３９】
実施の形態１５．
図３９は本発明の実施の形態１５による電磁界分布計算方法を示すフローチャート図、図４０〜図５０は、図３９の各ステップを説明する図である。
本実施の形態は金属構造体が自動車である場合の電磁界分布の最良の計算例を示すものであり、孔検出と閉導体検出とブロック分割と行列方程式の作成とを行うものである。【０１４０】
まず、ステップＳ１００で金属構造体の形状を入力し（図４０）、金属構造体を要素に分割する（図４１）。金属構造体の形状を入力すル方法としては、金属構造体の設計ＣＡＤ情報を変換する方法、金属構造体をＣＡＤソフトを用いて入力する方法、直接座標を手で入力方法などが具体例として考えられる。また、金属構造体を要素に分割する際、図４１では、自動車形状を三角要素に分割しているが、三角要素以外に、四角形要素・多角形要素を用いた要素分割も同じ定式化で扱うことが可能である。自動的に要素分割を行う手法では、デローニ法やペーブ法などが知られている。メッシュ編集ソフトウェアを用いて、手動で要素分割を行うことも可能である。【０１４１】
次に、ステップＳ８００で孔の検出を行う（図４２）。この例では自動車の窓を孔として検出している。この例で用いた形状ほど単純な形状であれば、目視によっても孔の検出は可能である。しかし、より複雑な形状では、実施の形態１２で述べたようなエッジ検出と点の内外判定によって自動的に孔を検出することが望ましい。【０１４２】
次に、ステップＳ８１０で孔要素の作成を行う（図４３）。この例では自動車の窓に三角要素を作成している。自動的に三角要素を作成するアルゴリズムは、デローニ法やペーブ法などが知られている。【０１４３】
次に、ステップＳ３１２で要素のブロック分割を行う（図４４）。この例では、自動車の車体金属構造を３つのブロックＡ、Ｂ、Ｃ（ボンネット、車両上部、車両下部）に分割している。ブロック間の相互作用が小さくなるようにブロック分割を行うことが望ましい。物理的距離が離れている構造体がある場合には、それぞれの構造体を別のブロックに割り当てるとブロック間の相互作用は小さくなる。孔要素と孔を構成する金属要素とは相互作用が強いため、同じブロックに割り当てることが望ましい。図４４では車両上部のブロックＢに、孔要素と孔を構成する金属要素とを割り当てている。仮に別のブロックに割り当てても計算時間が多少長くなるだけで計算は可能である。【０１４４】
次に、ステップＳ２０４で、各ブロック毎に、行列方程式を作成する。この場合、金属要素と孔要素を含むブロックにおいては、金属要素に対しては、金属表面の磁界の法線成分が０という境界条件を課し、孔要素に対しては、孔の閉曲面を貫く鎖交磁束の総和が０であるというファラデーの法則を課す。実際には、実施の形態１２の図３２で説明したように、金属要素に対しても孔要素に対しても、磁界の法線成分が０という境界条件を課し、同じ孔に属する孔要素のループ電流はすべて等しいという条件を課す。得られる行列方程式は、実施の形態１２で述べた式（３２）と同様の形で表すことができる。ただし、式（３２）は孔の数が２個であるが、図４４の例では、孔の数が８個あるので、孔に対する行列が６列、及び６行分追加される。また、この場合のＮは、各ブロックにおける金属要素の数に相当する。孔要素（窓）が合計８枚で６７６要素、車上部の金属要素がＮ＝８３２要素、車下部の金属要素がＮ＝１５３０要素、ボンネット部の金属要素がＮ＝１８９８要素である。１枚の窓には未知数が１個であるから、窓と車上部からなるブロックＢの未知数は、８＋８３２＝８４０個である。したがって、ブロックＢの行列方程式は、８４０×８４０の行列方程式となる。車下部からなるブロックＣの未知数は、１５３０個である。したがってブロックＣの行列方程式は、１５３０×１５３０の行列方程式となる。ボンネット部からなるブロックＡの未知数は、１８９８個である。したがってブロックＡの行列方程式は、１８９８×１８９８の行列方程式となる。【０１４５】
次に、ステップＳ７００では、各ブロックに対して、閉じた導体の検出を行う（図４５）。この例ではブロックＡにのみ閉じた導体構造が存在し、自動車のボンネット部を閉じた導体として検出した。この例では、「エッジを持つ辺を除去する」というアルゴリズムを繰り返し適用し最後に要素が残るか否かで、閉じた導体であるか開いた導体であるかを判定した。この例のように単純化した車モデルの場合には、目視によっても閉じた導体構造を検出できるが、実車規模の複雑な形状においては自動的に検出を行うほうが望ましい。【０１４６】
次に、ステップＳ７１０では、閉じた導体を含むブロックにおいて、閉じた導体を構成する１要素に流れるループ電流の値を固定値にする（図４６）。この例では、ボンネットの上面の真中の１要素１８に流れるループ電流を０に固定した。この段階で、ブロックＡは、１８９７×１８９７の行列方程式であり、ブロックＢは、８４０×８４０の行列方程式であり、ブロックＣは、１５３０×１５３０の行列方程式である。【０１４７】
次に、ステップＳ３３０では、各ブロックの行列方程式を独立に解くことにより、各ブロック毎に電流分布の初期値を求め、ブロック間の相互作用を無視した電流分布を求める。各ブロックの行列方程式は小規模であるため、従来行われている手法（ガウスの消去法、ＬＵ分解法、ＧＭＲＥＳ法など）によって解くことが可能である。【０１４８】
次に、ステップＳ４０１では、得られた各ブロックの電流分布より、ビオサバールの法則を用いて磁界分布を求める（図４７）。ここで求まった磁界分布は、ブロック間の相互作用を無視したものである。【０１４９】
次に、各ブロック間の相互作用を考慮して、得られた電流分布の補正を行う。
まず、ステップＳ２０５で、各ブロック間の相互作用の行列を作成する。各ブロック間の相互作用の行列は、実施の形態８で述べた方法により作成する。すなわち、ブロック間の相互作用を、ブロック番号（ここでは、例えばブロックＡをブロック１、ブロックＢをブロック２、ブロックＣをブロック３とする）を示す添え字をつけたＤで表す（例えばブロック１とブロック２との間の相互作用を表す行列はＤ12となる。前記自動車の場合、Ｄ12は１８９７×８４０の行列で表される。）。【０１５０】
次に、ステップＳ３４２では、ステップＳ２０５で得られた各ブロック間の相互作用の行列を用い、各ブロック間の相互作用を考慮して各ブロックの電流分布を補正する。補正の式は、各ブロック１〜３の電流分布を表す行列をＪ１、Ｊ２、Ｌ３、各ブロック１〜３の外部磁束密度の法線方向成分を表す行列をＥ1、Ｅ2、Ｅ3と表現すると、式（１４）となる。ｋは補正回数（反復回数）である。【０１５１】
【数２３】
0004665969.tif 000024
【０１５２】
ステップＳ４００では、得られた電流分布を用いて、金属構造体周辺の空間の磁場分布を計算する。すなわち、得られた電流分布よりＤｙａｄｉｃ−Ｇｒｅｅｎ関数を用いて、電界分布を求めることによって、電磁場分布が求まる。
ステップＳ５００では、得られた磁場分布と前回の補正計算で得られた磁場分布とを比較して収束判定を行い、収束していない場合は、ステップＳ３４２に戻り、さらに補正計算を繰り返す。収束している場合は、ステップＳ６００で、収束した磁場分布を出力する。【０１５３】
なお、前記実施の形態では、３ブロックに分割する例を示したが、ブロック分割の数は任意である。【０１５４】
また、図４８に１回補正を行った磁界分布、図４９に９回補正を行った磁界分布、図５０に１０回補正を行った磁界分布を示す。
９回補正を行った磁界分布と１０回補正を行った磁界分布がほぼ等しいので、補正を１０回目で終了した。【０１５５】
本実施の形態による電磁界分布計算方法を用いることにより、開いた導体のみならず、閉じた導体や孔を含む導体等、種々の形状の金属構造体の電磁界分布解析が可能となる。
さらに、実車規模の大規模問題でも、ブロックに分割して計算を行うため、必要なメモリサイズが小さくなる。また、計算時間が要素数の２乗に比例するため１万要素以上の計算では従来に比べ１０倍以上高速に計算できる。
また、一般に用いられる電界積分方程式に比べ、行列が実数行列となるため、必要なメモリサイズが半分になり計算が高速になる。【０１５６】
なお、前記実施の形態では、電磁界分布を計算する際に、孔検出→ブロック分割→行列方程式を作成→閉導体検出→電流分布の計算の順で計算を行ったが、「孔検出」、「ブロック分割」、「行列方程式を作成」、及び「閉導体検出」の順番を入れ替えても良い。
図５１〜図５４は前記順番を入れ替えた変形例である。各図において、計算順所は以下のとおりである。
・図５１：孔検出→行列方程式を作成→閉導体検出→ブロック分割→電流分布の計算
・図５２：ブロック分割→孔検出→行列方程式を作成→閉導体検出→電流分布の計算
・図５３：孔検出→閉導体検出→行列方程式を作成→ブロック分割→電流分布の計算
・図５４：閉導体検出→孔検出→行列方程式を作成→ブロック分割→電流分布の計算【０１５７】
なお、図５１〜図５４に示す例以外にも、「孔検出」、「ブロック分割」、「行列方程式を作成」、及び「閉導体検出」の順番を任意に入れ替えて計算しても良く、同様の効果がある。【０１５８】
また、電流分布の計算に際し、ブロック毎の初期値の求め方、反復計算による補正方法は前記方法（ステップＳ３３０、ステップＳ３４２）に限らず、実施の形態４〜７に示した方法を用いても良い。【０１５９】
また、図３９、及び図５１〜図５４に示すフローチャートでは、収束の判定に磁場分布を用いたが、電流分布を収束の判定に用いることもできる。反復計算において、ある反復における電流分布と次の反復における電流分布が十分近ければ収束したとみなすことが出来る。
【図１】
0004665969.tif 000025
【図２】
0004665969.tif 000026
【図３】
0004665969.tif 000027
【図４】
0004665969.tif 000028
【図５】
0004665969.tif 000029
【図６】
0004665969.tif 000030
【図７】
0004665969.tif 000031
【図８】
0004665969.tif 000032
【図９】
0004665969.tif 000033
【図１０】
0004665969.tif 000034
【図１１】
0004665969.tif 000035
【図１２】
0004665969.tif 000036
【図１３】
0004665969.tif 000037
【図１４】
0004665969.tif 000038
【図１５】
0004665969.tif 000039
【図１６】
0004665969.tif 000040
【図１７】
0004665969.tif 000041
【図１８】
0004665969.tif 000042
【図１９】
0004665969.tif 000043
【図２０】
0004665969.tif 000044
【図２１】
0004665969.tif 000045
【図２２】
0004665969.tif 000046
【図２３】
0004665969.tif 000047
【図２４】
0004665969.tif 000048
【図２５】
0004665969.tif 000049
【図２６】
0004665969.tif 000050
【図２７】
0004665969.tif 000051
【図２８】
0004665969.tif 000052
【図２９】
0004665969.tif 000053
【図３０】
0004665969.tif 000054
【図３１】
0004665969.tif 000055
【図３２】
0004665969.tif 000056
【図３３】
0004665969.tif 000057
【図３４】
0004665969.tif 000058
【図３５】
0004665969.tif 000059
【図３６】
0004665969.tif 000060
【図３７】
0004665969.tif 000061
【図３８】
0004665969.tif 000062
【図３９】
0004665969.tif 000063
【図４０】
0004665969.tif 000064
【図４１】
0004665969.tif 000065
【図４２】
0004665969.tif 000066
【図４３】
0004665969.tif 000067
【図４４】
0004665969.tif 000068
【図４５】
0004665969.tif 000069
【図４６】
0004665969.tif 000070
【図４７】
0004665969.tif 000071
【図４８】
0004665969.tif 000072
【図４９】
0004665969.tif 000073
【図５０】
0004665969.tif 000074
【図５１】
0004665969.tif 000075
【図５２】
0004665969.tif 000076
【図５３】
0004665969.tif 000077
【図５４】
0004665969.tif 000078