よくわからないOpenFOAM
By S.Negishi
プロローグ ― この記事を読む人たちへ
この記事では、OpenFOAMの経験者向けに、OpenFOAMのよくわからない機能や現象について紹介している。
今回は、OpenFOAMをよく知らない人にも雰囲気だけでも分かるようにどのような人向けに、どのようなものか少し具体的に説明しようと思う。
OpenFOAM (Open source Field Operation And Manipulation)は数値解析、 主にCFD (Computer Field Dynamics, 数値流体力学)用のC++で作ったツールボックス（アプリケーションの塊）である。
OpenFOAMはオープンソースソフトウェアであり、 通常、無料でソースコードが公開され、利用することができる。
無料で流体解析ができることはOpenFOAMを利用する上での長所である。
ただし、オープンソースソフトウェアの性質上、留意しなければならないことがある。
-
一つ目はこのソフトウェアは無保証であるということだ。 もしこのソフトに間違い(不具合)があったとしても、修正がすぐにされるとは限らないし、 基本的にはよほどのことがない限り損失の補償もされないだろう。
-
二つ目はGPL (General Public License)による制約である。 改変したソースコードを自分で利用する分には公開する必要はない。 しかし、改変したソースコードを利用して作ったソフトウェアは、原則提供先にソースコードを公開しなければならない。
-
三つ目はOpenFOAMソルバーの詳細な分かりやすいマニュアルがあるとは限らないことである。 マニュアルがない場合は自分でインターネットや参考書、ソースコードを調べるなどをして、自力で利用情報を分かりやすく整理する必要がある。
OpenFOAMは非圧縮性または圧縮性流体の解析をはじめ、混相流や遷音速流、超音速流、希薄流、粘弾性流体等の特殊な場合での流体解析も一応可能である。
ただし、ANSYS Fluentなどの商用ソルバーと比べると一般に計算の収束性や信頼性に問題を抱えている場合があるため注意が必要である。
OpenFOAMは無料で利用可能でCFDに関する研究や学習する上で有用である。 ただし、利用できるようにするまでのコストが商用ソフトよりかかってしまう場合もある。
そのことを踏まえたうえで使ってほしい。
実録！揚力係数の災難（その1）
車に似せたベンチマーク用のモデルにAhmed bodyというものがある。
Ahmedという学者が考案したからAhmed body というらしい。
（下図の寸法はmm。参照: Ahmed bodyとか S.R. Ahmed, G. Ramm, Some Salient Features of the Time-Averaged Ground Vehicle Wake, SAE-Paper 840300, 1984）
この全長1m強のAhmed bodyに40m/sの風を流した場合の流体計算をOpenFOAMでやってみる。
まず対象モデルは上図のφ=35°に設定、流体は非圧縮の空気でソルバーはsimpleFoamを使おう。
乱流モデルはSST k-ωにして、揚力係数はfunction objectsのforceCoeffsを使えばよいだろう。
さて、結果を見るか…
forceCoeffs forceCoeffs1 output:
Cm = -0.28220571
Cd = 0.33437225
Cl = -3.0134781
Cl(f) = -1.7889447
Cl(r) = -1.2245333
揚力がマイナスになってる…
これはいったいどういうことなのだろうか？
実録！揚力係数の災難（その2）
OpenFOAMの場合$0$という名のフォルダで初期値（時間若しくは反復回数$＝0$）の設定を行う。
その中のファイル名pを見てみると以下のような記述がある。
FoamFile
{
version 2.3;
format ascii;
class volScalarField;
location "0.org";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 1.0e05;
boundaryField
{
"inlet.*."
{ type zeroGradient; }
"outlet.*."
{ type fixedValue;
value $internalField; }
"frontBackTop.*."
{ type symmetry; }
"body.*."
{ type zeroGradient; }
"bottom.*."
{ type zeroGradient; }
}
dimensionsは圧力の次元である。OpenFOAMの場合、非圧縮流体は密度で圧力値を割っていることに注意してほしい。
internalFieldは解析領域内の初期値、boundaryFieldは各境界条件を記述している。
この中に不具合の原因があるのだが分かるだろうか？
実録！揚力係数の災難（その3）
それでは答えである。
実はinternalFieldの値が不適切だ。以下のように変更する。
変更前:
dimensions [0 2 -2 0 0 0 0];
internalField uniform 1.0e05;
boundaryField
変更後:
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
すると結果はClが正となりもっともらしい値となる。
forceCoeffs forceCoeffs1 output:
Cm = 0.028143459
Cd = 0.33437225
Cl = 0.060316112
Cl(f) = 0.058301515
Cl(r) = 0.0020145975
なぜこのようなことで値が変わるのだろうか？
じつはCd, Cm, Clといった値は各セル位置での(圧力×セル面素ベクトル)の値を定義式に従って計算後、
それをすべて足し合わせた値として出力している。
そのため、圧力の初期値に100000といった大きな値を設定すると計算時に桁落ちが発生し、正しい結果が得られない場合がある。
ANSYS Fluentのような商用ソルバーの場合は対策がなされている場合もあるが、 OpenFOAMのようなフリーソフトの場合、各ユーザが注意と対策を行うことが肝要である。
もし興味があればお問い合わせいただけると、嬉しい。
最後まで、読んでいただき、どうもありがとうございました。
弊社の解析事例
弊社の流体解析事例については、下記のリンクからご覧ください。