こんにちは（@t_kun_kamakiri）
最近DEMの見習い中で、DPM(Discrete Phase Model 分散モデル)と粒子間の衝突の取り扱いを簡略化したMP-PIC法（multiphase particle-in-cell method）を使った解析に取り組んでいます。
MP-PIC法については理解が浅いので参考記事を載せておきます。
- Multiphase particle-in-cell method
- multiphase particle-in-cell method（O’Rourke 1997）
- MP-PIC法による噴流層の解析
- https://www.mr-cfd.com/services/fluent-modules/discrete-phase-model-dpm/
理論的な解説は以下の論文が参考になりそうです。
粒子間相互作用をマクロで近似するMPPIC法。
・packing model
・collisional damping models
・collisional isotropy model
・blended acceleration model
これらがOpenFOAMに実装されていると。https://t.co/3Q7GgWQiM4— カマキリ🐲CAE頑張る (@t_kun_kamakiri) July 14, 2023
今回はOpenFOAMのMP-PIC法のソルバのMPPICFoamで流体の乱流モデルの$k$-$\omega$SSTの解析をしてみます。
設定は少々適当ですので、「やってみた」系の記事なります。
OpenFOAM v2212
チュートリアルのコピー
まずはチュートリアルのコピーをします。
1 | cp -r /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/MPPICFoam/injectionChannel MPPI_injectionChannel |
フォルダ移動します。
1 | cd MPPI_injectionChannel |
とりあえずチュートリアを計算してみましょう。
1 | blockMesh |
メッシュ作成できたらParaViewで確認してみましょう。
1 | MPPICFoam |
これで解析実行がされます。
結果はParaViewで確認することができますが、ここでは割愛します。
粒子の設定の確認
このチュートリアルの肝になるのが粒子の設定の部分です。
constant/kinematicCloudProperties
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 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 | solution { active true; coupled true; transient yes; cellValueSourceCorrection on; interpolationSchemes { rho.air cell; U.air cellPoint; mu.air cell; } averagingMethod basic; integrationSchemes { U Euler; } sourceTerms { schemes { U semiImplicit 1; } } } constantProperties { rho0 1000; alphaMax 0.9; } subModels { particleForces { WenYuDrag { alphac alpha.air; } } injectionModels { /* % matlab/octave calculation of flow rate alpha=0.1 U=20 A=0.01^2 Q=U*A d=650e-6 v=(4/3)*pi*(d/2)^3 n=1 rate=alpha*Q/(v*n) */ model1 { type patchInjection; parcelBasisType fixed; patch lowerInlet; U0 (18.7939 6.8404 0); nParticle 1; parcelsPerSecond 1390885; sizeDistribution { type normal; normalDistribution { mu 650e-6; sigma 25e-6; minValue 500e-6; maxValue 800e-6; } } flowRateProfile constant 1; massTotal 0; SOI 0; duration 60; } model2 { type patchInjection; parcelBasisType fixed; patch upperInlet; U0 (18.7939 -6.8404 0); nParticle 1; parcelsPerSecond 1390885; sizeDistribution { type normal; normalDistribution { mu 650e-6; sigma 25e-6; minValue 500e-6; maxValue 800e-6; } } flowRateProfile constant 1; massTotal 0; SOI 0; duration 60; } } dispersionModel none; patchInteractionModel localInteraction; localInteractionCoeffs { patches ( "(frontAndBack|walls)" { type rebound; e 1; mu 0; } "(lowerInlet|upperInlet|outlet)" { type escape; } ); } surfaceFilmModel none; packingModel none; dampingModel relaxation; relaxationCoeffs { timeScaleModel { type nonEquilibrium; alphaPacked 0.58; e 0.9; } } isotropyModel stochastic; stochasticCoeffs { timeScaleModel { type isotropic; alphaPacked 0.58; e 0.9; } } stochasticCollisionModel none; } cloudFunctions {} // ************************************************************************* // |
設定の主要な部分だけ解説をしておきます。
- particleForces：流体から粒子への力
- injectionModels：patchInjectionでは指定したパッチから粒子を投入
- patchInteractionModel：localInteractionでは指定したパッチ面との接触定義（反発係数や摩擦係数を指定）
- dampingModelとisotropyModel がMP-PICに関係している部分？よくわかっていない。ちなみに$FOAM_TUTORIALS/lagrangian/MPPICFoam/cyclone/constant/kinematicCloudPropertiesを見ると以下のようになっています。1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253subModels{(省略）packingModel implicit;explicitCoeffs{particleStressModel{type HarrisCrighton;alphaPacked 0.6;pSolid 10.0;beta 2.0;eps 1.0e-7;}correctionLimitingMethod{type absolute;e 0.9;}}implicitCoeffs{alphaMin 0.0001;rhoMin 1.0;applyLimiting true;applyGravity false;particleStressModel{type HarrisCrighton;alphaPacked 0.6;pSolid 5.0;beta 2.0;eps 1.0e-2;}}dampingModel none;isotropyModel stochastic;stochasticCoeffs{timeScaleModel{type isotropic;alphaPacked 0.6;e 0.9;}}stochasticCollisionModel none;}
HarrisCrightonは以下の応力の式です。
粒子追跡は以下でも関連記事を書いていますのでご参考ください。
乱流モデルの設定
$k$-$\omega$SST乱流モデルの設定ファイルをコピーします。
lagrangianのチュートリアルの中で乱流モデルの「k」ファイルを使っているケースはないかを探します。
1 | find /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/ -name "k.*" |
候補として1つだけありました。
1 | /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/MPPICFoam/cyclone/0.orig/k.air |
このチュートリアルの乱流モデルは何を使っているのか確認してみましょう。
1 | vi /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/MPPICFoam/cyclone/constant/turbulenceProperties.air |
$FOAM_TUTORIALS/lagrangian/MPPICFoam/cyclone/constant/turbulenceProperties.air
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | simulationType LES; LES { LESModel kEqn; turbulence on; printCoeffs on; delta cubeRootVol; cubeRootVolCoeffs { } } |
どうやらLESのようです。
今回はRANSの$k$-$\omega$SSTにしたいので、少し違いますが参考にはなりそうなので、こちらのチュートリアルを参考にさせてもらいます。
まずはconstant/turbulenceProperties.airを編集します。
constant/turbulenceProperties.air
1 2 3 4 5 6 7 8 | simulationType RAS; RAS { // RASModel kEpsilon; RASModel kOmegaSST; turbulence on; printCoeffs on; } |
境界条件
k.airとnut.airをコピーします。
omega.airは自分で作る必要があります。k.airを適当にコピーして編集することにします。
1 2 3 | cp -r /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/MPPICFoam/cyclone/0.orig/k.air 0/ cp -r /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/MPPICFoam/cyclone/0.orig/nut.air 0/ cp -r 0/k.air 0/omega.air |
では、各ファイルを編集していきます。
0/k.air
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 | FoamFile { version 2.0; format ascii; class volScalarField; object k.air; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 -2 0 0 0 0]; internalField uniform 1; boundaryField { walls { U U.air; type kqRWallFunction; value $internalField; } lowerInlet { type fixedValue; value $internalField; } upperInlet { type fixedValue; value $internalField; } outlet { type inletOutlet; U U.air; phi phi.air; inletValue $internalField; value $internalField; } frontAndBack { type symmetry; } } |
冒頭のヘッダー部分は「k.air」であることを確認します。
outletは「phi phi.air」が必要だそうです。
0/nut.air
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 | FoamFile { version 2.0; format ascii; class volScalarField; object nut.air; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 -1 0 0 0 0]; internalField uniform 0; boundaryField { walls { type nutkWallFunction; value $internalField; } lowerInlet { type calculated; value $internalField; } upperInlet { type calculated; value $internalField; } outlet { type calculated; value $internalField; } frontAndBack { type symmetry; } } |
冒頭のヘッダー部分は「nut.air」であることを確認します。
outletは「phi phi.air」が必要だそうです。
0/omega.air
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 | FoamFile { version 2.0; format ascii; class volScalarField; object omega.air; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 0 -1 0 0 0 0]; internalField uniform 1; boundaryField { walls { k k.air; type kqRWallFunction; value $internalField; } lowerInlet { type fixedValue; value $internalField; } upperInlet { type fixedValue; value $internalField; } outlet { type inletOutlet; k k.air; phi phi.air; inletValue $internalField; value $internalField; } frontAndBack { type symmetry; } } |
冒頭のヘッダー部分は「omega.air」であることを確認します。
outletは「phi phi.air」「k k.air」が必要だそうです
それぞれ.airの記述が必要な理由がわかっていませんが、計算に必要な物理量を書いているということでしょうか・・・0/Uは「phi phi.air」が必要です。
離散化スキームと代数ソルバ
使用する物理量が変わったので離散化スキームと代数ソルバの設定も変更が必要です。
1 2 | cp -r /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/MPPICFoam/cyclone/system/fvSchemes system/ cp -r /usr/lib/openfoam/openfoam2212/tutorials/lagrangian/MPPICFoam/cyclone/system/fvSolution system/ |
system/fvSchemes
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 | ddtSchemes { default Euler; } gradSchemes { default Gauss linear; } divSchemes { default none; div(alphaPhi.air,U.air) Gauss linearUpwindV unlimited; div(((alpha.air*nuEff.air)*dev2(T(grad(U.air))))) Gauss linear; div(phiGByA,kinematicCloud:alpha) Gauss linear; div(alphaPhi.air,omega.air) Gauss limitedLinear 1; div(alphaPhi.air,k.air) Gauss limitedLinear 1; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } wallDist { method meshWave; } |
fvSchemesには壁関数の計算で壁からの距離を必要とするため、「wallDist」を設定する必要があります。
system/fvSolution
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 | solvers { p { solver GAMG; tolerance 1e-06; relTol 0.01; smoother GaussSeidel; } pFinal { solver GAMG; tolerance 1e-06; relTol 0; smoother GaussSeidel; } "(U|k|epsilon|omega).air" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-05; relTol 0.1; } "(U|k|epsilon|omega).airFinal" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-05; relTol 0.1; } kinematicCloud:alpha { solver GAMG; tolerance 1e-06; relTol 0.1; smoother GaussSeidel; } } PIMPLE { nOuterCorrectors 1; nCorrectors 2; momentumPredictor yes; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; } |
計算実行
では設定が終わりましたので、MPPICFoamを実行します。
1 | MPPICFoam |
結果をParaViewで確認してみましょう。
omegaの初期設定が適当なので流速が少し不自然なですね。
今回のMPPICFoamは非定常の解析なので、omegaの値も見積もって適切に設定する必要があります。
まとめ
今回はOpenFOAMをMPPICFoamで乱流モデル$k$-$\omega$SSTを設定してみました。
どういったアルゴリズムで解かれているか詳しく見ていかないといけませんが、とりあえず解析はできました。
わかっていない点
- outletは「phi phi.air」の記述が必要であるという意味は理解できていません。
- 連続相（流体）の流速に乱流場の変動成分を加える分散モデル(dispersionModel)の設定がありますが、どういうときに設定する必要があるのかわかっていません。
・GradientDispersionRAS：乱流エネルギーkの方向に一様乱数で乱れ速度を加える
・StochasticDispersionRAS：連続相の速度に一様乱数で発生させた方向に乱れ速度を加える
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍が無い中唯一わかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの？」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ（NextPublishing）)
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
OpenFOAMコードの辞書的な扱いとしては以下の参考書が大変役に立ちます。
本記事ではESI版のOpenFOAMを使っているため本書のFoundation版で対応していない部分がありますが、その辺を考慮しても持っていて損はないでしょう。
本記事の粒子追跡については結構詳しく書かれていますので、とても参考になります。
メッシュ作成を極めたいって方はこちらもお勧めです。
OpenFOAMでメッシュ作成
Follow @t_kun_kamakiri