衝撃波管問題（Riemannソルバー）

衝撃波管問題（Sod問題に代表される1次元Riemann問題）は、圧縮性流れの数値スキームの精度と安定性を評価するための最も基本的なベンチマークだからだよ。厳密解が解析的に得られるため、数値解との比較で数値拡散やオーバーシュート、衝撃波解像度を定量的に評価できる。

1次元Euler方程式の初期値問題として定式化するんだ。管の中央にダイアフラム（隔膜）があって、左側が高圧 $(\rho_L, u_L, p_L)$、右側が低圧 $(\rho_R, u_R, p_R)$ の状態。ダイアフラムを破ると、3つの波が発生する。

1次元Euler方程式の保存形はこうなる。

衝撃波を跨ぐ保存則から導かれるのがRankine-Hugoniot関係式だ。衝撃波速度を $s$ とすると、

接触不連続面の左右の領域で、圧力と速度は等しい（$p^ = p_L^ = p_R^$, $u^ = u_L^ = u_R^$）けど、密度は不連続になる領域だよ。この $p^*$ を求める方程式は非線形なので、Newton-Raphson法で反復的に解く。

Godunovの考え方は、各セル界面でのフラックスを局所的なRiemann問題の解から求めるというものだ。セル $i$ と $i+1$ の境界では、左状態 $\mathbf{U}_i$ と右状態 $\mathbf{U}_{i+1}$ から成るRiemann問題を解き、界面フラックス $\mathbf{F}_{i+1/2}$ を決定する。

厳密Riemannソルバーは毎回Newton法の反復が必要で計算コストが高い。そこで実用的な近似Riemannソルバーが多数開発されてきた。代表的なものを比較しよう。

HLLは波を2本しか考慮しないから、接触不連続面（3本目の波）の情報が失われるんだ。HLLCのCはContact waveのCで、3本目の波を追加することで接触不連続面を解像できるようになった。Toroの教科書に詳しい導出がある。

Roeの考え方は、非線形のEuler方程式を界面で線形化するというものだ。Roe平均状態を使ってヤコビアン行列 $\hat{A}$ を構成し、

膨張衝撃波（expansion shock）という非物理的な解を生成する場合がある。これはエントロピー条件に違反するもので、Harten-Hyman entropy fixで修正する必要がある。また、低マッハ数領域では過度な数値拡散が生じる問題もあり、Low-Mach Roe修正（Rieper fix等）が提案されているよ。

MUSCL（Monotonic Upstream-centered Scheme for Conservation Laws）法で2次精度化するのが標準だ。セル界面の左右状態を線形再構成で外挿する。

Sod (1978) の標準的な初期条件は以下の通りだ。管の長さを1、ダイアフラム位置を $x = 0.5$ とする。

rhoCentralFoamを使う。0/ディレクトリの初期条件はsetFieldsユーティリティで設定するのが簡単だ。system/setFieldsDictに左右の状態を定義して実行する。メッシュはblockMeshDictで1次元に1000セル程度配置すれば十分だよ。

4つの観点で評価するのが標準的だ。

1. 衝撃波の解像度: 衝撃波が何セルで捕捉されているか。理想は2-3セル。5セル以上は数値拡散が過大。

いくつかの標準的な問題を紹介しよう。それぞれ異なる側面をテストできる。

その通り。衝撃波がサイン波状のエントロピー変動に突入するという設定で、衝撃波後方に高周波の密度変動が生じる。1次精度や散逸の大きいスキームではこの高周波成分が消えてしまう。WENO5やMP5といった高次精度スキームの性能評価にはうってつけだよ。

ソフトごとの対応状況をまとめるよ。

良い指摘だ。KTスキームはRiemann問題を解かずに、局所的な最大波速を使って数値粘性を制御する中心差分系のスキームだ。実装が簡単でロバストだけど、接触不連続面の解像度はRoeやHLLCに劣る。OpenFOAMでRoeソルバーを使いたい場合は、rhoPimpleFoamにカスタムフラックスを実装するか、blastFOAMのようなサードパーティライブラリを利用する手がある。

Fluent GUIでの手順はこうだ。

もちろん。Fluent Journalファイルでパラメータスタディを回すのは実務でよくやる。初期条件の圧力比やメッシュ密度を変えて一括実行し、格子収束性を確認するのに便利だよ。

SU2は2次精度のMUSCL再構成がデフォルトで使えるし、JST（Jameson-Schmidt-Turkel）スキームの人工散逸パラメータも調整可能だ。設定ファイル（.cfg）での指定はこんな感じだね。

その通り。0に近いほど1次精度に近づいて振動は消えるけど散逸が増す。0.05-0.1が衝撃波管問題ではバランスの良い値だよ。

圧縮性流れの高次精度スキームとして最も研究されているのがWENO（Weighted Essentially Non-Oscillatory）スキームだ。基本的な考え方は、複数のステンシル候補から非線形重みを使って最適な補間を選択するというもの。

WENO5（5次精度）の場合、3つの3点ステンシルの重み付き組み合わせで5次精度の精度を実現する。衝撃波近傍では振動しないステンシルに自動的に重みが集中し、滑らかな領域では全ステンシルが等しく寄与して高次精度が発揮される。

残念ながら主要な商用CFDソフト（Fluent、STAR-CCM+、CFX）にはWENOは標準搭載されていない。使えるのはオープンソースや研究コードだ。OpenFOAMにはblastFOAMプロジェクトでWENO実装がある。また、JAXAのFaSTARやDLRのTAUコードにも実装されているよ。

DG法はセル内を高次多項式で近似し、セル界面でRiemannソルバーを使ってフラックスを評価するという、有限要素法とGodunovスキームのハイブリッドだ。任意の高次精度を達成でき、非構造格子との相性も良い。

衝撃波管問題へのDG法の適用で注意すべきは、不連続近傍でのGibbs振動への対策だ。以下のアプローチが使われる。

Fluentの最新版（2024R2以降）にPoly-Hexcore格子と組み合わせたDG法ベースのソルバーが試験的に導入されている。オープンソースではNektar++、Flexi、HOPRなどがDG法による圧縮性流れソルバーを提供しているよ。

たくさんあるよ。いくつか紹介しよう。

MHDでは磁場の影響で速い磁気音波、遅い磁気音波、アルヴェン波、接触不連続面の7つの特性速度があるからね。HLLDソルバーはこれを効率的に近似した手法で、宇宙物理の数値シミュレーションで広く使われている。

不連続面近傍のオーバーシュート/アンダーシュートは高次精度スキームの宿命だ。対策の選択肢を整理しよう。

Fluentでは「Solution Methods」でSpatial DiscretizationをFirst Order Upwindに一旦落として収束させ、その解を初期値にしてSecond Orderに切り替えるという手順が有効だ。また、Gradient LimiterのMethodをStandardからMultidimensional Limiterに変更すると振動が抑制されることもある。

Carbuncle現象は、強い垂直衝撃波の数値解が格子に依存した非物理的な不安定構造を示す問題だ。Roe系スキームで特に発生しやすい。衝撃波面が格子に平行に配置されたとき、衝撃波面上に不規則な凹凸が成長してしまう。

対策はいくつかある。

それも理由の一つだよ。AUSM系は衝撃波安定性に優れるから、汎用ソルバーのデフォルトには適しているんだ。

Godunov系スキームは低マッハ数で過度な数値拡散を生じる。これはフラックスの散逸項が $O(1/M)$ でスケーリングするためだ。対策として以下がある。

All-speed schemeは活発な研究分野で、STAR-CCM+のCoupled Flowソルバーはそれに近い設計思想だよ。非圧縮極限と圧縮性を一つのフレームワークで扱える点が強みだ。