目次
第 1 章 序論 ............................................................................................. 1
1.1. 背景 ...................................................................................................... 1
1.2. イオンエンジンの原理 ........................................................................... 2
1.2.1. マイクロ波放電式イオンエンジン ...................................................... 3
1.2.2. 電子サイクロトロン共鳴 .................................................................... 4
1.2.3. マイクロ波放電式中和器 .................................................................... 5
1.3. 研究目的 ............................................................................................... 8

第 2 章 数値解析手法 ............................................................................... 92.1. 解析全体のフロー .................................................................................. 92.2. プラズマ粒子計算手法 (13) ..................................................................... 11
2.2.1. 基礎式 ............................................................................................... 11
2.2.2. メッシュサイズ・タイムステップ (14) ................................................ 122.2.3. 粒子計算 PIC 法 (15) ........................................................................... 13
2.2.4. 衝突計算 MCC 法 ............................................................................. 13
2.3. DSMC ................................................................................................. 19
2.3.1. 中性粒子解析手法 ............................................................................. 20
2.3.2. 中性粒子の運動 ................................................................................ 20
2.3.3. 中性粒子の初期条件と境界条件 ........................................................ 20
2.3.4. 中性粒子間の衝突 ............................................................................. 21
2.4. FDTD-PIC-MCC .................................................................................. 22
2.4.1. プラズマ粒子計算 ............................................................................. 22
2.4.2. 電磁波解析手法 (26)(27) ........................................................................ 25
2.4.3. FDTD-PIC-MCC のフローチャート .................................................. 32
2.5. Full-PIC-MCC ..................................................................................... 33
2.5.1. 運動方程式解法 ................................................................................ 34
2.5.2. PIC 法重みつけ計算 .......................................................................... 35
2.5.3. 電場解析手法 (31) ............................................................................... 38
2.5.4. Full-PIC-MCC のフローチャート ..................................................... 40
2.6. 計算体系および静磁場分布 .................................................................. 41
2.7. 弱連成 ................................................................................................. 43
2.8. 各種物理条件と計算条件 ...................................................................... 44
2.8.1. 物理条件 .......................................................................................... 44
2.8.2. 計算条件 .......................................................................................... 44

第 3 章 結果と考察 ................................................................................ 46
3.1. DSMC ................................................................................................. 46
3.2. 各物理値分布 ...................................................................................... 47
3.2.1. 電子エネルギー分布 ......................................................................... 47
3.2.2. 空間ポテンシャル分布 ...................................................................... 48
3.2.3. 空間密度分布.................................................................................... 51
3.3. FDTD-PIC-MCC の計算不安定性 ......................................................... 53
3.4. 電子・イオン電流 ................................................................................ 53
3.4.1. 各部材の電子・イオン電流 ............................................................... 53
3.4.2. 引き出し電子電流 ............................................................................. 56

第 4 章 結論 ........................................................................................... 58
謝辞 ....................................................................................................... 59
参考文献 ................................................................................................ 60

第1章
1.1.

序論

背景

2010 年，小惑星探査機「はやぶさ」が小惑星イトカワからサンプルを採取後，地球への帰還に成功した．この「はやぶさ」のミッションには，長時間の安定した動作をする動力源が必要不可欠であり，それを実現したのがマイクロ波放電式イオンエンジンμ10 である．
マイクロ波放電式イオンエンジンとは, 電気推進の一種であり, 電気エネル
ギーを用いて推進剤を電離させ，生成されたプラズマを様々な方法で加速し，その反作用で推力を得る (1)．従来の化学推進機に比べ，電気推進機は得られる
推力ははるかに小さいが，燃費を示す指標である比推力(Specific impulse：Isp)
は一桁以上大きくすることが可能である．つまり，電気推進機は加速に時間がかかるが，消費する燃料は化学推進機と比較し非常に尐なくて済むということである．電気推進機は，イオンエンジンやアークジェット，MPD スラスタ，
RF スラスタ，パルスプラズマスラスタ，ホールスラスタなどが開発・研究されている．その中でも，イオンエンジンは他の電気推進機と比較しても比推力が高く(Isp 数 1000s 以上)，長期間のミッションである深宇宙探査に最適である．
この中でもマイクロ波放電式イオンエンジンは，プラズマ生成にマイクロ波を使用し, 長寿命が特徴のイオンエンジンである．JAXA の宇宙科学研究所では18000 時間（1997 年開始）と 20000 時間（2000 年開始）の耐久試験に成功しており，はやぶさのミッションで用いられた 4 台を一つのシステムとして構成されたμ10 は，積算作動時間 40000 時間を達成し，高い信頼性を実証した (2)．
2014 年には，はやぶさの後継機「はやぶさ 2」が打ち上げられ，地球からの直線距離がイトカワの約 2 倍で，約 3 億 km の小惑星リュウグウからのサンプ
ルリターンの実現を目標としている (3)．はやぶさ 2 で使用されるイオンエンジ
ンμ10 は推力が，はやぶさの 8 mN から 10 mN と 25%アップしている．
この「はやぶさ」，「はやぶさ 2」の活躍を契機に，マイクロ波放電式イオンエンジンは信頼性や実用性が高く評価され，超小型衛星の市場活性化に伴うマ
イクロ波放電式イオンエンジンの需要は年々増加している (4)．その使用例の一つとして，宇宙デブリ問題が挙げられる．現在，宇宙デブリは 1cm 以下の小型のものを含めると数百万個と言われており，今後宇宙開発を行っていく上で解1

決する必要のある問題である．宇宙デブリ問題には，様々な解決策が提案されているが，その中の一つにイオンエンジンを用いた小型衛星によるものがある(5)．デブリ診断用機器を搭載した母艦と，デブリ除去機能を持った子機から構成されており，低コストかつ短期開発でデブリ除去衛星を打ち上げることが可能である．また，その他にも超低高度衛星等のミッションでのイオンエンジンの需要も高まっており，観測分解能の向上やセンサの小型化・低消費電力化が図られ，災害監視，地球環境観測，サイエンスなどの様々な分野への応用が期待されている (6)．上記のようなミッションの幅を広げるためにも，イオンエン
ジン自体の性能を向上させることは非常に重要となってくると言える．
さて，上述したマイクロ波放電式エンジンに対する需要に応えるためには，現状抱えているある問題点を解決する必要がある．それは中和器の性能务化である．中和器とは，電子を放出することでイオンを中和する装置である．中和器を使用しない場合，宇宙機自体が負に帯電し，排出したイオンが逆流するため，十分な推力が得られない．はやぶさにおけるマイクロ波放電式イオンエン
ジンμ10はミッション中，中和器に問題が発生した．μ10は単独で 14000 時間以上という最長作動実績を上げたが，中和器は务化の兆候を見せた．4 機搭載された内の 1 機はイオン源の点火不良，3 機は中和器が想定より早く寿命を迎え，運用停止の機器に瀕することになった．バックアップ回路により，残る 1 機の中和器とイオン源を組み合わせたクロス運転を行うことで，地球帰還を果たした．これによりマイクロ波放電式イオンエンジンシステムのロバスト性の高さが証明されたが，中和器の寿命が不安要素であるのは明らかとなった．イオン
エンジンが長期のミッションに適しているからこそ，マイクロ波放電式イオンエンジンの長寿命化は必須であり，そのためには中和器の性能向上が重要なものとなっている．

1.2.

イオンエンジンの原理

ここではまず，イオンエンジンの原理を通して，本研究で解析を行う中和器の位置づけおよび役割を説明する．イオンエンジンの概略図を図 1-1 に示す．
イオンエンジンはアーク放電やマイクロ波などで推進剤を加熱・電離させてプラズマを生成し，2 枚もしくは 3 枚からなる多孔状の電極に 1000 V 程度の高電圧を印加させ，イオンを加速させる静電加速型の推進装置である．イオンエン
ジンは以下の 3 つの領域から構成されている．

2

① 推進剤を電離させるイオン生成部（Ionization）② 生成されたイオンを静電的に加速して推力を得る加速部（Acceleration）③ 放出されたイオンビームを電気的に中和する中和部（Neutralization）
これらの各過程はそれぞれ，イオン源，加速電極，中和器によって行われる．

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

1.2.1.

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

プラズマの生成方法の違いにより，イオンエンジンは直流放電式，RF (Radio
Frequency) 誘導放電式，マイクロ波放電式等に大別される．本研究の解析対象であるマイクロ波放電式イオンエンジンの概念図を図 1-2 に示す．マイクロ波
放電式イオンエンジンは，マイクロ波帯域の交流電場により電子を加速させ中性粒子と衝突電離によりプラズマの生成を行う．

3

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

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

このときの円軌道の半径 rL はラーマ半径と呼ばれ，次式で与えられる．

rL 

mv
eB

(1.2.2)

ここで v⊥ は電子の B に対する垂直な速度成分の大きさである．プラズマを構4

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

 ce 

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

E

電子

B
E

B

電子

電子
E

B
B

電子
E

図 1-3 ECR の原理

1.2.3.

マイクロ波放電式中和器

マイクロ波放電式イオンエンジンには，そのイオンビームに同量の電子を供給
し電気的に中和する中和器という装置が搭載されている．これにより，マイク5

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

6

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

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

7

1.3.

研究目的

マイクロ波放電式中和器の性能向上に向け，本研究室ではこれまで数値解析
によるアプローチを試みてきた．その理由は，実験による内部物理計測ではマ
イクロ波放電式中和器の空間スケールが小さいため，プローブ等の測定機器を掃引するとプラズマの擾乱を引き起こしてしまい，正しい結果が得られにくい
こと．更に，何万時間に及ぶ地上試験を行うのは，時間的・金銭的コストが莫大になるということなどが挙げられる．数値解析であれば，これらの問題を解消し，中和器内部の物理現象を捉えることが可能である．以上より数値解析は有用であるといえる．
イオンエンジンの数値解析による研究は多くの研究機関で実施されているが，中性粒子解析，マイクロ波解析，プラズマ解析，損耗解析の異なる時間スケー
ルの物理現象を統合的に解析するツールの開発は行われていない．本研究室では，笹川により，中性粒子解析と損耗解析を行うツールである DSMC の開発が行われた (11)．また，楠田により後述する FDTD-PIC-MCC と Full-PIC-MCC を連成的に解析するツールの開発が行われ，マイクロ波解析，プラズマ解析を行うことが可能となった (12)．しかし，この連成解析ツールは不完全であり，中和器より引き出される電子電流の値も実験値の 1/100 程度と極めて小さい値となっている．
そ こ で ， 本 研 究 で は 連 成 解 析 ツ ー ル の 妥 当 性 向 上 の た め ， DSMC と
FDTD-PIC-MCC，Full-PIC-MCC それぞれの結合，つまり連成解析に中性粒子の分布を取り入れた連成解析ツールの開発を行う．また，それにより解析結果にどのような変化が見られるか考察を行う．

8

第2章
2.1.

数値解析手法

解析全体のフロー

この章では本研究で用いている数値解析手法について述べていく．各計算モ
デルについて詳しく述べる前に，まずは解析ツールの全体の流れを述べる．本研究室にて開発された数値解析コードは 3 つの物理現象に基づき, 3 種の計算モデルを用いている．表 2-1 に示す．
表 2-1
物理対象
モデル名
DSMC
中性粒子, スパッタ粒子の輸送および壁面損耗解析
電子とマイクロ波の相互作用による電子のエネルギー獲 FDTD-PIC-MCC
得解析
イオンと電子の挙動. 静電場解析
Full−PIC−MCC
これら 3 つの計算モデルを弱連成的に解析させることで，プラズマ挙動の解明
および寿命の評価ができる．図 2-1 に解析全体のフローチャートを示す．まず，前章で述べたようにマイクロ波放電式中和器の基本原理は永久磁石とマイクロ波の ECR 共鳴であるため，静磁場の情報が必要となる．中和器に永久磁石を配置した際の静磁場分布は，解析を始める前に市販ソフトで計算しておく(後述)．
さらに，中和器内部の中性粒子の密度分布を，DSMC を用いて計算する．これ
らの中性粒子密度データ，静磁場データを以降の FDTD-PIC-MCC および Full
− PIC-MCC の 入 力 デ ー タ と し て 使 用 す る ． こ の 中 性 粒 子 密 度 デ ー タ を
FDTD-PIC-MCC および Full−PIC-MCC に反映させたことが，本研究の要点である．
FDTD-PIC-MCC では，一様状態から電子とマイクロ波の吸収が計算され，電子エネルギー獲得確率が算出される．これは，単位時間単位体積あたりの電子のエネルギー増分を確率密度分布として算出後， 電子のマイクロ波吸収量の効果として Full−PIC-MCC に渡す．Full−PIC-MCC では，イオンと電子を一様状態から計算し，静電場の効果による粒子の輸送を解析する．プラズマ密度が変化するとカットオフ密度(電磁場がほとんど浸透不可能である密度)以上の分布が出現してきて，密度分布が変化するため，電子のマイクロ波吸収量が変化する． 実際に，Full-PIC の計算において，初期電子エネルギー獲得確率でエ9

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

10

図 2-1 解析全体のフローチャート

2.2.

プラズマ粒子計算手法 (13)FDTD-PIC-MCC と Full-PIC-MCC は，プラズマ粒子計算のモデルは同等のものを利用している．但し，Full−PIC-MCC は 2 次元円筒座標，FDTD-PIC-MCC
は 3 次元直交座標を用いているため差分式が異なる．また，衝突計算に関して
は，DSMC，FDTD-PIC-MCC，Full-PIC-MCC はすべて共通である．よって，共通項をこの節で述べ，共通項以外については，以降の節で各モデルそれぞれの詳細として述べる．

2.2.1.

基礎式

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

示す．
m

dv
 q E  v  B 
dt
dx
v
dt

(2.2.1)(2.2.2)

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

2.2.2.

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

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

 0 k bTeDe 

ne e 2(2.2.3)に対して，∆≤ 𝜆𝐷𝑒(2.2.4)であることが望ましい．また，タイムステップはプラズマ周波数 ω pe，

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

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

c 

eB
me

12

(2.2.8)に対して，

t  c1(2.2.9)である必要がある．
また，粒子はタイムステップ毎に１セル以上移動してはならない．特に荷電粒子は格子ごとに規定される電磁力が働くため，この条件を満たすことは重要となる．これを Courant の安定化条件という．3 次元の際，時間ステップ t は，
1
Dt £
1
1
1
(2.2.10)
v
+
+
2
2
2
(Dx ) (Dy ) (Dz )を満たさなければならない．この条件は極めて厳しく，わずかでも満足しなけ
れば計算が不安定となり，解析結果が発散してしまう要因となる．実際には，
さらにこの 10 分の 1 程度に設定することが望ましいとされている．ここにおいて，v は，粒子の最大熱速度または，電磁波速度つまり光速である．FDTD-PICと Full-PIC において，条件が異なるので各モデルの計算条件において説明する．

2.2.3.

粒子計算 PIC 法 (15)粒子ごとに基礎式を記述し解くものを粒子法という．本解析では 1012 個程の粒子計算を行うため，最新のスーパーコンピュータを用いても計算コストの面から実用的でない．そこで，粒子シミュレーションで取り扱う粒子として，実際の粒子の電荷質量比を一定に保ったまま多数の粒子の電荷と質量を一つにまとめた超粒子を使用する．１セルに含まれる超粒子数が 10 個以上になる事と取り扱う粒子数よりも超粒子の重みが尐ない事（統計性から）から決定する．
また，本数値解析では，電場量を空間格子点上で求める.そのため,各粒子の電荷を格子点上に振り分ける必要がある.また,粒子の運動方程式を解く際に,電場解析によって格子点状に割り振られた値を各粒子に内挿する必要がある.
このような計算を行うことを PIC 法（Particle In Cell method）と呼び，粒子シミュレーションにおいての標準的な手法である．

2.2.4.

衝突計算 MCC 法

粒子コードにて，衝突を取り扱う場合，粒子の衝突の解析には，
MCC(Monte-Carlo-Collision)法を用いるのが，計算速度の面から考えて最適な手法である．MCC 法では，粒子同士の衝突判定を各粒子間の衝突周波数から決定する確率と乱数を用いて判定する．詳細を以下に述べる．
13

2.2.4.1.

衝突種別

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

2.2.4.2.

平均自由行程と透過度 (17)面  ，厚さ l ，その物質の数密度 N の板を考える．エネルギー E の粒子ビームが，面  に垂直に当るものとする．衝突をせずに板を透過してくる粒子数
は，入射粒子数より減尐する.その減尐度は，  を衝突の断面積とすれば，dn  nNdl(2.2.11)と表わされ，これよりn  n0 exp  Nl (2.2.12)となる．ここで，n0 は面 （ l  0 ）でのビーム中の粒子である．それゆえ，（ l ,l  dl ）間で最初の衝突が起こる確率は，
Pinitial l dl  exp  Nl Ndl(2.2.13)となり， l 以下で最初に衝突が起こる確率は次式のようになる．
Pl    Pinitial l dl
l0 1  exp  Nl (2.2.14)

最初の衝突までに粒子が移動する平均距離  は Pl の一次のモーメントである．
   lP l dl
0(2.2.15)1N
これを平均自由行程という．よってモンテカルロ法によって，任意の点から出

発し最初の衝突が起こるまでの距離 l を決めるには，乱数 r を用いて，r  Pl  l
 1  exp    (2.2.16)で表わされる．ここで，乱数は(0,1)間に一様に分布する乱数であり，生成周期の短い fotran の組み込み乱数 ran()の使用は避け Mersenne Twister (19)を使用した．周期は(219937-1)である．式(2.2.16)を変形して，l  ln 1  r 
14(2.2.17)

となる．ここで， r は一様乱数であるから， 1  r を r とすれば，l  ln(r )(2.2.18)と表わされる．また，速さ v の粒子が気体分子と衝突する際の断面積を  とすれば，この粒子が１秒間に気体分子と衝突する回数  は，
v


 Nv

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

2.2.4.3.

衝突確率・衝突判定

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

 total  ngas totalvm

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

(2.2.21)(2.2.22)

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

2.2.4.4.

衝突種別判定乱数 r2 を用いて以下のように衝突の種類を決める．ここでは例として，電子と中性粒子の衝突について示す．弾性衝突断面積を s elastic ，励起衝突断面積を ex citatio n，電離衝突断面積を  io n izatio とする．
n
 0  r2  elas tic
 to tal elas tic
 r2  elas tic ex citatio n to tal to tal elas tic  ex citatio n




 r2  elas tic ex citatio n io n izatio n to tal to tal



であれば弾性衝突であれば励起衝突
であれば電離衝突

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

15

る．

2.2.4.5.

エネルギーと散乱角

粒子は衝突の前後でエネルギーの受け渡しを行い散乱される．以下にそれぞれの衝突におけるエネルギーと散乱角について述べる．
①弾性衝突
まず，分子同士の弾性衝突を考える．速度 v1 ,v 2  を持つ分子対の衝突後の速度を v1 ,v 2  とすると，
1
v1  v 2  u2
1
v 2  v1  v 2  u2v1 (2.2.23)(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 とすれば，次式が成立する (18)．
m  M cos M sin v e 
ve 
h(2.2.25)mMmM
h の直交座標成分はhx  ver cos hy  

hz = -(2.2.26)

vex vey cos   ve vez sin vervex vez cosj + ve vey sin jver(2.2.27)(2.2.28)

ここで，2
ver  vey
 vez2ve  v  v  v
2
ex

2
ey

(2.2.29)2ez(2.2.30)であり，  は散乱角，  は任意の角である．式(2.2.25)はエネルギー損失を含んでいる．
次に，衝突前後のエネルギー収支を考える ．衝突前の電子のエネルギーを
Eincident,e とし，衝突後の電子のエネルギーを Escattered,e とすれば，エネルギー収支

は以下のように表わされる．

Escattered,e  Eincident,e(2.2.31)電子は，エネルギー Eincident,e が小さいとき等方散乱され，大きいときは主に前方16散乱となる．これを考慮した散乱角  の確率密度を g   として以下に示す (33) ．
ここで， Ein cid enet, の単位を eV とする．
Eincident,e sin 
g   
21  Eincident,e sin 2 ln 1  Eincident,e 2
これを解けば，cos   1 

0     

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

(2.2.32)(2.2.33)となる．これを式(2.2.25)に用いて衝突後の速度を求める．このとき任意の角は  2 r4(2.2.34)

として与える．

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

式(2.2.33)を用いるが式(2.2.33)中の Ein cid enet, には式(2.2.35)から求められる
Escatterede, を用いる．励起エネルギー Eex citatio nは，キセノンの場合 8.34 eV を用いて
いる．
③電離衝突
入射電子のエネルギーを Eincident,e ，散乱電子のエネルギーを Escattered,e ，生成電子のエネルギーを Ecreatede, ，電離エネルギーを Eio n izatio とする．エネルギー収支はnEscattered,e  Ecreated,e  Eincident,e  Eionization(2.2.37)と表わされる．
電離後の散乱電子のエネルギーは以下のように求められる．

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

(2.2.38)

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

E  Eincident,e  Escattered,e(2.2.39)と表わされる．式(2.2.36)から電離後の電子の速度は，v~  v e 1 EEincident,e(2.2.40)となり，励起衝突の場合と同様に，この速度で弾性衝突をすると考える ．式(2.2.27)中 v e に v~ を用いて，衝突後の速度 v e を求める．散乱角には式(2.2.33)を用いるが，式(2.2.33)中の Ein cid enet, には式(2.2.38)から求められる Escattered,e を用いる．
電離エネルギー Eio n izatio は，キセノンの場合
12.13 eV である．
n
電離後の生成電子の持つエネルギーは，Ecreated,e  E  Eionization(2.2.41)となる．また，入射電子の速度 v e を用いて生成電子の速度はv created,e veve

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

2.2.4.6.

Null-collision 法 (21)(22)(23)

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

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

18

Null-Collision Process が選ばれたときには，粒子の運動は何ら変化しないものとする．図 2-2 に電子と中性粒子の衝突におけるキセノンガスの中性粒子密度 で 規 格 化 し た 衝 突 周 波 数 を 示 す ． 本 計 算 で は Null-Collision 法 で 用 い た co llisio n ng asの値は 6.5×10-13 m3/sec である．

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

2.2.4.7.

タイムステップの制約

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

P i2
1- P i(2.2.46)であるので，衝突ミスの回数を１％いかにするには𝑃𝑖 ≦ 0.1とする必要がある．
nullcollision を用いる場合は，𝑃𝑛𝑢𝑙𝑙 ≦ 0.1 すなわち
0.1
Dt £

n null

(2.2.47)

となる．

2.3.

DSMC

DSMC は DSMC 法 (24)(25)(Direct Simulation Monte Carlo method)を用いて，
19

中性粒子, スパッタ粒子の輸送および壁面損耗解析を行う．本研究ではスパッタ粒子解析は行っていないため，ここでは中性粒子解析についてのみ記述する．

2.3.1.

中性粒子解析手法

推進剤の中性粒子に対する粒子解析に用いた DSMC 法は，多くの量の粒子の位置座標と速度成分を記憶し，粒子の移動と粒子間衝突の計算を繰り返すこと
によって流れ場の状態を算出する手法である．この方法は多くの乱数を用いる
ため，ここでも Mersenne Twister (19)を使用している．

2.3.2.

中性粒子の運動中性粒子は電荷を持たないため，電場・磁場の影響を受けない．よってそれらを考慮する必要はない．計算では 1 ステップ経過するごとに全ての粒子を時間ステップ𝛥𝑡だけ動かす．動かす前の位置を𝒓𝟎 ，速度を𝒄𝟎 とすると𝛥𝑡後には𝒓 = 𝒓𝟎 + 𝒄𝟎 𝛥𝑡(2.3.1)になる．
また，全ての粒子の運動を追跡するのは計算コストの上で非効率であるため，ひとつの粒子が複数の粒子を表す超粒子として扱う．

2.3.3.

中性粒子の初期条件と境界条件

中性粒子の粒子数は任意の推進剤流量を設定し，それに対応した粒子数をス
テップごとに流入するようにし，徐々に中和器内部が粒子で満たされる．流入する粒子の初期速度 (gap)はマクスウェル分布より抽出した．気体定数を R とす
ると初期速度は x,y,z 方向に対しそれぞれ
𝑐𝑥 = 𝐴𝑐𝑜𝑠𝜃

(2.3.2)𝑐𝑦 = 𝐴𝑠𝑖𝑛𝜃

(2.3.3)𝑐𝑧 = 𝐴′ 𝑠𝑖𝑛𝜙′(2.3.4)となる．ここで𝜃 = 2π𝑈，𝜙 ′ = 2𝜋𝑈で𝑈は0 < 𝑈 < 1の一様乱数で𝐴 = (−2𝑅𝑇0 𝑙𝑛𝑈)1/2
𝐴′ = (−2𝑅𝑇0 𝑙𝑛𝑈)1/2
20(2.3.5)(2.3.6)である．
また，中性粒子解析における境界条件は中和器の全ての壁面で乱反射が行われるとし，オリフィスから中和器外部に出たものについては計算領域外となる
ので消失させる．

2.3.4.

中性粒子間の衝突

中性粒子解析においては中性粒子間の衝突については最大衝突数法を用いている．𝛥𝑡を平均自由時間以下にとると，分子の運動と衝突を分離してよいという，分離の原理が成り立つ．このため，全ての粒子を無衝突で𝛥𝑡走られた後，衝突して粒子は速度を衝突後のものに書き換えればよい．各セル内で起こる衝突は他のセルでの衝突に無関係に扱う．したがってどの分子も，衝突相手は自分の属するセルの粒子から選ばれる．ここで粒子は直径 d の剛体球であるとする．
あるセルに N 個の超粒子が含まれているとする．これらの粒子の平均速度を𝐶̅ と すると，１つの粒子が単位時間に衝突する回数ν(s−1 )は𝐶̅
(2.3.7)
ν = = 𝑛𝑔̅ 𝜍𝑇𝜆となる．ここで，σT は衝突断面積であり，λ = 1/√2𝑛𝜍𝑇 を用いた．また，２粒子間の平均相対速度は𝑔̅ = √2𝐶̅ で与えられることを用いた．セル内の超粒子 i に着目すると式(2.3.2)の𝑔̅ で 近似できる．
𝑁1𝑔̅ = ∑ 𝑔𝑖𝑗
𝑁(2.3.8)𝑗=1

ただし，𝑔𝑖𝑗 = |𝒄𝒋 − 𝒄𝒊 |で，𝑔𝑖𝑖 = 0である．しかしボルツマン方程式から，上式の1/Nを1/(N − 1)で置き換えてはいけないことが分かっている．𝑔̅ を 式(2.4.2)に代入し，νΔtがνΔt < 1のとき粒子 i の衝突確率𝑃𝑖 を表すことに注意すると𝑁

𝑃𝑖 = 𝑁−1𝑁

𝑛Δ𝑡 ∑ 𝑔𝑖𝑗 𝜍𝑇 = ∑ 𝑃𝑖𝑗
𝑗=1(2.3.9)𝑗=1

ただし𝑃𝑖 = 𝑁 −1 𝑛Δ𝑡𝑔𝑖𝑗 𝜍𝑇 である．𝑃𝑖𝑗 は粒子 i がΔ𝑡後に粒子 j と衝突する確率であるので相対速度𝑔𝑖𝑗 が大きいほど衝突を起こしやすいことが分かる．今着目しているセル内でΔ𝑡時間内に起こる衝突の回数𝑁𝑐𝑜𝑙 は
𝑁

𝑁

𝑁

1
1
𝑁𝑐𝑜𝑙 = ∑ 𝑃𝑖 =𝑛Δ𝑡 ∑ ∑ 𝑔𝑖𝑗
2
2𝑁
𝑖=1

(2.3.10)𝑖=1 𝑗=1となる．全ての𝑃𝑖 が１のとき，衝突する粒子の数は N であり，衝突回数は N/2
21である．このため，因子 1/2 が必要である．式(2.3.5)の二乗和を求めるためには(𝑁 2 − 𝑁)/2個の𝑔𝑖𝑗 を計算する必要があるが，計算量が N の二乗に比例するため計算効率を著しく低下させる．そこで𝑔𝑖𝑗 の最大値𝑔𝑚𝑎𝑥 を適当に見積もり，式(2.3.5)の𝑔𝑖𝑗 を全て𝑔𝑚𝑎𝑥 で置き換えた衝突数𝑁𝑚𝑎𝑥 を用いる．
1(2.3.11)
𝑁𝑚𝑎𝑥 = (𝑁 − 1)𝑛𝑔𝑚𝑎𝑥 𝜍𝑇 Δ𝑡2
また，数密度 n は次の式で表される．
𝑊𝑁
𝑛=(2.3.12)Δ𝑥Δ𝑦Δ𝑧
ここで，W は重みである．λは𝜍𝑇 に反比例するので𝑁𝑚𝑎𝑥 は𝜍𝑇 に依存しなくなる．
𝑁𝑚𝑎𝑥 は𝑁𝑐𝑜𝑙 を過大に見積もったものであり補正する必要がある．まず，着目しているセルに対して𝑁𝑚𝑎𝑥 回の衝突を計算する．その回ごとに二つの粒子を選
びその粒子対(i,j)に対する𝑔𝑖𝑗 /𝑔𝑚𝑎𝑥 の確率で本当の衝突を起こす．
対(i,j)が衝突するとき衝突方向は乱数により選ばれるようにした．衝突後の速度(𝑐𝑖 ′, 𝑐𝑗 ′)は
1(2.3.13)𝑐𝑖 ′ = 𝒘𝑖𝑗 − 𝑔𝑖𝑗 𝑹
2
1(2.3.14)𝑐𝑗 ′ = 𝒘𝑖𝑗 + 𝑔𝑖𝑗 𝑹
2で与えられる．ただし，𝒘𝑖𝑗 は重心の速度，R はランダムな方向を持つ単位ベクトルである．

2.4.

FDTD-PIC-MCC

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

2.4.1.

プラズマ粒子計算

2.4.1.1.

運動方程式解法

基礎式は，2.2.1 で述べたものであり，式(2.2.1)に leap-frog 法による時間中心差分を施すと，
x n 1  x n
 v n 1 2
t

(2.4.1)となる．ここで，各変数の上付き添え字は時間に関するもので， n はタイムス22

テップである．式(2.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.4.2)である． n-1 2 はある時刻 n から t 2 だけ以前の値であることを示す．同様に，n  1 2 はある時刻 n から t 2 だけ進んだ値であることを示す．
プラズマ粒子の運動を求めるために， 式(2.4.2)において未知量である v n1 2 について解き，求めた v n1 2 と式(2.4.1)より新たな粒子の位置 x n  1 を求めることができる．式(2.4.2)を変形すると，

Av n1 2  S(2.4.3)と書ける．ここで， a ・  ・  ・  を次のように定義する．
qt n

Bx
2m
qt n

By
2m
qt n
 
Bz
2m
  1  2   2   2(2.4.4)(2.4.5)(2.4.6)(2.4.7)

なお，ここでの 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.4.8)(2.4.9)

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

2.4.1.2.

(2.4.10)(2.4.11)

PIC 法重みつけ計算

3 次元における重みつけ方法を示す．まず，粒子が位置するセル座標 (i , j , k )を求める．次に，図 2-3 に示す dx1, dx 2, dy1, dy 2, dz1, dz 2 を求める．その値を用いて，図 2-3 にあるような直方体 V1〜V8 の体積を求める．この体積比から粒子が感じる電磁場を内挿することができる．
23

具体的には，座標 (x , y , z )上にある粒子が感じる電場 E (x , y , z )は，計算セルの体積を V ，図に示すような粒子によって仕切られた空間体積をそれぞれ V1～V8とすると，格子点上の重みは

An =

Vn
V

(n = 1 〜 8)

(2.4.12)となり，
E ( x, y, z ) = A1E (i, j, k ) + A2 E (i +1, j, k ) + A3 E (i +1, j, k +1) + A4 E (i, j, k +1)+ A5 E (i, j +1, k ) + A6 E (i +1, j +1, k ) + A7 E (i +1, j +1, k +1)(2.4.13)+ A8 E (i, j +1, k +1)と表わされる．

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

2.4.1.3.

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

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

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

r (i , j , k ) =24

また，電流密度はj  qnv(2.4.15)より，

j (i , j , k ) =2.4.1.4.

V1
qNv
V2

(2.4.16)

粒子の壁面境界条件

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

2.4.1.5.

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

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

2.4.2.

電磁波解析手法 (26)(27)

2.4.2.1.

基礎式

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

H(2.4.17)t
E  H  0
J
(2.4.18)t
ここで， E は電場， H は磁場，  0 は真空の透磁率，  0 は真空の誘電率， J は電
  E   0

流密度を表わす．

2.4.2.2.

差分式

三次元の FDTD 法で用いられる電場，磁場の空間配置を図 2-4 に示す．このように，電場・磁場ともに空間的に交互に配置され，電場はセルの各辺に沿って，磁場はセル面の中心で与えられる．これは Yee によって提案され (28)，電場の回転（   E ）が磁場を，磁場の回転（   H ）が電場の方程式を満たすように配置されている．
この電磁場の配置のもと，式(2.4.17), 式(2.4.18)を空間，時間とも中心差分した式を以下に示す．式((2.4.19)~(2.4.24))ここで， 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
ïî
ïþ
t
E xn 1 i  1 2 , j , k   E xn i  1 2 , j , k   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z26

(2.4.19)(2.4.20)(2.4.21)(2.4.22)
 J x i  1 2 , j , k 

E yn 1 i , j  1 2 , k   E yn i , j  1 2 , 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.4.23)

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



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.4.24)

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


1

n

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

27

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

2.4.2.3.

アンテナ解析

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

V は中心導体と外部導体間の電圧, 𝑟 は動径方向の距離, 𝑎 はアンテナの半径である.
また, 同軸ケーブルのインピーダンス Z は以下の式で表される.
28

𝜇 1 𝑅𝑍=√ln
𝜀 2𝜋 𝑟

(2.4.27)

ここで, R は外部導体の半径, r は中心導体の半径である. マイクロは電力を P
とすると，
𝑃=

𝑉2
2𝑍

(2.4.28)となり，式(2.4.27)，式(2.4.28)より， V を求め，式(2.4.25)，式(2.4.26)よりア
ン テ ナ に 与 え る 電 場 を 求 め る ． こ こ で ， Ex を 𝑟 = 𝛥𝑥/2 の 電 場 で ， Ey を
𝑟 = 𝛥𝑦/2 の電場で代表させると，
Ex (IFED, JFED, KFED) = -E x (IFED -1, JFED, KFED)

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

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

(2.4.30)

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

図 2-5 同軸線路給電

29図 2-6 アンテナモデル

2.4.2.4.

PEC(Perfect Electric Conductor )境界条件

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

2.4.2.5.

Mur の吸収境界条件 (29)

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

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

(2.4.31)と表わされ，次の微分方程式を満足する．
E z 1 E z
0(2.4.32)
x v tx = 0 の境界で完全に反射がないならば，電場は式(2.4.31)の形を保ったまま伝播するはずである．すなわち，境界でも式(2.4.32)が満足される．これを FDTD
30

法の形式に定式化するため，時間について差分化すると¶E z E zn - E zn -1
=
¶t
Dt

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

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.4.34)となる. 式(2.4.33)の Ezn 1 2, Ezn1 1 2 も FDTD 法では直接求めることはできな
いため，前後の値で平均して求める．そして， E z についてまとめるとEzn ( 0, j, k +1 2) = Ezn-1 (1, j, k +1 2)(2.4.35)vDt - Dx n
Ez (1, j, k +1 2) - Ezn-1 ( 0, j, k +1 2)}{
vDt + Dxとなる．これを Mur の 1 次吸収境界条件という．上式は吸収境界に平行な電場+成分について適用される．

2.4.2.6.

ポインティングベクトル (30)

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

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

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

2.4.2.7.

投入電場の補正 (37)

投入マイクロ波電力の制御因子は, 先に示したアンテナ根本に励振する電場強度である.時間発展するにつれ,電子 の状態は変化するため,アンテナから見た負荷インピーダンスは変化し,一定の電場を根本に与えた場合,投入 されるマ31

イクロ波電力も時間とともに変化する.本計算 では,一定のマイクロ波電力を放電室内に投入するために,以下の式にしたがって,根本に与える電場強度をマイ
クロ波 1 周期毎に更新した.
𝐸𝑛𝑒𝑤 = 𝐸𝑜𝑙𝑑 +𝑃𝑡𝑎𝑟𝑔𝑒𝑡𝜃( 𝑃− 1) 𝐸𝑜𝑙𝑑𝑜𝑙𝑑2(2.4.37)

ここで，𝐸𝑜𝑙𝑑 , 𝐸𝑛𝑒𝑤 は，更新前後の励振電場強度，𝑃𝑜𝑙𝑑 は更新前のマイクロ波 1周期で平均した投入マイクロ波電力，𝑃𝑡𝑎𝑟𝑔𝑒𝑡 は設定投入マイクロ波電力とする．

2.4.3.

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

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

32

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

2.5.

Full-PIC-MCC

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

2.5.1.

運動方程式解法

基礎式は，2.2.1 で述べた．式(2.2.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.5.1)

これらを四次精度 Runge-Kutta 法にて解いた．
tbr1  2br 2  2br 3  br 4 6tb1  2b 2  2b 3  b 4 vn 1  vn 
6tvzn 1  vzn  bz1  2bz 2  2bz 3  bz 4 6vrn 1  vrn (2.5.2)

但しここで定数 b はそれぞれ以下を用いる．
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 

br 4 

(2.5.3)qEr  v  tb 3 Bz  1 v  tb 3 v  tb 3 mrqE  vz Br  vr Bz   vr v
m
r
q

tb
tbr1   1 
tbz1 
tbr1 z1 b 2   E   v z  Br   vr 
 Bz    v z  vr m2 2   r
2 
2 



b1 

b 3 

q
tbz 2 tbr 2   1 
tbz 2 tbr 2  Br   vr 
 Bz    v z  vr  E   v z 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 m
r34(2.5.4)q
Ez  v Br mtb 1  q Br 
bz 2   E z   v m
2  
bz1 

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

(2.5.5)q
Ez  v  tb 3 Br m

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

𝑞
𝐸,
𝑚 𝑟

𝑣𝑧 =

𝑞
𝐸
𝑚 𝑧

(2.5.6)とした．

2.5.2.

PIC 法重みつけ計算

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

Sn
S

n  1 ～ 4となり，
Er, z   A1Ei, j   A2 Ei  1, j   A3 Ei  1, j  1  A4 Ei, j  1と表わされる．

35(2.5.7)(2.5.8)図 2-8 PIC 法におけるセルの構成

2.5.2.1.

電荷密度計算

電荷密度は図 2-8 に示されるセルの格子点上で定義され，セル中に N 個の粒子を代表した 1 個の超粒子が存在したとき，格子点 i, j  上に割り振られる電荷
密度は，
 i, j  

S1 qN
S V

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

2.5.2.2.

電子質量の変更 (31)

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

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

衝突断面積を修正する必要がある．表 2-2 に現実値と仮想値を示す．
表 2-2 質量比の変更による補正係数

2.5.2.3.

現実値

仮想値

𝑚𝑒

𝛼 𝑚𝑒

𝑣𝑚

𝑣𝑚
√𝛼

∆t

√𝛼 ∆t

B

√𝛼 𝐵

𝜍

√𝛼 𝜍

粒子の境界条件

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

2.5.2.4.

Digital Filter (35)

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

side
cornerMΦi,j,k + SΦi,j,k+ KΦi,j,k(2.5.11)
M + 6S + 12K
new
Φi,j,k を格子点上に定義される値とし、Φi,j,k
は Digital Filter を通した新しいnewΦi,j,k←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 + Φ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+ Φi−1,j−1,k + Φi−1,j+1,k + Φi+1,j+1,k + Φi+1,j−1,k(2.5.12)(2.5.13)

ここで、M、S、K の値をそれぞれ 12、4、1 とすると式(2.5.11)は次式となる。
1
1 side1 cornernewΦi,j,k← Φi,j,k + Φi,j,k+ Φi,j,k(2.5.14)4
1248
電子密度およびイオン密度，電位を計算した後に,それぞれの値にこの Digital
Filter を作用させる.

2.5.3.

電場解析手法 (31)2.5.3.1.

基礎式

ポアソン方程式および電場と電位の関係式を基礎式として計算を行った．
rÑ 2j = (2.5.15)
e0

E = -Ñj

(2.5.16)式(2.5.15)を差分化し,連立方程式として解いた．差分式を以下に示す．
i , jaa1  i 1, j  aa 2  i 1, j  aa3  i , j 1  aa 4  i , j 1  aa5  i , j 
01
aa1  aa 2  2
z1
1
aa3  2 
(2.5.17)r
2ri , j r
aa 4 

11
2
r
2ri , j r

æ 1
1 ö
aa5 = -2 ç 2 + 2 ÷è Dz Dr ø

ここで， i，j はそれぞれ z，r 方向の格子点番号に対応している．𝛷𝑖,𝑗 は i，j における電位を表し，．𝜌𝑖,𝑗 はその点での電荷密度を表す．𝛥 𝑧，𝛥 𝑟は各方向の格子サイズであり，𝑟𝑖,𝑗 はその格子点の r 方向における位置を表す．
また，r = 0 上での係数は，
38

aa1  aa 2 1
z 2

aa3 = 04
aa4 = 2Dr

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

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.5.19)2Dとなり，静電場を求める．

2.5.3.2.

PETSc

式(2.5.15)はすべての格子点に１つずつ与えられるので，連立方程式として解くことができる．その方程式を数多く反復計算する際に用いる計算手法として，
PETSc というライブラリを使用している．PETSc とはアメリカの Argonne
National Laboratory 内 に 発 足 し た Mathematics and Computer Science
Division project の メ ン バ ー で あ る Satish Balay, William Gropp, Lois
Curfman McInnes, Barry Smith の手によって開発されているプログラムライ
ブラリであり現在アメリカでは航空力学や音響学等の数値解析に利用されている. 数多くの反復解法と前処理と呼ばれる反復計算をする前に行列式を適切な形に加工する方法の組み合わせを容易に試行できる．
利用方法としては，式((2.5.17),(2.5.18))で示した係数行列と電荷密度分布を
PETSc の用意するデータ構造体に格納し，実行オプションとして，反復法と前処理を指定するものである．今回，選択した反復解法として，Bi-CGStab 法，
前処理法として，Geometric-Algebraic Multigrid 法を選択した．これは，数種の反復解法と前処理の組み合わせを現モデルで試行し，安定解が得られる組み合わせの中で，最も収束性の高かったからである．尚，計算のアルゴリズム
については，今回記載しない．

2.5.3.3.

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

2.5.4.

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

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

40

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

2.6.

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

図 2-10 に本研究室所有のマイクロ波放電式中和器を yz 断面で切り取った断面図を示す．ここにおいて，z 方向を図の左から右方向，y 方向を図の下から上方向としており、rz 断面図の際も r 方向が y 軸の向きと一致している．形状としては，円筒径の放電室(Wall)を箱状の永久磁石(サマリウムコバルト磁石)が z
軸方向に 3 個を 1 組として，棒磁石として 15 組の磁石が周方向に均等に配置されている．磁石を鉄で出来た FrontYork および BackYork が挟んでおり，こ41

の間に磁力線を形成する形状となっている．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-11 に示す．これは，Advanced Science Laboratory, Inc 社の 3 次元磁場解析ソフト Amaze を用いて解析した磁場強度分布である．今回使用するマイクロ波周波数は 2.45GHz であるので，それに対応する ECR 磁場強度である 85mT を赤線で示している．
赤丸で囲んだところに，磁気チューブと呼ばれる電子補足層があり，この領域
を電子は，E × Bドリフトまたは grad B ドリフトをしながら移動する．

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

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

2.7.

弱連成

弱連成では，FDTD−PIC-MCC から Full−PIC-MCC にはエネルギー獲得確率を Full−PIC-MCC から FDTD−PIC-MCC には電子密度及び電子エネルギーである．ここで，上述したように 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 のフィードバックはあまり意味をなさない．よって，あ43

る程度電子とイオンを動かす必要があり，この時間を経験的に 250ns 程度とし
た．

2.8.

各種物理条件と計算条件

2.8.1.

物理条件

表 2-3 に本計算に用いた物理条件を示す． 本研究室の実験体系と同等のものを扱っている．イオン初期エネルギーは室温として，初期密度初期電子エネ
ルギーに関しては，過去の実験データより設定した (38)．
表 2-3 物理条件
Xe
0.05 mg/s
2.45 GHz
8W
1.0×1017 m-3
2.0 eV
0.026eV(300K)
40 V

推進剤
推進剤流量
マイクロ波周波数
投入電力
初期密度
電子初期エネルギー
イオン初期エネルギー
出口電位 𝛷𝑝𝑙𝑢𝑚𝑒2.8.2.

計算条件

FDTD−PIC-MCC
メッシュサイズを 0.5 mm とし，x 方向に 40 格子 y 方向に 40 格子 z 方向に
36 格子とした．メッシュサイズ 0.5 mm より FDTD における計算制約条件は，−1
プラズマ周波数より， 0.2 𝜔𝑝𝑒= 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 となり，さらにその 1/10 程度とすると 7×10-12
sec となる．このなかで最も制約条件の厳しいのは，Courant 条件であるためこれを満たしさえすればよく．FDTD10 回につき PIC を 1 回解くとすると，2.5
×10-12 sec で十分計算の制約を満たす．

44

表 2-4 計算条件(FDTD−PIC-MCC)
0.5 mmメッシュサイズ
タイムステップ(FDTD)
2.5 ×10-13 sec
タイムステップ(PIC)
2.5 ×10-12 sec
超粒子の重み9.85 ×104
超粒子数
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 となる．一般的にはさらにその 1/10とすべきなので，4.0×10-12 sec，最小値は 4.0×10-13sec となる．しかし，タイムステップを尐しでも広くとるために，Courant の 1/8 程度として，5.0×10-13
sec とした．これを仮想値に変換すると，2.5×10-11 sec つまり 0.025ns とした．
表 2-5 計算条件(Fll−PIC-MCC)
0.025 mm
メッシュサイズ
タイムステップ(現実値)
5.0 ×10-13 sec
タイムステップ(仮想値)
2.5 ×10-11 sec
超粒子の重み3.90 ×104
超粒子数(イオン・電子)
1000 万個
1 セルあたりの粒子数
50 個

45

第3章
3.1.

結果と考察

DSMC

連成解析に用いる入力データとして，DSMC により中性粒子密度分布の計算をした．推進剤であるキセノンガスの流量は，0.05 mg/sとしている．このときの中性粒子の密度分布を図 3-1 に示す．中和器内部全体で密度の平均を取った
値は2.4 × 1020 𝑚−3 となった．ここで，アンテナ付近には他の領域に比べて高い
値の密度の分布が存在していることが分かる．このことが，連成解析にどのような影響を与えたかを，次節以降で述べる．

図 3-1 中性粒子密度分布

46

3.2.

各物理値分布

この節では連成解析ツールにより解析することができた 3250ns 時での結果を，中性粒子密度データを導入する前の解析結果と比較，さらには文献(39)と比較することにより，得られた各物理値に対しての考察を行う．

3.2.1.

電子エネルギー分布

図 3-2 に連成解析により得られた電子エネルギー分布を示す．まず，図 3-2(b)から見てみると，磁気チューブ内で電子温度が高い領域がみられ，アンテナ前方ではエネルギーの授受があまり起こっていない．それに対して中性粒子密度
データを導入した場合では，磁気チューブ内のエネルギーが低く，アンテナ前方からオリフィス出口付近のエネルギーが高くなっている．これは，図 3-1 から，中和器内部ではアンテナ前方付近に中性粒子密度が高い分布が多いため，
密度データを導入した連成解析ではアンテナ前方での衝突による新たなプラズ
マができやすくなっているからだと考えられる．その結果，エネルギーを受け取る電子がアンテナ前方に多く存在するようになり，分布に変化が生じている．
図 3-2(a)の電子エネルギー分布を文献(39)と比べてみても，アンテナ前方と磁気チューブでエネルギーが高くなっているといる傾向は一致しているため，
この解析結果は妥当性があるといえる．しかし，オリフィス出口付近のエネル
ギーが文献(39)よりも高くなっている．これは，オリフィス出口にかけた出口
電圧により電子が引き寄せられすぎているからである．実際の実験では中和器
外部からコレクタに電圧をかけて電子を引き出しており，文献(39)でも中和器
外部に出口電圧をかけている．よって，より現実値に近い結果を得るためには，
本解析コードの計算領域を広げる必要がある．

47

(a)

中性粒子密度データ導入後(3250 ns)(b)

中性粒子密度データ導入前(3250ns)
図 3-2 電子エネルギー分布

3.2.2.

空間ポテンシャル分布

図 3-3 に空間ポテンシャル分布を示す．ここで，各部材とプラズマ領域の境48

界面でシースが確認できる．図 3-3(a)より，シース電位は，オリフィスプレー
トで 20～35V 程度，フロントヨークでは 10～15V，壁面では 8～15V，アンテナに関してはアンテナ前方で 30V，それ以外のアンテナ周辺は 10～20V となっている．バックヨークにおいて 5V 程度ではあるが z=0～1mm と厚い幅でシースができている．これは，バックヨーク付近は次節の密度分布から分かるように電子密度が低く，デバイ長が3.8 × 10−2 mm と他の領域と比べて長くなるためである．これらのシース電位を文献(39)と比較すると，オリフィスとアンテナの値は近いものとなっている．しかし，それ以外の部材付近のシース電位は文献(39)よりも小さい値となっていることが分かる．これは，連成解析開始時にイオンは中和器内部に均一に配置するようにしているため，3250ns 時はまだアンテナより後方の領域にあるイオンが，定常状態より多く残っているからだと考えられる．定常状態とされる20μsまでの解析を行えるようにする必要がある．
また，密度データを導入した連成解析の出口電圧を，40V から変更して 30Vと 0V で行った場合の空間ポテンシャル分布を図 3-4 に示す．計算領域を広げた場合，z=18mm の領域には 30～40V のポテンシャルが分布すると予測されるため，図 3-3(a)と図 3-4(a)の間のような分布になると思われる．

(a)

中性粒子密度データ導入後(3250 ns)

49

(b)

中性粒子密度データ導入前(3250 ns)
図 3-3 空間ポテンシャル分布(a)

出口電圧 30V50

(b)出口電圧 0V

図 3-4 空間ポテンシャル分布

3.2.3.

空間密度分布図 3-5 に電子の空間密度分布を示す．ここで，プラズマとは準中性であり，電子密度とイオン密度においてはほとんど定性的，定量的差異はない．本解析の結果でも，電子とイオンの密度分布はほぼ同様なものとなった．そのため，
ここでは電子の密度分布のみを掲載している．図 3-5(a)より，アンテナ前方からオリフィス出口にかけての電子密度が高い値となっており，エネルギー獲得のモード変化の影響を受けていることが確認できる．しかし，文献(39)と比較すると，密度分布が中心軸により過ぎている．これは，本解析ではエネルギー獲得機構がアンテナ前方とオリフィス出口付近に偏っているからであり，これは後述する FDTD-PIC-MCC の計算不安定性によるものだと思われる．

51

(a)

中性粒子密度データ導入後(3250 ns)(b)

中性粒子密度データ導入前(3250 ns)
図 3-5 電子密度分布

52

3.3.

FDTD-PIC-MCC の計算不安定性

現在の解析コードでは，連成解析を進めて行くごとに FDTD-PIC の計算が不安定になり，連成解析の計算が発散してしまう．図 3-6 に FDTD-PIC-MCC より得られた電子エネルギー分布を示す．左が x=10mm での yz 断面，右が z=8mm
での xy 断面図である．図 3-6 より，アンテナ周辺でエネルギーが欠損するバ
グが起こっていることが確認できる．アンテナ周辺は電子のエネルギー獲得に非常に大きく関わってくる部分である．改良することにより，定常状態までの解析はもちろん，アンテナ周辺や磁気チューブ内でのより高いエネルギー授受の分布を得られることが見込まれる．よって，この部分のコード改良が急務だといえる．

図 3-6 電子エネルギー分布

3.4.

電子・イオン電流

3.4.1.

各部材の電子・イオン電流

図 3-7 に中性粒子密度データを導入した連成解析の各部材に対する電子・イオ
ン電流値を示す．また，図 3-8 に中性粒子密度データを導入する前の連成解析
での電子・イオン電流を示す．赤が電子電流，青がイオン電流を表す．
53

まず，中性粒子密度データを導入した連成解析の出口電位 40V と 30V の場合について考察する．イオン電流について考えると，図 3-7(a)からオリフィスプ
レートが最も値が高くなっていることが分かる．これは図 3-3(a)のシース電位からしても明らかである．それ以外では，壁面，アンテナの順で高い値となっているが，アンテナは表面積が狭く壁面は表面積が広いため，単位面積当たりのイオン電流量ではアンテナの方が壁面より高いといえる．図 3-7(b)と比較す
るとすべての部材において図 3-7(a)の方がイオン電流値は高い．これは出口電
圧 40V の方が，ポテンシャル勾配が広く分布しており，どの部材のシース電圧も高くなっているからである．これは電子電流についても同じことが言え，出口電位 40V の方がすべての部材において電子電流値が低くなっている．
次に図 3-7(c)を見ると，出口電位によるポテンシャル勾配が存在しないため，
イオン電流と電子電流の値が，どの部材でもほぼ等しい値となっていることが分かる．しかし，マイクロ波の影響により磁気チューブ付近の部材やアンテナ
ではわずかにシースが形成されるため，その部分の部材では若干であるがイオン電流が支配的となっている．
最後に中性粒子データ導入前について見てみると，イオン電流・電子電流ともにバックヨークが最も高い値となっている．これは図 3-5(b)から分かるように，プラズマがアンテナ後方に滞留してしまっているからである．

25
20

Current, mA

15
10
Ion
Electron

5
0
-5
-10
Antenna

FrontYork(a)

BackYork

Wall

OrificePlate

出口電圧 40V (3250ns)
54

10
8
6

Current, mA

4
2
Ion

0

Electron

-2
-4
-6
-8
-10
Antenna

FrontYork(b)

BackYork

Wall

OrificePlate

出口電圧 30V (3250ns)

8
6

Current, mA

4
2
Ion

0

Electron

-2
-4
-6
-8
Antenna

FrontYork(c)

BackYork

Wall

OrificePlate

出口電位 0V (3250ns)

図 3-7 壁面電流・および引き出し電流（中性粒子密度データ導入後）

55

30

20

Current, mA

10
Ion

0

Electron
-10

-20

-30
Antenna

FrontYorkBackYork

Wall

OrificePlate(a) 出口電位 40V (3250ns)
図 3-8 壁面電流・および引き出し電流（中性粒子密度データ導入前）

3.4.2.

引き出し電子電流表 3-1 にそれぞれの条件で連成解析を行った際の，オリフィス出口からプルーム領域へと排出される引き出し電子電流値を示す．中性粒子データ導入後の連成解析において，出口電圧が 0V の場合は引き出し電子電流値が 1.87mA であるのに対し，出口電圧 30V と 40V ではその 10 倍以上の値が得られている．
これは，図 3-3(a)と図 3-4(a)で示した出口電圧 30V と 40V のポテンシャル分布
が文献(39)と同様に，オリフィス出口から中和器内部に向けて滑らかなポテン
シャル勾配を形成しているからである．また，30V と 40V を比較すると約 3 倍の差が確認できる．文献 (39)ではオリフィス出口からアンテナ前方 にかけて
25V 以上でのポテンシャル勾配が広く分布しているため，オリフィス出口への
25V 以上のポテンシャル勾配の領域が広いほど引き出し電流値は高いということが伺える．それを踏まえ図 3-3(a)と図 3-4(a)をみると，出口電位 40V の方が25V 以上のポテンシャルとなっている領域が広いことは明らかであり，これにより引き出し電流に差が生じたことが分かる．
中性粒子データ導入の前後を比べると，中性粒子データを導入することにより引き出し電流値が 2 オーダー上がっていることが確認できる．これは上述した 25V 以上のポテンシャルの領域の広さはもちろん，図 3-5 の電子密度の分布も影響しているといえる． しかし，実験によりえられたデータから，引き出し56電流値は 120mA と見積もられており，今回解析により得られた値をもってしても実験値の 1/2 以下である．よって，前節にて述べた FDTD-PIC-MCC の不具合をはじめとするエネルギー獲得に関しての計算モデルの改善を行い，定常状態までの解析を行う必要がある．
表 3-1

引き出し電子電流値

出口電位 (V)
中性粒子データ導入後

中性粒子データ導入前

引き出し電子電流値 (mA)

40

51.4

30

18.6

0

1.87

40

0.42

57

第4章

結論

本研究ではマイクロ波放電式中和器の性能を統合的に解析できるツールの開発を行った．
DSMC による中性粒子密度を，FDTD-PIC-MCC と Full-PIC-MCC へ導入した連成解析を可能にした．これにより中和器内部のポテンシャル勾配がオリフ
ィス出口からアンテナ後方へと連続的に分布されるようになり，電子がプルー
ム領域へと引き出されやすい状態へ変化していく過程を捉えることができた．
プルーム領域への引き出し電子電流の解析をすると 51.4mA となり，中性粒
子密度データを導入する前と比べ，2 オーダーほど出力を向上させた．しかし，
これは実験値の 1/2 以下の値であるため，エネルギー獲得に関しての計算モデ
ルの改善する必要がある．また，解析結果の妥当性向上のために，計算領域の拡大を行う必要がある．

58

謝辞
まず，本研究に関して，終始ご指導ご鞭撻を頂きました山本直嗣准教授に心より感謝致します．数値解析の結果が出た際に共に考察をしていただけるなど，研究に対しての理解を深める手助けをしていただいた森田太智助教，研究室内で発表する際に誰よりも多くのコメントをくださいました中島秀紀名誉教授に深く感謝いたします．研究に関する事務手続きや，学生生活の様々なサポートをしていただいた馬渡隆子氏に深く感謝致します．
また，2 年間同じ苦楽を共にし，一緒に頑張って来た同期の皆様，研究室での楽しく充実した日々を与えてくださった M1 と B4 の皆様，D1，D2 および本研究室を卒業された諸先輩方に深く感謝致します．
最後に，ここまで学生生活を送らさせて頂き，私が進む進路を応援して下さった家族に深く感謝致します．

59

参考文献
(1)

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

細田聡史，國中均，”イオンエンジンによる小惑星探査機「はやぶさ」の帰還運用”，J.Plasma Fusion Res. Vol.86,No.5,(2010),pp.282-292(3)吉川真; はやぶさ 2 プロジェクト準備チーム”小惑星探査ミッション「はやぶさ 2」” 第 11 回宇宙科学シンポジウム: 講演集録 Vol 11 ,2011(4)中村 友哉, “超小型衛星ビジネスがもたらす宇宙イノベーショ
ン” ,JSASS-2014-4323(5)荒 木 友 太 , “ デ ブ リ 除 去 実 証 衛 星 「 ADRAS-1 」 の 開 発 状 況 ”,
JSASS-2015-4355.

(6)高山慎一郎, ”超低高度衛星技術試験機 (SLATS)の開発状況”, JSASS2014-4075(7)

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

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

(9)

大道渉，”マイクロ波放電式中和器のイオン電流測定”，STEP-2011-033(10)

N. Wallace, P. Jameson, IEPC, 2011, 327(11)

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

平成 25 年度九州大学修士論文．
(12)楠田将士，“マイクロ波放電式中和器の引き出し電流電圧特性解析”，

平成 26 年度九州大学修士論文．
(13)

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

(14)

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

ロ法”(海文堂,1998)(15)

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

Giovanni Lapenta, Jeremiah U. Brackbill, Journal of Computational

Physics, 115, (1994), 213.
(17)

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

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

Mersenne

Twister

Home

Page

:

(http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html)(20)

南部健一：“コンピュータアナリシスシリーズ 760

原子･分子モデルを用

いる数値シミュレーション(21)

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

V. Vahedi and M. Surendra, Computer Physics Communications, 87(1995), 179.
(22)

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

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

(24)

保原充，大宮司久明：”数値流体力学:基礎と応用”(東京大学出版会，1992)(25)

南部健一，渡辺安雄，米村茂，：”中間流・分子流のシミュレーション法-大面積微小ギャップの排気への応用-“，真空, 43(2000)(26)

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

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

K. S. Yee, IEEE Tarns. Antennas Propagation, 14 (1966), 302.

(29)

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

(30)

著=W. H Hayt，訳=山中惣之助，岡本考太郎，宇佐美興一：“工学系の基礎 電磁気学”(朝倉書店, 2001)(31)

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

田中守也，”電子・イオンビーム工学”(電気学会，1995)(33)

C. B. Opal, W. K. Peterson and E. C. Beaty, The Journal of Chemical

Physics, 55 (1971), 4100.
(34)

保原充，大宮司久明：”数値流体力学:基礎と応用”(東京大学出版会，1992)(35)

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

simulation." (1985).
(36)

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

Phys. B65B, (1989)(37)

渡邊裕樹, “マイクロ波放電式電子源内の電子加熱解析”.STEP-2013-014(38)

小野寺範義, 竹ヶ原春貴, 西山和孝, 船木一幸, 國中均: マイクロ波放電型中和器の電子放出機構, 日本航空宇 宙学会論文集, 第 49 巻, 第 564号, pp.27-31, 2001.
(39)

窪田

健 一 , “マ イ ク ロ 波 放 電 式 中 和 器 の 電 子 引 出 に 関 す る 3 次 元

Hybrid-PIC 解析”, 日本航空宇宙学会論文集, Vol. 63, No.5, pp. 197-203, 2015.

61