Transcript PPT - 理論天文学宇宙物理学懇談会

page 1 / 51
多次元 輻射輸送・輻射流体力学シミュレーション
中本泰史（筑波大学）
1. 輻射とは
2. 輻射輸送
3. 輻射流体力学
4. 最近の計算例
5. 今後
「輻射輸送」は多様
「多次元」輻射輸送・輻射流体力学は
line continuum
注：本講演では主に，２・３次元空間内の
未熟・未発達！
1D continuumの輻射輸送・輻射流体力学計算に
2D ついて概観する．
１次元についてはlineも含め長い歴史と多量の
始まったばかり・黎明期にある！
3D ノウハウがあるが，それらは別の機会に．
参考書: 「Foundations of Radiation Hydrodynamics」 --- RHD，continuum
Mihalas & Mihalas, 1999, Dover
「Radiative Transfer in Moving Media」
--- 1D球対称，line
Sen & Wilson, 1998, Springer
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 2 / 51
1. 輻射とは
宇宙物理学における「輻射」の役割
1. 観測手段
2. 力（輻射圧，輻射粘性）
3. エネルギー輸送（加熱・冷却，電離・解離，化学反応）
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 3 / 51
「輻射」の記述


輻射強度 (Intensity)
I(x,y,z, ,,; t)
dE I dt dA d d


「光子」の分布関数との関係
dA
d
(,)
(x,y,z)

６次元位相空間上の関数


h 3 2
I(x,n,; t)  ch  3 f (x,n,; t)
c
注：n, νは，光子の運動量 p
の極座標表示に対応する．
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 4 / 51
輻射輸送方程式の導出 (1)
Boltzmann 輸送方程式
f
f
f f 
 v  F    
t
x
p t coll
輻射輸送方程式

1 I
 n I     I
c t
(x,n,;t)

: emissivity
(x,n,;t) : opacity（単位長さあたり）

第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 5 / 51
輻射輸送方程式の導出 (2)



dA
n
d

x
dl x  x
 
I(x  x,n,; t  t)  I(x,n,; t)dAdtdd
1 I(x,n,
 ; 
t) 
I(x,n,; t)
 

dl dAdtdd

c

t
l
 (x,n,; t)  (x,n,; t)I(x,n,; t) dl dAdtdd
1 I(x,n,; t) I(x,n,; t)
 (x,n,; t)  (x,n,; t)I(x,n,; t)

c
t
l
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 6 / 51
輻射輸送方程式の導出 (3)

と
I
n'
の
具体的な内容


散乱
散乱
I

吸収

n
放射

x

1 I
 n I  absB  sca  I(n')(n,n')d  abs  sca I

c t



1 I
 n I    I
c t
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 7 / 51
物質分布の対称性と位相空間の次元
方向１次元
Pz
Pz
方向２次元
θ
θ
φ Py
Px
空間２次元
空間１次元
I(y,z,,, )
I(z,,)


I(r,,)
r

I(R,,, )
I(R,Z,,, )
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
「輻射」が関与する問題の類型化：「何」が知りたいのか
ρ, T: given
Flux
ρ: given
T, E, Flux
ρ: given
page 8 / 51
I, ni, T, E, Flux
ρ, V, T, E, Flux
Radiation第15回
Hydrodynamics
(RHD)
理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 9 / 51
2. 輻射輸送（定常問題）
1 I
 n I    I
c t
定常性：
多くの現象では，

thydro > tradiation
V
hydro
である．
 c


定常性を仮定して良い．
n I    I
注：
輻射流体力学計算の
アルゴリズムでも，
定常輻射輸送方程式を
解くことは基本となる．
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 10 / 51
輻射輸送方程式の形式解（積分型）
l
dI
   I
dl
x
0

1 dI

 I 
 dl


dI
 I  S
d

Source Function
n, :
fix
I(x)  I(0)e 


x
  0  dl

 S(t)e
( t )
0
dt
光学的厚さ (深さ)：
空間距離 x を光子の平均自由行程 1
を単位にして測った距離
 1
光学的に厚い，不透明

 1
光学的に薄い，透明
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 11 / 51
数値解法
反復解法が必要


n I  (S  I)  absB(T)  sca  I(n')(n,n')d  abs  sca I
・Λ-iteration
・ALI (Accelerated Lambda Iteration)
・Variable Eddington Factor 法
・Feautrier法
・Rybicki法
・Complete Linearization 法
・...
・Monte Carlo法
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 12 / 51
差分化（一例）


IX  I0e  S0
0  X

2
1e  e

 1 e
 SX

Λ-iteration
(1) S, χを既知とする
(2) I を更新
(3) S, χを更新
(4) (1)に戻る
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 13 / 51
数学的構造

1
I(x)  I(0)e 
abs  sca


NX NY NZE
 
N 
N NI
n

 

B(T)   sca  I(n')(n,n')d' e( t )dt
abs
0
NX NY NZ  N N N
 bn1
Feautrier 法
(Rybicki法, Complete Linearization法, ...)
IX1,1,1  
IX1,1,1  bX1,1,1(T) 
   
 


 n1  n1
 


n
Λ-iteration法
IX1,n1, 1 
 b
I  I

IX1,n1, 1  bX1,n1, 1(T) 

 
 


I
I
b
(T)
 X 2,1,1  

 X 2,1,1   X 2,1,1
n  
n1
n1

 


E 
I



I

b
*
 *
 
 


IX 2,n 2, 2  
IX 2,n 2, 2  bX 2,n 2, 2 (T)
Accelerated
IX 3,1,1  
 Iteration

ILambda
bX 3,1,1(T) (ALI)
X 3,1,1
(Operator Perturbation
Method)

 
 




 
 






I  I b
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 14 / 51
計算量を減らす工夫
計算量
~ Niter  N  NX NY NZ  N N N

~ Nit er  NX NY NZ  N N N
・Short Characteristics法
・Accelerated Ray Tracing法
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 15 / 51
Monte Carlo法について
・Λ-iteration と同等
・位置・方向・振動数をランダムに振る
・点光源は扱いにくい
計算量
~ Nit er  NX NY NZ  N N N
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 16 / 51
Adaptive Ray Tracing
(Abel & Wandelt 2002)
角分解能を維持するため，
適宜 ray を分割
点光源から周囲への
輻射輸送の際(のみ)に有効
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 17 / 51
並列化
n I    I
右辺を既知とみなす反復解法では，
輻射強度は方向・振動数毎に独立

並列化は容易
p
p
・計算順序に注意
・MWF法で解決


方向・振動数分割

x
x
配位空間分割

第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史

page 18 / 51
3. 輻射流体力学
M


R
一般的な基礎方程式

;
F
1 I
 n I    I
c t

輻射エネルギー密度
(1 e /c 2 )
0 


p

M0  
0
p 



p


 E0 c1F0 
R0   1

c F0 P0 
Ex,,t 
1
I(x,n,; t) d

c
輻射エネルギーフラックス
Fx,,t   n I(x,n,; t) d
輻射圧テンソル

Px,,t 
1
nn I(x,n,; t) d

c
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 19 / 51
仮定
注：「輻射流体力学」の基礎方程式は扱う問題に応じていろいろ
あり得るが，以下ではその一例を紹介する．
・非相対論的 O(1) --- (v/c の 0次まで)
・任意の光学的厚さで成り立つ
・散乱を無視 (η0 は等方的)
注：以下に述べる方法でも，
v = 0 ならば非等方散乱も
扱うことが出来る．
・主にcontinuumを扱う
V  c*
V  c*
static diffusion
dynanic diffusion
free streaming
relativistic flow
 1



 1

c*  c /max(1, )
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 20 / 51
Comoving Frame
O(v/c)
z
1 I0 (0, 0 )
I ( , )
O(1) の式を得る
 (0   ) 0 0 0
V
c
t
z
 V 
 c*
 
c*
2 0 V Va
c*  c /max(1, )

 2 I0 (0, 0 )
1 0 
0 
c z c 

 1

0 V static

0adiffusion  dynanic diffusion



I (0, 0 )
 


 0 
 0  0 0
c z c 2 

2  1

relativistic flow
free streaming
 1 0 V  20a 
I
(

,

)

c 2  0 0 0
 c z
1 I
 n I    I

0 (0, 0 )   0 (0, 0 )I0 (0, 0 )
c t
すべての場合で無視できるもの以外は残す
Lab Frame


DE0 F0 V
2a


E0  P0   2 F0 
Dt
z z
c
4 ( )  c ( )E ( )d
 

0
0
0
1 DF0 P0 2 V
a
1


F

E

P




c 2 Dt z c 2 z 0 c 2 0 0
c
0


0
0
0
0
0
0 ( 0 )F0 ( 0 )d 0
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 21 / 51
輻射流体力学 基礎方程式： O(1)
注：O(v/c) の項を一旦は残した
上でO(1)の式を作らないと，
エネルギー保存則や作用・反作用
が満たされない式が出来てしまう．
見かけ上 v/c の項が残っているのは
重要な意味がある．
D
   v  0
Dt
Dv
1

 p    F F
Dt
c
D egas 
   p  v  4  c E E 
Dt   
３次元
D E 
    F  v :P  4  c E E 
Dt  
D F 
1
    P   F F
Dt  
c
６次元

1 I
 n I    I
c t
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 22 / 51
解法
1. Flux Limited Diffusion (FLD) 法
2. Simple Method
3. Variable Eddington Factor (VEF) 法
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
1. Flux Limited Diffusion (FLD) 法
Diffusion 近似
 1
  L
c
D D
  v
 0 v  0
Dt Dt
Dv Dv
1 1
 
  p p

E  F F
Dt Dt
3 c
 c

D  D eEgas
E
egas  p
p vv 4
 vcEE E 
 
Dt  Dt   
3
3

V



page 23 / 51
E  aT 4
1

1 

P  E 1 

3 

1

c
F   E
3
D E 
    F  v :P  4  c E E 
Dt  
D F 
1
    P   F F
Dt  
c

第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史


page 24 / 51
Flux Limiter:
F  
c
E
3
6  3R
6  3R  R2
E
R
E
(R) 
FLD
D
   v  0
Dt
Dv
1

 p   E
Dt
3
 c

D 
E 
E
 egas   p  v    v     E 
Dt 
 
3
 3


 c
 E R  0


c
3
F   E  
3
cE E R  

E

・３次元空間の計算
・エネルギー式は implicit積分
・F // ∇E の制約がある
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 25 / 51
2. Simple Method
D
   v  0
Dt
Dv
1

 p     F F
Dt
c
D egas 
   p  v  4  c E E 
Dt   

n I    I
E
1
I d

c
F   n I d

・３次元空間＋
６次元位相空間
・流体は explicit，
エネルギー式は，
implicit積分
・輻射場と反復計算
(整合性をとるため)
・アルゴリズムは
比較的簡単


第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 26 / 51
3. Variable Eddington Factor (VEF) 法
D
   v  0
Dt
Dv
1

 p    F F
Dt
c
D egas 
   p  v  4  c E E 
Dt   
n I    I
D E 
    F  v :P  4  c E E 
Dt  
D F 
1

    P   F F
Dt  
c
PfE
nn I d

f
 I d
Eddington
Factor
1
nn I d

c
1
E   I d
c
P

第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 27 / 51
Variable Eddington Factor 法の特徴
・空間３次元，位相空間６次元問題
・流体は explicit，エネルギー式は implicit 積分
・エネルギー式の積分は，やや複雑
・Eddington Factor f は，輻射場の等方・非等方性を表している
→ E, F そのものよりも，誤差が小さい
・ Eddington Factor f は，時間変化は大きくない
→ 輻射場との反復計算は少なくてすむ
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 28 / 51
各計算法の比較
RHD
HD
Polytrope
FLD
VEF
Simple
計算量
小
<
中
<<
大
<
~
大
複雑さ
小
<
中
<<
大
<
大
精度
低
<<
中
<
高
応用例
多
少々
高
極少数
注：FLDは空間部分だけの計算量で輻射が扱えるので，
計算量が少ない割には「輻射」という「質」の良い計算が
出来る．その意味でコストパフォーマンスが良い．
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 29 / 51
4. 最近の計算例
1. 1D軸対称 RHD
フィラメント状星間雲の重力収縮
2. 2D軸対称 輻射輸送
T Tauri型星・原始星の構造推定
3. 2D軸対称 RHD
大質量星形成
4. 3D 輻射輸送
宇宙の再電離過程
5. 3D RHD
銀河形成
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 30 / 51
1D RHD: フィラメント状星間雲の重力収縮
(Ogochi & TN, in prep.)
・フィラメント状分子雲 ＝ 無限に長いシリンダー
・磁場，回転，乱流，化学反応は無い
・星間輻射による加熱（境界条件）
・1D 軸対称 RHD (VEF)
加熱 = 冷却
=1
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 31 / 51
2D RHD：大質量形成
(Yorke & Sonnhalter 2002)
・30, 60, 120 Msun の分子雲
・2D RHD (FLD) で重力収縮計算
・Nested Grid (64x64, 3階層)
・FLDは波長依存性を考慮
輻射圧のために
降着が抑止される
形成された星の質量
初期 波長依存
灰色
30
60
120
19.1
20.1
22.9
31.6
33.6
42.9
初期 60 Msun Clump
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 32 / 51
2D 輻射平衡計算
T Tauri型星
(Kikuchi, TN, & Ogochi 2002)
T Tauri型星 HL Tau
・2次元 軸対称
・輻射平衡，VEF
Kikuchi, TN, & Ogochi 2002
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 33 / 51
ディスク・ハロー モデル
ハロー：エンベロープの
内側100 AU程度の領域
中心星の放射を
エンベロープが散乱
加熱されたディスク
からの赤外放射
中心星からの放射を直接吸収
するよりも得をする
エンベロープはディスクからの
赤外放射に対して光学的薄い
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 34 / 51
密度・温度分布
ディスク
中心星
ハロー
Kikuchi, TN, & Ogochi 2002
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 35 / 51
近赤外散乱光イメージ
観測 (HL Tau)
Close et al. (1997)
モデル計算 (i = 60o)
Kikuchi (in prep.)
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 36 / 51
原始星の構造推定
(Nakazato, TN, & Umemura 2003)
bp
 2-D axisymmetric
 three components;
Observer
i
Envelope
• Central star:
Luminosity
L*
Temperature T
*
Circumstellar disk
Mass
M*
• Circumstellar disk:
1
p 15
.
Surface density at 1AU
Power law index
• Envelope:
100AU
Density at 1AU
1
Power law index
q 15
.
Opening angle
bp
Central star
Temperature: Radiative equilibrium


0

 B d  0  J d
abs
abs
Outflow region
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 37 / 51
parameters
Model SED
νLν [erg s-1]
1036
観測角度依存性 大
i の推定
原始星構造の推定
1034
1  1013 g cm3
L*  10
. L
T*  4000K
M*  05
. M
i  0
30
60
90
1  2000 g cm2
p  q  3/ 2
bp  25
(νLν)max
1032
1030

LSED   L d log
0
1028
1011
1012
1013
 [Hz]
1014
1015
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 38 / 51
 IRAS04365+2535（TMC1A）
Class I 天体
νLν [erg s-1]
L* 1.0L
1 1012.5 g cm-3
1 10000 g cm-2
BP  20o
i
 22o
 [Hz]

(Nakazato, TN, Umemura 2003)
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 39 / 51
球対称エンベロープの場合
r-z平面上の温度分布
VEFとFLDは良い一致を示す
(Kikuchi, in prep)
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 40 / 51
Diskの温度分布
VEF
FLD (non-grey)
FLD (grey)
(Kikuchi, in prep)
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史

3D 輻射輸送
(TN, Umemura, & Susa 2001)
N3＝1283
3,
page 41 / 51
in (8Mpc) Nangle =
1282 background
Isotropic
UV: I21z=0.1
Zel’dovich
approximation:
= 15
Radiative Transfer
Ionization Equilibrium
Neutral Fraction:
nHI
X HI 
nH
CP-PACS 256PU
(理論ピーク ~80 GFlops)
~ 100 hour
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
Reionization History of an Inhomogeneous Universe
page 42 / 51
I21=0.1
Z=15
Z=9
Z=7
Z=5
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
Shadowing Effect
Inhomogeneous
注：同じ物質量でも，物質の
分布の仕方の違いにより，透過
する輻射量が異なってくる
page 43 / 51
Homogeneous
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 44 / 51
Shadowing Effect
(1) 3D RT
Raises Neutral Fraction
in
Optically Thick Region
(2) Inhomogeneity
Raises Neutral Fraction
in
All Region
Shadowing Effect
is prominent at
XHI  103
(TN, Umemura, & Susa 2001)
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 45 / 51
3D RHD: Cosmic Reionization (Gnedin 2000)
XHI
z=9
Cosmological HD ~1010-12 Msun
RT (Local Optical Depth 近似)
Star Formation
H, He
H2 form/dest
ngas
T
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
Optically Thin Variable Eddington Tensor
page 46 / 51
(Gnedin & Abel 2001)
XHI
LOD
OTVET
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 47 / 51
3D RHD (SPH): 銀河形成
(Susa & Umemura, in prep.)
NSPH  NCDM  215
zc ~ 1.5, 6

89
M

2
10
M
t
ot



Initial Condition by GRAFIC in
COSMICS
+ Rigid rotaion with λ=0.05
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史

page 48 / 51
UVが(相対的に)弱く，
遮蔽効果も効いている時
(ガスが冷えて星が形成される)
zc ~ 6
9
M

2
10
M
t
ot

青色は星
その他はガス温度

第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
UVが(相対的に)強く，遮蔽効果も効いていない時
(ガスが冷えず星が形成されにくい)
page 49 / 51
Mt ot  2 108 M
zc ~ 1.5 
青色は星
その他はガス温度


第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 50 / 51
Heterogeneous Multi Computer System
GRAPE6
8boards~8Tflops
PCI bus
~200MB/s
×8
Alpha Cluster
(16×alpha21264)
CP‐PACS
2048PU~600Gflops (peak)
Internal Communication
1：1→280MB/s
16 parallel
40MB/s (peak)
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史
page 51 / 51
5. 今後
詳細な観測データ
↓
詳細な輻射輸送・RHDシミュレーションの意義が
より大きくなる
多次元輻射輸送・輻射流体力学シミュレーションは黎明期
今後：
輻射輸送
・Ray Tracing の工夫
・Iteration の改良
RHD
・近い将来
・その先
Flux Limited Diffusion (FLD)
Variable Eddington Factor (VEF)
第15回 理論懇シンポジウム 国立天文台 2002年12月26日 中本泰史