第1章

序論

1.1 背景1945 年、アーサー・C・クラークは” Extra-Terrestrial Relays — Can Rocket Stations
Give Worldwide Radio Coverage? “(1-1) の中で，静止衛星を用いた通信方法を提唱した．70 年近く経過した現在，彼の思い描いた”地球外の中継”は現実のものとなり，
彼のアイディアの中核を成した人工衛星は今や，我々の生活に欠かせないものとなっている．ここ数十年の間に着実に進歩を遂げてきた人工衛星の中で，近年注目を集めているのは小型人工衛星である．
従来の大型衛星は重要なミッションを担うため，信頼性の高さが望まれていた．
そのため，使用される部品の多くは過酷な宇宙環境に耐えうる専用品でなければな
らず，開発費は膨大なものとなった．衛星所有にかかる膨大な初期投資は企業の参入を妨げ，宇宙の商業利用の障害となっていた．
これに対して小型衛星は，開発費が安く済み，開発期間を短縮できるという利点がある．また、低質量であることを利用して一回の打ち上げで複数の衛星を宇宙に上げることも可能となる．低リスクの打ち上げが期待できる小型衛星は，企業参入の敷居を下げ，新サービスを提供する機会を生むのではないかと注目を受けている．
衛星は大気の存在や地球，月の引力の影響を受けるので，軌道・姿勢制御のための推進器を要する．現在，衛星の寿命を決定するのは推進剤の枯渇，すなわち推進器の燃料の枯渇であり，衛星の推進器は高効率であることを求められる．これは小型衛星も同様なのだが，小型衛星にはそのサイズゆえに小型且つ高効率の推進器が必要となってくる．この要求を満たすものとして使用されているのが電気推進器である。
電気推進(1-2)とは，電気エネルギーを推進剤に与えることで反力を得る推進機であり，従来の化学反応を利用したロケット推進と違って静電力や電磁力を利用できる．
このため比推力（推進剤流量あたりの推力）を一桁以上大きくすることが可能であり，大幅な推進剤の低減が望める．現在使用されている代表的な電気推進機はアー
クジェットスラスタやホールスラスタ，そしてイオンエンジンである．イオンエン
ジンは他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制御や惑星探査に適している．

1

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

図 1.1 マイクロ波放電式イオンエンジン概略図

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

2

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

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

(1.2.1)

ここで m は電子の質量，v は速度ベクトル，e は電荷量，B は磁束密度である．このときの円軌道の半径 rL はラーマ半径と呼ばれmv
rL (1.2.2)eB
で与えられる．ここで v⊥は電子の B に対する垂直な速度成分の大きさである．
プラズマを構成している荷電粒子は全て反磁性体である．そのためサイクロトロン運動の回転の方向は，外部磁場の向きに対して荷電粒子の回転によってできる磁場が常に逆を向く方向である．つまり電子は磁場に対して右回りの回転運動を行う．
この回転運動の角周波数はサイクロトロン角周波数と呼ばれ

 ce 

3eB
m(1.2.3)

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

E

B

電子

E

B

電子

電子
E

B
B
電子

E

図 1.3 ECR の原理

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

4

1. ヒーターを用いないため，予備加熱が不要となり素早い起動が出来る．
2. マイクロ波放電式イオンエンジン駆動用のマイクロ波電源が 1 台あれば中和器のプラズマも同時に生成することが可能なため，システム構成の簡素化や信頼性の向上，また衛星の重量軽減が図れる．
マイクロ波放電式中和器は深宇宙探査機「はやぶさ(1-3)」に搭載され十分な中和性能と寿命が得られたという実績がある．しかし，小型衛星に搭載するようなサイズの中和器は未だ実用化されておらず，その開発が急務となっている．
当研究グループは小型衛星への搭載を想定した小型マイクロ波放電型イオンエンジンを開発してきた．エンジン本体について優れた性能が得られた．一方中和器においては要求されるサイズ，消費電力でも中和性能の目標値をクリアしておらず，
現在改善に向けて研究中である．
図 1.4 はマイクロ波放電型小型中和器の外観写真である．放電室は内径 21 mmの円筒形状であり、オリフィス開口径 10 mm、放電室長 32 mm である．放電室外部と出口部分に Sm-Co 磁石を配置し構成された磁気回路によって ECR に必要な磁場を供給する。アンテナはモリブデン製である．

図 1.4

マイクロ波放電式小型中和器の写真

5

図 1.5

図 1.6マイクロ波放電式小型中和器の内部構造

マイクロ波放電式小型中和器運転時

61.3 研究目的
本研究室では，これまで小型イオンエンジンの開発，特に中和器の開発に注力してきた．中和器の形状や内部の磁場を変更して実験 (1-4)(1-5)(1-6)を行ってきたが，中和器が小型であるために内部のプラズマを直接測定できないという問題が生じている．また，より良い性能を得るために多くのパラメータを試行錯誤して変更していくのは，時間，費用共に限界がある．そこで，内部物理の解明及び中和器の性能を評価する数値計算コードを開発してきた．これまでの計算(1-7)(1-8)(1-9)の問題点は，計算コストが障壁となり，定常状態まで計算できていなかった点である．本研究では，マルチグリッド法を用いて計算コストを低減し，プラズマ電位分布を厳密に解析した場合のプラズマ密度分布及び密度分布のパラメータ依存性を算出し，それに基づく設計の最適化を目的とする．
参考文献(1-1) Arthur. C. Clarke ， ”EXTRA-TERRESTRIAL RELAYSCan Rocket
Stations Give World-wide Radio Coverage?”, 1945, Wireless World．
(1-2) 栗木恭一・荒川義博[編]：”電気推進ロケット入門”, (東京大学出版会,2003)(1-3) I. Funaki, H. Kuninaka and K. Toki, J. Propul. Power, Vol. 20, (2004), pp.718-726.
(1-4) N. Yamamoto, H. Kataharada, H. Masui, H. Ijiri and H. Nakashima,
AJCPP2005-22093, 2005.
(1-5) N. Yamamoto, H. Kataharada, T. Chikaoka, H.Masui and H. Nakashima,
IEPC-2005-036, 2005.
(1-6) T. Chikaoka, S. Kondo, N. Yamamoto, H. Nakashima and Y. Takao, Proceedingsof the 25th International Symposium on Space Technology and Science, 2006,
pp.254-259.
(1-7)
増井博一， “マイクロ波放電型宇宙推進用プラズマ源に関する数値解析”,
平成 17 年度九州大学博士論文．
(1-8)
金川隆保, “マイクロ波放電式小型イオンエンジンの数値解析”，平成 19
年度九州大学修士論文．
(1-9) 新谷將 “マイクロ波放電式小型イオンエンジンの数値解析”, 平成 22 年度
九州大学修士論文．

7

第2章

数値解析

2.1 粒子計算手法(2-1)
2.1.1 基礎式
プラズマの挙動を解析する際の基礎式として運動方程式を用いる．式を以下に示す．

mdv
 q E  v  B dt(2.3.1)dx
v
dt(2.3.2)

ここで， m は粒子質量， v は粒子の速度， x は粒子の位置， t は時間， q は粒子の持つ電荷，E は電場， B は磁束密度を表わす．(2.1)式は二次元円筒座標系においては次のように書き下
せる。

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

(2.3.3)

これらを四次精度 Runge-Kutta 法にて解いた。

vrn1  vrn tbr1  2br 2  2br 3  br 4 6vn1  vn tb1  2b 2  2b 3  b 4 
6vzn1  vzn tbz1  2bz 2  2bz 3  bz 4 
6(2.3.4)

但しここで定数 b はそれぞれ以下を用いる。
qEr  v Bz   1 v v
mr
tb1   1 tb1 tb1 q
 Bz    v  v br 2   Er   v m
2   r
2 
2 


br1 

br 3 

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

br 4 qEr  v  tb 3 Bz  1 v  tb 3 v  tb 3 mr

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

b1 

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

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

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

一般的な実験室プラズマは 1 cm3 当りに 1010～1014 個が含まれ，この粒子全部を取り扱うことは最新のスーパーコンピュータを用いても不可能であり，実用的でない．そこで粒子シミュレーションで取り扱う粒子として，実際の粒子の電荷質量比を一定に保ったまま多数の粒子の電荷と質量を１つにまとめた超粒子(2-2)を使用する．超粒子数は，１セルに含まれる粒子数が 10個以上になる事と取り扱う粒子数よりも超粒子の重みが尐ない事（統計性から）から決定する．

2.1.2 メッシュサイズ・タイムステップ(2-3)
計算に用いるメッシュサイズ  は，電子温度 Te，ボルツマン定数 kb，電子密度 ne から計算される．
デバイ長 λDe，De 

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

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

ne e 2me  0 pe 
に対して，

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

c 

eB
me

(2.3.13)に対して，t  c1(2.3.14)である必要がある．

2.1.3

PIC 法(2-4)

電場解析では，電場量を空間格子点上で求める.そのため,各粒子の電荷を格子点上に振り分ける必要がある.また,粒子の運動方程式を解く際に,電場解析によって格子点状に割り振られた値を内挿する必要がある.このような計算を行うため，粒子シミュレーション
において標準的な手法である PIC 法（Particle In Cell method）を用いる．
まず，粒子の位置から近傍のセル座標 i, j  を求める．次に，図 2.1 に示す dr1, dr 2, dz1, dz 2を求める．その値を用いて，図中にあるような長方形 S1～S 4 の面積を求める．この面積比から粒子が感じる電磁場を内挿することができる．
具体的には，座標 r, z  上にある粒子が感じる電場 E r , z  は，計算セルの体積を S ，図に示すような粒子によって仕切られた空間体積をそれぞれ S1～S8 とすると，格子点上の重
みは

An 

Sn
S

n  1 ～ 4

(2.3.15)となり，

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

10(2.3.16)図 2.1.1 PIC 法におけるセルの構成

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

 i, j  

S1 qN
S V

(2.3.17)となる．すべての粒子について和をとれば，格子点状での電荷密度が得られる．
2.1.5 荷電粒子と中性粒子の衝突
荷電粒子は電磁場の影響を受けて，中性粒子と衝突する．プラズマは電子，イオン，中性粒子からなるとする．電子とイオンの数密度は中性粒子の数密度より非常に小さいため，電子-電子，電子-イオン，イオン-イオンの各衝突は無視できるものとする．電子と中性粒子は弾性，励起，電離衝突をすると仮定する．なお，以下に用いる乱数は(0,1)間に一様に分布する乱数である．生成周期の短い組み込み乱数 ran()の使用は避け Mersenne
Twister(2-5)を使用した．周期は(219937-1)であり，生成速度も早い．

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

11

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

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



1
N

(2.3.22)

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



v


 Nv

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

12(2.3.26)

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

   E 

(2.3.27)

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

 total  ngas totalvm(2.3.28)

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

(2.3.29)

 1  exp   total t 

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

2.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 に示す．今回の計算では，多価電離は考慮していない．

13

図 2.1.2 電子と Xe の衝突断面積

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

v1 

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

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

v e 

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

h の直交座標成分は

14(2.3.32)hx  ver cos hy  

hz  

(2.3.33)

vex vey cos   ve vez sin (2.3.34)vervex vez cos   ve vey sin (2.3.35)

ver

ここで，

ver  vey2  vez2

(2.3.36)ve  vex2  vey2  vez2

(2.3.37)であり，  は散乱角，  は任意の角である．(2.3.32)式はエネルギー損失を含んでいる．

次に，衝突前後のエネルギー収支を考える．衝突前の電子のエネルギーを Eincident,e とし，衝突後の電子のエネルギーを Escattered,e とすれば，エネルギー収支は以下のように表わされる．

Escattered,e  Eincident,e(2.3.38)

電子は，エネルギー Eincident,e が小さいとき等方散乱され，大きいときは主に前方散乱となる．これを考慮した散乱角  の確率密度を g   として以下に示すエラー! 参照元が見つかりません。．ここで， Eincident,e の単位を eV とする．

g   

Eincident,e sin 0     21  Eincident,e sin 2 ln 1  Eincident,e 2(2.3.39)

これを解けば，

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

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

として与える．

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

15

Escattered,e  Eincident,e  Eexcitation(2.3.42)Ev~  v e 1  excitation
Eincident,e(2.3.43)と表わされ，励起後の速度はとなる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失った後に(2.3.43)式中 v~ の速度で弾性衝突をすると考えることができる．(2.3.32)式中 v e に v~ を用
いて，衝突後の速度 v e を求める．このときの散乱角には(2.3.40)式を用いるが(2.3.40)式中の Eincident,e には(2.3.42)式から求められる Escattered,e を用いる．励起エネルギー Eexcitation は，キセノンの場合 8.34 eV を用いている．
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを Escattered,e ，生成電子のエネル
ギーを Ecreated,e ，電離エネルギーを Eionization とする．エネルギー収支はEscattered,e  Ecreated,e  Eincident,e  Eionization(2.3.44)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる(2-12)．

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

(2.3.45)

ここで，B はガスに固有の値でキセノンでは 8.7 eV となる．
また，入射電子が電離によって失うエネルギーを E とすると，

E  Eincident,e  Escattered,e(2.3.46)と表わされる．(2.3.43)式から電離後の電子の速度は，v~  v e 1 EEincident,e(2.3.47)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．(2.3.34)式中 v e に v~ を用いて，衝突後の速度 v e を求める．散乱角には(2.3.40)式を用いるが，(2.3.40)式中の Eincident,e には(2.3.45)式から求められる Escattered,e を用いる．電離エネルギー Eionization は，キセノンの場合
12.13 eV である．

電離後の生成電子の持つエネルギーは，Ecreated,e  E  Eionization

16(2.3.48)

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

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

2.1.10 Null-collision 法(2-11)(2-12)(2-13)
計算時間短縮のため，Null-collision 法を用いる．これには，架空の衝突断面積  fake を導入する．すべての電子のエネルギー E に対して

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

 constant(2.3.52)を満たす  fake を導入すれば， collision が一定，すなわち Pcollision,m が一定となりエネルギーに依存しない．したがって，粒子毎に Pcollision,m を計算する必要が無くなり，全粒子に対して衝突計算を行わず，衝突を起こす粒子のみについて衝突の種類を決めればよい．これにより，計算時間を大幅に短縮できる．衝突の種類を決めるときは， 2.1.8 において，

 t o t a l  e l a s t i c  e x c i t a t i on  i o n i z a t i on  fa k eとし， elastic   excitation   ionization
 r2  1
 totalのとき Null-Collision Processとする．
Null-Collision Process が選ばれたときには，粒子の運動は何ら変化しないものとする．図2.1.3 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度で規格化した衝突周波数を示す．本計算では Null-Collision 法で用いた

17の値は 2×10-12 m3/sec である．

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

2.1.11 電子質量の変更
今回の計算において，電子の速度成分は最大で 105 m/s 程度，イオンにおいては最大 103 m/s
程度であり，電子の速度を基準に時間刻みを設定するとイオンはほとんど動かない．そこで，電子の質量を実際のものよりα倍大きくした．質量の変更により，電子の速度は√α分の 1 に低下，イオンとの速度差は小さくなり計算結果の収束を早めることができる．

me _ virtual  me _ real

(2.3.53)

また，電子の質量変更によって.ECR 層の位置が変化してしまうことを防ぐため，電子の感じる磁場に補正を与えた．Rarmor 半径は以下の式で与えられるので，電子の感じる磁場を√α倍することによって，Rarmor 半径を現実のものと等しくした．

rc 

mv
eB

(2.3.54)2.1.12 粒子の境界条件
放電室内のプラズマが実際に壁面に衝突した場合，放電室壁面は導体であるので，多くの場
合，電子は電流として導体内に流れ，イオンは導体から電子を奪い中性粒子に戻る．そこで，今回の計算では電子，イオン共に壁面に衝突した粒子は消滅し，衝突後の挙動は追わないものとした．また，中心軸を超えた粒子は，軸対称性から反射されて運動を続けるとした．

18

2.2 マイクロ波解析手法(2-14)(2-15)(2-16)マイクロ波の伝搬解析は，別途 FDTD-PIC コード用いて行った．FDTD-PIC コードとは，
FDTD 法（Finite Difference Time Domain method）と PIC 法を組み合わせたコードであり，電子の挙動とマイクロ波の伝搬を同時に解析できる。このコードを使用して，ECR 層を通過する電子のエネルギー増分を調査，ECR 層を通過する際に電子が得るエネルギーの確率分布を求め，今回のコードに導入した．すなわち，今回の作成したコード内で電子は ECR 層を通過する際，上述の確率分布に従って円周方向へ加速される．

19

2.3 電場解析手法(2-17)
2.3.1

基礎式

ポアソン方程式を基礎式として計算を行った．
．

 2  0(2.3.1)(2.2.1)式を差分化し,連立方程式として解いた．差分式を以下に示す．

aa1  i 1, j  aa 2  i 1, j  aa3  i , j 1  aa 4  i , j 1  aa5  i , j 

1
z 2
1
1
aa3  2 r
2ri , j ri , j0

aa1  aa 2 

aa 4 

(2.3.2)1
1
2
r
2ri , j r

1 
 1
aa5   2  2 
r 
 z
ここで， i，j はそれぞれ z，r 方向の格子点番号に対応している．φi,j は i，j における電位を表し，ρi，j はその点での電荷密度を表す．⊿z，⊿r は各方向の格子サイズであり，ri,j はその格子点の r 方向における位置を表す．

2.3.1 反復法
(2.2.2)式はすべての格子点に１つずつ与えられるので，連立方程式として解くことができる．
この連立方程式を解く方法として，反復法を用いた．
I 番目の方程式が以下のように表せる連立方程式を考える．
na x  b
j 1

i, j

j

i

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

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

20(2.3.4) bi   ai , j x j ( k 1) 


j i
(k )

xi  
ai , j(2.3.5)

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

xi

(k )

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

ωを反復解の収束率を加速させるような値に選択することにより，ガウス-ザイデル法より早い収束を得ることができる．ωには最適値が存在するが，最適値予め知ることは一般的には不可能である．そのため，しばしばω=2-O(h)という発見的な評価式が使用される．ここで h は元の物理領域を離散化するときのメッシュ間隔である．今回の計算では，この SOR 法に対し後述のマルチグリッド法を適用した．

2.3.2 マルチグリッド法
反復法の計算コストを低下させるテクニックとしてマルチグリッド法がある．Jacobi 法や
ガウス-ザイデル法のような単純な反復法は，誤差の高周波成分を最初に収束させることが知られている．この特徴を利用し，複数枚の粗さの異なるグリッド(空間格子)を用意して，それぞれのグリッドに対する高周波成分を収束させるのがマルチグリッド法である．
マルチグリッド法の計算の進め方は数多くあるが，今回はオーソドックスな V サイクル法を用いた(図 2.3.1)．図 2.3.2 に流れ図を示す．まず，細かいグリッドと粗いグリッドを用意し，細かいグリッドでほんの数回解く．得られた解から残差を計算し，その残差を粗いグリッドへ射影する．次に，射影された残差について粗いグリッドで完全に解く．得られた解を細かいグ
リッド上の解に補正として与え，細かいグリッド上で再度数回解く．得られた解から残差を計算し，十分な収束が見られるまで処理を繰り返し，最終的な解を得る．実際には上記の処理だ
けでは 2 グリッドしか用いていないのでマルチグリッド法とは呼称されない．修正方程式解法に再帰的にマルチグリッド処理を適用することで初めてマルチグリッド法と呼称される．

21

今回の計算では，SOR 法を smoother としてグリッドを 3 枚使用し，残差が 1×10-5 以下になるまで計算を繰り返した．

図 2.3.1 V サイクル方式
細かい格子から順に計算を行い，最終的に細かい格子で計算を終える

22

図 2.3.2 マルチグリッド法の流れ図

23

2.3.3 不等間隔格子
プラズマ密度分布の算出にはシースの解析が重要であるため，壁面付近の格子サイズはデバ
イ長以下を要求される．小型中和器の放電室を直径 20mm，長さ 30mm 程度としたとき，デバイ長は 23.5μm である(電子密度 1×1017m-3，電子温度 2eV)．このデバイ長は全体のスケー
ルに対して 1200 分の 1 であり，全体を一様にこの精度で解くことは大幅な計算速度の低下をもたらす．そこで，場所ごとに適切な格子サイズを設定できるように不等間隔格子を採用した．
使用した不等間隔格子は 2 つの 格子サイズを持ち，それらのサイズの比は 1:4 となる．これらの格子のうち細かい格子を壁面付近に優先して配置することで，計算コストは低減される．
不等間隔格子を用いた場合のポアソン方程式解法手順は以下のようになる．
① 粗い格子で全体についてマルチグリッド法で解く．
② 粗い格子で得られた解を元に，細かい格子で壁面付近を局地的に解く．
①，②をそれぞれ別個にマルチグリッド法で解いたため，使用した 3 枚の格子サイズは①，②で異なる．格子サイズをまとめたものを下表に示す。最も細かい格子は，計算時間の問題からデバイ長と同程度とした．

表 2.3 格子サイズ
粗い格子

細かい格子

100 μm

25 μm

マルチグリッド法で200 μm

50 μm

使用される仮想格子

400 μm

100 μm

実際に物理量を割り振る格子

2.3,4 電位の境界条件
壁面の電子を零，中心軸上では軸対称性から電位の傾きが零となるように境界条件を設定した．また，中和器のオリフィス開口部の電位も零とした．実際にイオンスラスターとして中和器を動作させた場合，オリフィスは壁面の電位より高い電位を持つと考えられる．今回の計算
では，オリフィス簡素化のために零と置いた．

2.4 全体の流れ

24

計算の全体の流れ図を図 2.4 に示す．

図 2.4 全体の流れ図

2.4.1 タイムステップ
レーザートムソン散乱法を用いた放電室内の電子温度測定実験より，電子の平均エネルギーは数 eV 程度になると予想されたのでタイムステップを 7×10-10 sec として計算を行った．
2.4.2 初期条件
初期条件として，電子，イオン共に放電室内に一様に分布させた．また，速度分布はマクス
ウェル分布に従うとし，電子温度は 2 eV，イオン温度は室温(300K)とした．

25

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

Mersenne Twister Home Page :(http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html)(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)

増井博一， “マイクロ波放電型宇宙推進用プラズマ源に関する数値解析”, 平成 17

年度九州大学博士論文．
(2-15)

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

新谷將 “マイクロ波放電式小型イオンエンジンの数値解析”, 平成 22 年度九州大学

修士論文．
(2-17)

Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June Donato, Jack
Dongarra, Victor Eijkhout, Roidan Pozo, Charies Romine, Henk va der Vorst, ：“反復法 Templates” (朝倉書店, 1996)

26

第3章
3.1

計算コード評価

マルチグリッド法の妥当性検証

計算速度の向上に対して，マルチグリッド法が有効であるかどうか検証を行った．今回使用
した smoother である SOR 法単独の解法とマルチグリッドを用いた解法の計算精度，速度について比較した．
計算体系として 1 辺の長さが 1 の正方形の領域を設定し，下記の境界条件を与えた．

 ( x,0)  sin(x) (0, y)   (1, y)   ( x,1)  0
この時，理論解は以下の式で与えられる．

 ( x, y) sin x sinh 1  y 
sinh  

計算格子は等間隔格子を用い，メッシュ数は 200 とした．また，SOR 法の緩和係数は 1.4 と設定した．マルチグリッド法の Smoother に使用した SOR 法についても同様である．解の最大値が 1 であることから，残差ベクトルの大きさが 1×10-5 以下となった場合を収束条件とした．
図 3.1 に，各格子点における SOR 法，マルチグリッド法それぞれの解と理論解との誤差率を示す．SOR 法，マルチグリッド法ともに 0.36%以下の誤差で計算を行うことができた．最
大誤差率に若干の違いが見られるが，これは収束条件を 1×10-5 以下としたために反復回数にずれが生じたことに由来する．この結果より，SOR 法，マルチグリッド法ともに必要な計算精度で計算を行うことが可能であることがわかる．

図 3.1 理論解との誤差率

27

表 3.1 に上記の計算を行った際の計算時間を示す．マルチグリッド法を用いたことにより，計算時間は約 9 分の 1 となり大きく減尐した．この結果から，SOR 法を単独で用いるのではなく，マルチグリッド法を適用することによって計算コストの低減が望めることがわかる．

表 3.1 計算時間
SOR 法

マルチグリッド法

19.96 sec

2.375 sec

28

第4章

結果と考察

4.1 計算体系
実験に使用した中和器について計算を行った．図 4.1 に計算体系を示す．図の中和器は直径21 mm，長さ 32 mm の放電室を有しており，直線状アンテナでマイクロ波を照射する．アンテナの直径は 1 mm，長さは 29 mm とした．マイクロ波周波数は 2.45 GHz とし，投入電力
は 8W を仮定した．放電室の周囲と内部には棒状の Sm-Co 磁石が設置してあり，これらの磁石で磁気回路を形成する．

図 4.1 計算体系
赤色部について計算を行った．

4.2 磁場
使用した計算領域は図 4.2 に示すような磁力線分布を持つ．

図 4.2 磁力線分布
緑色部は ECR 層29

磁石とヨークの手前にそれぞれ磁気ミラーがあり，電子は磁石-ヨーク間で往復運動すると予想される(以降,この領域を磁気チューブと呼ぶ)．今回使用したマイクロ波周波数は 2.45GHz なので，87.5mT の領域が ECR 層となる．図 4.2 の磁場強度分布では二つの ECR 層が存在し，ECR 層で加速された電子が往復運動をすることによって，磁気チューブ内でプラズマを生成することが期待できる．

4.3 電子のエネルギー獲得確率
今回 FDTD-PIC コードでは中和器内部の領域を 3 領域に分割し，それぞれの領域におけるECR 層を通過した電子のエネルギー増分を調査した．分割領域を図 4.3.1 に示す．z 軸 16mm
以下を領域 1 とし，16 mm 以上の領域さらに 2 分割，r 方向の 5.25 mm 以下を領域 2，残った領域を 3 とした．エネルギー獲得確率の妥当性を示すため，アンテナ長を 25mm，29mm，33mm と変更して比較計算を行った．FDTD-PIC コードから得られた電子のエネルギー獲得
確率を図 4.3.2-4 に示す．アンテナ長に関わらず，領域 3 におけるエネルギーの変化は見られなかったので図には載せていない．領域 1，領域 2 においてアンテナ長に関わらず，電子のエネルギーは増加する確率が高いことが確認でき，ECR によるエネルギーの受け渡しが起きていることがわかる．
25mm から 33mm へアンテナ長を変更することは，アンテナを領域 2 の ECR 層に近づけることと同義なので，領域 2 の ECR 層周辺の電場強度は増加するはずである．図より，アンテ
ナ長を長くするに従って，領域 2 で高いエネルギーを獲得する確率が上昇していることがわかる．このことから，FDTD-PIC コードによるこの計算結果は妥当なものといえる．

図 4.3.1 分割領域30

図 4.3.2 アンテナ長 25mm

図 4.3.3 アンテナ長 29mm

図 4.3.3 アンテナ長 33mm

31

4.4 電子密度分布
図 4.3.1 にマイクロ波投入電力 8 W，推進剤流量 0.02 mg/s として 10μsec まで計算を行った時の電子密度分布を示す．磁気チューブ内に電子密度の高い領域が発生し， 密度は 5×1017～1×1018 m-3 程度となったことが図 4.4.1 よりわかる．レーザートムソン散乱計測法による密度測定実験(1-1)(1-2)では，イオンエンジン内部の密度は 1×1018 m-3 程度 と計測されており，計算された密度は実験値に近い値であることがわかる．また，図 4.2 と図 4.4.1 を比較すると，磁力線に沿った密度分布が確認できる．このことから，ECR 層で加速された電子が磁気ミラー間を往復運動した結果，このような分布を形成したといえる．可視化イオン源を用いた実験(1-3)では(図 4.4.2)，図 4.4.3 のようなプラズマ発光強度分布を得ている．図 4.4.3 では磁力線に沿った領域に強い発光強度を得ており，今回得られた電子密度分布と合致している．

図 4.4.1 電子密度分布32図 4.4.2 可視化イオン源

図 4.4.3 プラズマ発光強度分布4.5 イオン密度分布
図 4.5 にイオン密度分布を示す．イオンも電子と同じように磁気チューブ内で高い密度分布を持ち，7×1017～2×1018 m-3 程度の密度となっている．磁場に閉じ込められた電子に引き寄せられる形でイオンもこのような分布をとったと考えられる．また，壁面から離れた領域ではイオン密度が粗く分布しているように見える．これは，不等間隔格子によって粗く計算された電位の影響が現れたものである．電子密度分布でこの粗さが現れないのは，イオンの運動が電位と初速のみに依存するのに対して，電子の運動はそれらに加えて中性粒子との衝突の影響があるからだと考えられる．

図 4.5 イオン密度分布33

4.6 電子のエネルギー空間分布
電子の平均エネルギーの空間分布を図 4.5 に示す．ECR 層から離れた点では精々10～20 eV程度だが，ECR 層付近では高エネルギーの電子が存在することが図より分かる．ECR によるエネルギーの受け渡しによって，ECR 層付近の電子は高いエネルギーを持つようになったといえる．また，全電子の平均運動エネルギーは 19.1 eV となった．
レーザートムソン散乱計測法を用いた実験によると，マイクロ波放電式小型イオンエンジン
内の電子のエネルギーは 5-6 eV 程度であるとされており，イオンエンジンに非常に近い体系
である中和器内においても同程度のエネルギーを持った電子が観測されるはずである．電子エ
ネルギーは分布を持っているので，単純に実験値と計算値を比較することはできないが，多くの電子は実験値よりかなり高いエネルギーを持っているといえる．全体的に実験値より高いエネルギーを持った原因として，FDTD-PIC コードによるエネルギー獲得確率分布に問題があった可能性が考えられる．今回，FDTD-PIC コードでは中和器全体を 3 領域に分割し，ECR 層付近の電子のエネルギー増分を観察した．この領域の分割数が尐なすぎたために，マイクロ波のエネルギーが実際には届いていない領域においても，電子にエネルギーが与えられていた可能性がある．解決策は，FDTD-PIC コードの領域分割数を増やすことだが，分割数を増やしたことによる計算コストの増加が問題となる．計算コードの並列化や計算機の変更による計算速度の向上を図った上で，領域分割数を増加させた計算を行う必要性がある．
.

図 4.6 電子のエネルギー空間分布34

4.7 電位分布
電位分布を図 4.7.1 に示す．図 4.4.1、図 4.5 と比較すると，プラズマ密度が高い領域において 15～20 V の電位となっていることがわかる．壁面付近の電位分布をピックアップしたものを図 4.7.2 に示す．図 4.7.2 より，壁面付近の電位は 0.1 mm 程の間に急激に上昇していることがわかる．これはシースだと考えられる．

図 4.6.1 電位分布

図 4.6.2 壁面付近の電位分布354.8 電子損失
図 4.8 に 10μsec 時の各部位に流れる電子電流密度を示す．図 4.8 より，磁気ミラー周辺において電子の損失が大きいことがわかる．内壁に入射する電子はほとんどなく，磁気ミラーに近い磁石とヨークで損失が大きいことから，電子損失の改善において，磁石，ヨークの改善を優先して行うべきであるといえる．例えば，今回の計算条件では磁石，ヨークの電位は 0V を設定しているが，もっと高い電位を付与することでプラズマの損失を低減できるだろう．
また，アンテナでの損失が内壁を超えている．これは，やはりアンテナが磁気チューブに近いためであり，今回の中和器で電子損失を改善するのであれば，内壁よりもアンテナに重きを置くべきであるといえる．

図 4.8 部位ごとの電子電流密度364.9 電子損失の時間変動
図 4.9 に時間ごとの損失電子電流値を示す．初期条件として一様なプラズマを仮定しているので，計算開始直後は多くの電子が壁面入射している．しかし、3 μsec 程経過すると，内壁に入射する電子を磁石とヨークが上回っている．8 μsec 以降に注目すると，磁石，ヨークに入射する電子が増加していることがわかる．このことから，本計算は収束しておらず，プラズ
マの定常状態を算出するには更なる計算が必要であることがわかる．

図 4.9 各部位に流れる電子損失電流374.10 壁面に入射するイオンのエネルギー分布「はやぶさ」に搭載された中和器は，イオンによってスパッタされた壁面の Fe 粒子が再堆積することで性能が低下したと考えられている(1-4)．そこで，壁面に入射するイオンのエネル
ギー分布を今回作成したコードを用いて調査した．図 4.10 に壁面に入射するイオンのエネル
ギー分布を示す．以下に得られた知見を示す．
① ピーク
アンテナ，内壁，バックヨークの分布は 15-20eV の区間で強いピークを示し，磁石はそれらよりも低いエネルギーでピークが出ていることがわかる．アンテナ，内壁，バックヨークに入射するイオンは先述の電位分布において 15-20V の領域で生成されたものであるということができる．対して，磁石に入射するイオンにはそれらの電位より低電位の領域で生成されたイオンが多く含まれている．磁石に非常に近い領域で高いイオン密度が見られることから，これは磁石付近のシースでプラズマが生成されているためだと考えられる．磁石の分布において，15-20eV のイオンが殆ど見られないのは，15-20V で生成されたイオンがまだ磁石に到達していないからであり，計算を進めていけばそれらのイオンが観察されるようになるだろう．
② 低エネルギーイオン
内壁，磁石，ヨークのエネルギー分布図において，わずかに低エネルギー(数 eV)のイオンが確認できる．また，低電位の領域と接していないアンテナでは，低エネルギーのイオンは見られない．これは初期条件により一様に分布していた低エネルギーのイオンが，低電位の領域から加速されずに入射したものだと考えられる．すなわち，この現象は計算の初期条件によるものであり，定常状態に近づくにつれてこれらのイオンは姿を消すだろう．

将来的には，これらのエネルギー分布を使用した Fe スパッタのシミュレーションコードを作成し，寿命評価を行いたいと考えている．

38

図 4.10 壁面に入射するイオンのエネルギー分布

参考文献
(1-1)

N Yamamoto, K Tomita, N Yamasaki, T Tsuru, T Ezaki, Y Kotani, K Uchino, H

Nakashima, “Measurements of electron density and temperature in a miniature
microwave discharge ion thruster using laser Thomson scattering technique”, Plasma
sources science and technology, 19(2010)．
(1-2)

杉田 健策， “小型マイクロ波放電式イオンスラスタのマイクロ波周波数依存性”，

平成 23 年度九州大学修士論文．
(1-3)

小谷 優介， “可視化イオン源を用いたイオンスラスタの内部物理診断”， 平成 21 年度九州大学修士論文．
(1-4)

大道 渉， 小泉 宏之， 国中 均, 細田 聡史， 西山 和孝, “マイクロ波放電式中和器の性能低下メカニズムとその実験的検証”，STEP-2010-058

39

第5章

結論

本研究では，マルチグリッド法によって計算コストを低減したプラズマシミュレーションコ
ードを作成し，プラズマ損失を厳密に解析した場合のプラズマ密度分布を示した．結果より，中和器内部のプラズマ密度生成は磁気チューブ内で盛んであり，チューブ内において高い密度
でプラズマが分布していることがわかった．また，磁気チューブに近い磁石，ヨークで多くのプラズマが失われていることが示された．また，各部位に入射するイオンのエネルギー分布を示すことができた．
今回開発したコードを用いて計算を進めていくことで，定常状態のプラズマ密度分布を算出することができる．それにより，最適な設計を計算機上で見出すことが可能となる．しかしながら，一様分布のプラズマから計算をスタートさせた場合，10μsec でも計算が収束しておらず，実時間における計算時間は 6.5 日程度となる．イオンが中和器の端から端まで移動した場
合，100μsec 程度の時間が必要であることから，プラズマの定常状態を算出するには今回の計算の十倍，すなわち 65 日の日数を要すると考えられる．この日数を実用的とするかどうかはユーザーによるが，さらなる計算コストの削減が可能であるならば，それは望ましい．実際，
今回作成したコードは並列計算の余地があり，更なる高速化の可能性を秘めている．今後は，更なる計算コストの削減を行い，定常状態を目指すことが重要となるだろう．

40

第6章

今後の展望

5.1 電子質量
今回電子の質量を実際の 2500 倍に増加させ計算を行った．電子の質量の増加によって，プラズマ電位は実際よりも低い値が算出されていると考えられるが，その程度は定かではない．
プラズマ電位が実際と異なる場合，シース，延いてはプラズマ密度分布も変化していると考えられるので，質量の変更が計算結果にどの程度影響を与えているのか調査が必要である．

5.2 FDTD-PIC コードの初期条件
FDTD-PIC コードの初期条件として，一様なプラズマ分布を用いた．従って，今回の計算結果である損失を考慮したプラズマ密度分布を初期条件とすることで，より厳密なエネルギー獲
得確率分布を算出することが可能となる．プラズマの定常状態を見るためには，計算のリスタートが必要であることがわかっているので，リスタートの際に FDTD-PIC コードのデータを更新していくと効率良く定常状態のプラズマ密度が計算することができる．

41

謝辞
本研究を進めるにあたって，多くの方々にお世話になりました．宇宙推進器の研究に携わる機会を与えてくださり，また，時として歩みの遅い研究内容に対して辛抱強くご指導してくださった中島秀紀教授に深く感謝いたします．また，研究内容について快く議論に応じ，数多くの助言を下さった山本直嗣准教授に深く感謝いたします．不慣れな事務処理についてサポートしていただいた大神めぐみ秘書，馬渡隆子秘書に深く感謝いたします．
研究室の皆様には公私ともに大変お世話になりました．皆様のおかげで充実した楽しい 2 年間を過ごすことができました．富士山や宝満山への登山は実に良い思い出です．また，研究の面でも多くの方々からサポートや助言をいただいたことに深く感謝しております．この場を借りてお礼申し上げます．

42