こんにちは（@t_kun_kamakiri）
今回は境界条件の設定ファイルに直接C++のコードを埋め込んで初期化設定を行う方法について解説をします。C++と言ってもOpenFOAMで用意されている関数を上手く使って書くため特別すごいアルゴリズムを組むわけではありません。
前回の記事では境界条件に関数を与えるためにC++のコードを書きました。
今回は境界条件ではなく内部の領域に関数で条件を与える方法を解説します。
以下のように$(x-x_{0})^2+(y-y_{0})^2\leq r^2$を液体領域（$\alpha=1$)とするにはどうするのかというのが本記事の主題です。
OpenFOAMv2012(WSL2)
チュートリアルのコピー
今回は気液混相流ソルバのinterFoamのdamBreakのチュートリルをコピーします。
1 | cp -r $FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreak/damBreak damBreak_alpha |
コピーしてきたチュートリアルを編集します。
具体的にはメッシュのスケールを「scale 1;」に変更します。
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 68 69 70 71 72 73 74 75 76 77 78 | scale 1; vertices ( (0 0 0) (1 0 0) (1 1 0) (0 1 0) (0 0 1) (1 0 1) (1 1 1) (0 1 1) ); blocks ( hex (0 1 2 3 4 5 6 7) (50 50 50) simpleGrading (1 1 1) ); edges ( ); boundary ( leftWall { type wall; faces ( (0 4 7 3) ); } rightWall { type wall; faces ( (5 6 2 1) ); } frontWall { type wall; faces ( (4 7 6 5) ); } backWall { type wall; faces ( (3 2 1 0) ); } bottomWall { type wall; faces ( (0 1 5 4) ); } top { type patch; faces ( (7 3 2 6) ); } ); mergePatchPairs ( ); |
これでメッシュを作ると以下のようになります。
1 | blockMesh |
blockMeshを実行して、ParaViewで確認します。
damBreakのチュートリアルは液体領域（$\alpha=1$)はsetFieldsで作るのですが、以下のように自身で好きなように配置させたい場合は難しかったりします。
$(x-x_{0})^2+(y-y_{0})^2\leq r^2$を液体領域（$\alpha=1$)とするにはどうするのか？
内部領域をcodeStreamで与える
では、「0/alpha」を以下のように編集します。
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 68 69 70 | dimensions [0 0 0 0 0 0 0]; //internalField uniform 0; internalField #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); const fvMesh& mesh = refCast<const fvMesh>(d.db()); scalarField alpha(mesh.nCells(), 0.); const scalar x0 = 2.0; const scalar y0 = 3.0; forAll(alpha, i) { const scalar x = mesh.C()[i][0]; const scalar y = mesh.C()[i][1]; const scalar r = sqrt(pow(x - x0,2) + pow(y - y0,2)); // 以下の関数でハートを作ることができる //const scalar rhart = sqrt(pow(x - x0,2) + pow((y-y0) - pow(x-x0, 2/3),2.0)); if (r <= 0.5) { alpha[i] = 1.0; } } alpha.writeEntry("", os); #}; }; boundaryField { leftWall { type zeroGradient; } rightWall { type zeroGradient; } lowerWall { type zeroGradient; } atmosphere { type inletOutlet; inletValue uniform 0; value uniform 0; } defaultFaces { type empty; } } |
一番初めに「fvCFD.H」をインクルードしています。
fvCFD.Hをインクルードすることで以下のクラスを全て継承していることになります。
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 | #ifndef fvCFD_H #define fvCFD_H #include "Time.H" #include "fvMesh.H" #include "fvc.H" #include "fvMatrices.H" #include "fvm.H" #include "linear.H" #include "uniformDimensionedFields.H" #include "calculatedFvPatchFields.H" #include "extrapolatedCalculatedFvPatchFields.H" #include "fixedValueFvPatchFields.H" #include "zeroGradientFvPatchFields.H" #include "fixedFluxPressureFvPatchScalarField.H" #include "constrainHbyA.H" #include "constrainPressure.H" #include "adjustPhi.H" #include "findRefCell.H" #include "IOMRFZoneList.H" #include "constants.H" #include "gravityMeshObject.H" #include "columnFvMesh.H" #include "OSspecific.H" #include "argList.H" #include "timeSelector.H" #ifndef namespaceFoam #define namespaceFoam using namespace Foam; #endif #endif |
「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つ目の[]がセルのidに対応し、2つ目の[]が座標(x->0,y->1,z->0)に対応します。
では、inerFoamを実行して計算しましょう。
計算をスタートさせると「0/alpha.water」が動的にコンパイルされて内部領域の初期状態を作ってくれます。
1 | interFoam |
ParaViewで確認します。
気体と液体の界面がぼやけることがinterFoam(VOF法）のデメリットですが非常に簡単に気液混相流の解析ができるのが特徴です。
今回は試していませんが、interFoamは界面がぼやけてしまうという問題があり、いくつかの改善点が必要です。
- 界面圧縮項を大きくする
- メッシュを細かくする
- Adaptive mesh refinement（界面のメッシュのみ細かく）
- interIsoFoamソルバを使う
interIsoFoam
気液混相流に対してisoAdvector位相分率ベースの界面捕捉アプローチを用いたinterFoamから派生したソルバーで、適応的再メッシュを含むメッシュの動きやメッシュトポロジー変更を任意に行うことができます。
interIsoFoamを使うと界面のぼやけるのがinterFoamと比べるとかなり改善されています。
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍として重宝しているわかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの？」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ（NextPublishing）)
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
OpenFOAMコードの辞書的な扱いとしては以下の参考書が大変役に立ちます。
本記事ではESI版のOpenFOAMを使っているため本書のFoundation版で対応していない部分がありますが、その辺を考慮しても持っていて損はないでしょう。
Follow @t_kun_kamakiri 