粒子を用いたプラズマ推進機の運動論的シミュレーション

概要

粒子法を用いたプラズマ推進機の運動論的シミュレーション
Full-kinetic particle simulation for plasma thruster
研究代表者：西山和孝（宇宙航空研究開発機構）
nishiyama@ep.isas.jaxa.jp
研究分担者：月崎竜童（宇宙航空研究開発機構）
tsukizaki.ryudo@jaxa.jp
担当：計算と比較用のレーザー計測実験
：張科寅（宇宙航空研究開発機構）
choh.shinatora@jaxa.jp
担当：計算コード開発
：山下裕介（東京大学大学院）
yamashita@ep.isas.jaxa.jp
担当：計算コード開発
：濃野歩（東京大学大学院）
nono-ayumu303@ecc.u-tokyo.ac.jp
担当：計算コード開発
研究目的 (Research Objective):
イオンエンジンをはじめとするプラズマ推進機では、推進剤を電離することでプラズマを発生
させ、イオンを高速後方排気することで推進力を得る。従来の化学推進と比べて燃費が良いこと
から、宇宙用推進機として利用される。宇宙機に搭載されたプラズマ推進機の寿命を律速する要
囚となるものの一つとして、電荷交換衝突イオン（ Charge exchange ion, CEX ion, 以下 CEX イ
オン）が挙げられる。 CEX 衝突は、高速イオンが推進機から漏れ出た推進剤ガス（中性粒子）の

間で起こる衝突であり、低速なイオンが発生する。この低速イオンは、宇宙機排気プラズマ間の
電位差によって宇宙機に向かって、逆流イオンとして輸送される。この逆流イオンは周辺の電子
機器の故障などの要因となることが懸念される。特に、当研究グループによる先行実験から、 CEX
イオンによって推進機周辺部の導電部位がスパッタリング損耗を受けると、推進機作動に影馨を
及ぼすことが明らかとなっている。以上のことから、 CEX イオンの物理過程を明らかにする必要
がある。
本研究では、 CEX イオンの物理過程を明らかにするべく、 CEX イオンの生成を考慮したシミ

Hybrid particle in cell
ュレーションモデルを構築した。宇宙機排気プラズマ間領域の電位構造を
(Hybrid-PIC) 法によって解いた。 Hybrid-PIC では、イオンは粒子として扱われ、電子は流体とし
て扱われる。

計算手法 (Computational Aspects):
本研究の手法に関する特徴は以下である。

1) 計算領域は 2 次元直交座標で行われた（ Fig. 1。
）

-65-

2) イオンビームは、実験的な測定値 2)をもとにした解析的な近似モデルで代替した。下式
は、今回新たに考案された混合ガウス分布によるビーム電流密度加の近似モデルであ
る
。
M

加＝午〗％
叫

l
bはビーム電流値、

m=1

1

［2下西

exp

（
長
2）

r
b=行＋ ztana、r
Tは推進機半径、 aはビーム発散角、

％ ,篇, w叫ますべてフィッティングパラメータである。なおここではM =3 とした。 2 次元空
間でも電流保存則が満たされるように設計した。
3)
3) 中性粒子分布は、解析的な近似にて表現した。下式は、推進機下流静電グリッドの各
孔から漏れ出す中性粒子の数密度分布を 3 次元空間でモデル化したものである。今回
のシミュレーションでは 2 次元であるが、簡単のためにこのモデルを用いた。 Ngriんま孔
数
、 0はグリッド面法線と孔中心からの位置ベクトルがなす角、 nnoはグリッド上流部での
中性粒子密度、 rは孔中心からの距離、 Aは孔の断面積である。式中の和記号は、各孔
からの漏れ出し効果を足し合わせることを意味する。
Ngrid

cos

四＝〗叫0冗r2%
4
4) 準中性電子流体近似下で、 CEX イオンを粒子として扱った。
5) 電子温度は、電子がポリトロープ流体であると仮定することでモデル化した。ポリトロピッ
クindex はシングルプローブで実験的に計測された電子温度分布に最も整合するように
定められた。
4)
6) 宇宙機表面は、シース効果を考慮した境界条件にてモデル化した。
7) 計算結果を別途実施された実験結果と比較検証した。
電子流体は、無磁場仮定下でドリフト拡散近似を用いて解かれ、 CEX イオンを含めたプラズ

マ密度が定常状態に遷移し、十分に時間が経てば、計算を打ち切った。時間刻み幅は

1.0 nsで、全体の計算継続時間は、物理時間で5.0 ms とした。定常状態に遷移するまでに、
1) 電位計算のみを行う段階、 2) 電位計算と粒子の運動の計算を行う段階、と 2 段階に計
算することで計算の発散を抑えた。
格子による空間分割数は 101 X 61であり、用いた超粒子は60,000個である。超粒子の挙動
を解く際には OpenMP を用いて並列化した一方で、電位計算時における楕円型方程式を
解く際には並列化はせず、プロセッサ間の通信量を抑えた。そのため、ほぼプロセッサ数に
比例した計算速度の上昇が実現できている。シミュレーションコードは、当研究グループに
よって今回新たに開発されたものである。

-66-

Fig.1 Left: 計算領域における座標の定義。 Right: 解析対象であるマイクロ波放電式イオン
エンジンの写真
研究成果（ Accomplishments ）:
401-ヽ
r

―
→

-+-y=O(m),experimental
y=O(m),numerical
y=O.l(m).expenmental
-- y=O.l(m).numerical
• • · y=0.2(m),expenmental
・
・
・
・y=0.2(m),numerical

-．

e-30

>

［
］ぇ
：：
：
ー
：
／
ご：ご7
:
-

： ご 零← ▼

C
l
.

.
.
・
;
1・ .
.
•.•.

'

。
鼻

I
-・:・．．．.・・
0
.
0

0
.
2

0
.
4

0.
6
z(m)

0
.
8

1
.
0

Fig. 2 プラズマ電位の計算結果と実験的計測結果の比較
Figure 2 にプラズマ電位の計算結果と、別途実施された実験的計測結果の比較図を示
す。ここで、実験計測値とはエミッシブプローブを用いた計測値である。また、 z=
0.0 m に推進機が設置されており、推進機中心軸上からの距離yでの位置におけるプラ
ズマ電位値が示されている。推進機中心軸から距離y =0.0, 0.2 m の位置では、計算値
と実験値は比較的良い一致が確認された。一方で、推進機中心軸から距離y =0.3 m の
位置では、 z軸方向で推進機に近い位置では、計算値と実験値で乖離が比較的大きい
ことが確認された。

-67-

X10-1

1.
0e+15

3

2
1.
0e+14

0 1 2︱ ︱

(E)A

1.
0e+13

1.0e+12

1.0e+11

0.
2

04

0.
6

0.
8

1.
0e+10

z(m)

Fig. 3 プラズマ密度の計算結果
Figure 3 にプラズマ密度の計算結果を示す。推進機下流直下で密度は最大であり、 np>
1015m -3であると確認された。また、推進機の両脇部分で密度は最小であり、 np<
1010m -3であると確認された。この領域は、推進機から排気されるイオンビームから
死角になる位置であることから、密度が最小であることに説明がつく。
以上の結果から、 Figure 2 で確認された y=0.3, zく 0.2 m の領域での計算値と実験
値の乖離は、実験値が妥当でないことに由来すると考えられる。なぜなら、先行研究
から、プラズマ密度が np<1011m -3では、エミッシブプローブによる電位計測は妥

1)

当性を欠くことが示されているからである。この領域は、後に議論するように、 CEX
イオンの影響が無視できない領域であり、本研究の重要な解析対象である。数値計算
による低密度領域での電位取得が、実験的計測と比較して優位であることを再確認し
た
。

Figure 4 に宇宙機表面導電部位におけるイオンエネルギー分布関数を示す。宇宙機表
面に衝突する CEX イオンの最大エネルギー値は、 40eV であり、計算領域内の最大電
位値と一致した。また、 35eV 付近にピークを持つ分布であることを確認した。
また、 Figure 5 に宇宙機表面導電部位における衝突角度分布関数を示す。ここで、衝突角
0
50゜から 90゜
度が゜で、粒子は壁面に対して垂直に衝突しているとした。衝突角度は概ね
に
分布しており、壁面に対して比較的大きな角度で衝突していることが確認された。
以上のことから、宇宙機表面導電部位に衝突するイオンは、推進機下流直下で生成され、
大きな角度を持って入射していると考えられる。

-68-

I
o
nenergyd
i
s
t
r
i
b
u
t
i
o
nf
u
n
c
t
i
o
n
1.
0
﹂ N
Slun03pazl-eEO

642
8
0000

゜

0
.
0

10

2
0
3
0
E
n
e
r
g
y(
e
V
)

40

50

Fig. 4 宇宙機表面導電部位におけるイオンエネルギー分布関数
I
n
c
i
d
e
n
tangled
i
s
t
r
i
b
u
t
i
o
nf
u
n
c
t
i
o
n
1.
0

0
﹂

0
-5

゜

e EO
﹂ N
S
J
U
n
0
:
:
Jp8Z!ー

42
86
0000
60

70

80

I
n
c
i
d
e
n
c
ea
n
g
l
ee
°
()

90

Fig. 5 宇宙機表面導電部位における衝突角度分布関数

公表状況（Publications ）:
（論文）
（口頭）
1.

Ayumu Nono, Yusuke Yamashita, Ryudo Tsukizaki, Kazutaka Nishiyama, Investigation
of ion back flow by Hybrid-PIC simulation considering experimental current density
distribution at the conductive surface for microwave discharge ion thruster, The 75th
Annual Gaseous Electronics Conference, 2022 Oct., Sendai

2.

, 下 裕 介 月, 崎 竜 童 西 山
, 和 孝 マ イ ,ク ロ 波 放 電 式 イ オ ン エ ン ジ ン に お け る バ ッ
濃野歩山
,
4
, 2023 年 1月 相
, 模原
クフローイオンに関する数値計算令和年度宇宙輸送シンポジウム

3.

, 和 孝 宇 宙, 機 プ ラ-ズ マ 間 相 互 作 用 評 価 に 向 け た
濃 野 歩 ,山 下 裕 介 月, 崎 竜 童 西 山
Hybrid-PIC 計 算, 第 36 回数値流体シンポジウム , 2022 年 12 月,オ ン ラ イ ン

-69-

1) Kemp RF, Sellen JM. Plasma Potential Measurements by Electron Emissive Probes. Rev
Sci Instrum. 1966;37: 455–461.

参考文献
2) Tani Y, Tsukizaki R, Koda D, Nishiyama K, Kuninaka H. Performance improvement of
the 10 microwave discharge ion thruster by expansion of the plasma production
volume. Acta Astronaut. 2019;157: 425–434.
^
U

3) Bird GA. Molecular gas dynamics. 1976.
4) Hara K, Mikellides IG. Characterization of low frequency ionization oscillations in Hall
thrusters using a one-dimensional fluid model. 2018 Joint Propulsion Conference.
American Institute of Aeronautics and Astronautics; 2018. ...

この論文で使われている画像