第1章

序論

1.1 背景1957 年にロシア（旧ソ連）が世界初の人工衛星「スプートニク 1 号」を打ち上げてから宇宙開発を取り巻く状況は大きく変化した．当初は，アメリカとソ連による軍事技術への転用
や国際的イメージ向上を目的としたものが中心であったが，70 年代に入ると欧州と日本が宇宙開発に加わり，国際協力の下で宇宙開発が行われるようになった．現在までに5000 個以上もの数の人工衛星が打ち上げられており，衛星通信システム，衛星放送シス
テム，気象衛星システム，地球観測衛星システム，化学観測衛星システム，カーナビゲー
ションやスマートフォンの GPS など身近にも多く利用されている(1-1)．このため，人工衛星
は私たちの生活に密接にかかわり欠かせないものとなっている．
現在，人工衛星の軌道制御や姿勢制御を行うために推進機が用いられている．推進機
は大きく分けて化学推進機と電気推進機とがある．電気推進機は，電気エネルギーを用いて推進剤を電離させ，生成されたプラズマを様々な方法で加速し，その反作用で推力を得る．これに対し，従来の化学推進機は化学反応を利用し作動流体を噴出し推力を得る．電気推進機は化学推進機と比較して得られる推力ははるかに小さいが，燃費を示す指標である比推力(Specific impulse：Isp)を一桁以上大きくすることが可能である．つまり，電気推進機は推力が小さいため加速に時間がかかるが，消費する燃料は化学推進機と比較し非常に尐なくて済む．
また，電気推進機の中にもいくつか分類があり，イオンスラスタやアークジェット，MPD
スラスタ，RF スラスタ，パルスプラズマスラスタ，ホールスラスタなどが開発・研究されている．特に，イオンスラスタは他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制御や惑星探査に適している．
この中でもマイクロ波放電式イオンスラスタは放電電極を使用せず推進剤にマイクロ波を放射することで放電を行うため，スラスタの構造が簡素化でき耐久性を向上させることができる．このマイクロ波放電式イオンスラスタは，小惑星探査機「はやぶさ」にも用いられ大きな注目を受けた．「はやぶさ」に搭載されたイオンエンジンμ10 は，一つのシステ
ムとしての積算作動時間は 4 台で 40,000 時間に迫り，世界一の宇宙実績と高い信頼性を示した(1-2)．しかしイトカワから地球への帰還途中，単独で 14,000 時間以上という最長作動実績を上げていたイオンエンジンの中和器が劣化の兆候を見せた (1-3)．中和器とは電子を放出することでイオンビームを中和する装置である．イオンエンジンは，中和器無しで用いるとイオンビームの逆流を招いてしまい，推進器としての役割を失ってしまう．現段階では中和器の性能低下がイオンエンジンの寿命の決定要因となっており，中和器を長寿命化することが求められている．

1

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

図 1.1 イオンエンジン概略図

1.2.1 マイクロ波放電式イオンエンジン
イオンスラスタはプラズマの生成方法によって直流放電式，RF (Radio Frequency) 誘導
放電式，マイクロ波放電式等に大別される．本研究の対象であるマイクロ波放電式イオンエンジンの概念図を図 1.2 に示す．マイクロ波放電とは，マイクロ波帯域の交流電場によ2

って電子を加速し，中性粒子と衝突電離により電子の数が増加し気体がプラズマ化され放電が維持されることである．

図 1.2

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

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

  ev  B

dt(1.2.1)

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

mv 
eB

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

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

E

B

電子

E

B

電子

電子
E

B
B
電子

E

図 1.3 ECR の原理

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

1.
2.

ヒーターを用いないため，予備加熱が不要となり素早い起動が出来る．
マイクロ波放電式イオンエンジン駆動用のマイクロ波電源が 1 台あれば中和器のプラズマも同時に生成することが可能なため，システム構成の簡素化や信頼性の向上，また衛星の重量軽減が図れる．

マイクロ波放電式中和器は深宇宙探査機「はやぶさ(1-5)」に搭載され十分な中和性能が得られたという実績がある．当研究グループでもマイクロ波放電型イオンエンジンを開発してきた． 図 1.4 はマイクロ波放電型中和器の外観写真である．

図 1.4 マイクロ波放電式中和器の写真

5

図 1.5 マイクロ波放電式中和器運転時

6

1. 3 マイクロ波放電式中和器の性能低下
1.3.1 スパッタリング(1-6)
イオンによるスパッタリングは金属ターゲットの表面にエネルギーを持ったイ.オンが衝突することで発生する．このイオンは入射イオンと呼ばれ，これと直接衝突したターゲット原子は一次反跳原子と呼ばれる．この一次反跳原子は通常元の格子位置から弾かれたら他のターゲット原子と衝突を起こす．この衝突の連鎖は衝突カスケードと呼ばれ二次，三次と反跳原子が他の原子を弾けるだけのエネルギーを持つ限り続く．最終定期にターゲット原子がスパッタされて外に飛び出すためには原子がターゲット表面付近で外部に向かうような軌道で弾かれ，且つその原子のエネルギーが表面結合エネルギーを超える必要がある．この衝突カスケードの様子を模式的に示した者を図 1.6 に示す．スパッタリング
は通常スパッタ率(Y )で評価される．これは一つの入射イオン当たりいくつのスパッタ粒子が発生するかを示すもので図 1.6 の場合は一つのイオンが入ってきたのに対し外側に出ている原子の数が三つなのでスパッタ率 Y＝3 となる．
スパッタ率はイオンの質量，エネルギー，入射角，ターゲット原子の質量，原子構造，表面結合エネルギーに依存する．入射イオンの質量が増えるにつれスパッタ率は増加する．
また，エネルギーに関していえばイオンのエネルギーがスパッタリング閾値以下であるとスパッタリングを起こさず，これを超えると増加する．しかしながら，0.1 eV あたり超えると入射イオンがターゲットの表面に対し深く入り込みすぎるため表面から粒子が出にくくなりス
パッタ率は減尐する．次に角度について述べる．スパッタ率は入射角がターゲット表面の法線に対する角度が 60°になるまで角度と共に増加しその後は 90°で 0 になるまで減り続ける．
Sigmund はスパッタリングを図 1.7 に示すようにシングルノックオン型，線形カスケード型，スパイク型の 3 タイプに分類した．これは入射イオンのエネルギーによって大まかに分類できる．まず 1 eV 以下程度の低エネルギーで起こるシングルノックオン型スパッタリングでは衝突カスケードは発生しない．この場合，最終的に外部に弾き出されるスパッタ粒子は元の入射イオンにより弾き出されるため極尐ない衝突しか起こさない．また，線形カスケ
ード型のスパッタリングは 1~300 eV 程度で起こり，その名の通り衝突カスケードを起こす．
この場合のスパッタ粒子は反跳原子によって弾き出されるが反跳原子と他のターゲット原子との間で衝突が起こるのは希である．100 eV 以上若しくは入射イオンのフラックスが大きい場合に起きるスパイク型のスパッタリングも衝突カスケードを起こし，またイオンのエネ
ルギーが高いため多くの衝突が起きる．特にこの場合の特徴として既に他の原子により弾かれ動いている原子に対して更に衝突するといったことが起きる．イオン推進機のような比較的低エネルギーのイオンを扱うものでは主にシングルノックオン型のスパッタリング
を起こす．

7

図 1.6 衝突カスケード(a)シングルノックオン型(b)衝突カスケード型(b)スパイク型< 1 keV
1 keV~300 keV>100 keV
図 1.7 Sigmund によるスパッタリングの分類

8

1.3.2 中和器電圧の上昇(1-3)
中和器の性能低下は中和器電圧の上昇という形で現れる．図 1.8 に宇宙機の電位の模式図を示す．通常は図 1.8 の(a)の様な電位構造を示す．電源のグランドは宇宙機構体に接続され，これが周辺のプラズマの電位と等しくなる．このプラズマの電位に対して，イオ
ンビームと中和電子がそれぞれ加速されて電流ループが閉じる．中和電流はイオンビー
ムと等しくなるように定電流制御され，その電流を引き出すのに必要な中和器電圧が自動的に調節される．しかし，劣化が進み電源の出力できるレンジを超えるとイオン電流と中和電流のバランスは崩れ中和器が機能しなくなる．その様になってしまうと，keV のエ
ネルギーを持つイオンビームが負に帯電した宇宙機の電荷に引き戻されグリッドや宇宙機表面の激しい損耗を引き起こす(図 1.8(b))．
この様な中和器電圧の上昇の原因の実験的検証は「はやぶさ」プロジェクト後，宇宙航空研究開発機構(JAXA)の宇宙科学研究所(ISAS)により行われた．これによるとスパッタリングが根本的な原因であるとされている (1-7)．すなわち以下のようなメカニズムを想定した．図 1.9 にその模式図を示す．中和器内部のプラズマはシースにより壁面に対して正の空間電位を持つ．このため，イオンは壁面へ向かって加速され衝突する．このイオンのエ
ネルギーがスパッタリングの閾値を超えるとスパッタリングを起こす．このとき発生した壁面
材料のスパッタ粒子は中和器内部を運動し最終的に再び壁面に付着する．この様なスパッタ粒子の一部はアンテナを絶縁している絶縁体上へと堆積する．これによりマイクロ波の伝搬モードの変化もしくは絶縁性能低下により，プラズマとマイクロ波のカップリングの低下をもたらし，結果としてプラズマ密度の低下が起こると考えられる．これに伴う電子ビームの減尐を防ぐため，中和器本体の中和器電圧を上昇させる (1-3)．しかしながら，この中和器電圧を上げると中和器内部の空間電位の上昇を招くため(1-8)，更なるスパッタリ
ングを誘発し悪循環となる．スパッタ率は電圧に対して指数的に増加することから，中和器電圧が高い状態で運転を続けると中和器の急激な劣化を招く(1-3)．スパッタリングのうち，特に二価のキセノンイオンのスパッタリングが中和器内部の損耗の主要因を占めている(1-7)．これは中和器内部の空間電位が 20~30 V 程度であるのに対し，ほとんどの壁面材
料のスパッタリング閾値が 30 eV 以上であるため一価イオンではほとんどスパッタリングの閾値を超えないためである．
また，図 1.10 に中和器を長時間作動させ続けたときの中和器電圧の変化を示す(1-7)．
中和器電圧の上昇が 1 万 5000 時間経過後に顕著に現れているは空間電位の上昇によるスパッタ率の指数関数的増加に起因すると考えられる．

9

(a) 通常の状態(b) 中和不良時の状態

図 1.8 中和器の電位の模式図

図 1.9 中和器の性能低下メカニズム10図 1.10 中和器電圧の時間変化 (大道 渉ほか, “マイクロ波放電式中和器の性能低下
メカニズムとその実験的検証”，STEP-2010-058)

11

1.4 研究目的
本研究室では，これまで実験や数値解析を用いて中和器の性能評価や内部の物理メ
カニズムの解明を目的として中和器の開発に注力してきた．しかし，中和器の寿命を直接的に評価することは行われておらずそのメカニズムはまだ解明されていない．中和器
損耗メカニズムを実験により確かめるためには，約 14,000 時間以上の時間と数千万円のコストがかかり実験により評価することは困難である．そこで，効率的に中和器の性能の変化を評価する数値解析コードを開発し，そのコードの妥当性を示し，本コードを用いて，小型イオンエンジン用中和器の寿命を評価することが本研究の目的である．
参考文献
(1-1)
(1-2)
(1-3)

茂原正道・鳥山芳夫, “衛星設計入門”, 2002.
大道渉，”マイクロ波放電式中和器のイオン電流測定”，STEP-2011-033
細田聡史，國中均，”イオンエンジンによる小惑星探査機「はやぶさ」の帰還運用”，J.Plasma Fusion Res. Vol.86,No.5,(2010),pp.282-292(1-4)
N. Wallace, P. Jameson, IEPC, 2011, 327(1-5)
I. Funaki, H. Kuninaka and K. Toki, J. Propul. Power, Vol. 20, (2004), pp.718-726.
(1-6)
Kirk A. Zoerb, “Differential sputtering yields of refractory by ion bombardment at
normal and oblique incidence”, Thesis, Department of Mechanical Engineering,
Colorado State University, Fort Collins, Colorado (2007)(1-7)
大道渉， 小泉宏之， 国中均, 細田聡史， 西山和孝, “マイクロ波放電式中和器の性能低下メカニズムとその実験的検証”，STEP-2010-058(1-8)
Noriyoshi Onodera,et al., “Electron Emission Mechanism of Microwave
Discharge Neutralizer”, IEPC-99-160,2002

12

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

 q E  v  B 

dtdx(2.1.1) v(2.1.2)dt

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



dt
dt



m

dv z
dt

r

z

r

r

z

v r v

(2.1.3)r

E  v B 

q



z

E  v B  v B  

q





r

m

dv 

2

v

E  v B  

q



z

m

r

これらを四次精度 Runge-Kutta 法にて解いた．
n 1

 vr 
nvrn 1

 v 
nv

n 1

 vz 
n

vz

t
6
t
6
t
6

b r 1  2 b r 2  2 b r 3  b r 4 b  2 b  2 b  b 
1

2

3

4

(2.1.4)

b z 1  2 b z 2  2 b z 3  b z 4 

但しここで定数 b はそれぞれ以下を用いる．
br1 

q
m

E  v B   1 v v
r



z



r



br 2 

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

br 3 

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

br 4 qm

E  v   tb B   1 v   tb v   tb 
r



3

z

r



3



13

3

(2.1.5)b 1 qm

E  v B  v B  


z

r

r

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

qm
qm

E  v   tb B  v   tb B   1 v   tb v   tb 


z

z3

r

r

r3

z

r

z

z3

r

z



r

 tb  1  q  Br 
 E z   v  m 
2  
qmr3

E  v B (2.1.7)

 tb  2  q  Br 
bz3 
 E z   v  m 
2  
bz 4 

(2.1.6)

E  v   tb B 
z



3

r

一般的な実験室プラズマは 1 cm3 当りに 1010～1014 個が含まれ，この粒子全部を取り扱う
ことは最新のスーパーコンピュータを用いても不可能であり，実用的でない．そこで粒子シミュレーションで取り扱う粒子として，実際の粒子の電荷質量比を一定に保ったまま多数の粒子の電荷と質量を１つにまとめた超粒子 (2-2)を使用する．超粒子数は，１セルに含まれる粒子数が 100 個以上になる事と取り扱う粒子数よりも超粒子の重みが尐ない事（統計性から）から決定する．
2.1.2 メッシュサイズ・タイムステップ (2-3)
計算に用いるメッシュサイズ  は，電子温度 Te，ボルツマン定数 kb，電子密度 ne から計算される．
デバイ長 λDe，
 De 

 0 k bTe
nee

2

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

nee

2m e 0(2.1.10)に対して，
1

 t   pe

14

(2.1.11)

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

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

eB

(2.1.13)meに対して，
 t   c

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

Sn
S

n  1 ～ 4 (2.1.15)となり，

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

15(2.1.16)図 2.1.1 PIC 法におけるセルの構成

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

S 1 qNS

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

16

2.1.6 平均自由行程と透過度(2-6)
面  ，厚さ l ，その物質の数密度 N の板を考える．エネルギー E の粒子ビームが，面

 に垂直に当るものとする．衝突をせずに板を透過してくる粒子数は，入射粒子数より減尐する.その減尐度は  を衝突の断面積とすれば，dn   nN  dl

(2.1.18)n  n 0 exp  N  l (2.1.19)と表わされ，これよりとなる．ここで， n は面  （ l  0 ）でのビーム中の粒子である．それゆえ，（ l , l  dl ）間で0

最初の衝突が起こる確率は，
Pinitial

l dl  exp  N  l N  dl

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

(2.1.22)1N

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

v


  Nvとなり，これを衝突周波数という．

17(2.1.26)

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

(2.1.27)

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

totalがわかる．したがって，ある粒子 m の衝突確率 Pcollision, m は以下の式で表わすことができる．
 total  n gas  total v m(2.1.28)

Pcollision, m  1  exp  n gas  total v m  t (2.1.29)

 1  exp   total  t 

ここで， n gas は中性粒子の数密度， v m は粒子 m の速度， 

total

は全衝突断面積である．ある乱数 r1 （ 0  r1  1 ）を用いて， Pcollision, m  r1 のときに粒子 m はタイムステップ  t の間に衝突が起こるとする．
2.1.8 衝突計算
乱数 r 2 を用いて以下のように衝突の種類を決める．ここでは例として，電子と中性粒子の衝突について示す．弾性衝突断面積を 
断面積を 
 elastic
 total

ここで， 

 elastic

 r2 

 totaltotal

excitation，電離衝突

であれば弾性衝突

 total elastic   excitation

 elastic   excitation

，励起衝突断面積を とする．

ionization

0  r2 elastic

であれば励起衝突 total r2   elastic   excitation

 elastic   excitation   ionization  ionization

 totalであれば電離衝突である．このように断面積の比と乱数を対応させ，衝突の種類を決める．計算では中性粒子としてキセノンを用い，その断面積データ(2-7)(2-8)(2-9)

を図 2.1.2 に示す．今回の計算では，多価電離は考慮していない．

18

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

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

1v 2 

12

2

v 1  v 2  u  (2.1.30)

v 1  v 2  u  (2.1.31)と表わすことができる．ここで， u   v 2  v 1 （相対速度）である．電子と中性粒子の衝突も同様に考える．衝突前の電子，中性粒子の速度をそれぞれ v ・ v n の衝突後の速度をそe

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

m  M cos m  Mh の直交座標成分は

19

ve 

M sin m  M

h(2.1.32)

h x  v er cos hy  
hz  

(2.1.33)v ex v ey cos   v e v ez sin (2.1.34)v er
v ex v ez cos   v e v ey sin (2.1.35)

v er

ここで，
v er 

ve v ey  v ez
22

(2.1.36)v ex  v ey  v ez
222

(2.1.37)であり，  は散乱角，  は任意の角である．(2.1.32)式はエネルギー損失を含んでいる．
次に，衝突前後のエネルギー収支を考える．衝突前の電子のエネルギーを E
衝突後の電子のエネルギーを E

scattered, eincident, eとし，とすれば，エネルギー収支は以下のように表わされる．
E scattered, e  E incident, e(2.1.38)

電子は，エネルギー E incident, e が小さいとき等方散乱され，大きいときは主に前方散乱となる．これを考慮した散乱角  の確率密度を g    として以下に示す (2-19)．ここで， E

incident, eの単位を eV とする．
g   

E incident ,e sin 0     
2 
2 1  E incident ,e sinln 1  E incident ,e 2

これを解けば，cos   1 (2.1.39)2 1  1  E incident, e  3r

(2.1.40)E incident, eとなる．この値を(2.1.32)式に用いて衝突後の速度を求める．このとき，任意の角は  2  r4

(2.1.41)

として与える．
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーをE excitationとし，入射電子のエネルギーを E

incident, e，励起後のエネルギーを E

scattered, e，励起後の速度を v~ とする．エネルギー収支は，
E scattered, e  E incident, e  E excitationと表わされ，励起後の速度は

20(2.1.42)v~  v e

1E excitation(2.1.43)

E incident, eとなる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失った
後に(2.1.43)式中 v~ の速度で弾性衝突をすると考えることができる．(2.1.32)式中 v に v~を用いて，衝突後の速度 v  を求める．このときの散乱角には(2.1.40)式を用いるが
ee(2.1.40)式中の E incident, e には(2.1.42)式から求められる E scattered, e を用いる．励起エネルギーEは，キセノンの場合 8.34 eV を用いている．
excitation③電離衝突
入射電子のエネルギーを E
ルギーを E

created, eincident, e，散乱電子のエネルギーを E，電離エネルギーを E

ionization

scattered, e，生成電子のエネとする．エネルギー収支はE scattered, e  E created, e  E incident, e  E ionization(2.1.44)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる(2-12)．
E incident, e  E ionization-1 
E scattered, e  B tan  r5 tan 
2B








(2.1.45)

ここで，B はガスに固有の値でキセノンでは 8.7 eV となる．
また，入射電子が電離によって失うエネルギーを  E とすると，
 E  E incident, e  E scattered, e(2.1.46)と表わされる．(2.1.43)式から電離後の電子の速度は，v~  v e1E(2.1.47)E incident, eとなり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．(2.1.34)式中 vに v~ を用いて，衝突後の速度 v  を求める．散乱角には(2.1.40)式を用いるが，(2.1.40)式
中の Eには(2.1.45)式から求められる Eを用いる．電離エネルギー E
は，キセノンの場合 12.13 eV である．
ee

inc ide nt, eionization

scattered, e

電離後の生成電子の持つエネルギーは，E created, e   E  E ionization(2.1.48)となる．また，入射電子の速度 v を用いて生成電子の速度はev created, e ve

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

を用いて(2.1.40)式から求める．

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

計算時間短縮のため，Null-collision 法を用いる．これには，架空の衝突断面積 

fake

を

導入する．すべての電子のエネルギー E に対して
 fake   total(2.1.50)であるとし，例えば電子と中性粒子の衝突においては total   elastic   excitation  ionization  fake(2.1.51)とする．ここで， total v を満たす 

fakeを導入すれば， collision collision

 constant(2.1.52)

n gas

が一定，すなわち Pcollision, m が一定となりエネルギーに依存しない．したがって，粒子毎に Pcollision, m を計算する必要が無くなり，全粒子に対して衝突計算を行わず，衝突を起こす粒子のみについて衝突の種類を決めればよい．これにより，計算時間を大幅に短縮できる．衝突の種類を決めるときは，2.1.8 において，
 total   elastic   excitation  ionization  fake とし， elastic   excitation   ionization
 total r2  1のとき Null-Collision Processとする．
Null-Collision Process が選ばれたときには，粒子の運動は何ら変化しないものとする．
図 2.1.3 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度で規格化した衝突周波数を示す．本計算では Null-Collision 法で用いた
m3/sec である．

22

 collisionn gasの値は 2×10-12図 2.1.3 Xe の数密度で規格化した電子の Xe の衝突断面積
2.1.11 電子質量の変更
今回の計算において，電子の速度成分は最大で 105 m/s 程度，イオンにおいては最大
103 m/s 程度であり，電子の速度を基準に時間刻みを設定するとイオンはほとんど動かない．そこで，電子の質量を実際のものよりα倍大きくした．質量の変更により，電子の速度は√α分の 1 に低下，イオンとの速度差は小さくなり計算結果の収束を早めることができる．
m e _ virtual   m e _ real

(2.1.53)

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

mv
eB

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

2.2 マイクロ波解析手法(2-14)(2-15)マイクロ波の伝搬解析は，別途 FDTD-PIC コード用いて行った．FDTD-PIC コードとは，
FDTD 法（Finite Difference Time Domain method）と PIC 法を組み合わせたカップリングコードであり，電子の挙動とマイクロ波の伝搬を同時に解析できる．このコードを使用して，ECR 層を通過する電子のエネルギー増分を調査，ECR 層を通過する際に電子が得るエ
ネルギーの確率分布を求め，プラズマ粒子解析に導入した．
2.2.1 基礎式
FDTD 法は数値電磁場解析の分野において，モーメント法や有限要素法と同様に，アンテナや平面回路など数多くの複雑な問題に適用され，主要な解析手法として多く用いられている．計算には，マクスウェル方程式を基礎式として用いる．以下に基礎式として用いる微分形式のマクスウェル方程式を示す．
  E  0  H  0

E
t

H

(2.2.1)t J(2.2.2)

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

0

2.2.2 差分式
三次元の FDTD 法で用いられる電場，磁場の空間配置を図 2.2.1 に示す．このように，
電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．これは Yee によって提案され(2-16)，電場の回転（   E ）が磁場を，磁場の回転（   H ）が電場の方程式を満たすように配置されている．
この電磁場の配置のもと， (2.2.1)式, (2.2.2)式を空間，時間とも中心差分した式を以下に示す．((2.2.3)~(2.2.8))式ここで，  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) +Dt
m0 (i, j +1 2, k +1 2)

ìï Eyn (i, j +1 2, k +1) - Eyn (i, j +1 2, k ) Ezn (i, j, k +1 2) - Ezn (i, j +1, k +1 2) üï
+
í
ý
Dz
Dy
ïî
ïþ
1
1
n+
nDt
H y 2 (i +1 2, j, k +1 2) = H y 2 (i +1 2, j, k +1 2) +m0 (i +1 2, j, k +1 2)
ìï Ezn (i +1, j, k +1 2) - Ezn (i, j, k +1 2) Exn (i +1 2, j, k ) - Exn (i +1 2, j, k +1) üï
+
í
ý
Dx
Dz
ïî
ïþ
24

(2.2.3)(2.2.4)

n+

1

n-

1

H z 2 (i +1 2, j +1 2, k ) = H z 2 (i +1 2, j +1 2, k ) +Dt
m 0 (i +1 2, j +1 2, k )(2.2.5)

ìï Exn (i +1 2, j +1, k ) - Exn (i +1 2, j, k ) Eyn (i, j +1 2, k ) - Eyn (i +1, j +1 2, k ) üï
+
í
ý
Dy
Dx
ïî
ïþ
n 1

Ex

i  1 2 , j , k   E xn i  1 2 , j , k  t 0 i  1 2 , j , k 1n
 n 1
 H z 2 i  1 2 , j  1 2 , k   H z 2 i  1 2 , j  1 2 , k 

y


n



n 1

Ey

1

n

1

H y 2  i  1 2 , j , k  1 2   H y 2 i  1 2 , j , k  1 2 zi , j  1 2 , k   E yn i , j  1 2 , k  (2.2.6) J x  i  1 2 , j , k t

 0 i , j  1 2 , k 1n
 n 1
 H x 2 i , j  1 2 , k  1 2   H x 2 i , j  1 2 , k  1 2 

z


n



n 1

Ez

1
2

Hz

n

1

i  1 2 , j  1 2 , k   H z 2 i  1 2 , j  1 2 , k xi , j , k  1 2   E zn i , j , k  1 2  (2.2.7) J y  i , j  1 2 , k t

 0 i , j , k  1 2 

1n
 n 12
2
 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
y(2.2.8)2
i , j  1 2 , k  1 2  J z  i , j , k  1 2 

1((2.2.3)~(2.2.8))式から，時間ステップがn  1 2 での磁場 H
1

n

1
2

は，その半ステップ前でのn

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

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

26

2.2.3 アンテナ解析
本計算ではアンテナを完全導体で構成されているものとしてモデル化する．給電方法
は同軸ケーブルを模擬し，図 2.2.2 に示すように動径方向の４つの電場を中心導体に与える同軸線路給電を用いる.また, アンテナ導入部には図 2.2.3 のように Mur の一次吸収境界条件を与え, 波長に比べて太さが無視できるくらい十分細い直線状のアンテナを解析する場合，その中心軸に沿って電界を 0 とすればよい．また, 今回は L 字型アンテナを模擬した．給電方式について詳しく述べると，中心導体に与える電場の x 方向成分を Ex，
y 方向成分を Ey は近似的に
𝐸𝑥 =

𝑉
𝛥𝑥
𝑟𝑙𝑛 ( 𝑎 )(2.2.9)𝐸𝑥 =𝑉𝛥𝑦𝑟𝑙𝑛 ( 𝑎 )(2.2.10)で与えられる. ただし, V は中心導体と外部導体間の電圧, 𝑟 は動径方向の距離, 𝑎
はアンテナの半径である.
また, 同軸ケーブルのインピーダンス Z は以下の式で表される.
𝜇 1 𝑅
ln
𝜀 2𝜋 𝑟𝑍=√(2.2.11)

ここで, R は外部導体の半径, r は中心導体の半径である. マイクロは電力を P とすると，
𝑉2(2.2.12)𝑃=
2𝑍となり，(2.2.11)式，(2.2.12)式より，V を求め，(2.2.9)式，(2.2.10)式よりアンテナに与える電場を求める．ここで，Ex を 𝑟 = 𝛥𝑥/2の電場で，Ey を 𝑟 = 𝛥𝑦/2 の電場で代表させると，
Ex (IFED, JFED, KFED) = -E x (IFED -1, JFED, KFED)
V (nDt)=Dx æ Dx ö
ln ç ÷
2 è a ø(2.2.13)Ey (IFED, JFED, KFED) = -E y (IFED, JFED -1, KFED)=-

V (nDt)Dy æ Dy ö
ln ç ÷
2 èaø

ここで，式中の(𝐼𝐹𝐸𝐷, 𝐽𝐹𝐸𝐷, 𝐾𝐹𝐸𝐷)は電場を励振する座標を与える．

27(2.2.14)図 2.2.2 同軸線路給電

図 2.2.3 アンテナモデル

28

2.3 電場解析手法(2-17)
電場解析により電場を解きそれと共に電子とイオンに対し運動方程式を解くことでプラズマ解析を行った．このコードを電子，イオンの両方のプラズマ粒子を扱っていることから
Full-PIC コードと呼ぶ．
2.3.1 基礎式
ポアソン方程式を基礎式として計算を行った．
  



2(2.3.1)0(2.2.1)式を差分化し,連立方程式として解いた．差分式を以下に示す．
aa 1   i 1 , j  aa 2   i 1 , j  aa 3   i , j 1  aa 4   i , j 1  aa 5   i , j 
aa 1  aa 2 

aa 3 

1r

2

 i, j
0

1
z

2

1
2 ri , j  r(2.3.2)
aa 4 

1
r

2



1
2 ri , j  r

1 
 1
aa 5   

2
2 
r 
 z

ここで， i，j はそれぞれ z，r 方向の格子点番号に対応している．φi,j は i，j における電位を表し，ρi，j はその点での電荷密度を表す．⊿z，⊿r は各方向の格子サイズであり，ri,jはその格子点の r 方向における位置を表す．
2.3.2 反復法
(2.2.2)式はすべての格子点に１つずつ与えられるので，連立方程式として解くことができる．この連立方程式を解く方法として，反復法を用いた．
I 番目の方程式が以下のように表せる連立方程式を考える．
na x b
i, j

j 1

j

i

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

次の式から xi を求める．
 bi   a i , j x j 


ji


xi 
ai, j

(2.3.4)

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

(2.3.5)

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

(k )

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

ωを反復解の収束率を加速させるような値に選択することにより，ガウス-ザイデル法より早い収束を得ることができる．ωには最適値が存在するが，最適値予め知ることは一般的には不可能である．そのため，しばしばω=2-O(h)という発見的な評価式が使用される．
ここで h は元の物理領域を離散化するときのメッシュ間隔である．今回の計算では，このSOR 法に対し後述のマルチグリッド法を適用した．
2.3.3 マルチグリッド法
反復法の計算コストを低下させるテクニックとしてマルチグリッド法がある．Jacobi 法や
ガウス-ザイデル法のような単純な反復法は，誤差の高周波成分を最初に収束させることが知られている．この特徴を利用し，複数枚の粗さの異なるグリッド(空間格子)を用意して，それぞれのグリッドに対する高周波成分を収束させるのがマルチグリッド法である．
マルチグリッド法の計算の進め方は数多くあるが，今回はオーソドックスな V サイクル法を用いた(図 2.3.1)．図 2.3.2 に流れ図を示す．まず，細かいグリッドと粗いグリッドを用意し，細かいグリッドでほんの数回解く．得られた解から残差を計算し，その残差を粗いグリッド30

へ射影する．次に，射影された残差について粗いグリッドで完全に解く．得られた解を細かいグリッド上の解に補正として与え，細かいグリッド上で再度数回解く．得られた解から残差を計算し，十分な収束が見られるまで処理を繰り返し，最終的な解を得る．実際には上記の処理だけでは 2 グリッドしか用いていないのでマルチグリッド法とは呼称されない．
修正方程式解法に再帰的にマルチグリッド処理を適用することで初めてマルチグリッド法と呼称される．
今回の計算では，SOR 法を smoother としてグリッドを 3 枚使用し，残差が 1×10-5 以下になるまで計算を繰り返した．

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

31

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

32

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

細かい格子

100 μm

25 μm

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

50 μm

使用される仮想格子

400 μm

100 μm

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

33

2.4 中性粒子・スパッタ粒子計算手法
推進剤の中性粒子とスパッタ粒子に対する粒子解析は各々DSMC 法 (2-20)(2-21)(Direct
Simulation Monte Carlo method)を用いた．この DSMC 法は，多くの量の粒子の位置座標と速度成分を記憶し，粒子の移動と粒子間衝突の計算を繰り返すことによって流れ場の状態を算出する手法である．この方法は多くの乱数を用いるためここでも Mersenne
Twister(2-5)を使用した．
2.4.1 中性粒子の運動
ここで扱う粒子は，全て電荷を持たない粒子である．このため，これらの粒子は電場・磁場の影響を受けない．計算では 1 ステップ経過するごとに全ての粒子を時間ステップ𝛥𝑡だけ動かす．動かす前の位置を𝒓𝟎 ，速度を𝒄𝟎 とすると𝛥𝑡後には𝒓 = 𝒓𝟎 + 𝒄𝟎 𝛥𝑡(2.4.1)になる．
また，全ての粒子の運動を追跡するのは計算コストの上で非効率であるため，ひとつの粒子が複数の粒子を表す超粒子として扱う．
2.4.2 中性粒子の初期条件と境界条件
中性粒子の粒子数は任意の推進剤流量を設定し，それに対応した粒子数をステップごとに流入するようにし，徐々に中和器内部が粒子で満たされるようにした．流入する粒子の初期速度(gap)はマクスウェル分布から抽出した．気体定数を R とすると初期速度は x,y,z方向に対しそれぞれ
𝑐𝑥 = 𝐴𝑐𝑜𝑠𝜃

(2.4.2)𝑐𝑦 = 𝐴𝑠𝑖𝑛𝜃

(2.4.3)𝑐𝑧 = 𝐴′𝑠𝑖𝑛𝜙′(2.4.4)となる．ここで𝜃 = 2π𝑈，𝜙 ′ = 2𝜋𝑈で𝑈は0 < 𝑈 < 1の一様乱数で𝐴 = (−2𝑅𝑇0𝑙𝑛𝑈)1/2
𝐴′ = (−2𝑅𝑇0 𝑙𝑛𝑈)1/2(2.4.5)(2.4.6)である．
また，中性粒子解析における境界条件は中和器の全ての壁面で乱反射が行われるとしオリフィス口を出たものについては計算領域外となるので消去した．
2.4.3 中性粒子間の衝突
中性粒子解析においては中性粒子間の衝突については最大衝突数法を用いた．𝛥𝑡を34

平均自由時間以下にとると，分子の運動と衝突を分離してよいという，分離の原理が成り立つ．このため，全ての粒子を無衝突で𝛥𝑡走られた後，衝突して粒子は速度を衝突後のものに書き換えればよい．各セル内で起こる衝突は他のセルでの衝突に無関係に扱う．
したがってどの分子も，衝突相手は自分の属するセルの粒子から選ばれる．ここで粒子は直径 d の剛体球であるとする．
あるセルに N 個の超粒子が含まれているとする．これらの粒子の平均速度を𝐶̅ とすると，１つの粒子が単位時間に衝突する回数ν(s −1 )は𝐶̅
(2.4.7)
ν = = 𝑛𝑔̅ 𝜎𝑇𝜆となる．ここで，σT は衝突断面積であり，λ = 1/√2𝑛𝜎𝑇 を用いた．また，２粒子間の平均相対速度は𝑔̅ = √2𝐶̅ で与えられることを用いた．セル内の超粒子 i に着目すると式(2.4.2)の
𝑔̅ は
𝑁

1
𝑔̅ = ∑ 𝑔𝑖𝑗
𝑁(2.4.8)𝑗=1で近似できる．ただし，𝑔𝑖𝑗 = |𝒄𝒋 − 𝒄𝒊 |で𝑔𝑖𝑖 = 0である．しかしボルツマン方程式から，上式の1/Nを1/(N − 1)で置き換えてはいけないことが分かっている．𝑔̅ を式(2.4.2)に代入し，
νΔtがνΔt < 1のとき粒子 i の衝突確率𝑃𝑖 を表すことに注意すると
𝑁

𝑁

𝑃𝑖 = 𝑁 −1 𝑛Δ𝑡 ∑ 𝑔𝑖𝑗 𝜎𝑇 = ∑ 𝑃𝑖𝑗
𝑗=1(2.4.9)𝑗=1

ただし𝑃𝑖 = 𝑁 −1 𝑛Δ𝑡𝑔𝑖𝑗 𝜎𝑇 である．𝑃𝑖𝑗 は粒子 i がΔ𝑡後に粒子 j と衝突する確率であるので相対速度𝑔𝑖𝑗 が大きいほど衝突を起こしやすいことが分かる．今着目しているセル内でΔ𝑡
時間内に起こる衝突の回数𝑁𝑐𝑜𝑙 は 35
𝑁

𝑁

𝑁

1
1
𝑁𝑐𝑜𝑙 = ∑ 𝑃𝑖 =𝑛Δ𝑡 ∑ ∑ 𝑔𝑖𝑗
2
2𝑁
𝑖=1

(2.4.10)𝑖=1 𝑗=1となる．全ての𝑃𝑖 が１のとき，衝突する粒子の数は N であり，衝突回数は N/2 である．このため，因子 1/2 が必要である．式(2.4.5)の二乗和を求めるためには(𝑁 2 − 𝑁)/2個の𝑔𝑖𝑗 を計算する必要がある．この計算量がほぼ N の二乗に比例すること計算効率を著しく低下させる．そこで𝑔𝑖𝑗 の最大値𝑔𝑚𝑎𝑥 を適当に見積もり式(2.4.5)の𝑔𝑖𝑗 を全て𝑔𝑚𝑎𝑥 で置き換え
た衝突数𝑁𝑚𝑎𝑥 を用いる．
1
𝑁𝑚𝑎𝑥 = (𝑁 − 1)𝑛𝑔𝑚𝑎𝑥 𝜎𝑇 Δ𝑡2(2.4.11)となる．数密度 n は𝑛=𝑊𝑁
Δ𝑥Δ𝑦Δ𝑧(2.4.12)で表される．W は重みである．λは𝜎𝑇 に反比例するので𝑁𝑚𝑎𝑥 は𝜎𝑇 に依存しなくなる．
35

𝑁𝑚𝑎𝑥 は𝑁𝑐𝑜𝑙 を過大に見積もったものであり補正する必要がある．まず，着目しているセルに対して𝑁𝑚𝑎𝑥 回の衝突を計算する．その回ごとに二つの粒子を選びその粒子対(i,j)に対する𝑔𝑖𝑗 /𝑔𝑚𝑎𝑥 の確率で本当の衝突を起こす．
対(i,j)が衝突するとき衝突方向は乱数により選ばれるようにした．衝突後の速度(𝑐𝑖 ′, 𝑐𝑗 ′)は
1
𝑐𝑖 ′ = 𝒘𝑖𝑗 − 𝑔𝑖𝑗 𝑹
2(2.4.13)1(2.4.14)𝑐𝑗 ′ = 𝒘𝑖𝑗 + 𝑔𝑖𝑗 𝑹
2で与えられる．ただし，𝒘𝑖𝑗 は重心の速度，R はランダムな方向を持つ単位ベクトルである．
2.4.4 スパッタ粒子の運動
スパッタ粒子の運動は中性粒子と同様に扱い式(2.4.1)を用いた．
2.4.5 スパッタ粒子の初期条件と境界条件
スパッタ粒子のステップ当たりの発生数はスパッタ率を計算することで求めた．スパッタ率の計算は，実験データ及びシミュレーションから得たデータをもとに，任意の入射イオ
ンと固体に対して成立する山村らによる以下の半経験式(2-22)(2-23)が提案されている．
Q𝛼 ∗
𝑆𝑛
𝐸𝑡ℎ
√
(
)
[1
]
𝑌 𝐸 = 0.042−
𝑈𝑠 1 + Γ𝑘𝑒 𝜀 0.3
𝐸s(2.4.15)

ここで，Q，s はフィッティングパラメータで𝑈𝑠 はターゲット原子の表面結合エネルギーである．その他の値は以下のように定義される．
𝛼 ∗ = 0.249 (M2 0.56
M2 1.5)+ 0.0035 ( )𝑀1𝑀1M1 ≤ 𝑀2

M2 −0.15
M2 1.5(2.4.16))+ 0.165 ( )M1 ≥ 𝑀2
𝑀1
𝑀1
ここで，M1 , 𝑀2 はそれぞれイオンとターゲット原子の質量である．また，核的阻止能𝑆𝑛 は= 0.088 (𝑆𝑛 = 84.78

𝑍1 𝑍2

𝑀1
𝑠𝑛TF
1/2
2/3
2/3
𝑀
+
𝑀
1
2
(𝑍1 + 𝑍2 )(2.4.17)となる．ここで，𝑍1 , 𝑍2 はそれぞれイオンとターゲット原子の原子番号である．また，換算核的阻止能は𝑠𝑛TF

36

𝑠𝑛TF = 3.441√𝜀 ln(ε + 2.718)
1 + 6.355√𝜀 + 𝜀(6.882√𝜀 − 1.708)(2.4.18)

ここで換算エネルギー𝜀は𝜀=0.3255

𝑀2
𝐸
1/2
2/3
2/3
𝑀
+
𝑀
1
2
𝑍1 𝑍2(𝑍1 + 𝑍2 )(2.4.19)W
1 + (𝑀1 /7)3(2.4.20)となる．また，ΓはΓ=となる．ここで，W はターゲット原子に依存するフィッティングパラメータである．Lindhardの電子的阻止能𝑘𝑒 は
ke = 0.079

3/2 3/2(𝑀1 + 𝑀2 )3/2
3/2

3/2

𝑀1 𝑀2

𝑍1 𝑍2

(𝑍13/2 +𝑍23/2 )(2.4.21)

3/4となる．最後に閾値のエネルギーE𝑡ℎ は𝐸𝑡ℎ ==1 + 5.7(𝑀1 /𝑀2 )𝛾

M1 ≤ 𝑀2

6.7
𝛾

M1 ≥ 𝑀2(2.4.22)

ここで，𝛾は弾性散乱におけるエネルギー移行因子で𝛾=

4𝑀1 𝑀2(𝑀1 + 𝑀2 )2(2.4.23)である．
今回ターゲットとして計算に用いたアルミニウム，鉄，銅，モリブデンに対するパラメータ
Us , Q, W, sの値は以下のようになる．
表 2.1 スパッタ率に関するパラメータtarget

Z

Us (

)

Q

W

s

Al
Fe

13
26

3.36
4.28

1.0
0.75

2.17
1.20

2.5
2.5

Cu
Mo

29
42

3.49
6.82

1.0
0.85

0.73
2.39

2.5
2.8

また発生した粒子の初期速度は初期エネルギー(2-24)を Thompson による式により求めることで得た．エネルギー𝐸1 を持った一次イオンによるスパッタリングの場合，粒子の放出
エネルギー𝐸の分布関数𝑓(𝐸)は37

𝑓 (𝐸) ∝ 𝐸1 − √(𝐸𝑏 + 𝐸)/Λ𝐸1(𝐸 + 𝐸𝑏 )3(2.4.24)

Λ = 4𝑚𝑖 𝑚𝑡 /(𝑚𝑖 + 𝑚𝑡 )2(2.4.25)で与えられる．ここで，𝐸𝑏 はターゲット材料の原子結合エネルギー，𝑚𝑖 と𝑚𝑡 はそれぞれイオンとターゲットの材料原子の質量である．Thompson の理論によれば，放出角度は余弦則に従う．すなわち，放出角度は法線から角度𝜃を取れば分布関数はcos𝜃に比例するためこれを用いた．
境界条件は C/C について再付着率 0.78 という実験結果((2-25)があるため，Fe に対しても同じであると仮定してこれを利用した．つまり壁面へと到達したスパッタ粒子の 78%が壁面へと付着したとし計算から消去し，22%は乱反射するとした．
2.4.6 スパッタ粒子・中性粒子間の衝突
中和器の放電室中には推進剤の中性粒子最も多く存在しスパッタ粒子の解析をする上でその存在は無視することができない．スパッタ粒子と中性粒子の間の衝突はまず，自由行程を計算しスパッタ粒子の変位が自由行程を越えた時点で衝突を行うようにした．スパ
ッタ粒子の平均自由行程𝜆は Tuner の論文(2-26)より次のように与えられる．
𝜆=√𝜋𝑚𝑔 𝑣𝑝2
2𝜎𝑛𝑔 𝑘𝑇𝑔

1
𝑚𝑔
ψ(𝑣𝑝 √2𝑘𝑇 )(2.4.26)𝑔𝑥

ψ(𝑥 ) = 𝑥 xp(−𝑥 2 ) + (2𝑥 2 + 1) ∫ xp(−𝑦 2 ) 𝑑𝑦

(2.4.27)0

ここで𝑚𝑔 , 𝑛𝑔 , 𝑇𝑔 はそれぞれ中性粒子原子の質量，密度，温度で𝑣𝑝 はスパッタ粒子の速度，𝜎は衝突断面積である．この平均自由行程を用い自由行程𝑙はλを平均値として分布関数𝑓(𝑙) ∝ xp(−𝑙/𝜆)となる．
衝突後の方向は乱数により与え，衝突後の速度は衝突の判定が行われるごとに中性粒子の速度をマクスウェル分布により抽出し，その粒子に対する相対速度として，(2.4.13),(2.4.14)と同様に考えた．

38

2.5 全体の流れ
本研究ではマイクロ波によるプラズマの生成からプラズマの運動，更にそれに伴うスパッタリングを評価しそれによりマイクロ波の性能が低下するという繰り返しを全て数値計算
により評価することを目標としている．これを行うにあたり，まず実験値等を初期条件としス
パッタリング解析を単独で行った．図 2.5.1 にこの計算の流れを示す．現段階ではスパッタリング解析からマイクロ波解析へのフィードバックをするようなループ計算は行っておらずマイクロ波解析からスパッタリング解析までで評価を行った．

図 2.5.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)
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)
39

(2-11) V. Vahedi and M. Surendra, Computer Physics Communications, 87 (1995), 179.
(2-12) 電気学会・マイクロ波プラズマ調査専門委員会：“マイクロ波プラズマの技術” (オーム社, 2003)(2-13) K. Koura, Phys. Fluids, 29(1986), 11.
(2-14) 橋本修, 安部琢美：“FDTD 時間領域差分法入門” (森北出版, 1996)(2-15)(2-16)(2-17)(2-18)

宇野亨：“FDTD 法による電磁界およびアンテナ解析” (コロナ社, 1998)
K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.
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)
廣池匠哉, “マルチグリッド法を用いたマイクロ波放電式中和器のプラズマ損失

機構解明”，平成 24 年度九州大学修士論文．
(2-19) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical Physics, 55(1971), 4100.
(2-20) 保原充，大宮司久明：”数値流体力学:基礎と応用”(東京大学出版会，1992)(2-21) 南部健一，渡辺安雄，米村茂，：”中間流・分子流のシミュレーション法-大面積微小ギャップの排気への応用-“，真空, 43(2000)(2-22) Yamamura, Y. and Tawara, H.,”Energy dependence of Ion-Induced Sputtering
yields from Monotonic Solids at Normal Incidence”, NIFS-DATA-23 (1995)(2-23) 百武徹，西田迪雄，剣持貴弘，村本哲也，"イオンエンジン加速グリッドにおけるスパッタリング解析に向けて”，宇宙航空研究開発機構特別資料，JAXA-SP-06019，(2007)(2-24) M. W. Thompson, Philos. Mag. 18, (1968)(2-25) Marker C, Clemons L, Banks B, Miller S, Snyder A, Hung C, Karniotis C, Waters
D, “Transport of Sputtered Carbon During Ground-Based Life Testing of Ion
Thrusters”, NASA TM-2005-213798(2-26) G. M. Tuner, I.S. Falconer, B. W. James, and D. R. McKenzie, J. Appl. Phys.
B65B, (1989)40

第 3 章 結果と考察
3.1 100 時間耐久試験との比較
3.1.1 計算体系
本計算を行うに伴い，コードを評価するためにμ10 用と同出力の引き出し電子電流
200 mA 級中和器で 100 時間の耐久試験を行った．この中和器の断面図を図 3.1 に示す．
放電室はより比較しやすくするようにスパッタリングの閾値が比較的低くスパッタリングを起こしやすい金属である銅を用いた．放電室の長さは 12 mm，直径は 18 mm でこれを挟み込むように軟鉄製のヨークを配置する．オリフィス径は 5 mm で長さは 8 mm である．実験ではアンテナは L 字型を用いているが計算では簡単のため棒状にした．また，アンテ
ナの周りには絶縁体を取り付けてありアンテナを絶縁している．さらに実験において空間電圧が高くなりスパッタリングをより起こしやすくなるように引き出し電圧を 50 V と高い値を設定した．また，計算は図の青色の部分に対して行った．

41

図 3.1 計算体系 (200 mA 級中和器)

3.1.2 中性粒子密度分布
スパッタ粒子と中性粒子間の衝突を評価するために DSMC 法を用いて中性粒子密度分布を計算した．推進剤はキセノンガスで流量は 0.049 mg/s として流入粒子数と流出粒子数が同程度となる 35 msec まで計算を行った．この時の中性粒子密度分布を図 3.2 に示す．推進剤の流入口付近は比較的高い5.0 × 1020 𝑚−3 程度であるが多くの領域では
2 × 1020 𝑚−3 程度であることが分かる．密度を中和器内部全体で平均を取った値は2.39 × 1020 𝑚−3 となった． コ ン ダ ク タ ン ス を 用 い た 計 算 に よ り 予 測 さ れ る 密 度 は
2.00 × 1020 𝑚−3 となり数値解析の結果は妥当であると考えられる．

(a) z=0.0035 m(b) x=0 m
図 3.2 中性粒子密度分布

3.1.3 スパッタリングによる質量変化
図 3.2 の中性粒子密度分布を用いスパッタリング解析を行いこれにより各部品の 1 時間での質量変化を計算した．この際，各壁面へのイオン電流値は大道による作動実験 (3-1)
で得られた表 3.1 のような値を用いた．また，イオンを壁面へ加速する空間電位は小野寺らによる研究(1-8)において 0.098 mg/s で 23 V とされている．今回流量が半分に絞られており，その場合空間電位は上昇するため計算においてプラス 5 V となる 28 V として計算を行った．それぞれの材料のスパッタリング閾値は表 3.2 に示すようになる．今回扱う粒子
は二価イオンであるため 28 V 電位で加速されるとそのエネルギーは二倍の 56 eV となる．
このため，今回の計算において全ての材料でスパッタリングを起こす．二価イオン存在比は前述の小野寺らによる研究から引き出し電圧 50 V で 0.02 とされているためこれを用い42た．図 3.3 に各部品の質量変化の実験値と計算結果を比較したものを示す．絶縁体，放電室，オリフィスに対しては何れも10−3 mg オーダーで実験と計算はよく一致している．
しかしながら，ヨークは実験のほうが質量が大きく，アンテナは計算では増加している質量が実験では減尐していることが分かる．この二つの部品はいずれも10−2 mg オーダーの比較的大きなずれが見られた．
この実験で実際に用いたバックヨークと絶縁体に対して蛍光 X 線分析により付着物の分析を行った．この結果のうち，数値計算に用いた鉄，銅，アルミニウム，モリブデンのみで割合を出したものと計算結果を図 3.4，図 3.5 に示す．図を見て分かるようにバックヨー
クにおいて実験ではアルミニウムが多く検出されていることが分かる．また，絶縁体については実験ではモリブデンが多くついたことが分かる．
以上のことから，質量を比較したときにヨーク部で差が出た原因は計算においてアルミニウムの付着が尐なくなっていた事による予測される．また，実験で絶縁体に多くのモリブデンが付着していたことはアンテナの質量が減っていたことに対応していると考えられる．すなわち，オリフィスとアンテナ付近で多くスパッタリングが起きており，このためこれらの場所に大きなエネルギーを持ったイオンが存在すると考えられる．
今回計算と実験に差が出たのはイオンのエネルギー分布が得られなかったためと考えられ今後，数値解析によりイオンのエネルギー分布を得ることでより正確な計算が行うことができると考えられる．
表 3.1 実験によるイオン電流値
バックヨーク
イオン電流密度

83.3

側壁

/mm2

154.2

フロントヨーク
/mm2

118.8

/mm2

オリフィス
77.1

表 3.2 各壁面材料のスパッタリング閾値

閾値

Fe

Cu

Al

Mo

34.24 eV26.60 eV40.15 eV46.82 eV43/mm2

図 3.3 各部品の質量変化(a) 計算(b) 実験
図 3.4 バックヨークの付着物の割合

44

(c) 計算(d) 実験
図 3.5 絶縁体の付着物の割合

3.2 マイクロ波・プラズマ解析とのカップリング
これまでのスパッタリング解析は純粋に数値解析のみで行ったものではなく壁面へ入射するイオンの電流やエネルギーは実験値を用いた．しかしながら，イオンのエネルギーは本来一定ではなく分布を持つことが予想される．また，本来目標としているように長時間行う実験を避けるためには数値解析のみで評価できる必要がある．このため，イオンの挙動をプラズマ解析により計算しこれをスパッタ粒子の発生に利用した．これにより実際の中和器を用いることなく数値計算のみによって，中和器内部のスパッタリングを予測することができる．
3.2.1 計算体系
プラズマ解析を用いたスパッタリング解析には体系として本研究室で試作した引き出し
電子電流 12 mA 級中和器に対し各解析を行った．この中和器の断面図を図 3.6 に示す．
放電室には銅よりスパッタリングを起こしにくいモリブデンを用いた．放電室の長さは 32mm，直径は 21 mm でこれを挟み込むように軟鉄製のヨークを配置する．また，放電室の周囲とフロントヨークの内側に径が 7 mm であるオリフィスに沿うように Sm-Co 磁石を配置している．前出の体系と同様アンテナの周りには絶縁体が取り付けてある．計算は図の水色の部分に対して行った．

45

図 3.6 計算体系 (12 mA 級中和器)

3.2.2 中性粒子密度分布
12 mA 用の中和器に対しても，中性粒子密度分布を計算した．推進剤はキセノンガス
で流量は 0.049 mg/s として流入粒子数と流出粒子数が同程度となる 60 msec まで計算を行った時の中性粒子密度分布を図 3.7 に示す．また密度を中和器内部全体で平均を取った値は2.02 × 1020 𝑚−3 となった．

46

(a) z=0.0035 m(b) x=0 m
図 3.7 中性粒子密度

3.2.3 壁面に入射するイオンのエネルギー分布
図 3.8 は廣池によって行われたプラズマ解析の結果(3-3)である．どの部位に対しても得られたエネルギーは 10~20 eV となった．すなわち，二価イオンの場合は 20~40 eV となると考えられる．また，0.1 nsec の間に各部位に流入するイオン数を図 3.9 に示す．今回の入射する二価イオンの量としてスパッタリング解析に用いた値はこのイオン数に二価イオン
存在比を掛けたものである．

図 3.8 壁面に入射するイオンのエネルギー分布 (廣池匠哉, “マルチグリッド法を用いたマイクロ波放電式中和器のプラズマ損失機構解明”，平成 24 年度九州大学
修士論文)

47

図 3.9 各部位に流入するイオン数

3.2.4 絶縁体上へのスパッタ粒子の堆積
絶縁体に堆積したスパッタ粒子による電気抵抗の変化で評価を行った．まず，寿命末期の抵抗を知るために，200 mA 級中和器で壁面を銅からモリブデンに取り替えたものに対してμ10 において性能が低下した 15,000 時間時点での抵抗を計算すると 50 kΩとなった．これに対し 12 mA 級中和器で絶縁部抵抗の時間変化を求めたものを図 3.10 に示す．図を見ると分かるように 100,000 時間を越えても 100 kΩ程度と比較的高い抵抗値を示した．このため，12 mA 級中和器は十分長時間利用が可能であるといえる．またこの様な差が生じた理由はそれぞれのマイクロ波投入電力の差によると考えられる．200mA 級
中和器では電力が 8 W であったのに対し，小型中和器では 2 W であるため放電室内のイオンのエネルギーに違いが出たと考えられる．

48

図 3.10 絶縁部抵抗の時間変化

参考文献
(3-1)
(3-2)

大道渉，”マイクロ波放電式中和器のイオン電流測定”，STEP-2011-033
Noriyoshi Onodera,et al., “Electron Emission Mechanism of Microwave

Discharge Neutralizer,” IEPC-99-160,2002(3-3)
廣池匠哉, “マルチグリッド法を用いたマイクロ波放電式中和器のプラズマ損失
機構解明”，平成 24 年度九州大学修士論文．

49

第 4 章 結論
本研究ではマイクロ波放電式中和器の損耗メカニズムを解明するため数値解析コードの開発を行った．また，それを用いて実際に小型イオンエンジン用の中和器の寿命評価を行った．得られた結論を以下に示す．
1) DSMC 法を用いたスパッタリング解析コードを作成し，イオンのエネルギーを一定として計算を行った．各部品の質量変化を比較することでコードが妥当であることが確認できた．しかし，付着物の比率を見ると絶縁体へのモリブデンの付着とヨークへのアル
ミニウムの付着が低く見積もられていることが分かった．この原因はエネルギーに分布を持たせていないことに起因すると考えられ，プラズマ解析によるイオンのエネルギー分布を用いることが必要であると考えられる．
2) マイクロ波解析，プラズマ解析の結果を利用し分布を持った形でイオンのエネルギーを与えて数値解析による小型中和器の寿命の予測を行った．「はやぶさ」に用いられたのと同じ体系のものと比較することで 5 万時間以上の作動時間が確認された．

50