令和元年度
九州大学大学院 総合理工学府
先端エネルギー理工学専攻

修 士 論 文論文名
数値解析によるマイクロ波放電式イオンスラスタのアンテナ形状依存性に関する研究

氏 名

西牟田 直

指導教員名

山本 直嗣

目次
第 1 章 序論 ......................................................... 1
1.1.

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

1.2.

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

1.2.1.

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

1.2.2.

マイクロ波放電式⼩型イオンエンジン ......................................................................4

1.3.

過去の研究 ....................................................................................................................6

1.4.

研究の⽬的 ....................................................................................................................6

第 2 章 数値解析手法 ................................................. 82.1.

解析全体の流れ .............................................................................................................8

2.2.

プラズマ粒⼦計算⼿法(2-1)..............................................................................................92.2.1.

基礎式.......................................................................................................................92.2.2.

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

2.2.3.

粒⼦計算 PIC 法 ...................................................................................................... 11

2.2.4.

衝突計算 MCC 法.................................................................................................... 12

2.3.

FDTD-PIC-MCC........................................................................................................ 19

2.3.1.

プラズマ粒⼦計算 ................................................................................................... 19

2.3.2.

電磁波解析⼿法(2-10)(2-11)........................................................................................... 23

2.3.3.

サブグリッド法(2-16)................................................................................................. 30

2.3.4.

FDTD-PIC-MCC のフローチャート ....................................................................... 32

2.4.

計算体系および静磁場分布.......................................................................................... 34

2.5.

各種物理条件と計算条件 ............................................................................................. 37

2.5.1.

物理条件 ................................................................................................................. 37

2.5.2.

計算条件 ................................................................................................................. 38

第 3 章 結果と考察 .................................................. 41
3.1.

サブグリッド適⽤コードの妥当性検証 ........................................................................ 41

3.2.

従来コードとの電界強度の⽐較................................................................................... 42

3.3.

FDTD-PIC コード ...................................................................................................... 45

3.3.1.

電界強度分布 .......................................................................................................... 45

3.3.2.

電⼦密度分布 .......................................................................................................... 47

3.3.3.

電⼦エネルギー分布 ................................................................................................ 49

3.3.4.

EEDF...................................................................................................................... 51

3.3.5.

電⼦損失 ................................................................................................................. 52

第 4 章 結論 ........................................................ 54

3

第1章 序論
1.1.

背景2014 年に打ち上げられた宇宙探査機“はやぶさ 2”は目標である小惑星“リュウグウ”への到着が成功し，岩石のサンプルを地球に持ち帰るべく，現在期間運用中
である．地球から“リュウグウ”までは約 3 億万 km もの距離があり，このミッシ
ョンに要求されるのは長時間の安定動作かつ高い速度増分である．この要求を満たすべく“はやぶさ 2”に搭載されているのが，マイクロ波放電式イオンエンジン“μ10”である．
イオンエンジンは上記のような深宇宙探査以外にも人工衛星のエンジンとしても応用されている．1957 年に旧ソ連により世界初の人工衛星である“スプートニ
ク 1 号”が打ち上げられたが，現在では GPS 機能や天気予報など身近なサービスに人工衛星が活用されている．これらの人工衛星が小型化が進んでいる．小型衛星
においてはサイズの制約から高燃費の推進器開発が必須となる．この高燃費の推進を実現する推進器として利用されているのが電気推進である．
電気推進(1-1)とは，電気エネルギーを推進剤に与えることで反力を得る推進機であり，従来の化学反応を利用したロケット推進と違って静電力や電磁力を利用できる．このため比推力（推進剤流量あたりの推力）を一桁以上大きくすることが可能であり，大幅な推進剤の低減が望める．現在使用されている代表的な電気推進機はアークジェットスラスタやホールスラスタ，そしてイオンエンジンである．イオン
エンジンは他の電気推進機と比較しても比推力が高く，長期間の人工衛星の姿勢制
御や惑星探査に適している．

1.2.

イオンエンジンの原理

イオンエンジンの概念図を図 1.1 に示す．イオンエンジンはアーク放電やマイク
ロ波などで推進剤を加熱・電離させてプラズマを生成し，2 枚ないし 3 枚からなる
多孔状の電極に 1000 V 程度の高電圧を印加させてイオンを加速するという静電加速型の推進装置である．イオンエンジンは主に 3 つの領域から構成されている．

① 推進剤を電離するイオン生成部（Ionization）② 生成されたイオンを静電的に加速して推力を得る加速部③ 放出されたイオンビームを電気的に中和する中和部
1（Acceleration）（Neutralization）

これらの各過程はそれぞれイオン源，加速電極，中和器によって行われる．

図 1-1

1.2.1.

イオンエンジンの概念図

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

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

①プラズマ生成に電極を必要とせず，それによるエンジン寿命の制約がない．
②ヒーターを用いないため，予備加熱が不要となりエンジンを即座に立ち上げることが出来る．
③1 台のマイクロ波電源で複数台のエンジンにプラズマを生成できるため，中和器

2

のプラズマも同時に生成することが可能となり，システム構成の簡素化・信頼性の向上・衛星の重量軽減をもたらす．

図 1.2

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

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

mdv= -ev ´ B
dt

(1.1)

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

rL =

mv ^
eB

(1.2)

3

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

w ce =

eB

(1.3)

m

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

E

B

電子

E

B

電子

電子
E

B
B

電子

E

図 1.3 ECR の原理

1.2.2.

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

本研究室では HAYABUSA(1-2)に搭載されたマイクロ波放電式中和器を元にマイクロ波放電式小型イオンエンジンを作成している．図 1.4 にスラスタの写真を，図1.5 にその構成図を示す．Al 製円筒放電室の周囲に 4×4×12 mm の長手方向磁化4磁石を複数設置し，それらを軟鉄製ヨークで挟み込んで放電室内に磁気回路を形成させている．磁石は永久磁石の Sm-Co (サマリウムコバルト) を用い，アンテナ・
グリッドはモリブデン (Mo) 製である．

図 1.4
Al

マイクロ波放電式小型イオンエンジンの写真
Mo

Soft iron

Sm-Co magnet

Propellant gas
Antenna

Front yoke

Microwave

Housing
Back yoke
Discharge chamber
Screen grid

図1.5

Acceleration grid

マイクロ波放電式小型イオンエンジンの構成図

5

1.3.

過去の研究

これまで本研究室でマイクロ波放電式小型イオンエンジンの磁場強度・アンテナ形状・放電室形状などを変化させて実験を行ってきた(1-3),(1-4),(1-5)．しかし，形状を変更させたことによって，放電室内部の物理メカニズムがどのように変化し，どのような理由で性能に影響をおよぼしているのかは解明されていない．その原因はエンジン自体が小型であるために静電プローブ等の装置を入れることで放電室内部のプラズマに擾乱を引き起こし，正確な測定ができないからである．
そこで数値解析は非常に有効な手段である．近年の計算機の発展のおかげで実験では知り得る事のできない情報を入手できるようになってきた．それにより，実験による時間的・金銭的コストを抑えた研究開発ができるようになってきている．本研究室では増井らによって，3 次元数値解析コードが開発されてきた(1-6)．
このコードは電子とマイクロ波との相互作用の解析を目的としているが，計算コストの兼ね合いからメッシュサイズが粗いため，アンテナ形状の尖りや丸みを綺麗に再現することが出来ない．また，メッシュサイズが粗いことにより，アンテナの表皮深さが過大評価されているとの報告もされている(1-7)．

1.4.

研究の目的

以上の背景より，本研究の目的は３次元数値解析コードに局所的にメッシュサイ
ズが細かい領域を設け，アンテナ形状による放電室内部への影響を調査することで
ある．

参考文献
(1-1) 栗木恭一・荒川義博[編]：”電気推進ロケット入門”, (東京大学出版会,2003)(1-2) I. Funaki, H. Kuninaka and K. Toki, J. Propul. Power, Vol. 20, (2004), pp.718-726.
(1-3) N. Yamamoto, H. Kataharada, H. Masui, H. Ijiri and H. Nakashima, AJCPP200522093, 2005.
(1-4) N. Yamamoto, H. Kataharada, T. Chikaoka, H.Masui and H. Nakashima, IEPC-2005036, 2005.
(1-5) T. Chikaoka, S. Kondo, N. Yamamoto, H. Nakashima and Y. Takao, Proceedings of the25th International Symposium on Space Technology and Science, 2006, pp.254-259.
6

(1-6) 増井博一， “マイクロ波放電型宇宙推進用プラズマ源に関する数値解析”, 平
成 17 年度九州大学博士論文．
(1-7) 牛尾康一， “マイクロ波放電型プラズマスラスタに関する研究”, 平成 30 年度九州大学博士論文．

7

第2章
2.1.

数値解析手法
解析全体の流れ

数値解析手法について, 述べていく.本研究室では表 2-1 に示す 2 つの数値解析
コードを組み合わせて計算を行うことを目的としてコード開発が進められている．
表 2-1 数値解析コード
物理対象

モデル名

電子とマイクロ波の相互作用による電子のエネルギー獲得解

FDTD-PIC-MCC

析Full−PIC−MCC

イオンと電子の挙動，静電場解析

これら 2 つの計算モデルを弱連成的に解析させることで，プラズマ挙動の解明および寿命の評価ができる．図 2-1 にこれらのコードの連成手法を示す．まず，前章で述べたようにマイクロ波放電式イオンエンジンの基本原理は永久磁石とマイク
ロ波の ECR 共鳴であるため，静磁場の情報が必要となる．スラスタに永久磁石を配置した際の静磁場分布は，解析を始める前に市販ソフトで計算しておく(後述)．
この静磁場データを持って，以降の FDTD-PIC-MCC および Full−PIC-MCC の計算が行われる． FDTD-PIC-MCC では，一様状態から電子とマイクロ波の吸収が計算され，電子エネルギー獲得確率が算出される．これは，単位時間単位体積あた
りの電子のエネルギー増分を確率密度分布として算出後，電子のマイクロ波吸収量の効果として Full−PIC-MCC に渡す．Full−PIC-MCC では，イオンと電子を一様状態から計算し，静電場の効果による粒子の輸送を解析する．プラズマ密度が変化するとカットオフ密度(電磁場がほとんど浸透不可能である密度)以上の分布が出現してきて，密度分布が変化するため，電子のマイクロ波吸収量が変化する．実際に，Full-PIC-MCC の計算において，初期電子エネルギー獲得確率でエネルギーを与え続けると，特定の領域のみでマイクロ波から電子へのエネルギー伝達が行われ続けた結果として電子温度が上昇し，計算が発散してしまう．このため，プラズマ密度の変化を FDTD-PIC-MCC に再びフィードバックすることで，カットオフの効果を取り入れ, 電子密度の空間分布によるマイクロ波の伝搬の変化を考慮している．また，エネルギー保存のために，Full−PIC-MCC 計算での電子エネルギーも同様にフィードバックを行う．Full−PIC-MCC の計算において，逐次プラズマ密度8の変化とともに，FDTD-PIC-MCC にフィードバックしながら繰り返し計算し，定常状態を目指す．このように，時間スケールの異なる現象過程をそれぞれ独立させて計算させ，それぞれの結果を受け渡しながら計算する事を弱連成解析という．現段階では Full-PIC-MCC コードが未完成なため，3D-FDTD-PIC-MCC のみを用
いて計算を行い，3 つのアンテナ形状を用いて電子加熱特性の解析を行った．

図 2-1 解析全体フロー

2.2.

プラズマ粒子計算手法(2-1)

Full-PIC-MCC も FDTD-PIC-MCC でもプラズマ粒子計算はモデルとしては．
同等のものを利用している．但し，Full−PIC-MCC は円筒座標，FDTD-PIC-MCC
は直交座標を用いている．そのため差分式が異なる．そのため，共通項のみをこの節で述べ，詳しい解法については，各モデルの節で述べる．

2.2.1.

基礎式

プラズマの挙動を解析する際の基礎式として運動方程式を用いる．式を以下に示す．
mdv= q(E + v ´ B )
dt

9

(2.1)

dx
=v
dt

(2.2)

ここで， m は粒子質量， v は粒子の速度， x は粒子の位置，t は時間， q は粒子の持
つ電荷， E は電場， B は磁束密度を表わす．

2.2.2.

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

粒子法におけるメッシュサイズは， CFD のような流体計算とはことなり格子形
成法自体が解析の精度に影響をあたえることはなく，対象とする物理現象を再現することができるかどうか次第である．まず，解析したい物理現象に対するメッシュ
サイズを以下のようなデバイ長などから決定し，メッシュサイズからくるスキームの制約からタイムステップが決定されることになる．
計算に用いるメッシュサイズ D は，電子温度 Te，ボルツマン定数 kb，電子密度

ne から計算される．
デバイ長 λDe，lDe =e 0 k bTene e 2

(2.3)

に対して，

D £ lDe(2.4)

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

ne e 2w pe =
me e 0(2.5)

Dt £ wpe-1

(2.6)

に対して，

とすれば，数値的な不安定性は生じない．時間積分の精度まで考えると，通常はDt £ 0.2wpe-1(2.7)

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

10

wc =

eB
me

(2.8)

に対して，

Dt << wc

-1

(2.9)

である必要がある．
また，粒子は時間刻みごとに１セル以上飛び越えてはならないという制約がある．
特に荷電粒子は格子ごとに規定される電磁力が働くため，この条件が必須となる．
これを Courant の安定化条件という．3 次元の際，時間ステップ Dt は，
1

Δt ≤
v

1
2

(Δx )

+

1
2

(Δy )

+

1

(Δz )2(2.10)

を満たさなければならない．この条件は極めて厳しく，わずかでも満足しなければ計算が不安定になる．実際には，この 10 分の 1 程度が望ましいとされている．ここにおいて，𝑣は，粒子の最大熱速度または，電磁波速度つまり光速である．FDTD-PIC と Full-PIC において，条件が異なるので各モデルの計算条件において説明する．

2.2.3.
2.2.3.1.

粒子計算 PIC 法
超粒子(2-3)粒子法では，粒子ごとに基礎式を記述し，解くものであるが，マイクロ波放電式
イオンエンジンでは，1 m3 当りに 1017～1018 個が含まれ，計算機全体でも 1012 個あり，この粒子全部を取り扱うことは最新のスーパーコンピュータを用いても不可能であり，実用的でない．そこで粒子シミュレーションで取り扱う粒子として，実際の粒子の電荷質量比を一定に保ったまま多数の粒子の電荷と質量を１つにまとめた超粒子を使用する．この超粒子数に含まれる粒子数が少なければ，少ないほど，ノイズが小さくなるが，計算コストも高くなるので，そのトレードオフになる．超粒子数は，１セルに含まれる粒子数が 10 個以上になる事と取り扱う粒子数よりも超粒子の重みが少ない事（統計性から）から決定する．

11

2.2.3.2.

重み付け

シミュレーションでは，電場量を空間格子点上で求める.そのため,各粒子の電荷を格子点上に振り分ける必要がある.また,粒子の運動方程式を解く際に,電場解析
によって格子点状に割り振られた値を内挿する必要がある.このような計算を行うことを PIC 法（Particle In Cell method）と呼び，粒子シミュレーションにおいて標準的な手法である．内挿の概念図は，各コードで異なるため，後で各モデルにつて述べる．

2.2.4.

衝突計算 MCC 法

粒子コードにて，衝突を取り扱う場合，粒子の衝突の解析には，MCC(MonteCarlo-Collision)法を用いるのが，最速の手法である．粒子同士の衝突判定を各粒子感の衝突周波数から決定する確率と乱数を用いて判定するというのがモンテカル
ロ衝突法である．以下に詳しく述べる．

2.2.4.1.

衝突種別

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

2.2.4.2.

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

面 a ，厚さ l ，その物質の数密度 N の板を考える．エネルギー E の粒子ビームが，面 a に垂直に当るものとする．衝突をせずに板を透過してくる粒子数は，入射粒子数より減少する.その減少度は， s を衝突の断面積とすれば，dn = -nNsdl(2.11)

n = n0 exp(- Nsl )(2.12)と表わされ，これよりとなる．ここで，n0 は面 a（ l = 0 ）でのビーム中の粒子である．それゆえ，（ l ,l + dl ）間で最初の衝突が起こる確率は，
Pinitial (l )dl = exp(- Nsl )Nsdl

12(2.13)

となり， l 以下で最初に衝突が起こる確率は，P(l ) = ò Pinitial (l )dll0= 1 - exp(- Nsl )(2.14)となる．最初の衝突までに粒子が走る平均距離 l は P(l )の一次のモーメントである．
¥

l = ò lP(l )dl
0

=

1
Ns

(2.15)

これを平均自由行程という．よってモンテカルロ法によって，任意の点から出発し最初の衝突が起こるまでの距離 l を決めるには，乱数 r を用いて，
r = P (l )æ lö
= 1 - expç - ÷è lø(2.16)

l = -lln(1 - r )(2.17)

で表わされ，これを変形して，となる．ここで， r は一様乱数であるから 1 - r を r とすれば，

l = -lln(r )(2.18)と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を s とすれば，この粒子が１秒間に気体分子と衝突する回数 u は，

u=

v

l
= sNv

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

衝突確率・衝突判定

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

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

u total = ngass total v m

13(2.21)

Pcollision, m = 1 - exp(- ngass total v m Dt )= 1 - exp(- u total Dt )(2.22)

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

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

* 0 £ r2 £s elastics total

であれば弾性衝突*

s elastics+ s excitation£ r2 £ elastics totals totalであれば励起衝突*s elastic + s excitations+ s excitation + s ionization£ r2 £ elastics totals totalであれば電離衝突

ここで，σ total = σ elastic + σ excitation + σ ionization である．このように断面積の比と乱数を対応させ，衝突の種類を決める．計算では中性粒子としてキセノンを用い，その断面積
データ(2-5)を図 2-2 に示す．今回の計算では，多価電離は考慮していない．

図 2-2 電子と Xe の衝突断面積14

2.2.4.5.

エネルギーと散乱角

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

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

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

(2.25)

h の直交座標成分は
hx = ver cos f

hy = -

vex vey cos f - ve vez sin fver(2.26)(2.27)v v cos ϕ + ve vey sin ϕhz = − ex ezver(2.28)

ver = vey2 + vez2

(2.29)

ve = vex2 + vey2 + vez2(2.30)

ここで，であり， c は散乱角， f は任意の角である．(2.25)式はエネルギー損失を含んでいる．

15

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

Escattered, e = Eincident, e(2.31)

電子は，エネルギー E incident, e が小さいとき等方散乱され，大きいときは主に前方散乱となる．これを考慮した散乱角 c の確率密度を g(c )として以下に示す(2-14)．ここで， E incident, e の単位を eV とする．
g (c ) =Eincident ,e sin c(0 £ c £ p )cìü2í1 + Eincident ,e sin 2 ln 1 + Eincident ,e ý2îþ(2.32)

これを解けば，{

2 1 - (1 + Eincident, e ) 3
cos c = 1 +Eincident, er}(2.33)となる．この値を(2.25)式に用いて衝突後の速度を求める．このとき，任意の角は
f = 2p r4(2.34)

として与える．
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エネルギーを Eexcitation とし，入射電子のエネルギーを E incident, e ，励起後のエネルギーを Escattered,e ，励起後の速度を v~ とする．エネルギー収支は，Escattered, e = Eincident, e - E excitation(2.35)と表わされ，励起後の速度はE
v~ = v e 1 - excitation
Eincident, e(2.36)となる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネルギーを失った後に(2.36)式中 v~ の速度で弾性衝突をすると考えることができる．(2.25)式
中 v e に v~ を用いて，衝突後の速度 v e¢ を求める．このときの散乱角には(2.33)式を用いるが(2.33)式中の E incident, e には(2.35)式から求められる Escattered,e を用いる．励起エネルギー E excitation は，キセノンの場合 8.34 eV を用いている．
③電離衝突
入射電子のエネルギーを E incident, e ，散乱電子のエネルギーを Escattered,e ，生成電子のエ16ネルギーを Ecreated,e ，電離エネルギーを Eionization とする．エネルギー収支はEscattered, e + E created,e = Eincident, e - Eionization(2.37)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる．

ìæ Eincident,e - Eionization öü÷÷ýEscattered,e = B tanír5 tan -1 çç
2B
è
øþ
î

(2.38)

ここで，B はガスに固有の値でキセノンでは 8.7 eV となる．
また，入射電子が電離によって失うエネルギーを DE とすると，
DE = Eincident, e - Escattered, e(2.39)と表わされる．(2.36)式から電離後の電子の速度は，

v~ = v e 1 -DEEincident,e(2.40)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．(2.27)式
中 v e に v~ を用いて，衝突後の速度 v e¢ を求める．散乱角には(2.33)式を用いるが，(2.33)
式中の E incident, e には(2.38)式から求められる Escattered,e を用いる．電離エネルギー Eionization
は，キセノンの場合 12.13 eV である．
電離後の生成電子の持つエネルギーは，E created,e = DE - Eionization(2.41)となる．また，入射電子の速度 v e を用いて生成電子の速度はv created,e =veve

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

Null-collision 法

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

s fake を導入する．すべての電子のエネルギー E に対して
s fake < s total(2.43)であるとし，例えば電子と中性粒子の衝突においてはs total = s elastic + s excitation + s ionization + s fakeとする．ここで，17(2.44)

s totalv =ucollisionn gas= constant(2.45)

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

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

2.2.4.7.

タイムステップの制約

モンテカルロ衝突法では，タイムステップの間にそれぞれの粒子についての衝突
を 1 回ずつしか処理できない．しかし，有限のタイムステップを用いる限りは，1
タイムステップの間に複数回の衝突が含まれうる．複数回の衝突による衝突のミスを減らすためにはタイムステップを十分小さくしなければらない．

18

ある粒子のタイムステップあたりの衝突ミス回数は∞P i2
r ~ ∑P =1− P ik=2k
i(2.46)であるので，衝突ミスの回数を１％いかにするには𝑃# ≦ 0.1とする必要がある．
Null-collision を用いる場合は，𝑃()** ≦ 0.1 すなわち

Δt ≤2.2.4.8.

0.1
ν null(2.47)

乱数生成

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

2.3.

FDTD-PIC-MCC

FDTD-PIC-MCC は，電子のみを粒子として扱い，イオンと中性粒子は背景としておいておき，電子の時間スケールにて３次元的に解く．後述する FDTD(Finite
Difference Time Domain)法によりマイクロ波解析を行い，電場 E の効果として，
PIC 法にて内挿し，電子の電磁運動方程式を解くものである．さらに， 前述の衝突計算を MCC 法にて実施したものである．以下に詳しく述べる．

2.3.1.

プラズマ粒子計算

2.3.1.1.

運動方程式解法

基礎式は，2.2.1 で述べた(2.1)式に leap-frog 法による時間中心差分を施すと，
x n +1 - x n= v n +1 2
Dt

(2.48)となる．ここで，各変数の上付き添え字は時間に関するもので， 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 ý
Dt
mî
2
þ

(

)

(2.49)

である． n-1 2 はある時刻 n から Dt 2 だけ以前の値であることを示す．同様に，

19

n + 1 2 はある時刻 n から Dt 2 だけ進んだ値であることを示す．
プラズマ粒子の運動を求めるために， (2.49)式において未知量である v n +1 2 について解き，求めた v n +1 2 と(2.48)

式より新たな粒子の位置 x n + 1 を求めることができる． (2.49)式を変形すると，

Av n +1 2 = S(2.50)と書ける．ここで， α ・ b ・ g ・ d を次のように定義する．

a=

qDt n
Bx
2m

(2.51)

b=

qDt n
By
2m

(2.52)

g=

qDt n
Bz
2m

(2.53)

d = 1+a 2 + b 2 + g 2(2.54)

なお，ここでの E ， B は粒子位置におけるそれぞれの値を表わしている．これらを用いて行列 A とベクトル S は，以下のように構成される．
æ 1
ç
A=ç g
ç- b
è

-g
1

a

b ö
÷
-a ÷
1 ÷ø

æ v xm + cE xn + gv ym - bvzm öç÷
S = ç v ym + cE yn + avzm - gv xm ÷ç v m + cE n + bv m - av m ÷
z
x
y ø
è z

(2.55)(2.56)

ここで，c = qDt m ，m = n - 1 2 である．さらに (2.50)式を次のようにして変形し，速度を求めることができる．

2.3.1.2.

v n +1 2 = A −1 S

(2.57)

æ 1 + a 2 ab + g ag - b ö÷1ç
A -1 = ç ab - g 1 + b 2 bg + a ÷dç
2 ÷èag + b bg - a 1 + g ø(2.58)

PIC 法による重み付け計算

3 次元における重みつけ方法を示す．粒子の位置から近傍のセル座標 (i , j , k )を求める．次に，図 2-4 に示す dx1, dx 2 , dy1, dy 2 , dz1, dz 2 を求める．その値を用いて，

20

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

Vn
V

(n = 1  8)

(2.59)となり，

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)と表わされる．

図 2-4 PIC 法におけるセルの構成

21(2.60)

2.3.1.3.

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

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

ρ (i , j , k ) =

V1
qN
V2

(2.61)となる．すべての粒子について和をとれば，格子点状での電荷密度が得られる．
また，電流密度は

j = qnv(2.62)

より，
j (i , j , k ) =2.3.1.4.

V1
qNv
V2

(2.63)

粒子の壁面境界条件

FDTD-PIC-MCC では，静電場の効果を考慮していないため，電子のシース電位による反射を仮定する必要がある．シース電位については，Full−PIC-MCC にて解析可能であるが，現在電位分布の影響を連成することはできていない．そのため，実験値と過去の数値解析の結果から 20V と設定し,壁面に達した場合，一律で20eV 以上は消失し，それ以上は反射するものした．また，Plume 領域 引き出し電流として，排出される出口に達した場合，無条件で消失するものとした．また,このシース電位の効果を弱錬成するときに渡すべきであるが，シース電位の判定が難しいため, 連成は行っていない.

2.3.1.5.

電子のエネルギー獲得確率計算方法

Full−PIC-MCC との連成で重要となる電子エネルギー獲得確率の計算方法を示す．1 つの観測時間をマイクロ波の１周期とし，タイムステップの間つまり運動方
程式を解いた前後で変化したエネルギー増分をカウントし，エネルギー増分のヒス
トグラムとして，確率密度分布とする．また，各エネルギー増分はセル毎で観測されており，粒子の位置情報をセルに割り当てている．

22

2.3.2.
2.3.2.1.

電磁波解析手法(2-10)(2-11)
基礎式

FDTD 法は数値電磁場解析の分野において，モーメント法や有限要素法と同様に，アンテナや平面回路など数多くの複雑な問題に適用され，主要な解析手法として多く用いられている．計算には，マクスウェル方程式を基礎式として用いる．以
下に基礎式として用いる微分形式のマクスウェル方程式を示す．

¶H¶t(2.64)¶E+J¶t(2.65)

Ñ ´ E = - µ0
Ñ ´ H = e0

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

差分式

三次元の FDTD 法で用いられる電場，磁場の空間配置を図 2-5 に示す．このように，電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．これは Yee によって提案され(2-13)，電場の回転（ Ñ ´ E ）が磁場を，磁場の回転（ Ñ ´ H ）が電場の方程式を満たすように配置されている．
この電磁場の配置のもと， (2.64)式, (2.65)式を空間，時間とも中心差分した式を以下に示す．((2.66)~(2.71))式ここで，Dx, Dy, Dz は空間幅， Dt は時間幅，添え字(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
&%
)%

23

(2.66)(2.67)

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 ) +Dt
e 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 )
í
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 )E yn +1 (i , j + 1 2 , k ) = E yn (i , j + 1 2 , k ) +(2.68)

Dz(2.69)üï- J x (i + 1 2 , j , k )ýï
þ

Dt
e 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)
í
Dz
ï
î

(2.70)

11n+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 )ýDxïþE zn +1 (i , j , k + 1 2) = E zn (i , j , k + 1 2) +Dt
e 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
í
Dx
ï
î

(2.71)

11n+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)ý
Dy
ï
þ
1

n+

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

n+

2

24

は E n から計算されるので，アンテナの図 2-5 FDTD 法における電場と磁場の配置

2.3.2.3.

アンテナの模擬および給電方法

本計算ではアンテナを完全導体で構成されているものとしてモデル化する．給電方法は同軸ケーブルを模擬し,図 2-6 にセントラルヨーク内部の給電面での断面を示す．ただし，メッシュサイズは 0.5 mm としている（局所グリッドについては後述）．アンテナは直径 1 mm，セントラルヨーク内径は 2mm とし，図中の灰色のセルが完全導体，白色のセルが電磁波が伝播する自由空間である．また，アンテナ導入部には図 2-7 のように Mur の一次吸収境界条件を与えている．自由空間に投入する電界は以下の式から決定した．
𝐸- =

ただし,

𝑉0
𝑟

(2.72)

𝐸2 = 𝐸- cosθ

(2.73)

𝐸7 = 𝐸- sinθ(2.74)

V は中心導体と外部導体間の電圧, 𝑟 はアンテナの中心軸からの距離，はアンテナ中心軸を原点として各給電電界と x 軸がなす各出である．V0 はセント25ラルヨーク内を通過するポインティングベクトルをもとに決定した．

図 2-6 同軸線路給電

図 2-7 アンテナモデル

26

2.3.2.4.

PEC(Perfect Electric Conductor)境界条件

境界面が金属のような完全導体である場合，境界に平行な電場成分を 0 とする
PEC 境界条件を用いる．放電室内壁は導体を仮定しているので，この境界条件を適用する．

2.3.2.5.

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

開放領域を取り扱う場合，計算領域を仮想的な境界で閉じなければならない．この境界を吸収境界といい，この吸収境界が完全でなければ，反射波が解析領域に戻ってしまい，誤差の原因となる．ここでは，定式化が容易で広く使用されている
Mur の吸収境界条件を用いる．

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

Ez = Ez (x + vt )(2.75)と表わされ，次の微分方程式を満足する．

¶E z 1 ¶E z
=0¶x v ¶t(2.76)

x = 0 の境界で完全に反射がないならば，電場は(2.75)式の形を保ったまま伝播するはずである．すなわち，境界でも(2.76)式が満足される．これを FDTD 法の形
式に定式化するため，時間について差分化すると∂E z E zn − E zn −1
=
∂t
Δt
n−

∂E
=v z
∂x

(2.77)

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

E (i, n - 1 2) - E
Dt
n
r

n-1
r

n-

1
2

n-

1
2

(i, n - 1 2) = v Er (i, n) - Er (i, n - 1)Dz

(2.78)となる. (2.77)式の Ezn (1 2), Ezn -1 (1 2)も FDTD 法では直接求めることはできないた
め，前後の値で平均して求める．そして， E z についてまとめると

27

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.79)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．
2.3.2.6.

投入エネルギーの計測

空間を流れる電力を求めるためには，ポインティングベクトルを使用する．ポイ
ンティングベクトルは以下の式で表され，電磁波のエネルギーの流れを表す．
𝑺 =𝑬×𝑯(2.80)

電磁波の基本的な性質として，伝播経路のインピーダンスが変化すると，その地
点で波の一部が反射してしまう．この波を反射波といい，スラスタにマイクロ波を投入した際にも投入したうちの一部がアンテナから放射され電子に吸収されるが，一部は反射波として返ってくる．従ってある一面のポインティングベクトルを計算しただけでは，反射波を差し引いた正味の投入電力しか計算できない s．しかし，実際の実験では投入は電力が制御パラメータであり，数値計算でもこの値を算出する必要がある．そこで図 2-8 に示すように，二つの面でポインティングベクトルを計算した．緑色の線で表した面に電場を投入すると，その両側に同じエネルギーを持った投入波が伝播していく．アンテナ先端側に伝播した波は放電室内へと放射され電子に吸収されるが一部は反射波として返ってくる．そして，反射波はアンテナ根元まで到達し，先述の通り Mur の一時吸収境界条件で吸収される．そのため，電場投入地点からアンテナ側では，投入波(FWD)と反射波(REF)の進行方向が逆で，ポインティングベクトルはこれらの波が打ち消し合った値 S1 となる．一方，電場投入地点よりもバックヨーク側では投入波と反射波の向きがそろっており，ポインティングベクトルはこれらが強め合った値 S2となる．そのため，投入波のポインティングベクトル Sfwd と反射波のポインティン
グベクトル Sref は以下のように表される．
𝑆FGH =𝑆I + 𝑆K
228(2.81)

𝑆-LF =−𝑆I + 𝑆K
2

(2.82)

この Sfwd の値を投入電力として扱い、この値が目標値となるようにアンテナ-セン
トラルヨーク間の電圧を決定した。そして、FDTD-PIC コードの計算の間、その電圧は一定とした。

図 2-8 マイクロ波給電の計算

2.3.2.7.

投入電圧の補正(2-15)

投入マイクロ波電力の制御因子は, 先に示したアンテナ根本に励振する電場強度である.時間発展するにつれ,電子 の状態は変化するため,アンテナから見た負荷
インピーダンスは変化し,一定の電場を根本に与えた場合,投入 されるマイクロ波
電力も時間とともに変化する.本計算 では,一定のマイクロ波電力を放電室内に投入するために,以下の式にしたがって,根本に与える電場強度をマイクロ波 1 周期毎に更新した.
𝑉(LG = 𝑉N*H +𝑃QR-SLQ𝜃P 𝑃− 1T 𝑉N*H
N*H(2.83)

2

ここで，𝑉N*H , 𝑉(LG は，更新前後の励振電場強度，𝑃N*H は更新前のマイクロ波 1 周期で平均した投入マイクロ波電力，𝑃QR-SLQ は設定投入マイクロ波電力，θは補正係数
29

とする．

2.3.3.

サブグリッド法(2-16)アンテナ形状をより明細に表現するためにはメッシュサイズをより細かくする必要がある．しかし，計算領域を一様に微細セルを用いてモデル化するとセル数が急激に増加し，膨大な計算コストが掛かる．今回，微細セルを用いる必要があるのはアンテナ周辺の空間のみであるため，局所的に微細構造を扱うことができる手法であるサブグリッド法を用いた．サブグリッド法を選択した理由は，解析対象物ごとに局所セルの形状，配置を変える必要がなく，異なるアンテナ形状に同一のコードで汎用的に使うことができるためである．
サブグリッド法とは，Yee セルの中にさらに細かな Yee セルを作り，微細構造をモデル化する方法である．周囲の大きなセル全体を主グリッド，主グリッド内の細かなセル全体を局所グリッド，それらの境界を MG-LG 境界という．
今回は電界を含むように MG-LG 境界をとった．例としてセルサイズ比を 1:3 としたときの主グリッドと局所グリッド内の電磁界の配置を図 2-9 に示す．セルサ
イズ比を奇数字にすることで主グリッドの電磁界 E, H と局所グリッドの電磁界 e,h が一致する場所を設けている．局所グリッドの電界 e はこれを囲む 4 つの主グリッドの電界 E を用いて空間的に補間している．また，時間的な補間に関しては局
所グリッドに属する電磁界の計算回数をセルサイズ比回行うことで補間している．
サブグリッド法のフローチャートを図 2-10 に示す．まず局所グリッド内も含めた主グリッドに属する電磁界を計算した後，境界面において局所グリッドの電界を主グリッドの電界から空間的に補間する．補間した電界から時間変化における局所グリッドの電磁界を解き，境界近傍の電界及び境界面の磁界を各グリッドで平均化する．これは境界面での発散を防ぐためである．これらの計算を繰り返すことで主グリッドと局所グリッドを補間している．
セルサイズ比については先述の理由により奇数字，また将来的に局所グリッド内にさらに局所グリッドを作りメッシュサイズをデバイ長(今回の計算条件では 20µm 程度)に近づける予定であるため，今回はセルサイズ比を 1:5 として計算を行った．

30

図 2-9 主グリッドと局所グリッド内の電磁界31図 2-10 サブグリッド法のフローチャート

2.3.4.

FDTD-PIC-MCC のフローチャート

図 2-11 に FDTD-PIC-MCC のフローチャートを示す．静磁場とプラズマ密度は事前に計算したものを使い，初期化処理の際，粒子の配置を密度分布に従って作成する．FDTD 部でマイクロ波の供給を行いながら，シームレスにマクスウェル方程式を解き，電磁場を求める．その後，PIC のタイムステップのタイミングになった
時に，PIC 部で運動方程式を解いた後，衝突計算を実施する．粒子の移動の更新後，
32

格子点に電流密度を割り当て，それを FDTD 計算に渡す．最終的に繰り返し回数が最大になれば計算を終了する．

図 2-11 FDTD-PIC-MCC のフローチャート

33

2.4.

計算体系および静磁場分布

図 2-12 に本研究室所有のマイクロ波放電式イオンスラスタと解析領域を示す．
これは yz 断面で切り取った断面図である．形状としては，円筒形の放電室を棒状の永久磁石（サマリウムコバルト磁石）10 個が周方向に均等に配置されている．
磁石を鉄でできたフロントヨークおよびバックヨークが挟んでおり，この間に磁力線を形成する形状となっている．図 2-13 に計算領域を示す．計算領域には主グリッドと局所グリッドを配置しており，局所グリッドの電磁波の計算にサブグリッド法を用いている．アンテナを囲う領域を局所グリッドとし，x:6-16 mm, y:6-16 mm,
z:9-10mm と設定した．図 2-14 各グリッドで描画されたアンテナ形状を示す．左側がメッシュサイズ 0.5 mm で描画されたアンテナ，右側が 0.25 mm で描画されたアンテナである．メッシュサイズを細かくすることで形状再現性が向上し，アン
テナ先端の尖り，丸みがより際立っている．
この放電室内部の静磁場強度分布 x 中心断面図を図 2-15 に示す．これは，Advanced Science Laboratory, Inc 社 の 3 次元磁場解析ソフト Amaze を用いて解析した磁場強度分布である．今回使用するマイクロ波周波数は 2.45GHz であるので，このときの ECR 層の磁場の強度は 87.5 mT である．
12

3

2

PEC
10.0

Antenna
1

21

9

Front yoke
Back yokeyxzSoft ironAl

Molybdenum

MagnetCalc Region

図 2-12 イオンスラスタの解析領域

34

図 2-13 計算領域(a) 円盤型

35

(b) 十字型(c) 星型
図 2-14 各メッシュサイズにおけるアンテナ形状36

図 2-15 放電室内の静磁場分布

2.5.

各種物理条件と計算条件

2.5.1.

物理条件

表 2-2 物理条件に本計算に用いた物理条件を示す． 本研究室の実験体系と同等のものを扱っている．推進剤流量は計算上では，一定の中性粒子密度として扱われている．これは，別途 DSMC 計算にて，0.2 mg/s を本スラスタに流した時，中性粒子密度がガス口を除き，スラスタ全体で約 4.3×1018 m-3 という結果が得られているため，衝突計算における中性粒子密度として，この値で一定とした．イオン初期エネルギーは室温として，初期密度初期電子エネルギーに関しては，過去の実験データより設定した．(2-17)37

表 2-2 物理条件

2.5.2.

推進剤

Xe

推進剤流量

0.2 mg/s

マイクロ波周波数

2.45 GHz

投入電力

8W

中性粒子初期密度4.3×1018 m-3

電子初期エネルギー

2.0 eV

計算条件

主グリッドのメッシュサイズを 0.5 mm とし，x 方向に 44 格子 y 方向に 44 格子 z
方向に 32 格子とした．局所グリッドはアンテナを含むように配置し，メッシュサ
イズを 0.1 mm とし，x 方向に 100 格子 y 方向に 100 格子 z 方向に 10 格子とした．主グリッドのメッシュサイズ 0.5 mm より FDTD における計算制約条件は，
XI
プラズマ周波数より， 0.2 𝜔WL= 8.0 ×10-12 sec , 電磁波の伝搬速度を光速と仮定

すると Courant 条件より，9.6 ×10-13 sec 以下であればよく計算誤差を考え，2.5
×10-13 sec とした．PIC における制約条件としては，サイクロトロン周波数より，2.0 ×10-12 sec , 衝突計算の制約より，7.3×10-10 sec，最大電子速度を 50eV とした時の 1 セルあたりを飛び越えないための Courant 条件は，6.9×10-11 sec となり，さらにその 10 分の 1 程度とすると 7×10-12 sec となる．このなかで最も制約条件の厳しいのは，Courant 条件であるためこれを満たしさえすればよい．
FDTD10 回につき PIC を 1 回解くとすると，2.5×10-12 sec で十分計算の制約を満たすが，このタイムステップで計算を行うと，アンテナ周囲で計算が不安定となることが報告されているため，本コードでは FDTD と PIC を同じタイムステップ

38で計算することとし，2.5×10-13 sec とした．各種計算条件を表 2-3 に示す．
表 2-3 計算条件(FDTD-PIC-MCC)

主グリッド : 0.5 mm

メッシュサイズ

局所グリッド : 0.1mm
主グリッド : 2.5 ×10-13 sec

タイムステップ

局所グリッド : 0.5 ×10-13 sec

超粒子の重み1.0 ×105

超粒子数

360 万個

1 セルあたりの粒子数

100 個

また，FDTD の境界条件であるが，アンテナの根本のみに Mur の１次吸収境界をかけている．

参考文献
(2-1)

内藤裕志, プラズマ・核融合学会誌, 74 (1998) ,470.

(2-2)

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

石黒静児, プラズマ・核融合学会誌, 74 (1998) , 591(2-4)

高村秀一：“プラズマ理工学入門”(森北出版, 1997)(2-5)

NIFS DATABASE 核融合科学研究所データベース(http://dbshino.nifs.ac.jp/index-j.html)(2-6)

Mersenne Twister Home Page : (http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html)(2-7)

V. Vahedi and M. Surendra, Computer Physics Communications, 87(1995), 179.
(2-8)

電気学会・マイクロ波プラズマ調査専門委員会：“マイクロ波プラズマの技術” (オーム社, 2003)(2-9)

K. Koura, Phys. Fluids, 29(1986), 11.

(2-10) 橋本修, 安部琢美：“FDTD 時間領域差分法入門” (森北出版, 1996)(2-11) 宇野亨：“FDTD 法による電磁界およびアンテナ解析” (コロナ社, 1998)
39

(2-12) G. Mur, IEEE Trans. Electromagnetic Compat., 40 (1998), 100.
(2-13) K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.
(2-14) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical
Physics, 55 (1971), 4100.
(2-15) 渡邊裕樹, “マイクロ波放電式電子源内の電子加熱解析”.STEP-2013-014(2-16) M.W.Chevalier,R.J.Luebbers, and V.P.Calle,”FDTD Local Grid with
Material Traverse”IEEE Trans. Antennas Propagat.,45,3,pp.411-421,March
1997(2-17) 小野寺範義, 竹ヶ原春貴, 西山和孝, 船木一幸, 國中均: マイクロ波放電型
中和器の電子放出機構, 日本航空宇 宙学会論文集, 第 49 巻, 第 564 号,
pp.27-31, 2001.

40

第3章
3.1.

結果と考察
サブグリッド適用コードの妥当性検証

新しく作成したサブグリッド法を適用した FDTD コードが正しく動作しているかどうか検証した．計算条件を表 3-1 に示す通りで，投入エネルギーは一定にした．
このとき図 3-1 に示す範囲のポインティングベクトルを各コードで計算し，比較を行った．図 3-2 に結果を示す．極値において差が見られるが，波形は概ね一致しており，サブグリッド法を用いたコードの妥当性が確かめられた．
表 3-1

計算条件

主グリッド : 1.0 mm

メッシュサイズ

局所グリッド : 0.2mm
主グリッド : 1.0 ×10-12 sec

タイムステップ

局所グリッド : 0.2 ×10-12 sec
2.45 GHz

マイクロ波周波数

図 3-1

計算体系

41

0.5
without Subgrid
with Subgrid

0.4
0.3
Energy [W]

0.2
0.1
0
-0.1
-0.2
-0.3
0

0.2

0.4

0.6

0.8
Time [nsec]

1

1.2

1.4

1.6

図 3-2 通常の FDTD とサブグリッド法を適用した FDTD の比較

3.2.

従来コードとの電界強度の比較

本研究室で開発された従来のコードとサブグリッド法を適用したコードで比較を行った．今回，局所グリッドの範囲は x=2 mm~20mm，y=2 mm~20 mm，z=9mm~10 mm とし， FDTD のみの結果を比較するため電子を配置せずに計算を行った．z=9 mm の xy 断面を比較した．各アンテナ形状における，マイクロ波の第2 周期から第 3 周期(0.82 ns~1.22 ns)までで平均した電界強度分布を図 3-3 に示
す．形状について比較すると星型，十字型はアンテナ先端の鋭角部分に強電界が発生していることがわかる．このため，円盤型よりも十字型，星型の方がアンテナ先端の電界強度は高い．この形状による強度の差は形状再現率が向上したサブグリッ
ド適用コードでは顕著に現れており，鋭角部分では電界強度はより強く，湾曲部分
ではより弱くなっている．
サブグリッド適用により形状再現性が向上し，鋭角部，湾曲部が強調された円盤型
アンテナ，星型アンテナにおける z=9 mm，x=11 mm における y 方向の電界強度分布を図 3-4 に示す．両コードにおいて電界強度はアンテナに近づくにつれて上昇
42

し，アンテナから離れると急激に上昇したのち，緩やかに減衰している．ここで，アンテナ最近傍の電界強度が極値でないのはアンテナ表面の電界が 0 の点を含めて空間的に平均された値を出力しているためである．アンテナ近傍で極値をとるまでの距離は，従来コードだと 1 mm であるのに対し，サブグリッド適用コードでは
0.2 mm と両コードとも 2 セルであり，この距離はメッシュサイズに依存していることがわかる．また，サブグリッド適用コードにおいて星型ではアンテナ近傍の電界強度が増幅，円盤型では減衰していることがわかる．
これらの結果の違いは，メッシュサイズが小さくなったことでデータ数が増え，境界部の計算を明細に行うことができるようになったためであり，サブグリッド適用
コードの有意性が示された．

(a)円盤型43

(b)十字型(c)星型
図 3-3 各コードで計算された電界強度分布の比較

44

1.4E+05

Antenna at local grid

Antenna at main grid

STAR with Subgrid

STAR without Subgrid

Round with Subgrid

Round without Subgrid

Electric field [V/m]

1.2E+05
1.0E+05
8.0E+04
6.0E+04
4.0E+04
2.0E+04
0.0E+00
4

6

8

10

12
y [mm]

1416

18

図 3-4 円盤型，星型における y 方向の電界強度分布

3.3.

FDTD-PIC コード

電子を配置した状態で，FDTD-PIC コードで計算を行った．各図は x=11 mm の
zy 断面，z=9 mm の xy 断面を示している．

3.3.1.

電界強度分布図 3-5 にマイクロ波第 9 周期~第 10 周期(3.7 ns~4.1 ns) までで平均した電界強度
分 布 を 示 す ．ア ン テ ナ 後方 およ び ア ン テナ 表 面と オ リ フ ィス の 間の 空 間 に
1.2×104V/m 程度の強電界が発生している．アンテナ表面積の大きい円盤型，星型
は強電界の範囲が広いのに対し，表面積の小さい十字型では範囲が狭く，アンテナ
表面積の大きさが電磁波の伝搬に有効であることがわかる．
しかし，この分布はアンテナ周辺の電界強度は弱くなっており，電子を配置せずに計算を行った時と異なっている．これは PIC により計算された電流密度がサブグ
リッド法を適用した FDTD に正しく反映されていないからだと考えられる．このため，局所グリッドに対応した PIC コードの開発が今後の課題としてあげられる．

45

(a)円盤型(b)十字型46

(c)星型
図 3-5 電界強度分布

3.3.2.

電子密度分布図 3-6 に各アンテナ形状における マイクロ波 10 周期後の 4.5ns 時の電子密度分布を示す．共通してアンテナ，壁面近くでは電子密度が低くなっている．これは，シース電位以上の電子がアンテナ，壁面に衝突しているからだと考えられる．しかし，各アンテナ形状による違いははっきりと見られない．これは解析時間が短く電離があまり起きていないことによるものであるため電子密度の差を見るためには
今後解析時間を伸ばす必要がある．

47

(a)円盤型(b)十字型48

(c)星型
図 3-6 電子密度分布

3.3.3.

電子エネルギー分布

図 3-7 に同じく 4.5 ns 時の電子エネルギー分布を示す．各アンテナ形状においてアンテナ先端で電子が加熱されている．円盤型では境界が直線で構成されている部分は電子エネルギーが低く，ジャギー部分は高くなっている．また，十字型，星型でもアンテナ先端の鋭角部分が加熱されていることがわかる．このことから，ア
ンテナ先端の鋭角部分が電子加熱に有効であることがわかる．

49

(a)円盤型(b)十字型50

(c)星型
図 3-7 電子エネルギー分布

3.3.4.

EEDF

アンテナ近傍の領域での電子エネルギー分布関数を磁場に垂直方向と平行方向に分けて算出した．バックヨーク側の ECR 層を含むよう算出領域は図 3-8 の範囲とした．EEDF を図 3-9 に示す．
磁場に垂直方向，平行方向ともに星型が最も高エネルギー電子数が多く，次いで十字型，円盤型となっている．これは実験におけるアンテナ表面の電子温度の結果と同じ傾向にある．

51

図 3-8 EEDF 算出領域

図 3-9 各アンテナ形状における電子エネルギー分布関数

3.3.5.

電子損失52アンテナに到達した電子を損失としてカウントした．図 3-10 に各アンテナ形状
における電子損失を示す．円盤型，十字形，星型の順で電子損失が多いことがわかる．これは高エネルギー電子数と逆の傾向を示している．このことから電子損失によってアンテナ周辺の電子エネルギーが低くなったと考えられる．各アンテナのメ
ッシュ数は表 3-2 に示す通りであり，円盤型，星型，十字形の順にメッシュ数が少ない．従って電子損失はアンテナ表面積に依存していない．しかし，ECR 層を結
ぶ磁気チューブと接触する直径 6 mm 以上の領域のメッシュ数は表に示す通り，
円盤型，十字型，星型の順に少ない．つまり，電子が往復する磁気チューブと接触するメッシュ数が電子損失の数に依存すると考えられる．

Number of particles

2400
2200
2000
1800
1600
star

round
Antenna shapecross

図 3-10 電子損失

表 3-1 各アンテナ形状のメッシュ数

アンテナ形状

円盤型

星型

十字型

全体メッシュ数

498

278

198

磁気チューブ接触メッシュ数

274

68

72

53

第4章

結論

本研究ではマイクロ波放電式イオンスラスタのアンテナ形状が放電室内部に与える影響を解明するためサブグリッド法を適用した数値解析コードを開発し，電子
加熱解析を行った．得られた結論を以下に示す．

1) サブグリッド法を適用すると従来のコードよりアンテナ形状再現性が向上し，
アンテナ先端の鋭角部分，湾曲部分がより強調されるため，アンテナ先端部分の電界強度の差が大きくなった．このことから，アンテナ形状依存性を解明するためにはサブグリッドを用いたコードが必要であることがわかった．
2) 円盤型，十字型，星型のアンテナ形状のうち，最も電子加熱に有効なのは星型であった，これは鋭角なアンテナ先端部分に強電界が発生すること，磁気チューブに接触するアンテナのメッシュ数が少なく電子損失が小さいことが寄与している．
しかし，本解析コードでは，局所グリッドを FDTD コードにしか適用できておらず，粒子計算には局所グリッドが対応していないため PIC コードにも局所グリッドを対応させる必要がある．

54