卒

題

業

論

文

目

磁場中レーザー生成プラズマの挙動解析

氏

名

富永

宙志

所属教育分野プラズマ理工学

指 導 教 員

中島
山本

秀紀
直嗣

教授
准教授

九州大学工学部エネルギー科学科
提出年月

平成２４年２月

目次
第1章

序論 ..................................................................................................................................... 1

1.1

背景 ......................................................................................................................................... 1

1.2

研究の目的 .............................................................................................................................. 2第2章

磁気スラストチャンバー .................................................................................................... 3

2.1

システムの背景....................................................................................................................... 32.2

推進原理.................................................................................................................................. 3

2.3

磁気スラストチャンバーの利点と欠点 .................................................................................. 5

第3章

計算コード .......................................................................................................................... 6

3.1

１次元放射流体コード ............................................................................................................ 6

3.2

１次元のコードを使う目的 .................................................................................................... 8

3.3

STAR1D による実際の結果 ................................................................................................... 9

3.4

カップリングコード ............................................................................................................. 12

3.5

カップリングコードの妥当性について ................................................................................ 14

3.6

3D ハイブリッドコード........................................................................................................ 15

3.6.1

計算方法 ........................................................................................................................ 17

3.6.2

プラズマ粒子の運動 ...................................................................................................... 173.6.3

イオン数密度、電流密度 ............................................................................................... 19

3.6.4

電磁場の計算方法 .......................................................................................................... 20

3.6.5

Digital Filter .................................................................................................................... 22

3.6.6

ハイブリッドコードに課せられる制限 ......................................................................... 23第4章

計算体系と結果 ................................................................................................................. 24

計算体系 1-1 ............................................................................................................................. 25
計算体系 1-2 ............................................................................................................................. 27
計算体系 1-3 ............................................................................................................................. 29
計算体系 2-1 ............................................................................................................................. 33
計算体系 2-2 ............................................................................................................................. 35
計算体系 2-3 ............................................................................................................................. 37
計算体系 3-1 ............................................................................................................................. 40
計算体系 4 ................................................................................................................................ 43
第5章

まとめ、今後の課題.......................................................................................................... 49

参考文献 ................................................................................................................................... 50
付録 1.......................................................................................................................................... A

第1章
1.1

序論
背景

現在多くの国が乗り出している宇宙探査。アメリカやロシアを筆頭に、ヨーロッパやアジアの国々でもその研究が進められている。そんな中注目されているのが磁気ノズルというプラズマの制御技術である。磁気ノズルを用いると、爆発的に膨張するプラズマを磁場との相互作用を利用し、固体壁との相互作用無く運動方向を変えて推力を得ることが出来る。従って、他の宇宙推進システムと比較して、高い排出速度（即ち高い推進材利用効率）と大きな推力を同時に達成可能である。この磁気ノズルの概念図を Fig1-1 に示す。
磁気ノズルは過去に設計、シミュレーションが行われてきた。Hyde は、超伝導(SCM)コイルを１つ用いたモデルを提案している[1]。また、２次元 MHD シミュレー
ションコードを用いてプラズマ挙動と推進効率を計算し、1 つのコイルで推進効率65％という結果を得ている[2]。長峯と中島[3]は３次元(3D) hybrid code を用いてプラズマ挙動と推進効率を計算し、推進効率 65％という結果を得ている。Vchivkov[4]は同様のコードで推進効率の最適化を行い、推進効率 70％という結果を得ている。また、実験との比較も行い、推進効率の一致を得ている[5]。

Figure 1-1 磁気ノズルの概念図

1

1.2

研究の目的

本研究室ではレーザー核融合ロケットについて、レーザー生成プラズマを用いて、大阪大学で実験を行っている。その実験と同じ、または近い体系でシミュレーショ
ンを行い、どの条件下で推力が向上するかを探り、またその結果を実験にフィードバックしていくことを目的とした。そこで、本研究ではレーザープラズマの生成過程については、大阪大学より提供を受けた１次元放射流体コード（STAR1D）を用いて計算した。プラズマの磁場中の挙動については、3D hybrid code にて計算を行った。
このように過程を分けて２つのコードをカップリングさせてシミュレーションすることで、より磁気ノズルにおける推力評価の精度向上を試みた。

この論文の構成は以下の通りである。
第 1 章 本研究の背景と目的を述べた。
第 2 章 磁気スラストチャンバーシステムについて述べる。
第 3 章 本研究に用いた計算コードについて述べる。
第 4 章 解析結果について述べる。
第 5 章 結論、まとめを示す。

2

第2章
2.1

磁気スラストチャンバー
システムの背景

現在、NASA などの宇宙機関において火星への有人探査の計画が進められている。
火星は、人類が地球圏外へ進出する際、比較的早い段階に移住できると考えられる唯一の天体である。火星の大気には二酸化炭素が豊富に存在し、大量の水が凍った状態で眠っているからである。また、火星は人類にとって次の居住地となるばかり
ではなく、火星での経験が、他の惑星の探査・開発にとって重要なものとなり、前哨基地としての役割も担うことになるだろう。しかし、近い将来、火星に恒久的な基地や居住施設を建設することを考えると、従来の宇宙輸送手段である化学推進のロケットでは建設に必要な膨大な量の資材を運搬することは容易ではない[6]。
そこで、その際に使用される宇宙輸送のシステムとして考えられているのが、磁気スラストチャンバーシステムを利用するレーザー核融合ロケットである。ロケットの推進システムにエネルギー密度が大きい核融合を用いると、化学推進に比べて高い比推力と高い比出力を同時に得ることができる。推進材の排出速度は107 [𝑚⁄𝑠𝑒𝑐 ]程度になると見られ、比推力（排出速度を重力加速度で割った値）は106 [𝑠𝑒𝑐]となる。
この値は、化学推進ロケットの 2000 倍程度であり、ペイロード（積載貨物）の割合が多く取れることになる[7]。さらに、ロケットの外部空間をエンジンとして利用するので、その分宇宙船に占めるエンジンの質量が小さいので高比出力が達成されると考えられる。これは、輸送時間の短縮にもつながり、宇宙線による被爆など人体への悪影響も抑えられる。
レーザー生成プラズマと磁場の相互作用を利用するプラズマ推進器の研究は 1970
年代に始まり，その目的はレーザー核融合による惑星間航行を目指した高推力な推進システムであった。この推進器の特徴は，高速航行によるミッション期間の短縮と，ミッションに応じた設計柔軟性，そして，長寿命性である．さらに、ロケットの外部空間をエンジンとして利用するので、その分宇宙船に占めるエンジンの質量が小さいので高比出力が達成されると考えられる。本研究ではこのようなレーザー
生成プラズマと磁場の相互作用を利用するプラズマ推進器を磁気スラストチャンバ
ーとよぶ。

2.2

推進原理

磁気ノズルでの推進原理を Fig.2-1 に示す。プラズマは反磁性体であり、磁場によってその運動方向が変化する。コイルによって作られた磁場中でターゲットにレーザーを照射し、プラズマを発生させたとき、このプラズマは大きな運動エネルギー3を持っている[a]。コイル方向に膨張するプラズマは Larmor 運動し、反磁性電流が流れるように動きながら磁場を圧縮する。磁場は圧縮されることでより強度が増し、プラズマに対して磁場に直角に強い力を及ぼすようになる。これによってプラズマは減速され、プラズマから磁場へエネルギーが輸送される[b]。さらに、磁場によって減速されたプラズマはやがて停止し、磁場のエネルギーの方がプラズマのエネル
ギーよりも強くなると磁場は押し戻す。それに伴って、プラズマは磁場により後方へ加速される[c]。このときに、磁場からプラズマへエネルギーが与えられる。このようにして、磁場がプラズマ粒子の運動方向を変化させ、押し出す反作用により推力を得ることができる[8]。

[a]

[b]

[c]

Figure 2-1 磁気ノズルによる推進の原理

4

2.3

磁気スラストチャンバーの利点と欠点

大きな推力を持つ推進器でなければ、加速力が不足し、必要な速度増分を得るには長い加速時間を有することになる。これはミッションの長期化を意味し、有人惑
星間飛行にとっては致命的な欠点となる [9]。大きな推力を得られる磁気スラストチ
ャンバーによるミッション期間の短縮は，有人飛行において大きな障害となる宇宙線による被爆期間を短くすることができ，将来の有人飛行にとって有益である。また、この推進器は惑星間航行だけではなく地球近傍におけるミッションにおいても活躍することができると期待されている。レーザー生成プラズマはレーザーパワーと燃料推進剤の質量は独立関係なので、プラズマの排出質量と排出速度を広範囲に制御することができ、ミッションに応じた最適な設計が行える[10]。そして、プラズ
マは磁場によって排出方向を制御されるので、ノズルなどの固体壁に接触損傷を与えず長寿命な推進器となる。つまり、この推進器は地球近傍から惑星間航行までのあらゆるミッションに対応できる万能な推進器である。
しかし、レーザー核融合ロケットの実現には多くの課題がある。レーザー核融合を実現するためのレーザードライバ、大型超伝導コイル、大量の核融合燃料ペレット生成などの開発が急務であるが、これらは現在活発に研究が行われている地上発
電用核融合技術がそのまま流用できる。さらに、2003 年に Orth が提案したレーザ
ー核融合ロケットの概念設計[11]では、宇宙船全体の質量が 6500ton もあり、2005 年に本研究グループが行った新しい概念設計でも 250ton の巨大な宇宙船になるという計算結果が導き出された。この巨大な構造物を宇宙で組み立て、さらに高速での運転を行うための技術も必要である。

5

第3章

計算コード

この章では、次章以降の実験解析に使用する、計算コードについて説明する。最初にプラズマのレーザーによる生成過程について計算する１次元放射流体コード、次にその後に使用する 3D ハイブリッドコードに、データを受け渡すカップリングコードについて説明する。最後に、磁場中に置いたプラズマの挙動について計算する 3D ハ
イブリッドコードについて説明する。

3.1

１次元放射流体コード

大阪大学より提供を受けたコードで、名前を STAR1D[12]と呼ばれている。特徴としては、高いエネルギー密度を持ったレーザー生成プラズマについて計算することができる。また、平面や球面幾何形状における１流体２温度の挙動について扱う。１次元、２次元コード共に一流体二温度モデルを採用しており、１つの運動の式と、イオン、
電子それぞれのエネルギー式を別々に解く。（1）式から（4）式に１次元の流体の基
礎式を示す。
基礎方程式
𝑑𝜌
𝜕𝜈
= −𝜌
𝑑𝑡
𝜕𝑥
・・・・・・(1)𝑑𝜈
1 𝜕(𝑝𝑖 + 𝑝𝑒 + 𝑞)
=−
𝑑𝑡
𝜌
𝜕𝑥
・・・・・・(2)
𝜌𝑐𝑉𝑖

𝑑𝑇𝑖
𝑑𝜈
𝜕
𝜕𝑇𝑖
= −𝑃𝑇𝐻𝑖+(𝜅𝑖) + 𝛼(𝑇𝑒 − 𝑇𝑖 ) + 𝑄𝑖
𝑑𝑡
𝑑𝑥 𝜕𝑥
𝜕𝑥
・・・・・・(3)

𝜌𝑐𝑉𝑒

𝑑𝑇𝑒
𝑑𝜈
𝜕
𝜕𝑇𝑒
= −𝑃𝑇𝐻𝑒+(𝜅𝑒) + 𝛼(𝑇𝑒 − 𝑇𝑖 ) + 𝑄𝑒
𝑑𝑡
𝑑𝑥 𝜕𝑥
𝜕𝑥

・・・・・・(4)
ここで(1)式は連続の式、(2)式は運動の式、(3)式はイオンのエネルギー式、(4)式は電子のエネルギー式であり 𝑝𝑇𝐻𝑖 = 𝑇𝑖 (𝜕𝑝𝑖 ⁄𝜕𝑇𝑖 ) 、𝑝𝑇𝐻𝑒 = 𝑇𝑒 (𝜕𝑝𝑒 ⁄𝜕𝑇𝑒 ) である。 は流体の速度、𝜌は密度、𝑝𝑖・𝑝𝑒 はそれぞれイオン、電子の圧力、𝑞は人工粘性である。
𝑐𝑉𝑖 ・𝑐𝑉𝑒 はそれぞれイオン、電子に対する比熱、𝑇𝑖 ・𝑇𝑒 はそれぞれイオン温度、電子温度である。イオンと電子のエネルギー式(3)、式(4)を解く際にはイオン、電子の熱伝導、およびイオンと電子間のエネルギー緩和も同時に解く。(3)式および(4)式の右
辺第２項はそれぞれ熱伝導項を表しており，熱伝導係数にはSpitzer の熱伝導係数[13]を採用し、熱束制限拡散近似モデル[14]により熱流計算を行う。ここで電子熱流の熱束制限係数は0.1 とした。また(3)式と(4)式の右辺第３項は電子とイオンのエネルギ
ー緩和を表しており，Spitzerによる熱緩和係数𝛼を用いる[15]。右辺第４項の𝑄𝑖・𝑄𝑒 は、
6それぞれイオン、電子のエネルギー式における生成項であり、人工粘性に伴う運動エネルギーからイオンの熱エネルギーへの変換分は𝑄𝑖 に代入される。また、後に述べるレーザー吸収とプラズマ中の輻射吸収により電子が受け取るエネルギーは𝑄𝑒 に代入される。次にレーザーエネルギー吸収について記す。レーザー光の吸収は逆制
動放射[16]を仮定し、計算領域端からターゲットに垂直に入射されるレーザー光線に沿って、局所的な電子密度 𝑒 、平均電離度 、電子温度𝑇𝑒 から求められる吸収係数を用いて各位置におけるレーザー光の吸収を計算した。
次に輻射輸送について記す。輻射輸送は熱束制限多群拡散近似によりX 線スペク
トルを多数のエネルギー群に分け、それぞれのエネルギーについて拡散近似により解く[17]。(5)式に輸送方程式を示す。
𝑑 𝐸𝜈
𝜕
𝜕𝐸𝜈
𝜌 ( )−(𝐷𝜈
) = 4𝜋𝜂𝜈 − 𝑐𝜒𝜈 𝐸𝜈
𝑑𝑡 𝜌
𝜕𝑥
𝜕𝑥
・・・・・・(5)輻射輸送方程式
ここで𝐸𝜈 は、振動数𝜈の光子に対する輻射エネルギー密度である。𝐷𝜈 は拡散係数を表す。拡散近似により𝐷𝜈 = 𝑐⁄( 𝜒𝜈 )と定義されるが、最大エネルギー束 𝜈 = 𝑐𝐸𝜈 を
超えないように熱束制限がかけられる．𝑐は光速である。𝜂𝜈 と𝜒𝜈 はそれぞれ放射率𝜂𝜈 、
吸収率𝜒𝜈 を表す。𝜒𝜈 = 𝜌𝜅𝜈 であり、𝜅𝜈 は振動数𝜈の光子に対するopacityである。
（５）
式の右辺は流体からの輻射によるエネルギー損失を表しており、スペクトル積分した値は、電子のエネルギー式の生成項𝑄𝑒 に代入される。光子のエネルギー群は１次元計算では0-1.5 keV を1500グループに不当間隔で分割し、EUV 領域では0.2 eV 刻
みで計算を行う。２次元計算では計算負荷が大きいため全エネルギー領域は40グループに分割され、EUV 領域では0.4 eV 刻みである。上記（１）式から（５）式において、𝑝𝑖 ・𝑝𝑒 ・𝑝𝑇𝐻𝑖 ・𝑝𝑇𝐻𝑒 ・𝑐𝑉𝑖 ・𝑐𝑉𝑒 ・𝜂𝜈 と𝜒𝜈 は密度と温度の関数として表され、状態方程式コードおよび原子過程コードで計算されたテーブルを引用することにより、実際の物性に即したシミュレーションを可能としている。
𝑞 ∶ 人工粘性𝜕𝑝𝑖
𝜕𝑝𝑒
𝑃𝑇𝐻𝑖 = 𝑇𝑖 ( ) , 𝑃𝑇𝐻𝑒 = 𝑇𝑒 (
)
𝜕𝑇𝑖
𝜕𝑇𝑒
𝑐𝑉𝑖 , 𝑐𝑉𝑖 ∶ 比熱
𝛼 ∶ 熱緩和係数
𝑄𝑖 , 𝑄𝑒 ∶ 生成項𝐸𝜈 ∶ 振動数𝜈の光子に対する輻射エネルギー密度𝜅𝜈 ∶ 振動数𝜈の光子に対する不透明度（𝑜𝑝𝑎𝑐𝑖𝑡𝑦）𝜂𝜈 ∶ 放射率𝜒𝜈 ∶ 吸収率 = 𝜌𝜅𝜈𝐷𝜈 ∶ 拡散係数 = 𝑐⁄( 𝜒𝜈 )7

イオン、電子のエネルギー保存式にそれぞれのエネルギーや熱の収支を考慮した生成項𝑄を加えている。イオンの位置、イオンの速度、イオン価数、イオン質量数を出し、これを 3D ハイブリッドコードにデータとして受け渡す。

エネルギーの投入について
ターゲットが平板の場合、投入したエネルギー値と出力した速度や密度、メッシ
ュの位置などから得たエネルギーの値が異なる。このことについては、可能性として付録 1 に載せている。

3.2

１次元のコードを使う目的

放射流体コードには、１次元（STAR1D）と２次元（STAR2D）とがある。もちろん、計算の精度や、ターゲットなどの条件、体系の幅などは２次元の方がより良いことは言うまでもない。では、なぜ１次元を使用するのか。その理由は１次元において計算できる体系においては２次元で計算するよりも計算時間が短いところにある。計算精度を高めて行う２次元放射流体コードの欠点は、何よりもその膨大なデータを結果として出さなければならないところである。データ量が多くなれば自然と計算時間も長くなる。確かに２次元にしかできない体系もあるが、平板のみ、球のみなど１次元でもできる体系に絞って考えると、やはり計算時間が短いことは強みである。
データの比較をしたり、実験にフィードバックしたりする上で、まずは計算時間の短い１次元で、より効果の高い条件の目途を立てることができ、その後へつなぐスパンが短くなる。以上が、１次元放射流体コードを使用する目的である。

8

3.3

STAR1D による実際の結果

ここで、STAR1D で実際に計算を行い、得た結果をここに示す。計算体系は Fig.3-1の通りである。

TargetLaser
Thickness of a target
7.0 × 10−5 𝑐𝑚

Figure 3-1 計算体系

Thickness of a target：0.7𝜇𝑚
Target：plastic
Laser wavelength：1.06𝜇𝑚
Focal position：−0.520𝑐𝑚
Intensity：0. 92 × 1012 𝑊/𝑐𝑚2
F number：20
Target mesh number：200
4.6𝐽
Table 1 各種パラメータ9𝑛𝑠
94.00E+00
レーザーピーク時（10ns）

3.00E+00
2.00E+00
レーザー方向1.00E+00
0.00E+00-1.00E+00

2.50E-09
5.00E-09
7.50E-09
1.00E-08
1.25E-08
1.50E-08
1.75E-08
2.00E-08
2.25E-08
2.50E-08
2.75E-08
3.00E-08
3.25E-08
3.50E-08
3.75E-08
4.00E-08
4.25E-08
4.50E-08
4.75E-08
5.00E-08

position of mesh (r)半値幅（4.5𝑛𝑠）-2.00E+00

time (t)Figure 3-2 STAR1D における r-t 線図
Fig.3-2 は、メッシュの位置の時間変化を表したグラフである。つまりプラズマの広がりを時間ごとに表している。
図を見て分かるように、レーザーが入射してきた方向の方が時間的に早く広がっており、背面の部分において、最初はあまり広がらず、レーザーピークのあたりから一気に広がっている。これは、背面への衝撃波によって時間差でプラズマ化していると考えられる。
次のページの Fig.3-3 は、それぞれの時間における圧力の最大値がどの位置にあるのかを表している。このグラフを見ると、一番大きな圧力は、レーザーピークの時（10ns）である。また 5ns まではレーザーの入射方向に最大圧力点があるが、レーザーピーク
時に近づくにつれて背面方向に移動している。前述の通り、メッシュの広がりがレー
ザーピークのあたりから背面方向にも急激に広がっていることと合わせて、背面に衝撃波が伝わるその分、プラズマの広がり方に時間のズレが生じていると考えられる。

10

9.00E+10

レーザーピ－ク

time
[sec]

highest pressure[Ba(=× 10−1 Pa)]

8.00E+10
7.00E+10

レーザー入射方向

6.00E+10
5.00E+10
4.00E+10
3.00E+10
2.00E+10
1.00E+10
0.00E+00

-0.1

-0.05

0

2.50E-09
5.00E-09
7.50E-09
1.00E-08
1.25E-08
1.50E-08
1.75E-08
2.00E-08

0.05

position (r) [cm]Figure 3-3 時間ごとの最大圧力の位置
1E+24
1E+23

初期におけるターゲットは、
0 から − 0.7μmまで広がっている

1E+22

2.5ns

density g⁄cm−3

1E+21
14.5ns

1E+20
1E+19
1E+18
10ns

1E+17

-2.00E-01-1.00E-01

0.00E+00

1.00E-01

2.00E-01

3.00E-01

4.00E-01

5.00E-01

6.00E-01

1E+16

position of mesh (r) [cm]Figure 3-4 メッシュの広がりと密度
Fig.3-4 は、計算が始まってすぐの時（2.5ns）、レーザーピークの時（10ns）と、後ろ側の半値幅を越えた時（14.5ns）の粒子密度の比較である。図を見れば明らかなように、
レーザーピーク時までには、まだレーザーの入射方向にしかあまり広がっていないが、半値幅を越えると、粒子が左右の両方に広がっているのが分かる。
11

カップリングコード

3.4

以下にはカップリングコードの説明をする。STAR1D の output ファイルから必要なデータを読み取って、それを 3D ハイブリッドコードに入力できるようにする。メ
ッシュごとに、プラズマ半径（プラズマが広がっている位置:leng）、価数、密度、圧力、温度のデータを読み取る。
さらにターゲットの形が球か薄膜平板かを選択し、照射半径（𝑟𝑠𝑝𝑜𝑡 = 250𝜇𝑚）、
カップリングさせるときの角度（𝛼）、メッシュ数を入力する。

粒子の位置に関して
x,y 方向
STAR1D においては、１次元であるため z 方向についてのみ計算している。そのため、3D ハイブリッドに渡すデータで、i 番目の粒子の x,y 方向の位置はランダムに、また円領域内に一様に振り分けて x(i)、y(i)を決める。
z 方向
ランダムな数 xx(0 から 1 の間)を決める。（i は粒子の番号、j はメッシュの番号）
𝑙𝑒 𝑔(𝑗 − 1) + {𝑙𝑒 𝑔(𝑗) − 𝑙𝑒 𝑔(𝑗 − 1)} × 𝑥𝑥 = 𝑧(𝑖)上記で決まった x(i)、y(i)、z(i)に関して zx 平面において最初に入力した角度𝛼だけ傾けるために、𝑧𝛼 (𝑖) = 𝑧(𝑖)𝑐𝑜𝑠(𝛼) − 𝑥(𝑖)𝑠𝑖 (𝛼)𝑥𝛼 (𝑖) = 𝑧(𝑖)𝑠𝑖 (𝛼) + 𝑥(𝑖)𝑐𝑜𝑠(𝛼)𝑦𝛼 (𝑖) = 𝑦(𝑖)このようにして座標変換する。

速度についてx、y 方向10−30 𝑚/𝑠𝑒𝑐だけ外側へ速度を持つ。
𝑥(𝑖) = 1.0 × 10−30 [𝑚/𝑠𝑒𝑐]𝑦(𝑖) = 1.0 × 10−30 [𝑚/𝑠𝑒𝑐]z 方向
最初に STAR1D から読み取ったデータを線形補完する。
（𝑖 は粒子の番号、𝑗 はメッシュの番号）
𝑧(𝑖) = 𝑧 (𝑗 − 1) + {𝑧(𝑖) − 𝑙𝑒 𝑔(𝑗 − 1)} ÷ {𝑙𝑒 𝑔(𝑗) − 𝑙𝑒 𝑔(𝑗 − 1)}× { 𝑧 (𝑗) − 𝑧 (𝑗 − 1)}
速度についても位置とおなじように座標変換する。
𝑧𝛼 (𝑖) = 𝑧(𝑖)𝑐𝑜𝑠(𝛼) − 𝑥(𝑖)𝑠𝑖 (𝛼)12𝑥𝛼 (𝑖) = 𝑧(𝑖)𝑠𝑖 (𝛼) + 𝑥(𝑖)𝑐𝑜𝑠(𝛼)𝑦𝛼 (𝑖) = 𝑦(𝑖)Figure 3-5 vz(i)の値の決め方

プラズマ全体の質量について
STAR1D の出力から次のように計算される
𝑚(𝑗) = 𝜌(𝑗) × 𝜋 × 𝑟𝑠𝑝𝑜𝑡 2 × (𝑙𝑒 𝑔(𝑗) − 𝑙𝑒 𝑔(𝑗 − 1))
𝑚𝑡𝑜𝑡𝑎𝑙 = 𝑚(1) + 𝑚(2) + ・・・ + 𝑚(𝑖) + ・・・ + 𝑚(メッシュ数)
3D ハイブリッドコードでは、超粒子を使用している。超粒子 1 個の重さは𝑚𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 = 𝑚𝑡𝑜𝑡𝑎𝑙 ⁄ 𝑝（ 𝑝は、超粒子の個数）

エネルギーについて𝑛𝑝

1
𝐸 = × 𝑚𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 × ∑( 𝑥𝛼 (𝑖)2 + 𝑦𝛼 (𝑖)2 + 𝑧𝛼 (𝑖)2 )
2
𝑖

13

𝑥
𝑙𝑎𝑠𝑒𝑟
𝛼

𝑧

𝑡𝑎𝑟𝑔𝑒𝑡𝑐𝑜𝑖𝑙

Figure 3-6 カップリングコードにおける概念図

3.5

カップリングコードの妥当性について

実際に計算した値から評価する。
プラズマ質量・カップリングコードにおいて
𝑚𝑡𝑜𝑡𝑎𝑙 = 1.45651 587 44 42 × 10−10 [𝑘𝑔]・厚さ0.7𝜇𝑚の薄膜ターゲット、初期状態において
𝑚𝑡𝑜𝑡𝑎𝑙 = 1.06 × 103 × 𝜋 × (2.5 × 10−4 )2 × 0.7 × 10−6
= 1.45651 587 44 42 × 10−10 [𝑘𝑔]但し、ターゲット（プラスチック）の密度 : 1.06 𝑔⁄𝑐𝑚3 、照射径 : 250𝜇𝑚である。
初期のプラズマ質量と、カップリングコードで計算した結果との誤差は、約2.81 × 10−7 %で、ほぼ無視できる誤差となっている。
エネルギー
STAR1D で、カップリングのタイミング時（50ns）における出力から、粒子の運動エネルギーを出す。

14

メッシュ数𝐸𝑆𝑇𝐴𝑅1𝐷 =

∑
𝑗

×(

2
1
{ × 𝜌(𝑗) × 𝜋 × 𝑟𝑠𝑝𝑜𝑡(照射径) × (𝑙𝑒 𝑔(𝑗) − 𝑙𝑒 𝑔(𝑗 − 1))2𝑧 (𝑗) + 𝑧 (𝑗 − 1)22) } = 4.6 90045 [𝐽]

次に前記のエネルギーの式より、3D ハイブリッドコードに入力する速度から、カップリングコードにおけるエネルギーを計算する。
（式は 3.3；エネルギーについて）𝐸 = 4.640246702 6 6 7 [𝐽]STAR1D の出力とカップリングコードにおけるエネルギーの計算結果の誤差は、
約2.68 × 10−2 %で、ほぼ無視できる。
妥当性の評価としては、質量についても、エネルギーについても問題なく、カップリングコードとしては妥当であると言える。

3.6

3D ハイブリッドコード

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

計算に使用する方程式を導出するための基礎方程式は次のように与えられる。以下の方程式は直交座標系および、MKSA 単位系で記述されている。まず、粒子として取り扱うイオンの運動方程式を示す。
midv i
 ZeE  v i  B 
dt
dx i
 vi
dt

(3.1)
(3.2)

15

ここで v i はイオンの速度、Z はイオンの価数、e は電気素量、x i はイオンの位置、miはイオンの質量、また、 E 、 B はそれぞれ電場、磁場を表している。
流体として扱う電子の運動方程式を電場の解法に用いる。電子の運動方程式は以下の式で表される。
ne medv e ene E  v e  B   Pe
dt

(3.3)

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

Pe  neTe(3.4)

ここで Te は電子温度である。(3.4)式を(3.3)式に代入して Pe を消去する。

0  ene E  v e  B  neTe (3.5)

荷電準中性条件および、ダーウィン近似を用いたアンペアの法則により以下の式が得られる。

  B p   0 J e  J i (2.6)

ここで  0 は真空中の誘電率であり、 J i 、 J e はそれぞれイオン電流密度、電子電流密度である。また、 B p はプラズマにより生成される磁場で以下のようになる。

B p  B  B ext  B pick

(3.7)

ここで、B は計算領域内部の合成磁場、B ext は外部磁場、および B pick はピックアップ
コイルにより生成される磁場である。ここで、ピックアップコイル磁場は、次式で示される。

  B pick   0 J pick

(3.8)

但し、 J pick はピックアップコイルに流れる電流密度である。

J i はイオンの数密度と電流密度から計算することが出来る。 J e は電子電流密度の定16

義により以下のようになる。

J e  ene v e(3.9)(3.6)、(3.9)式を用いて、 J e を消去し v e について求めると以下のようになる。

ve  
1  1
   B p  J i ene   0
(3.10)

これを(3.5)式に代入して電場 E について整理する。荷電準中性条件より ene  Zeni を考慮し次式を得る。

E


1 1
   B p  B  J i  B  ZTe ni 
Zeni   0


(3.11)

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

3.6.1(3.12)

計算方法(3.1)式から(3.12)式を数値計算によって解くためには、基礎式を差分式に直し、数値的に安定で、物理的に妥当な空間幅、時間幅で解かなければならない。
計算の大まかな手順を以下に示す。
[１]外部磁場および、イオンの粒子位置、速度の初期設定を行う。
[２]格子点上でのイオンの数密度、電流密度を計算する。
[３](3.11)式により、電場を求める。
[４][１]で設定した外部磁場と、[３]で求めたプラズマによる電場と磁場から、(3.1)式、(3.2)式を用いて新たなイオン粒子の位置、速度を求める。
[５](3.12)式により磁場を求める。
[２]～[５]を繰り返し、プラズマの時間的変化を追う。
３次元ハイブリッドコードに対する各方程式の差分式を以下の小節で示す。

3.6.2

プラズマ粒子の運動

本コードでは、いくつかのイオンを超粒子と呼ばれる集団として扱い、この超粒子17

10 万個を用いた。つまり、超粒子を個別に運動させることでプラズマ粒子の運動を追う。(3.2)式に、leap-frog 法[17]による時間中心差分を施すと次式となる。
x in 1  x in
 v in 1 / 2
t(3.13)

ここで各変数の上付き添え字は時間を表し、下付き添え字は i 番目のイオンを表す。
また、 t は時間幅を表す。
(3.1)式に時間中心差分を施すと以下の式となる。
v in1 / 2  v in1 / 2 Ze  n 1 n1 / 2
(3.14)
 v in1 / 2  B n 
E  v i
t
mi 
2






n  1 / 2 はある時刻 n から t / 2 だけ以前の値であることを示す。同様に、 n  1 / 2 は
t / 2 だけ進んだ値であることを示す。上付き添え字が n 、 n  1 / 2 である変数は既知量である。
プラズマ粒子の運動を求めるためには、まず(3.14)式において、v in 1 / 2 を左辺にまとめて、未知量である v in 1 / 2 を解く。次に、求めた v in 1 / 2 を(3.13)式の右辺に代入し、新たな粒子の位置 x in 1 を求めることが可能となる。
(3.14)式を変形すると次式となる。
Av n1 / 2  S
 、  、  、  を次のように定義する。
qt n

Bx
2 mi



qt n
By
2 mi

 

qt n
Bz
2 mi

(3.15)

  1 2   2   2
なお、ここでの E 、B は粒子位置におけるそれぞれの値を示している。これらを用いて、行列 A とベクトル S は以下のように構成される。
 1

A 
 



1



 

 

1 

 vixm  cE xn  v ym  v zm S   viym  cE yn  v zm  v xm  v m  cE n  v m  v m 
z
x
y 
 iz
18

(3.16)(3.17)

ここで、 c  qt mi 、 m  n  1 / 2 である。さらに、(2.15)式を次のように変形して、速度を求める式が得られる。
v n1 / 2  A 1S(3.18)

ここで、 A 1 は以下のように得られる。
 1 21A 1       3.6.3

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

(3.20)

イオン数密度、電流密度

イオン数密度を求めるためには、空間内に広がっているプラズマを格子点上に割り付
けなければならない。そのために、PIC（Particle In Cell）法を用いる。今、１つの超粒子の中に N 個のイオンが含まれている場合を考える。Fig.3-7 に示すように、格子点内に存在する超粒子の位置によって、8 個の格子点によって囲まれた体積 V の立方体を、

V1 、 V2 、～ V8 の体積を持つ直方体に分割する。格子点 i, j, k  に割り振られる超粒子に含まれるイオンの数は次式で与えられる。

ni (i, j, k ) V2
NV(3.21)

すべての粒子についての和を ni についてとれば、ある格子点上の ni を定義することが出来る。
イオン電流密度 J i は、数密度 ni を用いて次式により求めることが出来る。

J i  Zeni v i(3.22)

19

Figure 3-7

3.6.4

電磁場の計算方法(1) プラズマ領域での計算
電磁場を求めるためには、(3.11)式、(3.12)式を用いる。(3.12)式を片側差分して、B n1 / 2
だけを左辺に残すと次式となる。
B n 1 / 2  B n t  En2(3.23)

また、(2.11)式は以下の式になる。

En 


T
1 1
1 n
  B np  B n 
J i  B n  e nin ni   0 Ze
Ze
e






(3.24)

ここで、B n 、E n は既知であるとすると、 (3.23)式から B n1 / 2 を求めることが出来る。
次に、 B n1 / 2 、 J n1 / 2 、 v n1 / 2 、 nin 1 / 2 を(3.24)式に代入して E n1 / 2 を求め、線形外挿より電11
磁場の予測値 E npred、 B npredを次式より求める。

1
E npred E n  2E n1 / 2(3.25)
20

1
B npred B n 1 / 2 t
1  E npred2(3.26)

1
1
ここで、下付き添え字 pred は予測値であることを表している。 E npred、B npredを求めた3 / 2
3/ 2
3/ 2
3/ 2
後、それらによるイオンの運動を解き、 v in pred、nin pred、J inpredを求める。次に、B npredを

次式により求める。
3 / 2
1
B npred B npredt  E n1 / 2
2(3.27)

3 / 2
そして、(2.24)式で E n を求めた時と同様にして、 E npredを求める。最後に次式に示す

線形内挿から次のタイムステップの値、 E n1 、 B n 1 を求める。
E n1 1 n 1 / 2 1 n3 / 2E E pred22(3.28)t  E n 1
2

(3.29)

B n 1  B n 1 / 2 

これらの値 E n1 、 B n 1 より、 v in 3 / 2 、 nin 3 / 2 、 J in3 / 2 を求めることが出来る。
ここで、電場を求める際に用いた一連の手法は、予測子－修正子法と呼ばれる収束計算の一種である。この方法では、n タイムステップから n  1 タイムステップの値を求めるために、n  3 / 2 タイムステップの予測値を求め内挿計算を行う。そして、求めた n  1
タイムステップの値を再度予測値として用いる。具体的には、(3.29)式により求められ1
た B n 1 を次の予測値 B npredとし、(3.25)式から(3.29)式までを次式で定義する条件[18]を満足するまで、繰り返し計算を行う[19]。

max i , j ,k Bin, j ,1k( q 1)  Bin, j ,1k( q )max i , j ,k Bin, j ,1k( q )(3.30)

ここで q は収束計算の回数を表し、 i 、 j 、 k はそれぞれ x、y および z 方向の格子番号を表している。  は収束計算を打ち切るための値である。
電場を計算するために用いる(3.24)式は、分母に数密度を含んでいるため、 ni がゼロとなる領域ではこの式は発散してしまう。さらに、数値計算の特性上 ni があまり小さい
値になると、発散もしくは計算誤差が生じる。このことを避けるために、数密度がある
特定の値 nv 以上の場合においてのみこのような電場の計算を行う。

21

(2)

真空領域での電場の計算

真空領域では、ラプラス方程式  2 E  0 を SOR（Successive Over Relaxation）法を用いて解く。SOR における加速係数  は正方形領域の場合は次のように表される。
2(3.31)
1  sin  NGMX 
ここで NGMX は x、y または z 方向の格子数である。  の値は通常 1.8 から 1.9 付近である。また、   1 の場合には、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 ) E p (i, j , k  1)  E p (i, j , k  1) 6(3.32)

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

 p  E p (i  1, j , k )  E p (i  1, j, k ) E p (i, j  1, k )  E p (i, j  1, k ) E p (i, j , k  1)  E p (i, j , k  1)

(3.33)

 6 E p (i, j , k ) , ( p  x, y, z )3.6.5

Digital Filter

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

inew, j ,k 

cornerMi , j ,k  S iside, j , k  K i , j , k

M  6S  12 K

(3.34)

 i , j ,k を格子点上に定義される値とし、 inew
, j , k は Digital Filter を通した新しい  i , j , k を表
corner
す。また、  iside, j , k と  i , j , k は side term と corner term であり、それぞれ以下の式で表さ

れる。

22

 iside, 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(3.35)

 icorner, 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(3.36)

  i 1, j 1,k   i 1, j 1,k   i 1, j 1,k
ここで、 M 、 S 、 K の値をそれぞれ 12、4、1 とすると(2.34)式は次式となる。
1
4

inew
i , j ,k 
, j ,k 

1 side1
 i , j ,k   icorner, j ,k12
48(3.37)

数密度、イオン電流密度および電場を計算した後に、それぞれの値にこの Digital Filterを作用させる。

3.6.6

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

ハイブリッドコードにおいて、差分化の時間幅 t 、空間幅 x は次の条件を満たす必要がある。
CFL（Courant-Friedrichs-Lewy）条件
x(3.38)t 
VA
これは、t , x  0 で解が収束する条件である。V A はアルフベン速度であり次式で定
義される。

VA 

B0
 0 ni mi(3.39)

また、付加的な数値的安定条件は以下である。

 ci t  0.2(3.40)

ここで、  ci はイオンサイクロトロン周波数である。
メッシュ幅の最小値を求める条件として以下の式がある。

x  c  pi

(3.41)

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

第4章

計算体系と結果

計算体系の概略図を Fig.4-1 に示す。コイルの右端がターゲット（原点）から zcenter（25mm）だけ離れている。コイルの半径は外径が 30mm、内径が 25mm、X-Y 平面に平行である。ターゲットは薄膜で、計算領域は Z 軸方向に-50～50mm、X,Y 軸方
向にそれぞれ-50mm～50mm の立方体とした。最初にカップリング時の入射角度による比較、次にコイルとターゲットの距離による比較、そして、コイルの大きさによる比較を行った。計算は、100000 time step の 5000[ns]まで実行した。平板にした理由としては一つ目に、計算が容易であること。二つ目に、レーザーが当たっている
その反対側、背面におけるプラズマの広がりを考えたかったからである。
（直接あたった部分だと、膨張、粒子のスピードが速くて捉えにくく、すぐに磁場の影響が少ないところに行ってしまう。）
最後に１月に大阪大学で行った実験と近い体系にて計算を行い、実験結果との比較を目指した。
評価としては、推力（ibit : impulse bit）を比較して行った。推力の算出は、3D ハ
イブリッドコードを用いて次のような計算を行った。
𝑖𝑏𝑖𝑡(𝑖𝑚𝑝𝑢𝑙𝑠𝑒 𝑏𝑖𝑡) = ∑{𝑚 𝑧 − 𝑚 𝑧 (𝑡 = 50[ 𝑠])}Z 軸方向の運動量の和から初期（3D ハイブリッドコードにデータを入力した時点：
50[ns]）における運動量を引いて、初期状態の推力を 0 として計算した。

Figure 4-1

計算体系の概略図
24

1. カップリング時の入射角度（30°、45°、60°）による推力の比較を行った。

計算体系 1-1
Irradiation diameter(rspot)
500𝜇𝑚

25𝑚𝑚

30°
Angle of rotation
Coil
7.0 × 10−5 𝑐𝑚
Thickness of a target

Figure 4-2 計算体系 1-1（30°傾けて入射した場合）

Thickness of a target：0.7𝜇𝑚
Target：plastic
Laser wavelength：1.06𝜇𝑚
Focal position：−0.520𝑐𝑚

Intensity：0. 92 × 1012 𝑊/𝑐𝑚2
F number：20
Target mesh number：200
Energy4.6𝐽

Table 2 計算体系 1-1 における各種パラメータ25

9𝑛𝑠

Figure 4-3 計算体系 1-1 における粒子の動き
Fig.4-3 は、カップリング時において、z 軸方向に対して 30°傾けて計算したものである。
26

計算体系 1-2

Irradiation diameter(rspot)
25𝑚𝑚

Laser

500𝜇𝑚

45°
Angle of rotation
Coil

Thickness of a target
7.0 × 10−5 𝑐𝑚

Figure 4-4 計算体系 1-2（45°傾けて入射した場合）

ここでレーザー強度などのパラメータは、30°の時と同じものを使用して計算した。

27

Figure 4-5 計算体系 1-2 における粒子の動き
Fig.4-5 は、カップリング時において、z 軸方向に対して 45°傾けて計算したものである。
28

計算体系 1-3

Irradiation diameter(rspot)

Laser

25𝑚𝑚
Angle of rotation
60°

500𝜇𝑚

Coil

Thickness of a target

7.0 × 10−5 𝑐𝑚

Figure 4-6 計算体系 1-3（60°傾けて入射した場合）

ここでレーザー強度などのパラメータは、30°の時と同じものを使用して計算した。

29

Figure 4-7 計算体系 1-3 における粒子の動き
Fig.4-7 は、カップリング時において、z 軸方向に対して 60°傾けて計算したものである。
30

4.50E-06
4.00E-06
3.50E-06
3.00E-06
2.50E-06
2.00E-06
1.50E-06
1.00E-06
5.00E-07
0.00E+00-5.00E-07-1.00E-06

30°
45°
60°

0
3.50E-07
7.00E-07
1.05E-06
1.40E-06
1.75E-06
2.10E-06
2.45E-06
2.80E-06
3.15E-06
3.50E-06
3.85E-06
4.20E-06
4.55E-06
4.90E-06

ibit [Ns]impulse bittime [s]Figure 4-8 角度による推力の比較

Angle between V and B [degree]

100
90
80
70
60
50

30°

40

45°

30

60°

20
10
0

-2.00E-02 -1.50E-02 -1.00E-02 -5.00E-03 0.00E+00 5.00E-03

Z [m]Figure 4-9 速度と磁場ベクトルの間の角度比較31

結果としては、
45°の時に最も高い推力が出ていることが Fig.4-8 から見て取れる。
その理由として考えられるのは、45°の時、粒子の運動がコイル磁場に対してより垂直になっているからではないかと考えられる。そうすることでラーマー運動により発生するコイル磁場と逆向きのプラズマ磁場が強くなり、より磁場との相互作用が強くなって推力が高まったと推察できる。実際に 50ns 時（カップリング）の粒子の速度ベクトルとその粒子の位置における磁場ベクトルを Fig.4-9 に比較してみた。確かに 30°より 45°の方が、角度が垂直に近いが、60°よりは劣っていた。そうであるにも関わらず、45°の時に推力が高いのは、60°の時には、粒子の動きが磁場に沿って曲がっているため、推力があまり出なかったのだと考えられる[Fig.4-10 参照]。今後
はターゲット形状を変更したり、他の角度でカップリングさせたりすることで、様々な体系での計算を行いたい。それによりさらに良い条件を見つけていきたい。

プラズマが運動している場所は磁場が小さくなっている。
2 つの図を見比べると、粒子がコイル磁場に沿って運動していることが確認できる。

Figure 4-10 5050ns における速度と磁場ベクトルの比較（計算体系 1-3）
左：磁場ベクトル

右：粒子の速度ベクトル32

2. コイルとターゲットとの距離（5mm,15mm,25mm,35mm）による推力の比較を行った。

計算体系 2-1
Irradiation diameter(rspot)
5𝑚𝑚

Laser

500𝜇𝑚

45°
Angle of rotation
Coil

Thickness of a target
7.0 × 10−5 𝑐𝑚

Figure 4-11 コイルとターゲットとの距離が 5mm の場合

Thickness of a target：0.7𝜇𝑚
Target：plastic
Laser wavelength：1.06𝜇𝑚
Focal position：−0.520𝑐𝑚

Intensity：0. 92 × 1012 𝑊/𝑐𝑚2
F number：20
Target mesh number：200
Energy4.6𝐽

Table 3 計算体系 2-1 における各種パラメータ339𝑛𝑠

コイルからすぐに離れてしまい、磁場の影響の少ないところまで飛んで行ってしまう。

その分磁場によって戻ってくる
粒子が少ない。

Figure 4-12 計算体系 2-1 における粒子の動き
Fig.4-12 は、コイルとターゲットの距離が 5mm の時（比較している中で１番近い）の結果である。
34

計算体系 2-2

Irradiation diameter(rspot)
15𝑚𝑚

Laser

500𝜇𝑚

45°
Angle of rotation
Coil

Thickness of a target
0.7 × 10−4 𝑐𝑚

Figure 4-13 計算体系 2-2（コイルとターゲットとの距離が 15mm の場合）

ここでレーザー強度などのパラメータは、5mm の時と同じものを使用して計算した。

35

5mm の時よりは多少戻ってきている。

Figure 4-14 計算体系 2-2 における粒子の動き
Fig.4-14 は、コイルとターゲットの距離が 15mm の時（比較している中で 2 番目に近い）の結果である。
36

計算体系 2-3

Irradiation diameter(rspot)

Laser

500𝜇𝑚

5𝑚𝑚

45°
Angle of rotation
Coil

Thickness of a target
0.7 × 10−4 𝑐𝑚

Figure 4-15 計算体系 2-3（コイルとターゲットとの距離が 35mm の場合）

ここでレーザー強度などのパラメータは、5mm の時と同じものを使用して計算した。

37

この辺りの粒子
は磁場の影響が少ない。

粒子の量は多い。しかし、磁気ノズルも広がり、粒子がそれに沿って広がっているので、推力が低い。

Figure 4-16 計算体系 2-3 における粒子の動き
Fig.4-16 は、コイルとターゲットの距離が 35mm の時（比較している中で 1 番遠い）の結果である。

38impulse bit

5mm
15mm

ibit [Ns]4.50E-06
4.00E-06
3.50E-06
3.00E-06
2.50E-06
2.00E-06
1.50E-06
1.00E-06
5.00E-07
0.00E+00-5.00E-07-1.00E-06

25mm

04.00E-07
8.00E-07
1.20E-06
1.60E-06
2.00E-06
2.40E-06
2.80E-06
3.20E-06
3.60E-06
4.00E-06
4.40E-06
4.80E-06

35mm

time [s]Figure 4-17 距離による推力の比較

Fig.4-17 は、同じ 45°でカップリングさせた時、コイルとターゲットの距離が推力にどう関係しているか見るため、距離を 10mm ごとに変えて計算結果を比較したものである。25mm は計算体系 1-2 と同じである。
この比較により、25mm の時が、一番推力が高くなった。理由として考えられるのは、距離が近すぎたり、遠すぎたりすると、粒子に与える磁場の影響が少なくなり、推力につながらなかったことが挙げられる。近すぎれば、コイルを通り抜けた粒子がそのまま流れて行ってしまい、遠すぎればコイルより手前側（ｚ軸方向で正の方向）の粒子が、磁場の影響の少ないところにすぐに行ってしまう。さらにコイルが離れていることで磁気ノズルも広がり、それに伴って粒子が広がるので計算領域内において z 軸方向の速度成分が減ってしまう。25mm は、ちょうどコイルの内径と等しく、より多くの粒子を磁場の影響下におくことができる距離であると考えられる。実験においても、プラズマに対してより磁場の影響を及ぼせるように（粒子がコイル磁場の中に納まると予想して）、コイルとターゲットの距離はコイルの内径と等しい。つまり実験におけるコイル・ターゲット間の距離に関しての条件は、推力を高くするのに良い条件であることが計算によっても確認できた。
39

3.

コイルの内径（12.5mm,25mm）による推力の比較を行った。

計算体系 3-1
Irradiation diameter(rspot)

Laser

500𝜇𝑚
10𝑚𝑚

45°

Coil

Angle of rotation
Thickness of a target
0.7 × 10−4 𝑐𝑚

Figure 4-18 計算体系 3-1（コイルの大きさが 12.5mm の場合）Thickness of a target：0.7𝜇𝑚
Target：plastic
Laser wavelength：1.06𝜇𝑚
Focal position：−0.520𝑐𝑚

Intensity：0. 92 × 1012 𝑊/𝑐𝑚2
F number：20
Target mesh number：200
4.6𝐽
Table 4 計算体系 3-1 における各種パラメータ409𝑛𝑠

レーザー照射により広がった正面部分が戻ってきた粒子に押され、このように崩れた形になっている。

コイルが小さいため粒子が外側に広がっている。

コイルが小さいため、通り抜けずに返ってくる粒子が多い。

外側の粒子がコ
イルに巻きつくように運動する。

Figure 4-19 計算体系 3-1 における粒子の動き
Fig.4-19 は、計算体系 1-2 のコイルの半分の大きさ（12.5mm）のもので計算した結
果である。
41

1.40E-05

impulse bit

1.20E-05
1.00E-05

ibit [Ns]8.00E-06
6.00E-06

12.5mm

4.00E-06

25mm

2.00E-06-2.00E-06

04.00E-07
8.00E-07
1.20E-06
1.60E-06
2.00E-06
2.40E-06
2.80E-06
3.20E-06
3.60E-06
4.00E-06
4.40E-06
4.80E-06

0.00E+00

time [s]Figure 4-20 コイルの大きさによる推力の比較

この比較は、45°でカップリングした時のコイルの大きさを変えて推力を比較したものである。グラフからも明らかなように、コイルが小さい時の方が、推力が上がっている。この理由としては、コイルが大きい時に比べて小さい時、コイルを通り抜けていく粒子の量が少ないのだと考えられる。それは単純にコイルの中心部分が小さくなり、粒子がそこを通り抜けにくくなっていると考えられるからだ。
ただ、コイルを小さくすれば推力が必ずしも上がるとは考えていない。ターゲット
が平板の時、粒子は STAR1D を使っているという理由もあるが、極端にレーザーの照射方向に広がっていく。つまり角度の広がりが少なく、コイルの中心部分から通り抜けていきやすいということである。ターゲットに球を使うなどすれば、通り抜ける粒子の量も変わり、コイルを小さくするメリットはなくなる。さらに、コイルが小さくなれば、磁場の影響する範囲も小さくなる。
あくまでもこの体系においては、コイルを小さくした方が推力も高くなることを確
認できた。

42

4.

本研究室では、大阪大学と共同でレーザー生成プラズマの磁気中の挙動について実験を行っている。今年、2012 年の１月に大阪大学において激光Ⅻ号を用いた実験を行った。その実験結果との比較をするため、激光Ⅻ号での実験体系に近い体系での計算を行った。

計算体系 4

Coil

500𝜇𝑚

25𝑚𝑚
Target2 Target1

Laser
(2)

(1)

7.0 × 10−4 𝑐𝑚
Thickness of targets
0.5 , 2.0 × 10−4 𝑐𝑚
Figure 4-21 計算体系 4（激光Ⅻ号に近い体系）

実験において、実際にはターゲット（1）はシェル状であり、薄膜ではないが、STAR1D においてシェルは計算できないため、同じ厚さの薄膜で代用した。ターゲット（1）と（2）の間は、密度の低い（プラズマに影響を与えない程度：
1.06 × 10−8 𝑔⁄𝑐𝑚3 ）厚さ1𝑚𝑚のプラスチックを挟んでいる。

43

Coil

25𝑚𝑚
TargetLaser

7.0 × 10−4 𝑐𝑚

Thickness of targets
0.5 , 2.0 × 10−4 𝑐𝑚
Figure 4-22 激光Ⅻ号による実験の概念図

Thickness of a target1：7.0𝜇𝑚
Thickness of target2：0.5𝜇𝑚 , 2.0𝜇𝑚
Target：plastic
Laser wavelength：1.06𝜇𝑚
Focal position：−0.75𝑐𝑚
Intensity：1.76 × 1014 𝑊/𝑐𝑚2
F number：3

Mesh number：1300
4.6𝐽
Table 5 計算体系 4 における各種パラメータ9𝑛𝑠

44

Figure 4-23 計算体系 4 における粒子の動き（Target2 が厚さ𝟎. 𝟓𝛍𝐦の時）

45

Figure 4-24 計算体系 4 における粒子の動き（Target2 が厚さ𝟐. 𝟎𝛍𝐦の時）

46

1.00E-03
9.00E-04
8.00E-04
7.00E-04
6.00E-04
5.00E-04
4.00E-04
3.00E-04
2.00E-04
1.00E-04
0.00E+00

0.5
2

08.00E-07
1.60E-06
2.40E-06
3.20E-06
4.00E-06
4.80E-06
5.60E-06
6.40E-06
7.20E-06
8.00E-06
8.80E-06
9.60E-06

ibit(Ns)

impulse bittime(s)Figure 4-25 ターゲット（2）の厚さによる比較
Fig.4-25 の比較は、激光Ⅻ号での実験条件において、ターゲット（2）の厚さを変更して計算し、その結果を表したものである。計算時間は、前記の体系の倍である
200000 time step、10000[ns]まで実行した。それでも計算時間が若干足りず、推力が収束していない。このグラフを見る限りターゲット（2）が薄い方が、推力が出ていることが分かる。
Fig.4-23 や Fig.4-24 を見比べてもあまり違いが見えてこない。負の向きに移動する粒子による推力に関しては、Fig.4-26 から分かるように、厚さによる違いがあまり見られなかったため、正の向きに移動する粒子数の差を計算した。すると、薄い方が厚い方より約 3%多かった。z 軸に正方向の粒子の速度合計を比較してみると、薄い方が約 4.9%多かった。推力の差は薄い時の推力の方が約 5.5%多かった。つまり、粒子の量の違い（ターゲット（2）が薄い方が z 方向に正の向きの粒子を多く持っていること）よりも、戻ってくる粒子が持つ速度の違い（厚いターゲットよりも薄い方が速いということ）の方が、より推力の差に影響していると考えられる。
この結果を実験と比較しようと試みたが、実験の推力計測がうまくいかず、結果が芳しくないものであった。そのため比較ができなかった。
この結果についてはまだ計算時間が足りていない。今後はこの計算について、さらに計算時間を伸ばして推力が収束するまで計算し、比較していきたいと考えている。

47

正の方向に戻ってきている粒子数が薄い方が多いという結果よりも、戻ってきている粒子の速度が、薄い方が速いという結果の方が、より推力に影響していると考えられる。

粒子数に 1322 個の違いがある。薄い方が厚い方より約 3%多い。

Target2の厚さ z軸の正方向に動く粒子の数[個] 負の方向に動く粒子の数（個）
0.5µｍ
44606
55394
2.0µｍ
43284
56716
速度（正方向の粒子）の合計[m/s] 推力の比較（10050ns）[Ns]0.5µｍ1838662328
0.003869864
2.0µｍ
1752311956
0.003668876
Table 6 Target2 の厚さによる比較

4.00E-03
3.90E-03
3.80E-03
3.70E-03
3.60E-03
3.50E-03
3.40E-03
3.30E-03
3.20E-03
3.10E-03
3.00E-03
2.90E-03

推力には2.01 × 10−4 [Ns]の違いがある。薄い時の推力の方が約 5.5%多い。

0.5

ibit− (Ns)

-2.70E-03-2.80E-03

09.50E-07
1.90E-06
2.85E-06
3.80E-06
4.75E-06
5.70E-06
6.65E-06
7.60E-06
8.55E-06
9.50E-06-2.60E-03-2.90E-03-3.00E-03
0.5-3.10E-03

22-3.20E-03

08.50E-07
1.70E-06
2.55E-06
3.40E-06
4.25E-06
5.10E-06
5.95E-06
6.80E-06
7.65E-06
8.50E-06
9.35E-06

ibit+ (Ns)

速度には8.6 × 107 [m/s]の違いがある。
薄い方が約 4.9%多い。

-3.30E-03-3.40E-03

time(s)time(s)Figure 4-26 ターゲット（2）の厚さによる比較
左：z 軸に正の方向の推力

右：z 軸に負の方向の推力

48

第5章

まとめ、今後の課題

計算 1
ターゲットを 30°、45°、60°傾けてカップリングさせ比較した結果は 45°が、一番推力が高くなった。今後他の様々な角度での計算を行っていきたい。

計算 2
コイルとターゲットの距離の比較については、コイルの内径と等しい距離が推力が高くなることが分かった。これにより、現在行っている実験条件も、推力を出すのに良い条件であることが確認できた。今後は他の体系についてもこの距離の比較を行っていきたいと考えている。

計算 3
薄膜平板をターゲットとした体系においては、コイルを小さくした方が、推力が高くなることが確認できた。他の体系においてはどうなるのか今後検証していきたい。
まとめとしては、実験において傾けて入射するのは、粒子がコイルに当たるなどしてしまうため難しい。そのためコイルとターゲットの距離を考えつつ、他のターゲット形状においても計算を繰り返し、より最適な条件を見つけていきたい。

計算 4
激光Ⅻ号での実験結果との比較は、実験の推力計測がうまくいかず、比較ができなかった。今後実験グループとのコミュニケーションをとり、実験、計算の方針を検討していきたい。

49

参考文献[1] R. A. Hyde, L. L. Wood, Jr.,J. H. Nuckolls: AIAA Paper, No.72-1063, (1972)[2] R. A. Hyde: UCRL-88857, (1983)[3] Y. Nagamine, H. Nakashima: Analysis of Plasma Behavior in a Magnetic Thrust
Chamber of Laser Fusion Rocket,Fusion Technology, 35 (1999), 62[4] K. V. Vchickov, H. Nakashima, F. Ichikawa, Y. P. Zakharov: Vacuum, 73 (2003), 427[5] K. V. Vchickov, H. Nakashima, Y. P. Zakharov, T. Esaki, T. Kawano,
T. Muranaka: Jpn. J. Appl. Phys., 42 (2003), 6590
[6] 川渕亮,磁気ノズルを用いたレーザー核融合ロケットの推力ベクトルの制御技術,学士論文,九州大学,(2006)[7] 松田伸夫,整形ターゲットを用いる磁気ノズル内のプラズマ挙動解析,修士論文,九州大学大学院,(2009)[8] A.Maeno, N.Yamamoto, H.Nakashima: J. Plasma Fusion Res. Vol.86, No.10 (2010), 595[9] A.Ando: J. Plasma Fusion Res. Vol.83, No.11 (2007), 276[10] N.Uchida: J. Plasma Fusion Res. Vol.83, No.11 (2007), 271[11] C. D. Orth: VISTA - A Vehicle for Interplanetary Space Transport Application Powered by
Inertial Confinement Fusion, UCRL-LR-110500 (2003)[12] A.Sunahara, et al.: J. Plasma Fusion Res. Vol.83, No.11 (2007), 920-921[13] Lyman Spitzer Jr. and Richard Harm, Phys. Rev. 89, (1953), 977[14] R.C. Malone, R.L. McCrory and R.L. Morse, Phys. Rev.Lett. 34, 721 (1975).
[15] Lyman Spitzer Jr., Physics of Fully Ionized Gases (John
Wiley & Sons Inc., 1962) p. 115.
[16] G.N. Minerbo, J. Quant. Spectrosc. Radiant,20, (1978), 541[17] D. S. Harned: J. Comp. Phys., 47 (1982), 452[18] E. J. Horowiz, D. E. Schumaker, and D. V. Anderson: J. Comp. Phys., 84 (1989), 279[19] A. Mankofsky, R. N. Sudan, and J. Denavit, J. Comp. Phys., 70 (1987), 89[20] C. K. Bridsall and A. B. Langdon: "Plasma Physics via Computer Simulation", IOP
Publishing Ltd, (1997)50

謝辞
本研究において、あらゆることを御指導頂いた中島秀紀教授に深く感謝致します。
研究室生活において身近で御指導して頂いた山本直嗣准教授に深く感謝致します。快
適な研究環境を提供して下さった大神めぐみ秘書に深く感謝致します。
本研究にて使用させてもらった計算コードに携わっていらっしゃる砂原淳先生、城崎知至先生に深く感謝致します。
同じ LFR 組として計算についての多くのアドバイスを頂き、大阪大学での実験を共に行った安永将人氏、花屋倫生氏、秀平靖磨氏、日永智之氏に深く感謝いたします。研究についてのアドバイスをしていただいたり、楽しい研究室生活にしたりして頂いた
諸先輩方、下川床潤氏、杉田健策氏、平岡義章氏、中村祐輔氏、平野賢治氏、廣池匠哉
氏に深く感謝致します。遠いインドの地からこの研究室に来てくださった Indranuj Dey 氏に深く感謝いたします。同じ学年としてアドバイスし合い、共に日々頑張ってきた坂本雄樹氏、笹川裕太郎氏、三輪祐大氏に深く感謝致します。
特に前野旭弘氏には、自分の研究の合間を縫って、研究に不可欠なコードや計算機のことなど多大なる知識やアドバイスを頂き、普段の生活においても、車について面倒を見てくださり、私だけでなく、研究室全体を明るくしてもらい深く感謝致します。
最後に、これまで支えてくださった家族、友人に深く感謝致します。

付録 1
STAR1D において、投入したエネルギーは
𝑡𝑜𝑡𝑎𝑙𝐸 = 4.0005 629295792 [𝐽]カップリングコードから 3D ハイブリッドコードにデータとして渡す速度から得た粒子の運動エネルギーは、（3 章の通り）
𝐸 = 4.640200 71817542 [𝐽]となっており、投入したエネルギーとカップリング時のエネルギーに差がある。これは、STAR1D 中のエネルギー投入方法によるものではないかと考えられる。
Fig.にその投入方法の概念図を示す。

𝑙𝑎𝑠𝑒𝑟

𝑡𝑎𝑟𝑔𝑒𝑡Figure A
図のようにターゲットが薄膜の時、投入するにはレーザーを何本も入射するような形となる。それぞれが投入エネルギー（前記の例の場合、4J）を持っている。そのレーザー同士が図のように重なってしまい、投入したエネルギーよりもカップリング時のエネルギーの方が大きくなっているのではないかと考えられる。コードの性質として、ターゲットに薄膜を使う場合には、この点は課題として残る。

A