こんにちは(@t_kun_kamakiri)
今回はOpenFOAMで以下のような配管内の水の流れを解析し、境界と任意の断面での流量を計算します。
OpenFOAMで境界面での流量を計算する。
境界の「inlet」「outlet1」「outlet2」での流量はOpenFOAMのfunction objectの機能を使うことで簡単に出力することができます。
c1~c5の断面での流量を出力するにはどうしたらいいのかなと思って自分なりに試してみた内容は以下の記事に書いています。
【使用環境】
OpenFOAM-v2012(WSL2)
モデルを作成
まずはモデルを作成する必要があるので簡単に手順を解説します。
- モデル作成と面に名前を付ける:FreeCAD→modelフォルダ
- メッシュ作成:cfMesh→meshフォルダ
- 計算実行:simpleFoam→runフォルダ
それぞれを保存するフォルダを作ります。
1 2 3 |
mkdir model mkdir mesh mkdir run |
モデル作成と面に名前を付ける:FreeCAD
モデルを保存するmodelフォルダに移動します。
1 |
cd model |
今回FreeCADで作成したモデルは以下です。
形状に特に意味はありません。太い円柱と細い円柱があった方が面白いかなと思っただけです。
作成したmodel.stlをmodelフォルダに保存します。
FreeCADで出力されたモデルは寸法は数値しか拾っていないのでm単位に変換する必要があります。
1 |
surfaceConvert -scale 0.001 model.stl model_m.stl |
次に、メッシュ作成の際にcfMeshで使う特徴線を持ったモデルに変換しておきます。
1 |
surfaceFeatureEdges model_m.stl model_m.fms |
メッシュ作成時に使うのは.stlファイルではなくこちらの.fmsファイルです。
メッシュ作成:cfMesh→meshフォルダ
円柱のようなメッシュはsnappyHexMeshよりcfMeshの方が簡単・早い・きれいに作れるのでcfMeshを使います。
以下に過去に検討した記事を書きましたのでご参考ください。
snappyHexMeshは特に境界層がきれいに入らないという問題があります。
meshフォルダに移動します。
1 |
cd mesh |
cfMeshのチュートリアルをコピーしてきて編集します。
1 2 |
cp -r $FOAM_TUTORIALS/modules/cfmesh/tutorials/cartesianMesh/multipleOrifices cfMesh cd cfMesh |
cfMeshは「テトラメッシュ」「ヘキサメインのメッシュ」「ポリヘドラルメッシュ」とあり、チュートリアルが分かれていますが共通のファイルで良く、コマンドが変わるだけなのでチュートリアルのどれを持ってきても良いです。
自分はデフォルトで色々設定が書いている上記のチュートリアルをコピーして使っています。
model作成で作った「model_m.fms」をコピーします。
1 |
cp ../../model/model_m.fms . |
フォルダ構成は以下です。
1 2 3 4 5 6 7 8 9 10 |
. ├── Allclean ├── Allrun ├── model_m.fms ├── multipleOrifices.stl └── system ├── controlDict ├── fvSchemes ├── fvSolution └── meshDict |
system/meshDicr
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 |
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | cfMesh: A library for mesh generation | | \\ / O peration | | | \\ / A nd | Author: Franjo Juretic | | \\/ M anipulation | E-mail: franjo.juretic@c-fields.com | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object meshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // surfaceFile "model_m.fms"; //minCellSize 1.0; maxCellSize 0.01; // boundaryCellSize 2.0; localRefinement { "(sideWall1|sideWall2|inlet|outlet1).*" { cellSize 0.004; } "(sideWall3|outlet2).*" { cellSize 0.002; } } boundaryLayers { // nLayers 3; // thicknessRatio 1.2; // maxFirstLayerThickness 0.5; patchBoundaryLayers { "(sideWall1|sideWall2|inlet|outlet).*" { nLayers 3; thicknessRatio 1.2; maxFirstLayerThickness 0.2; allowDiscontinuity 0; } "(sideWall3|outlet2).*" { nLayers 3; thicknessRatio 1.2; maxFirstLayerThickness 0.2; allowDiscontinuity 0; } } optimiseLayer 1; } // ************************************************************************* // |
以下のコマンドでメッシュを作ります。
1 |
cartesianMesh |
メッシュができれば 「post.foam」など拡張子が.foamという空ファイルを作ってParaViewで読み込んで確認します。
計算実行:simpleFoam→runフォルダ
メッシュができたので計算をさせます。
フォルダ移動、チュートリアルのコピーをします。
1 2 3 |
cd ../../run cp -r $FOAM_TUTORIALS/tutorials/incompressible/simpleFoam/pitzDaily simpleFoamRun cd simpleFoamRun |
そして先ほどcfMeshで作成したメッシュ情報を一式コピーします。
1 |
cp -r ../../mesh/cfMesh/constant/polyMesh constant/ |
そして以下適当に解析条件を設定していきます。
constant/turbulenceProperties
乱流モデルは使いません。
1 |
simulationType laminar; |
constant/transportProperties
できるだけ層流での計算をさせたいのであえて粘性を少し大きくします。
デフォルトの10倍(レイノルズ数$Re=\frac{UL}{\nu}=\frac{0.1\times 0.04}{10^{-4}}=40$)。
1 2 3 |
transportModel Newtonian; nu 1e-04; |
あとは境界条件を設定します。
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 |
dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0.1); boundaryField { ".*" { type noSlip; } inlet { type fixedValue; value uniform (0 0 0.1); } "outlet.*" { // type zeroGradient; type pressureInletOutletVelocity; value uniform (0 0 0); } } |
0/p
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
boundaryField { ".*" { type zeroGradient; } inlet { type zeroGradient; } "outlet.*" { // type fixedValue; // value uniform 0; type totalPressure; p0 uniform 0; value uniform 0; } } |
非圧縮の層流だとこれだけですね。
スキームも何でもよかったのでほとんどデフォルトのままです。
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 40 41 42 |
ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) bounded Gauss limitedLinear 1; div(phi,k) bounded Gauss limitedLinear 1; div(phi,epsilon) bounded Gauss limitedLinear 1; div(phi,omega) bounded Gauss limitedLinear 1; div(phi,v2) bounded Gauss limitedLinear 1; div((nuEff*dev2(T(grad(U))))) Gauss linear; div(nonlinearStress) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } wallDist { method meshWave; } |
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 |
solvers { p { solver GAMG; tolerance 1e-06; relTol 0.1; smoother GaussSeidel; } "(U|k|epsilon|omega|f|v2)" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-05; relTol 0.1; } } SIMPLE { nNonOrthogonalCorrectors 0; consistent yes; residualControl { p 1e-2; U 1e-3; "(k|epsilon|omega|f|v2)" 1e-3; } } relaxationFactors { equations { U 0.9; // 0.9 is more stable but 0.95 more convergent ".*" 0.9; // 0.9 is more stable but 0.95 more convergent } } |
これでだいたい計算条件設定は終了です。
あとは「system/controlDict」に計算制御や出力設定などをします。
まずは境界の流量出力から見ていきます。
境界面の流量出力
system/controlDict
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 |
application simpleFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 2000; deltaT 1; writeControl timeStep; writeInterval 100; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; functions { inlet { type surfaceFieldValue; libs ( "libfieldFunctionObjects.so" ); writeControl timeStep; log yes; writeFields no; regionType patch; name inlet; operation sum; fields ( phi ); // Optional (inherited) entries writePrecision 8; writeToFile yes; useUserTime yes; region region0; enabled yes; timeStart 0; timeEnd 10000; executeControl timeStep; executeInterval 1; writeInterval 100; } outlet1 { $inlet; name outlet1; } outlet2 { $inlet; name outlet2; } residuals { type solverInfo; //libs("libutilityFunctionObjects.so"); libs (utilityFunctionObjects); fields (".*" ); } continuityError1 { type continuityError; libs ( fieldFunctionObjects ); // Mandatory entries (unmodifiable)// Optional entries (runtime modifiable) phi phi; // Optional (inherited) entries writePrecision 8; writeToFile yes; useUserTime yes; region region0; enabled yes; log yes; timeStart 0; timeEnd 10000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 50; } } |
では計算を実行します。
1 |
simpleFoam > log.simpleFoam & |
ParaViewで確認します。
流量、残差、質量保存のデータは「postProcessing」に出力されます。
これをPythonなどでグラフ化すると以下のようになります。
それだけではなく、連続の式(非圧縮であれば$\nabla \cdot \boldsymbol{v}=0$)を満たすのでlocalやglocalが0に近くなっていれば問題ないです。
流量はoutlet1とoutlet2が自然流出の条件にしているので定常解析とは言え揺らぎながら落ちついていきます。
inetl部分は$Q=vA=0.1\times (\pi0.2^2 )=0.01256$[m3/s]です。
(ちょっと理論値とずれてる?)
境界面流量はこのようにfunction objectに面指定で流量出力を設定すると得ることができます。
まとめ
本記事では境界面の流量を出力する方法を解説しました。
境界面の流量はOpenFOAMの既存のfunction objectを使えば簡単に出力することができます。
では、c1~c5のような中途半端なところでの流量を知りたい場合はどうすればいいのか?
参考書
有限体積法の勉強にはこの本で間違いないです。
対流項目が中心差分で安定しないことを理論的に評価しておりとても参考になります。
OpenFOAMの基本的な設定を学ぶことができる書籍として以下の参考書がお勧めです。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ(NextPublishing))
本書はESI版ではなくFoundation版でありOpenFOAMv2.xで書かれたやや古い内容です。その点に注意してESI版のお使いのversionとの違いを丁寧に調べることができる方であれば、十分参考になる書籍です。
OpenFOAM全般の設定などまとめられた書籍として以下のものがあります。
こちらもFoundation版なので、ESI版のお使いのversionとの違いを丁寧に調べることができる方であれば、十分参考になる書籍です。
こちらのオープンCAE学会から技術書典にOpenFOAMでメッシュ作成【PDF版】が出ています。難しい内容ではありますが、OpenFOAMのメッシュに関する内容が1000円で学べます。