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

修 士 論 文論文名
マイクロ波放電式小型中和器の二次元数値解析

氏

名

下川床

潤

指導教員名

中島

秀紀

目次
第 1 章 序論
1.1 背景1

1.2 イオンエンジンの原理

2

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

2

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

41.3 研究目的

6

第 2 章 数値解析
2.1 粒子計算
2.1.1 基礎式

8

2.1.2 メッシュサイズ

92.1.3 PIC 法

10

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

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

112.1.6 平均自由行程と透過度

122.1.7 衝突確率・衝突判定

13

2.1.8 衝突計算

13

2.1.9 エネルギーと散乱角142.1.10 Null-collision 法

172.2 電磁波解析手法
2.2.1 基礎式

192.2.2 差分式

192.2.3 アンテナ解析

212.2.4 安定条件

222.2.5 PEC 境界条件

22

2.2.6 Mur の吸収境界条件

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

23

2.2.8 補正ポテンシャル

24

2.3 カップリングコード

28

2.3.1 計算体系292.3.2 計算格子

30

2.3.3 時間刻み

30

2.3.4 電子初期条件

30

2.3.5 粒子数

302.3.6 質量比の変更

312.3.7 磁場強度

312.3.8 境界条件

32

2.3.9 計算コードの並列化

32

第 3 章 コード評価
3.1 エネルギーバランス

34

3.2 粒子軌道

35

3.3 計算速度36

第 4 章 計算結果と考察
4.1 電子密度分布374.2 イオン密度分布394.3 電界強度分布404.4 電子エネルギー分布

424.5 粒子損失45

4.6 電離衝突位置

46

4.7 実験との比較
第 5 章 結論
謝辞

第1章

序論

1.1 背景天気予報やカーナビゲーションシステム、携帯電話等は現在の我々の生活において欠かすことの出来ないものになっており、今後も更に宇宙の商業利用は活発化していくと考えられる。それに伴い、小型の商業用人工衛星に関する研究(1-1)(1-2)が進められている．
従来の大型衛星は重要なミッションを担うために失敗が許されず信頼性が重視されていた．そのため使用部品も宇宙環境に耐えうる高価な専用品でなければならず
膨大な開発費用が高く，さらに耐久テスト等で開発期間も必要であった．特に衛星所有にかかる膨大な初期投資は企業の参入を難しくし、宇宙の商業利用の障害となっていた．
これに対して小型衛星は信頼性よりも低コスト・短期開発を目標としている．宇宙用部品に民生部品を流用することによって費用を抑え，また，先端技術を使用した部品を投入することも可能である．更に小型衛星はコンパクトで重量が軽いため
一回打ち上げで複数の衛星を宇宙空間に輸送することができるという点でも有利である．このような衛星開発は宇宙ビジネスをさらに活性化し、我々の生活がより便利になるようなサービスが生まれることに繋がっていく．小型人工衛星の抱える問題とし、そのサイズの制約が挙げられる。省スペース化するために、多くは、モーメンタムホイール(1-3)や磁気トルカ(1-3)，コールドガスジェット(1-3)等の低性能のエンジンが搭載されているだけである．
人工衛星や宇宙ステーションなどの地球周回物は，地球が楕円球であることや大気抵抗の存在，また太陽や月の引力による摂動があるため，その姿勢制御・軌道制御のために推進機を搭載している．現在，衛星の寿命を決定するものは機器の耐年数よりもむしろ推進剤の枯渇である．特に小型衛星においてはサイズの制約上，搭載可能な推進剤の量も制限されてしまう．従って，少量の推進剤で長期間運用可能な燃費の良い推進機の開発が重要となる．そして，その要求に応えるものとして利用されているのが電気推進である．
電気推進(1-4)とは，電気エネルギーを推進剤に与えることで反力を得る推進機であり，従来の化学反応を利用したロケット推進と違って静電力や電磁力を利用できる．
このため比推力（推進剤流量あたりの推力）を一桁以上大きくすることが可能であり，大幅な推進剤の低減が望める．現在使用されている代表的な電気推進機はアー
クジェットスラスタやホールスラスタ，そしてイオンエンジンである．イオンエン
ジンは他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制御や惑星探査に適している．

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 =

3eBm(1.2.3)

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

E

B

電子

E

B

電子

電子
E

B
B
電子

E

図 1.3 ECR の原理

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

4

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

図 1.4

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

5

図 1.5

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

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

1.3 研究目的
これまで本研究室でマイクロ波放電式小型イオンエンジンの磁場強度・放電室形状などを変化させて実験を行ってきた(1-7)(1-8)(1-9)．しかし，形状を変更させたことによって，放電室内部の物理メカニズムがどのように変化し，どのような理由で性能に影響をおよぼしているのかは解明されていない．その原因の一つとして中和器自体が小型であるために従来のプローブ法等の装置での計測では放電室内部のプラ
ズマに擾乱を引き起こし，正確な測定が困難であることが挙げられる．レーザ計測
6等(1-11)を用いれば計測が可能であるが，放電室にレーザ掃引用の光路を確保しなければならず，一つのスラスタで複数の位置で計測することは難しい．
その様な場合に数値解析は非常に有効な手段となる．近年の計算機の発展のおかげで実験では知り得る事のできない情報を入手できるようになってきた．それにより，実験による時間的・金銭的コストを抑えた研究開発ができるようになってきている．本研究室では新谷らによって，数値解析コードが開発されてきた(1-11)(1-12)．
このコードでは計算時間の関係で定常状態に至っていないという問題があり，また各種パラメータに関する依存性の調査もまだまだこれからの課題である．
本研究の目的は，マイクロ波放電式小型イオンエンジンの内部プラズマの数値解
析コードに二次元化を含めた改良を施すこと，また，そのコードを用いて投入マイ
クロ波周波数による依存性の調査を行うことである．

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

H. Sahara, S. Nakasuka and C. Kobayashi, AIAA paper 2005-3956, 2005.
茂原正道・鳥山芳夫, “衛星設計入門”, 2002.
栗木恭一・荒川義博[編]：”電気推進ロケット入門”, (東京大学出版会,2003)(1-5) N. Wallace, P. Jameson, IEPC, 2011, 327(1-6) I. Funaki, H. Kuninaka and K. Toki, J. Propul. Power, Vol. 20, (2004), pp.718-726.
(1-7) N. Yamamoto, H. Kataharada, H. Masui, H. Ijiri and H. Nakashima,
AJCPP2005-22093, 2005.
(1-8) N. Yamamoto, H. Kataharada, T. Chikaoka, H.Masui and H. Nakashima,
IEPC-2005-036, 2005.
(1-9) 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-10) 鶴哲平， “マイクロ波放電式中和器の性能向上に関する研究” , 平成 20 年度九州大学修士論文．
(1-11) 増井博一， “マイクロ波放電型宇宙推進用プラズマ源に関する数値解析”,
平成 17 年度九州大学博士論文．
(1-12) 新谷將 “マイクロ波放電式小型イオンエンジンの数値解析”, 平成 22 年度
九州大学修士論文．

7

第 2 章 数値解析
2.1 粒子計算
粒子計算(2-1)
2.1.1

基礎式

プラズマの挙動を解析する際の基礎式として運動方程式を用いる．式を以下に示
す．

m

dv
= q(E + v × B )dt(2.1.1)dx=v
dt(2.1.2)

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

(2.1.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 )6(2.1.4)∆t(bz1 + 2bz 2 + 2bz3 + bz 4 )6
但しここで定数 b はそれぞれ以下を用いる。
vzn +1 = vzn +q(Er + vφ Bz ) + 1 vφ vφmr∆tbφ 1 ∆tbφ1 ∆tbφ1   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.1.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.1.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.1.7)

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

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

9(2.1.10)∆t ≤ ω pe−1

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

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

(2.1.13)に対して，∆t << ω c−1

(2.1.14)である必要がある．
2.1.3

PIC 法(2-4)

シミュレーションでは，電磁場量を空間格子点上で求める．そのため，計算粒子の電荷，電流を空間の格子点上に分配し，マクスウェル方程式を用いて電磁場を計算しなければならない．また，粒子が感じる電磁場を計算する際，格子点上に割り振られた電磁場を内挿する必要がある．このような計算を行うため，粒子シミュレーションにおいて標準的な手法である PIC 法（Particle In Cell method）を用いる．
まず，粒子の位置から近傍のセル座標 (i, j ) を求める．次に，図 2.1 に示す

dr1, dr 2, dz1, dz2 を求める．その値を用いて，図中にあるような長方形 S1～S 4 の面積を求める．この面積比から粒子が感じる電磁場を内挿することができる．
具体的には，座標 (r, z ) 上にある粒子が感じる電場 E (r , z ) は，計算セルの体積を S ，
図に示すような粒子によって仕切られた空間体積をそれぞれ S1～S8 とすると，格子点上の重みは

An =
となり，

Sn
S

(n = 1 ～ 4)

E (r , z ) = A1E (i, j ) + A2 E (i + 1, j ) + A3 E (i + 1, j + 1) + A4 E (i, j + 1)と表わされる．

10(2.1.15)(2.1.16)図 2.1.1

2.1.4

PIC 法におけるセルの構成

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

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

S1 qN(2.1.17)
S Vとなる．すべての粒子について和をとれば，格子点状での電荷密度が得られる．また，電流密度は
j = qnv(2.1.18)

ρ (i, j ) =より，

j (i, j ) =

2.1.5

S1 qN
v
S V

(2.1.19)

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

荷電粒子は電磁場の影響を受けて，中性粒子と衝突する．プラズマは電子，イオン，中性粒子からなるとする．電子とイオンの数密度は中性粒子の数密度より非常に小さいため，電子-電子，電子-イオン，イオン-イオンの各衝突は無視できるものとする．電子と中性粒子は弾性，励起，電離衝突をすると仮定する．なお，以下

11

に用いる乱数は(0,1)間に一様に分布する乱数である．生成周期の短い組み込み乱数
ran()の使用は避け Mersenne Twister(2-5)を使用した．周期は(219937-1)であり，生成速
度も早い．

2.1.6

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

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

λ = ∫ lP (l )dl
0

=

1
Nσ

(2.1.24)

これを平均自由行程という．よってモンテカルロ法によって，任意の点から出発し最初の衝突が起こるまでの距離 l を決めるには，乱数 r を用いて，
r = P (l ) l
= 1 − exp −  λ(2.1.25)で表わされ，これを変形して，

l = −λln(1 − r )となる．ここで， r は一様乱数であるから 1 − r を r とすれば，
l = −λln(r )(2.1.26)(2.1.27)と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を σ とすれば，この粒子が１秒間に気体分子と衝突する回数 υ は，

12

υ=

v

λ
= σNv

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

2.1.7

衝突確率・衝突判定
衝突確率・衝突判定

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

υ total = ngas σ total v m(2.1.30)

Pcollision,m = 1 − exp(− ngasσ total v m ∆t )= 1 − exp(− υ total ∆t )(2.1.31)

ここで，ngas は中性粒子の数密度，vm は粒子 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.1.32)v ′2 =

1(v1 + v 2 + u′)2(2.1.33)と表わすことができる．ここで， 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)．
m + M cos χM sin χv e′ =
ve +
h
m+M
m+M
h の直交座標成分は

14(2.1.34)hx = ver cos φ
hy = −
hz = −(2.1.35)

vex vey cos φ − ve vez sin φ(2.1.36)v er
vex vez cos φ − ve vey sin φ(2.1.37)ver

ここで，ver = vey2 + vez2

(2.1.38)ve = vex2 + v ey2 + vez2

(2.1.39)であり， χ は散乱角，φ は任意の角である．(2.1.34)式はエネルギー損失を含んでいる．
次に，衝突前後のエネルギー収支を考える．衝突前の電子のエネルギーを Eincident,e とし，衝突後の電子のエネルギーを Escattered,e とすれば，エネルギー収支は以下のように表わされる．
E scattered, e = E incident, e(2.1.40)

電子は，エネルギー Eincident, e が小さいとき等方散乱され，大きいときは主に前方散乱となる．これを考慮した散乱角 χ の確率密度を g (χ ) として以下に示す(2-25)．ここで， Eincident,e の単位を eV とする．
Eincident ,e sin χg (χ ) =χ
21 + E incident ,e sin 2 ln 1 + E incident ,e 2(0 ≤ χ ≤ π )(2.1.41)

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

2 1 − (1 + Eincident, e ) 3
Eincident, er}(2.1.42)となる．この値を(2.1.34)式に用いて衝突後の速度を求める．このとき，任意の角は
φ = 2π r4(2.1.43)
として与える．
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーを Eexcitation とし，入射電子のエネルギーを Eincident,e ，励起後のエネルギーを Escattered,e ，励起後の速度を v~ とする．エネルギー収支は，

15

E scattered, e = E incident, e − E excitation(2.1.44)と表わされ，励起後の速度はE
v~ = v e 1 − excitation
E incident, e(2.1.45)となる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失った後に(2.1.45)式中 v~ の速度で弾性衝突をすると考えることができる．(2.1.34)式
中 v e に v~ を用いて，衝突後の速度 v e′ を求める．このときの散乱角には(2.1.42)式を用いるが(2.1.42)式中の Eincident, e には(2.1.44)式から求められる E scattered, e を用いる．励起エネルギー E excitation は，キセノンの場合 8.34 eV を用いている．
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを Escattered,e ，生成電子のエネルギーを Ecreated,e ，電離エネルギーを Eionization とする．エネルギー収支はE scattered, e + E created, e = E incident, e − Eionization(2.1.46)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる(2-12)．
 E incident, e − E ionization E scattered,e = B tan r5 tan -1 
2B




(2.1.47)

ここで，B はガスに固有の値でキセノンでは 8.7 eV となる．
また，入射電子が電離によって失うエネルギーを ∆E とすると，∆E = E incident, e − E scattered, e(2.1.48)と表わされる．(2.1.45)式から電離後の電子の速度は，
v~ = v e 1 −∆E
E incident, e(2.1.49)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．(2.1.36)
式中 v e に v~ を用いて，衝突後の速度 v e′ を求める．散乱角には(2.1.42)式を用いるが，(2.1.42)式中の Eincident,e には(2.1.47)式から求められる Escattered,e を用いる．電離エネルギ

ー Eionization は，キセノンの場合 12.13 eV である．
電離後の生成電子の持つエネルギーは，E created, e = ∆E − E ionization

16(2.1.50)

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

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

ネルギー E created, e を用いて(2.1.42)式から求める．
2.1.10Null-collision 法(2-11)(2-12)(2-13)

計算時間短縮のため，Null-collision 法を用いる．これには，架空の衝突断面積 σ fakeを導入する．すべての電子のエネルギー E に対して
σ fake < σ total(2.1.52)であるとし，例えば電子と中性粒子の衝突においてはσ total = σ elastic + σ excitation + σ ionization + σ fake(2.1.53)とする．ここで，σ total v =υ collisionn gas= constant(2.1.54)を満たす σ 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 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度で規υn格化した衝突周波数を示す．本計算では Null-Collision 法で用いた collision gas の値
は 2×10-12 m3/sec である．

17

図 2.1.3 Xe の数密度で規格化した電子の Xe の衝突断面積182.2 電磁波解析手法(2-14)(2-15)
2.2.1

基礎式

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

∇ × H = ε0

∂E
+J
∂t

(2.2.1)(2.2.2)

ここで， E は電場， H は磁場， µ 0 は真空の透磁率， ε 0 は真空の誘電率， J は電流密度を表わす．
円筒形の放電室内を伝播するマイクロ波が周方向に m 個のモードをもつとき、マクスウェル方程式を円筒座標系のもとに展開すると以下のようになる。
∂Er1  ∂H φ m= − 
± H z − σEr + J r 
r∂tε  ∂z∂Eφ∂t=

1  ∂H r ∂H z
− σEφ − J φ −∂rε  ∂z(2.2.3)1  1 ∂ (rH z ) m∂E zm H r + σE z + J z 
=− 
∂t
r
ε  r ∂r
1  ∂Eφ m 
∂H r
= 
± E z µ  ∂z
∂t
r

∂H φ
∂t

=

1  ∂E z ∂Er −µ  ∂r∂z (2.2.4)∂H z
1  ∂ (rEφ ) m = − 
m H r ∂tµ  ∂r
r


2.2.2

差分式

三次元の FDTD 法で用いられる電場，磁場の空間配置を図 2.2.1 に示す．このように，電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．これは Yee によって提案され(2-16)，電場の回転

19（ ∇ × E ）が磁場を，磁場の回転（ ∇ × H ）が電場の方程式を満たすように配置されている．
この電磁場の配置のもと，(2.2.3)式，(2.2.4)式を空間，時間とも中心差分した式を以下に示す．
（(2.2.5)~(2.2.11)式）ここで， ∆r , ∆z は空間幅， ∆t は時間幅，添え字(i, j ) はセル中での座標，添え字 n は時間ステップを表わす．
E rn +1 (i, j ) = E rn (i, j ) +n +1 2
n +1 2∆t  H φ (i, j ) − H φ (i, j − 1)  m−
H zn +1 2 (i, j ) − J r (i, j )

ε 0ε r ∆z ri +1 2

Eφn +1 (i, j ) = Eφn (i, j ) +∆t  H rn +1 2 (i, j ) − H rn+1 2 (i, j − 1)   H zn +1 2 (i, j ) − H zn+1 2 (i − 1, j )  − 
 − J φ (i, j )
ε 0ε r 
∆z
∆r

 


E zn+1 (i, j ) = E zn (i, j ) +(2.2.5)n +1 2
n +1 2∆t  1  ri +1 2 H φ (i, j ) − ri −1 2 H φ (i − 1, j )  m n+1 2− H r (i, j ) + J z (i, j ) 
 r
ε 0ε r  ri 
∆r

 i

(2.2.6)(2.2.7)

H rn +1 2 (i, j ) = H rn −1 2 (i, j ) − Eφn (i, j + 1) − Eφn (i, j ) ∆t  m n


(
)
E
i
j
,
−



µ 0 µ r  ri z∆z

(2.2.8)Hφn +1 2 (i, j ) = Hφn −1 2 (i, j ) −∆t  Ern (i, j + 1) − Ern (i, j )   E zn (i + 1, j ) − E zn (i, j )  − 


∆z
∆r
µ0 µr  
(2.2.9)H zn +1 2 (i, j ) = H zn −1 2 (i, j ) + 1 ri +1 Eφn (i + 1, j ) − ri Eφn (i, j ) ∆t  m n


(
)
E
i
,
j
−

 ri +1 2µ 0 µ r  ri +1 2 r∆r
(2.2.10)

特に、r=0 の位置に配置された磁場 Hz は以下の差分式を使用する。
4∆t
H zn +1 2 (0, j ) = H zn −1 2 (0, j ) −Eφn (1, j )µ 0 µ r ∆r(2.2.5)~(2.2.11)式から，時間ステップが1n + 1 2 での磁場 H
n−

n+

1
2

(2.2.11)は，その半ステップ前

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

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

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

21

図 2.2.2
2.2.4

アンテナモデル

安定条件

計算結果の正確さを保つため，空間の増分値 ∆r , ∆z は波長 λ に比べて十分小さくなければならない．（ λ 10 以下）
時間ステップ ∆t は，Courant の安定化条件より
1
∆t ≤
1
1
v
+
2
(∆r ) (∆z )2(2.2.12)を満たさなければならない．この条件は極めて厳しく，わずかでも満足しなければ計算が不安定になる．
2.2.5 PEC(Perfect Electric Conductor )境界条件
境界条件
境界面が金属のような完全導体である場合，境界に平行な電場成分を 0 とする
PEC 境界条件を用いる．放電室内壁は導体を仮定しているので，この境界条件を適

用する．
2.2.6

Mur の吸収境界条件(2-17)

開放領域を取り扱う場合，計算領域を仮想的な境界で閉じなければならない．この境界を吸収境界といい，この吸収境界が完全でなければ，反射波が解析領域に戻ってしまい，誤差の原因となる．ここでは，定式化が容易で広く使用されている
Mur の吸収境界条件を用いる．
z = n∆z の吸収境界に E r 成分を持つ平面波が入射する場合を考える．伝播速度をc とすると， z 方向に進む波は

22

Er = Er (z + ct )(2.2.13)と表わされ，次の微分方程式を満足する．
∂Er 1 ∂Er−=0∂z c ∂t(2.2.14)z = n∆z の境界で完全に反射がないならば，電場は(2.2.13)式の形を保ったまま伝播

するはずである．すなわち，境界でも(2.2.14)式が満足される．これを FDTD 法の形式に定式化するため，時間について差分化すると∂Er Ern − Ern−1
=
∂t
∆t
n−

(2.2.15)

1

∂E 2
=v r
∂z

となる．右辺の z についての微分は，電場が z = (n − 1)∆z, n∆z に割り当てられていることから， z = (n − 1 2)∆z で行い，(2.2.14)式に代入すると
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.16)となる。(2.2.15)式の Ern (i, n − 1 2 ), Ern−1 (i, n − 1 2 ) も FDTD 法では直接求めることはできないため，前後の値で平均して求める．そして， Er についてまとめるとErn (i, n ) = Ezn−1 (i, n − 1) +c∆t + ∆z n
Er (i, n − 1) − Ern−1 (i, n )c∆t − ∆z{}(2.2.17)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．
2.2.7 ポインティングベクトル(2-18)
空間を流れる電力を求めるためには，ポインティングベクトルを使用する．
Maxwell の方程式を変形すると，以下の式となる．

− ∫ ( E × H ) ⋅ dS = ∫ J ⋅ Edv +svol∂  εE 2 µH 2 ∫  2 + 2 dv∂t vol

(2.2.18)

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

23

2.2.8 補正ポテンシャル式(2.2.6)～(2.2.8)で計算される電界 E にはプラズマ中の荷電粒子の移動によって発生する電流 J の項が反映されている．電流 J の算出方法は 2.1.4 節で述べた．しかし，この方法で求めた J は厳密にはガウスの法則を満足しておらず荷電粒子の位置から計算される電荷密度ρとの間に僅かな差が生じる（参考文献）．E を補正する
ために補正ポテンシャル∇∇δφを導入する．補正後の電界を E’とすればE ′ = E − ∇δφ(2.2.19)
ガウスの法則を満たすには∇ ⋅ E′ =ρε0∇ ⋅ ( E − ∇δφ ) =(2.2.20)ρ
ε0(2.2.21)補正ポテンシャルについて解けば以下の偏微分方程式を得る．
∇ 2δφ = ∇ ⋅ E −ρ
ε0(2.2.22)

円筒座標では以下のようになる．
∂ 2δφ 1 ∂  ∂δφ   ∂E r E r ∂E z  e+++
 − (ni − ne )
r
=
∂z 2
r ∂r  ∂r   ∂rr∂z  ε 0(2.2.23)φは空間電位，ε0 は真空中の誘電率，ni と ne はイオン・電子の数密度である．
差分化すれば
左辺＝

1 
1 
 2
 2
 2 +
δφr + ∆r , z +  2 −δφr − ∆r , z
r∆r 
r∆r 
 ∆r
 ∆r
+

右辺＝

2
2
2 
 2
δφr , z + ∆z + 2 δφr , z − ∆z −  2 + 2 δφr , z
2
∆z
∆z
∆z 
 ∆r

1
(Err+∆r ,z − Err−∆r ,z ) + 1 Err ,z + 1 (Ezr ,z+∆z − Ezr ,z−∆z ) − ρ r ,z
2∆r
r
2∆z

(2.2.24)(2.2.25)

連立方程式を Ax=b に帰着させる。
r 方向のメッシュ数を Nr 個、z 方向のメッシュ数を Nz 個とすると未知数 Nr×Nz
個を並べたベクトル x は以下のように作る。

24

 δφ1,1 
 M 


 δφ Nr ,1  δφ1, 2 
 M 
x=

 δφ Nr , 2  M  δφ 
 1, Nz 
 M 
δφ

 Nr , Nz 

(2.2.26)

同様に(2.2.23)の右辺からベクトル b を作る。ここで(2.2.25)を Ri , j と書くようにす
れば R1,1  M  RNr ,1 
 R1, 2 
 M 
b=

 RNr , 2 
 M 
R 
 1, Nz 
 M 
R

 Nr , Nz 

(2.2.27)

次に(2.2.24)から行列 A を作る。(2.2.24)を次のように書き換える。
(2.2.24)=a1δφ r + ∆r , z + a 2 δφ r − ∆r , z + a3δφ r , z + ∆z + a 4δφ r , z − ∆z + δφ r , z(2.2.28)

但し、
1 
 2
 2 +

r∆r 
∆r
a1 = − 2 
 2 2 + 2∆z  ∆r

(2.2.29)

1 
 2
 2 −

r ∆r ∆r
a2 = −2 
 2
 2 + 2∆z  ∆r

(2.2.30)

a3 = a 4 = −

2
∆z 2
2 
 2
 2 + 2∆z  ∆rである。（ a5 の係数を 1 にした。）

25(2.2.31)

Ax=b に帰着させるため、A は以下のような行列になる。

 a3


 1 a1O a2 O O

α =
,β = OO O a1
a2 1 
a3 


(2.2.32)

α β



β O OA=O O β
β α 


(2.2.33)として、場 の 方 程 式 を 解 く に あ た り ， GMRes(Generalized Minimal Residual method) 法(2-21)(2-22)(2-23)(2-24)を採用した．
この方法は非対称行列を有する線型方程式の反復解法であり，次のような考え方にもとづいている．解くべき方程式はAx = b
(2.2.34)
である．ここに A は(N×N 行列)であり，簡単のために実行列であるとする．b は既知の N ベクトル，x は未知の N ベクトルである．いま，適当に初期過程 x0 を選ぶと，(2.2.24)の残差 r0 は
r0 = b - Ax0(2.2.35)
る．ここに zn は，Krylov 空間
Kn(A,r0) = Span(r0,Ar0,…,An-1r0)(2.2.36)から||b - Axn|| = ||r0 - Azn||(2.2.37)を最小にするように選ぶものとする．記号 Span()は，括弧内のベクトルの線形結合で作られる空間を表す．
(2.2.27)の最小化を行うために，Kn(A,r0)の正規直交基底(v1,v2,…,vn)を次のようにして逐次構成する．まず||r0|| = γ とすると v1 = r0/γ となる．いま，第 n 番目の基底
vn までが求められているとする．このとき，
wn=Avn, hin=vi・wn(2.2.38)とおき，Schmidt の直交化により
Vn+1 = wn - Σhinvi(2.2.39)とおけば，Vn+1 ⊥ Kn(A,r0)である．さらに Vn+1 = Vn+1/||Vn+1||とすると，
hn+1 n := vn+1・wn = vn+1・Vn+1 = ||Vn+1||(2.2.40)を得る．
次に，(2.2.27)における最小化を実行する．第 n 回目の繰り返しにおいては zn∈

26

Kn(A,r0)であるから，zn = Σvjyj(2.2.41)と書ける．これを用いて(2.2.27)を v1,…,vn+1 座標系で書くと，(i = 1,2,…,n+1)

min || γδi1 – Σhijyj||(2.2.42)が得られる．上式の行列 hij は(n+1)×n の上 Hessenberg 行列である．

 h11

 h21
 0
Qn 
 M
 0

 0


h12
h22
h32

h13
h23
h33

L
L
L

O
L

O
0

O
hn−1 nL00h1n   h11
 
h2 n   0
h3n   0
=M   Mhnn   0
hn +1 n   0

h12
h22
0
O
L
L

L h1n 
h23 L h2 n 
h33 L h3n 
O O M 

0 0 hnn 
0 0 0 h13

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

 h11

 0
 0

 M

 0
 0


h12
h22
0
O
L
L

h13 L h1n  y1 
 Q11n  n
h23 L h2 n  
 Q21  y2 

h33 L h3n  
y 3 = γ  Q31n O O M  


M


M

0 0 hnn   n
 yn 
 Qn1 0 0 0 (2.2.44)を後退代入によって解くことにより構成され，(2.2.32)の残差 En は|γQnn+1 1|と求められる．以上のプロセスを残差が許容限度 ε に対して En<ε を満たすまで繰り返して，
近似解 xn を得る．以上が，GMRES の基本的なアルゴリズムである．
電場は空間電位の勾配より求めた．ただし，円周方向成分については無視した．

27

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

28

図 2.3.1

Coupling code のフローチャート

2.3.1 計算体系
図 2.3.2 は当研究室で開発されている小型マイクロ波放電型中和器の断面図であり, 計算では水色の領域を取り扱う.

29

図 2.3.2

計算体系

2.3.2 計算格子
PIC 法を用いるために格子を刻む必要があるが，そのために望ましい長さは，デバイ長である．そうすることによって，壁面でのシースが計算できる．しかし，計算コストの点からデバイ長で格子を刻むことは現時点では不可能である．今回は電子とマイクロ波との相互作用を調べることが目的であるからデバイ長で計算格子を作る必要はない．そこで本計算では計算格子の長さを 0.25 mm とし，z 方向に 176
格子 r 方向に 42 格子の計 7392 格子作成した．
2.3.3

時間刻み

PIC 法 の 時 間 刻 み は ∆tpic = 1.0 × 10−12∆tFDTD = 2.0 × 10

−14

sec ， FDTD 法 の 時 間 刻 み はsec とすれば、2.1.2 や 2.2.4 で述べた安定条件を満たす．FDTD

コードが 50 回計算された後，PIC コードを 1 回計算させた．
2.3.4 電子初期条件
Ｔe＝0.5eV となるマクスウェル分布の初期速度を持った電子を放電室内に一様に分布させた．
2.3.5 粒子数
1 ボクセル当たり数個しか粒子がない場合，ノイズが大きくなり望ましくない(2-12)．通常，100 個以上の粒子が存在しているとよいとされている．今回は 1 ボク
セルに平均 130 個程度の粒子が含まれるように粒子数を設定した．
30

2.3.6 質量比の変更(2-27)
今回の計算において，電子の速度成分は最大で 105 m/s 程度，イオンにおいては最大 103 m/s 程度であり，電子の速度を基準に時間刻みを設定するとイオンはほとんど動かない．そこでイオンの質量をα分の一まで軽くする．エネルギーは実際の値と同じとするため速度が α 倍となる．これによって領域内の物理現象を早く収束させ計算時間を短縮できる．

mi _ virtual =1αmi _ real(2.3.1)

ただしこれにより物理現象が影響を受けないようにするため以下のように修正を行った．
Rarmor 半径はmv(2.3.2)eB
で与えられるので、イオンの感じる磁束密度をα分の一にすることによって仮想値と現実値をおなじとする．これはイオンの移動量が修正されてしまうと，固定された幾何形状の計算領域に対して分布が変わってしまう可能性があるためである．
2.3.7 磁場強度
今回はマイクロ波周波数依存性を調べるために電源として 900MHz、1200MHz、1600MHz、2.45 GHz を用いる。マイクロ波周波数と電子サイクロトロン周波数が一致するときの磁場強度はそれぞれ 32mT、43mT、57mT、87.5mT である。その値を示す領域を ECR 層と呼ぶ。この ECR 層の位置が各周波数で同じになるように磁石の数を変えて磁場強度を調節し計算を行った．磁場計算には Advanced Science
Laboratory, Inc の 3 次元静磁場解析ソフトを使用した(2-28)．図 2.7 にそれぞれの周波数毎の ECR 位置を示す．
表 2.1
rc =

周波数

放電室外側の磁石数

フロントヨークの磁石数

Case1

900 MHz

6個

3個

Case2

1200 MHz

8個

4個

Case3

1600 MHz

13 個

3個

Case4

2.45 GHz

15 個

3個

31

(a) case1 のときの ECR 領域(b) case2 のときの ECR 領域(c) case3 のときの ECR 領域(d) case4 のときの ECR 領域

図 それぞれの磁場強度における ECR 領域
2.3.8 境界条件
粒子に関しては,軸に達したものはその反対側でも同様に粒子が存在するものと考え，鏡面反射させた．絶縁されているアンテナに衝突した粒子は計算体系から除外するが，電荷は残るとした．計算領域右端はグリッドがあるので，電子はグリッドの開口率で反射，イオンに関してはすべて吸収とした．ヨークに衝突したものと計算領域の上部，左端に達したものについては計算から除外した．場については軸
ではノイマン境界条件とし，そのほかではディリクレ境界条件を用い 0 V と設定した．
2.3.9 計算コードの並列化(2-24)
計算にかかる時間を短縮することはより多くの条件に対して検証するために非常に重要である。演算の高速化の為，並列計算が有効である（参考文献）．並列計算を実現されるためにはコンパイラに特別な指示文やサブルーチン，環境変数を与える必要があるが，本研究では OpenMP という規格に従い共有メモリシステム上でスレッド並列演算を行えるようにプログラムを書き加えた．OpenMP と共に普及している分散メモリシステム上で演算を行う為の標準規格に MPI がある．OpenMPには MPI よりも逐次計算プログラムから並列計算プログラムへの書き換えが簡易であるという利点をもつ．

参考文献
(2-1)
内藤裕志, プラズマ・核融合学会誌, 74 (1998) ,470.
(2-2)
Giovanni Lapenta, Jeremiah U. Brackbill, Journal of Computational Physics, 115,
(1994), 213.
32

田中實, 山本良一：“計算物理学と計算化学―分子動力学法とモンテカルロ法” (海文堂,1998)

(2-3)

(2-4)
(2-5)

石黒静児, プラズマ・核融合学会誌, 74 (1998) , 591.
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) 橋本修, 安部琢美：“FDTD 時間領域差分法入門” (森北出版, 1996)(2-15)(2-16)(2-17)(2-18)

宇野亨：“FDTD 法による電磁界およびアンテナ解析” (コロナ社, 1998)
K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.
G. Mur, IEEE Trans. Electromagnetic Compat., 40 (1998), 100.
著=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)
藤野清次，張紹良：”反復法の数理”, (朝倉書店, 1996)(2-23) 樫山和男，西村直志，牛島省：”並列計算法入門”, (丸善, 2003)(2-24) 牛島省：”OpenMP による並列プログラミングと数値計算法”, (丸善, 2006)(2-25) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical Physics, 55(1971), 4100.
(2-22)(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.

(2-28)

http://www.asl-i.com/index.html

33

第 3 章 計算コード評価
3.1 エネルギーバランス
エネルギー収支の確認を行った。放電室でのプラズマ生成でのエネルギー収支式はWmicrowave _ in = ∆U electron + ∆U ion + W microwave _ out(3.1.1)となる．ここで、 Wmicrowave _ in はマイクロ波として投入されたエネルギー， ∆U electron は単位時間あたりの電子のエネルギー増加量， ∆U ion は単位時間あたりのイオンのエネルギー増
加量， Wmicrowave _ out は流出したマイクロ波のエネルギーである． (3.1.1)を書き換え、
Wbalance を以下のように定義し監視した．
Wbalance = Wmicrowave _ in − ∆U electron − ∆U ion − Wmicrowave _ out(3.1.2)
図 3.1 に投入電力 2W のときのそれぞれのエネルギーとそのバランスを示す．それぞれの項は投入マイクロ波 5 周期分の平均値を取りプロットしている． Wbalance = 0 が常に成り立っておりエネルギーは系全体で収束していることが確認できた．

図 3.1 エネルギーバランス

34

3.1.2 粒子軌道
粒子の運動が正しく計算されているかどうかを確認した．図 3.2 は磁場中を運動する電子の軌道である．電子はラーマー運動をしながら磁場に巻き付いて移動し，放電室内を左右に往復運動していることが分かる．一方，イオンはラーマー半径が非常に大きいため，磁力線に巻き付いていない．両者とも運動中にエネルギーの変化はなく PIC 法によって運動方程式が正しく計算されていることが確認できた．

図 3.2 電子の軌道

図 3.3 イオンの軌道

図 3.4 電子とイオンのエネルギーの時間変化

35

3.1.3 計算速度
これまでの三次元カップリングコードでは計算コストが大きく，プラズマ生成の初期段階までしかシミュレーションすることが出来なかった．本研究では位置の上方を一次元減らした二次元カップリングコードに改良した．両コードを用いてを同じ条件の計算を行い，演算にかかった時間を比較した．二次元コードでは円筒座標を採用し，円周方向への物
理現象は対照的であるとの仮定を用いたためメッシュ数が大幅に少なくなっている．尚、計算では電子のみを扱った．計算条件を表 3.1 に示す．
表 3.2 に三次元のカップリングコードの計算時間を 100 としたときの二次元カップリング
コード，また二次元カップリングコードをメモリ共有システムを用いて 8 スレッド並列計算した場合の計算時間を比較した．特に FDTD 部分は二次元化によって 1%以下の時間で計算が完了しており，高速化に成功している．二次元化によってほとんど高速化できなかった粒子計算については並列化することで約 71%計算時間を短縮することに成功した．
表 3.1 計算速度を比較する際の計算条件
メッシュ数
総メッシュ数
メッシュサイズ
粒子数 [個]FDTD タイムステップ
PIC タイムステップ
計算時間

[個]
[個]
[mm]
[個]
[sec]
[sec]
[n sec]

三次元
x:60×y:60×z:32
115,200
x:0.25×y:0.25×z:0.25
1,000,000
1×10-14
5×10-130～0.1

表 3.2 三次元カップリングコードを 100 としたときの計算時間の比較
三次元
二次元
FDTD
54.1
0.31
PIC
45.9
42.1
100.0
42.5
全体

二次元
r:30×z:32960
r:0.25×z:0.25
1,000,000
1×10-14
5×10-130～0.1

二次元(8 並列)
0.27
12.2
12.5

120
PIC

Calculation time [-]

100

FDTD

80
60
40
20
0
3D

2D

図 3.5 計算時間の比較
36

2D(8Parallel)

37

第 4 章 結果と考察
4.1 電子密度分布
図 4.1 に各マイクロ波周波数における電子密度分布を示す。900 MHz，1200 MHz では磁気チューブの形成が確認できる．一方，1600 MHz，2450 MHz ではチューブの影は薄く，放電室中程に均等に広がっている．
バックヨークに比べフロントヨーク部の密度が濃いのは，磁力線が集中しており，磁束密度が高いため，磁気ミラー効果が大きくなった為だと考えられる（磁力線は図 3.2 に示した）．
アンテナ近傍を見ると，周波数が低いほど密度も低いことが分かる．2 章の表 2.1 に示したように低い周波数の計算条件ほど磁石数が少なく磁束密度が低いので，電子の束縛が弱く，半径方向へ移動しやすい．アンテナ付近の電子は放電室中央へ拡散したり，アンテナへ衝突して消滅し，結果として密度が薄くなったものと考えられる．

37

図4.1 電子エネルギーの空間分布（上から 900 MHz，1200 MHz，1600 MHz，2450 MHz，）

38

4.2 イオン密度分布
図 4.2 はイオン密度分布を示したものである．電子密度分布に近い分布をとっていることが分かる．アンテナ近傍の密度が低い．しかし，後述のイオンの損失数ではアンテナに到達して消滅したイオンは少ない．このことから，イオンはアンテナに衝突したのではなく，アンテナに吸収されなかった放電室中程の電子に引っ張られて移動したと考えられる．

図4.2 電子エネルギーの空間分布（上から 900 MHz，1200 MHz，1600 MHz，2450 MHz，）

39

4.3 電界強度分布
電界強度分布図 4.3 に電界の絶対値を時間平均し，その空間分布を 10000 V/m，6000 V/m，2000
V/m の等高線で表示した．900 MHz では 10000 V/m 以上の強電界領域が放電室中程まで広がっている．これに対し，2450 MHz では強電界領域は非常に狭くマイクロ波が浸透していないことが分かる．4.1 節で確認したように 900 MHz においてはアンテナ付近で電子密度が低く，逆に 2450 MHz では密度が濃い．そのためマイクロ波がカットオフされたと考えられる．

40

図4.3 電界強度の空間分布（上から 900 MHz，1200 MHz，1600 MHz，2450 MHz，）

41

4.4 電子エネルギー分布
電子エネルギー分布
図 4.4 に空間電子エネルギーの空間分布を示す。900 MHz ではフロント－バックヨーク
間のアンテナ近傍において電子エネルギーが高く，キセノンの電離エネルギーである
12.13eV 以上の電子が広範囲に分布していることが分かる．2450 MHz においても高エネルギーの電子がアンテナ近傍に分布していることが確認できるが，900 MHz の場合よりも全体に低温であり，キセノンの電離エネルギーに達している領域は少ない．4.3 節で示したようにマイクロ波が放電室内に届かず電子の加熱に使われなかったからだと考えられる．
図 4.5 に投入電力 Wmicrowave _ in と流出電力 Wmicrowave _ out との差 Wabsorbed の周波数毎の時間得推移を示す．1600 MHz，2450 MHz では系に吸収された電力が小さいく，マイクロ波
はアンテナ近傍のプラズマに反射されることを示唆している．

42

図4.4 電子エネルギーの空間分布（上から 900 MHz，1200 MHz，1600 MHz，2450 MHz，）

43

図4.5 系に吸収されたエネルギー44

4.5 粒子損失
フロントヨーク，バックヨーク，内壁，アンテナ，に到達したものと外部に排出された電子・
イオンについてカウントした．電子の損失数を図 4.6 に，イオンの損失数を図 4.7 に示す．
電子イオンともに周波数が高くなるにつれて損失は減少しており，特に電子にこの傾向が強い．また，電子の損失数に比べ，イオンの損失数のほうが二倍近くまで増えている．
電子は周波数増加とともに磁場閉じ込めの効果が強くなったため損失が減った．イオンについては磁場に閉じ込められた電子に引っ張られる効果によって電子と同傾向が出たが，その規模は小さかった．
35000
900 MHz
30000
25000
20000

1200 MHz
1600 MHz
2450 MHz

15000
10000
5000
0
Frontyork

Backyork

Antenna

Inner wallExit

図4.6 電子の損失数
70000
900 MHz
60000

1200 MHz
1600 MHz

50000

2450 MHz
40000
30000
20000
10000
0
Frontyork

Backyork

Antenna

Inner wall

図4.7 イオンの損失数

45

Exit

4.6 電離衝突位置
電離衝突位置
電離衝突位置の結果を図 4.8 に示す．900 MHz，1200MHz では放電室全体で電離が起きているのに対し，1600 MHz，2450 MHz ではフロント－バックヨーク間の領域での電離が集中している．プラズマ生成位置も磁場閉じ込めの度合いに依存し，磁場強度が強
いほど壁面から離れた位置で電離を引き起こすことが分かった．

図4.8 電離衝突位置（上から 900 MHz，1200 MHz，1600 MHz，2450 MHz，）

46

4.7 実験との比較
図 4.9 に実験によって得られた中和器の引き出し電子電流と，同体系で計算した電子
エネルギーの周波数依存関係を示す．実験では ECR 層をバックヨーク上部に設置するために表 4.1 のような磁石配置を用いた．そのときの ECR 層の位置を図 4.10 に示す．
投入電力は 2 W，推進剤流量は 0.05 sccm である．
実験では周波数が高いほど引き出し電子電流が高く，1600 MHz のときに 16.6 mA と最も性能がよかった．また，2W，2450 MHz ではプラズマが点火しなかった．
シミュレーションでは 900 MHz のときがもっともプラズマの加熱効率がよく，実験と逆の傾向になった．プラズマの生成効率と中和器の中和性能（引き出し電子電流）が直結していておらず，密度の分布や磁場閉じ込めの効果による電子の引き出しにくさ等が複雑に関係しているためと考えられる．今後の課題として，電子引き出しまでをモデルに組み込んだ計算コードに発展させることが必要である．

図 4.9 実験値との比較
表 4.1
周波数
900 MHz
1200 MHz
1600 MHz
2.45 GHz

放電室外側の磁石数
4個
6個
9個
13 個

47

フロントヨークの磁石数
3個
3個
3個
3個

4.10 実験で用いた体系の ECR 位置(a)900 MHz，(b)1200 MHz，(c)1600 MHz，(d)2450 MHz

48

第 5 章 結論
本研究ではマイクロ波放電式小型中和器の内部物理を解明するために，三次元数
値解析プログラムの二次元化に取り組んだ．また，それを用いてプラズマとマイクロ波の相互作用について調査した．得られた結論を以下に示す．
1) FDTD-PIC コードの二次元化及び並列化によって計算速度の向上に成功した．
FDTD は二次元化により約 100 倍の速度に，PIC は並列化によって約 4 倍の速度になり，全体ではこれまでの約 8 倍の速度が得られた．
2) FDTD-PIC コードを使用した二次元数値解析によってマイクロ波周波数と電子の加熱の関係について調査した．
1. 高い周波数では磁場強度が強いため電子が束縛の影響を強く受け，アンテ
ナ付近のプラズマ密度が高くなると分かった．
2. 高い周波数ではアンテナ近傍のプラズマによってマイクロ波が反射されていることが示された．
3. 高い周波数では放電室内部に強電界領域が形成されないため，電子は加熱されにくく，2 W においては低い周波数のほうが効率よく加熱できることが示された．
3) 磁場強度が強まることで閉じ効果が増し，プラズマの生成が壁面から離れた領域で行われることが分かった．

謝辞
本研究を進めるにあたって，数多くの方々にお世話になりました．
中島先生・山本さんには研究の方向性に迷ったり問題に直面した際には的確な助言によって解決に導いていただきました．充実した研究を行うことが出来ました．ありがとう御座いました．
研究は異なりましたが中島研究室の皆さんには公私共々いろいろとお世話になりました．この場を借りて感謝いたします．