人工粘性 — CAE用語解説

良い質問だ。衝撃波のような急峻な圧力・密度変化を伴う現象を数値的に解く時、離散化した格子点上では物理量が不連続に「飛び跳ねる」現象（オーバーシュート/アンダーシュート）が起きる。これは数値的不安定性を引き起こし、計算を破綻させる。人工粘性は、この飛び跳ねを「なだらかに拡散させる」ダンパーの役割を果たし、衝撃波面を数値的に捉えられるようにするんだ。

まさにそのトレードオフが核心だ。人工粘性が強すぎると衝撃波面がぼやけ、厚みが数セル以上になってしまう。逆に弱すぎると不安定になる。目標は、衝撃波の厚みをせいぜい1〜2セル幅に収めつつ、計算を安定させること。例えば、超音速流れの解析では、衝撃波の内部構造そのものには興味がなく、その前後のジャンプ条件（ランキン・ユゴニオ条件）が正しく満たされることが重要になる。

von Neumann と Richtmyer が提案した古典的で線形な形式はこれだ。圧力に追加する人工粘性項 \( q \) は、速度勾配の2乗に比例する。

有限体積法（FVM）ベースのソルバー、例えばOpenFOAMの`rhoCentralFoam`や、商用コードのFluentの密度ベースソルバーでは、セル界面での流束計算に組み込まれることが多い。具体的には、圧力項 \( p \) を \( p + q \) に置き換えた「有効圧力」を用いて流束を計算する。離散化では、速度勾配 \( \frac{\partial u}{\partial x} \) はセル中心の速度値から、1次風上差分や中心差分で評価される。

その通り。これは「メッシュ依存性」という大きな問題を引き起こす。メッシュを細かく（\( \Delta x \to 0 \)）すると理論上 \( q \to 0 \) となり、衝撃波面で再び振動が発生する。これを防ぐため、現代的な手法ではメッシュ幅に依存しない、あるいは弱く依存する非線形人工粘性や、Jameson型の圧力に基づく人工粘性が使われる。例えば、Aerojet Rocketdyneの内部コードでは、\( \Delta x \) の代わりに局所的な圧力勾配を用いた修正項を加えている。

関連しているが目的が異なる。LS-DYNAやAbaqus/Explicitで言う人工粘性（バルク粘性）は主に2つの役割がある。1つは衝撃波の制御（これは流体と同じ）。もう1つは、減積分要素（1点積分）で発生するゼロエネルギー変形モード、いわゆる「アワーグラスモード」の抑制だ。こちらは体積ひずみ速度に比例する線形項と、その2乗に比例する二次項の和で表されることが多い。定数のデフォルト値は、線形項係数が0.06、二次項係数が1.5程度だ。

デフォルト値は一般的なケースで「とりあえず動く」ように設定されている。しかし、特に極超音速（マッハ5以上）や爆轟問題では調整が必須だ。まずは係数を0にして計算を走らせ、衝撃波面の前後で解が激しく振動することを確認せよ。次に、係数を少しずつ（例えば0.1刻みで）上げていき、振動が収まる最小の係数を見つける。これが「必要最小限の人工粘性」を決める第一歩だ。

必ず既知の理論値や実験データと照合する。最も基本的な検証は、1次元衝撃波管問題（Sod問題など）だ。衝撃波の伝播速度、衝撃波後方の圧力・密度が、ランキン・ユゴニオの関係式から求まる理論値とどの程度一致するかを定量評価する。例えば、圧力の相対誤差が5%以内に収まっているか。商用ソフトでは、Ansys Fluentのチュートリアル「Supersonic Flow Over a Wedge」が、くさび角とマッハ数から決まる衝撃波角を検証する良い例だ。

同じ係数設定で、グリッド数を2倍、4倍と系統的に細かくした複数のメッシュで解析を実行する。衝撃波の厚み（圧力が急上昇する領域のセル数）と、衝撃波前後のジャンプ量（圧力比など）を記録する。理想は、メッシュを細かくしてもジャンプ量が一定値に収束し、かつ衝撃波厚みが1〜2セルに保たれること。もしメッシュを細かくするたびに振動が再発するなら、使用している人工粘性の形式がメッシュ依存性を強く持っている証拠だ。

目的と設定方法が根本的に異なる。Fluent（密度ベースソルバー）では、主に「Flux Type」として`Roe`を選択した場合の「Roe-FDS」スキームに、内在する人工粘性（行列拡散）が含まれる。ユーザーが直接いじるのは「Pressure-Based」なら「Skewness Correction」や「Multigrid」の設定、あるいは「Second Order」から「First Order」に落とすことが実質的な人工粘性の増加になる。一方、Abaqus/Explicitでは、ステップ定義内の「Bulk Viscosity」セクションで線形係数と二次係数を直接入力する。デフォルトはそれぞれ0.06と1.5で、金属成形などの大変形解析ではこの値を大きく（例えば0.12と3.0）して安定性を高めることがある。

OpenFOAMはユーザーがソースコードレベルで実装を確認できる好例だ。超音速ソルバー`rhoCentralFoam`では、`centralCourantNo`スキームが使われ、その中で`viscous`というスイッチで人工粘性のオンオフが制御される。具体的な計算は`TVD`（Total Variation Diminishing）の考え方に基づく「限界関数」の形で実装されており、`fvSchemes`ファイルで`divSchemes`を`Gauss limitedLinearV 1`などと設定することで、その強度（1が最大制限）を間接的にコントロールする。より直接的な古典的形式を試したければ、ユーザーが自分で`q`項をソース項として追加する必要がある。

その通りだ。NASAで開発された`OVERFLOW`や、日本で広く使われる`UPACS`（現`Cflow`）のような航空宇宙特化型ソルバーでは、人工粘性は単なる係数調整ではなく、スキーム選択そのものに深く関わる。例えば`JST`（Jameson-Schmidt-Turkel）スキームは、2次精度の中心差分に、4次差分の人工粘性項（滑らかな領域での影響を最小化）と2次差分の人工粘性項（不連続部での減衰）を巧みに組み合わせている。ユーザーは、この2次項と4次項それぞれに掛ける係数（σ2, σ4）を調整する。経験則では、σ4は1/256程度、σ2は0.5〜1.0の範囲で設定される。

それは最も典型的な対処法の一つだ。特にAbaqus/Explicitで大変形・接触問題を解く際に「負の体積」エラーが出たら、まずバルク粘性の二次係数をデフォルトの1.5から2.0や2.5に上げてみる。これは要素が急激に圧縮される際のエネルギー散逸を増やすことで、変形を現実的に「遅らせ」、極端な形状歪みを防ぐ。ただし、過度に強くすると（例えば係数5.0以上）、問題の物理的本質が失われ、エネルギー減衰が大きすぎる非現実的な結果になるので注意が必要だ。

それは「過剰な人工粘性」による衝撃波面の「ぼやけ」が原因の可能性が高い。粘性が強すぎると、衝撃波面が複数セルにわたって拡散し、その分、不可逆的なエントロピー生成（つまりエネルギー損失）が過大評価される。その結果、衝撃波通過後の圧力・温度が理論値より低く、密度が高く算出される。対策は、使用しているスキームの人工粘性係数を少しずつ下げて、衝撃波面の厚みが1〜2セルになる最小値を見つけること。Fluentで`Roe`スキームを使っている場合、`Solution Methods`の`Roe-FDS Flux Type`で`None`から`Entropy Fix`に変えるだけでも改善されることがある。

まさに古典的人工粘性のメッシュ依存性が原因だ。先述の式

ここが非常に重要で、見過ごされがちな点だ。バルク粘性はエネルギーを散逸させるため、部材の変形に必要な塑性仕事（エネルギー）の一部を「余計に」吸収してしまう。その結果、部材が実際より硬く、脆く見える可能性がある。評価方法は、感度解析を行うことだ。線形係数と二次係数をデフォルト値の±50%の範囲で変化させ、部材の吸収エネルギーや最大変位量がどの程度変化するかを調べる。もし係数を変えることで結果が大きく（例えば10%以上）変わるなら、その解析結果は人工粘性の設定に強く依存していると言え、結果の信頼性には注意が必要だ。LS-DYNAのマニュアルでは、このような感度解析を行うことを強く推奨している。