卒

題

業

論

文

目

数値解析によるマイクロ波中和器の寿命低下に関する研究

氏

名

江川

雄亮

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

指 導 教 員

中島

秀紀

教授

山本

直嗣

准教授

森田

太智

助教

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

平成２７年２月

目次
第1章

序論

1.1 背景1

1.2 イオンエンジン

1

1.2.1 イオンエンジンの作動原理

11.2.2 イオンエンジンの種類

2

1.2.3 電子サイクロトロン共鳴

41.2.4 磁気ミラー効果

5

1.2.5 マイクロ波放電式中和器

71.3 マイクロ波放電式中和器の性能低下

10

1.3.1 はやぶさにおける中和器性能

101.3.2 スパッリング

10

1.3.3 中和器の電圧上昇

131.4 研究目的

16

第2章

数値解析

2.1 プラズマ粒子計算手法

18

2.1.1 基礎式

18

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

19

2.1.3 PIC 法

20

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

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

222.1.6 平均自由行程と透過度

222.1.7 衝突確率・衝突判定

23

2.1.8 衝突計算

232.1.9 エネルギーと散乱角242.1.10 Null-collision 法

27

2.2 電磁波解析手法

29

2.2.1 基礎式

29

2.2.2 差分式

29

2.2.3 アンテナ解析

31

2.2.4 安定条件

342.2.5 PEC 境界条件

34

2.2.6 Mur の吸収境界条件

342.2.7 ポインティングベクトル

35

2.3 計算条件

36

2.3.1 計算体系

37

2.3.2 計算格子

372.3.3 タイムステップ

37

2.3.4 粒子数

372.3.5 推進剤

38

2.3.6 磁場強度

382.3.7 絶縁体・抵抗被膜

38

第3章

結果と考察

3.1 絶縁体の挿入

44

3.2 抵抗被膜の挿入

47

第4章
謝辞

結論

第１章
1.1

序論

背景

2010 年,マイクロ波放電式イオンエンジンμ10 を搭載した小惑星探査機「はやぶさ」
は,地球重力圏外にある天体の固体表面に着陸してのサンプルリターンに世界で初めて成功した. また、その後継機である「はやぶさ２」が 2014 年 12 月 3 日に打ち上げられたことはまだ記憶に新しいだろう.これによりマイクロ波放電式イオンエンジンは宇宙実証され,累積稼働時間は約 40000 時間を超え，その有用性の高さが実証された.マイク
ロ波放電式イオンエンジンとは電気推進と言われるものの一種で,電気エネルギーを用いて推力を得る推進機である.従来の化学反応を利用した化学推進と違って静電力や電磁力を利用できる．このため比推力（推進剤流量あたりの推力）を一非常に大きくすることが可能であり，大幅な推進剤の低減が望める．現在使用されている代表的な電気推進機はアークジェットスラスタやホールスラスタ，そしてイオンエンジンである．イオンエンジンは他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制御
や惑星探査に適している．惑星探査においては，イオンエンジンは比推力が高い一方で，推力は化学推進に比べて著しく小さい.そのため,十分な速度が得られるまで多くの時間を要する. 「はやぶさ」も帰還までに約 7 年間かかり, 長時間運用が必須となったた
め，長寿命化が求められている．

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

1

図 1.2.1

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

1.2.2 イオンエンジンの種類
プラズマ生成の方式によってイオンエンジンの分類を行う事が多い。これまでに研究されてきているイオンエンジンは大きく以下の 5 種類に分けることが出来る。
1) 直流放電型（カウフマン型）ex) UK-10（ESA）
2) 直流放電型（カスプ磁場型）ex) NSTAR（NASA），XIPS25（Boeing)3) 高周波誘導加熱型 ex) RIT10（ESA）
4) マイクロ波放電型 ex) μ10（JAXA），HiPEP（NASA）
5) 接触電離型1)および 2)は、放電プラズマの生成に熱陰極からの電子による電離衝突を用いるタイ
プで、1)では放電室外部磁石による軸方向発散性の磁場が、2)では放電室内部磁石によるカスプ状磁場が、プラズマの閉じこめと電離能率向上のために形成されている。最も早く宇宙作動に成功したのはカウフマン型であり、現在は三菱電機や英国 QinetiQ 社で製造されている。リングカスプ型のイオンエンジンは世界で最も多く搭載されており、米国の深宇宙探査機 DeepSpace1、Dawn の为推進機である NSTAR や、L3 社が商用衛星に使用している XIPS シリーズがある。NSTAR エンジンはイオンエンジン単体での最高作
動時間（16000 時間）の記録を持つ。このような圧倒的な宇宙実績がある直流放電型イ
オンエンジンであるが、放電室内で発生したイオンが高い放電電圧によって陰極をスパ
2

ッタリング浸食し、発生した金属粉が電極短絡を誘発するなどの故障が発生する。また
電子発生源となる陰極も、単純なフィラメントからホローカソードなどの耐久性に優れたものに変わったが、温度サイクルなどによるホローカソードのヒーターの断線や、ゴミの付着による放電電極の短絡などの場合、システムの全損故障へと発展する危険性がある。
3)は RF（Radio Frequency）型(1-1)とも呼ばれ、kHz - MHz 帯の高周波を放電室外部から加えることで電子に電離衝突に必要なエネルギーを与え、プラズマを生成する。
4)は ECR（Electron Cyclotron Resonance)型とも呼ばれ、周波数が GHz 程度のマイクロ波により静磁場中の電子を共鳴的に加速・加熱することで電子による電離衝突が為される。3)や 4)は「無電極プラズマ生成」を実現することで、大きなアドバンテージを持つ。無電極放電によって生成されたプラズマは放電質壁面との電位差が小さく、スパッタリングが発生しない。また１次電子の発生源も必要としない。
5)の接触電離型は推進剤よりも仕事関数の大きい物質に推進剤を接触させることでイオン生成を行う機構である。
本研究では 4)のマイクロ波放電型イオンエンジン(1-2)(1-3)を使い実験を行う。マイク
ロ波放電とは、マイクロ波帯域の交流電場によって電子を加速し、中性粒子と衝突電離により電子の数が増加し気体がプラズマ化され放電が維持されることである。マイクロ波放電型イオンエンジンのシステム全体図を図 1.2.2 に示す。プラズマ源を正に、加速電極と中和器を負にバイアスするための３つの個別の電源が接続されている。プラズマ
源と中和器の電源を１つの電源で置換されることが多いが、「はやぶさ」のイオンエン
ジンシステムでは中和器のヘルスモニタや宇宙機電位の固定のために個別の電源を置き、それぞれの電源のリターンラインを宇宙機構体に接続している。

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

3

1.2.3

電サイクロトロン共鳴

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

mdv
 ev  B
dt

(1.2.1)

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

rL 

mv
eB

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

 ce 

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

電子

電子
E

B
B

4

電子
図 1.2.3 ECR の原理

E

ロン共鳴現

1.2.4 磁気ミラー効果
ECR によって電子を加速することにおいて、電子を外に逃がさずに効率よく加速させることができる理由に磁気ミラーによる電子の閉じ込めがある。磁場に垂直な方向の荷電粒子の運動は, 磁力線の回りをまわるような軌道を描くことは先で述べた。ここでは、プラズマ中に図 1.2.4 に示すような磁力線が存在する状況を考える。この図の点 O に旋回の中心を持ち、速度 𝑣0 （磁場に垂直な成分, 磁場に平行な成分）を持つ粒子があるとする。荷電粒子の運動エネルギーの保存より
𝑣 2 = 𝑣∦2 + 𝑣⊥2
2
2= 𝑣∦,0+ 𝑣⊥,0= 𝑣02(1.2.4)1となる。また、磁場が緩やかに変化する場合磁気モーメント 𝜇 = (2) 𝑚𝑣⊥2 /𝐵 が一定に保たれることを考慮すると、磁場強度が増加するに従い、速度の平行成分が減小していき、ある磁場強度で荷電粒子は磁場の弱い方向へ反射される。このような磁場による荷電粒子の反射を磁気ミラーと呼ぶ。𝜇 の一定性より、
122( )m𝑣⊥0

𝐵0

=

1
2

( )m𝑣02𝐵(1.2.5)である。また、𝑣⊥ = 𝑣 sin 𝜃 = 𝑣0 sin 𝜃05(1.2.6)

であることを考慮し、式(1.2.6)を式(1.2.5)に代入すると、
sin2 𝜃0sin2 𝜃
=𝐵0
𝐵(1.2.7)となる。よって反射点での磁場強度は、𝑣∦ = 0 と θ = 90° であるので、𝐵𝑟𝑒𝑓 =

𝐵0
sin2 𝜃0となる。

図 1.2.4 磁場による電子の回転

6(1.2.8)図 1.2.5 磁気ミラーによる粒子の閉じ込め
1.2.5 マイクロ波放電式中和器
マイクロ放電式イオンエンジンには，そのイオンビームに同量の電子を供給し電気的に中和する中和器という装置が搭載される．これにより，エンジンや宇宙機本体が負に帯電し噴出したイオンを引き戻すことなく，運転を継続することが出来る．中和器が故障し衛星からイオンのみが放出されることになるとイオンビームの抽出が不可能となる．
一般的な直流放電型イオンエンジンでは中和器としてホローカソードが使用されており，スラスタ本体と同様に電極の損耗や電源構成の複雑化が故障の要因となっている．
ヨーロッパで盛んに研究が行われている高周波放電型イオンスラスタ(1-4)でも中和器にはホローカソードが使用されているが，この問題は解決されていないままである．
このマイクロ波放電は直流放電と比べた場合，以下のようなメリットがある．
1. ヒーターを用いないため，予備加熱が不要となり素早い起動が出来る．
2. マイクロ波放電式イオンエンジン駆動用のマイクロ波電源が 1 台あれば中和器のプラズマも同時に生成することが可能なため，システム構成の簡素化や信頼性の向上，また衛星の重量軽減が図れる．
先に述べたようにマイクロ波放電式中和器ははやぶさに搭載され十分な中和性能と寿命が得られたという実績がある．しかし，小型衛星に搭載するようなサイズの中和器は未だ実用化されておらず，その開発が急務となっている．
本研究室では小型衛星への搭載を想定した小型マイクロ波放電型イオンエンジンを開発してきた．エンジン本体について優れた性能が得られた．一方中和器においては要求されるサイズ，消費電力でも中和性能の目標値をクリアしておらず，現在改善に向けて力を入れている．
図 1.2.6 は本研究において開発中のマイクロ波放電型小型中和器の外観写真である．
放電室は内径 21 mm の円筒形状であり，オリフィス開口径 10 mm，放電室長 32 mm であ7

る．放電室外部と出口部分に Sm-Co 磁石を配置し構成された磁気回路によって ECR に必
要な磁場を供給する．アンテナはモリブデン製である．

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

図 1.2.7マイクロ波放電式小型中和器運転時
81.3 マイクロ波放電式中和器の性能低下
1.3.1 はやぶさにおける中和器性能
小惑星探査機「はやぶさ」に搭載されたマイクロ波放電型イオンエンジン「 µ10」はJAXA 宇宙科学研究所の電気推進工学部門にて研究開発され、そして世界で唯一の実用化された ECR 放電型イオンエンジンである。宇宙科学研究所のマイクロ波放電型イオンエ
ンジンの開発は１９８０年代から始められ、１９８９年に１号機（Y-I）が完成する。
その後改良が重ねられ、２００１年には「µ10」において、推進剤利用効率にして４倍、イオン生成コストにして 1/10 もの性能向上と２万時間の耐久性を達成した。「µ10」に搭載された中和器は、先にも述べたように、従来の中和器と異なり、プラズマ化の際に電極を必要としないマイクロ波放電を採用している。投入マイクロ波電力は 8 W、ガス
は Xe を使用し、その流量は 0.05 sccm(standard cubic centimeter per minute : 標準状態における 1 分あたりの体積流量)、中和するための電子電流が、実際に推力を得るために陽イオンを放出しているイオン源と同等の 135 mA で中和器電圧が 22 V 程度となる
1.3.2 スパッリング
イオンによるスパッタリングは金属ターゲットの表面にエネルギーを持ったイ.オンが衝突することで発生する．このイオンは入射イオンと呼ばれ，これと直接衝突したター
ゲット原子は一次反跳原子と呼ばれる．この一次反跳原子は通常元の格子位置から弾かれたら他のターゲット原子と衝突を起こす．この衝突の連鎖は衝突カスケードと呼ばれ二次，三次と反跳原子が他の原子を弾けるだけのエネルギーを持つ限り続く．最終定期にターゲット原子がスパッタされて外に飛び出すためには原子がターゲット表面付近
で外部に向かうような軌道で弾かれ，且つその原子のエネルギーが表面結合エネルギーを超える必要がある．この衝突カスケードの様子を模式的に示した者を図 1.3.1 に示す．
スパッタリングは通常スパッタ率(Y )で評価される．これは一つの入射イオン当たりい
くつのスパッタ粒子が発生するかを示すもので図 1.3.1 の場合は一つのイオンが入ってきたのに対し外側に出ている原子の数が三つなのでスパッタ率 Y＝3 となる．
スパッタ率はイオンの質量，エネルギー，入射角，ターゲット原子の質量，原子構造，
表面結合エネルギーに依存する．入射イオンの質量が増えるにつれスパッタ率は増加する．また，エネルギーに関していえばイオンのエネルギーがスパッタリング閾値以下であるとスパッタリングを起こさず，これを超えると増加する．しかしながら，0.1 eV あたり超えると入射イオンがターゲットの表面に対し深く入り込みすぎるため表面から粒子が出にくくなりスパッタ率は減尐する．次に角度について述べる．スパッタ率は入射角がターゲット表面の法線に対する角度が 60°になるまで角度と共に増加しその後は 90°で 0 になるまで減り続ける．
Sigmund はスパッタリングを図 1.3.2 に示すようにシングルノックオン型，線形カスケード型，スパイク型の 3 タイプに分類した．これは入射イオンのエネルギーによって9

大まかに分類できる．まず 1 eV 以下程度の低エネルギーで起こるシングルノックオン型スパッタリングでは衝突カスケードは発生しない．この場合，最終的に外部に弾き出されるスパッタ粒子は元の入射イオンにより弾き出されるため極尐ない衝突しか起こさない．また，線形カスケード型のスパッタリングは 1~300 eV 程度で起こり，その名の通り衝突カスケードを起こす．この場合のスパッタ粒子は反跳原子によって弾き出されるが反跳原子と他のターゲット原子との間で衝突が起こるのは希である．100 eV 以上
若しくは入射イオンのフラックスが大きい場合に起きるスパイク型のスパッタリング
も衝突カスケードを起こし，またイオンのエネルギーが高いため多くの衝突が起きる．
特にこの場合の特徴として既に他の原子により弾かれ動いている原子に対して更に衝突するといったことが起きる．イオン推進機のような比較的低エネルギーのイオンを扱うものでは为にシングルノックオン型のスパッタリングを起こす．

10

図 1.3.1(a)シングルノックオン型< 1 keV
図 1.3.2

衝突カスケード(b)衝突カスケード型(b)スパイク型

1 keV~300 keV>100 keV
Sigmund によるスパッタリングの分類

11

1.3.3中和器の電圧上昇(1-2)

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

12

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

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

図 1.3.4 中和器の性能低下のメカニズム13

図 1.3.5 中和器電圧の時間変化(1-5)（出展：大道 渉， 小泉 宏之， 国中 均, 細田 聡史， 西山 和孝, “マイクロ波放電式中和器の性能低下メカニズムとその実験的検証”，STEP-2010-058）

14

1.4 研究目的
本研究室では，これまで実験や数値解析を用いて中和器内部のスパッタ粒子の挙動を解析してきた．しかし，スパッタ粒子による中和器内部への影響の調査は行われていなかった．この評価を実験により繰り返し行うことは困難である．中和器の性能低下は
約 10,000 時間以上作動させ続けてから発生する．このため，実際に実験により様々な体系で中和器を作動させ，最適な形状を調査することは莫大な時間と費用を要することが問題となる．この様な場合，数値解析は非常に有効な手段である．
本研究の目的は，数値解析を用いることで時間と費用を削減し,スパッタ粒子による中和器内部への影響を調査し,中和器性能低下の要因を考えることで寿命評価へ貢献することである.
そのため,中和器内部のアンテナ周囲に絶縁体を挿入し,さらにスパッタ粒子を模擬した抵抗被膜を絶縁体表面に表現し,電界強度と電子温度の変化を解析した.

15

参考文献
(1-1) 大熊 拓児 : 誘導結合型プラズマ電子源のホールスラスタへの適用検討, 首都大学東京大学院, 修士論文, 2010(1-2) 細田聡史，國中均 : イオンエンジンによる小惑星探査機「はやぶさ」の帰還運用, J. Plasma Fusion Res. Vol.86, No.5 (2010)282-292(1-3) 荒川義博,國中均,他 : イオンエンジンによる動力航行 (宇宙工学シリーズ
8,2006)(1-4) N. Wallace, P. Jameson, IEPC, 2011, 327(1-5) 大道 渉， 小泉 宏之， 国中 均, 細田 聡史， 西山 和孝, “マイクロ波放電式中和器の性能低下メカニズムとその実験的検証”，STEP-2010-058(1-6) Noriyoshi Onodera,et al., “Electron Emission Mechanism of Microwave
Discharge Neutralizer,” IEPC-99-160,2002

16

第2章

数値解析

2.1 プラズマ粒子計算手法(2-1)
2.1.1 基礎式
プラズマの挙動を解析する際の基礎式として運動方程式を用いる．式を以下に示す．
mdv
 q E  v  B dt(2.1.1)dx
v
dt(2.1.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.1.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.1.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(2.1.5)と書ける．ここで， a ・  ・  ・  を次のように定義する．


qt n
Bx
2m

(2.1.6)



qt n
By
2m

(2.1.7)

17

 

qt n
Bz
2m

(2.1.8)  1  2   2   2(2.1.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.1.10)(2.1.11)

ここで， c  qt m ， m  n  1 2 である．さらに(2.5)式を次のようにして変形し，速度を求めることができる．
v n+1 2 = A -1 S 1  21
A 1              
1   2       1   2 

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

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

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

ne e 2
 pe me  0に対して，
18(2.1.16)t   pe1(2.1.17)とすれば，数値的な不安定性は生じない．時間積分の精度まで考えると，通常は

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

c 

eB
me

(2.1.19)に対して，

t  c1(2.1.20)である必要がある．廣池らによると尐なくとも 100 分の一にしないとエネルギーが保存されないことが明らかになったため，
Dt <<

wc

(2.1.21)100となるように設定した．
2.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)

となり，

19

(2.1.22)

E ( x, y, z ) = A1E (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)(2.1.23)+ A8 E (i, j +1, k +1)と表わされる．

図 2.1.1

PIC 法におけるセルの構成

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

V1
qN
V2

(2.1.24)となる．すべての粒子について和をとれば，格子点状での電荷密度が得られる．また，
電流密度は
j  qnv(2.1.25)より，
j (i , j , k ) =

V1
qNv
V2

2.1.5 荷電粒子と中性粒子の衝突
20(2.1.26)

荷電粒子は電磁場の影響を受けて，中性粒子と衝突する．プラズマは電子，イオン，
中性粒子からなるとする．電子とイオンの数密度は中性粒子の数密度より非常に小さい
ため，電子-電子，電子-イオン，イオン-イオンの各衝突は無視できるものとする．電子と中性粒子は弾性，励起，電離衝突をすると仮定する．なお，以下に用いる乱数 (2-5)
は(0,1)間に一様に分布する乱数である．
2.1.6 平均自由行程と透過度(2-5)
面  ，厚さ l ，その物質の数密度 N の板を考える．エネルギー E の粒子ビームが，面
 に垂直に当るものとする．衝突をせずに板を透過してくる粒子数は，入射粒子数より減尐する.その減尐度は  を衝突の断面積とすれば，dn  nNdl(2.1.27)と表わされ，これよりn  n0 exp  Nl (2.1.28)となる．ここで， n0 は面  （ l  0 ）でのビーム中の粒子である．それゆえ，（ l ,l  dl ）間で最初の衝突が起こる確率は，

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



1
N

(2.1.31)

これを平均自由行程という．よってモンテカルロ法によって，任意の点から出発し最初の衝突が起こるまでの距離 l を決めるには，乱数 r を用いて，r  Pl  l
 1  exp    (2.1.32)で表わされ，これを変形して，l  ln 1  r となる．ここで， r は一様乱数であるから 1  r を r とすれば，l  ln(r )21(2.1.33)(2.1.34)

と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を  とすれば，この粒子が１秒間に気体分子と衝突する回数  は，



v


 Nv

(2.1.35)となり，これを衝突周波数という．
2.1.7 衝突確率・衝突判定
すべての衝突の断面積は，荷電粒子の全エネルギーの関数として表わされる．
   E 

(2.1.36)

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

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

(2.1.37)(2.1.38)

ここで， ngas は中性粒子の数密度， v m は粒子 m の速度，  total は全衝突断面積である．
ある乱数 r1 （ 0  r1  1）を用いて， Pcollision,m  r1 のときに粒子 m はタイムステップ t の間に衝突が起こるとする．
2.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 である．このように断面積の比と乱数を対応させ，衝突の種類を決める．計算では中性粒子としてキセノンを用い，その断面積データ
22(2-7)(2-8)(2-9)

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

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

v1 ,v 2  とすると，
v1 1
v1  v 2  u2(2.1.39)v 2 

1
v1  v 2  u2(2.1.40)と表わすことができる．ここで， 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)．
23

v¢e =

m + M cos c
M sin c
ve +
h
m+M
m+M

(2.1.41)h の直交座標成分はhx  ver cos hy  
hz  

(2.1.42)vex vey cos   ve vez sin (2.1.43)vervex vez cos   ve vey sin (2.1.44)

ver

ここで，
ver 
ve 

2
2
vey
 vez

(2.1.45)

2
2
2
vex
 vey
 vez

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

Escattered,e  Eincident,e(2.1.47)

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

g( c ) =Eincident,e sin cc
2 1+ Eincident,e sin 2 ln 1+ Eincident,e2{}(0 £ c £ p )(2.1.48)

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

2 1- (1+ Eincident,e ) 3rEincident,e}(2.1.49)となる．この値を(2.1.41)式に用いて衝突後の速度を求める．このとき，任意の角は  2 r4(2.1.50)として与える．
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーを
Eexcitation とし，入射電子のエネルギーを Eincident,e ，励起後のエネルギーを Escattered,e ，励起後の24

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

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

(2.1.54)

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

E  Eincident,e  Escattered,e(2.1.55)と表わされる．(2.1.52)式から電離後の電子の速度は，v~  v e 1 EEincident,e(2.1.56)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．(2.1.41)式中v e に v~ を用いて，衝突後の速度 v e を求める．散乱角には(2.1.49)式を用いるが，(2.1.49)
式中の Eincident,e には(2.1.54)式から求められる Escattered,e を用いる．電離エネルギー Eionization は，
キセノンの場合 12.13 eV である．
電離後の生成電子の持つエネルギーは，Ecreated,e  E  Eionizationとなる．また，入射電子の速度 v e を用いて生成電子の速度は25(2.1.57)

v created,e 

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

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

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

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

 total   elastic   excitation   ionization   fake とし， elastic   excitation   ionization r2  1
 totalのとき
Null-Collision
Processとする．
Null-Collision Process が選ばれたときには，粒子の運動は何ら変化しないものとする．図 2.1.3 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度で規格化した衝突周波数を示す．本計算では Null-Collision 法で用いた 𝑣𝑐𝑜𝑙𝑙𝑖𝑠𝑖𝑜𝑛 /𝑛𝑔𝑎𝑠 の値は 2×10-12 m3/sec である．

26

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

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

  E   0
  H  0

ここで， E は電場， H は磁場，  0 は真空の透磁率， 0 は真空の誘電率，J は電流密度を表わす．
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
ïî
ïþ
n+

1

n-

1

H y 2 (i +1 2, j, k +1 2) = H y 2 (i +1 2, j, k +1 2) +Dt
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
ïî
ïþ
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 )ìï 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
ïî
ïþ

28

(2.2.3)(2.2.4)(2.2.5)

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 
H
i

1
2
,
j

1
2
,
k

H
 z
z

y


n



Hy

1
2

n

i  1 2 , j , k  1 2  H y

1
2

z

E yn 1 i , j  1 2 , k   E yn i , j  1 2 , k  (2.2.6)i  1 2 , j , k  1 2
 J x i  1 2 , j , k t
 0 i , j  1 2 , k 1n n  122i , j  1 2 , k  1 2
H
i
,
j

1
2
,
k

1
2

H
 x
x

z





H

1
n
2
z

i  1 2 , j  1 2 , k   H

1
n
2
z

i  1 2 , j  1 2 , k xE zn 1 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  122i  1 2 , j , k  1 2
H
i

1
2
,
j
,
k

1
2

H
 y
y

x





H

1
n
2
x

i , j  1 2 , k  1 2  H

1
n
2
x

i , j  1 2 , k  1 2y(2.2.8) J z i , j , k  1 2



n

1

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

29

図 2.2.1

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

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

𝑉
𝛥𝑥
𝑟𝑙𝑛 ( 𝑎 )(2.2.9)𝐸𝑥 =𝑉𝛥𝑦𝑟𝑙𝑛 ( 𝑎 )(2.2.10)で与えられる. ただし, V は中心導体と外部導体間の電圧, 𝑟 は動径方向の距離,
𝑎 はアンテナの半径である.
30

また, 同軸ケーブルのインピーダンス 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 øここで，式中の(𝐼𝐹𝐸𝐷, 𝐽𝐹𝐸𝐷, 𝐾𝐹𝐸𝐷)は電場を励振する座標を与える．

図 2.2.2

同軸線路給電31(2.2.14)図 2.2.3 アンテナモデル

32

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

1

Dt £
v

1
1
1
+
+
2
2
(Dx ) (Dy ) (Dz )2(2.2.15)を満たさなければならない．この条件は極めて厳しく，わずかでも満足しなければ計算が不安定になる．
2.2.5 PEC(Perfect Electric Conductor )境界条件
境界面が金属のような完全導体である場合，境界に平行な電場成分を 0 とする PEC 境
界条件を用いる．放電室内壁は導体を仮定しているので，この境界条件を適用する．
2.2.6 Mur の吸収境界条件(2-17)
開放領域を取り扱う場合，計算領域を仮想的な境界で閉じなければならない．この境
界を吸収境界といい，この吸収境界が完全でなければ，反射波が解析領域に戻ってしま
い，誤差の原因となる．ここでは，定式化が容易で広く使用されている Mur の吸収境界条件を用いる．

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

(2.2.16)と表わされ，次の微分方程式を満足する．
E z 1 E z0
x v t(2.2.17)x = 0 の境界で完全に反射がないならば，電場は(2.2.16)式の形を保ったまま伝播するはずである．すなわち，境界でも(2.2.17)式が満足される．これを FDTD 法の形式に定式化するため，時間について差分化すると¶E z E zn - E zn -1
=
¶t
Dt
n-

=v

¶E z
¶x

1
2

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

33

n

1

n

1

Ern i, n  1 2  Ern1 i, n  1 2
E 2 i, n   Er 2 i, n  1
v r
t
z
となる．(2.2.18)式の Ez 1 2, Eznn 1(2.2.19)

1 2 も FDTD 法では直接求めることはできないため，前後の値で平均して求める．そして， E z についてまとめるとEzn ( 0, j, k +1 2) = Ezn-1 (1, j, k +1 2)+

vDt - Dx n
{Ez (1, j, k +1 2) - Ezn-1 ( 0, j, k +1 2)}
vDt + Dx

(2.2.20)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．
2.2.7 ポインティングベクトル(2-18)
空間を流れる電力を求めるためには，ポインティングベクトルを使用する．Maxwellの方程式を変形すると，以下の式となる．

  E  H   dS   J  Edv svol  E 2 H 2   2  2 dv (2.2.21)
t vol

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

34

2.3 計算条件
本研究ではマイクロ波によるプラズマの生成からプラズマの運動，更にそれに伴うス
パッタリングを評価しそれによりマイクロ波の性能が低下するという繰り返しを全て数値計算により評価することを目標としている．これを行うにあたり,笹川の結果をもとに,スパッタ粒子の導電率を算出し,スパッタ粒子の堆積によるマイクロ波の性能低下の評価を行った.

図 2.3.1 計算の全体の流れ

35

2.3.1 計算体系
図 2.3.2 に本研究において使用した計算体系を示す．円筒形の放電室をヨークで囲み，オリフィス開口部より電子を放出する．スペーサーによってオリフィスとの位置を調節している．磁石の個数は 14 個とし，放電室を囲んだ．.

図 2.3.2 計算体系
2.3.2 計算格子
従来は計算格子の長さを 0.5 mm とし，x 方向に 40 格子 y 方向に 40 格子 z 方向に 36
格子を作成していたが,今回はアンテナとヨークの間にセラミックと抵抗被膜を挿入するため,より細かい格子の長さを 0.25 mm とし, x 方向に 80 格子 y 方向に 80 格子 z 方向に 72 格子,合計 460800 格子を作成した.
2.3.3 タイムステップ
PIC 法の時間刻みは t pic  5.0  10

13

sec，FDTD 法の時間刻みは t FDTD  2.0  10

12

sec

とすれば，2.1.2 や 2.2.4 で述べた安定条件を満たす．FDTD コードが 10 回計算された
後，PIC コードを 1 回計算させた．
2.3.4 粒子数
1 ボクセル当たり数個しか粒子がない場合，ノイズが大きくなり望ましくない (2-12)．
通常，10 個から 100 個以上の粒子が存在しているとよいとされている．今回は 1 ボクセ
ルに平均 25 個程度の粒子が含まれるように粒子数を設定した．
36

上記の計算体系における体積 5.46×10-6 m3，実験による粒子密度は 10-17~10-18 m-3 と分かっているので，平均化して 5×10-17m-3 とした．超粒子数は 100 万個として，超粒子 1 つが何個代表するかは，体積と粒子密度より超粒子 1 つは約 27 万個を代表すると決定した．
2.3.5 推進剤
推進剤はキセノンとし，中性粒子密度は 2.1×1020 /m3，推進剤流量は 0.05 mg/s とした．
2.3.6 磁場強度
今回は電源としてはやぶさにも採用されている 2.45GHz を用いる．マイクロ波周波数と電子サイクロトロン周波数が一致するときの磁場強度は 0.085 mT である．その値を示す領域を ECR 層と呼ぶ．磁場計算には Advanced Science Laboratory, Inc の 3 次元
静磁場解析ソフトを使用した(2-28)．
2.3.7 絶縁体・抵抗被膜(2-14)
図 2.3.3 のように,アンテナとヨークの間にセラミックの絶縁体を挿入した.このとき絶縁体はアンテナから 1 セル分離れていて,厚さ 0.25 mm, 長さ 7.75 mm, 比誘電率 7.4としている. また絶縁体の表面に堆積するスパッタ粒子を,抵抗被膜を用いて表現した.
抵抗被膜は非常に薄い膜を表現でき,セルを小さくすることにはメモリに限界がある
ため,今回のように解析対象として極めて薄い薄膜を解析する場合に有効である.FDTD
法では,一般的に物体をセル単位で構成するため,抵抗被膜の厚さを有限とする必要がある.つまり,非常に薄い抵抗被膜を表現する場合, 図 2.3.3(a)および(b)のように最低でも 1 セルで構成し,モデル化することになり,セルサイズ分の厚みを有する抵抗被膜として表現することになる.しかし,下記の近似式(2.3.1)を用いると,図 2.3.3(c)のように 1 つのセルの 1 表面にのみ導電率を入力することで,非常に薄い抵抗被膜を表現することができる.この理由として,図 2.3.5(a)のように,真空と抵抗被膜の境界におけるセ𝑛+1/2ルの磁界𝐻𝑦𝑛+1/2の計算は,一般的な物体のモデル化では,磁界𝐻𝑦として考え,境界の𝑛導電率を真空との平均値とし,導電率の影響は電界𝐸𝑧1までの位置として考えるためで
𝑛+1/2ある.しかし,磁界𝐻𝑦𝑛は,真空の導電率を用いて計算された電界𝐸𝑧2と抵抗被膜の導電𝑛率を用いて計算された電界𝐸𝑧1の位置における抵抗被膜部分の導電率を平均値としない

場合,図 2.3.5(b)のように,真空と抵抗被膜との境界部分の磁界の値になると考えられる.このように,磁界(𝐻𝑦 )の位置は,電界(𝐸𝑧1 )と電界(𝐸𝑧2 )の間,すなわち,セルの中心の位置に存在し,セル単位で考えると,抵抗被膜の導電率の影響は半セル分まで影響していることになる.また,1 つのセルの 1 表面のみに導電率を入力して抵抗被膜を表現した
37

場合,左右で半セル分の影響を受けているため,両側で 1 セル分の影響を受けていると考えることができる.そのため,抵抗被膜を構成するセルに対して入力する導電率の値は,抵抗被膜近似式(2.3.1)から導出し,その際に抵抗被膜の厚さ t は,仮想的に 1 セル分のサイズを用いて計算する.
以上の理由から,抵抗被膜は絶縁体の外側に厚さを持たず,セルの境界面上に挿入した.このとき長さは 7.75 mm とし,導電率は抵抗被膜近似式及び笹川の結果をもとに 2.58
×10-2,1.29×10-3S/m を用い,FDTD 法の真空の伝導率での補正因子の式に組み込んだ.このとき配列を用いて,アンテナの中心軸からの距離で絶縁体,抵抗被膜,ヨーク,プラズ
マを領域判定し,アンテナの中心軸から 5.0×10-4m の時,算出した導電率を代入するよう設定した.
抵抗被膜近似式を以下に示す.

Rs  1 / t

(2.3.1)tx (1, j, j, k )  (2.0  eps  sig ) /( 2.0  eps  dt  sig )(2.3.2)tx (1, j, j, k )  (2.0  eps  sig ) /( 2.0  eps  dt  sig )(2.3.3)tx (1, j, j, k )  (2.0  eps  sig ) /( 2.0  eps  dt  sig )(2.3.4)tx (2, j, j, k )  2.0  dt /( 2.0  eps  dt  sig )(2.3.5)tx (2, j, j, k )  2.0  dt /( 2.0  eps  dt  sig )(2.3.6)tx (2, j, j, k )  2.0  dt /( 2.0  eps  dt  sig )(2.3.7)厚さ t m,導電率σ S/m, 面抵抗値 Rs Ω
真空の導電率の補正因子の式

誘電率 eps, 導電率 sig
FDTD 法の電磁界の時間配置を考えたとき,式(2.2.1)及び(2.2.2)で記したマクスウェル方程式を,𝐁 = 𝜇0 𝐇, 𝐃 = 𝜀0 𝐄, 𝐉 = σ𝐄を用いて解くと,
𝜕𝑬
𝜎1= − 𝑬+ ∇×𝑯
𝜕𝑡
𝜀0
𝜀0

(2.3.8)𝜕𝑯
1= − ∇×𝑬𝜕𝑡𝜇0となる。前述のように,時間ステップが n  1 2 での磁場 H
電場 E n と１ステップ前の磁場 H
場 E n1 は，半ステップ前の磁場 H

n

1
2

n

1
2

(2.3.9)n

12は，その半ステップ前でのから求めることができ，時間ステップが n  1 での電と１ステップ前の電場 E n から求めることができるため,電界及び磁界の時間微分は,
𝜕𝑬
𝑬𝑛 − 𝑬𝑛−1|=
𝜕𝑡 𝑡=(𝑛−1∆𝑡)
∆𝑡
2

38

(2.3.10)1

1

𝜕𝑯
𝑯𝑛+2 − 𝑯𝑛−2|=
𝜕𝑡 𝑡=𝑛∆𝑡∆𝒕(2.3.11)となる.このとき式(2.3.10),(2.3.11)を式(2.3.8),(2.3.9)に代入すると,
1
1
𝑬𝑛 − 𝑬𝑛−1𝜎1= − 𝑬𝑛−2 + ∇ × 𝑯𝑛−2
(2.3.12)
∆𝑡
𝜀0
𝜀0
1

1

𝑯𝑛+2 − 𝑯𝑛−2
1= − ∇ × 𝑬𝑛∆𝑡
𝜇0(2.3.13)となる.さらに,式(2.3.12)の右辺第 1 項の電界はt = (n − 1⁄2)∆tにおける値となっているため,σ𝑬𝑛−1⁄2 を近似すると次のようになる.
σ𝑬𝑛−1⁄2 ⇒ 𝜎𝑬𝑛 − 𝑬𝑛−1
2

よって,𝑬𝑛 , 𝑯𝑛+1⁄2 についてまとめると,
1
2𝜀0 − 𝜎∆𝑡 𝑛−12∆𝑡
𝑬𝑛 = −
𝑬+∇ × 𝑯𝑛−22𝜀0 + 𝜎∆𝑡
2𝜀0 + 𝜎∆𝑡1

1

𝑯𝑛+2 = 𝑯𝑛−2 −1∇ × 𝑬𝑛𝜇0(2.3.14)(2.3.15)(2.3.16)となる.この時の式(2.3.15)の右辺第 1 項,第 2 項を真空の導電率の補正因子の式としている.

39

(a)

(b)

(c)

図 2.3.3 抵抗被膜の近似モデル

図 2.3.4 絶縁体の挿入イメージ図

40

(a)

(b)
図 2.3.5 導電率の影響範囲(2-14)

41

参考文献
(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)
高村秀一：“プラズマ理工学入門”(森北出版, 1997)(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) 橋本修, 安部琢美：“FDTD 時間領域差分法入門” (森北出版, 1996)(2-15) 宇野亨：“FDTD 法による電磁界およびアンテナ解析” (コロナ社, 1998)(2-16) K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.
(2-17) G. Mur, IEEE Trans. Electromagnetic Compat., 40 (1998), 100.
(2-18) 著=W. H Hayt，訳=山中惣之助，岡本考太郎，宇佐美興一：“工学系の基礎電磁気学”(朝倉書店, 2001)(2-19) V. P. Gropinath and T. A. Grotjohn, IEEE Trans. Plasma Sci., 23, (1995),602.
(2-20) T. A. Grotjohn, Rev. Sci. Instrum. 67, (1995), 921.
(2-21) Y. Saad, and M.H. Schultz , SIAM. J. Sci. Stat. Comput., 7, pp.856-869 (1986)(2-22) 藤野清次，張紹良：”反復法の数理”, (朝倉書店, 1996)(2-23) 樫山和男，西村直志，牛島省：”並列計算法入門”, (丸善, 2003)(2-24)
新谷將， “マイクロ波放電式小型イオンエンジンの数値解析”, 平成 22 年度九州大学修士論文．
(2-25) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical Physics,
55 (1971), 4100.
(2-26) M. Srendra, D. B. Graves and G. M. Jellum, Physical Review A, 41 (1990),
1112(2-27) M. Hirakawa, IEPC-97-021, (1997), 159.
42

(2-28) http://www.asl-i.com/index.html

43

第3章

結果と考察

3.1 絶縁体の挿入
図 3.1 にアンテナからマイクロ波が入射していく様子を表すために,電界強度,図 3.2に時間遷移する電子エネルギー分布を示す.図はそれぞれ,アンテナの中心の yz 面
(x=16.5 mm),xy 面(z=10.5 mm)の断面図を見たものであり,時間は 7.75 ns 時を採用した.
セラミックを入れたとき,アンテナ近傍の電界強度の高い範囲が小さくなった.FDTD 法
では出力が 8 W になるようにポインティングベクトルによって投入電力が調整されている.そのため,セラミックを挿入したことで電力の反射が抑えられ,結果として,投入電力が減尐し,その結果電界強度の高い範囲が小さくなったと考えられる.これは,電子がエネルギーの減尐につながり,その結果として,電子の平均エネルギーも同様にセラミックを入れる前はアンテナ近傍で 30 eV の電子温度が見られたが,セラミックを入れた
とき,アンテナ近傍では 18 eV と,低くなっている.
また一方で,セラミックを入れたときアンテナから離れたエリアでは電界強度と電子のエネルギーが 60000 V/m から 100000 V/m へと全体的に高くなった.これは電力の反射が減ったことで,今までアンテナ近傍でしか安定したプラズマが見られなかったが,より広い範囲で安定したプラズマが存在できるようになったと考えることができる.

44

図 3.1 (a) 電界強度

（セラミックなし）

図 3.1 (b) 電界強度

（セラミックあり）

45

図 3.2(a) 電子エネルギー分布 （セラミックなし）

図 3.2(b) 電子エネルギー分布

46（セラミックあり）

3.2

抵抗被膜の挿入

絶縁上にスパッタ粒子が堆積した時の電界強度を図 3.3 に,電子のエネルギー分布を図
3.4 に示す.それぞれには異なる導電率を代入していて,中和器を稼働してから(a)稼働直後(b)7500 時間後(c)150000 時間後の状態である.時間は 4.1 ns 時を採用した.
中和器の稼働時間が長くなるほど,アンテナ近傍及び中和器全体の電界強度は弱くなった.これは時間の経過とともにセラミック上にスパッタ粒子が堆積したことで,マイクロ波のエネルギーが堆積層によって消費され,結果的にプラズマに投入されるエネル
ギーが減尐し,電子の平均エネルギーも低くなったと考えられる. このとき,電界強度は中和器稼働直後では z=8～12 mm の範囲で 100000 V/m を示していたが,稼働後 150000
時間後には 100000 V/m を示す範囲は z=10～12 mm と範囲が狭まり,電子温度は中和器稼働直後ではアンテナ周りに 30 eV 程度の値が見られたが,稼働後 150000 時間後では 10 eV程度の値しか見られなかった.これにより,プラズマ生成率が下がり,中和器内部でのプラズマ生成密度が小さくなる.よって,電子を引き出すための電圧を上げる必要があり,
シースの拡大,さらなるスパッタ粒子の生成と,悪循環を引き起こし,寿命を早めてしまう.

47

図 3.3 (a)

電界強度(稼働直後)

図 3.3 (b)電界強度(7500 時間後)48

図 3.3 (c)電界強度(150000 時間後)49

図 3.4 (a)電子エネルギー分布(稼働直後)

図 3.3 (b)電子エネルギー分布(7500 時間後)50

図 3.4電子エネルギー分布

51(150000 時間後)

第4章

結論

本研究では,スパッタ粒子がマイクロ波伝送へ与える影響を解明することを目的とし,1 セルを従来の 2 倍の細かさ(0.25 mm×0.25 mm×0.25 mm)に細かくして,絶縁体の挿入及びスパッタ粒子を,抵抗被膜を用いて表現した数値解析コードの開発を行った.
FDTD 法と PIC 法を用いたカップリングコードを改良し,電界強度及び電子のエネルギ
ー分布の計算を行った.抵抗被膜の抵抗値で電界強度及び,電子温度を比較することで,スパッタ粒子が絶縁体に付着し,導電性が生まれ,電界強度は中和器稼働直後では z=8～
12 mm の範囲で 100000 V/m を示していたが,稼働後 150000 時間後には 100000 V/m を示
す範囲は z=10～12 mm と範囲が狭まり,電子温度は中和器稼働直後ではアンテナ周りに30 eV 程度の値が見られたが,稼働後 150000 時間後では 10 eV 程度の値しか見られなかった.これにより,プラズマ生成率が下がり,中和器内部でのプラズマ生成密度が小さくなる.よって,電子を引き出すための電圧を上げる必要があり,シースの拡大,さらなるスパッタ粒子の生成と悪循環を引き起こしてしまう.以上からスパッタ粒子が中和器の寿命を低下させている要因の一つと結論づけることができる.
今後は,現在の計算はプラズマの状態が一定である条件で計算しているため,プラズマ解析によるデータを用いることが必要である.また 1 セルを 2 倍にしたことで 1 セルに入る粒子数が尐なくなった.よって粒子数の数と解析時間を増やすために計算コストの削減が必要とされる.

52

謝辞
最後に執筆に当たって,適切なご指導と丁寧なご校閲をいただいた中島秀紀教授,山本直嗣准教授,森田太智助教授に深く感謝の意を表するとともに,本研究にあたり,公私にわたり助言をいただいた楠田将士氏,並びに中島研究室関係者へあわせて謝意を表す
る.

53