令和元年度
九州大学大学院 総合理工学府
先端エネルギー理工学専攻

修 士 論 文論文名
レーザー核融合ロケットにおける中性子遮蔽を考慮した円錐形推進剤の推進性能評価

氏

名

指導教員名

山村

徹

山本 直嗣

目次
第1章

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

1.1

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

1.2

先行研究 ............................................................................................................................. 4

1.3

目的 ..................................................................................................................................... 5

第2章

シミュレーションコード ...................................................................................................... 6

2.1

概要 ..................................................................................................................................... 6

2.2

二次元輻射流体コード (Star2D) について ..................................................................... 62.2.1

基礎方程式 .................................................................................................................. 7

2.2.2

計算領域，計算格子 .................................................................................................... 8

2.3

3D hybrid code について ................................................................................................ 11

2.3.1

基礎方程式 ................................................................................................................ 11

2.3.2

計算方法 .................................................................................................................... 14

2.4

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

計算体系 ........................................................................................................................... 27

2.5.1

Star2D における計算体系......................................................................................... 27

2.5.2

3D hybrid code における計算体系 ........................................................................... 31

第3章

計算結果 ........................................................................................................................ 32

3.1

評価方法 ........................................................................................................................... 32

3.2

球状ターゲット ................................................................................................................... 32

3.2.1

Star2D の計算結果 ................................................................................................... 32

3.2.2

3D hybrid code の計算結果...................................................................................... 36

3.3

円錐ターゲット ................................................................................................................... 38I

3.3.1

Star2D の計算結果 ................................................................................................... 38

3.3.2

3D hybrid code の計算結果...................................................................................... 42

第4章

考察 ................................................................................................................................ 44

4.1

Star2D ............................................................................................................................... 44

4.2

3D hybrid code ................................................................................................................. 50

第5章

結論 ................................................................................................................................ 53

参考文献 .......................................................................................................................................... 54
謝辞 .................................................................................................................................................. 56

II

第1章
1.1

序論

背景

宇宙開発の次なるステップとして有人惑星間飛行が挙げられる。特に，火星は人類が地球圏外へ進出する際，比較的早い段階に移住できると考えられる唯一の天体である．
火星の大気には二酸化炭素が豊富に存在し，大量の水が凍った状態で眠っているからである．また，火星は人類にとって次の居住区になるばかりでなく，火星での経験が，ほかの惑星の探査・開発にとって重要なものとなり，前哨基地としての役割も担うことになると考えられる．しかし，近い将来，火星に恒久的な基地や居住施設を建設することを考えると，従来の宇宙輸送手段である化学推進ロケットでは，推力は大きいが比推力が低く，火星まで移動するために大量の燃料を必要とするため，建設に必要な膨大な量の資材を運搬することが困難である．また，電気推進によるものは燃料利用効率が高く，長時間使用することができるが，推力が小さいため，目的地に到着するまでの時間が長くなる．
そのため，無重力による人体への影響，宇宙放射線による被曝を考慮してできるだけ短い時間でミッションを終わらせる必要がある有人飛行には向かない．そこで，その際に使用される宇宙輸送手段として考えられているのが，レーザー核融合ロケットである．レー
ザー核融合ロケットの概念設計はローレンスリバモア国立研究所によって提案された (1)．
VISTA(2)のフライトモデルを図 1.1 に示す．

図 1.1

レーザー核融合ロケット VISTA のフライトモデル
1

レーザー核融合ロケットは，核融合による莫大なエネルギーで，大推力・高比推力を同時に達成するため，従来の化学推進ロケットで火星往復 500 日以上かかるところ，レーザ
ー核融合ロケットなら 145 日で帰還可能と推定されている(2)．
核融合によって生じるエネルギーは，プラズマ・輻射・中性子の形で放出され，そのう
ちプラズマは荷電粒子の集合体であるため，磁場によってその運動方向が変化する．レーザー核融合ロケットは，磁気ノズルを用いることで，レーザー核融合により爆発的に膨張する核融合プラズマと外部磁場との相互作用を利用し，固体壁との相互作用なく運動方向を変えて推力を得る．この推進原理は，“磁気スラストチャンバー”と呼ばれており，
図 1.2 に VISTA における磁気スラストチャンバーのメカニズムを示す．超電導磁石によって作られた磁場中でレーザー核融合を起こし，膨張する核融合プラズマは磁場を圧縮す
る．磁場は圧縮されることで強度が増し，プラズマに対して磁場に直角に強い力を及ぼす
ようになる．これによってプラズマは減速され，プラズマから磁場へエネルギーが輸送される．さらに，磁場によって減速されたプラズマはやがて停止し，磁場のエネルギーの方がプラズマのエネルギーよりも大きくなると磁場はプラズマを押し戻す．それに伴って，プラズマは磁場により後方に加速される．このときに，磁場からプラズマへエネルギーが与えられる．このようにして，磁場がプラズマ粒子の運動方向を変化させ，プラズマを押し出す
反作用により推力を得ることができる．

図 1.2

VISTA における磁気スラストチャンバーのメカニズム2

磁気ノズルは過去に設計，シミュレーションが行われてきた．Hyde は，2 次元 MHD シ
ミュレーションを用いてプラズマ挙動と推進効率を計算し，1 つのコイルで推進効率 65%という結果を得ている(3)．長峰と中島は 3D hybrid code を用いてプラズマ挙動と推進効率を計算し，推進効率 65%という結果を得ている(4)．Vchikov は同様のコードで推進効率の最適化を行い，推進効率 70%という結果を得ている(5)．また，実験との比較も行い，推進効率の一致を得ている(6)．
レーザー核融合ロケットの課題の一つとして，核融合により莫大なエネルギーを生み出す反面，同時に高速中性子が生じ，宇宙機方向に飛んでいき，その中性子がコイルや宇宙機を加熱してしまうことが挙げられる．
この中性子からコイルや宇宙機を守るため，密度が大きく厚い中性子遮蔽体が必要となる．また，中性子遮蔽体に落ちた中性子のエネルギーを捨てるための重い排熱装置も必要であり，これらの重量は推進剤を除いた VISTA の全重量の 32%を占めるため，宇宙機の重量増加の原因となる．表 1.1 に VISTA における各システムの重量分布を示す(2)．これらを減量することは，積み荷を増やすことや比出力の向上，すなわち宇宙機の性能向上につながるため，宇宙機への中性子放射を削減することが重要である．

3

表 1.1VISTA における各システムの重量分布

Sub-system

1.2

Mass [ton]Total Mass [ton]Payload System

289

Propellant System

4458

Propellant

4140

Tritium refrigerator

50

Propellant tanks

207

Tritium tank

61

Driver System

268

Thrust Chamber System

805

Coil

216

Coil shield + structure

523

Radiators

66

Auxiliary Systems

180

Total DRY MASS [ton]1860

Total WET MASS [ton]

6000

先行研究

Orth による VISTA の概念設計(2)では，図 1.3(a)に示すように，ターゲットは D-T 燃料を覆うように 50 g の固体水素推進剤が球状に配置された形状となっている．それに対し，核融合により生じる宇宙機方向への中性子放射を削減するために，Şahin(7),(8)により，図1.3(b)のような推進剤形状にした場合の VISTA に対する中性子遮蔽計算が行われた．
この円錐ターゲットは，燃料ペレットに円錐形状の固体水素推進剤が取り付けられ，円錐の中心角は 100 °で VISTA の超電導コイル方向にある．この形状は，固体水素推進剤をコイル方向に集中配置することにより，推進剤中で入射中性子を減衰させる．球状ター
ゲットは，一次元中性子輸送計算コード ANISN(9)を使用し，円錐ターゲットは，r-z 体系の二次元中性子輸送計算コード DORT(10)を用いて， DT 核融合反応後，各形状における4

固体水素推進剤を通過する中性子計算が行われた．ANISN では 30 群の中性子エネル
ギー群，12 群の γ 線エネルギー群に分け，DORT では 5 群の中性子エネルギー群，2 群の γ 線エネルギー群に分けられ計算された．その結果，円錐ターゲットの場合，球状ター
ゲットに比べ，コイル方向への入射中性子束が約 1 桁減衰することがわかっている．また，これにより，中性子遮蔽体及び排熱装置の重量も軽減できることもわかっている．

図 1.3

1.3

ターゲット形状 (a)球状ターゲット，(b)円錐ターゲット

目的

前述の円錐ターゲットは，コイル方向への中性子放射を削減でき，中性子遮蔽体や排熱装置の質量を軽減させるメリットがあるが，推進効率などの性能は明らかでなく，具体的なプラズマ計算や実験は未だ報告されていない．
そこで，本研究の目的は，核融合によるコイル方向への中性子放射を推進剤中での削減が見込める円錐ターゲットの推進性能を，二次元輻射流体コード及び 3D hybrid
code を用いた数値計算により算出することとする．また，円錐ターゲットの推進性能を評価するため，球状ターゲットの推進性能と比較を行う．

5

第2章
2.1

シミュレーションコード

概要

この章では解析に使用した計算コードについて述べる．本研究では大きく二つの計算コードを用いたカップリングシミュレーションを行う．初めに，核融合反応直後を仮定した
高温の核融合プラズマと低温の推進剤との相互作用を二次元輻射流体コード(Star2D)で計算する．次に，カップリングコードにて二次元データを三次元データに拡張し，流体を粒子として変換する．最後に磁場中でのプラズマの挙動を 3D hybrid code で計算し，その結果より球状ターゲットと円錐ターゲットの推進性能を評価する．
本研究で取り扱うターゲットの核融合プラズマは，核融合反応直後を仮定した高温であるため輻射の影響を考慮する必要がある。そのため輻射流体計算によりプラズマの膨張過程を計算する。また，一般にプラズマの運動は，統計力学的な複雑さに加え，その運動がさらに電磁気力に支配されるため大変複雑である．本研究では輻射輸送，状態方程式，電磁気力など多くの物理モデルが重要となるが，すべての物理モデルを考慮した場合計算時間が膨大になる．そのため輻射輸送，状態方程式などの物理モデルを輻
射流体コードで取り扱い，3D hybrid code では電磁気力が関与した物理モデルを取り扱ったカップリングシミュレーションを行う．以降で各計算コードの詳細を説明する．

2.2

二次元輻射流体コード (Star2D) について

二次元輻射流体コード(Star2D)(11),(12)は，レーザー核融合のシミュレーションを行う爆縮コードとして，砂原 敦 氏が開発・提供しているものである．「アブレーション，衝撃波，
レーザー光線追跡，レーザー吸収，電子加熱，熱伝導，状態方程式，原子過程」などの物理モデルを取り扱っている．本研究では，核融合反応後，高温プラズマ状態にある核融合燃料が膨張し，推進剤を加熱して膨張する過程を計算するために使用させて頂いた．注意する点として，2 次元輻射流体コードでは実際の 3 次元での現象を 2 次元で計算シミュレーションするので，解析しようとする対象には物理量に対称性があることが必要条件となる．今回のシミュレーションの場合，ターゲット形状・ターゲット密度・磁場に同じ軸の対称性が仮定できるように条件設定した．出力された 2 次元データの対称性を仮定
6

した，3 次元への拡張は「2.4 カップリングコード」で説明する．

2.2.1

基礎方程式

以下の方程式は直交座標系および，CGS 単位系で記述されている．イオン，電子の 2温度を仮定した輻射流体の方程式を以下に示す．

連続の式

運動方程式

𝜌

𝑑𝜌
= −𝜌∇ ∙ 𝒖
𝑑𝑡

(2.1)

𝑑𝒖
= −∇(𝑝i + 𝑝e + 𝑞)
𝑑𝑡

(2.2)

𝑑𝑇i
= −𝑝i ∇ ∙ 𝒖 − ∇ ∙ (𝜅i ∇𝑇i ) + 𝛼(𝑇e − 𝑇i ) + 𝑄i
𝑑𝑡

(2.3)

𝑑𝑇e
= −𝑝e ∇ ∙ 𝒖 − ∇ ∙ (𝜅e ∇𝑇e ) − 𝛼(𝑇i − 𝑇e ) + 𝑄e
𝑑𝑡

(2.4)

イオンエネルギー式

𝜌𝑐𝑣i

電子エネルギー式

𝜌𝑐𝑣e

ここで，𝜌 は密度，𝒖 は流体速度， 𝑝i はイオン圧力，𝑝e は電子圧力，𝑞 は人工粘性，𝑇i はイオン温度，𝑇e は電子温度， 𝑐𝑣i はイオン比熱，𝑐𝑣e は電子比熱，𝜅i ・𝜅e はそれぞれイオ
ン・電子の Spitzer の熱伝導係数(13)，𝛼 は Spitzer の熱緩和係数(14)，𝑄i は“人工粘性に伴う運動エネルギーからイオンの熱エネルギーへの変換分”，𝑄e は“レーザー吸収とプラズ
マ中の輻射吸収により電子が受け取るエネルギー”である．(※添え字，はそれぞれイオ
ン，電子の物理量であることを示す．)

次に輻射輸送について記す．熱束制限多群拡散近似(15)により光子のエネルギーを多数のエネルギー群に分け，それぞれのエネルギーについて拡散近似(16)により解いた．式(2.5)に輻射輸送方程式を示す．

𝜌𝑑 𝐸𝑣( ) − ∇ ∙ (𝐷𝑣 ∇𝐸𝑣 ) = 4𝜋𝜂𝑣 − 𝑐𝜒𝑣 𝐸𝑣
𝑑𝑡 𝜌
7

(2.5)

ここで，𝐸𝑣 は周波数𝑣に対するエネルギー密度，𝐷𝑣 は拡散係数，𝜒𝑣 は吸収率，𝜂𝑣 は輻
射率，𝑐 は光速である．式(2.5)の右辺は流体からの輻射によるエネルギー損失を表しており，積分した値は電子エネルギー式の生成項𝑄e に一部として代入される．

2.2.2

計算領域，計算格子

Star2D の計算領域の概略を図 2.1 に示す．

図 2.1

Star2D の計算領域概略

イメージとして各物理量は，Z 軸周りの対称性を持つことになる．計算は二次元の r-θ
体系で行われる．計算領域は r 方向，θ 方向をそれぞれ 1 からn𝑟 ，1 からn𝜃 個の計算格子に分けている．計算格子の座標値は，格子境界(ib , jb )で定義され，下付き添え字 b を用いて 𝑟b (ib )，𝜃b (jb )で表される．また，格子中心(ic , jc )での座標値を式(2.6)，式(2.7)のように定義している．

𝑟c (ic ) =𝑟b (ib ) + 𝑟b (ib − 1)
2

(2.6)

𝜃c (jc ) =𝜃b (jb ) + 𝜃b (jb − 1)
2

(2.7)

8

物理量は格子中心で定義され，r 方向速度は𝒗𝑟c (ic )，θ 方向速度は𝒗𝜃c (ic )，質量密度は𝜌c (ic , jc )，イオン温度は𝑇ic (ic , jc )，電子温度は𝑇ec (ic , jc )，イオン圧力は𝑝ic (ic , jc )，電子圧力は𝑝ec (ic , jc )で表される．
境界条件の設定を容易にするために，上下左右の端に 1 ずつダミーメッシュを置いている．ダミーメッシュには，内側のメッシュで定義された物理データがコピーされる．反射
境界条件の場合には，境界の速度を 0 として，内側の格子での速度の符号を反転させたものを設定する．メッシュの位置については補外され，

𝜃b (nym2) = 2𝜃b (nym2 + 1) − 𝜃b (nym2 + 2)

(2.8)

𝜃b (nyp2 − 1) = 2𝜃b (nyp2 − 2) − 𝜃b (nyp2 − 3)

(2.9)

ここで，nym2は 0°側の境界部分のダミーメッシュ番号，nyp2は逆に 180°側の境界部分のダミーメッシュ番号を指す．境界をjb = 0 (0°側)，jb = ny (180°側)とすると，jb = −1 ，
jb = ny + 1 の位置である．図 2.2 に θ 方向の計算格子のイメージを示す．また，r 方向についても同じように書ける．

図 2.2

θ 方向の計算格子のイメージ

9

次に，物理量の更新について記す．物理量は運動方程式に従い更新される．
・速度の更新

n+

1
2

n−

𝒖ib − 𝒖ib
𝑡

n+

1
2

n
1 = 𝑆ib

1
n−
2−𝑡 2

(𝑝 + 𝑞)nic+1 − (𝑝 + 𝑞)nicn∆𝑚ib(2.10)

1
n
n
(∆𝑚ic
)
+ ∆𝑚ic+12(2.11)n∆𝑚ib=0=1(∆𝑚1n )
2

(2.12)n∆𝑚ib=ny
=

1
n
(∆𝑚ny
)
2

(2.13)∆𝑚ic∆𝑉icn+1(2.14)n∆𝑚ib=

n
ここで，𝑆ibはメッシュの面積を示す．

・密度の更新n𝜌ic=

ここで，𝑉icn はメッシュの体積である．

10

2.3

3D hybrid code について

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

2.3.1

基礎方程式

計算に使用する方程式を導出するための基礎方程式は次のように与えられる．以下の方程式は直交座標系および，MKSA 単位系で記述されている．

イオンの運動方程式

電子の運動方程式

アンペアの法則

𝑚i

𝑛e 𝑚e

𝑑𝒗i
= 𝑍𝑒(𝑬 + 𝒗i × 𝑩)𝑑𝑡𝑑𝒗e
= −𝑒𝑛e (𝑬 + 𝒗e × 𝑩) − ∇𝑃e𝑑𝑡∇ × 𝑩p = 𝜇0 (𝑱e + 𝑱i )(2.15)(2.16)(2.17)

ファラデーの電磁誘導の法則

𝜕𝑩= −∇ × 𝑬𝜕𝑡(2.18)

電子の状態方程式

𝑃e = 𝑛e 𝑇e(2.19)

11

以降，計算に使用する方程式の導出を説明する．まず粒子として取り扱うイオンの運動方程式を示す．

𝑑𝒗i= 𝑍𝑒(𝑬 + 𝒗i × 𝑩)
𝑑𝑡

(2.15)

𝑑𝒙i 𝑑𝒚i 𝑑𝒛i)𝒗i = (,,
𝑑𝑡 𝑑𝑡 𝑑𝑡

(2.20)

𝑑𝒙i
𝑑𝒚i
𝑑𝒛i, 𝒗𝑦i =
, 𝒗𝑧i =
𝑑𝑡
𝑑𝑡
𝑑𝑡

(2.21)𝑚i

𝒗𝑥i =

ここで𝑚i はイオンの質量 ，𝒗i はイオンの速度、 𝑍 はイオンの価数，𝑒 は電気素量，
𝒙i , 𝒚i , 𝒛i はイオンの位置，𝒗𝑥i , 𝒗𝑦i , 𝒗𝑧i はそれぞれ x，y，z 方向のイオン速度，また，
𝑬，𝑩はそれぞれ電場，磁場を表している．

流体として扱う電子の運動方程式を電場の解法に用いる．電子の運動方程式は以下の式で表される．

𝑛e 𝑚e𝑑𝒗e= −𝑒𝑛e (𝑬 + 𝒗e × 𝑩) − ∇𝑃e𝑑𝑡(2.16)

ここで𝒗e は電子の速度，𝑚e は電子の質量，𝑛e は電子数密度，𝑃e は電子圧力である．電子は慣性無しの流体としてみなすので，慣性による効果は無視し，式(2.16)の左辺は 0 にできる(𝑚e → 0)．さらに電子圧力𝑃e に関しては電子の状態方程式より以下のようになる．

𝑃e = 𝑛e 𝑇e(2.19)

ここで𝑇e は電子温度である．式(2.19)を式(2.16)に代入して𝑃e を消去する．

0 = −𝑒𝑛e (𝑬 + 𝒗e × 𝑩) − ∇(𝑛e 𝑇e )12

(2.22)

荷電準中性条件および，ダーウィン近似を用いたアンペアの法則により以下の式が得られる．
∇ × 𝑩p = 𝜇0 (𝑱e + 𝑱i )(2.17)

ここで𝜇0 は真空中の誘電率であり，𝑱i ，𝑱e はそれぞれイオン電流密度，電子電流密度である．また，𝑩p はプラズマにより生成される磁場で以下のようになる．
𝑩p = 𝑩 − 𝑩ext

(2.23)

ここで，𝑩 は計算領域内部の合成磁場，𝑩ext は外部磁場である．𝑱i はイオンの数密度と電流密度から計算することが出来る．𝑱e は電子電流密度の定義により以下のようになる．

𝑱e = −𝑒𝑛e 𝒗e(2.24)

式(2.17)，(2.24)を用いて𝑱e を消去し𝒗e について求めると以下のようになる．

𝒗e = −

1 1( ∇ × 𝑩p − 𝑱i )𝑒𝑛e 𝜇0(2.25)

これを式(2.22)に代入して電場𝑬 について整理する．荷電準中性条件より𝑒𝑛e = 𝑍𝑒𝑛i を考慮し次式を得る．

𝑬=

1 1{ (∇ × 𝑩p ) × 𝑩 − 𝑱i × 𝑩 − 𝑍𝑇e ∇𝑛i }𝑍𝑒𝑛i 𝜇0(2.26)

ただし，式(2.26)から得られる電場はプラズマが存在する領域内でのみ使用され，プラズ
マが存在しない真空領域での電場は後節で述べるラプラス方程式から得られる．
プラズマによる磁場を記述する方程式として，以下のファラデーの電磁誘導の法則を用いる．
𝜕𝑩= −∇ × 𝑬𝜕𝑡
13(2.18)

2.3.2

計算方法

式(2.15)から式(2.26)を数値計算によって解くためには，基礎式を差分式に直し，数値的に安定で，物理的に妥当な空間幅、時間幅で解かなければならない．
計算の大まかな手順を以下に示す．
[1]外部磁場および，イオンの粒子位置，速度の初期設定を行う．
[2]格子点上でのイオンの数密度，電流密度を計算する．
[3](2.25)式により，電場を求める．
[4][1]で設定した外部磁場と，[3]で求めたプラズマによる電場と磁場から，式(2.15)，式(2.20)を用いて新たなイオン粒子の位置、速度を求める。
[5] 式(2.26)により磁場を求める．
[2]～[5]を繰り返し，プラズマの時間的変化を追う．
3D hybrid code に対する各方程式の差分式を以下の小節で示す．

・プラズマ粒子の運動の計算
3D hybrid code では計算負荷を減らすために，いくつかのイオン粒子を 1 つの大きな
電荷，大きな質量をもった粒子と見なす。このような粒子を超粒子と呼ぶ．超粒子を個別に運動させることでプラズマ粒子の運動を追う．今回の計算では超粒子数を 100 万個としてシミュレーションを行った．つまり， 式(2.20)に，leap-frog 法(18)による時間中心差分を施すと次式となる．

1
𝒙𝑛+1
− 𝒙𝑛i
𝑛+
i
= 𝒗i 2∆𝑡(2.27)

ここで各変数の上付き添え字は時間を表し，下付き添え字はi 番目のイオンを表す．また，∆𝑡 は時間幅を表す．
式(2.15)に時間中心差分を施すと以下の式となる．

𝑛+

𝒗i

1
2

𝑛−

− 𝒗i
∆𝑡

1
2

1
𝑍𝑒 𝑛 1 𝑛+12
𝑛−=
{𝑬 + (𝒗i − 𝒗i 2 ) × 𝑩𝑛 }
𝑚i
2

14

(2.28)

1

∆𝑡

1

∆𝑡

𝑛 − 2 はある時刻𝑛 から 2 だけ以前の値であることを示す．同様に，𝑛 + 2 は 2 だけ進ん1だ値であることを示す．上付き添え字が𝑛 ，𝑛 − 2 である変数は既知量である．
𝑛+1/2

プラズマ粒子の運動を求めるためには，まず式(2.28)において，𝒗i
𝑛+1/2

未知量である𝒗i𝑛+1/2を解く．次に，求めた𝒗iを左辺にまとめて，を式(2.27)の右辺に代入し，新たな粒子の位置𝒙𝑛+1を求めることが可能となる．式(2.28)を変形すると次式となる．
i1

𝑨𝒗𝑛+2 = 𝑺

(2.29)

𝛼=

𝑞∆𝑡 𝑛
𝐵
2𝑚i 𝑥

(2.30)

𝛽=

𝑞∆𝑡 𝑛
𝐵
2𝑚i 𝑦

(2.31)

𝛾=

𝑞∆𝑡 𝑛
𝐵
2𝑚i 𝑧

(2.32)

𝛼 , 𝛽 , 𝛾 , 𝛿 を次のように定義する．

𝛿 = 1 + 𝛼 2 + 𝛽2 + 𝛾 2(2.33)

なお．ここでの𝑬 ，𝑩 は粒子位置におけるそれぞれの値を示している．これらを用いて，
行列𝑨 とベクトル𝑺 は以下のように構成される．

1
𝑨=( 𝛾
−𝛽

−𝛾
1
𝛼

15

𝛽
−α)
1

(2.34)

𝑚𝑣i𝑥+ 𝑐𝐸𝑥𝑛 + 𝛾𝑣𝑦𝑚 − 𝛽𝑣𝑧𝑚𝑚
𝑺 = (𝑣i𝑦 + 𝑐𝐸𝑦𝑛 + 𝛾𝑣𝑧𝑚 − 𝛽𝑣𝑥𝑚 )𝑚𝑣i𝑧+ 𝑐𝐸𝑧𝑛 + 𝛾𝑣𝑥𝑚 − 𝛽𝑣𝑦𝑚𝑞∆𝑡(2.35)

1

ここで，𝑐 = 𝑚 ，𝑚 = 𝑛 − 2 である．さらに，式(2.29)を次のように変形して，速度を求めiる式が得られる．

1

𝒗𝑛+2 = 𝑨−1 𝑺(2.36)

ここで，𝑨−1 は以下のように得られる．

1 + 𝛼21
𝑨−1 = (𝛼𝛽 − 𝛾𝛿𝛼𝛾 + 𝛽𝛼𝛽 + 𝛾1 + 𝛽2𝛽𝛾 − 𝛼𝛼𝛾 − 𝛽𝛽𝛾 + 𝛼 )1 + 𝛾2(2.37)

・イオン数密度，電流密度
イオン数密度を求めるためには，空間内に広がっているプラズマを格子点上に割り付
けなければならない．そのために，PIC（Particle In Cell）法を用いる．今，１つの超粒子の中に N 個のイオンが含まれている場合を考える．図 2.3 に示すように，格子点内に存在する超粒子の位置によって，8 個の格子点によって囲まれた体積 V の立方体を，𝑉1 ，𝑉2 ，～𝑉8 の体積を持つ直方体に分割する．格子点(i，j，k)に割り振られる超粒子に含まれるイオンの数は次式で与えられる．

𝛿𝑛i (i，j，k) =𝑉2𝑁𝑉(2.38)

すべての粒子についての和を𝛿𝑛i についてとれば，ある格子点上の𝑛i を定義することが出来る．また，イオン電流密度𝑱i は，数密度𝑛i を用いて次式により求めることが出来る．

𝑱i = 𝑍𝑒𝑛i 𝒗i
16(2.39)

図 2.3

PIC (Particle In Cell)法

・電磁場の計算方法(1) プラズマ領域での計算
電磁場を求めるためには，式(2.25)，式(2.26)を用いる．式(2.26)を片側差分して，𝑩𝑛+1/2 だけを左辺に残すと次式となる．

1

𝑩𝑛+2 = 𝑩𝑛 −∆𝑡
∇ × 𝑬𝑛2(2.40)

また，式(2.25)は以下の式になる．

𝑬𝑛 =

1 1
1 𝑛
𝑇e
{
(∇ × 𝑩𝑛p ) × 𝑩𝑛 −𝑱i × 𝑩𝑛 − ∇𝑛i𝑛 }𝑛i 𝜇0 𝑍𝑒𝑍𝑒𝑒

(2.41)

ここで，𝑩𝑛 ，𝑬𝑛 は既知であるとすると，式(2.40)から𝑩𝑛+1/2 を求めることが出来る．次に，
𝑛+1/2

𝑩𝑛+1/2 ，𝑱𝑛+1/2 ，𝒗𝑛+1/2 ，𝑛iを式(2.41)に代入して𝑬𝑛+1/2 を求め，線形外挿より電17𝑛+1
磁場の予測値𝑬𝑛+1
pred ，𝑩pred を次式より求める．

1

𝑛+
𝑛
2
𝑬𝑛+1
pred = −𝑬 + 2𝑬

1

𝑛+
2−
𝑩𝑛+1
pred = 𝑩∆𝑡∇ × 𝑬𝑛+1
pred2(2.42)(2.43)𝑛+1
ここで，下付き添え字pred は予測値であることを表している．𝑬𝑛+1
pred ，𝑩pred を求めた後，
𝑛+3/2

𝑛+3/2

𝑛+3/2

それらによるイオンの運動を解き𝒗pred ，𝑛i pred ，𝑱i 𝑛+3/2を求める．次に，𝑩pred を次pred式により求める．

𝑛+

3

2
𝑩pred
= 𝑩𝑛+1pred −

1∆𝑡∇ × 𝑬𝑛+22(2.44)

𝑛+3/2

そして，式(2.41)で𝑬𝑛 を求めた時と同様にして，𝑬pred を求める．最後に次式に示す線
形内挿から次のタイムステップの値𝑬𝑛+1 ，𝑩𝑛+1 を求める．

𝑬𝑛+1

1 𝑛+1 1 𝑛+32
= 𝑬 2 + 𝑬pred22

𝑩𝑛+1 = 𝑩𝑛+2 −∆𝑡∇ × 𝑬𝑛+1
2

𝑛+3/2

𝑛+3/2

1

これらの値𝑬𝑛+1 ，𝑩𝑛+1 より𝒗i

𝑛+3/2

，𝑛i

，𝑱i

(2.45)(2.46)を求めることが出来る．ここで，電場を求める際に用いた一連の手法は，予測子－修正子法と呼ばれる収束計算の一種であ3る．この方法では，𝑛 タイムステップから𝑛 + 1 タイムステップの値を求めるために，𝑛 + 2
タイムステップの予測値を求め内挿計算を行う．そして，求めた𝑛 + 1 タイムステップの値を再度予測値として用いる．具体的には，式(2.46)により求められた𝑩𝑛+1 を次の予測値
(19)
𝑩𝑛+1
を満足するまで繰り返し計
pred とし，式(2.42)から式(2.46)までを次式で定義する条件

算を行う．
18

𝑛+1(𝑞+1)

maxi，j，k |𝐵i，j，k𝑛+1(𝑞)− 𝐵i，j，k |𝑛+1(𝑞)
maxi，j，k |𝐵
|
i，j，k

≤𝜀

(2.47)

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

(2) 真空領域での電場の計算
真空領域では，ラプラス方程式∇2 𝑬 = 𝟎 を SOR（Successive Over Relaxation）法を用いて解く．SOR における加速係数𝛼 は正方形領域の場合は次のように表される．

𝛼=2π1 + sin (NGMX)(2.48)

ここでNGMX は x，y，または z 方向の格子数である．𝛼 の値は通常 1.8 から 1.9 付近である．また，𝛼 = 1 の場合には，SOR 法はガウス－ザイデル法に等しい．
SOR 法では，次式によって変数の収束を加速する．

𝐸p (i，j，k) = 𝛼{𝐸p (i − 1，j，k) + 𝐸p (i + 1，j，k) + 𝐸p (i，j − 1，k)+𝐸p (i，j + 1，k) + 𝐸p (i，j，k − 1) + 𝐸p (i，j，k + 1)}/6(2.49)+(1 − 𝛼)𝐸p (i，j，k) ，(p = 𝑥，𝑦，𝑧)
次式で定義される残差𝜉𝑥，𝑦，𝑧 が、あらかじめ決められた値以下になれば真空部分の電場の計算を終了する．

19

𝜉p = 𝐸p (i − 1，j，k) + 𝐸p (i + 1，j，k) + 𝐸p (i，j − 1，k)+𝐸p (i，j + 1，k) + 𝐸p (i，j，k − 1) + 𝐸p (i，j，k + 1)

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

new
𝜙i，j，k←side
corner𝑀𝜙i，j，k + 𝑆Φi，j，k+ 𝐾Φi，j，k

𝑀 + 6𝑆 + 12𝐾(2.51)

new
𝜙i，j，k を格子点上に定義される値とし，𝜙i，j，k
は Digital Filter を通した新しい𝜙i，j，k を表
side
cornerす．また，Φi，j，kとΦi，j，kは side term と corner term であり，それぞれ以下の式で表され

る．

side
Φi，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(2.52)

cornerΦi，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(2.53)+𝜙i−1，j−1，k + 𝜙i−1，j+1，k + 𝜙i+1，j+1，k + 𝜙i+1，j−1，k

ここで，M，S，Kの値をそれぞれ 12，4，1 とすると式(2.51)は次式となる．

new
𝜙i，j，k←1
1 side1 corner
𝜙i，j，k + Φi，j，k+ Φi，j，k
4
12
48

20

(2.54)

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

・計算制限-CFL 条件
Hybrid code において，差分化の時間幅Δ𝑡 ，空間幅∆𝑥 は次の条件を満たす必要がある．
CFL（Courant-Friedrichs-Lewy）条件

Δ𝑡 ≤

∆𝑥
𝑽A

(2.55)

これは，Δ𝑡 ，∆𝑥 → 0 で解が収束する条件である．𝑽A はアルフベン速度であり次式で定
義する．

𝑽A ≡

𝐵0
√𝜇0 𝑛i 𝑚i(2.56)

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

𝜔ci Δ𝑡 < 0.2(2.57)

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

∆𝑥 >

𝑐
𝜔pi

𝑐

(2.58)

この𝜔 より∆𝑥 が小さい場合にも Hybrid code は発散することなく計算を行えるが，その解piが正しいものである保証はない．ここで，𝜔pi はイオンプラズマ周波数である．

21

2.4

カップリングコード

このコードは輻射流体コードのデータを 3D hybrid code で扱えるようにするためにデータ変換をし，2 つのコード間を連結する役割を持つ．具体的には単位系，座標系，次元数，
またイオンの取り扱い方などの違いに対するデータ変換を行う．
本研究において，カップリングで 3D hybrid code に受け渡しが必要なデータは「イオン位置座標，イオン速度，イオン電荷，イオン質量」である．ここではイオンに関連する物理量のみで電子に関する物理量が含まれない.
式(2.26)からわかるように，3D hybrid code の計算で必要な電子の物理量は，電子温度
𝑇e のみである．しかし，過去の計算結果から電子温度𝑇e は 3D hybrid code の計算結果にあまり影響を与えないことがわかっているので，計算時間短縮のため𝑇e = 0 とした．そのため，電子に関する物理量の受け渡しは必要ない．

これらの必要なデータは 2 次元輻射流体コード(STAR2D)で求められた各計算格子の境界で定義される位置座標と，各計算格子の中心で定義される速度，質量密度，イオン
数密度，イオン電荷を変換して得る．

行ったデータ変換の流れは以下のようになる．

①2 次元輻射流体コードのデータの CGS 単位系から MKS 単位系への変換する．
②流体であったイオンを粒子に置き換え，イオン粒子を“超粒子”として扱う．
③座標系を変換する．

以下にその詳細を示す．

①2 次元輻射流体コードのデータの CGS 単位系から MKS 単位系への変換
表 2.1 に示すように単位系を変換する．

22

表 2.1

カップリングコードでの単位系変換
CGS

→

MKS

𝑟 [𝑐𝑚]→

𝑟 × 10−2 [𝑚]𝑣𝑟 [𝑐𝑚/𝑠]→

𝑣𝑟 × 10−2 [𝑚/𝑠]𝑣𝜃 [𝑐𝑚/𝑠]→𝑣𝜃 × 10−2 [𝑚/𝑠]Mass density𝜌 [𝑔/𝑐𝑚3 ]→𝜌 × 103 [𝑘𝑔/𝑚3 ]Ion density

𝑁i [𝑐𝑚−3 ]→𝑁i × 106 [𝑚−3 ]Position coordinate

Velocity②イオン粒子の“超粒子”としての取り扱い
二次元輻射流体コードにて計算したイオン流体より，超粒子へと変換する．前述したように，3D hybrid code では計算負荷を減らすために，いくつかのイオン粒子を超粒子とお𝑛く．実際の粒子数を𝑛 ，超粒子数を𝑛s としたとき，粒子数比ratio𝑛 = 𝑛s を用いて計算を行った．流体から粒子に変換するために，イオンの数密度を用いて粒子として配置していく．
また，超粒子という仮定をおいたため，超粒子数の電荷及び質量を定める必要がある．
二次元輻射流体コードのi , j 番目の格子の質量密度𝜌 (i, j)，イオン数密度𝑁i (i, j)，イオ
ンの平均電荷𝑣𝑎𝑙𝑒𝑛𝑐𝑒 (i, j)，体積𝑉 (i, j) [式(2.61)より]から式(2.59)及び式(2.60)のように超粒子の平均電荷𝑣𝑎𝑙𝑒𝑛𝑐𝑒s ，平均質量𝑚𝑎𝑠𝑠s を求めた．

𝑣𝑎𝑙𝑒𝑛𝑐𝑒s =∑𝑛i=1 ∑𝑛j=1 𝑁i (i, j) × 𝑉(i, j) × 𝑣𝑎𝑙𝑒𝑛𝑐𝑒(i, j)

𝑚𝑎𝑠𝑠s =𝑛s∑𝑛i=1 ∑𝑛j=1 𝜌(i, j) × 𝑉(i, j)𝑛s

23(2.59)(2.60)

③二次元座標系 r-θ から三次元座標系 x-y-z への変換
二次元輻射流体コードの計算データのある一つの計算格子(i, j)のイメージを図 2.4，
2.5 に示す．前述したように物理量は軸対称性(z 軸対称性)を念頭に設定してあるので，物理量の軸対称性を仮定すると，二次元輻射流体コードの計算データは𝜑 を0 → 2π で積分することで図 2.6 のように拡張できる．

図 2.4

二次元輻射流体コードの計算格子(i, j)のイメージ

図 2.5 計算格子(i, j)のイメージ

24

図 2.6

計算格子(i, j)拡張のイメージ

図 2.4 の要素の体積𝑉 (i, j)，イオン粒子数𝑛 (i, j)は，式(2.61)，(2.62)のようになる．

2
3) × (cos 𝜃j − cos 𝜃j−1 )𝑉(i, j) = − π(𝑟i3 − 𝑟i−1
3(2.61)𝑛(i, j) = 𝑁ic (i, j) × 𝑉(i, j) × ratio(2.62)

ここで，π は円周率，𝑟i ，𝜃j は𝑟 方向 i 番目，𝜃 方向 j 番目の格子の境界の位置座標(図2.4)，𝑁ic (i, j) はその格子中心で定義されるイオン数密度[m−3 ]，ratio は前述した実際の粒子数と超粒子数の粒子数比である．
3D hybrid code ではイオンを粒子として取り扱うので，各格子ごとに各格子の粒子数だ
けイオンの位置座標，イオン速度を設定する．ここで分布が均等になるように乱数(random number)を使って式(2.63)，(2.64)，(2.65)のように定義した．
𝑟∫𝑟

4π𝑟 2 𝑑𝑟∫𝑟

4π𝑟 2 𝑑𝑟

i−1
𝑟i
i−1

25

= 𝑈1

(2.63)

𝜃

∫𝜃

𝑗−1

sin 𝜃 𝑑𝜃𝜃∫𝜃 j sin 𝜃 𝑑𝜃(2.64)= 𝑈2

j−1(2.65)

𝜑 = 2π𝑈3式(2.63)，(2.64)より，𝑟 ，𝜃 は式(2.66)，式(2.67)のように表せる．
1

3
3
)}3
𝑟 = {𝑟i−1+ 𝑈1 (𝑟i3 − 𝑟i−1

(2.66)cos 𝜃 = cos 𝜃j−1 + 𝑈2 (cos 𝜃j − cos 𝜃j−1 )(2.67)式 (2.66) ， (2.67) ， (2.65) よ り ， イ オ ン の 位 置 座 標 (𝑥, 𝑦, 𝑧) は 式 (2.68) ， イ オ ン 速 度(𝑣𝑥 , 𝑣𝑦 , 𝑣𝑧 )は式(2.69)のように設定する．下付き添え字は超粒子数𝑛s の中のk 番目の粒子であることを示している．

𝑥k = 𝑟 sin 𝜃 cos 𝜑{ 𝑦k = 𝑟 sin 𝜃 sin 𝜑𝑧k = 𝑟 cos 𝜃(2.68)

v𝑥 = v𝑟 sin 𝜃 cos 𝜑 + v𝜃 cos 𝜃 cos 𝜑{ v𝑦 = v𝑟 sin 𝜃 sin 𝜑 + v𝜃 cos 𝜃 sin 𝜑v𝑧 = v𝑟 cos 𝜃 − v𝜃 sin 𝜃(2.69)

26

2.5

計算体系

各コードにおける計算モデル及び計算体系を以下に記す．

2.5.1

Star2D における計算体系

図 2.7，2.8 に Star2D における球状及び円錐ターゲットの計算モデルを示す．また，表2.2 に Star2D における球状及び円錐ターゲットの初期パラメータを示す．
各ターゲットの材質は H で，Star2D では流体近似をするため，本来真空である背景領
域には低密度かつ低温の H プラズマを背景ガスとして置いている．
円錐ターゲットの質量と体積は VISTA で考えられているターゲット質量(50 g)の 1/500とした．また，推力は質量に依存するため，比較対象である球状ターゲットの質量も同質量とし，この質量を基に各ターゲットの質量密度と寸法を決定した．
計算は前述の通り，r-θ 体系で，球状ターゲットの計算領域は 6 cm で，メッシュ分割数
は，半径方向 120 メッシュ，周方向 10 メッシュとした．また，円錐ターゲットは計算領域 12cm，半径方向 240 メッシュ，周方向 180 メッシュとした．
また本コードでは，核融合反応後，高温プラズマ状態にある核融合燃料が膨張し，低温状態にある推進剤を加熱して膨張する過程を計算するが，核融合プラズマの内部エネ
ルギーは，VISTA モデルの 1/500，レーザーエネルギー5 MJ，エネルギーゲイン 750，また核融合エネルギーの 40%が荷電粒子の寄与として考え，3 MJ とした．この内部エネル
ギーより，核融合プラズマのイオン温度𝑇i 及び電子温度𝑇e を求める．内部エネルギーU
は，式(2.70)のように表せる(21)．

𝑈=

𝑝𝑉
𝛾−1

(2.70)

ここで，𝑝 は圧力，𝑉 は体積，𝛾 は比熱比である．また，圧力 𝑝 は式(2.71)のように表せる．

𝑝 = 𝑝i + 𝑝e

𝑝i 及び 𝑝e はイオン圧力，電子圧力を表し，それぞれ式(2.72)のように表せる．
27(2.71)

k 1∙ 𝜌𝑇
𝑚H 𝑀 i
k 𝑍
𝑝e =
∙ 𝜌𝑇
{
𝑚H 𝑀 e𝑝i =(2.72)

ここで，𝑘 はボルツマン定数，𝑚H はプロトンの質量，𝑀 は質量数，𝜌 は質量密度，𝑍 は電荷数である．これらの式(2.71)及び式(2.72)により，式(2.70)は式(2.73)のように表せる．

𝑈=k𝜌(𝑇i + 𝑍𝑇e )𝑉𝑚H 𝑀(𝛾 − 1)

(2.73)

この理想気体の内部エネルギー式である式(2.73)を用いて核融合プラズマのイオン温度及び電子温度を𝑇i = 𝑇e として求めた．
また，円錐ターゲットの推進剤は球状ターゲットの推進剤に比べて長いため，中性子による加熱を過小評価しないために荷電粒子エネルギー(3MJ)の約 0.1%として考慮し推進剤温度を決定した．

図 2.7

Star2D における球状ターゲットの計算モデル28

図 2.8

Star2D における円錐ターゲットの計算モデル

表 2.2

Star2D における各ターゲットの初期パラメータ
Spherical

Conical1.044 × 105

Fusion plasma temperature [eV]
Propellant temperature [eV]

0.415

0.025

Background gas temperature [eV]0.025

Mass [mg]100

また，推進剤中で中性子を減衰させることが見込める円錐ターゲットについて，今回の計算体系に対して，中性子減衰率を評価する．中性子束の減衰率は式(2.74)のように表せる(22)．

𝑥𝑁= 𝑒 −𝜆𝑁0(2.74)

ここで， 𝑁0 は物質に入射する中性子束，𝑁 は𝑥 地点での中性子束，𝜆 は物質の平均自由行程を示す．また，平均自由行程 𝜆 は式(2.75)のように表せる．
29

𝜆=

1
𝜎𝑛

(2.75)

ここで，𝜎 は散乱過程を表す断面積，𝑛 は物質の数密度を表す．円錐推進剤の質量密度は7.655 × 10−2 g/cm3 ，陽子の質量は1.673 × 10−24 g より，数密度 n は4.576 ×1022 cm−3 である．また，式(2.75)と 14 MeV 中性子に対する水素の断面積 (23) 𝜎 =6.878 × 10−25 cm2 より，平均自由行程は 31.77 cm である．また，円錐推進剤の長さは3

今回の体系の 500 倍の体積であるため，𝑥 = 1.85 × √500 = 14.68 cm である．以上より，
63

円錐ターゲットの推進剤中での中性子の減衰率は100 となる．

30

2.5.2

3D hybrid code における計算体系

図 2.9 に 3D hybrid code の計算体系，表 2.3 に 3D hybrid code における初期パラメー
タを示す．円錐ターゲットの推進剤がコイル方向を向くように，コイル角度を決定し，それに応じコイル位置，コイル半径のパラメータを決定した．
球状ターゲットの計算領域は，x，y，z 各方向±1.5 m，円錐ターゲットの計算領域は，x，y 方向は±2 m，z 方向は±2.5 m とした．メッシュ分割数は各ターゲット，各方向ともに 5 cm
ずつ均等に区切った 120 メッシュとした．

図 2.9

表 2.3

3D hybrid code の計算体系

3D hybrid code における初期パラメータ

SCM current [MA]3

Magnetic field strength at the center of coil [T]10

Magnetic field strength at the plasma position [T]

1.4

31

第3章
3.1

計算結果

評価方法

本研究の目的は，球状ターゲットと円錐ターゲットの推進性能を比較することであるが，パルス駆動であるレーザー核融合ロケットでは推力が時間とともに変化するため，シミュレ
ーション結果の評価方法は推力の時間積分である力積(Impulse bit)を用いた．2.3 節で説明した 3D hybrid code で運動量保存が成り立つとして，式(3.1)で与えられる Impulse
bit 𝐼 によって評価した．

𝐼 = −∆ [∑ 𝑚𝑣𝑧 ] = − [∑ 𝑚𝑣𝑧 − ∑ 𝑚𝑣𝑧0 ](3.1)

Z 軸の負の方向が宇宙機の推進方向と考え，− ∑ 𝑚𝑣𝑧 は Z 方向の運動量総和，∑ 𝑚𝑣𝑧0 は Z 方向の初期運動量総和となる．
また，もう一つの評価方法として，式(3.2)で与えられる運動量変換効率 𝜂 を用いた．

𝜂=𝐼−∆[∑ 𝑚𝑣𝑧 ] −[∑ 𝑚𝑣𝑧 − ∑ 𝑚𝑣𝑧0 ]==∑ 2𝑚|𝑣0 |∑ 2𝑚|𝑣0 |∑ 2𝑚|𝑣0 |(3.2)

ここで，∑ 𝑚|𝑣0 | は初期運動量の絶対値の総和を示す．

3.2

球状ターゲット

以下に，球状ターゲットにおける Star2D 及び 3D hybrid code の結果を示す．

3.2.1

Star2D の計算結果

図 3.1，3.2 に Star2D における球状ターゲットのイオン数密度，イオン温度の分布を示す．図 3.1，3.2 より，イオン数密度が1029 m−3 であった低温状態の推進剤プラズマが，核融合プラズマとの衝突により加熱され，等方膨張していくことが確認できる．
32

また，図 3.3 に Star2D における球状ターゲットのエネルギー推移を示す．縦軸エネル
ギー [MJ]，横軸が時間 [ns]を表し，紫が運動エネルギー，緑が内部エネルギー，青が輻射エネルギー，黄色がこれらの合計エネルギーとなっている．合計エネルギーは，初期状態の 3 MJ に対し，2.4 MJ で一定となっている．
また， 3D hybrid code では内部エネルギー(プラズマ温度)を考慮せず，低温(電子・イ
オン温度をゼロ)を仮定するので，内部エネルギーが十分下がり，運動エネルギーが最大となるタイミングを 3D hybrid code の初期状態とする．よって，球状ターゲットで運動エネ
ルギーが 2.0 MJ で最大となる 175 ns 時点でのデータを 3D hybrid code での初期状態として渡した．

図 3.1

イオン数密度分布[(a)：初期状態，(b)：50 ns 後，(c)：125 ns 後，(d)：175 ns 後]33

図 3.2

イオン温度分布[(a): 初期状態，(b)50 ns 後，(c)125 ns 後，(d)175 ns 後]

図 3.3

球状ターゲットのエネルギー推移

34

図 3.4

Star2D による終状態におけるイオン数密度(左図)と流体速度(右図)の径方向分布

図 3.5

3D hybrid code での初期状態における全粒子分布(左図)と Y=0 面での粒子分布(右図)

図 3.4 に Star2D による終状態のデータの r 方向のそれぞれイオン数密度，絶対速度(√𝑣𝑟2 + 𝑣𝜃2 )を示す．このデータを拡張したイオン超粒子分布を図 3.5 に示す．拡張されたプラズマは半径 6 cm である．図 3.5 の右図は Y=0 の断面図であり，図 3.4 のイオン数
密度分布に示すように，プラズマ中のイオン超粒子は密度が高い r=3 cm 付近に多く配置され，密度が低いプラズマの内側と外側は超粒子が少なく配置されるため，このプラズ
マの形状は球殻となっている．
表 3.1 に 3D hybrid code の初期状態となる，カップリングしたプラズマの各パラメータを
35

示す．

表 3.1

球状ターゲットの 3D hybrid code でのプラズマの各初期パラメータ

Mass [kg]1.0 × 10−4

Valence

0.2

Number of super-particles

1.0 × 106

Number of ion particles per super-particle

6.0 × 1016

Kinetic energy [MJ]

3.2.2

2.0

3D hybrid code の計算結果

カップリングコードにより拡張されたデータを，3D hybrid code の初期状態として，磁場中でのプラズマの時間発展を計算した．図 3.6 に磁場中でのプラズマの挙動を示す．
Z=‐1.2 m にコイルがあり，初期状態としてプラズマは原点に配置される．プラズマが膨張するにつれて膨張方向が変わり，プラズマがキノコ型になり広がっているのが確認できる．

図 3.7 に Z 方向の運動量推移を示す．3.2 μs までは，磁力線に沿って‐Z 方向に進む
コイル側の粒子が多いため Z 方向の運動量は負の値となっている．3.2 μs 以降はプラズ
マが磁場に押され，+Z 方向の運動量が増加している．この結果から，6.7 μs で Impulse
bit は 12 Ns，運動量変換効率は 47%であった．

36

図 3.6球状ターゲットにおける磁場中でのプラズマの挙動[(a)：初期状態，(b)：3 μs 後，(c)：5 μs 後，(d)：6.7 μs 後]

図 3.7

球状ターゲットにおける Z 方向の運動量推移
37

3.3

円錐ターゲット

以下に，円錐ターゲットにおける Star2D 及び 3D hybrid code の結果を示す．

3.3.1

Star2D の計算結果

図 3.8，3.9 に Star2D における円錐ターゲットの各時間のそれぞれイオン数密度及びイ
オン温度の分布を示す．これらの図より，核融合プラズマとシェルは時間が経つにつれ，
イオン密度及び温度が小さくなり膨張するのに対し，推進剤は加熱されないことがわかる．
図 3.10 に Star2D における円錐ターゲットの各エネルギーの推移を示す．初期内部エ
ネルギーが 3 MJ に対し，合計エネルギーは 2.7 MJ で一定となっている．
運動エネルギーが最大となるのは，2.9 ns のときで 2.5 MJ である．この時点のデータを
3D hybrid code の初期状態とする．
ここで，図 3.11 及び 3.12 に Star2D による終状態の円錐推進剤の中心方向である
θ=133 °及び推進剤と反対側である θ=45 °でのイオン数密度，絶対速度の r 方向分布を示す．また図 3.13 に Star2D による終状態の円錐推進剤の中心方向である θ=133 °での電離度の r 方向分布を示す．図 3.11 及び 3.12 から，加熱されない推進剤部分は密度が高く，速度も持っていないことが確認できる．また図 3.13 から円錐推進剤の先端部以外は電離度が小さいことが確認できる．そのため，3D hybrid code に渡すときにはこのプラズ
マ化していない推進剤部分に粒子を配置しないようにイオン密度から粒子分布に変換する．具体的には，電離度が 0.9 以下，またはイオン数密度が5 × 1026 m−3 以上，絶対速度が 1 km/s 以下の部分を除いて超粒子に変換した．
図 3.13 に二次元データを三次元データに拡張したイオン超粒子分布を示す．円錐推進剤部分を取り除いたので，このプラズマは核融合プラズマとシェルによるものであり，
Y=0 の断面図のように粒子が偏って配置されている．プラズマ半径は 1.8 cm である．表
3.2 に 3D hybrid code で初期状態となる，カップリングしたプラズマの各パラメータを示す．

38

図 3.8

イオン数密度の時間変化図[(a)：初期状態，(b)：1.0 ns 後，(c)：2.0 ns 後，(d)：2.9 ns 後]

図 3.9

イオン温度の時間変化図[(a)：初期状態，(b)：1.0 ns 後，(c)：2.0 ns 後，(d)：2.9 ns 後]39

図 3.10

図 3.11

円錐ターゲットのエネルギー推移

Star2D による計算の終状態におけるイオン数密度の径方向分布[左：θ=133 °方向(円錐推進剤方向)，右：θ=45 °方向]40

図 3.12

Star2D による計算の終状態における流体速度の径方向分布[左：θ=133 °方向(円錐推進剤方向)，右：θ=45 °方向]

図 3.13

Star2D による計算の終状態における電離度の θ=133 °方向分布

図 3.143D hybrid code での初期状態における全粒子分布(左図)と Y=0 面での粒子分布(右図)41表 3.2

円錐ターゲットの 3D hybrid code でのプラズマの各初期パラメータ
Mass [kg]

5.8 × 10−7

Valence

0.99

Number of super-particles

1.0 × 106

Number of ion particles per super-particle

3.4 × 1014

Kinetic energy [MJ]

2.5

3.3.2

3D hybrid code の計算結果

図 3.15 に円錐ターゲットにおける磁場中でのプラズマの挙動を示す．コイルによる磁場によりコイル側に膨張するプラズマは減速，それ以外のところは膨張し，1 μs にはキノコ型になり広がっているのが確認できる．
図 3.16 に円錐ターゲットにおける Z 方向の運動量推移を示す．初期状態(0 μs)として，-Z 方向側のプラズマは本来推進剤があった部分に超粒子を配置していないため，-Z 方向側の超粒子が少ない分，質量及び-Z 方向の速度を持つ超粒子が+Z 方向側に比べ少ない．そのため初期状態に Z 方向の運動量は正の値をとる．
3D hybrid code の計算開始後，0.04 μs まで磁力線に沿って-Z 方向に膨張する粒子が多くなるが，その後プラズマが磁場に押し返されることによって+Z 方向の運動量が増える．
0.8 μs 以降は運動量の増減がなくなり，1 μs のとき，Impulse bit は 0.75 Ns，運動量変換効率は 55%となった．

42

図 3.15

円錐ターゲットにおける磁場中でのプラズマの挙動[(a)：初期状態，(b)：0.25 μs 後，(c)：0.50 μs 後，(d)：1.00 μs 後]

図 3.16

円錐ターゲットにおける Z 方向の運動量推移
43

第4章
4.1

考察

Star2D・円錐ターゲットにおける円錐推進剤が加熱しないことについて
Star2D 計算において，球状ターゲットは核融合プラズマが推進剤をプラズマ化させて膨張するのに対し，円錐ターゲットでは円錐推進剤部分が加熱されなかった．これは中心の核融合燃料から見る立体角が非常に小さいために，核融合プラズマ(流体)が推進剤部分とほとんど衝突しないため加熱されないと考えられる．
図 4.1，4.2，4.3 に θ=133 °(円錐推進剤方向)における初期状態，2.9 ns 及び 40 ns のときのイオン温度，絶対速度の分布を示す．なお，横軸は推進剤が配置されている 0.15
cm から 2 cm としている．

図 4.1

初期状態の円錐推進剤部分における r 方向の各分布[(a)イオン温度分布，(b)絶対速度分布]44図 4.2

2.9 ns 時の円錐推進剤部分における r 方向の各分布[(a)イオン温度分布，(b)絶対速度分布]

図 4.3

40 ns 時の円錐推進剤部分における r 方向の各分布[(a)イオン温度分布，(b)絶対速度分布]

図 4.4，4.5，4.6 より，初期状態に比べ，円錐推進剤部分の中心から見たときの先端部
分は，3D hybrid code の初期状態とする 2.9 ns には温度において 1 cm，絶対速度において 0.4 cm まで影響を及ぼしている．時間が経過した 40 ns のときにはイオン温度，絶対速度はともに 1 cm まで増加している．すなわち円錐推進剤全体のうち，推進剤先端部以外の 46%に変化はなかった．
また，表 4.1 に円錐推進剤部分のみの初期状態，2.9 ns 及び 40 ns のときの運動エネ
ルギー，内部エネルギー，その合計エネルギーを示す．

45

表 4.1

円錐推進剤部分の各時間における各エネルギー

Kinetic energy [kJ]Internal energy [kJ]

Total energy [kJ]

0 ns

0

6.0

6.0

2.9 ns

0.4

10

10.4

40 ns

0.1

5.3

5.4

表 4.1 より，2.9 ns のとき，計算体系全体の合計エネルギーに比べ，推進剤中の合計エネルギーは 0.33%であり，40 ns のときは 0.18%であった．このことから，円錐ターゲットの円錐推進剤は推力にほとんど寄与しないと考えられる．

・エネルギー保存について
図 3.3，及び図 4.1 より，Star2D において，初期に与えた 3 MJ と一定となるエネルギーが一致していない．初期に設定した 3 MJ に対し，球状ターゲットは 2.4 MJ，円錐ターゲ
ットは 2.7 MJ で一定となっている．初期状態からのエネルギー減少は，それぞれ 20%，
10%となった．この原因として，2 つの理由が考えられる．1 つ目は，コード上でイオン及び電子温度が上手く計算できていなかったことである．式(4.3)に Star2D で考えている保存型のオイラー方程式を示す．

𝜕𝑼 𝜕𝑭 𝜕𝑮++=𝑺
𝜕𝑡 𝜕𝑥 𝜕𝑦(4.3)

ここで，𝑼 は保存量，𝑭 ，𝑮 はそれぞれ x，y 方向の流束，𝑺 はソース項を表す．また，𝑼 ，𝑭 ，𝑮 は式(4.4)のように表せる．

𝜌𝜌𝑢
𝑼 = ( 𝜌𝑣 ) ,
𝜌𝑒t

𝜌𝑢
𝜌𝑣
𝜌𝑢𝑣
𝜌𝑢 + 𝑝T
)
(
𝑭=(
,
𝑮
=
𝜌𝑣𝑢
𝜌𝑣 2 + 𝑝T )(𝜌𝑒t + 𝑝T )𝑢(𝜌𝑒t + 𝑝T )𝑣
2

また，式(4.5)，(4.6)にイオン温度及び電子温度の時間発展式を示す．
46

(4.4)

𝜕𝑇i
+ (𝒗 ∙ ∇)𝑇i }𝜕𝑡𝜕𝑝i= − {𝑇i ( ) + 𝑞} (𝒗 ∙ ∇) + ∇ ∙ (𝜅i ∇𝑇i ) + 𝛼(𝑇e − 𝑇i )𝜕𝑇i 𝜌𝜌𝑐𝑣i {𝜕𝑇e+ (𝒗 ∙ ∇)𝑇e }𝜕𝑡𝜕𝑝e) } (∇ ∙ 𝒗) + ∇ ∙ (𝜅e ∇𝑇e ) + 𝛼(𝑇e − 𝑇i ) + 𝑆L + 𝑆rad
= − {𝑇e (
𝜕𝑇e 𝜌

(4.5)

𝜌𝑐𝑣e {

ここで，に各成分について表 4.2 に示す．

表 4.2

各成分の説明𝜌Mass density𝒗 = (𝑢, 𝑣)

Velocity1𝑒t = 𝑒int + |𝒗2 |2Specific total energy

𝑒int = 𝑒i int + 𝑒e int

Specific internal energy𝑒i = 𝑒i (𝜌, 𝑇i )(ion part)𝑒e = 𝑒e (𝜌, 𝑇e )(electron part)

𝑝T = 𝑝i + 𝑝e

Total pressure𝑝i = 𝑝i (𝜌, 𝑇i )Ion pressure𝑝e = 𝑝e (𝜌, 𝑇e )Electron pressure with cold part𝑇i

Ion temperature

𝑇e

Electron temperature

𝜅i

Ion conductivity

𝜅eElectron conductivity𝑞Artificial viscosity𝛼Thermal relaxation coefficient of Spitzer

𝑆L

Laser heating

𝑆rad

Radiation heating

47

(4.6)

式(4.3)より求まるのは，密度 𝜌，運動量 𝜌𝑢，合計エネルギー 𝜌𝑒t (運動エネルギー＋
内部エネルギー)であり，合計エネルギーから運動エネルギーを引くことで内部エネルギ
ーを求めることができるが，そこからイオン温度と電子温度の 2 つの温度を計算することができない．そのため，イオン温度と電子温度を求めるためには，別途式(4.5)及び式(4.6)に式(4.3)より求めた内部エネルギーから解く必要がある．そこから求めたイオンおよび電子温度をもとに，内部エネルギーを再計算して式(4.3)の内部エネルギーを置き換える，という計算を以前の Star2D では行われていた．その置き換えの際に内部エネルギーに誤差が生じ，初期の温度変化が激しい部分が上手く計算できず，エネルギーが保存されていなかった．
そこで，式(4.5)及び式(4.6)が移流項＋pdV 項からなっていることを利用し，式(4.3)の保存型から求まる単位体積当たりの内部エネルギー 𝜌𝑒t −

(𝜌𝑢)2
𝜌

と式(4.5)及び式(4.6)の時

間発展式の内部エネルギーが一致するように，pdV の圧縮 dV を調節することにより，双方の内部エネルギーが一致するときのイオン温度及び電子温度が求められるようになった．これによりエネルギーが保存されるようになった．
2 つ目は，メッシュの切り方である．どちらのターゲットも r 方向のメッシュ幅は均一に500 μm ずつとしていた．しかし，核融合プラズマ及び核融合プラズマと推進剤の間，推進剤内面近傍は，最初に核融合プラズマが膨張して推進剤に衝突し，そこから推進剤内に衝撃波が駆動されるとともに，反射衝撃波が再び中心に向かって進む．すなわち，今回の計算体系において，中心付近の領域の計算が最も激しい現象が起きるところである
ため，そこに十分な格子を置く必要があった．
そこで，改良した Star2D により，中心近傍は 50 μm，外側は 5000 μm と，不均一間隔にメッシュ幅を変更し，メッシュ数を 300 と以前よりも中心付近を細かくなるように変更した．
以上の変更を加え，新たに計算した球状及び円錐ターゲットの Star2D での計算結果
を図 4.4，4.5 に示す．これらの変更により，エネルギーが保存されるようになった．図 3.3，
3.9 の以前のデータと比較して，内部エネルギー及び運動エネルギーの増減の傾向もどちらのターゲットにおいても一致している．

48

図 4.4改良した Star2D における球状ターゲットのエネルギー推移

図 4.5

改良した Star2D における円錐ターゲットのエネルギー推移

49

4.2

3D hybrid code

・Z 方向の運動量について
図 3.6，3.15 より，球状ターゲットは 3.2 μs，円錐ターゲットは 0.04 μs まで-Z 方向の運動量が増えている．図 4.5 に球状ターゲットにおける，初期状態及び運動量がマイナスにピークの時である 2.2 μs 時の全超粒子の Z 方向の速度分布を示す．図 4.6 から，初期状
態と 2.2 μs 時を比較すると，2.2 μs 時-Z 方向の速度を持つ粒子が多くなり，+Z 方向に速度を持つ粒子の数が少なくなっている．さらに時間が経ち，3D hybrid code 終時間である
6.7 μs 時には+Z 方向の速度を持つ粒子の方が多くなっていることが確認できる．これは初期状態では等方プラズマであるが，2.2 μs 付近ではコイル側の粒子が先に磁場により速度方向が変えられ-Z 方向の速度を持つ粒子が多くなることが考えられる．その後プラズマは膨張していきコイル側以外の粒子も磁場に押し返され+Z 方向の速度を持つ粒子が多くなることで推力を得る．

図 4.6

球状ターゲットにおける全超粒子の Z 方向の速度分布[(a)初期状態，(b)2.2 μs，(c)6.7 μs]
50

・Impulse bit の比較
3D hybrid code により Impulse bit は，球状ターゲットは 12 Ns，円錐ターゲットは 0.75
Ns となり，円錐ターゲットの Impulse bit は球状ターゲットの 6.3%であることがわかった．こ
こで ，超粒子の質量𝑚 及び超粒子数𝑛s よりプラズマの質量は，𝑀p = 𝑚𝑛s とする．また，
1

プラズマの平均的速度を𝒗𝑧 = √𝑛 ∑ 𝑣i2 とすると，Z 方向の運動エネルギー𝐸 は式(4.7)
𝑠

で表せる．

𝐸=

1
𝑀 𝒗 2
2 p 𝑧

(4.7)

また，式(3.1)と式(4.7)より Impulse bit 𝐼 は式(4.8)のように表せる．

𝐼 = √2𝐸𝑀p(4.8)

式(4.8)より，Impulse bit はプラズマのエネルギーを固定した場合，プラズマの質量√𝑀p に比例する．そのため，球状ターゲットと円錐ターゲットのプラズマの質量はそれぞれ 100 mg，0.5 mg であるので，√𝑀p はそれぞれ 10 mg1/2，0.7 mg1/2 となり，円錐ターゲッ
トの√𝑀p は球状ターゲットの 7%である．そのため，3D hybrid code で計算した両者のImpulse bit の違いは妥当であると考えられる．
また hybrid 計算において，円錐ターゲットのプラズマの運動エネルギー及び質量は，
球状ターゲットに比べそれぞれ 1.25 倍，0.417 倍であった．そのため，式(4.8)より円錐ターゲットの Impulse bit は球状ターゲットの 6.4%と計算できるので，本シミュレーションで求めた Impulse bit は妥当であると考えられる．

51

・円錐ターゲット
円錐ターゲットの円錐推進剤部分の中性子減衰率は 63%であり，円錐推進剤中で約4 割の中性子束の減衰が見込め，それに伴い中性子遮蔽体及び排熱装置の重量の軽減も見込める．そのため，円錐推進剤を中性子遮蔽材とする分には有効であると考えられる．しかし，Impulse bit が球状ターゲットに比べ円錐ターゲットは 1/16 であるため，その分火星往復のためには，燃料や推進剤等の重量が増加する．そのため，推力と中性子遮蔽を同時に見込める推進剤形状に最適化する必要がある．

52

第5章

結論

・輻射流体(Star2D)とハイブリッド計算を組み合わせることで，核融合プラズマと推進剤プラズマの相互作用，加熱・膨張過程を検証するカップリングシミュレーションを作成できた．

・二次元輻射流体コードにおいて，球状ターゲットは高温の核融合プラズマが推進剤を加熱，プラズマ化し，膨張するのを確認できた．

・二次元輻射流体コードにおいて，円錐ターゲットの円錐推進剤は先端部分のみが加熱されることがわかった．円錐推進剤部分の 78%は加熱されず，カップリングタイミング時の円錐推進剤部分の合計エネルギーは全体の 0.33%だった．そのため，円錐ターゲッ
トの円錐推進剤は推力にほとんど寄与しない．

・3D hybrid code において，二次元輻射流体コードで計算した球状及び円錐ターゲットのデータを初期状態とした膨張プラズマが，コイルによる磁場と相互作用し，推力が得られることが確認できた．

・球状ターゲットの Impulse bit は，12 Ns で，運動量変換効率は 47%に対し，円錐ターゲ
ットの Impulse bit は 0.75 Ns で，運動量変換効率は 55%だった．

・円錐ターゲットの推進剤中で中性子束はおよそ 63%に減衰でき，中性子遮蔽体及び排
熱装置の重量軽減が見込めるが，球状ターゲットに比べ，Impulse bit が 1/16 であるため，推力が見込めない．そのため，中性子遮蔽と推力を同時に見込める推進剤形状に最適化する必要がある．

53

参考文献
(1) Hyde, R. A., Wood, L. L., and Nuckolls, J. H., “Prospects for Rocket Propulsion with
Laser-Induced Fusion Microexplosions,” AIAA Paper 72-1063, 1972.
(2) Orth, C. D., “VISTA - A Vehicle for Interplanetary Space Transport Application
Powered by Inertial Confinement Fusion,” Lawrence Livermore National Laboratory
UCRL-LR-110500, 2003.
(3) Hyde, R. A., “A Laser Fusion Rocket for Interplanetary Propulsion,” Lawrence
Livermore National Laboratory UCRL-88857, 1983.
(4) Nagamine, Y., and Nakashima, H., “Analysis of Plasma Behavior in a Magnetic Thrust
Chamber of a Laser Fusion Rocket,” Fusion Science and Technology, Vol. 35, Issue 1,
1999, pp. 62-70.
(5) Vchivokv, K. V., Nakashima, H., Ichikawa, F., Zakharov, Y. P., Vacuum, 73, 2003, 427.
(6) Vchivokv, K. V., Nakashima, H., Zakharov, Y. P., Esaki, T., Kawano, T., andMuranaka, T., “Laser-Produced Plasma Experiments and Particle in Cell Simulation to
Study Thrust Conversion Processes in a Laser Fusion Rocket,” Jpn. J. Appl. Phys., Vol.
42, Issue 10, 2003, pp. 6590-6597. doi: 10.1143/JJAP.42.6590(7) Şahin, S., Şahin, M. H., “Evaluation of the Neutron and Gamma-Ray Heating in the
Radiation Shielding and Magnet Coils of the VISTA Spacecraft,” Fusion Science andTechnology, Vol. 33, Issue 4, 1998, pp. 418-434.
(8) Şahin, S., Şahin, M. H., “Optimization of the radiation shielding mass for magnet coilsof the VISTA spacecraft,” Annlas of Nuclear Energy 28, 2001, pp.1413-1429.
(9) W.W. ENGLE, Jr., “ANISN, A One-Dimensional Discrete Ordinates Transport Codewith Anisotropic Scattering,” ORNL-K1693, Oak Ridge National Laboratory ,1970.
(10) “TORT-DORT, Two- and Three-Dimensional Discrete Ordinates Transport,” CCC543, Radiation Shielding Information Center, Oak Ridge National Laboratory, 1993.
(11) Sunahara, A., Nishihara, K, Sasaki, A., “Optimization of extreme ultraviolet emission
from laser-produced tin plasmas based on radiation hydrodynamics simulations,”
Plasma and Fusion Research, 3:043-043, 2008.

54

(12) Sunahara, A., Johzaki, T., Nagatomo, H., Mima, K., “Generation of pre-formed plasma
and its reduction for fast-ignition,” Laser and particle beams, 30(01):95-102, 2012.
(13) Spitzer Jr., L., and Härm, R., “Transport Phenomena in a Completely Ionized Gas,”
Phys. Rev., Vol. 89, Issue 5, 1953, pp. 977-981. doi: 10.1103/PhysRev.89.977(14) Spitzer Jr., L., Physics of fully ionized gases, second ed., pp. 30–40, Interscience
Publications, New York, 1962.
(15) Kruer. W. L., “The Physics of Laser Plasma Iteractions,” Addison-Wesley Pub, 1988,
p48(16) Malone, R. C., McCrory, R. L., and Morse, R. L., “Indications of Strongly FluxLimited Electron Thermal Conduction in Laser-Target Experiments,” Phys. Rev. Lett.,Vol. 34, 1975, pp. 721-724. doi: 10.1103/PhysRevLett.34.721(17) Harned, D. S., “Quasineutral Hybrid Simulation of Macroscopic Plasma Phenomena,”
J. Comput. Phys., Vol. 47, Issue 3, 1982, pp. 452-462. doi: 10.1016/00219991(82)90094-8(18) Horowitz, E. J., Shumaker, D. E., and Anderson, D. V., “QN3D: A ThreeDimensional Quasi-neutral Hybrid Particle-in-Cell Code with Applications to the Tilt
Mode Instability in Field Reversed Configurations,” J. Comput. Phys., Vol. 84, Issue
2, 1989, pp. 279-310. doi: 10.1016/0021-9991(89)90234-9(19) Mankofsky, A., Sudan, R. N., and Denavit, J., “Hybrid Simulation of Ion Beams in
Background Plasma,” J. Comput. Phys., Vol. 70, Issue 1, 1987, pp. 89-116. doi:
10.1016/0021-9991(87)90003-9(20) Birdsall, C. K., and Langdon, A. B., “Plasma Physics via Computer Simulation,”Francis & Taylor, New York, 2004, Appendix C.
(21) Christiansen, J. P., Ashby, D. E. T. F., and Roberts, K. V., “MEDUSA, A OneDimensional Laser Fusion Code,” Computer Physics Communication., 7, 1974, pp.271287.
(22)高村秀一, 「プラズマ理工学入門」, 森北出版株式会社, 1997, pp.20.
(23)

Cross

sectiontable,JENDL-4.0

https://wwwndc.jaea.go.jp/cgi-

bin/Tab80WWW.cgi?lib=J40&iso=H001 (最終閲覧日：2020 年 1 月 29 日)

55

謝辞
本研究を進めるにあたり，研究の場を提供していただき，あらゆることをご指導いただ
いた山本直嗣教授に深く感謝致します．プログラミングの基礎から研究内容の細部に至るまでたくさんの助言を下さった森田太智助教授に深く感謝致します．親身に研究のアドバイスをたくさん下さった中島秀紀教授に深く感謝致します．快適な研究環境を提供してくださった大内田映子秘書に深く感謝致します．
本研究にて使用させてもらい，また多くの助言など頂いた，輻射流体コードに携わっていらっしゃる砂原淳先生，城崎知至先生に深く感謝致します．
いつも優しく LFR の魅力や研究内容に助言を下さった枝本雅史氏に深く感謝致します．いつも親身になって多くのアドバイスを下さった児島富彦氏に深く感謝致します．研
究や進路についての助言をくださった諸先輩方，講義や研究面で苦楽を共にした黒江
雄貴氏，高木麻理子氏，安達美咲氏，西岡祐輝氏，齊木淳氏，西牟田直氏に深く感謝いたします．また修士 1 年，学部 4 年生の皆さんに心から感謝申し上げます．

最後に，今まで支えてくれた友人，家族に深く感謝致します．

56