第1章

序論

1.1 背景1957 年 10 月，人類初の人工衛星であるスプートニク 1 号がソビエト連邦から打ち上げられた．それ以降，宇宙開発が活発に行われるようになり，現在では，人工衛星は，情報通信，気象観測，地球観測など様々なミッションに利用されている．
今後も様々な分野で宇宙開発が活発に行われていくと考えられる．
近年，小型衛星が注目されている．小型衛星の定義は，必ずしも明確ではないが，大別すると質量が 1[ton]より小さい衛星のことを指す．小型衛星のメリットとして，
低コスト，開発期間の短縮，ミッションの多様化などが挙げられる (1-1)．一方，搭載できる推進剤が制限されるため，推進剤が衛星の運用に大きく影響するという欠点が挙げられる．このような問題を解決するために，小型の燃費の良いスラスタが必要不可欠であり，電気推進が注目されている．中でも，イオンエンジンは，推力密度は低いものの他の電気推進と比べ比推力が高いという特徴をもつことが知られている(1-2)．実用化に当たっては，小型のイオンエンジンの寿命の認定が必要不可欠になる．
イオンエンジンの寿命認定においては，実時間寿命試験が行われてきた (1-3)．まだ，我々の記憶にも新しい小惑星探査機「はやぶさ」に搭載されたイオンエンジン
も 18,000 時間に亘る寿命評価試験を行っている(1-4)．このように，イオンエンジンの開発期間の大半を寿命評価試験に費やさねばならず，改良を含めた新機種開発が進みにくい．更なる長時間作動が達成できれば，より広範囲に亘る深宇宙探査や衛星の長寿命化が期待できるが，作動時間の長期化に伴い必要な寿命評価試験も長期化することから各種ミッションへのタイムリーなイオンエンジンの投入が妨げられていると言える(1-5)．
現在，宇宙航空研究開発機構(JAXA)により，多大な時間とコストが発生する耐久性能試験の負担を軽減することを目的とし，寿命評価試験を数値解析によって支援 す る イ オ ン 加 速 グ リ ッ ド 耐 久 認 定 用 数 値 解 析 JIEDI(JAXA Ion Engine
Development Initiative)コードの研究開発が行われている(1-5)．JIEDI コードの詳細については，3.1 項にて述べる．この JIEDI コードの抱える課題として，求解速度の向上と妥当性の検討が挙げられる(1-6)．
また，グリッドの設計に関しては，長寿命かつ推力を確保したグリッドの設計に関する明確な設計指針はなく，これまで得られた経験や，実験を行い試行錯誤的に決定されてきた(1-7)．そのため，最適なグリッドの設計についても実験による多大な時間とコストが生じている．そこで，長寿命かつ一定量の推力が確保されたグリッドの設計を短期間で行うことができれば，深宇宙探査といった様々なミッションの促進に貢献できると考えられる．

1

1.2 目的
本研究では，長寿命かつ一定量の推力が確保されたグリッドの設計をすることのできるグリッドの最適化ツールの開発を目的として，以下のことを行った．
１）JIEDI コードの求解速度の向上
２）遺伝アルゴリズムを用いたグリッド最適化ツールの開発
１）の「JIEDI コードの求解速度の向上」については，目標として，1 日 1 ケース以上の解析ができる求解速度を目指す．現状の JIEDI コードの求解速度は，使用する計算機・設定により異なるため，一概に同じということはできないが，電子の逆流が生じ，寿命を迎えるまで 10 日～2 日程度の時間を要している．この JIEDI
コードの求解速度を向上させ，グリッド最適化ツールに適用し，開発期間の短縮に貢献できる求解速度を目指す．
２）の「遺伝アルゴリズムを用いたグリッド最適化ツールの開発」については，最適化の手法として，遺伝アルゴリズムを用いることにした．遺伝アルゴリズムの詳細については 3.2 節にて述べる．遺伝アルゴリズムを用いた理由は，ニューラル
ネットワークのように教師信号を必要とせずに，適応関数を決定するだけでよいた
め，容易に用いることができ，さらに適応関数を変更するだけで任意の最適化できるため，設計者の理想に沿った設計ができ，応用性が高いと考えたからである．この遺伝アルゴリズムに JIEDI コードを組み込み，グリッド最適化ツールの開発を行
った．

2

参考文献
(1-1) 茂原正道，鳥山芳夫：衛星設計入門，培風館(2002)(1-2) 荒川義博，國中均，中山宜典，西山和孝：イオンエンジンによる動力航行，コロナ社(2006)(1-3) 早川幸男，北村正治：キセノンイオンエンジンの地上設備における運転と評
価法，航空宇宙技術研究所資料 752 号(2000)(1-4) 國中均：「はやぶさ」小惑星探査機に搭載されたマイクロ波放電式イオンエンジン，
J. Plasma Fusion Res. Vol.82 No.5(2006)， pp.300-305(1-5) 中野正勝，梶村好宏，船木一幸：JIEDI-1 の開発（イオンエンジングリッド損
耗評価コードの改良），イオン加速グリッド耐久認定用数値解析 JIEDI(JAXA
Ion Engine Development Initiatives)ツールの研究開発ワークショップ論文集(2010)，pp.1-9(1-6) 船木一幸，中野正勝，中山宜典，梶村好宏：JIEDI ツー開発上の課題と 2007年度の活動方針，イオン加速グリッド耐久認定用数値解析 JIEDI(JAXA Ion
Engine Development Initiatives)ツールの研究開発ワークショップ 2006 年度報
告書(2007)，pp.75-81(1-7) 早川幸男，北村正治，宮崎勝弘：カスプ磁場型キセノンイオンエンジンの設
計論，航空宇宙技術研究所資料 751 号(2000)

3

第2章

イオンエンジン

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

図 2-1

イオンエンジンの概念図

ここで，イオンビームの加速部について注目する．イオンビーム加速部は，2 枚または 3 枚のグリッドと呼ばれる多数の小穴の開いた電極で構成される．これらのグリッドに電位差を与えることにより，プラズマ生成したイオンを引出して加速す
る．推力は，イオンビーム排出の反作用によって得られる．数枚のグリッドの内，プラズマ生成部に接するグリッドをスクリーングリッドといい，千数百[V]の正電
4圧が加えられる．スクリーングリッドの隣に設置されるグリッドはアクセルグリッドと呼ばれ，このグリッドはマイナス数百[V]の負電圧が加えられる．この負電圧はイオンエンジン外部からの電子が加速部への逆流を防ぐ働きがある．この 2 枚のグリッド間の電位差によりイオンビームを加速する．2 枚グリッドシステムでは，
このような構造となるが，3 枚グリッドシステムでは，アクセルグリッドの隣にディセルグリッドと呼ばれる電位を 0[V]としたグリッドが加えられる．このディセ
ルグリッドは，アクセルグリッドとの間にスクリーン‐アクセル間とは逆の電位勾配を形成し，電子の逆流を防ぐ役割を持つ(2-2)．図 2-2 に軸上電位分布の概略図を
示す．

図 2-2

軸上電位分布5

2.2 イオンエンジンの寿命
イオンエンジンの寿命を決める要因は，主にグリッドの損耗と中和器の劣化である．
グリッドの損耗については，電位の関係により負の電位であるアクセルグリッド，またはディセルグリッドが損耗を受ける．アクセルグリッドの損耗による孔の広がりは，中性粒子の閉じ込め性能低下と，孔中央部の空間電位の上昇を引き起こす．
空間電位がマイナス数[V]まで上昇すれば，中和器から放出された電子がイオン源
内部に逆流する．イオンと電子の移動度の違いから，イオン源内部からイオンを引き出すよりも電子がイオン源内に進入するほうが容易であるため，この逆流電子によるイオン源‐電源‐中和器の電流ループができてしまい，イオンを引き出しにく
くなってしまう．また，極端な損耗はグリッドの機械的な破壊を引き起こし，電極間の短絡に発展する(2-3)．この原因については，次項にて言及する．
中和器の劣化については，中和器電圧の上昇のことで，この電圧がある閾値を超えると，中和器内部は自身で生成したイオンからのスパッタリングにより激しい損
耗を受け，ついには放出されるイオンと中和電子のバランスが崩壊する．この状態になると，宇宙器の電位は大きく変化し，アクセル・ディセルグリッドの過電流，グリッドや宇宙器の表面の激しい損耗を引き起こす(2-3)．
このような要因によりイオンエンジン自体の性能低下・停止を引き起こし，寿命を迎える．本研究においては，グリッドの損耗を対象とする．

2.3 グリッド損耗のメカニズム
グリッドの損耗は，イオンまたは中性粒子によるスパッタリングによって生じる．
イオンエンジンの場合，ほとんどが電荷交換イオンによるスパッタリングである．
ここで，スパッタリングとは，高速の粒子が材料表面に衝突することで，その材料を構成する粒子が叩き出される現象であり，1 衝突粒子に対するスパッタ率つまり放出粒子数は，衝突条件によって大きく異なる(2-4)．グリッドの損耗に大きく影響を与えるものは以下の通りである．各衝突の概念図を図 2-3 に示す．

2.3.1

電荷交換衝突

電荷交換衝突は，イオンと中性粒子が電荷を交換する衝突である．イオン加速部+
では(2.3.1)式で示すようにグリッド間の電位差により加速されたイオン𝑋𝑓𝑎𝑠𝑡
とそれよりも圧倒的に遅い(熱速度程度)中性粒子𝑋𝑠𝑙𝑜𝑤 が衝突して電荷を交換し，遅いイ
オンと速い中性粒子を生成する．

6

+
+
𝑋𝑓𝑎𝑠𝑡
+ 𝑋𝑠𝑙𝑜𝑤 → 𝑋𝑓𝑎𝑠𝑡 + 𝑋𝑠𝑙𝑜𝑤(2.3.1)

この反応により中性粒子となった元イオンの𝑋𝑓𝑎𝑠𝑡 は衝突前のイオンの速度で直進する．直進先にグリッドが存在する場合，グリッドに衝突してグリッド損耗の原因+となる．一方新たに電荷を得た𝑋𝑠𝑙𝑜𝑤の速度は中性粒子の熱速度程度であり，負電位であるアクセルグリッドに向けて加速され，グリッドと衝突し，グリッド損耗を引き起こす(2-5)．

2.3.2

弾性衝突

弾性衝突では，イオンと中性粒子が衝突する際にエネルギーが保存されたまま衝
突後の速度ベクトルが変化する．速度ベクトルが変化したことにより，イオンまたは中性粒子の一部がグリッドに衝突する軌道を有することになる(2-5)．
+
𝑋𝑓𝑎𝑠𝑡+ 𝑋𝑠𝑙𝑜𝑤 → 𝑋 + + 𝑋

2.3.3(2.3.2)

直接衝突

イオン加速系の設計が不適当あるいは意図的であれば，ビームの周辺部や過収束したビームがグリッド端に衝突して損耗を与える．主流イオンビームの衝突は電荷交換衝突や弾性衝突とは異なり，衝突粒子が数ケタ多く，イオンの持つエネルギーも大きいのでスパッタ量として甚大である(2-5)．

図 2-3

グリッド損耗の概要(2-4)

7

2.4 電荷交換イオンによるグリッド損耗
図 2-2 に示したように，イオンビーム加速部の孔一対について見た場合の軸上の電位分布は，スクリーングリッドからアクセルグリッドにかけて電位が下がり，電位の底を向かえて電位が上がり周囲のプラズマ電位まで回復する．電位の最小値はアクセルグリッドであるため，電荷交換イオンがアクセルグリッドに集中しやすい．
また，グリッドの損耗状況は 2 枚グリッドと 3 枚グリッドで異なる．
2 枚グリッドシステムにおけるグリッドの損耗は孔と孔の間の六角形の頂点部
分がひどく損耗を受ける“Pits & Grooves”と呼ばれる形状となる．この損耗の主な原因はアクセルグリッド下流側で生じた電荷交換イオンや，アクセルグリッド上流で生じ，アクセルグリッド下流側に回り込んだ電荷交換イオンが衝突することにより生じる．Pits & Grooves 領域は主流ビームに囲まれた領域であるため，電位が周囲と比べて低くなり，低速の電荷交換イオンを集めやすい．図 2-4 に Deep space
1 に搭載された 2 枚グリッドシステムのイオンエンジンにおける実時間耐久試験のグリッド損耗の写真を示す．図 2-4(a)が作動から 125 時間経過したものであり，図2-4(b)が 30,352 時間経過したものである．時間が経過することによって，六角形の頂点部分が損耗受けていることと，孔径が拡大していることがわかる．このようなPits & Grooves の形状の損耗が続くことで孔と孔の間がつながってしまい，グリッドの構造が損なわれることになる．
3 枚グリッドシステムでは，ディセルグリッドがアクセルグリッド下流側で役割を果たすために Pits & Grooves 状の損耗は生じない．一方でグリッド孔の周辺部が特に損耗を受ける．

(a)125 時間経過(b)30,352 時間経過
図 2-4 時間耐久試験でのアクセルグリッド損耗写真(Deep Space 1)(2-6)

8

参考文献
(2-1)栗木恭一，荒川義博：電気推進ロケット入門，東京大学出版会(2003)(2-2)中野正勝：電気推進のシステム解析と性能評価，東京大学博士論文(1999)(2-3)細田聡史，國中均：イオンエンジンによる小惑星探査機「はやぶさ」の帰還運用，J. Plasma Fusion Res. Vol.86 No.5 (2010) pp.282-292(2-4)大川恭志，早川幸男，北村正治：グリッド損耗評価試験技術の現状，イオン加速グリッド耐久認定用数値解析 JIEDI(JAXA Ion Engine Development Initiatives)
ツールの研究開発ワークショップ 2006 年度報告書(2-5) 船木一幸，篠原育，中野正勝，梶村好宏，宮坂武志，中山宜典，百武徹，和田元，剣持貴弘，村本哲也，國中均：イオンエンジングリッド損耗解析ツール
JIEDI，宇宙航空研究開発機構研究開発資料(2012)(2-6) John R. Anderson，Ira Katzt，Dan Goebel：Numerical Simulation of Two-Grid Ion
Optics Using a 3D Code

9

第3章

数値解析

3.1 イオンエンジングリッド損耗評価コード(3-1)
現在，寿命評価試験を数値解析によって支援するイオン加速グリッド耐久認定用
数値解析 JIEDI(JAXA Ion Engine Development Initiative)コードの研究開発が行われている．JIEDI コードは多孔からなるイオンエンジンの加速部の一孔に注目し，そこで引き起こされるグリッド損耗と寿命の評価を行う．
図 3-1 に JIEDI コードのフローチャートを示す．初期メッシュを生成後，中性粒子密度計算を行い，損耗計算にて，電位計算・イオンビームの軌道計算を行い，その時刻のグリッド損耗率を計算する．ここで，グリッド損耗率とは，単位時間当たりの体積損耗率のことである．そして，算出されたグリッド損耗率に適切なタイ
ムステップを掛けることにより，次の時間のグリッド形状を求める．新しいグリッ
ド形状を用いて，中性粒子密度計算・損耗計算を繰り返すことで，徐々に進行していくグリッド損耗をシミュレートすることができる．本研究では，電子の逆流が生じるまでの時間を寿命とした．

図 3-1

JIEDI コードのフローチャート

10

JIEDI コードにおける解析領域は，図 3-2 に示すように多数の孔が空いているグ
リッドの一孔部分であり，対称性を利用し一孔の 1/12 とすることができる．

図 3-2 解析領域
JIEDI コードに用いられているイオンビームの計算方式は，Flux-Tube(FT)法あるいは Streamline PIC 法と呼ばれる手法である．この方法は，流体解析において流線を求めるのと同じように，荷電粒子の流れを，電流値を持つビームとして扱う手法である．定常状態のみに着目しているために，イオンがグリッド間を通過する数
10[ns]レベル以下の非定常な現象を解析することはできないが，グリッド損耗のよ
うに数時間から数万時間のオーダーで進む現象を解析する上ではイオンの流れを定常状態として扱って問題はない．一方，より厳密な物理的再現を行っているのがParticle In Cell(PIC)法である．PIC 法はあくまで荷電粒子を粒子として捉え，その振る舞いを逐次追跡する手法である．これら 2 つのアプローチは異なるが定常解は 2
つとも同じになる．また，FT 法の特徴として，荷電粒子をビームとして捉えるため，粒子を逐一シミュレートする PIC 法に比べ計算負荷が圧倒的に軽いという点が挙げられる．したがって，同精度ならば，FT 法のほうが，圧倒的に計算量が少な
くてすみ，JIEDI コードにおいて目標とするのは，過渡現象の解析などではなく，
加速グリッドの損耗の様子であるので，FT 法を採用している．
本研究では，JAXA より提供いただいた JIEDI コードを用いてグリッドの最適
化を行う．

3.1.1

中性粒子密度計算

グリッドのスパッタをもたらす高エネルギー粒子は，電荷交換衝突や弾性衝突で生成されるが，これらは放電室より漏れ出た中性粒子とイオンビームとの衝突によって生成される．そのため，中性粒子密度と速度の決定が必要である．
11

3.1.2

損耗計算

損耗計算のフローチャートを図 3-3 に示す．データの読み込みによりメッシュ，
境界条件，中性粒子密度のデータを読み込む．その後の電位計算では，有限要素法
を用い，電荷分布がない状態を仮定してポアソン方程式を解く．
𝜌𝑖 − 𝜌𝑒∇2 𝜙 = −(3.1.1)𝜀
ここで，𝜙は電位，𝜌𝑖 はイオン電荷密度，𝜌𝑒 は電子密度，𝜀は真空の誘電率である．
次に上流境界より主流イオンビームを走らせ，イオンが作る電荷密度を計算する．
イオンビームの軌道計算は以下の式より行う．
𝑀𝒗∇𝒗 = −𝑞∇𝜙(3.1.2)

ここで，𝑀は質量，𝒗は速度ベクトル，𝑞は電荷である．イオンビームの軌道追跡の際に，そのイオンビームが要素をどれだけの時間で通過したか求めることができるので，要素番号𝑚を通過するイオンの電荷密度𝜌𝑖,𝑚 は，イオンビームの滞在時間をΔ𝑡，要素の体積を𝑉とすると，通過する全てのイオンビームについて和をとることで，以下の式より計算できる．
𝑄𝑚 ∑𝑘 𝐼𝑘 Δ𝑡𝑘,𝑚(3.1.3)
𝜌𝑖,𝑚 ==𝑉𝑚
𝑉𝑚
また，プラズマとの電位差からボルツマンの関係式を用いることで電子密度を決定する．放電室においては，(3.1.4)式を用い，下流の中和プラズマにおいては，(3.1.5)式を用いる．
𝜙 − 𝜙0(3.1.4)𝜌𝑒 = 𝜌𝑒0 ()
𝑇𝑒
𝜌𝑒 = 𝜌𝑒∞ (𝜙 − 𝜙∞)
𝑇𝑒(3.1.5)

ここで，𝜌𝑒0 は放電室内の電子密度，𝜙0 は放電室内のプラズマ電位，𝜌𝑒∞ は下流の電子密度，𝜙∞ は下流プラズマ電位，𝑇𝑒 は電子温度である．反復計算により，電位と電子密度を決定し，その後，再びイオンビームを走らせて，イオン電荷密度を求め，再度繰り返し計算により電子密度を決定する．この作業をイオンビーム軌道と電位分布が自己無頓着に決定されるまで繰り返す．
主流イオンビームの軌道が確定した後に，電荷交換イオンと中性粒子の追跡，弾性散乱イオンと中性粒子の追跡を行い，それらの粒子の衝突によるグリッドの損耗率を評価する．

12損耗率は以下のようにして算出する．グリッドに入射するイオンビーム電流を𝐼とすれば，単位時間当たりのスパッタ量𝑚̇𝑖 は，スパッタ率を𝑌，グリッドの入射エ
ネルギを𝐸，入射角度を𝜃として，
𝐼(3.1.6)× 𝑌(𝐸, 𝜃)𝑒から得られる．また，高速中性粒子の流束をΓ𝑛 とすると単位時間当たりのスパッタ
量は同様に，ṁ𝑖 =ṁ𝑛 = Γ𝑛 × 𝑌(𝐸, 𝜃)(3.1.7)となる．(3.1.6)式と(3.1.7)式によりスパッタ量が計算できれば，グリッドから放出されるグリッド材原子の流束を以下の式により，求めることが可能となる．
Γ𝑔 = ∑ ṁ𝑖 + ∑ ṁ𝑛
𝐼𝑘1
Γ𝑔 = ∑× 𝑌(𝐸𝑘1 , 𝜃𝑘1 ) + ∑ Γ𝑛𝑘2 × 𝑌(𝐸𝑘2 , 𝜃𝑘2 )𝑒
𝑘1(3.1.8)𝑘2

ここで，Γ𝑔 はグリッド材原子の流束である．これらをグリッド材の表面から各方向に重み付けして放出し，軌道を追跡することでグリッド原子がそのまま損失するのかグリッドに付着して再付着するのか判明する．放出方向の重み付けには，スパッタ原子の放出方向成分まで含む微分スパッタ率のデータが必要であるが，放出方向のデータがない場合については，スパッタ後の原子の放出方向を全方向に一様として微分スパッタ率としている．JIEDI コードにおいて，2 つのモデルが用意されており，それらを表 3.1 に示す．
表 3.1 スパッタ率モデル
モデル名

グリッド材料データ

Matsunami

Mo

𝑌(𝐸, 𝜃)，Xe，全スパッタ収量

Williams

C/C

𝑌(𝐸, 𝜃)，Xe，微分スパッタ収量

本研究においては，J.Williams のモデルを用いて計算を行っている．
放出されたグリッド原子で再びグリッドに到達したものに再付着率𝑓𝑠 を掛けることで付着量を決定する．この際，到達量の(1 − 𝑓𝑠 )は再び放出されることになるのでこれも追跡して同様の操作を繰り返すことになるが，実験により求められた再体
積率は到着後の再放出の繰り返し後の最終的なものの値であるはずなので，(1 − 𝑓𝑠 )を放出させてはいない．以上により，グリッド表面における損耗率，再付
着率が求まったので，損耗率から再付着率を引いた正味の損耗量が以下の式により求まる．計算では，グリッド表面を面で分割しているため，𝑘番目のグリッド面の正味の損耗量は，𝑘面に出入りする全ての粒子を考慮して，以下の式によって求まる．
13ṁ𝑙𝑜𝑠𝑠,𝑘 ==

∑
入射イオン

𝐼𝑖
× 𝑌(𝐸𝑖 , 𝜃𝑖 ) +𝑒Γ𝑛 × 𝑌(𝐸𝑛 , 𝜃𝑛 )∑入射中性粒子(3.1.9)−𝑓𝑠∑

Γ𝑔

入射グリッド原子

図 3-3

3.1.3

損耗計算フローチャート

グリッド損耗

損耗計算によって得られた各グリッド面の損耗率に適切なタイムステップを掛けることで，損耗体積を求める．各グリッド面の損耗体積とグリッド構成する各メ
ッシュの体積を比較し，損耗体積がメッシュ体積よりも大きい場合，そのメッシュを除去する．除去したメッシュの損耗率は隣接するメッシュに受け渡す．この操作を入力した指定時間を迎えるまで繰り返す．

14

3.2 遺伝アルゴリズム
遺伝アルゴリズム(GA : Genetic Algorithms)は，生物の進化(選択淘汰，突然変異)の原理に着想を得たアルゴリズムである(3-2)．
GA のフローチャートを図 3-4 に示す．

図 3-4

遺伝アルゴリズムフローチャート

GA による解探索は，個体と呼ばれる解の集合である母集団に対して，選択・交叉・突然変異と呼ばれる遺伝的操作を繰り返し適用することによって行われる．ある世代を形成している個体間で交叉を行い，新しい個体を生成し，突然変異により，一定の確率で遺伝子を変化させる．こうして生成された個体を評価し，その中で適応度の高い個体ほど高い確率で生き残り，次の世代を形成する．このように世代の更新を繰り返し，最適解に近づいていく．
GA の特徴として，
１）適応関数値のみを用いる直接探索法である
２）一点探索ではなく，多点探索である
３）確率的オペレータを用いる探索であるという点が挙げられる(3-3)．このような特徴を持ち合わせているため，適応度を決める情報だけが必要とされ，定式化が困難な問題や性能が明確でないものにも有効であるという利点を持つ．一方，多点探索であることによる計算負荷の増大や，確率に多く依存しているので，場合によっては時間をかけても良質な解が得られない
場合があるという欠点も持つ．

15

ここで GA に用いられる用語の定義を表 3.2 に示す．
表 3.2

遺伝アルゴリズムで用いる用語

染色体

何かしらの情報を持った個体(図 3-5)遺伝子

染色体に書き込まれた情報の構成要素(図 3-5)個体

染色体で特徴づけられた個々(図 3-5)

集団

個体の集まり(図 3-5)

集団サイズ

集団内の個体数(図 3-5)世代

遺伝的操作の前または後の集団(図 3-5)

図 3-5

遺伝アルゴリズムで用いる用語

16

3.2.1

初期染色体の発生

初期集団をランダムに発生させる．ここで，発生させる初期集団のサイズについて考えなければならない．この値は直観的に決められることが多いが，難解な問題に対しては，大きな値を採用し，そうでなければ比較的小さな値を採用する．一般
的に集団サイズを小さくすると計算コストを軽減することが可能だが，確率的揺らぎの効果が相対的に大きくなり，結果を解釈する際には注意が必要となる (3-4)．また，集団サイズや後述の交叉率・突然変異率などのパラメータの設定によって，シミュレーションの結果が大きく左右されるので，予備実験などを通して慎重に設定する必要がある(3-5)．
本研究において，GA と JIEDI コードを組み合わせて用いるために MPI による並列化を行った．この計算を行うにあたり，九州大学情報基盤研究開発センターにある高性能アプリケーションサーバ SR16000，“hakozaki”と称されるモデルの計算機を用いた．この計算機で最大限効率的に並列計算を行うができるのが 8 並列であったため，本研究においては集団サイズを 8 と設定した．本来ならば，予備実験等
により慎重に設定しなければならないが，本研究においては計算の効率を重視した．
そのため，初期集団のサイズについて今後検討する必要がある．
また，本研究における遺伝子数は 6 つとし，各遺伝子と各遺伝子の変化幅を表3.3 に示す．この変化幅は，本研究室において開発が行われている小型イオンエン
ジンの値の近傍の範囲を設定した．
表 3.3 遺伝子の変化幅
遺伝子

変化幅

スクリーングリッド厚さ ts

0.20< ts[mm]<

1.0

スクリーングリッド孔径 ds

0.80< ds[mm]<1.6

グリッド間距離 lg

0.20< lg[mm]<

1.0

アクセルグリッド厚さ ta

0.20< ta[mm]<

1.0

アクセルグリッド孔径 da

0.60< da[mm]<

1.2

アクセルグリッド印可電圧 va

-200

<

<

-100

va[V]

本研究では，実数エンコードを用いた．エンコードとは，入力信号やデータをある規則に基づいて符号化することであり，ここでは，現実の問題に対する答えである遺伝子の表現型を GA で必要な部分だけを表す情報を示す遺伝子型に変換することである．エンコードの種類として，2 進数，順列，実数の 3 種類がある．これまでに，2 進数を用いた GA に比べ，実数型 GA のほうが良好な解が得られることが報告されている(3-6)．

17

3.2.2

適応度計算

自然界においては，環境に対して適応する度合いが高い個体が生き残り繁殖する．
逆に適合度の低い個体は死滅する．GA においてもそれは同様であり，求める解に近いほど，次世代へ残る確率が高くなる．
本研究では，JIEDI コードで得られた結果を用いて適応度の評価を行う．適応度の計算式については，5 章にて述べる．

3.2.3

選択

この操作は自然界における適者生存を模倣したもので，次世代の個体を作成する
両親となる個体を選ぶ．代表的な選択として 4 つの方法が挙げられる．
a. ルーレット選択
b. ランキング選択c. トーナメント選択
d. エリート選択
以下にそれぞれについて説明する．
a. ルーレット選択
遺伝的アルゴリズムの提唱者 J. H. Holland が提案した選択手法であり，適応度を基に図 3-6 に示すようなルーレット盤を作成し，そのルーレット盤を利用してランダムに個体を選択する手法である．個体𝑖のルーレットにおける割合P𝑖 は，(3.2.1)
式に示すように，適応度の合計値に対する各個体の適応度f𝑖 によって決定される．
この手法の特徴として，適合度の低い個体も選択される可能性を持つが，適合度の低いものほど選択されにくいという点が挙げられる．

表 3.4

ルーレット選択の例

個体

適応度f𝑖

割合P𝑖 [%]

A

36

12

B

135

45

C

9

3

D

120

40

合計

300

100
図 3-6

18

ルーレット盤

p𝑖 =

𝑓𝑖
∑ 𝑓𝑖

(3.2.1)

この選択手法はシンプルであるが，適応度が負の値をとる場合，このままでは使えない点に注意が必要である．
b. ランキング選択
適応度によって順位づけを行い，あらかじめ各順位によって定められた確率によって定められた確率によって次世代に残る確率が決まる選択手法である．例を表
3.5 に示す．この選択もルーレット選択と同様に適合度の低い個体も選択される可能性を持つが，いくら適応度が低いものであっても，一定確率で選択される可能性を持つという特徴を持つ．
表 3.5

ランキング選択の例

順位

個体

適応度f𝑖

割合P𝑖 [%]

1

B

135

50

2

D

120

25

3

A

36

15

4

C

9

10

この選択手法では，適応度が負の値をとる場合でも採用することができるが，新しい世代の評価の度にソートする手間がかかる．
c. トーナメント選択
全ての個体群の中から，一定の数だけランダムに個体を選び出し，その中で最も適応度の高い個体を選択する手法である．概念図を図 3-7 に示す．ここで，「一定の数」のことをトーナメントサイズという．このトーナメントサイズが小さいほど適合度の低い個体が次世代に残る可能性があるが，適応度が高い個体でもトーナ
メントに参加できなかった場合は，選択されない恐れもある．一方，トーナメント
サイズが大きすぎると，適応度が高い個体ばかりが選択される恐れがある．

表 3.6

トーナメント選択の例
個体

適応度f𝑖

A

36

B

135

C

9

D

120
図 3-7
19

トーナメント選択概念図

d. エリート選択
最も単純な選択手法であり，適応度が最も高い個体を一定数分だけ次世代に残す選択手法である．このエリート選択された個体は，交叉・突然変異の対象とせずに，
無条件で次世代に残す．また，エリート選択のみが単独で採用されることなく，通常は他の選択手法と組み合わせて行う．
本研究においては，ルーレット選択とエリート選択を組み合わせて用いた．これ
らの選択手法を用いた理由として，シンプルな手法でかつ計算の収束性を高めるために，この 2 つの手法を用いることとした．

3.2.4

交叉

交叉とは，親となる個体間で染色体の交換・組み替えを行い，新しい個体を生成する操作である．
交叉は，選択された個体すべてに対して行うのではなく，交叉率によって定められた確率で適用される．この交叉率は，図 3-8 に示すように，選択された個体のう
ち優秀な個体の適応度の値によって異なる．適応度と適応度の平均値を比較し，平均値以下であれば，交叉率は最大となり，平均値以上であれば，線形に減少していき，適応度が最大で，交叉率が最小となる．本研究においては，交叉率の最大確率を 75%，最小確率を 50%に設定した．

図 3-8

交叉率の変化20

交叉の手法として，用いるエンコードによって手法が異なる（表 3.7）．
表 3.7

交叉の種類

エンコード

交叉方法
1 点交叉

二進数

2 点交叉
多点交叉
一様交叉
中間組み替え実数

ブレンド交叉
単峰性正規分布交叉

本研究では，実数エンコードを用いたので，実数による交叉方法について述べる．
(1)中間組み替え
この交叉は，以下の式によって行われる．
𝑣𝑐 = 𝑣𝑝1 𝑎 + 𝑣𝑝2 (𝑎 − 1)

(3.2.2)(3.2.3)𝑎 = 𝑟 + 𝑑(2𝑟 − 1)
ここで，𝑣𝑐 は次世代の個体，𝑣𝑝1 ∙ 𝑣𝑝2 は現世代の選択された個体，𝑎は(3.2.3)式で定義される -d～1+d の乱数，𝑟が 0～1 の乱数であり，𝑑は図 3-9 に示すように次世代の染色体の発生幅を決定している(3-7)．

図 3-9染色体の発生幅(2)ブレンド交叉
ブレンド交叉は，2 つ親個体によって定まる各変数の区間をα倍拡張してできる
超直方体内に子個体を一様に乱数によって生成する交叉である(3-8,9)．
(3)単峰性正規分布交叉
単峰性正規分布交叉は変数間依存性への対処を目標に設計された交叉である．こ21

の交叉では，3 つの親個体によって決まる正規乱数を用いて子個体を生成する．子個体の 2 つの親個体を結ぶ線分の周辺に正規分布に従って生成され，3 番目の親は正規分布の標準偏差を用いるために用いられる(3-8,9)．
本研究では，最も単純な手法である中間組み替えを用い，𝑑 = 0.25として交叉を行った．しかし，この手法が最も適した手法であると限らないため，今後，他の交叉手法でも計算を行い，最適なものを選択する必要がある．

3.2.5

突然変異

突然変異は，遺伝子を一定の確率で変化させる操作である．交叉と突然変異のみ
では，初期母集団の組み合わせでしか探索が行われないため，局所解に陥ると脱出が不可能になる．
突然変異も交叉同様に，すべての個体に対して行われるのではなく，突然変異率によって定められた確率で適用される．この突然変異率は，図 3-10 に示すように，個体の適応度の値によって異なる．適応度と適応度の平均値を比較し，平均値以下
であれば，交叉率は最大となり，平均値以上であれば，線形に減少していき，適応度が最大値で，交叉率が最小となる．本研究においては，突然変異率の最大確率を40%程度，最小確率を 25%に設定した．

図 3-10

突然変異率の変化

22

突然変異の手法として，摂動に類似した手法を用いた．摂動とは，予め変異させる幅を決め，その幅だけ値を変化させるものである．本研究においては，変異させる幅もランダムに発生させ，突然変異を以下の式によって行った．
𝑉𝑎𝑙𝑢𝑒 𝑀𝑢𝑡𝑎𝑡𝑒𝑑 = 𝑉𝑎𝑙𝑢𝑒 + 𝑠・𝑟(3.2.2)

ここで，𝑉𝑎𝑙𝑢𝑒 𝑀𝑢𝑡𝑎𝑡𝑒𝑑 は突然変異後の値，𝑉𝑎𝑙𝑢𝑒は突然変異前の値，𝑠はランダムに-1 か 1 のどちらかの値をとり，𝑟が変異幅であり，0.001 ～0.5 の乱数である．
本研究においてはこのような手法を用いているが，パラメータ・手法については，
今後詳細な検討が必要となる．

3.2.6

終了条件

一般的に計算を終了する条件として，いくつの種類がある．
a) あらかじめ決めておいた回数分処理を繰り返した．
b) あらかじめ決めておいた閾値を超える適応度を持つ個体が現れたc) あらかじめ決めておいた閾値を超える適応度の平均値を持つ世代(集団)が現れた
d) 世代ごとの適応度の平均値の増加率が，ある閾値以下で，かつある一定期間続いた．
などの条件が挙げられる．しかし，それぞれ次のような特徴がある．
a)の場合，終了までの時間経過がある程度まで見積もれるが，必ずしも相応しい最適解が得られる保証はない．
b)または c)の場合，想定していた最適解を得ることができるが，それが真の最適解である保証はない．また，適応度の値を高く設定しすぎると，繰り返し処理が永久に終わらない恐れがある．
d)は世代を経ることによる成長が頭打ちになることを意味する．終了までの時間
経過や最適解の水準を事前に見積もることができない．
これらのことを踏まえ，扱う問題に応じて，適宜使い分ける必要がある(3-3, 5)．
本研究における計算の終了条件は，以下のように設定した．
１）100 世代以内に 95%の遺伝子がエリートの染色体と一致した場合
２）エリートの染色体との一致の条件として，エリートの染色体との差が±0.1[mm]もしくは，10[V]以内ならば同一と判定
上記の条件より，終了判定を行った．グリッドの形状の最適化においては，最適解を想定することが難しいため，100 世代という制限を設けて，適応度が一定になるまで計算を行うことにした．

23

参考文献
(3-1)JIEDI ツール開発チーム：JAXA Ion Engine Development Initiatives(JIEDI) JIEDI-1
Version βコードサブセット版 OPTJ 解説書(3-2)北野宏明：遺伝的アルゴリズム，人工知能学会誌 Vol.7 No.1(1992)，pp.26-37(3-3) 花 田 良 子 ： 遺 伝 的 ア ル ゴ リ ズ ム ， 関 西 大 学 理 工 学 会 理 工 学 と 技 術
Vol.15(2008)，pp.13-19(3-4)古谷博史：遺伝的アルゴリズムにおける確率的揺らぎの効果，社団法人 情報処理学会 研究報告(3-5)目時伸哉：遺伝的アルゴリズム GA 入門 ～遺伝的アルゴリズムを用いたテス
トケースの生成～，(http://www.zipc.com/cesl/info/04.pdf)(3-6)小野功， 山村雅幸， 喜多一：実数値 GA とその応用，人工知能学会誌 Vol. 15
No. 2(2000)(3-7) Cody Coffman Farnell ：PERFORMANCE AND LIFETIME SIMULATION OFION THRUSTER OPTICS，Doctoral Thesis，The Colorado State University(2007)(3-8) 白川真一，矢田紀子，長尾智晴： 遺伝的プログラミングによる実数値 GA の性能差を強調する探索空間の生成，進化計算学会論 Vol.1 No.1(2010)(3-9)小野 功，佐藤浩，小林重信：単峰性正規分布交叉 UNDX を用いた実数値 GA
による関数最適化，人工知能学会誌 Vol.14 No.6

24

第4章

JIEDI コードの求解速度の向上

4.1 計算モデル
本研究室で開発が行われている小型イオンエンジンに対して JIEDI コードを適用した．本研究室では，2cm 級のマイクロ波放電式イオンエンジンの開発が行われている．表 4.1 に本計算で適用したグリッドパラメータ，表 4.2 に計算条件を示す．
また，図 4.1 に当研究室で開発が行われている小型イオンエンジン，図 4.2 に計算モデルを示す．

表 4.1 グリッドパラメータ
スクリーン

アクセル

開口率[%]

58.0

19.7

孔直径[mm]1.2

0.7

厚さ[mm]

0.3

孔の中心間距離[mm]

1.5

材料

CC

グリッド間距離[mm]

0.2

孔の数

91

表 4.2

計算条件

推進剤

Xenon

スクリーングリッド印可電圧[V]1000

アクセルグリッド印可電圧[V]-150一孔当たりの電流[A]1.319×10-4

NP/H[A/V3/2]

9.393×10-10

上流領域長さ[mm]1.2

下流領域長さ[mm]

3.27

2 価イオン比0.1

25

50 mm

50 mm
図 4.1

2cm 級のマイクロ波放電式イオンエンジン

図 4.2

解析モデル

一孔当たりの電流は，本研究室で計測された小型イオンエンジンのイオン電流値
12[mA]と孔の数から決定し，イオンビームの引き出しに関するパラメータである
NP/H は，以下の式によって決定した．
NP/H=

𝐽∙𝑙𝑒 2
𝑉0 3/2 ∙𝐷 2

(4.1)

ここで，J は一つのグリッド孔から引き出されるイオン電流，𝑙𝑒は実行電極間隙，𝑉0は引出電圧，Dは電極孔径である．上流領域の長さは，スクリーングリッドの孔径26と同程度にした．また，下流領域長さは，中和がその範囲内で十分に行われる設定する必要があるため，電位勾配が 0 になる場所までの距離の 3 倍以上に設定する必要がある．本計算における下流領域の長さについて検討した結果，3 倍では中和が十分ではなく，4 倍では計算が収束しづらくなる可能性があったため，本計算では，下流領域の長さを電位勾配が 0 になる場所までの距離の 3.5 倍に設定した．2 価イ
オン比については，損耗解析において重要なパラメータであることが報告されている(4-1)が，正確な値が明らかではないため，μ10 と同じ値に設定した．
本章において，まず求解速度向上のために，適切なメッシュサイズやグリッド形
状更新のタイムステップを調査した結果について述べる．その後，電位分布計算の処理を見直すことにより求解速度の向上を行った結果を述べる．

4.2 メッシュサイズの検討
メッシュ分割の方向として，軸方向と半径方向の 2 つが挙げられる．軸方向の分割数は上流・下流域がプラズマ状態であるため，メッシュサイズがデバイ長程度，つまり，0.0183[mm]，下流領域で 0.133[mm]程度に制限される．一方，半径方向の分割数もデバイ長程度にすることが望ましいが，特に制限はない．これは，本計算
条件において半径方向の分割数を 8 以上とすることで，メッシュサイズがデバイ長以下になるためである．しかしこれまで，適切な半径方向の分割数について言及されてこなかった．そのため，半径方向のメッシュ分割数を 16，20，24，28，32 となるモデルを用意し，それぞれに対し，計算を行い，一条件で計算が収束するまでの計算時間と初期状態でのアクセルグリッドの単位時間当たりのスパッタ原子数の変化について比較を行った．図 4.3 にその結果を示す．

図 4.3

半径方向の分割数の変化に対する計算時間とスパッタ原子数の変化
27

図 4.3 よ り， ス パッ タ 原 子 数は 半 径 方向 の 分 割 数が 24 程度 か ら 2.20 ×
11

10 [atoms/sec]に収束することが明らかになった．一方，計算時間は，分割数の増
加に伴い増えている．これは，分割数が増えることによりメッシュ数が増加し，これに伴い演算サーバの負荷も増大するためである．これらのことを踏まえ，計算時
間・精度の観点から，適切な半径方向の分割数は 24 であると考えられる．

4.3 タイムステップの検討
図 3-1 で示したように，JIEDI コードでは，中性粒子密度計算，損耗計算を行い，グリッド表面の損耗分布を求め，得られた損耗率に適当なタイムステップを乗じてグリッド表面座標を更新した後，メッシュの再構築を行うことによりグリッド損耗を逐次シミュレートする．この JIEDI コードにおいて，タイムステップについては，解析を行う人の主観に基づいて決定され，短いタイムステップを用いたシミュレー
ションについては妥当性が検討されてきた．そこで，本研究では除去割合というパラメータを導入し，最適なタイムステップの検討を行った．除去割合𝜂は以下の式
によって定義した．

η=𝜂𝑒𝑟𝑜𝑑𝑒𝑑× 100
𝜂𝑡𝑜𝑡𝑎𝑙(4.2)

ここで，𝜂𝑒𝑟𝑜𝑑𝑒𝑑 は入力されるタイムステップを基に決定される除去されるメッシュ
数，𝜂𝑡𝑜𝑡𝑎𝑙 は半径方向の分割数を 24 として計算されたアクセルグリッドのメッシュ数(648 メッシュ)とする．本計算においては，除去割合を 1.0%，2.5%，5.0%，10%，
15%とし，それぞれに対応する除去メッシュ数を決定する．表 4.3 に各除去割合に対応する除去メッシュ数を示す．
表 4.3

除去割合と除去メッシュ数

除去割合

除去メッシュ数

1.0%

6

2.5%

16

5.0%

32

10%

64

15%

97

除去されるメッシュ数は損耗計算後に得られた損耗率と入力したタイムステッ
プを基に決定され，このメッシュ数が表 4.3 に示される数と同じになるように，タイムステップを変更した．同じ除去割合で，電子が逆流するまで繰り返し計算を行
28

い，計算時間，損耗状況，各パラメータの比較を行う．計算時間としては，本条件での寿命となる電子逆流が見られるまでの計算時間の比較を行った．図 4.4 に各除去割合における総計算時間を示す．図 4.4 より，除去割合が 1.0%の場合，電子の逆流までの計算時間は 281 時間とおよそ 10 日程度の時間が必要なのに対し，除去割合を 15%とした場合 32 時間と 2 日程度となり，1.0%の場合と比較して 1/5 程度になるなど，メッシュ除去割合を増やすに従い，減少していることがわかる．

図 4.4

削減割合と電子逆流が生じるまでの計算時間の関係

図 4.5

アクセルグリッドの流入する全イオン電流の変化

29図 4.5 にアクセルグリッドに流入する全イオン電流の変化を示す．ここで全イオ
ン電流は，電荷交換イオン電流と弾性散乱イオン電流の和を示しており，大半を電荷交換イオン電流が占めている．図 4.5 より，全イオン電流は，除去割合に依存せず，時間が経過するにつれて減少している．除去割合に依存しない理由として，図4.6 に示すように電位分布に大きな違いが見られないことが挙げられる．各除去割合における電位分布に大きな違いが見られないため，流入するイオン電流に大きな違い見られなかったと考えられる．また，全イオン電流が時間の経過につれて減少する理由として，アクセルグリッドの孔部が削れることによって，孔径が広がり，衝突するイオンが減少したことが原因と考えられる．電荷交換イオン・弾性散乱イ
オンは図 4.7 に示すように主流イオンの周囲に広がって流れている．作動直後はアクセルグリッドの孔部に衝突するが，時間が経過するにつれてアクセルグリッドの孔径が広がり，衝突するイオンが減少することにより，電流値も減少したと考えら
れる．

(a)

(b)

(c)

(d)

(e)
図 4.6

10,000 時間時における電位分布
除去割合 (a)1.0%，(b)5.0%，(c)10%，(d)15%

図 4.7

イオンビームの加速領域30図 4.8 にアクセルグリッドのスパッタ原子数の変化を示す．図 4.8 に見られるようにアクセルグリッドのスパッタ原子数において，除去割合が 10%・15%の場合が異なる変化を示している．この理由として，除去割合が大きいことによって，グリッドの損耗状況が異なったことが原因だと考えられる．図 4.9 に初期状態のグリッド，図 4.10 に除去割合 1%と 15%の場合におけるグリッドの損耗状況の変化を示す．
除去割合が 1%の場合，図 4.9 の初期状態から 10,000 時間にかけてスクリーングリ
ッド上流側が損耗し，20,000 時間にかけてアクセルグリッドの下流面の孔径が広がるように削れている．それ以降は，Pits & Grooves やアクセルグリッド表面が削れていき，最終的に寿命を迎える．一方，除去割合が 15%の場合，10,000 時間付近の損耗状況と電子逆流が生じる時間を除き，同様な変化をしている．10,000 時間付近
において，スクリーングリッド上流側が損耗することに加え，アクセルグリッドの下流面の孔径が広がるように削れていることがわかる．ここで 10,000 時間付近の

図 4.8

アクセルグリッドのスパッタ原子数の変化

図 4.9 初期状態のグリッド

31

10,000[hours]

20,000[hours]

30,000[hours]

電子逆流時(a)除去割合：1%(b)除去割合：15%

図 4.10 グリッドの損耗状況の変化

32削減状況に注目する．図 4.11 に 10,000 時間付近における除去割合が 1%，5%，10%，15%の損耗状況を示す．図 4.11 より，除去割合が 1%，5%ではスクリーングリッド
上流側が損耗しているだけだが，10%，15%ではこれに加えてアクセルグリッドの孔径の下流側が広がるように削れていることがわかる．これは，損耗量が時間と共に減少して，本来ならば削れるはずがない下流側が削れた結果である．これにより，主流イオンビームの発散角が表 4.4 に示すように 1.0%・2.5%に比べ大きくなり，
結果として，電荷交換するイオンの量自体も減少し，結果として，アクセルグリッ
ドに流入するイオンがわずかに減少し，スパッタ原子数が減少したと考えられる．
また，図 4.8 のアクセルグリッドのスパッタ原子数は，除去割合が 10%・15%の場合を除き，10,000 時間付近でスパッタ原子数は一定となり，その後減少している．
このスパッタ原子数が一定となる変化は，アクセルグリッドの損耗状況によって説明できる．10,000 時間付近の損耗状況は，アクセルグリッドの孔部の損耗による孔径の変化がほとんど生じなくなる．そのため，ビームの発散角がおよそ 16.1°程度になり，全イオン電流の変化が小さくなることで，アクセルグリッドのスパッタ原
子数の変化もあまり生じずに一定となったためと考えられる．

(a)除去割合 1%(b) 除去割合 5%(c) 除去割合 10%(d) 除去割合 15%
図 4.11 10,000 時間付近の損耗状況

33

表 4.4

ビーム発散角

除去割合

発散角

1.0%

15.644°

2.5%

15.696°15%

16.639°

図 4.12 は軸上の最少電位の変化を示している．この最小電位は，イオンエンジンの性能低下を生じさせる電子逆流の発生基準であり，寿命を決定する要因となる．
本計算においては，下流の電子温度を 2[eV]と仮定しているため，グリッド孔中心部の軸上電位が-2[V]より大きくなると電子逆流が生じる設定となっている．除去割合が 15%の場合を除き，本イオンエンジンの寿命は 38,000 時間程度と変化がな
いことがわかった．しかしながら，15%の場合，45,000 時間と 7,000 時間も寿命を過大評価していることが確認された．
以上の結果より，計算時間と最適な寿命評価の観点から，除去割合が 5.0%となるタイムステップで十分であると考えられる．

図 4.12

アクセルグリッドの流入する全イオン電流の変化

34

4.4

LU 分解の改善

JIEDI コードにおいて最も計算時間を必要とする部分は，損耗計算である．この損耗計算における計算処理には，電位計算，ビーム軌道計算，有限要素法を用いるための形状関数の計算などの処理があり，損耗計算における各計算処理の割合を図4.13 に示す．図 4.13 より，損耗計算における計算時間の 80%近くを電位計算の際に用いられている LU 分解が占めていることが明らかになった．この LU 分解の処理に注目し，求解速度の改善を行う．
JIEDI コードでは，ポアソン方程式を有限要素法で解くために，以下のような式に変形を行っている(式の変形については，参考文献(3-1)を参照のこと)．ここで，𝑁は形状関数，𝜙は要素内の電位，𝑅𝐻𝑆はポアソン方程式の右辺の項を示している．

𝜕𝑁1
𝜕𝑥
∭
⋮
𝜕𝑁8
( 𝜕𝑥

𝜕𝑁1
𝜕𝑦
⋮
𝜕𝑁8
𝜕𝑦

𝜕𝑁1
𝜕𝑁1
𝜕𝑥
𝜕𝑧
𝜕𝑁1
⋮
𝜕𝑦
𝜕𝑁8
𝜕𝑁1
𝜕𝑧 )
( 𝜕𝑧

𝜕𝑁8
𝜕𝑥
𝑁1
𝜙1
𝜕𝑁8
…
𝑑𝑥𝑑𝑦𝑑𝑧 ( ⋮ ) = ∭ ( ⋮ ) 𝑅𝐻𝑆 𝑑𝑥𝑑𝑦𝑑𝑧
𝜕𝑦
𝑁8
𝜙8
𝜕𝑁8
…
𝜕𝑧 )
…

(4.3)

ここで
𝜕𝑁1
𝜕𝑁8
𝜕𝑁1
…
𝜕𝑥
𝜕𝑥
𝜕𝑧
𝜕𝑁1
𝜕𝑁8
…
𝑑𝑥𝑑𝑦𝑑𝑧
⋮
𝜕𝑦
𝜕𝑦
𝜕𝑁8
𝜕𝑁1
𝜕𝑁8
𝜕𝑧 )
…
( 𝜕𝑧
𝜕𝑧 )
𝑁1
𝑭 = ∭ ( ⋮ ) 𝑅𝐻𝑆 𝑑𝑥𝑑𝑦𝑑𝑧
𝑁8

𝜕𝑁1
𝜕𝑥
𝑲=∭
⋮
𝜕𝑁8
( 𝜕𝑥

𝜕𝑁1
𝜕𝑦
⋮
𝜕𝑁8
𝜕𝑦

(4.4)

(4.5)

と置くと
𝜙1
𝑲( ⋮ ) = 𝑭
𝜙8

(4.6)

となる．この変形後の(4.6)式の係数行列 K 対し，LU 分解を行っている．ここで，
LU 分解とは，係数行列 A を下三角行列 L と上三角行列 U の積に分解することである(4-2,3,4)．LU 分解を行った後に共役勾配法(4-2,3,4,5)を用いて，電位分布を求めている．
JIEDI コードでは，3.1.2 節で述べたように電位計算とイオンビームの軌道計算を軌道と電位分布が自己無頓着に決定されるまで繰り返している．この繰り返しにおいて，改善前のコードでは，電位計算を行うたびに LU 分解を行っている．(4.5)式からわかるように，イオンビームの軌道計算により，行列 F は毎回更新されてい35

く．一方，LU 分解を行う係数行列 K は，グリッドの形状変化がない限り変化しな
い．そのため，繰り返し LU 分解を行うことは，無駄な処理であり，LU 分解を行
う回数を削減することで，求解速度の向上ができるのではないかと考えた．このLU 分解の処理を，初回のループのみ LU 分解を行い，LU 分解した行列を保存する
ことにより，それ以降のループでは保存した行列を読み込むことで，LU 分解の回数を削減することとした．変更前・変更後のコードで初期状態まで計算を行い，計算時間を比較した結果を表 4.5 に示す．また，図 4.22 に変更後の損耗計算における計算時間の割合を示す．

図 4.11 損耗計算における計算時間の割合(変更前)

図 4.12

損耗計算における計算時間の割合(変更後)
表 4.5

計算時間
[hours]

計算時間の比較

変更前

変更後

22.4

1.64

36

表 4.5 より，変更後の計算時間が変更前と比較して，およそ 1/10 程度まで削減されていることがわかる．また，図 4.12 より LU 分解の占める割合が 10%程度まで削減されていることがわかる．以上の結果より，LU 分解の改善により，求解速度を大幅に改善することができた．また，この結果は，目標とした 1 日 1 ケースの解析に十分に達成することのできる求解速度であると考えられる．この改善されたコードを用いて，GA によるグリッドの最適化を行う．

参考文献(4-1)中野正勝：イオン加速グリッド損耗における 2 価イオンの感度解析，日本航空
宇宙学会論文集 Vol.59 No.695(2011)，pp344-351(4-2)川上一郎：数値計算の基礎，http://www7.ocn.ne.jp/~kawa1/(4-3)戸川隼人：シリーズ新しい応用の数学 17 共役勾配法，教育出版株式会社(1977)(4-4)渡部善隆：連立一次方程式の基礎知識 ~および Gauss の消去法の安定性について~，http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/MANUSCRIPT/TUT
ORIAL/leq.pdf(4-5)長谷川里美，長谷川秀彦，藤野清次：応用数値計算ライブラリ
Templates，株式会社朝倉書店(1996)

37

反復法

第5章
5.1

グリッド形状の最適化

GA と JIEDI コードの設定

本研究で設定した GA のパラメータを表 5.1 に示す．
表 5.1

GA における設定パラメータ
8

集団サイズ

6
・スクリーングリッド厚さ ts
・スクリーングリッド孔径 ds
遺伝子数・グリッド間距離 lg
・アクセルグリッド厚さ ta
・アクセルグリッド孔径 da・アクセルグリッド印可電圧 va

エンコード

実数エンコード

選択手法

エリート選択，ルーレット選択

交叉率

最大：75% 最少 50%

交叉手法

中間組み替え突然変異率

最大：40% 最少 25%

突然変異手法摂動の類似系

終了条件

100 世代以内に遺伝子の 90%がエリー
トの遺伝子と一致した場合

また，初期染色体を表 3.3 に示す範囲内で，表 5.2 に示すように以下のように発生
させた．
表 5.2
1

2

3

4

5

6

7

8

ts[mm]

0.30

0.25

0.30

0.40

0.30

0.40

0.30

0.25

ds[mm]

1.20

1.00

1.10

1.10

1.40

1.00

1.20

1.10

lg[mm]

0.15

0.16

0.14

0.13

0.18

0.13

0.20

0.16

ta[mm]

0.30

0.40

0.50

0.30

0.30

0.40

0.30

0.60

da[mm]

0.70

0.80

0.60

0.70

0.80

0.70

0.70

0.80

va[V]

-169

-158

-126

-139

-169

-169

-150

-139

38

本研究では，JIEDI コードによって得られた結果を用いて，適応度の計算を行った．JIEDI コードで計算を行うためには，遺伝子の値によって決められるスクリー
ングリッド厚さ，スクリーングリッド孔径，グリッド間距離，アクセルグリッド厚
さ，アクセルグリッド孔径，アクセルグリッド印可電圧といった 6 個のパラメータだけでは計算を実行することはできず，その他のパラメータの設定が必要となる．
そのため，表 5.3 に示すように必要なパラメータの設定を行い，すべての個体に対し，同様の設定を用いた．
表 5.3

JIEDI コードにおける設定パラメータ

スクリーングリッド印可電圧[V]1000

イオンビーム電流値[mA]

12グリッドの孔の数

91

スクリーングリッドの開口率0.58

2 価イオン比0.1

これらの値は本研究室のグリッドで設定された，または計測されたパラメータである．本来ならば，汎用性を持たせるために個体ごとに設定・検討する必要がある
が，計算が複雑になる可能性があるため，本研究では，表 5.3 のように固定し，本研究室で開発されているグリッドの最適化を行うこととする．このパラメータの設定に関しては，今後の課題とする．

5.2 初期状態の値を用いた最適化
初期状態とは，イオンエンジンの作動直後の状態のことである．この状態まで計算を行い，得られた結果を用いて最適化を行った．適応度関数の与え方は，さまざ
まな方法が考えられるが，長寿命かつ一定量の推力が確保されたグリッドを設計するために，アクセルグリッドのスパッタ原子数𝐼𝑎𝑑𝑝 と，イオンビームの下流方向の速度成分の割合を示す推力係数𝑓を用いて，以下のように設定した．

𝑓𝑖𝑡𝑛𝑒𝑠𝑠 =
𝑓=

𝑓
𝐼𝑎𝑑𝑝

∑ 𝐽𝑖 𝑣𝑖,𝑧∑ 𝐽𝑖 𝑣𝑖

(5.1)
(5.2)

ここで，𝑓𝑖𝑡𝑛𝑒𝑠𝑠が適応度，𝐽𝑖 がイオンビーム電流，𝑣𝑖 がイオンの速度を示す．(5.1)
式より，推力係数を最大で，アクセルグリッドのスパッタ原子数が最少のもの，つまり推力があり，電子逆流の要因となるアクセルグリッドの損耗が最も少ないもの
39

が優秀な個体となるように設定した．
このような設定のもとで計算を行った結果，7 世代で収束した．各世代における遺伝子の値とエリートの変化を図 5.1 に示し，収束した結果を表 5.4 に示す．また，比較として，当研究室で開発が行われているグリッドパラメータも表 5.4 に示す．

(a)スクリーングリッド厚さ(b)スクリーングリッド孔径(c)グリッド間距離(d)アクセルグリッド厚さ(e)アクセルグリッド孔径(f)アクセルグリッド印可電圧
図 5.1 各世代におけるパラメータの変化

40図 5.1 より，すべてのパラメータにおいて，世代を経るにつれて徐々に収束していっていることがわかる．また，本計算においては，非常に近い範囲内で初期値を発生させているため，早い段階で収束したと考えられる．
表 5.4 より，GA による計算結果はアクセルグリッドの厚さと印可電圧を除き，当研究室のモデルと同じパラメータとなった．アクセルグリッドの厚さは，当研究室のグリッドより薄くなり，印可電圧はより大きな負の電圧が印可される結果となった．この違いが，グリッドの寿命にどのような影響を与えたのか，それぞれのパラメータで電子逆流が生じるまで計算を行った．表 5.5 にそれぞれのパラメータで電子逆流が生じるまでの時間を示す．また，図 5.2 に電子逆流までのアクセルグリ
ッドのスパッタ原子数の変化と，軸上最小電位の変化，図 5.3 にグリッドの損耗形状の比較を示す．
表 5.4

グリッドパラメータの比較
Laboratory Model

GA Model(BOL)

スクリーングリッド厚さ[mm]

0.30

0.30

スクリーングリッド孔径[mm]

1.20

1.20

グリッド間距離[mm]

0.20

0.20

アクセルグリッド厚さ[mm]

0.30

0.20

アクセルグリッド孔径[mm]

0.70

0.70

アクセルグリッド印可電圧[V]

-150

-170

表 5.5
電子逆流が生じる時間[hour]

電子逆流が生じる時間の比較
Laboratory Model

GA Model38,00018,000(a)スパッタ原子数の変化(b)軸上最小電位の変化
図 5.2 電子逆流時までのパラメータの変化

410[hours]

8,000[hours]

12,000[hours]

18,000[hours](a)Laboratory Model
図 5.3(b)GA Model
グリッドの損耗状況の変化

42表 5.4 より，当研究室のモデルに対し，GA によるモデルは，およそ半分程度の時間で電子の逆流が生じる結果となった．図 5.2 (a)より，GA によるモデルのスパ
ッタ原子数は，当研究室のモデルに比べ小さいため，適応度は GA によるモデルのほうが高いことがわかる．図 5.2(b)より，動作開始時の軸上の最小電位に注目すると，アクセルグリッドの厚さ薄くなっている分，より小さな負の電位を形成している．また，GA によるモデルはスパッタ原子数が少ないにも関わらず，電位の上昇は当研究室のモデルと同程度の割合で生じていることがわかる．これは，当研究室のモデルよりも大きな負電圧が印加されることにより，アクセルグリッドの損耗が早い段階で生じたためと考えられる．図 5.3 より，8,000 時間程度までは両モデル
とも，スクリーングリッド上流側とアクセルの孔径が全体的に広がるように削れている．しかし，12,000 時間に達すると，本研究室のモデルではアクセルの下流側の孔径が広がるように削れているのに対し，GA によるモデルでは，すでにアクセル
グリッドの下流側の表面が削れていることがわかる．このように早い段階でアクセ
ルグリッドの損耗が生じたために，寿命が短くなったと考えられる．
以上の結果より，長寿命という観点において，GA を用いたグリッドの最適化は行われなかった．このようになった原因として，初期状態の値を用いて評価したことによって，動作後のグリッド形状の変化やパラメータの変化などが考慮されていなかったためと考えられる．表 5.5 に初期状態と損耗時における適応度を示す．推力係数に関しては，初期状態・電子逆流時ともに同程度の値となっている．一方，スパッタ原子数に関しては，初期状態・電子逆流時において大小の関係が逆転している．初期状態では，GA によるもののほうが当研究室によるものに比べ値が小さくなっているため，適応度も GA のモデルのほうが高くなっている．電子逆流時においては，当研究室のものほうが値が低くなっているため，適応度もこちらのほうが高い．この表のように，初期状態・電子逆流時で得られる値が異なるため，初期状態で評価する場合，長寿命となるグリッドを設計することが難しいことが明らかになった．そのため，長寿命グリッドを設計するために，適応度関数に電子逆流時の値を用いる必要があると考えられる．
表 5.5 適応度の比較

初期状態

電子逆流時

Laboratory Model

GA Model

推力係数

0.995

0.994

アクセルグリッド
スパッタ原子数

0.218

0.156

適応度

4.56

6.37

推力係数

0.995

0.995

アクセルグリッド
スパッタ原子数0.0793

0.105

適応度

12.5

9.48

43

5.3 電子逆流時の値を用いた最適化
5.2 節で示したように，初期状態の値を用いて評価を行った場合，長寿命という観点では最適化が行われなかった．そのため，長寿命グリッドの設計を行うために，
電子逆流が生じる時間を用いてグリッドの評価を行った．以下に本計算で設定した適応度関数を示す．
(5.3)
𝑓𝑖𝑡𝑛𝑒𝑠𝑠 = 𝐿𝑖𝑓𝑒
ここで𝐿𝑖𝑓𝑒は電子逆流が生じる時間である．つまり，電子逆流が生じる時間が最も遅いものが優秀な個体となるように設定した．
このような設定のもとで計算を行った結果，15 世代で収束した．各世代における遺伝子の値とエリートの変化を図 5.4 に示し，収束した結果を表 5.6 に示す．また，比較として，当研究室で開発が行われているグリッドパラメータも表 5.6 に示す．図 5.1 と同様に世代を経るにつれて徐々に収束していることがわかる．また，
図 5.1 においてはエリートの変化が表れているのに対し，図 5.4 では，どのパラメータにおいてもエリートの個体変化がないことがわかる．表 5.6 より，電逆流が生じる時間で評価した場合，スクリーングリッド・アクセルグリッドの孔径が大きくなり，グリッド間距離が広がり，さらに動作開始時で評価した場合と同様に印加電圧は-170[V]とより大きな負電圧が印加される結果となった．この違いがグリッドの寿命にどのように影響を与えたのか，電子逆流が生じるまで計算を行った．表
5.7 に電子逆流が生じるまでの時間を示す．また，図 5.5，図 5.6 にグリッドの損耗形状の比較を示す．

表 5.6

グリッドパラメータの比較
Laboratory Model

GA Model(EOL)

スクリーングリッド厚さ[mm]

0.30

0.30

スクリーングリッド孔径[mm]

1.20

1.40

グリッド間距離[mm]

0.20

0.40

アクセルグリッド厚さ[mm]

0.30

0.30

アクセルグリッド孔径[mm]

0.70

0.80

アクセルグリッド印可電圧[V]

-150

-170

44

(a)スクリーングリッド厚さ(b)スクリーングリッド孔径(c)グリッド間距離(d)アクセルグリッド厚さ(e)アクセルグリッド孔径(f)アクセルグリッド印可電圧
図 5.4 各世代におけるパラメータの変化

表 5.5
電子逆流が生じる時間[hour]

電子逆流が生じる時間の比較
Laboratory Model

GA Model38,00086,000450[hours]

10,000[hours]

20,000[hours]

30,000[hours](a)Laboratory Model
図 5.5(b)GA Model
グリッドの損耗状況の変化

46

(a)Laboratory Model
図 5.6(b)GA Model
電逆流時における損耗状況(a)スパッタ原子数の変化(b)軸上最小電位の変化
図 5.7 電子逆流時までのパラメータの変化
表 5.5 に示すように，当研究室のモデルに対し，GA によるモデルは 2 倍程度の時間で電子の逆流が生じる結果となり，長寿命という観点でグリッドを設計するこことができた．図 5.5，図 5.6 に示した損耗状況を比較すると，大きな違いはあまり見られないが，寿命は 2 倍程度伸びている．この原因として，アクセルグリッドの厚さは変化せず，より大きな負の電圧が印加されたことが原因であると考えられる．図 5.7 に電子逆流までのアクセルグリッドのスパッタ原子数の変化と，軸上最小電位の変化を示す．図 5.7(a)より，スパッタ原子数は動作開始時においては，GAによるモデルのほうが小さくなっているが，電子逆流時においては同程度の値となっている．また，スパッタ原子数の変化は，当研究室のグリッドに比べ，緩やかに減少している．これは，スクリーングリッドの孔径が増加することにより，中性粒
子密度が減少し，電荷交換イオンなどのグリッドの損耗に影響する要因が減少した
ため，スパッタ原子数が動作開始時において小さく，その後，緩やかに減少したと考えられる．図 5.7(b)より，動作開始時に注目すると，アクセルグリッドの印加電圧が大きい分，当研究室のモデルと比べ 2 倍程度，電位の谷は深くなっている．以上の結果より，寿命が 2 倍程度延びた原因として，スクリーングリッドの孔径が広がることによるスパッタ原子数が減少したこと，そしてアクセルグリッドにより大きな負の電圧が印加されたためと考えられる．
47

以上の結果より，電子逆流までの時間を用いて評価することにより，長寿命グリッド設計をすることができた．しかし，表 5.3 に示すようにパラメータを固定して計算しているため，正確な計算は行われていない．そのため，今後，パラメータを正しく設定し，得られたパラメータが最適解であるかなども検討する必要がある．

48

第6章

結論

本研究では，長寿命かつ一定量の推力が確保されたグリッドを設計するためのグ
リッド最適化ツールの開発を行った．最適化ツールとして，遺伝アルゴリズムを採用し，この遺伝アルゴリズムに JIEDI コード組み込み，グリッドの評価を行った．
1) 遺伝アルゴリズムに JIEDI コードを組み込むために，求解速度の向上を行った．
・計算時間と精度の観点から，半径方向の分割数は 24 程度，除去割合が 5.0%となるタイムステップが適切であると考えられる．また，本計算では，グリッドの計算自体は規格化して行っているため，形状によらずこのタイムステ
ップとメッシュサイズの分割数は汎用性を持つ．
・LU 分解の処理を改善することにより，計算時間を従来の 1/10 程度に削減することができた．この求解速度の向上により，1 日 1 ケース程度の解析が可能となった．
2) 遺伝アルゴリズムを用いたグリッド最適化ツールの開発を行った．
・JIEDI コードで初期状態まで計算を行った結果を用いて，GA によるグリッドの最適化を行った．その結果，BOL では確かに最適化できた．しかしながら，長寿命という観点からグリッドは最適化されなかった．この原因として，初期状態の結果を用いて評価したために，時間経過によるグリッドの形状の変化やパラメータの変化を考慮できなかったからである．
・JIEDI コードで電子逆流まで計算を行った結果を用いて，GA によるグリッド
最適化を行った．電子逆流が生じる時間を用いて評価を行い，当研究室で開発が行われているグリッドの 2 倍程度の寿命を持つグリッドのパラメータを決定することができた．しかし，この得られたパラメータが最適解なのか検討する必要がある．
また，その他の課題として，・収束性を向上させるために，GA における交叉手法・突然変異手法・パラメー
タの検討
・適応関数の設定の違いによる最適化の違いの検討
・GA に用いる際に設定した JIEDI コードのパラメータの検討・得られた解の妥当性検討が挙げられる．

49謝辞
まず始めに，宇宙推進という分野の研究に携わる機会を与えて頂いた中島秀紀教
授に心より感謝致します．研究を進めていく中で，的確なアドバイスや考え方等を，いつも親しみやすく丁寧にご指導くださいました山本直嗣准教授に深く感謝致します．
研究を行うにあたり，JIEDI コードに関し様々なアドバイスをしていただいた都
立産業技術高等専門学校の中野正勝准教授に深く感謝いたします．
事務手続きなどでお世話になった大神めぐみ氏，馬渡隆子氏に感謝致します．そして，ともに研究に励んだ秀平靖磨氏，日永智之氏，平野賢治氏，廣池匠哉氏に深く感謝いたします．研究室をいつも楽しく笑いのある空間にしてくれた前野旭弘氏，
笹川裕太郎氏，富永宙志氏，豊田裕司氏，木邊厚視氏，楠田将司氏，湯上尚之氏に感謝致します．また，卒業された先輩方やこの場所で研究を行うにあたり，お世話になった多くの方々に感謝致します．

50