OpenFOAM後処理
理論：数値を取り出す仕組み
概要
OpenFOAMの後処理はParaViewで見るだけではないのですか？
見るだけならParaViewで十分だが、設計に使う数値は「関数オブジェクト」で取り出すのが基本だ。system/controlDictのfunctionsに、forces（力とモーメント）、forceCoeffs（抗力・揚力係数）、probes（点の時刻歴）、sets・surfaces（線や面上の分布）、fieldAverage（時間平均）、yPlus、wallShearStressなどを書いておくと、計算と同時にpostProcessingフォルダへ書き出される。計算後に postProcess -func で同じ処理を既存の結果に適用することもできる。数値を取り出すときに間違えやすいのは、基準となる面積や速度、密度の与え方、時間平均を始める時刻、出力の間隔で、どれも結果の桁や意味を変えてしまう。
係数の定義と単位
F_d は抗力、ρ∞・U∞ は基準の密度と速度、A_ref は基準面積、D は円柱の直径、f は渦放出の周波数、St はストローハル数、p_kin は非圧縮ソルバーの運動学的圧力。
paraFoamとParaView
OpenFOAMの結果を開くparaFoamは、実はParaViewを起動するスクリプトで、OpenFOAMの形式を読むリーダーを組み込んでいる。ParaView自体は米国の国立研究所とKitware社が開発した可視化ソフトで、OpenFOAMとは別のプロジェクトだ。ParaView付属のリーダー（.foamファイル）を使えばparaFoamを使わずに開けるので、Windowsのパソコンで計算機上の結果を見るときによく使われる。
計算例
Cd = 1.33・St = 0.165は設定を検算するための例示入力で、本サイトの計算実績ではない。基準値との精密比較には同じ閉塞率・境界条件のデータが必要。Foundation v13ではfoamPostProcess等、利用版に合うコマンドを使う。
例：Re = 100の円柱（直径0.01 m、流速1 m/s、動粘度1×10⁻⁴ m²/s、メッシュ厚さ0.002 m）
| 項目 | 値 | 備考 |
|---|---|---|
| 抗力係数（例示値） | ≈ 1.33 | 時間平均 |
| ストローハル数（例示値） | ≈ 0.165 | |
| メッシュ厚さ分の抗力 | 1.33×10⁻⁵ N | rhoInf = 1 |
| Aref = D と誤ったときの表示 | 0.00266 | 500分の1 |
| 渦放出の周波数・周期 | 16.5 Hz / 0.061 s | |
| 1周期20点のサンプリング間隔 | 0.003 s | |
| 10周期の平均に必要な時間 | 0.61 s | 初期の過渡を除いた後 |
例示値は代表的な数値。forceCoeffsの基準面積は直径×メッシュ厚さ、渦放出の周波数は揚力係数の時刻歴をFFTして求める。
抗力係数の時刻歴を平均したら1.45になり、文献の1.33と9％ずれました。メッシュを細かくするべきですか？
その前に、平均のとり方を確かめよう。円柱の後ろの渦放出は、流れ始めてからしばらくは成長途中で、周期的な状態になるまでに数十周期かかることがある。最初から平均すると、立ち上がりの過渡が混ざってずれる。fieldAverageやforceCoeffsの平均は、揚力係数の振幅が一定になってから始め、整数個の周期（例えば10周期、この例では0.61秒）で平均する。周期の途中で切ると、振動の分だけ平均がぶれる。そのうえで差が残るなら、計算領域の大きさ（側面の壁が近いと抗力が増える。円柱の直径の20倍以上の幅が目安）、出口の位置、時間刻み（クーラン数）、そしてメッシュを順に確かめる。
実務での進め方
- 必要な数値（力・係数・点の時刻歴・断面の分布）を決め、controlDictの関数オブジェクトに書く。
- 基準の面積・長さ・速度・密度を、2次元ならメッシュ厚さを含めて設定する。
- 非定常なら、現象の周期の1/20以下の間隔で出力し、過渡が終わってから整数周期で平均する。
- yPlusやwallShearStressで壁の扱いが妥当かを確かめる。
- 例示値や実験と比べ、差の原因を領域・時間刻み・メッシュの順に調べる。
「並列計算の結果はそのままでは開けない」
decomposeParで分割して並列計算した結果は、processor0、processor1…のフォルダに分かれて保存される。そのままParaViewで開くと、最初の時刻しか見えない、と戸惑うことがある。reconstructParで結合するか、.foamファイルを開くときに「Decomposed Case」を選べば全体が見える。関数オブジェクトの出力はpostProcessingフォルダに1つにまとまって書かれるので、数値の取り出しはこちらを使うほうが確実だ。
よくある間違い
間違いと対策
| 間違い | 起きること | 対策 |
|---|---|---|
| 2次元でArefに厚さを入れない | 係数が桁違い | 直径×メッシュ厚さ |
| 運動学的圧力をPaと誤解 | 圧力・力が密度倍ずれる | rhoInfを指定 |
| 過渡を含めて平均 | 平均値がずれる | 周期定常後に整数周期 |
| 出力間隔が粗い | 周波数を誤る | 周期の1/20以下 |
| 力の方向ベクトルの誤り | 揚力と抗力が入れ替わる | liftDir・dragDirを確認 |
関連する内容も知りたいです。
関連するページはOpenFOAM入門、OpenFOAMメッシュ生成、OpenFOAM並列計算、OpenFOAM乱流モデルが参考になる。
よくある質問
OpenFOAMで力や係数はどう求めますか？
system/controlDictのfunctions（関数オブジェクト）にforcesやforceCoeffsを書き、計算と同時に出力させます。forceCoeffsでは基準速度magUInf、基準長さlRef、基準面積Aref、流体密度rhoInf、力の方向（liftDir・dragDir）を指定します。計算後にpostProcess -funcで既存の結果から求めることもできます。
2次元計算で抗力係数が極端に小さく出ます。
参照面積Arefにメッシュの厚さを入れていない可能性が高いです。OpenFOAMの2次元計算は厚さ1セルの3次元メッシュなので、力はその厚さ分しか働きません。直径0.01 mの円柱を厚さ0.002 mのメッシュで解いたなら、Aref = 0.01 × 0.002 m² です。Aref = 0.01と入れると、抗力係数1.33が0.00266と表示されます。
非圧縮ソルバーの圧力の値がおかしいのですが？
simpleFoamやpimpleFoamなどの非圧縮ソルバーが解く p は、圧力を密度で割った「運動学的圧力」（単位 m²/s²）です。空気（密度1.2 kg/m³）なら、p = 50 m²/s² は60 Paに相当します。forcesでもrho rhoInfとrhoInfの値を与えないと、力が密度1として計算されます。
参考資料・検算の位置づけ
掲載値は本文の仮定に基づく教育用の参考計算。資料は理論・設定の参照先であり、本記事の数値を各ソフトで再現した証明ではない。設計・製造の判断には対象条件での検証が必要。
関連トピック
なった
詳しく
報告