OpenFOAM入門：平面ポアズイユ流れのシミュレーション
1. はじめに
OpenFOAMは、流体などの数値解析のシミュレーションおよび前処理・後処理を行うためのOSSのツールボックスです。
この記事では、理論解が求められる平面ポアズイユ流れをOpenFOAMでシミュレーションして、結果を可視化・比較するまでのチュートリアルを紹介します。
流体力学の数値解析だけでも門外漢には非常に難しいのに加え、OpenFOAMには入門者用の親切なドキュメントがほとんどありません。筆者は現在も初心者ですが、これからOpenFOAMで解析を始める方が少しでも簡単に入門できるようにこの記事を作成しました。読者の方が解析を少し自分で修正できるレベルくらいにまでなれるよう、丁寧に説明していきたいと思います。
対象読者
- 流体（動）力学の入門書レベルの知識は何となくはわかる方
- 流体の数値解析を行う必要があるものの、ほとんど何一つわからない方
- OpenFOAMを触ってみた・触ろうとしてみたけどほとんどわからなかった方
- Linux・Pythonの基礎は問題ない方
なお、筆者自身が上記の対象読者と同じレベルですので、内容に何か間違いがあればご指摘いただければ幸いです。
検証環境
この記事で紹介するチュートリアルは、OpenFOAM 10で動作検証を行いました。同じ意味の設定でもバージョンによって書式が違う場合がありますので、このチュートリアルの実行にはOpenFOAM 10を利用してください。
なお、OpenFOAMのバージョン表記は、OpenFOAM 10などのメジャーバージョンと、v2212などといったYYMM形式のバージョンが混在しています。履歴については以下のリストの一番上のページなどをご覧ください。なお、2,3番目の公式系のサイトのバージョン履歴ページには片方ずつしか載っていません。
また、筆者はUbuntu 22.04 LTSでOpenFOAMを利用しています。Ubuntuでのインストールは以下の通りです。
$ sudo sh -c "wget -O - https://dl.openfoam.org/gpg.key > /etc/apt/trusted.gpg.d/openfoam.asc"
$ sudo add-apt-repository http://dl.openfoam.org/ubuntu
$ sudo apt update
$ sudo apt install openfoam10
APTでGPGキーを登録してupdate,installでインストールできます。
その他の環境での利用については、OpenFOAM 10のリリースページの指示に従ってインストールをお願いします。環境依存のインストールや不具合については流石にこの記事のカバー外ですが、他の環境ではDocker上のUbuntuまたはWSLで実行されるようになっていますので、インストール以後の操作はどの環境でも問題なく実行できるのではないかと思います。
2. 数値解析の流れ
大まかに言えば、数値解析の流れは以下の通りです。
- 解析対象の形状（＝メッシュ）を定義する。
- シミュレーションの初期条件・境界条件やソルバのパラメータを設定する。
- シミュレーションを実行する。
- シミュレーションの結果を確認する。
実際の解析ではデバッグや調整で1から4を繰り返したり行ったり来たりすることになりますが、最小限としてはこの四つを上から順に実行するのみです。以降の節で一つずつ説明していきます。
なお、この記事の以降の節では、OpenFOAMの一般的なチュートリアルのように公式のサンプルの全ファイルをコピーしてきて一部を紹介・改変するのではなく、空のディレクトリからスタートして、ユーザが作成する必要のあるファイルの全文を一つずつ定義していきます。また、コマンドによって自動生成されるファイルについては、それとわかるように説明していきます。
そのため、まずは以下のように、作業用のディレクトリを作成してください。
$ mkdir openfoam_sample1
$ cd openfoam_sample1
3. メッシュの作成
まず最初の作業ステップとしては、解析対象のメッシュを定義します。
メッシュとは、解析対象の物体や空間をポリゴンによってモデリングしたものです。
コンピュータでは（任意の形状の）連続空間を取り扱うことができないため、数値解析では多くの場合、解析対象をポリゴンによって近似的に表現します。圧力や流速などの数値はメッシュの頂点でのみ定義されます。また、ナビエ–ストークス方程式などの微分方程式を計算する際は、ある頂点の値の更新については辺で接続された隣接頂点の値（＋自身の値や全体設定のパラメータなど）のみが使用されます。
blockMeshによる格子メッシュの定義
メッシュを作成する方法は色々ありますが、今回はblockMeshを使用します。その他の方法についてはまた別の機会に紹介するかもしれません。
blockMeshはOpenFOAMに同梱されているツールの一つです。格子状の単純なメッシュを定義するためのツールで、テキストファイルでメッシュの形状を定義します。
まず、systemディレクトリを作成してください。主にシミュレーション全体の設定ファイルを格納する固定名のディレクトリです。blockMeshの定義ファイルもこの中に格納します。
$ mkdir system
次に、以下のファイルをcontrolDictというファイル名でsystem以下に保存してください。
// ファイルの種類・概要。
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// 各種定義。
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.5;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
controlDictはシミュレーション全体の設定を記述するためのパス・名称固定のファイルです。
このファイルの中身については詳細には説明しませんが、一般的なシミュレーションでよく出てくる概念のパラメータがほとんどなので、何となくイメージできるのではないかと思います。
blockMeshコマンドの実行時には、このcontrolDictが存在しないとエラーになります。deltaTなどのシミュレーションの関連の値は無関係ですが、writeCompressionなどのファイル形式の設定が参照されているようです。
最後に、以下のファイルをblockMeshDictという名称でsystem以下に保存してください。
// ファイルの種類・概要。
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// 定義ファイル内の座標をメッシュに変換する際のスケール。
convertToMeters 0.001; // 頂点をmm単位で定義。
// 格子の頂点の定義。
// 横300mm, 縦100mm, 厚さ10mm。
// 一次元流れなので厚さは不要だが、blockMeshの構造上必要。
vertices
(
// 右手系の(x y z)
(0 0 0)
(300 0 0)
(300 100 0)
(0 100 0)
(0 0 10)
(300 0 10)
(300 100 10)
(0 100 10)
);
// 格子の形状の定義。
// hex (六面体の頂点リスト) (x y zの分割数) simpleGrading (分割の重みづけ)
blocks
(
hex (0 1 2 3 4 5 6 7) (60 20 1) simpleGrading (1 1 1)
);
// 境界面の定義
boundary
(
// 流入口（x正を右、y正を上、z正から負に向かって見下ろした時の、左面）
inlet
{
type inlet; // 面の種類
faces // どの面をこの境界として定義するか。
(
(0 4 7 3)
);
}
// 流出口（右面）
outlet
{
type outlet;
faces
(
(2 6 5 1)
);
}
// 管路壁面（上下）
fixedWalls
{
type wall;
faces
(
(3 7 6 2)
(1 5 4 0)
);
}
// 厚み部分（前後）
frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(4 5 6 7)
);
}
);
blockMeshDictはblockMeshによって生成する格子メッシュを定義するための、パス・名称固定のファイルです。blockMeshを使用するためには、固定の辞書構造のデータを上記のように定義します。各項目の説明はファイル内のコメントをご覧ください。
この状態で、openfoam_sample1をカレントディレクトリとしてblockMeshコマンドを実行すると、メッシュが生成されます。生成されたメッシュデータはconstant/polyMesh以下に複数のファイルに分けて保存されます。
$ blockMesh
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 10-c4cf895ad8fa
Exec : blockMesh
(途中省略)
Writing polyMesh
----------------
Mesh Information
----------------
boundingBox: (0 0 0) (0.3 0.1 0.01)
nPoints: 2562
nCells: 1200
nFaces: 4880
nInternalFaces: 2320
----------------
Patches
----------------
patch 0 (start: 2320 size: 20) name: inlet
patch 1 (start: 2340 size: 20) name: outlet
patch 2 (start: 2360 size: 120) name: fixedWalls
patch 3 (start: 2480 size: 2400) name: frontAndBack
End
$ ls constant/polyMesh
boundary faces neighbour owner points
エラーになった場合は、エラーメッセージの内容を確認してみてください。
メッシュの可視化・確認
メッシュの可視化にはparaViewを使用します。paraViewもまたOpenFOAMに同梱されているツールの一つです。
paraViewは、正確には、メッシュおよびシミュレーション結果の可視化ツールです。そのため、メッシュデータ単体としては問題がなくても、シミュレーションデータや設定ファイルに不整合があるとうまく動きません。
今回のチュートリアルでも、可視化の前にまず、system/fvSchemesおよびsystem/fvSolutionを作成します。
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSolution;
}
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
fvSchemeは離散化スキームの定義、fvSolutionはシミュレーションのソルバの定義を行うファイルです。このファイルの中身については、この記事では省略させてください。
この状態で、paraFoamコマンドを実行すると、paraViewのGUIアプリケーションが立ち上がります。paraFoamとparaViewがややこしいのでご注意ください。
$ paraFoam
最初は何も表示されていない状態で正常です。左メニューの"Apply"を押下するとメッシュが表示されます。
画面上部メニューの"Surface"を"Wireframe"など他の項目に切り替えると、表示の仕方が切り替わります。
一般的な3DCGソフトと同様に、メッシュはマウスで回転や拡大縮小ができます。
なお、paraViewが正常に動作しない原因としてよくあるものは、blockMeshDictのboundaryと、後で紹介する0/Uや0/pのboundaryFieldが整合していない場合です。これらのファイルはシミュレーションの境界条件を定義するためのファイルです。この問題は、既存のプロジェクトをコピーしてきて、メッシュの改変だけが完了した段階で可視化する場合に起こりがちです。
解決の方法は、0/Uなどを一旦削除してしまうか、または左メニューの"Fields"のチェックを外してから"Apply"を押下するかのどちらかです。
4. シミュレーションの設定
次に、シミュレーションを実行するための設定を作っていきます。
まず前提として、OpenFOAMには複数のソルバ（シミュレータ）が用意されています。
| ソルバ | 圧縮・非圧縮 | 定常・非定常 | 層流・乱流 |
|---|---|---|---|
icoFoam |
非圧縮性 | 非定常 | 層流 |
simpleFoam |
非圧縮性 | 定常 | 乱流 |
pisoFoam |
非圧縮性 | 非定常 | 乱流 |
pimpleFoam |
非圧縮性 | 非定常 | 乱流 |
その他については以下をご覧ください。
今回はicoFoamを使用します。ポアズイユ流れは定常流ですが、時間発展で安定した状態を見てみたいと思います。
境界条件の設定
ソルバごとに必要な設定は異なりますが、まず、多くのソルバで共通となる、流速と圧力について設定を作成します。
初期条件は、初期時刻のディレクトリを作成し、その中に定義していきます。具体的にはsystem/controlDictのstartTime、すなわち0です。
$ mkdir 0
流速は0/U、圧力は0/pに定義します。以下をファイルに保存してください。
// ファイルの種類・概要。
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// [kg m s K mol A Cd] で何の成分を定義しているかを記述。
// [m/s]
dimensions [0 1 -1 0 0 0 0];
// 内点の初期値。
internalField uniform (0 0 0);
// 境界の初期値。
boundaryField
{
inlet
{
type zeroGradient; // 勾配0。
}
outlet
{
type zeroGradient; // 勾配0。
}
fixedWalls
{
type noSlip; // 滑らない。
}
frontAndBack
{
type empty;
}
}
// ファイルの種類・概要。
FoamFile
{
format ascii;
class volScalarField;
object p;
}
// [kg m s K mol A Cd] で何の成分を定義しているかを記述。
// [m^2/s^2]。非圧縮性の場合、圧力は密度で割った単位になるとのこと。
dimensions [0 2 -2 0 0 0 0];
// 内点の初期値。
internalField uniform 0;
// 境界の初期値。
boundaryField
{
inlet
{
type fixedValue; // 常に定数値
value uniform 1; // 1 [m^2/s^2]
}
outlet
{
type fixedValue; // 常に定数値
value uniform 0; // 1 [m^2/s^2]
}
fixedWalls
{
type zeroGradient; // 勾配0。
}
frontAndBack
{
type empty; // 勾配0。
}
}
流速は初め全て0.0にしておきます。圧力は入口と出口に差をつけて、密度あたりの圧力勾配が（値は幾つでも良いのですが）1になるように設定しておきます。
その他物理パラメータの設定
その他、icoFoamの実行には、constant/physicalPropertiesが必要です。以下をファイルに保存してください。
FoamFile
{
format ascii;
class dictionary;
location "constant";
object physicalProperties;
}
nu [0 2 -1 0 0 0 0] 0.01; // 動粘度 0.01 [m^2/s]
rho [ 1 -3 0 0 0 0 0 ] 1; // 密度 1 [kg/m^3]
内容としては、流体の動粘度と密度を設定しているだけです。なお、ファイルがないとicoFoamの実行時にエラーになりますが、設定については未定義の場合にはデフォルト値が用いられる（ものがある）ようです。
5. シミュレーションの実行
ここまでの設定がうまくいっていれば、あとはコマンド一発でシミュレーションが実行できます。
$ icoFoam
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 10-c4cf895ad8fa
Exec : icoFoam
...
Time = 0.5s
Courant Number mean: 0.276748 max: 0.412884
smoothSolver: Solving for Ux, Initial residual = 0.000271161, Final residual = 6.34569e-06, No Iterations 8
smoothSolver: Solving for Uy, Initial residual = 0.516283, Final residual = 8.18554e-06, No Iterations 23
DICPCG: Solving for p, Initial residual = 4.1799e-06, Final residual = 8.56715e-07, No Iterations 7
time step continuity errors : sum local = 7.98172e-09, global = 2.80603e-09, cumulative = -2.33437e-05
DICPCG: Solving for p, Initial residual = 3.57872e-06, Final residual = 9.83001e-07, No Iterations 6
time step continuity errors : sum local = 9.15828e-09, global = 1.89267e-09, cumulative = -2.33418e-05
ExecutionTime = 0.12316 s ClockTime = 1 s
End
system/controlDictに設定したendTimeまでの間、deltaTとwriteIntervalの間隔で中間状態が保存されていることが確認できます。
$ ls
0 0.1 0.2 0.3 0.4 0.5 constant system
$ ls 0.1
U p phi uniform
6. 結果の可視化
結果の可視化は前述の通りparaViewを使用します。
$ paraFoam
今回は、"Apply"を押下する前に"Last Frame"ボタンを押下して、全てのシミュレーション結果を読み込んでください。
可視化の前に結果を読み込むことにより、カラーマップの範囲などが適切に自動調整されます。
また、シミュレーション結果に何か問題がある場合は、この時点でparaViewが落ちたりします。
"Apply"を押下したあとは、左上メニューで"vtkBlockColors"を"U"や"p"に変更してみてください。
流速および圧力の分布が確認できます。
7. 理論解との比較
それでは最後に、シミュレーションの結果と理論解を比較してみます。
平面ポアズイユ流れの理論的な説明については、例えば以下のページを参照してください。
まず、理論式としては、密度を
となります。
この式に、physicalPropertiesで設定したblockMeshDictで設定した
となります。
一方、paraViewの"U"の結果を見てみると、管の中央で最大値4.1e-1になっていることがわかります。
OpenFOAMの理解を深める練習として、各種設定ファイルの数値を修正して再実行した場合に、理論解と結果が一致するかどうかを確かめてみるのも良いかもしれません。
8. その他入門向け資料の紹介
日本語の入門向け資料としては、以下のサイトのシリーズが参考になると思います。それぞれ丁寧に書かれていて、筆者もお世話になっています。
Discussion