こんにちは（@t_kun_kamakiri）
今回は境界条件の設定ファイルに直接C++のコードを埋め込んで初期化設定を行う方法について解説をします。C++と言ってもOpenFOAMで用意されている関数を上手く使って書くため特別すごいアルゴリズムを組むわけではありません。
OpenFOAMのカスタマイズにはC++の知識に加えてOpenFOAMのAPIを理解する必要があり、初心者にとってはハードルが高いでしょう。
そこで、OpenFOAMのカスタマイズはしてみたいけど難しそうと感じている人は、まずは本記事のような直接コードを書きこむcodeStreamを試してみると良いです。
- 境界条件をcodeStremを使ってコードを直接埋め込む
- 円筒内流れの流入境界条件を関数で与える
前回の記事で「円筒内流れ」をOpenFOAMでやってみました。
以下、blockMeshでメッシュを作成して円筒内流れを解析した結果です。
blockMeshでの結果
【OpenFOAM(円筒内の流れ)】blockMeshでメッシュ作成(1)
層流条件においては、流入境界での流速は1.0m/sの一定値ですが助走区間を経て徐々に2次関数で近似できるようになり、OpenFOAMの結果が２次関数と一致していることがわかります。
u=2u_{b}\bigg(1-\big(\frac{r}{R}\big)^2\bigg)\tag{1}
\end{align*}
では、はじめから境界条件を(1)で与えるようにするにはどうすれば良いのか？というのが本記事の主題です。
OpenFOAMv2012(WSL2)
モデルファイルの入手
ベースとなるモデルファイルは以下からダウンロードできます。
こちらのモデルを使って境界条件のファイル設定を書き換えます。
具体的には「0/U」の境界面の名前「inlet」の条件をcodeStreamでコードを書いて設定を行います。
境界条件の設定
codeStream
0/Uの境界条件を2次関数で与えます。
0/U
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 | dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 1.0); boundaryField { inlet { type fixedValue; // value uniform (0 0 1.0); value #codeStream { codeInclude #{ #include "fvCFD.H" #}; codeOptions #{ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude #}; code #{ const IOdictionary& d = static_cast<const IOdictionary&>( dict.parent().parent() ); const fvMesh& mesh = refCast<const fvMesh>(d.db()); const label id = mesh.boundary().findPatchID("inlet");// パッチ面に対応したidを取得 const fvPatch& patch = mesh.boundary()[id]; // idからパッチ情報を取得 vectorField Uz_inlet(patch.size(), vector(0, 0, 0)); //Uの初期化(0 0 0) const scalar ub = 1.0; // 平均速度 const scalar Radius = 0.01; //半径 forAll(Uz_inlet, i) { /* 1目の[]はid 2目の[]は座標0->x, 1->y, 2->z const scalar x = mesh.C()[i][0]; const scalar y = mesh.C()[i][1]; でも可 */ const scalar x = patch.Cf()[i][0]; const scalar y = patch.Cf()[i][1]; const scalar r = sqrt(pow(x,2) + pow(y,2)); Uz_inlet[i] = vector(0.0, 0.0, 2.0*ub*(1-pow(r/Radius,2) )); } Uz_inlet.writeEntry("", os); #}; }; } outlet { type zeroGradient; } sideWall { type noSlip; } } |
fvMeshがメッシュを扱うクラスでメッシュに関する情報である「セル、面、エッジ、頂点」などを保持しています。
具体的にはfvMeshは、
- polyMesh
- lduMesh,
- fvSchemes,
- surfaceInterpolation, // needs input from fvSchemes
- fvSolution,
- data
などを基底クラスに持つ多重に継承したクラスです。
fvMeshはtypedefで大文字から始まるMeshに定義されていますが、ややこしいですが同じクラスのことですね。
1 | typedef fvMesh Mesh |
fvMeshが以下のようにboundary()という関数を持っています。
1 | Foam::fvBoundaryMesh & boundary( )const |
なので、
1 | mesh.boundary() |
とするとfvBoundaryMeshクラスが返ってきて、fvBoundaryMeshクラスの中のfindPatchID()関数を呼び出すことでパッチ面のidを取得することができます。
1 | Foam::label findPatchID (const word & patchName) const |
labelはOpenFOAMで定義された整数型のことです。
今回は境界面としてinletを探すようにしているため、findPacthID(“inlet”)としています。
あとはvecterFieldのUz_inletがリストになっているのでfor文で回して、以下の2次関数ををUz_inletに与えれば良いです。
u=2u_{b}\bigg(1-\big(\frac{r}{R}\big)^2\bigg)\tag{1}
\end{align*}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | forAll(Uz_inlet, i) { /* 1目の[]はid 2目の[]は座標0->x, 1->y, 2->z const scalar x = mesh.C()[i][0]; const scalar y = mesh.C()[i][1]; でも可 */ const scalar x = patch.Cf()[i][0]; const scalar y = patch.Cf()[i][1]; const scalar r = sqrt(pow(x,2) + pow(y,2)); Uz_inlet[i] = vector(0.0, 0.0, 2.0*ub*(1-pow(r/Radius,2) )); } |
「fvMesh」のクラスを使うことでメッシュ情報はだいたい取得できます。
以下の表を参考に関数を呼び出し必要な情報を取得します。
| Class | Description | Symbol | Access function |
| volScalarField | Cell volumes | $V$ | V() |
| surfaceVectorField | Face area vector | $\boldsymbol{S_{f}}$ | Sf() |
| surfaceScalarField | Face area magnitude | |$\boldsymbol{S_{f}}$ | | magSf() |
| volVectorField | Cell centres | $\boldsymbol{C}$ | C() |
| surfaceVectorField | Face centres | $\boldsymbol{C_{f}}$ | Cf() |
| surfaceScalarField | Face fluxes | $\boldsymbol{\phi_{f}}$ | Phi() |
今回のように、
1 2 | const scalar x = mesh.C()[i][0]; const scalar y = mesh.C()[i][1]; |
とすればセル中心の座標を取得することがでいます。
以上が全体の流れです。
では、計算を実行させてみましょう。
1 | simpleFoam > log.simpleFoam & |
結果はParaViewで確認します。流入境界条件が思った通りの速度分布になっているかを確認しましょう。
※まず、いきなりこのようなコードを書くのはしんどいのでcodeStreamが使わているチュートリアルを探します。
1 | find $FOAM_TUTORIALS -type f |xargs grep "codeStream" |
いくつか候補が出てきますので適当なファイルを開いて眺めてみると勉強になります。
1 2 3 4 5 6 7 8 | /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interFoam/RAS/floatingObject/constant/dynamicMeshDict.sixDoF: momentOfInertia #codeStream /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/overInterDyMFoam/floatingBody/background/constant/dynamicMeshDict: momentOfInertia #codeStream /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/overInterDyMFoam/floatingBodyWithSpring/background/constant/dynamicMeshDict: momentOfInertia #codeStream /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/pimpleFoam/LES/periodicHill/steadyState/system/blockMeshDict:edges #codeStream /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/pimpleFoam/LES/periodicHill/transient/system/blockMeshDict:edges #codeStream /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/blockMeshDict.main:vertices #codeStream /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/simpleFoam/bump2D/system/blockMeshDict:edges #codeStream /opt/OpenFOAM/OpenFOAM-v2012/tutorials/basic/potentialFoam/cylinder/system/blockMeshDict:vertices #codeStream |
codedFixedValue
境界条件の設定は、
で行うこともできます。
0/U
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 1.0); boundaryField { inlet { type codedFixedValue; value uniform (0 0 0); name inlet_BC; code #{ /* fvPatchField<vector> f //fvPatchVectorFieldでも良い ( patch().lookupPatchField<volVectorField, vector>("U") ); */ const fvPatch& boundaryPatch = patch(); //パッチ情報の取得 const vectorField& Cf = boundaryPatch.Cf(); //パッチ面の中心情報を取得 vectorField& field = *this; // 流速valueとnameを保持 const scalar ub = 1.0; // 平均速度 const scalar Radius = 0.01; //半径 forAll(Cf, i) { /* 1目の[]はid 2目の[]は座標0->x, 1->y, 2->z */ const scalar x = Cf[i][0]; const scalar y = Cf[i][1]; const scalar r = sqrt(pow(x,2) + pow(y,2)); field[i] = vector(0.0, 0.0, 2.0*ub*(1-pow(r/Radius,2) )); } #}; } outlet { type zeroGradient; } sideWall { type noSlip; } } |
今回はいきなりfvPathクラスを「boundaryPatch = patch()」と置き換えて以下のCf()関数を呼び出して、パッチ面の情報を取得しています。
1 | const Foam::vectorField & Cf() const |
あとの流れはcodeStreamと同じです。
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍として重宝しているわかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの？」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ（NextPublishing）)
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
OpenFOAMコードの辞書的な扱いとしては以下の参考書が大変役に立ちます。
本記事ではESI版のOpenFOAMを使っているため本書のFoundation版で対応していない部分がありますが、その辺を考慮しても持っていて損はないでしょう。
Follow @t_kun_kamakiri