こんにちは（@t_kun_kamakiri）
本記事では、OpenFOAMを初めての方を対象に、チュートリアルを使ってメッシュ作成、計算の実行まで分かりやすく解説します。
👇自宅で流体解析をしてみませんか？
- OpenFOAMを始めたいけど、どこから手をつければいいのか分からない…
- チュートリアルの流れを理解して、自分の解析に応用したい！
そんな方に向けて、基本操作を丁寧に説明するので、ぜひ一緒に手を動かしてみてください！
本記事で扱うチュートリアルはOpenFOAMがベンチマークとして下記の実験データとの比較のために行っており、より実践的な内容がまとめられています。
OpenFOAM v2412(WSL Ubuntu 22.04)
チュートリアルのコピー
OpenFOAMは一から解析設定をすることは少なく、既存のチュートリアルから自分が解析したいものに近いものを探し出し、編集して使うことを前提としています。
そのため数多くのチュートリアルが存在します。
今回はこの中で熱の問題も解くことができるソルバbuoyantSimpleFoamを使います。
その中には、buoyantCavity というチュートリアルがありますので、そちらを作業フォルダにコピーします。
※作業フォルダンはどこでも良いです。
|
1 |
cp -r $FOAM_TUTORIALS/heatTransfer/buoyantSimpleFoam/buoyantCavity . |
$FOAM_TUTORIALS = /usr/lib/openfoam/openfoam2412/tutorials
では、フォルダを移動します。
|
1 |
cd buoyantCavity |
フォルダの構成は以下のようになっています。
|
1 2 3 4 5 6 7 8 9 |
$ tree -L 1 . ├── 0.orig ├── Allclean ├── Allrun ├── README.md ├── constant ├── system └── validation |
validationフォルダに実験データが保存されていますが、本記事では扱いません。
メッシュ作成
まずはメッシュ作成を行います。
メッシュ作成はOpenFOAMにblockMeshを用います。
blockMeshの設定はsystem/blockMeshDictで起こないます。
|
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 |
vertices ( ( 0 0 -260) // 0 (76 0 -260) // 1 (76 2180 -260) // 2 ( 0 2180 -260) // 3 ( 0 0 260) // 4 (76 0 260) // 5 (76 2180 260) // 6 ( 0 2180 260) // 7 ); edges ( ); blocks ( hex (0 1 2 3 4 5 6 7) (35 150 15) simpleGrading (1 1 1) ); boundary ( frontAndBack { type wall; faces ( (0 1 5 4) (2 3 7 6) ); } topAndBottom { type wall; faces ( (4 5 6 7) (3 2 1 0) ); } hot { type wall; faces ( (6 5 1 2) ); } cold { type wall; faces ( (4 7 3 0) ); } ); mergePatchPairs ( ); |
テキストでの設定ですので、対応する節点の座標と面（名前を付ける）を設定することで、六面体メッシュを作成しています。
$L=0.076\,\text{m}$、$D=0.026\,\text{m}$、$H=2.18\,\text{m}$
では、blockMeshを実行してメッシュを行います。
|
1 |
blockMesh |
メッシュが終わりましたら、一度ParaViewでメッシュ状態の確認をします。
ParaViewでOpenFOAMの解析結果を可視化するには、拡張子が.foamの空ファイルを作成し、それをParaViewで読み込むだけで簡単に結果を表示できます。
|
1 |
touch post.foam |
ParaViewで確認します。
メッシュ状態が良ければ、次に解析設定に移ります。
解析の設定
ここからは解析の設定を行います。
今回はチュートリアルの設定のまま変更をしませんので、ファイルの中を確認するだけになります。
物性値の設定
まずは物性値の設定。
constant/thermophysicalProperties
|
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 |
thermoType { type heRhoThermo; mixture pureMixture; transport const; thermo hConst; equationOfState perfectGas; specie specie; energy sensibleEnthalpy; } mixture { specie { molWeight 28.96; } thermodynamics { Cp 1004.4; Hf 0; } transport { mu 1.831e-05; Pr 0.705; } } |
ここでは、以下の値の設定をしています。
- モル質量：$28.96$ [g/mol]
- 定圧比熱：$1004.4$ [J/kg K]
- 粘性係数：$1.831e-05$ [Pa・s]
- プラントル数：$0.705$ [-]
乱流モデルの設定
ここでは乱流モデルの設定をしています。
constant/turbulenceProperties
|
1 2 3 4 5 6 7 8 9 10 |
simulationType RAS; RAS { RASModel kOmegaSST; turbulence on; printCoeffs on; } |
乱流モデルは$k$-$\omega$SSTにしています。
境界条件の設定
0.orig/U、0.orig/p_rgh、0.orig/Tから設定します。
こちらは乱流モデルの有無に限らず必ず設定するファイルです。
p_rgh は、圧力 p から重力ポテンシャルを引いたもの
0.orig/pは0.orig/p_rghから計算される設定にしています。
0.orig/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 |
dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { frontAndBack { type noSlip; } topAndBottom { type noSlip; } hot { type noSlip; } cold { type noSlip; } } |
流速は全て滑り無し条件です。
0.orig/p_rgh
|
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 |
dimensions [1 -1 -2 0 0 0 0]; internalField uniform 1e5; boundaryField { frontAndBack { type fixedFluxPressure; value uniform 1e5; } topAndBottom { type fixedFluxPressure; value uniform 1e5; } hot { type fixedFluxPressure; value uniform 1e5; } cold { type fixedFluxPressure; value uniform 1e5; } } |
fixedFluxPressure は、壁面（壁境界）における圧力勾配を速度と粘性応力のバランスから決定します。
0.orig/p
|
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 |
dimensions [1 -1 -2 0 0 0 0]; internalField uniform 1e5; boundaryField { frontAndBack { type calculated; value $internalField; } topAndBottom { type calculated; value $internalField; } hot { type calculated; value $internalField; } cold { type calculated; value $internalField; } } |
圧力pはp_rghから計算されています。
0.orig/T
|
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 |
dimensions [0 0 0 1 0 0 0]; internalField uniform 293; boundaryField { frontAndBack { type zeroGradient; } topAndBottom { type zeroGradient; } hot { type fixedValue; value uniform 307.75; // 34.6 degC } cold { type fixedValue; value uniform 288.15; // 15 degC } } |
type zeroGradient;は温度勾配0ですので、断熱条件$\dot{q}=-\lambda\frac{\partial T}{\partial n}=0$を意味しています。
続いて、乱流モデルに関係する境界条件の設定です。
乱流モデルに壁関数の記事もご参照ください。
0.orig/k
|
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 |
dimensions [0 2 -2 0 0 0 0]; internalField uniform 3.75e-04; boundaryField { frontAndBack { type kqRWallFunction; value uniform 3.75e-04; } topAndBottom { type kqRWallFunction; value uniform 3.75e-04; } hot { type kqRWallFunction; value uniform 3.75e-04; } cold { type kqRWallFunction; value uniform 3.75e-04; } } |
0.orig/omega
|
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 |
dimensions [0 0 -1 0 0 0 0]; internalField uniform 0.12; boundaryField { frontAndBack { type omegaWallFunction; value uniform 0.12; } topAndBottom { type omegaWallFunction; value uniform 0.12; } hot { type omegaWallFunction; value uniform 0.12; } cold { type omegaWallFunction; value uniform 0.12; } } |
0.orig/nut
|
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 |
dimensions [0 2 -1 0 0 0 0]; internalField uniform 0; boundaryField { frontAndBack { type nutUWallFunction; value uniform 0; } topAndBottom { type nutUWallFunction; value uniform 0; } hot { type nutUWallFunction; value uniform 0; } cold { type nutUWallFunction; value uniform 0; } } |
0.orig/alphat
|
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 |
dimensions [1 -1 -1 0 0 0 0]; internalField uniform 0; boundaryField { frontAndBack { type compressible::alphatWallFunction; Prt 0.85; value uniform 0; } topAndBottom { type compressible::alphatWallFunction; Prt 0.85; value uniform 0; } hot { type compressible::alphatWallFunction; Prt 0.85; value uniform 0; } cold { type compressible::alphatWallFunction; Prt 0.85; 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 43 44 45 46 47 |
ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) bounded Gauss limitedLinear 0.2; energy bounded Gauss limitedLinear 0.2; div(phi,K) $energy; div(phi,h) $energy; turbulence bounded Gauss limitedLinear 0.2; div(phi,k) $turbulence; div(phi,epsilon) $turbulence; div(phi,omega) $turbulence; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear orthogonal; } interpolationSchemes { default linear; } snGradSchemes { default orthogonal; } 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 42 43 44 45 46 47 48 49 50 51 |
solvers { p_rgh { solver GAMG; tolerance 1e-7; relTol 0.01; smoother DICGaussSeidel; } "(U|h|k|epsilon|omega)" { solver PBiCGStab; preconditioner DILU; tolerance 1e-8; relTol 0.1; } } SIMPLE { momentumPredictor no; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; residualControl { p_rgh 1e-4; U 1e-4; h 1e-4; // possibly check turbulence fields "(k|epsilon|omega)" 1e-3; } } relaxationFactors { fields { rho 1.0; p_rgh 0.7; } equations { U 0.3; h 0.3; "(k|epsilon|omega)" 0.7; } } |
計算が発散したり収束が悪かったりする場合は、relaxationFactorsの各物理量の値を小さくするなどして調整します。
計算制御の設定
最後に計算時間（ステップ数）の指定を行います。
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 |
application buoyantSimpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 1000; deltaT 1; writeControl timeStep; writeInterval 50; purgeWrite 3; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; |
計算の実行
Allrunスクリプトには、このチュートリアルを実行するためのコマンドが書かれています。
Allrun
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
#!/bin/sh cd "${0%/*}" || exit # Run from this directory . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions #------------------------------------------------------------------------------ restore0Dir runApplication blockMesh runApplication $(getApplication) runApplication postProcess -func sample -latestTime if notTest "$@" then runApplication validation/plot fi |
Allrunスクリプトを実行するとメッシュ作成から計算実行までの一連の流れを行ってくれますが、はじめは勉強だと思ってひとつずつ実行してどのような挙動をするのかを確認する方が良いでしょう。
まずは、0.origのオリジナルフォルダ0フォルダとしてコピーします。
Allrunスクリプトではrestore0Dirとなっている箇所です。
|
1 |
cp -r 0.orig 0 |
０フォルダがコピーできたら、次にメッシュ作成ですが、メッシュ作成は既にblockMeshを実行してできているので飛ばします。
以下のコマンドで計算実行が行われます。
runApplication $(getApplication)の箇所です。
|
1 |
buoyantSimpleFoam |
計算が終われば結果をParaViewで確認してみましょう。
ParaViewでは$y=0$の高さの温度データをグラフにしています。
まとめ
本記事では、OpenFOAM を用いたキャビティ内の乱流自然対流解析について、チュートリアルを通じて基本的な流れを解説しました。
計算実行の流れとしては、チュートリアルをコピーし、必要な設定を行った後、buoyantSimpleFoam を実行し、結果を ParaView で可視化することで解析内容を確認しました。
今回の内容を通じて、OpenFOAM における基本的な解析の進め方や設定ファイルの構成について理解が深まったのではないでしょうか。
以下、余裕のある方向けにAllrunスクリプトに関する補足を付け加えています。
Allrunスクリプトの説明
Allrunスクリプトの内容を知らなくても良いですが、少し調べてみるのも勉強になります。
Allrunの中は3行目の. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctionsでさらにスクリプトを読むようにしており、このスクリプトがAllrunの重要な役割を担っています。
$WM_PROJECT_DIRが何なのかを確認してみましょう。
|
1 2 |
$ echo $WM_PROJECT_DIR /usr/lib/openfoam/openfoam2412 |
つまり、ここのパスに従ってRunFunctionsを読み込んでいるということになります。
vi /usr/lib/openfoam/openfoam2412/bin/tools/RunFunctionsなので中身を確認してみると良いでしょう。
もしくは、以下のコマンドで作業フォルダにコピーして確認しても良いです。
|
1 |
cp -r /usr/lib/openfoam/openfoam2412/bin/tools/RunFunctions . |
AllrunスクリプトはrunApplication $(getApplication)としているため、runApplication とgetApplicationを確認します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
getApplication() { # Re-use positional parameters for automatic whitespace elimination set -- $(foamDictionary -entry application -value system/controlDict 2>/dev/null) if [ "$#" -eq 1 ] then echo "$1" else echo "Error getting 'application' from system/controlDict" 1>&2 echo false # Fallback return 1 fi } |
foamDictionary -entry application -value 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 |
$ foamDictionary -value system/controlDict FoamFile { version 2; format ascii; class dictionary; object controlDict; } application buoyantSimpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 1000; deltaT 1; writeControl timeStep; writeInterval 50; purgeWrite 3; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; |
つまりsystem/controlDictの中を辞書型で出力しています。
例えば、
- キー：application、値： buoyantSimpleFoam;
- キー：startFrom 、値： startTime;
のように、キーと値が対になっています。
|
1 2 |
$ foamDictionary -value system/controlDict -entry application buoyantSimpleFoam |
つまり、runApplication $(getApplication)の部分は、runApplication buoyantSimpleFoamとなっているわけです。
そうするとrunApplicationに対してbuoyantSimpleFoamは第一引数になっているので、runApplicationの中の$1をbuoyantSimpleFoamに置き換えて読んでいくわけです。
|
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 |
runApplication() { local appName appRun optValue logFile logMode # Any additional parsed arguments (eg, decomposeParDict) local appArgs # Parse options until executable is encountered while [ "$#" -gt 0 ] && [ -z "$appRun" ] do case "$1" in ('') ;; # Ignore junk (-a | -append) logMode=append ;; (-o | -overwrite) logMode=overwrite ;; (-s | -suffix) logFile=".$2" shift ;; (-decompose-dict=*) optValue="${1#*=}" case "$optValue" in ('' | none | false) ;; ## Ignore (*) appArgs="$appArgs -decomposeParDict $optValue" ;; esac ;; (-decomposeParDict) optValue="$2" shift case "$optValue" in ('' | none | false) ;; ## Ignore (*) appArgs="$appArgs -decomposeParDict $optValue" ;; esac ;; (*) appRun="$1" ;; esac shift done appName="${appRun##*/}" logFile="log.$appName$logFile" if [ -f "$logFile" ] && [ -z "$logMode" ] then echo "$appName already run on $PWD:" \ "remove log file '$logFile' to re-run" else echo "Running $appRun on $PWD" if [ "$logMode" = append ] then $appRun $appArgs "$@" >> $logFile 2>&1 else $appRun $appArgs "$@" > $logFile 2>&1 fi fi } |
細かいところまで見ていかなくても、appRun="$1"がappRun="buoyantSimpleFoam"になっているので、$appRun $appArgs "$@" > $logFile 2>&1でbuoyantSimpleFoam > log.buoyantSimpleFoamとなっているのがわかるでしょう。
実際には既にlog.buoyantSimpleFoamが存在すると計算が実行されているものとして、計算が実行されないという仕組みになっています。
計算力学技術者のための問題アプリ
計算力学技術者熱流体２級、1級対策アプリをリリースしました。
- 下記をクリックしてホームページでダウンロードできます。
- LINE公式に登録すると無料で問題の一部を閲覧できます。
※LINEの仕様で数式がずれていますが、アプリでは問題ありません。
- 計算力学技術者の熱流体2級問題アプリ作成
リリース後も試行錯誤をしながら改善に努め日々アップデートしています。
備忘録として作成の過程を綴っています。
OpenFOAMに関する技術書を販売中！
OpenFOAMを自宅で学べるシリーズを販売中です。
OpenFOAM初学者から中級者向けの技術書となっていますので、ぜひよろしくお願いいたします。
次回の技術書典18に向けて内容を考え中です。
乞うご期待！！
お勧めの参考書
本記事の内容について触れている書籍がこちらです。
CFD（流体解析）のガイドブックというタイトルだけあって、熱流体に関する内容は網羅的に書かれています。
乱流モデルの数式の展開が非常に丁寧なのはこちらの参考書です。
今まで読んだ本の中で途中式もしっかり書いてあって一番丁寧でした。
乱流モデルの話だけでなく、混相流（気液、固液）や粒子法、浅水方程式の話も乗っているので重宝しています。
乱流モデルはこちらもお勧めです。
前半は数値シミュレーションの離散化の話で、後半に乱流モデルの話が出てきます。
乱流モデルのざっくりした解説と流体全般の基礎知識にはこちらがちょうど良いでしょう。
Follow @t_kun_kamakiri
