ヘリウム原子ビームを用いた逆磁場ピンチプラズマの磁場方向分布計測

目次 第１章 序論 -1- １-１ 核融合とプラズマ -1- １-２ 逆磁場ピンチプラズマ -4- １-３ 磁場方向分布計測方法 -6- １-３-１ ピックアップコイルやループを用いた磁場計測 -6- １-３-２ シュタルク効果分光計測 -6- １-３-３ 原子ビームを用いた磁場計測 -7- １-４ 研究目的 -9- 第２章 RFP プラズマの平衡 -10- ２-１ MHD（電磁流体力学）方程式系 -10- ２-２ Grad-Shafranov 方程式の導出 -10- ２-３ 数値計算法 -18- ２-３-１ 有限差分法 -18- ２-３-２ SOR 法 -20- ２-４ 圧力分布及びポロイダル電流束関数 -23- ２-５ Grad-Shafranov 方程式の規格化 -24- ２-６ Grad-Shafranov 方程式の差分化 -26- ２-７ 平衡計算モデル -29- ２-８ 装置参考装置 -32- ２-９ RFP プラズマにおける平衡解計算結果 -33- 第３章 電子衝突シミュレーション -37- ３-１ 中性粒子の電離 -37- ３-２ モンテカルロ法 -38- ３-３ ビームの分散 -39- ３-４ 電子衝突シミュレーション結果 -39- 第４章 粒子軌道シミュレーション -42- ４-１ プラズマ中のイオンビームの運動 -42-

４-２ クーロン衝突 -43- ４-２-１ Slowing-Down Collision -43- ４-２-２ ピッチ角散乱 -46- ４-３ 粒子軌道集計法 -48- ４-４ 粒子軌道シミュレーション結果 -50- 第５章 結論 -53- ５-１ まとめ -53- ５-２ 今後の課題 -53- 参考文献 -54- 謝辞 -56-

1 第1 章 序論 1-1 核融合とプラズマ 現在，我々の生活を担っている発電方式として原子力発電，火力発電，水力 発電，新エネルギー発電（太陽光発電、風力発電など）などが挙げられる。し かし現在，可動している発電方式は問題を抱えている。たとえば，水力発電， 波力発電などの新エネルギー発電方式に関しては，エネルギー使用の増加に伴 い上に挙げた発電方式のみでは充分なエネルギーを供給できていないのが現状 である。原子力発電，火力発電に関しては，燃料となる資源の枯渇が挙げられ る。石油で考えると 1 世紀半で約 90 倍も消費量が増加している。現在公表され ている世界のエネルギー資源確認埋蔵量は、石油が 1 兆 2,378 億バーレルで約 42 年分（BP 統計 2008），原子力発電のウランが 547 万トンで約 100 年分 （OECD/NEA-IAEA URANIUM 2007）であり，今尚，世界のエネルギー消費量 は着実に増えつづけている。また，このエネルギー資源枯渇の問題は電力需要 だけでなく日本を含む国々，特に発展途上国の人口増加と産業の発展によるエ ネルギー消費量の大幅な増加も後押ししているのが現状である。このエネルギ ー消費の約 90％が化石燃料を使用しており，それに伴い二酸化炭素の大量排出 による環境問題も深刻化している。二酸化炭素は，地球温暖化の原因とされる 温室効果ガスであり，大量に放出されると地球全体の平均気温が上昇し，海水 の膨張や氷河などの融解により，気候メカニズムに変化が起こり，異常気象が 頻発化する恐れがある。それにより自然生態系や生活環境，農業などへの影響 が懸念されている。 そこで，化石燃料に変わる大規模エネルギー源として注目されているのが， 核融合反応により発生させたエネルギーの利用による発電である。 核融合によるエネルギー利用の利点として，まず無尽蔵なエネルギー資源と いうことが挙げられる。最終的に核融合反応に使用する資源は主に水素の同位 核種の重水素であり，重水素は海水より容易に入手可能でき，ほぼ無限のエネ ルギー資源であると言われている。さらに炉心は超高温だが，蓄えられている エネルギーはわずかであるので暴走の危険が無いこと，化石燃料で問題になっ ている二酸化炭素を発生しないので地球の温室効果は無く，多少の中性子が放 出されるだけで放射性廃棄物は少ない等の利点がある。そのため、実現すれば エネルギー問題を解決することができる。 核融合を起こすには比較的に軽い原子を高温･高圧において物質の第四状態 であるプラズマ状態にしなければならない。プラズマ状態となった原子は正の 電荷を持つイオンと負の電荷を持つ電子の 2 つに分かれる。このとき，2 つの正

2 の電荷を持った原子核同士を核力（原子核の間に働く引力）の働く距離まで十 分に近づけクーロン力（静電的な斥力）に打ち勝ち，1 つに融合し別の原子核を つくり出す。この際，衝突前の総質量よりも衝突後の総質量の方が軽くなる。 これを質量欠損といい，この質量欠損分は，アインシュタインの導き出した式 2 mc E  から生成粒子の運動エネルギーとなる。このエネルギーを電気エネル ギーとして利用するのが核融合エネルギーである。利用可能な核融合反応とし て考えられているものは以下の通りである。 ) MeV 03 . 3 ( p ) MeV 01 . 1 ( T D D   ) MeV 45 . 2 ( n ) MeV 82 . 0 ( He D D  3  ) MeV 06 . 14 ( n ) MeV 52 . 3 ( He T D  4  ) MeV 67 . 14 ( p ) MeV 67 . 3 ( He He D 3  4  4.8MeV He T n Li6   4 n He T 2.5MeV) n( Li7    4 ここで、Dは重水素、Tは三重水素、pは陽子、nは中性子を表す。 Fig. 1-1 核融合反応 以上のことより，核融合反応を定常的に起こすには，プラズマを高温かつ高 圧で，反応が起こるまで閉じ込める必要がある。現在では核融合を実現するた めに 2 種類の閉じ込め方式が存在する。

1 つは慣性閉じ込め方式（ICF：Inertial Confinement Fusion）である。慣性閉じ 込め方式は主に高出力のレーザーを用いて核融合を起こすレーザー核融合が主

3 流となっている。原理は燃料ペレットという球状の燃料を考える。球の殻部を 重水素や三重水素の固体でつくり，球の内部には重水素や三重水素の気体で満 たしておく。この燃料ペレットに外部より均一に高出力のレーザーを照射する と，殻部の燃料は高圧のプラズマとなり中心に向けて加速され，燃料は高温高 密度状態となり，核融合が起こる。慣性閉じ込め方式の利点としては，低エネ ルギーでの発電が可能であり，磁場を発生するための巨大なコイルなどを必要 としないので，施設がコンパクトとなり，低コストであると言うことがいえる。 Fig. 1-2 慣性閉じ込め

2 つ目は磁場閉じ込め方式（MCF：Magnetic Confinement Fusion）である。核 融合を行う際，プラズマ状態にし，核融合を行う。プラズマ状態となった物質 は電荷を帯びているため電磁気学的な作用を受ける。つまり，磁力線の周りを 回転しながら磁力線方向に沿って自由に運動することができる。核融合を実現 するためには，高温プラズマを一定時間閉じ込めておかなければならないが， プラズマが放電容器に直接触れると，プラズマの温度は急激に下がってしまう。 そこで，プラズマが放電容器に直接触れないようにし，かつ高圧力にプラズマ を閉じ込めるために，磁場によりプラズマを閉じ込める方式を磁場閉じ込め方 式という。その代表的なものに「ミラー型」「ヘリカル型」「トカマク型」の 3 つのタイプがある。本研究は磁場閉じ込め方式の中でも「トカマク型」に属す る逆磁場ピンチと呼ばれるプラズマ核融合に関する研究である。

4

Fig. 1-3 磁場閉じ込め 1-2 逆磁場ピンチプラズマ

前節で述べた磁場閉じ込め方式の一つとして，逆磁場ピンチプラズマ（RFP : Field Reversed Pinch）[1] がある。逆磁場ピンチでは，ドーナツ状の放電容器（ト ーラス）の中心を貫く磁束を急激に変化させることによりトーラス内のプラズ マにトロイダル方向の電流が流れ，ポロイダル磁場（Bp）が作られる。また外 部のポロイダル方向に設置されたコイルに流れる電流によって作られるトロイ ダル磁場（Bt）が作られ、ポロイダル磁場とトロイダル磁場のらせん状の合成 磁場により，プラズマを閉じ込めている。らせん状の磁場にすることで磁力線 の網目が細かくなり，プラズマを閉じ込めやすくなる。よって逆磁場ピンチで はトカマクに比べ小さな外部トロイダル磁場で，プラズマを閉じ込めることが でき，磁場の利用効率が高くなり，弱い磁場で高温，高密度のプラズマを閉じ 込めることが出来るという利点を持っている。 さらに，強力なプラズマ電流の抵抗加熱だけで核融合反応の点火を実現でき るので，高価で複雑な補助加熱装置が不要となる。また，閉じ込め特性は直接 アスペクト比（トーラスの大半径と小半径の比）に依存しないため，炉の設計 の自由度が増し，エネルギーを取り出しやすくなる。 同じトカマク型に属するトカマクプラズマと大きく異なるところはプラズマ中 のトロイダル磁場分布にある。Figure 1-4 に示す様に，RFP では，トロイダル 磁場の向きがプラズマ中心部と境界部分で逆転しており，プラズマ内部ではト ロイダル磁場とポロイダル磁場が同程度の大きさを持っている。これがRFP の 由来である。また，このような磁場配位を取ることにより，弱いトロイダル磁 場で大きなプラズマ電流を流すにもかかわらず，プラズマの電磁流体的な安全 性を確保することが可能となっている。これらの点から、逆磁場ピンチは、磁

5 場閉じ込め方式として有効なものであるといえる。 しかし，RFPでは，解決が必要な課題が多く残っている。プラズマの安定性， 閉じ込め性能の大幅な改善や逆磁場ピンチ配位を定常的に維持するための技術 開発等である。特にプラズマの閉じ込め性能においてプラズマ中の磁場配位を 維持するために必要なダイナモ効果により磁場揺動が起こる。このため磁力線 がストキャスティックになり閉じ込めが悪化することが観測されている［2］。 磁力線がストキャスティックになると、そこでの熱や粒子の拡散が異常に大き くなるという現象が知られている。このことからプラズマ装置内での磁力線の 方向分布を解明することがRFPの実用化を進めるにあたり重要な要素の1つだ といえる。 Fig. 1-4 磁場の径方向分布

6 1-3 磁場方向分布計測方法 プラズマ実験において強力な磁場とプラズマの両者が存在する装置内の磁場の 空間分布を測定する必要のある分野において用いられる磁場計測方法をいくつ か紹介する。 1-3-1 ピックアップコイルやループを用いた磁場計測 プラズマ中で磁場を計測する方法として使われているものがいくつかある。 その中でも最も多く使用されている方法が，磁気ピックアップコイルプローブ を使用した計測方法である[3]。この方法はプラズマ中に挿入し，誘導磁場の変 動を検出し，低温プラズマの測定に広く使用される。そのコイル（またはルー プ）には，その場所でのコイルを貫く磁束の変化（コイルへの鎖交磁束の変化） に比例した起電力が生じる。この起電力を計測・解析（積分）することにより， 磁界の時間変化を求めることができる。 しかしこの計測法には欠点がある。RFP のような高温プラズマで使用すると， プローブ自体とプローブの材料が熔けてプラズマへ混入してしまい，純粋なプ ラズマを測定できなくなってしまう。また，プラズマの温度が高くなるとプロ ーブ自体がプラズマ粒子による損傷を受け破壊され計測ができなくなることも ある。したがって，プローブは通常の磁気閉じ込め実験において高温，高密度 プラズマでの計測には適していないといえる。 1-3-2 シュタルク効果分光計測 他に，運動シュタルク効果（MSE）に基づく計測法が挙げられる[4]。シュタ ルク効果とは外部電場が発生した場合に原子のエネルギー準位にずれが起こり， それらが出す光のスペクトルが分裂する現象である。電場のシュタルク効果に よる波長の変化を計測する。ここで，磁場をB，原子ビームの速度をｖとすると， 運動論的シュタルク効果による電場EはEvBとなり，この電場の強度とシュ タルク効果による波長変化には一定の関係があることから，波長変化から電場 の強度が求まり、ビーム速度は入射条件で与えられているので電場強度から磁 場を求めることができる。 シュタルク効果は高電子密度が高ければ高いほど効果がある。そのため，高 温プラズマにおけるトカマクやヘリカル装置のような強い磁場を持つ高温プラ ズマの閉じ込めのシステムで広く使用されている。また，この計測法では高エ ネルギーの大容量原子ビームを必要とするため，特殊な大型装置でないと利用

7 が難しい。さらに，装置が大掛かりなものとなるため計測装置が非常に高価と なることが問題である。 1-3-3 原子ビームを用いた磁場計測 この節で紹介する計測方法は，日本大学の平野氏より提案された手法である [2]。これは，原子ビームを用いた磁場計測であり，高温プラズマ中での局所範 囲における磁場の方向分布を求めることができる。新しい計測方法では高温プ ラズマにも適応しており，また磁場の強弱に影響されることはない。この計測 法では計測用原子ビームとしてヘリウム(He)を採用する。He を採用した理由と して， ①長く 1 価のイオンで滞在している。 ②水素と比べて 1 価及び 2 価への電離エネルギーが高い。 ③低エネルギーで高速のビーム及び浸透距離が得られる。 ④低 Z で不純物としての影響が少ない というような利点が挙げられる。 では計測方法の原理について説明する。システムの概略図をFig. 1-5 に示す。 He 原子ビームは装置の赤道面に，磁場に垂直に入射する。入射された He ビー ム原子はプラズマ中の電子との衝突により，He+にイオン化される。その後， He+イオンは，磁場中で旋回運動をおこし，旋回運動中に励起し基底状態に戻 るとき468.6 nm の発光スペクトルを放出する。また，He+イオンが完全にイオ ン化するか，様々なドリフト運動によって検出領域から離れるまで発光スペク トルを放射する。この発光スペクトルをドップラー計測器で波長幅を計測する ことで磁場方向分布を求めることができる。 He+は一定の速さv で旋回運動をし，その回転軸は磁力線のなす角つまり磁_b 場の方向に一致する。磁場と回転軸の関係は Fig. 1-6 のように表される。装置 内ではポロイダル磁場のみが存在する場所(A)，トロイダル磁場のみが存在する 場所(C)，及びその両方が存在する場所(B)が存在する。これらにより磁場分布に よ っ て 磁 力 線 の 作 る 角 度 に 変 化 が 起 こ る の が わ か る 。 磁 場 の 傾 き 角 は p t t B B   　/ 　B  より求まる。また，ドップラー計測器によりドップラー幅を求 めることでき，それは計測場所による磁場分布(A)，(B)，(C)により異なり，Fig.1-7 のような結果になる。波長の変化幅は磁場の傾き角度が90°の場合に最も大き く、角度が小さくなるに従って小さくなり、角度が0°，即ち観測点から見て回 転面が垂直となる場合には 0 となる。これよりドップラー幅は旋回運動の回転

8 軸の傾きと連動しているのがわかる。また，原子ビームのエネルギーv はビー_b ムの入射時に決めているので、Fig. 1-7 で示した波長の最大変化幅は計算で容易 に求めることができる。ここで，₀ 468.6nm，c は光速を表す。この波長の最 大変化幅と計測した波長の変化幅との比の大きさから，発光点における磁力線 の傾斜角度を求めることができる。 そこで，He+イオンからの線スペクトルを計測し，場所によってスペクトル 形状に変化があることから，旋回運動（磁力線の方向）の軸の傾きが決まり， 最終的に磁場の方向分布を計測することができる。 Fig. 1-5 計測原理

9

Fig. 1-6 磁場と回転軸のなす角の関係

10 1-4 研究目的 前節で本研究において使用する計測方法を紹介した。この計測方法は理論上 可能であるということだけで，実際のところ正確に計測できるのかはまだ実証 されていない。 そこで本研究ではヘリウム原子ビームを用いた新しい磁場方向計測の実現可 能性をシミュレーションで検討する。実際にヘリウム原子ビームの電離位置シ ミュレーション，粒子軌道シミュレーションを行い，旋回運動の回転軸のなす 角を求めることで，磁場方向計測の可能性を評価することが本研究の目的であ る。

11 第2 章 RFP プラズマの平衡 本章では，プラズマの平衡状態とはプラズマの巨視量が時間と共に変化しな い状態のことである。ここでは平衡状態にある RFP プラズマを扱う際に用いら れる方程式について説明する。 平衡状態にあるRFP プラズマが圧力，磁場，電流密度などについてどのよう な値を持つのかを2 次元円筒座標系rz空間で数値計算する。 2-1 MHD（電磁流体力学）方程式系 プラズマの性質を調べるにはプラズマ全体の動きを正確に捉えることが必要 であり，その巨視的な動きを調べるために有効なモデルが MHD（電磁流体力学） モデルである。このモデルではプラズマを全体で一つの流体と見なしており， それを取り扱う方程式を MHD 方程式という。本研究では高温プラズマを仮定し ているため，以下のような電気伝導度を無視した理想 MHD 方程式を用いている。

 

0      u   t (2.2.1) u



u



u  jB       _{ }   p t  (2.2.2)





0       _{ }   _ p t u (2.2.3) 0   u B E (2.2.4) t       E B (2.2.5) j B



₀   (2.2.6) 0    B (2.2.7) ここで，E，B，j はそれぞれ電場強度，磁場強度，電流密度を表している。また，



，u，p，



，



0はプラズマの質量密度，速度，圧力，比熱比，真空の透磁 率を表している。 2-2 Grad-Shafranov 方程式の導出 プラズマの平衡は一般的にGrad-Shafranov 方程式を解くことによって求め

12 られる。プラズマは平衡状態にあるとし，時間に関する項は0 となるので t0 となる。次に，RFP は軸対称性を持った円筒座標系でモデル化される。 また，この平衡状態は静的であるため流速u0（CT 全体は移動していない） とする。さらに，電荷準中性条件が成り立つとする。 プラズマ中のイオン及び電子の運動方程式（流体方程式）は，力学的及び熱 的平衡により圧力テンサを0（等方的）として，



i



i i i



i



i ie i i i qn p t n m u u u  Eu B  R     _{ }   (2.2.1)



e



e e e



e



e ei e e e q n p t n m u u u  Eu B  R     _{ }   (2.2.2) 先程の熱的平衡より，プラズマ中のイオンはほぼ静止している（u_i 0）と仮定 し，プラズマ中の電子は質量が小さいので0 と近似すると， ie i i in p q E R  0 (2.2.3)



e



e ei e en p q Eu B  R  0 (2.3.4) ここで，プラズマ電流は定義式(1.5)より， e e e i i in qn q u u j  (2.3.5) e e en q u j (2.3.6) 流体方程式に戻り，式(2.2.3)と(2.2.4)を加算すると，



qini qene



Eqene



ueB

 

 pi pe



RieRei  0 (2.2.7) ここで，衝突項（摩擦項）は作用反作用による完全弾性衝突より，R_ieR_ei 0 となり，プラズマ全体の圧力を



p_i p_e



pとおくと，



q_in_i q_en_e



Eq_en_e



u_eB



p (2.2.8)

13 最初に電荷準中性条件が成り立つと仮定した事により電場の項は消え，



u_eB

 

 _{e e}u_e



B e en q n q であるので，式(2.2.6)より式(2.2.8)は次のように書き換 えられる。 p   B j (2.2.9) これがプラズマ平衡の式である。 次に，式(2.2.9)の両辺に

B

を内積すると，







 0   p B j B B (2.2.10) 同様に，(2.2.9)式の両辺にjを内積すると，







 0   p j j B j (2.2.11) 式(2.2.10)と(2.2.11)より，磁場と電流はpに対して垂直かつ等圧面上にあるこ とが分かる。 よって，この磁場による磁気面（磁力線が沿う等圧面）を関数で定義し，



 rrBzdr 0 ` `  (2.2.12) とする。 円筒断面積dS 2rdrを利用して，円筒座標系に書き換えると式(2.2.12)は，





    S SB ndS  B dS   2 1 2 1 　 (2.2.13) となる。式(2.2.23)はポロイダル磁束関数と言い，磁気面上の磁束の強さを表す。 ここで，ベクトルポテンシャルAより磁場Bは， A B (2.2.14) なので，式(2.2.13)は，







   S A dS   2 1 (2.2.15) となる。 ストークスの定理より，式(2.2.15)は，



  C dl A   2 1 (2.2.16) となる。この線積分の積分路は，半径rの円周に沿った閉曲線である。

14 円の中心は円筒座標系の中心（r 0）に相当し，積分路は円の正の方向であ るとする。 よって， r A rd A d C     2 2 0   





A l (2.2.17) であり， r A_



 (2.2.18) となる。 また，式(2.2.14)を円筒座標系で計算すると，軸対称性より 0    なので，

 

z z r r z z r r rA r r r A z A z A r A r A r A z A z A e e e e e e A B                                       _                                    1 　　　 　　　 　 　 (2.2.19) であり，式(2.2.18)より， z r z A B_r           1  _(2.2.20)

 

r r rA r r Bz        1 _ 1  (2.2.21) となる。 次に，アンペールの式を利用してポロイダル電流関数との関係を求める。 アンペールの式(2.2.23)はマクスウェル方程式の内，アンペール-マクスウェル式 (2.3.22)の電場（電束密度）の時間変化（変位電流）を無視し，磁場を磁束密度 に書き換えたものである。 t       H D j (2.2.22) B j 0



(2.2.23) ここで，アンペールの式(2.2.23)を半径

r

の円で面積分すると， S を円の面積と して，

15









    S S0j ndS B ndS (2.2.24)









    S S0j dS B dS (2.2.25) となり，右辺はストークスの定理より，





_



    C S B dS B dl (2.2.26) となる．この線積分の積分路は，半径

r

の円周に沿った閉曲線である。 円の中心は円筒座標系の中心（r 0）に相当し，積分路は円の正の方向であ るとする。 すると，





  2 _  0 B rd d C l B (2.2.27) である． ここで，式(2.2.25)の左辺を使用して式(2.2.28)を定義する。



  S d F 0j S 2 1 _  (2.2.28) 式(2.2.28)はポロイダル電流関数，またはトロイダル磁場関数と言い， 平曲面中にどれだけ電流が流れるのかを表す。 Fig. 2-1 トロイダル方向とポロイダル方向 式(2.2.26)及び(2.2.27)より，





d B rd B r S     2 2 0     





B S (2.2.29) となり，式(2.2.29)を式(2.2.25)に代入し，定義式(2.2.28)を利用すると， r B d S     2 2 1 2 ₀       _



j S (2.2.30)

16 となり，F B_rであるので r F B  _(2.2.31) ここで，式(2.2.23)中のBを円筒座標系で計算すると，





z r z r r z z z r r z r r B B r r B r B z B z B B r B B B z r r e e e e e e e e e B       _                                                             1 1 1 　　　 (2.2.32) となり，軸対称性より 0    なので，

 

z z r r z z r r rB r r r B z B z B r B r B r B z B z B e e e e e e B                                       _                                   1 　　　 (2.2.33) である。 よって，式(2.2.23)より電流の各方向成分を求めると，式(2.2.31)より，





z F r z B j r r _                  0 0 0 1 1     B (2.2.34)





               r B z B j r z 0 0 1     B (2.2.35)





_{ }

r F r rB r r jz z _                 0 0 0 1 1 1     B (2.2.36) となる。 ここで，式(2.2.20)及び(2.2.21)を利用して式(2.2.35)を書き換えると，

17                                                                   2 2 0 0 0 1 1 1 1 1 1 z r r r r r r r r z z r r B z B j r z         (2.2.37) である。 式(2.2.10)にいったん戻り，右辺を展開すると軸対称性より 0    なので，





z p B r p B B B B z p p r p p z r z z r r z r                          　　　 e e e e e e B _{  }  _(2.2.38) となり，式(2.2.20)及び(2.2.21)より， z p r r r p z r z p B r p B p _{r z}                          1 1 　　　 B (2.2.39) である。 ここで，圧力pが,の関数であると仮定すると式(3.39)は，

















_                                                       z r r z p r p z r r z r z p r r r p z r                  , 1 , 1 , 1 , 1 (2.2.40) となり，式(2.2.40)及び(2.2.10)が成り立つためには，

 

0 , _      p _(2.2.41) とならなければならない。よって，pはによらないので，p p

 



となる。

18 同様に，式(2.2.11)の右辺を展開すると軸対称性より 0    なので，





z p j r p j j j j z p p r p p z r z z r r z r                          　　　 e e e e e e j _{  }  _(2.2.42) となり，式(2.2.34)及び(2.2.36)より， z p r F r r p z F r z p j r p j p _{r z}                        0 0 1 1   j (2.2.43) である。 ここで，Fが,の関数であると仮定すると

















_                                                       z r p r z p F r F z r p r z p r z p r F r r p z F r                , 1 , 1 , 1 , 1 0 0 0 0 (2.2.44) となり，式(2.2.44)及び(2.2.11)が成り立つためには，

 

0 , _      F _(2.2.45) とならなければならない。よって，Fもによらないので，F F

 



となる。 ここで式(2.2.9)に戻り，

r

方向成分に注目すると，  B j B j r p z Z     (2.2.46) なので，式(2.2.46)に式 (2.2.21) ,(2.2.31), (2.2.35), (2.2.36)を代入すると， r F r F r r r z r r r r r r p                                            0 2 2 0 1 1 1 1      (2.2.47) となり， p p

 



及びF F

 



より，

19 r d dF r F         (2.2.48) r d dp r p         (2.2.49) なので，式(2.2.47)は， F r d dF r z r r r r r d dp r _                                                        ₀ 2₂ 2 1 _(2.2.50) となり，移項して を r   取り出すと， r z r r r r r d dF F r d dp r                                        _        ₀ 2₂ 2 1 _(2.2.51) である。 よって， 2 2 0 2 1 z r r r r d dF F d dp r                       (2.2.52) となり，円筒座標系におけるGrad-Shafranov 方程式が導出された。 2-3 数値計算法 ここでは本研究で用いた主な計算法について述べる。数値計算で微分量を求 めるのに離散化された値で近似する方法として差分法を用いた。また，RFP プ ラズマの平衡解を求めるのにSOR 法を用いる。 2-3-1 有限差分法 差分法は，時間や空間などの微分量を離散化された値を用いて近似する方法 である。例えば，x の一次元空間において x の関数 f

 

x があるとする。エラー! 参 照元が見つかりません。2-1 のような，ある任意の点x₀での微分量 0 x x x f    を，離 散点x₀Δx，x₀ Δxでの値 f



x₀Δx



， f



x₀Δx



を用いて近似すると，

20









Δx Δx x f Δx x f x f x x 2 0 0 0        (2.3.1) と表すことができる。但し，Δx は十分小さいとする。このように，微分量を離 散点で近似することを差分化するという。また，この差分式の真の値との誤差 は

 

_{Δx オーダーなので２次精度の差分式という。また求めたい点を中心に同じ}2 だけ離れた２点を使った差分を中央差分という。 Fig. 2-2 差分法 (2.3.1)の差分式は Taylor 展開から導くことができる。Taylor 展開を用いると



x Δx



f ₀ ， f



x₀Δx



は，それぞれ，



_



_

 

_  

_{ }

_

_{ }

 

_{ }

_

_{ }

 

_{ }

_ x f Δx x f Δx x Δxf x f Δx x f _{0 0} 1 2 2 3 3 ! 3 1 ! 2 1 _(2.3.2)





 

 

_{ }

_{ }

 

_{ }

_{ }

 

_{ }

      Δx f x Δxf x Δx f x Δx f x x f 0 0 1 2 2 3 3 ! 3 1 ! 2 1 _(2.3.3) と表すことができる。 f n

 

x の添え字の n は x の n 階微分を表す。 (2.3.2)式と(2.3.3)式の差をとると，  

_{ }

_















_

_{ }

 

_{ }

_ x f Δx Δx Δx x f Δx x f x f 1 0 0 2 3 6 1 2 (2.3.4) となり，Δx が十分小さいという仮定から右辺第二項目は第一項目に比べれば十 分小さいので切り捨てると，まさに(2.3.1)式と等しくなる。 ２次精度の中央差分は，２点を使って２次精度の近似を行っているが，使う 点を増やすことによって，精度を上げることが可能である。ただ，単に精度の 上げても解の収束性や計算スキーム全体の安定性によって正しく計算できると 0

x

Δx

 

x

₀

f

Δx

Δx

x

₀



Δx

x

₀





x

Δx



f

₀



f



x

₀



Δx



x

21 は限らないこともあり，注意が必要となります。また使う点が増えることによ って計算領域の境界条件が厳しくなるため本研究では，２次精度または４次精 度を用いている。 ４次精度の中央差分も同様にしてTaylor 展開から導くことができる。  

_{ }





 









 





Δx Δx x f Δx x f Δx x f Δx x f x f 12 2 2 8 _{0 0 0 0} 1         (2.3.5) また二階微分についても本研究ではよく使うので，２次精度，４次精度につい てそれぞれ示す。  

_{ }





  



 

2 0 0 0 2 2 Δx Δx x f x f Δx x f x f      (2.5.6)  

_{ }

 





 









 





 

2 0 0 0 0 0 2 12 2 2 16 30 Δx Δx x f Δx x f Δx x f Δx x f x f x f           (2.3.7) 2-3-2 SOR 法 次のような n 元連立一次方程式を解くとする。                                      n n nn n n n n b b b x x x a a a a a a a a a          2 1 2 1 2 1 2 22 21 1 12 11 (2.3.8) もしガウス・ジョルダン法のような直接法を用いた場合，その計算回数は 3 n に 比例するため n が大きい程数値計算にかかる時間は大きくなってしまう。しかし 収束法を用いた場合，解の収束条件によっては直接法を用いた場合より格段に 早く計算することが可能である。SOR(Successive Over Relaxation)法は、その 収束法の一つである。

連立方程式

b

Ax (2.3.9) を満たすxを数値計算により求めるとする。ここで真の解をxとする。

22 そしてある計算を m 回繰り返し計算した解をxmとしたとき、 x x    m m lim (2.3.10) が成り立つとする。このような計算回数を増やすことで真の解に近づけていく 方法を収束法と言う。 このある計算方法というのは簡単な手順で作ることができる。行列Aを， T S A  (2.3.11) と分解し， b Tx Sx_{m 1}  _m (2.3.12) とすればよい。x_mがα に収束したとすれば(2.3.9)と(2.3.12)から，αと x が等し いことは自明である。つまりx_mが収束すれば真の解 x が求まることに等しいと いえる。 ここで誤差について考えると，真の解では(2.3.12)式は， b Tx Sx  (2.3.13) であるので差をとると，



x xm



T



x xm



S  _1   (2.3.14) ここでxx_m1とxx_mは真の解からの誤差を表す。xx_m ε_mとおくと， m m   S T  1 (2.3.14) 最初の計算であらわれた誤差が₁とすると， 1  m m       S T _(2.3.15) と表すことができる。 (2.3.11)のS を対角行列とすると， xm1を求めやすくなる。つまり，行列Aを 次のように分解する。

23                                       0 0 0 0 0 0 0 0 0 2 1 2 21 1 12 22 11 2 1 2 22 21 1 12 11                      n n n n nn nn n n n n a a a a a a a a a a a a a a a a a a (2.3.16) 右辺の第1 項が S ，第２項がTである．また S の逆行列S1は，                     1 1 22 1 11 1 0 0 0 0 0 0 nn a a a        S (2.3.17) である。つまりm1回計算したときの近似解は，xm S



Txmb



1 1 なので，  





       





 





       





 





       





 





       m





n nn m n m n m n n nn m n m n n m m m m m n n m m m m m n n m m m m x a x a x a x a b a x x a x a x a x a b a x x a x a x a x a b a x x a x a x a x a b a x                                  4 4 3 3 1 1 1 1 3 4 34 3 33 1 31 3 1 33 1 3 2 4 24 3 23 1 21 2 1 22 1 2 1 4 14 3 13 2 12 1 1 11 1 1 (2.3.18) と計算することができる。これをヤコビ法という。 そして計算効率をよくする為に，  1 2  m x を計算するときにx₁ m の部分に計算で求 め ら れ た  1 1  m x を 代 入 し て 計 算 す る と い う よ う に ， x_nm1を 計 算 す る 時 に    m n m x x₁  ₁の部分に計算で求められた    11 1 1    m n m x x  を代入する方法をガウス・ザ イデル法という。(2.3.18)式は，  





       





 





       





 





       





 





     1  1





4 4 1 3 3 1 1 1 1 1 3 4 34 1 3 33 1 1 31 3 1 33 1 3 2 4 24 3 23 1 1 21 2 1 22 1 2 1 4 14 3 13 2 12 1 1 11 1 1                                    m n nn m n m n m n n nn m n m n n m m m m m n n m m m m m n n m m m m x a x a x a x a b a x x a x a x a x a b a x x a x a x a x a b a x x a x a x a x a b a x      (2.3.19) と書き直すことができる。 ここで，ガウス・ザイデル法では 1 回の計算にx_m1x_mだけ解を修正する。 これを 1 度に



x_m1x_m



だけ修正するというのが SOR 法である。は加速係 数を表す。加速係数は1 だとガウス・ザイデル法に等しくなり，1 以下だと遅く

24 なる。1 以上だとより早く収束するが 2 以上だと発散してしまうため 1～1.9 の 間で選ばれるのが一般的である。SOR 法を用いた場合(2.3.19)は，  





       





   



   



 





       





   



   



 





       





   



   



 





       





   



   m



n m n m n m n m n nn m n m n m n n nn m n m m m m m n n m m m m m m m m m n n m m m m m m m m m n n m m m m x x x x x a x a x a x a b a x x x x x x a x a x a x a b a x x x x x x a x a x a x a b a x x x x x x a x a x a x a b a x                                                        1 1 1 1 4 4 1 3 3 1 1 1 1 1 3 1 3 3 1 3 3 4 34 1 3 33 1 1 31 3 1 33 1 3 2 1 2 2 1 2 2 4 24 3 23 1 1 21 2 1 22 1 2 1 1 1 1 1 1 1 4 14 3 13 2 12 1 1 11 1 1          (2.3.20) と計算できる。 2-4 圧力分布及びポロイダル電流束関数 (2.2.52)式を磁束関数ψについて解くためには，圧力p

 

ψ 及びポロイダル電流 束関数F

 

ψ の分布を知らなければならない。ここでは，次のように与えて計算 することにする。 圧力p

 

ψ は，

 



b_p



ap p p ψ  01ˆ (2.4.1) ただし，



0

 

/ lim 0



ˆ











   ,_lim 0.0 (2.4.2) また，ポロイダル電流束関数F

 

ψ は，

 

         2 . 0 ψ 2 . 0 ψ 2 2 1 0 1 0 　　　 　　　　　　     B B B A A F (2.4.3) と設定した。ここで，ˆは規格化した磁束関数で，₀，_limはそれぞれ磁気軸 及び境界における磁束関数であり，p₀は磁気軸におけるプラズマ圧力値を表す。これ により各乗数



a ，_p b_p



や定数A₀,A₁,B₀,B₁,B₂用いて圧力分布及びポロイダル電

25 流分布を制御することで、様々な平衡配位が得られる。 2-5 Grad-Shafranov 方程式の規格化 数値計算で物理量を含む式を扱う場合，その中に出てくる変数や定数に勝手 な単位を使い，数値の間に成り立つ関係を示しても，単位を変えると意味がな くなってしまうため，規格化（無次元化）する必要がある。また，規格化を行 うことにより，ある物理量と別の物理量との相対的な大きさなどが明確になる など利点がある。 ここでは Grad-Shafranov 方程式の規格化を行う。規格化変数は以下の通り である。 a r r  ， a z z  ， max ψ ψ ψ  ， max p p p  ， max F F F  (2.5.1)

ここで， a はトーラス小半径，ψ_max,P_max, F_maxは磁気軸での磁束関数及び圧力関 数を表す。また，





0 2 2 max max / / ~    a p  ， max max/ ~ F a q  (2.5.2) とおくと，規格化基準量は



_



  I I F  0 _tfc  max 2 ，



max aF ~maxq， 0 2 2 max max / ~  _         a p (2.5.3) とかける。I_tfc，I_はトロイダル磁場コイル電流，ポロイダル電流を表す。 これに従い，(2.2.52)式は次のように規格化される。      d d ~ 1 d d ~ 1 2 2 2 2 F F q p r z r r r r                _(2.5.4) 変数の上についている『－』は規格化された物理量を表し，その次元は無次 元である。 次に規格化した磁場B_r,B_,B_z及び電流密度 j_r,j_,j_zを求める。 先ほど同様に，以下のように規格化変数を設定する。

26 2 max/ a B B   ， ) /( ₀ 3 max a j j    (2.5.5) (2.5.5)式を(2.2.20)，(2.2.21)，(2.2.31)の磁場の式へ，また(2.2.34)～(2.2.36)の 電流密度の式へ適応すると， z r B_r     1  ， r F q B_  _~1 ， r r B_z    1  ， (2.5.6) z F r q j_r     _~1 1 ，     d d ~ 1 d d ~ 2 F r F q p r j   ， r F r q j_z    _~1 1 (2.5.7) 次に，RFP の平衡を考えるうえで重要なパラメータとしてピンチパラメータ ，逆転パラメータFがある。これらは 断面積ポロイダル磁場

 

r,z r z B a B a a z a R z a R d d 1 2 2 2 2 2

 

_      _  _ (2.5.8) 表面平均トロイダル磁場

 

     d 2 1 2 0 B r,z B 



(2.5.9) 表面平均トロイダル磁場

 

   d 2 1 2 0 p p B r,z B 



(2.5.10) を用いて  B Bp   ，   B B F  (2.5.11) で表すことができる。各規格化基準量を変更することにより，1.450.05， 05 . 0 1 . 0    F になるよう最適な値を取るパラメータを求める。このパラメータ の決定はこの章の参考対象装置で行う。 2-6 Grad-Shafranov 方程式の差分化 数値計算するために，Grad-Shafranov 方程式を差分化する。差分精度は 2 次 精度の中央差分を用いた。(2.5.4)の第一項目は次のように差分化できる。

27









 

2 , 1 , , , 1 2 1 , 2 1 , 2 1 2 1 , 2 1 , 2 1 , 2 1 , 2 1 ψ ψ 2 1 ψ ψ 2 1 ψ ψ 1 ψ ψ 1 ψ ψ 1 ψ 1 r Δ r Δ r r Δ r r r Δ r Δ r r Δ r r r Δ r r r r r r r j i j i j i j i i i j i j i i i j i j i j i j i                                                                                           　









 

                       2 2 ψ ψ 2 ψ ψ 2 ψ 2 , 1 , 1 , 1 , , 1 r Δ r r Δ r r Δ r Δ r r j i j i j i j i j i (2.6.1) 同様にして第二項目は，

 

2 1 , , 1 , 2 2_ψ ψ ₂ψ ψ r Δ z j i j i j i        (2.6.2) である。よってψ についてまとめると， _i,j









                                      _{ }     d d ~ 1 d d ~ ψ ψ 2 2 ψ ψ 2 ψ ψ 1 ψ 2 ₂ 1 , 1 , , 1 , 1 , 1 , 1 , F F q p r r Δ r r Δ r r Δ r r i j ij j i j i j i j i j i (2.6.3) となる。ただし，

28                               1 2 2 2 2 r r r r r  (2.6.4) (2.6.4)式を 2 章で説明した SOR 法を用いて数値計算することによって磁束 関数 を求めることができ，磁束関数 から続いて磁場B_r,B_,B_z及び電流密度 を求めることができる。

29 2-7 平衡計算モデル 従来のRFP プラズマにおけるシミュレーションでは円筒近似モデルを用いて 計算していたが，より実験に近い数値計算を行うために，円筒座標系を用いた トーラス座標で計算を行う。 Figure 2-3 のようにトーラス座標を用いれば装置を円筒近似する必要がない がトーラス座標を用いると，境界条件の設定が非常に困難になる。またr，，z の全成分の数値計算が必要になり，計算式が非常に複雑になってしまい計算時 間が大幅にかかってしまう。 そこで，本研究ではFig. 2-4 の装置を Fig. 2-5 のような円筒座標系に配置す る。このような円筒座標系を用いたトーラス座標モデルを用いることにより， 以下のような利点が挙げられる。 ① 実際のプラズマの駆動に近い分布を再現 ② 境界条件を Fig. 2-5 のような簡単な曲線境界に設定できる ③ 軸対称により r, z成分のみに依存した数値計算が可能 この座標系ではトロイダル磁場，ポロイダル磁場などのプラズマ実験におけ る物理量は Fig. 2-6 のように表すことができる。また，以前の円筒近似での解 析は r方向のみの 1 次元の分布を求めれば，対称性より 3 次元の分布を導出で きたが，トーラス座標では r, z方向の２次元の分布を求め， 方向の対称性よ り，３次元の分布を導出することができる。

30

Fig. 2-3 トーラス座標(b,,)の概念図

31

Fig. 2-5 円筒座標系を用いたトーラス座標

32 2-8 参考対象装置 円筒近似計算での数値解析の参考対象装置としてイタリアにある RFX_mod 装置を計算の際の各パラメータ，初期値等で使用した。主要装置パラメータを Table 1 に示す。また，平衡計算で使用した計算用パラメータを Table 2 に示す。 Table 1. 装置参考パラメータ

装置パラメータ

トーラス大半径 r 2.00〔m〕 トーラス小半径 a 0.46〔m〕 プラズマ電流 IP 2.0〔MA〕 電子温度 T _e 250.0〔eV〕 イオン温度 T i 250.0〔eV〕 電子密度 N _e 19〔 -3〕 m 10 0 . 1  Table 2. 平衡計算パラメータ

平衡計算パラメータ

ピンチパラメータ  1.47 逆転パラメータ F -0.09 圧力乗数



a ，p bp



（4 , 2） ベータ値  _0.002 安全係数 q _0.4 全ポロイダル電流  I 11 [MA] トロイダルコイル電流 tfc I -0.3 [MA]

33

2-9 RFP プラズマにおける平衡解計算結果

RFP プラズマの磁束関数



の計算結果をFig. 2-7 に示す。また磁束，磁場， 圧力，密度及び温度の赤道面における分布を，それぞれFig. 2-8，Fig. 2-9， Fig. 2-10，Fig. 2-11 に示した。A とはアスペクト比のことを指し，トーラス大 半径をトーラス小半径で割ったものである。この研究ではアスペクト比を 4.35 に設定している。磁束関数は磁気軸が外側にシフトしているのがわかる。磁場 に関しては，トロイダル磁場B_tが中心と装置壁で反転し、ポロイダル磁場B の_p 分布が第 1 章 2 節で述べた RFP の磁場分布に近い値となっていることが見て取 れる。電流密度も磁場分布に沿って同様の分布を取っている。また,圧力pが中 心付近で平らになっており，温度T，密度 n などの各種物理量も磁気軸で最大と なりフラットな値を取っていることがわかる。以上のことから求めた分布は実 験プラズマを再現できていると言える。 Fig. 2-7 プラズマの磁束関数

34

Fig. 2-8 赤道面上における磁束分布

35

Fig. 2-10 赤道面上トロイダル電流及びポロイダル電流

36

37 第3 章 電子衝突シミュレーション 3-1 中性粒子の電離 原子ビームとしてHe 原子をプラズマに入射する。装置内のプラズマは水素 プラズマを想定する。入射したビームがプラズマに進入し，その中性粒子はプ ラズマ中の粒子と衝突して電離し，イオンビームに変換される。その際の反応 過程はいくつかあり，主に ①荷電交換反応：He+ H+ _{→ He}+ _{+ H} ②イオン衝撃による電離：He + H+ _{→ He}+ _{+ H}+ _{+ e}- ③電子衝撃による電離：H + e- _{→ H}+ _{+ e}- _{+ e}- がある。本研究ではこれらの反応のうち，③の電子電離のみを考慮する。また， He は一価のイオンとしての反応を仮定する。 電離判定の起こりやすさを表すものとして断面積がある。これが大きければ 反応が起こりやすいことを表している。この断面積は実験で調べられており， Fig. 3-1 が He の電離断面積を示している。 Fig. 3-1 He の電離断面積

38 3-2 モンテカルロ法 モンテカルロ法は,フォン・ノイマンとユラムによって1945 年頃から用いら れはじめたものであり,当初は「決定論的な数学の問題を乱数を用いて解くこ と」が目的であった。しかし,最近では精度の高い乱数が得られることにより, それを用いたシミュレーションが盛んとなり,現在ではその手法をモンテカル ロ法と呼ぶ。モンテカルロ法は,確率論的モデルに基づいた確率過程を問題とす るのである。モンテカルロ法では,現象のモデルを構築し,その様子を調べよう というものである。例としては,形状の複雑な領域の面積の計算などがある。シ ミュレーションでは,乱数を発生しながら確率過程を計算機上で試行していく。 その最も基本的なアルゴリズムを以下に示す。 Fig. 3-2 モンテカルロ法のアルゴリズム[6] コンピュータ上で発生させた乱数は先にも述べたように,一様に理想的な乱 数ではなく,疑似乱数と呼ばれるものである。疑似乱数というのは,理想的な一 様乱数を模倣したものであるので,やはり完全に一様なものではない。そこで, 性質の良い乱数をいかに多数発生させることができるかどうかが,モンテカル ロ・シミュレーションにおいて良質の結果を得るためのポイントとなる。 あるイベントが起こる断面積を とすると,そのイベントが t の間に起こっ たとする判定方法は,0~1までの一様な擬似乱数を用いて

n



 

v

v



t



ξ

(3.2.1) と計算することができる。ここで n はプラズマ密度, v は中性粒子とプラズマの 相対速度である。ここで注意するのは t のとり方である。これを大きくすると (3.2.1)の左辺は 1 を超えてしまう。このとき t は長く取りすぎており,1 を超え た分だけ起こる回数を過小評価することになる。そのため, (3.2.1)の左辺は 1 よ

39 りも小さくなるように注意しなければならない。 3-3 ビームの分散 RFP へ NB が入射された際に，Fig. 4-2 に示すようなビームの分散が実験に おいて確認されている。本研究では，より現実的なシミュレーションを行なう ためにこの分散を考慮し計算している。分散による粒子数の振り分けは正規乱 数によって行っており，ビーム軸でもっとも粒子が多く，そこから離れるにつ れて粒子が少なくなるようにしている。本研究では原子ビーム入射は垂直方向 入射におけるビーム分散のみ考慮した。 Fig. 3-3 ビーム分散モデル 3-4 電子衝突シミュレーション結果 まず，He 原子ビームを入射による電離シミュレーションを行った。RFP プラ ズマでのビームエネルギーE の違いによる電離位置分布を Fig. 3-4 ，_b Fig. 3-5， Fig. 3-6 に示す。また，この時ビームの拡がりは 2°に設定し，10000 個の粒子 を入射させた。Fig. 3-4 ，Fig. 3-5，Fig. 3-6 を見ればわかるようにビームエネ

ルギーE が高くなるにつれ，より装置内部で浸透しているのがわかる。どのグ_b ラフも浸透長が 0.2 付近でピークと取っているがこれはプラズマ内の電子密度 勾配に比例しているからである。 また，Fig. 3-7 にE の変化におけるシャインスルー率と平均浸透長の関係性_b

1









1 tan 





       2₂ 2 exp  

40 を示した。シャインスルー率とは全粒子数における電離せずに装置壁付近まで 到達した粒子の割合を示したものである。E が高くなるにつれ平均浸透長，シ_b ャインスルー率ともに上昇することが確認できた。 Fig. 3-4 E_b 100　[eV] における電離位置 Fig. 3-5 E_b 1[keV] における電離位置

41

Fig. 3-6 E_b 10[keV] における電離位置