卒

題

業

論

文

目

2 次元輻射磁気流体コードを用いた
磁気スラストチャンバ―におけるプラズマ挙動の解析

氏

名

児島

富彦

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

指 導 教 員

山本

直嗣

教授

森田

太智

助教

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

平成２９年２月

目次
第1章

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

1.1

研究背景 ...............................................................................................................................1

1.2

磁気スラストチャンバーの推力発生原理 ..............................................................................2

1.3

先行研究 ...............................................................................................................................3

1.4

研究目的 ...............................................................................................................................8

第2章

輻射磁気流体計算コード.....................................................................................................9

2.1

PINOCO の適応例 ................................................................................................................92.2

計算体系 .............................................................................................................................10

2.3

輻射流体方程式 ................................................................................................................... 11

2.4

一般極座標変換 ...................................................................................................................122.5

イオン・電子の熱伝導式 .....................................................................................................132.6

磁場の式 .............................................................................................................................142.7

初期磁場の設定 ...................................................................................................................15

2.8

境界条件 .............................................................................................................................17

第3章

計算手法 ...........................................................................................................................18

3.1

CIP 法.................................................................................................................................18

3.1.1 CIP 法の説明 ..................................................................................................................18
3.1.2
3.2

PINOCO での表記........................................................................................................20

ALE 法................................................................................................................................21

3.2.1

ALE 法の概要...............................................................................................................21

3.2.2

PINOCO における ALE ...............................................................................................21

第4章

計算結果 ...........................................................................................................................25

4.1

100 J×6 ビーム ...................................................................................................................25

4.2

レーザーエネルギー依存性 .................................................................................................30

4.2.1

100 J ............................................................................................................................30

4.2.2

50 J ..............................................................................................................................33

4.2.3

10 J ..............................................................................................................................354.3

磁気圧の効果の確認 ............................................................................................................38

第5章

考察 ..................................................................................................................................415.1

実験結果との比較 ...............................................................................................................41

5.2

計算で生じた問題点のまとめ ..............................................................................................45

第6章

結論・今後の展望 .............................................................................................................46

参考文献 .........................................................................................................................................47
謝辞 ................................................................................................................................................49

第 1 章 序論
1.1 研究背景2016 年 9 月に行われた IAC(国際宇宙会議)にて、米国の宇宙企業 SpaceX の創業者イーロン・マスク氏が「Making Humans a Multiplanetary Species」と題して壮大な火星への移住計画を発表した。同氏の計画では地上から地球周回軌道上まで移民を乗せたロケットを打ち上げ、軌道上で推進剤を補給し火星を目指す。他惑星への移住は１つの惑星に縛られて生活することからの脱却という面と、いずれ人類が直面するエネルギー
資源の枯渇問題や隕石の衝突による急激な気候変動で生活困難となる際の解決策 となる面など多くの利点がある。計画でのロケットの推進剤にはメタンを用いており、火星資源からメタンを生成することも予定されている。しかし、そのような化学推進ロケットの推力では火星到達までに 8 ヶ月、往復で 500 日以上かかるとされている(1)。また比推力、即ち単位推進剤あたりに得られる推力が低いため大量の推進剤が必要となる。ミッションの長期化に伴う宇宙飛行士の身体的負担の軽減の面でも、目標達成のためには大推力・高比推力を同時に達成する宇宙船が必要である。推進剤量と身体的負担を減らす 宇 宙 船 と し て レ ー ザ ー核 融合 ロ ケッ ト が提 案 され た (2) 。 そ して 米国 Lawrence
Livermore 国立研究所の Orth らのグループがレーザー核融合ロケットの概念例として
火星往復をわずか 145 日程度で実現できる「VISTA」を提案した(3)。図 1.1 の概念図で示すように、レーザー間接照射法で生成した核融合プラズマと超伝導コイルで生じる磁場との相互作用で図面左方向に推進する。なお、核融合によって生じる中性子線による宇宙船の損耗を減らすため中抜けの円錐の形状をしてあり、外縁部に居住区・ペイロー
ドを設ける。しかし、その推進原理の実証はされていなかった。そこで本研究室では推進機構を磁気スラストチャンバーと呼び、模擬実験による検証と数値的検証の双方を行ってきた。模擬実験では大阪大学レーザーエネルギー学研究センターの共同研究施設を利用して、レーザー生成プラズマと磁場との相互作用を確認する実験を行い、数値的検証では推進理論の検証、模擬実験の数値的再現を試みてきた。

1図 1.1 「VISTA」概念図(4)1.2 磁気スラストチャンバーの推力発生原理
図 1.2 で磁気スラストチャンバーの推力発生原理を説明する。まず、コイルに電流を流し外部磁場を発生させる。その外部磁場中でターゲットにレーザーを照射し慣性核融合を起こし、核融合プラズ
マを生成する。膨張するプラズマによって磁場は押し広げられ、圧縮される。そして、磁場が圧縮される領域では磁場のエネルギー（磁気圧）が増え、やがて膨張するプラズマの運動エネルギーとつり合う。その後、磁場はバネのような働きをし、プラズマを押し返す。押し返されたプラズマは図面右方向に加速され、排出される。その反作用で推力を得る。

図 1.2 磁気スラストチャンバーの推力発生原理

21.3 先行研究
本研究グループはこれまで、MHD コードを用いたプラズマの挙動解析(5)、さらにイオンを粒子、電子を流体と仮定した 3 次元ハイブリッドコードを用いた数値計算(6)で磁気スラストチャンバーの挙動の解析を行ってきた。また核融合プラズマの代わりにレー
ザーアブレーションプラズマを利用した模擬実験による検証を行ってきた。そして、実験と数値計算におけるプラズマの分布でほぼ一致が見られた(7)。さらに実験においてコイル磁場とプラズマの相互作用により、推力が発生することが確認された(8)。しかし、それまでの数値計算ではプラズマがすでに広がった状態から計算していたため、模擬実
験の再現性は不十分であった。そこでグループはレーザー入射過程も計算に考慮するため、1 次元輻射流体コード(Star1D)を導入した。即ちレーザーアブレーションプラズマ
生成過程を Star1D で計算し、その後、磁場とプラズマの相互作用を 3 次元ハイブリッ
ドコードで計算した。その複数のコードを組み合わせた Star1D+3 次元ハイブリッド計
算で模擬実験との比較を行った(9)が、実験との相違が生じた。その解決のためグループ
は 2 次元輻射流体コード(Star2D)を導入した。それまで Star1D の結果を 3 次元ハイブ
リッドコードに受け渡す際に、球対称を仮定して 1 次元から 3 次元にデータを拡張してきたが、Star2D ではすでに 2 次元の結果を持っているため、入射レーザーに対して軸対称で拡張することができる。多本数のレーザーを用いてターゲット表面全体にレーザーを照射する実験では、アブレーションプラズマは等方膨張するため、1 次元から球対称を仮定して拡張していい。しかし、1 本のレーザーでプラズマを生成し、プラズマの空間分布を調べる際には、レーザー照射面と裏面でアブレーションされる量が異なるた
め 2 次元で計算する必要がある。図 1.3 が大阪大学レーザーエネルギー学研究センターの EUV レーザーを用いて行われた実験であり、レーザーアブレーションプラズマの角度分布が検証された。図 1.4 が実験結果と 1 次元輻射流体コード Star1D と、当時新たに導入した Star2D との比較である(10)。なお、数値計算ではある角度の位置にイオンがいくつ存在しているかを求めイオン電流に換算している。即ち、ここでのイオン電流はプラズマ中のイオン分布である。実験値では 0°から角度が大きくなるほどイオン電流が大きくなり、45°で極大を取っている。Star1D ではその傾向が見えなか った が、
Star2D では 60°で極大が見られ、角度依存性はある程度示す事ができた。

3

図 1.3 EUV チャンバー実験体系図(10)レーザー入射線上を 0°としチャージコレクタを A=－14°,
A’=14°,E=45°,F=60°,G=75°の位置でイオン電流を測定

図 1.4 実験結果(黒)と Star1D(青),Star2D(赤)の計算結果(10)4図 1.5 は 15°における実験で計測したイオン電流の時間発展
（破線、
y 軸左）
と Star2D+3
次元ハイブリッド計算で得られた、同じ計測位置でのイオン粒子数の時間発展（ 実線、y 軸右）である。実験値と計算値をイオン電流密度に換算して比較したところ、計算値の方が小さかった。その定量的相違が生じた原因には(1)遠方で計測しているため、イオ
ン―電子再結合によりイオン密度が減少していること、(2)数値計算コードの再現性がよくないことが考えられた。

図 1.5 15°における実験で計測したイオン電流(破線,左軸)と数値計算においてその位置に存在するイオン超粒子の個数(実線,右軸)の時間変化
次にシミュレーションフローについて説明する。図 1.6 に示すように、従来の数値計
算では複数の数値計算コードを用いて 1.2 節で述べたプラズマと磁場の相互作用を検証してきた。まず 1 次元または 2 次元輻射流体計算を行いレーザーによるプラズマ生成過程を計算する。そして、生成されたプラズマの内部エネルギー、運動エネルギー、放射エネルギーのうち運動エネルギーが支配的になる時点（レーザー入射から 75 ns より後(9) ）でのイオンの位置・速度と外部磁場データを初期値として 3 次元ハイブリッド計算

によりプラズマと磁場の相互作用（～数 µs）を計算する。ここで運動エネルギーが支配的になる時点を初期値とする理由は 3 次元ハイブリッドコードではイオン温度、電子温
度を 0 としているためである。

図 1.6 従来の数値計算の流れ
5しかし、輻射流体計算において、アブレーションプラズマの膨張を捉えきれていなかった。図 1.7 が Star2D でプラズマの生成過程を計算したもので、レーザー照射から 30 ns
でのイオン数密度の空間分布である。速度 150 km/s 以上のアブレーションプラズマが計算領域の境界まで達している。そのため 3 次元ハイブリッド計算に初期値を渡す時間
では、速いプラズマは計算領域から抜けており、実験の再現性として不十分である。その原因の 1 つには、本来は生成プラズマが広がる時点で外部磁場が圧縮され、磁場によるプラズマの減速が起きるが、Star2D では磁場を導入していないため考慮できなかったことが挙げられた。本来ならば、３次元ハイブリッド計算をする際にプラズマ生成時の相互作用で変化した磁場とプラズマの密度・速度を初期値として入れるべきで ある。

図 1.7 Star2D の計算によるイオン数密度の空間分布(11)

数値計算との定量的比較を目的としたパラメータ取得のため、ターゲット付近でのプラズマ計測が行われた。図 1.7 は 2015 年 1 月に行われた大阪大学レーザーエネルギー学
研究センターの激光Ⅻ号レーザーを利用して行われた Mach-Zehnder 干渉計測のセットアップである。レーザー条件は表 1 で行われた。図 1.8 がレーザーエネルギー100 J×
6 本で行った、レーザー入射から 200 ns 後の電子密度分布であり、1017 ―1018 個/cm3のアブレーションプラズマが計測された(12)。図 1.8 においてコイルは図の左側に位置し、そのコイル磁場によってアブレーションプラズマの先端が減速され、かつ遅れて広がった内部のアブレーションプラズマが堆積するという、特徴的シェル構造が確認された

6

図 1.8 干渉計測のセットアップ(12)表1

レーザー条件(12)入射レーザー本数6本

波長1053 nm

パルス幅

1.3 ns

レーザーエネルギー

10, 50, 100 J/本

図 1.9 干渉計測された電子数密度のグラフ(12)ターゲット中心は(0,0)の位置であり、プラズマのシェル構造が確認できる。
7

1.4 研究目的
1.3 節で述べた数値計算の問題点の解決に向け、従来行っていた輻射流体計算の代わりに輻射磁気流体計算を行い、一連の数値計算の再現性の向上を目指す。また、ターゲット近傍でのプラズマのパラメータの計測が既に行われており、実験値と直接比較を行うことで定量的な一致を目指す。実験を再現するシミュレーションコードを開発することで、プラズマと磁場のエネルギースケーリング則の構築につなげる。スケール則は将来のレーザー核融合ロケット開発に向け、高エネルギーにおける磁気スラストチャンバーの設計に役立つと考えられる。しかし、輻射磁気流体計算を磁気スラストチャンバーの数値計算に導入することは新たな試みであり問題が生じた。本研究ではその問題解決を行い、再現性の向上を目指す。
本論文の構成は以下の通りである。第 2 章では本研究に用いた輻射磁気流体 コード(PINOCO)の概要を述べる。第 3 章では PINOCO の計算手法を説明する。第 4 章では計算結果を述べる。第 5 章に考察、第 6 章に結論、今後の展望を述べる。

8

第 2 章 輻射磁気流体計算コード
2.1 PINOCO の適応例
本 研 究 に 用 い る 輻 射 磁 気 流 体 コ ー ド (precision integrated implosion numericalobservation code: PINOCO)は大阪大学レーザーエネルギー学研究センターの長友氏が開発・提供しているもので、レーザー核融合におけるレーザーとターゲットの輻射計算と、生成されるプラズマの挙動を流体計算するものである(13)。レーザー核融合の方法には大きく 2 つがある。1 つは中心点火方式、もう 1 つは高速点火方式である。大まかに説明すると、前者はレーザーによってターゲット表面をアブレーションさせ、その反作用にて等方的に爆縮させ核融合を起こす方法である。後者は前者と同様にレーザーによって爆縮され中心密度が上昇したターゲットに、高速電子を打ち込んで加熱し、核融合を起こすものである。さらに最近の高速点火方式では、磁場を加えることで電子の広がりを抑え、爆縮中心に高速電子を打ち込むという研究がされている。図 2.1 は長友氏が高速点火方式において PINOCO を用いて計算した例であり、これまで輻射流体計算していたコードに磁場を導入することで、磁場構造が爆縮の際にどの様に影響するかを示したものである(14)。PINOCO は輻射流体計算と磁気流体計算において成果を上げており、磁気スラストチャンバーのシミュレーションに応用して成果を上げると期待できる。

図 2.1 (a)(b)1.6 ns における磁場と質量密度、(c)(d)1.88 ns における磁場と質量密度(14)92.2 計算体系
計算は図 2.2 に示すように 2 次元極座標で行う。ターゲット中心から計算領域の外周までを 1.7 cm としている。これは 1.3 節の先行研究で述べた輻射流体計算における計算領域より十分大きい。メッシュ分割数は半径方向に 720、周方向に 240 としている。
また、ターゲットの内部と、ターゲットと背景領域の境界付近を細かく、外部を粗くして全体として不均一に取っている。
ターゲットの材質は実験と同じ CH とし、半径は 250µm としている。実験ではシェルターゲットが用いられたが、計算では簡略のため、中実球とする。また流体近似をするため、実験では真空である背景領域に低密度かつ低温の CH プラズマを背景ガスとして置いている。レーザーはメッシュに沿ってレイトレースを用いて照射される。実験のレーザービーム本数は有限であるが、本研究では 6 ビー
ム分のレーザーエネルギーが均等にターゲット表面に照射されると仮定している。コイ
ル半径は 1.8 cm、ターゲット中心からの位置を 1.8 cm としている。これは先行実験と同条件の磁場強度、即ちターゲット位置において 1.1 T となるように設定している。

図 2.2 計算体系10

2.3 輻射流体方程式
イオン、電子の 2 温度を仮定した輻射流体の方程式系は、プラズマ密度ρ、速度 u、イオンと電子の単位質量当たりのエネルギーεi、εe 、全圧力 p=pi+pe 、イオン圧力 pi、
電子圧力 pe などを用いて次のように書ける。なお計算には cgs 単位系を用いている。
質量保存の式𝑑𝜌= −𝜌𝛻 ∙ 𝒖
𝑑𝑡

(2.1)

𝑑𝒖
= −𝛻𝑝
𝑑𝑡

(2.2)

𝑑𝜀𝑖
= −𝑝𝑖 𝛻 ∙ 𝒖 − 𝛻 ∙ 𝑞𝑖 + 𝑄𝑒𝑖
𝑑𝑡

(2.3)

運動方程式
𝜌
イオンのエネルギー式
𝜌

電子のエネルギー式
𝜌𝑑𝜀𝑒= −𝑝𝑒 𝛻 ∙ 𝒖 − 𝛻 ∙ 𝑞𝑒 − 𝑄𝑒𝑖 + 𝑆 𝐿
𝑑𝑡

(2.4)

ただし、qi,e 、Qei、SL はそれぞれ、熱伝導による熱流束、電子―イオンの熱緩和、レーザー吸収項を表す。なお、添え字の i,e はそれぞれイオン、電子を示す。
輻射輸送方程式は周波数νに対する強度 Iν の放射係数ην 、吸収係数をχν として、光子に関する 6 次のボルツマン方程式
1 𝜕𝐼 𝜈+ 𝛺 ∙ 𝛻𝐼 𝜈 = 𝜂𝜈 − 𝜒 𝜈 𝐼 𝜈
𝐶 𝜕𝑡

(2.5)

で表される。計算負荷を考慮すると直接このボルツマン方程式系を解くことは困難であるから、等散乱を仮定し、0 次、および 1 次のモーメント式から次式の多群拡散近似
した輻射輸送エネルギー密度 E ν と輻射エネルギー流速 Fν の方程式を解く(15) 。
𝜕𝐸 𝜈+ 𝛻 ∙ 𝑭𝜈 = 4𝜋𝜂𝜈 − 𝜒 𝜈 𝑐𝐸 𝜈
𝜕𝑡

(2.6)

1 𝜕𝑭𝜈
+ 𝑐 ∙ 𝛻𝑷𝜈 = 𝜒 𝜈 𝑭𝜈
𝐶 𝜕𝑡

(2.7)
11

ただし、輻射エネルギー流束は、エディントン因子 f ν を用いて次式のようになる。
𝑐𝑭𝜈 = − 𝜒𝜈 ∇ ∙ (𝒇𝜈 𝐸 𝜈 )(2.8)

2.4 一般極座標変換
任意の(x,y)空間を扱えるよう座標系は一般曲線座標系(ξ,η)を用いる。
(2.9)

x=x(ξ(x,y),η(x,y)), y=y(ξ(x,y),η(x,y))とし、それぞれ x,y で微分すると
𝜉𝑥 𝑥 𝜉 + 𝜂𝑥 𝑥 𝜂 = 1 ← 𝜕𝑥/𝜕𝑥
𝜉𝑦 𝑥 𝜉 + 𝜂𝑦 𝑥 𝜂 = 0 ← 𝜕𝑥/𝜕𝑦𝜉𝑥 𝑦𝜉 + 𝜂𝑥 𝑦𝜂 = 0 ← 𝜕𝑦/𝜕𝑥{𝜉𝑦 𝑦𝜉 + 𝜂𝑦 𝑦𝜂 = 1 ← 𝜕𝑦/𝜕𝑦

(2.10)

すなわち、行列表示では𝑥𝜉[𝑦𝜉

𝑥 𝜂 𝜉𝑥𝑦𝜂 ] [𝜂𝑥𝜉𝑦1𝜂𝑦 ] = [00]1(2.11)

𝜉[𝜉𝑥𝑦

𝜂𝑥 𝑥 𝜉
𝜂𝑦 ] [𝑥 𝜂𝑦𝜉1𝑦𝜂 ] = [00]1(2.12)

あるいは

よって変換のマトリックスは𝑥𝜉
𝑑𝑥[ ] = [𝑦𝑑𝑦𝜉𝑥 𝜂 𝑑𝜉
𝑦𝜂 ] [𝑑𝜂 ](2.13)𝑦𝜂𝑑𝜉[ ] = 𝐽 [−𝑦𝑑𝜂𝜉−𝑥 𝜂 𝑑𝑥𝑥 𝜉 ] [𝑑𝑦 ](2.14)

この逆行列を作ると𝐽 −1 = 𝑥 𝜉 𝑦𝜂 − 𝑥 𝜂 𝑦𝜉

(2.15)

ここで、J はヤコビアンと呼ばれる。
𝐽≡

𝜕(𝜉, 𝜂) 計算空間のセル面積≡
𝜕(𝑥, 𝑦) 物理空間のセル面積(2.16)

ヤコビアンは上式で表されるように座標(x,y)で定義される物理空間のセル面積と一般座
標(ξ,η)で定義される計算空間のセル面積の比で、各座標間の単位面積比に相当する。
よって変換のメトリックスは12𝜉𝑥 = 𝐽 ∙ 𝑦𝜂𝜂𝑥 = −𝐽 ∙ 𝑦𝜉𝜉𝑦 = −𝐽 ∙ 𝑥 𝜂
𝜂𝑦 = 𝐽 ∙ 𝑥 𝜉(2.17)

なお、コーシー・リーマンの関係が成り立つ。
𝜉𝑥 = 𝜂𝑦 ,𝜉𝑦 = −𝜂𝑥

2.5 イオン・電子の熱伝導式
式(2.3)、式(2.4)における熱伝導項について説明する。フーリエの法則により x 方向は
𝜕𝑇
𝜕𝑥

(2.18)

𝜕𝑇
𝜕𝑇+ 𝜂𝑥 )
𝜕𝜉𝜕𝜂

(2.19)

𝜕𝑇
𝜕𝑇+ 𝑥𝜉 )
𝜕𝜉𝜕𝜂

(2.20)

𝑞𝑥 = 𝜅
座標変換をすると= 𝜅 (𝜉𝑥
ヤコビアンを用いて= 𝜅𝐽 (−𝑥 𝜂となる。y 方向も同様にして

𝜕𝑇
𝜕𝑦(2.21)

𝜕𝑇
𝜕𝑇+ 𝜂𝑦 )
𝜕𝜉𝜕𝜂

(2.22)

𝜕𝑇
𝜕𝑇+ 𝑥𝜉 )
𝜕𝜉𝜕𝜂

(2.23)

𝑞𝑦 = 𝜅= 𝜅 (𝜉𝑦= 𝜅𝐽 (−𝑥 𝜂よって∇∙ q =∂𝑞𝑥 1 𝑦𝜕𝑞𝑦+
∂x 𝑦 𝜕𝑦= 𝜉𝑥(2.24)𝜕𝑞𝑥𝜕𝑞𝜉 𝜕(𝑦𝑞𝑦 ) 𝜂𝑦 𝜕(𝑦𝑞𝑦 )+ 𝜂𝑥 𝑥 + 𝑦
+
𝜕𝜉
𝜕𝜂
𝑦 𝜕𝜉
𝑦 𝜕𝜂

(2.25)𝐽𝑥 𝜂 𝜕(𝑦𝑞𝑦 ) 𝐽𝑥 𝜉 𝜕(𝑦𝑞𝑦 )𝜕𝑞𝑥𝜕𝑞− 𝐽𝑦𝜉 𝑥 −
+
𝜕𝜉
𝜕𝜂
𝑦
𝜕𝜉
𝑦
𝜕𝜂

(2.26)= 𝐽𝑦𝜂

13熱伝導は古典的な Spitzer-Hӓrm モデルを用いており、電子熱伝導テンソルは𝛋 ∙ ∇𝑇𝑒 = 𝜅∥ ∇ ∥ 𝑇𝑒 + 𝜅⊥ ∇⊥ 𝑇𝑒 + 𝜅∧ 𝐡 × ∇𝑇𝑒(2.27)𝑛 𝑇𝜏𝜅∥ = 𝛾0 ( 𝑒 𝑒 𝑒 ) ≡ 𝛾0 𝜅0𝑚𝑒(2.28)(𝛾1′ 𝜒 2 + 𝛾0′ )𝜅⊥ = 𝜅0∆(2.29)(𝛾1" 𝜒 2 + 𝛾0" )∆(2.30)

𝜅∧ = 𝜅0を満たす。ここでの∇∥ 𝑇𝑒、∇ ⊥ 𝑇𝑒 はそれぞれ𝐡(𝐡 ∙ ∇𝑇𝑒 )、𝐡 × (∇𝑇𝑒 × 𝐡)を表す。また式中の
𝛾0 、𝛾0′ 、𝛾0" 、𝛾1 、𝛾1′ 、𝛾1" 、𝜒 は熱伝導係数である(16)。
レーザーによる熱伝導は、主にイオンと比べて質量が小さい電子に働き、吸収によって熱伝導が急勾配になり熱流束が急に増加してしまうが、電子は各々の平均自由行程以上には広がれないはずであるため、電子熱伝導においてフラックスリミッターを設定している。

2.6 磁場の式
磁場の時間変化は以下のように表される∂𝑩= ∇ × {𝑽 × 𝑩}∂t𝑐𝑚 𝑒 𝑐 𝛼0− [∇ × (𝒉{(∇ × 𝑩) ∙ 𝒉} )𝑒4𝜋𝑒 𝑛𝑒 𝜏𝑒(2.31)

𝑚 𝑐 1 − 𝛼′+∇ × ( 𝑒𝒉 × {(∇ × 𝑩) × 𝒉})]4𝜋𝑒 𝑛𝑒 𝜏𝑒
右辺第 2 項は拡散項である。B、c、e、me、ne、τe、h はそれぞれ磁場、光速、電気素
量、電子質量、電子数密度、イオン=電子衝突緩和時間、磁場方向の単位ベクトルを表す。また、α0、α’は Braginskii の拡散係数(17)である。
また、ローレンツ力を考慮した運動量保存則は
𝜌

𝑑𝒖
𝑩2
1
( 𝑩 ∙ 𝛻) 𝑩
= −𝛻 (𝑃𝑡ℎ + ) +
𝑑𝑡
8𝜋
4𝜋

となる。第 1 項は熱圧力と磁気圧、第 2 項は磁場の張力を表す。
14(2.32)

2.7 初期磁場の設定
ビオ・サバールの法則
𝑩(𝒙) =

𝜇0 𝐼 𝑑𝒙′ × (𝒙 − 𝒙′ )∫|𝒙 − 𝒙′ |34𝜋(2.33)

を利用し単ループのコイルが作る磁場を考える。まず、コイル軸上の磁場を計算 する。
図 2.3(左)のように円の中心を z 軸として、原点から磁場を求める点へ向かうベクトルを𝒙、コイル上の点へ向かうベクトルを𝒙′ とする。𝒙 − 𝒙′ は図のようになり、𝐼𝑑𝒙′ と𝒙 − 𝒙′ の外積は図のように斜め上向きとなる。これを一周分積分する。図 2.3(右)のようにコイルの半径を R、磁場求める点までの距離を Z とすると、𝑑𝒙′ は大きさ𝑅𝑑𝜙で円周向きのベク
トルとなる。𝑑𝒙′ と𝒙 − 𝒙′ は直交しているため外積は単なる掛け算となり、
𝑑𝑩(𝒙) =

𝜇0 𝐼 𝑑𝒙′ × (𝒙 − 𝒙′ ) 𝜇0 𝐼 𝑅𝑑𝜙=|𝒙
4𝜋⃗ − 𝒙′ |34𝜋 (𝑅 2 + 𝑍 2 )(2.34)

これを一周積分すると、xy 方向のベクトルは打ち消され、z 方向のみになるため、三角形の相似より z 軸方向の微小磁束密度は=𝑅

𝜇0 𝐼 𝑅𝑑𝜙√𝑅 2 + 𝑍 2 4𝜋 (𝑅 2 + 𝑍 2 )(2.35)となり積分すると
2𝜋

𝑩(𝒙) = ∫0𝑅

𝜇0 𝐼 𝑅𝑑𝜙√𝑅 2 + 𝑍 2 4𝜋 (𝑅 2 + 𝑍 2 )

𝜇0 𝐼
𝑅2
=
2𝜋 (𝑅 2 + 𝑍 2 )32となる。

図 2.3 円電流の軸上の磁場(18)15(2.36)

次にコイル軸上以外の磁場を求める。図 2.4 のように単位ベクトルを用いて(2.37)

𝒙 = 𝑥𝒆𝒙 + 𝑧𝒆𝒛
また、𝒙′ は角度𝜙を用いて
𝒙′ = 𝑅 cos 𝜙 𝒆𝒙 + 𝑅 sin 𝜙 𝒆𝒚 = 𝑅𝒆𝝆(2.38)である。ここで𝒆𝝆は円筒座標で z から遠ざかる距離である ρ が増加する方向を向いた単位ベクトルである。𝑑𝒙′ は式(2.38)を微分して
𝑑𝒙′ = 𝑅𝑑𝜙(−sin 𝜙 𝒆𝒙 + cos 𝜙 𝒆𝒚 ) = 𝑅𝑑𝜙𝒆𝝓(2.39)となる。ここで𝒆𝝓 は円筒座標の角度方向を向いた単位ベクトルである。
また、𝒙 − 𝒙′ は
𝒙 − 𝒙′ = 𝑥𝒆𝒙 + 𝑧𝒆𝒛 − 𝑅(cos 𝜙 𝒆𝒙 + sin 𝜙 𝒆𝒚 )(2.40)

式(2.40)の自分自身との内積を取ることにより|𝒙 − 𝒙′ | = √𝑧 2 + 𝑅 2 + 𝑥 2 − 2𝑅𝑥 cos 𝜙(2.41)

が得られる。
以上を組み合わせて𝑩(𝒙) ==

𝜇0 𝐼 𝑑𝒙′ × (𝒙 − 𝒙′ )∫|𝒙 − 𝒙′ |34𝜋𝜇0 𝐼 𝑅𝑑𝜙𝒆𝝓 × (𝑥𝒆𝒙 + 𝑧𝒆𝒛 − 𝑅𝒆𝝆 )∫4𝜋 (𝑧 2 + 𝑅 2 + 𝑥 2 − 2𝑅𝑥 cos 𝜙)32

(2.42)

ここでの外積は𝒆𝝓 × 𝒆𝒛 = 𝒆𝝆𝒆𝝓 × 𝒆𝝆 = −𝒆𝒛(2.43)

𝒆𝝓 × 𝒆𝒙 = (−sin 𝜙 𝒆𝒙 + cos 𝜙 𝒆𝒚 ) × 𝒆𝒙 = −cos 𝜙 𝒆𝒛

であるので𝑩(𝒙) =

𝜇0 𝐼 𝑅𝑑𝜙(−𝑥 cos 𝜙 𝒆𝒛 + 𝑧𝒆𝝆 + 𝑅𝒆𝒛 )∫34𝜋(𝑧 2 + 𝑅 2 + 𝑥 2 − 2𝑅𝑥 cos 𝜙 )2(2.44)

計算の簡略化のために、𝑧 2 + 𝑅 2 + 𝑥 2 ≫ 2𝑅𝑥 cos 𝜙として𝑩(𝒙) =

𝜇0 𝐼 𝑅𝑑𝜙(−𝑥 cos 𝜙 𝒆𝒛 + 𝑧𝒆𝝆 + 𝑅𝒆𝒛 )∫34𝜋(𝑧 2 + 𝑅 2 + 𝑥 2 )2として計算している。

16(2.45)

図 2.4 円電流の軸上以外での磁場(18)
図 2.5 が実際の計算での磁場形状である。矢印は磁場ベクトルを表し、矢印の大きさは磁場の絶対値の 10-5 倍に比例している。また、メッシュは不均一でありターゲット付近を特に細かく取っているため、塗りつぶしたように 見えている。

図 2.5 初期磁場の形状

2.8 境界条件
境界条件は磁場のみ境界で固定している。その他の物理量は境界での勾配を 0 とし、自由境界としている。
17

第3章

計算手法

3.1 CIP 法
3.1.1 CIP 法の説明
流体方程式中の移流項の数値解析は難しいため、CIP 法(19)を用いる。代表的な例として、速度 u（一定）の 1 次元の波の伝播を考える。その方程式は
∂𝑓
𝜕𝑓
+𝑢
=0
∂𝑡
𝜕𝑥

(3.1)

であり、これを移流方程式と呼ぶ。これを 1 次の風上差分で解くと、図 3.1 のような問題点が生じる。(a)の実線が初期の波形であり、連続体では破線のように移動するだけである。このとき、次の時刻でのメッシュ上の値は白丸となる。しかし、(b)のようにメッ
シュの上の値のみで次のステップに進むと、コンピュータは(c)のように解釈してしまう。
このように波形が変わったまま、次の初期値として移流させると(d)、さらに波形が変わってしまう(e)。これはメッシュ間の波形を線形補間したためである。

図 3.1 移流計算における問題点(19)

この移流計算の問題を解決する手法として CIP 法が提案された。式 3.1 を微分すると、
∂𝑔
𝜕𝑔
𝜕𝑢
+𝑢
=− 𝑔
∂𝑡
𝜕𝑥
𝜕𝑥

(3.2)
18

ただし、𝑔 = ∂𝑓/ ∂𝑥 と定義した。伝播速度 u が一定（∂𝑢/ ∂𝑥 = 0 ）の場合、式(3.2)は式(3.1)と一致し、微分 𝑔 が u の速度で伝播することを表す。図 3.2 のように CIP 法ではメッシュ上の値 f だけではなく、その微分（傾き）も速度 u で伝播させる。そうする
ことで、移動後の波形に図中の矢印のような制限が加わる。よって、傾きの情報を加えることで、移動後の波形が移動前の波形と非常に近くなる。

図 3.2 CIP 法の原理(19)

その具体的方法として、2 つのメッシュ i、i-1 の間のプロファイルを𝐹𝑖 (𝑥 ) = 𝑎𝑖 (𝑥 − 𝑥 𝑖 )3 + 𝑏𝑖 (𝑥 − 𝑥 𝑖 )3 + 𝑐𝑖 (𝑥 − 𝑥 𝑖 ) + 𝑑𝑖(3.3)

のように 3 次多項式で表せば、隣り合う 2 つの格子点上で与えられた 4 つの情報 𝑓𝑖 、
𝑔𝑖 、𝑓𝑖−1 、𝑔𝑖 −1 から以下の関係式を用いて 4 つの未知数𝑎、𝑏、𝑐、𝑑 を決定することができる。
𝐹𝑖 (𝑥 𝑖 ) = 𝑑𝑖 = 𝑓𝑖𝑑𝐹𝑖 (𝑥 𝑖 )𝑑𝑥 = 𝑐𝑖 = 𝑔𝑖
𝐹𝑖 (𝑥 𝑖−1 ) = −𝑎𝑖 ∆𝑥 3 + 𝑏𝑖 ∆𝑥 2 − 𝑐𝑖 ∆𝑥 + 𝑑𝑖 = 𝑓𝑖−1
𝑑𝐹𝑖 (𝑥 𝑖−1 )= 3𝑎𝑖 ∆𝑥 2 − 2𝑏𝑖 ∆𝑥 + 𝑐𝑖 = 𝑔𝑖 −1
𝑑𝑥

(3.4)

格子点間[i-1,1]の内部のプロファイルを表す関数𝐹𝑖 (𝑥 ) が端点でこれらの値とその微分値を取るとすると、式(3.4)より各係数が求められ
𝑔𝑖 + 𝑔𝑖−1 2(𝑓𝑖 − 𝑓𝑖−1 )(∆𝑥 )2 −(∆𝑥)33(𝑓 − 𝑓 ) 2𝑔 + 𝑔
𝑏𝑖 = − (𝑖−1 )2 𝑖 + 𝑖 ∆𝑥 𝑖−1
∆𝑥
𝑎𝑖 =

(3.5)

となる。こうして次の時刻 n+1 での値はこのプロファイルを𝑢∆t だけ移流させたもの、
19

すなわち
𝑓𝑖𝑛+1 = 𝐹𝑖 (𝑥 𝑖 − 𝑢∆𝑡)
= 𝑎𝑖 𝜉3 + 𝑏𝑖 𝜉2 + 𝑔𝑖𝑛 𝜉 + 𝑓𝑖𝑛
𝑑𝐹 (𝑥 − 𝑢Δ𝑡)
𝑔𝑖𝑛+1 = 𝑖 𝑖𝑑𝑥= 3𝑎𝑖 𝜉2 + 2𝑏𝑖 𝜉 + 𝑔𝑖𝑛(3.6)となる。式(3.6)を逐次計算することで、移流方程式の解が求まる。

3.1.2 PINOCO での表記
PINOCO は 2 次元であるため、半径方向と角度方向の 2 速度を扱う。よって微分の式は∂𝒇 ∂𝒇𝑢 ∂𝒇𝑣
+
+
=𝑔
∂𝑡
∂x
∂y

(3.7)

となる。全微分系で表すと
𝐷𝒇 ∂𝒇
∂𝒇
∂𝒇
∂𝑢
∂𝑣
≡
+𝑢
+𝑣
= 𝑔 −𝒇
−𝒇
𝐷𝑡 ∂𝑡
∂x
∂y
∂x
∂y

(3.8)

ここで、式(3.8)の右辺を
𝑔−𝒇

∂𝑢
∂𝑣
−𝒇
= 𝑔∗
∂x
∂y

(3.9)

とする。ここでの 𝑔∗ は反復計算における中間子である。式(3.9)の第 2 項と第 3 項は数値誤差を抑えるために置いている。係数は𝑎𝑖 = 𝒇𝑖𝑏𝑖 = 𝒇′𝑖
3𝑆 − (2𝑓𝑖′ + 𝑓𝑖′+1 )𝑐𝑖 = 𝑖∆𝑖(3.10)としている。また、一般曲線座標系での速度を U、V とすると、
𝑈
𝜉[ ] = [𝜂𝑥𝑉𝑥𝜉𝑦 𝑢𝜂𝑦 ] [𝑣](3.11)と表すことができ、変換を行うと
𝑥𝜉𝑢[ ] = [𝑦𝜉𝑣𝑥𝜂 𝑈𝑦𝜂 ] [𝑉 ](3.12)とできるので、式(3.8)、式(3.9)より
𝐷𝒇 ∂𝒇
∂𝒇
∂𝒇
≡
+𝑈
+𝑉
= 𝑔∗
𝐷𝑡 ∂𝑡
∂ξ
∂η
と表せる。したがって 𝒇 を成分に分けると
20(3.13)

𝐷𝒇𝝃

≡

𝐷𝑡

𝐷𝒇𝜼
𝐷𝑡

≡

∂𝒇𝝃
∂𝑡

∂𝒇𝜼
∂𝑡

+𝑈

+𝑈

∂𝒇𝝃
∂ξ

∂𝒇𝜼
∂ξ

+𝑉

+𝑉

∂𝒇𝝃
∂η

∂𝒇𝜼
∂η

= 𝑔𝜉 ∗(3.14)= 𝑔𝜂 ∗(3.15)となる。

3.2 ALE 法
3.2.1 ALE 法の概要
ALE（Arbitrary Lagrangian-Eulerian）法とは、流体の運動方程式の計算においてオ
イラー法とラグランジュ法の中間に位置する解析法である。オイラー法とは、座標の値に注目し、各時刻で流れによって座標の値がどれほど変化したのか、言い換えればメッ
シュにどれ程の流量が起きているかを見ることで流体の運動を計算する手法である。一方、ラグランジュ法とは流体を細かくメッシュで切り分け、それぞれの動き（メッシュの位置変化）を追うことで、流体の運動を計算する手法である。オイラー法では流れを見るのに適し、ラグランジュ法は物体間の境界を追うのに適している。
それに対し ALE 法では、流体の運動を見る点が空間（メッシュ）に固定されたオイ
ラー法の立場と、流体の運動とは独立 に物体との境界を任意に設定されるラグランジ
ュ法の立場の双方を混在させる(20)。

3.2.2 PINOCO における ALE
レーザー核融合において Rayleigh-Taylor 不安定性は重要な課題の一つである。これは本研究で取り扱うレーザーアブレーションプラズマ生成過程でも重要である。不安定性は初期条件の微小な擾乱が指数関数的に成長するため，ターゲットにおける初期の微小擾乱を追跡できる程度の細かい計算格子が必要となる。しかし、計算領域全てにそのように細かく格子を張ることは、実行に必要な計算機の容量、解析時間が膨大になり非現実的である。そこで，必要な物理モデルをできる限り考慮しつつ、実行する上で適切なモデルで近似することが求められる。1 次元ではラグランジュ法を使うことが有効であるが，2 次元シミュレーションでは流体力学的不安定性などによって格子に歪みが生じ計算継続が困難になることから，PINOCO では CIP 法に格子の移動速度を考慮したALE 法を用いている(13)。
21次に ALE 表記の移流項を含んだ計算を説明する。上記の CIP 法によって求める移流速度を利用する。
静止座標系における反変速度を𝑈0 、𝑉0 とする
𝑈0 = 𝜉𝑥 𝑢 + 𝜉𝑦 𝑣
𝑉0 = 𝜂𝑥 𝑢 + 𝜂𝑦 𝑣

(3.16)

また ALE による計算格子の移動速度を𝑈𝐺 、𝑉𝐺 とする
𝑈𝐺 = 𝜉𝑡
𝑉𝐺 = 𝜂𝑡

(3.17)

式(3.16)、式(3.17)より移動座標系における反変速度𝑈、𝑉 は𝑈 = 𝜉𝑡 + 𝜉𝑥 𝑢 + 𝜉𝑦 𝑣 = 𝑈𝐺 + 𝑈0
𝑉 = 𝜂𝑡 + 𝜂𝑥 𝑢 + 𝜂𝑦 𝑣 = 𝑉𝐺 + 𝑉0

(3.18)

𝜉𝑡 = −𝑥 𝑡 𝜉𝑥 − 𝑦𝑡 𝜉𝑦𝜂𝑡 = −𝑥 𝑡 𝜂𝑥 − 𝑦𝑡 𝜂𝑦(3.19)

ただし、であり、変換のメトリックス
𝑔11 = 𝜉𝑥2 + 𝜉𝑦2
𝑔12 = 𝜉𝑥 𝜂𝑥 + 𝜉𝑦 𝜂𝑦𝑔21 = 𝜂𝑥 𝜉𝑥 + 𝜂𝑦 𝜉𝑦𝑔22 = 𝜂2𝑥 + 𝜂2𝑦(3.20)とする。
基礎方程式（質量保存式、運動方程式）は以下のようになる。
𝜕𝜌
𝜕𝜌
𝜕𝜌
+𝑢
+𝑣
=0
𝜕𝑡
𝜕𝑥
𝜕𝑦

(3.21)

𝜕𝑢
𝜕𝑢
𝜕𝑢
1 𝜕𝑝
+𝑢
+𝑣
=−
𝜕𝑡
𝜕𝑥
𝜕𝑦
𝜌 𝜕𝑥

(3.22)

𝜕𝑣
𝜕𝑣
𝜕𝑣
1 𝜕𝑝
+𝑢
+𝑣
=−
𝜕𝑡
𝜕𝑥
𝜕𝑦
𝜌 𝜕𝑦

(3.23)

式(3.21)を書き換えると𝜕𝜌𝜕𝜉 𝜕𝜌 𝜕𝜂 𝜕𝜌𝜕𝜉 𝜕𝜌 𝜕𝜂 𝜕𝜌
+ 𝑢(
+
) +𝑣(
+
)=0
𝜕𝑡
𝜕𝑥 𝜕𝜉 𝜕𝑥 𝜕𝜂𝜕𝑦 𝜕𝜉 𝜕𝑦 𝜕𝜂

(3.24)

式(3.18)の U、V を使ってまとめると
𝜌𝑡 + 𝑈𝜌𝜉 + 𝑉𝜌𝜂 = 𝜉𝑡 𝜌𝜉 + 𝜂𝑡 𝜌𝜂

(3.25)∂𝑢 𝜕𝜏𝜕𝜉𝜕𝜂≡ 𝑢𝜏 + 𝑢𝜉 +
𝑢 ⋯
∂t 𝜕𝑡
𝜕𝑡
𝜕𝑡 𝜂

(3.26)

また22

を利用して、式(3.22)、式(3.23)を書き換えると𝜕𝑢𝜕𝜉 𝜕𝑢 𝜕𝜂 𝜕𝑢𝜕𝜉 𝜕𝑢 𝜕𝜂 𝜕𝑢
1 𝜕𝑝
+𝑢(
+
)+𝑣(
+
)=−
𝜕𝑡
𝜕𝑥 𝜕𝜉 𝜕𝑥 𝜕𝜂𝜕𝑦 𝜕𝜉 𝜕𝑦 𝜕𝜂𝜌 𝜕𝑥(3.27)𝜕𝑢𝜕𝜉 𝜕𝑣 𝜕𝜂 𝜕𝑣
𝜕𝜉 𝜕𝑣 𝜕𝜂 𝜕𝑣
1 𝜕𝑝
+𝑢(
+
)+𝑣(
+
)=−
𝜕𝑡
𝜕𝑥 𝜕𝜉 𝜕𝑥 𝜕𝜂𝜕𝑦 𝜕𝜉 𝜕𝑦 𝜕𝜂𝜌 𝜕𝑦(3.28)

式(3.24)と同様にして
𝑢 𝑡 + 𝑈𝑢 𝜉 + 𝑉𝑢 𝜂 = −

1 𝜕𝑝
𝜌 𝜕𝑥

(3.29)

𝑣𝑡 + 𝑈𝑣𝜉 + 𝑉𝑣𝜂 = −

1 𝜕𝑝
𝜌 𝜕𝑦

(3.30)

式(3.29)× 𝜉𝑥 + 式(3.30)× 𝜉𝑦 とすると
𝜉𝑥 𝑢 𝑡 + 𝜉𝑥 𝑈𝑢 𝜉 + 𝜉𝑥 𝑉𝑢 𝜂 + 𝜉𝑦 𝑣𝑡 + 𝜉𝑦 𝑈𝑣𝜉 + 𝜉𝑦 𝑉𝑣𝜂
=−(3.31)𝜉𝑥𝜉𝑦(𝜉𝑥 𝑝𝜉 + 𝜂𝑥 𝑝𝜂 ) − (𝜉𝑦 𝑝𝜉 + 𝜂𝑦 𝑝𝜂 )𝜌𝜌書き換えると𝜉𝑥 𝑢 𝑡 + 𝜉𝑥 𝑈𝑢 𝜉 + 𝜉𝑥 𝑉𝑢 𝜂 + 𝜉𝑦 𝑣𝑡 + 𝜉𝑦 𝑈𝑣𝜉 + 𝜉𝑦 𝑉𝑣𝜂
(3.32)𝜉𝜉= − 𝑥 (𝜉𝑥 𝑝𝜉 + 𝜂𝑥 𝑝𝜂 ) − 𝑦 (𝜉𝑦 𝑝𝜉 + 𝜂𝑦 𝑝𝜂 )𝜌𝜌
上式は(𝜉𝑥 𝑢)𝑡 − 𝑢 (𝜉𝑥 )𝑡 + 𝑈{(𝜉𝑥 𝑢)𝜉 − 𝑢 (𝜉𝑥 )𝜉 } + 𝑉{(𝜉𝑥 𝑢)𝜂 − 𝑢(𝜉𝑥 )𝜂 }+(𝜉𝑦 𝑣) − 𝑣(𝜉𝑦 ) + 𝑈 {(𝜉𝑦 𝑣) − 𝑣(𝜉𝑦 ) } + 𝑉 {(𝜉𝑦 𝑣) − 𝑣(𝜉𝑦 ) }
𝑡

𝑡

𝜉

𝜉

𝜂

𝜂

(3.33)𝜉𝜉= − 𝜌𝑥 (𝜉𝑥 𝑝𝜉 + 𝜂𝑥 𝑝𝜂 ) − 𝜌𝑦 (𝜉𝑦 𝑝𝜉 + 𝜂𝑦 𝑝𝜂 )とも書けるので、式(3.16)、式(3.18)、式(3.19)を使ってまとめると(𝑈0 )𝑡 + 𝑈 (𝑈0 )𝜉 + 𝑉(𝑈0 )𝜂
1= − 𝜌 (𝑔11 𝑝𝜉 + 𝑔12 𝑝𝜂 ) + 𝑢(𝜉𝑥 )𝑡 + 𝑣(𝜉𝑦 )(3.34)𝑡+𝑈 {𝑢(𝜉𝑥 )𝜉 + 𝑣(𝜉𝑦 ) } + 𝑉 {𝑢(𝜉𝑥 )𝜂 + 𝑣(𝜉𝑦 ) }
𝜉

𝜂

となる。
同様にして、式(3.29)× 𝜂𝑥 + 式(3.30)× 𝜂𝑦 とすると

23

(𝑉0 )𝑡 + 𝑈 (𝑉0 )𝜉 + 𝑉(𝑉0 )𝜂
1= − 𝜌 (𝑔21 𝑝𝜉 + 𝑔22 𝑝𝜂 ) + 𝑢(𝜂𝑥 )𝑡 + 𝑣(𝜂𝑦 )(3.35)𝑡+𝑈 {𝑢 (𝜂𝑥 )𝜉 + 𝑣(𝜂𝑦 ) } + 𝑉 {𝑢(𝜂𝑥 )𝜂 + 𝑣(𝜂𝑦 ) }𝜉𝜂

したがって、既知の静止座標系の反変速度、式(3.16）を代入し、式(3.34)、式(3.35)を連立して解くことで、未知の移動座標系の反変速度𝑈、𝑉を求める。

24

第 4 章 計算結果
4.1 100 J×6 ビーム
背景プラズマの質量密度を 1.0×10-8 g/cm3 にして計算をした。なお電子数密度で換算すると 1.0×1015 個/cm3 に値する。レーザーの条件は 100 J ビームの 6 本分の 600 J がターゲットに均等に照射されると仮定している。レーザー波長は 1053 nm、パルス幅は1.3 ns とし、先行研究で行われた実験と等しくしている。図 4.1、4.2、4.3 はそれぞれ
t=1.0，5.0，10 ns における電子数密度と磁場ベクトルを重ねた図（上）と x=0 における y 軸上での電子数密度(ne )とイオン温度(Ti)、電子温度(Te )の空間分布（下）である。
まず電子数密度に注目すると、5.0 ns において ne =1017 個/cm3 のプラズマが 0.5 cm
程まで膨張しているのがわかる。ここで 1017 cm-3 の電子数密度は、先行研究で行われた干渉計測によって見積もられた電子数密度の下限値に相当する。またその先に急峻な密度変化があるが、これは背景に置いた低密度プラズマを伝播する衝撃波とアブレーショ
ンプラズマの先端（アブレーションフロント）の重なりである。このうち、衝撃波は本来真空である領域にも背景プラズマを置いているため生じるものであり、実験では存在しない。ターゲット中心位置からアブレーションフロントに向かうにつれて数密度は小さくなっており、これは先行研究での数値解析、即ち Star1D を用いた輻射流体計算の結果と同じ傾向である。
次にイオン温度・電子温度に注目すると、1.0 ns では電子温度の方がイオン温度より高くなっている。これはレーザーパルス幅が 1.3 ns であるため、レーザーが入射されていることと、電子の方がイオンに比べて 2000 倍以上軽いためレーザーエネルギーを吸収しやすく、電子温度の方が上昇しやすいことによるものである。また 5.0 ns、10.0 ns
ではアブレーションフロント付近においてイオン温度が電子温度より高くなって いる。
レーザーエネルギーを吸収して膨張するアブレーションプラズマが背景プラズマ に衝突すると、背景プラズマ中には順行衝撃波が、アブレーションプラズマには逆行衝撃波が励起されると考えられる。衝撃波の下流では、上流の運動エネルギーが熱エネルギーへと変わるため、アブレーションプラズマの先端では逆行衝撃波によって、背景プラズ
マは順行衝撃波によって加熱される。そのためイオン温度が高くなっている。通常プラ
ズマではイオンと電子の衝突により、イオンと電子はほぼ等温となるが、背景領域のプラズマ密度が低いのとアブレーションプラズマの先端の密度が低いため、衝突による温度緩和が起こりにくく、イオン温度が高いままになっている。一方、アブレーションプ
ラズマの内側では密度が高いため温度緩和が生じ、イオンと電子がほぼ等温にとなっている。
25

図 4.1 レーザー照射から 1.0 ns 後における電子数密度 ne 個/cm3 と磁場ベク
トル（上）y 軸上の電子数密度 ne 個/cm3(黒)、イオン温度 Ti keV(赤)、電子温度 Te keV(青)の空間分布（下）26図 4.2 レーザー照射から 5.0 ns における電子数密度 ne と磁場ベクトル（上）x=0 における y 軸上における電子数密度 ne 、イオン温度 Ti、電子温度 Te の空間
分布（下）

27

図 4.3 レーザー照射から 10.0 ns 後における電子数密度 ne と磁場ベクトル（上）y 軸上における 電子数密度 ne 、イオン温度 Ti、電子温度 Te の空間分布（下）また、磁場の絶対値の空間分布を見ると(図 4.4.1、4.4.2、4.4.3)、プラズマの膨張によりアブレーションフロントで磁場が圧縮され値が大きくなっている。10.0 ns において、
x=0 での y 軸の空間分布を見ると(図 4.4.4)、アブレーションフロントにおいて磁場強度が初期値 1 T の 4 倍程になっており、磁場が圧縮されているのが確認できる。また、アブレーションプラズマ内の磁場はほぼ 0 であり、反磁性キャビティができている。
28

図 4.4.1

t=1.0 ns における磁場空間分布

図 4.4.2

t=5.0 ns における磁場空間分布

図 4.4.3

t=10.0 ns における磁場空間分布29

図 4.4.4 電子密度（実線）と磁場（点線）の y 軸上空間分布

4.2 レーザーエネルギー依存性
次に入射レーザーエネルギーを 100 J×6、50 J×6、10 J×6 ビームと変化させてシミュレーションを行い、レーザーエネルギーによる依存性を検証した。

4.2.1

100 J

図 4.5.1-4.5.6 が 4.1 節の計算を継続したものである。t=10 ns では、図 4.5.1 のように衝撃波と共にアブレーションプラズマが膨張し、それにより外部磁場の形状が変化するのが確認できる。しかし、t=30 ns では、図 4.5.2 のようにアブレーションプラズマが計算領域の外周まで到達し、外部磁場は領域外へ掃き出されている。境界は自由境界にしているため、領域内で減速しなかったプラズマは領域外に抜けていく。t=50 ns では
外周部の磁場が強くなり始め(図 4.5.3)、時間がさらに経過し、t=80 ns ではその磁場が周囲のプラズマを帯状に押し戻しているのが見えている(図 4.5.5)。しかし、この磁場変化とプラズマの相互作用は目的に沿ったものではない。50 ns から磁場が生じたのは境界条件として外周で初期磁場の値を固定しており、初期磁場で抑えきれなかったプラズ
マが領域外に抜けた後、固定した磁場がアブレーションプラズマ内に拡散したからである。しかし、t=100 ns では磁場拡散が進み、初期磁場に近い状態になる。そして、比較的遅いアブレーションプラズマの膨張が磁場によって押さえつけられている(図 4.5.6)。
t<50 ns において磁場が領域外に出てしまうため実験の再現には至っていないが、t=80
ns 以降、密度の低いアブレーションプラズマに拡散することで再び戻った磁場と、速度が遅いアブレーションプラズマの相互作用は確認できている。
30

図 4.5.1 t=10 ns における ne [個/cm3]と磁場の空間分布

図 4.5.2 t=30 ns における ne と磁場の空間分布

図 4.5.3 t=50 ns における ne と磁場の空間分布31図 4.5.4 t=60 ns における ne と磁場の空間分布

図 4.5.5 t=80 ns における ne と磁場の空間分布

図 4.5.6 t=100 ns における ne と磁場の空間分布

32

4.2.250 J

次にレーザーエネルギーを半分の 50 J×6 ビームにして検証をした結果を図 4.6.14.6.6 に示す。
図 4.6.1 のように 10 ns では図 4.5.1 に比べて膨張プラズマの半径が小さくなってい
るのがわかる。これはプラズマが吸収するレーザーエネルギーが小さくなった効果によるものである。しかし、30 ns では 100J の時と同様に衝撃波と共にアブレーションプ
ラズマが計算領域の外周まで到達している(図 4.6.2)。そして、50 ns では 100 J の時と同様に外周から磁場が入り込んできている(図 4.6.3)。t=60 ns(図 4.6.4)、80 ns(図 4.6.5)と時間がたつにつれ、100 J の時と同様に磁場が帯状にプラズマを押し戻していく。
しかし、その帯状のプラズマの形状が 100 J の時のプラズマと異なっている。レーザ
ーエネルギーを減らした影響で衝撃波とアブレーションプラズマが持つ運動エネ ルギ
ーが小さくなる。運動エネルギーの減少により、プラズマの膨張の速度が小さく なる。
そのため、計算領域外に抜けることのできる速度を持ったプラズマの量が減り、外周付
近に滞在するアブレーションプラズマの量が増える。したがって 100 J の時程、容易に磁場が拡散できないのである。そして、磁場による押し返しの帯状プラズマがy ≤ 0.5 cm
では、より中心に近く位置するのは初期磁場形状によるものである。図 2.2 で示したよ
うにコイルは計算領域外にあり、その初期形状は図 2.5 の通りである。そのため、x ≤0 cmにおける外周上では、磁場の強さは y に反比例する。磁場の強さが大きい程、磁気圧は大きくなるため、図のような形状になっている。
t=100 ns(図 4.6.6）では磁場拡散が進み、100 J の時と同様に磁場によるプラズマの押さえつけが見えている。100 J のとき(図 4.5.6)と比較すると、磁場がプラズマを抑えつける位置が x 軸正の向きに移動している。また、アブレーションプラズマは小さくなっている。これもレーザーエネルギーを小さくした影響である。

図 4.6.1 t=10 ns における ne [個/cm3]と磁場の空間分布33

図 4.6.2 t=30 ns における ne と磁場の空間分布

図 4.6.3 t=50 ns における ne と磁場の空間分布

図 4.6.4 t=60 ns における ne と磁場の空間分布34図 4.6.5 t=80 ns における ne と磁場の空間分布

図 4.6.6 t=100 ns における ne と磁場の空間分布

4.2.3

10 J

次にレーザーエネルギーを 1/10 の 10 J×6 ビームにして検証をした。図 4.7.1-4.7.6 がその結果である。t=10 ns(図 4.7.1)では 50 J の場合(図 4.6.1)に比べてさらに膨張プラズ
マの半径が小さくなっているのがわかる。これはレーザーエネルギーによる依存性を表している。t=30 ns(図 4.7.2)では 100 J、50 J の時とは異なり、衝撃波とアブレーショ
ンプラズマが計算領域の外周まで到達していない。そして外部磁場がまだ計算領域内に残っている。t=50 ns(図 4.7.3)ではコイル側の磁場は存在しているが、x ≥ 0の領域ではプラズマが計算領域外に抜け、それに押し出される形で磁場もなくなっている)。t=60 ns、
80 ns（それぞれ図 4.7.4、図 4.7.5）と時間がたつにつれ、50 J の時と同様に磁場が帯状にプラズマを押し戻していく。しかし、その帯状のプラズマの形状は 50 J の時のプラズマと少し異なっている。これは 50 J の時よりも更に領域内に留まるプラズマの量が増えたためである。t=100 ns(図 4.7.6) ではこれまでと同様に磁場の押さえつけが見えており、境界の位置、プラズマの大きさはさらに小さくなっている。
35

図 4.7.1 t=10 ns における ne [個/cm3]と磁場の空間分布

図 4.7.2 t=30 ns における ne と磁場の空間分布

図 4.7.3 t=50 ns における ne と磁場の空間分布36図 4.7.4 t=60 ns における ne と磁場の空間分布

図 4.7.5 t=80 ns における ne と磁場の空間分布

図 4.7.6 t=100 ns における ne と磁場の空間分布374.3 磁気圧の効果の確認
次に、計算領域を半径 3.2 cm まで広げ、コイル周りの磁場を計算領域内に入れた。また磁場拡散項を無視することで磁気圧による効果を調べた。
コイル付近では磁場が強いため、磁気圧、つまり式(2.32)に示す右辺第 1 項の影響が大きく、背景に置いた静止しているプラズマに流れが生まれ、思惑と異なった磁場の移流が生じていることが分かった。図 4.8.1 と図 4.8.2 を見比べてわかるように、外部磁場の形状がコイル位置から次第に変化し、アブレーションフロント（衝撃波）が到達するまで保たれていない。しかし、磁気圧がプラズマに及ぼす効果が表れていることは確認できた。

図 4.8.1 0.1 ns における ne と磁場の空間分布(100 J×6 ビーム)

図 4.8.2 5 ns における ne と磁場の空間分布
38そこで磁気圧が働く条件を変更した。電子数密度でしきい値を 1017 cm-3 と設定し、その値未満の領域では磁気圧を働かせないようにした。このしきい値は実験結果で計測できた電子数密度を参考にした。その結果図 4.9.1、図 4.9.2 の通りコイル付近での初期磁場を保つことに成功した。しかし、計算時間が進むにつれ、アブレーションプラズマの膨張が進み、コイル位置まで達するとコイル付近の磁場の形状が変わってしまった(図4.9.3)。実験ではコイルに電流を流す時間、言い換えると磁場継続時間が 200 µs である
ため、コイル付近の磁場が変化することは起こりえない。また、コイル付近の磁場は元々強く、アブレーションプラズマの膨張による圧縮で更に強くなるが、アブレーションフ
ロントを抑えていない。これはアブレーションプラズマの膨張と衝撃波により圧縮された磁場ごと流されているとわかった。

図 4.9.1 しきい値を設けた結果(t=0.1 ns)

図 4.9.2 しきい値を設けた結果(t=5 ns)

39

図 4.9.3 しきい値を設けた結果(t=40 ns)

40

第 5 章 考察
5.1 実験結果との比較
図 5.1(a)に 100 J×6 ビームの条件におけるシミュレーション結果を、図 5.1(b)-図5.1(c)には実験から得られた干渉計測画像とそれを解析した電子密度の x-y 分布を示す。
図 5.1(a)では初期ターゲット位置において非常に高い密度となっているが、これは計算の簡略化のためターゲットを中実球としたためである。また、図 5.1(c)では電子がx ≤ 0
側に偏って分布している。しかし、図 5.1(a)ではこの傾向はない。相違が生じた原因は背景プラズマ中の磁場が移流したため、実験における磁場が再現できていないことにある。磁場の移流が改善できれば、減速されるプラズマが多くなり、実験の再現に近づくと考えられる。

(a)

(b)
(c)
図 5.1 100J×6 ビーム, レーザー照射から 50 ns 後における結果(a)PINOCO,(b)干渉計測画像(12),(c)電子数密度の解析結果(12)41図 5.2(a)に 50 J×6 ビームの条件におけるシミュレーション結果を、図 5.2(b)-図 5.2(c)には実験から得られた干渉計測画像とそれを解析した電子密度の x-y 分布を示す。実験
結果では 100 J のときと同様にプラズマが偏って分布している。しかし PINOCO では
その特徴は見られない。これもアブレーションプラズマの一部が計算領域から抜けているため、密度を過小評価しているからと考えられる。

(a)

(b)
(c)
図 5.2 50J×6 ビーム, レーザー照射から 50 ns 後における結果(a)PINOCO,(b)干渉計測画像(12),(c)電子数密度の解析結果(12)42図 5.3(a)に 10 J×6 ビームの条件におけるシミュレーション結果を、図 5.3(b)-図 5.3(c)には実験から得られた干渉計測画像とそれを解析した電子密度の x-y 分布を示す。(a)、(c)のどちらにも x=2 mm 付近で 1017 cm-3 の電子数密度が見えている。しかし、(a)のプラズマの方が大きく広がっている。この違いが生じた原因として、PINOCO では衝撃波の伝播により磁場が移流したため、実験よりプラズマが大きく広がっているということが考えられる。

(a)

(c)

(b)

図 5.3 10J×6 ビーム, レーザー照射から 50 ns 後における結果(a)PINOCO,(b)干渉計測画像(12),(c)電子数密度の解析結果(12)43図 5.4(a)に 10 J×6 ビームの条件でレーザー照射から 100 ns のシミュレーション結果を、図 5.4(b)-図 5.4(c)には実験から得られた干渉計測画像とそれを解析した電子密度の x-y 分布を示す。(a)、(c)のどちらにも x=－3 mm 付近で電子数密度 1017 cm-3 程度のプラズマシェル構造が見えている。PINOCO において磁気圧の効果が働き、磁場とプラズマの相互作用が起きていることが確認できる。しかし、この磁気圧は再び拡散した磁場によるものであり、初期磁場によるものとは異なる。その原因である、衝撃波の伝播
で背景プラズマ中の磁場が移流する問題を改善することで、実験の再現に近づくと考え
られる。

(a)

(c)

(b)

図 5.4 10J×6 ビーム, レーザー照射から 100 ns 後における結果(a)PINOCO,(b)干渉計測画像(12),(c)電子数密度の解析結果(12)445.2 計算で生じた問題点のまとめ
4.1 節で述べたイオン温度が電子温度より上がる点について、改良すべきであると考える。先にも述べたがレーザーエネルギー吸収は電子に影響するので、レーザー入射時に電子温度が高く、イオン温度の方が低くなるべきであり、その点はよく再現されている。しかし、低密度・低温の背景プラズマを置くことによって生じるイオン温度の加熱
は、実験では起きない現象であり、実験の再現性に影響し得る。そのため、改良が必要とされる。
4.2 節で述べた衝撃波とアブレーションプラズマの膨張により、磁場が移流し領域外に流されてしまうのは問題である。また、境界に磁場を固定したことで、後の磁場がプラズマを抑え込むのが確認できたが、実験結果が示すようにプラズマの膨張によって圧縮された磁場による押し戻しが起きるべきである。しかし、それが再現できていない理
由として磁気圧によるアブレーションプラズマの抑え込みができていないのだと 考えられた。
4.3 節ではその 4.2 節での問題解決に向け、磁場拡散項を無視して磁気圧の影響を検証したところ、磁気圧の効果は確認できた。しかしアブレーションプラズマの減速は不十分であり、磁場拡散と磁気圧の影響ではなく、移流によって磁場は外側に流され、ア
ブレーションプラズマの膨張により圧縮されつつも、留まっていないことが原因とわかった。今後、実験の再現に向けて背景プラズマ内の磁場の移流の取り扱いにおける改良が必要である。その改良として、アブレーションプラズマと背景プラズマの境界を追跡し、背景プラズマ内での磁場の移流を調整する。

45

第6章

結論・今後の展望

本研究では、従来の磁気スラストチャンバーの模擬実験に対する再現の数値解析結果に相違が生じた問題において、レーザープラズマ生成過程の計算において磁場変化を十分に考慮できていないことが原因であると指摘した。また、その問題解決に向け、従来の数値計算の流れの中で輻射流体コード(Star1D、Star2D)を用いていた箇所を、大阪大
学 レ ー ザ ー エ ネ ル ギ ー 学 研 究 セ ン タ ー の 長 友 氏 が 開 発 し た 輻 射 磁気 流 体 コード(PINOCO)を用いて計算することを試みた。
これまでレーザー核融合の爆縮過程の計算用に開発された PINOCO を、我々の実験体系に導入した。その際、レーザー入射において等方照射を仮定し、また実験ではチャ
ンバー内を真空にしているが、本研究の数値計算では背景領域に低温・低密度のプラズマをおいた。これは背景領域を MHD 近似するためである。
計算から、磁気スラストチャンバーの推進原理に通ずる、磁場の圧縮を確認した。しかし、圧縮された磁場がアブレーションプラズマを抑える構造の再現には至っていない。
レーザーエネルギーを変えて計算することで、アブレーションプラズマの半径の変化が確認でき、レーザーエネルギー依存性を示した。計算において磁場が計算領域外には掃きだされる問題が生じた。また、レーザー入射から 50 ns 以降では、境界に固定され磁場が再び背景プラズマ内を中心に向かって拡散した結果、初期磁場に近い磁場が戻り、遅いアブレーションプラズマと磁場との相互作用が見えた。しかし、100 ns では磁場が押さえつける構造ができ、実験の状態に近い結果が得られた。
また、磁気圧の効果を調べたところ、コイル中心から磁気圧による流れが生じている問題が見つかったが、計算上、磁気圧にしきい値を設けることで解決はできた。そして磁場が流される原因が背景においたプラズマ内の移流が大きいことにあるとわか った。
今後のイオン温度の問題、磁場の移流の問題の解決に取り組み、実験の再現の追求を行う。そして、レーザー核融合ロケットの実現に向けた、スケーリング則の構築に役立
つコードとなるよう改良を進めていく。

46

参考文献
(1) 中島秀紀, 原子力ロケットの現状と未来 (2014).
(2) R. A. Hyde, L. L. Wood, and J. H. Nuckolls, AIAA Pap. 72, 1063 (1972).
(3)C. D. Orth, VISTA-A Vehicle for Interplanetary Space Transport Application
Powered by Inertial Confinement Fusion, UCRL-TR-110500 (2003)(4)中島秀紀，梶村好宏，核融合ロケットの開発に向けてレーザー核融合ロケット，

Journal of Plasma and Fusion Research Vol. 83，No. 3 (2007)(5)T. Nagamine et al., "Development of a MHD code Using CIP Method and its
Application to ICF Rocket ," Proceedings of the 12th International Conference on
Magnetohydrodynamic Electrical Power Generation 1, 559 (1996)(6)Konstantin V. VCHIVKOV et al., "Laser-Produced Plasma Experiments andParticle in Cell Simulation to Study Thrust Conversion Process in a Laser Fusion
Rocket ," Japanese Journal of Applied Physics 42, 6590 (2003)(7)花屋倫生，永久磁石を用いた磁気ノズル中のプラズマ挙動に関する実験及び解析，学士論文，九州大学，(2008)(8)前野旭弘, レーザー生成プラズマと磁場の相互作用による運動量変化に関する研究 ,
修士論文, 九州大学大学院, (2010).
(9)前野旭弘, 磁気スラストチャンバーシステムの実験的実証と数値シミュレーションに関する研究,博士論文,九州大学大学院,(2014)(10)富永宙志, 2 次元放射流体コードを用いた磁場中プラズマの挙動解析,修士論文,九州大学大学院（2014）(11)枝本雅史 et al. Development of a magnetic thrust chamber for a laser fusion
rocket，American Institute of Aeronautics and Astronautics, (2016)(12)Ryosuke Kawashima, Measurement of plasma structure in a Magnetic Thrust
Chamber,修士論文,九州大学大学院（2015）(13)長友英夫，城_知至，砂原淳，三間圀興,高速点火核融合の統合シミュレーション－多階層プラズマシミュレーションシステム“FI 3”－,J. Plasma Fusion Res. Vol.82, No.3
(2006)(14)H. Nagatomo et al., Computational study of magnetic field compression by laserdriven implosion, Nucl. Fusion 55 (2015) 093028
(15)

D.M.

Mihalas,andB.

Weibel-Mihalas,

Hydrodynamics",Oxford University Press, (1984).
(16)Robinson A.P.L. et al 2014 Nucl. Fusion 54 054003
47

FundationsofRadiation

(17)Braginskii S.I. 1965 Review of Plasma Physics vol 1 (NewYork: Consultants

Bureau)(18)前野昌弘: “よくわかる電磁気学” (東京図書,2010)(19)矢部孝・内海隆行・尾形陽一:“CIP 法 原子から宇宙までを解くマルチスケール解法”(森北出版,2003)(20)野村卓史・ 飯島政義, ALE 法に基づく粘性流体と構造との相関問題の有限要素解析
手法, 土木学会論文集 第 416 号/1-13 (1990)48

謝辞
本研究を進めるにあたり、研究の場を提供していただき、また、研究に関する数多くの助言をいただきました山本直嗣教授に心から感謝いたします。そして、研究について親身に見てくださり、幅広く丁寧にご指導いただきました森田太智助教に心より感謝いたします。また、研究についてのアドバイスや議論をしてくださった中島秀紀名誉教授に心より感謝いたします。
数値解析を行うにあたり、計算コードの提供、コードの改良についての指導をして下さった、大阪大学レーザーエネルギー学研究センターの長友英夫准教授に深く感謝いた
します。
事務手続きなどお世話になった馬渡隆子氏に感謝いたします。そして、ともに研究やスポーツに励んだ研究室同期生の竹末一平氏に深く感謝いたします。また、公私共に関わった李後毅氏、齋藤直哉氏、牛尾康一氏、市丸智裕氏、飯島健介氏、上野文輔氏、江川雄亮氏。枝本雅史氏、三浦智之氏、Ali Yousefian 氏、板谷佑太朗氏、岩本政隆氏、川原友太郎氏、小山智樹氏、福本湧氏、淵上太貴氏、Huang Yan 氏ら研究室の先輩に感謝いたします。
最後に、22 年間その成長を支え続けてくれた家族に深く感謝いたします。

49