令 和 二 年 度 修 士 論 文
非断熱トラップを用いた
D-3He 核融合炉心プラズマに
おける
14.7MeV 陽子閉じ込めと構造材への影響
指導教員 高橋 俊樹 教授
群馬大学大学院理工学府 理工学専攻
電子情報・数理教育プログラム
角掛 東
目次
第
1 章 序論 ... 1
1.1 研究背景 ... 1 1.2 非断熱トラップ... 5 1.3 研究目的 ... 7第
2 章 非断熱トラップの磁場構造 ... 8
2.1 Accessible region ... 8 2.2 ヘルムホルツコイル半径 ... 8 2.3 装置モデル... 10 2.4 磁場計算 ... 11 2.5 磁場構造 ... 11 2.6 核融合反応確率... 12 2.7 Simpson 則... 13 2.8 プラズマの平衡計算 ... 14 2.9 密度計算 ... 16 2.10 核融合発生率分布 ... 16第
3 章 陽子軌道計算 ... 17
3.1 軌道計算パラメータ ... 17 3.2 速度の等方性... 18 3.3 陽子軌道計算に用いる運動方程式 ... 19 3.4 Slowing-Down collision ... 19 3.5 Runge-Kutta 法 ... 20第
4 章 支持構造物命中粒子 ... 22
4.1 支持構造物に命中した粒子割合 ... 22第
5 章 ローレンツ力により生じる支持構造物へ加わる応力解析 ... 25
5.1 ヘルムホルツコイルに加わるローレンツ力 ... 25 5.2 Von Mises 応力 ... 26 5.3 CAE 解析 ... 27第
6 章 まとめと今後の課題 ... 30
6.1 まとめ ... 30 6.2 今後の課題... 30謝辞 ... 31
参考文献 ... 32
1
第
1 章 序論
1.1 研究背景
現在，人口の爆発的な増加などからエネルギーについての問題が深刻化している．人類 の多くは石油や石炭などの化石燃料に依存しているものの，しかしながらこれらは枯渇や 二酸化炭素による温暖化といった問題がつきまとい，技術の発達に伴って更に資源の需要 増加が見込まれる近い未来において，化石燃料を始めとしたエネルギー資源の追求は避け られない課題となる． その背景から，エネルギーの生成方法としては数多くの研究がなされた．特に水力や火 力，原子力発電が2019 年 3 月現在の主な方式であるが，水力には立地条件上の限界が存在 し，発電量の増大には広大な面積を発電のために用いなければならない．一方火力発電は 石油や石炭などの化石燃料を燃焼させるために化石燃料の枯渇後に使用できず，また二酸 化炭素などの温室効果ガス排出量も増加させる．原子力発電はこれらの温室効果ガスの心 配はないものの，しかしながらウランなど化石燃料の問題に加え，東日本大震災の福島第 一原発事故以後は放射能汚染の影響が危惧されている． 風力や太陽光発電などではこれらを代替できるほどのエネルギー生成能力はなく，水力 以上に膨大な設置面積を必要とするため，化石燃料枯渇後の発電を考慮すると大出力かつ 放射能汚染の恐れが少ないエネルギー取得方法が求められている． 以上の条件をある程度満たす候補として，核融合発電が挙げられる．この発電方式は原 子の核融合反応によって生じた質量の欠損からエネルギーを得る方式で，現在ではフラン スのカダラッシュにおいてトカマク式核融合炉の一つである国際熱核融合実験炉(ITER)計 画が進められており，実用化に向けた動きが行われている． Fig. 1-1. 国際熱核融合実験炉（ITER）模式図[1]．2 しかし，核融合反応の実用化に関しては未だ多くの課題が残されており，ITER で想定さ れているD-T 反応などは以下の Table 1 に示す通り 14.1 MeV の高速中性子を放出するが， この中性子は放射化の原因となるため核融合炉の放射化に対する危惧が生じ，クリーンな エネルギー源を求める場合には必ずしも適当ではない． Table 1-1. 主な核融合反応の例． D – D 50% → T(1.01 MeV) + p(3.02 MeV)
50% → 3He(0.82 MeV) + n(2.45 MeV)
D – T → He(3.5 MeV) + n(14.1 MeV) D – 3He → He(3.6 MeV) + p(14.7 MeV)
そのため本研究ではD-3He 反応について考える．この反応は D-T 反応と比べ放射化の原 因となる中性子が発生しないため，放射能汚染の原因が少なく，また陽子が荷電粒子であ ることから直接的にエネルギーを得ることができる可能性が生じる．実際にはTable 1-1 か らD-D 反応が発生するため放射化の危惧はあるものの D-T 反応に比べ少ない．
しかしD-3He 反応について，Fig. 1-2 に示すように D-T 反応などと比べると要求される温 度が非常に高い． Fig. 1-2. 核融合の反応率．
3 これは従来と同様の方法を取って反応を行った場合，シンクロトロン放射損失が発生す ることでエネルギー効率が悪化してしまい発電の目的では不適格となることを示している． そのため，プラズマの閉じ込め効率の検討が行われてきた．このプラズマ閉じ込め効率を 表す値として以下に示すベータ値が存在し，D-3He 反応による発電には D-T を目的とした炉 よりもベータ値の高いプラズマ閉じ込め方式が必要となる．ベータ値は次の式で表される． 𝛽 = プラズマ圧 外部閉じ込め磁気圧 (1.1) この式はプラズマ圧力と外部磁場による磁気圧の比で，1 に近いほど効率的なプラズマ閉 じ込めである．この値はITER を始めとしたトカマク式核融合炉では，およそ 0.01 オーダー である．また日本の核融合研究所にあるLHD 装置はヘリカル方式を取っているが，これは 目標としても0.041 である[2]．しかし，これらはD-3He 反応の実用化に足るほどの値ではな く，D-3He 核融合炉の実用化を目指す場合はより高ベータが望まれる． そこで，1 に近い高ベータ値を実現できる方式の一つとして，磁気ミラー方式が挙げられ る．磁気ミラー方式はFig. 1-3 のように作られる磁場構造であり，(1.2)式で表される断熱不 変量である磁気モーメントにより，弱磁場領域と強磁場領域において(1.3)式のエネルギー保 存により装置端の磁場を十分に強くし𝑣∥= 0となれば粒子は反射する． 𝜇 ≡1 2 𝑚𝑣⊥2 𝐵 = const. (1.2) d d𝑡( 1 2𝑚𝑣⊥ 2+1 2𝑚𝑣∥ 2) (1.3) Fig. 1-3. 磁気ミラー[2]． しかしこの方式では閉じ込められない粒子も存在する．磁場 B0領域内で𝑣⊥= 𝑣⊥0，𝑣∥= 𝑣∥0である粒子を考えると，この粒子がBmの地点で反射されるとすると，速度を𝑣⊥= 𝑣⊥𝑚， 𝑣∥= 0とすれば，磁気モーメントの保存から 1 2 𝑚𝑣⊥02 𝐵0 =1 2 𝑚𝑣⊥m2 𝐵m (1.4) ここで，(1.3)式より
4 𝑣⊥m2= 𝑣⊥02+ 𝑣∥02≡ 𝑣02 (1.5) であるため， 𝐵0 𝐵m =𝑣⊥0 2 𝑣02 ≡ sin2𝜃 (1.6) とできる．このときの𝜃は粒子が閉じ込められる最小の角度であり，これ以下の角度を持 つ粒子は閉じ込めることができない．これをロスコーンといい，これにより磁気ミラー方 式は安定性に欠く．
5
1.2 非断熱トラップ
ロスコーンを抑制する方式として，百田らの提唱する非断熱トラップがある[3]．これは (1.2)式によって示される磁気モーメントの保存を磁場が 0 となる領域を作ることによって 破り，粒子挙動を確率的なものとすることでロスコーンを抑制し安定したプラズマの閉じ 込めが期待でき，副次的に極小磁場となるため磁気ミラー以上の高ベータが期待できる． 磁気モーメントが破れ跳躍する例をFig. 1-3 に示す．確率的な挙動を行う粒子は開放端を持 つ装置内においてFig. 1-4 のように，左右のいずれかに等確率で移動する．このとき粒子の 閉じ込め時間は装置の連結数の自乗に比例すると見積もられており，効率の良い閉じ込め を実現できると予想されている[4]． また，装置が開放端を持つため，Fig.1-5 に示すように陽子を直接エネルギーに変換する エネルギー変換器と装置との接続も可能である． Fig. 1-3. 磁気モーメント跳躍[5]． Fig. 1-4. 装置内の確率的挙動[4]．6 Fig. 1-5. 装置とエネルギー変換器の接続． 非断熱トラップを形成するには，磁場が0 となる領域を作らなければならない．そこで， Fig.1-6 に示されるような装置を用いる．この装置はソレノイドコイルとミラーコイルによ って作られた装置内の磁場をヘルムホルツコイルによってキャンセルすることにより極端 な弱磁場領域を生成する事ができる． しかし，ヘルムホルツコイルが磁場中に存在するため，ヘルムホルツコイル間にはロー レンツ力が働く．また，核融合により発生した14.7 MeV のエネルギーを持つ陽子が支持構 造物に命中する．そのため，これらにより支持構造物が破断されると，装置は非断熱トラ ップを維持することができず，またヘルムホルツコイルが衝突することが予想される． Fig. 1-6 非断熱トラップ生成を考慮した装置． また，ソレノイドコイル並びにミラーコイルについては14.7 MeV の陽子を閉じ込めるた め大電流が必要となり，伴ってヘルムホルツコイルに加わる電流も大きい値でなければ磁 場を打ち消すことができない．これより常伝導帯を用いて装置を作成した場合，抵抗によ る電力損失が膨大となり発電を目的とする装置としてふさわしくない．そのため，コイル には超伝導体を扱う必要がある．
7
1.3 研究目的
以上の要因から，非断熱トラップを形作る核融合炉の実用化のためには，ヘルムホルツ コイル支持構造物に対して加わるローレンツ力並びに衝突する陽子が与える影響について 考慮し，支持構造物に要求される強度について求める必要がある． 本研究では支持構造物に要求される強度をローレンツ力やヘルムホルツコイルの臨界電 流密度から計算して支持構造物が必要な太さを決定し，粒子軌道計算を行って支持構造物 にどれほど陽子が命中したとき破断されるかを導出することを目的とする．8
第
2 章 非断熱トラップの磁場構造
2.1 Accessible region
装置の詳細な概要を決定するにあたり，まずはソレノイドコイルとミラーコイルについ てのパラメータをそれぞれ求める必要があるが，この時考慮するべき事項は装置内に粒子 が閉じ込めることができるかどうかである．そこでまずは粒子が動くことのできる範囲に ついて考える．定常かつ軸対称な系における単一粒子の全エネルギー𝐻は 𝐻 =1 2𝑚𝑣 2+ 𝑞𝜙 (2.1) ここで，電場は静電場であるため， 𝐸 = −∇𝜙 (2.2) また，正準角運動量は 𝑃θ= 𝑚𝑣θ𝑟 + 𝑞𝜓 (2.3) ここで，静電ポテンシャル𝜙 = 0とすると 𝐻 =1 2𝑚𝑣 2 =1 2𝑚(𝑣r 2+ 𝑣 z2+ 𝑣θ2) ≥ 1 2𝑚𝑣θ 2 (2.4) (5.10)式を変形すれば 𝑚𝑣θ𝑟 = 𝑃θ− 𝑞𝜓 (2.5) これより， 1 2𝑚𝑣θ 2= 1 2𝑚𝑟2(𝑚𝑣θ𝑟)2 = 1 2𝑚𝑟2(𝑃θ− 𝑞𝜓)2 (2.6) (5.11)式と(5.13)式から 𝐻 ≥ 1 2𝑚𝑟2(𝑃θ− 𝑞𝜓) 2= 𝑈(𝑟, 𝑧; 𝑃 θ) (2.7) この𝑈は擬似的なポテンシャルと見なすことができ，Störmer potential と呼ばれる．(2.7) 式は𝑣2⁄2𝑚のエネルギーを持つ粒子が運動可能な領域を表す式となり，この運動可能な領域をAccessible region（または Störmer region）と呼ぶ．
2.2 ヘルムホルツコイル半径
次にヘルムホルツコイルの半径について決定する方法を考える．当研究室の過去研究[4] より，ある一定以上の半径を取るとプラズマが装置壁と接触してしまうことがわかってい る。プラズマが壁と接触した場合，装置壁との相互作用について危惧が生じてしまうため， 相互作用が起きない程度の半径を設定する必要がある．そこで，Fig. 2-1 にヘルムホルツコ イル半径の変化に対する閉じ込め熱エネルギー量を当該の過去論文より抜粋し示す．9 Fig. 2-1. ヘルムホルツコイル半径の変化に対する閉じ込め熱エネルギー量[4] Fig.2-1 より，閉じ込め熱エネルギー量はヘルムホルツコイル半径が 1.4 [m] 程で一度ピー クを持ち，1.8 [m] になると再度急激に増加することがわかる，しかし，1.8 [m] 以降のもの はそれ以前とは明らかに状態の異なる変化であり，これはプラズマが壁と接触しているた めに閉じ込められる熱エネルギー量が多くなったことが考えられる．そのため，ヘルムホ ルツコイル半径は1.4 [m] 付近が好ましいと考えることができ，本実験ではヘルムホルツコ イル半径 𝑟a＝1.4 [m] として以降のシミュレーションを行う．
10
2.3 装置モデル
以上に示したAccessible region が装置内に核融合生成粒子を閉じ込めることができるよう に設定し，ヘルムホルツコイル半径 𝑟a＝1.4 [m]として設定した計算モデルから作成した， 非断熱トラップにおけるコイル配置と流れる電流値についてFig. 2-2 に示す．本研究ではシ ミュレーションに使用されているパラメータは全て以下の値に固定されている．また，ヘ ルムホルツコイルに流れる電流値についてはソレノイドコイルとミラーコイルが装置中央 に作る磁場を打ち消すように設定し，次の式で与えられる． 𝐼h= ( 5 4) 3 2 𝑟a[{ 1 2𝑛s𝐼s( 𝑙s 𝑙s2+ 𝑟s2 )} + {𝑛m𝐼m( 1 2 𝑙s+ 𝑙m (12 𝑙s+ 𝑙m) 2 + 𝑟m2 − 𝑙s 𝑙s2+ 𝑟m2 )}] (2.8) Fig. 2-2. 装置モデル． また，ヘルムホルツコイルの支持構造物については本数を四本かつ高さを0.6[m]，底面積 を0.15×0.15[m2]とし，間隔が一定となるように配置する．考慮する素材についてはデータ のあるもので最も適したと思われるものを採用し，ヘルムホルツコイルについては超伝導 体YBCO SuperPower SCS を，支持構造物については純タングステンを考慮し，シミュレー ションはFortran f90 にて行う．11
2.4 磁場計算
このコイル配置による磁場の生成について，磁場計算を行う．磁場計算は軸対称として， 各コイル（コイル半径𝑎，コイル電流𝐼）によるベクトルポテンシャルを 𝐴φ(𝑟, 𝑧) = 𝜇0𝐼 4𝜋∫ 𝑎 cos 𝜒 d𝜒 √𝑧2+ 𝑟2+ 𝑎2− 2𝑎𝑟 cos 𝜒 2𝜋 0 (2.9) で求め，磁束関数 𝜓(𝑟, 𝑧) = 𝑟𝐴𝜑(𝑟, 𝑧) (2.10) から磁場を求めると，次のように磁場𝐁が求められる． 𝐵𝑟(𝑟, 𝑧) = − 1 𝑟 𝜕𝜓(𝑟, 𝑧) 𝜕𝑧 (2.11) 𝐵𝑧(𝑟, 𝑧) = 1 𝑟 𝜕𝜓(𝑟, 𝑧) 𝜕𝑟 (2.12) 数値計算においては規格化をした上で，(2.9)式は台形公式により導出し，(2.11)並びに(2.12) 式は4 次精度の有限差分法によって計算した．2.5 磁場構造
計算により得られた磁場構造をFig. 2-3 に示す．示された磁場は無限長ソレノイドコイル の作る磁場を基準として規格化されたもので，r, z 座標が 0 に近い付近において値が 0 に近 い弱磁場領域が形成されていることが見て取れ，これにより非断熱トラップが形成される ことがわかる． Fig. 2-3. 規格化磁場構造．12
2.6 核融合反応率
磁場を求めることができたため、この装置内で発生した核融合の単位面積毎秒当たりに 発生する核融合反応の回数である核融合反応率について求めることができる。 核融合反応率𝑃 [m−3s−1]は次の式によって示される． 𝑃 = 𝑛1𝑛2〈𝜎𝑣〉 (2.13) 𝑛1, 𝑛2は反応する粒子の密度 [m-3]，〈𝜎𝑣〉は核反応率 [m3s-1] であり，これらを求める必要 がある． ここで，核反応率〈𝜎𝑣〉は核反応断面積𝜎 [barn] とMaxwell-Boltzmann 分布により求められ， D-3He 反応では，以下の(2.14)式によって核反応断面積が示される[6]． 𝜎 =647 + 25900{(1.297 − 3.98 × 10 −3𝐸)2+ 1}−1 𝐸{exp(89.27𝐸−1 2⁄ ) − 1} (2.14) また，Maxwell-Boltzmann 分布により核反応率〈𝜎𝑣〉は 〈𝜎v〉 = ( 8 𝜋𝜇12(𝑘𝑇)3 ) −12 ∫ 𝜎(𝐸)𝐸exp(− 𝐸 𝑘𝑇) (2.15) 𝐸はエネルギー[keV]，𝜇12はD と3He の換算質量 [kg], 𝑘𝑇は温度[keV]である． (2.15)式を解くためには，(2.14)式を代入した上で，積分区間を設定し，積分計算を行わな ければならない．そこで，今回はSimpson 則により計算を行う． 10-6 10-5 0.0001 0.001 0.01 0.1 1 1 10 100 1000 σ [b ar n ] Energy[keV] Fig. 2-4. D-3He 反応の反応断面積．13
2.7 Simpson 則
Simpson 則は定積分の近似解を求める場合に用いられ，同じく定積分を求める台形公式な どよりも誤差が小さい．この公式は ∫ 𝑓(𝑥)d𝑥 𝑏 𝑎 (2.16) となる定積分において，Fig. 2-5 のように等間隔の x 座標を持つ三点の関数値を 𝑦(𝑥) = 𝑎(𝑥 − 𝑥2)2+ 𝑏(𝑥 − 𝑥2) + 𝑐 (2.17) で補完することを考える．ここで，各点間の幅ℎから 𝑦1 = 𝑎ℎ2− 𝑏ℎ + 𝑐, 𝑦2 = 𝑐, 𝑦3= 𝑎ℎ2+ 𝑏ℎ + 𝑐 (2.18) であるため， a =𝑦1+ 𝑦3− 2𝑦2 2ℎ2 , 𝑏 = 𝑦3− 𝑦1 2ℎ , 𝑐 = 𝑦2 (2.19) 以上より，被積分関数𝑓(𝑥)は𝑥3から𝑥1までの範囲で ∫ 𝑓(𝑥)d𝑥 𝑥3 𝑥1 ≅ ∫ y(𝑥)d𝑥 𝑥3 𝑥1 =ℎ 3(𝑦1+ 4𝑦2+ 𝑦3) (2.20) (2.20)式を任意の範囲 a から b までに一般化すると ∫ 𝑓(𝑥)d𝑥 𝑏 𝑎 ≅ℎ 3∑[𝑓(𝑥2𝑗−1) + 𝑓(𝑏𝑥2𝑗+1) + 𝑓(𝑏𝑥2𝑗+2)] 𝑁−1 𝑗=0 =ℎ 3[𝑓(𝑎) + 𝑓(𝑏) + 4 ∑ 𝑓(𝑥2𝑗−1) 𝑁 𝑗=1 + 2 ∑ 𝑓(𝑥2𝑗) 𝑁−1 𝑗=1 ] (2.21) 以上の式を用いることで，シミュレーション上で積分計算を行うことができる．14 Fig. 2-5. Simpson 則の導出．
2.8 プラズマの平衡計算
以上の計算によって核反応率が求められたため，次にD, 3He の密度について求める．こ れらを求めるためには，プラズマを電磁流体的に近似し，圧力平衡状態を求めることがで きる．電磁流体方程式は 𝜌m{ 𝜕𝐕 𝜕𝑡 + (𝐕・∇)𝐕} = −∇𝑝 + 𝐣 × 𝐁 (2.22) ただし，𝜌mは電荷密度，𝐕は流速ベクトルである．𝐕 = 𝟎 ならびに d𝐕 d𝑡⁄ = 0 の場合， ∇𝑝 = 𝐣 × 𝐁 (2.23) また，定常状態となるため ∇ × 𝐁 = 𝜇0𝐣 (2.24) ∇・𝐁 = 0 (2.25) ∇・𝐣 = 0 (2.26) これより，(2.23)式から 𝐁・∇𝑝 = 0 (2.27) となり，これは等圧面と磁気面の一致を示す． ここで，(2.4), (2.5)式により，磁場の各成分と(2-27)式から −𝜕𝜓 𝜕𝑧 𝜕𝑝 𝜕𝑟+ 𝜕𝜓 𝜕𝑟 𝜕𝑝 𝜕𝑧 = 0 (2.28) これを満たす𝑝を考えた場合，𝑝は磁束関数のみに依存する事がわかる．つまり 𝑝 = 𝑝(𝜓) (2.29) また，(2.24)式と軸対称性を考えると，15 𝜇0𝑗θ= 𝜕 𝜕𝑧(− 1 𝑟 𝜕𝜓 𝜕𝑧) − 𝜕 𝜕𝑟( 1 𝑟 𝜕𝜓 𝜕𝑟) = − 1 𝑟 𝜕2𝜓 𝜕𝑧2 − 𝜕 𝜕𝑟( 1 𝑟 𝜕𝜓 𝜕𝑟) = −𝑟𝜕𝜓 𝜕𝑟( 1 𝑟 𝜕𝜓 𝜕𝑟) − 𝜕2𝜓 𝜕𝑧2 (2.30) (2.23)式から 𝜕𝑝 𝜕𝑟 = 𝑗θ𝐵𝑧 (2.31) 𝜕𝑝 𝜕𝑧 = −𝑗θ𝐵𝑟 (2.32) (2.5)式並びに(2.30)式より 𝑗θ= 𝑟 d𝑝(𝜓) d𝜓 (2.33) 以上をまとめると， 𝑟 𝜕 𝜕𝑟( 1 𝑟 𝜕𝜓 𝜕𝑧) + 𝜕2𝜓 𝜕𝑧2 = −𝜇0𝑟2 d𝑝(𝜓) d𝜓 (2.34) (2.34)式は Grad-Shafranov 方程式と呼ばれ，圧力勾配を与えて解くことで圧力平衡を求め， 分布を導くことができる．
16
2.9 密度計算
Grad-Shafranov 方程式に，以下で示された圧力勾配を与え，圧力平衡を求めた． d𝑝(𝜓) d𝜓 = 2𝑎𝜓 (6.23) ここで，𝑎 = −100とし，反復法を用いて Grad- Shafranov 方程式を解き，装置内のプラズ マ圧力分布を求めた．プラズマ圧力は密度𝑛 [m-3] と温度𝑘𝑇 [keV] により次の(6.24)式のよ うに表すことができるため，これにより密度を求める事ができる． 𝑝 = 𝑛𝑘𝑇 (6.24) 今回は計算の簡単化のため，温度を一様に100 [keV] として装置全体に与え，D と3He の 密度比は1 となるように考える．2.10 核融合発生率分布
以上の結果を元に導いた核融合発生率分布をFig. 2-6 に示す． Fig. 2-6. 核融合発生率分布． 核融合発生率は最大値が1 となるよう規格化され，規格化するのに用いた最大の反応率𝑃0 は𝑃0 = 2.41 × 1023[m−3s−1] である．17
第
3 章 陽子軌道計算
3.1 軌道計算パラメータ
シミュレーションに用いた軌道計算のパラメータを，Table 3-1 に示す．装置や磁場構造 については第2 章に示したものと同様のものを用いる． Table 3-1. 軌道計算モデル． パラメータ 粒子数 [個] 100000 初速度 [m/s] 等方的に与える ステップ数 100000 与えたステップ数は，実時間にして0.1[us]と等価である．また，今回は粒子の初期座標 についてFig.2-6 に示した核融合反応率に応じて設定する．装置は軸対象であり，体積効果 を考慮した場合には粒子の初期配置は以下のFig.3-1 に示すような形を取る． Fig. 3-1. 粒子の初期配置． 14.7 [MeV] の陽子について，速度をエネルギーが全て運動エネルギーに変換されている として古典論から求めた．また，等方的な速度については後述する． ここで，𝐸𝑘は運動エネルギー[J]，𝑚pは陽子質量[kg]である． 数値計算においてはr, z 座標について 1 で規格化し，微小時間はサイクロトロン周波数の 逆数から𝑚 𝑒𝐵⁄ 0で規格化，速度はz 座標の最大値を規格化微小時間で割ることで規格化して いる．ただし，𝐵0は磁場の規格化定数．18
3.2 速度の等方性
Table 3-1 に示すように，軌道計算を行う粒子の初期速度についてはこれらの分布などが 方向に依存しないよう等方的に与えている． 粒子の初期速度に等方的な分布を与えるには単位球面上に一様に分布する(x,y,z)の組を 考える必要がある．それぞれ[-1,1], [0,2π]の範囲で一様に分布する確率変数 X1,X2を用いる と，等方的な分布を与えるための値は次のように示される[7]． 𝜃 = arccos(−𝑋1) (3.1) 𝜑 = 𝑋2 (3.2) 𝑥 = √1 − 𝑋12cos𝑋2 (3.3) 𝑦 = √1 − 𝑋12sin𝑋2 (3.4) 𝑧 = 𝑋1 (3.5) この値を用いて粒子に初期速度を与える．粒子の初期速度の大きさ|𝐯𝟎|を用いると初期 速度𝑣𝑟0, 𝑣𝜃0, 𝑣𝑧0は 𝑣𝑟0= √𝑥2+ 𝑦2|𝐯𝟎| (3.6) 𝑣𝜃0= 𝜃|𝐯𝟎| (3.7) 𝑣𝑧0= 𝑧|𝐯𝟎| (3.8) と与えられる．これらの式によって与えられる初期速度の分布はFig.3-2 に示すように 球状となる． (a)x-y 平面分布 (b)y-z 平面分布 (c)z-x 平面分布 Fig. 3-2. 粒子の初期速度分布．19
3.3 陽子軌道計算に用いる運動方程式
陽子の軌道計算には運動方程式 𝑚α d𝐯α d𝑡 = 𝑞α(𝐄 + 𝐯α× 𝐁) − ∑ 𝑚α𝜐αβ(𝐯α− 𝐮β) 𝛽 (3.9) を用いて計算を行う．ここで，αはテスト粒子で，βはバックグランド粒子を示している． 𝜐αβはテスト粒子とバックグランド粒子との衝突周波数であり， 𝜐αβ = ∑ 𝑛β𝑞α2𝑞β2𝜆 4𝜋𝜀02𝑚α2 (1 +𝑚α 𝑚β ) 𝛽 𝜂(𝑥) 𝐯α3 (3.10) 𝜂(𝑥) = 2 √𝜋∫ 𝑒 𝑡 √𝑡d𝑡 𝑥 0 (3.11) により与えられる．本研究ではテスト粒子に核融合により発生した陽子を，バックグラ ウンド粒子に電子を用いている．3.4 Slowing-Down collision
(3.9)式の第二項は，プラズマ中の粒子が運動中にバックグラウンド粒子とのクーロン衝突 により速度が変化することを表しており，この摩擦による速度の減速現象をSlowing-Down Collision という． Slowing-Down Collision の衝突周波数は， 𝜐sl= ∑ 𝑛S𝑒2𝑒S2𝜆 4𝜋𝜀02𝑚2 (1 + 𝑚 𝑚S ) 𝑆 𝜂(𝑥) 𝑣3 = ∑𝑛S𝑒 2𝑒 S2𝜆 4𝜋𝜀02𝑚2 (1 + 𝑚 𝑚S ) 𝑆 ( 𝑚S 2𝑘𝐵𝑇𝑆 ) 3 2 ・𝜂(𝑥) 𝑥32 (3.12) 𝑥 ≡ 𝑚S𝑣 2 2𝑘𝐵𝑇𝑆 (3.13) で表される．この𝜐slは，温度𝑇𝑆は Maxwell 分布したプラズマ中に存在するテスト粒子の 速度であるe-folding time の逆数で，slowing-down collision frequency と呼ぶ．また，(3.12)式で与えられる slowing-down collision frequency は𝑣2 ≫ 𝟐𝑘𝐵𝑇𝑆⁄𝑚Sの場合は
𝜂(𝑥) → 1であり，反対に𝑣2 ≪ 𝟐𝑘 𝐵𝑇𝑆⁄𝑚Sの場合は𝜂(𝑥) ≅ 4 3√𝜋𝑥 3 2である．よって，(3.12)式は 2 つの極限を持ち， 𝜐sl= { ∑ 𝑛S𝑒 2𝑒 S2𝜆 8√2𝜋𝜀02𝑚2 (1 + 𝑚 𝑚S ) ・𝜀−32 𝑆 (∵ 𝑣2 ≫ 2𝑇 𝑆⁄𝑚S, 𝜀 = 1/2𝑚𝑣2) ∑𝑛S𝑒 2𝑒 S2𝜆 3𝜋𝜀02𝑚2 (1 + 𝑚 𝑚S ) 𝑆 (𝑚S 2𝑇S ) (∵ 𝑣2≪ 2𝑇 𝑆⁄𝑚S) (3.14) となる．すなわち，プラズマの熱速度(𝑣th=√2𝑇𝑆⁄𝑚S)より充分に大きな速度をもつテス ト粒子の𝜐slは，自身の速度の 3 乗に逆比例する．反対に熱速度より充分に小さいテスト粒
20 子は熱速度の3 乗に逆比例した一定値に近づく． ところで，(3.9)式の右辺第二項について，電子全体の流速は 0 であるため𝐮βは無視する ことができる．これより陽子の減速量について 𝑚p d𝐯p d𝑡 = −𝑚p𝜐pe𝐯p (3.15) となるが、この両辺に陽子の速度を掛けることで 𝑚p𝐯p・ d𝐯p d𝑡 = −𝑚p𝜐pe|𝐯p| 2 (3.16) 整理すると d d𝑡( 1 2𝑚p|𝐯p| 2 ) = −1 2𝑚p|𝐯p| 2 ・2𝜐pe (3.17) このとき，1 2𝑚p|𝐯p| 2 はエネルギーとなるため， d𝐸 d𝑡 = −2𝜐pe𝐸 (3.18) と表すことができる．(3.10)式を𝐸について解くとエネルギーの初期値𝐸0を用いて 𝐸 = 𝐸0exp(−2𝜐pe𝑡) (3.19) が得られ、𝑡秒後のエネルギーを求めることができる．ところで今回の計算時間は実時間にし ておおよそ0.1[us]に等しいため与えられたパラメータを用いて(3.10)式を解いた場合，Slowing-Down Collision の 衝 突 周 波 数 の 値 は お お よ そ𝜐pe =7.5 × 104 ほ ど が 得 ら れ ， こ の 値 と 𝐸0 = 14.7 [MeV],𝑡 = 0.1 × 10−6 [s]を(3.19)式に代入することで減衰後の陽子のエネルギーを求め ると𝐸 ≈ 14.7 [MeV]となり，本計算ではほとんど減衰していないことが分かる。
3.5
Runge-Kutta 法
(3.9)式の運動方程式をシミュレーションを用いて解くためには，一次の微分方程式につい て考慮する必要がある．そこで，本シミュレーションでは四次の Runge-Kutta 法を使用し， 結果を導出する． シミュレーションで用いられる微分方程式の近似的解法には，Euler 法と Runge-Kutta 法な どがある．Euler 法は精度に欠くため，本実験では四次の Runge-Kutta 法を用いる．Runge-Kutta 法は d𝑦 d𝑡 = 𝑓(𝑡, 𝑦) (3.20) で表される常微分方程式について，解が𝑦𝑖であるときから微小時間ℎ後の微分方程式の解 𝑦𝑖+1は 𝑦𝑖+1= 𝑦𝑖+ 1 6(𝑘0+ 2𝑘1+ 2𝑘2+ 𝑘3) (3.21)21 ここで 𝑘0= ℎ𝑓(𝑡𝑖, 𝑦𝑖) (3.22) 𝑘1 = ℎ𝑓(𝑡𝑖+ ℎ 2, 𝑦𝑖+ 𝑘0 2) (3.23) 𝑘2 = ℎ𝑓(𝑡𝑖+ ℎ 2, 𝑦𝑖+ 𝑘1 2) (3.24) 𝑘3= ℎ𝑓(𝑡𝑖+ ℎ, 𝑦𝑖+ 𝑘2) (3.25) これらはFig. 3-1 のように，微小時間幅の半分において傾き𝑘0を計算し，さらにこうして 導出した傾きによって𝑦𝑖+1 2 の値を仮定してから，それを初期値として同様に傾き𝑘1を求め， 𝑦𝑖に戻って再度同様の手順を取って𝑘2, 𝑘3を求める．最終的な微分方程式の解𝑦𝑖+1は𝑘1, 𝑘2 に重みを持たせた加重平均を計算して導出することができる． Fig. 3-3. Runge-Kutta 法の導出．
22
第
4 章 支持構造物命中粒子
4.1 支持構造物に命中した粒子割合
支持構造物が破壊されるかどうかを調べるためには，第 3 章で行った粒子の軌道計算の 結果から，計算粒子のうち支持構造物へ命中した粒子の割合を導く必要がある．ここで， 軌道計算の結果計算が終了した粒子の割合について以下のTable4-1 に示す． Table 4-1. 計算が終了した粒子割合． 装置左側に流出[%] 装置右側に流出[%] 支持構造物に命中[%] 合計[%] 13.9 14.2 4.08 32.2 この結果から，計算が終了した粒子のうちおよそ[%]が支持構造物へと命中することがわ かる．このとき，それぞれの支持構造物に命中する粒子の割合は Table4-2 に示すようにな る． Table 4-2. 支持構造物へ命中した粒子割合． 支持構造物の角度[rad] 左コイルへの命中割合[%] 右コイルへの命中割合[%] 𝜋 4⁄ 0.568 0.490 3𝜋 4⁄ 0.501 0.510 5𝜋 4⁄ 0.461 0.552 7𝜋 4⁄ 0.487 0.507 支持構造物への粒子の影響を論じるには，これらが支持構造物へと命中した箇所につい て調べる必要がある．粒子の命中箇所については正負の座標に配置された二本のヘルムホ ルツコイルについて支持構造物が各コイルに 4 本ずつ配置されており，それぞれ 4 つの面 に命中したものを考えなければならないため合計で32 パターンの画像を出すことができる が，今回は負の座標にあるコイルについて角度𝜋 4⁄ に位置する支持構造物の正面面へ命中し た粒子の命中分布を代表して図示する．23 Fig. 4-1. 粒子の支持構造物命中分布． 次に，これらの粒子が支持構造物へと加えるエネルギーを求める．支持構造物毎に加わる エネルギーの大きさ𝑊 [W]については次に示す式から求めることができる． 𝑊 = 𝑃𝐸𝑝∬ 𝑛D𝑛3He〈𝜎𝑣〉d𝑉 (4.1) 𝑃は命中の確率であり，𝐸𝑝は陽子のエネルギー[J]である．それぞれの値は以下の式により 与えられる． 𝑃 = (各支持構造物への命中粒子数) (総粒子数) (4.2) 𝐸𝑝= 1 2𝑚𝑝|𝐯𝑝| 2 (4.3) ところで陽子のエネルギーはほとんど減衰しないため，全ての陽子のエネルギーを 𝐸𝑝 ≈14.7 × 106⁄(1.60 × 10−19)として計算することができる．装置内全体の核融合反応率の 総和は∬ 𝑛D𝑛3He〈𝜎𝑣〉d𝑉 = 4.76 × 1027であるため，これらを(4.1)式に代入すると 𝑊 = 4.37 × 1015𝑃 (4.4) となる．これにより求められる支持構造物へ加わるエネルギー量をTable 4-3 に示す．
24 Table 4-3. 支持構造物へ加わるエネルギー． 支持構造物の角度[rad] 左コイル支持構造物へ 加わるエネルギー[W] 右コイル支持構造物へ 加わるエネルギー[W] 𝜋 4⁄ 2.48 × 1013 2.14 × 1013 3𝜋 4⁄ 2.19 × 1013 2.23 × 1013 5𝜋 4⁄ 2.01 × 1013 2.41 × 1013 7𝜋 4⁄ 2.13 × 1013 2.21 × 1013 以上より，支持構造物には0.1[us]に最大で一本あたり総粒子数の 0.568％の粒子が命中し， 2.48 × 1013[W]のエネルギーが加えられることがわかる． 将来的にはこれらの結果を用いて，支持構造物が粒子により破断されるかどうかを判定す る必要がある．
25
第
5 章 ローレンツ力により生じる支持構造物へ加わる応力解析
5.1 ヘルムホルツコイルに加わるローレンツ力
装置内部には磁場が存在するため，装置内部に存在し電流が流れているヘルムホルツ コイルにはローレンツ力が働く．ローレンツ力は次の式で示される． 𝐅 = 𝐼d𝐬 × 𝐁 (5.1) ここで，B はヘルムホルツコイルにかかる磁場であるが，この磁場が自身を含まない全て のコイルが生成することを考慮する必要がある．そこで，片側のヘルムホルツコイルを 取り除いた状態で磁場計算を行い，取り除いたヘルムホルツコイルの座標に生成された 磁場の強さを確認することでヘルムホルツコイルにかかる磁場を求める．このとき，磁場 は次のようになることが確認できた． ( 𝐵r 𝐵θ 𝐵z ) = ( −0.0964 [T] 0 0.291 [T] ) (5.2) (3.1)式について，d𝐬 × 𝐁を先んじて解くと，次のようなテンソル量で表される． d𝐬 × 𝐁 = 𝑅d𝜃𝐞θ× (𝐵r𝐞r+ 𝐵z𝐞z) = ( 𝑅𝐵zd𝜃 0 −𝑅𝐵rd𝜃 ) (5.3) r 方向の力について，r 方向の力を一周積分すると打ち消し合い，理論上 0 となる．これ を確かめる．ヘルムホルツコイルは半径𝑅の円であるため，線素ベクトル𝑅d𝜃𝐞θをデカルト 座標系に変換すると， 𝑟 = √𝑥2+ 𝑦2 𝐞θ= −sin𝜃𝐞x+ cos𝜃𝐞y (5.4) (5.5) これより，x 軸方向並びに y 軸方向の線素ベクトルは (d𝑥d𝑦) = (−𝑅sin𝜃d𝜃 𝑅cos𝜃d𝜃) (5.6) となる．これを𝜃について一周分積分すれば，三角関数の性質から 0 となることがわかり， 次の式よりヘルムホルツコイルに加わるr 方向の力は理論上 0 となることが確認できる． (𝐹𝐹x y) = ( −𝐼h𝐵z∮ 𝑅sin𝜃d𝜃 𝐼h𝐵z∮ 𝑅cos𝜃d𝜃 ) = 𝟎 𝐹r = √𝐹x2+ 𝐹y2= 0 (5.7) (5.8) このため，以後はz 軸方向の力のみについて考えればよい．ヘルムホルツコイルのある 微小な一点に対してz 軸方向にかかる力は，(5.1)，(5.2)，(5.3)式より d𝐹z= −𝐼h𝑅𝐵rd𝜃 (5.9) 一周積分して26 𝐹z= −2𝜋𝐼h𝑅𝐵r (5.10) この式に求められた磁場の値と装置パラメー夕を代入することで，ヘルムホルツコイル 間に加わるローレンツ力の値を求めることができる．符号から，ヘルムホルツコイルに 加わる力は引力であることがわかり，その大きさは(5.2)，(5.10)式より 𝐹z= 3.18 [MN] (5.11) となることがわかった．
5.2 von Mises 応力
以上の手順で求められたヘルムホルツコイル間に加わるローレンツ力からヘルムホルツ コイルの支持構造物に加わる応力を求めることで，破断に耐えるためのパラメータを導出 する必要がある．ここで，計算する応力としてvon Mises 応力を考慮する．von Mises 応力は本来 6 軸ある 応力の値を一つのスカラー値として代表して示された値であり，von Mises 応力𝜎Misesは 主応力𝜎1, 𝜎2, 𝜎3を用いて以下のように計算される． 𝜎Mises= √ (𝜎1− 𝜎2)2+ (𝜎2− 𝜎3)2+ (𝜎3− 𝜎1)2 2 (5.12) このvon Mises 応力から材料に必要な耐性を求める際には，材料が塑性変形しないような 閾値について考える必要がある．ところで炭素鋼以外の鋼材が塑性変形する場合，材料の 変形度合いは0.2[％]以上のひずみが生じていることになるため，逆説的にそれ以下の ひずみに抑えた場合は弾性変形であると言える．このひずみが 0.2[%]時点での応力の値を 材料の0.2％耐力または単に耐力と言い，この値が von Mises 応力よりも下であれば 塑性変形が起こらない閾値となっている。以下のFig.5-1 に炭素鋼以外の一般的な応力 ひずみ線図と0.2％耐力のモデルを示す。 Fig. 5-1. 応力ひずみ線図と降伏点[8]．
27
5.3 CAE 解析
von Mises 応力を求めるためには，支持構造物に対して(5.11)式で示されたローレンツ力の 値を用いて応力計算を行う必要がある．そこでAltair 社製の CAE ソフトウェアである HyperWorks を用い，計算を行った．このソフトウェアで計算を行う場合、以下の Fig.5-1 に 示した手順を取る必要がある． Fig. 5-1. 標準的な解析の工程[9]． 以上の図に示した工程から，CAD によるモデリングを行う必要がある．今回は HyperMesh によるモデリングを2D で行い，以下に計算に使用したモデルを示す． Fig. 5-2. 支持構造物の CAE 計算モデル．28 ここでCAE の計算に際し，有限要素法を用いるためメッシュの切り分けが必要となるが， 本研究では計算モデルが単純であるため正方形要素のみを考慮する．また，入力した材料 に関する情報は純タングステンに準拠し，以下のTable5-1 に示したものを使用する． Table 5-1. CAE 解析に用いる材料特性[10]． パラメータ ヤング率 [N/mm2] 345000 ポワソン比 0.284 この条件で解析を行った場合の支持構造物へ加わる応力を，HyperWorks に同梱されて いる構造最適設計ソフトウェアOptiStruct を用いて計算する．例として，軌道計算に用いた 支持構造物の断面積0.15×0.15[m2]に加わる von Mises 応力のコンタープロットを Fig.5-3 に
示す．
Fig. 5-3. 断面積 0.15×0.15[m2]の支持構造物に加わる von Mises 応力．
この結果から，固定端に接続されている部位で最大の値である41.2[kN/mm2]が現れる ことが分かり，考慮した素材である純タングステンにはこれ以上の耐力が必要となる． しかし純タングステンの耐力は80-150[kg/mm2]であり[11]，これを[kN/mm2]に変換すると 0.785-1.47 [kN/mm2]となるため，要求された耐力を満たしてはいない．そこで，この耐力を 満たすことができるような条件を探す．そこで断面積が0.15×0.15[m2]の支持構造物に対し， ヘルムホルツコイルの半径を変化させた場合に支持構造物へ加わる最大の応力とヘルムホ ルツコイル半径との相関をFig.5-4 に示す．
29
Fig. 5-4. ヘルムホルツコイル半径と支持構造物へ加わる応力の相関．
この結果と準タングステンの耐力を比較することで、ヘルムホルツコイルがおおよそ 0.275[m]以下の半径を持っていればローレンツ力による破断の影響を抑えることができる ことがわかる．この時に陽子の閉じ込めが成立するかを確かめるために(2.7)式で示される Accessible region を求めると，以下の Fig.5-5 に示す分布を得ることができる．
Fig. 5-5 ヘルムホルツコイル半径が 0.275[m]の時の陽子の Accessible region．
青の部分がAccessible region の範囲を示しており，この範囲を陽子が自由に運動できる． これより，この条件下では装置内の下部にAccessible region が形成されており，範囲は狭い ながらも陽子の閉じ込めができていることが分かる．
30
第
6 章 まとめと今後の課題
6.1 まとめ
今回考慮した装置内で発生した陽子は，支持構造物一本に対し最大で0.568 [%]が命中し， 2.48×1013[W]のエネルギーを与えることが分かる． またヘルムホルツコイルに加わるローレンツ力についてヘルムホルツコイル半径が 1.4[m]としたときには 3.18[MN]が加わり，41.2[kN/mm2]の応力が支持構造物に加わるが，純 タングステンではこれを耐えることができない．そこでこれを耐えるためにヘルムホルツ コイルの半径を変化させると，0.275[m]以下であれば破断せず，かつ閉じ込め領域を形成し うることがわかった．6.2 今後の課題
本研究にて行ったシミュレーションでは，プラズマの影響をGrad-Shafranov 方程式で計算 しているが，これが正しく論じられる手法であるかどうかは議論の余地がある．そのため 将来的にはKinetic な計算を行い現在の計算と結果を比較することが望ましい． また，支持構造物に与えられたエネルギーと命中分布から，核融合生成粒子による破断 判定を行う必要がある．そのためには支持構造物へと命中した粒子の角度を求め，原子の はじき出しの影響などを計算する必要がある． 他にも，ローレンツ力による支持構造物の問題について，引き続き計算を繰り返し行い 最適なパラメータを求めることやより強い材料で比較するなどの課題があり，これらに関 してアプローチを行うべきである．31
謝辞
本研究を進めるにあたり，研究の理論や方向性など多くのご指導いただいた髙橋俊樹准 教授に深く感謝致します． また，研究テーマの違いはありましたが，私のシミュレーションに対する不慣れや理論 に対する勉強不足に対し的確な指摘をいただいた本研究室の先輩方や，研究室での生活中 で多くの気付きを与えてくれた後輩一同には大変お世話になりました．この場を借りて御 礼申し上げます． 最後に，6 年間の大学生活を陰ながら支えてくださった群馬大学関係者の皆様並びに家族 一同に対し感謝の意を示し，謝辞とさせていただきます．32
参考文献
[1] 国際熱核融合実験炉 ITER ウェブサイト; http://www.naka.jaea.go.jp/ITER/iter/index.html. [2] 研究最前線/平成 27 年度核融合科学研究所プロジェクト成果報告, NIFS NEWS No. 229,
2(2016).
[3] Francis F. Chen : INTRODUCTION TO PLASMA PHYSICS AND CONTROLLED FUSION, p.10, Plenum Press, New York(1984).
[4] H. Momota et al., J. Fusion Energy 27, 77(2008).
[5] 柳雅裕, 非断熱性トラップによるプラズマ閉じ込めに関する研究, 修士論文(2017). [6] Thomas A. Mehlhorn, 2013, “NRL Plasma Formulary”, Naval Research Laboratory. [7] 小林信之，等方的ビーム/計算ノート， http://be.nucl.ap.titech.ac.jp/~koba/cgi-bin/moin.cgi/等方的ビーム [8] 明治大学，8.応力ひずみ線図｜材料力学， http://www.messe.meiji.ac.jp/~experiments/f09.pdf(閲覧日：2021.1.21)． [9] Altair Engineering，有限要素シミュレーションの実践(2019). [10] 岳石電気株式会社 ，タングステン基本物性， http://www.takeishi.co.jp/technology/leading01.html(閲覧日：2021.1.24)． [11] 日本ウエルディング・ロッド株式会社 ，実用金属および合金の物理的性質， http://www.nihonwel.co.jp/pdf_data/Capter17/alloy%20property.pdf (閲覧日：2021.1.24)．
![Fig. 5-3. 断面積 0.15×0.15[m 2 ]の支持構造物に加わる von Mises 応力．](https://thumb-ap.123doks.com/thumbv2/123deta/6324759.1097792/31.892.186.705.286.369/Fig53断面積15×15m2の支持構造物に加わるvonMises応力.webp)
