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

修 士 論 文論文名
マイクロ波放電式
小型イオンエンジンの数値解析

氏

名

金川

隆保

指導教員名

中島

秀紀

教授

目次
第1章

序論

1

1.1

研究背景

1

1.2

イオンエンジンの原理

1

1.2.1

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

2

1.2.2

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

4

1.3

過去の研究

6

1.4

研究目的

6

第2章

数値解析

9

2.1

3D simulation

9

2.1.1

粒子シミュレーション

9

2.1.2

電磁波解析手法

19

2.1.3

coupling コード

25

2.1.4

計算条件

27

2D3V コード

30

2.2

2.2.1

PIC 法

30

2.2.2

静電場計算

30

2.2.3

計算条件

32

第3章

数値解析妥当性検証

36

3.1

磁場データ変換36

3.2

PIC コード検証

37

3.3

GMRES 法コード検証

38

計算結果と考察40

3D simulation

40

第4章
4.1

4.1.1

電子軌道

40

4.1.2

∆E の磁場依存性

42

4.1.3

空間電子温度分布

45

4.1.4

EEDF (Electron Energy Distribution Function)

47

4.2

48

4.2.1

電子密度分布48

4.2.2

プラズマポテンシャル分布

49

4.2.3

電離衝突位置50

第5章
謝辞

2D simulation

結論

53
54

第1章
1.1

序論

背景

GPS機能を搭載したカーナビゲーションは以前から存在していたが，最近，GPS機能付き携帯電話も耳にするようになった．GPSとは複数の人工衛星から受信した信号の位相差によって自分の現在の位置を知るためのシステムである．このように多くの人が個人レベルでも人工衛星を使用するようになってきており，宇宙利用が身近なものになってきている．この傾向は加速し，人々の要求を満たすために今後ますます多くの衛星が地球近傍を周回することになるだろう．しかし，大型衛星では開発期間が 5 ~ 10 年以上，打上げコストが数百億 ~ 数千億円であり開発は容易ではない．そのため開発コストの低減・打ち上げ失敗リスクの分散などが可能となる小型衛星の研究が進められている (1-1),(1-2) ．大型衛星の打ち上げは国家規模の大型プロジェク
トとなるが，民間企業・大学レベルで大きさ数十センチメートル程度の小型衛星を作り，それを大型衛星の打ち上げに便乗させる計画が予定されるなど，
小さな団体でも宇宙利用ができるようになってきた (1-3) ．それによりさらに挑戦的なプロジェクトやタイムリーな宇宙開発が行われることが期待される．
人工衛星や宇宙ステーションなどの地球周回物は，地球が楕円球であることや大気抵抗の存在，また太陽や月の引力による摂動があるため，その姿勢
制御・軌道制御に推進機が必要である．現在，衛星の寿命を決定するものは機器の耐年数よりもむしろ推進剤の枯渇である．特に小型衛星においてはサ
イズの制約上，推進剤を大量に搭載することが不可能であり，少量の推進剤
で長期間運用可能な燃費の良い推進機の開発が急がれる．その課題を解決す
る打開策として電気推進の採用が挙げられる．
電気推進

(1-4)

とは，電気エネルギーを推進剤に与えることで反力を得る推進機であり，従来の化学反応を利用したロケット推進と違って静電力や電磁力を利用できる．このため比推力（推進剤流量あたりの推力）を一桁以上大きくすることが可能であり，大幅な推進剤の低減が望める．現在使用されている代表的な電気推進機はアークジェットスラスタやホールスラスタ，そしてイオンエンジンである．イオンエンジンは他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制御や惑星探査に適している．

1.2

イオンエンジンの原理

イオンエンジンの概念図を図 1.1 に示す．イオンエンジンはアーク放電やマイクロ波などで推進剤を加熱・電離させてプラズマを生成し，2 枚ないし3 枚からなる多孔状の電極に 1000 V 程度の高電圧を印加させてイオンを加速するという静電加速型の推進装置である．イオンエンジンは主に 3 つの領域から構成されている．
①

推進剤を電離するイオン生成部（Ionization）②

生成されたイオンを静電的に加速して推力を得る加速部（Acceleration）③

放出されたイオンビームを電気的に中和する中和部 （Neutralization）

これらの各過程はそれぞれイオン源，加速電極，中和器によって行われる．

図 1.1

イオンエンジンの概念図

1.2.1 マイクロ波放電式イオンエンジン
イオンエンジンはプラズマ生成方法により直流放電式やマイクロ波放電式などに分類される．本研究の対象であるマイクロ波放電式イオンエンジンの概念図を図 1.2 に示す．マイクロ波放電とは，マイクロ波帯域の交流電場によって電子を加速し，中性粒子と衝突電離により電子の数が増加し気体がプラズマ化され放電が維持されることである．このマイクロ波放電は直流放電と比べた場合，以下のようなメリットがある．

①プラズマ生成に電極を必要とせず，それによるエンジン寿命の制約がない．
②ヒーターを用いないため，予備加熱が不要となりエンジンを即座に立ち上げることが出来る．
③1 台のマイクロ波電源で複数台のエンジンにプラズマを生成できるため，中和器のプラズマも同時に生成することが可能となり，システム構成の簡素化・信頼性の向上・衛星の重量軽減をもたらす．

図 1.2

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

マイクロ波放電式イオンエンジンはプラズマの生成に電子サイクロトロ
ン共鳴 (Electron Cyclotron Resonance : ECR) を利用し，電子の加熱効率を上げている．
ECR の原理を図 1.3 に示す．真空中に存在する荷電粒子は磁場中でローレンツ力を受ける．この向心力のため磁力線に巻きつくようなサイクロトロン運動と呼ばれる回転運動が現れる．磁界における電子の運動方程式は次式で
表される．

m

dv
= −ev × B
dt

(1.1)

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

mv ⊥
eB

(1.2)

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

ω ce =

eB

(1.3)

m

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

E

B

電子

E

B

電子

電子
E

B
B

電子

図 1.3 ECR の原理

E

1.2.2 マイクロ波放電式小型イオンエンジン
本研究室ではHAYABUSA (1-5) に搭載されたマイクロ波放電式中和器を元にマイクロ波放電式小型イオンエンジンを作成している．図 1.4 にスラスタの写真を，図 1.5 にその構成図を示す．Al製円筒放電室の周囲に 4×4×12mmの長手方向磁化磁石を複数設置し，それらを軟鉄製ヨークで挟み込んで放電室内に磁気回路を形成させている．磁石は永久磁石のSm-Co (サマリウ
ムコバルト) を用い，アンテナ・グリッドはモリブデン (Mo) 製である．図1.6 に実験結果を示す．磁場強度によってイオンビーム電流が変化している．
つまり，最適な磁場強度が存在していることが分かる．しかし，その原因は完全には解明していない．

図 1.4

マイクロ波放電式小型イオンエンジンの写真

Mo

Al

Soft iron

Sm-Co magnet

Propellant gas
Antenna

Front yoke

Microwave

Housing
Back yoke
Discharge chamber
Screen grid

Acceleration grid

図1.5 マイクロ波放電式小型イオンエンジンの構成図

Ion beam current, mA

15
10
Not Ignited.

50

m = 0.018 mg/s
Pi = 8 W
Vb =1500 V
910

11 12 13 14 15
Number of magnets

16

図 1.6 マイクロ波放電式イオンエンジンの実験結果

17

1.3

過去の研究

これまで本研究室でマイクロ波放電式小型イオンエンジンの磁場強度・ア
ンテナ形状・放電室形状などを変化させて実験を行ってきた (1-6),(1-7),(1-8) ．しかし，形状を変更させたことによって，放電室内部の物理メカニズムがどのように変化し，どのような理由で性能に影響をおよぼしているのかは解明されていない．その原因はエンジン自体が小型であるために静電プローブ等の装置を入れることで放電室内部のプラズマに擾乱を引き起こし，正確な測定ができないからである．
そこで数値解析は非常に有効な手段である．近年の計算機の発展のおかげ
で実験では知り得る事のできない情報を入手できるようになってきた．それにより，実験による時間的・金銭的コストを抑えた研究開発ができるようになってきている．本研究室では増井らによって，3 次元数値解析コードが開発されてきた (1-9) ．このコードは電子とマイクロ波との相互作用の解析を目的としており，イオンがほとんど動かないタイムスケールであるためイオンを背景として扱っている．そのため定常状態の解明には至っていない．そこでタイムスケールを大きくとり，定常状態までの計算を行うことが次の課題として挙げられる．

1.4

研究の目的

以上の背景より，本研究の目的は，マイクロ波放電式小型イオンエンジン
内部のプラズマ状態を模擬できる数値解析コードの開発及びそのコードを用いたプラズマの観察である．
参考文献(1-1) M.Kato, S. Takayama, U. Nakamura, K. Yoshihara and H. Hashimoto, 56 thInternational Astronautics Congress paper IAC-05.B5.6.B.01, 2005.
(1-2) H. Sahara, S. Nakasuka and C. Kobayashi, AIAA paper 2005-3956, 2005.
(1-3)http://www.astroarts.co.jp/news/2007/05/17small_satellite/index-j.shtml.
(1-4) 栗木恭一・荒川義博[編]：”電気推進ロケット入門”, (東京大学出版会,2003)(1-5) I. Funaki, H. Kuninaka and K. Toki, J. Propul. Power, Vol. 20, (2004), pp.718-726.
(1-6) N. Yamamoto, H. Kataharada, H. Masui, H. Ijiri and H. Nakashima,
AJCPP2005-22093, 2005.
(1-7) N. Yamamoto, H. Kataharada, T. Chikaoka, H.Masui and H. Nakashima,
IEPC-2005-036, 2005.
(1-8) T. Chikaoka, S. Kondo, N. Yamamoto, H. Nakashima and Y. Takao, Proceedings ofthe 25th International Symposium on Space Technology and Science, 2006,
pp.254-259.
(1-9) 増井博一， “マイクロ波放電型宇宙推進用プラズマ源に関する数値解析”, 平
成 17 年度九州大学博士論文．

第2章

数値解析

2.1 3D simulation
2.1.1 粒子シミュレーション (2-1)
2.1.1.1 基礎式
プラズマの挙動を解析する際の基礎式として運動方程式を用いる．式を以
下に示す．

m

dv
= q(E + v × B )
dt

(2.1)

dx
=v
dt

(2.2)

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

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

(2.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 ⎬
m⎩
2
∆t
⎭

(

)

(2.4)

である．n- 1 2 はある時刻 n から ∆t 2 だけ以前の値であることを示す．同様に，n + 1 2
はある時刻 n から ∆t 2 だけ進んだ値であることを示す．
プラズマ粒子の運動を求めるために，(2.4)式において未知量である v
き，求めた vn +1 2と(2.3)式より新たな粒子の位置 x

n +1

n +1 2について解を求めることができる．(2.4)式

を変形すると，

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

(2.5)

α=

q∆t n
Bx
2m

(2.6)

β=

q∆t n
By
2m

(2.7)

γ =

q∆t n
Bz
2m

(2.8)

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

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

v n +1 2 = A −1 S

⎛ 1+α 21⎜
A −1 = ⎜ αβ − γ
δ⎜⎝ αγ + β2.1.1.2

αβ + γ αγ − β ⎞⎟1 + β 2 βγ + α ⎟βγ − α 1 + γ 2 ⎟⎠(2.12)(2.13)

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

計算に用いるメッシュサイズ ∆ は，電子温度 T e ，ボルツマン定数 k b ，電子密度 n eから計算される．
デバイ長 λ De ，λ De =

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

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

(2.17)とすれば，数値的な不安定性は生じない．時間積分の精度まで考えると，通常は∆t ≤ 0.2ω pe−1

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

(2.19)

に対して，∆t << ω c

−1

(2.20)

である必要がある．

2.1.1.3

PIC法 (2-4)

シミュレーションでは，電磁場量を空間格子点上で求める．そのため，計算粒子の電荷，電流を空間の格子点上に分配し，マクスウェル方程式を用いて電磁場を計算しなければならない．また，粒子が感じる電磁場を計算する
際，格子点上に割り振られた電磁場を内挿する必要がある．このような計算を行うため，粒子シミュレーションにおいて標準的な手法である 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.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)
(2.22)と表わされる．

図 2.1

2.1.1.4

PIC 法におけるセルの構成

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

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

V1
qN
V2

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

j = qnvより，

j (i , j , k ) =2.1.1.5

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

V1
qNv
V2

(2.24)(2.25)

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

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

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

λ = ∫ lP(l )dl0(2.32)

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

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

l = −λln (1 − r )となる．ここで， r は一様乱数であるから 1 − r を r とすれば，(2.33)(2.34)

l = −λln(r )(2.35)と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を σ とすれば，この粒子が１秒間に気体分子と衝突する回数 υ は，

υ=

v

λ
= σNv

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

2.1.1.7

衝突確率，衝突判定

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

σ = σ (E )(2.37)

荷電粒子速度がわかれば，荷電粒子と中性粒子の各衝突における衝突周波数がわかり，全衝突周波数 υ total がわかる．したがって，ある粒子 m の衝突確率

Pcollision,m は以下の式で表わすことができる．
υ total = ngasσ total v m

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

2.1.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 である．このように断面積の比と乱数を対応させ，衝突の種類を決める．計算では中性粒子としてキセノンを用い，その断面積データ (2-7), (2-8), (2-9) を図 2.2 に示す．今回の計算では，多価電離は考慮していない．

図 2.2 電子と Xe の衝突断面積
2.1.1.9

エネルギーと散乱角

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

v1′ =

1
(v1 + v 2 − u′)2(2.38)

v ′2 =

1(v1 + v 2 + u′)2(2.39)と表わすことができる．ここで， u′ = v ′2 − v1′ （相対速度）である．電子と中性粒子の衝突も同様に考える．衝突前の電子，中性粒子の速度をそれぞれ

v e ・ v n の衝突後の速度をそれぞれ v e′ ・ v ′n とする．電子と中性粒子の衝突では，電子は中性粒子に比べ質量が非常に小さいため，電子のみが散乱される．
中性粒子の速度 v n を無視すれば，相対速度は u = v n − v e ≅ −v e としてよい．よって，電子の質量を m ，中性粒子の質量を M とすれば，次式が成立する (10) ．

v e′ =

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

( 2.40)h の直交座標成分は

hx = ver cos φ
hy = −
hz = −(2.41)

vex vey cos φ − ve vez sin φ(2.42)vervex vez cos φ − ve vey sin φ(2.43)ver

ここで，ver = vey2 + vez2

(2.44)

ve = vex2 + vey2 + vez2(2.45)

であり， χ は散乱角，φ は任意の角である．(2.40)式はエネルギー損失を含んでいる．

次に，衝突前後のエネルギー収支を考える．衝突前の電子のエネルギーを
Eincident,e とし，衝突後の電子のエネルギーを Escattered,e とすれば，エネルギー収支

は以下のように表わされる．

Escattered,e = Eincident,e(2.46)

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

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

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

2 1 − (1 + Eincident, e ) 3
Eincident, er}(2.48)となる．この値を (2.40) 式に用いて衝突後の速度を求める．このとき，任意の角は

φ = 2π r4(2.49)

として与える．

②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネ

ルギーを E excitation とし，入射電子のエネルギーを Eincident,e ，励起後のエネルギーをE，励起後の速度を v~ とする．エネルギー収支は，
scattered,eEscattered,e = Eincident,e − E excitation(2.50)と表わされ，励起後の速度はE
v~ = v e 1 − excitation
Eincident,e(2.51)となる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失った後に(2.51)式中 v~ の速度で弾性衝突をすると考えることができる．(2.40)式中 v e に v~を用いて，衝突後の速度 v e′ を求める．このときの散乱角には(2.48)式を用いるが(2.48)
式 中 の Eincident,e に は (2.50) 式 か ら 求 め ら れ る E scattered,e を 用 い る ． 励 起 エ ネ ル ギ ー

E excitation は，キセノンの場合 8.34 eV を用いている．
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを Escattered,e ，生成電子のエネルギーを E created,e ，電離エネルギーを Eionization とする．エネルギー収支
はEscattered,e + E created,e = Eincident,e − Eionization(2.52)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる (2-12) ．

⎧⎛ Eincident,e − Eionization ⎞⎫⎟⎟⎬Escattered,e = B tan ⎨r5 tan -1 ⎜⎜
2B
⎝
⎠⎭
⎩

(2.53)

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

v~ = v e 1 −∆E(2.55)

Eincident,eとなり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．(2.40)式中 v eに v~ を用いて，衝突後の速度 v e′ を求める．散乱角には(2.48)式を用いるが，(2.48)式中の Eincident,e には(2.53)式から求められる E scattered, e を用いる．電離エネルギー E ionization は，キセノンの場合 12.13 eV である．

電離後の生成電子の持つエネルギーは，E created,e = ∆E − Eionization(2.56)

となる．また，入射電子の速度 v e を用いて生成電子の速度はv created,e =veve

2 E created,e(2.57)mと表わされる．この速度で弾性衝突をすると考えれば， v created,e を v e の代わりに (2.40) 式に代入して，衝突後の速度を求める．このときの散乱角は生成電子のエネルギー E created,e を用いて (2.48) 式から求める．
Null-collision法 (2-13),(2-14),(2-15)

2.1.1.10

計算時間短縮のため，Null-collision 法を用いる．これには，架空の衝突断面積 σ fake を導入する．すべての電子のエネルギー E に対して
σ fake < σ total(2.58)であるとし，例えば電子と中性粒子の衝突においてはσ total = σ elastic + σ excitation + σ ionization + σ fake(2.59)とする．ここで，σ total v =υ collisionn gas(2.60)= constantを満たす σ fake を導入すれば，υ collision が一定，すなわち Pcollision,m が一定となりエネルギーに依存しない．したがって，粒子毎に Pcollision,m を計算する必要が無
くなり，全粒子に対して衝突計算を行わず，衝突を起こす粒子のみについて
衝突の種類を決めればよい．これにより，計算時間を大幅に短縮できる．衝突 の 種 類 を 決 め る と き は ， 2.1.1.8に お い て ，

σ total = σ elastic + σ excitation + σ ionization + σ fake とし，∗

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

υ collision ngasの 値 は 2×

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

2.1.2 電磁波解析手法 (2-16), (2-17)
2.1.2.1

基礎式

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

∂H
∂t
∂E
∇× H = ε0+J∂t∇ × E = −µ 0(2.61)(2.62)

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

2.1.2.2

差分式

三次元のFDTD法で用いられる電場，磁場の空間配置を図 2.4 に示す．このように，
電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．これはYeeによって提案され (2-18) ，電場の回転（ ∇ × E ）が磁場を，磁場の回転（ ∇ × H ）が電場の方程式を満たすように配置されている．

この電磁場の配置のもと，(2.61)式，(2.62)式を空間，時間とも中心差分した式を以下に示す．（(2.63)~(2.68)式）ここで，∆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) +∆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) ⎫⎪
+
⎬
⎨
∆y
∆z
⎪⎭
⎪⎩
(2.63)
n+

1

n−

1

H y 2 (i + 1 2 , j , k + 1 2) = H y 2 (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) ⎫
+
⎨
⎬
∆z
∆x
⎩
⎭
n+

1

n−

1

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

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

E yn +1 (i , j + 1 2 , k ) = E yn (i , j + 1 2 , 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
⎪
⎭

(2.67)

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

+

Hx

1
2

n+

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

nの電場 E と１ステップ前の磁場 H
の電場 E

n−

1
2

だけを与える．

1
n+
2

n+

1
2

⎬
⎪
⎭
(2.68)

は，その半ステップ前でから求めることができ，時間ステップが n + 1 では，半ステップ前の磁場 Hできる．初期値については，Hz∆y(2.63)~(2.68)式から，時間ステップが n + 1 2 での磁場 H
n +1

⎫

1
2

n+

1
2

と１ステップ前の電場 E n から求めることがは E n から計算されるので，アンテナの動径方向に E n

図 2.42.1.2.3

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

アンテナ解析

本計算ではアンテナを完全導体で構成されているものとしてモデル化する．給電方
法は同軸ケーブルを模擬し，図 2.5 に示すように動径方向の４つの電場を中心導体に与える同軸線路給電を用い，アンテナ根元には後述する Mur の一次吸収境界条件を与える．波長に比べて太さが無視できるくらい十分細い直線状のアンテナを解析する場
合，図 2.6 に示すようにその中心軸に沿って電界を 0 とすればよい．また，今回は円盤アンテナを模擬し，円盤の表面の電界を 0 とした．

図 2.5

同軸線路給電図 2.62.1.2.4

太さ零のアンテナモデル

安定条件

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

∆t ≤

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

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

PEC(Perfect Electric Conductor )境界条件

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

2.1.2.6

Mur吸収境界条件 (2-19)

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

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

E z = E z ( x + vt )と表わされ，次の微分方程式を満足する．

∂E z 1 ∂E z
=0−∂x v ∂t(2.71)

x = 0 の境界で完全に反射がないならば，電場は(2.70)式の形を保ったまま伝播するはずである．すなわち，境界でも(2.71)式が満足される．これを FDTD 法の形式に定式化するため，時間について差分化すると∂E z E zn − E zn −1
=
∂t
∆t
n−

(2.72)

1

∂E 2
=v z
∂x

となる．右辺の x についての微分は，電場が x = 0, ∆x に割り当てられていることから，

x = ∆x 2 で行い，(2.71)式に代入すると

E (1 2) − E∆tnzとなる。(2.72)式の E z (1 2 ), E z
n

n −1

n −1
z

n−

1
2

n−

1
2

(1 2) = v E z (1) − E z (0)(2.72)∆x(1 2) も FDTD 法では直接求めることはできないため，前後の値で平均して求める．そして， E z についてまとめると

E zn (0, j , k + 1 2) = E zn −1 (1, j , k + 1 2) +v∆t − ∆x n
E z (1, j , k + 1 2 ) − E zn −1 (0 , j , k + 1 2)v∆t + ∆x{}(2.73)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．
2.1.2.7

ポインティングベクトル (2-20)

空間を流れる電力を求めるためには，ポインティングベクトルを使用する．Maxwellの方程式を変形すると，以下の式となる．

∂ ⎛ εE 2 µH 2 ⎞⎟⎟dv+− ∫ ( E × H ) ⋅ dS = ∫ J ⋅ Edv + ∫ ⎜⎜
t
∂
2
2
⎠
s
vol
vol⎝

(2.74)

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

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

図 2.7

Coupling code のフローチャート

2.1.4 計算条件
2.1.4.1 計算体系
図 2.6 は当研究室で開発されている小型マイクロ波放電型イオンエンジンの断面図であり，赤の網掛けの部分が計算領域である．

図 2.6
2.1.4.2

計算体系

計算格子

PIC法を用いるために格子を刻む必要があるが，そのために望ましい長さは，デバ
イ長である．そうすることによって，壁面でのシースが計算できる．しかし，3 次元座標系では計算コストの点からデバイ長で格子を刻むことは現時点では不可能である．
今回は電子とマイクロ波との相互作用を調べることが目的であるからデバイ長で計算格子を作る必要はない．そこで本計算では計算格子の長さを 5.0×10 -4 mとした．
2.1.4.3

時間刻み

PIC 法 の 時 間 刻 み は ∆t pic = 1.0 × 10

−12sec 、 FDTD 法 の 時 間 刻 み は∆t FDTD = 5.0 × 10 −13 sec とすれば、2.1.1.2 や 2.1.2.3 で述べた安定条件を満たす．FDTD
コードが 2 回計算された後，PIC コードを 1 回計算させた．
2.1.4.4 電子初期条件
電子は放電室内に 0 eV で一様に分布させた。

2.1.4.5 粒子数
1 ボクセル当たり数個しか粒子がない場合，ノイズが大きくなり望ましくない (2-14) ．
通常，10 から 100 の間の粒子が存在しているとよいとされている．今回は 1 ボクセ
ルに平均 20 個の粒子が含まれるように粒子数を設定した．
2.1.4.6 磁場強度
放電室外側の磁石の数を 8 個・9 個･10 個･11 個と変更し計算を行った．その放電
室 内 の 磁 場 強 度 の 変 化 の 様 子 を 図 2.7 に 示 す ． 磁 場 計 算 に は Advanced Science
Laboratory, Incの 3 次元静磁場解析ソフトを使用した (2-23) ．マイクロ波電源として周
波数が 2.45 GHzのものを使用している．マイクロ波周波数と電子サイクロトロン周波数が一致するときの磁場強度は 87.5 mTであり，その値を線でつなぎ，それをECR層と呼ぶ．図より，磁石数が増えるにしたがってECR層とアンテナとの距離が離れていくことがわかる．

(a)

(b)

(c)

(d)

図 2.7

磁場強度分布(a), N mag =8, (b) N mag =9, (c) N mag =10, (d) N mag =11.

2.2

2D3V コード

放電室内のプラズマを定常状態に達するまで計算するためには，プラズマの生成と損失について考慮しなければならない．生成については上述の coupling code を用いて解析ができる．損失についてはシース近辺の静電場を解く必要があり，そのためにはイオンを粒子として扱わなければならない．3 次元では計算コストが大きすぎるた
め，計算が不可能だった．そこで計算体系の次元を 2 次元とし，マイクロ波伝播の計算は行わず，マイクロ波による電子の ECR 層での加熱は 2.1 で述べた 3 次元数値解
析から得られた値を元にモデル化した．そのコードを 2D3V コードと呼ぶ．2D3V コ
ードとは各粒子は Z ・ R ・ θ の三つの速度成分を持つが，位置の情報は θ 成分を無視
し Z ・ R の二成分としたものである．すなわち粒子は全方向に動くが，位置の θ 成分
は常に無視される． θ 方向一周分の粒子をひとつの ZR 平面に投影するイメージであ
る．

2.2.1

PIC 法

すでに 2.1.1 に述べたように，粒子の計算では PIC 法を用いる．円筒座標系の電磁場中の運動方程式は(2.75)式の通りであり，ここでは四次精度 Runge-Kutta 法にて解
いた．
2
qj
dvr
vθ
=
( Er + vθ Bz ) +
dt
mj
r

q
dvθ
vv
= j ( Eθ + vz Br − vr Bz ) − r θdtmjr(2.75)

qjdvz=( Ez − rθ Br )dtmj
jは粒子種をあらわし，j

=

e,

i

であり， q e =-e， q i =e である．また，衝突計算には同じくNull-Collision法を採用している．電子―中性粒子の弾性・励起・電離
衝突のみ考慮している．

2.2.2

静電場計算

静電場を解析するためにポアソン方程式を解いた．円筒座標での基礎方程式を以下に示す．

∂ 2φ 1 ∂ ⎡ ∂φ ⎤e+r
= − (ni − ne )
2
⎢
⎥
∂z
ε0
r ∂r ⎣ ∂r ⎦(2.76)

φは空間電位，ε 0 は真空中の誘電率，n i と n e はイオン・電子の数密度である．場の方
程 式 を 解 く に あ た り ， GMRes(Generalized Minimal Residual method) 法(2-24),(2-25),(2-26),(2-27) を採用した．

この方法は非対称行列を有する線型方程式の反復解法であり，次のような考え方にもとづいている．解くべき方程式はAx = b

(2.77)

であるとする．ここにAは(N×N 行列)であり，簡単のために実行列であるとする．b は既知の N ベクトル， x は未知の N ベクトルである．いま，適当に初期過程 x 0 を選ぶと，(2.77)の残差 r 0 は

r 0 = b - Ax 0

(2.78)と書ける．いま，第 n 回の繰り返しにおける解 x n を x n = x 0 + z n と書くことにする．ここに z n は，Krylov空間

K n (A,r 0 )=Span(r 0 ,Ar 0 ,,A n-1 r 0 )…

(2.79)から||b - Ax n || = ||r 0 - A z n ||(2.80)

を最小にするように選ぶものとする．記号 Span()は，括弧内のベクトルの線形結合で作られる空間を表す．
(2.80)の最小化を行うために， K n (A, r 0 )の正規直交基底(v 1 , v 2 ,…, v n )を次のようにして逐次構成する．まず||r 0 || = γとすると v 1 = r 0 /γとなる．いま，第 n 番目の基底 v nまでが求められているとする．このとき，

w n =Av n ,h in =v iwn・(2.81)とおき，Schmidt の直交化により

V n+1 = w n - Σ h in v i(2.82)とおけば， V n+1 ⊥ K n (A, r 0 )である．さらに V n+1 = V n+1 /||V n+1 ||とすると，

h n+1

n

:=

v n+1 ・ w n=

v n+1 ・ V n+1=||V n+1 ||(2.83)
を得る．
次に，(2.80)における最小化を実行する．第 n 回目の繰り返しにおいては z n ∈ K n (A, r 0 )であるから，

z n = Σv j y j(2.84)と書ける．これを用いて(2.80)を v 1 ,…, v n+1 座標系で書くと，min || γδ i1 – Σ h ij y j ||(i = 1,2,…, n+1)が得られる．上式の行列 h ij は( n+1)×n の上Hessenberg行列である．

(2.85)

⎛ h11
⎜
⎜ h21
⎜ 0
n⎜
Q
⎜ M
⎜ 0
⎜
⎜ 0
⎝

h12

h13

L

h22

h23

L

h32

h33

L

O
L

O
0

O
hn −1 nL00h1n ⎞ ⎛ h11⎟ ⎜h2 n ⎟ ⎜ 0
h3n ⎟ ⎜ 0⎟=⎜
M ⎟ ⎜ Mhnn ⎟⎟ ⎜⎜ 0
hn +1 n ⎟⎠ ⎜⎝ 0

h12
h22
0
O
L
L

L h1n ⎞⎟
h23 L h2 n ⎟h33 L h3n ⎟⎟O O M ⎟⎟
0 0 hnn ⎟
0 0
0 ⎟⎠
h13

(2.86)とするような( n +1)×(n +1)直交行列 Q n が簡単に構成できるので，(2.85)の解は線形方程
式

⎛ h11
⎜
⎜ 0
⎜ 0
⎜
⎜ M
⎜
⎜ 0
⎜ 0
⎝

h12
h22
0
O
L
L

L h1n ⎞⎟⎛ y1 ⎞⎛ Q11n ⎞⎜ n⎟
h23 L h2 n ⎟⎜ ⎟⎜ Q21 ⎟
⎜ y2 ⎟⎟
h33 L h3n ⎟⎜ ⎟y 3 = γ ⎜ Q31n ⎟⎟⎜O O M ⎟⎜ ⎟
⎜
M
M ⎟
⎟⎜ ⎟
0 0 hnn ⎟⎜ ⎟⎜ n⎟
⎝ yn ⎠⎝ Qn1 ⎠
0 0
0 ⎟⎠
h13

(2.87)

を後退代入によって解くことにより構成され，(2.85)の残差 E n は|γQ n n+1 1 |と求められる．以上のプロセスを残差が許容限度εに対して E n <εを満たすまで繰り返して，近似解 x n を得る．以上が，GMRESの基本的なアルゴリズムである．
電場は空間電位の勾配より求めた．ただし，円周方向成分については無視した．

E = −∇φ

2.2.3(2.88)

計算条件

2.2.3.1 計算領域
フロントヨークからバックヨークにかけての磁気チューブ内でのプラズマ挙動，およびその両端のシースポテンシャルが放電室内部の物理現象に深く関与するものと考え，図 2.8 での赤枠内を計算領域とした．格子数は全体で 92×78 にとった．格子が粗いとプラズマシースにおける非中性の現象を正しく追うことができないため，格子間隔を Debye 長と同程度にする必要がある．等間隔の格子点が望ましいが，計算コストが大きくなりすぎるため，式(2.89)を用いて不等間隔に格子を作成した．

xi exp(aI ) − exp(a ( I − i ))=exp(aI ) − 1
L(i = 0,1,L , I )（2.89）

Lは体系長さ，Iは分割数、x 0 , x 1 , …, x I はセルの境界の座標，aは格子の粗さを調節する係数である．今回は特に重要な磁気チューブ中のアンテナ前面・上部・フロント
ヨークの付近とグリッド付近ではDebye長と同程度にした．式(2.89)を用いて作成した計算格子を図 2.9 に示す．

図 2.8

計算体系概念図

図 2.9

計算格子

2.2.3.2 境界条件
粒子に関しては,軸に達したものはその反対側でも同様に粒子が存在するものと考え，鏡面反射させた．絶縁されているアンテナに衝突した粒子は計算体系から除外す
るが，電荷は残るとした．計算領域右端はグリッドがあるので，電子はグリッドの開口率で反射，イオンに関してはすべて吸収とした．ヨークに衝突したものと計算領域の上部，左端に達したものについては計算から除外した．場については軸ではノイマン境界条件とし，そのほかではディリクレ境界条件を用い 0 V と設定した．
2.2.3.3 質量比の変更 (2-28)
今回の計算において，電子の速度成分は最大で 10 5 m/s程度，イオンにおいては最大 10 3 m/s程度であり，電子の速度を基準に時間刻みを設定するとイオンはほとんど動かない．それを避けるため，電子の質量を実際よりも大きくし，電子とイオンの速度差を小さくした．一般的には逆に重粒子の質量を小さくする方法がとられることが多いが，その場合プラズマ密度・平均自由行程・推進剤流量・衝突断面積など補正すべき物理量が非常に多くなってしまう．これに対し今回用いた電子を重くする方法では，プラズマ密度は速度の遅いイオンによって支配されると考えられるため補正の必要は無く，補正が必要なパラメタはLarmor運動及び衝突周波数のみであり扱いが簡単である．電子を実際よりどれだけ重くするかを表すために．重み係数 α を導入する．

me _ verchal = α me _ real
重み係数は，電子の速度には平方根(2.90)

α で効いてくるので，時間刻みを α 倍粗くすることができる．今回の計算では α =2500 とした．これにより実際は 1:2.4×10 5である電子とXeイオンとの質量比は 1:96 にまで近づく．重み係数の導入により変化する物理量を，以下のように補正する．
①Larmor 運動
Larmor 半径、Larmor 周波数はそれぞれ次のように書ける．

rc =

mv
eB

(2.91)eBm(2.92)

ωc =イオンについては重み係数の影響は無い．電子の Larmor 半径を補正するため，電子の感じる磁場の大きさを

α 倍大きくする．このとき Larmor 周波数は α 倍小さ

くなるが，Larmor 半径に比べその影響は小さいと考えられるのでここでは無視する．
②衝突断面積
衝突断面積は電子速度に依存する．電子の速度はる影響を，衝突断面積を1α倍となるので，その小さくな

α 倍大きくすることで補正する．

参考文献(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)(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) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical
Physics, 55 (1971), 4100.
(2-12) M. Srendra, D. B. Graves and G. M. Jellum, Physical Review A, 41 (1990),
1112.
(2-13) V. Vahedi and M. Surendra, Computer Physics Communications, 87 (1995),
179.
(2-14) 電気学会・マイクロ波プラズマ調査専門委員会：“マイクロ波プラズマの技術”

(オーム社, 2003)(2-15) K. Koura, Phys. Fluids, 29(1986), 11.
(2-16) 橋本修, 安部琢美：“FDTD 時間領域差分法入門” (森北出版, 1996)(2-17) 宇野亨：“FDTD 法による電磁界およびアンテナ解析” (コロナ社, 1998)(2-18) K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.
(2-19) G. Mur, IEEE Trans. Electromagnetic Compat., 40 (1998), 100.
(2-20) 著=W. H Hayt，訳=山中惣之助，岡本考太郎，宇佐美興一：“工学系の基礎電磁気学”(朝倉書店, 2001)(2-21) V. P. Gropinath and T. A. Grotjohn, IEEE Trans. Plasma Sci., 23, (1995),602.
(2-22) T. A. Grotjohn, Rev. Sci. Instrum. 67, (1995), 921.
(2-23) http://www.asl-i.com/index.html(2-24) Y. Saad, and M.H. Schultz , SIAM. J. Sci. Stat. Comput., 7, pp.856-869(1986)(2-25) 藤野清次，張紹良：”反復法の数理”, (朝倉書店, 1996)(2-26) 樫山和男，西村直志，牛島省：”並列計算法入門”, (丸善, 2003)(2-27) 牛島省：”OpenMP による並列プログラミングと数値計算法”, (丸善, 2006)(2-28) M. Hirakawa, IEPC-97-021, (1997), 159.

第3章

数値解析妥当性検証

2D3V-DSMC-PIC コードは今年度，新しく開発したコードである．新しく作ったプ
ログラムにはバグはつきものであるので，正しくプログラムできているかを確認した．
3.1

磁場データ変換

磁場計算ソフトでは等間隔格子上でしか磁場データを出すことができない．しかし，計算では図 2.9 に示すような格子点上に磁場のデータが必要となる．そこで不等間隔格子に応じてデータを線形補間する．その結果を図 3.1 に示す．二つを比較し等高線
は概ねあっており，等間隔格子上のデータを不等間隔格子上に配置することができた．

(a)

(b)
図 3.1
(a) 等間隔格子上磁場データ,

磁場データ(b) 不等間隔格子上磁場データ

3.2

PIC コード検証

不等間隔格子・2D3V という条件で PIC が正しく動いているかを確認した．図 3.2に電子軌道，図 3.3 にエネルギーの履歴の結果を示す．電子はラーマ運動をしながら磁力線に沿って移動していることがわかる．また電子のエネルギーも計算中，一定に保たれており電子の運動方程式は正しくプログラムされていることが確認された．

図 3.2 磁場中の電子軌道

図 3.3

電子エネルギー履歴3.3 GMRES 法コード検証
GMRES 法の妥当性を検証するため，理論解と比較した．計算体系として図 3.4 のような一辺の長さが 1 の正方形の領域を扱い，下記の境界条件を与えた．
φ(x,0) = sin (πx), φ(x,1) = φ(0,y) = φ(1,y) = 0
このときの理論解は次式で与えられる．

ϕ ( x, y ) =sin(πx) sinh[π (1 − y )]sinh(π )

Ｙ
Ｘ
図 3.4

計算体系

図 3.5(a)が理論解，(b)が GMRES 法を用いて解いたときの計算結果である．両者に見た目で差はない．そこで理論解と GMRES 法の差をとったものが図 3.6 である．誤差はあるものの反復法である以上これは避けられない．結果として，良い一致を示しており，これでポアソン方程式を GMRES 法を用いて正しく解けていることが確認され
た．

(a)

(b)

図 3.5

理論解と GMRES 法の比較

(a) 理論解, (b) GMRES 法図 3.6

誤差

第4章
4.1

計算結果と考察3D simulation

4.1.1 電子軌道
放電室内の一つの電子に注目した．その軌道の結果を図 4.1 に，電子のエネルギー，
感じる磁場を図 4.2 に示す．図 4.1 より電子は磁力線に沿って移動していることがわ
かる．図 4.2 より電子はフロントヨーク側に移動し，そこで磁場により反射されている．バックヨーク側に移動し，電子の感じる磁場が 87.5 mT 前後の時，ECR によって加熱されていることがわかる．ここで電子の ECR を通過する前と後のエネルギーの変化をΔE とする．図 4.3 に示すようなイメージである．
フロントヨーク側の ECR ではほとんど加熱されていないのは，アンテナからの電界強度が十分でないためと考えられる．

図 4.1 電子軌道

図 4.2 電子のエネルギーと感じる磁場の時間履歴E

電子
図 4.3 ΔE の概念図

4.1.2 ⊿E の磁場依存性
ΔE の磁場依存性を調べた．解析領域を図 4.4 に示すように 3 つに分けた．この理由として，ECR とアンテナとの距離だけでなく，位置関係によってもΔE が変化すると考えたからである．

図 4.4

解析領域

結果を図 4.5〜4.8 に示す．ヒストグラムの幅は 0.1 eV となっている．図 4.5 より 1と 2 に比べて 3 の領域ではエネルギーがあまり変化していないことがわかる．これは図 4.9 に示すとおり，アンテナの端の方が真ん中よりも電界強度が強くなるためである．そのためアンテナ正面の 3 の領域では電界強度が弱く，エネルギーの変化が少な
い．1 と 2 ではほぼ同じような結果となった．また磁気チューブを含む 1 と 2 で最大のエネルギー変化でも 10 eV 程度である．これは電子は 1 回の ECR の通過ではキセ
ノンを電離させるエネルギー(12.13 eV)を得ることができないということである．電子がキセノンを電離させるためには複数回磁気チューブ間を往復する必要がある．そのためフロントヨーク・バックヨーク付近での電子の振る舞いが重要となってくる．
磁石数が増えるにしたがってΔE の絶対値は小さくなっている．これはアンテナから離れているため，アンテナからの電界の供給が弱くなっているため当然の結果といえる．
今回はシースを境界条件として扱っているが，ヨーク付近での物理を正確に再現する必要があるため，ポアソン方程式によりシースを自己矛盾なく解くコードの開発が望まれる．今回の計算結果である電子とマイクロ波との相互作用ΔE と 2.2 で述べた
コードとをカップリングさせることによって電磁波解析をする必要がなく，計算コストを抑えることができる．

(a)

(b)

(a)

(b)

(c)

(c)
図 4.5

8個

(a)領域 1, (b)領域 2, (c)領域 3

図 4.6

9個

(a)領域 1, (b)領域 2, (c)領域 3

(a)

(a)

(b)

(b)

(a)
(c)
図 4.7

(c)
10 個

(a)領域 1, (b)領域 2, (c)領域 3

図 4.8

11 個

(a)領域 1, (b)領域 2, (c)領域 3

図 4.9 円盤アンテナ正面電界強度

4.1.3空間電子温度分布

図 4.10 に 50 nsでの空間電子温度分布を示す．等高線は外側から 1 eV，3 eV，5 eVとなっている．この結果よりアンテナとECRとの距離が一番近い N mag =8 が高エネル
ギー領域が大きい．これはやはりECRとアンテナからの距離が近いためである．また加熱された電子は磁気チューブに沿って移動していることがわかる．4.1.2 よりアンテ
ナ正面ではエネルギーの変化が少なかったが，図 4.10 からもアンテナ正面はほとんど加熱されないことがわかる．

(a)

(b)

(c)

(d)

図 4.10空間電子温度分布(a) N mag =8, (b) N mag =9, (c) N mag =10, (d)4.1.4

EEDF(Electron Energy Distribution Function)

図 4.11 に示す赤枠の中の電子のエネルギー分布を調べた．図 4.12 に EEDF の結果を示す．今回の結果ではわずかであるが，磁石数が増えるに従い，高エネルギー電子の数が減っていることがわかる．今回は電子の初期条件 0 eV から始めたので，初期の速度が低く,電子が往復運動する回数が少なかったので,もう少し計算時間を延ばせ
ば各磁石の結果の差は開くと予想される． 30 nsec まで本計算は行っているが，これは実時間では 1 日以上かかっていることと計算誤差が大きくなり正しくない値を出すことからその時間までしか計算を行わなかった．

図 4.11 EEDF 計算領域

図 4.12 EEDF 計算結果

4.2

2D simulation

4.2.1 電子密度分布
図 4.13 に磁石数を変えたときの電子密度分布を示す．3 つに共通していることは磁気チューブ内では周りに比べて電子密度が高いことが分かる．これは磁気チューブ内
で電子が閉じ込められているためと考えられる．また磁場強度の増加に伴い，磁気チ
ューブ内の電子密度がわずかではあるが，高くなっていることがわかる．これは磁場閉じ込めの効果が増大したためである．

磁気チューブ

ヨーク
アンテナ
(a)

(b)

(c)

図 4.13 電子密度分布(a)N mag =9, (b)N mag =10, (c) N mag =11

4.2.2 イオン密度分布
図 4.14 に磁石数を変えたときのイオン密度分布を示す．3 つ共通して電子密度と同じく磁気チューブ内で密度が高くなっていることが分かる．これはイオンは磁場にほとんど捕らわれていないが，この領域では電子が多いために中性を保とうとするため，イオンが多いと考えられる．

磁気チューブヨーク アンテナ
(a)

(b)

(c)
図 4.14 イオン密度分布(a)N mag =9, (b)N mag =10, (c) N mag =11

4.2.3 電離衝突位置
図 4.15 に電離衝突位置の結果を示す． N mag =9, N mag =10 では磁気チューブ中で盛んに電子衝突が起こっていることがわかる．実験での光学測定において磁気チューブ
で は プ ラ ズ マ 密 度 が 高 い と 報 告 (3-1) さ れ て い る こ と か ら ， 実 験 と 良 い 一 致 を 示 す ．

N mag =11 ではECR層で電離衝突が頻繁に起こっていることがわかる．

図 4.15 電離衝突位置

4.2.4電子・イオン損失

電子はヨーク，アンテナ，グリッドに到達したもの，イオンはヨーク，アンテナに到達したものを損失としてカウントした．図 4.16 に電子の損失を，図 4.17 にイオンの損失を示す．電子・イオン共にわずかではあるが，磁石数増加に伴い，損失は減少していることが分かる．これは ∆E が小さくなりプラズマ生成が減ったことと，磁場閉じ込めの効果が増大したことためである．

Number of particles

94000
92000
90000
88000
86000
8

9

10

11

Number of magnets
図 4.16

電子損失

36000

Number of particles

35500
35000
34500
34000
33500
33000
32500
32000

8

9

10

Number of magnets
図 4.17

イオン損失

11

参考文献
(3-1) 近藤慎哉， “30 W 級小型イオンスラスタの内部測定”, 平成 19 年度九州大学修
士論文．

第5章

結論

本研究ではマイクロ波放電式小型イオンエンジン内を，FDTD-PIC コー
ドを用い，電子とマイクロ波の相互作用について調査した．その後，プラズ
マの生成だけでなく，損失についても調べるため 2D3V-DSMC-PIC コードを開発し計算を行った．以下に得られた結論を示す．
(1) FDTD-PIC コードにより，電子がマイクロ波から得るエネルギー増分（ ∆E）と磁場強度の関係を調査した．その結果，磁場強度の増加に伴い，∆Eは減少することが確認できた．これは ECR 層での電場強度の減少によるものである．また円盤アンテナではアンテナ中央よりもアンテナの端の方が電界強度が強く，電子加熱が効率的に行われることがわかった．
(2) イ オ ン を 粒 子 と し て 扱 う こ と が 可 能 で シ ー ス を 計 算 で き る
2D3V-DSMC-PIC コードを開発した．磁石数の増加に伴い，磁場による閉
じ込め効果が増大し，壁面での損失が減少することが確認できた．しかし，今回の計算では実験で測定されたよりも 1000 分の 1 程度の密度にしか到達することができなかった．その対策として，今回検証しなかった質量比・超粒子数・粒子数・衝突計算について考える必要がある．

謝辞
本論文の作成にあたり，丁寧なご指導を頂き，また遅々としてなかなか進まない時でも辛抱強く我慢していただきました中島秀紀教授に深く感謝致します．渡辺幸信准教授には違った視点から貴重なアドバイスを頂き，梶村好宏講師には計算について基本的なことや貴重な助言を頂き，山本直嗣助教には日本にいる時でもアメリカにいる
時でも時間や曜日を問わず様々な疑問に答えてくださったことに深く感謝致します。
いつも明るく研究室を華やかにして頂いた岡安悦子秘書に深く感謝致します。

研究室の皆様には、この 3 年間本当にお世話になりました。特に 3 年間共に研究したコンピュータに強い川渕亮氏，細かいミスも見逃さない近藤慎哉氏，追い込み型の西嶋康太氏，逆境に負けない三好誠氏には本当にお世話になりました．2 年間共に過ごした永田英隆氏，大宅将史氏，松田伸夫氏，本岡親英氏のおかげで楽しく過ごせました．1 年間だけでしたが，鶴哲平氏，七村和彰氏，安部晋一郎氏，小谷優介氏，林田憲治氏にもお世話になりました．半年間だけでしたが，高橋奈菜子氏にはとても刺激を与えていただきました．研究室は違いましたが，Kadrev Dimtre Nikolov 氏，林真照氏，叶涛氏にも刺激を受けました．卒業されました増井博一氏，小川真司氏、加来大輔氏、笠岳幸氏、近岡貴行氏、野田賢治氏、宗貞伸孝氏に深く感謝します。先輩方の研究に取り組む姿勢などとても参考になりました．皆様とこの研究室で出会い、共に研究できたことを幸せに感じています。