卒

題

業

論

文

目

小型プラズマ源における作動ガスが性能へ及ぼす影
響

氏

名

松尾

拓哉

所属教育分野プラズマ理工学

指 導 教 員

中島

秀紀

教授

山本

直嗣

准教授

九州大学工学部エネルギー科学科
提出年月

平成２７年２月

目次
第1章

序論 .............................................................................................................. 1

1.1

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

1.2

電気推進 .............................................................................................................. 1

1.3

小型マイクロ波プラズマ ..................................................................................... 2

1.3.1

マイクロ波放電............................................................................................. 2

1.3.2

ECR プラズマ ............................................................................................... 2

1.4

磁気チューブ....................................................................................................... 4

1.5

小型マイクロ波プラズマ推進器の原理 ............................................................... 6

1.6

研究目的 .............................................................................................................. 8第２章

数値解析................................................................................................. 102.1

粒子計算(2-1) ............................................................................................. 10

2.1.1

基礎式 .................................................................................................... 10

2.1.2

メッシュサイズ・タイムステップ(2-3) .................................................... 11

2.1.3 PIC 法(2-4) ............................................................................................... 122.1.4

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

荷電粒子と中性粒子の衝突.................................................................... 14

2.1.6

平均自由行程と透過度(2-5) ...................................................................... 142.1.7

衝突確率・衝突判定 .............................................................................. 15

2.1.8

衝突計算 ................................................................................................ 15

2.1.9

エネルギーと散乱角 .............................................................................. 17

2.1.10 Null-collision 法(2-11)(2-12)(2-13) ................................................................ 192.2

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

2.2.1

基礎式 .................................................................................................... 22

2.2.2

差分式 .................................................................................................... 23

2.2.3

アンテナ解析 ......................................................................................... 24

2.2.4

安定条件 ................................................................................................ 27

2.2.5 PEC(Perfect Electric Conductor )境界条件 .......................................... 27
2.2.6 Mur の吸収境界条件(2-17) ........................................................................ 27
2.2.7
2.3

ポインティングベクトル(2-18) ................................................................. 28

カップリングコード ................................................................................. 29

2.3.1

計算体系 ................................................................................................ 31

2.3.2

計算格子 ................................................................................................ 31

2.3.3

時間刻み .......................... エラー! ブックマークが定義されていません。

2.3.4

電子初期条件......................................................................................... 31

2.3.5

粒子数 ................................................................................................... 32

2.3.6

推進剤 ................................................................................................... 32

2.3.7

磁場強度 ................................................................................................ 32

2.3.8

計算コードの並列化(2-24)........................................................................ 32

第 3 章 計算結果と考察 .................................................................................. 36
3.1

アルゴンとキセノンの電子の平均エネルギー .................................................. 36

3.2

アルゴンとキセノンのイオンの電離数 ............................................................. 39

第 4 章 結論 ..................................................................................................... 40
謝辞

第1章

序論

1.1 背景昨年,宇宙飛行士の若田光一氏が宇宙国際ステーションに長期滞在したことが話題になったが,そのミッションのひとつの超小型衛星の放出が見事に成功し超小型衛星が大きく取り上げられた.(1-1)今回の超小型衛星の放出の为な目的は宇宙より地球を撮影することだったが,小型衛星は天気予報やカーナビゲーションシステムなど私たちの生活に欠かせないものとなっている.小型衛星とは,厳格な定義はないが,多くの場合５００ｋｇ
より軽い衛星のことを指す.(1-2)小型衛星の最大の特徴は大型・中型衛星に比べ,製作コス
トが低く,製作時間が短いことである.(1-3)従来までの大型・中型衛星では製作コストが数百億円かかり,製作期間も長かったので民間企業の宇宙開発への参加を妨げ,失敗した時のリスクが計り知れなかった.しかし,小型衛星の開発が進んだことにより,宇宙開発への参入リスクが大幅に低減され,より多くの民間企業の新規参入や宇宙開発のビジネス化など宇宙事業が活性化している.また学生においても,九州工業大学が宇宙航空研究開発
機構の人工衛星「アストロ H」との相乗りで開発中の超小型衛星「風龍四号」の打ち上げが決まった.(1-4)このように小型衛星の開発による製作費の低コスト化と製作期間の短
縮が宇宙産業に競争力を生み,さらなる宇宙の開発につながると思われる.
人工衛星や宇宙ステーションなどの地球周回物は，地球が楕円球であることや大気抵
抗の存在，また太陽や月の引力による摂動があるため，その姿勢制御・軌道制御のために推進機を搭載している．現在，衛星の寿命を決定するものは機器の耐年数よりもむしろ推進剤の枯渇である．特に小型衛星においてはサイズの制約上，搭載可能な推進剤の量も制限されてしまう．従って，尐量の推進剤で長期間運用可能な燃費の良い推進機の開発が重要となる．そして，その要求に応えるものとして利用されているのが電気推進である．

1.2 電気推進(1-5)
宇宙機の推進方式としては化学推進と電気推進に大別される.化学推進とはロケットエンジンなどにみられるように推進剤を燃焼室で燃やし,高温高圧のガスを排出することによって推力を得るものである.特徴は,高い推力が得られ,比推力（単位推進剤重量によって 9.8 N の推力を持続させることができる時間,すなわち燃費）は小さい点である.
比推力が小さいと大量の推進剤を搭載せねばならなくなるため,小型衛星の推進方式としては適していない.
一方,電気推進とは,太陽光エネルギーや原子力エネルギーを一旦電気エネルギーに変換し,推進剤を加熱・電離させ,電磁気力を用いて推進剤を排出し,その反作用によって推力を得る推進方式である.電気推進の特徴としては,ロケットエンジンなどの化学推進に比べて推力は小さいものの,比推力が約 10 倍である点であり,搭載しなければならない燃料を大幅に減らすことができる点である.代表的な電気推進機には,イオンスラスタ,アー1

クジェットスラスタ,MPD（Magneto-Plasma-Dynamic）スラスタ及びホールスラスタなどがあり,ミッションよって使い分けられている.
故に,本研究では小型衛星に適用可能な小型かつ小電力な電気推進の開発を目指す.研究で用いる推進器としては 2012 年度より同研究室の豊田により研究が始められた新しいタイプの推進器を用いる.（以下,本研究で用いる推進器を小型マイクロ波プラズマ推進器,もしくは本推進器と呼ぶ）この推進器の最大の特徴は小型かつ小電力という目標を達成することが理論的に可能な点である.これらの特徴はこの推進器において必要な構成要素がマイクロ波源と推進剤タンクと流量制御系(流量制御器,減圧器など)のみである
ことによる.
他の電気推進器においては,プラズマ生成部のほかにビーム中和部が必要であり,そのための電力やスペースが必要となるため目標達成にあまり向いていない.例えば本研究室でも小型化の研究が行われてきたイオンスラスタは本推進器と同様マイクロ波でプ
ラズマを生成したのちに,グリッドに電圧をかけることによってイオンを加速するイオ
ン加速部とビームを中和する中和部が必要となる.また,現在盛んに研究が行われている
ホールスラスタはプラズマ生成部と加速部が１つになっていることや,推力密度が高いという特徴を持っており小型化に向いているが,プラズマ生成とビームの中和のための電子を放出するホローカソードと呼ばれる部分が必要となる(1-6).しかし,部品や要素が多くなると故障のリスクが増える等のデメリットが発生する.よって,もしイオンスラスタ
やホールスラスタで同じ目標を達成できたとしても小型マイクロ波プラズマ推進器を用いるほうがよいと言える.

1.3 小型マイクロ波プラズマ
1.3.1 マイクロ波放電
マイクロ波放電型とはマイクロ波帯域の交流電界で電子を加速し,この電子が中性粒
子と電離衝突してプラズマを生成,放電を維持するというものである.マイクロ波放電方
式は「はやぶさ」にも用いられ,高い信頼を得ている(1-7)(1-8).このマイクロ波放電を推進器に採用することで以下のようなメリットが得られる.
① 放電用電極を必要としないのでスラスタの長寿命化と構造の簡略化が可能② 予備加熱が不要となるので,スラスタの迅速なスタートが可能
以上のような理由により本研究ではマイクロ波放電方式を採用している.

1.3.2

ECR プラズマ(1-9)

ECRの原理と概念図を図1.1と図1.2 にそれぞれ示す.真空中に存在する荷電粒子は磁場中でローレンツ力を受ける.この向心力のため電子は磁力線に巻きつくようなサイク
ロトロン運動と呼ばれる回転運動を行う.磁界における電子の運動方程式は次式で表さ
れる.
2

ｍ

d𝒗
= 𝑒𝒗 × 𝑩
d𝑡(2-1)

ここでｍ は電子の質量、 𝒗 は速度ベクトル、𝑒 は電荷量、𝑩 は磁束密度である。このときの円軌道の半径𝑟𝐿 はラーマー半径と呼ばれ、
𝑚𝒗⊥
𝑟𝐿 =
(2-2)
|𝑒|𝐵
で与えられる。ここで𝒗⊥ は電子の𝑩 に対する垂直な速度成分の大きさである。
プラズマを構成している荷電粒子は全て反磁性体である。そのためサイクロトロン運動の回転の方向は、外部磁場の向きに対して荷電粒子の回転によってできる磁場が常に外部磁場と逆を向く方向である。つまり、電子は磁場に対して右回りの回転運動を行う。
この回転運動の角周波数はサイクロトロン角周波数と呼ばれ
|𝑒|𝑩
(2-3)
𝑚
で与えられる。電子は電場と逆向きに加速されるため、磁場中の電子の回転方向と逆向
きに、回転する速さが等しい電場をかけると電子は連続的に加速され、電場から効率的にエネルギーを受けとることができる。これが電子サイクロトロン共鳴現象である。
ECRプラズマ発生装置では、導波管を通して角周波数𝜔 のマイクロ波を入射し、プラズマを生成する。生成されたプラズマ中の電子は、サイクロトロン角周波数に従って磁力線方向に向かって右回りの旋回運動を行う。一方、強磁場側から入射したマイクロ波はプラズマ中を浸透し、電子サイクロトロン波と呼ばれる右回りの円偏波を励起する。この波は弱磁場側へ伝播し角周波数が電子サイクロトロン角周波数と一致する層で急速に減尐してマイクロ波から電子にエネルギーが吸収される。つまり、𝜔𝑐 = 𝜔 の𝜔𝑐 =関係が成り立つときにECRが生じる。
また逆に、ECRを起こすときに必要な磁束密度の大きさを𝐵
𝑚𝜔
𝐵 =
| |

とすれば (2-3)式より(2-4)となる。ECRによってエネルギーを増大させた電子は，磁力線に拘束されながら次々と効率よく周辺の中性粒子と電離衝突を繰り返す。このようにして生成されたプラズマがECRプラズマである。

3

1.4

図1.1

ECRの原理図1.2

ECRの概念図

磁気チューブ

ECR プラズマにおいて、電子は電場によって中性粒子を電離させるのに十分なエネル
ギーを持つまで加速されなければならない。そこで、本推進器では電子を十分加速する
ために、電子が往復運動を行う磁気チューブと呼ばれる磁束管を利用し、電子の動きを制限している。以下、磁気チューブの原理について、詳しく説明する。
磁場に垂直な方向の荷電粒子の運動は、磁力線の回りをまわるような軌道を描くこと
は 1.3 節で述べた。ここでは、プラズマ中に図 1.3 に示すような磁力線が存在する状況を考える。この図の点 O に旋回の中心を持ち、速度 ν0 (磁場に垂直成分ν⊥ 、平行成分ν// )を持つ粒子があるとする。荷電粒子の運動エネルギーの保存より
4

v 2  v||2  v 22
 v||,20  v ,0

v

(2-5)

2
0

となる。また、磁場が緩やかに変化する場合、磁気モーメント   (1/2)mv2 /B が一定に保たれることを考慮すると、磁場強度が増加するに従い、速度の平行成分が減尐してい
き、ある磁場強度で荷電粒子は磁場の弱い方向へ反射される。このような磁場による荷電粒子の反射を磁気ミラーと呼ぶ。μ の一定性より、 

1 2
1 2
mv 0
mv 
2
2

B0
B

(2-6)

である。また、v  v sin  v 0 sin0であることを考慮し、式（2-7）を式（2-6）に代入すると、
sin2  0 sin 2 


B0
B

(2-7)

(2-8)

となる。よって、反射点での磁場強度は、 v||  0 と   90 であるのでB
Bref  20(2-9)

sin  0

となる。



図 1.3 のように両端に磁気ミラーを持ち、荷電粒子が往復する領域を磁気チューブと呼ぶ。小型マイクロ波プラズマ推進器では電子がこの磁気チューブを往復しアンテナ近傍を通過するたびにエネルギーを受け取り、加速される。

v0
B

v0

O

図 1.3

磁気ミラー（チューブ）による粒子の閉じ込め

5

1.5

小型マイクロ波プラズマ推進器の原理

本研究で扱う推進器について、推進器の概略図を図 1.4 に示す.放電室外部には Sm-Co
磁石が配置されており、磁石から出た磁力線が軟鉄製のヨークを通って放電室内に導かれている。このときセントラルヨークとフロントヨーク間で磁気チューブが形成されている。そして、この磁気チューブを往復する電子をマイクロ波によって加速し、その電子に中性粒子が衝突、電離することによってプラズマが生成される。このとき生成されたイオンはプラズマ電位と宇宙空間電位(即ち 0 V)の電位差により加速、排出される。また、電子については磁気チューブ内で高いエネルギーを持ったものが電位差を乗り越え、
イオンと共に排出される。これにより、推進器から排出されるビームは中性に保たれる。
このとき、排出される電子が高いエネルギーのもののみだと、イオンとのバランスがとれないように思われるかもしれないが、電子は速度が速く順次補充されるため問題はない。また、このバランスは自然かつ自動的に行われる。本推進器において推力が発生していることは同研究グループの豊田によって行われたイオンエネルギーアナライザー
を用いて行われた実験から確かめられている。

6

図 1.4推進器の概略図と軸上電位分布71.6 研究目的
本研究室ではマイクロ波小型推進器の磁場強度・アンテナ形状・放電室形状などを変化させて,実験を行い，中和器の性能評価及び内部プラズマの計測を実施してきた．000
しかし,放電室内部の物理メカニズムがどのように変化し，どのような理由で性能に影響をおよぼしているのかは解明されていない．その原因はエンジン自体が小型であるために静電プローブ等の装置を入れることで放電室内部のプラズマに擾乱を引き起こし，
正確な測定ができないからである． レーザ計測等を用いれば計測が可能であるが，放電室にレーザ掃引用の光路を確保しなければならず，一つのスラスタで複数の位置で計測することは難しい．そこで数値解析は非常に有効な手段である.近年の計算機の発展もあり実験では知り得る事のできない情報を入手できるようになったそれにより，実験による時間的・金銭的コストを抑えた研究開発ができるようになってきている．本研究室では，下川床、楠田らによっての内部物理の解明および，中和器の性能を評価する数
値解析コードの開発を行なってきた 00 エラー! 参照元が見つかりません。エラー! 参照元が見つかりません。．
これまでの課題として，これまで電子のエネルギー授受は ECR 層で効率的に行われていると考えられてきたが，実験結果より，ECR 層以外でもエネルギーのやり取りが行われている可能性が示唆された．そこで電子の平均エネルギーの分布を調べることによりエネルギーが授受されている場所を特定する.
また,イオンエンジンでは推進剤としてキセノンが为流である.その理由としてはイオン化エネルギーが低く,推進効率が高いことが挙げられる(2).しかし,本研究ではより安価で比推力のあるアルゴンを用いている.そこでプラズマ源における作動ガスであるキセ
ノンとアルゴンの性能に及ぼす影響の違いを数値解析することでアルゴンの推進剤としての適正を確認する.

8

参考文献
(1-1)日本経済新聞(Web 版)“ 宇宙ステーションから衛星放出日本の技術米社活用”、2014/01/27(1-2) 飯田尚志：小型衛星の現状と将来動向、電子情報通信学会技術研究報告. SANE, 宇宙・航行エレクトロニクス 93(398), 49-54, 1993-12-17(1-3) 趙孟佑、増井博一、入江大樹、今里昂史：超小型衛星の検証コストに関する試算、第 54 回宇宙科学技術連合講演会講演集、JSASS-2010-4049、(2013)(1-4)日本経済新聞(Web 版)“ 九工大の小型衛星、１５年度打ち上げへ”、2014/0９/06(1-5) 荒川義博、國中均、中山宜典、西山和孝：”イオンエンジンによる動力航行”(コロナ社 2006)(1-6) 江崎徹：低電力ホールスラスタの開発、平成 22 年度九州大学修士論文(1-7) 國中均：
「はやぶさ」小惑星探査機に搭載されたマイクロ波放電式イオンエンジン、
プラズマ・核融合学会誌 82(5), 300-305, 2006-05-25、(2006)(1-8) 細田聡史，國中 均：イオンエンジンによる小惑星探査機「はやぶさ」の帰還運用、
J. Plasma Fusion Res. Vol.86, No.5 (2010)282-292(1-9) 高村秀一：プラズマ理工学入門、森北出版(1997)(1-10) N. Yamamoto, H. Kataharada, H. Masui, H. Ijiri and H. Nakashima, AJCPP2005-22093,
2005.

(1-11) N. Yamamoto, H. Kataharada, T. Chikaoka, H.Masui and H. Nakashima, IEPC-2005-036,
2005.

(1-12) T. Chikaoka, S. Kondo, N. Yamamoto, H. Nakashima and Y. Takao, Proceedings of the25th International Symposium on Space Technology and Science, 2006, pp.254-259.

(1-13) 鶴哲平， “マイクロ波放電式中和器の性能向上に関する研究” , 平成 20 年度九州大学修士論
文．

(1-14) 増井博一， “マイクロ波放電型宇宙推進用プラズマ源に関する数値解析”, 平成 17 年度九州大学博士論文．

(1-15) 金川隆保，“マイクロ波放電式小型イオンエンジンの数値解析”， 平成 19 年度九
州大学修士論文．
(1-16) 新谷將， “マイクロ波放電式小型イオンエンジンの数値解析”, 平成 22 年度九州大学修士論
文．

9

第２章

数値解析

2.1 粒子計算エラー! 参照元が見つかりません。
2.1.1 基礎式
プラズマの挙動を解析する際の基礎式として運動方程式を用いる．式を以下に示す．

mdv
 q E  v  B dt(2.1.1)dx
v
dt(2.1.2)

ここで， m は粒子質量， v は粒子の速度， x は粒子の位置， t は時間， q は粒子の持つ電荷， E は電場， B は磁束密度を表わす．一般的な実験室プラズマは 1 cm3 当りに 1010～1014 個が含まれ，この粒子全部を取り扱うことは最新のスーパーコンピュータを用いても不可能であり，実用的でない．そこで粒子シミュレーションで取り扱う粒子として，実際の粒子の電荷質量比を一定に保ったまま多数の粒子の電荷と質量を１つにまとめた超粒子エラー! 参照元が見つかりません。を使用する．超粒子数は，１セルに含まれる粒子数が 10 個以上になる事と取り扱う粒子数よりも超粒子の重みが尐ない事（統計性から）から決定する．
（2.2）式に leap-frog 法による時間中心差分を施すと，

x n1  x n
 v n1 2
t(2.1.3)となる．ここで，各変数の上付き添え字は時間に関するもので， 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.1.4)である． n-1 2 はある時刻 n から t 2 だけ以前の値であることを示す．同様に， n  1 2 はある時刻
n から t

2 だけ進んだ値であることを示す．

プラズマ粒子の運動を求めるために，(2.4)式において未知量である v n 1 2 について解き，求めたv n 1 2 と(2.3)式より新たな粒子の位置 x n  1 を求めることができる．(2.4)式を変形すると，
Av n1 2  S(2.1.5)と書ける．ここで， a ・  ・  ・  を次のように定義する．



qt n
Bx
2m

(2.1.6)



qt n
By
2m

(2.1.7)

10

 

qt n
Bz
2m

(2.1.8)  1  2   2   2(2.1.9)

なお，ここでの 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.1.10)(2.1.11)

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

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

(2.1.12)(2.1.13)2.1.2 メッシュサイズ・タイムステップエラー! 参照元が見つかりません。
計算に用いるメッシュサイズ  は，電子温度 Te，ボルツマン定数 kb，電子密度 ne から計算される．
デバイ長 λDe，De 

 0 k bTene e 2(2.1.14)に対して，

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

ne e 2me  0 pe 

(2.1.16)に対して，

t   pe111(2.1.17)

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

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

c 

eB
me

(2.1.19)に対して，

t  c1(2.1.20)である必要がある．廣池らによると尐なくとも 100 分の一にしないとエネルギーが保存されないことが明らかになったため，

Dt <<

wc

(2.1.21)100となるように設定した．

2.1.3PIC 法エラー! 参照元が見つかりません。
シミュレーションでは，電磁場量を空間格子点上で求める．そのため，計算粒子の電荷，電流を空間の格子点上に分配し，マクスウェル方程式を用いて電磁場を計算しなければならない．また，粒子が感じる電磁場を計算する際，格子点上に割り振られた電磁場を内挿する必要がある．このような計算を行うため，粒子シミュレーションにおいて標準的な手法である PIC 法（Particle In Cell method）
を用いる．

ま ず ， 粒 子 の 位 置 か ら 近 傍 の セ ル 座 標 (i , j , k )を 求 め る ． 次 に ， 図 2.1 に 示 す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.1.22)となり，

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)
12(2.1.23)

と表わされる．

図 2.1.1 PIC 法におけるセルの構成

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

V1
qN
V2

(2.1.24)となる．すべての粒子について和をとれば，格子点状での電荷密度が得られる．また，電流密度はj  qnv(2.1.25)より，

j (i , j , k ) =

V1
qNv
V2

(2.1.26)

2.1.5 荷電粒子と中性粒子の衝突
荷電粒子は電磁場の影響を受けて，中性粒子と衝突する．プラズマは電子，イオン，
中性粒子からなるとする．電子とイオンの数密度は中性粒子の数密度より非常に小さい
ため，電子-電子，電子-イオン，イオン-イオンの各衝突は無視できるものとする．電子と中性粒子は弾性，励起，電離衝突をすると仮定する．なお，以下に用いる乱数エラー!
参照元が見つかりません。

は(0,1)間に一様に分布する乱数である．
13

2.1.6 平均自由行程と透過度エラー! 参照元が見つかりません。
面  ，厚さ l ，その物質の数密度 N の板を考える．エネルギー E の粒子ビームが，面  に垂直に当るものとする．衝突をせずに板を透過してくる粒子数は，入射粒子数より減尐する.その減尐度は
 を衝突の断面積とすれば，dn  nNdl(2.1.27)n  n0 exp  Nl (2.1.28)と表わされ，これよりとなる．ここで， n0 は面  （ l  0 ）でのビーム中の粒子である．それゆえ，（ l ,l  dl ）間で最初の衝突が起こる確率は，

Pinitial l dl  exp  Nl Ndl(2.1.29)となり， l 以下で最初に衝突が起こる確率は，P ( l ) = ò Pinitial ( l )dll0= 1- exp ( -Ns l )(2.1.30)となる．最初の衝突までに粒子が走る平均距離  は Pl の一次のモーメントである．
   lP l dl
0



1
N

(2.1.31)

これを平均自由行程という．よってモンテカルロ法によって，任意の点から出発し最初の衝突が起こるまでの距離 l を決めるには，乱数 r を用いて，r  Pl  l
 1  exp    (2.1.32)で表わされ，これを変形して，l  ln 1  r となる．ここで， r は一様乱数であるから 1  r を r とすれば，l  ln(r )(2.1.33)(2.1.34)と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を  とすれば，この粒子が１秒間に気体分子と衝突する回数 は，



v


 Nv
14

(2.1.35)

となり，これを衝突周波数という．

2.1.7 衝突確率・衝突判定
すべての衝突の断面積は，荷電粒子の全エネルギーの関数として表わされる．

   E 

(2.1.36)

荷電粒子速度がわかれば，荷電粒子と中性粒子の各衝突における衝突周波数がわかり，全衝突周波数

 total がわかる．したがって，ある粒子 m の衝突確率 Pcollision,m は以下の式で表わすことができる．
 total  ngas totalvm
Pcollision,m  1  exp  ngas totalvm t 
 1  exp   total t 

(2.1.37)(2.1.38)

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

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

 0  r2  elastic
 total

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

ここで， total   elastic   excitation   ionization である．このように断面積の比と乱数を対応させ，衝突の種類を決める．計算では中性粒子としてキセノンとアルゴンを用い，その断面積データエラー! 参照元が見つかりません。エラー! 参照元が見つかりません。エラー! 参照元が見つかりません。

の計算では，多価電離は考慮していない．

15

はそれぞれ図 2.1.2 と図 2.1.3 に示す．今回

elastic
excitation
ionization

4.5
4.0

σ (×10-19m2)

σ(×10-19m2)

3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

2

4

6

8

10

Energy1/2(eV)1/2)
1/2(eV
Energy

図 2.1.2 電子とキセノンの衝突断面積

elastic
excitation
ionization

4.5
4.0

σ (×10-19m2)

σ(×10-19m2)

3.5
3.0
2.5
2.0
1.5
1.0
0.5
0.0
0

2

4

6

8

Energy1/2
1/2(eV(eV)
1/2)
Energy

図 2.1.3 電子とアルゴンの衝突断面積1610

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

v1 

1
v1  v 2  u2(2.1.39)v 2 

1
v1  v 2  u2(2.1.40)と表わすことができる．ここで， u  v 2  v1 （相対速度）である．電子と中性粒子の衝突も同様に考える．衝突前の電子，中性粒子の速度をそれぞれ v e ・ v n の衝突後の速度をそれぞれ v e ・ v n とす
る．電子と中性粒子の衝突では，電子は中性粒子に比べ質量が非常に小さいため，電子のみが散乱される．中性粒子の速度 v n を無視すれば，相対速度は u  v n  v e  v e としてよい．よって，電子の質量を m ，中性粒子の質量を M とすれば，次式が成立するエラー! 参照元が見つかりません。．

v¢e =m + M cos c
M sin c
ve +
h
m+M
m+M

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

vex vey cos   ve vez sin vervex vez cos   ve vey sin ver(2.1.42)(2.1.43)(2.1.44)

ここで，

ver 

ve 

2
2
vey
 vez

(2.1.45)

2
2
vex
 vey
 vez2(2.1.46)であり，  は散乱角， は任意の角である. エラー! 参照元が見つかりません。式はエネルギー損失を含んでいる．
次に，衝突前後のエネルギー収支を考える．衝突前の電子のエネルギーを Eincident,e とし，衝突後の電子のエネルギーを Escattered,e とすれば，エネルギー収支は以下のように表わされる．

Escattered,e  Eincident,e(2.1.47)

電子は，エネルギー Eincident,e が小さいとき等方散乱され，大きいときは为に前方散乱となる．これ
17

を考慮した散乱角  の確率密度を g   として以下に示すエラー! 参照元が見つかりません。．ここで， Eincident,e の単位を eV とする．

g( c ) =Eincident,e sin cc
2 1+ Eincident,e sin 2 ln 1+ Eincident,e2{(0 £ c £ p )}(2.1.48)

これを解けば，{cos c =1+

2 1- (1+ Eincident,e ) 3r}Eincident,e(2.1.49)となる．この値をエラー! 参照元が見つかりません。式に用いて衝突後の速度を求める．このとき，任意の角は

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

Escattered,e  Eincident,e  Eexcitation(2.1.51)Ev~  v e 1  excitation
Eincident,e(2.1.52)と表わされ，励起後の速度はとなる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失った後にエラー!
参照元が見つかりません。式中 v~ の速度で弾性衝突をすると考えることができる．エラー! 参照元が
見つかりません。式中 v e に v~ を用いて，衝突後の速度 v e を求める．このときの散乱角にはエラー! 参
照元が見つかりません。式を用いるがエラー! 参照元が見つかりません。式中の Eincident,e にはエラー! 参照元が見つかりません。式から求められる Escattered,e を用いる．
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを Escattered,e ，生成電子のエネルギーを
Ecreated,e ，電離エネルギーを Eionization とする．エネルギー収支はEscattered,e  Ecreated,e  Eincident,e  Eionizationと表わされる．
電離後の散乱電子のエネルギーは以下のように求められる(2-12)．

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

(2.1.54)

ここで，B はガスに固有の値である.
また，入射電子が電離によって失うエネルギーを E とすると，

E  Eincident,e  Escattered,e(2.1.55)と表わされる．エラー! 参照元が見つかりません。式から電離後の電子の速度は，v~  v e 1 EEincident,e(2.1.56)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．エラー! 参照元が見つかりません。式中 v e に v~ を用いて，衝突後の速度 v e を求める．散乱角にはエラー! 参照元が見つかりません。式を用いるが，エラー! 参照元が見つかりません。式中の Eincident,e にはエラー! 参照元が見つかりません。式から求められる Escattered,e を用いる．電離エネルギー Eionization は，キセノン,アルゴンの場合それぞれ 12.13,15.76 eV である．
電離後の生成電子の持つエネルギーは，Ecreated,e  E  Eionization(2.1.57)となる．また，入射電子の速度 v e を用いて生成電子の速度はv created,e 

2 Ecreated,e

ve
ve

m

(2.1.58)と表わされる．この速度で弾性衝突をすると考えれば， v created,e を v e の代わりにエラー! 参照元が見つかりません。式に代入して，衝突後の速度を求める．このときの散乱角は生成電子のエネルギー

Ecreated,e を用いてエラー! 参照元が見つかりません。式から求める．
2.1.10Null-collision 法エラー! 参照元が見つかりません。エラー! 参照元が見つかりません。エラー! 参照元が見つかりません。
計算時間短縮のため，Null-collision 法を用いる．これには，架空の衝突断面積  fake を導入する．
すべての電子のエネルギー E に対して

 fake   total(2.1.59)であるとし，例えば電子と中性粒子の衝突においては total   elastic   excitation   ionization   fake(2.1.60)とする．ここで， totalv  collisionn gas

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

null
elastic
excitation
ionization

null
elastic (×10-11m3/s)

ν/n(×10-13m3/s)

1

0.1

0.01

1E-3

1E-4
0.01

0.1

1

10

100

A (eV)
ENERGY
Energy(eV)

図 2.1.4 キセノンの数密度で規格化した電子のキセノンの衝突断面積

20

null
elastic
excitation
ionization

excitation (×10-11m3/s)

ν/n(×10-13m3/s)

1

0.1

0.01

1E-3

1E-4

1E-5
0.01

0.1

1

10

100

ENERGY (eV)Energy(eV)

図 2.1.5 アルゴンの数密度で規格化した電子のアルゴンの衝突断面積

2.2 電磁波解析手法エラー! 参照元が見つかりません。エラー! 参照元が見つかりません。
2.2.1 基礎式
マイクロ波伝播解析には，FDTD 法（Finite Difference Time Domain method）を用いる．FDTD
法は数値電磁場解析の分野において，モーメント法や有限要素法と同様に，アンテナや平面回路など
数多くの複雑な問題に適用され，为要な解析手法として多く用いられている．計算には，マクスウェ
ル方程式を基礎式として用いる．以下に基礎式として用いる微分形式のマクスウェル方程式を示す．

H
t(2.2.1)E
Jt(2.2.2)

  E   0
  H  0

ここで， E は電場， H は磁場，  0 は真空の透磁率，  0 は真空の誘電率， J は電流密度を表わす．

21

2.2.2.差分式
三次元の FDTD 法で用いられる電場，磁場の空間配置を図 2.2.1 に示す．このように，電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．
これは Yee によって提案されエラー! 参照元が見つかりません。，電場の回転（   E ）が磁場を，磁場の回転（   H ）が電場の方程式を満たすように配置されている．
この電磁場の配置のもと， エラー! 参照元が見つかりません。式, エラー! 参照元が見つかりません。式を空間，時間とも中心差分した式を以下に示す．(エラー! 参照元が見つかりません。~エラー!
参照元が見つかりません。)式ここで， x, y, z は空間幅， t は時間幅，添え字 i, j , k  はセル中での座標，添え字 n は時間ステップを表わす．
n+

Hx

1
2

n-

1
2

(i, j +1 2, k +1 2) = H x (i, j +1 2, k +1 2) +Dt
m0 (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
ïî
ïþ
1
1
n+
nDt
H y 2 (i +1 2, j, k +1 2) = H y 2 (i +1 2, j, k +1 2) +m0 (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
ïî
ïþ
1
1
n+
nDt
H z 2 (i +1 2, j +1 2, k ) = H z 2 (i +1 2, j +1 2, k ) +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
ïî
ïþ
E xn 1 i  1 2 , j , k   E xn i  1 2 , j , k  



Hy

1
2

n

1
2

i  1 2 , j , k  1 2  H y i  1 2 , j , k  1 2E yn 1 i , j  1 2 , k   E yn i , j  1 2 , k  (2.2.4)(2.2.5)t 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

(2.2.3)z(2.2.6) 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


1
1
n
n

H z 2 i  1 2 , j  1 2 , k   H z 2 i  1 2 , j  1 2 , k  J y i , j  1 2 , k 
x



22

(2.2.7)

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
 H y i  1 2 , j , k  1 2  H y i  1 2 , j , k  1 2

x



(2.2.8)

1
1
n
n

H x 2 i , j  1 2 , k  1 2  H x 2 i , j  1 2 , k  1 2
 J z i , j , k  1 2y(エラー! 参照元が見つかりません。~エラー! 参照元が見つかりません。)式から，時間ステップがn  1 2 での磁場 H

n

1
2

は，その半ステップ前での電場 E n と１ステップ前の磁場 Hとができ，時間ステップが n  1 での電場 E

n 1

は，半ステップ前の磁場 HnE から求めることができる．初期値については， H

n

1
2

n

FDTD 法における電場と磁場の配置

23

1
2

1
2

から求めること１ステップ前の電場

は E から計算されるので，アンテナの動径方向に E rn だけを与える．

図 2.2.1nn

2.2.3 アンテナ解析
本計算ではアンテナを完全導体で構成されているものとしてモデル化する．給電方法は同軸ケーブ
ルを模擬し，図 2.2.2 に示すように動径方向の４つの電場を中心導体に与える同軸線路給電を用いる.
また, アンテナ導入部には図 2.2.3 のように後述する Mur の一次吸収境界条件を与え, 波長に比べて太さが無視できるくらい十分細い直線状のアンテナを解析する場合，その中心軸に沿って電界を 0とすればよい．また, 今回は L 字型アンテナを模擬した．給電方式について詳しく述べると，中心導体に与える電場の x 方向成分を Ex，y 方向成分を Ey は近似的に
𝑉
𝐸𝑥 =
𝛥𝑥
(2.2.9)𝑟𝑙𝑛 ( )𝑎𝑉
𝐸𝑥 =𝛥𝑦
(2.2.10)𝑟𝑙𝑛 ( )𝑎で与えられる. ただし, V は中心導体と外部導体間の電圧, 𝑟 は動径方向の距離, 𝑎 はアンテナの半径である.
また, 同軸ケーブルのインピーダンス Z は以下の式で表される.
𝜇 1 𝑅
ln
𝜀 2𝜋 𝑟𝑍=√(2.2.11)

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

𝑉2
2𝑍

(2.2.12)となり，エラー! 参照元が見つかりません。式，エラー! 参照元が見つかりません。式より，V を求め，エラー! 参照元が見つかりません。式，エラー! 参照元が見つかりません。式よりアンテナに与える電場を求める．ここで，Ex を 𝑟 = 𝛥𝑥/2の電場で，Ey を 𝑟 = 𝛥𝑦/2 の電場で代表させると，

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

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

24(2.2.14)図 2.2.2

同軸線路給電図 2.2.3 アンテナモデル

25

2.2.4 安定条件
計算結果の正確さを保つため，空間の増分値 r , z は波長  に比べて十分小さくなければならな
い．
（  10 以下）
時間ステップ t は，Courant の安定化条件より

1

Dt £
v

1
1
1
+
+
(Dx )2 (Dy )2 (Dz )2(2.2.15)を満たさなければならない．この条件は極めて厳しく，わずかでも満足しなければ計算が不安定になる．

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

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

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

(2.2.16)と表わされ，次の微分方程式を満足する．

E z 1 E z

0
x
v t

(2.2.17)x = 0 の境界で完全に反射がないならば，電場はエラー! 参照元が見つかりません。式の形を保ったまま伝播するはずである．すなわち，境界でもエラー! 参照元が見つかりません。式が満足される．
これを FDTD 法の形式に定式化するため，時間について差分化すると¶E z E zn - E zn -1
=
¶t
Dt
n-

¶E
=v z
¶x

1
2

(2.2.18)となる．右辺の x についての微分は，電場が x = 0, Dx に割り当てられていることから， x = Dx 2 で行
い，エラー! 参照元が見つかりません。式に代入すると

26

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zn

(2.2.19)  n1 1 2 も FDTD 法では直接求めることとなる．エラー! 参照元が見つかりません。式の Ez 1 2 , Ez

はできないため，前後の値で平均して求める．そして， 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.2.20)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．

2.2.7 ポインティングベクトルエラー! 参照元が見つかりません。
空間を流れる電力を求めるためには，ポインティングベクトルを使用する．Maxwell の方程式を変形すると，以下の式となる．

  E 2 H 2 dv  E  H   dS   J  Edv   t vol 2
2 
s
vol

(2.2.21)

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

27

2.3 カップリングコード
本研究では，2.1 で述べたプラズマ解析の PIC 法と 2.2 で述べたマイクロ波解析の FDTD 法を組み合わせることエラー! 参照元が見つかりません。エラー! 参照元が見つかりません。によりプラズマの挙動とマイクロ波の伝播を同時に解析する．図 2.3.1 にカップリングコードのフローチャートを示す．FDTD コードで計算されたマイクロ波の電場は PIC コードに渡され，その電場の効果が粒子の運動に影響する．また，電子の挙動により計算される電流密度が FDTD コードに渡され，電磁波に粒子の効果が加えられる．
このように PIC コード・FDTD コードの間で物理的に相互作用を及ぼすデータをやり取りすること
で，プラズマおよびマイクロ波の挙動を解析する．なお，タイムステップについては PIC 法ではエラー! 参照元が見つかりません。式により決定され，FDTD 法ではエラー! 参照元が見つかりません。
により決定され計算の安定化条件が異なる．必ずしも両者のタイムステップを一致させる必要は無い
ので，本研究では PIC 法の計算負荷を軽減するために FDTD コードを 10 回，連続で計算させた後に PIC コードを 1 回計算することとした．

28

図 2.3.1 Coupling code のフローチャート

29

2.3.1 計算体系

図 2.3.2 計算体系
図 2.3.2 に本研究において使用した計算体系を示す．円筒形の放電室をヨークで囲み，オリフィス
開口部より電子を放出する．スペーサーによってオリフィスとの位置を調節している．磁石の個数は14 個とし，放電室を囲んだ．.

2.3.2 計算格子
PIC 法を用いるために格子を刻む必要があるが，そのために望ましい長さは，デバイ長である．
そうすることによって，壁面でのシースが計算できる．しかし，計算コストの点からデバイ長で格子を刻むことは現時点では不可能である．今回は電子とマイクロ波との相互作用を調べることが目的であるからデバイ長で計算格子を作る必要はない．そこで本計算では計算格子の長さを 0.5 mm とし，
x 方向に 40 格子 y 方向に 40 格子 z 方向に 36 格子合計 57600 格子作成した．そのうちの電子が動く
自由領域は 3832 格子である．

2.3.3 時間刻み
PIC 法の時間刻みは t pic  2.5  10

14

sec，FDTD 法の時間刻みは t FDTD  2.5  10

13

sec とすれば，2.1.2 や 2.2.4 で述べた安定条件を満たす．FDTD コードが 10 回計算された後，PIC コードを1 回計算させた．

30

2.3.4 電子初期条件
Te＝2.0eV となるマクスウェル分布の初期速度を持った電子を放電室内に一様に分布させた．本来であれば，FDTD−PIC コードとは別に，プラズマの挙動を観測するコードの結果を元に電子密度が決定されるが，今回においては，まだ導入に至っていない．

2.3.5 粒子数
1 ボクセル当たり数個しか粒子がない場合，ノイズが大きくなり望ましくないエラー! 参照元が見つかりませ．通常，10 個から 100 個以上の粒子が存在しているとよいとされている．今回は 1 ボクセルに平ん。

均 25 個程度の粒子が含まれるように粒子数を設定した．
上記の計算体系における体積 4.79×10-7m3，実験による粒子密度は 10-17~10-18 m-3 と分かっているので，平均化して 5×10-17m-3 とした．超粒子数は 400 万個とした.

2.3.6 推進剤
推進剤はキセノンとアルゴンとし，中性粒子密度はそれぞれ 2.1×1020 /m3，2.04×1019 /m3 とし,推進剤流量はそれぞれ 0.015mg/s とした．

2.3.7 磁場強度
今回は電源としてはやぶさにも採用されている 2.65GHz を用いる．マイクロ波周波数と電子サイ
クロトロン周波数が一致するときの磁場強度は 0.085T である．その値を示す領域を ECR 層と呼ぶ．
磁場計算には Advanced Science Laboratory, Inc の 3 次元静磁場解析ソフトを使用したエラー! 参照元が見．図 2.3.2 に上述の計算体系における磁場強度及び ECR 位置を示す．

つかりません。

図 2.3.2 磁場強度

31

2.3.8 計算コードの並列化エラー! 参照元が見つかりません。
計算にかかる時間を短縮することはより多くの条件に対して検証するためには非常に重要である．
1 台のコンピュータ(1 つの CPU)で行っていた計算をネットワークでつながれた複数のコンピュータ(複数の CPU)で計算できるように改良している.用いたのは MPI(Message Passing Interface) ，分散メモリ型の並列計算機で 2 つのプロセス 間でのデータをやりとりするために用いるメッセージ通
信操作の仕様標準である．

32

参考文献
(2-1)

内藤裕志, プラズマ・核融合学会誌, 74 (1998) ,470.

(2-2)

Giovanni Lapenta, Jeremiah U. Brackbill, Journal of Computational Physics, 115,

(1994), 213.
(2-3)

田中實, 山本良一：“計算物理学と計算化学―分子動力学法とモンテカルロ法”(海文堂,1998)(2-4)

石黒静児, プラズマ・核融合学会誌, 74 (1998) , 591.

(2-5)

丹慶勝市, 他：“Numerical Recipes in C (日本語版)” (技術評論社, 1993)高村秀一：“プラズマ理工学入門”(森北出版, 1997)(2-6)

高村秀一：“プラズマ理工学入門”(森北出版, 1997)(2-7)

M. Hayashi, J. Phys. D, Appl. Phys., 16 (1983), 581.

(2-8)

M.S. Dababneh, et al., Phys. Rev. A., 22 (1980), 1872.

(2-9)

J. J. Szabo, et al., AIAA 00-3653, (2000)(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)

著=W. H Hayt，訳=山中惣之助，岡本考太郎，宇佐美興一：“工学系の基礎電磁気学”(朝倉書店, 2001)(2-19)

V. P. Gropinath and T. A. Grotjohn, IEEE Trans. Plasma Sci., 23, (1995), 602.

(2-20)T. A. Grotjohn, Rev. Sci. Instrum. 67, (1995), 921.

(2-21)

Y. Saad, and M.H. Schultz , SIAM. J. Sci. Stat. Comput., 7, pp.856-869 (1986)(2-22)

藤野清次，張紹良：”反復法の数理”, (朝倉書店, 1996)(2-23)

樫山和男，西村直志，牛島省：”並列計算法入門”, (丸善, 2003)(2-24)

新谷將， “マイクロ波放電式小型イオンエンジンの数値解析”, 平成 22 年度九州大学修士論文．
(2-25)

C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical Physics, 55 (1971),

4100.
(2-26)

M. Srendra, D. B. Graves and G. M. Jellum, Physical Review A, 41 (1990), 1112

(2-27)

M. Hirakawa, IEPC-97-021, (1997), 159.

(2-28)

http://www.asl-i.com/index.html
33

第3章

計算結果と考察

3.1 アルゴンとキセノンの電界強度
図 3.1.1 にアンテナからマイクロ波が入射していく様子を表すために作動ガスをキセノンとしたときの電界強度の空間分布を,図 3.1.2 に作動ガスをアルゴンとした時の電界強度の空間分布を示す.図はそれぞれアンテナの中心における xz の断面図であり、それぞれの時間は(a)1 ns, (b)25 ns, (c)50 ns である．(a)はマイクロ波を投入して初期の段階であり、アンテナの軸付近で電界強度が高くなっている. (b)ではアンテナの軸に最も電界が強い点がみられ、その周辺に 50000V/m 程度の電界が広がっている. (c)では,アンテナの軸の周辺と同程度の強度をもった電界が ECR 層に存在している.このことから ECR
が上手く作用していることがわかる.また,(c)におけるキセノンとアルゴンの電界強度を比較すると ECR 層においてキセノンは 40000V/m 程度の電界強度の電界が見られ,アルゴンでは 60000V/m 程度の電界強度の電界が見られる.

34

(a)

(b)

(c)
図 3.1 Xe の電子の平均エネルギー分布(a) 1ns (b) 25ns (c) 50ns

35

(a)

(b)

(c)
図 3.1 Xe の電子の平均エネルギー分布(a) 1ns (b) 25ns (c) 50ns

36

3.2 キセノンとアルゴンの電子の平均エネルギー
図 3.1 に作動ガスをキセノンとしたときの電子の平均エネルギーの空間分布を,図 3.2に作動ガスをアルゴンとした時の電子の平均エネルギーの空間分布を示す.図はそれぞれアンテナの中心における yz の断面図であり、それぞれの時間は(a)1 ns, (b)25 ns, (c)50
ns である．まず (a)はマイクロ波を投入して初期の段階であり、(a)において最も平均エネルギーが高いのは、アンテナ先端付近で 12eV 程度である.電界強度においてもアンテ
ナの先端に高い値が集中していたことより,初期段階においては電子の平均エネルギー
は ECR 層よりもアンテナが作り出す電界強度に依存していることがわかる. 次に(b)を見ると，アンテナの先端の下方，つまり BackYork 傍に約 20 eV の高エネルギー値が見られる．ECR 層の位置がこの位置にあるため，ECR 層における共鳴により，高エネル
ギーが集中している事が分かる．また，FrontYork と BackYork との間に磁気チューブが見られ，磁気ミラー効果により，電子の閉じ込めが起きていると確認できる．最後に(c)では，(a),(b)と同じく，高エネルギーの電子は ECR 層よりも，アンテナ近傍に存在
し，これよりアンテナが作り出す強電界による電子の加熱が ECR 加熱よりも支配的であることがわかる．アンテナの軸近傍において，アンテナの先端の下方と同程度の約 25
eV のエネルギー値が分布しており，アンテナの強電界による影響が考えられる.
電子の平均エネルギーはどの時間においても Ar の方が高い.50ns では,アルゴンの電子の平均エネルギーはアンテナ先端付近の点では 30eV となり,全体的には 20eV 程度で広がっている.キセノンでは,このようになった理由として,Ar は衝突断面積が小さいため,
電子との衝突頻度が小さくなり電子の平均エネルギーが効率的に上昇したためであると考えられる.キセノンの電子の平均エネルギーはアンテナの先端付近の点では 20eV となり,全体的には 12eV 程度で広がっている.また,キセノンはアルゴンに比べて電子の平均エネルギーが広がっている範囲が狭いことがわかる.このようになった理由として,アルゴンは衝突断面積が小さいため,電子が中性粒子との衝突によりエネルギーを失う頻度が小さくなり,結果として,電子の平均エネルギーが効率的に上昇したためであると考
えられる.

37

(a)

(b)

(c)
図 3.1 Xe の電子の平均エネルギー分布(a) 1ns (b) 25ns (c) 50ns

38

(a)

(b)

(c)
図 3.2 Ar の電子の平均エネルギー分布(a) 1ns (b) 25ns (c) 50ns

39

3.3 キセノンとアルゴンのイオンの電離数
プラズマの生成量を測るうえで重要な指標となるのが,イオンの電離数である.計算領域でのイオン電離数はキセノン,アルゴンそれぞれ 255077 個,3162 個であった.アルゴンの電離したイオン数はキセノンの 0.0124 倍であり,プラズマの生成量がキセノンよりも大幅に尐なかった.したがって,アルゴンの流量を 80.7 倍にすることによりキセノンと同程度の性能が得られる.しかし,これは現実的ではないので,磁場形状の変更により,電離衝突数を増やす必要があることが示唆された.

40

第４章

結論

電子エネルギーの分布より ECR 層で電子がエネルギーを授受していることが確認された.また, ECR 層よりもアンテナ先端の方がプラズマが生成されやすいことが示唆された.また作動ガスが性能に及ぼす違いにおいて,推進剤がアルゴンである場合プラズマ
生成量という面ではキセノンより务ることがわかった.しかし,アルゴンは衝突断面積が小さいことにより電子の平均エネルギーが効率よく上がるという良い点もある.アルゴンのプラズマ生成量の尐なさを改善するためには今後磁場形状を改善し,電離衝突数を増やす必要がある.

41

謝辞
本論文の研究に当たって１年間研究の場を与えてくださり,また, 様々なご指導を頂いた中島秀紀教授に深く感謝いたします.研究や論文制作において遅くまで指導をしてくださり,また,日頃より親切に話しかけてくださった山本直嗣准教授に深く感謝致します.研究についてご指導頂き,また親切に話しかけてくださった森田事務の事で支えていただき,また明るく接してくださった馬渡隆子事務補佐員に深く感謝いたします.
また,研究室での生活を明るく充実したものにしてくださった研究室の皆様に感謝致します.研究についての疑問点についていつも親身になって答えてくださり,また,明るく接してくださった楠田将

士氏,牛尾康一氏には深く感謝いたします. 研究・論文等でアドバイスをしていただくと共に常に、明るく接してくださった李後毅氏、伊藤匠氏、齋藤直哉氏、坂井寿光氏、木邊厚視氏、田中啓介氏,奥田雄氏,川島諒祐,高瀬紘平氏,中野和彦氏,山口敦氏に感謝致します.
同じ B4 として研究室の生活を楽しいものにしてくれた飯島健介氏,江川雄亮に深く感謝致します。

42