MWE 2018 FR6A-1

モーメント法入門
Introduction to the Method of Moments
大貫 進一郎Shinichiro Ohnuki
日本大学 理工学部
概要
モーメント法（Method of Moments）は，アンテナ設計や電磁波散乱問題の解析法として広く利用される．
基本的には電磁界の積分方程式に対する数値解析法であり，解析対象物のみを離散化すれば良いこと，吸収境界条件を用いることなく散乱・放射問題を解析できること，などが特徴として挙げられる． 工学系の学部生にとって，マックスウェル方程式を直接的に解析する手法に比べ，モーメント法の物理的および数学的概念は掴みにくいように感じている．また計算コードを作成する際，プログラム上のテクニックや数学的な予備知識が必要となるため，具体的な問題が解析できるようになるまで，ある程度の時間を要している．本基礎講座では，筆者が大学および大学院で行っている電磁気学，電磁波工学，およびコンピュータ演習科目の内容を基に，モーメント法を平易に説明し，初歩的な電磁界解析が実行できることを目的とする．

zθ in
H (i )x

図φin

E (i )y

モーメント法による電磁波散乱解析のイメージ

Abstract
The method of moments (MoM) has been widely used for antenna design and electromagnetic scattering problems.
The MoM is one of the computational techniques for solving integral equations. The advantages of this techniqueare (1) electromagnetic objects are discretized using meshes and (2) absorbing boundary conditions do not need tobe considered for scattering and radiation problems. In this talk, the author explains fundamental concepts of theMoM using basic knowledge of electromagnetics, electrometric waves, and computer algorithms.

1. はじめに電磁界数値解法のひとつとして，モーメント法(Method of Moments) はアンテナ設計から電磁波散乱問題の解析まで，幅広く利用されてきた[1]-[4]．現在
では多くの商用シミュレータに解析ソルバとして搭載されている．また，高速解法の発展とともに，従来は大型コンピュータでしか扱えなかった大規模問題が，デスクトップ PC でも解析できるようになっている[5]-[7]．
マックスウェル方程式を差分法などで直接解く手法に比べ，モーメント法の物理的および数学的概念は，工学系の学部生にとって掴みにくいように感じている．また計算コードを作成する際，プログラム上のテクニックや数学的な予備知識も必要となるため，具体的な問題が解析できるようになるまで，ある程度の時間を要している．しかし，大学学部レベルの電磁気学および線形代数の知識だけで十分に理解できる問題も存在し，科学計算言語 MATLAB などを利用するとプログラムも比較的容易に作成できる．
本基礎講座では筆者の大学および大学院で行われている電磁気学，電磁波工学，コンピュータ演習科目の内容を基に，モーメント法を出来るだけ平易に解説し，初歩的な電磁界解析が実行できることを目的とする．
2. 電磁気学の復習から大学の電磁気学において，科学史的な展開により講義を行う場合，クーロンの法則から説明を始めることが多い[8][9]．
空間座標の原点に置かれた点電荷Q が距離 𝑅𝑅 の点に作る電界 E は次式により与えられる．

E＝Q

4𝜋𝜋𝜀𝜀0 𝑅𝑅 2𝒂𝒂𝑹𝑹(1)ここで，𝜀𝜀0 は真空の誘電率，𝒂𝒂𝑹𝑹 は 𝑅𝑅 方向の単位ベクトルを表す．
点電荷が k 個存在する場合の電界は，それぞれの電荷が作る電界をベクトル的に足し合わせることで求まる．
𝑘𝑘E＝ �𝑖𝑖=1Q 𝑖𝑖

4𝜋𝜋𝜀𝜀0 𝑅𝑅𝑖𝑖2𝒂𝒂𝑹𝑹𝒊𝒊(2)

電荷が領域内で連続的に分布している場合，電荷の代わりに電荷密度 ρ を用いると定式化は便利である．この時，観測点における電界は，電荷が密度 ρで存在する領域 𝑣𝑣 全体に渡る積分で与えられる．

E＝ �𝑣𝑣

ρ 𝒂𝒂𝑹𝑹4𝜋𝜋𝜀𝜀0 𝑅𝑅2𝑑𝑑𝑑𝑑(3)

一方，点電荷 Q が距離 𝑅𝑅 の点に作る電位 𝑉𝑉 については，式(1)の電界を球座標系の半径方向に対して積分することで与えられる．
𝑄𝑄
𝑉𝑉 =(4)4𝜋𝜋𝜀𝜀0 𝑅𝑅

ここで，無限遠点（𝑅𝑅 → ∞）における電位をゼロと仮定した．
有限の体積内に電荷が密度ρで分布する場合，観測点での電位 𝑉𝑉 は電荷が存在する領域を積分することにより求められる．

𝑉𝑉 = �𝑣𝑣

𝜌𝜌𝜌𝜌𝜌𝜌
4𝜋𝜋𝜀𝜀0 𝑅𝑅(5)電荷がある軸方向に対して一様に分布している二次元問題では，軸方向に対する積分を数学的に評価することにより，線電荷に対する表現式が得られる．
密度ρで線電荷が分布している場合の電位 𝑉𝑉 は，平面領域 𝑠𝑠 に対する以下の積分により与えられる．

𝑉𝑉 =1
� 𝜌𝜌 ln 𝑅𝑅 𝑑𝑑𝑑𝑑
2𝜋𝜋𝜀𝜀0 𝑠𝑠(6)3. 電磁気学における積分方程式
入門的な電磁気学の教科書における演習問題では，式(3)や(5)を用いて電界や電位を求める．通常これらの問題では電荷密度 𝜌𝜌 が与えられ，ある領域に対する積分を行う．
【教科書の問題例】
半径 a の円板上に電荷が密度 𝜌𝜌 で一ので様に分布している．このとき円板中心軸上の点 d における電位を求めなさい．
また関連の演習として，図１に示す平行平板コンデンサの問題がある．平板の両端に電圧 𝑉𝑉 の直流電源を繋ぎ，平板間にのみ一様な電界を仮定する．この時の静電容量や電極間に蓄えられるエネルギーを求める．
一方，学生実験などで実施している，より現実的な問題の設定条件は次の通りである．平行平板の両端は電源に接続されているため，平板上での電位は 𝑉𝑉である．しかし，平板端部での電界の乱れは無視できず，平板間だけでなく外部にも電界は存在する．この

yΔi

平板1VdRijyz

𝑑𝑑2x𝑑𝑑2

平板2jx

図 1 平行平板コンデンサの座標系

図 2 平行平板コンデンサのモデル化

条件下においても，電位を求める式(6)は正しいが，積分に含まれる電荷密度 𝜌𝜌 は未知となることが教科
書の問題と異なる．このように未知数が積分に含まれることから，積分方程式と呼ばれる．この電荷分布を数値的に求める手法の一つとしてモーメント法がある．

ここで 𝑅𝑅𝑖𝑖𝑖𝑖 は，区間 i 番目と j 番目の中点間の距離を表す．電位 𝑉𝑉𝑖𝑖 をすべての微小区間(i = 1～2n)に対して書き表すと，以下に示す連立一次方程式が得られる．

4. モーメント法の基礎：静電容量解析
図 1 の平行平板コンデンサについて，その静電容量をモーメント法により求める [10]．平板は z 軸方向に一様な二次元構造を仮定すると，任意の点 (x, y)における電位 V (x, y) は次式により与えられる．
𝑉𝑉(𝑥𝑥, 𝑦𝑦) =1� 𝜌𝜌(𝑥𝑥′, 𝑦𝑦′) ln 𝑅𝑅 𝑑𝑑𝑑𝑑′𝑑𝑑𝑑𝑑′
2𝜋𝜋𝜀𝜀0(7)ここで，𝜌𝜌(𝑥𝑥′, 𝑦𝑦′)は平板上の点(𝑥𝑥′, 𝑦𝑦′)における電荷密度を表し，距離 R は以下により定義する．
𝑅𝑅 = �(𝑥𝑥 − 𝑥𝑥 ′ )2 + (𝑦𝑦 − 𝑦𝑦 ′ )2(8)次に，図 2 に示すような微小区間Δで，平板１を n個，同様に平板 2 を n 個に分割する．平板 1 の左側から順番に 1 から n，平板 2 の左側から順番に n +1から 2n まで番号を割り当てる．また，区間 𝑗𝑗 (: = 1～2n) 内における電荷密度 𝜌𝜌𝑗𝑗 は一定と仮定する．この時，i 番目の区間における電位 𝑉𝑉𝑖𝑖 は，平板上に存在する電荷の重ね合わせにより，次式で与えられる．
2𝑛𝑛

𝑉𝑉𝑖𝑖 = � 𝐴𝐴𝑖𝑖𝑖𝑖 𝜌𝜌𝑗𝑗(9)𝑗𝑗=1𝐴𝐴𝑖𝑖𝑖𝑖 =1� ln 𝑅𝑅𝑖𝑖𝑖𝑖 𝑑𝑑𝑑𝑑′
2𝜋𝜋𝜀𝜀0 𝛥𝛥𝑗𝑗(10)𝐴𝐴⎡ 11𝐴𝐴21⎢⎢ ⋮⎣𝐴𝐴2𝑛𝑛,1

ただし，𝐴𝐴𝑖𝑖𝑖𝑖 =𝐴𝐴12𝐴𝐴22⋮
𝐴𝐴2𝑛𝑛,2𝑉𝑉1𝜌𝜌1⋯ 𝐴𝐴1,2𝑛𝑛⎡ 𝑉𝑉 ⎤⎤
2⋯ 𝐴𝐴2,2𝑛𝑛 𝜌𝜌 2⎥⎥��=⎢⋮ ⎥ ⋮⎢ ⋮ ⎥⋯ 𝐴𝐴2𝑛𝑛,2𝑛𝑛 ⎦ 𝜌𝜌2𝑛𝑛⎣ 𝑉𝑉2𝑛𝑛 ⎦⎧⎪∆𝑗𝑗
ln 𝑅𝑅𝑖𝑖𝑖𝑖 ,2𝜋𝜋𝜀𝜀0𝑖𝑖 ≠ 𝑗𝑗

(11)

⎨ ∆𝑗𝑗 �ln ∆ −1.5�,𝑖𝑖 = 𝑗𝑗
𝑗𝑗
⎪2𝜋𝜋𝜀𝜀0⎩
また行列表現を用いると次式のように表せる．

(12)[A][ρ] = [𝑉𝑉](13)

[ρ] = [A]−1 [𝑉𝑉](14)

ここで，平行平板は電源に接続されているため，右辺の列ベクトル [𝑉𝑉] は既知である．ゆえに，電荷密度[ρ]は，両辺から[A]の逆行列を乗算することで求まる．

各区間の電荷密度 𝜌𝜌𝑗𝑗 が求まると，平行平板コンデンサの静電容量 Cl は式(15)で与えられる．
𝑛𝑛𝐶𝐶𝑙𝑙 = � 𝜌𝜌𝑗𝑗 ∆𝑗𝑗 /𝑉𝑉𝑗𝑗=1(15)図 3 はモーメント法により求めた平板上の電荷密度分布を示す．平板上の電位を一定とした場合，電荷密度は平板上で一様に分布しないことが確認できる．
端部に近づくにつれ密度は高くなり，端部において最大となる．

x 10-11412y[m]100-1-0.5-2

電荷密度[C/m]

30.5-3-1-401234

モーメント法は対象物を微小区間に離散化して解析を行うため，任意形状物体も容易に扱える．図 4 に示す同軸円筒コンデンサの静電容量をモーメント法により求める．この静電容量は数学的に厳密解が求まるため，微小区間の分割数に対するモーメント法の信頼性を検討した．図 5 の横軸は分割数の逆数，縦軸は厳密解との相対誤差を表す．分割数を増やしていくことで，計算結果は真値に近づくことが確認できる．
5. モーメント法の応用：電磁波散乱解析
本学の大学院では，モーメント法の応用例として二次元柱状物体の電磁波散乱問題を解析している．
ここでは，電磁波工学 [11] の講義内容と併せ，電磁波散乱解析に対するモーメント法の定式化を概説す
る．
電磁界が時間因子 exp(𝑗𝑗𝑗𝑗t) で正弦的に変化する
場合，マックスウェル方程式の時間微分項を 𝑗𝑗𝑗𝑗 で置き換えることにより次式が得られる．

5x[m]図 3 平行平板上の電荷密度分布

H 2 = 2[m]H1 = 1[m]𝛁𝛁 × 𝐄𝐄 = −𝑗𝑗𝑗𝑗𝑗𝑗𝐇𝐇z(16)𝛁𝛁 × 𝐇𝐇 = 𝑗𝑗𝑗𝑗𝑗𝑗𝐄𝐄 + 𝐉𝐉(17)𝛁𝛁 2 𝐸𝐸𝑧𝑧 + 𝑘𝑘𝐸𝐸𝑧𝑧 = 𝑗𝑗ωµ𝐽𝐽𝑧𝑧(18)

ここで，𝐉𝐉 は電流密度を表す．電界が z 方向成分のみを持つ E 波の場合，𝐸𝐸𝑧𝑧 に対する波動方程式は以下で与えられる．

ここで波数は 𝑘𝑘 = 𝜔𝜔√𝜀𝜀𝜀𝜀 で定義する．この波動方程式は z 方向に対して一様な線電流 I に対して数学的な解が求まる．𝝆𝝆′点に置かれた線電流が 𝝆𝝆 点に作る電界は次式により与えられる．

図 4 同軸円筒コンデンサの座標系

10-1

𝐸𝐸𝑧𝑧 =

相対誤差

10-2(19)(2)

ここで 𝜂𝜂 は電波インピーダンス，𝐻𝐻0 は外向き放射条件を満足する 0 次の第二種ハンケル関数である．
電流が密度 𝐽𝐽𝑧𝑧 で存在する場合，観測点 𝝆𝝆 における電界は，線電流がつくる電界の重ね合わせで表現でき，式(20)により求められる．

10-310-4
10-5 -510−𝑘𝑘𝑘𝑘 (2)𝐼𝐼𝐻𝐻0 (𝑘𝑘|𝝆𝝆−𝝆𝝆′ |)
410-4

10-310-2
分割数の逆数

図 5 分割数に対する静電容量の変化

10-1𝐸𝐸𝑧𝑧 (𝝆𝝆) =−𝑘𝑘𝑘𝑘(2)
� 𝐽𝐽𝑧𝑧 (𝝆𝝆′ )𝐻𝐻0 (𝑘𝑘|𝝆𝝆−𝝆𝝆′ |)𝑑𝑑𝑠𝑠 ′4(20)

この積分は，線電流が存在する柱状断面に対して実行する．

yEˆ z(i )rEincφinHinczzxyx

図 8 完全導体円柱の座標系

図 6 柱状物体の座標系

まず，導体表面を微小区間 ∆𝐶𝐶 に N 分割し，各区間での電流分布を関数により近似する．ここでは簡単のため，自身の区間においてのみ値が１，それ以外ではゼロとなるようなパルス関数を用いる．

: 領域の中点∆C( X n ,Yn )1 on Δ𝐶𝐶𝑛𝑛
𝑓𝑓𝑛𝑛 (𝝆𝝆) = �0 otherwise

導体表面電流をパルス関数により展開するとRmn𝐽𝐽𝑧𝑧 = � 𝛼𝛼𝑛𝑛 𝑓𝑓𝑛𝑛(24)𝒏𝒏

ここで，𝛼𝛼𝑛𝑛 は重み係数である．
式(24)を(22)に代入すると，積分方程式を以下の連立一次方程式に帰着できる．

( X m ,Ym )図 7 柱状物体のモデル化

モーメント法による数値解析例として，完全導体からなる円柱に，平面 E 波が入射する場合の散乱問題を検討する．完全導体では，表面 𝐶𝐶 において全電界の接線成分がゼロとなる．
𝐸𝐸𝑧𝑧 = 𝐸𝐸𝑧𝑧𝑖𝑖 + 𝐸𝐸𝑧𝑧𝑠𝑠 = 0

(23)

on 𝐶𝐶

(21)

ここで，𝐸𝐸𝑧𝑧𝑖𝑖 は入射電界, 𝐸𝐸𝑧𝑧𝑠𝑠 は散乱電界を表す．円柱
表面に誘起される電流 𝐽𝐽𝑧𝑧 により散乱電界𝐸𝐸𝑧𝑧𝑠𝑠 は作られるとすれば，式(20)を(21)に代入することで次式が得られる．
𝑘𝑘𝑘𝑘(2)
𝐸𝐸𝑧𝑧𝑖𝑖 (𝝆𝝆) =� 𝐽𝐽 (𝝆𝝆′ )𝐻𝐻0 (𝑘𝑘|𝝆𝝆−𝝆𝝆′ |)𝑑𝑑𝑙𝑙 ′4 𝑐𝑐 𝑧𝑧(22)

ここで，入射電界 𝐸𝐸𝑧𝑧𝑖𝑖 は既知である．導体上に誘起される電流 𝐽𝐽𝑧𝑧 は，境界条件を満足するようにモーメント法を用いて決定する．

[𝑙𝑙𝒎𝒎𝒎𝒎 ][𝛼𝛼𝒏𝒏 ] = [𝑔𝑔𝒎𝒎 ](25)

𝑔𝑔𝒎𝒎 = 𝐸𝐸𝑧𝑧𝑖𝑖 (𝑥𝑥𝒎𝒎 , 𝑦𝑦𝒎𝒎 )(26)

ここで(𝑥𝑥𝒎𝒎 , 𝑦𝑦𝒎𝒎 )は微小区間Δ𝐶𝐶𝑚𝑚 の中点を表し，[𝑙𝑙𝒎𝒎𝒎𝒎 ]は次式により定義する．
𝑙𝑙𝑚𝑚𝑚𝑚𝑘𝑘𝑘𝑘
(27)
(2)
=
� 𝐻𝐻 �𝑘𝑘�(𝑥𝑥 − 𝑥𝑥𝒎𝒎 )2 + (𝑦𝑦 − 𝑦𝑦𝒎𝒎 )2 � 𝑑𝑑𝑑𝑑
4 Δ𝐶𝐶𝑛𝑛 0

求めた電流密度を用いて，導体遠方における観測角𝜙𝜙での散乱電界は次式により与えられる．
𝐸𝐸𝑧𝑧𝑠𝑠 (𝜙𝜙)′′= 𝜂𝜂𝜂𝜂𝜂𝜂 � 𝐽𝐽𝑧𝑧 (𝑥𝑥 ′ , 𝑦𝑦 ′ ) 𝑒𝑒 𝑗𝑗𝑗𝑗(𝑥𝑥 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐+𝑦𝑦 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠) d𝑙𝑙′

ただし𝐶𝐶𝐾𝐾(𝜌𝜌) =1

�8𝜋𝜋𝜋𝜋𝜋𝜋𝑒𝑒 −𝑗𝑗(𝑘𝑘𝑘𝑘+3𝜋𝜋� )4(28)(29)図 9 は完全導体円柱の遠方散乱界を規格化周波数ka について計算したものである．ここで k は波数，a は円柱の半径を表す．ka が大きくなるにつれて，円柱後方での散乱波の値が大きくなること，サイドローブの影響が強く表れることを確認できる．
6. まとめ
本基礎講座では，著者が大学および大学院で行っている電磁気学，電磁波工学，コンピュータ演習科目の内容を基に，アンテナ解析や電磁波散乱解析などに広く用いられているモーメント法の基礎を紹介した．工学系の大学においては，電磁界解析を積分方程式により扱う機会が少なくなっているため，モーメント法の物理的および数学的なイメージは掴みにくくなっているように感じている．本講座は，工学系の大学で学ぶ電磁気学の基礎から展開し，モーメント法を出来るだけ平易に解説することを心掛けた．受講頂く方の電磁界解析の理解促進に少しでも貢献出来たら幸いである．

(a) ka = 1

文(b) ka = 590

1560120103015050180330210240300

献[1] R, F. Harrington, Field Computation by Moment Methods, NJ,
USA: IEEE Press, 1992.
[2] 熊谷信昭，森田長吉, 電磁波と境界要素法 (POD 版), 森
北出版， 2014 年．
[3] 徳丸仁, 基礎電磁波，森北出版，1992 年．
[4] C. A. Balanis, Advanced Engineering Electromagnetics, John
Wiley & Sons, 1989.
[5] W. C. Chew, J. M. Jin, E. Michielssen, J. Song, Fast andEfficient Algorithm in Computational Electromagnetics,
Norwood, MA, USA: Artech House, 2001.
[6] W. C. Chew, H. Y. Chao, T. J. Cui, C. C. Lu, S. Ohnuki, Y. C.
Pan, J. M. Song, S. Velamparambil, and J. S. Zhao, “FastIntegral Equation Solvers in Computational Electromagneticsof Complex Structures,” Engineering Analysis with Boundary
Elements, vol.27, no.8, pp.803-823, Sept. 2003.
[7] 大貫進一郎，千葉英利, “アルゴリズムによるシミュレーションの高速化技術,” 電子情報通信学会誌, vol. 96, no. 6,pp. 391-395, 2013.
[8] J. A. Edminister, マグローヒル大学演習 電磁気学（改訂2 版）, オーム社, 2007.
[9] 浜松芳夫，伊藤洋一，山﨑恆樹，大貫進一郎, 一番わかる電磁気学演習, オーム社, 2013.
[10] M. N. O. Sadiku, Numerical Techniques in Electromagnetics
with MATLAB, Third ed., NY, USA: CRC Press, 2009.
[11] 細野敏夫, 電磁波工学の基礎 (POD 版), 森北出版，2015.

270

著者紹介
(c) ka =10
図 9 完全導体円柱の遠方散乱界：平面 E 波入射

大貫 進一郎
日本大学理工学部, 教授, ohnuki.shinichiro@nihon-u.ac.jp