目次
第1章

序論

1

1.1 研究背景

1

1.2 イオンエンジンの原理

2

1.2.1

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

2

1.2.2

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

4

1.3 研究目的
第2章

数値解析

2.1 3D simulation

6
7
7

2.1.1

粒子シミュレーション

7

2.1.2

電磁波解析手法

17

2.1.3

coupling コード

23

2.1.4

計算条件

25

2.2 2D3V コード

28

2.2.1 PIC 法

28

2.2.2

静電場計算

28

2.2.3

計算条件

30

2.3 並列化
第3章
3.1

計算結果と考察32
34

3D simulation

34

3.1.1

投入電力

34

3.1.2

ΔE のマイクロ波周波数依存性

35

3.1.3

空間電子温度分布39

3.2

2D simulation

40

3.2.1

電子密度分布40

3.2.2

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

41

3.2.3

電離衝突位置

42

3.2.4

粒子損失

43

3.3 並列化
第4章
謝辞

結論

44
45
46

第1章

序論

1.1 背景1970 年に日本初の人工衛星「おおすみ」が打ち上げに成功してから 40 年，日本の宇宙開発は進歩を続け多くの衛星や探査機が宇宙で活躍してきた．中でも小惑星探査機はやぶさが小惑星イトカワに到着し，サンプル採集を試みたのち昨年地球に帰還したニュースは記憶に新しい．マイクロ波放電式イオンエンジンμ10 をはじめとした様々な新技術の導入により世界初のミッションに成功したのである．
JAXA では現在後継機としてはやぶさ 2 の開発が進められているが，民主党政権による事業仕分けの対象として扱われた．財政難の日本において大型ミッションは資金面での問題は今後もついて回るのではないではないかと思う．また，ロケットの打ち上げ数では宇宙の軍事利用を行っている欧米や中国に溝を開けられているのが現状である．
このような情勢を踏まえて日本が宇宙ビジネスで生き残る道の一つが，衛星その他機器の小型化である．その利点は，開発期間が短くなるためより新しい技術を導入しやすい．単純に小さいことによるコストの削減とそれに伴う制作回数の増加による打ち上げ失敗等のリスクの分散できる．さらに，大型プロジェクトの打ち上げに便乗することで打ち上げコストも大幅に削減可能となることである．実際日本でも 2009 年の温室効果ガス観測技術衛星いぶきの打ち上げや 2010 年の金星探査機
あかつきの打ち上げの際に民間企業や大学等の開発した人工衛星が相乗りする形
で宇宙に上がっている．
人工衛星や宇宙ステーションなどの地球周回物は，地球が楕円球であることや大気抵抗の存在，また太陽や月の引力による摂動があるため，その姿勢制御・軌道制御のために推進機を搭載している．現在，衛星の寿命を決定するものは機器の耐年数よりもむしろ推進剤の枯渇である．特に小型衛星においてはサイズの制約上，搭載可能な推進剤の量も制限されてしまう．従って，尐量の推進剤で長期間運用可能な燃費の良い推進機の開発が重要となる．そして，その要求に応えるものとして利用されているのが電気推進である．
電気推進とは，電気エネルギーを推進剤に与えることで反力を得る推進機であり，従来の化学反応を利用したロケット推進と違って静電力や電磁力を利用できる．このため比推力（推進剤流量あたりの推力）を一桁以上大きくすることが可能であり，大幅な推進剤の低減が望める．現在使用されている代表的な電気推進機はアークジェットスラスタやホールスラスタ，そしてイオンエンジンである．イオンエンジン
は他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制御や惑星探査に適している．

1

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

図 1.1

イオンエンジンの概念図

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

2

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

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

(1.1)

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

3

rL 

mv
eB

(1.2)

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

 ce 

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

E

B

電子

E

B

電子

電子
E

B
B
電子

E

図 1.3 ECR の原理
1.2.2 マイクロ波放電式小型イオンエンジン
本研究室では HAYABUSA(1-5)に搭載されたマイクロ波放電式中和器を元にマイクロ波放電式小型イオンエンジンを作成している．図 1.4 にスラスタの写真を，図1.5 にその構成図を示す．Al 製円筒放電室の周囲に 4×4×12 mm の長手方向磁化
磁石を複数設置し，それらを軟鉄製ヨークで挟み込んで放電室内に磁気回路を形成させている．磁石は永久磁石の Sm-Co (サマリウムコバルト) を用い，アンテナ・
グリッドはモリブデン (Mo) 製である．
4

図 1.4Al

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

Mo

Soft iron

Sm-Co magnet

Propellant gas

Antenna

Front yoke

Microwave

Housing
Back yoke
Discharge chamber

5

図1.5

Screen grid

Acceleration grid

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

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

本研究の目的は，マイクロ波放電式小型イオンエンジン内部のプラズマ状態を模擬するための数値解析コードの改良及びそのコードを用いて投入マイクロ波周波数による依存性の調査を行うことである．

参考文献(1-1) M.Kato, S. Takayama, U. Nakamura, K. Yoshihara and H. Hashimoto, 56th
International 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 年度九州大学博士論文．

6

第2章

数値解析

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

(2.1)

dx
v
dt

(2.2)

ここで， m は粒子質量， v は粒子の速度， x は粒子の位置， t は時間， q は粒子の持つ電荷， E は電場， B は磁束密度を表わす．一般的な実験室プラズマは 1 cm3
当りに 1010～1014 個が含まれ，この粒子全部を取り扱うことは最新のスーパーコン
ピュータを用いても不可能であり，実用的でない．そこで粒子シミュレーションで取り扱う粒子として，実際の粒子の電荷質量比を一定に保ったまま多数の粒子の電荷と質量を１つにまとめた超粒子(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 
t
m
2






(2.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と書ける．ここで，  ・  ・  ・  を次のように定義する．
qt n

Bx
2m

7

(2.5)

(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

 1

A 
 


 

 



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)

計算に用いるメッシュサイズ  は，電子温度 Te，ボルツマン定数 kb，電子密度
ne から計算される．
デバイ長 λDe，De 

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

  De(2.15)であることが望ましい．また，タイムステップはプラズマ周波数 ωpe，
ne e 2me  0 pe 

8(2.16)

に対して，
1
t   pe(2.17)とすれば，数値的な不安定性は生じない．時間積分の精度まで考えると，通常は
1
t  0.2 pe(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, dx2, dy1, dy 2, dz1, dz 2 を求める．その値を用いて，図中にあるような直方体 V1～V8の体積を求める．この体積比から粒子が感じる電磁場を内挿することができる．
具体的には，座標 x , y , z  上にある粒子が感じる電場 E x , y , z  は，計算セルの体積
を V ，図に示すような粒子によって仕切られた空間体積をそれぞれ V1～V8 とすると，格子点上の重みは
An 

Vn
V

n  1 ～ 8となり，
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

(2.21)

 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)と表わされる．

9図 2.12.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)間に一様に分布する乱数である．

10

2.1.1.6

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

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

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



v


 Nv

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

11(2.36)

2.1.1.7

衝突確率，衝突判定

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

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

(2.38)(2.39)

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

衝突計算

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

 elasticのとき弾性衝突 total elastic  excitation r2  elasticのとき励起衝突 total total elastic   excitation  excitation   ionization r2  elasticのとき電離衝突 total total 0  r2 

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

12

図 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 

13(2.41)

hy  hz  

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 とすれば，エネルギー収支は以下のように表わされる．
(2.46)

Escattered,e  Eincident,e

電子は，エネルギー 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)
として与える．
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーを Eexcitation とし，入射電子のエネルギーを Eincident,e ，励起後のエネルギーを Escattered,e ，励起後の速度を v~ とする．エネルギー収支は，Escattered,e  Eincident,e  Eexcitationと表わされ，励起後の速度は

14(2.50)

Ev~  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)式から求められる Escattered,e を用いる．励起エネルギ
ー Eexcitation は，キセノンの場合 8.34 eV を用いている．
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを Escattered,e ，生成電子のエネルギーを Ecreated,e ，電離エネルギーを Eionization とする．エネルギー収支はEscattered,e  Ecreated,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)式から求められる Escattered,e を用いる．電離エネルギー Eionization
は，キセノンの場合 12.13 eV である．
電離後の生成電子の持つエネルギーは，Ecreated,e  E  Eionization(2.56)となる．また，入射電子の速度 v e を用いて生成電子の速度はv created,e 

ve
ve

15

2 Ecreated,em(2.57)

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

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

 constant(2.60)

を満たす  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 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度で規格n化した衝突周波数を示す．本計算では Null-Collision 法で用いた collision gas の値は 2×10-12 m3/sec である．

16

図 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 ）が電場の方程式を満たすように配置さ17

れている．
この電磁場の配置のもと，(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



z
y


(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



x
z


n

Hz

1
2

n

i  1 2 , j  1 2 , k   H z

12

ti  1 2 , j  1 2 , k  
 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 



y
x


(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  12
 H z i  1 2 , j  1 2 , k   H z 2 i  1 2 , j  1 2 , k 

y


n



Hy

1
2

n

1
2

i  1 2 , j , k  1 2  H y i  1 2 , j , k  1 2z
 J x i  1 2 , j , k 


(2.66)

18

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


n



Hz

1
2

n

1
2

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

(2.68)(2.63)~(2.68)式から，時間ステップがn  1 2 での磁場 H は，その半ステップ前で
1n
nの電場 E と１ステップ前の磁場 H 2 から求めることができ，時間ステップが
n 1
1
n

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

19

図 2.4

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

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

図 2.5

同軸線路給電20図 2.6 太さ零のアンテナモデル
2.1.2.4

安定条件

計算結果の正確さを保つため，空間の増分値 x , y , z は波長  に比べて十分小さ
くなければならない．（  10 以下）
時間ステップ t は，Courant の安定化条件より
1
t 
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 の負の方向に進む波はEz  Ez x  vt

21(2.70)

と表わされ，次の微分方程式を満足する．
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

v

E z
x

(2.72)

12となる．右辺の x についての微分は，電場が x  0, x に割り当てられていることから， x  x 2 で行い，(2.71)式に代入すると

E 1 2  E
t
n
z

n 1
z

1 2

n

v

Ez

1
2

n

1
2

1  Ez 0(2.72)

xとなる。
(2.72)式の Ezn 1 2, Ezn1 1 2 も FDTD 法では直接求めることはできないため，前後の値で平均して求める．そして， E z についてまとめるとEzn 0, j , k  1 2  Ezn1 1, j , k  1 2 vt  x n
Ez 1, j , k  1 2  Ezn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/m2 の単位を持つ瞬時電力密度である．この式をア

22

ンテナ上で計算することにより投入電力を測定した．

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 回計算することとした．

23

図 2.7

Coupling code のフローチャート

24

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 12 sec ， FDTD 法 の 時 間 刻 み は
2.1.1.2 や 2.1.2.3 で述べた安定条件を満たす．FDTD
t FDTD  5.0  10 13 sec とすれば、コードが 2 回計算された後，PIC コードを 1 回計算させた．
2.1.4.4 電子初期条件
Ｔe＝0 .5eV となるマクスウェル分布の初期速度を持った電子を放電室内に一様に分布させた．

25

2.1.4.5 粒子数
1 ボクセル当たり数個しか粒子がない場合，ノイズが大きくなり望ましくない．通常，10 から 100 の間の粒子が存在しているとよいとされている．今回は 1
ボクセルに平均 20 個の粒子が含まれるように粒子数を設定した．
(2-14)

2.1.4.6 磁場強度
今回はマイクロ波周波数依存性を調べるために電源として 900MHz、2.45 GHz、
4.2GHz を用いる。マイクロ波周波数と電子サイクロトロン周波数が一致するときの磁場強度はそれぞれ 32mT、87.5mT、150mT である。その値を示す領域を ECR層と呼ぶ。この ECR 層の位置が各周波数で同じになるように放電室外側の磁石の数を 900MHz では 4 個、2.45GHｚでは 12 個、4.2GHz は 19 個に設定し計算を行った．磁場計算には Advanced Science Laboratory, Inc の 3 次元静磁場解析ソフトを使用した(2-23)．図 2.7 に例として 2.45GHz 用の解析結果をプログラムで読み込んだ後の放電室内磁場強度を示す．

図 2.7

放電室内磁場強度（2.45GHz 用）

26

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

qj
dv
vv

( E  vz Br  vr Bz )  r dt
mjr(2.75)qjdvz
( Ez  r Br )dtmj
j は粒子種をあらわし，j = e, i であり，qe = -e，qi = e である．また，衝突計算には同
じく Null-Collision 法を採用している．電子―中性粒子の弾性・励起・電離衝突のみ考慮している．

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

 2 1    er    ni  ne 2zr r  r 
0(2.76)

φは空間電位，ε0 は真空中の誘電率，ni と ne はイオン・電子の数密度である．場
27

の 方 程 式 を 解 く に あ た り ， GMRes(Generalized Minimal Residual method) 法(2-24),(2-25),(2-26),(2-27)を採用した．
この方法は非対称行列を有する線型方程式の反復解法であり，次のような考え方にもとづいている．解くべき方程式はAx = b
(2.77)
であるとする．ここに A は(N×N 行列)であり，簡単のために実行列であるとする．
b は既知の N ベクトル，x は未知の N ベクトルである．いま，適当に初期過程 x0を選ぶと，(2.77)の残差 r0 は
r0 = b - Ax0(2.78)と書ける．いま，第 n 回の繰り返しにおける解 xn を xn = x0 + zn と書くことにする．
ここに zn は，Krylov 空間
Kn(A,r0) = Span(r0,Ar0,…,An-1r0)(2.79)||b - Axn|| = ||r0 - Azn||(2.80)

からを最小にするように選ぶものとする．記号 Span()は，括弧内のベクトルの線形結合で作られる空間を表す．
(2.80)の最小化を行うために，Kn(A,r0)の正規直交基底(v1,v2,…,vn)を次のようにして逐次構成する．まず||r0|| = γ とすると v1 = r0/γ となる．いま，第 n 番目の基底 vnまでが求められているとする．このとき，
wn=Avn, hin=vi・wn(2.81)とおき，Schmidt の直交化により
Vn+1 = wn - Σhinvi(2.82)とおけば，Vn+1  Kn(A,r0)である．さらに Vn+1 = Vn+1/||Vn+1||とすると，
hn+1 n := vn+1・wn = vn+1・Vn+1 = ||Vn+1||(2.83)

を得る．
次に，(2.80)における最小化を実行する．第 n 回目の繰り返しにおいてはzn  Kn(A,r0)であるから，zn = Σvjyj(2.84)と書ける．これを用いて(2.80)を v1,…,vn+1 座標系で書くと，min || γδi1 – Σhijyj|| (i = 1,2,…,n+1)が得られる．上式の行列 hij は(n+1)×n の上 Hessenberg 行列である．

 h11 h12 h21 h22
 0 h
32
n
Q
  
 0 

 0 


h13
h23
h33

0
0





hn 1 n0h1n   h11 h12
 
h2 n   0 h22
h3n   00    
hnn   0 
hn 1 n   0 

28

h13
h23
h33

0
0





0
0

h1n 

h2 n 
h3n 
 

hnn 
0 

(2.85)(2.86)

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

 h11 h12
 0 h22
 0
0

  

 0 
 0 


h13
h23
h33

0
0





0
0

h1n 
 y1 
 Q11n 
 n
h2 n  y Q21 
2h3n  n 
y 3    Q31




 







hnn   n
 yn 
 Qn1 0 

(2.87)

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

2.2.3 計算条件
2.2.3.1 計算領域
フロントヨークからバックヨークにかけての磁気チューブ内でのプラズマ挙動，
およびその両端のシースポテンシャルが放電室内部の物理現象に深く関与するものと考え，図 2.8 での赤枠内を計算領域とした．格子数は全体で 92×78 にとった．
格子が粗いとプラズマシースにおける非中性の現象を正しく追うことができない
ため，格子間隔を Debye 長と同程度にする必要がある．等間隔の格子点が望ましい
が，計算コストが大きくなりすぎるため，式(2.89)を用いて不等間隔に格子を作成した．
xi exp(aI )  exp(a( I  i))（2.89）
(i  0,1, , I )L
exp(aI )  1
L は体系長さ，I は分割数、x0, x1, …, xI はセルの境界の座標，a は格子の粗さを調節する係数である．今回は特に重要な磁気チューブ中のアンテナ前面・上部・フ
ロントヨークの付近とグリッド付近では Debye 長と同程度にした．式(2.89)を用いて作成した計算格子を図 2.9 に示す．

29

図 2.8

計算体系概念図

図 2.9

30

計算格子

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

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

mv
eB

(2.91)eBm(2.92)

c 

イオンについては重み係数の影響は無い．電子の Larmor 半径を補正するため，電子の感じる磁場の大きさを  倍大きくする．このとき Larmor 周波数は  倍小
31

さくなるが，Larmor 半径に比べその影響は小さいと考えられるのでここでは無視する．
②衝突断面積
衝突断面積は電子速度に依存する．電子の速度は1倍となるので，その小さくなる影響を，衝突断面積を  倍大きくすることで補正する．
2.3 計算コードの並列化
計算にかかる時間を短縮することはより多くの条件に対して検証するためには非常に重要である。1 台のコンピュータ(1 つの CPU)で行っていた計算をネットワ
ークでつながれた複数のコンピュータ(複数の CPU)で計算できるように改良することにした。用いたのは MPI(Message Passing Interface) 、分散メモリ型の並列計算機で 2 つのプロセス 間でのデータをやりとりするために用いるメッセージ通
信操作の仕様標準である。

参考文献(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.
32

(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.

33

第 3 章 結果と考察
3.1 3Dsimulation
3.1.1 投入電力
投入電力の確認を行う。先に述べたとおり投入電力はアンテナヨークの放電質側開口部から流入するエネルギーをポインティングベクトルを計算することによって求める。図3.1 は縦軸にポインティングベクトルを時間積分した値，横軸に時間を取ったものである．
この傾きが投入電力にあたる。今回の計算に用いた 900MHz、2.45GHz、4.2GHz の投入
電力が同程度であることが確認した。

図 3.1 投入エネルギーの時間変化

34

3.1.2 ΔE のマイクロ波周波数依存性
ΔE の依存性を図 3.2 に示すような 3 つの領域に分けて調べた。これはアンテナとECR の距離だけでなく位置関係によって ΔE に差異が生じるためである。

図 3.2 解析領域
結果を図 3.3，3.4 に示す。ヒストグラムの幅は 0.1eV になっている。領域 1、2 に比べ領域 3 でエネルギー変化が小さいのはアンテナ側面に比べ正面では電界強度が低いことが原因である。周波数による違いを見ると、周波数増加に伴いわずかではあるが ΔE の減少が確認できる。前述の領域による違い同様、ECR 層近傍における電界強度によるものではないかという予測から各領域の ECR 層の 1 点を抽出しそれぞれ比較を行った結果が図 3.5 である。周波数が低いほうが電界強度が高いことが確認できた。

35

図 3.3（左）領域 1，（右）領域 236

図 3.4 領域 337

図 3.5 上）電界強度分布，下）ECR 層における電界強度比較

38

3.1.3 空間電子温度分布
図 3.6 に空間電子温度分布を示す。900MHz で確認できるアンテナ近傍の高温領域
が周波数増加とともに減少し、4.2GHz においてはほぼ消失している。また 4.2GHz では
壁面近傍で電子温度が高くなっている。これは先に示した電場強度を反映しているが、壁面近くでは計算における粒子数が少ないことが考えられるので正しい値を示しているか確認が必要であると思われる。

図 3.6 空間電子温度分布 左）900ＭＨｚ，中）2.45GHz，右）4.2GHz

39

3.2 2Dsimularation
3.2.1 電離衝突位置
電離衝突位置の結果を図 3.7 に示す。各周波数に共通して見られるのは，磁気チューブの中でもフロントヨークとアンテナヨークを結ぶ領域において衝突頻度が多いということである。これは磁力線と ECR 層の関係によりエネルギーの高い電子が多いという
FDTD の結果を反映したものと考えられる．周波数の増加とともに衝突位置がアンテナ前方とフロントヨークを結ぶ領域にも広がりを見せている．これは磁場強度の増加により電子の閉じ込め効果が高まり、ECR を通過する回数が伸びることで電子がより高いエネル
ギーに達したことを示すものである．

図 3.7 電離衝突位置上左）900MHz，上右）2.45GHz，下）4.2GHｚ

40

3.2.2 電子密度分布
図 3.8 に各マイクロ波周波数における電子密度分布を示す。共通して磁気チューブ内での密度が周囲に比べて高いことが確認できる。また、磁場強度が強く閉じ込め効果の高い高周波数になるほど電子密度が高くなっている。FDTD の方での電子温度分布とは異なりアンテナ前方の磁気チューブで高い値を示しているのもアンテナ後方とを結ぶ部
分よりも磁気チューブが短く磁場が強いためと考えられる。

図 3.8 電子密度分布上左）900MHz，上右）2.45GHz，下）4.2GHｚ

41

3.2.3 イオン密度分布
図 3.9 はイオン密度分布を示したものである．2.45GHｚと 4.2GHｚでは顕著に 900MHｚ
では微弱ではあるが磁気チューブで電子密度が高かった領域のイオン密度が高くなっていることが分かる．これはイオンには閉じ込め効果はあまり働いていないが中性を保とうとするものだと考えられる．

図 3.9 イオン密度分布上左）900MHz，上右）2.45GHz，下）4.2GHｚ

42

3.2.4 粒子損失
電子はヨーク，アンテナ，グリッドに到達したもの，イオンはヨーク，アンテナに到達したものを損失としてカウントした．電子イオンともに周波数が高くなるにつれて損失は減少した．これは磁場が強くなったことによる閉じ込め効果の増大を示す妥当な結果といえる．
300000250000200000
グリッド損失150000アンテナ損失ヨーク損失

10000050000
0900MHz

2.45GHz

4.2GHz

14000012000010000080000
アンテナ損失
60000

ヨーク損失

40000
20000
0
900MHz

2.45GHz

4.2GHz

図 3.10 粒子損失の周波数依存性，上）電子損失，下）イオン損失433.3 並列化
まずはコードの中でより多くの計算時間を要している場所がどこであるかを調べた．結果を示したのが図 3.11 である．粒子の運動を解く部分である maincalc とその結果から各種密度を算出する density の 2 つのコードが際立って計算に時間がかかっていることがわかる．粒子の運動はそれぞれ独立しているため、各 CPU に分割し並列して計算し、各種密度については各 CPU の小計を通信を行って合算することにした．この改良による計算時間と CPU 数の関係をまとめたのが図 3.12 である．最大 16 並列において計算時間
は約 6 分の 1 になった。

図 3.11 計算コード消費時間内訳
120

計算時間，％

100
80
60
40
20
0
1

2

3

4

8

12

並列数

図 3.12 並列数による計算時間の変化

44

16

第4章

結論

本研究ではマイクロ波放電式小型イオンエンジン内を，電子のみを粒子として扱う FDTD-PIC コードを用い電子とマイクロ波の相互作用について調査した．その後、イオンも含めたプラズマ全体の生成と損失を調べるため 2D3V-DSMC-PIC コードにて計算を行った．以上の結果から得られた結論を以下に示す．
1）FDTD-PIC コードにより，電子が ECR 加熱によってマイクロ波から得るエネ
ルギーの増分（ΔE）とマイクロ波周波数の関係を調査した．その結果，マイクロ波周波数の増加に伴いΔE の減尐を確認した．これは ECR 層における電界強度が高周波数になるにつれて減尐していたことが一因である．
2）2D3V-DSMC-PIC コードによって，マイクロ波周波数が高くなるのに伴う磁石数の増加から，磁場閉じ込め効果の増大を確認した．それに伴いΔE の低いアンテ
ナ前方の領域へ電離衝突位置が拡大している様子や，壁面への粒子損失の減尐が認められた．
3）メモリ分散が並列計算用通信メッセージライブラリである MPI を用いてコードの計算時間短縮に取り組んだ．もっとも比重の重かった粒子計算部分を並列化することにより最大 16 並列で計算時間を６分の１に短縮できた．

謝辞
3 年間ご指導いただきました，中島秀紀教授，山本直嗣助教には感謝の言葉もありません．無事修了できた暁にはどんな形であれ恩返しをしていければと思います．
研究では直接関係はありませんでしたが渡辺幸信准教授に頂いた研究にとどまらない様々なお話は人生の糧になるものと思います．
また，私が携わった研究を代々引き継いできた諸先輩，ともに学んだ同級の友にもこの場を借りてお礼を申し上げます．