こんにちは（@t_kun_kamakiri）
今回は境界条件の設定ファイルに直接C++のコードを埋め込んで初期化設定を行う方法について解説をします。C++と言ってもOpenFOAMで用意されている関数を上手く使って書くため特別すごいアルゴリズムを組むわけではありません。
前回の記事では境界条件に関数を与えるためにC++のコードを書きました。
今回は境界条件の特定領域に液相率を与え、その領域に流速$u_{x}$を与える方法を解説します。
液体領域（$\alpha=1$)とするにはどうするのかというのが本記事の主題です。
OpenFOAMv2012(WSL2)
チュートリアルのコピー
今回は気液混相流ソルバのinterFoamのdamBreakのチュートリルをコピーします。
1 | cp -r $FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreak/damBreak damBreak_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 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で作るのですが、以下のように自身で好きなように配置させたい場合は難しかったりします。
$\alpha =1, \boldsymbol{u}=(1.0,0,0)\,(0.5\leq y\leq 0.7\,\, , 0.4\leq z\leq 0.6) $
境界の特定領域をcodeStreamで与える
では、「0/alpha」を以下のように編集します。
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 71 72 73 74 75 76 77 78 79 80 | dimensions [0 0 0 0 0 0 0]; internalField uniform 0; boundaryField { leftWall { type codedFixedValue; value uniform 0; name inletProfile2; code #{ const fvPatch& boundaryPatch = patch(); const vectorField& Cf = boundaryPatch.Cf(); scalarField& field = *this; /* //take value from initialization at the given boundary surface */ field = patchInternalField(); scalar minz = 0.4; scalar maxz = 0.6; scalar miny = 0.5; scalar maxy = 0.7; scalar t = this->db().time().value(); forAll(Cf, faceI) { if ( (Cf[faceI].z() > minz) && (Cf[faceI].z() < maxz) && (Cf[faceI].y() > miny) && (Cf[faceI].y() < maxy) ) { if ( t < 1.) { field[faceI] = 1.; } else { field[faceI] = 0.; } } } #}; } rightWall { type zeroGradient; } frontWall { type zeroGradient; } backWall { type zeroGradient; } bottomWall { type zeroGradient; } top { type inletOutlet; inletValue uniform 0; value uniform 0; } } |
※field = patchInternalField(); の記述いるかな？ってあとで思った・・・
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 68 69 70 71 72 73 74 75 76 77 78 79 | dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { leftWall { type codedFixedValue; value uniform (0 0 0); name inletProfile1; code #{ const fvPatch& boundaryPatch = patch(); const vectorField& Cf = boundaryPatch.Cf(); vectorField& field = *this; scalar minz = 0.4; scalar maxz = 0.6; scalar miny = 0.5; scalar maxy = 0.7; scalar t = this->db().time().value(); forAll(Cf, faceI) { if ( (Cf[faceI].z() > minz) && (Cf[faceI].z() < maxz) && (Cf[faceI].y() > miny) && (Cf[faceI].y() < maxy) ) { if ( t < 1.) { field[faceI] = vector(1,0,0); } else { field[faceI] = vector(0,0,0); } } } #}; } rightWall { type fixedValue; value uniform (0 0 0); } bottomWall { type fixedValue; value uniform (0 0 0); } frontWall { type fixedValue; value uniform (0 0 0); } backWall { type fixedValue; value uniform (0 0 0); } top { type pressureInletOutletVelocity; value uniform (0 0 0); } } |
以下の表を参考に関数を呼び出し必要な情報を取得します。
| 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 3 4 5 6 7 | if ( (Cf[faceI].z() > minz) && (Cf[faceI].z() < maxz) && (Cf[faceI].y() > miny) && (Cf[faceI].y() < maxy) ) {処理} |
とすればセル中心の座標によって条件を変えることができます。
では、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 