J. Plasma Fusion Res. Vol.88, No.9 (2012) 467‐474

講座 輻射流体シミュレーション

１．輻射流体シミュレーションの基礎
大西直文
東北大学大学院工学研究科航空宇宙工学専攻（原稿受付：２
０
１
２年５月２
５日）

輻射流体シミュレーションの基本的な構造と，解くべき方程式系の性質について概説する．また，現在利用されている代表的な輻射輸送計算手法について紹介し，それらの特徴と利点，欠点についていくつかの観点から比較を行う．現在の輻射流体シミュレーションにおいては，輻射モーメント方程式に基づく数値解法が最もよく輸送計算に用いられているため，その背景を説明した後に，古典的な流束制限拡散法に始まり非等方性を考慮できる手法までエディントン因子の評価法という側面からまとめて紹介する．
Keywords:
radiation transfer, laser-produced plasma, astrophysical flow, flux-limited diffusion, Eddington factor

１．
１ はじめに体中の原子過程という３つの要素で構成されるが，本講座

高温プラズマでは，電磁波や中性子，ニュートリノなど，

ではそのうち後者２つに焦点を当て，それらの計算手法の輻射と呼ばれる比較的平均自由行程の長いエネルギーキャ

概要を紹介することで，これから輻射流体コード開発を試リアが，しばしばプラズマの流体力学的挙動を大きく左右

みる読者の一助としたい．

する．特に，エネルギー密度の高い爆発的現象において輻

輻射輸送方程式の解法や流体計算との結合方法は，プラ射の寄与が顕著となることが多く，それゆえ衝撃波の物理

ズマの光学的厚さや輻射強度に応じて適切に選択する必要とともに研究されてきた［１］．そのような輻射流体力学的がある．次節に続く本章では，代表的な輻射輸送計算手法

現象は，宇宙物理においては古くから研究されており，実を挙げながらどのような手法を選択すべきかを述べる．典験室においてもレーザー生成プラズマの基本的挙動を決定

型的な輻射流体シミュレーションでは，輻射輸送計算に全する．例えば重力崩壊型の超新星爆発では，原始中性子星

計算時間の９割以上を費やすため，なるべく計算負荷の低から放射されるニュートリノがバウンス時に生成される衝い手法を選択することが求められる．また，レーザープラ

撃波を加熱し，爆発現象に大きく寄与すると考えられていズマのように輻射場が等方に近い場合に用いられる輻射

る［２］．また，慣性核融合ではそもそも輻射によって爆縮モーメント方程式の解法について詳説し，そこに非等方効を駆動する方式もあるため，現象の予測に輻射流体力学の果を取り入れることのできる可変エディントン因子法を紹

理解が不可欠である［３］．

介する．第２章では，超新星爆発におけるニュートリノ輸

ところが，そこでは連続体と飛程の異なる粒子（光子を送を例に，非等方性の強い輻射場における輸送計算手法に含む）が混在した輸送現象や，さらにその粒子の生成，消ついて紹介する．角度方向の情報を保持する輸送計算手法

失過程が複雑に絡んでいる．光と流体という異なる時間・は一般に計算負荷が高くなるが，大型計算機の発展に伴

空間スケールの現象を含む物理の解析的な取り扱いは難し

い，そのような手法を用いた輻射流体計算が実現可能にく，数値シミュレーションの担う役割が極めて重要になっなってきている．代表的手法である !! 法について，ニューてくる．現在，世界的に見ると，輻射流体コードと呼ばれトリノ輸送計算への応用を紹介する．最後に第３章では，るものは宇宙物理と慣性核融合のコミュニティを中心とし輸送計算に必要となる放射係数，吸収係数がどのように算て多数存在する．しかし，輻射流体力学は，流体力学，熱出されるかについて述べる．流体の輻射特性は，それを構

力学，統計力学，運動論，原子核物理など，専門性の異な成する原子が輻射の放出・吸収によって状態遷移すること

るいくつかの要素から構成されるため，コード開発には分で決まる．この過程を原子過程と呼び，適切なモデル化が野を横断する複数の研究者によるプロジェクト的な枠組み必要となる．代表的なモデル化とその計算手法，および輻

が必要である．残念ながらそのような組織的研究が不得意

射スペクトルの表現方法について解説する．

な我が国では，成熟した輻射流体コードが存在していると

１．
２ 輻射流体コードの構造は言い難く，したがって輻射流体コードに関する指南書も十分でないのが現状である．輻射流体コードは大きく分け

輻射流体力学の方程式は，輻射と相互作用する物質て流体，輻射輸送，およびそれらの相互作用を決定する流（以降，流体と表現する）を支配する流体方程式をベースに

1. Fundamentals of Radiation Hydrodynamics Simulation
OHNISHI Naofumi

author’s e-mail: ohnishi@rhd.mech.tohoku.ac.jp

４６７!2012 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012

し，輻射との相互作用項を付加するという形で表すことがまま輻射流体コードの構造を表す．輻射流体コードを構成できる．例えば流体が非粘性圧縮性である場合，輻射流体

する３つの要素に対しては，それぞれ独立に数値計算手法

方程式は以下のようになる［４］．

が発展している．多くの輻射流体コードでは，各要素に対応した支配方程式を別々の手法で解き，要素間で物理量の

(
'
'#(
!&"
'
#!
(
'

（１）

('
'#!!"
$#(
#!
'##!&!P(
!&"
'
(
'

（２）

している．

('
$!"(
P(
#!!*
$!&(
!&"
)
'
##"
'
&"
(
'

（３）

性質，特に解を得るにあたっての困難さの要因となるもの

やりとりをすることで結合計算を行う方式を採用している．図にはそのやりとりに用いられる代表的な物理量を記数値計算法の紹介に入る前に，もう少し支配方程式系のについて述べておく必要がある．まずはこの方程式系が極

ここで '，#，&，$はそれぞれ密度，流速，圧力，流体のめて多次元の問題であることに注意されたい．それは，後全エネルギーを表す．流体に対する付加項は，運動方程式

述するように輻射輸送方程式が輻射を担う粒子に対するボに輻射エネルギー流束 !の時間微分項と輻射圧力テンソル

ルツマン方程式であること，さらには，原子過程が多電子

P の発散項，エネルギー保存式には輻射エネルギー密度 "系の量子力学に支配されることに由来する．したがって，の時間微分項と輻射エネルギー流束の発散項，および輻射

輻射流体力学の方程式系を「真面目に」解くことは断念せ圧による仕事の項 #"P(が入る．ただし，流れ場が相対'&"ざるを得ない．さらに方程式系には硬直性がある．例えば論的である場合，相互作用を流体の共動系で評価する必要

多くの場合，流体の伝搬速度は音速であり，輻射の伝搬速があるため，方程式の形が若干異なってくる．また，宇宙度は光速であるため，それらの違いが硬直性を招くことは流体現象ではしばしば流体の圧力に比べて輻射圧が卓越す

容易に想像できる．このため，収束解を得るには多大な演るが，レーザー生成プラズマでは無視できる場合が多いた

算量を費やすことになる．

め，式（２）の輻射に関連する項，および式（３）の右辺は省以上のことから，輻射流体シミュレーションには如何に略されることが多い．輻射流体力学方程式の詳しい導出に「サボる」かという観点が必要となる．実のところ，輻射流

ついては，参考文献を参照されたい［４‐６］．

体力学の解を正確に知ることなど，ほとんどの場合不可能

輻射に関連する量 "，!，P は，周波数 %の輻射強度なのである．つまり，輻射流体力学の妥当な解を得ること（specific intensity）#%に対する伝播方向 Ω のモーメントとは，「適当な」近似をして演算量を減らし，必要な計算精して定義される．

度と使用できる計算資源のトレードオフの中でベストな解%%
Ω%%#! "%%%"$! ! "#
%%
! $&$
!
%

%

Ω%%#! !%%%!$! ! #
%Ω%
!

$&

%

!

%

P $! ! "#Ω%%#! P %%%
%ΩΩ%
! $&$
!

を探す作業になる．したがって輻射流体シミュレーション（４）の研究者には，様々な近似解法とその適用範囲を知り，場合によって解法を使い分けることが求められる．次節では（５）

輻射輸送方程式について概要を説明した後，代表的な近似解法をいくつか紹介する．

（６）

１．
３ 輻射輸送方程式

ここで，$，"%，!%，P %はそれぞれ光速，周波数 %におけ以降本章では輻射は光子（電磁波）であるとし，また，相る輻射エネルギー密度，輻射エネルギー流束，輻射圧力テ対論的な効果は無視する．輻射輸送方程式は，光子の分布

ンソルである．そして位置 "と方向ベクトル Ω の関数であ関数 # についてのボルツマン方程式から導出でき，輻射強

Ω!
る#
"
!
'
(は，後述する輻射輸送方程式によって支
%##
%'

度#
%%
# について以下のように記述される．
%#$

配される．ただし，輻射流体力学と呼ばれる現象においては，流体と輻射の方程式が独立に成立するのではなく，流体中の原子過程に従って輻射が流体から発せられ，また同時に流体に吸収，散乱されるため，それらは相互に依存している．さらに，流体の内部エネルギーは電子とイオンの熱運動だけでは決まらず，特に部分電離プラズマではイオ
ン化ポテンシャルの部分が比熱 !( の大部分を占めることになるため，原子過程を適切に考慮することで状態方程式という形でこの効果を流体へ取り入れる必要がある．すな
わち輻射流体力学は複合系の力学であり，方程式系は高度に非線形となる．
図１

長い距離を伝搬する輻射，中距離の輸送を担う流体，そして微視的スケールの原子過程と，輻射流体力学は多階層構造の問題とみることもできる．このような輻射流体力学の現象論的構造を示したのが，図１である．この図はその４６８

輻射流体コードの構造．輻射流体コードは大きく３つの要
素に分けられ，それぞれの要素間で物理量のやり取りが行われる．図中の括弧で囲った物理量は，レーザープラズマ
では通常考慮されないが，状況によっては（慣性系の変換や輻射による原子の励起および電離）必要になる情報を表す．

Lecture Note

1. Fundamentals of Radiation Hydrodynamics Simulation

N. Ohnishiが必要ないことがわかる．この事実が，後に述べるモーメ

#
!,
'
"Ω#
'#
'$&
'!#
'#
'
%,
*

ント方程式を用いた輻射流体方程式の解法に繋がる．
輻射輸送方程式の解は，方程式中の相互作用項がどの程

Ω!
Ω!
Ω ,)
Ω ,!
Ω)
Ω ,!
"
!
*
%(
*
%(
#Ω ,
)
+
)
!#
"
#
!
"*'! *
'(
'(
"(
（７）

度，あるいはどのように寄与するかに大きく依存する．本

ここで，&
'，#
'はそれぞれ放射係数，吸収係数を表す．まを基にした解の性質の分類について述べる．また，最後にΩ)は方向 Ω ,から方向 Ω への
Ω ,!
た，*'は散乱係数で，%(輻射輸送方程式の近似解法を数例紹介し，その特徴につい節の残りでは，いくつかの重要となる概念を挙げ，それら

散乱カーネルを表す．#'は定義より明らかなようにエネルてまとめる．

ギー流束の単位を持つが，独立変数として位置３次元と方１．
３．
１ 光学厚さ向２次元を持つため，ある点におけるある方向のみの，つ輻射輸送方程式の解の性質を分類するために最もよく使まり単位立体角中のエネルギー流束を表しており，ベクト

われる指標が光学厚さである．光学厚さ +'は以下のようにル量としてのエネルギー流束 !'とは異なる．

定義される．

光子はその性質から伝搬速度が変化しない．したがっ#'"*')+#)('$!

て，式（７）の左辺には通常のボルツマン方程式にみられる（１４）

ような加速度による速度空間の移流項がない．右辺は流体との相互作用項であり，第１項と第２項はそれぞれ輻射のここで，)は輻射の伝搬距離を表す．ほぼ同じ概念に光学

放出，吸収を表す．これらの項によって，流体から放出さ深さがあるが，光学深さは特定の区間に対して使われるのれたエネルギーは輻射に，輻射が吸収された場合は流体へ

に対し，光学厚さは系全体の値として語られる場合が多とエネルギーが移動する．式（７）の右辺第３項は散乱項でい．ここでは，特定の区間について議論するわけではないある．ここでは，コンプトン散乱等のエネルギー変化を生ので，光学厚さという言葉を使う．

じる散乱は考えていない．そのため，輻射は流体によってさて，この光学厚さであるが，吸収係数や散乱係数がそ散乱するが，流体へのエネルギーフィードバックはない．

れぞれに対応する平均自由行程 &$!
"'$の逆数で与えられ散乱によって光子の伝搬方向は変化するため，積分中の第ることから，積分区間を通して光子がどれだけ流体と衝突

１項 は 別 の 方 向 か ら 考 え て い る 方 向 に 入 っ て く る

するかを表す．ここで，'は衝突相手の数密度，$は衝突断（scattered-in）光子，第２項は今の方向から別の方向に移面積である．+'&!であれば光子にとって流体場は不透明

動する（scattered-out）光子の寄与を表している．

（opaque）であり光 学 的 に 厚 い（optically thick），+'%!

ここで後の議論のために，輻射輸送方程式に対し方向ベであれば流体場は透明（transparent）で光学的に薄い（op-クトル Ω のモーメントを取ってみる．簡単のため散乱の寄tically thin）という．

与を無視すると，０次および１次のモーメント方程式はそ

まず，光学的に厚い場合を考えよう．光学的に厚い場合，れぞれ以下のように表される．

光子は流体中ですぐに吸収あるいは散乱されてしまう．とこ ろ で，局 所 熱 平 衡（local thermodynamic equilibrium:

,
"'
!'$"
(&
#'"'
"'#
'!%
,
*

（８）

!,!'
P'$!#'!'
"%'#
%,
*

（９）

LTE）状態にある系ではキルヒホッフの法則から放射係数が吸収係数に比例する．
&
'$#
'!'

これらの両辺をさらに周波数 'で積分すると，（１５）

ここで，!'はプランク関数で温度の関数である．そのため，吸収が顕著な場合は放射係数も高くなっており，多数

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

（１０）

!,!
P $!!#'!'#'
"%'#
%,
*

（１１）

の光子が放出されている．放出された多数の光子がその周辺で吸収されると，そこで無視できない加熱が起こり再び放出が起こる．結果として光子は吸収，放出を繰り返し，放出は等方的に起こるため，光子は実効的にランダムとなり，式（２），（３）に現れる輻射に関連した項は，これウォークする状況になる．散乱が卓越する場合も多数の散らの方程式を用いて書き換えることができる．

乱を繰り返すため，やはりランダムウォークする．した,()#))##"()"'#($!!#'!'#',*

がって，光学的に厚い場合，輻射輸送は拡散的な現象にな（１２）り，後述するように輻射輸送方程式は方向２次元を落として輻射エネルギー密度の拡散方程式へと帰着することがで,$P （１３）
#+
$"()
(&
#'"')
"
#'"##
"'#
*
(
$!
(
(
'#)
'!%
,
*

きる．
では，光学的に薄い場合はどうであろうか．光学的に薄例えば，式（１３）の右辺第２項を無視すると，流体のエネル

い場合，光子はほとんど流体と衝突することなく系をすりギー方程式のソース項がその時刻での &'，#'，"'の積分だ

抜けていく．したがって，光子がただ移流する自由流（freeけで表されることを意味している．つまり，式（８），（９）-stream）の状況となる．この場合，問題は非常に簡単にから "'や !'が求まるのであれば，輻射流体計算には #'なったように感じられる．確かに，まったく放出以外の相
４６９

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012

互作用がなければ，単に輻射冷却が起こるだけあり，輻射

ルは量子力学的な準位間の遷移過程によって決まるので，輸送方程式を解く必要すらなくなる．つまり，図１の上と電子の束縛状態間の遷移はいわゆる線スペクトルとなっ右の要素だけで系を表現することができる．ところが，流て，エネルギー空間における変動が激しい．エネルギー空体の挙動に無視できないぐらいの吸収が起こるとなると，間の粗視化を行わず，特定の光子エネルギーへのすべてのやはり何らかの形で輸送方程式を解く必要が生じる．光学輻射遷移過程を逐次評価する手法を line-by-line 法とよぶ的に薄い場合，流体の時間スケールでは伝搬速度の大きいが，この手法は輻射スペクトルの詳細を知りたいときには輻射はほぼ定常になるため，輸送方程式の時間微分項を省有効であるものの，膨大な計算時間を要するため輻射流体

略できる場合が多い．ただし，非等方性が重要となるため，計算には不向きである．
そこで，エネルギー空間をある幅を持った領域で分割方向の次元を落とすことはできない．
１．
３．
２ LTE と non-LTE

し，その領域ごとに輻射輸送方程式を解くといった方針が図１に示したように，輻射流体計算において，流体と輻取られる．このようにエネルギー空間を有限のエネルギー

射輸送の橋渡しをするのが原子過程で決まる物質の輻射特群（bin ともいう）に分割することを多群近似，あるいはマ性である．物質が局所的に熱平衡にあるとする LTE の場

ルチバンドモデルとよぶ［８］．分割数が多いほど実際のス

合，物質の熱力学的状態は密度と温度に従い，その状態にペクトルに近づけることができるが，それに比例して計算

よって輻射特性も決定づけられる．これは，プラズマが比量も増大するため，不必要に分割すべきではない．また，

較的低温で高密度のときに成立する．ところが，レーザー

輻射流体シミュレーションは，通常，流体力学的な挙動にプラズマのように高温で低密度な領域が含まれる場合や，注目している場合が多く，輻射スペクトルの詳細まで再現

物質の状態が瞬時に変化するような 状 況 で は 必 ずしもすることは求められない．実際，エネルギー空間の解像度

LTE が成立しない．このような状態を非局所熱平衡（nonが低くても流体の挙動は十分予測できるため，分割数は通-LTE）という．

常数十のオーダーで，多くても百前後であることが多い．

LTE でなければ non-LTE な の で，一 口 に non-LTE と

ただし，特定のスペクトルがエネルギー輸送に大きく貢献

いっても様々な状態が考えられる．最もLTEから遠い状況しており，かつそのスペクトル構造が輸送に大きく影響す下では輻射特性に寄与する束縛電子の量子力学的状態が統る場合はその限りではない．

計的分布に従わないため，時々刻々状態遷移過程を追跡す

ちなみに，まったくエネルギー空間の分割を行わず，平る必要がある．しかし，多くの場合，ある熱力学的条件下均化された放射係数や吸収係数を用いて輸送計算を行うこでの定常状態を仮定できるため，non-LTE といっても時間とを灰色（gray）近似とよぶ．特に，光学的に薄い場合に依存のないある種の平衡状態を表す．これを準定常状態

はプランク平均，光学的に厚い場合にはロッセランド平均（quasi-steady state: QSS）という．最もよく仮定されるQSSとよばれる平均化が施される［１］．

は衝突輻射平衡（collisional radiative equilibrium: CRE）あ１．
３．
４ 様々な解法
さて，上述のような輻射輸送方程式の性質を考えながるいは衝突輻射定常状態（collisional radiative steady state:ら，いくつかの代表的な解法について概観してみる．

CRSS）と呼ばれる状態で，電子衝突過程として電離・再結合および励起・脱励起が起こると同時に，輻射過程とし

まず，モンテカルロ法はおそらく最も古くから輻射輸送ては再結合・脱励起だけが起こると仮定し，それらが全体方程式に用いられている解法の一つであろう．モンテカル

としてつり合っている状態である［７］．比較的高温である

ロ法では，サンプル粒子となる光子をランダムに発生さ程度密度が高いときにこの仮定は成り立つ．一方，アブせ，それに確率論的な衝突イベントを繰り返し与えること

レーションプラズマのように高温かつ低密度な場合は衝突で，輻射場の発展を追う［９］．光子発生の試行回数が多い

過程がほとんど起きないので，基底状態からの電子衝突とほど正確な解を得ることができるが，逆に少ないと統計的輻射遷移のみによって QSS になっており，これをコロナ平な誤差が生じるため，ノイズの多い解となる．モンテカル

衡とよぶ．

ロ法は，多次元問題の解の概形を得るには適しているが，流体計算や輻射輸送計算に使用する放射係数や吸収係数

高精度の解を得るには膨大な数のサンプル粒子が必要となは，密度や温度といった熱力学的諸量のみによって決定さる．また，光学的に厚い場合，衝突が頻繁に起こるため，自れなければ，膨大な数の状態遷移を記述する大規模な連立

由行程を衝突区間にとると光子がほとんど空間を移動でき

レート方程式をインラインで解く必要がある．そのため，ないという問題も生じる．そのため，どちらかというと光

通常は上記のような QSS を仮定し，輻射流体計算に必要な学的に薄い場合に適した解法である．また，散乱が容易に物理量を密度と温度の関数として予めテーブル化してお導入できるのも特徴の一つである．
非等方性の強い輻射場，つまり光学的に薄い系の解法とく．それらの具体的な計算手法については，本講座の第３章を参照されたい．

して，現在最もよく用いられるのは通称 "! 法とよばれる

１．
３．
３ 多群近似離散座標（discrete ordinate）法であろう［１０］．これは輻，つまり光子エネルギー #!の関数
輻射強度 !
!は周波数 !

射輸送方程式を伝搬方向も含めて離散化して解く手法であである．したがって !
!は光子エネルギー空間にも次元をる．ただし，当然ながら他の手法に比べて演算量が多い．

もっており，数値計算するためには，自然な発想としてこまた，方向次元の離散化が不十分であると，遠くの発光源の空間についても離散化することを考える．輻射スペクトが捉え難くなるいわゆる ray effect と呼ばれる問題が生じ
４７０

Lecture Note

1. Fundamentals of Radiation Hydrodynamics Simulation

N. Ohnishi

る［１１］．ただし，散乱も考慮することができ，モンテカル
ロ法のような統計ノイズもないため，計算資源さえ許せば

f %%広範な問題に適用することができる．この手法の応用例に

P%
$
!%

ついては，第２章で詳しく説明する．

%"
%ΩΩdΩ

（１６）

%"
dΩ
%

もしこの f %を何らかの方法で決めることができるなら，P %解の定常性を仮定することができれば，最も簡単な手法
は光線追跡法であろう．この手法では，解を必要とする場を!%で表すことができるので，式（８），（９）からなる方程所に向かって，外部境界から輸送方程式を光路に沿って積式系を閉じることができる．これは，流体の方程式におけ分するだけである．ただし，この手法でも #% 法と同様にる状態方程式の役割に似ており，モーメント方程式に基づ

ray effect が生じるので，なるべく多くの方向に対して積分

く数値解の信頼性はエディントン因子の決定方法に大きくを実行する必要がある．局所的な発光が重要になる場合，依存する．

それを遠くで捉えるためには膨大な光線数を必要とする．

エディントン因子の決定方法の詳細については後の議論

また，そのままでは散乱を扱うことが難しいため，定常かに譲ることとして，まずは f %が何らかの方法で得られた場

つ散乱の無視できる系で有効である．

合に，モーメント方程式を如何に離散化するかを考える．
P %＝f %!%を式（９）に代入すると，以下の式を得る．

最後に，本章の後半で詳しく説明するモーメント方程式に基づく解法は，光学的に厚い系で用いられる手法であ

!'!%
!
"$&#
'
$!#%!%%!%($'(

る．輻射輸送方程式の方向モーメントを取った方程式を扱うため，方向の次元がない分，演算量を大幅に削減できる．

（１７）

こ の と き，!%の 時 間 微 分 が 無 視 で き る と す る と，

ただし，非等方性を捉えることが苦手であり，それを克服するためには後述する付加的な手続きを必要とする．輻射

（８）から以下の
!
!%$!'
$!
#%(
&#
'
%!%(と表せるので，式と流体との相互作用が強い領域は光学的に厚い場合が多い

式を得る．

ため，輻射流体シミュレーションでは安定で演算量も少な

$
'
!%
!
&$
#%!%
!&#
'
! &#
"$$
%!%(
%!$
#%
'
(

いモーメント方程式に基づく手法が最もよく用いられる．
以上，これらの手法をまとめると表１のようになる．光（１８）

学厚さの程度，散乱の重要性に加え，何らかの解を得よう

この式は!%に対する拡散方程式になっており，一般的な拡とするならば，計算精度の重要度，計算資源の限界等を考

散方程式の解法を用いて数値解を得ることが可能となる．

慮した上で計算手法を選定する必要がある．

また，前節で説明したように光学的に厚い系では輻射場が拡散的になるため，上式はそのような状況下で成り立つ式

１．
４ モーメント方程式の解法

で あ る と 考 え ら れ る．特 に，光 学 的 に 十 分 厚 け れ ば，

ここから本章の残りを使って，モーメント方程式（８），

I（I は単位テンソル）となって，問題は最も簡単化
f %$'
!
!
#(

（９）の解法について解説する．前節で述べたように，輻射される．例えば，平板一次元であれば有限体積法とオイ

流体シミュレーションにおける輸送計算で最もよく用いら

ラー陰解法によって次のように離散化できるであろう．

れる手法である．モーメント方程式は元の輻射輸送方程式'' '"!
!'"!!%!
!&'"!!!&' $ %& &! #&"! &"!
'")"(#&"!!"")

よりも方向ベクトルの分だけ次元が少ないが，その代わり
未知変数が "
%から !%，!%，P %の３つに増えている．ところがモーメント方程式は２つであるから，このままでは解

' '"!
'
%
! !%
!'"!
'
' '"!
!& （１９）#&&$!$$$! & & ' &!! &!! $&#&!!!"")を得るための式が足りない．式の数を増やすためにさらに高次のモーメントを取ったとしても，また新しく未知変数が追加されるため，方程式系を閉じることができない．そここでは添字の %は省略し，等間隔格子を仮定した．簡単

こでいわゆるクロージャーとして，次のエディントン因子のために $，#，%
はすべて 'ステップでの値を用いたので，（テンソル）を導入する．

上式は !&'"!に関する連立一次方程式になっており，それを SOR や ICCG などの連立一次方程式の解法を用いて解く

表１

ことで !&の時間発展を得ることができる［１２］．ただし，

輻射輸送計算手法の比較．これらの評価はあくまで筆者の経験によるものであり，数値的裏付けはない．それぞれに提案されている改良法やプログラミング技法によっても相対的評価は異なってくるが，各々の手法の基本的なアイデ
アにそうような評価にした．

採用する解法によって解を得るまでの計算時間が大きく変わるため，適切な解法の選択が必要である．
式（１６）からわかるとおり，エディントン因子は輻射強度"%で定義されるため，定義どおりに求めようとすると "%moment eqs.
ray trace Monte Carlo#%computational cost
low
high
high
mediumrequired memory
small
large
large
mediumanisotopic field
poor
fine
fine
fine
opaque fieldfine
medium medium
poor
scattering field
N/A
available complicated available
time-dependence available available
N/Aavailableのために元の輻射輸送方程式を解かなければならない．したがって，何らかの仮定の元に近似的に決定すべきである．特に，流れ場の状態に従ってエディントン因子を変化させる手法を可変エディントン因子（variable Eddington
factor: VEF）法という．エディントン因子の決定方法は大きく分けて２つある．１つは局所的な物理量のみによって
４７１

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012

決める方法であり，もう１つは非局所的な影響も取り入れ

光学的に薄い場合であり，*!'*"'となっていると考え%'る方法である．光学的に厚い状況では，光子が遠くまで衝られる．したがって，式（１７）は実質的に移流方程式になる．

突なしに進むことができないため，輻射場は局所的な情報

このときの移流速度は光速で決まり，陽解法で安定に計算で決まり，前者の取り扱いが適用できる．局所的な取り扱するためには #+を小さく設定する必要がある．そのため

いでは，あまり多くの演算を必要としないため，効率的な離散化は陰的に行われるべきであり，例えば平板一次元で計算が期待できる．しかし，光学的に薄い場合は輻射場があれば以下のような離散式が考えられる［１６］．

遠くの光源によって決まる可能性があるため，後者のよう

""
*""
!#)*
*
* *"" #
#
#
"""
!""
)
"
(%
$
")*""#")*"#+!%
!'
!
"
)
) )
#,
（
２１）
#
#)*
'
#
+
*""
#
"""
#)"""##!
**$)$)#+ ""'
#+
""'
#
#
"""
"""

な取り扱いが必要となる．
１．
４．
１ エディントン因子の局所近似法
エディントン因子を局所的に決定する方法は，これまでいくつも試みられている．基本的には，光学的に厚い限界

*
* *""
(
"*""!("!)"" )"" ) ) "（２２）#,で拡散近似に，薄い限界では自由流になるように構築される．ただし，光学的に薄い場合にエディントン因子を局所的に決定可能な根拠は原理的に存在しないため，決定的な

ここで添字 'は省略し，")，#)"""#はスタッガード格子で定手法はない［１３］．したがって，光学的に薄い領域での解が義している．このような離散化であれ ば，結 果 と し てそれほど重要でないが，厚い領域と混在しているような場")*""の連立一次方程式に帰着させることができる．

合に，安定かつ効率的に両方の領域をカバーするための手法とみるべきであり，以下に挙げる様々な方法から実際に

ただし，%" 法では光学的に薄い限界で伝搬速度が '"$+に漸近することが指摘されており，そのような領域では 'どれを使用するかは，そのような観点で選択されるべきでに正しく漸近する %""$法が好ましいとされている［１７］．

ある．

最大エントロピー法

流束制限拡散法

ボーズ−アインシュタイン統計に従う光子がエントロ

流束制限拡散（flux-limited diffusion: FLD）法は，エディ

ピーを最大にするように分布関数を求め，そこから決まる

I としたときのントン因子の決定法というよりは f '#(""$)

エディントン因子を用いる方法が Minerbo の最大エントロ式（１８）の近似解法である．エディントン因子が一定の値をピー法である［１８］．輻射場の軸対称性を仮定すると，エ

取 る と き，式（１８）に お け る エ ネ ル ギ ー 流 束 はディントン因子は以下のように表される．

!'#!(
'"
$')
$
'"'で表されるため，光学的に薄くなって$'が小さくなったり，"'の勾配が急峻になると，流束が光

""
"$
""
(
!(
)
!")
f '# (
I" (##速から決まる限界値（free-stream limit）である '"'を超え

ただし，"#!'"!'*である．Minerbo の原論文では， (*る可能性がある．その場合，非物理的な流束となるため，

+#*
は#
!'"
'
"'*の関数として与えられており，そのようなこれを超えないように流束を制限する必要がある．

表現は後述の $" 法にも使用することができる．

制限方法はいくつか考えられるが，例えば以下のような

一方，一種のFLDとして使用されることも多く，式（２０）ものが挙げられる．
'!'#!
'"'#!!)*('"'
$'"*
"
"'
$
'"'*

（２３）

と 同 様 に !'が "'の 勾 配 で 表 現 さ れ る と す る と， (は"$'"')をパラメータとして
&#*
(
'"'*

（２０）

(
#&"&#&#

ここで，!)*( は流束制限された拡散係数である．一見，こ（２４）

れはエディントン因子を操作しているようにもみえるが，と表される［１９］．このとき，"#'"'"*'"'*である．また，それとは若干思想が異なる．本来，一次元で考えた場合，&は流束制限関数であり，最大エントロピー法では以下のエディントン因子は光学的に厚い領域から薄い領域に向けようになる．
!
"&#) if !$&$$#$"#'""#
(
"
#
&(
&)
#$
"
!"
&
"
"
)
(
"&
"
"#
+
if
"
!&!&
$
#
$
#

て 1/3 か ら１に 変 化 す る よ うに 構 築 さ れ る べ き で あ る［１４］．拡散近似では輻射場はあくまで等方なままなので，いくら流束を制限したところで，非等方な自由流を表現す

Levermore-Pomraning 法と M1 法

ることは難しい．また，上記の表現では光学的に厚い領域
でも０次精度しか達成できないという指摘もあり，他の流輻射場の角度依存性を示す規格化した輻射強度を定義

束制限法が提案されている［１５］．

P 1 法と P 1/3 法%" および %""$法では，エディントン因子を一定とする
が，代わりにエネルギー流束の時間発展方程式も解くこと
で，エネルギー流束が流束限界を超えないようにする．例
I，%""$法では f '#I とする［１３］．
えば，%" 法では f '#(""$)
エネルギー流束は式（１７）の時間微分を残したまま解くことによって求めるが，時間微分項が重要になってくるのは（２５）

し，それが時間的，空間的に十分緩やかに変化するという仮定に基づいてエディントン因子を決定することもできる［２０］．この Levermore-Pomraning（L-P）法も FLD として使用されることが多く，&#*"$'"')をパラメータと('"'*すると，以下のような &の近似式を使って式（２３），（２４）からエディントン因子を求めることができる［１４］．
&(
&)
#
４７２

#"&
&"&#
&"$

（２６）

Lecture Note

1. Fundamentals of Radiation Hydrodynamics Simulation

N. Ohnishi

リッドと捉えることができる．
最大エントロピー法と L-P 法で評価した )を図２にそれ

光線追跡法によるエディントン因子の評価

ぞれ示した．'は光学厚さの程度を示しており，光学的にエディントン因子は流体の時間スケールにおいて定常で厚いとき 'は小さくなり， )は 1/3 に近づく．光学的に薄あると仮定すると，エディントン因子を評価するための輸い場合は 'が大きくなり，)が１に漸近することがわかる．

送方程式の時間微分項を無視できる．そこで，以下の定常

ここで，)は式（２３）中の変数であるが，一次元のときにはな輸送方程式を光線追跡法で解くことを考える．

これがそのままエディントン因子になる．

&#
Ω#
%#
$$
$$#
$!"
$#
$
&*

Levermore はこれとは別に，現在 $! 法のクロージャー
としてよく用いられるエディントン因子の決定法も提案し（２８）ている［１９］．$! 法は主に宇宙物理における輻射流体計算

ただし，*は光路に沿った距離を表し，簡単のため散乱項で採用されており，式（８），（１７）を !$と !$から求まるエは省略した．光線追跡のアルゴリズムは単純であり，#$

デ ィ ン ト ン 因 子 を 用 い て 解 く 手 法 で あ る．こ の 場 合，
'$&
!$!
"
(
!$&を用いて，は場所と方向の関数であるから，エディントン因子を得た)$

い場所 &からある方向に遡る光路と外部境界との交点 %を求め，点 % から点 &に向かって上式を積分すれば点 &

'
"
#"$

"

'
"
%""!$!#

"

での #
$や "
$はセルで定義される値$が求まる．このとき，#（２７）なので，通過するセルをどれだけ横切るかを逐次評価する必要がある．

という式からエディントン因子が求められ，光学的に薄い

上 記 の よ う な 計 算 は long characteristics 法 と よ ば れ

状況下の検証も行われている［２１］．
１．
４．
２ 非局所的に決定するエディントン因子［２２］，手続きとしては単純であるが，エディントン因子が局所的にエディントン因子を決定する場合，遠くの光源

必要なすべての格子点で，ある程度の方向離散点の数だけの影響をうまく取り入れることができない．モーメント方

光線追跡を行わなければならず，演算量はたちまち膨大な

程式に基づく解法は本来そのような状況に適していないものとなる．そこで，図３に示すように，光線追跡は隣りが，これを克服するために非局所的なエディントン因子ののセル境界からのみ行い，初期値となる境界での輻射強度評価が試みられている．

は周辺の値を補間することによって評価する short characteristics 法が提案されている
［１６，
２３］．

式（１６）の定義に立ち戻ると，エディントン因子を評価するためには #$が求まっていればよい．ところが，そもそもshort characteristics 法によって演算量を抑えることは#$は輻射輸送方程式を解かなければ決められないので，こできるが，光源から遠くにいくほど非等方性が失われてしれを用いてエディントン因子を定義するのは本末転倒であまう問題がある．一方 long characteristics 法では，演算量ろう．しかし，計算負荷を小さく抑えながらある程度の輸を厭わなければ精度よくエディントン因子を得ることがで

送方程式の解が得られるのならば，その解からエディント

きるが，有限の方向離散点ではいずれ遠くの光源を捉えらン因子だけを評価し，時間発展はモーメント方程式を用いれなくなるという光線追跡法の根本的な問題は避けることて評価するという手続きが可能となる．これは，エディン

ができない
［２４］．

トン因子を評価するための輻射輸送計算にはそれほど精度を求めず，精度が重要となるエネルギー輸送の評価に関してはモーメント方程式の解法に委ねるという発想から導かれており，連続体的なアプローチと運動論的手法のハイブ

図３

図２

エディントン因子の近似法の比較．実線は最大エントロ
ピー法で破線は L-P 法．一次元のときのエディントン因子となる )は光学厚さを表すパラメータ R = | %E | / ( "E ) に対して 1/3 から 1 まで変化する．

４７３

光線追跡法によるエディントン因子評価法の概念図．long
characteristics 法では，ある方向の輻射強度が点 P で必要な場合，そこから光路を遡って外部境界との交点 O を求める．そして，点 O から P に向かって，輸送方程式を積分する．一方，short characteristics 法で点 P´の輻射強度を評価する場合，隣のセル境界との交点 O´から P´の間のみ積分する．積分の初期値となる点 O ’における輻射強度は，周辺の輻射強度から補間して求める．

Journal of Plasma and Fusion Research Vol.88, No.9 September 2012［７］D. Salzmann, Atomic Physics in Hot Plasmas (Oxford University Press, Oxford, 1998).
［８］D.E. Cullen and G.C. Pomraning, J. Quant. Spectrosc. Radiat. Transfer 24, 97 (1980).
［９］Y. Wu, M.F. Modest, and D.C. Haworth, J. Comput. Phys.
223, 898 (2007).
［１０］小林啓祐：原子炉物理（コロナ社，東京，１
９
９
６）
．
［１
１］P.J. Coelho, J. Quant. Spectrosc. Radiat. Transfer 73, 231(2002).
［１
２］小国力 編著，村田健郎，三好敏郎，J.J.ドンガラ，長谷川秀彦著：行列計算ソフトウェア−WS，スーパーコ
ン，並列計算機（丸善，東京，１
９
９
１）
．
［１
３］G.L. Olson, L.H. Auer and M.L. Hall, J. Quant. Spectrosc.
Radiat. Transfer 64, 619 (2000).
［１４］N.J. Turner and J.M. Stone, Astrophys. J. Suppl. Ser. 135,
95 (2001).
［１５］J.E. Morel, J. Quant. Spectrosc. Radiat. Transfer 65, 769(2000).
［１６］J.C. Hayes and M.L. Norman, Astrophys. J. Suppl. Ser. 147,
197 (2003).

モンテカルロ法とのハイブリッド計算
連続体に対する方程式と運動論的方程式のハイブリッド
解法を考えるとき，計算格子ベースの手法と粒子的手法の結合計算が最も現実的な選択である．実際，光線追跡法も一種の粒子的手法であり，そのような発想の範疇を脱しない．ところが，光線追跡法は粒子の飛行経路を固定した計算手法であるため，予め用意した飛行経路にない情報を取り入れることができないという意味では，計算格子ベースの手法と解釈することもできる．すなわち，計算空間中の点光源を隈なく捉えるためには網羅的に飛行経路，つまり離散方向を設定する必要があり，それは計算資源の面から考えて到底現実的ではない．
そこで，もう一つの粒子的手法としてモンテカルロ法を考えることは有益であろう［２５］．前節で述べたように，モンテカルロ法は解の統計誤差を除去するために多くのサンプル粒子を必要とするが，解の概形を得るだけであれば，比較的少ない計算負荷で目的を達成することができる．エディントン因子を評価するだけであれば，多くの計算負荷［１７］K.H. Simmons and D.Mihalas, J. Quant. Spectrosc. Radiat.
Transfer 66, 263 (2000).
［１８］G.N. Minerbo, J. Quant. Spectrosc. Radiat. Transfer 20,
541 (1978).
［１９］C.D. Levermore, J. Quant. Spectrosc. Radiat. Transfer 31,
149 (1984).
［２０］C.D. Levermore and G.C. Pomraning, Astrophys. J. 248,
321 (1981).
［２１］D. Aubert and R. Teyssier, Mon. Not. R. Astron. Soc. 387,
295 (2008).
［２２］梅村雅之，中本泰史：ながれ 15, 457 (1996).
［２３］J.M. Stone, D. Mihalas and M.L. Norman, Astrophys. J.
Suppl. Ser. 80, 819 (1992).
［２４］N. Ohnishi, K. Sugai and Y. Ogino, J. Phys. Conf. Ser. 244,
022078 (2010).
［２５］N. Ohnishi, Toward an accurate numerical simulation ofradiation hydrodynamics in laser ablation plasmas,
Proc.9th International Conference on High Energy Density
Laboratory Astrophysics, 2012.

を必要としない可能性が高い．発生確率を放射係数で重み付けすることで，光源からより多くのサンプル粒子を生成し，点光源の寄与をなるべく全域に行き渡らせることも可能である［９，２５］．

参考文献［１］Y.B. Zel'dovich and Y.P. Raizer, Physics of Shock Waves andHigh-Temperature Hydrodynamic Phenomena (Dover, New
York, 2002).
［２］K. Kotake, K. Sato and K. Takahashi, Rev. Prog. Phys. 69,
971 (2006).
［３］S. Atzeni and J. Meyer-ter-Vehn, The Physics of Inertial Fusion: Beam Plasma Interaction, Hydrodynamics, Hot Dense
Matter (Oxford University Press, Oxford, 2004).
［４］D. Mihalas and B.W. Mihalas, Foundations of Radiation Hydrodynamics (Dover, New York, 1999).
［５］G.C. Pomraning, The Equations of Radiation Hydrodynamics(Pergamon Press, Oxford, 1973).
［６］J.I. Castor, Radiation Hydrodynamics (Cambridge University Press, Cambridge, 2004).

!!!!!!!!!!!!!!!!!

おお

にし

なお

ふみ

大 西 直 文東北大学大学院工学研究科航空宇宙工学専攻准教授．２０
０
１年３月大阪大学大学院工学研究科電気工学専攻にて学位を取得（博士（工学））．専門は輻射流体力学，高温気体!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

力学．出身は島根県で大学は大阪，就職は関東を飛び越して
仙台，と東にシフトし，おかげで大きな地震に２度遭いました．帰省が大変で，毎年盆には家族を乗せて仙台から島根を走破します．関西より西の道路で関東ナンバーを見ると，勝ったと思ってしまいます．

!!!!!!!!!!!!!!!!!!!!!!!!!!

４７４