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

修 士 論 文論文名
マイクロ波放電式中和器の引き出し電流電圧特性解析

氏

名

指導教員名

楠田 将士

中島

秀紀

山本

直嗣

森田

太智

目次
第 1 章 序論 ............................................................................................. 1
1.1. 背景 ...................................................................................................... 1
1.2. イオンエンジンの原理 ........................................................................... 3
1.2.1. マイクロ波放電式イオンエンジン ...................................................... 4
1.2.2. マイクロ波放電式中和器 .................................................................... 7
1.3. マイクロ波放電式中和器の性能低下 ....................................................... 9
1.3.1. マイクロ波放電式中和器の性能低下 ................................................... 9
1.3.2. スパッタリング (1-13) ............................................................................ 9
1.3.3. 中和器電圧の上昇 (1-2) ........................................................................ 121.4. 研究目的 ............................................................................................. 15

第 2 章 数値解析手法 ............................................................................. 18
2.1. 解析全体流れ ...................................................................................... 18
2.2. プラズマ粒子計算手法 (2-2) .................................................................... 21
2.2.1. 基礎式 .............................................................................................. 21
2.2.2. メッシュサイズ・タイムステップ (2-3) ............................................... 21
2.2.3. 粒子計算 PIC 法 (2-4) .......................................................................... 23
2.2.4. 衝突計算 MCC 法 ............................................................................. 23
2.3. FDTD-PIC-MCC .................................................................................. 31
2.3.1. プラズマ粒子計算 ............................................................................. 31
2.3.2. 電磁波解析手法 (2-13)(2-14) .................................................................... 34
2.3.3. FDTD-PIC-MCC のフローチャート .................................................. 42
2.4. Full-PIC-MCC ..................................................................................... 44
2.4.1. 運動方程式解法 ................................................................................ 44
2.4.2. PIC 法重みつけ計算 .......................................................................... 46
2.4.3. 電場解析手法 (2-18) ............................................................................. 49
2.4.4. Full-PIC-MCC のフローチャート ..................................................... 53
2.5. 弱連成 ................................................................................................. 55
2.6. 計算体系および静磁場分布 .................................................................. 55
2.7. 各種物理条件と計算条件 ...................................................................... 58
2.7.1. 物理条件 .......................................................................................... 58
2.7.2. 計算条件 .......................................................................................... 58

第 3 章 結果と考察 ................................................................................ 62
3.1. FDTD−PIC-MCC ................................................................................. 623.1.1. 電子軌道 .......................................................................................... 62
3.1.2. 空間電子密度分布・エネルギー分布・電界強度分布 ........................ 64
3.2. Full−PIC-MCC .................................................................................... 69
3.2.1. 空間エネルギー分布 ......................................................................... 69
3.2.2. 空間密度分布.................................................................................... 70
3.2.3. 空間ポテンシャル分布および静電場分布 .......................................... 71
3.2.4. 壁面電流および引き出し電流 ........................................................... 73
3.2.5. IEDF ................................................................................................ 75

第 4 章 結論 ........................................................................................... 77

第1章
1.1.

序論

背景

2010 年に宇宙探査機“はやぶさ”の帰還により, 日本での宇宙への関心度は大きく高まった. はやぶさのミッション目的は, 小惑星“イトカワ”に到達し, その岩石のサンプルを地球に持ち帰るというものであった . ミッションに要求されていたのは , 長時間の安定動作かつ高いΔ V(速度増分)であった. それを実現したのが,マイクロ波放電式イオンエンジン”μ10”である. なぜミッション実現足り得たのかを以下に示す.
マイクロ波放電式イオンエンジンとは, 電気推進の一種であり, 電気エ
ネルギーを用いて推進剤を電離させ，生成されたプラズマを様々な方法で加速し，その反作用で推力を得る (1-1)．これに対し，従来の化学推進機は化学反応を利用し作動流体を噴出し推力を得る．電気推進機は化学推進機と比較して得られる推力ははるかに小さいが，燃費を示す指標である比推力(Specific impulse：Isp )を一桁以上大きくすることが可能である．つまり，電気推進機は推力が小さいため加速に時間がかかるが，消費する燃料は化学
推進機と比較し非常に尐なくて済む．電気推進機イオンスラスタやアーク
ジェット，MPD スラスタ，RF スラスタ，パルスプラズマスラスタ，ホールスラスタなどが開発・研究されている．特に，イオンスラスタは他の電気推進機と比較しても比推力が高く(Isp 数 1000s 以上)，長期間のミッションである深宇宙探査に適している．この中でもマイクロ波放電式イオンエンジンとは, 詳しい原理は後述するが, プラズマ生成に電子レンジでお馴
染みのマイクロ波を使用し, 長寿命に特徴があったイオンエンジンであり,
1997 年からの 18000 時間と 2000 年からの 20000 時間の地上耐久試験を経て, 2003 年打ち上げの「はやぶさ」に搭載された推力が 8 mN(ミリニュートン), 比推力が 3000 秒の μ10 である. はやぶさでのミッションにおいて,マイクロ波放電式イオンエンジンμ10 は, 一つのシステムとしての積算作動時間は 4 台で 40,000 時間に迫り，世界一の宇宙実績と高い信頼性を示した (1-2)
そして, 2014 年 12 月 3 日, はやぶさ 1 のミッションの成功に湧いた 4 年後,“はやぶさ 2” が打ち上げられ,地球を飛び立ちイオンエンジン μ10 に火をつけた. 今回のミッションは, 地球からの直線距離がイトカワの約 2 倍で1

約 3 億 km の小惑星”1999JU3” からのサンプルリターンである (1-3). 改良されたμ10 は推力が 8mN から 10mN という前回の 25%アップされている.
話は変わって, はやぶさ 2 の相乗り衛星として, “PROCYON”という超小型衛星がはやぶさ 2 と打ち上げと同時に打ち上げられた. この衛星のミッション目的は, 50kg 級の低コストで短期開発可能な超小型宇宙機で深宇宙探査ができることを実証することである. この衛星の背景には, 近年，
数 10kg 以内かつ数億円以内の地球周回の超小型人工衛星が，超小型化・低コスト化技術として大学を中心に研究・開発・軌道上実証され，実用的なレベルに達しつつあり (1-4) ,この技術を持ってして, 超小型衛星が深宇
宙探査に利用されることは, 低コストかつ短期間開発ができれば,ミッシ
ョ ンの 頻 度 を 上げ , 宇 宙開 発 の 加 速が 促 さ れる も の と され て い る . こ の“PROCYON”のミッションを支えるものが, I-COUPS (Ion thruster and COld
gas thruster UnifiedPropulsion System)と呼ばれる低消費電力かつ小型を実現したイオンエンジンシステムであり,このイオンエンジンにもマイクロ波
放電式を採用しており,超小型衛星に載せられる程度に低電力と小型化を実現している. マイクロ波放電式イオンエンジンは, 小型推進機としての信頼性や実用性が高まっているのである. この I-COUPS に代表されるように, 超小型衛星の市場活性化に伴う超小型衛星への推進系の必要性は年々増している (1-5), . その理由が, 超低高度衛星等のミッションの需要の高まりであり, 観測分解能の向上やセンサの小型化・低消費電力化が図られ,
災害監視，地球環境観測，サイエンスなどの様々な分野への応用が期待されている (1-6). 更に, 宇宙デブリつまり宇宙空間にある衛星の残骸などのゴミの観点から言うと, 現在, 衛星には, デオービット機構と呼ばれる地球近傍の衛星において, ミッションが終了または, 何らかのトラブルによりミッション続行が不可能と判断された場合, 地球に向けて軌道をはずれ,
大気圏で燃やしてしまう機能が付いているものもある. しかし, デオービ
ット機構の不具合が原因でデオービットできないという状況も考えられる.
その一方, 推進機を搭載していると, 衛星は, 地球との微重力を推進機にて軌道修正, 姿勢制御を行い, ミッション期間は, 推進剤の総量で決定される. そのために推進剤が枯渇した途端に, 地球との微重力に耐えられずに自然原理的にデオービットされるため, 衛星自身がデブリにならないた
めのシステム的な特徴を持っているといえよう. 国際的な超小型衛星市場においても, 既に市場は激化の一途をたどっており (1-7), ミッションの幅を広げるためにも, 小型推進機の役割は重要となってくると言える. 以上より, マイクロ波放電式イオンエンジンシステムは, 超小型衛星用の小型推2進機として, 市場からの高い需要と信頼性や実用性に応えることが可能である.
さて, はやぶさにおけるマイクロ波放電式イオンエンジンμ10 にも問題点があった. その問題点とは, 中和器である. 中和器とは, 詳しい原理は,次項で説明するが, 電子を放出することでイオンビームを中和する装置である. イオンビームを中和しないと, エンジンが負に帯電することで, 排出したイオンが逆流し, 性能が得られない為である. μ10 は単独で 14,000
時間以上という最長作動実績を上げていた中で , 中和器が务化の兆候を見せた (1-1) ． 4 機搭載されたうちの 1 機がイオン源（イオン噴射装置）の点火不良, 3 機が中和器（電子噴射装置）が想定より早く寿命を迎え運用停止の危機に瀕することとなった.バックアップ回路により, 残る 1 機の中和器とイオン源を組み合わせた「クロス運転」を行うことで危機を脱出
し, 地球帰還を果たし, マイクロ波放電式イオンエンジンシステムのロバ
スト性の高さが証明されたのだが, やはり中和器の寿命が今後の不安要素であるのは, 明らかである. 長期のミッションに適しているからこそ, マ
イクロ波放電式イオンエンジンシステムの長寿命と言うのは, 必須であり,信頼性向上には, 中和器の性能向上が鍵となっている.

1.2.

イオンエンジンの原理

マイクロ波放電式中和器の説明の前に, イオンエンジンの原理を通して,
中和器の位置づけおよび役割を詳しく説明する .

イオンエンジンの概略

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

3

図 1-1 イオンエンジン概略図

1.2.1.

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

イオンエンジンの中でのマイクロ波放電式を説明していく. イオンスラ
スタはプラズマの生成方法によって直流放電式，RF (Radio Frequency) 誘導放電式，マイクロ波放電式等に大別される．本研究の対象であるマイク
ロ波放電式イオンエンジンの概念図を図 1.2 に示す．マイクロ波放電とは，マイクロ波帯域の交流電場によって電子を加速し，中性粒子と衝突電離に
より電子の数が増加し気体がプラズマ化され放電が維持されることである．

4

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

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

(1.2.1)

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

rL 

mv
eB

(1.2.2)で与えられる．ここで v ⊥ は電子の B に対する垂直な速度成分の大きさである．
5

プラズマを構成している荷電粒子は全て反磁性体である．そのためサイクロトロン運動の回転の方向は，外部磁場の向きに対して荷電粒子の回転によってできる磁場が常に逆を向く方向である．つまり電子は磁場に対して右回りの回転運動を行う．この回転運動の角周波数はサイクロトロン角
周波数と呼ばれ ce 

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

E

電子

B
E

B

電子

電子
E

B
B

電子
E

図 1-3 ECR の原理

6

1.2.2.

マイクロ波放電式中和器

マイクロ放電式イオンエンジンには，そのイオンビームに同量の電子を供給し電気的に中和する中和器という装置が搭載される．これにより，エンジンや宇宙機本体が負に帯電し噴出したイオンを引き戻すことなく，運転を継続することが出来る．中和器が故障し衛星からイオンのみが放出されることになるとイオンビームの抽出が不可能となる．
一般的な直流放電型イオンエンジンでは中和器としてホローカソードが使用されており，スラスタ本体と同様に電極の損耗や電源構成の複雑化が故障の要因となっている．ヨーロッパで盛んに研究が行われている高周波放電型イオンスラスタ (1-11) でも中和器にはホローカソードが使用されているが，この問題は解決されていないままである．
このマイクロ波放電は直流放電と比べた場合，以下のようなメリットがある．
1. ヒーターを用いないため，予備加熱が不要となり素早い起動が出来る．
2. マイクロ波放電式イオンエンジン駆動用のマイクロ波電源が 1 台
あれば中和器のプラズマも同時に生成することが可能なため，システム構成の簡素化や信頼性の向上，また衛星の重量軽減が図れる．
マイクロ波放電式中和器は深宇宙探査機「はやぶさ (1-12)」に搭載され十分な中和性能が得られたという実績がある．当研究グループでもマイクロ波
放電型イオンエンジンを開発してきた． 図 1.4 はマイクロ波放電型中和器の外観写真である．

7

図 1-4 マイクロ波放電式中和器の写真

図 1-5 マイクロ波放電式中和器運転時

81.3. マイクロ波放電式中和器の性能低下
1.3.1.

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

小惑星探査機「はやぶさ」に搭載されたマイクロ波放電型イオンエンジン”µ10”は JAXA 宇宙科学研究所の電気推進工学部門にて研究開発され,そして世界で唯一の実用化された ECR 放電型イオンエンジンである. 宇宙科学研究所のマイクロ波放電型イオンエンジンの開発は１９８０年代から始められ, １９８９年に１号機（Y-I）が完成する. その後改良が重ねられ, ２００１年には「µ10」において, 推進剤利用効率にして４倍, イオン
生 成 コ ス ト に し て 1/10 も の 性 能 向 上 と ２ 万 時 間 の 耐 久 性 を 達 成 し た .
“µ10”に搭載された中和器は, 先にも述べたように, 従来の中和器と異なり, プラズマ化の際に電極を必要としないマイクロ波放電を採用している.
投 入 マ イ ク ロ 波 電 力 は 8 W, ガ ス は Xe を 使 用 し , そ の 流 量 は 0.5
sccm(standard cubic centimeter per minute : 標準状態における 1 分あたりの体積流量), 中和するための電子電流が, 実際に推力を得るために陽イオンを放出しているイオン源と同等の 135 mA で中和器電圧が 22 V 程度となる.

1.3.2.

スパッタリング (1-13)

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

10

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

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

1 keV~300 keV>100 keV

Sigmund によるスパッタリングの分類

11

1.3.3.

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

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

(a) 通常の状態
図 NV(b) 中和不良時の状態
中和器の電位の模式図

13

図 X

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

図 Y 中和器電圧の時間変化

14

1.4.

研究目的

マイクロ波放電式中和器の現在の課題をまとめると,① プラズマのエネルギー獲得機構② シース付近のプラズマの挙動③ 損耗メカニズム
これら 3 つの解明が急務である. これを解析するために, 本研究室では,数値解析によるアプローチを試みてきた. これは, マイクロ波放電式中和器のような小さい空間スケールの物体には, プローブ等の測定機器を掃引すると, プラズマの擾乱を引き起こし, 正しい結果が得られにくい. 更に,損耗過程のような何万時間に及ぶ地上試験を行うのは, 時間と金銭コストが莫大になる. 数値解析であれば, 時間的金銭コストを抑えつつ, 実験上
では, 見られないダイナミクスを解析することが可能である. 以上の理由で, 数値解析は有用であるといえる.
数値解析の側面からイオンエンジンや, 他静電スラスタ等の解析は内外で, 多く実施されているが, マイクロ波解析, イオン挙動, 損耗解析の異なる時間スケールの物理現象を連成的に解析するまでは行われていない .
そのために,
よって, 本研究目的は, マイクロ波放電式中和器内部の物理現象を統合的に解析できるツールの開発であり, 本論では, マイクロ波中和器の代表的な特徴のひとつである引き出し電圧よって, 電子電流が高電流モードと低電流モードとをモード遷移する()という (1-17) 代表的な特徴をとらえることを可能にすることを目的とする.

15

図 1-6 マイクロ波放電式中和器の IV 特性

参考文献
(1-1)

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

細田聡史，國中均，”イオンエンジンによる小惑星探査機「はやぶさ」の帰還運用”，J.Plasma Fusion Res. Vol.86,No.5,(2010),pp.282-292(1-3)

吉川, 真; はやぶさ 2 プロジェクト準備チーム”小惑星探査ミッシ

ョン「はやぶさ 2」” 第 11 回宇宙科学シンポジウム: 講演集録 Vol
11 ,2011(1-4)

Shinichi Nakasuka, "From Education to Practical Use of Nano-satellites- Japanese University Challenge towards Low Cost Space Utilization –", 8th
IAA Symposium on Small Satellite for Earth Observation, Berlin, Germany,
2011.
(1-5)

小泉宏之, “小型衛星用推進機－マイクロスラスタ－の開発と今後の課題”, JSASS-2014-4333(1-6)

高 山 慎 一 郎 , “ 超 低 高 度 衛 星 技 術 試 験 機 （ SLATS ） の 開 発 状 況 ”,

JSASS-2014-4075(1-7)

中 村 友 哉 , “超 小 型 衛 星 ビ ジ ネ ス が も た ら す 宇 宙 イ ノ ベ ー シ ョン” ,JSASS-2014-4323(1-8)

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

2005.
(1-9)

茂原正道・鳥山芳夫, “衛星設計入門”, 2002.

(1-10)

大道渉，”マイクロ波放電式中和器のイオン電流測定”，

STEP-2011-033(1-11)

N. Wallace, P. Jameson, IEPC, 2011, 327(1-12)

I. Funaki, H. Kuninaka and K. Toki, J. Propul. Power, Vol. 20, (2004),

pp.718-726.
(1-13)

Kirk A. Zoerb, “Differential sputtering yields of refractory by ion

bombardment at normal and oblique incidence”, Thesis, Department of
Mechanical Engineering, Colorado State University, Fort Collins, Colorado(2007)(1-14)

大道 渉， 小泉 宏之， 国中 均, 細田 聡史， 西山 和孝, “マイク ロ 波 放 電 式 中 和 器 の 性 能 低 下 メ カ ニ ズ ム と そ の 実 験 的 検 証 ”，
STEP-2010-058(1-15)

大道渉“マイクロ波中和器の内部务化機構の解明と長寿命化に関する研究”,STEP-2013-008(1-16)Noriyoshi Onodera,et al., “Electron Emission Mechanism of Microwave

Discharge Neutralizer,” IEPC-99-160,2002(1-17)

I. Funaki., H. Kuninaka, "Overdense Plasma Production in a Low-power

Microwave Discharge Electron Source," Jpn. J. Appl. Phys., Vol. 40, 2001,
pp. 2495-2500.

17

第2章
2.1.

数値解析手法

解析全体流れ数値解析手法について, 述べていく. 各計算モデルについて, 詳しく述べる前に, 解析ツールの全体図を述べる. 本研究室にて開発された数値解
析コードは, 3 つの物理現象に基づき, 3 種の計算モデルを用いている. 表
2-1 に示す.
表 2-1
物理対象

モデル名

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

FDTD-PIC-MCC

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

Full−PIC−MCC

中性粒子, スパッタ粒子の輸送および壁面損耗解析

DSMC

これら 3 つの計算モデルを弱連成的に解析させることで,プラズマ挙動の解明および寿命の評価ができる. 図 2-1 に解析全体のフローチャート（もしくは解析全体の流れ）を示す. 1 章で述べたが, マイクロ波放電式中和器の基本原理は, 永久磁石とマイクロ波の ECR 共鳴であるため, 静磁場の情報が必要である．中和器に永久磁石を配置した際の静磁場分布は先に市販ソフトで計算しておく(後述)． この静磁場データを持って，以降のFDTD-PIC-MCC お よ び Full − PIC-MCC の 計 算 が 行 わ れ る ．
FDTD-PIC-MCC では，一様状態から電子とマイクロ波の吸収が計算され，電子エネルギー獲得確率が算出される．これは， 単位時間単位体積あたりの電子のエネルギー増分を確率密度分布として算出後， 電子のマイクロ波
吸収量の効果として，Full−PIC-MCC に渡す． Full−PIC-MCC では，イオンと電子を一様状態から計算させ，静電場の効果による粒子の輸送を解析する．プラズマ密度が変化するとカットオフ密度(電磁場がほとんど浸透
不可能である密度)以上の分布が出現してきて，密度分布が変化するため，
電子のマイクロ波吸収量が変化する． 実際に，Full-PIC での計算において，初期電子エネルギー獲得確率でエネルギーを与え続けると，特定の領域のみでマイクロ波から電子へのエネルギー伝達が行われ続け多結果として,電子温度が上昇し,結果として,計算が発散してしまう．このため，プラズマ密度の変化を FDTD-PIC-MCC に Full−PIC-MCC に再びフィードバ
ックすることで，カットオフの効果を取り入れ, 電子密度の空間分布によ18

るマイクロ波の伝搬の変化も考慮することが可能となった エネルギー保存のために，Full−PIC 計算での電子エネルギーも同様にフィードバックしている．Full−PIC の計算において， 逐次プラズマ密度の変化とともに，
FDTD-PIC にフィードバックしながら繰り返し計算し，定常状態を目指す．
このように，時間スケールの異なる現象過程をそれぞれ独立させて計算さ
せ，それぞれの結果を受け渡しながら計算する事を弱連成解析という．本来は，タイムステップ毎に，マイクロ波と静電場を計算すれば，
Self-Consistent に数値解析を行うことができるが，時間スケールの異なる現象を計算することは，非常に高い計算コストが生じてしまうため，計算時間短縮しながら，解析物理対象のモデリングを行えるため，弱錬成解析は有用と言える．Full−PIC-MCC にて定常状態まで解析後，Full−PIC で算出した各壁面へのイオンのエネルギー分布関数(IEDF)を算出する． この結果を用いて，壁面でのスパッタ粒子の発生確率を算出でき，DSMC にてスパッタ粒子および中性粒子の輸送を解析する．1 章にてアンテナの周りの絶縁体へのスパッタ粒子の堆積が絶縁の抵抗値を下げ，その影響により中和器性能低下をもたらすと述べたが，これを模擬し，DSMC によって得られたスパッタ粒子の堆積量によって，寿命を算出するものとする． また，アンテナ部へのスパッタ粒子の堆積がマイクロ波の伝搬モードへどのように影響をおよぼすかを与える効果を取り入れるために，
FDTD-PIC-MCC にアンテナの周りに抵抗皮膜を挿入する効果を取り入れることで，従来の定常状態解析から，損耗状態解析を行えるものとして，寿命評価だけでなく，損耗状態時の中和器内部のプラズマダイナミクスを解析することを可能にしている．尚，本論では，FDTD-PIC-MCC と Full
−PIC-MCC の弱錬成解析による解析のみを扱う．DSMC 解析は笹川によって，解析がなされた (2-1)．
また，DSMC から FDTD-PIC-MCC への連成は，江川らが，別途開発中である19

図 2-1 解析全体フロー

20

2.2.

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

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

2.2.1.

基礎式

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

dv
 q E  v  B 
dt

(2.1)

dx
v
dt

(2.2)

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

2.2.2.

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

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

 0 k bTene e 2

に対して，

21

(2.3)

  De(2.4)

であることが望ましい．また，タイムステップはプラズマ周波数 ω pe ，ne e 2
 pe me  0(2.5)

t   pe1

(2.6)

に対して，

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

t  0.2 pe1(2.7)

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

c 

eB
me

(2.8)

に対して，

t  c

1

(2.9)

である必要がある．
また，粒子は時間刻みごとに１セル以上飛び越えてはならないという制約がある．特に荷電粒子は格子ごとに規定される電磁力が働くため，この条件が必須となる．
これを Courant の安定化条件という．3 次元の際，時間ステップ t は，
1
Dt £
1
1
1
(2.10)v++
2
2
2(Dx ) (Dy ) (Dz )を満たさなければならない．この条件は極めて厳しく，わずかでも満足しなければ計算が不安定になる．実際には，この 10 分の 1 程度が望ましいとされている．ここにおいて，v は，粒子の最大熱速度または，電磁波速度つまり光速である．FDTD-PIC と Full-PIC において，条件が異なるので各モデルの計算条件において説明する．

22

2.2.3.

粒子計算 PIC 法 (2-4)

2.2.3.1.

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

2.2.3.2.

重み付け

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

2.2.4.

衝突計算 MCC 法

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

2.2.4.1.

衝突種別

本モデルにおいて荷電粒子は電磁場の影響を受けて，中性粒子と衝突する．プラズマは電子，イオン，中性粒子からなるとする．電子とイオンの数密度は中性粒子の数密度より非常に小さいため，電子-電子，電子-イオ

23ン，イオン-イオンの各衝突は無視できるものとし，電子-中性粒子のみの衝突を取り扱う．衝突種別は弾性，励起，電離衝突のみを仮定する

2.2.4.2.

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

面  ，厚さ l ，その物質の数密度 N の板を考える．エネルギー E の粒
子ビームが，面  に垂直に当るものとする．衝突をせずに板を透過してくる粒子数は，入射粒子数より減尐する.その減尐度は  を衝突の断面積とすれば，dn  nNdl(2.11)と表わされ，これよりn  n0 exp  Nl (2.12)となる．ここで，n0 は面  （ l  0 ）でのビーム中の粒子である．それゆえ，（ l ,l  dl ）間で最初の衝突が起こる確率は，
Pinitial l dl  exp  Nl Ndl(2.13)となり， l 以下で最初に衝突が起こる確率は，Pl    Pinitial l dl
l0 1  exp  Nl (2.14)となる．最初の衝突までに粒子が走る平均距離  は Pl の一次のモーメン
トである．
   lP l dl
0



1
N

(2.15)

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

で表わされ，これを変形して，となる．ここで， r は一様乱数であるから 1  r を r とすれば，l  ln(r )(2.18)と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を  と

24

すれば，この粒子が１秒間に気体分子と衝突する回数  は，



v


 Nv

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

2.2.4.3.

衝突確率・衝突判定

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

(2.20)

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

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

(2.21)(2.22)

ここで， ng asは中性粒子の数密度， v m は粒子 m の速度，  t o t a lは全衝突断面積である．ある乱数 r1（ 0  r1  1）を用いて， Pcollision,m  r1 のときに粒子 m は
タイムステップ t の間に衝突が起こるとする．

2.2.4.4.

衝突種別判定乱数 r2 を用いて以下のように衝突の種類を決める．ここでは例として，

電子と中性粒子の衝突について示す．弾性衝突断面積を s elastic ，励起衝突断面積を  excitation ，電離衝突断面積を  ionization とする．

 elastic
 total

であれば弾性衝突 elastic  excitation
 r2  elastic
 total totalであれば励起衝突 elastic   excitation  excitation   ionization
 r2  elastic
 total totalであれば電離衝突 0  r2 

ここで，s total = s elastic + s excitation + s ionization である．このように断面積の比と乱数を対応させ，衝突の種類を決める．計算では中性粒子としてキセノンを用25

い，その断面積データ (2-7) を図 2-2 に示す．今回の計算では，多価電離は考慮していない．

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

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 として26

よい．よって，電子の質量を m ，中性粒子の質量を M とすれば，次式が成立する (2-7) ．
v e 

m  M cos M sin 
ve 
h
mM
mM

(2.25)

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

(2.26)

vex vey cos   ve vez sin verv v cosj + ve vey sin j
hz = - ex ezver(2.27)(2.28)

ここで，
ver 
ve 

2
vey
 vez2(2.29)

2
2
vex
 vey
 vez2(2.30)

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

Escattered,e  Eincident,e(2.31)

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

Eincident,e sin 0     21  Eincident,e sin 2 ln 1  Eincident,e 2

これを解けば，

2 1  1  Eincident,e  3
cos   1 
Eincident,er

(2.32)(2.33)となる．この値を(2.25)式に用いて衝突後の速度を求める．このとき，任意の角は  2 r4

27(2.34)

として与える．
②励起衝突
基底状態にある中性粒子をある準位に励起したとき，電子が失う励起エ
ネルギーを Ee x c i t at i oとし，入射電子のエネルギーを
Ei n c i d e net ，励起後のエネル
n
,
~
ギーを
，励起後の速度を v とする．エネルギー収支は，Escattered,eEscattered,e  Eincident,e  Eexcitation(2.35)と表わされ，励起後の速度はEv~  v e 1  excitation
Eincident,e(2.36)となる．励起衝突を励起と弾性衝突に分けて考えれば，励起衝突でエネル~ の速度で弾性衝突をすると考えることがでギーを失った後に(2.36)式中 v~ を用いて，衝突後の速度  を求める．このときのきる．(2.25)式中 に vveve

散乱角には(2.33)式を用いるが(2.33)式中の Ein cid enet,には(2.35)式から求められる Escattered,e を用いる．励起エネルギー Eexcitation は，キセノンの場合 8.34
eV を用いている．
③電離衝突
入射電子のエネルギーを Ei n ci d enet ，散乱電子のエネルギーを
Es cat t erede ，生成,,
電子のエネルギーを Ecreated,e ，電離エネルギーを Eionization とする．エネルギー収
支はEscattered,e  Ecreated,e  Eincident,e  Eionization(2.37)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる．

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

(2.38)

ここで， B はガスに固有の値でキセノンでは 8.7 eV となる．
また，入射電子が電離によって失うエネルギーを E とすると，
E  Eincident,e  Escattered,e(2.39)と表わされる．(2.36)式から電離後の電子の速度は，v~  v e 1 EEincident,e(2.40)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える．

28

(2.27)式中 v e に v~ を用いて，衝突後の速度 v e を求める．散乱角には(2.33)式を用いるが，(2.33)式中の Eincident,e には(2.38)式から求められる Escattered,e を用いる．電離エネルギー Eionization は，キセノンの場合 12.13 eV である．
電離後の生成電子の持つエネルギーは，Ecreated,e  E  Eionization(2.41)となる．また，入射電子の速度 v e を用いて生成電子の速度はv created,e veve

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

2.2.4.6.

Null-collision 法 (2-10)(2-11)(2-12)

計算時間短縮のため，Null-collision 法を用いる．これには，架空の衝突
断面積  fake を導入する．すべての電子のエネルギー E に対して
 fake   total(2.43)であるとし，例えば電子と中性粒子の衝突においては total   elastic   excitation   ionization   fake(2.44)とする．ここで， totalv collisionn gas

 constant(2.45)

を満たす  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-3 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度で規格化した衝突周波数を示す．本計算では Null-Collision 法
nで用いた collision gas の値は 6.5×10 -13 m 3/sec である．
29

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

2.2.4.7.

タイムステップの制約

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

P i2
1- P i(2.46)であるので，衝突ミスの回数を１％いかにするには𝑃 ≦ 0.1とする必要があ≦ 0.1 すなわち

る．nullcollision を用いる場合は，𝑃
Dt £2.2.4.8.

0.1

n null(2.47)

乱数生成

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

30

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
t

(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 
t
m
2






(2.49)

である． n-1 2 はある時刻 n から t 2 だけ以前の値であることを示す．同様に， n  1 2 はある時刻 n から t 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)と書ける．ここで， a ・  ・  ・  を次のように定義する．



qt n
Bx
2m

(2.51)



qt n
By
2m

(2.52)

 

qt n
Bz
2m

(2.53)

  1  2   2   2(2.54)

なお，ここでの E ， B は粒子位置におけるそれぞれの値を表わしている．
31

これらを用いて行列 A とベクトル S は，以下のように構成される．
 1   


A 
1  
(2.55)
  
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.56)

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

v n+1 2 = A -1 S 1  21
A 1        

2.3.1.2.

      
1   2       1   2 

(2.57)(2.58)

PIC 法重みつけ計算

3 次元における重みつけ方法を示す．粒子の位置から近傍のセル座標(i , j ,k )を求める．次に，図 2-4 に示す dx1,dx2,dy1,dy 2,dz1,dz 2 を求める．その値を用いて，図中にあるような直方体 V1〜V8 の体積を求める．この体積比から粒子が感じる電磁場を内挿することができる．
具体的には，座標 (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)と表わされる．

32(2.60)

図 2-4 PIC 法におけるセルの構成(3 次元)

2.3.1.3.

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

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

V1
qN
V2

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

33(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 つの観測時間をマイクロ波の１周期とし，タイムステップの間つまり運動方程式を解いた前後で変化したエネルギー増分をカウントし，エネルギー増分のヒストグラムとして，確率密度分布とする．マイクロ波
解析は計算コストの兼ね合いとある程度平均化された結果がでる 80 周期分を計算しているが，最初の 40 周期分は計算上の誤差が大きいため，40周期から 80 周期分の時間平均をとっている．また，各エネルギー増分はセル毎で観測されており，粒子の位置情報をセルに割り当てている． ヒス
トグラムのステップ数は，0.1eV，エネルギー範囲としては，-50〜50eVとした．

2.3.2.

電磁波解析手法 (2-13)(2-14)

2.3.2.1.

基礎式

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

34

H
t

(2.64)

E
J
t(2.65)

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

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

2.3.2.2.

差分式

三次元の FDTD 法で用いられる電場，磁場の空間配置を図 2-5 に示す．
このように，電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．これは Yee によって提案され (2-17) ，電場の回転（   E ）が磁場を，磁場の回転（   H ）が電場の方程式を満たすように配置されている．
この電磁場の配置のもと， (2.64)式, (2.65)式を空間，時間とも中心差分した式を以下に示す．((2.66)~(2.71))式ここで， x, y, z は空間幅， t は
時間幅，添え字 i, j, k  はセル中での座標，添え字 n は時間ステップを表わ
す．
n+

1

n-

1

H x 2 (i, j +1 2, k +1 2) = H x 2 (i, j +1 2, k +1 2) +Dt
m 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) üï
+
í
ý
Dz
Dy
ïî
ïþ
n+

1

n-

1

H y 2 (i +1 2, j, k +1 2) = H y 2 (i +1 2, j, k +1 2) +Dt
m 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) üï
+
í
ý
Dx
Dz
ïî
ïþ
n+

1

n-

1

H z 2 (i +1 2, j +1 2, k ) = H z 2 (i +1 2, j +1 2, k ) +Dt
m 0 (i +1 2, j +1 2, k )ìï Exn (i +1 2, j +1, k ) - Exn (i +1 2, j, k ) Eyn (i, j +1 2, k ) - Eyn (i +1, j +1 2, k ) üï
+
í
ý
Dy
Dx
ïî
ïþ

35

(2.66)(2.67)(2.68)

E xn 1 i  1 2 , j , k   E xn i  1 2 , j , k  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



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  z(2.69) 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



(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 
x



t
E zn 1 i , j , k  1 2   E zn i , j , k  1 2  
 0 i , j , k  1 2

1n
 n  122
 H y i  1 2 , j , k  1 2  H y i  1 2 , j , k  1 2

x



(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
y


1

n

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

36

図 2-5 FDTD 法における電場と磁場の配置

2.3.2.3.

アンテナ解析

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

𝑉
𝛥𝑥
𝑟𝑙𝑛 ( 𝑎 )37

(2.72)

𝐸 =で与えられる. ただし,𝑉𝛥𝑦𝑟𝑙𝑛 ( 𝑎 )(2.73)

V は中心導体と外部導体間の電圧, 𝑟 は動径方向の距離, 𝑎 はアンテナの半径である.
また, 同軸ケーブルのインピーダンス Z は以下の式で表される.
𝜇 1 𝑅
ln
𝜀 2𝜋 𝑟𝑍=√(2.74)

ここで, R は外部導体の半径, r は中心導体の半径である. マイクロは電力を P とすると，
𝑃=𝑉
2𝑍(2.75)となり，(2.74)式，(2.75)式より， V を求め，(2.72)式，(2.73)式よりアン
テ ナ に 与 え る 電 場 を 求 め る ． こ こ で ， Ex を 𝑟 = 𝛥𝑥/2の 電 場 で ， Ey を
𝑟 = 𝛥𝑦/2 の電場で代表させると，

Ex (IFED, JFED, KFED) = -E x (IFED -1, JFED, KFED)
V (nDt)=Dx æ Dx ö
ln ç ÷
2 è a ø(2.76)Ey (IFED, JFED, KFED) = -E y (IFED, JFED -1, KFED)=-

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

(2.77)

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

38

図 2-6 同軸線路給電

図 2-7 アンテナモデル

2.3.2.4.

PEC(Perfect Electric Conductor )境界条件

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

この境界条件を適用する．

2.3.2.5.

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

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

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

(2.78)と表わされ，次の微分方程式を満足する．

E z 1 E z0
x v t(2.79)

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

1
(2.80)n¶E z 2=v¶xとなる．右辺の x についての微分は，電場が x = 0, Dx に割り当てられていることから， x = Dx 2 で行い，(2.79)式に代入すると

E i, n  1 2  E
t
n
r

n1
r

n

1
2

n

1
2

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

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

(2.82)となる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場成分について適用される．

40

2.3.2.6.

ポインティングベクトル

空間を流れる電力を求めるためには，ポインティングベクトルを使用する．Maxwell の方程式を変形すると，以下の式となる．

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

(2.83)

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

2.3.2.7.

投入電場の補正 (2-26)

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

𝐸
ここで，𝐸

,𝐸

=𝐸

+

𝑃
𝜃( 𝑃

− 1) 𝐸(2.84)

2

は，更新前後の励振電場強度，𝑃

波 1 周期で平均した投入マイクロ波電力，𝑃
力とする．

41

は更新前のマイクロは設定投入マイクロ波電2.3.3.

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

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

42

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

43

2.4.

Full-PIC-MCC

Full-PIC-MCC は，解析領域を 2 次元円筒座標系とし，イオンと電子を共に粒子として扱う．時間スケールの異なるイオンと電子を共に粒子として扱うのは，計算コストが非常に高くあまり主流ではないが，その代わりにポアソン方程式を解くことが可能になり，静電場解析を可能にしている．
このモデルの対象とする物理現象は，シースの解析であるため，メッシュ
サイズがデバイ長に近い非常に細かいものとなっている．本来であれば，マイクロ波放電式中和器は，3 次元の非対称性をもつものであるが，モデル簡易化のために，2 次元円筒座標系とした．メッシュからくるタイムス
テップの制約の緩和のために，後述するが電子の質量比の変更を行い，タ
イムステップをイオンスケールに近づける工夫をしている．また，衝突計
算には，MCC 法を用いている．以下に詳しく述べる．

2.4.1.

運動方程式解法

基礎式は，2.2.1 で述べた．(2.1.1)式を 2 次元円筒座標系で，書き下すと以下のようになる．

v2
dvr q Er  v Bz   
dt
m
r
dv q
vv
 E  v z Br  vr Bz   r dtmr
dv z q E z  v Br 
dtm

(2.85)

これらを四次精度 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 
6vzn 1  vzn tbz1  2bz 2  2bz 3  bz 4 
6

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

44(2.86)qEr  v Bz   1 v v
mr
tb1   1 tb1 tb1 q
 Bz    v  v br 2   Er   v 
m
2
r
2
2

 




br1 

tb 2   1 
tb 2 tb 2 q Bz    v  v br 3   Er   v m
2   r
2 
2 

q
1
br 4  Er  v  tb 3 Bz  v  tb 3 v  tb 3 mr(2.87)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 

(2.88)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 

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 

(2.89)

tb 2  q Br 
bz 3   E z   v m
2  
bz 4 q
Ez  v  tb 3 Br m

また，本計算において，モデル簡易化のために，𝐸 = 0, 𝐵 = 0としている．
尚，r=0 での扱いは，𝑣 =0 として，

𝑣 =

𝐸,

𝑣 =

とした．

45

𝐸

(2.90)

2.4.2.

PIC 法重みつけ計算

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

An 

Sn
S

n  1 ～ 4

(2.91)となり，

Er, z   A1Ei, j   A2 Ei  1, j   A3Ei  1, j  1  A4 Ei, j  1(2.92)と表わされる．

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

2.4.2.1.

電荷密度計算

電荷密度は図 2-9 に示されるセルの格子点上で定義され，セル中に N 個

46

の粒子を代表した 1 個の超粒子が存在したとき，格子点 i, j  上に割り振られる電荷密度は， i, j  

S1 qN
S V

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

2.4.2.2.

電子質量の変更 (2-18)

今回の計算において，電子の速度成分は最大で 10 5 m/s 程度，イオンにおいては最大 10 3 m/s 程度であり，電子の速度を基準に時間刻みを設定するとイオンはほとんど動かない．そこで，電子の質量を実際のものよりα倍大きくした．質量の変更により，電子の速度は √𝛼分の 1 に低下，タイム
ステップを √𝛼 倍に伸ばすことが出来る．イオンとの速度差は小さくなり計算結果の収束を早めることができる．
me _ virtual  me _ real

(2.94)

質量を変更後の値を仮想値，変更前を現実値と呼ぶ．電子の質量変更によって各物理現象に矛盾を引き起こさないためにも各物理値を仮想値に修正する必要がある．考慮すべきはラーマー運動への影響から磁場，衝突計算への影響から衝突断面積を修正する必要がある．以下の表に現実値と仮想値をしめす．
表 2-2 質量比の変更による補正係数
現実値

仮想値

𝛼
𝑣

𝑣

√𝛼

Δt

√𝛼 Δt

B

√𝛼 𝐵

𝜍

√𝛼 𝜍

47

2.4.2.3.

粒子の境界条件

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

2.4.2.4.

Digital Filter

PIC 法で格子点上の量を決める場合, 例えば粒子の位置を格子上に割り振る場合には, メッシュ内に存在する粒子の数が尐ない領域ではその付近の格子に比べて数密度が小さくなり, 数密度の変化が不連続になる可能性がある. このような現象を避けるために Digital Filter と呼ばれる一種の補間を施すことがある. Digital Filter は, 格子上に定義される値, 例え
ば電場, 数密度等を計算する際に発生する高周波の雑音を除去する目的にも用いられる. (2-24)Digital Filter の式は以下のように表される.
new
Φi,j,k←side
cornerMΦi,j,k + SΦi,j,k+ KΦi,j,k

M + 6S + 12K(2.95)

new
Φi,j,k を格子点上に定義される値とし,Φi,j,k
は Digital Filter を通した新side
cornerし い Φi,j,k を表す. またΦi,j,kとΦi,j,kは side term と corner term であり, そ

れぞれ以下の式で表される.

side
Φi,j,k= Φi−1,j,k + Φi+1,j,k + Φi,j−1,k + Φi,j+1,k(2.96)+ Φi,j,k−1 + Φi,j,k+1
cornerΦi,j,k= Φi,j−1,k−1 + Φi,j+1,k−1 + Φi,j+1,k+1 + Φi,j−1,k+1+ Φi−1,j,k−1 + Φi−1,j,k+1 + Φi+1,j,k+1 + Φi+1,j,k−1(2.97)+ Φi−1,j−1,k + Φi−1,j+1,k + Φi+1,j+1,k + Φi+1,j−1,k
ここで,M,S,K の値をそれぞれ 12,4,1 とすると(3.44)式は次式となる.
new
Φi,j,k←1
1 side 1 corner
Φi,j,k +Φ+ Φ
412 i,j,k 48 i,j,k(2.98)

電子密度およびイオン密度，電位を計算した後に,それぞれの値にこのDigital Filter を作用させる.

48

2.4.3.

電場解析手法 (2-18)

2.4.3.1.

基礎式

ポアソン方程式および電場と電位の関係式を基礎式として計算を行った．

re0(2.99)

E = -Ñj

(2.100)Ñ 2j = -(2.99)式を差分化し,連立方程式として解いた．差分式を以下に示す．

aa1  i 1, j  aa 2  i 1, j  aa3  i , j 1  aa 4  i , j 1  aa5  i , j 
aa1  aa 2 

i , j
0

1
z 2

aa3 

1
1

2
r
2ri , j r

aa 4 

11
2
r
2ri , j r

(2.101)æ 1
1 ö
aa5 = -2 ç 2 + 2 ÷è Dz Dr ø
ここで， i，j はそれぞれ z，r 方向の格子点番号に対応している．𝛷 ,𝑗 は i，j における電位を表し，．𝜌 ,𝑗 はその点での電荷密度を表す．𝛥 𝑧，𝛥 𝑟は各方向の格子サイズであり，𝑟 ,𝑗 はその格子点の r 方向における位置を表す．
また，r = 0 上での係数は，
aa1  aa 2 1
z 2

aa3 = 0
aa4 =

4
Dr 2

æ 1
1 ö
aa5 = -3ç 2 + 2 ÷è Dz Dr øとする．(2-19)(2.100)式より静中心差分すると，

49(2.102)

E 1 =i+ , j
2

E

i, j+

1
2

=-

(j i+1, j - j i-1, j )2D(j i, j+1 - j i, j-1 )(2.103)

2Dとなり，静電場を求める．

2.4.3.2.

反復法と PETSc(2.99)式はすべての格子点に１つずつ与えられるので，連立方程式として解くことができる．I 番目の方程式が以下のように表せる連立方程式を考える．
na x  b
j 1

i, j

j

i

(2.104)x は変数，a，b は定数を表す．この時，x i 以外の変数が固定されているという仮定のもと，次の式から x i を求める．

 bi   ai , j x j 


j i


xi 
ai , j(2.105)

これは，次の式によって定義される反復法を示している．

 bi   ai , j x j ( k 1) 


j i
(k )

xi  
ai , j(2.106)

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

xi

(k )

 xi( k )  1   xi( k 1)(2.107)

ωを反復解の収束率を加速させるような値に選択することにより，ガウ
ス-ザイデル法より早い収束を得ることができる．ωには最適値が存在する
が，最適値予め知ることは一般的には不可能である．そのため，しばしば
ω=2-O(h)という発見的な評価式が使用される．ここで h は元の物理領域を離散化するときのメッシュ間隔である．
従来は，SOR 法にマルチグリッド法 (2-21) を適用したポアソンソルバを使用
してコードの高速化を図っていたが，並列化に不利である難点を抱えていた．今回，粒子コードの並列化を適用し，粒子計算における計算コストの削減を実現出来たが，ポアソンソルバの計算時間全体に占める割合が高く(全体の 7 割程度)，定常状態までの解析を不可能にしていた．このポアソ
ンソルバを最適化するためには，数多くの反復法を試行するまたは，並列化に対応させる必要があった．しかしながら，学習コストと実装コストの観点から 1 つ 1 つ手法を試してくのは，時間コストに対する成果が見られないと考えた．そこで，数多くの反復解法と前処理と呼ばれる反復計算をする前に行列式を適切な形に加工する方法(SOR も前処理の一つである)の組み合わせを容易に試行できるライブラリの利用を決定した．
今回利用する反復解法のライブラリは PETSc(The Portable Extensible
Toolkit for Scientific Computation)である．アメリカの Argonne National
Laboratory 内に発足した Mathematics and Computer Science Divisionproject のメンバーである Satish Balay, William Gropp, Lois Curfman
McInnes, Barry Smith の手によって開発されているプログラムライブラ
リであり現在アメリカでは航空力学や音響学等の数値解析に利用されている.

PETSc の特徴として以下に示す．

A) PETSc はオープンソースでありソースコードは, 以下の URL から入手できる.
http://www-fp.mcs.anl.gov/PETSc/
B) C, C++, FORTRAN の 3 種類のプログラミング言語をサポートしている.（それぞれの PETSc 関数モジュール自体は C 言語でプログラミングされているが, その関数モジュールをつなげて作成するプログ
ラムは上記 3 種類の言語のいずれでもよい. ）C) MPI（Message Passing Interface）を用いて, ベクトルの分散及び収集等といった並列化ベクトル計算や. 行列の組み合わせを簡単かつ効51

果的に行うことが可能である.
D)

疎行列のデータの格納をいくつかのフォーマットで行うことが可能である.

E) LU 前処理や ILU 前処理, block Jacobi 前処理などといった様々な
前処理プログラムを含んでいる. (前処理とは，前述のとおりはある問
題を特定の解法で解くとき, 予めより効果的に解けるような同値な問
題に変形することであり Krylov 部分空間法と併用することによって,収束性が著しく向上することが多く, 超大規模な問題の有効な解法として注目を浴びている.)F) GMRES, CG 法等といった様々な種類の Krylov 部分空間法（反復法）
プログラムを含んでいる.
G) 線形方程式・非線形方程式両方の solver を所有している.
利用方法としては，(2.101) (2.102)で示した係数行列と電荷密度分布を
PETSc の用意するデータ構造体に格納し，実行オプションとして，反復法と前処理を指定するものである．今回，選択した反復解法として，
Bi-CGStab 法，前処理法として，Geometric-Algebraic Multigrid 法を選択した．尚，計算アルゴリズムについては，今回記載しない．数種の反復解法と前処理の組み合わせを現モデルで試行し，安定解が得られる組み合
わせの中で，最も収束性の高いものを選択しただけである．

2.4.3.3.

電位の境界条件

反 復 法 の 境 界 上 兼 と し て , 大 別 し て 2 種 類 あ る ． Dirichlet 条 件 とNeumann 条件である．Dirichlet 条件とは，境界上の点の値を直接与える
ものである．一方，Neumann 条件とは，計算領域は有限であるが，その外にも広がっていると考えるのが妥当である．この境界外との連続性を保つために，その境界上の点における微分を直接与えるものである．今回，壁面に対して，Dirichlet 条件として電位 0 とし，r=0 の中心軸上ではノイマン条件として軸対称性から電位の傾きが零となるように境界条件を設定した．中和器は，イオンスラスタから排出されるイオンビームとのポテン
シャル差によって電子電流として引き出される．実験においては，中和器の排出口にコレクターと呼ばれるアノードをおき， これに数 10V の正電位をかけ，イオンビームとのポテンシャル差を模擬している．この排出出52

口，後述の計算体系図では，Plume に相当する部分に Dirichlet 条件として，𝛷𝑝をかけることで，数値解析により，中和器性能を評価できるものとしている．

2.4.4.

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

図 2-10 に Full-PIC-MCC のフローチャートを示す．静磁場およびエネルギー獲得確率は事前に計算しておく．初期化処理の際，粒子のデータは前回の情報を引き継ぎ，そのままリスタートさせる．ポアソンソルバにおける行列式を事前に作成しておき，各種 PETSc の初期化処理を行う．各格子点に電子およびイオンを重みつけにより，割り当て．電子密度，イオン密度を更新．電子密度，イオン密度より，電荷密度を計算し，その情報を PETSc に認識できる形に加工して，PETSc によりポアソン方程式を解
く．収束後，電位分布より静電場を算出後，その結果を用いて，PIC 計算を行う．電子に関しては，運動方程式を解く前に，エネルギー獲得確率よりエネルギーを供給しておく．PIC 計算後，電子の衝突計算を行う．

53

図 2-10 Full-PIC-MCC のフローチャート

54

2.5.

弱連成

FDTD-PIC-MCC と Full−PIC-MCC の弱連成は，繰り返し述べるが，
FDTD−PIC-MCC から Full−PIC-MCC にはエネルギー獲得確率を Full−
PIC-MCC から FDTD−PIC-MCC には電子密度及び電子エネルギーである．
弱連成のタイミングであるが，FDTD−PIC-MCC の計算時間は，マイクロ波周期によって決定し，約 50〜80 周期分を取れば， 数値的なばらつきはないものとした． 本計算では，約 75 周期分の 30ns を解析している．Full−PIC-MCC の 解 析 時 間 で あ る が ， 密 度 が 大 き く 変 わ ら な い 限 り Full−
PIC-MCC から FDTD−PIC-MCC のフィードバックは意味がないとしている．よって，ある程度電子とイオンを動かす必要があり，この時間を経験
的に 250ns 程度とした．連成のタイミングは，Full-PIC の計算の不安定性が高いという原因も有り，最適化されていないため，今後物理的な意味を持つ解析時間を選択する必要がある．

2.6.

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

図 2-11 に本研究室所有のマイクロ波放電式中和器と解析領域を示す．
これは， yz 断面で切り取った断面図である．ここにおいて，z 方向を図の左から右方向，y 方向を図の下から上方向としており,rz 断面図の際も r 方
向が y 軸の向きと一致している．形状としては，円筒径の放電室(Wall)を
箱状の永久磁石(サマリウムコバルト磁石)が z 軸方向に 3 個を 1 組として，
棒磁石として 15 組の磁石が周方向に均等に配置されている．磁石を鉄で出来た FrontYork および BackYork が挟んでおり，この間に磁力線を形成する形状となっている．BackYork から伸びているがアンテナであり，これは L 字型を取っており，座標系において，y 軸方向にアンテナの先が向いている．FDTD−PIC-MCC では 3 次元直交座標を取り，この図を周方向に全てとった形状となっており，(x，y，z ) = (20, 20, 18)[mm]である．一方，Full-PIC-MCC は，rz 円筒座標系を取り，(r, z) = (10 18)[mm]である．
この中和器内部の静磁場強度分布 x 中心断面図を図 2-12 に示す．これは ， Advanced Science Laboratory, Inc 社 の 3 次 元 磁 場 解 析 ソ フ ト
Amaze を用いて解析した磁場強度分布である． 今回使用するマイクロ波
周波数は 2.45GHz であるので，それに対応する ECR 磁場強度である 85mTを赤線で示している．赤丸で囲んだところに，磁気チューブと呼ばれる電55

子補足層があり，この領域を電子は，EｘB ドリフトまたは grad B ドリフ
トをしながら移動する.

図 2-11 本研究室所有のマイクロ波放電式中和器と解析領域

56

図 2-12 静磁場強度分布(赤線が ECR 層)

57

2.7.

各種物理条件と計算条件

2.7.1.物理条件
表 2-3 に本計算に用いた物理条件を示す． 本研究室の実験体系と同等のものを扱っている．コントロールパラメータとしては， 出口電位を 2 種類用意し，この違いにより，引き出し電流特性の違いや各種プラズマへの影響をみる．推進剤流量は計算上では，一定の中性粒子密度として扱われている．これは，別途 DSMC 計算にて，0.05 mg/s を本中和器に流した時，中性粒子密度がガス口を除き，中和器全体で約 2.0×1020 m-3 という結果が得られているため，衝突計算における中性粒子密度として，この値で一定とした．イオン初期エネルギーは室温として，初期密度初期電子エネル
ギーに関しては，過去の実験データより設定した． (2-27)
表 2-3 物理条件
推進剤

Xe

推進剤流量

0.05 mg/s

マイクロ波周波数

2.45 GHz

投入電力

8W

初期密度

1.0×1017 m -3

電子初期エネルギー

2.0 eV

イオン初期エネルギー

0.026eV(300K)
10. 40 V

出口電位 𝛷𝑝

2.7.2.

計算条件

FDTD−PIC-MCC
メッシュサイズを 0.5 mm とし，x 方向に 40 格子 y 方向に 40 格子 z 方
向に 36 格子とした．FDTD−PIC-MCC でみる物理対象には，シースの解析は不要である・マイクロ波波長は約 7cm 程度であり，本メッシュサイズ
で問題はないとされる．メッシュサイズ 0.5 mm より FDTD における計算
制約条件は，プラズマ周波数より， 0.2 𝜔𝑝−1 = 8.0 ×10-12

sec , 電磁波の伝搬速度を光速と仮定すると Courant 条件より，9.6 ×10-13

58

sec 以

下であればよく計算誤差を考え，2.5 ×10-13

sec とした．PIC におけるsec , 衝突

制約条件としては，サイクロトロン周波数より，2.0 ×10-12
計算の制約より，7.3×10-10

sec，最大電子速度を 50eV とした時の 1 セ

ルあたりを飛び越えないための Courant 条件は，6.9×10-11
さらにその 10 分の 1 程度とすると 7×10-12

sec となり，sec となる．このなかで最も制約条件の厳しいのは，Courant 条件であるためこれを満たしさえすればよく．FDTD10 回につき PIC を 1 回解くとすると，2.5×10-12

sec で十分計算の制約を満たす．
表 2-4 計算条件(FDTD−PIC-MCC)
メッシュサイズ

0.5 mmタイムステップ(FDTD)

2.5 ×10 -13 sec

タイムステップ(PIC)

2.5 ×10 -12 sec

超粒子の重み9.85 ×10 4

超粒子数

400 万個

1 セルあたりの粒子数

130 個

また，FDTD の境界条件であるが，体系図の Plume 面及びアンテナの根本のみを Mur の１次吸収境界をかけている．
Full-PIC-MCC
Full−PIC-MCC の解析対象は，シース解析であるために，メッシュサイ
ズがデバイ長以下であるのは，必須である．デバイ長は電子密度 2.0×1017
m-3 平均電子エネルギー8eV 程度とすると，0.047 mm 程度である．そのために計算精度を考え，0.025mm とした．この際の Courant 条件は，最大電子エネルギーを 50 eV と仮定すると，4.0×10-12 sec となる．さらに
その１０分の１を考えると，

4.0×10-12 sec，最小値は 4.0×10-13sec

であるが，タイムステップを尐しでも広く取るのために，Courant の８分の１程度として，5.0×10-13 sec とした．これを仮想値に変換すると，2.5
×10-11 sec つまり 0.025ns とした．

59

表 2-5 計算条件(Fll−PIC-MCC)
メッシュサイズ0.025 mm

タイムステップ(現実値)

5.0 ×10 -13 sec

タイムステップ(仮想値)

2.5 ×10 -11 sec

超粒子の重み3.90 ×10 4

超粒子数(イオン・電子)

1000 万個

1 セルあたりの粒子数

50 個

参考文献
(2-1)

笹川裕太郎，”数値解析によるマイクロ波放電式中和器の性能評

価”，平成 25 年度九州大学修士論文(2-2)

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

(2-3)

田中實, 山本良一：“計算物理学と計算化学―分子動力学法とモン

テカルロ法”(海文堂,1998)(2-4)

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

Giovanni

Lapenta,

JeremiahU.

Brackbill,

JournalofComputational Physics, 115, (1994), 213.
(2-6)

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

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

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

南部健一：“コンピュータアナリシスシリーズ 7を用いる数値シミュレーション(2-10)

V.

Vahediand

原子･分子モデル

日本機会学会編” (コロナ社, 1996)

M.

Surendra,

ComputerPhysics

Communications, 87 (1995), 179.
(2-11)

電気学会・マイクロ波プラズマ調査専門委員会：“マイクロ波プラ

ズマの技術” (オーム社, 2003)(2-12)

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

(2-13)

橋本修, 安部琢美：“FDTD 時間領域差分法入門” (森北出版, 1996)(2-14)

宇野亨：“FDTD 法による電磁界およびアンテナ解析” (コロナ社,

1998)(2-15)K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.

60

(2-16)

G. Mur, IEEE Trans. Electromagnetic Compat., 40 (1998), 100.

(2-17)K. S. Lee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.

(2-18)

Szabo Jr, James Joseph. “Fully kinetic numerical modeling of aplasma thruster.” Diss. Massachusetts Institute of Technology,
2001.
(2-19)

田中守也，”電子・イオンビーム工学”(電気学会，1995)(2-20) Richard Barrett, Michael Berry, Tony F. Chan, James Demmel,
June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
Charies Romine, Henk va der Vorst, ：“反復法 Templates” (朝倉書店, 1996)(2-21)

廣池匠哉,“マルチグリッド法を用いたマイクロ波放電式中和器のプラズマ損失機構解明”，平成 24 年度九州大学修士論文．
(2-22) C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal ofChemical Physics, 55 (1971), 4100.
(2-23) 保原充，大宮司久明：”数値流体力学:基礎と応用”(東京大学出版会，
1992)(2-24)

Birdsall, C. K., and A. B. Langdon. "Plasma physics via

computer simulation." (1985).
(2-25)

G. M. Tuner, I.S. Falconer, B. W. James, and D. R. McKenzie, J.

Appl. Phys. B65B, (1989)(2-26)

渡 邊 裕 樹 , “ マ イ ク ロ 波 放 電 式 電 子 源 内 の 電 子 加 熱 解

析”. STEP-2013-014(2-27))小野寺範義, 竹ヶ原春貴, 西山和孝, 船木一幸, 國中均: マイク

ロ波放電型中和器の電子放出機構, 日本航空宇 宙学会論文集, 第 49
巻, 第 564 号, pp.27-31, 2001.

61

第3章
1.1.

結果と考察

FDTD−PIC-MCC

FDTD−PIC-MCC の節では，主に妥当性の検証と連成解析の前後での各物理値の変化の様子を考察していく．

3.1.1.

電子軌道

電子のラーマー運動の確認の為に，電子軌道をしめす．コード単体の検証のために弱連成は行わず，Cold 状態での観測を実施した．図 3-1 に初期位置，30ns 後の電子流線をしめす．左が x=20mm の x 中心 zy 断面，右が z=18mm の z 中心 xy 断面図である．以降空間分布においては，特筆しない限り，同じ断面を示す．
まず，流線全体図から磁場との相互作用により，ラーマー運動が確認され，
ExB ドリフトの効果もみることが出来る．アンテナからフロントヨーク(計算体系図参照のこと)にかけて形成されている磁気チューブににかけて，電子の補足が多くあると仮定されているが，電子が磁力線に沿って，捉えられるのを良く表している．フロントヨークにおいて，電子の反射がなされているのがわかるが，これは，粒子計算特有のもので，流体計算等では，
このような詳細な物理は表現できない．赤線の粒子は，初期位置がアンテ
ナの先端付近としているが，これは，磁気チューブでの運動よりアンテナの先端で捉えられ，強電界の効果を与えられている事がわかる．また，初期位置をアンテナの根本に配置した粒子は，まっすぐ Plume に流れていっており，これは，発散磁場に沿って，電子が流されていることを示している．アンテナ前面の密度及び磁場形状によって，中和器の性能が変わるということが言えるため，性能評価をする際に，1 つの指針であることは間違いないだろう．
以上より，ラーマー運動と磁気チューブでのトラップが表現出来たので，
これも解析モデルの妥当性の証の 1 つと言える．

62

(a) 電子軌道(初期位置)(b) 電子軌道(30ns)
図 3-1 電子軌道63

3.1.2.

空間電子密度分布・エネルギー分布・電界強度

分布
電子のエネルギー吸収は，密度変化にてどのように変化するかをみるた
めに，以降の空間電子密度分布・エネルギー分布更に，電界強度分布の連
成解析前の比較を行う．尚，連成解析に使用したのは，𝛷𝑝= 40 V の密度変化を導入している．
まず，連成前の FDTD−PIC-MCC の１サイクル計算にて図 3-2 の(a)に示す空間電子密度分布を見ると，電子の密度は初期の 1.0×10 17 m -3 からほとんど変化していないのがわかる．なぜならば，FDTD−PIC-MCC の時間
スケールは，電子スケールであるため，ほとんど密度の変化にはつながらない．密度が変化するまで計算することも可能であるが，Full−PIC-MCC
である程度動かして反映させるほうがより効率の良い計算サイクルになる．
(b)の連成後の密度分布を見てみると電子密度が空間的に増加している．全体的に，約 4.0×10 17 m -3 ほどでまで増加し，Plume 部分に 6.0×10 17 m -3 と比較的高い密度を持つように変化している．さて，この密度変化がどのよ
うにエネルギー獲得に影響をおよぼすかを実際に，電子エネルギー分布を見て考察してく．図 3-3 に電子エネルギー分布を示す．まず，1 サイクル
分の Cold 状態スタートした(a)の図を見てみると，磁気チューブを中心に電子が加速されて行き平均エネルギー8 eV〜12 eV が分布している．それは，アンテナの先端側の反対側つまり y=0〜10mm の領域にも同様に磁気チューブが形成されているが，アンテナ先端近傍よりは低いが，磁気チューブにて加速が起こっている事がわかる．また，平均エネルギーの最も高い領域は，アンテナの周辺に 15eV 程度が分布しており，特にバックヨー
クとの境界において，約 20eV もの高いエネルギーの電子が存在している．
これは，アンテナの強電界による物で，ECR 加速とは別の加速機構を保持している可能性を示している．その証拠に(b)を見てみる．密度変化があっ
た後，全体的にエネルギーが増加しているが，これは，エネルギー保存のために Full−PIC-MCC の電子エネルギーを Update したためであり，平均
して 8 eV 程度のエネルギーが広く分布している．特筆すべきは，アンテ
ナの先端の逆方向の磁気チューブにおいては，エネルギー吸収が行われておらず，密度上昇に従って,マイクロ波の伝搬モードが変化したことを示している. そして,肝心のアンテナ周辺であるが, 密度変化が起こった後もアンテナ周辺の加熱は変わらずに起こっており, 定常状態になり, 密度が増64

加した後も, 表皮深さの進入がおこり, プラズマ生成に重要であることを示している. 逆にいうと,もしアンテナの伝搬性能が落ちた場合, プラズマ生成に影響を及ぼし, 性能低下に繋がることは間違いない．アンテナの性能低下を評価できるシステムが中和器の性能評価につながると考えられる．
マイクロ波性能を示す指標として， 電界強度分布は重要であると考え，
図 3-4 に電界強度分布をしめす．まず，アンテナ近傍に強電界が形成されていることが, はっきりと分かる. マイクロ波の伝搬が正しく行われているのかを見てみると，yz 断面でも xy 断面でも波が現れており，アンテナには完全導体条件を課しているため，本来であればマイクロ波が漏れることはないのだが，波のようにマイクロ波の反射が多く起こっている事がわ
かる． この強電界により, 電子の加速に寄与していることが電子エネルギ
ー分布からも明らかであり, (b)を見てみると，壁面付近に，5×10 4 V/m 程
度の電場が形成されている．この原因として，壁面付近の電流密度項の増加によるものと考えられる．時間変化後の電子の動きとのカップリングがしっかり反映されていることがわかる.

65

(a) 連成前(b) 連成後
図 3-2 空間電子密度分布 連成後と連成前

66

(a) 連成前(b) 連成後
図 3-3 電子エネルギー分布67

(a) 連成前(b) 連成後
図 3-4 電界強度分布

68

3.2.

Full−PIC-MCC

Full−PIC—MCC の節では，Plume 領域にかける電位(10 V，40 V )によって，各物理値がどう変化するかを検証していく．

3.2.1.

空間エネルギー分布

図 3-5 に電子とイオンそれぞれの空間エネルギー分布を示す．まず，電子について見ると，両者ともに平均して，6 eV の電子温度が全体に広がっているのがわかる．本来であれば，アンテナの付近や磁気チューブ近傍に高エネルギーが存在しているはずであるが，局所的な高まりは見て取れな
い．電子のエネルギー吸収がうまくやり取りできていない事がわかる．電子エネルギー吸収のモデル化が出来てない．今後，この機構の見直しが必要である． 電子エネルギー獲得確率の算出方法の改善や，ヒストグラムのエネルギー幅をさらに細かくする等,数値的不安定性を極力抑えたモデル化が必用である．
電子のエネルギー吸収を抜きに考えると，𝛷𝑝=40 V においては，

Plume 領域にかけ，電子 10eV 程度のエネルギーを持っている．これは，静電場による加速がなされているものと考えられる．𝛷𝑝

=10 V では，

そのような静電加速は行われていない．静電加速の影響を解析できていることを表している．
イオンに関していうと，シース電位により加速された高エネルギーイオ
ンが壁面付近に存在している．そのエネルギーは，4〜5 eV である． 放電
室内部のエネルギーは初期温度である室温程度に保たれたままである．シース加速を再現出来ているといえる．さらに𝛷𝑝=40 V では，Plume

出口付近のオリフィスへ 15eV 以上のイオンが存在し，Plume 径への高エネルギーイオンの衝突が示唆され，この部位のスパッタが多くなされる可能性があると考えられる． IEDF を後述するので，その際に，詳しく述べる．

69

(a)𝛷𝑝

=40 V

(b)𝛷𝑝
=10 V
図 3-5 電子・イオンエネルギー分布(左図：電子，右図：イオン)

3.2.2.

空間密度分布電子・イオンの空間密度分布を図 3-6 に示す．まず，イオンと電子においては，ほとんど定性的，定量的差異は，ないものと考えられうる．軸上で，高密度が存在するが，これは，粒子法ゆえのセルあたりの粒子数が尐ないためにでる影響であり，物理的意味は持たないので取り上げる必用はない．𝛷𝑝

=40 V においては，Plume 領域のみに 5.0×10 17 m -3 ほどの高い密度を示しているが，𝛷𝑝=10 V においては, Plume 領域から中空に向けて広く 6.0×1017 m -3 ほどの密度が分付していることがわかる．プ
ラズマ生成よりも Plume への排出が多いために，40 V の方が軸上付近で低密度になると推測されるが，エネルギー授受の課題であるとも言えるので，今後調査が必用である．

70

(a)𝛷𝑝

=40 V

(b)𝛷𝑝
=10 V
図 3-6 電子・イオン密度分布(左図：電子，右図：イオン)

3.2.3.

空間ポテンシャル分布および静電場分布

図 3-7 に空間ポテンシャルおよび， 静電場のベクトル図を示す．𝛷𝑝
=40 V においては Plume に 40 V がかかっており，そこから放電室内部へ
ポテンシャル勾配が形成されている．静電場のベクトル図においても，
Plume から放電室内部へ電場のベクトルが向いており，この電場に加速されて，Plume 領域へ電子が排出されるはずである． また，Plume 出口の
r 軸方向の強い電場は，解析上壁面を 0V としているため，この壁面への電場が高くなると考えられる．この電場の加速によってイオンが加速されているものとされる． また，シース電位は，アンテナ付近およびオリフィス
プレートで，20V 程度，その他壁面に対して 10〜13V である．
𝛷𝑝

=10 V では，放電室内部が 15〜20V 程度であるために，ポテン

71シャル勾配が Plume 側に向いている．そのため，電子は放電室内部側に静電加速され，電子は排出されにくくなることが見て取れる．これにより内部プラズマのポテンシャルが引き出し電子電流の性能に影響を及ぼし，これを解析的に表現する事が可能になった．シースの厚みが 40V の時より，厚くなっているが，これは，引き出し電圧と内部の電位差が尐ないために，勾配が緩やかになったと考えられる．壁面電流・引き出し電子電流性能評
価の節で実際に引き出し電流がどの程度際があるかを述べるが，簡易ながら高電流モード，低電流モードのメカニズムを解析可能になったことは，中和器解析ツールとしての実用性を証明したことになるであろう．

(a)𝛷𝑝

=40 V

(b)𝛷𝑝
=10 V
図 3-7 空間ポテンシャル(左図)・静電場分布(右図)72

3.2.4.

壁面電流および引き出し電流

図 3-8 に電子・イオンの各壁面電流値および引き出し電流を示す．青が電子電流，赤がイオン電流を示す．まず，各壁面の電流バランスとして，
Wall と Orifice Plate では，イオン電流が支配的であり，これは実験定性的とに一致する．しかし，BackYork においては，イオン電流が支配的でAntenna では，40V の際に，イオン電流が支配的となっている．これは，実験とは一致しない．この理由として，プラズマ密度が未だ低く，1×10 18
m -3 オーダーまでは，密度の上昇は考えられるので，今後更なる解析を進める必要がある．また，アンテナでのエネルギー授受が過小評価されているため，アンテナのシースを乗り越えて，入っていくような電子が尐なく，アンテナでの電子電流の増加につながらないと考えられる．アンテナのシース電位の低い 10V の方は若干であるが，電子電流の方が支配的となっており，実験との一致を示すが，エネルギー授受が過小評価されていることは変わりがないので，評価できないであろう．
引き出し電子電流については，40V の時約７mA ，10V の時 2 mA となった． モードジャンプといえるほどの電流増加とは言いがたいが， ，壁面電位によって，引き出される電子電流が増加するという定性的な性質を
本解析モデルで初めて実現することができた．定量的に見ると，各電流値
は，実験値とワンオーダー低い．引き出し電流で言うと，120mA 程度と見積られている． それ故に，電流評価方法の改善とエネルギーの授受モデルの改善が必要である．

73

(a)𝛷𝑝

=40 V

(b)𝛷𝑝
=10 V
図 3-8 壁面電流・および引き出し電流743.2.5.

IEDF

最後に各壁面に対するイオンのエネルギー分布関数(IEDF)を示す．40Vの時は，オリフィスへのエネルギー確率が高くなっている．これは，シー
ス加速によるものが大きい．スパッタ粒子のしきい値は，30eV 以上であるため，今回の解析よりオリフィス部において，スパッタが起こる確率が高くなるものと推測される．しかし，1 価イオンであれば，各壁面に対するエネルギー確率はスパッタを起こすほどの粒子存在しないため，今後 2価イオンの考慮を入れた計算モデルの改善が必要となってくる．壁面において，十分スパッタを起こすほどの粒子がいると解析できれば，この結果
を DSMC によるスパッタ粒子輸送解析に使用でき，寿命評価に関して，
より信頼性のある結果を得られるであろう

75

(a)𝛷𝑝

=40 V

(b)𝛷𝑝

=10 V

76

第4章

結論

本研究ではマイクロ波放電式中和器の性能を統合的に解析出来るツールの開発を行った．
FDTD−PIC-MCC と Full-PIC-MCC の連成解析を可能にした.
密度変化によるマイクロ波吸収機構の変化を取り入れることが可能になり,
より定性的な議論が出来るようになった.
引き出し電子電流電圧特性の評価を行うことを可能にした．その結果として， 引き出し電圧 40V の時に 7 mA，引き出し電圧 10V の時に 2 mA と定性的にモード遷移を解析することができた．しかし，実験値とはワンオーダーほど異なるので，今後計算モデルの改善が必要である．

77

謝辞
先ず始めに, 学部時代より宇宙への興味を持たせてくれ, ロマンと夢を教えてくれたと同時に, 貴重な研究の場を提供していただきました中島秀紀教 授に 深 く感 謝 致 しま す . 時に 白 熱す る議 論に よ って , 研 究者 とし て度々成長をさせて頂き, 多くの勉強の機会を頂きました山本直嗣准教授に深く感謝致します. 連日の研究で疲れた時に, 明るく労って頂き, 研究室を明るく支えてくださるとともに学生生活の様々なサポートをしていただ
いた馬渡隆子氏に深く感謝致します.3 年間同じ研究室で苦楽を共にしてきた木邊厚視氏, しっかり者でいつも余裕のある行動を取っていた伊藤匠氏,
チャレンジ精神やバイタリティにあふれていた田中啓介氏, いつも笑い溢れる話題を提供してくれた坂井寿光氏, 冷静と情熱を併せ持った所が魅力的な斎藤直哉氏という尊敬できる同輩の皆様と一緒に研究生活を遅れたことを深く感謝致します. また, 研究室では楽しく充実した日々を与えてくださった M1 と B4 の皆様, D1 および中島研究室を卒業された諸先輩方に深く感謝致します.
最後に, ここまで学生生活を送らさせて頂き, 私が進む進路を応援して下さった家族に深く感謝致します.

78