平成２９年度
九州大学大学院 総合理工学府
先端エネルギー理工学専攻

修 士 論 文論文名
アンテナ形状を考慮した
マイクロ波放電式中和器内部の解析コード開発

氏

名

小山

智樹

指導教員名

山本

直嗣

森田

太智

中島

秀紀

目次
第 1 章 序論 ........................................................................................................ 1
1.1.

背景 ............................................................................................................ 1

1.2.

イオンエンジンの原理 ............................................................................... 2

1.2.1. マイクロ波放電式イオンエンジン ......................................................... 3
1.2.2. マイクロ波放電式中和器 ....................................................................... 5
1.3.

研究目的 ..................................................................................................... 7

参考文献.............................................................................................................. 8
第 2 章 数値解析手法 ........................................................................................... 92.1.

解析全体の流れ .......................................................................................... 92.2.

旧解析コードの問題点 ............................................................................. 11

2.3.

本研究の解析の流れ ................................................................................. 11

2.4.

プラズマ粒子計算手法 ............................................................................. 122.4.1. 基礎式 .................................................................................................. 122.4.2. メッシュサイズ・タイムステップ ....................................................... 13
2.4.3. 粒子計算 PIC 法 ................................................................................... 14
2.4.4. 衝突計算 MCC 法 ................................................................................. 14
2.5.

FDTD-PIC-MCC ..................................................................................... 21

2.5.1. プラズマ粒子計算 ................................................................................ 21
2.5.2. 電磁波解析手法 .................................................................................... 24
2.5.3. FDTD-PIC-MCC のフローチャート .................................................... 30
2.6.

FULL -PIC-MCC ....................................................................................... 32

2.6.1. 運動方程式解法 .................................................................................... 32
2.6.2. MPI による領域分割と粒子分割 ........................................................... 34
2.6.3. PIC 法による重みつけ計算 ................................................................... 35
2.6.4. 中心軸近傍における座標変換 .............................................................. 38
2.6.5. 電場解析手法 ....................................................................................... 40
2.7.

FULL -PIC-MCC のフローチャート .......................................................... 44

2.8.

弱連成 ...................................................................................................... 46

2.9.

計算体系および静磁場分布 ...................................................................... 46

2.10.

各種物理条件と計算条件 ...................................................................... 48

2.10.1. 物理条件 ............................................................................................. 48
2.10.2. 計算条件 ............................................................................................. 48

第 3 章 結果と考察 ............................................................................................. 50
3.1.

エネルギー獲得確率 ................................................................................. 50

3.2.

各物理値分布 ............................................................................................ 53

3.2.1. 電子エネルギー分布 ............................................................................. 53
3.2.2. 電子密度分布 ....................................................................................... 543.3.

過剰電離による計算の不安定性 ............................................................... 55

第 4 章 結論 ....................................................................................................... 58
参考文献............................................................................................................. 59
謝辞 ................................................................................................................... 61

第1章 序論
1.1. 背景2010 年に小惑星探査機「はやぶさ」が小惑星イトカワからのサンプルリターンに成功した．続いて 2014 年に後継機「はやぶさ 2」が打ち上げられ，小惑星リュウグウからのサンプルリターンに向け飛行中である．
この「はやぶさ」ミッションにて成功の鍵を握っていたのが，直径10 𝑐𝑚のマイ
クロ波放電式イオンエンジンμ10 である．宇宙にて，長時間安定して動作する動力源として地球から数億 km 離れた小惑星の往復飛行を可能にした．
マイクロ波放電式イオンエンジンとは, 電気推進の一種であり, 電気エネルギーを用いて推進剤を電離させ，生成されたプラズマを様々な方法で加速し，その反作用で推力を得る(1-1)．これに対し，打上げロケットなどに代表される化学推進機は，推進剤の化学反応により生じる熱エネルギーを運動エネルギーに変換し大きな推進力を得る．電気推進機は化学推進機と比較して得られる推力は小さいが，燃費を示す指標である比推力(Specific impulse：𝐼𝑠𝑝 )を一桁以上大きくすることが可能である．
そのため電気推進機は，化学推進機と同じ量の燃料で，より大きな速度に達することができる．電気推進機には，イオンエンジンやアークジェット，MPD スラス
タ，RF スラスタ，パルスプラズマスラスタ，ホールスラスタなどが開発・研究されていが，これらの中でも特に，イオンエンジンは他の電気推進機と比較しても比推力が高く(𝐼𝑠𝑝 数 1000 s 以上)，長期間のミッションである深宇宙探査に適している．その一種であるマイクロ波放電式イオンエンジンは，プラズマ生成に電子レンジと同じ周波数のマイクロ波を使用する長寿命を特徴とした推進機である．「はやぶさ」では JAXA の宇宙科学研究所にて，1997 年からの 18000 時間と 2000 年からの 20000 時間の地上耐久試験に成功しており，マイクロ波放電式イオンエン
ジンμ10 は，一つのシステムとしての積算作動時間 4 台で 40,000 時間に迫り，高い信頼性を示した(1-2)．
しかし，長寿命を特徴とするマイクロ波放電式イオンエンジンにも課題がある．
それは中和器の劣化によってエンジンの寿命が制限されることである．詳しい原理は，次項で説明するが，中和器とは,電子を放出することでイオンビームを中和する装置である．イオンビームを中和しないと, エンジンが負に帯電することで, 排出したイオンが逆流し, 性能が得られないためである．中和器も，内挿されたアンテナによるマイクロ波を用いてプラズマを生成するが，そのアンテナ周りの絶縁体部へ内壁面からのスパッタ粒子が堆積し，絶縁の抵抗値を下げ，その影響により中和器性能が低下する現象が起こる(1-3) (1-4)．その後，中和機能を失ったエンジンは推1

力を長時間発生することは困難となる．そのため中和器の寿命がエンジンの寿命と言っても過言ではない．
以上のような問題を抱えているが，それでもやはり「はやぶさ」，「はやぶさ 2」の活躍により，マイクロ波放電式イオンエンジンの信頼性や実用性の評価は高く，
近年では小型化・大型化，高比推力化や商用化といったμ10 に継ぐ新たなμシリーズの展開が行われている(1-5)．例えば，μ10 の口径を 1/10 にして汎用性を高めた
μ1 の開発である．マイクロ波放電式イオンエンジンμ1 は 20-50kg 級の衛星に分散して搭載し精密な推力制御を行い，軌道修正や姿勢制御を目的としている (1-6)．
さらに，μ10 の口径を 2 倍にしたμ20 は，大推力化を図ることで小天体探査機への応用を目指し開発が行われている(1-7)．このように，既に高い性能が実証された製品を応用し汎用性を高めることで，総合的な開発コストの削減が期待できる．しかし，これらμ1，μ20 は，μ10 を単純にサイズダウン・アップしただけでなく，各部の形状や条件の最適化を必要とされ様々な工夫が施されている(1-6)．さらに，何万時間にも及ぶ耐久試験による寿命評価を行うとなると，結果的に開発までに要する時間や金銭的コストが莫大に掛かってしまう．そこで，新たな製品開発を進めるにあたり，寿命制限の一因とされる“中和器”の性能劣化の遷移を短期間・低コ
ストで評価する寿命評価システムの構築が求められている．

1.2. イオンエンジンの原理
マイクロ波放電式中和器の説明の前に, イオンエンジンの原理を通して, 中和器の位置づけおよび役割を詳しく説明する. イオンエンジンの概略図の図 1-1 に示
す．イオンエンジンはアーク放電やマイクロ波などで推進剤を加熱・電離させてプラズマを生成し，2 枚ないしは 3 枚からなる多孔状の電極に 1000 V 程度の高電圧を印加させてイオンを加速するという静電加速型の推進装置である．イオンエンジ
ンは主に 3 つの領域から構成されている．
① 推進剤を電離するイオン生成部 （Ionization）② 生成されたイオンを静電的に加速して推力を得る加速部 （Acceleration）③ 放出されたイオンビームを電気的に中和する中和部 （Neutralization）
これらの各過程はそれぞれイオン源，加速電極，中和器によって行われる．

2

図 1-1 イオンエンジン概略図

1.2.1.マイクロ波放電式イオンエンジン
イオンエンジンの中でのマイクロ波放電式を説明していく．イオンスラスタはプ
ラズマの生成方法によって直流放電式，RF (Radio Frequency) 誘導放電式，マイク
ロ波放電式等に大別される．本研究の対象であるマイクロ波放電式イオンエンジンの概念図を図 1-2 に示す．マイクロ波放電とは，マイクロ波帯域の交流電場によって電子を加速し，中性粒子と衝突電離により電子の数が増加し気体がプラズマ化
され放電が維持されることである．

図 1-2 マイクロ波放電式イオンエンジン概念図

マイクロ波放電式イオンエンジンはプラズマの生成に電子サイクロトロン共鳴3(Electron Cyclotron Resonance : ECR) を利用し，電子の加熱効率を上げている．
ECR の原理を図 1-3 に示す．真空中に存在する荷電粒子は磁場中でローレンツ力を受ける．この向心力のため磁力線に巻きつくようなサイクロトロン運動と呼ばれる回転運動が現れる．磁界における電子の運動方程式は次式で表される．
dvm
 ev  B
(1.2.1)dt
ここで m は電子の質量，v は速度ベクトル，e は電荷量，B は磁束密度ベクトルである．このときの円軌道の半径𝑟𝐿 はラーマ半径と呼ばれ
mv 
rL 
(1.2.2)eB
で与えられる．ここで v⊥は電子の B に対する垂直な速度成分の大きさである．
プラズマを構成している荷電粒子は全て反磁性体である．そのためサイクロトロン運動の回転の方向は，外部磁場の向きに対して荷電粒子の回転によってできる磁場が常に逆を向く方向である．つまり電子は磁場に対して右回りの回転運動を行う．
この回転運動の角周波数はサイクロトロン角周波数と呼ばれ
eB(1.2.3)
 ce 
mで与えられる．電子は電場と逆向きに加速されるため，磁場中の電子の回転方向と逆向きに，回転する速さが等しい電場をかけると電子は連続的に加速され，電場から効率的にエネルギーを受けとることができる．これが電子サイクロトロン共鳴現
象である．

E

B

電子

E

B

電子

電子
E

B
B
電子

図 1-3 ECR の原理

4E

1.2.2.マイクロ波放電式中和器
マイクロ放電式イオンエンジンには，そのイオンビームに同量の電子を供給し電気的に中和する中和器という装置が搭載される．これにより，エンジンや宇宙機本
体が負に帯電し噴出したイオンを引き戻すことなく，運転を継続することが出来る．
中和器が故障し衛星からイオンのみが放出されることになるとイオンビームの抽出が不可能となる．
一般的な直流放電型イオンエンジンでは中和器としてホローカソードが使用されており，スラスタ本体と同様に電極の損耗や電源構成の複雑化が故障の要因となっている．ヨーロッパで盛んに研究が行われている高周波放電型イオンスラスタ (1-8)
でも中和器にはホローカソードが使用されているが，この問題は解決されていないままである．
このマイクロ波放電は直流放電と比べた場合，以下のようなメリットがある．
1.
2.

ヒーターを用いないため，予備加熱が不要となり素早い起動が出来る．
マイクロ波放電式イオンエンジン駆動用のマイクロ波電源が 1 台あれば中和器のプラズマも同時に生成することが可能なため，システム構成の簡素化や信頼性の向上，また衛星の重量軽減が図れる．

マイクロ波放電式中和器は深宇宙探査機「はやぶさ(1-9)」に搭載され十分な中和性能が得られたという実績がある．当研究グループでもマイクロ波放電型イオンエンジンを開発してきた．図 1-4 はマイクロ波放電型中和器の外観写真である．

5

図 1-4 マイクロ波放電式中和器の写真

図 1-5 マイクロ波放電式中和器運転時

61.3. 研究目的
本研究室では，数値解析によるマイクロ波放電式中和器の寿命評価システムの構築を試みている．数値解析を用いる理由は，正確に短時間で低コストの評価ができるためである．実験による内部計測では，マイクロ波放電式中和器の空間スケールが小さいため，プローブ等の計測機器を掃引すると内部のプラズマの擾乱を引き起こしてしまい，正しい結果が得られにくい．さらに，数万時間にも及ぶ地上試験は，
時間的・金銭的コストが莫大になり，μシリーズのような応用開発の進展に影響を与える．これらの問題を解消し，中和器内部の物理現象を捉えることが可能である
数値解析は有用であると言える．
イオンエンジンの数値解析による研究は多くの研究機関で実施されているが，中性粒子解析，マイクロ波解析，プラズマ解析，損耗解析の異なる時間スケールでの物理現象を統合的に解析するツールの開発はまだ行われていない．本研究室では，笹川により，中性粒子解析と損耗解析を行うツールの解析が行われた．また，楠田らにより後述するマイクロ波解析とプラズマ解析を連成的に解析するツールの開発がなされた．しかし，この連成解析ツールは不完全であり，中和器より引き出される電子電流値が実験値の 1/2 以下となっている．その一因として，旧プラズマ解析コードの計算体系が中和器を正しく模擬できていないことが挙げられる．
図 1-6 に示すように，μ10 や本研究室で使用された中和器は L 字型アンテナを持つ非対称軸の円筒系である．しかし旧コードでは，それを対称として扱い，2 次
元円筒座標系(𝑟, 𝑧)で解いているため，L 字型アンテナからのエネルギー伝搬によるプラズマ挙動を解析できていない．
そこで，本研究では連成解析ツールのプラズマ解析コードを 2 次元(𝑟, 𝑧)から 3 次
元(𝑟, 𝜙, 𝑧)へ拡張し，本解析ツールの妥当性を検討する．

図 1-6 中和器のモデル概略図7参考文献(1-1) 栗木恭一・荒川義博[編]：”電気推進ロケット入門”, (東京大学出版会,2003)(1-2)
細田聡史，國中均，”イオンエンジンによる小惑星探査機「はやぶさ」の帰還運用”，J.Plasma Fusion Res. Vol.86,No.5,(2010),pp.282-292(1-3)
大道 渉， 小泉 宏之， 国中 均, 細田 聡史， 西山 和孝, “マイクロ波放電式中和器の性能低下メカニズムとその実験的検証”，STEP-2010-058(1-4)
大道渉“マイクロ波中和器の内部劣化機構の解明と長寿命化に関する研究”,STEP-2013-008(1-5)
西山和孝，細田聡史 宇宙航空研究開発機構“P4-102 ソーラー電力セイル用イオンエンジンの開発状況”(http://www.isas.jaxa.jp/j/researchers/symp/sss13/paper/P2-127.pdf)(1-6)
杉田裕人，小泉 宏之，国中 均，山極芳樹，松井信，“小型マイクロ波放電
式イオンエンジンにおける中性粒子密度測定”, STEP-2012-042(1-7)
西山和孝，細田聡史，國中均 宇宙航空研究開発機構“P2-178 DESTINY 用イオンエンジンμ20”(http://www.isas.jaxa.jp/j/researchers/symp/sss14/paper/P2-178.pdf)(1-8)
N. Wallace, P. Jameson, IEPC, 2011, 327(1-9)
I. Funaki, H. Kuninaka and K. Toki, J. Propul. Power, Vol. 20, (2004),
pp.718-726.

8

第2章 数値解析手法
2.1.

解析全体の流れ

数値解析手法について, 述べていく. 各計算モデルについて, 詳しく述べる前に,
解析ツールの全体図を述べる. 本研究室にて開発された数値解析コードは, 3 つの物理現象に基づき, 3 種の計算モデルを用いている. 表 2-1 に示す.
表 2-1

物理対象
電子とマイクロ波の相互作用による電子のエネルギー獲得解析
イオンと電子の挙動，静電場解析
中性粒子, スパッタ粒子の輸送および壁面損耗解析

モデル名
FDTD-PIC-MCC
Full-PIC-MCC
DSMC

これら 3 つの計算モデルを弱連成的に解析させることで，プラズマ挙動の解明および寿命の評価ができる．図 2-1 に解析全体のフローチャートを示す．まず，前章で述べたようにマイクロ波放電式中和器の基本原理はマイクロ波と電子のサイク
ロトロン周波数の ECR 共鳴であるため，静磁場の情報が必要となる．中和器に永久磁石を配置した際の静磁場分布は，解析を始める前に市販ソフトで計算しておく(後述)．この静磁場データを持って，以降の FDTD-PIC-MCC および Full−
PIC-MCC の計算が行われる． FDTD-PIC-MCC では，一様状態から電子とマイクロ波の吸収が計算され，電子エネルギー獲得確率が算出される．これは，単位時間単位体積あたりの電子のエネルギー増分を確率密度分布として算出後，電子のマイクロ波吸収量の効果として Full-PIC-MCC に渡される．Full-PIC-MCC では，
イオンと電子を一様状態から計算し，静電場の効果による粒子の輸送を解析する．
プラズマ密度が変化するとカットオフ密度(電磁場がほとんど浸透不可能である密度)以上の分布が出現してきて，密度分布が変化するため，電子のマイクロ波吸収量が変化する．実際に，Full-PIC-MCC の計算において，初期電子エネルギー獲得
確率でエネルギーを与え続けると，特定の領域のみでマイクロ波から電子へのエネ
ルギー伝達が行われ続けた結果として電子温度が上昇し，計算が発散してしまう．
このため，プラズマ密度の変化を FDTD-PIC-MCC に再びフィードバックすることで，カットオフの効果を取り入れ, 電子密度の空間分布によるマイクロ波の伝搬の変化を考慮している．また，エネルギー保存のために，Full−PIC-MCC 計算での電子エネルギーも同様にフィードバックを行う．Full−PIC-MCC の計算において，逐次プラズマ密度の変化とともに，FDTD-PIC-MCC にフィードバックしながら繰り返し計算し，定常状態を目指す．このように，時間スケールの異なる現象過9程をそれぞれ独立させて計算させ，それぞれの結果を受け渡しながら計算する事を
弱連成解析という．本来は，タイムステップ毎に，マイクロ波と静電場を計算すれ
ば，Self-Consistent に数値解析を行うことができるが，時間スケールの異なる現象を計算することは，非常に高い計算コストが生じてしまうため，計算時間を短縮しながら解析物理対象のモデリングを行えるため，弱錬成解析は有用と言える．
Full−PIC-MCC にて定常状態まで解析後，Full−PIC-MCC 算出した各壁面へのイオンのエネルギー分布関数(IEDF)を算出する．この結果を用いて，壁面でのス
パッタ粒子の発生確率を算出でき，DSMC にてスパッタ粒子および中性粒子の輸送を解析する．これにより中和器性能低下の遷移を模擬し，DSMC によって得られたスパッタ粒子の堆積量によって，寿命を算出するものとする．また，アンテナ部へのスパッタ粒子の堆積がマイクロ波の伝搬におよぼす影響を再現するために，
FDTD-PIC-MCC にアンテナの周りに抵抗皮膜を挿入する効果を取り入れることで，従来の定常状態解析から損耗状態解析を行えるものとして，寿命評価だけでなく損耗状態時の中和器内部のプラズマダイナミクスを解析することを可能にしている．

図 2-1 解析全体フロー

10

2.2.

旧解析コードの問題点

本研究室では，これまで各解析コードの開発が行われてきた．さらに，楠田により後述する FDTD-PIC-MCC と Full-PIC-MCC を連成的に解析するツールの開発が行われ，マイクロ波解析，プラズマ解析を行うことが可能となった(2-1)．
しかし，旧 Full-PIC-MCC コードでは，解析領域を 2 次元円筒座標系(𝑟, 𝑧)としているため，FDTD-PIC-MCC と受授する，電子のエネルギー獲得確率，プラズマ
密度や粒子エネルギーが周方向𝜙に平均されてしまう．そのため，L 字型のアンテナを特徴とするマイクロ波放電式中和器の 3 次元非対称性における粒子の挙動解
析が正確に行われていない．さらに，中心軸(z 軸)付近，𝑟𝑐ℎ 以内に存在する粒子の扱いは，𝜙方向の速度成分を”0”として扱い，中心軸𝑟 = 0に到達した粒子の𝑟 方向の速度成分を反転している．これらは，円筒座標系における運動方程式の慣性項に
1 ∕ 𝑟 があり，中心軸付近での計算が発散してしまうためである．この処置よって，𝜙方向の粒子移動が制限されることが懸念される．通常，粒子計算では 1 タイムス
テップで 1 セル以上超えないようタイムステップを設定するため，3 次元における𝜙方向の粒子移動は，1 タイムステップで 1 セル以内である．しかし，中心軸付近においては，中心軸を横切る粒子が存在するため，移動が 1 セル以内であっても，𝜙方向の移動は無視できない．

2.3.

本研究の解析の流れ

本研究の解析フローについて説明する．2.2 で述べたように，アンテナ形状を考慮するため，Full-PIC-MCC を MPI（Message Passing Interface）を用いて𝜙方向の領域分割と中心軸(z 軸)付近の座標変換を行う．それにより，FDTD-PIC-MCC
による電子エネルギー獲得確率を分割した領域数に応じて，より正確に
Full-PIC-MCC に渡すことが可能となる．さらに，Full-PIC-MCC による算出された 3 次元の 電子密度，電子エネルギーを用いて FDTD-PIC-MCC がマイクロ波による電子のエネルギー獲得を解析することで，アンテナ形状を考慮したイオン，電子の挙動の弱錬成解析ができる．本研究の解析フローを図 2-1 に示す．尚，本研
究では，FDTD-PIC-MCC と Full-PIC-MCCC の弱錬成解析のみを扱う．ただし，
両解析コードにおいて用いられる中性粒子密度は DSMC による解析結果を導入している．この DSMC 解析結果の導入に関しては，奥田により，中和器内の中性粒
子密度データを DSMC にて解析し，その結果を両コードに反映させることで，解析精度の向上が図られたためである(2-2)．Full-PIC-MCC から DSMC によるスパッ
タ粒子の解析は，今後開発を行う予定である．

11

図 2-2 本研究解析フロー

2.4.

プラズマ粒子計算手法(2-3)

Full-PIC-MCC も FDTD-PIC-MCC でもプラズマ粒子計算はモデルとしては．
同等のものを利用している．但し，Full−PIC-MCC は円筒座標，FDTD-PIC-MCC
は直交座標を用いている．そのため差分式が異なる．そのため，共通項のみをこの節で述べ，詳しい解法については，各モデルの節で述べる．

2.4.1. 基礎式
プラズマの挙動を解析する際の基礎式として運動方程式を用いる．式を以下に示す．
mdv
 qE  v  B 
dt
dx
v
dt

(2.1)
(2.2)

ここで， m は粒子質量， v は粒子の速度， x は粒子の位置， t は時間， q は粒子の持つ電荷， E は電場， B は磁束密度を表わす．

12

2.4.2. メッシュサイズ・タイムステップ(2-4)
粒子法におけるメッシュサイズは， CFD のような流体計算とはことなり格子形
成法自体が解析の精度に影響をあたえることはなく，対象とする物理現象を再現することができるかどうか次第である．まず，解析したい物理現象に対するメッシュ
サイズを以下のようなデバイ長などから決定し，メッシュサイズからくるスキームの制約からタイムステップが決定されることになる．
計算に用いるメッシュサイズ  は，電子温度 Te，ボルツマン定数 kb，電子密度

ne から計算される．
デバイ長 λDe，
 0 k bTeDe 

ne e 2

(2.3)

に対して，

  De(2.4)

であることが望ましい．また，タイムステップはプラズマ周波数 ωpe，

ne e 2me  0 pe 

(2.5)

に対して，

t   pe1(2.6)

とすれば，数値的な不安定性は生じない．時間積分の精度まで考えると，通常は

t  0.2 pe1(2.7)

が用いられる．
また，磁場中の荷電粒子のサイクロトロン運動にともなう現象をシミュレートするためには，電子サイクロトロン周波数 ωc

c 

eB
me

(2.8)

に対して，

t   c

1

(2.9)

である必要がある．
また，粒子は時間刻みごとに１セル以上飛び越えてはならないという制約がある．
特に荷電粒子は格子ごとに規定される電磁力が働くため，この条件が必須となる．
これを Courant の安定化条件という．3 次元の際，時間ステップ t は，

13

1

Dt £
v

1
1
1
+
+
2
2
(Dx ) (Dy ) (Dz )2(2.10)

を満たさなければならない．この条件は極めて厳しく，わずかでも満足しなければ計算が不安定になる．実際には，この 10 分の 1 程度が望ましいとされている．ここにおいて，𝑣は，粒子の最大熱速度または，電磁波速度つまり光速である．
FDTD-PIC と Full-PIC において，条件が異なるので各モデルの計算条件において説明する．

2.4.3. 粒子計算 PIC 法(2-5)
2.4.3.1.

超粒子(2-5)

粒子法では，粒子ごとに基礎式を記述し，解くものであるが，マイクロ波放電式
中和器では、1 m3 当りに 1017～1018 個が含まれ，計算機全体でも 1012 個あり、この粒子全部を取り扱うことは最新のスーパーコンピュータを用いても不可能であり，実用的でない．そこで粒子シミュレーションで取り扱う粒子として，実際の粒子の電荷質量比を一定に保ったまま多数の粒子の電荷と質量を１つにまとめた超粒子を使用する．この超粒子数に含まれる粒子数が少なければ，少ないほど，ノイズが小さくなるが，計算コストも高くなるので，そのトレードオフになる．超粒子数は，１セルに含まれる粒子数が 10 個以上になる事と取り扱う粒子数よりも超粒子の重みが少ない事（統計性から）から決定する．
2.4.3.2.

重み付け

シミュレーションでは，電場量を空間格子点上で求める.そのため,各粒子の電荷を格子点上に振り分ける必要がある.また,粒子の運動方程式を解く際に,電場解析
によって格子点状に割り振られた値を内挿する必要がある.このような計算を行うことを PIC 法（Particle In Cell method）と呼び，粒子シミュレーションにおいて標準的な手法である．内挿の概念図は，各コードで異なるため，後で各モデルにつて述べる．

2.4.4. 衝突計算 MCC 法
粒子コードにて，衝突を取り扱う場合，粒子の衝突の解析には，
MCC(Monte-Carlo-Collision)法を用いるのが，最速の手法である．粒子同士の衝突判定を各粒子感の衝突周波数から決定する確率と乱数を用いて判定するというのがモンテカルロ衝突法である．以下に詳しく述べる．

14

2.4.4.1. 衝突種別
本モデルにおいて荷電粒子は電磁場の影響を受けて，中性粒子と衝突する．プラ
ズマは電子，イオン，中性粒子からなるとする．電子とイオンの数密度は中性粒子の数密度より非常に小さいため，電子-電子，電子-イオン，イオン-イオンの各衝
突は無視できるものとし，電子-中性粒子のみの衝突を取り扱う．衝突種別は弾性，励起，電離衝突のみを仮定する
2.4.4.2.

平均自由行程と透過度(2-7)

面  ，厚さ l ，その物質の数密度 N の板を考える．エネルギー E の粒子ビームが，面  に垂直に当るものとする．衝突をせずに板を透過してくる粒子数は，入射粒子数より減少する.その減少度は，  を衝突の断面積とすれば，dn  nNdlと表わされ，これよりn  n0 exp Nl (2.11)(2.12)となる．ここで，n 0 は面 （ l  0 ）でのビーム中の粒子である．それゆえ，（ l ,l  dl ）間で最初の衝突が起こる確率は，
Pinitial l dl  exp Nl Ndl(2.13)となり， l 以下で最初に衝突が起こる確率は，
P l    Pinitial l dl
l0 1  exp Nl (2.14)となる．最初の衝突までに粒子が走る平均距離  は Pl の一次のモーメントである．
   lPl dl
0(2.15)
1N
これを平均自由行程という．よってモンテカルロ法によって，任意の点から出発し最初の衝突が起こるまでの距離 l を決めるには，乱数 r を用いて，

r  Pl  l
 1  exp   で表わされ，これを変形して，l  ln1  r (2.16)(2.17)となる．ここで， r は一様乱数であるから 1  r を r とすれば，l  ln(r )(2.18)と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を  とすれば，この粒子が１秒間に気体分子と衝突する回数  は，
15



v


 Nv

(2.19)となり，これを衝突周波数という．
2.4.4.3.

衝突確率・衝突判定

すべての衝突の断面積は，荷電粒子の全エネルギーの関数として表わされる．
   E 
(2.20)
荷電粒子速度がわかれば，荷電粒子と中性粒子の各衝突における衝突周波数がわかり，全衝突周波数  total がわかる．したがって，ある粒子 m の衝突確率 Pcollision,m は以下の式で表わすことができる．

 total  ngas total vm(2.21)

 1  exp  total t 

(2.22)

Pcollision,m  1  exp ngas total v m t 

ここで， ngas は中性粒子の数密度， vm は粒子 m の速度，  total は全衝突断面積である．ある乱数 r1 （ 0  r1  1 ）を用いて， Pcollision,m  r1 のときに粒子 m はタイムステップ t の間に衝突が起こるとする．
2.4.4.4.

衝突種別判定乱数 r2 を用いて以下のように衝突の種類を決める．ここでは例として，電子と中性粒子の衝突について示す．弾性衝突断面積を s elastic ，励起衝突断面積を  excitation ，
電離衝突断面積を  ionization とする．

 0  r2  elastic
 total

であれば弾性衝突

 elastic  excitation
 r2  elastic
 total total elastic   excitation  excitation   ionization r2  elastic
 total totalであれば励起衝突
であれば電離衝突

ここで，s total = s elastic + s excitation + s ionization である．このように断面積の比と乱数を対応させ，衝突の種類を決める．計算では中性粒子としてキセノンを用い，その断面積
データ(2-8)を図 2-3 に示す．今回の計算では，多価電離は考慮していない．

16

図 2-3 電子と Xe の衝突断面積

2.4.4.5.

エネルギーと散乱角

粒子は衝突の前後でエネルギーの受け渡しを行い散乱される．以下にそれぞれの衝突におけるエネルギーと散乱角について述べる．
①弾性衝突
まず，分子同士の弾性衝突を考える．速度 v1 ,v 2  を持つ分子対の衝突後の速度を

v1 ,v 2  とすると，
1
v1  v 2  u2
1
v 2  v1  v 2  u2v1 (2.23)(2.24)と表わすことができる．ここで， u  v 2  v1 （相対速度）である．電子と中性粒子の衝突も同様に考える．衝突前の電子，中性粒子の速度をそれぞれ v e ・ v n の衝突後の速度をそれぞれ v e ・ v n とする．電子と中性粒子の衝突では，電子は中性粒子に比べ質量が非常に小さいため，電子のみが散乱される．中性粒子の速度 v n を無視すれば，相対速度は u  v n  v e  v e としてよい．よって，電子の質量を m ，中性粒子の質量を M とすれば，次式が成立する(2-8)．
v e 

m  M cos M sin 
ve 
h
mM
mM

(2.25)

h の直交座標成分はhx  ver cos hy  

vex vey cos   ve vez sinv er
17

(2.26)(2.27)

v v cos j + ve vey sin j
hz = - ex ezver(2.28)

v er  v ey2  v ez2(2.29)v e  v ex2  v ey2  v ez2(2.30)

ここで，であり，  は散乱角，  は任意の角である．(2.25)式はエネルギー損失を含んでいる．
次に，衝突前後のエネルギー収支を考える．衝突前の電子のエネルギーを Eincident,e とし，衝突後の電子のエネルギーを E scattered, e とすれば，エネルギー収支は以下のように表わされる．

Escattered, e  Eincident,e(2.31)

電子は，エネルギー Eincident,e が小さいとき等方散乱され，大きいときは主に前方散乱となる．これを考慮した散乱角  の確率密度を g   として以下に示す(2-24)．ここで， Eincident,e の単位を eV とする．
Eincident ,e sin 
g   
21  Eincident ,e sin 2 ln 1  Eincident ,e 2
これを解けば，0     

2 1  1  Eincident,e  3
cos   1 
Eincident,er

(2.32)(2.33)となる．この値を(2.25)式に用いて衝突後の速度を求める．このとき，任意の角は  2 r4(2.34)

として与える．
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーを Eexcitation とし，入射電子のエネルギーを Eincident,e ，励起後のエネルギーを E scattered, e ，励
~ とする．エネルギー収支は，起後の速度を v

Escattered, e  Eincident,e  Eexcitation(2.35)と表わされ，励起後の速度はEv~  v e 1  excitation
E incident,e(2.36)となる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失~ の速度で弾性衝突をすると考えることができる．(2.25)式中った後に(2.36)式中 v~
v e に v を用いて，衝突後の速度 v e を求める．このときの散乱角には(2.33)式を用いるが(2.33)式中の Eincident,e には(2.35)式から求められる E scattered, e を用いる．励起エネ

18ルギー Eexcitation は，キセノンの場合 8.34 eV を用いている．
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを E scattered, e ，生成電子のエネルギーを E created, e ，電離エネルギーを Eionization とする．エネルギー収支はEscattered, e  Ecreated, e  Eincident,e  Eionization(2.37)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる．

 Eincident,e  Eionization Escattered, e  B tanr5 tan-1 
2
B




(2.38)

ここで，B はガスに固有の値でキセノンでは 8.7 eV となる．
また，入射電子が電離によって失うエネルギーを E とすると，
E  Eincident,e  Escattered, e(2.39)と表わされる．(2.36)式から電離後の電子の速度は，v~  v e 1 EEincident,e(2.40)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．(2.27)式
中 v e に v~ を用いて，衝突後の速度 v e を求める．散乱角には(2.33)式を用いるが，(2.33)式中の Eincident,e には(2.38)式から求められる E scattered,e を用いる．電離エネルギー
Eionization は，キセノンの場合 12.13 eV である．
電離後の生成電子の持つエネルギーは，E created, e  E  E ionization(2.41)となる．また，入射電子の速度 v e を用いて生成電子の速度はv created, e veve

2 Ecreated, em(2.42)と表わされる．この速度で弾性衝突をすると考えれば，v created, e を v e の代わりに(2.27)
式に代入して，衝突後の速度を求める．このときの散乱角は生成電子のエネルギー
E created, e を用いて(2.33)式から求める．

2.4.4.6.

Null-collision 法(2-11)(2-12)(2-13)

計算時間短縮のため，Null-collision 法を用いる．これには，架空の衝突断面積
 fake を導入する．すべての電子のエネルギー E に対して
 fake   total(2.43)であるとし，例えば電子と中性粒子の衝突においては total   elastic   excitation   ionization   fake(2.44)

19

とする．ここで， totalv collisionn gas

 constant(2.45)

を満たす  fake を導入すれば，  collision が一定，すなわち Pcollision,m が一定となりエネル
ギーに依存しない．したがって，粒子毎に Pcollision,m を計算する必要が無くなり，全粒子に対して衝突計算を行わず，衝突を起こす粒子のみについて衝突の種類を決めればよい．これにより，計算時間を大幅に短縮できる．衝突の種類を決めるときは，
2.1.8 節において，  total   elastic   excitation   ionization   fake とし， elastic   excitation   ionization r2  1のとき Null-Collision Process totalとする．
Null-Collision Process が選ばれたときには，粒子の運動は何ら変化しないものとする．図 2-4 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度でn
規格化した衝突周波数を示す．本計算では Null-Collision 法で用いた collision gas の値は 6.5×10-13 m3/sec である．

図 2-4 Xe の数密度で規格化した電子の Xe の衝突断面積

2.4.4.7. タイムステップの制約
モンテカルロ衝突法では，タイムステップの間にそれぞれの粒子についての衝突
を 1 回ずつしか処理できない．しかし，有限のタイムステップを用いる限りは，1
タイムステップの間に複数回の衝突が含まれうる．複数回の衝突による衝突のミスを減らすためにはタイムステップを十分小さくしなければならない．
ある粒子のタイムステップあたりの衝突ミス回数は

20

¥

P i2
r ~ åP =(2.46)
1- P ik=2であるので，衝突ミスの回数を１％いかにするには 𝑃𝑖 ≦ 0.1 とする必要がある．
nullcollision を用いる場合は，𝑃𝑛𝑢𝑙𝑙 ≦ 0.1 すなわち
0.1
Dt £(2.47)n nullki

2.4.4.8. 乱数生成
乱数は(0,1)間に一様に分布する乱数である．生成周期の短い fotran の組み込み
乱数 ran()の使用は避け Mersenne Twister(2-9)を使用した．周期は(219937-1)であり，生成速度も早い．

2.5.

FDTD-PIC-MCC

FDTD-PIC-MCC は，電子のみを粒子として扱い，イオンと中性粒子は背景としておいておき，電子の時間スケールにて３次元的に解く．後述する FDTD(Finite
Difference Time Domain)法によりマイクロ波解析を行い，電場 E の効果として，
PIC 法にて内挿し，電子の電磁運動方程式を解くものである．さらに， 前述の衝突計算を MCC 法にて実施したものである．以下に詳しく述べる．

2.5.1. プラズマ粒子計算
2.5.1.1.

運動方程式解法

基礎式は，2.4.1 で述べた(2.1)式に leap-frog 法による時間中心差分を施すと，
x n 1  x n
 v n 1 2
(2.48)
tとなる．ここで，各変数の上付き添え字は時間に関するもので， n はタイムステッ
プである．(2.1)式に時間中心差分を施すと，
v n 1 2  v n-1 2 q  n 1 n 1 2
  E  v
 v n-1 2  B n 
t
m
2


(2.49)

である． n-1 2 はある時刻 n から t 2 だけ以前の値であることを示す．同様に，n  1 2 はある時刻 n から t 2 だけ進んだ値であることを示す．

プラズマ粒子の運動を求めるために， (2.49)式において未知量である v n 1 2 について解き，求めた v n 1 2 と(2.48)

式より新たな粒子の位置 x n  1 を求めることができる． (2.49)式を変形すると，
Av n 1 2  Sと書ける．ここで， a ・  ・  ・  を次のように定義する．
qt n

Bx
2m

21

(2.50)(2.51)

q t nBy
2mq t n Bz
2m  1 2   2   2

(2.52)(2.53)(2.54)

なお，ここでの E ， B は粒子位置におけるそれぞれの値を表わしている．これらを用いて行列 A とベクトル S は，以下のように構成される．

 1

A 
 



1



 

 
1 

 v xm  cE xn  v ym  v zm S   v ym  cE yn  v zm  v xm  v m  cE n  v m  v m 
z
x
y 
 z

(2.55)(2.56)

ここで， c  qt m ， m  n  1 2 である．さらに (2.50)式を次のようにして変形し，速度を求めることができる．

v n+1 2 = A -1 S
 1   2       1
A 1      1   2    2        1   
2.5.1.2.

(2.57)(2.58)

PIC 法による重みつけ計算

3 次元における重みつけ方法を示す．粒子の位置から近傍のセル座標 (i , j , k )を求める．次に，図 2-5 に示す dx1, dx 2, dy1, dy 2, dz1, dz 2 を求める．その値を用いて，図中にあるような直方体 V1〜V8 の体積を求める．この体積比から粒子が感じる電磁場を内挿することができる．
具体的には，座標 (x , y , z )上にある粒子が感じる電場 E (x , y , z )は，計算セルの体積
を V ，図に示すような粒子によって仕切られた空間体積をそれぞれ V1～V8 とすると，格子点上の重みは

An =

Vn
V

(n = 1 〜 8)

(2.59)となり，
E ( x, y, z ) = A1E (i, j, k ) + A2 E (i +1, j, k ) + A3 E (i +1, j, k +1) + A4 E (i, j, k +1)+ A5 E (i, j +1, k ) + A6 E (i +1, j +1, k ) + A7 E (i +1, j +1, k +1)+ A8 E (i, j +1, k +1)と表わされる．

22(2.60)

図 2-5 PIC 法におけるセルの構成(3 次元)

2.5.1.3.

電荷密度および電流密度計算

電荷密度，電流密度はそれぞれ図 2-5 に示されるセルの格子点上で定義され，セル中に N 個の粒子を代表した 1 個の超粒子が存在したとき，格子点 (i, j, k )上に割り振られる電荷密度は，

V1qN(2.61)
V2となる．すべての粒子について和をとれば，格子点状での電荷密度が得られる．また，電流密度はj  qnv(2.62)r (i , j , k ) =より，j (i , j , k ) =

V1
qNv
V2

23

(2.63)

2.5.1.4. 粒子の壁面境界条件
FDTD-PIC-MCC では，静電場の効果を考慮していないため，電子のシース電位による反射を仮定する必要がある．シース電位については，Full−PIC-MCC にて解析可能であるが，現在電位分布の影響を連成することはできていない．そのため，実験値と過去の数値解析の結果から 20 V と設定し, 壁面に達した場合，一律で20 eV 以上は消失し，それ以上は反射するものした．また，Plume 領域 引き出し電流として，排出される出口に達した場合，無条件で消失するものとした．また,このシース電位の効果を弱錬成するときに渡すべきであるが、シース電位の判定が難しいため, 連成は行っていない.
2.5.1.5. 電子のエネルギー獲得確率計算方法
Full−PIC-MCC との連成で重要となる電子エネルギー獲得確率の計算方法を示す．1 つの観測時間をマイクロ波の１周期とし，タイムステップの間つまり運動方
程式を解いた前後で変化したエネルギー増分をカウントし，エネルギー増分のヒス
トグラムとして，確率密度分布とする．マイクロ波解析は計算コストの兼ね合いとある程度平均化された結果がでる 80 周期分を計算しているが，最初の 40 周期分は計算上の誤差が大きいため，40 周期から 80 周期分の時間平均をとっている．また，各エネルギー増分はセル毎で観測されており，粒子の位置情報をセルに割り当てている． ヒストグラムのステップ数は，0.1 eV，エネルギー範囲としては，-50〜50 eV とした．

2.5.2. 電磁波解析手法(2-14)(2-15)
2.5.2.1.

基礎式

FDTD 法は数値電磁場解析の分野において，モーメント法や有限要素法と同様に，アンテナや平面回路など数多くの複雑な問題に適用され，主要な解析手法として多く用いられている．計算には，マクスウェル方程式を基礎式として用いる．以下に基礎式として用いる微分形式のマクスウェル方程式を示す．
H  E   0
(2.64)tE  H  0
J
(2.65)t
ここで， E は電場， H は磁場，  0 は真空の透磁率，  0 は真空の誘電率， J は電流密度を表わす．

24

2.5.2.2.

差分式

三次元の FDTD 法で用いられる電場，磁場の空間配置を図 2-6 に示す．このように，電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．これは Yee によって提案され(2-16)，電場の回転（   E ）が磁場を，磁場の回転（   H ）が電場の方程式を満たすように配置されている．
この電磁場の配置のもと， (2.64)式, (2.65)式を空間，時間とも中心差分した式を以下に示す．((2.66)~(2.71))式ここで， x, y, z は空間幅， t は時間幅，添え字

i, j, k  はセル中での座標，添え字 n は時間ステップを表わす．
n+

1

n-

1

H x 2 (i, j +1 2, k +1 2) = H x 2 (i, j +1 2, k +1 2) +Dt
m 0 (i, j +1 2, k +1 2)

ìï Eyn (i, j +1 2, k +1) - Eyn (i, j +1 2, k ) Ezn (i, j, k +1 2) - Ezn (i, j +1, k +1 2) üï
+
í
ý
Dz
Dy
îï
þï
n+

1

n-

Dt

1

H y 2 (i +1 2, j, k +1 2) = H y 2 (i +1 2, j, k +1 2) +m 0 (i +1 2, j, k +1 2)

ïì Ezn (i +1, j, k +1 2) - Ezn (i, j, k +1 2) Exn (i +1 2, j, k ) - Exn (i +1 2, j, k +1) ïü
+
í
ý
Dx
Dz
îï
þï
n+

1

n-

1

H z 2 (i +1 2, j +1 2, k ) = H z 2 (i +1 2, j +1 2, k ) +Dt
m 0 (i +1 2, j +1 2, k )ïì Exn (i +1 2, j +1, k ) - Exn (i +1 2, j, k ) Eyn (i, j +1 2, k ) - Eyn (i +1, j +1 2, k ) ïü
+
í
ý
Dy
Dx
ïî
ïþ
t
E xn 1 i  1 2 , j , k   E xn i  1 2 , j , k  
 0 i  1 2 , j , k 1n n  12
 H z i  1 2 , j  1 2 , k   H z 2 i  1 2 , j  1 2 , k 

y


n



1

n

1

H y 2 i  1 2 , j , k  1 2   H y 2 i  1 2 , j , k  1 2 z

E yn 1 i , j  1 2 , k   E yn i , j  1 2 , k  



Hz

1
2

(2.67)(2.68)(2.69) J x i  1 2 , j , k t

 0 i , j  1 2 , k 1n n  12
 H x i , j  1 2 , k  1 2   H x 2 i , j  1 2 , k  1 2 

z


n

(2.66)n12

i  1 2 , j  1 2 , k   H z i  1 2 , j  1 2 , k x25(2.70) J y i , j  1 2 , k 

E zn 1 i , j , k  1 2   E zn i , j , k  1 2  t 0 i , j , k  1 2

1n
 n  122i  1 2 , j , k  1 2
H
i

1
2
,
j
,
k

1
2

H
 y
y

x





H

1
n
2
x

i , j  1 2 , k  1 2  H

1
n
2
x

i , j  1 2 , k  1 2y(2.71)
 J z i , j , k  1 2 ((2.66)~(2.71))式から，時間ステップが n1  1 2 での磁場 H

n

1
2

は，その半ステップ前n

での電場 E n と１ステップ前の磁場 H 2 から求めることができ，時間ステップが1n
n  1 での電場 E n1 は，半ステップ前の磁場1 H 2 と１ステップ前の電場 E n から求めnることができる．初期値については， H 2 は E n から計算されるので，アンテナの動径方向に E rn だけを与える．

図 2-6 FDTD 法における電場と磁場の配置

26

2.5.2.3. アンテナ解析
本計算ではアンテナを完全導体で構成されているものとしてモデル化する．給電方法は同軸ケーブルを模擬し, 図 2-7 に示すように動径方向の４つの電場を中心導体に与える同軸線路給電を用いる.また, アンテナ導入部には図 2-8 のように
Mur の一次吸収境界条件を与え, 波長に比べて太さが無視できるくらい十分細い
直線状のアンテナを解析する場合，その中心軸に沿って電界を 0 とすればよい．また, 今回は L 字型アンテナを模擬した．給電方式について詳しく述べると，中心導体に与える電場の𝑥方向成分を Ex，𝑦方向成分を Ey は近似的に
𝑉
𝛥𝑥𝑟𝑙𝑛 ( 𝑎 )𝑉𝐸𝑥 =𝛥𝑦𝑟𝑙𝑛 ( 𝑎 )𝐸𝑥 =(2.72)(2.73)

で与えられる. ただし, V は中心導体と外部導体間の電圧, 𝑟 は動径方向の距
離, 𝑎 はアンテナの半径である.
また, 同軸ケーブルのインピーダンス Z は以下の式で表される.
𝑍=√

𝜇 1 𝑅
ln
𝜀 2𝜋 𝑟(2.74)

ここで, R は外部導体の半径, r は中心導体の半径である. マイクロは電力を P とす
ると，
𝑃=

𝑉2
2𝑍

(2.75)となり，(2.74)式，(2.75)式より，V を求め，(2.72)式，(2.73)式よりアンテナに与える電場を求める．ここで，Ex を 𝑟 = 𝛥𝑥/2の電場で，Ey を 𝑟 = 𝛥𝑦/2 の電場で代表させると，Ex (IFED, JFED, KFED) = -E x (IFED -1, JFED, KFED)

V (nDt)Dx æ Dx ö
ln ç ÷
2 è aø
Ey (IFED, JFED, KFED) = -E y (IFED, JFED -1, KFED)=-

V (nDt)Dy æ Dy ö
ln ç ÷
2 èaø
ここで，式中の(𝐼𝐹𝐸𝐷, 𝐽𝐹𝐸𝐷, 𝐾𝐹𝐸𝐷)は電場を励振する座標を与える．
=-

27(2.76)(2.77)

図 2-7 同軸線路給電

図 2-8 アンテナモデル

2.5.2.4. PEC(Perfect Electric Conductor )境界条件
境界面が金属のような完全導体である場合，境界に平行な電場成分を 0 とする
PEC 境界条件を用いる．放電室内壁は導体を仮定しているので，この境界条件を
適用する．

28

2.5.2.5.

Mur の吸収境界条件(2-17)

開放領域を取り扱う場合，計算領域を仮想的な境界で閉じなければならない．この境界を吸収境界といい，この吸収境界が完全でなければ，反射波が解析領域に戻ってしまい，誤差の原因となる．ここでは，定式化が容易で広く使用されている
Mur の吸収境界条件を用いる．

x = 0 の吸収境界に x の正の方向から E z 成分を持つ平面波が入射する場合を考える．伝播速度を v とすると， x の負の方向に進む波はEz  Ez x  vt 

(2.78)と表わされ，次の微分方程式を満足する．
Ez 1 Ez
0(2.79)
x v tx = 0 の境界で完全に反射がないならば，電場は(2.78)式の形を保ったまま伝播するはずである．すなわち，境界でも(2.79)式が満足される．これを FDTD 法の形式に定式化するため，時間について差分化すると¶E z E zn - E zn -1
=
¶t
Dt

1
(2.80)n¶E z 2=v¶xとなる．右辺の x についての微分は，電場が x = 0, Dx に割り当てられていることから， x = Dx 2 で行い，(2.79)式に代入すると
E i, n  1 2  E
t
n
r

n1
r

n

1
2

n

1
2

i, n  1 2  v Er i, n   Er i, n  1z(2.81)となる. (2.80)式の E zn 1 2, E zn 1 1 2 も FDTD 法では直接求めることはできないため，前後の値で平均して求める．そして， E z についてまとめるとEzn ( 0, j, k +1 2) = Ezn-1 (1, j, k +1 2)+

vDt - Dx n
{Ez (1, j, k +1 2) - Ezn-1 ( 0, j, k +1 2)}
vDt + Dx

(2.82)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．
2.5.2.6. ポインティングベクトル
空間を流れる電力を求めるためには，ポインティングベクトルを使用する．
Maxwell の方程式を変形すると，以下の式となる．

  E  H   dS   J  Edv 
s

vol

29

  E 2 H 2   2  2 dv
t vol(2.83)

右辺の第 1 項は体積内で消費されるオーム性の全電力，右辺第 2 項は電界と磁界に蓄えられている全エネルギーであり，時間に関する偏導関数をとることによって，この項は体積内に蓄えられたエネルギーの時間増加率，すなわち体積内の蓄積エネ
ルギーを増加させる瞬時電力となる．よって左辺は，この体積に流れ込む全電力でなければならないから，その体積から流れ出る全電力となる．外積 E×H はポイン
ティングベクトルと呼ばれ，W/m2 の単位を持つ瞬時電力密度である．この式をアンテナ上で計算することにより投入電力を測定した．
2.5.2.7.

投入電場の補正(2-28)

投入マイクロ波電力の制御因子は, 先に示したアンテナ根本に励振する電場強度である.時間発展するにつれ,電子 の状態は変化するため,アンテナから見た負荷
インピーダンスは変化し,一定の電場を根本に与えた場合,投入 されるマイクロ波
電力も時間とともに変化する.本計算 では,一定のマイクロ波電力を放電室内に投入するために,以下の式にしたがって,根本に与える電場強度をマイクロ波 1 周期毎に更新した.
𝑃𝑡𝑎𝑟𝑔𝑒𝑡𝜃(− 1) 𝐸𝑜𝑙𝑑𝑃𝑜𝑙𝑑(2.84)
𝐸𝑛𝑒𝑤 = 𝐸𝑜𝑙𝑑 +2
ここで，𝐸𝑜𝑙𝑑 , 𝐸𝑛𝑒𝑤 は，更新前後の励振電場強度，𝑃𝑜𝑙𝑑 は更新前のマイクロ波 1 周期で平均した投入マイクロ波電力，𝑃𝑡𝑎𝑟𝑔𝑒𝑡 は設定投入マイクロ波電力とする．

2.5.3. FDTD-PIC-MCC のフローチャート
図 2-9 に FDTD-PIC-MCC のフローチャートを示す．静磁場とプラズマ密度は事前に計算したものを使い，初期化処理の際，粒子の配置を密度分布に従って作成す
る．FDTD 部でマイクロ波の供給を行いながら，シームレスにマクスウェル方程式を解き，電磁場を求める．その後，PIC のタイムステップのタイミングになった
時に，PIC 部で運動方程式を解いた後，衝突計算を実施する．粒子の移動の更新後，格子点に電流密度を割り当て，それを FDTD 計算に渡す．最終的に繰り返し回数が最大になれば計算を終了する．

30

図 2-9 FDTD−PIC-MCC のフローチャート

31

2.6.

Full-PIC-MCC

本 Full-PIC-MCC コードは，旧コードの問題点を解決するために，MPI を用いて𝜙方向に領域分割することで，解析領域を 3 次元円筒座標系(𝑟, 𝜙, 𝑧)とし，イオンと電子を共に粒子として扱い計算している．さらに中心座標付近においては，直交座標系に座標変換を行い，運動方程式を解いている．
粒子計算において，時間スケールの異なるイオンと電子を共に粒子として扱うのは，計算コストが非常に高くあまり主流ではないが，その代わりにポアソン方程式
を解くことが可能になり，静電場解析を可能にしている．このモデルの対象とする物理現象は，シースの解析であるため，メッシュサイズがデバイ長に近い非常に細かいものとなっている．メッシュからくるタイムステップの制約の緩和のために，後述するが電子の質量比の変更を行い，タイムステップをイオンスケールに近づける工夫をしている．また，衝突計算には，MCC 法を用いている．以下に詳しく述べる．

2.6.1. 運動方程式解法
基礎式は，2.4.1 で述べた．(2.1)式を円筒座標系で書き下すと以下のようになる．
ただし，本計算において，モデル簡易化のために， 𝐸𝜙 = 0, 𝐵𝜙 = 0としている．尚，
中心軸近傍である r <= 𝑟𝑐ℎ については後述する．

v2
dvr q
 Er  v Bz  
dt m
r
dv q
vv
 E  v z Br  vr Bz   r 
dt
m
r
dvz q
 E z  v Br 
dt m

(2.85)

これらを四次精度 Runge-Kutta 法にて解く．
tbr1  2br 2  2br 3  br 4 6tb1  2b 2  2b 3  b 4 vn 1  vn 
6tbz1  2bz 2  2bz 3  bz 4 v zn 1  v zn 6
但しここで定数 b はそれぞれ以下を用いる．
vrn 1  vrn 

32(2.86)qEr  v Bz   1 v v
mr
tb1   1 tb1 tb1 q
 Bz    v  v br 2   Er   v 
m
2
r
2
2

 



br1 

tb 2   1 
tb 2 tb 2 q Bz    v  v br 3   Er   v m
2   r
2 
2 


(2.87)qEr  v  tb 3 Bz  1 v  tb 3 v  tb 3 mr

br 4 

q
E  vz Br  vr Bz   vr vmrqtbz1 
tbr1   1 
tbz1 
tbr1 b 2   E   v z  Br   vr 
 Bz    v z  vr m2 2   r
2 
2 


b 1 b 3 

q
tbz 2 tbr 2   1 
tbz 2 tbr 2 (2.88) Br   vr 
 Bz    v z  vr  E   v z m2 
2   r
2 
2 



b 4 

q
E  vz  tbz 3 Br  vr  tbr 3 Bz   1 vz  tbz 3 vr  tbr 3 mrqEz  v Br mtb 1  q Br bz 2   E z   v 
m
2

 
bz1 

tb 2  q Br bz 3   E z   v 
m
2

 
q
bz 4  E z  v  tb 3 Br m(2.89)

以上の解から，時間発展による位置の更新を行うと，

𝑟 𝑛+1 = 𝑟 𝑛 + 𝑣𝑟𝑛+1 ∆𝑡
𝑣𝜙𝑛+1 ∆𝑡𝑛+1𝑛𝜙
= 𝜙 + 𝑛+1
𝑟
𝑧 𝑛+1 = 𝑧 𝑛 + 𝑣𝑧𝑛+1 ∆𝑡(2.90)となる．ここで，𝜙 𝑛+1 が領域を超える，つまり領域移動する粒子が存在する場合，それらの粒子を対象に MPI を用いて，粒子情報を次の領域に渡す．

33

2.6.2.

MPI による領域分割と粒子分割

MPI(Message Passing Interface)とは，並列計算機間の通信ライブラリである．
MPI を用いることで，複数の計算機でプログラムを並列に実行できるため，主に計算時間の短縮を目的として利用される．本コードでは，この MPI を用いて領域
分割，粒子分割を行なっている．図 2-10 は中和器の上面(𝑟𝜙面)の概略図である．
図 2-10 に示されるように，𝜙方向に領域分割し，各領域が，本コードの主となるプログラムを並列実行する．ここで各領域が扱う電子のエネルギー獲得確率を変化させ，𝜙方向の粒子運動を解くことで，アンテナ形状を考慮した粒子の挙動を解析することができる．さらに，各領域では粒子分割がなされ，プログラムの実行単位であるプロセス（図 2-10:process）数に応じた粒子数が振り分けられ実行される．
これにより計算時間の短縮を図っている．
しかし，本コードの領域分割・粒子分割において，MPI 通信の対象となるプロセ
スが異なるため，コミュニケータを使い分ける必要がある．ここで，コミュニケー
タとは MPI における操作の対象となるプロセスの集合体を意味する．領域分割における MPI の用途は，主に𝜙方向の粒子移動による粒子情報の通信であるため，
全プロセスが MPI の対象となる．一方で領域内の粒子分割における MPI の対象となるプロセスは，領域に依存するプロセスのみとなる．そのため，通常 MPI ではプログラム起動後に MPI_COMM_WORLD というコミュニケータが用意されるが，このコミュニケータは全プロセスを集合体として扱うため，各領域内のプロセスを操作するための新たなコミュニケータを作成した．
そこで，本研究で用いたのが MPI_COMM_SPLIT 関数である(2-18)．この関数を使用し，任意の新たなコミュニケータ NEW_COMM を作成し，各領域内のプロセス間の MPI 通信を可能とした．概念図を図 2-11 に示す．この関数を用いる際に，
color と key という変数の値を定義する．両者とも rank と呼ばれる各プロセスに割り当てられる ID 番号に与える値である．color は領域数を決定する値であり，同じ値を持つプロセスの集合体を扱う新たなコミュニケータを生成する．本コード
では，ここで与えた color 値が領域番号となる．key は，生成されたコミュニケー
タに割り当てられる新たな rank 番号を決定する値である．
これらより，本コードでは，領域間は MPI_COMM_WORLD を，領域内は
NEW_COMM を使用することで領域分割と粒子分割を行う．

34

図 2-10 領域分割モデル

図 2-11 コミュニケータ分割の概念図

2.6.3. PIC 法による重みつけ計算
まず，粒子の位置から近傍のセル座標 i, j  を求める．次に，図 2-12 に示す
dr1, dr 2, dz1, dz2 を求める．その値を用いて，図中にあるような長方形 S1～S 4 の面積を求める．この面積比から粒子が感じる電磁場を内挿することができる．
具体的には，座標 r, z  上にある粒子が感じる電場 E r , z  は，計算セルの体積を S ，
図に示すような粒子によって仕切られた空間体積をそれぞれ S1～S 4 とすると，格子点上の重みはAn SnSn  1 ～ 4

35

(2.91)

となり，

E r, z   A1E i, j   A2 E i  1, j   A3 E i  1, j  1  A4 E i, j  1(2.92)と表わされる．

図 2-12 PIC 法におけるセルの構成

2.6.3.1.

電荷密度計算

電荷密度は図 2-12 に示されるセルの格子点上で定義され，セル中に N 個の粒子を代表した 1 個の超粒子が存在したとき，格子点 i, j  上に割り振られる電荷密度は，
 i, j  

S1 qN
S V

(2.93)となる．すべての粒子について和をとれば，格子点状での電荷密度が得られる．
2.6.3.2. 電子質量の変更(2-19)
今回の計算において，電子の速度成分は最大で 105 m/s 程度，イオンにおいては最大 103 m/s 程度であり，電子の速度を基準に時間刻みを設定するとイオンはほとんど動かない．そこで，電子の質量を実際のものよりα倍大きくした．質量の変更
により，電子の速度は√𝛼分の 1 に低下，タイムステップを√𝛼 倍に伸ばすことが出来る．イオンとの速度差は小さくなり計算結果の収束を早めることができる．
me _ virtual  me _ real
(2.94)
質量を変更後の値を仮想値，変更前を現実値と呼ぶ．電子の質量変更によって各物
理現象に矛盾を引き起こさないためにも各物理値を仮想値に修正する必要がある．
考慮すべきはラーマー運動への影響から磁場，衝突計算への影響から衝突断面積を36

修正する必要がある．以下の表に現実値と仮想値をしめす．
表 2-2 質量比の変更による補正係数

現実値

仮想値

𝑚𝑒

𝛼 𝑚𝑒

𝑣𝑚

𝑣𝑚
√𝛼

2.6.3.3.

Δt

√𝛼 Δt

B

√𝛼 𝐵

𝜎

√𝛼 𝜎

粒子の境界条件

放電室内のプラズマが実際に壁面に衝突した場合，放電室壁面は導体であるので，多くの場合，電子は電流として導体内に流れ，イオンは導体から電子を奪い中性粒子に戻る．そこで，今回の計算では電子，イオン共に壁面に衝突した粒子は消滅し，衝突後の挙動は追わないものとした．
2.6.3.4. Digital Filter
PIC 法で格子点上の量を決める場合, 例えば粒子の位置を格子上に割り振る場
合には, メッシュ内に存在する粒子の数が少ない領域ではその付近の格子に比べて数密度が小さくなり, 数密度の変化が不連続になる可能性がある. このような現象を避けるために Digital Filter と呼ばれる一種の補間を施すことがある. Digital
Filter は, 格子上に定義される値, 例えば電場, 数密度等を計算する際に発生する
高周波の雑音を除去する目的にも用いられる. (2-26)Digital Filter の式は以下のように表される.
newΦi,j,k←side
cornerMΦi,j,k + SΦi,j,k+ KΦi,j,k

M + 6S + 12K(2.95)

new
Φi,j,k を格子点上に定義される値とし、Φi,j,k
は Digital Filter を通した新しいsidecornerΦi,j,k を表す。またΦi,j,kとΦi,j,kは side term と corner term であり、それぞれ以下の式で

表される。

37

side
Φi,j,k= Φi−1,j,k + Φi+1,j,k + Φi,j−1,k + Φi,j+1,k+ Φi,j,k−1 + Φi,j,k+1

cornerΦi,j,k= Φi,j−1,k−1 + Φi,j+1,k−1 + Φi,j+1,k+1 + Φi,j−1,k+1 + Φi−1,j,k−1+ Φi−1,j,k+1 + Φi+1,j,k+1 + Φi+1,j,k−1 + Φi−1,j−1,k+ Φi−1,j+1,k + Φi+1,j+1,k + Φi+1,j−1,k(2.96)(2.97)

ここで、M、S、K の値をそれぞれ 12、4、1 とすると(3.44)式は次式となる。
1
1 side 1 cornernewΦi,j,k← Φi,j,k + Φi,j,k+ Φi,j,k4
12
48(2.98)

電子密度およびイオン密度，電位を計算した後に,それぞれの値にこの Digital
Filter を作用させる.

2.6.4. 中心軸近傍における座標変換(2-20)式(2.85)により，2.2 で述べた問題が生じる．そのため本コードでは，中心軸（𝑧軸）付近における座標変換を行った．この座標変換では，計算中に𝑟𝑐ℎ 以内に到達した粒子を対象に行っている．変換フローは図 2-13 に示すように，円筒座標系から直交座標系に変換し，2.4.1 で述べた．(2.1)式の運動方程式を解き，再度，円筒座標系に戻す流れとなっている．ここで，𝑟𝑐ℎ は，座標変換後の PIC 法を解きやすくする
ため，メッシュサイズとしている．
以下に，変換の手順について説明する．
先ず，位置と速度の円筒座標系から直交座標系の座標変換を行うと，

𝑥 𝑛 = 𝑟 𝑛 cos 𝜙 𝑛
𝑦 𝑛 = 𝑟 𝑛 sin 𝜙 𝑛𝑧𝑛 = 𝑧𝑛
𝑣𝑥𝑛 = 𝑣𝑟𝑛 cos 𝜙 𝑛 − 𝑟 𝑛 𝑣𝜙𝑛 sin 𝜙 𝑛𝑣𝑦𝑛 = 𝑣𝑟𝑛 sin 𝜙 𝑛 + 𝑟 𝑛 𝑣𝜙𝑛 cos 𝜙 𝑛
𝑣𝑧𝑛 = 𝑣𝑧𝑛

(2.99)(2.100)となる．ここでの n は繰り返し回数を示す．また，電場𝐸 (𝐸𝑟 , 𝐸𝑧 )と磁場𝐵 (𝐵𝑟 , 𝐵𝑧 )については，それぞれ以下のように変換する．

𝐸𝑥 = 𝐸𝑟 cos 𝜙
𝐸𝑦 = 𝐸𝑟 sin 𝜙
𝐸𝑧 = 𝐸𝑧
𝐵𝑥 = 𝐵𝑟 cos 𝜙
𝐵𝑦 = 𝐵𝑟 sin 𝜙
𝐵𝑧 = 𝐵𝑧(2.101)(2.102)

また，直交座標系における運動方程式は，2.4.1 の(2.1)式 より

𝑑𝑣𝑥 𝑞= (𝐸𝑥 + 𝑣𝑦 𝐵𝑧 − 𝑣𝑧 𝐵𝑦 )
𝑑𝑡
𝑚
𝑑𝑣𝑦 𝑞
= (𝐸𝑦 + 𝑣𝑧 𝐵𝑥 − 𝑣𝑥 𝐵𝑧 )𝑑𝑡𝑚
38

(2.103)

𝑑𝑣𝑧 𝑞= (𝐸𝑧 + 𝑣𝑥 𝐵𝑦 − 𝑣𝑦 𝐵𝑥 )𝑑𝑡𝑚となり，これを 4 次精度 Runge-Kutta 法にて解く．𝑞は電荷，𝑚は質量である．

∆𝑡(𝑏 + 2𝑏𝑥2 + 2𝑏𝑥3 + 𝑏𝑥4 )6 𝑥1∆𝑡𝑣𝑦𝑛+1 = 𝑣𝑦𝑛 + (𝑏𝑦1 + 2𝑏𝑦2 + 2𝑏𝑦3 + 𝑏𝑦4 )6∆𝑡𝑣𝑧𝑛+1 = 𝑣𝑧𝑛 + (𝑏𝑧1 + 2𝑏𝑧2 + 2𝑏𝑧3 + 𝑏𝑧4 )6𝑣𝑥𝑛+1 = 𝑣𝑥𝑛 +(2.104)

但し, Δt はタイムステップである．定数 b については，それぞれ以下を用いる．

𝑞(𝐸 + 𝑣𝑦 𝐵𝑧 − 𝑣𝑧 𝐵𝑦 )
𝑚 𝑥
∆𝑡𝑏𝑦1
𝑞
∆𝑡𝑏𝑧1
𝑏𝑥2 = {𝐸𝑥 + (𝑣𝑦 +)𝐵𝑧 − (𝑣𝑧 +
)𝐵𝑦 }
𝑚
2
2
∆𝑡𝑏𝑦2
𝑞
∆𝑡𝑏𝑧2
𝑏𝑥3 = {𝐸𝑥 + (𝑣𝑦 +)𝐵𝑧 − (𝑣𝑧 +
)𝐵𝑦 }
𝑚
2
2
𝑞
𝑏𝑥4 = {𝐸𝑥 + (𝑣𝑦 + ∆𝑡𝑏𝑦3 )𝐵𝑧 − (𝑣𝑧 + ∆𝑡𝑏𝑧3 )𝐵𝑦 }
𝑚
𝑏𝑥1 =

𝑞
(𝐸 + 𝑣𝑧 𝐵𝑥 − 𝑣𝑥 𝐵𝑧 )
𝑚 𝑦
𝑞
∆𝑡𝑏𝑧1
∆𝑡𝑏𝑥1
𝑏𝑦2 = {𝐸𝑦 + (𝑣𝑧 +)𝐵𝑥 − (𝑣𝑥 +
)𝐵𝑧 }
𝑚
2
2
𝑞
∆𝑡𝑏𝑧2
∆𝑡𝑏𝑥2
𝑏𝑦3 = {𝐸𝑦 + (𝑣𝑧 +)𝐵𝑥 − (𝑣𝑥 +
)𝐵𝑧 }
𝑚
2
2
𝑞
𝑏𝑦4 = {𝐸𝑦 + (𝑣𝑧 + ∆𝑡𝑏𝑧3 )𝐵𝑥 − (𝑣𝑥 + ∆𝑡𝑏𝑥3 )𝐵𝑧 }𝑚𝑞𝑏𝑧1 = (𝐸𝑧 + 𝑣𝑥 𝐵𝑦 − 𝑣𝑦 𝐵𝑥 )
𝑚
∆𝑡𝑏𝑦1
𝑞
∆𝑡𝑏𝑥1
𝑏𝑧2 = {𝐸𝑧 + (𝑣𝑥 +)𝐵𝑦 − (𝑣𝑦 +
)𝐵𝑥 }
𝑚
2
2
∆𝑡𝑏𝑦2
𝑞
∆𝑡𝑏𝑥2
𝑏𝑧3 = {𝐸𝑧 + (𝑣𝑥 +)𝐵𝑦 − (𝑣𝑦 +
)𝐵𝑥 }
𝑚
2
2
𝑞
𝑏𝑧4 = {𝐸𝑧 + (𝑣𝑥 + ∆𝑡𝑏𝑥3 )𝐵𝑦 − (𝑣𝑦 + ∆𝑡𝑏𝑦3 )𝐵𝑥 }𝑚(2.105)𝑏𝑦1 =(2.106)(2.107)

ここでの電場𝐸，磁場𝐵には，2.5.1.2 で説明した FDTD-PIC 同様に 3 次元直交座標における PIC 法による重み付けを行う．
運動方程式の解から，∆𝑡後の時間発展による位置を算出すると，

𝑥 𝑛+1 = 𝑥 𝑛 + 𝑣𝑥𝑛+1 ∆𝑡
𝑦 𝑛+1 = 𝑦 𝑛 + 𝑣𝑦𝑛+1 ∆𝑡𝑧 𝑛+1 = 𝑧 𝑛 + 𝑣𝑧𝑛+1 ∆𝑡
39(2.108)

となり，更新した位置，速度を以下のように直交座標系から円筒座標系に座標変換する．

𝑟 𝑛+1 = √(𝑥 𝑛+1 )2 + (𝑦 𝑛+1 )2
𝑦 𝑛+1
𝑛+1
−1
𝜙
= tan
𝑥 𝑛+1
𝑛+1
𝑛+1
𝑧
=𝑧

(2.109)𝑣𝑟𝑛+1 = 𝑣𝑥𝑛+1 cos 𝜙 𝑛+1 + 𝑣𝑦𝑛+1 sin 𝜙 𝑛+1
𝑣𝜙𝑛+1 = −𝑣𝑥𝑛+1 sin 𝜙 𝑛+1 + 𝑣𝑦𝑛+1 cos 𝜙 𝑛+1
𝑣𝑧𝑛+1 = 𝑣𝑧𝑛+1(2.110)2.6.5. 電場解析手法(2-19)
2.6.5.1.

基礎式

ポアソン方程式および電場と電位の関係式を基礎式として計算を行った．

Ñ 2j = -re0(2.111)E = -Ñj

(2.112)(2.111)式を差分化し,連立方程式として解いた．差分式を以下に示す．

aa1  i 1, j  aa2  i 1, j  aa3  i , j 1  aa4  i , j 1  aa5  i , j 
aa1  aa2 

aa3 

i, j
0

1
z 2

1
1

2
r
2ri , j r

(2.113)1
1
2
r
2ri , j r
æ 1
1 ö
aa5 = -2 ç 2 + 2 ÷è Dz Dr ø
aa4 

ここで， 𝑖 ，𝑗はそれぞれ𝑧，𝑟方向の格子点番号に対応している．𝛷𝑖,𝑗 は𝑖 ，𝑗における電位を表し，．𝜌𝑖,𝑗 はその点での電荷密度を表す．𝛥 𝑧，𝛥 𝑟は各方向の格子サイズであり，𝑟𝑖,𝑗 はその格子点の𝑟方向における位置を表す．
また，𝑟 = 0 上での係数は，

aa1  aa2 

aa3 = 0

40

1
z 2

(2.114)

4
Dr 2
æ 1
1 ö
aa5 = -3ç 2 + 2 ÷è Dz Dr ø
aa4 =とする．(2-21)(2.112)式より静中心差分すると，

E 1 =i+ , j
2

E

i, j+

1
2

=-

(j i+1, j - j i-1, j )2D(j i, j+1 - j i, j-1 )(2.115)2Dとなり，静電場を求める．
2.6.5.2. 反復法と PETSc(2.111)式はすべての格子点に１つずつ与えられるので，連立方程式として解くことができる．I 番目の方程式が以下のように表せる連立方程式を考える．
na x  b
j 1

i, j

j

i

(2.116)

𝑥は変数，a，b は定数を表す．この時，𝑥𝑖 以外の変数が固定されているという仮定のもと，次の式から𝑥𝑖 を求める．

 bi   ai , j x j 


j i


xi 
ai , j
これは，次の式によって定義される反復法を示している．
 bi   ai , j x j ( k 1) 


j i
(k )

xi  
ai , j(2.117)(2.118)

この手法は Jacobi 法と呼ばれる．注目すべき点は，方程式を解く順番を気にする必要が無いという点である．これは，各方程式を別々に解いているからであり，
Jacobi 法は各変数を一斉に置き換える方法として知られる．これに対し，各変数を一斉に置き換えるのではなく，方程式をひとつずつ解いていき，更新した変数を次の方程式の計算ですぐさま利用可能なものとするのがガウス-ザイデル法である．
ガウス-ザイデル法は一般的に Jacobi 法よりも高速であるとされる．
逐次的過剰緩和法(SOR 法)は，ガウス-ザイデル法に補外法を適用することによって発明された．この補外法は直前の反復解と各要素に対して逐次的に計算された
ガウス-ザイデル法の反復解との次式で表されるような加重平均の形をとる(ωは41

補外因子)．

xi(k ) xi( k )  1   xi( k 1)(2.119)

ωを反復解の収束率を加速させるような値に選択することにより，ガウス-ザイ
デル法より早い収束を得ることができる．ωには最適値が存在するが，最適値予め知ることは一般的には不可能である．そのため，しばしばω=2-O(h)という発見的な評価式が使用される．ここで h は元の物理領域を離散化するときのメッシュ間隔
である．
従来は，SOR 法にマルチグリッド法(2-23)を適用したポアソンソルバを使用してコードの高速化を図っていたが，並列化に不利である難点を抱えていた．粒子コードの並列化を適用し，粒子計算における計算コストの削減をしているが，ポアソン
ソルバの計算時間全体に占める割合が高く(全体の 7 割程度)，定常状態までの解析を不可能にしていた．このポアソンソルバを最適化するためには，数多くの反復法を試行するまたは，並列化に対応させる必要があった．しかしながら，学習コストと実装コストの観点から 1 つ 1 つ手法を試してくのは，時間コストに対する成果が見られないと考えた．そこで，数多くの反復解法と前処理と呼ばれる反復計算をする前に行列式を適切な形に加工する方法(SOR も前処理の一つである)の組み合わせを容易に試行できるライブラリの利用を決定した．
今回利用する反復解法のライブラリは PETSc(The Portable Extensible Toolkit
for Scientific Computation)である．アメリカの Argonne National Laboratory 内に発足した Mathematics and Computer Science Division project のメンバーである Satish Balay, William Gropp, Lois Curfman McInnes, Barry Smith の手によって開発されているプログラムライブラリであり現在アメリカでは航空力学やPETSc の特徴として以下に示す．

音響学等の数値解析に利用されている.

A) PETSc はオープンソースでありソースコードは, 以下の URL から入手できる.
http://www-fp.mcs.anl.gov/PETSc/
B) C, C++, FORTRAN の 3 種類のプログラミング言語をサポートしている.
（それぞれの PETSc 関数モジュール自体は C 言語でプログラミングされているが, その関数モジュールをつなげて作成するプログラムは上記 3 種類の言語のいずれでもよい. ）C) MPI（Message Passing Interface）を用いて, ベクトルの分散及び収集等といった並列化ベクトル計算や. 行列の組み合わせを簡単かつ効果的に行うことが可能である.
D)

疎行列のデータの格納をいくつかのフォーマットで行うことが可能である.

42

E) LU 前処理や ILU 前処理, block Jacobi 前処理などといった様々な前処理プ
ログラムを含んでいる. (前処理とは，前述のとおりはある問題を特定の解法で解くとき, 予めより効果的に解けるような同値な問題に変形することであり
Krylov 部分空間法と併用することによって, 収束性が著しく向上することが多く, 超大規模な問題の有効な解法として注目を浴びている.)F) GMRES, CG 法等といった様々な種類の Krylov 部分空間法（反復法）プロ
グラムを含んでいる.
G) 線形方程式・非線形方程式両方の solver を所有している.
利用方法としては，(2.113) (2.114)で示した係数行列と電荷密度分布を PETScの用意するデータ構造体に格納し，実行オプションとして，反復法と前処理を指定するものである．今回，選択した反復解法として，Bi-CGStab 法，前処理法として，Geometric-Algebraic Multigrid 法を選択した．尚，計算アルゴリズムについては，今回記載しない．数種の反復解法と前処理の組み合わせを現モデルで試行し，安定解が得られる組み合わせの中で，最も収束性の高いものを選択しただけである．
2.6.5.3.

電位の境界条件反復法の境界上兼として, 大別して 2 種類ある．Dirichlet 条件と Neumann 条件である．Dirichlet 条件とは，境界上の点の値を直接与えるものである．一方，Neumann 条件とは，計算領域は有限であるが，その外にも広がっていると考えるのが妥当である．この境界外との連続性を保つために，その境界上の点における微分を直接与えるものである．今回，壁面に対して，Dirichlet 条件として電位 0 としている．中和器は，イオンスラスタから排出されるイオンビームとのポテンシャ
ル差によって電子電流として引き出される．実験においては，中和器の排出口にコ
レクターと呼ばれるアノードをおき， これに数 10V の正電位をかけ，イオンビー
ムとのポテンシャル差を模擬している．この排出出口，後述の計算体系図では，
Plume に相当する部分に Dirichlet 条件として，𝛷𝑝𝑙𝑢𝑚𝑒 をかけることで，数値解析により，中和器性能を評価できるものとしている．

43

2.7. Full-PIC-MCC のフローチャート
図 2-13 に Full-PIC-MCC のフローチャートを示す．静磁場およびエネルギー獲
得確率は事前に計算しておく．初期化処理の際，粒子のデータは前回の情報を引き継ぎ，そのままリスタートさせる．ポアソンソルバにおける行列式を事前に作成しておき，各種 PETSc の初期化処理を行う．各格子点に電子およびイオンを重みつけにより，割り当て．電子密度，イオン密度を更新．電子密度，イオン密度より，
電荷密度を計算し，その情報を PETSc に認識できる形に加工して，PETSc により
ポアソン方程式を解く．収束後，電位分布より静電場を算出後，その結果を用いて，PIC 計算を行う．電子に関しては，運動方程式を解く前に，エネルギー獲得確率よりエネルギーを供給しておく．PIC 計算後，電子の衝突計算を行う．

44

図 2-13 Full-PIC-MCC のフローチャート

45

2.8.

弱連成

FDTD-PIC-MCC と Full−PIC-MCC の弱連成は，FDTD−PIC-MCC から Full−
PIC-MCC にはエネルギー獲得確率を Full−PIC-MCC から FDTD−PIC-MCC には
電子密度及び電子エネルギーである．弱連成のタイミングは，FDTD−PIC-MCCの計算時間については，マイクロ波周期によって決定し，約 50〜80 周期分を取れ
ば， 数値的なばらつきはないものとした． 本計算では，約 75 周期分の 30ns を解析している．Full−PIC-MCC の解析時間であるが，密度が大きく変わらない限
り Full−PIC-MCC から FDTD−PIC-MCC のフィードバックは意味がないとしている．よって，ある程度電子とイオンを動かす必要があり，この時間を経験的に 250ns程度とした．

2.9.

計算体系および静磁場分布

図 2-14 に本研究室所有のマイクロ波放電式中和器と解析領域を示す．これは，𝑦𝑧断面で切り取った断面図である．ここにおいて，𝑧方向を図の左から右方向，𝑦方向を図の下から上方向としており，𝑟𝑧断面図の際も𝑟方向が𝑦軸の向きと一致している．形状としては，円筒径の放電室(Wall)を箱状の永久磁石(サマリウムコバルト磁石)が𝑧軸方向に 3 個を 1 組として，棒磁石として 15 組の磁石が周方向に均等に配置されている．磁石を鉄で出来た FrontYork および BackYork が挟んでおり，この間に磁力線を形成する形状となっている．BackYork から伸びているがアンテナであり，これは L 字型を取っており，座標系において，𝑦軸方向にアンテナの先が向いている．FDTD−PIC-MCC では 3 次元直交座標を取り，この図を周方向に全てとった形状となっており，(𝑥, 𝑦, 𝑧 ) = (20,20,18) mm である．一方，Full-PIC-MCC
は，𝑟𝜙𝑧円筒座標系を取り，(𝑟, 𝜙, 𝑧) = (10,2𝜋,18) mm である．
この中和器内部の静磁場強度分布 𝑥 中心断面図を図 2-15 に示す．これは，Advanced Science Laboratory, Inc 社 の 3 次元磁場解析ソフト Amaze を用いて解析した磁場強度分布である． 今回使用するマイクロ波周波数は 2.45 GHz であるので，それに対応する ECR 磁場強度である 85 mT を赤線で示している．赤丸で囲んだところに，磁気チューブと呼ばれる電子補足層があり，この領域を電子
は，𝐸ｘ𝐵ドリフトまたは𝑔𝑟𝑎𝑑 𝐵ドリフトをしながら移動する.

46

図 2-14 本研究室所有のマイクロ波放電式中和器と解析領域

図 2-15 静磁場強度分布(赤線が ECR 層)

47

2.10.

各種物理条件と計算条件

2.10.1.物理条件
表 2-3 に本計算に用いた物理条件を示す． 本研究室の実験体系と同等のものを扱っている．DSMC により中性粒子密度分布を連成解析に用いる入力データとして用いている．推進剤であるキセノンガスの流量は，0.05 mg/sとしている．イオン
初期エネルギーは室温として，初期密度初期電子エネルギーに関しては，過去の実験データより設定した．(2-29)
表 2-3 物理条件

Xe
0.05 mg/s
2.45 GHz
8W
1.0×1017 m-3
2.0 eV
0.026 eV(300K)
40 V

推進剤
推進剤流量
マイクロ波周波数
投入電力
初期密度
電子初期エネルギー
イオン初期エネルギー
出口電位 𝛷𝑝𝑙𝑢𝑚𝑒2.10.2. 計算条件
FDTD−PIC-MCC
メッシュサイズを 0.5 mm とし，𝑥方向に 40 格子𝑦方向に 40 格子𝑧方向に 36 格子とした．FDTD-PIC-MCC でみる物理対象には，シースの解析は不要である．マ
イクロ波波長は約 7cm 程度であり，本メッシュサイズで問題はないとされる．メ
ッシュサイズ 0.5 mm より FDTD における計算制約条件は，プラズマ周波数より，−10.2 𝜔𝑝𝑒= 8.0 ×10-12 sec，電磁波の伝搬速度を光速と仮定すると Courant 条件より，9.6 ×10-13 sec 以下であればよく計算誤差を考え，2.5 ×10-13 sec とした．PIC
における制約条件としては，サイクロトロン周波数より，2.0 ×10-12 sec，衝突計
算の制約より，7.3×10-10 sec，最大電子速度を 50 eV とした時の 1 セルあたりを飛び越えないための Courant 条件は，6.9×10-11 sec となり，さらにその 10 分の 1程度とすると 7×10-12 sec となる．このなかで最も制約条件の厳しいのは，Courant 条件であるためこれを満たしさえすればよく．FDTD10 回につき PIC を
1 回解くとすると，2.5×10-12 sec で十分計算の制約を満たす．

48

表 2-4 計算条件(FDTD−PIC-MCC)

0.5 mm
2.5 ×10-13 sec
2.5 ×10-12 sec
9.85 ×104
400 万個130 個

メッシュサイズ
タイムステップ(FDTD)
タイムステップ(PIC)
超粒子の重み超粒子数
1 セルあたりの粒子数

また，FDTD の境界条件であるが，体系図の Plume 面及びアンテナの根本のみを
Mur の１次吸収境界をかけている．
Full-PIC-MCC
Full−PIC-MCC の解析対象は，シース解析であるために，メッシュサイズがデバイ長以下であるのは必須である．デバイ長は電子密度 2.0×1017 m-3 平均電子エ
ネルギー8 eV 程度とすると，0.047 mm 程度である．そのために計算精度を考え，0.025 mm とした．この際の Courant 条件は，最大電子エネルギーを 50 eV と仮定すると，4.0×10-12 sec となる．さらにその 10 分の 1 を考えると，4.0×10-13 sec
であるが，タイムステップを少しでも広く取るために，Courant の 8 分の 1 程度として，5.0×10-13 sec とした．これを仮想値に変換すると，2.5×10-11 sec つまり0.025 ns とした．
表 2-5 計算条件(Full−PIC-MCC)

0.025 mm
5.0 ×10-13 sec
2.5 ×10-11 sec
3.90 ×104
1000 万個
50 個

メッシュサイズ
タイムステップ(現実値)
タイムステップ(仮想値)
超粒子の重み超粒子数(イオン・電子)
1 セルあたりの粒子数49

第3章 結果と考察
本研究では，図 3-1 に示すように，アンテナを中心に領域を 8 つに分割して計算を行った．以下に結果を表す．

図 3-1 中和器モデル（8 分割）

3.1.

エネルギー獲得確率

図 3-2 に FDTD-PIC-MCC による 30 ns 後の電子のエネルギー獲得確率の結果を表す．図 3-2 の⓪で示す獲得確率が，図 3-1 のアンテナを含む領域⓪で周方向に平均をとった値である．一方，旧コードで示す獲得確率は，2 次元のFull-PIC-MCC で使用していた値であり，周方向 2𝝅に平均をとった値である．また，エネルギー獲得確率はセル毎に算出されるため，図 3-3 で示す r=3 mm，z=8
mm のセルにおけるエネルギー獲得確率をそれぞれ示している．各エネルギー獲得
確率の平均獲得エネルギーは，アンテナ領域⓪が 0.1 eV，旧コードは 0.3 eV である．

図 3-2 電子エネルギー獲得確率50図 3-3 エネルギー獲得確率の算出位置（r = 3 mm, z = 8 mm）

これらの平均獲得エネルギーを比較すると，周方向にエネルギー獲得確率が異なっていることがわかる．そこで，領域毎に図 3-3 に示す位置のエネルギー獲得確率を算出し，平均獲得エネルギーをプロットした結果を図 3-4 に表す．

図 3-4 平均獲得エネルギー（30 ns , 8 分割）

図 3-4 より，やはり周方向に獲得確率が異なっており，アンテナ周りの領域ではエネルギーが低いことがわかる[1]．さらに，アンテナから反時計回りの領域②，③，④でエネルギーが高く，アンテナを中心として左右に分けて観ても非対称であることが確認できる[2]．[1]については，アンテナの強電場によりエネルギー損失が起こりエネルギー獲得確率の平均値が低くなるためだと考えられる．1.2.1 で述べたように，マイクロ波放電式中和器はアンテナによる交流電場で電子にエネルギーを与えプラズマを生成する．そのため，正の電場が掛かると電子がアンテナに引51

き寄せられ衝突し，エネルギー損失となる．本コードのエネルギー獲得確率の計算においては，アンテナに衝突した粒子は消失と捉え，負のエネルギーとしてカウントする．ここで，図 3-5 はアンテナ領域⓪と領域③のエネルギー獲得確率を表している．各平均獲得エネルギーは，領域⓪が 0.1 eV，領域③が 0.9 eV である．図 3-5
からも，アンテナ領域⓪は領域③と比べて負のエネルギー獲得確率が高くなってい
ることがわかる．以上から，アンテナ周辺の領域では平均獲得エネルギーが低くなると考える．

図 3-5 平均獲得エネルギー（領域⓪・領域③）[2]については，𝑔𝑟𝑎𝑑 𝐵ドリフトによってエネルギーが伝搬したためだと考えられる．本研究で扱う中和器モデルにおいて図 3-3 に示す位置では，電子は反時計方向にドリフトするため，30 ns 後の平均獲得エネルギーにおいて左右で非対称性が生じると考えられる．ここで，図 3-6 に 40 ns 後における各領域の平均獲得エネルギーを表す．図 3-4 の 30 ns 後の平均獲得エネルギーにおいては領域②，③，④で高いエネルギーが観られたが，図 3-6 の 40 ns 後においては，領域②，③，④，⑤で高いエネルギーが観られる．このことから時間の経過とともにドリフトに
より反時計方向にエネルギーの伝播が行われていることがわかる．さらに，30 ns
後の最大の平均獲得エネルギーは 0.9 eV であったが 40 ns 後においては，最も高い平均獲得エネルギーでも 0.4 eV と小さく，図 3-4 と図 3-6 の同じ領域を比べても 40 ns 後は平均獲得エネルギーが小さい値となることがわかる．これは，反時計方向に隣接する領域がアンテナ領域である領域⑦の 30 ns と 40 ns 後の平均獲得
エネルギーを比べると，30 ns 後が-0.2 eV，40 ns 後が-0.3 eV で 40 ns 後の方の損失が大きくなっていることから，30 ns と 40 ns では時間の経過とともに反時計方
向にドリフトによってエネルギーの伝播が行われるが，アンテナの損失も大きくなるため，40 ns 後の方の平均獲得エネルギーが小さくなると考えられる．

52

0
1

2

平均獲得エネルギー [eV]

1.0
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4

7

6

3

5
4

図 3-6 平均獲得エネルギー

3.2.

(40 ns , 8 分割)

各物理値分布

領域毎のエネルギー獲得確率を用いて Full-PIC-MCC により計算を行なった結果を示す．

3.2.1. 電子エネルギー分布
図 3-7 と図 3-8 に Full-PIC-MCC によるアンテナ領域⓪と領域③の電子エネルギー分布の結果を表す．r=3 mm，z=8 mm のセルにおける平均獲得エネルギーは
領域⓪が 0.1 eV，領域③が 0.9 eV である．図 3-7 と図 3-8 より，領域③の方が高エネルギーの電子が広範囲に分布していることがわかる．これは各領域でアンテナから得るエネルギーが異なるため異なる分布になったと考えられる．

図 3-7 電子エネルギー分布（領域⓪ 平均獲得エネルギー 0.1 eV）

53

図 3-8 電子エネルギー分布（領域③ 平均獲得エネルギー 0.9 eV）

3.2.2. 電子密度分布
図 3-9 と図 3-10 に Full-PIC-MCC によるアンテナ領域⓪と領域③の電子密度分布を示す．r=3 mm，z=8 mm のセルにおける平均獲得エネルギーは領域⓪が 0.1
eV，領域③が 0.9 eV である．図 3-9 と図 3-10 を比較すると，密度分布が周方向𝜙に不均一になることがわかる．これは 3.2.1 で述べたように，各領域でアンテナから得るエネルギーが異なることで電子エネルギー分布も各領域で不均一になるため密度分布に差が生じると言える．

図 3-9 電子密度分布（領域⓪ 平均獲得エネルギー 0.1 eV）

54

図 3-10 電子密度分布（領域③ 平均獲得エネルギー 0.9 eV）

3.3.

過剰電離による計算の不安定性

本 Full-PIC-MCC コードの計算において計算の不安定性が観られた．3.2.1 で示した電子エネルギー分布の結果では，25 ns 時で 30 eV を超える高エネルギー分布が観られたが，文献(2-30)と比較すると，電子エネルギーは数 eV であるのに対し，
本コードでは非常に高い値となっている．これは，3.1 で述べたように各領域でエネルギー獲得確率を算出し Full-PIC-MCC で領域分割を行なって計算することで，領域③のように，領域によっては旧コードで用いていたエネルギー獲得確率より大きいエネルギーを Full-PIC-MCC で与えることが一因として挙げられる．これによって過剰電離が起こり，計算が不安定になると考えている．そこで 2.8 で
Full-PIC-MCC の計算時間について述べたが，各領域でエネルギーが異なり粒子が十分移動する時間も異なるため，計算時間の再検討が必要である．また，
FDTD-PIC-MCC の計算時間について，過剰なエネルギー供給によって過剰電離が起こることも一因として考えられる．3.1 で述べたように計算時間によって電子の獲得エネルギー値は変化するため，計算時間の再検討を行い定常状態までの計算を目指す．
次に，アンテナを中心に領域を 2 つに分割して弱連成解析によって算出されたエネルギー獲得確率の結果を表す．図 3-11 は連成 1 回目の FDTD-PIC-MCC による
図 3-3 に示す位置の電子エネルギー獲得確率を表し，図 3-12 は連成 2 回目のFDTD-PIC-MCC の図 3-3 に示す位置の電子エネルギー獲得確率を表している．
平均獲得エネルギーは 1 回目のアンテナを含む領域が 0.2 eV，アンテナを含まない領域が 0.4 eV であり，2 回目ではアンテナを含む領域が 0.4 eV，アンテナを含まない領域が 0.2 eV であった．1 回目の Full-PIC-MCC の計算結果を弱連成によ55

り電子密度を更新した 2 回目の FDTD-PIC-MCC の電子密度分布を図 3-13 に表す．

図 3-11 電子エネルギー獲得確率（2 分割 連成 1 回目）

図 3-12 電子エネルギー獲得確率（2 分割 連成 2 回目）

図 3-13 FDTD-PIC-MCC 電子密度分布（練成後）56

3.2.2 で述べたように，周方向に密度分布が異なることがわかった．図 3-13 の領域を 2 つに分割した結果も粒子がアンテナに消失することによりアンテナを含む領域のアンテナ周りでは低い密度分布となっていることがわかる．本
Full-PIC-MCC の 3 次元の密度を更新することで，図 3-11 と図 3-12 に表すように FDTD-PIC-MCC の電子エネルギー獲得確率が変化することが確認できた．これは周方向に密度が不均一になることから，密度更新後の FDTD-PIC-MCC では，密度の低い領域でエネルギーを多く受けやすくなるためだと考えられる．
これ以降の解析は，過剰電離による不安定性の改善が求められるが，
Full-PIC-MCC の 3 次元化を行うことで周方向によるプラズマの不均一性や弱練成によるエネルギー獲得確率の変化が確認できた．

57

第4章 結論
本研究ではマイクロ波放電式中和器のアンテナ形状の効果を取り入れた解析ツールの開発を行った．
アンテナ周辺ではアンテナの強電場によりエネルギーの損失が大きく，各領域でエネルギー獲得確率が異なることがわかった．そこで，MPI によって𝜙方向の領域分割を行い Full-PIC-MCC を 3 次元(𝑟, 𝜙, 𝑧)に拡張しプラズマ解析を行なった．計算
結果では，異なるエネルギー獲得確率を用いることで𝜙方向に不均一性が確認でき，
このことからアンテナの再設計には 3 次元の解析ツールが必要であることがわかった．
しかし，本解析コードは計算の不安定性があるため本コードの計算時間の見直しが必要である．これを改善することで定常状態まで解くことを可能にし，3 次元化による解析コードの妥当性を実験値との比較から検証する必要がある．

58

参考文献
(2-1) 楠田将士，“マイクロ波放電式中和器の引き出し電流電圧特性解析”，平成
26 年度九州大学修士論文(2-2) 奥田雄，”中性粒子密度分布考慮したマイクロ波放電式中和器内部物理現象
解析コードの開発”，平成 27 年度九州大学修士論文(2-3) 内藤裕志, プラズマ・核融合学会誌, 74 (1998) ,470.
(2-4) 田中實, 山本良一：“計算物理学と計算化学―分子動力学法とモンテカルロ法” (海文堂,1998)(2-5) 石黒静児, プラズマ・核融合学会誌, 74 (1998) , 591(2-6)
Giovanni Lapenta, Jeremiah U. Brackbill, Journal of Computational
Physics, 115, (1994), 213.
(2-7) 高村秀一：“プラズマ理工学入門”(森北出版, 1997)(2-8) NIFS DATABASE 核融合科学研究所データベース(http://dbshino.nifs.ac.jp/index-j.html)(2-9) Mersenne Twister Home Page :(http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html)(2-10) 南部健一：“コンピュータアナリシスシリーズ 7 原子･分子モデルを用いる数値シミュレーション 日本機会学会編” (コロナ社, 1996)(2-11) V. Vahedi and M. Surendra, Computer Physics Communications, 87(1995), 179.
(2-12) 電気学会・マイクロ波プラズマ調査専門委員会：“マイクロ波プラズマの技術” (オーム社, 2003)(2-13) K. Koura, Phys. Fluids, 29(1986), 11.
(2-14) 橋本修, 安部琢美：“FDTD 時間領域差分法入門” (森北出版, 1996)(2-15) 宇野亨：“FDTD 法による電磁界およびアンテナ解析” (コロナ社, 1998)(2-16) K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.
(2-17) G. Mur, IEEE Trans. Electromagnetic Compat., 40 (1998), 100.
(2-18) (https://www.mpich.org/static/docs/v3.2/www3/MPI_Comm_split.html)(2-19) Szabo Jr, James Joseph. “Fully kinetic numerical modeling of a plasma
thruster.” Diss. Massachusetts Institute of Technology, 2001.
(2-20) 吉見直文, “磁場を利用した ICF システムの二次元電磁粒子コードによる解析”，平成 10 年度九州大学修士論文．
(2-21) 田中守也，”電子・イオンビーム工学”(電気学会，1995)(2-22) Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June
Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo, Charies Romine,
Henk va der Vorst, ：“反復法 Templates” (朝倉書店, 1996)
59

(2-23) 廣池匠哉, “マルチグリッド法を用いたマイクロ波放電式中和器のプラズ
マ損失機構解明”，平成 24 年度九州大学修士論文．
(2-24) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical
Physics, 55 (1971), 4100.
(2-25) 保原充，大宮司久明：”数値流体力学:基礎と応用”(東京大学出版会，1992)(2-26) Birdsall, C. K., and A. B. Langdon. "Plasma physics via computer
simulation." (1985).
(2-27) G. M. Tuner, I.S. Falconer, B. W. James, and D. R. McKenzie, J. Appl.
Phys. B65B, (1989)(2-28) 渡邊裕樹, “マイクロ波放電式電子源内の電子加熱解析”.STEP-2013-014(2-29) )小野寺範義, 竹ヶ原春貴, 西山和孝, 船木一幸, 國中均: マイクロ波放電型中和器の電子放出機構, 日本航空宇 宙学会論文集, 第 49 巻, 第 564 号,
pp.27-31, 2001.
(2-30) 窪田健一，渡邊裕樹，山本直嗣，中島秀紀 ：“マイクロ波放電式中和器の電子引出に関する 3 次元 Hybrid-PIC 解析”（日本航空宇宙学会論文集，2015）60

謝辞
本項を作成するにあたり，ご指導を頂き，自ら考え学ぶことの大切さを教えてくださった山本直嗣教授に深謝いたします．私に丁寧かつ的確な助言をくださった森田太智助教に感謝の意を表します．中島秀紀名誉教授には日頃から私の進捗具合を気にかけて下さり，貴重なアドバイスを頂きました．心
より感謝致します．研究活動に関する事務手続きや学生生活の様々なサポートをして下さった馬渡隆子氏，大内田映子氏に感謝致します．
研究のみならず公私ともにお世話してくださった齋藤直哉氏，江川雄亮氏に感謝致します．研究と共通の趣味について議論を重ね合った牛尾康一氏に感謝いたします．研究で行き詰まった時に助けて下さった枝本雅史氏に感謝します．研究室で充実した時間を共に過ごした，卒業された先輩方，M1，B4のみなさんに感謝します．
研究や就職活動，2 年間苦楽を共に過ごした板谷佑太朗氏，岩本政隆氏，川原友太郎氏，福本湧氏，淵上太貴氏に心から感謝します．
最後になりますが，ここまで温かい目で見守ってくれた家族に心から感謝します．

61