流体力学 (格子法)・・・宇宙気体力学への応用
松田卓也
神戸大学・理学部・地球惑星科学科
平成 15 年 8 月 19 日
概要
数値流体力学の宇宙気体力学への応用について述べる

1

宇宙気体力学1.1

宇宙はほとんど気体でできている物質には一般に固体、液体、気体の三態がある。この地球上では、地球は固体、海の水は液体、空気は気体である。しかし宇宙においては、固体は惑星以外は星間ダストとして存在する程度で、液体も中性子星のような特殊な状態以外ほとんど存在しない。宇宙に存在する物質はほとん
どが気体 (=ガス) と考えてよい。気体は中性気体と電離気体 (プラズマ)に分類できる。宇宙にはどちらもあるが、プラズマの方がメインであろう。しかし磁場の影響を無視できる場合は、中性気体もプラズマも、流体力学的振る舞いには差異はなく、統一して取り扱うことができる。

1.2

宇宙流体に現れる流れの分類

宇宙流体を考える場合に、現象の時間依存性によって分類することができる。

1. 静的: 時間依存性がなく、かつ速度が 0 の場合。静的な重力平衡にある星などの場合。
1

2. 定常: 時間依存性がないが、速度が 0 でない場合。太陽風のパーカー解などがその例である。
3. 動的: 時間依存性がある場合。超新星爆発や、星の形成過程などの場合。
そのほかにも様々な分類が考えられる。

• 磁場のありなし。磁場がある場合は電磁流体力学方程式 (MHD) を解く必要がある。
• 粘性のありなし。粘性を無視できる場合は Euler 方程式、無視できない場合は Navier-Stokes 方程式を解く。宇宙現象においては分子粘性はふつう無視できる。乱流粘性が重要であるが、これを統一的に取り扱う方法はない。
• 一般相対論的効果のありなし。ブラックホールや中性子星など、重力が強い場合は、一般相対論的効果を考慮する必要がある。
• 輻射輸送効果のありなし。輻射が重要な役割を果たす場合は、輻射流体力学の手法を用いなければならない。

2

数値計算法の歴史と現状

2.1

基礎方程式系

ガスの運動を微視的に記述する方程式として Boltzman 方程式がある。
気体分子の平均自由行程と系の大きさの比を Knudsen 数とよぶ。Knudsen数が 1 程度以上の流れは、希薄気体と呼ばれ、Boltzman 方程式を解析的または数値的に解く必要がある。Knudsen 数が 1 より十分小さい場合は、気体は連続体とみなしてよい。その場合、気体分子の運動を記述する分布関数は、マクスウエルの平衡分布に近くなる。その場合、分布関数をマクスウエル分布からのずれが小さいとして、Chapman-Enskog 展開を行うと、Navier-Stokes 方程式が導かれる。そのとき粘性や熱伝導のような輸送現象が自動的に導かれる。Navier-Stokes 方程式において粘性を無視した方程式を Euler 方程式と呼ぶ。それは次のような式で表される。

2∂Q ∂E ∂F∂G++++H =0∂t∂x∂y∂zρwρvρuρ ρuw  ρuv  ρu2 + p  ρu 

2Q =  ρv  , E =  ρuv ,  F =  ρv + p  , G =  ρvw  , ρρw2 + p  ρvw  ρuw
 ρw (e + p)w(e + p)u(e + p)ve0ρκx
H=ρκyρκzρ(uκx + vκy + wκz )p1
e=+ (u2 + v 2 + w2 )(γ − 1)ρ 2
ここで ρ, u, v, w, e, p, κx , κy , κz は密度、速度の x,y,z 成分、運動エネルギーを含む全エネルギー、圧力、外力の x,y,z 成分を表す。Q, E, F, G, H
は式に表されたような縦ベクトルで、Q は保存変数、E, F, G は流束ベク
トル、H はソース項である。このように、流速の発散として表された、方程式の表示を保存形とよぶ。
Euler 方程式は非粘性の圧縮性気体の運動を記述する方程式で、双曲型の方程式である。一般に二階の偏微分方程式は、楕円型、放物型、双曲型に分かれ、その解法は非常に異なる。楕円型の代表としては重力ポテンシャルを記述するポアソン方程式、放物型の代表としては熱伝導方程式がある。Euler 方程式は非線形の双曲型偏微分方程式であり、時間発展とともに、衝撃波や接触不連続面のような不連続を発生させる。そのような不連続を含む解を弱解とよぶ。

2.2

Euler 方程式の数値解法の諸問題

Euler 方程式を数値的に解く場合、上記のような保存形の方程式を離散化して解く方法と、非保存形の方程式を解く方法がある。昔の教科書では、非保存形の方程式が用いられた場合もあったが、近年の数値流体力学
では保存形を用いる場合が多い。保存形の特徴は、質量、運動量、エネルギーなどの保存量が、差分スキームなどを用いて偏微分方程式を離散化し3た場合でも保存することである。衝撃波の位置を決める Rankin-Hugoniot
条件が、数値解においても厳密に満たされるので、衝撃波の位置が正確に決まる。一方、非保存形を用いる利点は、保存量の保存性をモニターすることにより、計算精度をチェックすることができることであろう。しかし保存量が保存したからといって、計算が必ずしも正しいとは限らないことは注意する必要がある。
Euler 方程式の時間発展を計算する場合に、陽的解法と陰的解法がある。
陽的解法では、時間微分項以外を右辺においた場合、右辺はすべて古い時間の量、つまり既知量で表される。したがって時間積分は簡単である。
一方、陰的解法では右辺に未知量が現れるので、未知量を求めるには行列反転やニュートン・ラプソン法などの緩和法を用いなければならず、アルゴリズム的にも、計算時間的にも大変である。
陽的解法では、時間刻みをいくらでも大きくとることはできず、クーラン条件を満たさなければならない。そうでないと計算は不安定になり、解は発散する。陰的解法ではこのような制限はなく、原理的には時間刻みはいくら大きくてもとれる。このことがアルゴリズムの難しさや、計算時間の長さを相殺するので、定常問題では多く使われてきた。
しかし、時間変動が大きく、動的に変化するシステムでは、精度の関係から時間刻みは小さくとらねばならないので、この場合は陰的解法を採用する利点は薄れる。

2.3

数値計算法の歴史・・・前史

流体力学の基礎方程式を数値的に解くというアイデアは戦前からあり、Richardson により天気予報などが考えられたが、成功しなかった。コンピュータが無かったことと、安定性条件が知られていなかったからである。戦後になり von Neuman たちが電子式コンピュータを開発して初めて、数値流体力学が現実のものとなった。
von Neuman と Richtmyer (1950) は、圧縮性流体の計算途中に発生する衝撃波を安定にとらえるために、人工粘性項を導入した。Lax と Friedrichs(1954) は、差分法を工夫することにより、人工粘性項なしで衝撃波を安定にとらえることに成功した。Lax-Friedrichs 法が１次精度であり、解が非常になまるという欠点があり、Lax と Wendroff (1964) は、それを改良するために２次精度の手法を開発した。これが有名な Lax-Wendroff 法で、現在でも MHD 業界では、その簡便性のために使われている。MacCormack4

(1969) は、一見奇妙な差分スキームを採用したが、これは Lax-Wendroff法と同じく２次精度の手法であり、一部の業界で愛用された。
Lax-Wendroff 法、MacCormack 法はいずれも簡便であるので、実際問題にも大いに利用されたが、衝撃波や接触不連続面で、解が振動するという欠点を持っていた。厳密解が単調な場合に、数値解も単調であることの要請を単調性の保存という。後に単調性を保存する手法が模索されるようになる。
これらの手法はすべて非保存形で陽的解法であった。Beam と Warming(1976) は、４次の人工粘性を持つ陰的解法を開発した。Beam-Warming 法
は 70-80 年代に、数値流体業界で大流行となり、猫も杓子も Beam-Warming法であり、Beam-Warming 法でなければ計算法でないといった状態であった。Beam-Warming 法に基づく専用コンピュータを開発しようという意見すらあった。

2.4

一般曲線座標の導入

もうひとつ重要な発展は、一般曲線座標の導入である。宇宙気体力学
ではあまり問題にならないが、数値流体力学の工学的な応用では、飛行機や自動車、タービンなどの複雑な形状をした物体表面で境界条件を与える必要がある。これらをデカルト座標で表現すると、でこぼこな階段状になりおもわしくない。そこで物体表面に沿うような一般曲線座標が導入され、方程式がデカルト座標から一般曲線座標に変換された。
一般相対性理論を知っているものならば、一般曲線座標は良く知られたものである。しかし、数値流体業界にとっては初めてのものであり、当初は様々な誤解や無理解があり、方程式形が変な形になった。その影響は現在にまで及んでいる。それぞれの分野が他の分野の常識を知らないために、学問が変な方向に発展することは往々にしてあることだ。
一般曲線座標の導入に伴って、一般的な格子を能率良く生成する、格子生成法がひとつの大きな分野になった。デカルト座標にしろ一般曲線座標にしろ、それは構造格子と呼ばれている。規則正しい形をしているからである。それにたいして２次元なら三角形や多角形、３次元なら多面体を基礎とする非構造格子が現れた。現在の最先端は非構造格子である。
しかしコンピュータの発達に伴い、メモリ量が飛躍的に増えたため、ふたたび計算が簡単なデカルト座標が復活の兆しを見せている。格子数を十分大きくとれば、物体表面もよく近似することができるからである。デ5

カルト座標への先祖帰りは、まだ試行的段階ではあるが、宇宙気体力学では、十分に有用であろう。

2.5

中心差分から風上差分法へ、そして TVD へ

つぎの流行は風上差分法であった。先に述べた手法はいずれも中心差分法かそれに近いものであったが、風上差分法は考えている格子に影響を及ぼす格子からの影響を重視する手法である。
風上差分法自体は新しいものではなく Courant-Isaacson-Rees 法 (1952)として知られていた。この手法では、考えている格子に流体が流れ込んでくる方向 (風上) の影響のみを考える。しかも１次精度であった。圧縮性流体においては、影響は流速だけではなく、音速でも伝播する。音波の影響をきちんと考える手法として、Godunov (1959) は格子間境界でRiemann 問題を解く手法を開発した。Riemann 問題については後で述べる。Godunov の手法は当初はそれほど顧みられなかったが、後に大きな影響を及ぼすことになる。
ここで後の発展を制約する Godunov の定理というものを述べておく。
次の条件の内、二つしか同時に満足することはできない。

1. 方程式の非線形性
2. 単調性の保存
3. 2 次精度以上
線形のモデル方程式の場合、条件 2 と 3 を同時に満たす、つまり単調性を保存した 2 次精度以上の高次精度の手法を開発することができる。しかし Euler 方程式は非線形であるので、2 と３を同時に満たす数値解法は存在しないことになる。この困難を回避するには、衝撃波や接触不連続面のように単調性の保存が重要な局面では１次精度になり、他の場所では高次精度になるような手法が望まれる。それを手でやるのではなく、自動的に行う手法が望ましい。
それが 80 年代以降に爆発的に発展した風上差分法と、それからさらに発展した TVD 法である。これらの手法は、以前の手法と異なり、きわめて数学的であり、行列の固有値や Riemann 問題などの知識を必要とする。
TVD とは Total Variation Diminishing(全変動減少) の略である。全変動とは解の最大値と最小値の差の絶対値の和のことである。たとえば１6次元空間を衝撃波がその形を変えることなく伝播している場合を考える。
この場合、全変動は時間的に一定である。しかしこの問題をたとえば LaxWendroff 法で解くと、衝撃波の部分に振動が現れて、全変動は時間的に増大して行くであろう。そこで全変動が時間的に一定か、減少するような計算方法を考えれば、そのような振動は現れないであろう。つまり単調性が保存される。アルゴリズムに TVD 性を要求すると安定でかつ不要な振動が現れないようなスキームができるであろう。このような考えの基に、さまざまな TVD スキームが考案された。
もっとも TVD 性がはっきりしているのは１変数の方程式か、線形の方程式系に限られるので、TVD 法といえど万能であるという証明はない。
現在の数値流体力学業界では、アルゴリズム研究は下火になっている。
より高度なアルゴリズムを考えるよりは、コンピュータの発展を待って、格子数を増やす方が手っ取り早いからである。先に述べたように、高度な一般曲線座標や非構造格子を用いるよりは、格子数の多いデカルト座標を用いる方が、計算は簡単だし、コンピュータの発展とともに、いくらでも精度は増していくのである。(もっともそれは一部の研究者の意見で、商業用流体力学パッケージなどでは、複雑な非構造格子が採用されている。)
そこで一部の数値流体力学研究者は、輻射流体やその他の応用分野に研究の矛先を転じているのが現状である。つまり磁場や輻射を考慮しない、単純な数値流体力学はすでに成熟段階に達したと言っても良いであろう。もちろんさらなる研究の余地はあるが、収穫逓減の法則の段階に入っているのである。

2.6

天文業界はガラパゴス島であった

以前、筆者は、数値流体力学の手法に関して、天文業界は世の中の進化に取り残されたガラパゴス島のようなものだと述べたことがある。数値流体力学業界は上に述べたように、アルゴリズムが画期的に発達して変遷を繰り返したのだが、天文業界ではいつまでも von Neumann の人工粘性や Lax-Wendroff 法が幅を利かせていた。天文業界では外の世界ではほとんど使用されない SPH 法などという新種まで発生していたのである。
筆者は 90 年代の始めまで、幸か不幸か京都大学の航空工学科に所属しており、半ば義務的な要請もあり、数値流体業界とつきあいがあった。そこで天文業界のガラパゴス的現状を認識したのである。このことは日本

7だけではなく、世界的に見てもそうであった。しかし例外もあり、たとえば天文業界出身の van Leer や Woodward は数値流体業界に進出して成功を収めている。もっとも Woodward たちが開発し、西欧の天文数値流体業界では標準的手法となっている PPM 法は、数値流体業界ではあまり顧みられていない。この手法は Godunov 法に起源を有し、van Leer により発展された手法を元にしたものであり、出自は正しい。
ともかく異分野で交流することは意義がある。これは言ってみれば、貿易のようなものである。売買自体は何の価値も生み出さないが、あるところで産出する産物を、産出しないところに持っていくだけで儲かるのである。

2.7

私のアルゴリズム遍歴

ここで数値流体力学手法に関する、筆者の個人的な遍歴を披露する。
筆者は 1965 年に京都大学理学部物理学第二教室の天体核物理学研究室に入った。選んだテーマは超巨大質量星の一般相対論的重力崩壊の数値シミュレーションであった。選んだ手法は、当時の天文業界の標準的手法である von Neumann の人工粘性を用いた、１次元ラグランジュ法であった。格子数は３０でしかなかった。人工粘性の人工的効果に悩まされた。
７０年代初頭には Lax-Wendroff 法を用いたが、衝撃波面での振動に悩まされた。当時、数値流体力学研究の最先端はロスアラモス研究所であった。そこでは PIC, FLIC などさまざまな手法が試されていた。筆者が格子数 30 で呻吟しているときに、彼らはすでに２次元計算を行い、美麗なコンピュータ・グラフィックスを駆使していたのだ。
1974 年に S.-A.Sorensen が筆者の研究室に来たときに、FLIC 法を用いて、初めて近接連星系の降着円盤の２次元数値シミュレーションを行った (Sorensen, Matsuda & Sakutai, 1975)。これはかなり成功したと思うが、今から思えば、人工粘性の大きさの故に、重要な特徴 (渦状衝撃波)が、消えてしまっていた。このことは van Leer から直接に批判された。
1975 年に筆者が渡英したときに、ロスアラモスで研究されていた一つの手法、オイラー・ラグランジュ混合法を開発した。これは２次元空間の一方をオイラー的に、他方をラグランジュ的に解く方法である。正直言って、この手法は成功しなかった。
帰国後は MacCormack 法に少し手を染めた。８０年代に入って、現在は東北大学の教授をしている沢田恵介君が研究室に入ってきた。彼は文8献調査からまずは Beam-Warming 法と一般曲線座標を導入した。その後、
風上差分法の一つである Osher 法を学び、それを近接連星系の降着円盤の流れに応用した。始め１次精度の手法を採用して Sorensen たちと同様な結果を得ていたが、２次精度に拡張すると、降着円盤の上に渦状の衝撃波が存在することを発見した (1986, 1987)。これは筆者にとって衝撃的な発見であった。Sorensen たちが見逃していた構造を、より高次のアルゴリズムを採用することにより発見したのである。もう一つの要因は、当時、ベクトル型スーパーコンピュータが日本でも使えるようになったことが挙げられる。
筆者はその後、Osher 法をさまざまな問題に適用して成功した。しかし
太陽風と星間ガスの MHD 的相互作用の研究の時に、Osher 法では MHDへの拡張が困難なことを知り、van Leer の開発した Flux-Splitting 法をMHD に適用する城之内たちの手法を採用した。Flux-Splitting 法は後にAUSM 法へと発展し、筆者はやはり城之内、嶋たちが開発した AUSM 法の一変種である SFS 法を以後採用している。もっともこの手法を採用しているのは、今では筆者たちだけのようである。

2.8

風上差分法と Riemann 問題

圧縮性流体における風上差分法を理解するには、Riemann 問題について理解しなければならない。Riemannm 問題とは次のようなものだ。図に示されたように一つのガス容器がダイヤフラムで分離されていて、一方の部屋 (図では左側) には高圧のガス、他方は低圧のガスで満たされているとする。ある時刻 t = 0 にダイヤフラムが破られたとする。このような容器は、実験的には衝撃波管と呼ばれている。1 次元衝撃波管問題は、圧縮性数値流体力学のテスト問題である。この問題を数学的には Riemann問題と呼ぶ。
ダイヤフラムが破られると、衝撃波管の右の方向に衝撃波、左の方向に希薄波 (膨張波) が伝わっていく。元の高圧ガスと低圧ガスの境界は接触不連続面として存在する。圧力差が小さい場合を線形 Riemann 問題とよび、衝撃波と希薄波は下の図に示されるように、音波となる。下の図は横軸に空間座標、縦軸に時間座標を示した時空図であり、いろいろの波の伝播を表すのに便利である。図に u + c, u − c と表されたものは、流速 u に乗って右と左に伝播していく音波を表す。図に u と示されたものは、流れとともに流される接触不連続面である。これらの線を特性線と9

High pressureu-c

Low pressure

tuu+cx
図 1: 線形 Riemann 問題。

よび、u + c, u − c, u を固有値とよぶ。
Riemann 問題を数値的に解くことが風上差分法の要点である。その方法はおおざっぱにいって次の 3 つに分類できる。

1. Godunov の方法: Riemann 問題を数値的に正確に解く。正確に解くとは、高次方程式とか超越方程式を逐次近似法を用いて解くことをいう。
2. Osher 法: Riemann 問題を Riemann 不変量を使って近似的に解く。
3. Roe 法: 線形 Riemann 問題を正確に解く。
もちろんこの他の近似解法も存在する。どの方法を採用するかがアルゴリズムを特徴づけている。圧縮性数値流体力学の詳細に関しては、いくらもすぐれた教科書があるので、それに譲る。

10

3宇宙気体力学

3.1

近接連星系内のガス流

本稿では宇宙流体力学一般についてはふれない。とくに筆者が研究してきた、近接連星系内のガス流に関連して、その数値流体力学的研究について紹介する。

3.2

降着円盤

近接連星系を構成する星の周囲に存在する降着円盤は、宇宙物理学の重要なテーマである。降着円盤が存在する天体としては、激変星、新星、X 線星などがある。
これらの降着円盤内のガスの流れを記述する標準理論としては、Shakuraと Sunyaev によるものがある。これは、乱流粘性の現象論的パラメターからアルファ円盤モデルとも呼ばれている。活動銀河中心核に存在するとされる巨大ブラックホールまわりにも降着円盤が存在する。その降着円盤を表す解としては近年 ADAF モデルが注目を浴びている。これは一種の自己相似解である。
標準円盤モデルも ADAF も、近接連星系に存在する降着円盤を表す解としては、粗い近似でしかない。近接連星系では降着円盤にガスを供給する伴星が存在するが、その重力的影響が先の解では取り入れられていない。それを考慮するには、数値計算を行う必要がある。
以下では近接連星系のガス流の数値シミュレーションについて、その歴史と筆者の研究について述べる。

3.3

近接連星系内のガス流の数値シミュレーション的研究：
私的歴史

近接連星系内のガス流の数値シミュレーション的研究は Prendergast(1960) により始められた。もちろん当時はコンピュータも数値計算法も十分に発展していなかったので、見るべき成果はなかった。Prendergast
は同時に、棒状銀河内の流れについても先駆者的研究を行っていることは注目すべきである。Prendergast & Taam (1974) は、Prendergast の発明した Beam 法を用いて、はじめて近接連星系内のガス流のもっともら

11図 2: FLIC 法によるシミュレーション。

しい計算を行った。Beam 法は現在の格子ボルツマン法の先駆ともよぶべきものである。流れを計算するのに Euller 方程式ではなく、Boltzmann方程式を用いるのであるが、速度分布関数を静止粒子と上下左右方向への一定速度を持つ粒子に限定して表現する。それらの粒子の平均速度で流速が決まるのである。この手法はそのままでは、非常に人工的な拡散が大きいという欠点を持っている。
当時、筆者は Sorensen たちと同様な研究を進めていたので、彼らの論文はかなりショックであった。しかし、彼らの計算では、質量降着星のサイズが大きく、これは X 線星や激変星ではなく、むしろアルゴルに対応することを知って、筆者は少し安心した。
Sorensen, Matsuda & Sakurai (1974, 1975) は質量降着星をコンパクト星とするために、小さく設定した。といっても現実のコンパクト星よりは非常に大きいのであるが、L1 点から流れ出すガスの作るリング半径よりも、小さくとることにより、降着円盤の形成を可能にした。計算は公転面のみに限定した 2 次元計算であり、1 次精度の FLIC 法を採用した。
図 2 に、当時のわれわれの結果を示す。図は質量比 1 の近接連星系におけるロッシュローブあふれ流の計算で、密度と速度ベクトルを示す。この結果では、L1 点からコンパクト星に向かって伸びる、「象の鼻」とよばれる L1 流は現れているが、降着円盤はあまりきれいには出ていない。

12

ちょうど同じ頃 Lin & Pringle (1976) は、粘着粒子法とでもよぶべき手法を用いて同じ計算を行った。この手法は粒子法であるが、補助的に格子も採用する。同一格子内に入った粒子は衝突すると仮定される。衝突後の各粒子の軌道は、衝突時の運動量、角運動量の保存則から決められる。この手法は SPH 法の先駆と考えられなくもないが、むしろ Bird により開発された直接シミュレーション・モンテカルロ法 (Direct Simulation
Monte Carlo Method = DSMC) の先祖と考えた方がよいだろう。DSMCにおいては、粒子数を多くとり、一格子ないの粒子の衝突数は衝突確率から計算される。また衝突後の粒子の運動方向は、運動量保存則を満たす範囲内でランダムに選ばれる。筆者は現在、DSMC の宇宙流体力学への応用について研究中である。

3.4

降着流の現代的計算

Sawada, Matsuda & Hachisu (1986, 1987) は、Osher 法と富士通 VP200,400 スーパーコンピュータを用いて、再度、降着円盤の 2 次元数値シミュレーションを行った。その結果を図 3 に示す。図は質量比が 1 の近接連星系での計算で、密度分布を示している。はじめ 1 次精度で計算を行ったが、そのときは見えなかった渦状の衝撃波が、精度を 2 次にあげると出現した。一定程度以上の計算精度は必要という教訓である。図 4 には、さらに精度を上げた最近の 2 次元計算の結果を示す。
降着円盤上に渦状衝撃波があることは長い間、信じられなかった。3 次
元計算を行えばなくなるという意見もあった。それにたいして Sawada &Matsuda (1992) は 3 次元計算においても、渦状衝撃波は存在することを示した。図 5 には最近のわれわれによる 3 次元計算の結果を示す (Fujiewara
et al. 2001)。図はコンパクト星まわりにある降着円盤の、ある等密度面とそのうえの流れ模様、公転面での流れ模様を LIC 法で可視化したものである。あきらかに渦状の衝撃波が見られる。

3.5

渦状衝撃波の観測による発見

渦状衝撃波は、降着円盤理論の一つのネックとなっている角運動量輸送の問題を解決する可能性も秘めている。たとえそうでなくても、渦状衝撃波が存在すれば、なんらかの観測的な影響を及ぼすはずである。しかし一部の観測者の熱心な探索にも関わらず、渦状衝撃波は 1997 年まで13見つからなかった。
それが見つかったのは Steeghs, Harlaftis & Horne (1997) による、激変星 IP Peg. のドプラー・トモグラフィーによる観測であった。トモグラ
フィーとは、間接的な可視化手法である。たとえば X 線トモグラフィー
では、人体の周囲から X 線を照射して、その吸収量から人体の断面図を計算で求めることができる。
ドプラー・トモグラフィーでは、降着円盤などの高温のガスが発する水素やヘリウムの輝線の時間変動を観測する。人体では X 線の照射光は人体の周りを回るのだが、連星系の場合は相手が回ってくれる。こうして観測した輝線スペクトルの時間変化から、速度空間における輝線の分布が決まる。位置空間の分布ではないので、ドプラー図からだけでは、位置空間の構造はユニークには決まらない。しかしいろいろ興味ある情報が引き出せる。たとえば降着円盤の上に渦状に光る分布が存在すると、ドプラー図にもそれは存在する。実際 Steeghs たちが発見したのはそのようなものである。
図 6 の左にはドプラー図上でのヘリウムの輝線の明るさを示した。輪のようなものが降着円盤を表している。それが軸対称なものであれば、降着円盤も軸対称である。しかし図から明らかに、降着円盤は軸対称ではなく、渦状の構造が見られる。また伴星の表面が降着円盤の中心部に照らされている様子も見える。それが少し左に偏っているのは、伴星表面の流れによるものと推測される。このようにドプラー図からは、直接的ではないが、様々なことが読みとれるのである。
図 5 には Yukawa, Matsuda & Boffin (1997) による、SPH 法を用いた 2次元計算を元にして作った、ドプラー図も示す。観測との一致はよいことが分かる。とはいえ、このドプラー図の作り方は、単に密度の大小と輝線の強弱を比例させただけであり、本当は、輻射輸送を含むもっと複雑な計算を行う必要がある。

3.6

まとめ本稿では数値流体力学の歴史、およびその宇宙物理学への応用について述べた。そして筆者がかかわった近接連星系の降着円盤の数値シミュレーションにしぼって解説した。数値シミュレーション、数値実験の重要な効用として、単に現象を再現するだけではなく、観測的には知られていない現象を予言することであろう。

14

図 3: 2 次精度の Osher 法による計算

図 4: さらに精度をあげた計算15図 5: 最近の 3 次元計算

図 6: 観測によるドプラー図 (左) と 2 次元 SPH 計算による理論的ドプラー
図 (右)16