卒	
  業	
  論	
  文

題	
  	
  目	
  	
 

マイクロ波イオンエンジンの中和器における電磁場解析

氏 	
  	
  	
  	
  名

楠田	
  将士

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

指 導 教 員

中島	
  秀紀	
  教授
山本	
  直嗣	
  准教

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

目次
第１章

序論 .................................................................................... 1

1.1

背 景 .............................................................................................. 1

1.2

イ オ ン エ ン ジ ン の 原 理 ...................................................................... 21.2.1

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

1.2.2

マイクロ波放電式小型中和器 ............................................................................. 4

1.3

研 究 目 的 ........................................................................................ 7

第２章

数値解析 .............................................................................. 9

2.1

粒 子 計 算 (2-1) ................................................................................... 92.1.1

基礎式 ................................................................................................................. 9

2.1.2

メッシュサイズ・タイムステップ(2-3) ............................................................... 10

2.1.3

PIC 法(2-4) ........................................................................................................... 112.1.4

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

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

2.1.6

平均自由行程と透過度(2-5) ................................................................................. 13

2.1.7

衝突確率・衝突判定 ......................................................................................... 142.1.8

衝突計算 ........................................................................................................... 14

2.1.9

エネルギーと散乱角 ......................................................................................... 15

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

電 磁 波 解 析 手 法 (2-14)(2-15) ................................................................. 192.2.1

基礎式 ............................................................................................................... 19

2.2.2

差分式 ............................................................................................................... 19

2.2.3

アンテナ解析 .................................................................................................... 22

2.2.4

安定条件 ........................................................................................................... 24

2.2.5

PEC(Perfect Electric Conductor )境界条件 ..................................................... 24

2.2.6

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

2.2.7

ポインティングベクトル(2-18) ............................................................................ 25

2.3

カ ッ プ リ ン グ コ ー ド ....................................................................... 26

2.3.1

計算体系 ........................................................................................................... 28

2.3.2

計算格子 ........................................................................................................... 28

2.3.3

時間刻み ........................................................................................................... 29

2.3.4

電子初期条件 .................................................................................................... 29

2.3.5

粒子数 ............................................................................................................... 29

2.3.6

推進剤 ............................................................................................................... 29

2.3.7

磁場強度 ........................................................................................................... 29

2.3.8

計算コードの並列化(2-24) ................................................................................... 30

第３章

計算結果と考察 ................................................................... 32

3.1

投 入 電 力 ...................................................................................... 32

3.2

電 界 強 度 分 布 ................................................................................ 33

3.3

電 子 エ ネ ル ギ ー 分 布 ....................................................................... 35

3.4

電 子 エ ネ ル ギ ー 獲 得 確 率 ................................................................. 37

第４章

結論 .................................................................................. 42

謝辞 ............................................................................................. 43

第１章	
 
1.1

序論

背景

	
  昨今，天文ブームが高まっている．昨年においては，金環日食，金星の太陽面通過，
その他流星群に関する情報がメディアから多く発信されている．そのきっかけとも言えるのが，2010 年 6 月 13 日人類初の小惑星から物質を持ち帰った探査機通称「はやぶさ」の帰還である．このはやぶさの帰還によって，小惑星イトカワから物質を持ち帰ったという世界初の偉業に加えて，世界で初めてμ10 と呼ばれるマイクロ波放電式イオンエン
ジンが宇宙実証され，累積稼働時間は約 40000 時間を超え，その有用性の高さが実証された．JAXA を始めとした日本の技術がこのような実績を残すことに成功した要因として，日本のコスト管理と省エネ技術の高さがある．この技術の背景には衛星の小型化の研究の推進(1-1)(1-2)がある．
	
  はやぶさのような大規模なミッションとは別に，我々の身近なところでは，携帯やカーナビシステムの GPS 機能や天気予報等我々の生活において欠かすことのできないものになっており，その機能の役目を担っているのが，地球を周回している人工衛星である．このような宇宙の商業利用は活発化していき，より多くの衛星が求められていくであろう．しかし，従来の大型衛星の開発コストは高く，企業が易く参入できない状況であった．これに対して，小型衛星は開発期間が短くなるためより新しい技術を導入しや
すい．単純に小さいことによるコストの削減とそれに伴う制作回数の増加による打ち上げ失敗等のリスクの分散できるという利点がある．
	  人工衛星や宇宙ステーションなどの地球周回物は，地球が楕円球であることや大気抵
抗の存在，また太陽や月の引力による摂動があるため，その姿勢制御・軌道制御のために推進機を搭載している．現在，衛星の寿命を決定するものは機器の耐年数よりもむしろ推進剤の枯渇である．特に小型衛星においてはサイズの制約上，搭載可能な推進剤の量も制限されてしまう．従って，少量の推進剤で長期間運用可能な燃費の良い推進機の開発が重要となる．そして，その要求に応えるものとして利用されているのが電気推進である．
	  電気推進(1-3)とは，電気エネルギーを推進剤に与えることで反力を得る推進機であり，従来の化学反応を利用したロケット推進と違って静電力や電磁力を利用できる．このため比推力（推進剤流量あたりの推力）を一桁以上大きくすることが可能であり，大幅な推進剤の低減が望める．現在使用されている代表的な電気推進機はアークジェットスラ
スタやホールスラスタ，そしてイオンエンジンである．イオンエンジンは他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制御や惑星探査に適している．
惑星探査においては，イオンエンジンは比推力が高い一方で，推力密度が低いというデメリットを持つために，十分な速度が得られるまで多くの時間が求められる．したがって，長時間運用が必須となり，長寿命化が求められている．

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 はラーマ半径と呼ばれ

rL =

mv ⊥
eB

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

ω ce =
3eBm(1.2.3)

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

E

B

電子

E

B

電子

電子
E

B
B

電子

E

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

ヒーターを用いないため，予備加熱が不要となり素早い起動が出来る．

2.

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

4

先に述べたようにマイクロ波放電式中和器ははやぶさに搭載され十分な中和性能と寿命が得られたという実績がある．しかし，小型衛星に搭載するようなサイズの中和器は未だ実用化されておらず，その開発が急務となっている．
本研究室では小型衛星への搭載を想定した小型マイクロ波放電型イオンエンジンを開発してきた．エンジン本体について優れた性能が得られた．一方中和器においては要求されるサイズ，消費電力でも中和性能の目標値をクリアしておらず，現在改善に向けて力を入れている．
図 1.4 は本研究において開発中のマイクロ波放電型小型中和器の外観写真である．放電室は内径 21 mm の円筒形状であり，オリフィス開口径 10 mm，放電室長 32 mm である．放電室外部と出口部分に Sm-Co 磁石を配置し構成された磁気回路によって ECR に必要な磁場を供給する．アンテナはモリブデン製である．

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

5

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

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

6

1.3

研究目的

これまで本研究室でマイクロ波放電式小型イオンエンジンの磁場強度・放電室形状などを変化させて実験を行い，(1-6)(1-7)(1-8) 中和器の性能評価及び内部プラズマの計測を実施してきた．しかし，形状を変更させたことによって，放電室内部の物理メカニズムがどのように変化し，どのような理由で性能に影響をおよぼしているのかは解明されていな
い．その原因の一つとして中和器自体が小型であるために従来のプローブ法等の装置での計測では放電室内部のプラズマに擾乱を引き起こし，正確な測定が困難であることが挙げられる．レーザ計測等(1-10)を用いれば計測が可能であるが，放電室にレーザ掃引用の光路を確保しなければならず，一つのスラスタで複数の位置で計測することは難しい．
そのような場合に数値解析は非常に有効な手段となる．近年の計算機の発展のおかげ
で実験では知り得る事のできない情報を入手できるようになってきた．それにより，実験による時間的・金銭的コストを抑えた研究開発ができるようになってきている．本研究室では，下川床らによって中和器の内部物理の解明および，中和器の性能を評価する
数値解析コードの開発を行なってきた(1-10)(1-11)(1-12)(1-13)．これまでの課題として，これま
で電子のエネルギー授受は ECR 層で効率的に行われていると考えられてきたが，実験結果より，ECR 層以外でもエネルギーのやり取りが行われている可能性が示唆された．これを解析するために，従来のコードでは，後述する電子のエネルギー獲得確率は広い領域でしか解析することができなかった．そのために電子がどの領域でエネルギーをやり取りしているかが大まかにしか解析できなかった．そのためにより細かい領域ごとの解析が必要となった．解析領域の細分化には粒子数を増加する必要性，及び観測物理量の増大に伴う計算コストの増大が課題であった．更に，計算限界によって定常状態まで計算できていない問題や各種パラメータに関する依存性はマイクロ波周波数やアンテナ
位置については調査されてきたが，引き続きコードの改良とともに実施していく必要がある問題がある．
	  本研究の目的は，マイクロ波放電式イオンエンジンの内部プラズマの数値解析コードの改良及び，電子のエネルギー授受に関しての細かな領域解析を目的とする．

7

参考文献
(1-1)

M.Kato, S. Takayama, U. Nakamura, K. Yoshihara and H. Hashimoto, 56th

International Astronautics Congress paper IAC-05.B5.6.B.01, 2005.
(1-2)

H. Sahara, S. Nakasuka and C. Kobayashi, AIAA paper 2005-3956, 2005.

(1-3)

栗木恭一・荒川義博[編]：”電気推進ロケット入門”, (東京大学出版会,2003)(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)

N.

Yamamoto,

H.

Kataharada,

H.

Masui,

H.

Ijiri

and

H.

Nakashima,

AJCPP2005-22093, 2005.
(1-7)

N. Yamamoto, H. Kataharada, T. Chikaoka, H.Masui and H. Nakashima,

IEPC-2005-036, 2005.
(1-8)

T. Chikaoka, S. Kondo, N. Yamamoto, H. Nakashima and Y. Takao, Proceedings ofthe 25th International Symposium on Space Technology and Science, 2006, pp.254-259.
(1-9)

鶴哲平， “マイクロ波放電式中和器の性能向上に関する研究” , 平成 20 年度九

州大学修士論文．
(1-10)

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

成 17 年度九州大学博士論文．
(1-11)

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

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

州大学修士論文．
(1-13)

下川床潤，“マイクロ波放電式小型中和器の二次元数値解析”，平成 23 年度九

州大学修士論文．

8

第２章	
 

数値解析

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)と書ける．ここで， α ・ β ・ γ ・ δ を次のように定義する．	
 

α=

qΔt n
Bx
2m

(2.1.6)

β=

qΔt n
By
2m

(2.1.7)

γ =

qΔt n
Bz
2m

(2.1.8)9δ = 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)式を次のようにして変形し，速度を求めることができる．	
 

2.1.2

v n +1 2 = A −1 S

(2.1.12)⎛ 1 + α 2 αβ + γ αγ − β ⎞⎟1 ⎜A −1 = ⎜ αβ − γ 1 + β 2 βγ + α ⎟
δ ⎜
2 ⎟
⎝ αγ + β βγ − α 1 + γ ⎠(2.1.13)

メッシュサイズ・タイムステップ (2-3)

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

λDe =

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

Δ ≤ λDe(2.1.15)であることが望ましい．また，タイムステップはプラズマ周波数 ωpe，	
 ne e 2me ε 0ω pe =(2.1.16)に対して，	
 Δt ≤ ωpe−1(2.1.17)とすれば，数値的な不安定性は生じない．時間積分の精度まで考えると，通常は	
 10

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

(2.1.19)に対して，	
 Δt << ωc−1

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

ωc
100

(2.1.21)となるように設定した．
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 を求める．その値を用いて，図中にあるような直方体 V1V8 の体積を求める．この体積比から粒子が感じる電磁場を内挿することができる．	
 	  具体的には，座標 (x , y , z )上にある粒子が感じる電場 E (x , y , z )は，計算セルの体積を V ，
図に示すような粒子によって仕切られた空間体積をそれぞれ V1～V8 とすると，格子点上
の重みは
An =

Vn
V

(n = 1  8)

(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)+ A8 E (i, j +1, k +1)と表わされる．	
 11(2.1.23)	 図 2.1.1	
  PIC 法におけるセルの構成
2.1.4

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

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

V1
qN
V2

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

V1
qNv
V2

12

(2.1.26)

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

	  荷電粒子は電磁場の影響を受けて，中性粒子と衝突する．プラズマは電子，イオン，
中性粒子からなるとする．電子とイオンの数密度は中性粒子の数密度より非常に小さい
ため，電子-電子，電子-イオン，イオン-イオンの各衝突は無視できるものとする．電子と中性粒子は弾性，励起，電離衝突をすると仮定する．なお，以下に用いる乱数(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 以下で最初に衝突が起こる確率は，l

P (l ) = ∫ Pinitial (l )dl0= 1− exp (−Nσ 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 )(2.1.33)で表わされ，これを変形して，となる．ここで， r は一様乱数であるから 1 − r を r とすれば，
l = −λln(r )13(2.1.34)

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

υ=

v

λ
= σNv

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

衝突確率・衝突判定

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

σ = σ (E )(2.1.36)

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

υ total = ngasσ total v m(2.1.37)

Pcollision, m = 1 − exp(− ngasσ total v m Δt )(2.1.38)= 1 − exp(− υ total Δt )ここで， ngas は中性粒子の数密度， v m は粒子 m の速度， σ total は全衝突断面積である．ある乱数 r1（ 0 < r1 < 1）を用いて， Pcollision, m > r1 のときに粒子 m はタイムステップ Δt の間に衝突が起こるとする．
2.1.8

衝突計算

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

∗ 0 ≤ r2 ≤σ elastic
σ total

であれば弾性衝突∗

σ elasticσ+ σ excitation≤ r2 ≤ elastic
σ totalσ totalであれば励起衝突∗

σ elastic + σ excitationσ+ σ excitation + σ ionization≤ r2 ≤ elastic
σ totalσ totalであれば電離衝突

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

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

14

図 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)．
v!e =m + M cos χM sin χ
ve +
h
m+M
m+M
15

(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 = vey2 + vez2

(2.1.45)ve = vex2 + vey2 + vez2

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

Escattered, e = Eincident, e(2.1.47)

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

g(χ ) =

Eincident,e sin χχ
2 1+ Eincident,e sin 2 ln 1+ Eincident,e2{(0 ≤ χ ≤ π )}(2.1.48)

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

2 1− (1+ Eincident,e ) 3
Eincident,e}(2.1.49)となる．この値を(2.1.41)式に用いて衝突後の速度を求める．このとき，任意の角は

φ = 2π r4(2.1.50)として与える．
②励起衝突	
  基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーを
Eexcitation とし，入射電子のエネルギーを E incident, e ，励起後のエネルギーを Escattered,e ，励起後の速度を v~ とする．エネルギー収支は，Escattered, e = Eincident, e − E excitation
16(2.1.51)

と表わされ，励起後の速度はE
v~ = 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)
式中の E incident, e には(2.1.51)式から求められる Escattered,e を用いる．励起エネルギー E excitation は，
キセノンの場合 8.34 eV を用いている．
③電離衝突
	  入射電子のエネルギーを E incident, 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)式中の E incident, e には(2.1.54)式から求められる Escattered,e を用いる．電離エネルギー Eionization は，キセ
ノンの場合 12.13 eV である．
電離後の生成電子の持つエネルギーは，E created,e = ΔE − Eionizationとなる．また，入射電子の速度 v e を用いて生成電子の速度は

17(2.1.57)

v created,e =veve

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

E created,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)とする．ここで，σ total v =υ 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 である．

18

図 2.1.3 Xe の数密度で規格化した電子の Xe の衝突断面積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 ）が磁場を，19

磁場の回転（ ∇ × 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 ) +Δtµ 0 (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 ) '%
+
$
(
Δz
Δy
&%
)%
n+

1

n−

1

H y 2 (i +1 2, j, k +1 2 ) = H y 2 (i +1 2, j, k +1 2 ) +Δtµ 0 (i +1 2, j, k +1 2 )%# 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) %'
+
$
(
Δx
Δz
&%
)%
n+

1

n−

1

H z 2 (i +1 2, j +1 2, k ) = H z 2 (i +1 2, j +1 2, k ) +Δtµ 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 ) %'
+
$
(
Δy
Δx
%&
%)
E xn +1 (i + 1 2 , j , k ) = E xn (i + 1 2 , j , k ) +

+

1

n+

1

H y 2 (i + 1 2 , j , k − 1 2 ) − H y 2 (i + 1 2 , j , k + 1 2 )E yn +1 (i , j + 1 2 , k ) = E yn (i , j + 1 2 , k ) +(2.2.4)	 (2.2.5)	 (2.2.6)	 (2.2.7)	 Δtε 0 (i + 1 2 , j , k )1n+⎧ n + 12⎪ H z (i + 1 2 , j + 1 2 , k ) − H z 2 (i + 1 2 , j − 1 2 , k )
⎨
Δy
⎪
⎩
n+

(2.2.3)Δz⎫⎪− J x (i + 1 2 , j , k )⎬
⎪
⎭

Δt
ε 0 (i , j + 1 2 , k )1n+⎧ n + 12⎪ H x (i , j + 1 2 , k + 1 2) − H x 2 (i , j + 1 2 , k − 1 2)
⎨
Δz
⎪
⎩
1
1
n+
n+
⎫
H z 2 (i − 1 2 , j + 1 2 , k ) − H z 2 (i + 1 2 , j + 1 2 , k )⎪+− J y (i , j + 1 2 , k )⎬
Δx
⎪
⎭

20

E zn +1 (i , j , k + 1 2) = E zn (i , j , k + 1 2) +Δt
ε 0 (i , j , k + 1 2)

1n+
⎧ n + 12
2()(i − 1 2 , j , k + 1 2)
H
i
+
1
2
,
j
,
k
+
1
2
−
H
⎪ y
y
⎨
Δx
⎪
⎩

(2.2.8)

	
 

1
1
n+
n+
⎫
H x 2 (i , j − 1 2 , k + 1 2) − H x 2 (i , j + 1 2 , k + 1 2)⎪+
− J z (i , j , k + 1 2)⎬
Δy
⎪
⎭

n+

1

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

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

21

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.11)式，(2.2.12)式より，V を求め，(2.2.9)式，(2.2.10)式よりアンテナに与える電場を求める．ここで，Ex を 𝑟 = 𝛥𝑥/2の電場で，Ey を 𝑟 = 𝛥𝑦/2 の電場で代表させると，

E x (IFED, JFED, KFED) = −E x (IFED −1, JFED, KFED)=−

V (nΔt)Δx # Δx &ln % (
2 $ a '

(2.2.13)E y (IFED, JFED, KFED) = −E y (IFED, JFED −1, KFED)=−

V (nΔt)Δy # Δy &ln % (
2 $ a '  ここで，式中の(𝐼𝐹𝐸𝐷, 𝐽𝐹𝐸𝐷, 𝐾𝐹𝐸𝐷)は電場を励振する座標を与える．

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

図 2.2.3	  アンテナモデル23

2.2.4安定条件

	  計算結果の正確さを保つため，空間の増分値 Δr , Δz は波長 λ に比べて十分小さくなけ
ればならない．（ λ 10 以下）
	  時間ステップ Δt は，Courant の安定化条件より
1
Δt ≤
1
1
1
v
+
+
2
2
(Δx ) (Δy ) (Δz )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
Δt
n−

1

∂E 2
=v z
∂x

(2.2.18)となる．右辺の x についての微分は，電場が x = 0, Δx に割り当てられていることから，

x = Δx 2 で行い，(2.2.17)式に代入すると

24

n
r

E (i, n − 1 2) − E
Δt

n−1
r

n−

1
2

n−

1
2

(i, n − 1 2) = v Er (i, n) − Er (i, n − 1)Δz(2.2.19)となる．(2.2.18)式の Ezn (1 2), Ezn −1 (1 2)も FDTD 法では直接求めることはできないため，前後の値で平均して求める．そして， E z についてまとめるとEzn ( 0, j, k +1 2 ) = Ezn−1 (1, j, k +1 2 )+vΔt − Δx n
{Ez (1, j, k +1 2) − Ezn−1 (0, j, k +1 2)}
vΔt + Δx(2.2.20)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．
2.2.7 ポインティングベクトル (2-18)
	  空間を流れる電力を求めるためには，ポインティングベクトルを使用する．Maxwellの方程式を変形すると，以下の式となる．

∂ ⎛ εE 2 µH 2 ⎞⎟⎟dv− ∫ (E × H ) ⋅ dS = ∫ J ⋅ Edv + ∫ ⎜⎜
+
∂
t
2
2
⎠
s
vol
vol⎝

(2.2.21)

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

25

2.3

カップリングコード

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

26

図 2.3.1

Coupling code のフローチャート

27

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

図 2.3.2	
  計算体系

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

28

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 コードが 4 回計算された
後，PIC コードを 1 回計算させた．
2.3.4 電子初期条件
	  Te＝0.3eV となるマクスウェル分布の初期速度を持った電子を放電室内に一様に分布させた．本来であれば，FDTD−PIC コードとは別に，プラズマの挙動を観測するコードの結果を元に電子密度が決定されるが，今回においては，まだ導入に至っていない．
2.3.5

粒子数

	
  1 ボクセル当たり数個しか粒子がない場合，ノイズが大きくなり望ましくない(2-12)．
通常，10 個から 100 個以上の粒子が存在しているとよいとされている．今回は 1 ボクセ
ルに平均 25 個程度の粒子が含まれるように粒子数を設定した．
上記の計算体系における体積 5.46×10-6 m3，実験による粒子密度は 10-17~10-18 m-3 と分かっているので，平均化して 5×10-17m-3 とした．超粒子数は 100 万個として，超粒子 1 つが何個代表するかは，体積と粒子密度より超粒子 1 つは約 27 万個を代表すると決定した．
2.3.6 推進剤
推進剤はキセノンとし，中性粒子密度は 3.86×1019 /m3， 推進剤流量は 0.018mg/sとした．
2.3.7 磁場強度	  今回は電源としてはやぶさにも採用されている 4.2GHz を用いる．マイクロ波周波数と電子サイクロトロン周波数が一致するときの磁場強度は 0.15T である．その値を示す
領域を ECR 層と呼ぶ．磁場計算には Advanced Science Laboratory, Inc の 3 次元静磁場解
析ソフトを使用した(2-28)．図 2.3.2 に上述の計算体系における磁場強度及び ECR 位置を示す．

29

ECR layor

図 2.3.2 磁場強度分布

2.3.8 計算コードの並列化 (2-24)
	  計算にかかる時間を短縮することはより多くの条件に対して検証するためには非常に重要である．1 台のコンピュータ(1 つの CPU)で行っていた計算をネットワークでつながれた複数のコンピュータ(複数の CPU)で計算できるように改良することにした．用いたのは MPI(Message Passing Interface) ，分散メモリ型の並列計算機で 2 つのプロセス 間でのデータをやりとりするために用いるメッセージ通信操作の仕様標準である．

30

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

(2-28)

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

第３章	
 
3.1

計算結果と考察

投入電力

投入電力の確認を行う．先に述べたとおり投入電力はアンテナヨークの放電質側開口部から流入するエネルギーをポインティングベクトルを計算することによって求める．
図 3.1 に投入エネルギーの変化を示す．この傾きが投入電力にあたる．これを計算した
ところ，投入電場 66000 V/m に対して 0.9 W となった．これは増井らの計算と比較するとかなり低い値であり，今後見直す必要がある．

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

32

3.2

電界強度分布	  図 3.2 にアンテナからマイクロ波が入射していく様子を表すために電界強度を示す．
これは 125ns の時間平均をとり，アンテナの中心の(a)yz 面(x=16.5mm)，(b)xy 面(z=10.5mm)のそれぞれの断面図を見たものである．L 字型アンテナの先端付近に約240000 V/m ，つまり投入電場 66000 V/m の約 3.5 倍ほどが加わっていることがわかる．
また，アンテナの軸近傍でも先端程ではないが高い電界が分布して降り，これがアンテ
ナの強電界が存在しているといえる．

33

アンテナ
(a)

(b)
	
  図 3.2 電界強度(a) x=16.5mm ,(b) z=10.5mm

34

3.3

電子エネルギー分布	  図 3.3 に時間遷移する電子エネルギー分布を示す．図 3.2(a)と同じくアンテナの中心(x=16.5 mm)における yz 面の断面図であり，それぞれの時間は(a) 0.5 ns, (b) 75 ns, (c) 125
ns である．まず，(a)はマイクロ波を投入して初期の段階であり，(a)において最も平均
エネルギーが高いのは，アンテナが作り出す電界強度の一番強いアンテナ先端付近で 12eV 程度である．電界強度においてもアンテナの先端に集中していたことより，この場所において，初期段階においては，電子エネルギーは ECR 層よりもアンテナが作り出す電界強度に依存していることが分かる．次に(b)を見ると，アンテナの先端の下方，つまりBackYork 傍に約 30 eV の高エネルギー値が見られる．ECR 層の位置がこの位置にある
ため，ECR 層における共鳴により，高エネルギーが集中している事が分かる．また，FrontYork と BackYork との間に磁気チューブが見られ，磁気ミラー効果により，電子の閉じ込めが起きていると確認できる．最後に(c)では，(a),(b)と同じく，高エネルギーの
電子は ECR 層よりも，アンテナ近傍に存在し，これよりアンテナが作り出す強電界による電子の加熱が ECR 加熱よりも支配的であることがわかる．アンテナの軸近傍において，アンテナの先端の下方と同程度の約 25 eV のエネルギー値が分布しており，アンテナの強電界による影響が考えられる．
このように時間に沿って，アンテナの電界による影響から磁気チューブの形成及びアンテナの強電界が出来る様子を見ることができ，電子のエネルギー変化を見ることがで
きる．

35

(b)

(a)

(c)

図 3.3 電子エネルギー分布(a) 0.5ns (b) 75ns (c) 125ns

36

3.4

電子エネルギー獲得確率ここで，ΔE をある領域を通過する前と後のエネルギー変化量とする．電子エネルギ
ー獲得確率とは，ある領域を通過する電子のΔE 毎の個数を数え，その各ΔE の個数を観測した回数で割ったもので，それをその領域における電子のエネルギー変化量毎の確率とした．
	  解析領域を図 3.4 に示す．例によって，アンテナ中心(x=16.5mm)における yz 面の断面図である．各々の領域はそれぞれ表 3.1, 表 3.2 に示す．
表 3.1 解析領域範囲
各軸の範囲[mm]

領域

x

Y

z

1

16.0~16.5

17.5 ~24.0

5.0~15.0

2

16.0~16.5

22.5~23.0

10.5~11.0

3

16.0~16.5

17.5~18.0

11.5~12.0

4

16.0~16.5

12.0~12.5

10.5~11.0

516.0~16.5

20.0~20.5

7.5~8.0

表 3.2 解析領域セル数
領域各軸のセル数

総セル数

x

y

z

1

1

13

20

260

2~5

1

1

1

1

従来の FDTD-PIC コードでは，領域 1 のように広い解析領域しか見ることができなかった．しかし，本研究において，領域 2〜5 のように 1 セル(0.5mm×0.5mm×0.5mm)の立法格子において解析することができるようになった．その中でも ECR 層やアンテナの位置によってどのようにエネルギーがやり取りされているかを見るために 4 つのセルに注目した．
	  図 3.5 に電子エネルギー獲得確率を示す．ヒストグラムは 0.5eV，解析時間は 125nsとした．まず，領域 1 において，これは従来のコードにおけるエネルギー獲得確率であるが，領域が広いため，ECR 層を通る数が多く，収束した確率密度関数として見ることができる．平均化すると、＋1eV 程度であり，この領域全体を見ると，エネルギー授受が行われているということができる．しかし，この領域において細かい場所ごとのエネ
ルギーの受け取りを見ることはできていない．そこで，領域 2〜5 の細かい領域において，エネルギーの授受を見ていく．
	  領域 2 において，この領域は，ECR 層から離れており，アンテナの先端付近である点である．母数が少ないので，十分なデータ数が足りていない．エネルギー的には，高エネルギーを受け取っている電子が存在しているが，平均化すると，その領域のエネルギー収支が求められ，エネルギー収支は 0 になると思われる．よって，アンテナから強い
37

電界が出ていようとも ECR 層の位置によって，その効率は変化することが分かる．
	  ここで，領域 4 はアンテナから遠い位置にあり，ECR 層に近い点である．ECR 層に近いため母数が多く，綺麗な確率密度関数となっていることが分かる．しかし，エネルギ
ー授受は平均化すると約 0eV となり、ほぼ受けとっていない．ECR 層に近くともアンテ
ナから遠ければ，エネルギーを受け取らないことが分かる．
	  領域 3 はアンテナの根本付近で ECR 層に近い領域である．比較的高いエネルギーのやり取りがこの領域において行われている．この理由として，この領域はアンテナによる強電界が存在している領域であり，この影響により，多くのエネルギーが電子に与えられている事がわかる．この領域がマイクロ波から電子へのエネルギー授受において重要な役割となっていることが明らかになった．
	
  最後に、領域 5 を見ると，この領域は，アンテナ先端下方のバックヨークに近い側でECR 層に近い点である．領域 2 と比べて，母数が多く，綺麗な確率密度関数となっていることが分かる．エネルギー値としては，10eV~20eV のような高エネルギーを電子が受けとっている．この領域はアンテナからのエネルギーを最も効率よく受け取る領域と知られており，他の領域と比べてマイクロ波の電力を効率よく受けとっていることが確認できる．すなわち ECR 加熱よりもアンテナ近傍の強電界でマイクロ波から電子へのエネ
ルギー伝達が行われていることが分かった．これは，従来と比較して細かく見ることにより，初めて分かったことである．

38

図 3.4 解析領域

図 3.5 (a)領域 139

図 3.5 (b)領域 2

図 3.5 (c)領域 340

図 3.5 (d)領域 4

図 3.5 (e)領域 5

41

第４章	
 

結論

	
  解析領域を 1 セルつまり(0.5mm×0.5mm×0.5mm)毎の領域に細かくすることによって、エネルギーのやり取りを詳細に見ることができるようになった。領域を細かくすること
により、一つの領域に入る粒子数が少なくなるので、母数が減り確率の精度が荒くなる。
故に粒子数の数と解析時間を増やす必要があり、さらなる計算コストの削減が必要であった．そこで体系中の粒子数も 10 倍の 100 万個まで増やし，それに伴う計算コストの削減のために，大型計算機の使用と MPI による並列化を行った．今回計算機として使用したのは新高性能演算サーバシステム tatara である．これにより従来の計算機より約 10倍の速度で計算することが可能となった．MPI の並列数は最適値を探した結果，tatara
上で 4 並列にすることで，約 2 倍の計算速度になった．
これにより，モデル化せずに，各セル毎のエネルギー授受が明らかになり，ECR 層以
外でもアンテナの強電界によるプラズマ生成などを調査できるようになった。
具体的に述べると，各点での平均エネルギー変化量はアンテナ近傍の ECR 層に近い点
では，3.3eV， アンテナの強電界の影響が強い点では，8.3eV であることが明らかになったように，アンテナ近傍の強電界がエネルギー授受に重要な役割を果たしていることが分かった．このことから，アンテナ形状がマイクロ波放電型中和器における性能に大きく寄与することが示唆された．

今後の課題	  今回，計算機の性能向上により，より長い解析時間まで解析することが出来るようになった．それによって，現体系においてエネルギー値が発散することが判明した．計算誤差を抑えるために更なる精度が求められる．したがって，高精度のためのさらなる計算コストの削減が必要であるし，粒子の損失機構や境界条件の精度の問題，あるいはこの FDTD−PIC コードの根本的なバグが存在している可能性もあり，コード自体の見直しが求められる．
	  また，この FDTD−PIC コードの電磁場解析の結果つまり電子エネルギー獲得確率を本研究室で開発中であるプラズマの損失まで含めた数値解析コードに導入し，より正確な
プラズマの挙動を解析出来る．今回は初期条件として電子は一様に分布しているものとしたが，このプラズマの挙動解析の結果である電子密度分布を導入して計算することによりさらなる正確性の向上を可能となる．そして，この結果を持って，本研究において開発中である壁面の損耗も考慮したコードと統合し，中和器の寿命評価を行う予定であ
る．

42

謝辞
	
  本論文の作成に当たって 1 年間研究の場を与えてくれ，貴重なご意見またご指導をして頂いた中島秀紀教授に深く感謝します．研究について様々な疑問点に根気強く答えてくださり，忙しい中夜遅くまで添削をしてくださった山本直嗣准教授に深く感謝します．
疲れている時に声をかけてくたり，研究室を明るくしてくださった大神めぐみ秘書，馬
渡隆子秘書に深く感謝します．
	
  研究室の皆様には，1 年間本当にお世話になりました．特に同じ学年として時には活発な議論をしたり，研究で辛い時は励ましてくれたりなど 1 年間苦楽を共にしてきた木邊厚視氏と湯上尚之氏には本当にお世話になりました．また，同じ研究グループとして，
計算コードに関して細かい点の疑問に答えてくださったことや，研究で行き詰まっているところを相談に乗って頂いた廣池匠哉氏と笹川裕太郎氏に深く感謝します．数値解析において貴重なノウハウを教えて頂いた中村祐輔氏と富永宙志氏に深く感謝します．夜遅くまで研究室に残った際，話し相手になっていただいた秀平靖磨氏に深く感謝します．
明るく接してくれた平野賢治氏，日永智之氏，豊田裕司氏のおかげで研究室の生活がとても楽しいものとなりました．皆様とこの研究室で共に過ごせたことを心から感謝します．1 年間研究を支えてくれた家族に深く感謝します．最後に今まで私に関わってくれた全ての人にこの場を借りて感謝の意を表します．

43