J. Plasma Fusion Res. Vol.95, No.9 (2019)４
２３‐４２６

小特集

統合コードによる磁場閉じ込め核融合
プラズマシミュレーションの現状と今後の展望
Simulation of Magnetically Confined Fusion Plasma
by Integrated Code and Its Prospect

１．核融合プラズマシミュレーションに向けた統合コード開発
1. Development of Integrated Code for Fusion Plasma Simulation
林

伸 彦，福 山

淳１），村 上 定 義１），横 山 雅 之２），藤 田 隆 明３）

HAYASHI Nobuhiko, FUKUYAMA Atsushi１）, MURAKAMI Sadayoshi１）,
YOKOYAMA Masayuki２）and FUJITA Takaaki３）
量子科学技術研究開発機構，１）京都大学大学院工学研究科，２）自然科学研究機構核融合科学研究所，
３）

名古屋大学大学院工学研究科（原稿受付：２
０
１
９年５月２
０日）

核融合プラズマの様々な物理現象が連携する自律的，複雑な現象を扱う有効な手段の１つとして，統合コードが開発，利用されている．本小特集では，主に磁場閉じ込め核融合プラズマの統合コードとそれを構成する物
理モデルの概要，国内外の開発状況，物理解明や実験解析，シナリオ構築への適用例，今後の展望，等について体系的にまとめる．その中で本章は，小特集の最初の章として，統合コードの必要性・構想・方向性について概略する．また，日本の原型炉開発に向けた統合コードの開発工程と各工程の目標，将来の JT-60SA と国際熱核融合実験炉 ITER の実験に対する役割とそれぞれを利用したモデル検証のスケジュール，等について述べる．そして，本小特集で統合コードの詳細を説明するのに必要な第２章から第５章までの全体構成について述べ，次章に繋げる．
Keywords:
Integrated code, magnetic confinement, modeling, simulation

１．
１ はじめにコードは，物理解明や実験解析だけでなく，運転シナリオ

核融合プラズマの物理解明や実験解析には，数値シミュ

構築や炉設計にも活用されている．統合コードは，現在もレーションが重要な役割を果たす．核融合プラズマは，空日々，物理モデル等の改良・追加が進められていると共間的にはプラズマ全体にわたって，時間的には放電時間にに，近年では用途に応じて，第一原理コードそのものや機わたって，様々な時空間スケールの物理現象が相互に影響

械学習を用いたモデルも開発され組み合わされる等，予測を与えて自律的，複雑な挙動を示す．単一の第一原理シ

精度・信頼性の向上が図られ開発が進められている．

ミュレーションコードで異なった時空間スケールの多階層

磁場閉じ込め核融合プラズマの統合コードに関する本学の現象を扱う試みが進みつつあるが，現在でも非常に挑戦

会誌記事としては，２００５年の解説記事［１］で国内における的な課題である．一方，このようなプラズマを扱う有効な

核燃焼プラズマ統合コード構想（BPSI: Burning Plasma

手段の１つとして，磁場閉じ込め核融合プラズマを対象

Simulation Initiative）を中心に，京都大学で開発されていに，プラズマの輸送コードを中心として，MHD（Magneto

る統合コード TASK のモデルや結果について紹介されてHydrodynamics）平 衡 コ ー ド，加 熱・電 流 駆 動 コ ー ド，いる．また，小特集及びプロジェクトレビューの一部としMHD 安定性解析コード，ダイバータコード等，時空間スて，核融合科学研究所（核融合研）で開発されている統合

ケールの特定の階層を扱う各種コードを組み合わせて結合

コード TASK3D について計画の概要や開発の進展，炉設した統合コードが開発・利用されている．この様な統合計ツールとして見た場合の統合コードの役割が述べられて

National Institutes for Quantum and Radiological Science and Technology, Naka, IBARAKI 311-0193, Japancorresponding author’s e-mail: hayashi.nobuhiko@qst.go.jp

４２３!2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

ジュールを結合したものである．統合コードを構成する物

いる［２，
３］．
本小特集では，より広く，統合コードとそれを構成する

理モジュール群の詳細については，第２章で述べる．工学

物理モデルの概要，国内外の開発状況，物理解明や実験解

モジュールについては，紙面の都合上割愛するが，物理モ

析，シナリオ構築への適用例，今後の展望，等について体ジュールに並び重要である．特に，核融合プラズマでは，系的にまとめる．その中で本章は，小特集の最初の章とし外部からの加熱に比べ核融合反応による自己加熱が支配的て，統合コードの必要性・構想・方向性について概略す

であり，プラズマの自律性が強い状況で制御する必要がある．また，関連した日本の原型炉開発に向けた統合コード

り，従来の PID 制御以上の先進制御が必要になることが予の開発工程と各工程の目標，将来の JT-60SA や国際熱核融想される．また，核融合プラズマを確実に制御するための合実験炉 ITER 等の実験に対する役割とそれぞれを利用し計測器や制御コントローラの選定を行う計装制御の検討にたモデル検証のスケジュールについて述べる．そして，本も，統合コードの果たす役割は重要である．
上記の統合コードの構想は国内では BPSI でまとめられ，小特集で統合コードの詳細を説明するのに必要な第２章か

現状主に３つの統合コードの開発が進められている．量子ら第５章までの全体構成について述べ，次章に繋げる．

科学技術研究開発機構が中心に開発している TOPICS［２］，

１．
２ 統合コード開発の現状

京都大学と核融合研と九州大学が中心に開発している

様々な物理・工学モデルを統合したシミュレーション

TASK［３］，名古屋大学が中心に開発している TOTAL［４］

コードは，核融合プラズマ実験の解析ツールとして役立つである．各々の統合コードには他にはないモデルがあり，だけでなく，核融合プラズマで現れる様々な相互結合現象

必要に応じて他の統合コードがもつモデルを導入して協力の理解，核融合プラズマ挙動の予測にも有効である．

した研究を行っている．これらの国内の統合コードの詳細

統合コードに必要なモデルは，以下に例として挙げるよと国外の開発状況は第３章で述べる．さらに，これらの統

うに物理課題だけでなく工学課題も包含しなければならな

合コードを用いた物理解明の実例を第４章で紹介する．

い．

１．
３ 核融合原型炉に向けた統合コード開発の課題−物理課題：
・熱／粒子／運動量それぞれの輸送と閉じ込め原型炉では燃料増倍率を確保するために計測や加熱装置

・新古典輸送物理用のポートの数が制限されると考えられることから，核融・乱流物理

合プラズマの動的挙動を確実かつ正確に予測できるシミュ

・MHD 不安定性

レータ開発が必要になる．全ての重要な物理・工学課題を・高エネルギー粒子物理

網羅した統合コードはシミュレータになりえる．しかし，・粒子補給／加熱・電流駆動

統合コードによる予測の信頼性を向上させるためには，用・スクレイプオフ層（SOL）／ダイバータ，プラズマ−壁相いているモデル・コードの実験による検証が必須である．

互作用

大型ヘリカル装置（LHD）等の既存装置の実験および将来・運転シナリオ（着火からシャットダウンまで）の JT-60SA や ITER の実験を用いたモデル検証を効率良く−工学課題：

行うために，実験家を含めてより多くのユーザーが容易に・外部コイル電流／電圧

運転シナリオ検討に利用できる等の実験との連携を考えた

・電磁力・熱負荷開発・整備が必要である．未だ実験が行われていないプラ・計測モデリングズマの挙動を予測し，その予測を確認するための実験（運・モデルに基づいた実時間制御

転シナリオ）を策定・提案し，データを取得して検証する

・統合制御システム（プラズマ位置・形状，温度・密度等ということを，統合コードの開発プロセスの１つとして行う必要がある．

の分布，MHD 不安定等の同時制御）上記の課題の各々を取り扱えるコードをまず開発・改良

日本の核融合原型炉開発に向けたアクションプランのして，各々の課題に関連する物理機構を理解し予測する必チェックアンドレビュー（C&R）［４］には，「原型炉を見据要がある．世界には既に理論モデルから第一原理コードまえた高ベータ定常プラズマ運転技術の確立」の項目に統合

で様々なレベルのモデル・コードがあり，それらを活用でコードに関する記述があり，図１に示す様に，２０２５年からきる．そして，複数の物理機構が結合した現象を理解して数年以内の第２回中間 C&R に「ダイバータを含む統合シプラズマの自律的，複雑な挙動を予測するためには，上記

ミュレーションの JT-60SA 等による検証」，２０３０年代の原課題を記述するモデル・コードをモジュール化し，それら型炉段階への移行判断で「ITER 燃焼制御の知見を踏まえを統合したコードを開発する必要がある．ただし，各モ

た統合シミュレーションにより，非誘導定常運転の見通し

ジュールには扱える時空間スケールに限りがある等，利用を得る」となっている．この C&R と移行判断をクリアできに制限があることが考えられるので，モジュールの統合にるように統合コードの開発を進めねばならない．そのたあたっては第一原理シミュレーションと実験と比較して注め，原型炉設計合同特別チームの理論・シミュレーション

意深く検証する必要がある．従来の統合コードは基本的にワーキンググループのもとで，統合コードや第一原理コー

は，MHD 平衡における磁気面を横切る１次元のプラズマドを含めた開発計画が議論され，まとめられた［５］．
ITER と原型炉におけるプラズマの予測・制御や運転の輸 送 を 解 く 輸 送コ ー ド を 中 心 に，各 種 物 理・工 学 モ
４２４

Special Topic Article

1. Development of Integrated Code for Fusion Plasma Simulation

N. Hayashi et al.

と簡易化によるシミュレータ構築の２つの要求に適応した
統合コード開発が必要であろう．統合コードで使用するモ
デル開発に関しては，当然，第一原理シミュレーションと実験の知見が必須である．そのため，アクションプランにある第一原理シミュレーションコード開発計画と連携して，第一原理コード自体を統合コードに結合するだけでな
く，第一原理シミュレーションの結果を模擬できる簡易化したモデルの開発も欠かせない．また，実験データに基づ
図１

く経験（比例）則も簡易化したモデルとして活用できるが，統合コードに関わる核融合炉原型炉開発に向けたアクショ
ンプランのチェックアンドレビュー（C&R）の達成目標と，原型炉段階への移行判断，JT-60SA と ITER の実 験 ス ケ
ジュール．

別途，物理解明が必要となる．物理解明やシミュレータの目的に応じて様々なレベルのモデルを開発し使えるようにすべきである．現状でモデルがないものに関しては，まずシ ナ リ オ の 設 計 を め ざ し た モ デ ル・コ ー ド の 検 証 は，

簡易的なもので良いので最低限のモデルを開発し，改良し

JT-60SA 計画［６，７］の主要な目的の１つである．図２に示ていくことを考える必要がある．上記で述べたモデルの高す様に，JT-60SA 実験で，高ベータ・高自発電流割合のプ

度化や簡易化の最近の進展は，第５章で紹介する．

ラズマでの検証を行う．一連の予測・実験策定・検証方法

特定の統合コードの枠にとらわれずに統合コードの開発のサイクルを構築・確立して ITER を用いた研究の効率を協力して進めるためには，モジュールの相互利用が容易

化，原型炉の予測と運転シナリオ策定に活用する．そして，にできるようにする必要があり，そのためにモジュール間

ITER 実験が開始されれば，燃焼プラズマでモデル・コーのデータ交換等のインターフェースを定義する共通のフドの検証を行い，JT-60SA と ITER の両方の検証から，原レームワークが欠かせない．このフレームワークは，モ型炉の予測・シミュレータの信頼性向上を図ることが考えジュールの相互利用だけでなく，複数のモデル同士の比較られている．

検証や実験との比較検証にも役立つ．国内的には BPSI の原型炉に向けて統合コードに求められる役割を担うた

基 に 開 発 さ れ た BPSD（Burning Plasma Simulation

め，現状不足している物理・工学モデルの開発や既存モデ

Data）というデータインターフェースに基づくフレーム

ルの改良等を行い，その統合コードを基に運転シナリオのワークがあり，TASK 系の統合コードを中心に活用されて作成や制御検討等に使うシミュレータを開発する必要があいる．一方，ITER 機構には IMAS（Integrated Modelling

る．物理解明のためにモデルを追加することやモデル自体

and Analysis Suite）［８］というフレームワークの開発が進めの改良を行って，統合コードのモデルを高度化することられており，各国のコードのモジュール化および統合，国は，高い計算機資源を要求し長時間の計算を要するおそれ

際トカマク物理活動（ITPA : International Tokamak Phys-がある．モジュールによっては，スーパーコンピュータのics Activity）におけるシミュレーション・実験データの交

大規模並列計算や GPGPU 等をもつ特定の計算機での実行換に利用されている．これらのフレームワークを活用したを必要とする場合もある．そのために，統合コードの全て

協力が，統合コードの効率的な開発に求められる．上記ののモジュールを同一計算機で動かすことを想定するだけでフレームワークに関しては，第３章で紹介する．

なく，ネットワーク上にある複数の計算機を使って計算資

コードの検証には，精度確認（同じ現象を扱う複数の源の有効利用と分散並列処理による高速化を図る．これにコード間ベンチマーク等）と妥当性検証（コードの結果とは，最新の計算科学の知見の導入が必要であろう．一方，

実験結果との比較）の２つがある．上で述べたフレーム

シミュレータには，運転シナリオ作成や計装制御検討，実ワークを用いることができれば確認・検証を効率的に行う

時間制御等，様々なレベルのものがあるが，概して少ないことができ，コードの不確実性の特定とその定量化を容易

計算機資源で短時間の計算を求められる．そのため，デーに実現できる．そして，実験の測定データを活用してモデ

タベースや機械学習等も用いたモデルの簡易化が求められ

ルとコードを改良することにより，予測の不確実性を最小る．したがって，これからはモデル高度化による物理解明化しつつ，その不確実性を定量化することで予測の信頼性を向上させることができる．

１．
４ まとめ
本章は，小特集の最初の章として，統合コードの必要
性・構想・方向性を概略し，関連した日本の原型炉開発に向けた統合コードの開発工程と目標，JT-60SA や ITER等の実験に対する役割とそれぞれを利用したモデル検証のスケジュールについて述べた．第２章から第５章では，次の構成で統合コードの詳細を説明する．第２章で統合コー
図２

理論モデルとシミュレーションコードによる予測と実験との比較による検証の計画．

ドを構成する物理モジュール群について，第３章で国内外の統合コードの開発状況を紹介する．そして統合コードを
４２５

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019［４］文部科学省核融合科学技術委員会：
http : / / www. mext. go. jp / b ̲ menu / shingi / gijyutu /gijyutu2/074/index.htm［５］原型炉設計合同特別チーム，理論・シミュレーション
ワ ー キ ン グ グ ル ー プ 報 告 書：http://id.nii.ac.jp/1657/
00074926/［６］JT-60SA Research Plan v4.0, http://www.jt60sa.org/
pdfs/JT-60SA̲Res̲Plan.pdf［７］林 伸 彦，福 山 淳：プラズ マ･核 融 合 学 会 誌 88, 678(2012).
［８］F. Imbeaux et al., Nucl. Fusion 55, 123006 (2015).

使った実験解析や物理解明，運転シナリオ構築の実例を第４章で示す．第５章で統合コード開発の今後の展望として，進展著しい情報科学等の最新の知見を取り入れ，従来の枠に囚われない開発について紹介する．最後の第６章で
まとめる．

参考文献
［１］福山 淳，矢木雅敏：プラズマ・核融合学会誌 81, 747(2005).
［２］横山雅之：プラズマ・核融合学会誌 92, 601 (2016).
［３］横山雅之：プラズマ・核融合学会誌 92, 814 (2016).

４２６

J. Plasma Fusion Res. Vol.95, No.9 (2019)４
２７‐４３６

小特集 統合コードによる磁場閉じ込め核融合プラズマシミュレーションの現状と今後の展望

２．統合コードを構成する物理モジュール群
2. Physics Modules Constituting Integrated Code
村 上 定 義，本 多充１），相 羽 信 行１），松 山 顕 之１），林 伸 彦１），星 野 一 生２），
藤 田 隆 明３），福 山 淳，横 山 雅 之４，５）

MURAKAMI Sadayoshi, HONDA Mitsuru１）, AIBA Nobuyuki１）, MATSUYAMA Akinobu１）,
HAYASHI Nobuhiko１）, HOSHINO Kazuo２）, FUJITA Takaaki３）,
FUKUYAMA Atsushi and YOKOYAMA Masayuki４，５）
京都大学大学院工学研究科，１）量子科学技術研究開発機構，２）慶應義塾大学理工学部
３）

名古屋大学大学院工学研究科，４）自然科学研究機構核融合科学研究所，５）総合研究大学院大学（原稿受付：２
０
１
９年５月２
０日）

統合コードは，プラズマ熱・粒子輸送や MHD 平衡・安定性など各要素物理を解く様々な物理モジュールで構成されている．本章では，それら物理モジュール群について，仮定している物理モデルや基礎となる方程式等について紹介する．具体的には，プラズマ熱・粒子輸送，MHD 平衡・安定性，ディスラプション・ペレット，ダイバータプラズマ，波動加熱，NBI 加熱などについてである．また，ヘリカルなど３次元系プラズマに対する物
理モジュールとして，３次元 MHD 平衡，新古典輸送，加熱などについても紹介する．
Keywords:
integrated simulation, plasma transport, impurity transport, MHD (Magnet hydrodynamics) equilibrium, MHD stability,edge plasma, divertor plasma, plasma heating, 3D MHD

２．
１ はじめに２．
２ 巨視的なプラズマの発展を導く輸送方程式群

統合コードにおいては，幅広い時間・空間スケールの物

統合コードで表現されるプラズマの巨視的な時間発展は理現象を考慮して，核融合プラズマの時間発展をシミュ

輸送方程式群によって記述される．歴史的に，輸送方程式

レーションする必要がある．このため，スケールの異なるはプラズマの閉じ込め性能を解析，評価，模擬するために様々な物理現象をそれぞれモデル化した物理モジュール群

導出されてきたため，プラズマの閉じ込め時間程度の長いを結合し，シミュレーションを行う．物理モジュールの多時間スケールにおいて変化する現象に着目している．そのくは，それ単体でも物理解析を行うことができるが，統合ため，トーラスプラズマ中の衝突性輸送を記述する新古典

コードとして結合されることにより，他の物理モジュール理論と密接に関連しながら発展してきた［１］．本節では輸との連携を通して，より現実に近い現象の解析が可能とな送方程式群の構成と，その導出について記載する．

ると考えられる．

プラズマの振る舞いを一般的に記述できる実空間３次元

本章では，統合コードを構成している各要素物理を解く

速度空間３次元の運動論方程式の速度モーメントを取る

様々な物理モジュールについて紹介する．プラズマの発展と，速度 !の次数毎に流体方程式が階層的に得られる．偶を導く輸送方程式群，不純物輸送，MHD 平衡と安定性，過パリティの速度モーメントを取ることで，低い次数から順度現象のモデル，ダイバータモデル，加熱モデルなどにつに連続の式とエネルギー保存方程式が得られる．これらはいて，仮定している物理モデルや基礎となる方程式や具体それぞれ密度と圧力の時間発展を記述する．一方，奇パリ的なコード等について紹介する．また，ヘリカルなど３次ティの速度モーメントを取るとベクトル量の発展方程式，

元系 プ ラ ズ マ に 対 す る 物 理 モ ジ ュ ー ル と して，３次元

すなわち運動量方程式と熱運動量方程式が得られる．これMHD 平衡，新古典輸送，加熱などについても紹介する．紙

らはそれぞれプラズマ流と熱流の時間発展を記述する．厳

面の都合によりそれぞれ紹介は簡単なものとなっている

密にはモーメント方程式は有限のモーメントで打ち切るこが，その詳細については参考文献を参照していただきたとは出来ないが，高次の効果を低次の量で表現するクローい．

ジャを導入することで，適切な次数で打ち切った閉じた方

また，それぞれのモジュールについて実際の適用例につ程式系を構成している．得られた方程式群は速度空間で積いては，４章において紹介されている．

分されたため，空間３次元の方程式となっている．
プラズマの巨視的な時間発展に影響を及ぼす現象のみをKyoto University, KYOTO 615-8540, Japancorresponding author’s e-mail: murakami@nucleng.kyoto-u.ac.jp

４２７!2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

るわけだが，輸送オーダリングから径方向流束は $"3(程度取り出し，それ以外の速い現象を捨象するために，輸送"［２］．ここで，オーダリング *(&$'(，*&$
3(を導入するの大きさしか持たない，という結果が必然的に導かれる．

$はラーモア半径と巨視的な特長長さ L の比で表される十例えば粒子源項を無視した（１）式を見ると，左辺は $"'(/

分小さな値であり，'(はイオン周回周波数，3(はイオン熱であり，右辺は "!
$なので，"&/$"3(が導かれる．すなわ

速度である．輸送オーダリングを適用するということは，

ち径方向速度は熱速度より２オーダも小さく，実際のトー

端的に言えば時間とプラズマ流の速度に関して上記のオーラス装置では m/sを単位としてせいぜいオーダ１程度の大ダを持った量や項のみを方程式に残すということである．

きさしか持たないことがわかる．熱輸送についても同様の輸送オーダリングの適用によって，連続の式とエネルギー帰結となる．

保存方程式において主要項は " ($")である一方，運動量方

先ほど無視した運動量方程式について考察する［３］．運程式などの奇パリティの速度モーメントを持った方程式は動量方程式はベクトル量であるため，各方向に対して射影

一般に無視することができる．詳細は後述する．輸送オーを取った方程式を考える必要がある．径方向の射影を取っダリングによって磁力線方向の運動は磁力線を横切る運動

た場合の主要項から，径電場を決定する際に用いられる径よりも遙かに速いことがわかるため，密度や温度といった方向の力の釣り合い式が導かれる．輸送コードにおいて物理量は磁気面を横切る運動が変化する時間スケールより

は，径電場を評価したいときのみ，この方程式を解くことも遙かに速く磁気面内で一様になる．そのため，磁気面内になる．磁力線方向の運動方程式の主要項から，磁力線方で物理量の平均を取る磁気面平均を行うことは妥当である

向の力の釣り合い式が導かれ，熱流束に対しても同様の方と見ることができる上，ある局所量とその磁気面平均量は

程式を導き連立することで，新古典輸送に関連する諸量を同等であると見なすことができる．ゆえに，得られた空間導くことができる［１，４，５］．これらは輸送オーダリングよ３次元の連続の式とエネルギー保存方程式を磁気面平均すり１オーダ速い現象を扱うドリフトオーダリングで導出さることで，空間１次元となった以下の密度発展方程式と圧れているため，この時間スケールでは密度や圧力分布は時

力発展方程式が得られる．

間変化しておらず，それゆえ輸送オーダリングに基づく時間発展を扱う際には，各時刻で磁力線方向の力の釣り合い! * *./! *( ,)$! . *.","(,
*.*
* *
1
&

（１）

# ! * *.$!#0! *($! . *.&,"+"%,
,)
* *
1
&
"*.$!# *（２）が瞬時に成立しているとしている．このため，新古典輸送
モジュールは各時刻において磁力線方向の連立方程式を解
いて新古典諸量を計算し，輸送コードに渡して時間発展を解いている．

ここで，添え付き下字 ,は粒子種を表し，* は磁気面内のトロイダル方向に射影した運動量方程式を考える．

.

体積（* $*，/は密度，&は小径方向磁気面座標，*!
*&）! * *.. / '2
( , ,* ,)+)1*.*(は粒子源（消滅項も含む．熱も同じ），0は圧力，+ は熱交換項，%は加熱源を表す．"，& はそれぞれ粒子束，熱流

-

'"#
%,+
'#,)+
'")+
$!*
'#
"*
",/
,*
&
（５）
!,#
'(, +
"'(+
"*
,*

束を表し，一般に以下のように書かれる．
/,
"
", %,
/,3,!*
!, *
'&,
,
'&,
+
*
&

-

ここで，'は主半径，2)はトロイダル速度，%は応力テン
ソル，#)は摩擦力，")はトロイダル電場，(はポロイダル（３）%'), $) ""&, %*/,3, ),!*/,#, *
,
'&,
+
,
'&,
+
" , , （４）*& "磁束，(& はトロイダル運動量源である．運動量源が無いと仮定してこの方程式に輸送オーダリングを適用すると，主"ここで，) は温度，*
,
'&,
+と *
,
'&,
+ は磁気面の形状効果を

要項は右辺第２，３，４項となり，これは径方向の新古典粒

現すメトリック，3，3%'は粒子と熱のピンチ速度，!，#子束を表現していることとなる［２，４］．ゆえに，このままでは粒子と熱の拡散係数である．本来であれば，"，&は（１）は輸送オーダリングの時間スケールにおけるトロイダル角（２）式に対応する奇パリティのモーメント方程式，すなわ

運動量密度の時間発展を記述できないことになる．しかち運動量方程式と熱運動量方程式の従属変数であり，それし，文献［６］において解決策が示された．熱プラズマを構成

らを含めて連立して解き決定する必要がある．しかし輸送する全粒子種に対して（５）式の和を取ることで，粒子間摩

オーダリングによってそれらの方程式を無視しているた擦項である右辺第２項が消え，準中性条件から右辺第３項

め，フラックスを移流拡散形式で仮定して表現することでが消え，新古典粒子束の両極性から右辺第４項が消えるた

方程式を閉じている．ピンチ速度や拡散係数は，輸送コー

め，先ほどの主要項は全て消える．その結果より高次の効ドに結合された外部の新古典輸送モジュールや乱流輸送モ

果を持つ項間の釣り合いが成立し，全トロイダル角運動量

ジュールで評価された各係数の線形和を取ることで表現さ密度を!$! , .,/,'2,)と定義すると，
（５）
式は以下の

れるのが一般的である．

ようになる．

（１）（２）式は密度と圧力の巨視的な時間発展方程式であ! * *.
'"#%,+'(,& +
( !)
$!! *
'#
"! *
*.*
1
,
,

るが，右辺が零と仮定すれば，両式はそれぞれ局所毎の粒子数の保存と，圧力の断熱変化を表している．実際には右辺の拡散項などによって密度や圧力がゆっくり時間変化す（６）（６）式は，初期に回転速度を与えるかもしくは運動量入力
４２８

Special Topic Article

2. Physics Modules Constituting Integrated Code

S. Murakami et al.

を与えない限り，トロイダル回転は生じないことを示してより行われる．その際，対象とするすべてのイオン種と電いる．右辺第１項は（３）（４）式と同様，移流拡散型のモデ子についての連立方程式を解くことになり高原子番号の核ル表記で記述するが，（６）式の詳細も含めた具体的な表式種が存在する場合は計算コストが高くなるため，その対策は４章に譲る．

として異なる価数の同一核種イオンをまとめて取り扱う"reduced charge state"モデルが用いられている
［１０，
１２］．

輸送オーダリングの時間スケールで発展するもう一つの重要な量は，ポロイダル磁束 %である．ポロイダル磁束の

トロイダル回転するプラズマにおいては，プラズマとと時間発展は，プラズマが持つ抵抗によって電磁的に緩和しもに回転する座標系で観測すると遠心力が発生することにていく現象であり，電流のしみ込みなどとも関連する重要なる．遠心力は質量に比例するため質量の大きい不純物イな量である．新古典的な扱いのみで各係数が決定できるこ

オンでは遠心力の影響を大きく受ける．静止している座標ともあって，通常の輸送コードはポロイダル磁束の時間発系ではプラズマ回転に伴う径電場が遠心力と等価な作用を展を司る磁場拡散方程式を必ず解いている．紙面の都合する．この遠心力の効果により重いイオンは弱磁場側で高上，詳細な導出は［７］に譲る．

い密度をもつようになり密度分布にポロイダル非一様性が生じる．イオン密度分布のポロイダル非一様性の効果を新

２．
３ 不純物輸送

古典輸送に取り入れる試みがなされている［１３‐１６］．ドリ

不純物イオンの輸送解析において燃料（水素同位体）イ

フト運動論方程式を数値的に解いて新古典輸送を求める

オンの輸送解析と異なる点として以下の３つを挙げること

NEO コード［１７］では，遠心力の効果も含められている．以ができる：（!）イオンの電荷数が変化する，（"）電荷数が上の理論あるいは数値計算ではイオン軌道の磁気面からの大きいため燃料イオンに比べて衝突周波数が大きく新古典

ずれが考慮されていない．不純物イオンの旋回中心の軌道

輸送がより重要となる，（#）質量数が大きいため燃料イオ

解析に基づき，軌道上でのイオンの電離・再結合にともな

ンに比べてトロイダル回転の影響（遠心力の影響）が大きうピンチ（PHZ）および径電場により変形した軌道上でのい．これらの特徴は特にタングステンなど高原子番号の元

クーロン衝突に伴うピンチ（Er ピンチ）のモデルが提唱さ素について顕著となる．

れ［１８］，TOTAL コードに組み入れられている［１９］．

ある核種の #価のイオンの密度を $#とする（簡単のため

２．
４ ２次元 MHD 平衡核種を示す添字は省略する）と，その時間発展は以下の（７）式で与えられる．
! & "&$! & "&!
#$$$#$ #%#! & $
!$
#%
#""
#%
"&&
" &
%
$

磁場閉じ込め核融合プラズマのシミュレーション研究を行う際には，プラズマが電磁流体 力 学（MHD: Magneto
Hydrodynamics）的に平衡状態にある（MHD 平衡）ことを（７）

前提としており，この MHD 平衡の安定性や輸送特性など

"#
$$
$$
#!!$
#!!""
#"!$
#"!"!
#

をシミュレーションで調べることになる．トカマクプラズ
マの場合，軸対称であるという特徴を用いることで理想

ここで，!# は磁気面を横切る #価のイオンの粒子束，$$MHD 方程式から Grad-Shafranov 方程式と呼ばれる平衡方は電子密度，#"!価のイオンへの電#は#価のイオンから#

程式を導出することが可能であり，この方程式を数値的に

離速度係数，"# は #価のイオンから #!!価のイオンへの解くことでポロイダル断面の２次元 MHD 平衡を求める

再結合速度係数，!# はソース項である．電離・再結合過程［２０］．この方程式はプラズマの圧力（勾配）分布およびトとして電子衝突によるもののみを考えている．

ロイダル磁場関数あるいは電流分布を与えることで解くこ

不純物イオンの存在により，線スペクトル放射，再結合とができ，これらの分布を持つ平衡状態を求めることが可

放射，制動放射などで放射損失が発生する．それらの計算能である．また，輸送方程式との親和性の観点から，トロに必要な放射損失係数および上述の電離・再結合速度係数

イダル磁場関数あるいは電流分布の代わりに安全係数分布を求める方法としては，ADPAK コード［８］，ADAS コー

を与えて FCT（Flux Conserving Tokamak）平衡と呼ばれド／データベース［９］などがよく利用されている．

るものを求めることもある．量研にて開発が進められてい%#粒子束 !#は，新古典輸送によるもの !#に異常輸送によるトカマク統合シミュレーションコード TOPICS では必要"%るもの !#
を加えて評価される．不純物イオンは燃料イオに応じていずれの手法でも平衡状態を計算することが可能

ンよりも衝突周波数が高いため，新古典輸送が大きく異常である．２次元 MHD 平衡は，上記２．
１，２．
２節の１次元輸送輸送の寄与は相対的に小さい．衝突周波数の違いのため，方程式の磁気面座標を与えるとともに，輸送方程式と組み

燃料イオンがバナナ領域にある場合でも不純物イオンは合わせて解くことでプラズマの時間発展を評価する 1.5 次

Pfirsch Schlüter 領域にあることが多い．

元輸送コードとなる．

圧力勾配による反磁性ドリフト速度が電荷数に反比例す

２．
５ MHD 安定性るため，燃料イオンと不純物イオンのポロイダル流速が異なってポロイダル方向の摩擦力が生じ，不純物イオンの内

本節では，TOPICS で利用可能な MHD 安定性解析コー

向きピンチが駆動される．トーラスプラズマにおいては新ド・モジュールについて紹介する．TOPICS には，長〜中古典理論により粒子束が求められる．標準的な新古典輸送

間 波 長 理 想 MHD モ ー ド 安 定 性 解 析 コ ー ド MARG2D

理論［４］に基づく輸送係数の計算は NCLASS
［１０，
１１］等に

［２１，
２２］，短波長理想バルーニングモード安定性解析コー
４２９

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

ド BETA［２３］，抵抗性壁モード（RWM）安定性解析コー

組み合わせることで，トカマクプラズマにおける理想

ド MINERVA/RWMaC［２４，２５］，および新古典テアリング

MHD 安定性は実用上すべて考慮されたシミュレーション

モード（NTM）による磁気島幅評価モジュール［２６］が実装が可能である［２７］．

されている（表１）．これらは相互補完的な関係にあり，

理想 MHD モードの発生を回避することで，統合コード

TOPICS の利用者の目的に応じて必要なコード・モジューの計算で得られる放電シナリオが放電停止に至るかどうか

ルのみを有効にして計算を行うことが可能である．本小特を事前にある程度検討することが可能になったが，放電停集では紙面の都合上，統合コードへの実装として最も特徴

止の原因は理想 MHD モードの発生に限られない．MHD的である，MARG2D および MINERVA/RWMaC について

モードに原因を限ったとしても，RWM はモードの成長は紹介する．

遅いものの大域的なプラズマ変形を引き起こすため，放電

トカマクプラズマにおける理想 MHD モードは，大域的停止の原因になる．RWM は，プラズマを囲む真空容器なな MHD モードの不安定化による放電停止や高閉じ込め運

どに有限の電気抵抗が存在するために，プラズマ変形を妨転モード（H-mode）におけるエッジローカライズドモード

げるように真空容器などに流れる渦電流が失われることで（ELM）の発生などの主要な原因であり，プラズマが到達発生する．そのため，外部制御コイルを用いて失われた渦できる最大ベータ値やプラズマ電流量を制限する．そのた電流が作るはずだった磁場を補うことや，プラズマを回転

め，運転シナリオの検討や H-mode 実験・高ベータプラズ

させて渦電流が失われることを妨ぐことで，RWM は安定

マ実験の解析などを行う際には，同モードの安定性解析を化できることが知られている［２８，２９］．このうち，後者の回行うことは非常に重要である．

転 に よ る 安 定 化 を 考 慮 し た RWM の 安 定 性 解 析 は

統合コードでの安定性解析を実現する際に求められるの

MINERVA/RWMaC コードを用いることで可能である．

は．プラズマの安定・不安定の正確な判定を短時間で行う

このコードは，プラズマ回転を考慮した理想 MHD 安定性ことである．トカマクプラズマの線形理想 MHD 安定性解

解析コード MINERVA［３０］と，真空容器の電気抵抗の影

析では，プラズマが不安定である場合には数値計算で不安響を考慮した真空磁場計算コードRWMaCを組み合わせた定モードを同定することは比較的容易であるが，安定であ物である．前述のMARG2Dとは異なり，プラズマが安定でることを陽に示すことは容易ではなく，通常は不安定ある場合にはある程度の計算時間をかけて RWM 安定性評

モードが見つからない＝安定であるという間接的な方法

価を行う必要がある．しかし，RWM はプラズマを囲む真で示すにとどまる（詳細な理由は［２１］などを参照）．その空容器が存在しない（壁無し）と仮定し た 場 合 に 理 想ため，様々な計算条件で安定性を評価した上で不安定モー

MHD モードが不安定である状態ではじめて発生する（こドが見つからないことを示さなければならず，短時間で判の理想モードの発生は真空容器を設置した場合には渦電流定を行うことは困難である．さらに，運転シナリオを検討で抑えられるが，容器の電気抵抗によって完全には安定化する上で通常想定するプラズマは理想 MHD モードに対しできずにRWMが発生する）．そのため，この壁無し状態でて安定であるため，このような比較的長時間の計算が必要

理想 MHD モードが不安定になるまでは MARG2D で高速となる安定性解析を複数回繰り返さなければならず，統合に安定性を評価し，不安定になる状態と判定された場合のコードに理想 MHD 安定性解析を組み込む上で大きな障害

み MINERVA/RWMaC を使った RWM 解析を行うことで，となっていた．MARG2D では，数値解析の目的をプラズマ

計算時間の大幅な削減が可能である．

の安定・不安定を 判 定 す る こ と に 制 限 す る こ と で 理想

２．
６ 過渡現象のモデル（ペレット入射，ディスラ
プションなど）MHD モデルの一部を解析的に解き進めて得られる Newcomb 方程式の随伴固有値問題を，MPI ライブラリを用いた並列計算で数値的に解くことで上記の問題を解決してい２．
２節で記述されたように統合コードではプラズマの巨る［２１］．これにより，統合コードによるプラズマの時間発

視的な時間発展を記述する目的で輸送オーダリングを採用

展シミュレーションにおいて，大域的な MHD モードからし，それ以外の速い現象は輸送方程式の境界条件やソース

ELM までを解析対象とした理想MHD安定性解析が可能と項，もしくは輸送係数の急速な変化としてモデル化されなった［２２］．この MARG2D の解析 結 果，お よ び 前述のる．本節ではそのような例として，将来の核融合プラズマ

BETAコードを用いて短波長理想バルーニングモード解析

において燃料となる重水素と三重水素をプラズマ中に供給する手法として現時点で最も有望な手法である固体水素ペ表１TOPICS に実装され て い る MHD 安 定 性 解 析 コ ー ド・モ
ジュールの一覧．

レット入射［３１］を考える．水素ペレットはプラズマ周辺部
から中心に向かって侵入し，入射経路に沿って各磁気面にコード・モジュール名
MARG2D役割
長〜中間波長理想 MHD モード安定性解析
BETA
短波長極限理想バルーニングモード
安定性解析
MINERVA/RWMaC
抵抗性壁モード（RWM）安定性解析
NTM 磁 気 島 幅 評 価 モ 新古典テアリングモード（NTM）にジュール
よって生じる磁気島幅の評価

中性水素粒子の雲を生じる．中性粒子雲は電離して磁場に捕捉され，磁力線方向に音速程度で膨張すると同時に背景プラズマの磁場勾配によって磁場を横切る向きのドリフトを生じながら，背景プラズマに均質化する．トカマク装置
では溶発雲のドリフトはバルーニングモード的な振る舞いで書けるのでトーラス外側向きのドリフトが支配的であるが，ヘリカル装置ではインターチェンジ的な振る舞いも加４３０

Special Topic Article

2. Physics Modules Constituting Integrated Code

S. Murakami et al.

わり［３２］，密度デポジションのピーク位置は単純な描像でえば，鋸歯状振動のモデルは，不安定性発生時に，不安定は決まらず詳細な計算が必要となる．溶発雲のドリフトを性に影響される領域の圧力・密度分布を平均化して中心か計算するコードとしては HPI2 コード［３３］が最も多くの物

らのプラズマの吐き出しを模擬している［４２］．もう１つ

理効果を考慮しており，トカマク実験とヘリカル実験でのは，輸送コードの枠組みをそのまま使って輸送を増加して検証が進んでいる．統合コードにおいては固体水素ペレッ分布を変化させるもので，テアリングモードや ELM のモトの溶発分布は中性ガス遮蔽（Neutral Gas Shielding）モデ

デルで用いられている．テアリングモードは磁気島を形成

ル［３４］を 用 い る と 高速 に計 算 で き る の で，こ の 分 布 にしてその内部の分布を平坦化させることから，１次元の輸HPI2 コード等で計算したデポジションのシフトと拡散的送コードで磁気島のある有理面近傍の領域の拡散係数を増な広がりを考慮することで粒子デポジション分布を決定

幅させて分布を平坦化して炉心のプラズマ閉じ込め劣化をし，輸送コードのソース項に組み込むことができる．核融

模擬する［２６］．ELM モデルは，拡散係数か対流速度を増合原型炉設計では HPI２による溶発雲のドリフト計算やペ加してペデスタルの分布崩壊を模擬し，ペデスタルに蓄積

レット入射の粒子デポジション分布を用いた核燃焼制御のされたエネルギー・粒子を周辺プラズマに吐き出す．ELM

シミュレーションが行われ，燃料供給や粒子循環の検討がによる電流分布の変化もプラズマ抵抗を増幅して模擬す行われている［３５］．

る．ELM を引き起こす不安定性（ピーリング・バルーニン

他方，統合コードの過渡現象の扱いで注意しなければな

グモード）の判定は計算時間の短い簡易的な手法が多く採らないのは MHD 平衡の扱いである．普段は意識すること用されていたが，２．
５節で述べたように線形 MHD 安定性は少ないが，本来の意味での 1.5 次元輸送モデルでは MHD

コード MARG2D の高速化により統合コードの分布時間発オーダリングの現象は無視されているわけではなく，時間

展の時刻毎の判定が可能となった．図１は，TOPICS でペ

微分項，すなわち慣性項を落とした MHD 平衡方程式を境

デスタル領域の輸送を新古典レベルに低減してペデスタル

界値問題として解くことで方程式の体系に巧妙に取り込まを形成し，不安定性が発生した時点で，不安定モードによれて い る［３６］．こ の た め，MHD 現 象 や L-H 遷 移，ELM

る径方向摂動の固有関数分布に基づいて拡散係数を増幅しなど平衡分布を大きく変える現象のモデル化にあたってはた一例である［４３］．この ELM モデルに関しては，MHD

統合コードを第一原理的な方程式系に立ち返って考え，時の速い現象をそれより長い時間スケールで構築した輸送

間スケールの分離を注意深く考察することが重要である．

コードの枠組みで扱っているので，実験と非線形 MHD シ

このような問題の例はトカマク放電のディスラプション時

ミュレーション等との比較から注意して検証する必要があにおいて観測される垂直移動現象（Vertical Displacement

る．さらに，この ELM モデルでは増幅する拡散係数の最大

Event: VDE）である．ITER の VDE シミュレーションに用値か，その維持時間がフリーパラメータであり改良の余地いられている DINA コード［３７］では 1.5 次元輸送モデルでがあるが，第４章で記述する様に他の物理（ブートスト

VDE を扱うため，輸送方程式の全時間ステップで MHD

ラップ電流や SOL の磁力線方向輸送）で実験結果を再現

平衡を解き直し，輸送のメトリックと矛盾がないよう収束

し，その物理の一部（磁力線方向輸送）は非線形 MHD シ

計算を行うことでプラズマ形状の時間発展を追跡しておミュレーションでも実験を再現する上での重要性が確認さり，量研においても同様の思想に基づく ETA コードを開

れている．

発している
［３８］．このようなアプローチが成り立つのは，

２．
８ ダイバータモデル

垂直位置不安定性に対し壁安定化が効いたトカマク放電では VDE の成長率がアルフヴェン時間よりも十分小さいた

スクレイプオフ層（SOL）及びダイバータ領域（以下でめである．DINA と並ぶ VDE コードとして有 名 な TSC

特に断りがなければ，これらの領域をまとめてダイバータではこの性質を逆手にとってプラズマの実効質量を人為的と呼ぶ）では，磁気面は開いており，磁力線は固体壁（主に増加させた MHD 方程式を解くことで比較的小さい計算にダイバータ板）に当たる．ダイバータにおけるプラズマ，

量で VDE を計算することに成功している［３９］．トカマク

中性粒子，不純物，壁の間で起こる強い相互作用の結果，ディスラプションの最近の研究では，このようなアプロー
チをキンク変形した３次元のプラズマに適用できるか，という問題が検討されている［４０］．ヘリカル配位の３次元理想 MHD 平衡［４１］
とも関わりが深く，統合コードをより高度に発展させるための興味深い研究課題である．

２．
７ MHD 不安定性による輸送モデル
前節２．
６で述べられた過渡現象を含めて MHD 不安定性

図１の影響をプラズマの輸送に反映するため，様々なモデルが考案され使われてきた．大きく分けて２種類あり，１つは，輸送の時間スケールより MHD 不安定性の時間スケールが短いことから MHD 不安定性が発生した時刻で不安定性が影響する領域のプラズマ分布を変化させるものである．例
４３１

MHD 安定性コードの不安定モードに基づく ELM 輸送モデ
ルの説明．ペデスタル領域で輸送を新古典レベルに低減（（c）の“pre ELM”）し て ペ デ ス タ ル 圧 力 分 布 を 形 成（a）
し，その形成の時間発展の時刻毎に線形 MHD 安定性コードで安定性を解析して不安定な場合，不安定モードの径方
向摂動の固有関数分布（b）に基づいて拡散輸送を増幅する（（c）の“post ELM”）．

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

プラズマは径方向だけでなく，磁力線方向にも物理量が大イバータ板ではシース入り口で境界条件を与える．シース

きく変化する．また，中性粒子は磁力線にとらわれないた

理論に基づきイオン流速についてはボーム条件を課し，エめ，中性粒子の平均自由行程と，幾何形状（磁力線が壁に

ネルギー輸送についてはシースの熱伝達係数を用いて対流入射する角度やダイバータの形状・大きさ等）との相互関

フラックスとして与える．ダイバータ以外の壁側境界につ係が重要となるため，定量評価のためには本質的に２次元

いては，基本的には内外ダイバータ板が接続される磁気面

もしくは３次元で取り扱うことが必要になる．本節では，に境界を設置する（第１壁によって切断されない磁気面．

ダイバータモデルの中でも，特に統合ダイバータコードと外側赤道面の SOL でセパラトリクスから 3〜5 cm，プライ呼ばれる SONIC［４４］や SOLPS［４５］等を例に，そのモデル

ベート領域で数 cm）．密度・温度については勾配長を，流概要を説明する．詳細については，各コードの文献や文献

速についてはゼロを与える場合が多い．

［４６‐４８］を参照していただきたい．

中性粒子については，様々な原子分子過程や，壁との相統合ダイバータコードでは，トロイダル軸対称の仮定の互作用，幾何形状効果を考慮しやすいモンテカルロ法を用

もと，沿磁力線方向の輸送をポロイダル断面に投影し，磁いたテスト粒子モデル（MC モデル）が採用され，プラズマ

気面に沿った方向とそれに垂直な径方向の２次元空間でモ流体モデルと結合されている．MC モデルでは，プラズマデルを構築する．ポロイダル断面上での磁場配位や壁形状

流体モデルで計算された２次元プラズマ分布のもとで，プも考慮して，図２のような数値メッシュで解析を行う．最ラズマリサイクリング，体積再結合，ガスパフ等によって

近では，３次元磁場構造を考慮できる３次元モデルの開発生した中性粒子を，電離もしくは排気されるまで軌道追発・利用も進んでいる［４９］．

跡する．この際，電離，荷電交換，弾性衝突，解離等の原プラズマについては，Braginskii の流体方程式［５０］を基

子分子過程をモンテカルロ法により模擬する．衝突に伴うにした，イオン密度，イオンの運動量，イオンエネルギー，

運動量やエネルギー変化は各空間メッシュでスコアリング

電子エネルギーの保存式を解くことで，イオン密度，沿磁され，プラズマの粒子・運動量・エネルギーのソース項と力線方向のイオン流速，イオン温度，電子温度を求める．

してプラズマ流体モデルへ引き渡される．

SOLPS コードの様に，"!! ドリフトや反磁性ドリフトを

不純物は固体壁の損耗により発生する．またダイバータ

考慮する場合には，電流連続の式（電子の運動方程式）を

熱負荷低減のために入射される希ガスも不純物となる．こ

基礎方程式に加え，電位も同時に求めている．沿磁力線方のような不純物の輸送は，SOLPSをはじめとする多くのダ向の輸送についてはクーロン衝突に起因する古典的な輸送

イバータコードでは流体モデルによりプラズマと同様に扱係数を用いる．一方，磁力線を横切る輸送については異常われる．一方で，SONIC では運動論効果や固体壁との相互輸送を仮定し，入力パラメータとして粒子及び熱拡散係数

作用を考慮しやすいMCモデルIMPMCが採用されている．

を与える．現状では，経験的な値もしくは赤道面においてIMPMC では炭素等の比較的質量の軽い不純物を主な対象

特定の実験データを再現できるように空間分布を含め調整としており，不純物イオンのラーモア半径は十分に小さいされたパラメータが用いられている．境界条件は，炉心プと仮定し，旋回中心近似を用いている．摩擦力を含むクー

ラズマ側，壁側（SOL 端及びプライベート領域端），ダイ

ロン衝突，及びクーロン衝突と温度勾配に起因する熱力，バータ板で与えられる．炉心プラズマ側では，測定データ

電離・再結合過程を考慮して，磁力線に沿った旋回中心のを基にした境界値の設定のしやすさや，中性粒子や不純物運動を追跡する．しかし，ITER や原型炉では固体壁にタによるポロイダル分布の非一様性を考慮するために，セパ

ングステンが用いられる計画であり，タングステンは質量#"!
!
#"前後）に設定される．ダラトリクスのやや内側（$"が 大 き い た め ラ ー モ ア 半 径 が 大 き く な る．そ の た め，
IMPMC のような旋回中心近似ではなく，ラーモア旋回運動まで含め全軌道を追跡する IMPGYRO［５１］のようなコードが国内で開発されている．
以上のプラズマ，中性粒子，不純物のモデルは，プラズ
マ分布とそれに対応するソース分布を交換する形で結合され，相互矛盾のない解を求める．流体モデルの進行に合わせて MC モデルの計算を都度行うため，計算コストが膨大になる．条件やモデルによっては１ケ ー ス の 収 束 に
数ヶ月〜半年近い計算時間がかかることもある．SONIC
では開発初期から MPI による並 列 化 を 導 入 し て お り，
500〜1000 PE を用いて数日〜数週間程度で収束解を得る
ことができる．MPI-OpenMPによるハイブリッド並列化やMPMD を用いた並列化効率の改善等，さらなる高速化も進められている．しかし，依然として，運転シナリオ開発
や炉設計研究のような多数のパラメータサーベイには不向

図２きである．そのため，後述するような簡易モデルや，統合

統合ダイバータコードのメッシュ例．ここでは，SONICで用いられる JT-60U のメッシュ［４８］を示す．

ダイバータコードによる解析結果から作成した比例則等，
４３２

Special Topic Article

2. Physics Modules Constituting Integrated Code

S. Murakami et al.

ある程度精度を犠牲にして計算コストを低減させたモデしている．中性粒子による粒子ソースも D5PM のモデルとル・方法を併用して，研究を進めていくことが必要であTOPICS のモンテカルロコードを結合して粒子バランスを考慮したモデリングを行っている［４３］．一方，TOPICS

る．

にはダイバータコードSONICとも結合されており，図４の２．
９ ダイバータの簡易モデル

ように１点で結合するのではなく炉心と周辺の各々のコー

炉心の周辺領域にあるダイバータプラズマは，炉心からドで解く領域をある程度重複させている［５５］．TOPICSの熱・粒子流を受けて変化するとともに，炉心端のセパラ

はセパラトリックスまで１次元分布を解き，炉心端の重複

トリックスの温度・密度の境界条件として炉心プラズマへ

領 域 炉 心 側 の 熱・粒 子 束 を 計 算 し て SONIC に 与 え，

影響を与える．特に，炉心端のセパラトリックス内側のSONIC は重複領域を含んだ２次元分布を解いてセパラト

H-mode ペデスタルの形成やそこで発生するELMとの相互リックスの温度・密度を TOPICS の境界条件として与え作用を考慮するためには，ダイバータプラズマを含めた動る．これにより炉心端の２次元効果をSONIC側で考慮する的現象の統合モデリングが必要である．前節２．
８のダイことが可能となる．重複領域で炉心の輸送係 数 や NBI

バータコードを炉心輸送コードに結合するのは，定量評価（Neutral Beam Injection）等の粒子ソースは TOPICS が評のために必要だが，計算負荷の観点から多数の計算を行う

価し，セパラトリックス外からの中性粒子ソースと放射損パラメータランは難しい．ダイバータコードで求めた定常

失は SONIC が評価して情報交換する．ただし，TOPICS-解のデータベースを構築し，比例則を作成して炉心輸送

SONIC 結合は，D5PM のように重複領域を設けないでセパ

コードの計算で用いる手法［５２］もあるが，動的現象には適ラトリックスで結合したこともあり，数値的安定性と物理用できない．そのために，定量性をある程度犠牲にして，の両観点から結合手法の検討が必要である．炉心と周辺の定性的にでも物理機構を解明し，パラメータランで運転領

両者のモデルが流体モデルであれば上記の方法が使える

域の見通しを得るための簡易モデルが必要となる．従来，が，片方が粒子モデルで領域境界において粒子の出入りがダイバータコードの基礎方程式を磁力線方向に積分し，上ある場合は，粒子モデルでは境界から出て行く粒子束しか流の SOL と下流にあるダイバータプラズマの関係性を調評価できないので工夫が必要である．TOPICS では，炉心べる２点モデルが，物理解釈や実験解析に用いられてきたの不純物輸送コード IMPACT と SONIC の不純物粒子モン［５３］．それを磁力線の両端にある２つダイバータを考慮しテカルロコード IMPMC とを結合するためのモデルを検討て拡張し，図３のような周辺プラズマの静的平衡の多重解中である．

を見つけるために５点モデルが開発され，後に，初期値問

２．
１１ 中性粒子ビーム入射加熱，波動加熱，速度分布解析

題を解く動的５点モデル（D5PM: Dynamic 5-Point Model）
が開発された．動的現象を取り扱うために，イオンと電子の温度を分離し，粒子コードに基づいた運動論効果のモデ

核融合プラズマのさまざまな加熱法の中で，ジュール加ル（heat flux limit）を導入している．中性粒子はイオン化熱と断熱圧縮加熱は通常輸送方程式に含まれる．それに対の平均自由行程とダイバータ長の比を考えたリサイクリン

して，中性粒子ビームや核融合反応生成粒子等の高速粒子

グモデルを用い，イオン化によるエネルギー損失を経験的による加熱や外部から励起された波動による加熱は，加熱に増倍して不純物の放射損失を考慮する［５４］．

手法に応じて高速粒子の生成あるいは波動の吸収が評価される．簡易的には，流体的取り扱いやマクスウェル速度分

２．
１０ 炉心−ダイバータ結合モデル

布に対する波動の吸収等から，加熱パワー密度，駆動電流

炉心とダイバータのプラズマの相互作用を矛盾なく考慮

密度，回転トルク等が求められるが，より正確には衝突にするために両者を結合するモデルが必要である．炉心とダよる緩和や粒子損失等との釣り合いによって定まる速度分

イバータの両者のモデルが流体モデルであれば，一方は

布関数の変形から求めることが必要となる．そのため，加熱・粒子束，もう片方は温度・密度を相互に与えて解け

熱・電流駆動モジュールの機能は，高速イオンの生成，波る．上記の D5PM は統合コード TOPICS に結合され，セパ動の励起・伝播・吸収，速度分布関数の３つの解析に分け

ラトリックスで TOPICS が熱・粒子束を与えてD5PMがセパラトリックス温度・密度を解き，それを境界条件としてTOPICS が炉心内部を解くことを時間発展で陽的に繰り返

図４
図３

ダイバータ５点モデル（D5PM）の概念図．

４３３

統合コード TOPICS におけるダイバータ統合コード SONICとの炉心−ダイバータ結合モデル．

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

る位置やエネルギーを評価することもできる．波動との相ることができる．
中性粒子ビーム入射による高速イオンの生成は，入射中互作用は，i）サイクロトロン共鳴を通過する際に受ける速性粒子ビームの軌跡に沿って，場合によっては粒子のエネ

度変化を乱数として与える手法と ii）準線形理論から求め

ルギー構成を考慮に入れて，電離（励起を伴う多段電離をられる速度拡散と減速を満たす乱数を速度変化として与え含む）と荷電交換の衝突断面積から求めることができる．

る手法があり，後者の方が必要計算量は少ないとされてい高速イオンの生成の解析は，生成後の高速イオンの振る舞る．フォッカープランク方程式を用いた解析では，トカマ

いを記述する OFMC［５６］や NUBEAM［５７］等のコードに含クの軸対称性と軌道平均を用いた３次元解析，さらに空間

まれることが多いが，FIT3D［５８］のように独立して，ドリ輸送を無視した速度空間２次元の局所解析が行われていフト軌道平均された生成分布を求めるコードもある．

る．代 表 的 な 軌 道 追 跡 コ ー ド と し て は OFMC［５６］，核融合反応による高速イオンの生成は，簡易的にはマクGNET［６８］，NUBEAM［５７］等があり，フォッカープラン

スウェル速度分布を仮定して評価できるが，燃料イオンがク コ ー ド と し て は TASK/FP［６９］，BAFP［７０］，CQL3D

高速成分をもつ場合には速度分布関数を用いてビーム・バ［７１］等がある．

ルク成分による反応率を計算する必要がある．特に重水素

２．
１２ ３次元平衡・安定性の中性粒子ビーム入射加熱と三重水素のイオンサイクロト

平衡は，統合解析のみならず，様々な解析の根幹をなす

ロン共鳴加熱が共存する場合はビーム・ビーム成分の反応重要な基盤であるが，本節では，LHD 実験解析のために開も考慮に入れる必要がある．

発，運用されている統合輸送解析スイートTASK3D-a［７２］波動の励起・伝播・吸収の解析には，幾何光学近似を用
いて波束を追跡する光線追跡法やビーム追跡法と，マクス（４章参照）に話題を限って概説する．

ウェル方程式を周波数領域の境界値問題として解き，波動TASK3D-a では，トムソン散乱計測による電子温度分布の空間構造を求める波動光学的手法がある．幾何光学的手［７３］が得られているタイミングに対して，座標マッピング

法は波長の短い電子サイクロトロン波や低域混成波に用いシステムである TSMAP［７４］を介して３次元平衡を規定しられ，波動光学的手法は波長がプラズマサイズと同程度にている．TSMAP は，事前に構築された VMEC［７５］平衡なるイオンサイクロトロン波に用いられてきた．電子サイ

データベースから，計測電子温度分布に最も適した（電子クロトロン波の吸収の解析においては，相対論効果を取り

温度分布のピーク位置を磁気軸と考え，その磁気軸をはさ入れると速度空間におけるサイクロトロン共鳴曲線が楕円

んだトーラス内側と外側の分布対称性を満たす）平衡を選になる効果を取り入れることが必要である．代表的な幾何

び出し，その平衡データを用いて，実座標から実効小半径

光 学 的 解 析 コ ー ド と し て，EC-HAMAMATSU［５９］，へとマッピングを行うものである．これにより，実座標

TASK/WR［６０］，TORAY［６１］，GENRAY［６２］，TRAVIS（視線）に依存した計測データを，実効小半径の関数として［６３］等がある．低域混成共鳴周波数帯の遅波および速波に表すことが可能となるが，データベースで仮定したプラズ

ついて，導波管列アンテナや進行波アンテナによる励起の解析に通常は解析的モデルが用いられるが，アンテナ前面

マ小半径と，実験における蓄積エネルギーの９９％を含むプ
ラズマ小半径位置（"##）との差異を考慮して，"##に対応すの密度分布や密度揺動の影響を取り入れるため，波動光学

る磁気面を最外殻とした VMEC 計算を改めて行っている．

的手法との比較が進められている．イオンサイクロトロン

この再計算された VMEC 平衡は，TSMAP で選び出された周波数帯の波動光学的伝播解析にプラズマの有限ラーモア

VMEC 平衡とは完全には一致しないが，その差異は，統合

半径効果を取り入れると，一般には空間的な畳み込み積分

輸送解析における１次元輸送計算に実質的な影響を及ぼさを含む微積分方程式を解く必要がある．そのため，!）冷たないことが確かめられている［７６］．

いプラズマ近似で波数を評価（##"!!），"）空間の２階

しかし，高ベータプラズマや摂動磁場印加実験など，高微分までの近似（##"""），#）局所的な一様性（"!
!""）

精度の平衡が必要な解析を行う際には，入れ子状の磁気面を仮定してフーリエ展開，$）相互作用の空間局在性を用を仮定しない HINT2［７７］で平衡を構築しているのが現状

い て 積 分 等 の 手 法 が 導 入 さ れ て い る．具 体 的 に は!）である．また，上記の平衡規定方法では，回転変換や磁気TASK/WM［６４］，"）TORIC［６５］，#）AORSA［６６］等 がシアといった平衡量が必ずしも実際を再現している保証は用いられており，iv）のTASK/WF［６７］への適用が進めらないため，これらに大きく依存する MHD 安定性解析につれている．

いては，TASK3D-a の適用範囲を安定性解析に拡げていく

生成された高速粒子は，衝突による減速や空間拡散に上での大きな開発課題となっている．現状でも，VMEC

よって緩和する．その速度分布関数の緩和過程は，!）モン

計算の出力として，メルシエ安定性［７８］の評価指標は参照

テカルロ法による旋回中心の軌道追跡あるいは"）フォッ

できるが，それを用いた議論は行わないように喚起していカープランク型の偏微分方程式によって記述される．旋回る．

中心の軌道追跡においては，クーロン衝突や中性粒子とのこのように規定された VMEC 平衡データは，Boozer 座

衝突による減速，エネルギー拡散，ピッチ角散乱が乱数で標系［７９］にも変換される．TASK3D-a では，VMEC 平衡や記述され，速度空間２次元実空間３次元の分布関数が求め

Boozer 平衡のデータが，その後の加熱計算や新古典輸送計

られるとともに，背景粒子に与えられるパワー，運動量等算等に提供されていく流れとなっている．また，磁場スペが評価される．さらに損失粒子については，壁面に衝突す

クトルや実効ヘリシティ［８０］（リップル輸送の大きさの指４３４

Special Topic Article

2. Physics Modules Constituting Integrated Code

S. Murakami et al.

標）といった平衡量が，GIOTA コード［８１］を経て，LHD

生する高エネルギー粒子の軌道が複雑であり，加熱分布の解析データサーバ［８２］に登録され，磁場配位依存性の検討

解析には有限の粒子軌道効果を含むシミュレーションが必などの共同研究に広く活用されている．

要となる．このため，TASK3D では，２種類の NBI 解析
コードを開発している．１つは，簡易的に粒子軌道の効果

２．
１３ ３次元新古典輸送・径電場・加熱を含む FIT3D［５８］コードである．FIT3D は，また３つの３次元磁場配位であるヘリカル系装置のプラズマ輸送に

コ ー ド HFREYA，MCNBI，FIT で 構 成 さ れ て い る．

おいて新古典輸送およびそれに密接に関連する径電場が重HFREYA コードにより，イオン源のビーム粒子の発生か要な働きをする．統合シミュレーションコード TASK3D

らプラズマ中でのイオン化までをシミュレーションし，［８３，８４］においては，新古典輸送解析モジュール と して

MCNBI コードにより有限の軌道幅の効果を取り入れた

LHD/DGN［８５］モジュールがある．一般的な新古典輸送解

ビーム分布を計算している．次に，FIT コードによりビー

析には，各エネルギー成分に対する運動論方程式を解き，ムのエネルギー減衰，ピッチ角散乱をフォッカープランク

エネルギー積分を行うことによって輸送係数を得る．しか方程式の解析的近似解を用いて計算している．比較的短時しながら，ヘリカル系の新古典輸送においては，イオンお間に計算できることから多くのショットについて解析が可よび電子の輸送係数が径電場に強く依存するため，径電場能となる．また，最近では近似的に時間発展プラズマに対はイオンおよび電子流束の両極性条件から決定される．こ応した FIT̲CONV［９０］が開発されている．これについてのため，新古典輸送と径電場を同時に解く必要があり，多は，４．
９節において解説されている．

くの新古典輸送計算が必要となる．一方，ドリフト運動論

もう一方は，HFREYA コードの結果を用いて，GNET

方程式の磁力線方向の運動と "!!ドリフトのみを考慮す［６８］コードによりビーム減衰分布を粒子軌道の効果を完全ることにより，新古典輸送係数は，規格化することができ，に取り入れて解法する手法である．FIT3D の場合に比べ解規格化された拡散係数 ""#"!
"&%#'$#(は，規格化された衝析に長い CPU 時間が必要となり限定されたショットの解

突周波数 #"# %#，径電場 $"# %#& の関数として与えら
"
'
"
&
'
!

析のみ行うことができる．また，GNET コードでは，ICRF

れる．

ICRF加熱による加熱分布やNBI加熱とICRF加熱の重畳の加熱による RF 波との相互作用を取り入れることができ，"したがって，これら規格化された変数と " の関係を基場合において解析を行うことができる．さらに，ペレット

礎データとしてデータベース化することにより，異なる粒入射などプラズマの大きな時間変化を伴う場合の加熱解析

子種や温度・密度に対応した輸送係数を求めることができのため，プラズマ時間発展および複数のイオン粒子種に対る（図５）．

応 し た GNET-TD コ ー ド が 開 発 さ れ て い る．図６は，GNET-TD［９１］によって解析された炭素ペレット入射によ

LHD/DGN モジュールでは，規格化されたデータとしる加熱量の変化を示している．

て，DCOM［８６，
８７］および GSRAKE［８８］コードによる結果
を用いている．DCOM コードはモンテカルロ法を用いて新古典輸送係数を評価するコードであり，約２０００点の計算結

参考文献

果を用いている．また，GSRAKE は解析的近似解を用いて［１］F.L. Hinton and R.D. Hazeltine, Rev. Mod. Phys. 48, 239(1976). VII 節.
［２］P. Helander and D.J. Sigmar, Collisional Transport in Magnetized Plasmas (Cambridge University Press, Cambridge, 2002).
［３］M. Honda and A. Fukuyama, Comput. Phys. Commun.
208, 117 (2016).
［４］S.P. Hirshman and D.J. Sigmar, Nucl. Fusion 21, 1079
(1981).

新古典輸送係数を評価するコードであり，LHD/DGN においては特に無衝突領域における拡散係数を評価するために用いており，約２００点のデータを用いて，LHD/DGN が構築されている．得られた結果は DKES コードを用いたデータ
ベースとベンチマークを行っており，良い一致が得られている［８９］．
３次元磁場配位における加熱においては，加熱により発図６
図５

規格化拡散係数 D＊の DCOM+GSRAKE（NEW）と DCOM
（OLD）によるデータ点の比較［８
５］
．

４３５

GENT-TD による炭素ペレット入射放電の全加熱入力の時間発展（H プラズマと多種イオンを考慮した場合の比較）［９１］．

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019［５］M. Kikuchi and M. Azumi, Plasma Phys. Control. Fusion
37, 1215 (1995).
［６］M. Honda et al., Nucl. Fusion 53, 073050 (2013).
［７］M. Honda, Comput. Phys. Commun. 181, 1490 (2010). Appendix A.
［８］R.A. Hulse, Nucl. Technol. Fusion 3, 259 (1983).
［９］H.P. Summers, The ADAS User Manual, http://adas.
phys.strath.ac.uk［１０］W.A. Houlberg et al., Phys. Plasmas 4, 3230 (1997).
［１１］https://w3.pppl.gov/ntcc/NCLASS/readme̲nclass̲pt.
html［１２］A.G. Peers, Phys. Plasmas 7, 268 (2010).
［１３］F.L. Hinton and S.K. Wong, Phys. Fluids 28, 3082 (1985).
［１４］M. Romanelli and M. Ottaviani, Plasma Phys. Control.
Fusion 40, 1767 (1998).
［１５］T. Fülöp and P. Hender, Phys. Plasmas 6, 3066 (1999).
［１６］E.A. Belli et al., Plasma Phys. Control. Fusion 56, 124002(2014).
［１７］E.A. Belli and J. Candy, Plasma Phys. Control. Fusion 54,
015015 (2012).
［１８］K. Hoshino et al., Nucl. Fusion 51, 083027 (2011).
［１９］Y. Shimizu et al., Plasma Fusion Res. 10, 3403062 (2015).
［２０］T. Takeda and S. Tokuda, J. Comput. Phys. 93, 1 (1991).
［２１］S. Tokuda and T. Watanabe, Phys. Plasmas 6, 3012 (1999).
［２２］N. Aiba et al., Comput. Phys. Commun. 175, 269 (2006).
［２３］M. Azumi et al., Proc. 6th Int. Conf. on Plasma Phys. 1, 200(1984).
［２４］N. Aiba et al., Phys. Plasmas 18, 022503 (2011).
［２５］J. Shiraishi et al., Nucl. Fusion 54, 083008 (2014).
［２６］N. Hayashi et al., プラズマ･核融合学会誌 80, 605 (2004).
［２７］T. Ozeki et al., Fusion Sci. Tech. 50, 68 (2006).
［２８］武智 学 他：プラズマ・核融合学会誌 85, 147 (2009).
［２９］白石淳也：プラズマ・核融合学会誌 94, 183 (2018).
［３０］N. Aiba et al., Comput. Phys. Commun. 180, 1282 (2009).
［３１］B. Pégourié, Plasma Phys. Control. Fusion 49, R87 (2007).
［３２］A. Matsuyama et al., Nucl. Fusion 52, 123017 (2012).
［３３］F. Koechl et al., EFDA-JET-PR(12) 57.
［３４］P.B. Parks and R.J. Turnbull, Phys. Fluids 21, 1735 (1978).
［３５］S. Tokunaga et al., Fusion Eng. Des. 123, 620 (2017).
［３６］H. Grad and J. Hogan, Phys. Rev. Lett. 24, 1337 (1970).
［３７］R.R. Khayrutdinov and V.E. Lukash, J. Comput. Phys. 109,
193 (1993).
［３８］A. Matsuyama and S. Miyamoto, private commun.
［３９］S.C. Jardin et al., J. Comput. Phys. 66, 481 (1986).
［４０］L.E. Zakharov, et al., Phys. Plasmas 19, 055703 (2012).
［４１］K. Harafuji et al., J. Comput. Phys. 81, 169 (1989).
［４２］H. Natsume et al., Plasma Fusion Res. 9, 3403048 (2014).
［４３］N. Hayashi et al., Nucl. Fusion 49, 095015 (2009).
［４４］H. Kawashima et al., Plasma Fusion Res. 1, 031 (2006).
［４５］R. Schneider et al., Contrib. Plasma Phys. 46, 3 (2006).
［４６］畑山明聖：プラズマ・核融合学会誌 77, 420 (2001).
［４７］星野一生 他：プラズマ・核融合学会誌 86, 681 (2010).
［４８］藤間光徳 他：プラズマ・核融合学会誌 86, 685 (2010).
［４９］Y. Feng et al., Contrib. Plasma Phys. 44, 57 (2004).
［５０］S. Braginskii, Rev. Plasma Physics, 1, 205, Consultants Bureau, NewYork, 1965.
［５１］A. Fukano et al., J. Nucl. Mater. 363-365, 211 (2007).
［５２］H.D. Pacher et al., J. Nucl. Mater. 415, S492 (2011).

［５３］清水勝宏，滝塚知典：プラズマ・核融合学会誌 80, 183(2004).
［５４］N. Hayashi et al., J. Nucl. Mater. 363-365, 1044 (2007).
［５５］M. Yagi et al., Contrib. Plasma Phys. 52, 372 (2012).
［５６］K. Tani et al., J. Phys. Soc. Japan 50, 1726 (1981).
［５７］R.J. Goldston et al., J. Comp. Phys. 43, 61 (1981).
［５８］S. Murakami et al., Trans. Fusion Technol. 27, 256 (1995).
［５９］K. Hamamatsu and A. Fukuyama, Fusion Eng. Des. 53,
53 (2001).
［６０］A. Fukuyama, Fusion Eng. Des. 53, 71 (2001).
［６１］Y.R. Lin-Liu et al., Phys. Plasmas 10, 4064 (2003).
２］R.W. Harvey et al., Phys. Fluids B 5, 446 (1993).
［６［６３］N.B. Marushchenko et al., Plasma Fusion Res. 2, S1129(2007).
［６４］A. Fukuyama et al., Proc. 18th IAEA Fusion Energy Conf.
(Sorrento, 2000) IAEA-CN-77/THP2/26.
［６５］M. Brambilla, Plasma Phys. Control. Fusion 41, 1 (1999).
［６６］E.F. Jaeger et al., Phys. Plasmas 8, 1573 (2001).
［６７］A. Fukuyama and Y. Ichida, Proc. 1996 Int. Conf. on
Plasma Phys. 2, 1342 (1997).
［６８］S. Murakami et al., Nucl. Fusion 46, S425 (2006).
［６９］H. Nuga and A. Fukuyama, Plasma Fusion Res. 5, S2068(2010)［７０］M. Toma et al., Plasma Phys. Control. Fusion 57, 095007(2015).
［７１］R.W. Harvey et al., Phys. Rev. Lett. 62, 426 (1989).
［７２］M. Yokoyama et al., Nucl. Fusion 57, 126016 (2017).
［７３］K. Narihara et al., Rev. Sci. Instrum. 72, 1122 (2001).
［７４］C. Suzuki et al., Plasma Phys. Control. Fusion 55, 014016(2013).
［７５］S.P. Hirshman and J.C. Whiston, Phys. Fluids 26, 3553
(1983).
［７６］M. Yokoyama et al., Plasma Fusion Res. 8, 2403016 (2013).
［７７］Y. Suzuki et al., Nucl. Fusion 46, L19 (2006).
［７８］C. Mercier, Nucl. Fusion 1, 47 (1960).
［７９］A.H. Boozer, Phys. Fluids 25, 520 (1982).
［８０］C.D. Beidler and W.N.G. Hitchon, Plasma Phys. Control.
Fusion 36, 317 (1994).
［８１］M. Yokoyama et al., Research Report NIFS-810, National
Institute for Fusion Science, Japan (2005) on the numerical code originally developed by L. Hedrick (retired, Oak
Ridge National Laboratory).
［８２］M. Emoto et al., Fusion Eng. Des. 81, 2019 (2006).
［８３］A. Sakai et al., Plasma Fusion Res. 10, 3403048 (2015).
［８４］S. Murakami et al., Plasma Phys. Control. Fusion 57,
054009 (2015).
［８５］A. Wakasa et al., Contrib. Plasma Phys. 50, 582 (2010).
［８６］A. Wakasa et al., J. Plasma Fusion Res. SERIES 4 (2001)
408.
［８７］A. Wakasa et al., Jpn. J. Appl. Phys. 46, 1157 (2007).
［８８］C.D. Beidler and W. D. D'haeseleer, Plasma Phys. Control.
Fisoin 37, 463 (1995).
［８９］Y. Turkin et al., Phy. Plasmas 18, 022505 (2011).
［９０］H. Nuga et al., Plasma Fusion Res., accepted for publication(2019).
［９１］H. Yamaguchi and S. Murakami, Nucl. Fusion 56, 026003(2016).

４３６

J. Plasma Fusion Res. Vol.95, No.9 (2019)４
３７‐４４３

小特集 統合コードによる磁場閉じ込め核融合プラズマシミュレーションの現状と今後の展望

３．国内外の統合コード開発状況
3. Status of Integrated Code Development in Japan and Overseas福山淳，林伸 彦１），村 上 定 義，横 山 雅 之２，４），藤 田 隆 明３）

FUKUYAMA Atsushi, HAYASHI Nobuhiko1）, MURAKAMI Sadayoshi,
YOKOYAMA Masayuki2，４）and FUJITA Takaaki3）
京都大学大学院工学研究科，１）量子科学技術研究開発機構，２）自然科学研究機構核融合科学研究所，
３）

名古屋大学大学院工学研究科，４）統合研究大学院大学（原稿受付：２
０
１
９年５月２
０日）

磁気閉じ込め核融合プラズマの包括的シミュレーションに向けて，世界各国で統合コードの開発が進められるとともに，多くのモジュールを連携させるためのコードインターフェースの標準化を含む統合コードフレーム
ワークの試みも進展しつつある．本章では，まず国内で開発が進められている TOPICS（量研），TASK（京大），
TOTAL（名大），TASK3D（核融合研，京大）について，その構成と主な特徴を説明する．引き続いて，欧米の統合コード開発の現状を紹介するとともに，ITER プラズマ実験に向けて，コード間連携とデータの標準化をめざして開発・普及が進められている IMAS を説明し，日米におけるコード間連携の試みについても報告する．
Keywords:
integrated code, comprehensive simulation, integrated modeling framework, data structure, data interface

３．
１ はじめに準化の試みとともにその概要を説明する．

磁気核融合炉心プラズマにおいては，様々な物理現象が３．
２ 統合コード TOPICS

相互に影響を及ぼし合いつつ時間発展している．その振る舞いを記述するためには，第２章で紹介された個々の物理

TOPICS（TOkamak Prediction and Interpretation Code

現象を記述する解析コード（モジュール）を組み合わせて，System）は，元々は日本原子力研究所で安積正史氏が開発

炉心プラズマ全体の時間発展を自己無撞着にシミュレート

した 1.5 次元輸送コードであり，統合コードという概念がする必要がある．そのために多数の解析モジュールを内包

ない時代でも，内包する物理モジュールは小規模ながら作する多くの統合コードが開発されてきた．主要な磁気閉じりそのものは既に統合コード的なものであった．TOPICS

込め装置においてはその運転性能を予測し，実験結果を解は，JT-60 実験のプラズマ分布等を評価する実験解析ツー

析するための統合コードが開発され，新しい実験装置の性

ルと連携するインターフェースが整備され［１，２］，組織名能予測や解析モジュールの機能向上を図る上でも統合コーが量研に変わった現在に至るまで JT-60 の標準の輸送解析

ドが利用されている．国内で開発されている代表的な統合

コードとして使われてきた．プラズマの時間発展を予測す

コードには TOPICS（量子科学技術研究開発機構（量研）），るシミュレーションコードとしても開発が進められ，統合

TASK（京 都 大 学（京 大）），TOTAL（名 古 屋 大 学），

コ ー ド の 概 念 が 現 れ た 時 に，TOPICS-IB（TOPICS ex-

TASK3D（核融合科学研究所，京大）があり，以下ではそ

tended for Integrated simulations for Burning plasma）に呼れらの特徴を紹介するとともに，主要な機能を比較し，国称を一時変更していたが［３，４］，統合コードが一般的にな外の統合コードについても簡単に説明する．

るとともに呼称を元の TOPICS に戻した
［５］．ただし，当

統合コードとして多数の解析モジュールの連携を系統的

初 TOPICS という名称であった 1.5 次元輸送ソルバー部分に実現するためには，連携のフレームワークを構築するこは，図１に示すように統合コード TOPICS のコア部分になとが重要である．そのフレームワークには，モジュール間り，他の研究機関の研究者とも協力して様々な規模の物でやりとりするデータのデータ構造，データ交換インター

理・工学モジュールが導入・結合されている．コア部分

フェース，そしてモジュールの実行を制御する仕組み等がは，１次元の熱／粒子／運動量／電流輸送と２次元の自由含まれる．現在国際熱核融合実験炉 ITER プラズマ実験の

境界 MHD 平衡を，プラズマ電流立ち上げ／下げ時の体積モデリングと実験解析に向けた統合フレームワーク IMAS

変化も考慮して解き［６］，新古典理論と整合したトロイダ（Integrated Modelling and Analysis Suite）の開発が欧州を

ル・ポロイダル回転と径電場を評価できる［７］．

中心に進められており，日米におけるモジュール間連携標

TOPICS の物理モジュールには，燃料ペレットの入射とKyoto University, KYOTO 615-8540, Japancorresponding author’s e-mail: fukuyama@afportal.net

４３７!2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

し，終了して出力ファイルが現れたら TOPICS 本体に転送する．この機能は，グリッドミドルウェア AEGIS（Atomic
Energy Grid Infrastructure）に基づいた統合計算フレーム
ワ ー ク SOAF（Simple Orchestration Application Framework）［１
６］で も 実 装 さ れ TOPICS で 利 用 さ れ た．一 方，SONIC ではファイルによるデータ交換でなく並列計算のMPI（Message Passing Interface）に 基 づ い た MPMD
（Multiple-Program Multiple-Data）システムの統合計算機能を実装しており［１７］，両者を融合した統合計算機能の開発を進めている．

図１

３．
３ 統合コード TASK

量研で開発されている統合コード TOPICS の構成図．

統 合 コ ー ド TASK（Transport Analyzing System for
溶発および溶発プラズマの大半径外側へのドリフトを考慮

tokamaK）［１８］の出発点は，当時 JT-60 で実験が始まろうとした APLEX［８］，新 古 典 トロイダル 粘 性 によるトルクを 評していた低域混成波加熱を解析するために１９７７年に開発さ価する第一原理３次元新古典コード FORTEC-3D
［９，
１０］，線

れた光線追跡法による波動伝播解析コード［１９］である．そ形 MHD 安定性コード（回転なし：MARG2D，回転あり：の後，イオンサイクロトロン波を記述する波動光学的伝播

MINERVA），RWM 安 定 性 解 析 コ ー ド MINERVA/

解析［２０］，高速イオンの生成と電流駆動を記述する速度分

RWMaC，新古典輸送係数やブートストラップ電流を評価

布関数解析［２１］，プラズマの時間発展を記述する輸送・平する Matrix Inversion モデル［１，１１］，NBI（Neutral Beam

衡解析が加わり，波動による加熱と電流駆動の自己無撞着

Injection）や核融合反応で生成される高エネルギー粒子のなシミュレーションを実現する統合コードへと発展

加熱・電流駆動やトルクを評価するモンテカルロ粒子コー

し，１９９２年頃から TASK と呼ばれるようになった．

ド OFMC と２次元フォッカープランクコード BAFP，電磁統合コード TASK の特徴の一つは，当初からモジュール

波 に よ る 加 熱・電 流 駆 動 を 評 価 す る コ ー ド（EC: EC-

構成で開発されてきたことである，現在のモジュール構成

Hamamatsu，LH: Bonoli コ ー ド［１２］，IC: TASK/WM），

を図２に示す．各モジュールは単体でも動作可能である

炉心の不純物１次元輸送コードIMPACT［１３］，SOL／ダイ

が，他のモジュール内からの実行も可能であり，全体統括バータプラズマの統合コード SONIC と動的５点モデル D

モジュ ー ル TOT か ら 時 間 シ ー ク エ ン ス に 従 っ て 各 モ

5PM，等が結合されている．SONIC はそれ自体が統合コー

ジュールを実行させることも可能である．各種物理量の空ドであるが，別の小特集［１４］で既に紹介されているので詳間分布はデータ交換モジュール PL に保存されてデータ交

細はそちらをご覧いただきたい．上記の内のいくつかのモ換を行うが，後述のデータ交換フレームワーク BPSD を利

ジュールは第２章で（上記で参考文献がないものは第２章

用することも可能である．

を参照），それらを利用した研究は第４章で紹介されてい平衡，輸送，加熱等の基本的な解析モジュールに加えて，る．一方，TOPICS の工学モジュールには，加熱・電流駆先進的なモデリングの試みが行われている．輸送について動による規格化ベータ値や周回電圧等の複数パラメータのは，従来の拡散型輸送方程式による解析（TR，TI）に加え統合 PID 制御モジュール［５］等があり，JT-60SA で使用すて，磁気面平均された流体方程式とマクスウェル方程式をる MHD 平衡制御シミュレータ MECS
［１５］とも結合する予解くTASK/TX［２２］はより短い時間スケールの輸送現象を定である．
TOPICS の計算は，使用するモジュールが小規模であれ
ば単一の PC クラスタ計算機で行われ，モジュールが大規模並列計算を必要とする場合はそのモジュール部分だけ
スーパーコンピュータ（スパコン）で実行 し，TOPICS
本体は PC クラスタで実行してモジュールの実行をリモー
トで自動制御することが可能である．もちろん全ての計算をスパコン上で実行することも可能であるが，現状のシ
ミュレーションでは数日以上かかる場合もあり，スパコンのジョブクラスの時間制限によってリスタート作業が繁雑になったり，スパ コ ン の 障 害 に 直 接 影 響 さ れ て しまう（TOPICS 本体の計算を PC クラスタで継続していれば，他の正常なスパコンでの計算に切り替えることができる）．
モジュールのリモート制御のためスパコン上にインストー
ルした監視プログラムは，TOPICS 本体から入力ファイル

図２が転送されるとそのファイルを認識してモジュールを実行
４３８

統合コード TASK の構成図．

Special Topic Article

3. Status of Integrated Code Development in Japan and Overseas

A. Fukuyama et al.

記述することができる．さらに３次元速度分布関数の時間

ラズマにおける新古典輸送の軸対称成分は，文献
［３１，
３２］

発展を解析する TASK/FP［２３］は，波動加熱や NBI 加熱にの式に基づいて評価される．ヘリカルプラズマにおける新よる速度分布関数の変形を含めた運動論的輸送解析が可能

古典輸送の非軸対称成分は GIOTIA コード［３３］によりマルである（５．
５節参照）．統合コード TASK はその起源から波チ成分のヘリカルリップルの効果を含めて求められ，イオ

動伝播解析は充実しているが，MHD 安定性解析は手薄でンと電子の粒子束が等しくなる条件から径電場が決定されある．しかしながら，ICRF 波加熱のために開発された

る．異 常 輸 送 と し て は，GLF23［３４］，Bohm
［３
５］
，Gy-

TASK/WM を複素周波数に拡張して固有モードを計算す

roBohm［２７］，Mixed Bohm/GyroBohm［３６］，CDBM などることにより，高速イオンによって励起されるアルヴェンのモデルが導入されている．"!! シアによって異常輸送

固有モードの成長率を計算することができる［２４］．

を低減させる式により内部輸送障壁の形成をモデル化してTASK のその他の特徴としては．ソースコードが公開さいる［２７，
３６］．H モードについては，プラズマ境界部で人為れていて git サーバーから入手可能であること［２５］，計算的に異常輸送を低減させて境界輸送障壁を模擬している．

資源を要するモジュールを MPI により並列処理できるこ

不純物輸送は NCLASS［３７，３８］による新古典輸送に基づいと，同一コードでラップトップパソコンからクラスター，て IMPDYN［３９］で解析される．
MHD 不 安 定 性 と し て は，鋸 歯 状 振 動（sawtooth）

スーパーコンピュータまで処理が可能であること等があげ［３５］，新古典テアリングモード（NTM）［４０］，メルシエ条

られる．

件，短波長バルーニングモード［３０］を評価するモジュール

今後は運動論的波動伝播解析の拡張，運動論的輸送解析が組みこまれている．

の整備，周辺プラズマ輸送との結合等に加えて，ITER 物

ペレットについて中性ガス遮蔽モデル（NGS モデル）と

理データモデルとのインターフェースの開発を進める予定［４１］，ガスパ"!! ドリフトによるリロケーションモデル

である．

フおよびリサイクリング粒子についてモンテカルロ法によ

３．
４ 統合コード TOTAL

る AURORA コード［４２］が組み込まれている．SOL および

TOTAL（TOroidal Transport Analysis Linkage）コード

ダイバータプラズマに関しては２点モデル［４３］が用いられは，ヘリカルプラズマ用のコード（HSTR）［２６］をベースとている．中性粒子ビーム（NB）入射についてフォッカープ

して，ヘリカル，トカマクのいずれにも対応できるコード

ランクコード FIFPC［４４］での解析，電子サイクロトロン波として開発が進められてきた．これまでに，LHD プラズマ

入射について光線追跡による解析が可能である．加熱分の性能予測［２６］，LHD プラズマにおける電子熱輸送障壁

布，駆動電流分布を関数形で与えることもある［４５］．

の解析［２７］，ITER におけるタングステン不純物輸送の予核燃焼プラズマの予測計算においては，核融合出力とし測［２８］などに用いられた．

てはサーマル成分によるもののみが考慮され，プラズマパ

TOTAL コードの構成をフローチャート形式で図３に示

ラメータに応じてアルファ粒子のイオン加熱パワー・電子

す．ヘリカルプラズマの解析においてはVMECコード［２９］

加熱パワーが評価される．核融合出力の制御として，ペにより３次元平衡を計算し，トカマクプラズマの解析におレット入射周波数あるいは外部加熱パワーによる帰還制御

いては APOLLO コード［３０］で自由境界２次元平衡を計算が用いられる．

する．得られた平衡における１次元輸送が HTRANS コー

近年，プラズマのトロイダル回転に伴う高価数不純物のドにより解析される．

ピンチモデル（PHZ，Er ピンチ）
［４６］の式が導入され実験

トカマクプラズマにおける新古典輸送およびヘリカルプ

データの解析が行われた（第４章に記述）．それらのピンチ速度を案内中心軌道の解析により求めるコードも最近組
み込まれている．H モードについて，LH 遷移・HL 逆遷移の導入や EPED1 比例則に基づくペデスタル圧力限界の導入も行われて い る．さ ら に NB 電 流 駆 動 解 析 の た め のフォッカープランクコードの作成・導入も進行中である．
TOTAL コードの特徴の一つとして，比較的簡易なモデルを用いて高速に計算を行える点があり，それを活かした原型炉等におけるプラズマの実時間制御への適用などが開発の方向として考えられている．

３．
５ 統合コード TASK3D
LHD をはじめとしたヘリカル系プラズマを主な対象とした統合輸送解析スイート TASK3D の開発・運用が進んでいる．トカマクプラズマを対象とした統合解析コード
TASK［２５］を基盤として，３次元磁場構造に起因・関連する物理を記述するモジュールの追加，改良などが進められ
図３

TOTAL コードの構成（フローチャート）．

ている．TASK3D の開発戦略は以下の２通りである．LHD
４３９

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

実験データに基づいた【LHD 実験データ解析型（TASK3よるNBI加熱計算への必要データ提供を挙げることができ

D-Analysis，TASK3D-a シリーズ）】は，輸送現象に関するる（４．
１０節参照）．これにより，大規模シミュレーションに

LHD 実験解析の基盤となっており，各種モジュールの実験

よる LHD プラズマ解析の実施や国際コードベンチマーク

検証や高精度化などが行われている．一方，【到達パラ

活動［５２］が円滑にできるようになったことも TASK3D-a

メータ予測型】は，輸送モデルを導入して到達パラメータの成果であると言える．

（温度など）を予測するものである．以下の節では，それぞ複数イオン種から成る LHD プラズマや，LHD 重水素実れについて記述する．

験への対応をさらに進めて，軽水素，重水素，ヘリウム存

LHD 実験データ解析型統合輸送解析コード：TASK3D-a

在下での NBI 加熱計算
［５３］に対応した第３版 a03，第４版

３次元平衡，NBI 加熱，熱輸送過程といった各種要素解

a04 へと次々と機能拡張が行われた［５４］．この機能拡張に

析コード（モジュール）の統合作業を経て，２０１２年９月に，より，同位体効果研究の基盤となる解析データベースの作

初 版 で あ る TASK3D-a01 が リ リ ー ス さ れ た［４７］．

成が大きく進展している［５５，５６］．また，LHD 重水素実験

TASK3D-a01 の構築によって，熱輸送解析の自動化が進展における中性子計測，高エネルギー粒子挙動の定量評価のした．温度分布や密度分布の時間変化も考慮したダイナ基盤としての役割も果たしている（４．
９節参照）．

ミック輸送解析の結果も提供することができるようになっ

このように TASK3D-a は，要素モジュールの追加や高精

た．

度化が次々と進められ，多様な研究課題や系統的検証研究

よ り 広 範 な LHD 実 験 条 件 へ の 適 用 を 図 る た め，モのための統合解析基盤として運用が行われている．今後のジュールの追加や高精度化を行い，TASK3D-a02 としてリ

課題は枚挙に暇がないが，喫緊のものとしては，ICH モ

リースした．新古典熱・粒子拡散フラックスを評価する

ジュールの導入，モジュールとしての導入が可能な不純物

GSRAKE コード［４８］を組み込んだ．これにより，a01 での輸送コードの検討と導入，周辺プラズマ解析との連携など

実験的なエネルギーバランス解析とともに，新古典エネルが挙げられる．

ギーフラックス解析が同時に行われることになった．まITER 統合モデリング専門家会合において，大型実験にた，ECH 加熱分布を評価する LHDGauss コード［４９］および

おいて実験データと密に連携して運用されている統合解析

TRAVISコード［５０］も導入したことで，ECH吸収パワーをスイートとして TASK3D-a が認知されていることを最後取り込んだエネルギーバランス解析が行えるようになり，に述べておきたい

a01 段階での NBI プラズマ解析から解析対象の格段の拡張

到達パラメータ予測型統合輸送解析コード：TASK3Dが実現した．

トカマク輸送解析を対象とする TASK/TR を ベ ー ス

さらに，計算規模や必要資源の違いから統合が困難であに，３次元 MHD 平衡，新古典輸送，径電場効果，３次元磁る大規模シミュレーションコードについても，TASK3D-a

場加熱解析等を取り入れてヘリカル系輸送解析に拡張されから LHD プラズマの平衡や温度・密度分布などを提供す

た統合コードであり，TAS3D-aと区別するためTASK3D-p

るという連携を確立した．その一例として，GNET［５１］にと呼ばれることもある［５７，５８］．主に，加熱によって生成された高速粒子の振る舞いを記述する分布関数解析コードと組み合わせて，LHD，HSX
GNET［５９］，GNET-TD［６０］
等のヘリカル系プラズマの時間発展の予測に用いられている．また，LHD 実験結果との比較により熱・粒子輸送モデルの構築やヘリカル型核融合炉の検証が進められている．

３．
６ 統合コードの機能
各統合コードに含まれる解析モジュールの主な機能の一
覧を表１に示す．TOPICS はトカマクプラズマにおける物理現象を幅広くカバーしており，TASK は波動加熱・電流駆動の解析に，TOTAL はヘリカル系を含めた輸送解析にウェイトがある．TASK3D（TASK3D-a，TASK3D-p）もヘリカル系プラズマの主要な物理現象をカバーしており，
TASK3D-a は実験解析との連携に特徴がある．

３．
７ 海外の統合コード
図４

LHD 実験データ解析型（TASK3D-a シリーズ）の機能拡張の進展．中央部横一列で示された a01 版から，下段に記された機能拡張（a02 版），NBI 加熱モジュールを中心とした
複数イオン種（軽水素，ヘリウム，重水素）対応（a03，a04），さらに，複数イオン種の新古典輸送解析に関する大規模シミュレーションコードとの連携や高エネルギー粒子
挙動に関する機能拡張（上段部）などの進展を遂げている．

欧米においても，多数の統合コードが開発されている．
JINTRAC（JET Integrated Transport Code）［６１］は JETで開発された輸送コード JETTO［６２］をベースにしており，
最近は Core-SOL 結合を含めた解析を行っている．ASTRA（Automated System for TRansport Analysis）［６３］は IPPGarching で開発された輸送コードをベースにした統合コー
４４０Special Topic Article

表１

3. Status of Integrated Code Development in Japan and Overseas

各統合コードの機能一覧．

MHD 平衡

固定境界平衡自由境界平衡
平衡データ EFIT 形式
読込
MECS 形式
VMEC 形式
MHD 安定性 鋸歯状振動
NTM
RWM
ELM
AE
Mercier/n＝∞BM
輸送
バルクイオン
高速イオン
不純物イオン
中性粒子
トロイダル回転径方向電界
新古典輸送（2D）
新古典輸送（3D）
乱流輸送
高エネルギー 軌道追跡
粒子
FP 方程式NBI 加熱高速イオン生成
波動加熱
光線追跡法
ビーム追跡法
波動光学的解析
ペレット入射NGS モデル
ドリフト効果ダイバータ 簡易モデル
ダイバータコード

利用されている．

TOPICS TASK TOTAL TASK3D
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
△
○
○
○
△
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○
○

○
○
○
○
○
○

A. Fukuyama et al.

○
○
○

３．
８ 統合モデリング解析フレームワーク
統合コードは多数の解析モジュールから構成され，それらの間でデータをやり取りする．従来の多くの統合コード
は，特定の平衡・輸送コードをベースとし，それが MHDや加熱等の個々の解析モジュールとデータをやり取りする形で構成されていた．そのため，ある解析モジュールを別の統合コードに移植するたびに，新しいデータインター
フェースを開発する必要があった．また，統合コードもその利用目的に応じて多様な解析レベルがあり，それに合わせて適切な解析レベルと必要計算負荷の解析モジュールを使用する必要がある．これらの多様性に対応するために，共通の統合化フレームワークを構築する試みが行われてきた．
統合コードを構成するためのフレームワークは，１）解析モジュール間でやり取りするデータの種類とその構
造，２）データをやり取りするためのプログラムインター
フェース（API），３）解析モジュールの実行を制御する
ワークフロー管理等から構成される．

○
○
○
○

国内では，データ構造とデータ交換 API を定義する
BPSD（Burning Plasma Simulation Data）［７４］が開発され，
統合コード TASK で利用されている．その開発方針は
１）データ量を抑えた必要十分なデータセットであること，２）並列処理に対応したデータ構造であること，３）単

○
○
○

○
○

純かつ適応性の高いデータ交換 API であることである．
ITER 機構では，新しい統合コードを開発するのではな
く，既存の解析モジュールを組み合わせる統合コードパッケ ー ジ IMAS（Integrated Modelling and Analysis Suite）［７５］を構築している．IMAS に基づいて開発される統合ドである．CRONOS［６４］はCEA-Saint-Paul-lez-Duranceで開コードにはさまざまなレベルがあり，その対象は即時のプ発された統合コードであるが，粒子・熱輸送モデリングを

ラズマ再構成，制御システム開発のためのプラズマシミュ

簡略化して高速化を図った METIS［６５］も開発されている．

レータ，運転シナリオ開発等の運転支援から，実験結果と一方，統合モデリングのフレームワークを開発するための比較による物理モデルの検証や物理機構の解明等の研究に，EFDA ITM-TF（Integrated tokamak modeling - task

支援まで幅広く，解析モジュール間のデータのやりとりだ

force）が 設 置 さ れ，European ITM［６６］が 開 発 さ れ た．

けではなく，実験データの取り込みや蓄積を視野に入れてEuropean ITM は，モジュール間でやり取りするデータをいる．そのため，IMASのデータモデルは，密度や温度など定義する CPO（Consistent Physical Object）とそのデータの物理量だけでなく，実験装置や計測装置の装置パラメーをやり取りすプログラムインターフ ェ ー ス を 定 義するタ，生データに近い計測データ等を記述する５０種類以上の

UAL（Universal Access Layer）から構成され，モジュール

IDS（Interface Data Structure）を定義し，さらに増えつつの連携を定義するワークフローマネージャーとして Ke-

ある．

pler
［６７］が採用された．そして，ITM-TF の活動として，既IMAS のデータモデルや API は EFDA ITM-TF が開発し存の輸送コードとは独立に新しい輸送コード ETS（Euro-た CPO と UAL をベースとして開発され，改良が加えられ

pean Transport Solver）［６８］も開発されている．

ている．そのため，EU が開発しているコードの多くに

米国 PPPL において開発された輸送解析用統合コード

IMAS API が 実 装 さ れ，相 互 接 続 の 検 証 が 行 わ れ て い

TRANSP［６
９，７０］は世界的に広く利用されているが，時間

る．また，ITER 物理研究活動（ITPA: International Toka-

発展を予測するための統合コード PTRANSP［７０］にも拡張

mak Physics Activity）におけるコード間ベンチマークテされている．GA では輸送コード ONETWO［７１］をベース

スト等で，シミュレーションデータや実験データの交換にとして統合シミュレーションが行われていたが，最近は後IMAS のデータ形式 IDS が用いられることも多くなってき

述の統合化フレームワーク OMFIT［７２］に基づいて統合ている．また，EU における実験装置のデータを IMAS 形式

コードが開発されている．また，LLNL で開発された統合に変換するプロジェクトも進行している．ただ，データモ

コード CORSICA［７３］は，ITER の制御プログラム開発にデル仕様がかなり大きく，一つの IDS の全てのデータを設
４４１

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

EPS2015PAP/pdf/P5.145.pdf ).
［６］M. Honda, Comput. Phys. Commun. 92, 181 (2010).
［７］M. Honda et al., Nucl. Fusion 53, 073050 (2013).
［８］N. Hayashi et al., Nucl. Fusion 53, 123009 (2013).
［９］S. Satake et al., Plasma Fusion Res. 3, S1062 (2008).
［１０］M. Honda et al., Nucl. Fusion 55, 073033 (2015).
［１１］M. Honda, Phys. Plasmas 21, 092508 (2014).
［１２］P.T. Bonoli and R.C. Englade, Phys. Fluids 29, 2937 (1986).
［１３］N. Hayashi et al., Nucl. Fusion 58, 066001 (2018).
［１４］K. Hoshino, J. Plasma Fusion Res. 92, 882 (2016).
［１５］Y. Miyata et al., Plasma Fusion Res. 9, 3403045 (2014).
［１６］T. Tatekawa et al., 3rd Int. Joint Conf. on Comput. Sci. Opt.,
IEEE Computer Society, DOI 10.1109/CSO.2010.162(2010).
［１７］K. Hoshino et al., Plasma Fusion Res. 9, 3403070 (2014).
［１８］A. Fukuyama et al., Proc. 20th IAEA FEC (Villamoura,
Portugal, 2004) IAEA-CSP-25/CD/TH/P2-3.
［１９］百田 弘 他：JAERI-M 6964 (1977).
［２０］A. Fukuyama et al., Comp. Phys. Rep, 4, 137 (1986).
［２１］T. Morishita et al., J. Phys. Soc. Jpn. 57, 1238 (1988).
［２２］M. Honda and A. Fukuyama, J. Comp. Phys. 227, 2808
(2008); M. Honda and A. Fukuyama, Comp. Phys. Comm.
208, 117 (2016).
［２３］H. Nuga and A. Fukuyama, Prog. Nucl. Sci. and Tech. 2,
78 (2011).
［２４］A. Fukuyama and T. Akutsu, Proc. 19th IAEA FEC (Lyon,
France, 2002) IAEA-CN-94/TH/P3-14.
［２５］TASK, https://bpsi.nucleng.kyoto-u.ac.jp/task/［２６］K. Yamazaki and T. Amano, Nucl. Fusion 32, 633 (1992).
［２７］J. Garcia et al., Phys. Rev. Lett. 96, 105007 (2006).
［２８］Y. Murakami et al., J. Nucl. Mater. 313-316, 1161 (2003).
［２９］S.P. Hirshman et al., Comput. Phys. Commun. 43, 143(1986).
［３０］K. Yamazaki et al., Nucl. Fusion 25, 1543 (1985).
［３１］F.L. Hinton and R.D. Hazeltine, Rev. Mod. Phys. 48, 239(1976).
［３２］C.S. Chang and J.D. Callen, Phys. Fluids 26, 3315 (1982).
［３３］M. Yokoyama et al., NIFS-810 (2005).
［３４］R.E. Waltz et al., Phys. Plasmas 4, 2482 (1997).
［３５］H. Natsume et al., Plasma Fusion Res. 9, 3403048 (2014).
［３６］Y. Hori et al., Plasma Fusion Res. 5, S2034 (2010).
［３７］W.A. Houlberg et al., Phys. Plasmas 4, 3230 (1997).
［３８］https://w3.pppl.gov/ntcc/NCLASS/readme̲nclass̲pt.
html［３９］T. Amano et al., Internal report IPPJ-616, Institute of
Plasma Physics, Nagoya Univ. (1982).
［４０］S. Taniguchi et al., Plasma Fusion Res. 5, S2035 (2010).
［４１］T. Yamakami et al., Plasma Fusion Res. 8, 2403079 (2013).
［４２］M.H. Huges and D.E. Post, J. Comput. Phys. 28, 43 (1978).
［４３］K. Borrass et al ., Nucl. Fusion 31, 1035 (1991).
［４４］R.H. Fowler et al., Comput. Phys. Commun. 13, 323 (1978).
［４５］H. Natsume et al., Plasma Fusion Res. 10, 3403051 (2015).
６］K. Hoshino et al., Nucl. Fusion 51, 083027 (2011).
［４［４７］M. Yokoyama et al., Plasma Fusion Res. Special Issue 7,2403011 (2012).
［４８］C.D. Beidler and W. D. D'haeseleer, Plasma Phys. Contr.
Fusion 37, 463 (1995).
［４９］T. Ii. Tsujimura et al., Nucl. Fusion 55, 123019 (2015).
［５０］N. Marushchenko et al., Plasma Fusion Res. 2, S1129

定するには手間がかかるため，シミュレーションにおける頻繁なデータ交換には適さない面がある．部分的なデータ設定も許容されているが，必須項目が定義されていないた
め，その場合には相手モジュールのデータ設定・利用仕様を十分に理解する必要がある．
IMAS におけるワークフローマネージャーには，欧州と同じく Kepler［６７］が採用されている．Fortran から Kepler
へのインターフェースも開発され，そのインターフェース
を通してのデータ交換を，GUI（Graphic User Interface）を用いて設定することができる．EU を中心に広く使われるようになりつつあるが，今後の課題としては，並列処理への対応，多様な計算環境（バッチ処理，遠隔処理，多数ユーザー等）への適応性向上等が挙げられている．
米国では，GA を中心に統合コードフレームワーク OMFIT［７２］が開発されている．OMFITではPythonを使用したファイルベースのデータ交換を行っており，IMAS をそのまま OMFIT に組み入れるには多くの問題点があるた
め，IDS のデータ構造に基づく Python ライブラリ OMASを開 発 し，性 能 向 上 な ら び に 機 能 向 上（デ ー タ 一 貫 性チェック，エラーハンドリング，動的ツリー生成，可変データ切り出し，データベースの多様化）を実現している．
国 内 で は BPSD
［７４］と IMAS の 間 の デ ー タ イ ン タ ー
フェースを開発し，モジュール間のデータ交換は BPSDで行い，必要に応じて IMAS IDS を通してデータを交換・蓄積する方向で作業が進められている．

３．
９ おわりに統合コード開発は欧米と日本で進められており，ITER
運転に向けた国際連携も進みつつある．国内では４つの統合コードが開発されており，今後 ITER の運転や DEMOの設計に向けて，それぞれのコードの機能向上とともに，限られた人的資源を有効に利用するため，コード間連携を系統的に進める必要がある．また開発に必要な人材の育成に向けて，大学レベルでの解析モジュール開発をどのように進めるかが課題である．一方 ITER 機構で開発されている IMAS への国内の対応はやはり人的資源不足のためにかなり遅れている．国内で開発した解析モジュールの IMASへの提供や JT-60SA 実験環境における IMAS の利用をどのように進めるかは今後の課題である．また，周辺プラズマを含めた統合シミュレーションは一層重要な課題となりつつあり，高い信頼性をもつ周辺‐統合コードの開発が望ま
れる．

参考文献
［１］M. Kikuchi and M. Azumi, Plasma Phys. Control. Fusion 37, 1215 (1995).
［２］H. Shirai et al., Plasma Phys. Control. Fusion 42, 1193
(2000).
［３］N. Hayashi and JT-60 Team, Phys. Plasmas 17, 056112(2010).
［４］N. Hayashi et al., Plasma Fusion Res. 6, 2403065 (2011).
［５］N. Hayashi et al., 42nd EPS Conf. on Plasma Physics (Lisbon,
Portugal, 22-26 June 2015) P5.145 ( http://ocs.ciemat.es/

４４２

Special Topic Article

3. Status of Integrated Code Development in Japan and Overseas(2007).
［５１］H. Yamaguchi et al., Nucl. Fusion 56, 026003 (2016).
［５２］佐藤雅彦 他：プラズマ・核融合学会誌 93, 67 (2017).
［５３］P. Vincenzi et al., Plasma Phys. Control. Fusion 58, 125008(2016).
［５４］M. Yokoyama et al., Nucl. Fusion 57, 126016 (2017).
［５５］H. Yamada et al., 27th IAEA Fusion Energy Conference,
EX/P3-5 (Ahmedabad, Oct. 2018).
［５６］K. Tanaka et al., submitted to Nucl. Fusion (2019).
［５７］A. Sakai et al., Plasma Fus. Res. 10, 3403048 (2015).
［５８］S. Murakami et al., Plasma Phys. Control. Fusion 57,
054009 (2015).
［５９］S. Murakami et al., Nucl. Fusion 46, S425 (2006)［６０］H. Yamaguchi et al., Nucl. Fusion 56, 026003 (2016).
［６１］M. Romanelli et al., Plasma Fus. Res. 9, 3403023 (2014).
［６２］G. Cenacchi and A. Taroni, JET-IR (88) 03 (1988).
［６３］G.V. Pereverzev, P.N. Yushmanov, IPP-Report, IPP 5/98

A. Fukuyama et al.

(2002).
［６４］J.F. Artaud et al., Nucl. Fusion 50, 043001 (2010).
［６５］J.F. Artaud et al., Nucl. Fusion 58, 105001 (2018).
［６６］G.L. Falchetto et al., Nucl. Fusion 54, 043018 (2014).
［６７］https://kepler-project.org［６８］D. Coster et al., IEEE Trans. Plasma Sci. 38, 2085 (2010).
［６９］R.J. Hawryluk, Physics of Plasmas Close to Thermonuclear Conditions, ed. by B. Coppi, et al., (CEC, Brussels)
1, 19 (1980).
［７０］https://transp.pppl.gov［７１］W. Pfeiffer et al., GA-Report GA-A-16178 (1980).
［７２］O. Meneghini et al., Nucl. Fusion 55, 083008 (2015).
［７３］J.A. Crotinger et al., LLNL Report UCRL-ID-126284;
NTIS #PB2005-102154 (1997).
［７４］福山 淳：プラズマ・核融合学会誌 85, 597 (2009).
［７５］F. Imbeaux et al., Nucl. Fusion 55, 123006 (2015)４４３

J. Plasma Fusion Res. Vol.95, No.9 (2019)４
４４‐４５２

小特集 統合コードによる磁場閉じ込め核融合プラズマシミュレーションの現状と今後の展望

４．統合コードを使った実験解析・物理解明・運転シナリオ構築
4. Application of Integrated Code: Experiment Analysis,
Physics Understanding and Development of Operation Scenario
横 山 雅 之１，２），佐 竹 真 介１，２），本 多 充３），相 羽 信 行３），林 伸 彦３），矢 木 雅 敏３），
藤 田 隆 明４），村 上 定 義５），關 良 輔１，２），山 口 裕 之１），奴 賀 秀 男１）
YOKOYAMA Masayuki1,2）, SATAKE Shinsuke1,2）, HONDA Mitsuru3）, AIBA Nobuyuki3）,
HAYASHI Nobuhiko3）, YAGI Masatoshi3）, FUJITA Takaaki4）, MURAKAMI Sadayoshi5）,
SEKI Ryosuke1,2）, YAMAGUCHI Hiroyuki1）and NUGA Hideo1）
１）

自然科学研究機構核融合科学研究所，２）総合研究大学院大学，３）量子科学技術研究開発機構，
４）

名古屋大学大学院工学研究科，５）京都大学大学院工学研究科（原稿受付：２
０
１
９年５月２
０日）

前章までに紹介された統合コードによって，トカマクやヘリカル系研究の様々な課題に関する実験解析，物
理解明，運転シナリオ構築が行われている．本章では，それらの代表例をまとめて紹介する．得られた結果や物理的知見の重要性はそれぞれの参考文献を引用することとし，コードの統合がそれらの解析にいかに効果を発揮したか（統合の恩恵）を強調する形でまとめる．
Keywords:
integrated code, tokamak and helical system, experiment analysis, physics understanding, operation scenario

４．
１ 統合コードを用いたトロイダル回転分布予
測式は以下の通りである．
! %'
%
"
,
$!
,
(
% !, % (,!
)
&#*
+
&#+
*
# $'*-*!" ! ")
! !)
"
+
$ %
$,%
#
%
# * *

トカマクプラズマのトロイダル回転は，MHD 不安定性やペデスタル部における Edge Localized Mode（ELM）の

$
'
"#*&
"!"#*
!)
&$*
"# )
*
!# &
*)
*)
%$
"

抑制効果などと関係しているため，核融合炉の定常運転シ*

ナリオを考える上で重要な要素である．これを数値シミュ

%
!
!$
!&
"# *)
* *
"$

レーションによって予測するためには，様々な要因によっ

*%+
)
.

（１）

て生じるトルク分布を評価し，そのバランスによって決ま

ここで，!%# *(*)*)"$*"*は磁気面平均したトロイダルる回転分布を計算する必要がある．量子科学技術研究開発

角運動量密度であり，右辺は左から順に移流項，微視的乱機構（QST）と核融合科学研究所（NIFS）との共同研究に流による拡散とレイノルズ応力，,#$トルク，NBIトルク，

おいて，個々の物理機構で生じるトルクを個別のシミュ

誘導電場の効果であり，最終項が非軸対称性によって生じレーションコードで計算し，それらを結合してトロイダルる新古典トロイダル粘性（NTV）である．NTV の計算には

角運動量の輸送方程式を解く計算モデルを，統合輸送コー

３次 元 磁 場 配 位 用 に 開 発 さ れ た 新 古 典 輸 送 コ ー ド

ド TOPICS［１］を中心としたコード群によって構築した．

FORTEC-3D を適用した［２］．また，NBI トルクの評価にこの研究で特に注目したのは，現実のトカマク装置の磁場は OFMC コ ー ド［３］，誤 差 磁 場 の 影 響 を 含 め た３次 元に存在する微弱な非軸対称性によって生じるトルクの影響

MHD平衡の計算にはVMECコード［４］が使われている．乱である．摂動磁場の要因としては，トロイダルコイルによ流によるピンチ項はスケーリングによる半経験的なモデ

るリップル磁場，真空容器のポートや機器類による影響のル，乱流運動量拡散係数はモデル評価した熱拡散係数に仮他，ELM 制御のために外部コイルによって与えられる場

定したプラントル数を掛けて算出している．海外の先行研合もある．摂動磁場によるトロイダル回転への影響は多く

究では，NTV を近似的な計算法で評価しているのに対しのトカマク実験で確認されており，国際熱核融合実験炉て，FORTEC-3D ではドリフト運動論方程式を第一原理的

ITER や将来の核融合炉のようにプラズマ体積の大きい装に解いて NTV を高精度に評価するため，スーパーコン

置ではNBIによる外部トルク入射よりも重要になると予想

ピュータを用いた長時間計算が必要となったが，このようされている．

なマルチ・フィジクスの統合シミュレーションにおいて第一原理的な計算を結合した先駆的な研究例となった．図１

TOPICS が解くトロイダル角運動量密度の時間発展方程

National Institutes of Natural Sciences, National Institute for Fusion Science, SOKENDAI, 322-6 Oroshi, Toki 509-5292, Japancorresponding author’s e-mail: yokoyama@nifs.ac.jp

４４４!2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

図１

4. Application of Integrated Code: Experiment Analysis, Physics Understanding and Development of Operation Scenario

M. Yokoyama et al.

図２JT-60U における FST 挿入あり（右），なし（左）の場合のトロイダル回転分布の TOPICS による予測（破線が NTV
無視，実線が NTV の影響を入れた結果）と観測（点）の比較（文献［５］の図４（e）
，４（f）
から転載）
．

図３

ITER における径電場（右軸），様々な非軸対称磁場成分（TF：トロイダルリップル，FI：機器の磁性体，port：真空
容器のポート，TBM：テストブランケットモジュール）の影響による NTV（左軸）の径方向分布．“NBI”は NBI トル
ク入力（文献［６］の図３（a）
から転載）．

TOPICS を中心としたトロイダル回転分布予測のための統合輸送計算モデル（文献［５］の図１から転載）．

に，これらコードの統合概念図を示す．
JT-60U 実験では，トロイダルコイル間にフェライト鋼（FST）を挿入してリップル磁場を低減する改造を行ったが，その前後でトロイダル回転分布の変化が観測されていた．TOPICS を用いたシミュレーションにおいて，FSTを入れる前の回転分布は，NTV の影響を入れた方がそれを無視した計算に比べ観測結果をよく再現しており（図２），この統合モデルの妥当性が示された［５］．
ITER では 部 分 的 に テ ス ト ブ ラ ン ケ ッ ト モ ジ ュ ー ル（TBM）を導入した実験が検討されているが，そのような機器内の磁性体や真空容器のポートによる軸対称性の破れがプラズマ回転や閉じ込め性能に与える影響の評価が求められている．そこで，こうした誤差磁場を設計データから具体的に評価し，３次元 MHD 平衡磁場を VMEC で求めた
［６］．TBM による誤差磁場は設置位置に局在化しているが，その振幅はトロイダル磁場の１％程度に達した．この誤差磁 場 に よ る NBI 高 速 イ オ ン や DT 反 応 に よ る高速$粒子の軌道損失によるトルクを OFMC で評価したが，解生しないことを MHD 安定性解析に基づいて確認する必要

析した時刻ではLモードであったため，TBMの影響は無視がある．本原稿執筆時点で，ペデスタルの圧力分布を予測できるほど小さかった．一方，摂動磁場による NTV の大きする際に国内外で最も利用されているモデルは EPED

さは図３に示されるように，TBM による寄与がその他の［８，９］と呼ばれるものである．このモデルは，ペデスタルの要因のものと同程度に強く，かつプラズマ中心にまで影響

幅は輸送現象で，高さは type-I ELM の発生原因である理想が及ぶことがわかった．この NTV 分布を TOPICS に取り

MHD モード（ピーリングバルーニングモード）の安定性

込んで計算したところ，トロイダル回転は TBM によってで決定されるという発想に基づいて圧力ペデスタル分布を２割弱減速されると予測された．

決定する．EPED には数種類が存在するが，最初に提案さこの他にも，JT-60SA で予定されている摂動磁場印加実れた EPED1 モデル［８］では，圧力ペデスタルの幅 #はペデ

験において，プラズマ内部に共鳴有理面を持つ共鳴モード

スタルポロイダルベータ値 %&!&%$の平方根に比例する，とをかけた場合と非共鳴モードをかけた場合でのプラズマ応"!
#
いう実験的に知られている傾向を用いて #!!%&
（ここ
!
&%$

答を含めたMHD平衡計算とNTVの予測にも用いられていで ! は定数）として求め，理想 MHD モードが臨界安定とる［７］．また，ペデスタル部の MHD 安定性においてプラなるまで圧力ペデスタルの高さを上げていく（同時にズマ回転 の 影 響 を 評 価 す る た め に，TOPICS コ ー ド と%&"&%$が大きくなるため幅も広がる）ことで，type-I ELMMHD 安定性解析コードを結合した研究も進められている

発生直前の圧力ペデスタル分布を求めることができる．

（４．
２節参照）．

JT-60SA のペデスタル性能も当初はこの EPED1 モデルを

（佐竹，本多）用いて予測がされていたが，近年の研究成果によって，ペ４．
２ トカマクにおける H モード放電におけるペデスタル性能予測

デスタル領域のプラズマ回転およびイオン反磁性ドリフト

統合コードを用いたJT-60SA運転シナリオ構築の１例とれた［１０，
１１］ため，これらを考慮したペデスタル性能予測

効果を自己無撞着に考慮することが重要であることが示さ

して，H-mode 放電におけるペデスタル性能予測についてを行うことが求められた．

紹介する．ペデスタル性能は，コア領域のように乱流など

このような予測を実現するには，まずプラズマ回転分布による輸送現象だけで評価することができず，ELM が発を数値計算によって評価する必要がある．しかし，回転分
４４５

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

布を決定する物理機構は非常に多岐にわたるため，単一の昇する．上昇した SOL 温度・密度は，セパラトリックス端

物理モデルに基づく第一原理シミュレーションで評価するの径方向勾配を平坦化して拡散輸送を減らし，ELM エネ

こ と は 非 現 実 的 で あ る．そ こ で，前 述 の 統 合 コ ー ド

ルギー・粒子損失を低下させる．

TOPICS を用いて運動量のピンチ効果や拡散効果，中性粒

実験では ELM によるエネルギー損失は，ペデスタルプ

子ビームによるトルク入力，さらには必要に応じて３次元

ラズマの衝突周波数に依存する．図６は，シミュレーショ

平衡計算コード VMEC および新古典輸送シミュレーショ

ンで得られた依存性である．ここで，シミュレーション結

ンコード FORTEC-3D との連携計算による NTV を考慮し果を詳細に解析した結果，電子のエネルギー損失の衝突周た運動量バランスの式を解くことでプラズマ回転分布を予

波数依存性がブートストラップ電流と SOL 輸送に起因す

測することとした（４．
１節参照）．そして，EPED1 で採用さることがわかった．ブートストラップ電流は，衝突周波数れているペデスタル幅のスケーリング則はそのまま用い，の増加とともに減少し，ペデスタル領域の磁気シアを強め予測された回転分布を持つプラズマの MHD 安定性を拡張

る．増加した磁気シアは，不安定モードの固有関数の幅を

MHD 安定性解析コード MINERVA-DI を用いて評価する

狭め，その結果，ELM により増幅される輸送の範囲とセパことで，回転及びイオン反磁性ドリフト効果を考慮した

ラトリックス近傍の値を小さくする．一方，磁力線平行方

MHD 安定性で決定される圧力ペデスタル性能の予測を実向の電子の熱伝導が SOL 温度の上昇を決めるので，高衝突現した［１２］．この統合コードフレームワークを図４に示

周波数では熱伝導が低下，SOL の電子温度がさらに上昇

す．そ し て，こ の 拡 張 さ れ た EPED1 モ デ ル を 用 い てし，径方向勾配を平坦化させる．以上２つの物理機構によJT-60SA の運転シナリオ＃4.1 と呼ばれるもののペデスタ

り，電子のエネルギー損失が，衝突周波数の増加とともにル性能を再評価した結果，回転およびイオン反磁性ドリフ減少する．これらの機構は 1.5 次元輸送コードに付随するト効果を考慮しない場合に比べて圧力ペデスタル性能は１０％ほど低下したものの，同シナリオで設定されているプ
ラズマ全体の性能目標値を達成するために必要なペデスタ
ル性能は確保できることが示された［１２］．

（相羽）

４．
３ トカマクにおける H モードペデスタルとELM の統合モデリングH モードにおけるペデスタル形成と ELM の物理機構を解明するため，統合化シミュレーション研究を進めている．プ ラ ズ マ 性 能 に 対 す る ELM と ス ク レ ー プ オ フ 層（SOL）‐ダイバータプラズマの自己矛盾のない効果を明らかにするために，統合コード TOPICS で，図５に示すように 1.5 次元輸送コードを基に，ピーリングバルーニング
モード安定性コード（MARG2D）と SOL‐ダイバータモデ
ル（D5PM），中性粒子モデルとを統合化した．統合した
コードのシミュレーションにより，ペデスタルの密度・温度をガスパフにより変化させ，ELM の繰り返しサイクル

図５

TOPICS におけるペデスタル統合モデル．下部の囲みの領域はペデスタル形成と ELM のためのモデル．密度変化を考慮するために，炉心とダイバータ各々の中性粒子モデルを結合している．上部のモデルは ELM 緩和制御のためのペレットモデルである．

図６

H モードペデスタル部に蓄積されたエネルギーに対する，
ELM による損失割合（電子とイオンの寄与とその合計）のペデスタルプラズマの衝突周波数依存性．

におけるペデスタル成長と ELM 崩壊を模擬した．シミュレーションで ELM 崩壊が起きると，ペデスタルのエネル
ギー・粒子が SOL に流れ，SOL プラズマの温度・密度が上

図４

プラズマ回転・イオン反磁性ドリフト効果を考慮したペデ
スタル性能予測を実現する統合コードフレームワーク．

４４６

Special Topic Article

4. Application of Integrated Code: Experiment Analysis, Physics Understanding and Development of Operation Scenario

M. Yokoyama et al.

ブートストラップ電流モデルと SOL‐ダイバータモデルの衝突周波数を人為的に増倍しても図６と同様の傾向が得られることからも確認できた．このように，統合モデルのどの物理要素が得られた結果に大きく影響するのかモデルを人為的に修正操作して調べることは，物理要素が明確に分離されているが故にできる統合モデルの利点であり，物理機構解明の有効な手段である．上記の他，ELM によるイオンのエネルギー損失の衝突周波数依存性や粒子損失の衝突
周波数非依存性の物理機構を解明し，衝突周波数以外にもペデスタル肩内側の急峻な圧力勾配がエネルギー損失を増倍することをシミュレーションで明らかにし，実験データを解析して確認した［１３，
１４］．

図７

（林）

４．
４ トカマクにおける固体ペレット入射による
ELM 緩和のモデリング
ELMによる過渡的な熱負荷の低減がITERの最重要課題

ペデスタル蓄積エネルギーに対するELMによるエネルギー
損失割合と，ELM 誘起時のペレット位置（○□：JT-60U，●：ITER），ELM のモード幅（横線）のシミュレーション
結果．○●は ELM 制御に適した条件でペレットを入射した
場合，□はそれ以外の条件で入射した一例，影の領域はその他の JT-60U 結果．水平の点線は制御なしの自然発生
ELM の値，縦の破線はペデスタル肩位置．

の１つで あ る．ELM 緩 和 制 御 手 法 の１つ と し て 燃料ペレット入射が考えられている．この制御では，ELM によるプラズマ性能に影響を与えず粒子供給を少なくすることが４．
５ トカマクにおける LH 遷移の炉心・SOL／ダ
イバータ統合モデリング必要である．ペレットによる ELM の制御性の解明をめざ

出力や粒子制御は ITER や原型炉における重要な研究課

して，図５の下部に示されるペデスタル・ELM モデルに，題である［１８］．スクレープオフ層やダイバータ領域では燃上部のペレットモデル（APLEX）を結合した．ペレットモ

料イオン，不純物，中性粒子，及びプラズマ対向材料の相

デルは，ペレットの背景プラズマ中での溶発と形成された互作用過程が存在する．このような複雑なスクレープオフ

溶発プラズマ雲の運動，背景プラズマとの均質化を模擬す層やダイバータプラズマを理解するためには，シミュレー

る．ペレットの溶発プラズマは背景プラズマから加熱されション研究が不可欠であり，２次元のダイバータコードがてエネルギーを吸収し，電子熱速度の速い時間スケールで開発されてきた．SONIC は２次元プラズマ流体 コー ド

磁気面の電子温度を低下させ径方向に電子圧力勾配を作る

SOLDOR，中性粒子モンテカルロコード NEUT2D，不純物ことを想定し，上記のペデスタル・ELM モデルに結合しモンテカルロコード IMPMC から構成される統合コードでている．そのため，ペレットによる３次元の局所的摂動をある［１９，
２０］．この SONIC と 1.5 次元トカマク輸送コード考慮していないが，上記の径方向摂動による ELM 誘起機

TOPICS を結合し，プラズマコアと周辺部を矛盾無く統合

構は ASDEX Upgrade の実験観測［１５］から示唆されたもの

して LH 遷移のシミュレーション研究を行った．セパラト

である．この機構を実証するため統合モデルによるシミュ

リックス近傍磁気面を接合面として，SONIC 側境界条件にレーションを行い，ペレットのエネルギー吸収効果によりTOPICS から磁気面平均した粒子と熱フラッ ク ス を 渡

圧力分布に摂動が生じて ELM を誘起させることがわかっし，TOPICS 側境界条件に SONIC で評価した密度，温度，

た．ペレットで誘起される ELM によるエネルギー損失は，

粒子ソース等を返す．図８及び９に JT-60SA プラズマにおペレットの入射条件（ペレットの入射位置，大きさ，速度，

け る LH 遷 移 の シ ミ ュ レ ー シ ョ ン 結 果 を 示 す［２１］．

成長するペデスタルに対する入射タイミング）に依存する

TOPICS の輸送モデルには，径電場シア流による乱流抑制

ので，様々な条件で JT-60U と ITER 標準運転のプラズマの機構を取り入れている．図８に，TOPICS で計算した吸収シミュレーションを行い，ELM 緩和制御に適したペレッ

パワー，閉じ込め改善度 HH 因子，蓄積エネルギー，中心ト入射条件を調べた．その結果，図７に示されるように，

電子温度，線平均電子密度の時間発展を示す．３秒の時点

エネルギー損失を大きく低減するだけでなく，できるだけペレットがペデスタルに深く侵入してペデスタル肩近くで局在化した ELM を誘起させれば，エネルギー損失を大きく低減できることを明らかにした．この ELM 制御に適したペレット入射条件は，ペデスタル圧力が自然発生 ELMを起こす高さから５％程度低い時点で，ペデスタル蓄積粒
子数の数％程度をもつペレットを弱磁場からペデスタル頂上近くに届く速度で入射することである．ITER で必要なペレット粒子数は概ね現在の設計範囲内だが，その設計値
は単にペレットがペデスタル肩に届くように経験的に評価されており，今回の結果により，その設計に物理的背景を与えることができた［１６，
１７］．

図８

（林）
４４７

吸収パワー，HH 因子，蓄積エネルギー，中心電子温度，線平均電子密度の時間発展．

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019から追加熱を印加し，L モードから H モードへの遷移を自ランクコード EC-hamamatsu で評価する．異常熱輸送モデ

己無装着に実現したシミュレーション結果である．図９はルとして JT-60U や JET 実験で検証した CDBM モデルを用SONIC で計算した外側ダイバータストライク点における

いて温度分布を予測する．CDBM は，実験とほぼ同じ温度

電子密度，電子温度の時間発展である．TOPICS と SONIC

か，より低い温度を予測するため，予測はプラズマ性能のとの結合により，LH 遷移後の低温高密度ダイバータから下限を予測する保守的なものである．以上の統合モデルに

高温低密度ダイバータへ変遷する応答解析が初めて可能とよるシミュレーションで，NB と EC 加熱を合わせて合計 26なった．

MW程度のパワーで，温度分布に急峻な勾配がある内部輸送障壁をもつ高圧力（規格化圧力 "! !"!
!）定常運転プラ

プラズマコアと周辺部の統合化は欧米でも精力的に行わズマが得られた［２４，２５］．

れ て い る（FACETS［２２］，JINTRAC［２３］）．今 後 の 課 題としては，FACETS や JINTRAC とのベンチマークや不純

この運転シナリオではプラズマの密度が低く，不純物を物輸送までを考慮した統合化（例えば不純物輸送コードで入射して放射損失によりダイバータ板への熱負荷を低減す

ある IMPACT や TASK/TI と IMPMC や GYROIMP とのる必要がある．しかし，入射した不純物が炉心に蓄積しプ結合）を進めていく必要があろう．これらにより，非接触

ラズマ性能を低下させるおそれがある．これを調べるため

ダイバータ制御，燃料供給，ヘリウム灰廃棄等のシミュに，図１０の下部の統合モデルを用いる．しかし，統合ダイレーション研究へと展開していくことが期待される．

バータコード SONIC は 1.5 次元輸送コードに結合されたが，SONIC で計算時間がかかり連携計算で多くのパラメー（矢木）タランができる状態ではない．また，SONIC の不純物モン

４．
６ JT-60SA の高圧力定常運転プラズマの性能評価

テカルロコードIMPMCと炉心不純物輸送コードIMPACT

JT-60SA の高圧力定常運転シナリオのプラズマ予測を統授受だけができる状態である．このため，次に述べる方法

合コード TOPICS のシミュレーションで行っている．図１０で現状の統合モデルを保守的予測に用いる．上記の炉心性の上部に示されるように，中性粒子入射（NB）加熱・電流能評価時の加熱パワーを用いて SONIC でアルゴン（Ar）

駆動の評価は，1.5次元輸送コードに結合した粒子軌道追跡

不純物を入射してダイバータ板の熱負荷を低減し，その時の結合が済んでおらず，IMPMC から IMPACT へのデータ

モンテカルロコード OFMC か，より簡約化した モ デルの Ar のプラズマ最外殻磁気面での密度等の境界条件を評（フォッカープランクコード等）で行い，電子サイクロトロ

価した．そして，この境界条件を TOPICS に結合されたン（EC）波加熱・電流駆動は，光線追跡法・フォッカープ

IMPACT で用いて，炉心の Ar 蓄積量とプラズマ性能低下
が最大となる条件を設定しシナリオの成立性を評価した．
そのために，IMPACT での不純物粒子輸送は NCLASS で評価し，異常拡散輸送は新古典レベルの大きさで与えている．
この評価方法で得られた結果を図１１に示す．Ar を外側
バッフルから入射し，ダイバータ板の熱負荷が 10 MW/m2
以下になる SOL ダイバータプラズマが得られた．新古典の対流輸送により炉心には Ar の１６価から完全電離の１８価までが蓄積し，その他の低価数の Ar は境界領域に留まる．
元々壁材料で炉心に混入していた炭素の放射損失に加えてAr 蓄積で放射損失が増大したため電子温度が低下して，

図９

外側ダイバータストライク点における電子密度，電子温度の時間発展（文献［２１］の図４から転載）．

図１
１ JT-60SA の高ベータ定常運転シナリオにおける Ar 入射時の，（a）
ダイバータ領域の Ar 密度分布，（b）
外側ダイバー
タ板熱負荷分布，（c）
炉心での価数毎の Ar 密度分布，（d）
Ar を入射しない炭素 C のみの場合と比較した電子温度と放射損失分布．

図１
０ TOPICS における炉心プラズマ性能とダイバータ熱負荷評
価のための統合モデリング．矢印はデータ交換を表し，点線部分は未完成．

４４８

Special Topic Article

4. Application of Integrated Code: Experiment Analysis, Physics Understanding and Development of Operation Scenario

M. Yokoyama et al.

炉心プラズマ性能は下がる．しかし，この程度の放射損失析として，原型炉を想定した条件での内部輸送障壁プラズの増大は装置の NB の全パワーの範囲内で補えるレベル

マにおける不純物入射の影響［２９］や H モードプラズマにおで，炉心プラズマ性能を十分回復できることから運転シナ

ける不純物入射時の運転密度の電子密度分布形状依存性がリオ の 成 立 性 を 確 認 で き た［２６］．壁 材 料 が 金 属 と な る

調べられている．

（藤田）

ITER や原型炉では不純物入射は必須であり，用いた解析

４．
８ TASK3D を用いた LHD における高イオン温度プラズマの統合熱輸送シミュレーション

手法は ITER や原型炉のシナリオ構築にも適用できる．この手法は，もちろん統合モデルの結合を完成させ予測の精度を上げる必要はあるが，結合を完了して計算負荷の高い

LHD においては，垂直NBI加熱の増強により高イオン温統合シミュレーションをする前に行える統合モデルの有効度実験が進展している．高イオン温度プラズマの輸送特性な使い方の一例であろう．

を解明するため，ヘリカル型装置に対応した統合シミュ（林）

レーションコードTASK3Dの予測バージョン（TASK3D-a

４．
７ トカマクにおける 不 純 物 蓄 積 のトロイダル 回
転・電子密度分布形状に対する依存性［３０，３１］を用いた研究が進められている．TASK3D の輸送

不純物輸送解析に関して，トロイダル方向に回転するト計算モジュール TR には，輸送モデルとして新古典輸送おカマクプラズマにおける高価数不純物の内向きピンチのモ

よび異常輸送のモデルが組み込まれている．新古典輸送に

デル（PHZ［２７］）が TOTAL コードへ導入され，JT-60U

関 し て は DGN/LHD［３２］を 用 い，乱 流 熱 輸 送 に 関 し て実験の解析が行われた．PHZ は，粒子軌道の磁気面からのは，電子については gyro-Bohm モデルとし，イオンについずれにより不純物の価数が軌道に沿って変動することによては ITG 乱流を想定した gyro-Bohm-gradTi モデル（イオと 区 別 す る た め TASK3D-p と 呼 ば れ る こ と も あ る）

り，価数に反比例するトロイダルドリフトの小径方向成分

ン温度勾配 gyro-Bohm モデルと呼ぶ）を導入して い るの軌道平均が０とならずピンチを引き起こすものである．

［３３］．

この解析では，実験での密度分布，温度分布，トロイダル

!$
#
" #
"!
#$ " "#'#% """&!%& " & %，"'!%'
（２）
! '#"
!#
"$$
! # "'"$$

回転分布，径電場分布を固定としてタングステンイオンの輸送のみが解かれたが，PHZ がイオンの価数に依存するため，径方向輸送とピンチ速度の評価を並行して（反復して）これらの熱輸送モデルでは，最適な比例係数 C は一般的な行う必要がある．図１２に示すように，反プラズマ電流方向

実験データを用いて決定した．結果として，イオン温度勾のプラズマ回転の増大に伴ってタングステンの蓄積が増加

配 gyro-Bohm モデルでは，高イオン温度放電に対してイオ

するという実験結果と同様の傾向が得られた［２８］．ただン・電子温度ともよい一致が得られた［３１］．

し，本解析では文献［２７］で示されたもう一つのピンチモデ

LHD においては，高イオン温度実験において 10 keV にル（Er ピンチ）の仮定が実験で満たされなかったため，評到達する高イオン温度を達成している［３４］．この LHD 高

価されなかった．旋回中心軌道を数値的に計算することにイオン温度モードは，将来のヘリカル型核融合炉設計におよ り PHZ，Er ピ ン チ を 評 価 す る コ ー ド が 開 発 さ れ ていて重要な論点となると考えられ，その物理的解明が必要

TOTAL へ組み込まれ，解析が行われている．一方，電子となっている．一方，LHD における高イオン温度モードへ

密度勾配による不純物イオンの内向きピンチに着目した解の遷移は，炭素ペレット入射による急激なプラズマの変化を伴うため，その解明には，短時間での密度変化や不純物の加熱・輸送への影響など複雑な物理量の変化を総合的に考慮する必要がある．このため，統合熱輸送シミュレー
ション（TASK3D+GNET-TD）を行うことにより，LHD
高イオン温度プラズマの炭素ペレット入射による輸送低減の物理的な解明を行った．
まず，高エネルギー粒子閉じ込め解析コード GNET
［３１］
を基に時間発展プラズマに対する正確なNBI加熱解析を行うことができる GNET-TD コード［３５］を開発し，NBI ビー
ムイオンの速度空間分布および加熱吸収分布の時間発展を求めた．結果として，ペレット入射によるプラズマ変化に対応するNBIビーム分布の変化および加熱分布の変化を高

図１
２ PHZ モデルにより評価されたタングステン蓄積のトロイダ
ル回転依存性．横軸は磁気軸におけるトロイダル回転速度，縦軸は磁気軸における電子密度に対するタングステン
イオン密度の比．白抜き菱形は実験値．解析は A~E の５ケースの実験データについて下記の４通りの仮定で行われた：プラス印は PHZ なしで異常拡散係数が 0.01 m2/s の場
合，丸，三角，白抜き四角は PHZ ありで異常拡散係数がそ
れぞれ 0.01 m2/s，0.001 m2/s，0.0001 m2/s の場合（文献［２８］の図７から転載）．

精度に求めることが可能となり，ペレット入射後のビーム
粒子の急激な減衰や加熱分布の変化を明らかにすることができた．また，ペレット入射後に中心部のイオン加熱が増加していることが示された．
次に，得られたNBI加熱吸収分布を用いてTASK3Dにより統合熱輸送シミュレーションを行った．ここでは，電子，軽水素イオン，ヘリウムイオン，炭素イオンの４種の粒子
４４９

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

種を 仮 定 し，シ ミ ュ レ ー シ ョ ン を 行 っ た［３６］．図１３は，LHD 重水素実験では，プラズマの温度に対して中性粒子

TASK3D により得られた炭素イオン温度の時間発展をペビームのエネルギーが高く，かつ高電力であるため，核融レット入射後の時間に対して示したものである．ここで合反応によって生じる中性子線計測を用いたプラズマ内部

は，乱流輸送モデルに，不純物による改善度（実効電荷数の高速イオン密度の推定が可能である．TASK3D-a におけ!!""依存性）を含むモデルを導入している（詳細は文献［３８］る中性粒子ビーム入射加熱モジュールの一つであるを参照されたい）．結果として実験結果と同様な比較的速

FIT3D
［４０，
４１］においても，計測結果との比較を通じた

い時間スケールでのイオン温度上昇とその後の降下が得ら

コードの検証が行われている．FIT3D は，高速中性粒子のれた．これらの結果は，実効電荷数に関連して，乱流輸送

プラズマ中でのイオン化位置の評価（HFREYA），短時間が改善される物理機構の存在が示唆される．

の軌道追跡による初期イオン損失評価（MCNBI）と，定常（村上）的な高速イオン密度および加熱パワーの評価を行う解析的

４．
９ TASK3D-a を用いた LHD プラズマの高エネ
ルギー粒子解析なフォッカープランク計算（FIT）を連結したコードであ

京都大学と核融合科学研究所の共同研究を通じて，LHD

取り入れることが難しい高速イオンの軌道幅効果を中性粒実験解析のための統合輸送解析ツールTASK3D-a［３７］が構子ビーム加熱解析に取り入れることができ，かつ，フォッ

築されている．TASK3D-a によって，これまでに約２，０００

カープランクによる高速な加熱解析が可能となっている．

ショットに対して輸送解析が行われており，LHD 重水素実

FIT3D は，高速イオンの減速過程における粒子損失をほと

験での同位体効果の解明［３８］や，得られた多数の輸送係数

んど考慮していないため，加熱パワーを過大評価しているから実験的な輸送モデルの構築［３９］が進められている．

ことが想定されている．これまでに，不純物を考慮しない

る．この連結によって，一般にフォッカープランク解析に

特に，中性粒子ビーム入射加熱モジュールは，LHD 重水

場合，LHD の標準的な磁場配位における FIT3D の計算結素実験によって定量的な評価が可能となったことから，実果は中性子発生率の計測結果と比べて２倍程度の過大評価

験検証が精力的に進められている．図１４に TASK3D-a におをしていることがわかっている［４１］．

ける中性粒子ビーム入射加熱モジュールの計算フローを示

ペレット入射のように密度が大幅に変化する場合など

す．各モジュールの役割については，次の段落に示す．

の，フォッカープランクの定常解を使用するFIT3Dでは加熱パワーが評価できない場合に対応して，TASK3D-a では
加熱の時間発展を求めることのできる CONV̲FIT3D［４２］の実装や検証も行われている．CONV̲FIT3D は FIT3Dのフォッカープランク部分を置き換えることを目的としたコードであり，高速イオンの減速方程式を解くことによって，中性粒子ビームによる高速イオン密度および加熱パ
ワーの時間発展を得ることができる．CONV̲FIT3D は，
高速イオンの減 速 過 程 を 評 価 す る こ と が 可 能 で あ り，
HFREYA，MCNBI と連結することで，中性子発生率の減
衰時間を LHD 実験と検証することが可能になっている．
図１５は，短パルスの中性粒子ビーム入射時における中性

図１
３ TASK3D により得られた炭素ペレット入射後の温度（４つの規格化小半径位置）の時間発展（文献［３６］の図に追記改訂）．

子発生率の計測結果と CONV̲FIT3D による計算結果を比較したものである．ここでは，プラズマ中の不純物の存在を考慮していないこと，高速イオン減速過程における粒子損失を考慮していないことなどの要因によって，計算結果

図１
５ 中性粒子ビームの入射パワー（上図）と中性子発生率（下図）の時間変化．下図の実線が計測結果で点線が CONV_FIT3Dによる計算結果である．

図１
４ TASK3D-a における中性粒子ビーム入射加熱モジュールの計算フロー．

４５０

Special Topic Article

4. Application of Integrated Code: Experiment Analysis, Physics Understanding and Development of Operation Scenario

M. Yokoyama et al.

は過大評価になっている．この計測結果と計算結果との差

不純物イオンを含む各粒子種の温度・密度分布といったから，高速イオンの実効的な粒子閉じ込め時間とプラズマデータが必要となる．統合輸送解析スイート TASK3D-aの有効電荷数を推定することが可能であり［４２］，この推定の整備により，これらの解析・実験データはシミュレー

値を用いて中性子発生率の再評価を行った結果が図１６であション担当者がすぐ利用できる形式で自動的に準備される

る．その結果，計測結果と近い計算結果が得られたことかようになってきている．NBI 加熱に関しては，FIT3D モ

ら，もっともらしい計算結果が得られていると考えられジュールによる半自動的な加熱解析が可能となっている．

る．CONV̲FIT3D では，このような高速イオンの減速過FIT3D は案内中心軌道の取り扱いが簡略化された簡易的程の検証を通して，実効的な高速イオン閉じ込め時間のモなモジュールであるが，その主要な入力データは，完全なデ ル の 開 発 も 行 わ れ て お り，開 発 さ れ た モ デ ル は

案内中心軌道を取り扱う GNET と共通している．このたCONV̲FIT3D に取り入れられる予定である．

め，FIT3D による簡易的な解析の後，重要な放電・タイミ

また，３次元フォッカープランクコードである TASK/

ングのデータを個別に取り出し，GNET を用いてより正確

FP［４３］は実空間１次元，速度空間２次元空間において，中な加熱解析を個別に行うという一連のプロセスが，比較的性粒子ビーム由来の高速イオン速度分布の時間発展を計算スムーズに実行できるようになっている．

することが可能であり，TASK3D-a での実装が進められて

軽水素実験と重水素実験で得られた高イオン温度プラズ

いる．TASK/FP では，CONV̲FIT3D や FIT3D と異なり，マのイオン熱輸送の比較において，統合輸送解析スイート

高速イオン速度のピッチ角を考慮していることにより，高が可能としたこのプロセスを活用した．実験データに基づ

速イオンどうしの衝突による核融合反応の寄与を求めるこき，GNET を用いて複数イオン種を考慮した加熱解析をとができる．TASK/FP では速度空間２次元での高速イオ行った．その結果，軽水素で 8 keV，重水素で 10 keV をそン速度分布を評価できるため，荷電交換分光などの高速イ

れぞれ達成したプラズマにおいて，NBI 加熱によるイオン

オンの速度空間の情報も含めた計測との比較検証が可能でへの加熱吸収量は，軽水素と重水素とで大きな差がないこある．TASK/FP においても，FIT3D と同様に HFREYA，とがわかった．この結果から，重水素プラズマではイオンMCNBI と連結し，高速イオンの軌道幅効果を取り入れたの熱輸送が小さくなっていることが示唆された．重水素実実験との比較検証が可能になっている．

験では，イオン加熱にとって重要な垂直入射ビームの入射（關，奴賀）

エネルギーが増加したことにより，高速イオンの軌道幅が４．
１０ LHD 重水素実験における NBI 加熱の高精度解析と，TASK3D-a との連携大きくなり，捕捉粒子の損失が増えたことで，加熱吸収とLHD 重水素実験の重要課題として位置づけられている

重水素プラズマそれぞれに対して磁場配位依存性を比較し同位体効果の解明には，熱輸送解析において加熱吸収を正た一連の放電の解析では，磁気軸を内寄せするほど加熱吸確に評価することが重要である．ヘリカル系においては，収が増加するという結果が得られている．これは，磁気軸磁場配位の３次元性により，高エネルギー粒子の軌道は非位置の内寄せによる捕捉粒子軌道の改善から予想される結常に複雑な振る舞いを示す．高速イオンの大域的な案内中

果と矛盾しない．

しては軽水素時と同程度になった可能性がある．軽水素，

心運動を考慮して正確なNBI加熱解析が可能なシミュレー

前述のように，GNET による加熱解析に必要な入力デー

ションコードとしては，５次元位相空間における分布関数

タ群の多くは，統合輸送解析スイートによって準備されるを 求 め る こ と の で き る モ ン テ カ ル ロ コ ー ド（GNET

ようになっている．一方で，計算の実行自体は，大型並列［３１，３５］，MORH［４４］）が整備されてきている．GNET は，計算機等において担当者が手動で行っている．今後は，重水素実験へ向けて複数イオン種やプラズマ時間発展に対GNET を解析システムに組み込み，その実行まで含めて自応する拡張と，統合輸送コード TASK3D との連結も進め動化することで，ヘリカル系における輸送物理の解明をさられてきた（４．
８節参照）．

らに加速することが期待される．

（山口）

LHD 実験における特定の放電・タイミングに合わせた

４．
１１ おわりに

シミュレーションを行うためには，対応する３次元平衡，本章では，トカマクやヘリカル系研究における統合コードを使った実験解析・物理解明・運転シナリオ構築の具体例について，代表例をまとめて紹介した．トカマクプラズ
マ解析への３次元効果の導入，ペデスタルと ELM の統合
モデリング，炉心・SOL／ダイバータの統合モデリング，
JT-60SA や ITER に向けた運転シナリオ構築，不純物蓄積とトロイダル回転，密度分布形状の統合解析，ヘリカル系
においては，LHD の高イオン温度プラズマの輸送解析，高エネルギー粒子挙動の解析，熱バランスの自動解析システ

図１
６ 中性子発生率の時間変化．実線が計測結果で，点線が計算
結果．計算には図１５に示した計測結果と計算結果の差から推定される実効的な高速イオン閉じ込め時間と有効電荷数を用いた．

ムの構築など，非常に多岐にわたる課題に対して，統合
コードの開発と運用が進展していることがご理解いただけたものと期待する．さらに，階層間結合，大規模シミュ
４５１

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019［１
７］N. Hayashi et al., Contrib. Plasma Phys. 54, 599 (2014).
［１８］A. Loarte et al., Nucl. Fusion 47, 203 (2007).
［１９］K. Shimizu et al., Nucl. Fusion 49, 065028 (2009).
［２０］T. Takizuka, AIP Conf. Proc. 1237, 138 (2010).
［２１］M. Yagi et al., Contrib. Plasma Phys. 52, 372 (2012).
［２２］J.R. Cary et al., J. Phys. Conf. Series 125, 012040 (2008).
［２３］M. Romanelli et al., Plasma Fusion Res. 9, 3403023 (2014).
［２４］N. Hayashi et al., Nucl. Fusion 57, 126037 (2017).
［２５］N. Hayashi et al., 27th IAEA Fusion Energy Conf. (Gujarat, India, 2018), PPC/1-1.
［２６］N. Hayashi et al., Nucl. Fusion 58, 066001 (2018).
［２７］K. Hoshino et al., Nucl. Fusion 51, 083027 (2011).
［２８］Y. Shimizu et al., Plasma Fusion Res. 10, 3403062 (2015).
９］T. Yamakami, et al., Plasma Fusion Res. 9, 3403091 (2014).
［２［３０］A. Sakai et al., Plasma Fusion Res. 10, 3403048 (2015).
［３１］S. Murakami et al., Plasma Phys. Control. Fusion 57,
054009 (2015).
［３２］A. Wakasa et al., Jpn. J. Appl. Phys. 46, 1157 (2007).
［３３］A. Wakasa et al., Proc. 39th EPS Conf. and 16th Int. Cong.
Plasma Phys. (2012) P2. 028.
［３４］H. Takahashi et al., Nucl. Fusion 58, 106028 (2018).
［３５］H. Yamaguchi et al., Nucl. Fusion 56, 026003 (2016).
［３６］S. Murakami et al., Proc. IAEA Fusion Energy Conference 2014, St Petersburg, Russia (2014) TH/P6-38.
［３７］M. Yokoyama et al., Nucl. Fusion 57, 126016 (2017).
［３８］K. Tanaka et al., submitted to Nucl. Fusion (2019).
［３９］M. Yokoyama, Nucl. Fusion 59, 094004 (2019).
［４０］S. Murakami et al., Trans. Fusion Technol. 27, 256 (1995).
［４１］R. Seki et al., Plasma Fusion Res. 14, 3402126 (2019).
［４２］H. Nuga et al., Plasma Fusion Res. 14, 3402075 (2019).
［４３］H. Nuga et al., Nucl. Fusion 59, 016007 (2019).
［４４］R. Seki et al., Plasma Fusion Res. 5, 027 (2010).

レーションとの連携（簡約化モデルや機械学習（第５章参照）の導入など）などを通じて，統合コードの機能拡張や解析の高精度化，迅速化を図っていくことが必要である．
そのために，BPSI（Burning Plasma Simulation Initiative）［http://bpsi.nucleng.kyoto-u.ac.jp/bpsi/］，統合コード研究会（九州大学応用力学研究所共同研究），核融合科学研究
所一般共同研究，核融合エネルギーフォーラムのクラス
ター会合など，多様な機会における共同研究や情報交換が不可欠である．

参考文献［１］M. Honda et al., Nucl. Fusion 54, 114005 (2014).
［２］S. Satake et al., Plasma Phys. Control. Fusion 53, 054018(2011).
［３］K. Tani et al., J. Phys. Soc. Jpn. 50, 1726 (1981).
［４］S.P. Hirshman et al., Comput. Phys. Commun. 43, 143(1986).
［５］M. Honda et al., Nucl. Fusion 55, 073033 (2015).
［６］M. Honda et al., Nucl. Fusion 57, 116050 (2017).
［７］M. Honda et al., Nucl. Fusion 58, 112012 (2018).
［８］P.B. Snyder et al., Nucl. Fusion 49, 085035 (2009).
［９］P.B. Snyder et al., Phys. Plasmas 19, 056115 (2012).
［１０］N. Aiba et al., Nucl. Fusion 57, 022011 (2017).
［１１］N. Aiba et al., Nucl. Fusion 57, 126001 (2017).
［１２］N. Aiba et al., Plasma Phys. Control. Fusion 60, 014032(2018).
［１３］N. Hayashi et al., Nucl. Fusion 49, 095015 (2009).
［１４］N. Hayashi et al., Nucl. Fusion 51, 073015 (2011).
［１５］T. Szepesi et al., Plasma Phys. Control. Fusion 51, 125002(2009).
６］N. Hayashi et al., Nucl. Fusion 53, 123009 (2013).
［１

４５２

J. Plasma Fusion Res. Vol.95, No.9 (2019)４
５３‐４５７

小特集 統合コードによる磁場閉じ込め核融合プラズマシミュレーションの現状と今後の展望

５．統合コード開発の今後の展望
5. Future Prospects of Integrated Code Development
本 多 充，成 田 絵 美，林 伸 彦，矢 木 雅 敏，
福 山 淳１），村 上 定 義１），横 山 雅 之２，４），藤 田 隆 明３）
HONDA Mitsuru, NARITA Emi, HAYASHI Nobuhiko, YAGI Masatoshi, FUKUYAMA Atsushi１）,
MURAKAMI Sadayoshi１）, YOKOYAMA Masayuki２，４）and FUJITA Takaaki３）
量子科学技術研究開発機構，１）京都大学大学院工学研究科，２）自然科学研究機構核融合科学研究所，
３）

名古屋大学大学院工学研究科，４）総合研究大学院大学（原稿受付：２
０
１
９年５月２
０日）

これまでの各章で，統合コード開発の意義，構成要素の物理的背景，国内外の開発状況と，実験解析や運転
シナリオ構築への応用例などを見てきた．統合コード開発は巨視的な時間スケールにおけるプラズマの振る舞いを包括的に記述することを目的に，プラズマの持つ幅広い時空間スケールに渡る現象を階層ごとに分離したの
ち，それらを適切に再度組み上げる手法が採られてきた．現在でもプラズマ放電の全時間予測を現実的な時間で行う唯一の方法である一方，統合化されるモデル・コード数が増大し，各要素の物理の高度化に伴い計算時間も大きく増加する傾向にある．その問題の克服も念頭に，現在各極では伝統的な手法に囚われないアプローチに基づくモデル・コード開発が進展している．これらは従来型統合コードを完全に代替，置換するものではなく，取り組む問題に応じて使い分けされるものであり，両者は車の両輪として開発が強力に進められると考えられる．
Keywords:
integrated code, speed-up technique, real-time simulation, stiff transport model, gyrokinetic simulation, deep learning,kinetic modeling５．
１ はじめに算終了までに膨大な時間が掛かるか，悪くすれば物理的に第３章では，国内外で開発されている統合コードとその妥当な解が得られないという結果となってしまう．従来の数値実装について詳細に見てきた．それらの基本的な設計

枠組みを超えた時空間離散化の手法や新しいアイデアに基思想は，輸送ソルバを核として様々な物理モジュールが結づく実装によって，これらの問題を克服しようとする試み合しているモジュール方式とするものである．統合コードが最近活発に行われている．また，近年著しく興隆していによって得られる結果のるデータ駆動科学技術の一つである機械学習を用いて，高洗練度は各物理モジュールのそれに強く依存しているため，各物理モジュールの物理を度化した輸送モデルの応答を人工ニューラルネットワーク

高度化する方向に開発が進み，それに伴い計算時間の増大（ANN）に模擬させることで，輸送モデルを直接用いたとも招いた．統合コード開発の重要な目的の一つに，プラズ

きと比べ計算を格段に高速化させる試みも進んでいる．本

マ放電の全時間シミュレーションがあり，一つの放電の全章ではこのような国内外の取り組みについて，とりわけ数時間予測に膨大な時間を費やしていては，実際のプラズマ

値スキームに着目して紹介する．

放電実験前の放電予測はおろか，事前の多様な運転シナリ

加えて，第２章で紹介された波動加熱や中性粒子ビーム

オ開発もままならない．

加熱などを扱う加熱コードとの結合を目的とした新しい輸モデル高度化に伴う弊害の一つに，硬い輸送モデルの取送コードの開発も進んでいる．従来の輸送ソルバはプラズ扱いに関する問題がある．硬い輸送モデルとは，分布硬さ

マを流体と見なした輸送方程式を解いてきた．しかし，波と呼ばれる，ある臨界勾配値を超えて分布勾配が急峻化し動加熱や中性粒子ビーム加熱などは速度空間を局所的に歪た際に乱流流束が非線形的に急増する現象を再現する輸送めるため，マクスウェル分布を仮定した流体近似とは本質モデルである．ジャイロ運動論シミュレーションなどで現的に相容れない．従来の体系では，フォッカー・プランクれる現象を忠実に再現しようとした結果，第一原理計算とコードなどを用いて加熱された分布関数を速度空間上で衝同様の分布硬さが生じ，この乱流流束の勾配に対する非線突緩和させ，流体近似された輸送コードに取り込めるよう

形応答が輸送シミュレーションを収束させることを困難ににしてから結合していた．しかし，そもそも速度分布関数していた．言い換えると，高度な輸送モデルを用いると計そのものを従属変数としてプラズマの時間発展を解いていNational Institutes for Quantum and Radiological Science and Technology, Naka, IBARAKI 311-0193, Japancorresponding author’s e-mail: honda.mitsuru@qst.go.jp

４５３!2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

れば，加熱コードと直接結合することが可能であり，また温度分布を再構成している．密度分布やトロイダル回転分

非マクスウェル分布的な現象も扱うことができる．このよ

布 も 類 似 の 手 法 で 計 算 し て い る．数 値 実 装 の 大 半 はうなアイデアに基づく運動論的統合モデリングが現在日本

MATLAB で行われているが，高速化のために一部ルーチで進んでおり，その紹介も行う．

ンは C や Fortran で書かれている．これらの工夫によって，
ITER 標準運転シナリオの１
０００秒に渡るシミュレーション

５．
２ 実時間計算実現に向けたコード開発を約１分で終えることができる［４］．

２分程度のプラズマ放電実験の放電間に統合コードで次５．
３ 高度な乱流モデルや乱流シミュレーション
コードとの結合に向けたコード開発の放電の予測を行いたい場合，長くとも３０分程度で計算を終える必要がある．また，入力パラメータを少しずつ変えながら特定の物理量に対して運転シナリオを最適化したい

従来型の輸送ソルバに QuaLiKiz［５］や TGLF［６］のよう場合には多くの試行が必要となり，計算時間が短いに越しな硬い輸送モデルを組み込むと，密度や温度の勾配に対すたことはない．従来型の大規模な統合コード計算ではそのる乱流熱流束の非線形的な強い依存性のために定常解を得ようなニーズに対応できないため，従来型統合コードからることが困難になる．以下，本節では熱輸送に関して説明簡素化できるものは簡素化し，実時間を超える高速化をめするが，粒子や運動量に関しても同様である．加熱によっざした統合モデルが欧州で開発されている．

て温度勾配が増加したときに熱流束が急増すると，次の時

RAPTOR（RApid Plasma Transport simulatOR）［１］は間ステップで温度勾配が急減する．これを繰り返すことに

MATLAB で実装された１次元の輸送コードであり，簡易

よって振動解に陥ってしまう．この問題を解消する最も簡的な加熱モデル，鋸歯状振動や新古典テアリングモードの単な手法は時間刻み幅を１０−６秒程度まで小さくすることで物理モデルなどが結合されている．平衡は読み込んではいあり，そうすることで時間ステップあたりの温度勾配と熱るものの，発展は原則として解いていない．TOPICS 同様，

流束の変動幅を小さく抑えることができる．しかし，時間

解釈モード［２］と予測モードがあるが，本稿では予測モー刻み幅を１０分の１にすると単純計算で計算時間は１０倍になドのみ扱う．輸送コードの支配方程式は ASTRA に完全にるため，この手法は有効な解決策とは言えない．従来型の倣っているが，電流制御とプラズマ分布制御を矛盾なく行

輸送ソルバの枠組みを維持したまま硬い輸送モデルによるうことを当初の目的とした経緯もあって，初めはポロイダ

数値振動を取り除く方法が二つほど提案されているが

ル磁束と電子温度のみの分布発展を解くソルバのみが導入［７，８］，現実の体系で適用すると期待した結果が得られなされていた．現在では，イオン温度や密度にまで拡張されいこともしばしばある．本節では，分布硬さを有した輸送ている［３］．輸送方程式は時間方向に完全陰解法，空間方向

モデルを用いる際に生じる問題を克服するため開発されたには有限要素法で離散化されており，基底関数には３次 B

輸送コードや，局所ジャイロ運動論コードとの直接結合をスプラインが用いられている．３次 B スプラインを用いる前提に開発されたものについて概観する．

ことで２階微分まで連続となるため，ポロイダル磁束の２FASTRAN（FASt TRANsport）［９］は米国で開発されて階微分で記述される電流分布の径方向連続性を保証していいる１次元輸送コードであり，運転シナリオ開発時にはる．完全陰解法と B スプライン基底関数の採用に加え，関IPS-FASTRAN［１０］という統合モデルの核として運用され数形で書かれた輸送モデルを使用することで，大きな時間ている．FASTRAN では硬い輸送モデルを安定に解くため

刻み幅（〜0.1 s）や少ない径方向グリッド点数（〜１１点）なに，空間離散化に局所補間微分オペレータ（IDO）法を採

どの工夫が可能になり，実時間レベルの高速化が実現して用している．変数となる温度分布を，エルミート補間関数

いる．近年では，後述する機械学習を利用して高速化されを用いた５次の多項式で近似しており，１階微分値は多項

た輸送モデルを用いたシミュレーションも活発に行われて式の中にそのまま係数として残す一方，高階微分の係数はいる．

１階微分値と温度そのものから算出している．温度と温度

METIS（Minute Embedded Tokamak Integrated Simula-

勾配を定めなくてはいけないため，熱拡散方程式とその空tor）は０次元のスケーリング則と１次元の電流拡散を組

間微分式を有限差分化した式を連立して解いている．硬いみ合わせた輸送コードである［４］．METIS では RAPTOR

輸送モデルは勾配値に鋭敏に依存するため，勾配値（１階と同様，電流拡散を最も重要な物理量として捉えており，

微分値）を従属変数として一意に非線形収束で決定してしそのため電流拡散と軸対称平衡のカップリングは通常のまうことで安定に解くという試みである．一方，熱拡散方

1.5 次元の枠組みで扱っている．但し，高速化のために電流程式を微分した式も解かなくてはならないことで，必然的拡散方程式の径方向グリッド点 数 を２
１点 に 抑 え，Grad-に式中に拡散係数の２階微分項が存在する．FASTRAN

Shafranov 方程式のモーメントを取って得られる３つの常では拡散係数やソース項などは３次スプライン補間してお微分方程式を解くことで平衡を構築している．乱流熱輸送

り，それによってそれらの微係数を算出している．しかし，予測にはまだ大きな不確実性が含まれていることや輸送モ

硬い輸送モデルの出力する熱流束・拡散係数は実空間におデルの持つ強い分布硬さの取り扱いの困難さから，熱輸送

いて必ずしも連続的で滑らかな振る舞いを示すものではなに関しては大きく割り切って，時間発展については閉じ込

いため，これをスプライン補間によって連続かつ滑らかにめ時間スケーリングに基づく熱蓄積エネルギーのみを解してしまうことが妥当なのか，またその効果によって安定

き，それに合うような形で拡散係数を決め，電子・イオンに解けるように な っ た の で は な い か と の 疑 問 も 残 る．
４５４

Special Topic Article

5. Future Prospects of Integrated Code Development

M. Honda et al.

FASTRAN は温度のみならず密度やトロイダル回転，ポロ間方向に有限要素法で離散化されており，基底関数にエルイダル磁束も扱うことができ，既に DIII-D の運転シナリオ

ミート補間関数を用いている．そのため空間３次精度を有開発に応用が進んでいる［１０］．

している他，従属変数とその微係数も同時に定まるという

TGYRO［１１］は元々は局所ジャイロ運動論コード GYRO

特長を持つ．通常，ジャイロ運動論コードとの直接結合を［１２］とドリフト運動論コード NEO［１３］を直接結合してプめざす輸送コードは，ジャイロ運動論計算にかかる数値的

ラズマ分布を予測するために開発された枠組みであり，１コストからせいぜい１０点程度しか径方向グリッドを取れな次元の定常輸送方程式を解くソルバでもある．局所ジャイい．その程度の空間点から作られた温度分布に基づく温度

ロ運動論コード一般に言えることだが，GYRO や NEO は

勾配は極めて粗いものになり，温度の少しの変化で温度勾与えられたプラズマ分布に従う輸送流束を評価することが配が大きく変化する．GKVなどの局所ジャイロ運動論コーできる一方で，その輸送流束を受けてプラズマ分布がどの

ドは温度などの勾配に鋭敏な依存性を持つため，計算の安ように変化するかを予測することはできない．自己無撞着

定性に大きな問題が生じる．TRESS では温度勾配の算出

な分布予測のためには輸送流束に基づく輸送を解く繰り返

手法に任意性が含まれず，温度と無撞着に定まるため，少し計算の必要があり，そのため TGYRO が開発 された．

ない空間点数しか持たないジャイロ運動論コードとの直接TGYRO がコード開発のコンセプトとしてこれまで紹介し結合に適した実装と言える．時間方向には陰解法か２次のた他の輸送コードと異なるのは，ジャイロ運動論コードと後退微分が用いられる．安定化手法としてニュートンイタの直接の結合を念頭においていることと，もう一つは定常

レーション［８］も取り入れられている．乱流輸送レベルを

輸送方程式を解いていることである．輸送係数が温度等の決める飽和モデルが簡易的なものであるため，まだ実験と非線形の関数となっている定常輸送方程式を解くのは困難の比較を行っておらずコードの開発途上にあるが，捕捉電である．そこで，TGYRO ではニュートン法による繰り返子モードも取り入れ二粒子種の定常に至るまでの計算に成し計算で収束させている．解く方程式は，定常輸送方程式

功するなど［１７］，コードの動作面における実証は終えていを体積積分した方程式であり，ソース項などから事前に定る．

まっている目標流束に対して，乱流と新古典による流束が硬い輸送モデルを用いて定常温度分布を求めるために開合致するように温度分布を調整する．この作業を一般にフ発された定常輸送コード GOTRESS（Global Optimization

ラックスマッチングと呼ぶ．ニュートン法の過程でヤコビ

version of TRESS）［１９］は，定常輸送方程式の体積積分方程行列式を評価する必要があり，その際に熱流束に対する温式を解いているという点で TGYRO と類似しているが，方度勾配の微分を計算しなくてはならない．それを前進差分

程式の解法に大域的最適化手法を用いている点で大きく異で評価しているため，追加の熱流束計算が必要となる．論なっている．大域的最適化とは，ある関数の制約条件を満文［１１］で示された例では，DIII-D の 放 電 を 再 現 するシ

たす範囲内において最大値を見つける手法であり，代表的ミュレーションにおいて，GYRO と NEO をそれぞれ３
８４回な手法に遺伝的 ア ル ゴ リ ズ ム や 焼 き な ま し 法 が あ る．

呼び出し，１，５３６コアで５８時間の計算ののち収束解が得ら

GOTRESS は大域的最適化手法として遺伝的アルゴリズムれている．TGYRO は通常単体ではなく，OMFIT（Oneとネルダー・ミード法の両方を用いており，それらを適宜Modeling Framework for Integrated Tasks）［１４］という統使い分けながら フ ラ ッ ク ス マ ッ チ ン グ を 行 っ て い る．

合モデルフレームワークの中で運用されている．OMFITGOTRESS が解を見つけるアルゴリズムは，端的に言えば

は温度や平衡の予測計算のみならず，MSE計測に基づく平目標流束"と合致する輸送流束#"$!
"!$ #を作り出す温度!

衡再構成などの実験解析など，物理コードを用いるあらゆ$と温度勾配!
"!$ の組み合わせを，"$!
"!"#のパラメータ!

る作業工程を一元化した枠組みであり，実験データベース

空間において!
を最大化させることで求め
"
"!#"
$!
"
!$ #
!
$
$

や他の統合モデルとのインターフェースも備えている．

ており，最大化のために大域的最適化手法を用いている．

局所ジャイロ運動論コードとの直接結合をめざしたコー

TGYRO では規格化温度勾配のみを従属変数として解き，

ド開発は TGYRO 以外にも行われており，よく知られてい勾配を積分する事で温度分布を求めているが，GOTRESSるものとして Trinity［１５］がある．Trinity は GS2［１６］とのでは温度と温度勾配を同時に定めている．GOTRESS は通結合を目的として開発された．TGYRO と異なり非定常輸

常の空間離散化に伴う隣接グリッド間での変数の関連がな

送方程式を解いて時間緩和させ定常解を求めている点で通いため，径方向グリッドの数を極端に減らしてもその点に常の輸送コードと同じであるが，数値的安定に硬い輸送モ

おける解の精度に影響が出ず，また不等間隔に取っても問デルを取り扱うためにニュートン法 を 用 い て い る点で題がないなどの長所も有する．GOTRESS に実装されていTGYRO と類似している．GS2 の非線形計算から熱流束をる遺伝的アルゴリズムは並列化されているため計算を高速評価しているため，温度の定常分布を得るために，１万コ化しているほか，他のコードとファイル I/O を用いず MPI

ア前後を用いて数時間から１日程度を要している［１５］．類を介して情報の授受が可能になっている．TGLF のような似のコンセプトに基づくコードに TRESS+GKV［１７］があ並列化された輸送モデルとは，MPMD（Multiple Program

る．１次元輸送コード TRESS（Transport Equation Stable

Multiple Data）型で並列実行されている．GOTRESS は

Solver）に GKV［１８］が結合されており，GKV で計算した線

TOPICS とのベンチマークの結果，極めて良好な一致を示

形成長率から簡易的な飽和モデルで熱流束を評価し，非定している．また，従来型輸送コードでは困難であった，

常計算によって定常分布を見いだしている．TRESS は空

TGLF を用いた JT-60U 放電の再現シミュレーションにも４５５

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019

成功している．最近で は，平 衡・電 流 コ ー ド ACCOME

子軌道追跡によって記述されているが，高速のテイル成分と加熱コード OFMCを組み合わせた統合モデルGOTRESS＋だけの解析であったり，加熱機構毎に別々に解析されていの開発が進んでおり，JT-60SA シナリオにおいて TOPICSたり，径方向輸送が考慮に入れられていない等の問題点がとのベンチマークが成功裏に行われたほか，シナリオ開発

あった．これらを解決するために，プラズマを流体ではなにも適用が進みつつある［２０］．

く速度分布関数によって記述する運動論的輸送ソルバを中核とする統合コードの開発が進められている．

５．
４ 機械学習を利用した高速化

TASK/FP は三次元軌道平均多成分フォッカープランク

機械学習モデルはデータが持つ特徴量を入力として用モジュールであり，速度空間二次元，径方向一次元の分布

い，対象としたデータに適した種類や連続値を出力して予関数の時間発展を記述することができる．バルク成分も含測するもので あ り，多 種 あ る 機 械 学 習 モ デ ル の 中でもまれているため，径方向拡散モデルを導入して空間輸送を

ANN は入力と出力の複雑な関係を再現できるため，非線記述し，クーロン衝突項も運動量だけでなくエネルギーも形な乱流輸送モデルを高速で模擬するモデルとしての利用保存するように拡張されている．解析例として，ITER プ

が進みつつある．ANN のうち，ディープニューラルネット

ラズマにおける電子，D，T，He の４成分の速度分布関数

ワーク（DNN）と呼ばれるモデルは，特に複雑な現象を再の時間発展を記述し，電子サイクロトロン波電流駆動，D

現することに適しているが，多くの変数を訓練用のデータの中性粒子ビーム加熱，T のイオンサイクロトロン高調波から最適化してモデルを構築する必要がある．近年では

加熱，アルファ粒子生成による燃焼立ち上げのシミュレー

Python で書かれた機械学習ライブラリが非常に充実してションが行われた［２５，２６］．その結果，分布関数の変形を考きているため，あらゆる分野に機械学習の応用が進みつつ慮に入れることにより，核融合反応率がより正確に計算である．前述のQuaLiKizやTGLFなど輸送モデルは高度化にきるようになり，高速イオンの径方向拡散係数のエネル

よって計算量の増大が問題になっていたが，これらの挙動

ギー依存性が加熱分布および核融合出力パワーに大きな影を 模 擬 す る DNN モ デ ル で あ る QuaLiKiz-NN［２１］や響を与えることが示された．

TGLF-NN［２２］が開発され，熱流束などの評価に掛かる計

しかしながら，燃焼立ち上げの現実的なシミュレーショ

算時間は１０５倍ほど高速化されており，実時間計算が可能

ンに向けては多くの課題が残されており．位相空間におけとなるレベルとなっている．DNN モデルはモデルの訓練

る径方向拡散の適切なモデル化，輻射損失やジュール加熱に用いたデータが適用対象のパラメータから大きく離れての導入化，準中性条件の確保等を解決する必要がある．

いた場合全く用をなさないこともあるが，それを克服する

５．
６ おわりにため，データ生成から DNN モデル訓練，分布予測までを一貫して行うことで，特定のケースに高度に適合した DNN

従来型の統合コードとは異なる近年の新たな統合輸送モ

モデルによる分布予測を行う試みも行われている［２３］．まデリングの動きとその進展について，実時間を超えた計算

た，局所ジャイロ 運 動 論 コ ー ド と 実 験 デ ー タ に 基づく

速度をめざす超高速化モデリングと，硬い輸送モデルによDNN モデルも開発されており，正味の流束を予測する他る困難を克服するためのモデリングに着目して概観した．

の DNN モデルと異なり，複数の輸送過程に起因する流束

多様な運転シナリオ開発のための試行錯誤やプラズマ放電を個別に示すことにも DNN モデルが用いられている［２４］．

実験の合間での次の放電計画の策定には，精度は粗くても高速な計算が必要となる．一方で，開発したシナリオは高５．
５ 運動論的統合モデリング

度な輸送モデルを用いた精度の高いシミュレーションによ

統合コードの中核としてプラズマの密度，温度，電流，る検証が欠かせない．高速化と高精度化をめざしたモデリ

回転等の時間発展を記述する輸送ソルバには，これまでプングはどちらか一方で十分というものではなく，目的と検ラズマを流体として取り扱う輸送方程式が用いられてき証の段階に応じて使い分け，全体としては両方を組み合わた．すなわち，プラズマ各成分の速度分布関数はマクスせて用いるような運用がなされていくものと思われる．矛

ウェル分布で近似されていた．しかしながら，実際には電盾する高速化と高精度化の命題は同時に達成することは不磁波や NBI による加熱や核融合反応によって，高速粒子が可能であると考えられてきたが，それを可能とする技術と生成され，速度分布関数はマクスウェル分布から変形す

して DNN が注目を集めており，ますます適用が進んでいる．高速粒子の存在は，衝突を通してバルク粒子を加熱す

くのは間違いないであろう．ここで紹介したモデリングはるだけでなく，加熱の物理過程を通して加熱効率を変化さ

既に従来型統合モデルの核となっていた輸送コードを代替

せ，輸送現象や不安定性にも影響を与える．特に核燃焼プする形での運用も始まっており［１４，２０］，さらに複数の

ラズマの立ち上げ時は，外部加熱により生成された高速燃

DNN モデルを組み合わせた統合シミュレーションも進行

料イオンが核融合反応率に大きな影響を与える．したがっしている［２７］．従来型統合モデルの能力を拡張するようにて，核燃焼プラズマの挙動を正確に予測し，信頼できる制

今後も開発が進んでいくものと思われる．

御手法を確立するためには，速度分布関数の変形を考慮に参考文献

入れた輸送解析が必要となる．
従来の加熱モデリングにおいても，２．
１１節で述べたよう［１］F. Felici and O. Sauter, Plasma Phys. Control. Fusion 54,
025002 (2012).

に，速度分布関数の変形はフォッカープランク方程式や粒４５６Special Topic Article

5. Future Prospects of Integrated Code Development［２］P. Geelen et al., Plasma Phys. Control. Fusion 57, 125008(2015).

M. Honda et al.

［１６］M. Kotschenreuther et al., Comput. Phys. Commun. 88,
128 (1995).
１
６年 秋 季 大 会，
［１７］本 多 充，仲 田 資 季：日 本 物 理 学 会２０14aKA-7 (2016).
［１８］T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
［１９］M. Honda, Comput. Phys. Commun. 231, 94 (2018).
［２０］本多 充，成田絵美，林 伸彦：第３５回プラズマ・核融合学会年会 4Bp07 (2018).
［２１］J. Citrin et al., Nucl. Fusion 55, 092001 (2015).
［２２］O. Meneghini et al., Phys. Plasmas 21, 060702 (2014).
［２３］M. Honda and E. Narita, 2nd Int. Conf. Data Driven
Plasma Science (2019).
［２４］E. Narita et al., submitted to Nucl. Fusion (2019).
［２５］H. Nuga and A. Fukuyama, Plasma Fusion Res. 5, S2068(2010).
［２６］H. Nuga and A. Fukuyama, Prog. Nucl. Sci. and Tech. 2,
78 (2011).
［２７］O. Meneghini et al., Nucl. Fusion 57, 086034 (2017).

［３］F. Felici et al., Nucl. Fusion 58, 096006 (2018).
［４］J.F. Artaud et al., Nucl. Fusion 58, 105001 (2018).
［５］C. Bourdelle et al., Plasma Phys. Control. Fusion 58,
014036 (2016).
［６］G.M. Staebler et al., Phys. Plasmas 12, 102508 (2005).
［７］G.V. Pereverzev and G. Corrigan, Comput. Phys. Commun. 179, 579 (2008).
［８］S.C. Jardin et al., J. Comput. Phys. 227, 8769 (2008).
［９］J.M. Park et al., Comput. Phys. Commun. 214, 1 (2017).
［１０］J.M. Park et al., Phys. Plasmas 25, 012506 (2018).
［１１］J. Candy et al., Phys. Plasmas 16, 060704 (2009).
［１２］J. Candy and R. Waltz, J. Comput. Phys. 186, 545 (2003).
［１３］E. Belli and J. Candy, Plasma Phys. Control. Fusion 50,
095010 (2008).
［１４］O. Meneghini et al., Nucl. Fusion 55, 083008 (2015).
［１５］M. Barnes et al., Phys. Plasmas 17, 056109 (2010).

４５７

J. Plasma Fusion Res. Vol.95, No.9 (2019)４
５８

小特集 統合コードによる磁場閉じ込め核融合プラズマシミュレーションの現状と今後の展望

６．おわりに
6. Conclusion and Acknowledgements林

伸 彦，福 山

淳１），村 上 定 義１），横 山 雅 之２），藤 田 隆 明３），本 多
１）

充

１）

HAYASHI Nobuhiko, FUKUYAMA Atsushi , MURAKAMI Sadayoshi , YOKOYAMA Masayuki２）,
FUJITA Takaaki３）and HONDA Mitsuru
量子科学技術研究開発機構，１）京都大学大学院工学研究科，２）自然科学研究機構核融合科学研究所，
３）

名古屋大学大学院工学研究科（原稿受付：２
０
１
９年５月２
０日）

統合コードの開発は，複数のコードの連携以上に，各々19360415 ，20226017 ，22760665 ，23360416 ，25420895 ，のコードを開発している研究者・技術者の緊密な連携，す25820442，２６２８９３５６，26420862．また，核融合科学研究所

なわち「人の統合」が必要不可欠である．本小特集の各章

一般共同研究 NIFS11KNTT008，NIFS14KNTT025，NIFSで紹介されたモデル・コードの開発や開発したコードでシ

16KNST092，数値実験炉研究プロジェクト予算 NIFS11

ミュレーションを行って結果を得るために，様々な分野のUNTT006，NIFS17UNTT008 の支援も受けた．

専門家や学生の活動，シミュレーションを行うためのスパ

統合コードによるシミュレーション結果の一部は，核融コン等の計算機の利用があり，さらに，国内における核燃

合科学研究所のプラズマシミュレータ，及び，国際核融合

焼プラズマ統合コード構想（BPSI: Burning Plasma Simula-エ ネ ル ギ ー 研 究 セ ン タ ー の 計 算 機 セ ン タ ー（IFERC-

tion Initiative）［１］を中心に，各々の活動を束ねて組織の垣CSC）の HELIOS 等を使用して得られたものである．第一

根を超えて協力して統合コード開発を推進してきた．これ

原理コードを含む統合コードのシミュレーションや，第一らのことを継続・発展させることが第１章で述べた原型炉

原理コードの結果に基づいた統合コードに導入できるモデに向けた統合コード開発に求められる．本章では，小特集

ルの開発には，スパコンは必須のインフラとなっており，の第１章から５章までの各々に貢献して来た方々に謝意を

統合コード開発の進展に合わせてスパコンの導入整備も欠表すとともに，統合コード開発の今後のさらなる発展を期

かせない．

待し，まとめとする．

第１章で紹介した原型炉に向けた統合コード開発計画

ま ず，国 内 の 主 要 な 統 合 コ ー ド TOPICS，TASK，［２］に関して，有意義な議論を行っていただいた原型炉設

TOTAL の開発に貢献していただいた国内外研究機関の多計合同特別チームの理論・シミュレーションワーキンググ

くの関係者の，BPSI 活動や共同研究を通じた長きにわた

ループのメンバーに謝意を表します．長期の開発計画を着る 多 大 な ご 貢 献 と ご 支 援 に 感 謝 申 し 上 げ ま す．特 に，実に進めるためには，上記のスパコン等のインフラだけでTOPICS に関しては，滝塚知典 招聘教授（阪大），小関隆久

なく，最初に述べた「人の統合」が重要であり，今後ます博士，清水勝宏 博士，矢木雅敏 博士，井手俊介 博士（量ます次代を担う若手研究者や学生の育成が求められる．

研）に，TOTAL に関しては，山"耕造 名古屋大学名誉教

本小特集が，若手研究者や学生にとって統合コードを授に，ご指導いただき開発してきました，ご貢献に感謝し使った研究の入門書となり，２０２０年に運転 を 開 始 す る

ます．

JT-60SA や２
０２５年に運転を開始する ITER に向けた研究が筆者一同，JT-60 チーム，及び LHD 実験グループの協力

広がっていくこと，原型炉実現に向けて研究が加速すること有意義な議論に謝意を表します．統合コードで用いるモとを期待する．

デルの検証には実験との比較検証が欠かせず，今後より緊

参考文献

密に実験と連携して検証活動を行うことが，統合コード開発の加速に欠かせないと考えます．

［１］Burning Plasma Simulation Initiative, http://bpsi.nucleng.
kyoto-u.ac.jp/bpsi/［２］原型炉設計合同特別チーム，理論・シミュレーション
ワ ー キ ン グ グ ル ー プ 報 告 書，http://id.nii.ac.jp/1657/
00074926/

本小特集で紹介された研究内容については，以下（番号順）の科学研究費補助金の支援により推進されたものである．15K06656，16K06948，17K07001, 18360448，18560790，

National Institutes for Quantum and Radiological Science and Technology, Naka, IBARAKI 311-0193, Japancorresponding author’s e-mail: hayashi.nobuhiko@qst.go.jp

４５８!2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

!!
はやし

のぶ

林

伸 彦

小特集執筆者紹介

ひこ

""
ほん

量子科学技術研究開発機構，先進プラズマモ
デ リ ン グ グ ル ー プ リ ー ダ ー・上 席 研 究員．１
９９
９年慶應義塾大学工学博士取得，日本原子力研究所を経て，２０
１
６年より現職．専門は様々な物理がひしめくダイバータプラズマのモデリングを始めて以来，いろんなモデリング研究一筋．難聴がひどくなって聞きづらいので反応がイマイチな時はもう一声下さい．眼鏡に補聴器を内蔵するか，安くて使い勝手が良い首掛け型の集音器をもう少し補聴器レベルに誰か改良してもらえないでしょうか．売れると思います．

だ

みつる

本 多

充

量子科学技術研究開発機構 那珂核融合研究所 先進プラズマ研究部 上席研究員．２００７
年京都大学大学院工学研究科・博士（工学）．
主にトカマク中の輸送現象や輸送シミュレー
ションを研究対象としている．最近は機械学習やベイズ推定にも興味を持って，研究分野への適用を進めている．家庭菜園とピアノの練習（弾けるとは言っていない）を趣味としています．ピアノ教室に通う息子に腕前で抜かされないよう独習に励む日々です．
あい

ふく

やま

あつし

福 山

淳

まつ
かみ

さだ

まさ

ほし

ゆき

たか

のぶ

の

かず

お

慶應義塾大学理工学部・准教授．２
００
８年に慶應義塾大学大学院理工学研究科・後期博士課
程を修了（博士（工学））．日本原子力研究開発機構，量子科学技術 研 究 開 発 機 構 を 経 て
２
０
１
８年４月より現職．周辺・ダイバータプラズマに加え，イオン源プラズマにも手を広げ，モデリング・シミュレーション
研究を継続中．１年ほど前から飼い始めた猫が，ようやく「お手」を覚えてくれました．

自然科学研究機構 核融合科学研究所 ヘリカ
ル研究部 核融合理論シミュレーション研究系
教授．工学博士（京都大学大学院原子核工学専攻）．LHD 実験解析型統合輸送解析スイート
TASK3D-a の運用で蓄積した熱輸送データベースを用いて，統計手法やデータ駆動手法との連携を実践，模索しています．数年のもがきの時代を越えて，科学研究費補助金にも採択され，少し風が吹いてきている印象があります．
たあき

星 野 一 生横 山 雅 之ふじやま

量子科学技術研究開発機構，プラズマ理論シ
ミュレーショングループ・主幹研究員．京都大学博士（エネルギー科学）．磁場閉じ込め核
融合プラズマ中の MHD・輸送現象や高エネル
ギー粒子の物理について理論やシミュレーションを用いた研究に携わっています．ここ数年はトカマク放電のディスラプション現象とそれに関わる物理と炉工学の関係が興味の対象です．

京都大学大学院工学研究科 教授，理学博士（広島大）．ヘリカル系の加熱シミュレーショ
ンからスタートし，新古典輸送や統合輸送シ
ミュレーションなどの研究を進めてきま した．現在は，非軸対称トカマクを含めた広い３次元系トーラス
プラズマにおける運動論的な問題などを中心に研究を行っている．ここ数年，健康維持のため毎週スイミングクラブに通っている．月１回の記録会で，前年の自身の記録を更新することを目標にしているが．
．
．
やま

ゆき

松 山 顕 之よし

村 上 定 義よこのぶ

量子科学技術研究開発機構 先進プラズマモデ
リンググループ 上席研究員．２
００
５年総合研究大学院大学 博士（学術）．専門はトカマクプ
ラズマにおける MHD 安定性の理論・シミュ
レーション研究．半年ほど前の解説記事で執筆者紹介を書いてから所属と娘の年齢が変わりました（２歳）．しゃべるようになった娘からの命令をこなしながら，新たな環境での仕事を早く覚えられるよう努力中です．

京都大学名誉教授．岡山大学工学部助教授，京都大学大学院工学研究科教授を経て２０
１
７年３
月定年退職．プラズマ中の波動現象と輸送現象の理論・モデリング・シミュレーションに基づいて，炉心プラズマ統合モデリングの研究に従事．退職後も炉心プラズマ統合コードTASKの整備，拡張，国内の利用者サポートに加えて，タイ，オーストラリア，中国等の利用者とも共同研究を進めている．昨年１２月，娘に初孫誕生．
むらば

相 羽 信 行

さ

たけ

しん

すけ

佐 竹 真 介あき自然科学研究機構 核融合科学研究所 ヘリカ
ル研究部 核融合理論シミュレーション研究系
准教授，２
０
０
３年 総合研究 大 学 院 大 学 博 士（学 術）．３次 元 磁 場 配 位 中 の 新 古 典 輸 送 現象，新古典粘性のシミュレーションが主な研究テーマ．最適化
配位の研究をやりたいとこの数年言い続けながら，現在メインで取り組んでいるのは不純物の新古典輸送現象の研 究．
１歳になる娘を溺愛しつつも，子育てと仕事の両立に頭を悩ます日々です．

藤 田 隆 明名古屋大学大学院工学研究科 教授，工学博士（東京大学）．日本原子力研究所（日本原子力研究開発機構）での JT-60U 実験研究を経て，２
０
１
３年４月 よ り 現 職．小 型 ト ー ラ ス
TOKASTAR-2装置による実験研究と統合コードTOTALによるシミュレーション研究を行っています．キャンパス内の高校に通っている藤井聡太七段の活躍に刺激され，大学時代以来久しぶりに趣味と言えるくらいに将棋（中継観戦やネット
対局など）を楽しんでいます．

４５９

Journal of Plasma and Fusion Research Vol.95, No.9 September 2019
や

ぎ

まさ

とし

やま

矢 木 雅 敏量子科学技術研究開発機構 六ケ所核融合研究
所 核融合炉システム研究開発部部長，工学博士（京都大学大学院原子核工学専攻），日本原子力研究所研究員，九州大学応用力学研究所
教授，日本原子力研究開発機構プラズマ理論シミュレーショ
ングループリーダー，量子科学技術研究開発機構核融合炉システム研究開発部次長を経て２
０
１
９年４月より現職．九州大学
在任中に，福山さんと核燃焼プラズマ統合コード研究会を立ち上げる（今年で第１７回目を迎える）．これまで津々浦々を渡り歩く（岡山県，京都府，茨城県，福岡県，青森県）．
せき

りょう すけ

關

良 輔

ぐち

ひろ

ゆき

山 口 裕 之自然科学研究機構 核融合科学研究所 ヘリカル研究部 助
教．京都大学大学院工学研究科 博士（工学）．磁場閉じ込めプラズマにおける高速イオンの挙動とプラズマ中の輸送現象のシミュレーション・モデリングを中心に研究を行っている．
なり

た

え

み

成 田 絵 美量子科学技術研究開発機構 那珂核融合研究所
先進プラズマ研究部 研究員．２
０１
５年大阪大学
大学院・博士（工学）．トカマクプラズマにおける乱流輸送を主な研究対象として，実験解析やモデル開発を行っている．２
０
１
７年に機械学習の適用を始めたことをきっかけに，天体や材料など他分野の研究者と交流する機会を得て，刺激を受けている．趣味はグラフィックデ
ザイナー・絵本作家のディック・ブルーナ氏の作品鑑賞です．

自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 プラズマ加熱物理研究系 助教．工学博士（北海道大学大学院量子理工学専攻）ヘリカル系における NBI や ICRF の加熱シ
ミュレーションや，高速イオンの閉じ込めや高速イオン由来の不安定性のハイブリッドシミュレーションが主な研究テー
マ．２匹の猫を飼い始めました．まだ，子猫なので毎日走り
回っています．長期の出張にでると忘れられそうで心配です．

４６０