卒 業 論 文

題目

数値解析を利用したマイクロ波放電型中和器における磁場形状の最適化

氏

名

金川 隆保

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

指 導 教 員中島 秀紀 教授

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

平成 18 年 2 月

目次
第1章

序論

1

第2章

イオンエンジンの原理

42.1

イオンエンジン

42.2

マイクロ波放電型イオンエンジン

4

2.3

電子サイクロトロン共鳴

42.4

ECR プラズマの原理

5

第3章

計算手法

8

3.1

緒言

8

3.2

プラズマ解析手法

8

3.2.1

8

3.2.1.1

基礎式

8

3.2.1.2

計算方法

8

3.2.1.3

メッシュサイズ，タイムステップ

9

3.2.2

PIC 法

10

3.2.2.1

内挿計算

10

3.2.2.2

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

11

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

113.2.3.1

平均自由行程と透過度

113.2.3.2

衝突確率，衝突判定

12

3.2.3.3

衝突計算

133.2.3.4

エネルギーと散乱角133.2.3.5

Null-collision 法

16

3.2.3

3.3

粒子シミュレーション

マイクロ波解析手法

19

FDTD 法19

3.3.1.1

基礎式

19

3.3.1.2

差分式

19

3.3.2

安定条件

21

3.3.3

境界条件

21

3.3.1

3.3.3.1

PEC(Perfect Electric Conductor )境界条件

21

3.3.3.2

Mur 吸収境界条件

21

3.3.4
3.4

給電方法

PIC 法と FDTD 法のカップリング

22
25

第4章
4.1

計算結果と考察

小型中和器における電子加速機構

4.1.1 計算体系

27
27
27

4.1.2

計算パラメータと境界条件

29

4.1.3

電場強度分布

29

4.1.4

少数の電子

31

4.1.5 多粒子での解析結果

34

磁場強度依存性

38

4.2

4.2.1

ECR と 2nd Harmonic ECR との比較(少数の粒子)

38

4.2.2

ECR と 2nd Harmonic ECR との比較(多粒子)

41

4.3

放電室径依存性

第5章

結論

46
51

参考文献

52

謝辞

54

第1章

序論

1955 年 4 月、東京大学生産研究所の糸川博士が開発したペンシルロケットの水平発射に成功し、同年 8 月、開設された秋田の実験場でペンシルロケットの打上げが行なわれた。その後の日本のロケット開発の進歩は目覚ましく、糸川博士はその功績から日本のロケット開発の父と呼ばれる事となる。
現在、糸川博士の名を持つ小惑星「イトカワ」へのサンプルリターン計画「MUSES-C」
が宇宙航空研究開発機構により行われている。これまで人類がサンプルを持ち帰った天
体は月だけであるが、月は変成してしまったため、太陽系初期のころの物質について知ることができない。小惑星「イトカワ」は惑星が誕生するころの記録を比較的よくとどめている化石のような天体で、この小惑星からサンプルを持ち帰る技術（サンプル・リターン）が確立されれば、「惑星を作るもとになった材料がどんなものか」、「惑星が誕生するころの太陽系星雲内の様子はどうか」についての手がかりが得られる。また地球上でサンプルの分析が行えるため、回収される量が少量であってもその科学的意義は極めて大きいといえる。
このサンプルリターン計画では小惑星への航行に電気推進機が用いられている (1) 。
電気推進機では太陽電池からの電気エネルギーをガスに与え、ガスをプラズマ状態にし、高速のプラズマジェットを排気して推力を得る。電気推進は推力が小さいため加速には
多くの時間を必要とするが、少ない推進剤で大きな速度変化を生み出す事が可能である。
現在、「イトカワ」から地球に向かって航行中の惑星探査機「はやぶさ」では主推進機
として、電気推進機の一種であるマイクロ波放電型イオンエンジンが使用されている。
イオンエンジンではプラズマ化したガスからイオンを静電的に加速噴出することによって推力を得る。従来のイオンエンジンでは直流放電型が主流であったが、「はやぶさ」
では電子サイクロトロン共鳴（Electron Cyclotron Resonance : ECR）を利用したマイクロ波放電型が採用されている。
MUSES-C 計画のような大規模プロジェクトが行われている一方で、人工衛星を小型化することにより、大学、民間が参入できる小規模のプロジェクトも計画されている。
小型衛星の利用のメリットはプロジェクトが小規模になり実行し易くなるだけではな
く、複数打ち上げる事によるリスクの分散、ネットワーク的運用を可能にする事も大きなメリットである。こうした動向に対して、本研究室では図 1.1 に示す小型衛星用のマ
イクロ波放電型イオンエンジンの開発を行ってきた(2)。イオンエンジンでは衛星の電気的中性を保つために、イオンスラスタから放出されるイオンと同量の電子を放出する中和器と呼ばれる電子源が必要である。中和器は直接推力の発生に寄与しないために低電1力で動作する事，軽量である事が望ましい。しかし、中和器を小型化することにより性能の低下が本研究室の実験グループにより報告されている。またその原因を探るための実験も中和器自体が小型であるため非常に困難であり、内部のプラズマ生成機構は十分解明されていない。
中和器を小型化することによる性能低下の一つの原因として、ECR 層で加熱された電子が壁面に衝突し、プラズマ密度が低下するためと考えられている。損失を抑えるための手段として、2nd Harmonic ECR を用いることにより電子の軌道を閉じ込めることができるのではないかいう報告がある(3)。そこで本研究では数値解析を利用し、内部の電子の軌道に着目した。プラズマを粒子として扱うためにプラズマ解析には PIC 法(Particle In
Cell Method)、マイクロ波伝播には FDTD 法（Finite Difference Time Domain Method）を用
い、マイクロ波とプラズマの挙動を自己無撞着に計算する PIC+FDTD コードを使用した。
この論文の構成は以下の通りである。
第１章では研究の背景と目的について述べた。
第２章では ECR プラズマ推進システムについて述べる。
第３章ではマイクロ波による電磁波の解析にもちいる用いる FDTD 法、プラズマ中の粒子の挙動解析に用いる PIC 法、マイクロ波の伝播およびプラズマの挙動を自己無撞着に計算するための FDTD+PIC カップリングコードついて述べる。
第４章では計算結果について述べる。
第５章では結論を述べ，本研究の総括を行う。

2

図 1.1

小型衛星用のマイクロ波放電型イオンエンジン

3

第 2 章 イオンエンジンの原理

2.1

イオンエンジン

イオンエンジンとは燃料となる推進剤をプラズマ化し、3 枚からなる静電グリッドによりイオンを後方に噴射し、その反作用により推力を得る装置である。図 2.1 に概念図を示す。イオンエンジンの物理機構は主に以下の 3 つに分けられる。
1.推進剤を電離するイオン生成部 (Ionization)2.生成されたイオンを静電的に加速して推力を得る加速部 (Acceleration)3.放出されたイオンビームを電気的に中和する中和部(Neutralization)
これらの各過程はそれぞれイオン源、加速電極、中和器によって行われる。

2.2

マイクロ波放電型イオンエンジン

イオンエンジンはプラズマ生成方法により直流放電型、マイクロ波放電型などに分類される。本研究の対象であるマイクロ波放電型イオンエンジンの概念図を図 2.2 に示す。
マイクロ波放電とは、マイクロ波帯域の交流電場によって電子を加速し、中性粒子と衝突電離により電子の数が増加し気体がプラズマ化され放電が維持されることである。このマイクロ波放電は直流放電型と比べた場合、以下のようなメリットがある。
・プラズマ生成に電極を必要とせず、それによるエンジン寿命の制約がない
・ヒーターを用いないため、予備加熱が不要となりエンジンを即座に立ち上げることが出来る。
・1 台のマイクロ波電源で複数台のエンジンにプラズマを生成できるため、中和器のプラズマも同時に生成することが可能となり、システム構成の簡素化、信頼性の向上、衛星の重量軽減をもたらす。
マイクロ波放電型イオンエンジンはプラズマの生成に電子サイクロトロン共鳴(Electron Cyclotron Resonance : ECR) を利用し、電子の加熱効率を上げている。

2.3

電子サイクロトロン共鳴

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

4

m

dv
= −ev × B
dt

(2.1)

ここで m は電子の質量、v は速度ベクトル、e は電荷量、B は磁束密度である。このときの円軌道の半径 rL はラーマ半径と呼ばれ

rL =

mv ⊥
eB

(2.2)

で与えられる。ここで v⊥は電子の B に対する垂直な速度成分の大きさである。
プラズマを構成している荷電粒子は全て反磁性体である。そのためサイクロトロン運動の回転の方向は、外部磁場の向きに対して荷電粒子の回転によってできる磁場が常に逆を向く方向である。つまり電子は磁場に対して右回りの回転運動を行う。この回転運
動の角周波数はサイクロトロン角周波数と呼ばれ

ωce =

eB
(2.3)

m

で与えられる。電子は電場と逆向きに加速されるため、磁場中の電子の回転方向と逆向
きに、回転する速さが等しい電場をかけると電子は連続的に加速され、電場から効率的にエネルギーを受けとることができる。これが電子サイクロトロン共鳴現象である。

2.4

ECR プラズマの原理

ECR プラズマ発生装置では、導波管やアンテナを用いて角周波数ωrf のマイクロ波を入射し、プラズマを生成する。生成されたプラズマ中の電子はサイクロトロン角周波数
ωce に従って、磁力線方向に向かって右回りの旋回運動を行う。一方で、強磁場側から入射したマイクロ波はプラズマ中を浸透し、電子サイクロトロン波と呼ばれる右回りの円偏波を励起する。この波は弱磁場側へ伝播し、角周波数が電子サイクロトロン角周波数と一致する層で急速に減少してマイクロ波から電子にエネルギーが吸収される。
(2.4)
ωce = ωrf
この関係が成り立つときに ECR が生じる。
また逆に、ECR を起こすときに必要な磁束密度の大きさを Bres とすれば(2.3)，(2.4) 式よりBres =

mω rf

(2.5)

q

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

図 2.1 イオンエンジンの構成図

図 2.2 マイクロ波放電型イオンエンジンシステムの概念図
6

E

B

電子

E

B

電子

電子
B

E
B
電子

図 2.3

ECR の原理

図 2.4

ECR の概念図

7

E

第3章

3.1

計算手法

緒言

この章では，本計算で使用したカップリングコードの詳細について述べる。プラズマを解析する手法としては，プラズマを電磁流体として取り扱う方法とプラズマを粒子として取り扱う方法に大別される。本計算では特に電子の放電室内での軌道，中性粒子との衝突などの電子の粒子としての挙動を詳しく解析する為に，プラズマを粒子として取り扱う粒子法を採用した。その軌道計算には粒子法の中でも一般的に良く使われている
PIC(particle-in-cell) 法 を 採 用 し た 。 ま た ， マ イ ク ロ 波 の 解 析 に は FDTD
（finite-difference-time-domain）法と呼ばれる電磁波の解析に最適な手法を採用した。
FDTD 法はアルゴリズムが単純である事，境界条件の取り扱いが容易などの特徴を持つ。
以下にプラズマ解析手法，マイクロ波解析手法の順に説明を行い，両手法のカップリン
グ方式を述べる。

3.2

プラズマ解析手法

3.2.1

粒子シミュレーション

3.2.1.1

基礎式(4)

プラズマの挙動を解析する際の基礎式として運動方程式を用いる。式を以下に示す。

mdv= q(E + v × B )
dt

(3.1)

dx
=v
dt

(3.2)

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

3.2.1.2

計算方法（3.2）式に leap-frog 法による時間中心差分を施すと，
8

x n +1 − x n= v n +1 2
∆t

(3.3)

となる。ここで，各変数の上付き添え字は時間に関するもので， n はタイムステップである。（3.1）式に時間中心差分を施すと，
v n +1 2 − v n-1 2 q ⎧ n 1 n +1 2⎫
= ⎨E + v+ v n-1 2 × B n ⎬
m⎩
2
∆t
⎭

(

)

(3.4)

である。n- 1 2 はある時刻 n から ∆t 2 だけ以前の値であることを示す。同様に，n + 1 2 はある時刻 n から ∆t 2 だけ進んだ値であることを示す。
プラズマ粒子の運動を求めるために，（3.4）式において未知量である v n +1 2 について解き，求めた v n +1 2 と（3.3）式より新たな粒子の位置 x n + 1 を求めることができる。
（3.4）式を変形すると，

Av n +1 2 = S(3.5)と書ける。ここで， α ， β ， γ ， δ を次のように定義する。

α=

q∆t n
Bx
2m

(3.6)

β=

q∆t n
By
2m

(3.7)

γ =

q∆t n
Bz
2m

(3.8)

δ = 1+α 2 + β 2 + γ 2(3.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

(3.10)(3.11)

ここで， c = q∆t m ， m = n − 1 2 である。さらに（3.5）式を次のようにして変形し，速度を求めることができる。

v n +1 2 = A −1 S
⎛ 1+α 21⎜
A −1 = ⎜ αβ − γ
δ⎜⎝ αγ + βαβ + γ αγ − β ⎞⎟1 + β 2 βγ + α ⎟βγ − α 1 + γ 2 ⎟⎠9

(3.12)(3.13)

3.2.1.3

メッシュサイズ，タイムステップ(5)計算に用いるメッシュサイズ ∆ は，電子温度 Te，ボルツマン定数 kb，電子密度 ne から計算される。
デバイ長λDe，

ε 0 k bTeλ De =ne e 2(3.14)に対して，∆ ≤ λ De(3.15)であることが望ましい。また，タイムステップはプラズマ周波数 ωpe，

ne e 2me ε 0ω pe =(3.16)に対して，∆t ≤ ω pe−1

(3.17)とすれば，数値的な不安定性は生じない。時間積分の精度まで考えると，通常は∆t ≤ 0.2ω pe−1
(3.18)が用いられる。また，磁場中の荷電粒子のサイクロトロン運動にともなう現象をシミュレートするためには，電子サイクロトロン周波数 ωc

ωc =

eB
me

(3.19)

に対して，∆t << ω c

−1

(3.20)

である必要がある。
3.2.2
3.2.2.1

PIC 法(6)内挿計算

シミュレーションでは，電磁場量を空間格子点上で求める。そのため，計算粒子の電荷，電流を空間の格子点上に分配し，マクスウェル方程式を用いて電磁場を計算しなければならない。また，粒子が感じる電磁場を計算する際，格子点上に割り振られた電磁場を内挿する必要がある。このような計算を行うため，粒子シミュレーションにおいて標準的な手法である PIC 法（Particle In Cell method）を用いる。
ま ず ， 粒 子 の 位 置 か ら 近 傍 の セ ル 座 標 (i , j , k ) を 求 め る 。 次 に ， 図 3.1 に 示 すdx1, dx 2, dy1, dy 2, dz1, dz 2 を求める。その値を用いて，図中にあるような直方体 V1〜V8 の体積を求める。この体積比から粒子が感じる電磁場を内挿することができる。
具体的には，座標 (x , y , z ) 上にある粒子が感じる電場 E ( x , y , z ) は，計算セルの体積を V ，
図に示すような粒子によって仕切られた空間体積をそれぞれ V1〜V8 とすると，格子点上

10

の重みは
An =

Vn
V

(n = 1 〜 8)

(3.21)となり，
E ( x , y , z ) = A1 E (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)
(3.22)と表わされる。

3.2.2.2

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

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

ρ (i , j , k ) =

V1
qN
V2

(3.23)となる。すべての粒子について和をとれば，格子点状での電荷密度が得られる。また，
電流密度は

j = qnvより，
j (i , j , k ) =

V1
qNv
V2

(3.24)(3.25)

PIC 法では(3.23)式から求まる電荷密度を元にポアソン方程式を解くことにより，静電場の効果をシミュレート可能となる。静電場を考慮する為には(3.15)式を満たす事が必要である。(3.15)式の意味は，実際のプラズマ中での電場はデバイ長で遮蔽される為，この効果を PIC 法において再現する場合は，デバイ長以内のセルの電場の影響だけを考慮す
ればよい事を示している。一般のプラズマの場合，デバイ長は λDe ~ 10 −4 m である。
3.2.3

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

荷電粒子は電磁場の影響を受けて，中性粒子と衝突する。プラズマは電子，イオン，
中性粒子からなるとする。電子とイオンの数密度は中性粒子の数密度より非常に小さい
ため，電子-電子，電子-イオン，イオン-イオンの各衝突は無視できるものとする。電子と中性粒子は弾性，励起，電離衝突をすると仮定する。なお，以下に用いる乱数(7)は（0,1）間に一様に分布する乱数である。

3.2.3.1

平均自由行程と透過度

面 α ，厚さ l ，その物質の数密度 N の板を考える。エネルギー E の粒子ビームが，面 α11

に垂直に当るものとする。衝突をせずに板を透過してくる粒子数は，入射粒子数より減少する。その減少度は σ を衝突の断面積とすれば，dn = − nNσdlと表わされ，これよりn = n0 exp(− Nσl )(3.26)(3.27)となる。ここで， n0 は面 α （ l = 0 ）でのビーム中の粒子である。それゆえ，（ l ,l + dl ）間で最初の衝突が起こる確率は，
Pinitial (l )dl = exp(− Nσl )Nσdl(3.28)となり， l 以下で最初に衝突が起こる確率は，P(l ) = ∫ Pinitial (l )dll0= 1 − exp(− Nσl )(3.29)となる。最初の衝突までに粒子が走る平均距離 λ は P(l ) の一次のモーメントである。
∞

λ = ∫ lP(l )dl
0

=

1
Nσ

(3.30)

これを平均自由行程という。よってモンテカルロ法によって，任意の点から出発し最初の衝突が起こるまでの距離 l を決めるには，乱数 r を用いて，

r = P(l )⎛ l⎞= 1 − exp⎜ − ⎟⎝ λ⎠で表わされ，これを変形して，

l = −λln(1 − r )(3.31)(3.32)となる。ここで， r は一様乱数であるから 1 − r を r とすれば，
l = −λln(r )(3.33)と表わされる。また，速さ v の粒子が気体分子と衝突する際の断面積を σ とすれば，この粒子が１秒間に気体分子と衝突する回数 υ は，

υ=

v

λ
= σNv

(3.34)となり，これを衝突周波数という。

3.2.3.2

衝突確率，衝突判定

すべての衝突の断面積は，荷電粒子の全エネルギーの関数として表わされる。
σ = σ (E )(3.35)

12

荷電粒子速度がわかれば，荷電粒子と中性粒子の各衝突における衝突周波数がわかり，
全衝突周波数 υ total がわかる。したがって，ある粒子 m の衝突確率 Pcollision,m は以下の式で表わすことができる。

υ total = ngasσ total v m

Pcollision,m = 1 − exp(− ngasσ total v m ∆t )= 1 − exp(− υ total ∆t )(3.36)(3.37)

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

3.2.3.3

衝突計算

乱数 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 である。このように断面積の比と乱数を対応させ，衝突の種類を決める。計算では中性粒子としてキセノンを用い，その断面積
データ(8) (9) (10)を図 3.2 に示す。
3.2.3.4

エネルギーと散乱角

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

v1′ =

1
(v1 + v 2 − u′)213(3.38)

v ′2 =

1
(v1 + v 2 + u′)2(3.39)と表わすことができる。ここで， u′ = v ′2 − v1′ （相対速度）である。電子と中性粒子の衝突も同様に考える。衝突前の電子，中性粒子の速度をそれぞれ v e ， v n ，衝突後の速度をそれぞれ v e′ ， v ′n とする。電子と中性粒子の衝突では，電子は中性粒子に比べ質量が非常に小さいため，電子のみが散乱される。中性粒子の速度 v n を無視すれば，相対速度は

u = v n − v e ≅ −v e としてよい。よって，電子の質量を m ，中性粒子の質量を M とすれば，次式が成立する(11)。

v e′ =

m + M cos χM sin χ
ve +
h
m+M
m+M

(3.40)

h の直交座標成分は

hy = −
hz = −

hx = v er cos φ(3.41)

vex vey cos φ − v e vez sin φ(3.42)vervex vez cos φ − ve vey sin φ(3.43)ver

ここで，ver = vey2 + v ez2(3.44)

ve = v ex2 + v ey2 + v ez2(3.45)

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

Escattered,e = Eincident,e(3.46)

電子は，エネルギー Eincident,e が小さいとき等方散乱され，大きいときは主に前方散乱となる。これを考慮した散乱角 χ の確率密度を g (χ ) として以下に示す(12)。ここで， Eincident,e の単位を eV とする。

g (χ ) =これを解けば，

Eincident ,e sin χ(0 ≤ χ ≤ π )χ
⎫⎧2⎨1 + Eincident ,e sin 2 ln 1 + Eincident ,e ⎬
2
⎭
⎩

{

2 1 − (1 + Eincident, e ) 3
cos χ = 1 +Eincident, er

14(3.47)}(3.48)

となる。この値を(3.40)式に用いて衝突後の速度を求める。このとき，任意の角は
φ = 2π r4(3.49)
として与える。
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーをE excitation とし，入射電子のエネルギーを Eincident,e ，励起後のエネルギーを Escattered,e ，励起後の速度を v~ とする。エネルギー収支は，Escattered,e = Eincident,e − E excitation(3.50)と表わされ，励起後の速度はE
v~ = v e 1 − excitation
Eincident,e(3.51)となる。励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失った
後に(3.51)式中 v~ の速度で弾性衝突をすると考えることができる。(3.40)式中 v に v~ を用
e

いて，衝突後の速度 v e′ を求める。このときの散乱角には(3.48)式を用いるが(3.48)式中のEincident,e には(3.50)式から求められる Escattered,e を用いる。励起エネルギー E excitation は，キセ
ノンの場合 8.34 eV である。
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを Escattered,e ，生成電子のエネル
ギーを Ecreated,e ，電離エネルギーを E ionization とする。エネルギー収支はEscattered,e + E created,e = Eincident,e − Eionization(3.52)と表わされる。
電離後の散乱電子のエネルギーは以下のように求められる(13)。
⎧⎛ Eincident,e − Eionization ⎞⎫⎟⎟⎬Escattered,e = B tan ⎨r5 tan -1 ⎜⎜
2B
⎝
⎠⎭
⎩

(3.53)

ここで，B はガスに固有の値でキセノンでは 8.7 eV となる。
また，入射電子が電離によって失うエネルギーを ∆E とすると，∆E = Eincident,e − Escattered,e(3.54)と表わされる。(3.51)式から電離後の電子の速度は，

v~ = v e 1 −∆EEincident,e(3.55)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える。(3.40)式中 v e に
v~ を用いて，衝突後の速度 v ′ を求める。散乱角には (3.48) 式を用いるが， (3.48) 式中のe

Eincident,e には(3.53)式から求められる Escattered,e を用いる。電離エネルギー E ionization は，キセノ

15

ンの場合 12.13 eV である。
電離後の生成電子の持つエネルギーは，E created,e = ∆E − Eionizationとなる。また，入射電子の速度 v e を用いて生成電子の速度はv created,e =2 E created,eveve(3.56)(3.57)mと表わされる。この速度で弾性衝突をすると考えれば，v created,e を v e の代わりに(3.40)式に代入して，衝突後の速度を求める。このときの散乱角は生成電子のエネルギー E created,e を用いて(3.48)式から求める。
3.2.3.5

Null-collision 法(14) (15)

計算時間短縮のため，Null-collision 法を用いる。これには，架空の衝突断面積 σ fake を導入する。すべての電子のエネルギー E に対して
σ fake < σ total(3.58)であるとし，例えば電子と中性粒子の衝突においてはσ total = σ elastic + σ excitation + σ ionization + σ fake(3.59)とする。ここで，σ total v =υ collisionn gas= constant(3.60)

を満たす σ fake を導入すれば，υ collision が一定，すなわち Pcollision,m が一定となりエネルギーに依存しない。したがって，粒子毎に Pcollision,m を計算する必要が無くなり，全粒子に対して衝突計算を行わず，衝突を起こす粒子のみについて衝突の種類を決めればよい。これ
により，計算時間を大幅に短縮できる。衝突の種類を決めるときは，3. 2. 3. 3 節において， σ total = σ elastic + σ excitation + σ ionization + σ fake とし，∗

σ elastic + σ excitation + σ ionization≤ r2 ≤ 1
σ totalのとき Null-Collision Processとする。Null-Collision Process が選ばれたときには，粒子の運動は何ら変化しないものとする。図 3.3 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度で規格化した衝突周波数を示す。本計算では Null-Collision 法で用いた υ collision ngas の値は 2×10-12

m3/sec である。

16

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

17図 3.2 電子と Xe の衝突断面積

図 3.3 Xe の数密度で規格化した、電子の Xe の衝突断面積

18

3.3

マイクロ波解析手法(16) (17)3.3.1

FDTD 法

3.3.1.1

基礎式

マイクロ波伝播解析には，FDTD 法（Finite Difference Time Domain method）を用いる。

FDTD 法は数値電磁場解析の分野において，モーメント法や有限要素法と同様に，アンテナや平面回路など数多くの複雑な問題に適用され，主要な解析手法として多く用いられている。計算には，マクスウェル方程式を基礎式として用いる。以下に基礎式として用いる微分形式のマクスウェル方程式を示す。

∂H
∂t
∂E
∇× H = ε0+J∂t∇ × E = −µ 0(3.60)(3.61)

ここで， E は電場， H は磁場， µ 0 は真空の透磁率， ε 0 は真空の誘電率， J は電流密度
を表わす。
差分式

3.3.1.2

三次元の FDTD 法で用いられる電場，磁場の空間配置を図 3.4 に示す。このように，
電場，磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる。これは Yee によって提案され(18)，電場の回転（ ∇ × E ）が磁場を，磁場の回転（ ∇ × H ）が電場の方程式を満たすように配置されている。
この電磁場の配置のもと，(3.60)式，(3.61)式を空間，時間とも中心差分した式を以下に示す（(3.62)~(3.67)式）。ここで， ∆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) +∆t
µ 0 (i , j + 1 2 , k + 1 2)⎧⎪ E yn (i , j + 1 2 , k + 1) − E yn (i , j + 1 2 , k ) E zn (i , j , k + 1 2) − E zn (i , j + 1, k + 1 2) ⎫⎪
+
⎬
⎨
∆z
∆y
⎪⎭
⎪⎩
(3.62)
H1n+2y(i + 1 2 , j , k + 1 2) = H

1
n−
2
y

(i + 1 2 , j , k + 1 2) +∆t
µ 0 (i + 1 2 , j , k + 1 2)

⎧ E zn (i + 1, j , k + 1 2) − E zn (i , j , k + 1 2 ) E xn (i + 1 2 , j , k ) − E xn (i + 1 2 , j , k + 1) ⎫
+
⎨
⎬
∆x
∆z
⎩
⎭
(3.63)

19

n+

Hz

1
2

n−

1
2

(i + 1 2 , j + 1 2 , k ) = H z (i + 1 2 , j + 1 2 , k ) +∆tµ 0 (i + 1 2 , j + 1 2 , k )⎧⎪ E xn (i + 1 2 , j + 1, k ) − E xn (i + 1 2 , j , k ) E yn (i , j + 1 2 , k ) − E yn (i + 1, j + 1 2 , k ) ⎫⎪
+
⎨
⎬
∆y
∆x
⎪⎩
⎪⎭
(3.64)
E xn +1 (i + 1 2 , j , k ) = E xn (i + 1 2 , j , k ) +∆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+

+

1

n+

1

H y 2 (i + 1 2 , j , k − 1 2) − H y 2 (i + 1 2 , j , k + 1 2)∆z⎫⎪
− J x (i + 1 2 , j , k )⎬⎪⎭(3.65)

E yn +1 (i , j + 1 2 , k ) = E yn (i , j + 1 2 , k ) +∆t
ε 0 (i , j + 1 2 , k )1n+⎧ n + 122()
(i , j + 1 2 , k − 1 2)
+
+
−
H
i
,
j
1
2
,
k
1
2
H
⎪ x
x
⎨
∆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⎪⎭

(3.66)
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
⎪
⎩
n+

+

Hx

1
2

n+

⎫

1
2

(i , j − 1 2 , k + 1 2) − H x (i , j + 1 2 , k + 1 2) − J (i , j , k + 1 2)⎪
∆y

z

⎬
⎪
⎭

(3.67)(3.62)~(3.67)式から，時間ステップが n + 1 2 での磁場 H
場 E n と１ステップ前の磁場 H

1
n−
2

1
n+
2

は，その半ステップ前での電から求めることができ，時間ステップが n + 1 での電場

20

E n +1 は，半ステップ前の磁場 H

期値については， H

3.3.2

n+

1
2

n+

1
2

と１ステップ前の電場 E n から求めることができる。初は E n から計算されるので，アンテナの根元に E n だけを与える。

安定条件

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

∆t ≤

1
1
1
1
v
+
+
2
2
(∆x ) (∆y ) (∆z )2(3.68)

を満たさなければならない。この条件は極めて厳しく，わずかでも満足しなければ計算が不安定になる。

3.3.3
3.3.3.1

境界条件
PEC(Perfect Electric Conductor )境界条件

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

3.3.3.2

Mur 吸収境界条件(19)

開放領域を取り扱う場合，計算領域を仮想的な境界で閉じなければならない。この境
界を吸収境界といい，この吸収境界が完全でなければ，反射波が解析領域に戻ってしま
い，誤差の原因となる。ここでは，定式化が容易で広く使用されている Mur の吸収境界条件を用いる。
x = 0 の吸収境界に x の正の方向から E z 成分を持つ平面波が入射する場合を考える。伝

播速度を v とすると， x の負の方向に進む波は
E z = E z ( x + vt )(3.69)と表わされ，次の微分方程式を満足する。
∂E z 1 ∂E z−
=0∂x v ∂t(3.70)

x = 0 の境界で完全に反射がないならば，電場は(3.69)式の形を保ったまま伝播するはずである。すなわち，境界でも(3.70)式が満足される。これを FDTD 法の形式に定式化するため，時間について差分化すると

21

∂E z E zn − E zn −1
=
∂t
∆t
n−

=v

∂E z
∂x

(3.71)

12となる。右辺の x についての微分は，電場が x = 0, ∆x に割り当てられていることから，
x = ∆x 2 で行い，(3.70)式に代入すると
n−

1

n−

1

E zn (1 2) − E zn −1 (1 2)E 2 (1) − E z 2 (0)
=v z
∆t
∆x

(3.72)となる。(3.72)式の E zn (1 2), E zn −1 (1 2) も FDTD 法では直接求めることはできないため，前後の値で平均して求める。そして， E z についてまとめるとv∆t − ∆x n
E z (1, j , k + 1 2) − E zn −1 (0 , j , k + 1 2)v∆t + ∆x(3.73)となる。これを Mur の１次吸収境界条件という。上式は吸収境界に平行な電場成分につ{E zn (0, j , k + 1 2 ) = E zn −1 (1, j , k + 1 2 ) +}いて適用される。

3.3.4

給電方法

本計算における給電方法は同軸ケーブルを模擬し，図 3.5 に示すように動径方向の４つの電場を中心導体に与える同軸線路給電を用いる。また，放電室出口には，Mur の一次吸収境界条件を与え，中心導体には PEC 境界条件を与えた。中心導体に与える電場のx 方向成分 E x ，y 方向成分 E y は近似的に
Ex =V⎛ ∆x ⎞
rln⎜ ⎟
⎝ a ⎠

Ey =

V
⎛ ∆y ⎞
rln⎜ ⎟⎝ a ⎠(3.74)(3.75)

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

µ 1
R
ln
ε 2π r

(3.76)

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

22

P=

V2
2Z

(3.77)となり，(3.76)式，(3.77)式より V を求め，(3.74)式，(3.75)式よりアンテナに与える電場を求める。ここで， E x を r = ∆x 2 の電場で， E y を r = ∆y 2 の電場で代表させると，
E x (IFED , JFED , KFED ) = − E x (IFED − 1, JFED , KFED )=−V (n∆t )∆x ⎛ ∆x ⎞ln⎜ ⎟2 ⎝ a ⎠

(3.78)E y (IFED , JFED , KFED ) = − E y (IFED , JFED − 1, KFED )=−V (n∆t )∆y ⎛ ∆y ⎞ln⎜ ⎟2 ⎝ a ⎠

ここで，式中の (IFED , JFED , KFED ) は電場を励振する座標を与える。

23(3.79)

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

図 3.5 同軸線路給電

24

PIC 法と FDTD 法のカップリング

3.4

本研究では，3.2 で述べたプラズマ解析の PIC 法と 3.3 で述べたマイクロ波解析の FDTD
法を組み合わせる事(20) (21)によりプラズマの挙動とマイクロ波の伝播を同時に解析する。
図 3.6 にカップリングコードのフローチャートを示す。FDTD コードで計算されたマイクロ波の電磁場は PIC コードに渡され，その電磁場の効果が粒子の運動に影響する。また，電子の挙動により計算される電流密度が FDTD コードに渡され，電磁波に粒子の効果が加えられる。このように PIC コード，FDTD コードの間で物理的に相互作用を及ぼすデータをやり取りすることで，プラズマおよびマイクロ波の挙動を解析する。なお，
タイムステップについては PIC 法では(3.17)式により決定され，FDTD 法では(3.68)により決定され計算の安定化条件が異なる。必ずしも両者のタイムステップを一致させる必要は無いので，本研究では PIC 法の計算負荷を軽減するために FDTD コードを連続（2~40回）で計算させた後に PIC コードを 1 回計算する事とした。カップリングコードは 3 次
元 X – Y - Z 座標系を採用する事により，3 次元的構造を持つアンテナの解析を可能とした。
データの出力方法について，特に電子エネルギーの空間分布 E e についてはエネルギ
ー分布関数の平均から算出するのではなく，PIC 法のアルゴリズムを応用した方式を採用した。エネルギー分布関数から平均を算出する場合は出力する領域ごとに，全粒子に対して位置に関しての判定を行う必要がある。この方式で E e を出力する領域が多い場合には計算速度の低下につながる。これに対して PIC 法を利用した方式の場合の E e の計算は次式より求まり，E e (m ) = ∑ E ei Ai
i

∑A

i

(3.80)iとなる。ここで Eei は i 番目の電子のエネルギー，ｍは格子点の位置，Ai は(3.21)式で計算される格子点 m への重みを表している。 この方式では全粒子に対して１回だけ，どのセルに含まれるかを判定すれば全ての領域の E e が得られるという利点がある。

25

図 3.6 カップリングフローチャート

26

第4章

4.1

計算結果と考察

小型中和器における電子加速機構

小型中和器では放電室自体が小型なため加速された電子が壁面に衝突し、プラズマ密度の低下を招いていると考えられている。そこで、まず放電室内の少数の電子の軌道を解析した。その後、多粒子での解析を行い、平均電子エネルギー分布、電子エネルギー
分布関数、電離衝突の位置、壁面損失について述べる。
4.1.1 計算体系
図 4.1 に本研究室で使用されている小型中和器を示す。本計算ではこの体系を模擬し
図 4.2 に示すような体系を作成した。実験ではコイルにより磁場を発生させているので磁力線は+Z 方向に一様である。３つの電子の初期エネルギーを 1eV とし、それぞれ初期位置を図 4.2 に示す。また 4.1.5 で述べる多粒子については図 4.2 に示してある領域
1,領域 2( 8 × 8 × 6 mesh)において電子エネルギー分布関数を出力した。

図 4.1

小型中和器

27

図 4.2計算体系

28

計算パラメータと境界条件

4.1.2

計算パラメータを表 4.1 に示す。通常は放電室内には、電子、イオン、中性粒子が存在するが、イオンは電子に比べて非常に重いため、本計算時間内ではイオンはほとんど動かない。よって、本計算では電子のみの挙動を追うものとする。しかし、イオンによって放電室壁面にシースが形成され、電子が放電室内に閉じ込められる。この影響を考慮するため、壁面において高エネルギー電子がシースに入射したとき、シースポテンシャル V を超えたエネルギーであれば壁面損失となり、エネルギーが低ければ反射されるものとした。そのシースポテンシャル V を 20 eV とした。また放電室出口では本来、粒子は外に逃げていくが、今回はプラズマの生成および損失を調べるためであるので出口
では反射とした。
表 4.1 計算パラメータ--------------------------------------------------------------------------------------------------------マイクロ波周波数：2.45 GHz

タイムステップサイズ（マイクロ波）： 5 × 10 −13 sec（プラズマ粒子）： 2 × 10 −11 sec

初期電子エネルギー

：1 eV

推進剤流量

：0.05 sccm

推進剤

：Xe

入力電力

：2 W

メッシュサイズ： 0.408 mm

磁場強度：0.0875 T

超粒子数：21480 個--------------------------------------------------------------------------------------------------------4.1.3

電場強度分布図 4.3(a)~(c)に放電室内の電場強度の時間変化を示す。アンテナ近傍では電場強度は強いがアンテナから離れていくと急激に弱くなっている。放電室内では(a)~(c)の変化が繰
り返し起こっている。

29

(a)

5 × 10 −11 s 後(b)

1 × 10 −10 s 後(c)

1.5 × 10 −10 s 後

図 4.2 放電室内の電場強度の時間変化(a) 5 × 10 −11 s 後, (b) 1 × 10 −10 s 後, (c) 1.5 × 10 −10 s 後30

4.1.4少数の電子の解析結果

図 4.3 に電子 1~3 の YZ 断面の電子軌道、図 4.4 に最もラーマー半径が最大となる電子

1 の XY 断面の電子軌道、図 4.5 にそれぞれのエネルギーの変化を示す。
図 4.3 では電子は+Z 軸方向に移動している。放電室内の電子のラーマー半径は時間によって増大と減少を繰り返している事が分かる。この原因はマイクロ波と電子のサイク
ロトロン運動の位相に起因していると考えられる。全体が ECR 層なので当然、電子のサイクロトロン周波数とマイクロ波の周波数は一致する。そして電子の運動方向と電場の方向が逆であるとき電子は大きなエネルギーを電場から得る。電子は+Z 軸方向に移動するので、サイクロトロン運動とマイクロ波との間で位相のずれが生じる。そして、電子の運動方向と電場の向きが同じ場合に電子は急激にエネルギーを失うことになる。その後、電子が Z 方向に移動し続けることにより位相がずれて加速に適している位相になり電子はマイクロ波の電場によりエネルギーを得ることになる。したがって電子は加速、減速を繰り返すことにより旋回半径は振幅していると言える。電子 1~3 を比べるとアン
テナから遠くなるにつれてラーマー半径の最大値は小さくなっている。また図 4.4 より
電子 1 の最大のラーマー半径は 0.5 mm であり、放電室径に比べて非常に小さい。図 4.5よりも電子エネルギーが上昇と減少を繰り返していることがわかり、その振幅はアンテ
ナからの距離に依存して小さくなっている。

31

図 4.3 電子 1~3 の YZ 断面の電子軌道

図 4.4 電子 1 の XY 断面の電子軌道

32図 4.5 電子 1~3 の電子エネルギーの変化

33

4.1.5 多粒子での解析結果
次に放電室内に多粒子(100 万個)を配置した。図 4.6 に平均電子エネルギー分布、図 4.7に図 4.2 に示した領域 1,2 の電子エネルギー分布関数、図 4.8 に電離衝突した位置、図4.9 に壁面衝突した電子数を示す。
図 4.6 よりアンテナ近傍の電子エネルギーは高いが放電室壁面では電子エネルギーは低い。多粒子の場合でも図 4.5 のように電子のエネルギーはアンテナとの距離に依存していることが分かる。
図 4.7 の電子エネルギー分布関数を見ると領域 1 のアンテナ近傍には領域 2 に比べて電子エネルギーが高い電子が多い。また領域 2 のアンテナから離れた領域であってもキセノンの電離エネルギー12.13 eV 以上のエネルギーをもった電子が存在していることが分かる。
電離衝突は図 4.8 よりアンテナ近傍で多く起こっていることがわかる。電離衝突の位
置は、平均電子エネルギーが高く、高エネルギー電子の存在確率が高い領域と一致する。
しかし、アンテナから離れても数は少ないが電離衝突は起こっていることが分かる。
図 4.9 より壁面にはほとんど衝突していない。これは壁面付近の電子はアンテナからの電場強度が弱いため、エネルギーの上昇が少なく、シースポテンシャルを越えることができないためと考えられる。放電室径 9.8 mm では壁面でエネルギーの損失はなく、放電室径は十分広いと言える。

34

図 4.6 15 nsec 後の電子エネルギー分布

図 4.7 電子エネルギー分布関数35

(a)XY 断面(b)YZ 断面

図 4.8 電離衝突位置(a) XY 断面, (b) YZ 断面36

図 4.9 壁面に衝突した電子数

37

4.2

磁場強度依存性

4.2.1

ECR と 2nd Harmonic ECR との比較(少数の粒子)

図 4.10(a)に ECR での、図 4.10(b)に 2nd Harmonic ECR での電子１の軌道を示す。
図 4.10(a), (b)より ECR では位相のずれによりラーマー半径は時間により増大、減少を繰り返しているが、2nd Harmonic ECR ではラーマー半径は増大し続けている。これは 2nd

Harmonic ECR の方が位相のずれが少ないためであり、本計算では位相が 180 度ずれる前に計算が終了していると考えられる。
図 4.11(a),(b)より 2nd Harmonic ECR に比べ ECR の方がエネルギーの上昇は急である
ことが分かる。この結果より ECR が有効であると考えられる。その理由としては電離衝
突により発生した出口付近の電子は放電室から排出されるまでの時間が短いのでエネ
ルギーが上昇する時間が短い。電離衝突を引き起こすためには、その発生した電子を短時間でエネルギーを上昇させる必要があり、その点において ECR が有効である。

38

(a) ECR(b) 2nd Harmonic ECR

図 4.10 電子１の軌道(a) ECR, (b) 2nd Harmonic ECR

39

(a) ECR(b) 2nd Harmonic ECR

図 4.11 電子１のエネルギー変化(a) ECR, (b) 2nd Harmonic ECR
40

4.2.2

ECR と 2nd Harmonic ECR との比較(多粒子)図 4.12(a),(b)に平均電子エネルギーの分布、図 4.13(a),(b)に電子エネルギー分布関数、
図 4.14(a),(b)に電離衝突発生位置、表 4.1 に電離衝突発生数、図 4.14(a),(b)に壁面に衝突した電子数を示す。
多粒子の場合でも電子エネルギーの上昇は図 4.12 の(a)と(b)を比較すると ECR の方が有効であることが分かる。
図 4.13(a),(b)より領域 1,領域 2 ともに ECR の方が高エネルギーの電子の存在確率が高い。しかし、2nd Harmonic ECR でも領域 1 では電離衝突に必要なエネルギーを持っている電子は存在している。
電離発生位置を示した図 4.14(a),(b)を見ると ECR の方が電離衝突の発生領域が広いことが分かる。この結果より電離衝突には同じ電場強度であっても、ECR にすることで電離発生領域が広がることを確認できた。表 4.1 より電離衝突の発生数を比較すると ECRの方が多く電離衝突が発生している。
図 4.15(a),(b)より 2nd Harmonic ECR の方が電子は壁面に衝突している。2nd Harmonic

ECR では電子エネルギーは低いが、磁場が弱いため式(2.2)よりラーマー半径が 2 倍大きくなるので衝突数が多くなったと考えられる。
損失を抑えるためには軌道を閉じ込めることが重要であったが、この体系ではECR,2nd Harmonic ECR ともに電子軌道は放電室に比べ非常に小さく壁面への衝突による損失は少ないことが分かった。そこで重要なのは生成であり、2nd Harmonic ECR でもプラズマは生成されることが確認できたが、ECR に比べその効果は少ないことが分かった。以上のことより最適な磁場は ECR であることが分かった。

41

(a) ECR(b) 2nd Harmonic ECR

図 4.12 15 nsec 後の平均電子エネルギー分布(a) ECR, (b) 2nd Harmonic ECR

42(a) ECR(b) 2nd Harmonic ECR

図 4.13 15 nsec 後の電子エネルギー分布関数(a) ECR, (b) 2nd Harmonic ECR
43(a) ECR(b) 2nd Harmonic ECR

図 4.14 電離衝突発生位置(a) ECR, (b) 2nd Harmonic ECR

表 4.1 電離衝突発生数
ECR

19770 回

2ndECR

12087 回44

(a) ECR(b) 2nd Harmonic ECR

図 4.15 壁面に衝突した電子数(a) ECR, (b) 2nd Harmonic ECR

45

4.3

放電室径依存性

放電室径による影響を調べるために放電室径をΦ2.4mm~Φ10mm まで変化させて計算を行った。なお、推進剤流量を 0.05 sccm と一定にしている。放電室径を変化させる
ことによる中性粒子密度の変化はコンダクタンスの計算により設定した(22)。また初期粒子数は初期プラズマ密度を 1 × 1016 m-3 とするために各計算で調整している。
図 4.15 に横軸に放電室径、縦軸に壁面に衝突した電子の割合を示す。表 4.2 に 150 nsec
後の放電室内のプラズマ密度、平均の電子エネルギーを示す。図 4.16 に表 4.2 の値より概算した電流値を示す。電流値は式(4.1)を用いた(23)。
J = en p S

kTe1
exp(− )
2
me

(4.1)

e は電子の電荷、 n p は電子密度、 Te は電子エネルギー、 me は電子の質量、k はボルツマ
ン係数、S は放電室径の面積である。
図 4.15 より放電室径が 6.5 mm 以下で損失の割合が増えてはじめていることがわかる。
この結果から 6.5 mm より放電室径が小さくなると、20 eV 以上の電子の壁面付近での存在確率が増加していることが言える。よって壁面付近まで壁面に電子が衝突してしまう
電場強度が行き届いていることが分かる。
表 4.1 より径が小さくなるほどプラズマ密度は高くなっていることが分かる。これは中性粒子密度が高いため電離衝突が多く起こっているためと考えられる。電子エネルギ
ーは 9.8 mm, 2.4 mm では低くなっている。9.8 mm では電場強度が弱く壁面付近の電子エネルギーが低いためであり、2.4 mm では壁面に衝突し、エネルギーを失っているためであると考えられる。
図 4.16 より放電室径が 2.4~6.5 mm までは 4.1 mm がピークになっており、この範囲では 4.1 mm が最適な値であるといえる。放電室径を小さくすることで電流値の上昇が見られるのはプラズマ密度が上昇するためである。しかし、径を小さくすると壁面損失が増え、電子を引き出す面積が減るので、その結果、引き出すことができる電流値が下がる。これらの関係からどこかに最適値があるはずだが、9.8 mm では電流値が大きく上がっている。この計算は電流値に関しては正しく計算されていないといえる。その原因としては計算時間の問題が考えられる。径が小さいほど中性粒子密度が高く、中性粒子密度が高いとプラズマ密度の増加は著しく増加し、2.4 mm では初期条件と比べると約 7.2
倍、9.8 mm では約 1.02 倍となっており、計算時間が増えるとその差はさらに拡大する。
そのため、どの時間まで計算するのかというのが大きな問題となる。もうひとつの問題
として電流を引き出すために電子を加速していないということが考えられる。実際は、電圧をかけて電子を加速しているのでそれによって電子エネルギーが増加する。結果として、引き出される過程で電離衝突が起き、電子電流は増加する。本計算ではその効果を考慮していない。以上の理由より最適な放電室径を求めるために、電流値を見積もる

46

にはさらなる検討が必要であると考えられる。

47

図 4.16 放電室径と衝突した電子数の割合

48表 4.2 150 nsec 後のパラメータ
放電室径（mm）

中性粒子密度(m-3)

プラズマ密度(m-3)

電子エネルギー（eV）

9.8

9.31E+18

1.0229E+16

9.6

6.5

2.55E+19

1.1234E+16

22.22

4.9

5.83E+19

1.71254E+16

31.65

4.1

1.02E+20

3.02768E+16

31.05

2.4

2.06E+20

7.20995E+1618.53

推定電流値(mA)

1
0.8
0.6
0.4
0.2
0
0

2

4
6
放電室径(mm)

図 4.17 推定電流値

49

8

10

第5章

結論

カップリングコードを用いて、マイクロ波放電型小型中和器内におけるプラズマ生成
機構を解析した。その結果を以下に示す。
放電室径 10 mm の小型中和器で計算を行った。放電室全域には十分な電場強度がなく、主にプラズマ生成が起こっている領域はアンテナ近傍であった。ECR ではサイクロトロ
ン運動と電磁波の位相のずれの問題、電場強度が壁面付近では十分ではないため電子の壁面での衝突は非常に少ないことが分かった。その損失は 2nd Harmonic ECR と比べて少なかった。生成の面において、ECR の方がエネルギーの上昇が急であり、有効であることが分かった。以上の結果より最適な磁場は ECR であることが示された。
放電室を小さくしていくことで電子の壁面への衝突が増えていった。しかし、9.8 mm
から 6.5 mm までは電子の壁面への衝突の割合が少ないまま、ほとんど変化しなかった。
この結果よりアンテナ近傍以外ではエネルギーの高い電子の存在確率は少ないことが示された。得られた計算結果から電流値を概算し検討したが、今回の計算では最適値を見つけ出すことができなかった。計算条件の見直し、計算コードの改良などが今後の課題として挙げられる。

51

参考文献
(1) 栗木恭一, 荒川義博: “電気推進ロケット入門”（東京大学出版会, 2003）7.
(2) 片原田弘: 小型マイクロ波放電式イオンスラスタにおける推進性能評価, 九大修士論文, (2002)(3) 重水哲郎、他: 短ギャップ同軸形マイクロ波放電の開始条件に及ぼす磁界の影響, 電学論 A, 111(1991), 727.
(4) 内藤裕志: プラズマ・核融合学会誌, 74 (1998) ,470.
(5) 田中實, 山本良一: “計算物理と計算化学” (海文堂,1998)(6) 石黒静児: プラズマ・核融合学会誌, 74 (1998) , 591.
(7) 丹慶勝市, 他: “Numerical Recipes in C (日本語版)” (技術評論社,1993)(8) M. Hayashi, J. Phys. D: Appl. Phys., 16 (1983) , 581.
(9) M.S. Dababneh, et al.: Phys. Rev. A, 22 (1980) , 1872.
(10) J. J. Szabo, et al.: AIAA 00-3653, (2000)(11) 南部健一: “コンピュータアナリシスシリーズ 7 原子･分子モデルを用いる
数値シミュレーション 日本機会学会編” (コロナ社, 1996)(12) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical Physics, 55(1971) , 4100.
(13) M. Srendra, D. B. Graves and G. M. Jellum, Physical Review A, 41 (1990) ,1112.
(14) V. Vahedi and M. Surendra, Computer Physics Communications ,87 (1995) ,179.
(15) 電気学会・マイクロ波プラズマ調査専門委員会: “マイクロ波プラズマの技術” (オーム社, 2003)(16) 橋本修, 安部琢美: ”FDTD 時間領域差分法入門” (森北出版, 1996)(17) 宇野亨: “FDTD 法による電磁界およびアンテナ解析” (コロナ社, 1998)(18) K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966) , 302.
(19) G. Mur, IEEE Trans. Electromagnetic Compat., 40 (1998) , 100.
(20) V. P. Gropinath and T. A. Grotjohn, IEEE Trans. Plasma Sci., 23 (1995) , 602.
(21) T. A. Grotjohn, Rev. Sci. Instrum. 67 (1995) , 921.
(22) G.F.Weston: “超高真空技術の実際”（共立出版株式会社, 1991）(23) 石川順三: ”荷電粒子ビｰム工学” (コロナ社, 2001)

52