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

修 士 論 文論文名MPS における
磁気インフレーションに関する研究

氏

名

指導教員名

野田

中島

賢治

秀紀

教授

目次

第1章

序論

1.1

太陽風を利用した深宇宙探査…………………………

1

1.2

磁気セイル(Magnetic Sail)…………………………

3

1.3

磁気プラズマセイル(Magneto Plasma Sail)…………………………

3

1.3.1

太陽風と磁場の相互作用…………………………

4

1.3.2

磁気インフレーション…………………………

5

1.4

これまでの研究…………………………

6

1.5

研究目的

…………………………

8

第2章

ハイブリットコード

2.1

基礎方程式

………………………… 10

2.2

計算方法

………………………… 12

2.2.1

CAM-CL (Current Advance Method and Cyclic Leapfrog)法

2.2.2

ラプラス方程式………………………… 142.2.3

プラズマ粒子の運動………………………… 15

2.2.4

PIC(Particle In Cell)法………………………… 19

2.2.5

Digital Filter………………………… 21

…… 12

2.3

規格化

………………………… 22

2.4

ハイブリッドコードに課せられる制限………………………… 22

第3章

磁気インフレーション解析

3.1

2 次元球座標格子………………………… 24

3.2

全方位プラズマ噴射………………………… 26

3.2.1

解析モデル

………………………… 26

3.2.2

結果

………………………… 27

3.2.3

考察

………………………… 30

極方向プラズマ噴射………………………… 43

3.3

3.3.1

解析モデル………………………… 44

3.3.2

結果と考察

………………………… 44

第4章

磁気セイル推力解析

4.1

解析モデル

………………………… 59

4.2

抗力係数

………………………… 60

4.3

解析結果

………………………… 61

結論………………………… 68

第5章
参考文献

………………………… 70

謝辞………………………… 72

第1章

序論

1.1 太陽風を利用した深宇宙探査
宇宙探査によって得られる情報は地球や太陽系、さらには宇宙の誕生と進化を解明するために非常に有益なものである。宇宙探査はアメリカ・ロシア・ヨーロッパ、
近年では中国・インドなどによって活発に行われている。日本でも、地球と火星の間を円軌道で太陽を周回している小惑星 ITOKAWA に向け、世界初の試料採取を目的とした探査機「はやぶさ(図 1-1)」が 2005 年 5 月に打ち上げられた[1]。その後、
この「はやぶさ」は 2005 年 9 月に ITOKAWA に到着し、現在は帰還ミッションの過程にある。
日本の「はやぶさ」やアメリカの「Deep Space1(図 1-2)」など数多くの探査機が運用されているが、深宇宙探査において探査機は,巨大なペイロード比(初期質量と燃料を使い果たしたロケットの終期質量の比)や探査期間の短期化が求められるた
め,現在の推進システムでの深宇宙探査は困難であるといえる。

図 1-1

図 1-2

はやぶさ

(Copyright JAXA)

Deep Space1(Copyright NASA)

現在、深宇宙探査のため多くの新たな推進システムが考案されている。そのような次世代の推進システムの候補として挙げられるのが本研究の研究課題でもある
磁気セイルである。磁気セイル(Magnetic Sail)は 20 世紀の終盤に Zubrin によって考案された太陽風を利用する推進システムである。探査機から巨大な超伝導コイルを展開し電流を流すことで磁場を発生させ、その磁場で太陽風を受けることで推力を1

得て宇宙を航行する [2]。磁気セイルは太陽風を利用するために原理的に推進剤を必要としない。また、磁場を広範囲に展開させることができれば大きな推力を得ることができるために深宇宙探査に最適であると考えられる。しかし、磁気セイルを実現するためには、工学的に困難な巨大な超伝導コイルを必要とする。よって磁気セイルは実用化にいたるまでの研究はなされなかった。この磁気セイルの考えを更に発展させ、欠点を補った推進システムが 2000 年に Washington Univ.の Winglee によって考案された磁気プラズマセイル(MPS)である[3]。磁気プラズマセイルは巨大な磁場を機械的に生成するのではなく、探査機からのプラズマ噴射によって小規模な磁場を広範囲に展開させるという点で磁気セイルより現実的な推進システムであるといえる。図 1-3 に比推力を横軸に、推力電力比を縦軸にとり各電気推進とMPS の想定される性能領域をプロットしたグラフを示す[4]。ただしここに示す
MPS の性能は、理想的な条件のもと簡単な仮定をおいて概算して求められたものである。これらから、磁気プラズマセイルは深宇宙探査に有効な推進システムとして期待されている。次節では磁気セイルと磁気プラズマセイルについて詳しい説明
を行う。
.

1000

Electric Propulsion
Efficiency=100%

Resist Jet

推力電力比（ mN/kW)

MPS

100

DC arcjetHall thruster
MPD arcjet
Ion Thruster

101000 2000 3000

4000

比推力(s)
図 1-3

MPS 予測性能[4]

2

5000

1.2 磁気セイル(Magnetic Sail)
磁気セイルの推力発生原理は太陽風と地球磁気圏の相互作用に基づいて考えられている。
太陽コロナのプラズマは 100 万度という高温になるため、静的な釣り合い状態をとることができず、絶えず外向きに加速されて超音速の流れとなり、惑星間空間に向かって吹き出される[5]。これが太陽風と呼ばれる。地球近傍における太陽風のプラズマ密度は 5.0 10 6 m 3 、プラズマ速度は 4.0 10 5 m s であると観測されている[6]。コロナや太陽風のプラズマの電気伝導度は非常に高く、多くの場合、電気伝導度は無限大として考えられる。そのため、コロナ上空まで延びた磁力線は、太陽風に運ばれて惑星間空間へと広がる。これより、地球近傍にもこのような惑星間磁場が存在するがその磁場強度は平均して数 nT ほどと極めて小さい。コロナの磁力線が開いたところから吹き出してくる太陽風は定常的な惑星間空間の構造をつくるが、磁力線が閉じた領域では、太陽フレアやコロナガス噴出と呼ばれる爆発的な活動現象が発生する。この場合、太陽風は通常よりも高い密度と速度を持つこともある。しかし、この点に関しては理論的にシミュレーションで予測することは困難なため、本研究では太陽風を定常的なプラズマ流れとして取り扱う。
地球の磁場は、太陽風に対してあたかも壁のごとく立ちはだかり、太陽に面した側の地球磁場は圧縮され、流れの圧力とつり合う。地球磁場は、ほぼ球形の空間に閉じ込められて磁気圏を形成する。太陽風は進路を曲げられ、上下や脇にそれて夜側へと流れていく。このように地球磁気圏は太陽風から粒子の運動量およびエネル
ギーの注入を受けることになり、これが探査機への推力に利用できるとして考えられるようになった。具体的には、探査機に設置したコイルにより磁場を生成し、その磁場と太陽風の相互作用によって得られる太陽風の運動量変化を推力としようと考えたのである。このような推進システムが 1990 年に Zubrin によって考案された磁気セイル(Magnetic Sail)である[2]。磁気セイルの推力は単純に、(推力)= C d ×(太陽風エネルギー)×(磁気圏投影面積)で考えられる。この式から、数 N の推力を得るのは半径数 km 級の超伝導コイルが必要であるという概算がなされ、磁気セイルは実現不可能な推進システムとされた。

1.3 磁気プラズマセイル(Magneto Plasma Sail)
巨大な超伝導コイルが必要であるという磁気セイルの欠点を磁気インフレーシ
ョンによって解決しようという推進システムが 2000 年に Winglee によって考案さ3

れた M2P2(Mini Magnetospheric Plasma Propulsion)である[3]。M2P2 では推力発生機
構は磁気セイルと同様だが、磁場がプラズマ流に凍結して運ばれる性質を利用して探査機から噴射したプラズマで磁場を遠くまで展開する(磁気インフレーション)という点に関して磁気セイルと異なる。磁気インフレーションを用いることにより、初期磁場は小規模でよく、巨大な超伝導コイルは必要としなくなった。Winglee ら
によって M2P2 と名付けられたこのシステムだが、我々はプラズマを伴う磁気セイルであることから MPS (Magneto Plasma Sail)と呼ぶ。
MPS は主に大きな 2 つの研究内容に分けられる。推力発生機構である太陽風と磁場の相互作用と推力増大機構である磁気インフレーションである。次節ではこの
2 つの研究内容について詳しく述べる。

1.3.1 太陽風と磁場の相互作用
一般に、プラズマの質量密度を 、プラズマ速度 v を、磁場を B とすると、プラズマの動圧と磁場の磁気圧の比を表す無次元のパラメータ は、
v2 2
B2 2 0である。磁気セイル及び MPS では(1.1)

1 付近で磁場との相互作用により、その方向を大きく曲げられる。そのときの運動量変化が磁場を介して探査機へと伝わり、推進力となる。プラズマを太陽風とし前述の太陽風のパラメータを代入すると、となるとき B

1

41 nT である。

推力発生メカニズム
磁気セイル及び MPS では、1 となる磁気圏に誘導電流 j1 が発生する(図 1-4)。
B j1 0(1.2)

x 、s をそれぞれコイル、電流の位置ベクトル、0を真空の透磁率とする。この誘導電流を j s とすると、ビオ・サバールの法則からコイル位置 x に次のような磁界
B x がつくられる。
Bx

0

4

js

x s
x s

3

ds

(1.3)

この磁界 B x と初期のダイポール磁場をつくる際にコイルに流れている電流 j 2 によりローレンツ力 F x を得る(図 1-5)。
F x j2 B x
これをコイルにそって積分した値が推進力 F となる。
4

(1.4)

F

(1.5)

F x dxここで、磁場境界では太陽風は曲げられて運動量が変化する。境界は太陽風の衝突によってできるので誘導電流 j1 は太陽風の運動量の変化に依存する。そのため太陽風の運動量変化によっても、推進力を見積もることができる。このことは MHD コ
ードを用いた Nishida による解析によって確認されている[7]。

図 1-4

図 1-5

磁気圏境界の誘導電流

コイルに働くローレンツ力

1.3.2 磁気インフレーション
磁気インフレーションは磁場の凍結と呼ばれる物理現象を利用して行われる。磁
場の凍結について詳しく述べる。
磁場の凍結[8]一般化されたオームの法則とマクスウェル方程式より、次のような磁場に関する
誘導方程式が得られる。

dB
dt

2

v B

B

(1.6)

0

ここで、 は電気抵抗率を表し、 0 は真空透磁率を表す。右辺第一項は凍結に関する項であり、右辺第二項は拡散に関する項である。磁場の空間変化を特徴づける長さである L を用いると、
2

B

B
L2

(1.7)

と近似できる。(1.7)式を(1.6)式に代入すると、(1.6)式の第１項と第２項の大きさの
比は、

5

VB

v B
2

B

B
0

L2

VL

L

0

1

Rm

(1.8)

0

V は流体速度である。(1.8)式の R m を磁気レイノルズ数という。 Rm は、磁力線凍結の程度を見積もるための無次元量として用いられ、 R m成立し、 R m

1 のとき、磁場の凍結が1 のとき、磁場の変化は拡散方程式によって決定される。

1.4 これまでの研究
磁気セイルは 1990 年に Zubrin によって考案された[2]。その後、工学的な理由によって実現が困難とされたが、2000 年に Winglee によって磁気インフレーションを用いた M2P2 が考案された[3]。その後、Winglee によって M2P2 の地上実験が行われた(2001 年、2003 年)[9]。これに対し Khazanov は得られる推力が非常に微弱であるとして M2P2 の優位性を否定した[10]。日本では 2003 年から宇宙航空研究開発機
構(JAXA)を中心として Winglee が考案した M2P2 を MPS として再設計を行っている。
磁気セイルの研究に関しては、MHD・ハイブリット PIC・Full PIC を用いた解析と実験の双方から研究が行われている。Akita による Full PIC コードを用いた解析
では、磁気圏境界においてコイルを想定したループ電流周りに空間的な荷電分離が生じ、静電ポテンシャル場が形成されたことを確認した[11]。また、MHD を用いた解析では、Nishida が磁気セイルの推力発生機構の解明を、Asahi が磁気セイルの推力評価と設計をそれぞれ行った[12][13]。Fujita はハイブリットコードを用いて磁気セイルにおける抗力係数のラーマ半径代表長比 rLs / L に対する依存性を解明し定式化を行った[14]。ここで抗力係数 C d とは太陽風-磁気圏の相互作用の強さそのものを表す無次元パラメータである。抗力係数 C d の詳しい説明は第 4 章で行う。
Kojima は磁気セイルを模擬した地上実験を行い、抗力係数を用いて解析との比較を行った[15]。MHD・ハイブリット(PIC)・実験で得られた抗力係数を図 1-7 にまとめる。図 1-7 よりラーマ半径代表長比 rLs / L が小さい領域では MHD とハイブリッ
ト(PIC)の結果がよく一致している。ラーマ半径代表長比 rLs / L が大きくなるとハイ
ブリット(PIC)では有限ラーマ半径効果が顕著になり、MHD と比較して抗力係数 C dが小さくなることが分かる。また、ラーマ半径代表長比 rLs / L が大きい場合の実験結果はハイブリット(PIC)の解析結果とよく一致した。これらから、ハイブリット(PIC)の結果が実際の磁気セイルの抗力係数をよく近似できていると考えられる。
Fujita が定式化した式を以下に示す。

6

Cd

ここで、 R L

3.6 exp

0.28 R L2

RL

1

3 .4
exp
R L2

0.22
R L2

RL

1

(1.9)

rLs / L であり rLs は太陽風のラーマ半径を表す。Fujita は代表長 L として磁気圏半径を用いた。磁気圏半径とは太陽風の動圧と磁場の圧力が電磁流体的に釣り合う位置である。つまりコイル中心から B

41nT となる位置までの距離を表す。当研究室でもハイブリットコードを用いた解析を行い、(1.9)式とよく一致する結果を得た[16]。
磁気インフレーションに関しては、Nishida が MHD コードを用いて磁気インフレーションの解析を行い、r-3 で減衰していた極方向の初期磁場が磁気インフレーシ
ョン後は r-2.0 で減衰することを確認した[17]。また、Winglee によって磁気インフレーションについての実験が行われた。しかし、そこで確認された磁気インフレーシ
ョンは非常に微弱なものであった[9]。
Nagata は磁気インフレーションをともなった磁気セイルである MPS そのものについて解析を行った。Nagata は MPS の推力と磁気インフレーションで使用するプ
ラズマを直接後方噴射した場合に得られる推力とを比較し、MPS の噴射プラズマのβ値が 10-6 以下でなければ MPS は他の推進システムに対して優位性を得ることができないと結論付けた[18]。

図 1-7

効力係数[19]7

1.5 研究目的
前述したように、これまで Full PIC・ハイブリット PIC・MHD を用いた解析が行われてきた。Full PIC に関しては、メッシュサイズをプラズマのデバイ長以下にしなければいけないという制限があり、現在のところ局所的な解析しかできない。そのため、磁気セイル、磁気インフレーションの双方の研究課題について評価することが難しい。また、MHD はラーマ半径が代表長に対して無視できる場合についてのみ有効である rLs / L

0 。しかし、太陽風と磁場の相互作用に関しては、目標とされる磁気圏が数 10km なのに対して太陽風のラーマ半径は約 100km である。また、磁気インフレーションに関しては、噴射後のプラズマはコイルから離れるにしたがってラーマ半径が大きくなる。よって、MPS の解析において Full PIC・MHD
による解析は MPS の実現可能性を検証する為には不十分である。そのため本研究では、ハイブリッドコードを用いて広範囲を解析領域とし、かつ粒子性を考慮した解析を行う。
図 1-8 に(1.9)式から導出した磁気圏半径に対する推力を示す。

10000

Thrust[N]

1000
100
10
1
0.1
0.01
0.001
0

200

400

600

800

1000

Magnetosphere radius[km]
図 1-8

磁気セイル推力

この図から MPS において推力は磁気圏半径のみに依存することが分かる。MPS
では磁気圏半径が大きければ大きいほど大きな推力が得られる。特に 100km 以下の磁気圏半径の領域に関しては、磁気圏半径が少し広がるだけで急激に推力が大きくなる。このことから、磁気インフレーションにより磁気圏を少しでも拡大することは非常に重要であることが分かる。本研究では、初めに探査機を球と仮定し全方位に向けてプラズマを噴射する。前述したように Winglee が行った実験では微弱な

8磁気インフレーションしか確認されなかった[9]。これは噴射プラズマのラーマ半
径が大きく磁気インフレーションを効率的に行えなかったためと考え、今回の解析
では小さなラーマ半径が予測される条件で解析を行う。これによって、プラズマが磁場に巻きつき強い磁気インフレーションが起こることを期待する。次に、実際のMPS に装置におけるプラズマ噴射を仮定し、磁気インフレーションの解析を行う。
その際噴射プラズマエネルギーを一定とし密度と速度を変化させることによって、効率的な磁気インフレーションを如何に行うかについて議論する。
また、前述した Fujita によるハイブリットコードを用いた磁気セイルの解析は磁気圏半径が 10km から 300km における解析である[14]。磁気セイルが数 N の推力を得るのは磁気圏半径が数 10km 必要と考えられる。しかし、磁気圏半径が 10km 以
下になる可能性も考慮すると、10km 以下の磁気圏半径に関しても推力を定量的に評価することは重要なことである。したがって、本研究では Fujita が行った解析より更に小さい磁気圏半径での磁気セイルの推力の算出を行う。それによって、磁気圏半径が 10km 以下で(1.9)式が妥当であるのかについて議論を行う。

9

第2章

ハイブリットコード

この章では、次章以降の解析に使用する、ハイブリッドコードについて述べる[20]。ハイブリッドと呼ばれる理由は、本手法がプラズマ中のイオンを粒子とし、電子を慣性を無視した流体として取り扱うことに起因する。この計算手法は、系の物理現象の時間スケールが、イオンサイクロトロン周波数で特徴付けられる時間スケールよりも小さい現象を解析するために開発されたものである[21]。
このコードの特徴は、電子の運動（プラズマ振動、サイクロトロン運動）のタイ
ムスケールで起こる現象の影響を除去し、時間幅や空間幅を大きく出来る点にある。
具体的には荷電準中性条件を仮定し、イオンと電子の数密度があらゆるところで等しいと置き、プラズマ振動を除去する。さらに電磁放射等の高周波現象を無視する
ために、ダーウィン近似を用いる。また、電磁場はプラズマの運動も含めて、自己無撞着に解かれる。
以後に、本コードの解法の詳細を示す。

2.1 基礎方程式計算に使用する方程式を導出するための基礎方程式は次のように与えられる。以下の方程式は MKSA 単位系で記述されている。まず、粒子として取り扱うイオンの運動方程式は、
mi

dv i
dt

Ze E v i B

(2.1)

dx i
dt

(2.2)

vi

である。ここで v i はイオンの速度、Z はイオンの価数、e は電気素量、 x i はイオンの位置、 mi はイオンの質量、また、 E 、 B はそれぞれ電場、磁場を表している。
添え字 i はイオンを後述する e は電子を表す。
流体として扱う電子の運動方程式を電場の解法に用いる。電子の運動方程式は以下の式で表される。
n e medv e
dt

ene E v e B

Pe

(2.3)

ここで v e は電子の速度、 me は電子の質量、 ne は電子数密度、 Pe は電子圧力である。電子は慣性無しの流体としてみなすので、慣性による効果は無視し、(2.3)式の左辺は 0 にできる( me0 )。さらに電子圧力 Pe に関しては電子の状態方程式より、
Pe ne Te
(2.4)

10である。ここで Te は電子温度である。(2.4)式を(2.3)式に代入して Pe を消去する。
0
ene E v e B
ne Te
(2.5)
荷電準中性条件および、ダーウィン近似を用いたアンペアの法則により
Bp

0

Je

Ji

(2.6)

が得られる。但し B p は、プラズマにより生成される磁場で、以後特に断りのない限り、磁場 B は外部磁場 B ext とプラズマにより生成される磁場とのベクトル和であるとする( B p

B B ext )。また J i はイオン電流密度、J e は電子電流密度、 0 は真空中の透磁率である。J i はイオンの数密度と電流密度から計算することができる。J eは電子電流密度の定義より、
Je

ene v e(2.7)

である。(2.6)、(2.7)式を用いて J e を消去し v e について求めると、
1
ene

ve

1

Bp

Ji

(2.8)

0

となる。これを(2.5)式に代入して電場 E について整理する。荷電準中性条件より
ene

Zeni を考慮し、
E

1
Zeni

1

Bp

B J i B ZTe ni

(2.9)

0

を得る。但し、(2.9)式から得られる電場はプラズマが存在する領域内でのみ使用さ
れ、プラズマが存在しない真空領域での電場は後節で述べるラプラス方程式から得られる。
プラズマによる磁場を記述する方程式として、ファラデーの電磁誘導の法則
BtE(2.10)

を用いる。(2.10)式より(2. 9)式は、
B
t

1
Zeni

1

Bp

0

となる。

11

B

Ji B

(2.11)

2.2 計算方法(2.1)式から(2.10)式を数値計算によって解くためには、基礎式を差分式に直し、数値的に安定で、物理的に妥当な空間幅、時間幅で解かなければならない。
本コードでは電磁場のアルゴリズムとして、CAM-CL(Current Advanced Method
Cyclic Leapfrog)法[22]を用いた。また、座標系はカーテシアン座標系と球座標系を用いた。
この節では本コードの計算方法として、電磁場の解析手法(CAM-CL 法)、ラプラス方程式(SOR 法)、プラズマ粒子の運動、PIC(Particle In Cell)を示し、最後に Digital
Filter について説明を行う。このうち、プラズマ粒子の運動、PIC についてはカー
テシアン座標系と球座標系についてそれぞれ記述する。

2.2.1 CAM-CL (Current Advance Method and Cyclic Leapfrog)法
CAM-CL 法とは電磁場の解析手法として用いられるアルゴリズムである。まず、Current Advance Method を用いて電場計算を行い、次に Cyclic Leapfrog 法を用いて磁場計算を行う。図 2-1 に CAM-CL 法のフローチャートを示す。CAM-CL 法の利
点は電磁場の時間刻みを粒子の時間刻みに比べて小さくとれる点にある。
初期設定(磁場B0 ・粒子X 0,V0 )

n=1

・PIC法
イオン密度・電流計算
00(密度N ・電流J )・Current Advance Method
電場計算
n+1/2
(磁場E
)

・Cyclic Leapfrog法
磁場計算(磁場B1/2)

粒子計算
・PIC法
イオン密度・電流計算(粒子Vn+1,Xn+3/2・
密度Nn+3/2・電流J-,J+）

粒子計算(粒子X1/2)

密度・電流計算(密度Nn+1・電流Jn+1)

・PIC法
イオン密度・電流計算(密度N1/2・電流J+ )・Cyclic Leapfrog法
磁場計算(磁場Bn+3/2)
データ出力

図 2-1

CAM-CL 法フローチャート
12

n+1=n
n+3/2=n+1/2

Current Advance Method
Current Advance Method は電場の計算を行う手法である。まず、 J n i , nin 1 2 と後述する Cyclic Leapfrog 法から求めた B n 1 2 を用いて E* を求める。ここで、n 1 / 2 はある時刻 n から t / 2 だけ進んだ値であることを示す。
E*

1
Zenin 1 2

1

B np 1 2

Bn 1 2Jni Bn 1 2

ZTe nin 1 2(2.12)

0

この E* と J i から J in 1 2 を求める。ここで、J i とは PIC 法を用いてイオンの位置 x n 1 / 2 ,速度 v n から求めたイオン電流である。同様に、後述する J i とは PIC 法を用いてイオンの位置 x n 1 / 2 ,速度 v n から求めたイオン電流である。
Ji

xn 1 2 vn(2.13)

Ji

xn 1 2 vn(2.14)

イオン電流 J in 1 2 は以下の式から求められる。
J in 1 2

t q2 n 1 2 *n i (E
2 m

Jiv in B n 1 2 )(2.15)

この J in 1 2 から最終的に E n 1 2 を求める。
En 1 2

1
Zenin 1 2

1

B np 1 2

Bn 1 2

nn 1 2nn 3 2J in 1 2 B n 1 2

ZTe nin 1 2(2.16)0

また、 n n 1 , J n 1 はそれぞれ、
n

n 1

2
Jn 1

J

J
2

(2.17)(2.18)

で求められる。
Cyclic Leapfrog 法
Cyclic Leapfrog 法は磁場の計算を行う手法である。粒子の時間刻み幅 t を l 分割
しより細かな時間刻みで電磁場の計算を行い最終的に B n 1 2 を得る。Cyclic Leapfrog法の流れを図 2-2 に示す。

13

図 2-2Cyclic Leapfrog 法 l10

ここで、電磁場を求める式は(2.9)式, (2.10)式より、
Ek

1
Zenin

1

B pk

Bk

J in B k

ZTe nin(2.19)0

B0

Bn 1 2

B1

B0

Bk

B k -2

*
B10

B9

E0

h
2h

E k -1k

2〜10(2.20)E10hとなる。最終的に B n 1 2 は B10 と B *10 により求められる。
Bn 1 2

*
B10

B10
2

(2.21)

2.2.2 ラプラス方程式
電場を計算するために用いる式は、分母に数密度を含んでいるため、ni がゼロとなる領域ではこの式は発散してしまう。さらに、数値計算の特性上 ni があまり小さい値になると、発散もしくは計算誤差が生じる。このことを避けるために、数密度がある特定の値 nv 以上の領域をプラズマ領域、 nv 以下の領域を真空領域とする。
これまで示した電場の式(2.12),(2.16),(2.19)式はいずれもプラズマ領域におけるものであり、真空領域ではラプラス方程式を使用する。
真空領域では、ラプラス方程式

2

E

0 を SOR（Successive Over Relaxation）法を14

用いて解く。SOR における加速係数

は正方形領域の場合は次のように表される。
2(2.22)
1 sin NGMX

ここで NGMX は x、y または z 方向の格子数である。 の値は通常 1.8 から 1.9 付
近である。また、

1 の場合には、SOR 法はガウス−ザイデル法に等しい。カー

テシアン座標系における SOR 法を以下に示す。
SOR 法では、次式によって変数の収束を加速する。
E p (i, j , k )E p (i 1, j , k ) E p (i 1, j , k )E p (i, j 1, k ) E p (i, j 1, k )(2.23)E p (i, j , k 1) E p (i, j , k 1) 6
(1

) E p (i, j , k ) , ( px, y, z )次式で定義される残差 x , y , zpE p (i 1, j , k ) E p (i 1, j , k )E p (i, j 1, k ) E p (i, j 1, k )(2.24)E p (i, j , k 1) E p (i, j , k 1)
6 E p (i, j , k ) , ( px, y, z )が、あらかじめ決められた値以下になれば真空部分の電場の計算を終了する。

2.2.3 プラズマ粒子の運動
本コードでは、いくつかのイオンを超粒子と呼ばれる集団として扱う。この超粒子を個別に運動させることでプラズマ粒子の運動を追う。
（１） カーテシアン座標系(2.2)式に、leap-frog 法[23]による時間中心差分を施すと、
x in 1x in
tv in 1 / 2(2.25)となる。ここで各変数の上付き添え字は時間を表し、下付き添え字は i 番目のイオンを表す。また、 t は時間幅を表す。
(2.1)式に時間中心差分を施すと、
v in 1 / 2 v in 1 / 2
t

Ze n
E
mi

1 n 1/ 2vi
2v in 1 / 2

Bn(2.26)となる。上付き添え字が n 、 n 1 / 2 である変数は既知量である。
プラズマ粒子の運動を求めるためには、まず(2.26)式において、 v in 1 / 2 を左辺にま15とめて、未知量である v in 1 / 2 を解く。次に、求めた v in 1 / 2 を(2.24)式の右辺に代入し、新たな粒子の位置 x in 1 を求めることが可能となる。
(2.26)式を変形すると、Av n 1 / 2
となる。

、

S

(2.27)、 、 を次のように定義する。
q t n
Bx
2 mi
q t n
By
2 mi

(2.28)

q t n
Bz
2 mi
1

2

2

2

なお、ここでの E 、B は粒子位置におけるそれぞれの値を示している。これらを用いて、行列 A とベクトル S は以下のように構成される。
1A(2.29)

1
1

S

ここで、c

q t mi 、m

vixm

cE xn

v my

v zm

viym
vizm

cE yn
cE zn

v zm
v xm

v xm
v my

(2.30)

n 1 / 2 である。さらに、(2.25)式を次のように変形して、速度を求める式が得られる。
v n 1/ 2

A 1S(2.31)

1

ここで、 A は以下のように得られる。
A1

1

1

2

1

2

(2.32)
1

2（２） 球座標系
超粒子の位置を数値的に求めるアルゴリズムとして 4 次の Runge-Kutta 法[24]を用いた。
ここでは、まず極座標系における運動方程式を示し、その後 4 次の Runge-Kutta
法についての説明を行う。

16

カーテシアン座標系 x, y, z と球座標系 r , ,x r sin cos
y r sin sinの関係は、(2.33)r coszとなる。またベクトル A Ax , Ay , Az を極座標系 Ar , A , A に変換する行列は、
Ar

sin cossin sin

cos

Ax

A
A

cos cossin

cos sin
cos

sin
0

Ay
Az

(2.34)

ddrr cos cos
dt
dt
d
dr
sin sinr cos sin
dt
dt
d
dr
cos
r sin
dt
dt

d
dt
d
(2.35)
r sin cosdt

(2.33)式を時間 t で微分することによって、
dx
dt
dy
dt
dz
dt

x dr
r dt
y dr
r dt
z dr
r dt

x d
dt
y d
dt
z d
dt

xd
dt
y d
dt
z d
dt

sin cosr sin sindx dy dzから求める。
, ,dt dt dtを得る。ここで、(2.34)式から vr , v , v を v x , v y , v z

dx
dy
dz
sin sin
cos
dt
dt
dt
dx
dy
dz
cos cos
cos sin
sin
dt
dt
dt
dx
dy
sin
cos
dt
dt

sin cosvrv
v(2.36)(2.36)式に(2.35)式を代入することによって、球座標系における速度の式が導出され
る。
vr
v
v
同様にして、加速度

r

,

,

dr
r
dt
d
r
r
dt
d
r sin
dt

(2.37)
r sinを求めると、
r

r
2r

r

2

r 2 sin 2
r 2 sin cos

r

2r sin

r sin

17

2r

(2.38)cos

を得る。
イオンの運動方程式は、
mαqE v B(2.39)である。これより、
r

r

q
Er
m
q
E
m
q
E
m

v B

v B

v Br

vr B

vr B

v Br

(2.40)

ここで、(2.37)式を時間 t で微分すると、
dv r
dt
dv
dt
dv
dt

dr
r
r 2 r 2 sin 2
r
dt
dr
r
r
r
r 2 sin cos
dt
d r sin
r sin
r sin
r cos
dt

(2.41)を得る。ここで、(2.37),(2.39)式を(2.40)式に代入することにより以下の式を得る。
dv r
dt

q
Er
m

dv
dt
dv

q
E
m
q
E
m

dt

v B

2

v

v B

2

v
r

vr v

v Br

vr B

vr B

v Br

cot

(2.42)

r
v cot
r

vr

v

2

v

イオン粒子の運動方程式は(2.42)式を使って解かれる。具体的には n タイムステnップの電磁場を用いて、イオン粒子の位置 ri と速度 vi
4 次の Runge-Kutta 法を用いて求め、 vi

n 12から位置 rin 12

n 1

から、 t 後の速度 vi

n 12

を

を求める。

4 次の Runge-Kutta 法は(2.42)式から、とし、速度 vif r (v r , v , v )

q
n
Er
m

f (v r , v , v )

q
n
E
m

v Br

f (v r , v , v )

q
n
E
m

vr B

n 12

n

v B

n

n

n

n

v B

n

n

vr B

n

n

n

v

n2

v
rn

n

n

n

を求める。

18

v Br

n2

vr v

v

n

vr

n

n

v

n2

cot

rn
n
v cot
rn

n

(2.43)n

t f i (v r

n 12

ki2

t f i (v r

n 12

ki3

t f i (v rki4

t f i (v r

k i1

vi

n 12

vi

n 12

,v

n 12

,v

)

k r1 2 , v

n 12

k 1 2,v

n 12

k 1 2)n 12

kr2 2, v

n 12

k 2 2,v

n 12

k 2 2)n 12

kr3 , v
2k i 2

k i1

n 12

n 12

k 3,v

2k i 3

n 12

k 3)

ki4

n 1

r, ,

i

6

となる。イオン粒子の位置 ri(2.44)

は(2.37)式を時間中心差分することによって求められる。(2.37)式を時間中心差分すると、
vr
v
v

n 12

rn 1

rn
t

n 12

n 12

(2.45)式より、イオン粒子の位置 ri

n 1

rn 12

n 1

n

(2.45)tr n 1 2 sin

n 12

n 1

n

t

は、

rn 1

rn

n 1

n

n 1

n

tv r

n 12
n 12

tv
rn 12
n 12
tv
r n 1 2 sinとなる。ここで、 r n 1 2

rn 1 rn 2,

n 12

n 1

(2.46)n 12

n2 で求められる。

2.2.4 PIC（Particle In Cell）法
イオン数密度を求めるためには、空間内に広がっているプラズマを格子点上に割り付けなければならない。そのために、PIC（Particle In Cell）法を用いる。
今、1 つの超粒子の中に N 個のイオンが含まれている場合を考える。図 2-3,2-4に示すように、格子点内に存在する超粒子の位置によって、8 個の格子点によって囲まれた体積 V0 の立方体を、V1 、V2 、〜 V8 に分割する。カーテシアン座標系(図 2-3)において体積 V0 、 V1 、 V2 、〜 V8 は簡単に求めることができる。球座標系(図 2-4)については以下の式を用いて V0 を求める。
r dr

d

d

r 2 sin dr d d

V0
r

19

1
cos(
3

d ) cos( ) d(2.47)

V1 、 V2 、〜 V8 も同様にして求められる。
格子点 i, j , k に割り振られるイオンの数は、その格子点の対角線上にある体積 V（図 2-3 の場合は V2）を用いて、
ni (i, j , k )VN
V0(2.48)

で与えられる。すべての粒子についての和を ni についてとれば、ある格子点上の niを定義することが出来る。
N i, j , kni (i, j , k )(2.49)

イオン電流密度 J i i, j , k も同様にして求める。超粒子を速度 v i とすると、
ji (i, j , k )J i i, j , k

Ze

V
Nv i
V0
ji (i, j , k )で求めることが出来る。

図 2-3

PIC 法(カーテシアン座標系)

20(2.50)

図 2-4

PIC 法(球座標系)

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

S

M

side
i , j ,kK

corneri , j ,k6 S 12 K(2.51)

を格子点上に定義される値とし、 inew
, j , k は Digital Filter を通した新しい i , j , k を表

す。また、

side
i , j ,kとcorneri , j ,kは side term と corner term であり、それぞれ以下の式で表

される。
side
i , j ,kcorneri , j ,ki 1, j , k

i 1, j , k

i , j 1, k

i , j 1, k

i , j ,k 1

i , j ,k 1

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, ki 1, j 1, ki 1, j 1, ki 1, j 1, k

21

(2.52)(2.53)

ここで、 M 、 S 、 K の値をそれぞれ 12、4、1 とすると(2.49)式は次式となる。
1
1 side1 corner
new(2.54)
i , j ,ki , j ,ki , j ,ki , j ,k4
1248
数密度、イオン電流密度および電場を計算した後に、それぞれの値にこの Digital
Filter を作用させる。

2.3 規格化
コード内で扱われる物理量は全て、規格化された値となっている。ここでは、前項で示された電場について規格化を行う。規格化のために用いられるパラメータは、
初期磁場 B0 [T]、サイクロトロン周波数ci[1/s]、光速 c [m/s]、粒子のタイムステッ

プとなる微少時間幅 t [s]である。ここで、サイクロトロン周波数、
ZeB0
ci
mi

(2.55)

は、イオン質量 mi [kg]、初期磁場 B0 [T]、電荷数 Z 、電荷素量 e [C]からなる。
規格化する物理量として距離、速度、電磁場、電流密度、および粒子密度が挙げられる。これらは、以下のように規格化される。
r
c t
v
c
B

t ci
B0 c
c2 t3
J
Ze
3
ni c t

r~

E

v~
t ci
B0

B

tZe
mi

~
B

E

tZe
mi c

~
E
~
J
n~i

ここで、文字の頭についている波線(ティルディ)は規格化された量であることを意味する。

2.4 ハイブリッドコードに課せられる制限
ハイブリッドコードにおいて、差分化の時間幅 t 、空間幅 x は次の条件を満たす必要がある。
CFL（Courant-Friedrichs-Lewy）条件t
これは、 t , xx
VA(2.56)

0 で解が収束する条件である。 V A はアルフベン速度であり、B0(2.57)
VA
0 ni mi22

で定義される。また、付加的な数値的安定条件は
0 .2
ci t
である。ここで、

ci(2.58)

はイオンサイクロトロン周波数である。

メッシュ幅の最小値を求める条件としてxがある。この cpic(2.59)

pi

より x が小さい場合にもハイブリッドコードは発散することなく計算を行えるが、その解が正しいものである保証は無い。ここで、プラズマ周波数である。

23

pi

はイオン

第3章

磁気インフレーション解析

MPS の推力は、(推力)= C d × (太陽風のエネルギー)×(磁気圏投影面積)(3.1)

で表される。ここで C d とは抗力係数であり太陽風-磁気圏の相互作用の強さそのものを表す無次元パラメータである。抗力係数 C d の詳しい説明は第 4 章で行う。
(3.1)式の右辺の第 1 項・第 2 項に関しては人工的にはどうすることもできない。よって、MPS で大きな推力を得ようとする際は、磁気圏の面積を大きくするしかな
い。図 1-8 から分かるように MPS では磁気圏が大きければ大きいほど推力を得ることが可能となる。言い換えると、磁気インフレーションにより磁場を展開できればできるほど推力が大きくなるのである。したがってこの章では、まず探査機(コイル)に対して全方位にプラズマを噴射して磁気インフレーションの可能性についての議論を行い、その後実際の MPS を模擬したより現実的な、極方向のみのプラ
ズマ噴射について磁気インフレーションを行う。その際、噴射プラズマエネルギーを一定とし、同じエネルギーを用いた場合、密度と速度を変化させることによって、
いかに効率的な磁気インフレーションを行うかについて検討する。

3.1 2 次元球座標格子
磁気インフレーションの解析において、計算コストを軽減するために 2 次元球座
標格子を用いて解析を行った。図 3-1 に 2 次元球座標格子を示す。球座標系 r , ,において、φ方向対称と仮定したプラズマ噴射を行った。よってφ方向の対称性が
生まれ、
0

(3.2)

として 2 次元球座標格子とした。具体的な方程式は第 2 章で示した方程式を(3.2)
式を考慮して用いる。球座標では極(0, )上に座標特異点を持つ。今回の解析において、極軸上では、θ成分やφ成分の値は持たず、r 方向のみの電磁場や電流の値を持つ。したがって極上では r 方向の計算のみを行う。このような極上の座標特
異点をうまく解析するためにロピタルの定理を用いる。(2.10)式を例にとると、(2.10)式の r 方向の計算は、
Br
t

Er

cos
1
E
r sin

となる。ここでロピタルの定理、

24

1 E
r

(3.3)

1
E
0 sin

1
0 cos

E

1 E
r

2 E
r

lim

lim

E

(3.4)

を用いると、
Er

0

1 E
r

(3.5)

となる。このようにロピタルの定理を用いることによって座標特異点の問題を解消
する。

Z
θ

r

Y
φ
X

図 3-1

2 次元球座標格子

今回の解析では更に計算コストを軽減するため

0〜

2と

〜

性も考慮できるようなプラズマ噴射を行った。この対称性を考慮すると

2 の対称
0〜

2

での磁場、電場、イオン電流を
Br , B , B , Er , E , E , J r , J , J

とした際に

〜

(3.6)

2 の磁場、電場、イオン電流はBr , B , B , Er , E , E , J r , J , J(3.7)とできる。これによって計算を行う体系が 1/2 となり計算コストの軽減が可能である。
また、本解析では電磁場の境界条件として一次外挿(電磁場の勾配 0)を用いた。

25

3.2 全方位プラズマ噴射
まずは、プラズマを全方位に噴射することによって磁気インフレーションが可能であるのかを調べる。

3.2.1 解析モデル
計算体系を図 3-2 に示す。探査機を半径 1m と仮定し、探査機表面から全方位に向けて噴射する。計算パラメータを表 3-1 に示す。イオンラーマ半径を小さくとる
ことによって第 1 章で述べたような Winglee の実験[9]と異なる結果を得ることを期2)待する。また、イオン密度を変更させることによって、噴射面の赤道方向(におけるβ値を変化させている。今回はβ値を 10,1,0.1,0.01 の 4 ケースにおける解析を行った。β値はプラズマの動圧と磁場の圧力の比であり、(1.1)式で表される。
これより、磁場の圧力が強い低β(β<1)ではプラズマは広がらず磁気インフレーシ
ョンは起きないと考えられる。しかし、大きなエネルギーのプラズマでのみ磁気インフレーションが行えるとするならば、その高エネルギーのプラズマをイオンエン
ジン等のように直接利用したほうが有効である。したがって、MPS の優位性を示すためには低エネルギーのプラズマで磁気インフレーションを起こすことが必要不可欠である。よって、本解析では低βにおける解析を行い、低エネルギーのプラズマでの磁気インフレーションの可能性を調査する。また、様々なβ値で解析を行うことによって、磁気インフレーションのβ値依存性を調査する。さらに減衰率 nを定義する。減衰率 n は磁束密度が B

r n で変化することを表す。この減衰率を

磁気インフレーションにより、どの程度小さくできるかを確認する。

Ion
Ion mumber density [m-3]
Ion tnitial velocity [m/s]
Ion temperature [eV]Electron temperature [eV]
Magnetic flux density [T]Larmor radius [m]Time step (particle) [s]Time step (field) [s]Total time [s]
図 3-2

表 3-1

計算体系

26

Argon
3×1018～21
4×103
1
4
0.02
0.084
5×10-7=0.024ωci
2.5×10-8
2.5×10-3

計算パラメータ3.2.2 φ方向の成分の扱い
本解析では、プラズマを等方的に噴射した後、プラズマの密度に局所的な分布が生じる結果となった。図 3-2-3 にその分布の様子を示す（β= 10）。赤い点は、噴射プラズマ（イオン）を示している。この原因については、いまだ特定できていな
いが、φ方向の磁場成分を 0 とすることで、局所的な分布はなく、等方噴射したプラズマが定常状態になるまで解析できる結果を得ている。図 3-2-2 にφ方向磁場成分を 0 にした際の同時刻の結果を示す。本 3.2 章の解析結果はすべて、φ方向成分を 0 として解析を行った。

図 3-2-2

図 3-2-3φ方向成分＝0φ方向成分≠03.2.3 結果
図 3-3、図 3-4 に t=2.5×10-3[s]の磁力線と磁束密度を示す。図 3-3 より、磁力線
は全ケースにおいて展開されており、β値が強ければ強いほど遠くまで広がっていることが確認できる。また、磁気インフレーション後の磁力線は極付近では密に、赤道付近では疎になっていることが確認できる。これより図 3-4 のように磁束密度が磁気インフレーション前に比べ極付近では強く、赤道付近では弱くなっている。

27

初期磁力線

β=10

β=1

β=0.1

β=0.01

図 3-3磁力線(t=2.5×10-3[s])28初期磁束密度

β=10

β=1

β=0.1

β=0.01
図 3-4

磁束密度(t=2.5×10-3[s])293.2.4 考察
今回の磁気インフレーションの解析では低βを含む全ケースにおいて磁気イン
フレーション現象が確認できた。そこで、低βで磁気インフレーションが起こる原因について考察を行う。まず、βは(1.1)式を変形させて、
v2
B2

(3.8)

0

で表される。磁場に関しては、コイルからの距離 r として B

r 3 に従って減衰する。

また、プラズマを噴射後、プラズマが等速でコイルから離れると仮定すると質量密度はr 2 に従って減衰する。これより、(3.8)式の右辺の一部は、
B

2

r4

(3.9)

となる。また、噴射プラズマの電子は電子圧力によって加速される。それによってイオンと電子との間に電場ができ、その電場でイオンも加速される。したがって、イオンは噴射後加速され、速度 v は増加する。このプラズマの加速と(3.9)式によって、噴射後β値が急激に増加し、βが 1 を超えた領域から磁気インフレーションが起こると考えられる。しかし、β=0.01 のケースにおいて赤道付近のプラズマは磁場の圧力が強いため、磁場を横切って噴出していない。図 3-5 に赤道付近のイオン
粒子の速度ベクトルを示す。図 3-5 から分かるように、赤道付近で噴射されたイオ
ン粒子は低βのために磁場を横切って噴出していないが、極付近のプラズマが磁場に巻きついて赤道付近に侵入してきている。それらの赤道付近に侵入してきたプラ
ズマによって磁気インフレーションが起こったと考えられる[26]。
また、表 3-2 に各ケースでの極方向における減衰率を示す。表 3-2 よりβ値が大きい方が磁場がより拡大していることが分かる。

図 3-5

赤道付近の速度ベクトル、t=2.5e-04[s] (β=0.01)30βDecay index

10
2
表 3-2

1
2.2

0.1
2.3

0.01
2.8

Initial3

極方向の減衰率

図 3-6 a)〜図 3-6 v)にβ=1 における結果、図 3-7 a)〜図 3-7 v)にβ=0.01 のケース
における結果を示す。
図 3-6 a)〜図 3-6 v)ならびに、図 3-7 a)〜図 3-7 v)について、以下に考察する。
a) 極方向、ならびに赤道方向の磁束密度プロファイルを示す。β=1、β=0.01 のそれぞれのケースにおいて、赤道方向は初期値よりも小さくなり、極方向については、初期値よりも大きくなっている。これについては、先に述べたように、
極方向では、磁力線がインフレーションによって密になるため、また赤道方向
では、磁力線が疎になるためである。
b) 極方向、ならびに赤道方向のイオン密度のプロファイルを示す。β=1 の場合と比較してβ=0.01 の場合では赤道方向の密度は大きくなっている。つまり、低いβでは、赤道方向に噴射粒子が溜まっていることを意味している。
c) 極方向、ならびに赤道方向のイオンの平均速度のプロファイルを示す。β=1の場合と比較してβ=0.01 の場合では赤道方向の速度は小さい。
d) イオンの粒子軌跡を示す。β=0.01 の場合では、極方向へ噴射した粒子は磁力線に沿って移動する。磁力線を横切ることが出来ない赤道付近の粒子についても、磁力線に沿って移動しており、赤道方向の磁場の展開度合いが小さくなっている。
e) イオン密度のコンター図を示す。
f)βのコンター図を示す。

g) 磁場の絶対値|B|のコンター図を示す。
h) Br のコンター図を示す。
i)

Bθのコンター図を示す。

j)

磁場のベクトル図を示す。
X 軸上に、θ方向に向いた磁場が確認できるか、この磁場の絶対値は非常に小さい（すべてのベクトルは大きさによらず長さ一定のため）。

k) 電場の絶対値|E|のコンター図を示す。
l)

Er のコンター図を示す。

m) Eθのコンター図を示す。

31

電場が正負の縞模様になっている分布が見られる。この様子は、ｐ）のベクトル図でも確認できる。このような構造が現れる原因は不明である。このような電場が生じているにも関わらず、粒子軌道が影響を受けずｒ方向を向いている
理由として、電場の大きさとしては、粒子の軌道に影響を与えるほど大きくな
いということが考えられる。
n) Eφのコンター図を示す。
o) Eφ/|E|のコンター図を示す。
p) 電場のベクトル図を示す。
q) Jr のコンター図を示す。
r)

Jθのコンター図を示す。

s)

Jφのコンター図を示す。

t)

電流のベクトル図を示す。

u) イオン電流と磁場のベクトル図を示す。
v) divB/|B|のコンター図を示す。
低ベータの噴射においては、コイルから遠方へ行くに従い、ベータ値が大きくな
り、βが１を超え、展開が開始されると、噴射方向と、磁場のベクトルは同じ方向になり、結局噴射プラズマは磁場に沿う形で流れていくことになる。そのような状
況では、ラーマ半径はほぼ 0 となり、有限ラーマ半径効果については定常状態、遠方では、その影響は小さいものと考えられる。
MHD コードを用いた解析結果との違いに関して、赤道方向の磁場が異なるのであるが、それが有限ラーマ半径効果かどうかは今後検討が必要である。

32

図 3-6 a) 磁束密度

図 3-6 b) イオン密度(t=2.5e-03[s])33図 3-6 c) イオン速度(t=2.5e-03[s])図 3-6 d) イオン粒子軌跡34

図 3-6 e) イオン密度 (t=2.5e-03[s])図 3-6 f) β 値 (t=2.5e-03[s])35図 3-6 g)|B| (t=2.5e-03[s])図 3-6 h) Br (t=2.5e-03[s])図 3-6 i) Bθ (t=2.5e-03[s])図 3-6 j) r-θ 面の磁場方向 (t=2.5e-03[s])図 3-6 k) |E| (t=2.5e-03[s])図 3-6 l) Er (t=2.5e-03[s])36図 3-6 m) Eθ (t=2.5e-03[s])図 3-6 n) Eφ (t=2.5e-03[s])図 3-6 o) Eφ/|E| (t=2.5e-03[s])図 3-6 p) r-θ 面の電場方向 (t=2.5e-03[s])図 3-6 q) Jr (t=2.5e-03[s])図 3-6 r) Jθ (t=2.5e-03[s])37図 3-6 s) Jφ (t=2.5e-03[s])図 3-6 t) r-θ 面のイオン電流方向 (t=2.5e-03[s])図 3-6 u) イオン電流と磁場 (t=2.5e-03[s])38図 3-6 v) divB/|B| (t=2.5e-03[s])図 3-7 a) 磁束密度39

図 3-7 b) イオン密度 (t=2.5e-03[s])図 3-7 c) イオン速度 (t=2.5e-03[s])40図 3-7 d) イオン粒子軌跡

図 3-7 e) イオン密度 (t=2.5e-03[s])41図 3-7 f) β 値 (t=2.5e-03[s])図 3-7 g) |B| (t=2.5e-03[s])図 3-7 h) Br (t=2.5e-03[s])図 3-7 i) Bθ (t=2.5e-03[s])図 3-7 j) r-θ 面の磁場方向 (t=2.5e-03[s])42図 3-7 k) |E| (t=2.5e-03[s])図 3-7 l) Er (t=2.5e-03[s])図 3-7 m) Eθ (t=2.5e-03[s])図 3-7 n) Eφ (t=2.5e-03[s])図 3-7 o) Eφ/|E| (t=2.5e-03[s])図 3-7 p) r−θ 面の電場方向(t=2.5e-03[s])43図 3-7 q) Jr (t=2.5e-03[s])図 3-7 r) Jθ (t=2.5e-03[s])図 3-7 s) Jφ(t=2.5e-03[s])図 3-7 t) r−θ 面の電流方向 (t=2.5e-03[s])図 3-7 u) イオン電流と磁場 (t=2.5e-03[s])44図 3-7 v) divB/|B| (t=2.5e-03[s])3.3 極方向プラズマ噴射
これまでは磁気インフレーションの可能性を探るために全方位における噴射を仮定した解析を行い、磁気インフレーション現象が起こることを確認した。しかし、全方位噴射は現実的には困難である。よって今回は実際の装置を想定(図 3-8)してプラズマ噴射を極方向のみとして解析を行う。その際、噴射プラズマエネルギーを一定とし、同じ噴射エネルギーを用いた際に、どのような密度と速度で噴射すれば、効率的な磁気インフレーションを行うことができるかについて検討する。

図 3-8

MPS イメージ図

45

3.3.1 解析モデル
計算体系を図 3-9 に示す。体系の中心から 1m の位置から 100m の範囲までの解析を行う。体系の中心にコイルを置き初期磁場を生成する。プラズマ源はアルゴンを使用し、θ=0~30°から磁力線に沿って噴射する。それによって、噴射面での電場を 0 に設定する。噴射面ではプラズマが磁力線と同方向に進むため、プラズマと磁場は相互作用を起こさず、結果として噴射面付近では磁気インフレーションが起こらない。
本解析ではプラズマ噴射口における磁束密度を 0.01[T]とし、プラズマの密度と速度を変化させる。表 3-3 に計算パラメータを示す。全ケースにおいてプラズマの動圧と磁場の圧力の比であるβ値を等しく設定し同じプラズマエネルギー(β=5)
を噴射した。このようなパラメータでの磁気インフレーションの比較を行い議論する。それによって、効率的な磁気インフレーションを如何に行うかについて議論する。
なお、本章における解析では、前章で実施していたφ方向の磁場成分については
0 としていない。

Case
1
2
3
4
5
図 3-9

Velocity [m/s]
4000
8000
12000
16000
20000
表 3-3

計算体系

Density [m-3]3.75×1020
9.38×1019
4.17×1019
2.34×1019
1.50×1019

計算パラメータ

3.3.2 結果と考察
図 3-10、図 3-11 に t=0.01[s]におけるイオン粒子の分布図と磁力線を示す。図 3-10より、全ケースにおいて極方向のみに噴射していたイオン粒子は+θ方向にも移動していることが分かる。これは極付近で噴射されたプラズマの電子はその電子圧力46によって+θ方向に移動する。その電子とイオンとの間に生じた電場によってイオンが+θ方向に力を加えられるためである。図 3-10 から各ケースにおける結果を比較すると低初期速度であるほどイオン粒子がより多く赤道付近に分布していることが分かる。これは、低初期速度の方が長時間この力を受けるためである。また、
図 3-11 から全ケースにおいて磁気インフレーションが生じて磁場が拡大されていることが分かる。更に初期速度が速い方が磁場がより拡大されている。これは、磁場がプラズマに凍結しているためであり、各ケースにおける結果はそれぞれのイオ
ン粒子の動きに対応している。初期速度が速い方がより速度の r 成分を保つためそれに伴って磁力線もより拡大されている。
図 3-12 に t=0.01[s]における磁束密度分布を示す。図 3-11 から全ケースにおいて磁力線は極では密に、赤道では疎になっていることが分かる。したがって、図 3-12で磁束密度が極では強く、赤道では弱くなっている。図 3-13 に全ケースにおける極上での磁束密度の比較を示す。全ケースにおいて r=2[m]付近から磁気インフレーションが始まり磁場が拡大していることが分かる。各ケースにおける減衰率 nを図 3-14 に示す。また、各ケース間を補間した解析結果を図 3-14 に追記する。図3-14 より初期磁場では n=3 であったが磁気インフレーションによって n=2.1〜2.17に変化したことが分かる。極付近の磁力線は、初期には r 方向とθ方向に成分を持
つが、磁気インフレーション開始後、時間が経過するに従って噴射プラズマの速度
u と磁場 B が並行になる、つまり磁場成分が r 方向のみを持つ。このような場合、理論的には、極での磁束密度は中心からの距離 r の 2 乗に反比例する。よって本解析のような噴射方法を用いた場合には、磁気インフレーション後、極の磁束密度の減衰率は n=2 に近づくと考えられる。また、本解析ではプラズマ速度が速くなるにしたがって減衰率が n=2 に近づいた。これは、前述したように磁場がプラズマに凍結しているためと考えられる。
図 3-15a)〜3-15y)に Case5 における解析結果を示す。
図 3-15c) の極方向、ならびに赤道方向のイオンの平均速度のプロファイルに関して、赤道方向のイオンの速度は、極から噴射されたイオンが赤道へ達してからの速度の変化となる。コイルから遠方へ行くに従い、加速されていることが分かる。また、図 3-15d)のイオンの粒子軌跡について、本解析は、すべてβ= 5 のケースであり、磁力線にトラップされる粒子は確認できない。

47

Case1

Case2

Case3

Case4

Case5

図 3-10イオン粒子分布 (t=0.01[s])48Initial

Case1

Case2

Case3

Case4

Case5
図 3-11 磁力線 (t=0.01[s])49Initial

Case1

Case2

Case3

Case4

Case5
図 3-12磁束密度 (t=0.01[s])

50

r=2m

図 3-13極における磁束密度 (t=0.01[s])2.17

Decay index

2.16

Case1

2.15

Case2

Case3
Case4

2.14
2.13
2.12
2.11
2.1

Case5

2.09
0

5000

10000

15000

Velocity [m/s]図 3-14

51

減衰率

20000

25000

図 3-15 a) 磁束密度

図 3-15 b) イオン密度 (t=0.01[s])52図 3-15 c) イオン速度 (t=0.01[s])図 3-15 d) イオン粒子軌跡53

図 3-15 e) イオン密度 (t=0.01[s])図 3-15 f) β 値 (t=0.01[s])54図 3-15 g) |B| (t=0.01[s])図 3-15 h) Br (t=0.01[s])図 3-15 i) Bθ (t=0.01[s])図 3-15 j) Bφ(t=0.01[s])図 3-15 k) Bφ/|B| (t=0.01[s])図 3-15 l) r-θ 面の磁場方向 (t=0.01[s])55図 3-15 m) |E| (t=0.01[s])図 3-15 n) Er (t=0.01[s])図 3-15 o) Eθ (t=0.01[s])図 3-15 p) Eφ(t=0.01[s])図 3-15 q) Eφ/|E| (t=0.01[s])図 3-15 r) r-θ 面の電場方向(t=0.01[s])56図 3-15 s) Jr (t=0.01[s])図 3-15 t) Jθ (t=0.01[s])図 3-15 u) Jφ(t=0.01[s])図 3-15 v) r-θ 面のイオン電流方向(t=0.01[s])図 3-15 w) イオン電流と磁場 (t=0.01[s])57図 3-15 x) divB/|B| (t=0.01[s])図 3-15 y) コイル近傍での divB/|B| (t=0.01[s])

58

第4章

磁気セイル推力解析

前述したように Fujita による磁気セイルの解析[14]は磁気圏半径が 10km から300km における解析であるが、それ以下の磁気圏半径に関しても推力を定量的に評価することは重要なことである。したがって、本研究では Fujita が行った解析より更に小さい磁気圏半径での磁気セイルの推力の算出を行う。それによって、磁気圏半径が 10km 以下で(1.9)式との比較を行い妥当性を確かめる。

4.1

解析モデル

磁気セイル解析の計算体系を図 4-1 に示す。計算体系はカーテシアン座標系とし、YZ 面にコイルを置き初期磁場を生成する。太陽風を仮定した水素プラズマを-Z の面から Z 軸正の方向に入射した。電磁場の境界条件として境界では電磁場の勾配は0 とした。表 4-1 に計算に用いたパラメータを示す。太陽風の密度、速度、電子温
度は地球近傍の値を用いた。

ci

はイオンサイクロトロン周波数で、ここに示しているのは磁場強度が最も大きい位置における値である。また今回の解析では惑星間磁場(IMF)は極めて小さいために無視できるものとし、考慮していない。

L
図 4-1

計算体系

59

Io n

H ydro g en

-3

Io n N um ber D ensity [m ]Io n Initial V elo city [km /s]

5 ×10 6400

N um ber o f S uper P article

3000000

W eight o f S uper P article

1 0 1 7 ～ 10 2 1
10E lectro n T em perature [eV ]Δt
T o tal T ime S tep [step]

表 4-1＜ 0.01ω ci -1
50000step

計算パラメータ

4.2 抗力係数
第 1 章で述べたように磁気セイルの推力は、(推力)= C d ×(太陽風エネルギー)×(磁気圏投影面積)(4.1)と考えられていた。ここで、抗力係数(無次元推力値) C d の説明を行う。抗力係数（無次元推力値） C d は単純に(4.1)式の左辺を右辺で割った値であり、
F
Cd
1 2 v2S(4.2)

で表される。ここで 1 2 v 2 は太陽風のエネルギー、S は磁気圏の面積である磁気圏
投影面積を示し S

L2 である。抗力係数は、推力を規格化（無次元化）すること

により太陽風-磁気圏の相互作用の強さそのものを評価することができる値である。
MPS における多くの研究で代表長 L は磁気圏半径、つまりコイル中心から太陽風の動圧と磁場の圧力が電磁流体的に釣り合う位置までの距離 rC とする。Fujita によって定式化された(1.9)式をもう一度示す。

Cd

3.6 exp

0.28 R L2

RL

1

3.4
exp
R L2

0.22
R L2

RL

1

RL

rLs / L

(4.3)

ここで、 rLs は磁気圏付近における太陽風のラーマ半径を表す。ラーマ半径 rL はこの磁場強度と太陽風代表速度 400km/s から、磁気圏周辺でおよそ 100km と計算することができ、この値を参考に rLs / L は求められる。
前述したように Fujita は 10km から 3000km までの磁気圏半径についての解析を行った。よって、本解析では 10km 以下について解析を行い(4.3)式との比較を行う。
各ケースのパラメータを表 4-2 に示す。Case1〜Case3 までは Fujita と同じ領域での60

解析でありこの領域でコードの妥当性を確認する。Case4,Case5 に関しては Fujitaが行ったものより小さい磁気圏半径での解析である。
Case

1
2
3
4
5

L[km]

200
80
40
4
0.4

Coil Radius[km]Coil Current[kA]

125
50
25
2.5
0.25

表 4-2

5.25
2.1
1.05
0.105
0.01

rLs/L

0.5
1.25
2.5
25
250

System Size[km]

1000
400
200
20
2

磁気セイル解析ケース

推力 F の算出方法についての説明を行う。第 1 章で示したように、本来、磁気セ
イルは電磁気学的相互作用によってコイルがローレンツ力を受けて加速する。しかし、Nishida によってそのローレンツ力を磁気圏周辺の太陽風の運動量変化によって近似できるということが示されている[7]。本解析においては、コイル半径に対して十分なメッシュ数を用いていないため、上記のような電磁気学的な力を直接算出することは困難である。そこで、計算体系内の太陽風の運動量変化を推力として見積もることとしている。

4.3 解析結果
図 4-2 から図 4-4 に Case1 から Case3 の結果を示す。それぞれ X-Z 面における図であり、添え字の a,b,c はそれぞれ定常状態でのイオン粒子分布、イオン密度分布、
磁束密度分布を表す。イオン粒子に関しては、コイルの中心から±Y 方向に体系の1/4 の領域を切り取った結果である。Case4,Case5 についてはイオン粒子、イオン密度ともにほぼ変化は見られないほど相互作用は小さい。よって図 4-5、図 4-6 に Case4と Case5 の磁束密度分布のみ示す。図 4-2〜図 4-4 から分かるように、ラーマ半径
代表長比 rLs / L が小さくなるにつれてイオン粒子が磁場によって偏向されにくくなっていることが分かる。それにより、磁気圏境界のイオン密度も低くなっている。
図 4-2〜図 4-6 から磁束密度はラーマ半径代表長比 rLs / L が小さくなるにつれ変化
が小さくなる。
図 4-7 に本解析で得られた抗力係数 C d を(4.3)式と比較した図を示す。まず Case1から Case3 の Fujita も解析を行った領域において本解析の結果は(4.3)式とよく一致した。このことから今回の計算は妥当であると考えられる。次に、Fujita が解析を行っていない範囲で(4.3)式との比較を行う。この領域でも本解析は(4.3)式とよい一致を示した。これより、(4.3)式は磁気圏半径が 10km 以下でも適用可能であるとい61

うことが言える。 本解析において得られた推力を図 4-8 に示す。図 4-8 から、磁気セイルが数 N の推力を得るためには磁気インフレーションによって B
磁場を数 10km まで展開する必要があることが分かる。

6241nT の

(a)イオン粒子分布(b)イオン密度分布(c)磁束密度分布
図 4-2

Case1(t=3.0[s])63

(a)イオン粒子分布(b)イオン密度分布(c)磁束密度分布
図 4-3

Case2(t=1.2[s])64

(a)イオン粒子分布(b)イオン密度分布(c)磁束密度分布
図 4-4

Case3(t=0.6[s])65図 4-5

Case4 磁束密度分布(t=1.25×10-2[s])図 4-6

Case5 磁束密度分布(t=2.5×10-4[s])66Drag Coefficient Cd

10
Case1
Case2

1
Case3

0.1
Case4

0.01

FujitaCase5

Our results

0.001
0.1

1

10

100

1000

100

1000

Magnetosphere radius [km]図 4-7(4.3)式との比較

Thrust [N]

10000
100
1
0.01
0.00010.000001
0.1

1

10

Magnetosphere radius [km]
図 4-8

67

推力

第5章

結論

本研究では、MPS の実現のためのキーテーマのひとつである磁気インフレーシ
ョンに着目し、粒子性を考慮したハイブリットコードを用いて解析を行った。
まず、噴射プラズマのラーマ半径を小さくすることが予想できる条件化で磁気インフレーションの解析を行った。ラーマ半径を小さくすることによって微弱な磁気インフレーションしか起きなかった Winglee[9]とは異なる結果を得ることを期待した。その結果、
1) 磁気インフレーションが確認された。
2) 低βの場合でも磁場が密度に比べて減衰勾配が大きいこと、またイオンの速度 v が電場によって加速されることによって磁気インフレーションが起こった。
3) β=0.01 の場合、赤道付近のプラズマは強い磁場の圧力のために噴射後磁場に押し返される。しかし、極付近のプラズマが磁力線に巻きつき赤道に侵入し、そのプラズマによって赤道でも磁気インフレーションが起きた。
次に、MPS の実機を想定して極方向からのプラズマ噴射を行った。その際に噴射プラズマエネルギーを一定にし、速度と密度を変化させることによって効率的な磁気インフレーションを如何に行うかについての議論を行った。その結果、
1) 同じ噴射プラズマエネルギーでも、初期速度が速いプラズマを用いた方が効率よく磁気インフレーションを起こせる。
2) 極付近で噴射されたプラズマの電子はその電子圧力によって+θ方向に移動する。その電子とイオンとの間に生じた電場によってイオンは+θ方向に力を加えられるが、低初期速度の方が長時間この力を受けるため初期速度の遅い方がイオン粒子は+θ方向により多く移動した。
3) 同じ噴射プラズマエネルギーでも初期速度が速い方がより遠くまで磁力線を展開可能であることを確認した。これは、磁場がプラズマに凍結しているためであり磁力線はイオン粒子の動きに対応する。初期速度が速い方がより速度の r 成分を保つためそれに伴って磁力線もより拡大される。
4) 極方向の減衰率は初期速度が小さくなるほど小さくなった。これは 2)と同じ理由からである。
また、Fujita の結果[14]を補うものとして小さい磁気圏半径における磁気セイルの解析を行った。その結果、
1) 小さい磁気圏半径においても Fujita が定式化した(4.3)式は適応可能である。
2) 磁気圏半径が 4km、400m の場合においてそれぞれ 2[mＮ],2.7[μN]の推力68を得られることを確認した。
最後に、MPS のおおよその推力を見積もる。1[T]の超伝導コイルを用いて、磁気インフレーションを行い極方向の磁束密度を 2 乗で減衰させたとすると中心から約 5km の位置で太陽風エネルギーとつりあう。これを磁気圏半径とし、MPS の推力を求めると 8.6[ｍN]となる。この値は他の推進システムと同等程度であるが、第3 章第 3 節のようなβ値を仮定するとプラズマ源に莫大なエネルギーを投入する必要がある。したがって、今後は MPS の優位性を示すためにより低いエネルギーで磁気インフレーションを行い優れた推力電力比を示す必要がある。

69

参考文献
[1] Internet Homepage, http://www.jaxa.jp/[2] R. M. Zubrin, D. G. Andrews, “Magnetic Sails and Interplanetary Travel,”

Journal ofSpacecraft and Rockets, Vol.28, No.2, p.3118, 1991.
[3] R. M. Winglee, J.Slough, T. Ziemba, A. Goodson, “Mini-Magnetospheric Plasma
Propulsion: Tapping the energy of the solar wind for spacecraft propulsion”, Journal ofGeophysical Research, Vol.105, No.21, pp.67-77, 2000.
[4] 山川宏，小川博之，藤田和央，野中聡， 磁気プラズマセイルによる外惑星探査,日本航空宇宙学会論文集, Vol.52, No.603, 2004

[5] 恩藤忠典，丸橋克英「宇宙環境化学」，オーム社，2000.
[6] Internet Homepage, http://www2.nict.go.jp．
[7] 西田浩之，稲谷芳文，小川博之，船木一幸， 磁気セイルの推進特性解析 , ISAS
Research Note, ISAS RN 791, 2005.
[8] 高村秀一 ; プラズマ理工学入門 , 森北出版株式会社，1997[9] R.M. Winglee, J. Slough, T. Ziemba, P. Euripides and D. Gallagher, “Laboratory
Testing of Mini-magnetospheric Plasma Propulsion Prototype”, in Space Technology
and Applications International Forum-2001, edited by EL-Genk, M.S., CP552, AIP
Press, pp. 407-412, 2001.
[10] G. Khazanov, P. Delamere and K. KAbin, “FUNDAMENTALS OF THE PLASMASAIL CONCEPT: MHD AND KINETIC STUDIES” 39th AIAA / ASME / SAE /ASEE Joint Propulsion Conference and Exhibit 20-23 July 2003, Huntsville, Alabama
AIAA 2003 -5225,2003[11] Daisuke Akita, Kojiro Suzuki, “Kinetic Analysis on Plasma Flow of Solar Wind
Around Magnetic Sail” 36th AIAA Plasmadynamics and Lasers Conference 6 – 9 June
2005, Toronto, Ontario Canada

AIAA 2005-4791,2005[12] H. Nishida, H. Ogawa, I. Funaki, Y. Inatani, ”MHD analysis of Magnetosphere around
Magneto Plasma Sail”, ISTS 2006-b-47, 2006[13] R. Asahi, I. Funaki, K. Fujita, H. Nishida, H. Yamakawa, and H. Ogawa, “Performance
Evaluation of a Magneto Plasma Sail”, Proceedings of 8th International70

School/Symposium for Space Simulations,pp.359-360,2005[14] K. Fujita, “Particle Simulation of Moderately-Sized Magnetic Sails”, Journal of SpaceTechnology and Science, Vol.20, No.2, pp.26-31, 2005.
[15] H. Kojima, I. Funaki, Y. Shimizu, H. Yamakawa, Y. Nakayama, and S. Shinohara,“Experimental Simulation of a Plasma flow around Magnetic Sail” , Proceedings ofthe 31st International Electric Propulsion Conference(CD-ROM), IEPC – 2005 – 107 ,
2005.
[16] 篠原大介 ; 九州大学大学院総合理工学府修士論文 , 2005[17] H. Nishida, H. Ogawa, I. Funaki, and Y. Inatani, “Numerical Simulation of Propulsion
Systems Making Use of the Solar Wind “, 44th AIAA Aerospace Sciences Meeting andExhibit, Reno, Nevada,

AIAA-2006-769,2006[18] 永田靖典，大津広敬，小川博之，船木一幸，山川宏，藤田和央 磁気プラズマ
セイルの 3 次元 MHD 解析 ，題 49 回宇宙科学技術連合講演会, 2005[19] I. Funaki, H. Yamakawa, H. Ogawa, H. Nishida, I. Shinohara, K. Fujita, Y. Nakayama
and H. Otsu, ”Feasibility Study of Magnetoplasma Sail.”, 29th International Electric
Propulsion Conference, Princeton University, October 31-November 4, 2005[20] 江崎智規 ; 九州大学大学院総合理工学府修士論文 , 2001[21] 中原大作 ; 九州大学大学院総合理工学府修士論文 , 2003[22] A. P. Matthews, “Current Advance Method and Cyclic Leapfrog for 2D Multispecies
Hybrid Plasma Simulations”, Journal of Computational Physics, Vol. 112, pp.102-116,
1994.
[23] E. J. Horowiz, D. E. Schumaker, D. V. Anderson, J. Comp. Phys.,Vol. 84, p.279,

1989.
[24] 森正武

Fortran77 数値解析プログラミング ,岩波書店 1986[25] C. K. Bridsall, A.B. Langdon, “Plasma Physics via Computer Simulation”, IOP
Publishing Ltd, 1997.
[26] 野田賢治, 梶村 好宏, 中島 秀紀,

磁気プラズマセイルにおける磁場展開に関する研究 、プラズマ・核融合学会九州・沖縄・山口支部第 10 回支部大会, 200671

謝辞
初めに、このような貴重な研究の場を与えてくださった中島秀紀教授に心から感謝します。また、貴重なご意見、助言を頂きました渡辺幸信助教授に心から感謝します。本研究を直接指導して頂き、さまざまな助言を頂いきました梶村好宏講師に心から感謝します。プラズマに関する基本的は質問にも快く答えてくださり様々なアドバイスをしていただいた山本直嗣助手に心から感謝します。いつも研究室の雰
囲気を明るくしてくださった岡安悦子秘書に心から感謝します。
これまで MPS の研究に関して様々なアドバイスを頂いた宇宙航空研究開発機構(JAXA)の船木一幸先生を初めとする MPS 研究グループの方々に心より感謝します。
また、地球シミュレーションセンターの陰山聡先生には私がコードを製作する際に様々なアドバイスを頂きました。心より感謝します。
研究について方向性などさまざまなアドバイスを頂いた D3 の桑野寛久氏に心から感謝します。D2 の林真照氏にはパソコンのトラブルの際に助けて頂きました・
心から感謝します。M1 の時から互いに切磋琢磨し、公私共に助け合ってきた M2の小川真司氏・加来大輔氏・笠岳幸氏・近岡貴行氏・宗貞信孝氏に深く感謝いたします。また、M1 の金川隆保氏・川渕亮氏・近藤慎哉氏・永田英隆氏・西嶋康太氏・
三好誠氏、B4 の大宅将史氏・松田伸夫氏・本岡親英氏には公私両面で大変お世話になりました。心より感謝します。
B4 で当研究室に配属され、ここまで様々な指導をしていただいた当研究室の卒
業生の方々に心から感謝します。
最後に、私を支え暖かく見守ってくださった家族に心から感謝します。

72