こんにちは（@t_kun_kamakiri）
本記事では、OpenFOAMを初めての方を対象に、チュートリアルを使ってメッシュ作成、計算の実行まで分かりやすく解説します。
👇自宅で流体解析をしてみませんか？
前回の記事では、OpenFOAMのチュートリアルから熱流体解析の解析設定を行いました。
上のPC画面に映っているような2枚の温度差のある平板内での流れをシミュレーションしています。
本記事では前回の記事で作成したモデルを用いて結果処理を行います。
チュートリアルには実験データがありますので、OpenFOAMで得られた結果と比較することで、シミュレーションの妥当性の評価と、実験と合わなかった際の着眼点をこの記事で示します。
- OpenFOAMを始めたいけど、どこから手をつければいいのか分からない…
- チュートリアルの流れを理解して、自分の解析に応用したい！
そんな方に向けて、基本操作を丁寧に説明するので、ぜひ一緒に手を動かしてみてください！
本記事で扱うチュートリアルはOpenFOAMがベンチマークとして下記の実験データとの比較のために行っており、より実践的な内容がまとめられています。
OpenFOAM v2412(WSL Ubuntu 22.04)
値の出力設定（sampleの設定）
すでに前回の記事で計算が終わっていることを前提に話をしていきます。
ここでは、$y$方向違いの温度分布や速度分布の結果を出力します。
sampleを使用します。
計算が終わったあとでも値を出力することができるため便利なツールです。
system/sample
|
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 |
type sets; libs (sampling); interpolationScheme cellPointFace; setFormat raw; fields ( T U ); sets { y0.1 { type face; axis x; start (-1 0.218 0); end (1 0.218 0); } y0.2 { type face; axis x; start (-1 0.436 0); end (1 0.436 0); } y0.3 { type face; axis x; start (-1 0.654 0); end (1 0.654 0); } y0.4 { type face; axis x; start (-1 0.872 0); end (1 0.872 0); } y0.5 { type face; axis x; start (-1 1.09 0); end (1 1.09 0); } y0.6 { type face; axis x; start (-1 1.308 0); end (1 1.308 0); } y0.7 { type face; axis x; start (-1 1.526 0); end (1 1.526 0); } y0.8 { type face; axis x; start (-1 1.744 0); end (1 1.744 0); } y0.9 { type face; axis x; start (-1 1.962 0); end (1 1.962 0); } } |
例えば、$L=0.076\,\text{m}$、$D=0.026\,\text{m}$、$H=2.18\,\text{m}$ですので、$y/H=0.5$の高さで、値を取得したい場合を考えます。
その場合は、以下のようにstart位置とend位置を指定した線上の値を取得することになります。
|
1 2 3 4 5 6 7 |
y0.5 { type face; axis x; start (-1 1.09 0); end (1 1.09 0); } |
$y=0.5*H=0.5*2.18=1.09$ということですね。
値の出力
値の出力方法は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 |
以下の部分ですね。
|
1 |
postProcess -func sample -latestTime |
-latestTimeというオプションを付けることで最後の結果だけを使ってデータを出力してくれます。
グラフ化（画像の出力）
グラフ化には以下のコマンドを使用します。
|
1 |
validation/plot |
こちらはスクリプトになっておりgnuplotを実行してグラフを作成しています。
|
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 |
plot_x_vs_T_U() { index="$1" sampleFile="$2" benchmarkDir="$3" gnuplot<<PLT_T_U set terminal pngcairo font "helvetica,20" size 1000, 1000 set grid set key left top set xrange [0:0.08] set xlabel "Channel width, x [m]" # Benchmark - Experimental benchmarkT="$benchmarkDir/mt_z0_${index}0_lo.dat" benchmarkU="$benchmarkDir/mv_z0_${index}0_lo.dat" # OpenFOAM samples="$sampleFile" # plot temperature profiles set yrange [285:310] set ylabel "Temperature [K]" set output "OF_vs_EXPT_T$index.png" plot \ benchmarkT u (\$1/1000):(\$2+273.15) t "Expt 0.$index" w p lt 1 pt 6, \ samples u 1:2 t "OpenFOAM 0.$index" w l lt -1 # plot velocity profiles set yrange [-0.2:0.2] set ylabel "Vertical velocity component, Uy [m/s]" set output "OF_vs_EXPT_U$index.png" plot \ benchmarkU u (\$1/1000):(\$2) t "Expt 0.$index" w p lt 1 pt 6, \ samples u 1:4 title "OpenFOAM 0.$index" w l lt -1 PLT_T_U } #------------------------------------------------------------------------------ # Requires "gnuplot" command -v gnuplot >/dev/null || { echo "FOAM FATAL ERROR: gnuplot not found - skipping graph creation" 1>&2 exit 1 } SETSDIR="../postProcessing/sample" [ -d "$SETSDIR" ] || { echo "FOAM FATAL ERROR: result sets not available in directory $SETSDIR" 1>&2 exit 1 } #------------------------------------------------------------------------------ echo "" echo "# Plot:" echo "" # paths to data LATESTTIME=$(ls $SETSDIR) OFDATAROOT="$SETSDIR/$LATESTTIME" EXPTDATAROOT=./exptData # generate temperature and velocity profiles set="1 3 4 5 6 7 9" for i in $set do echo " processing temperature and velocity profiles at y/yMax of 0.$i" OF="$OFDATAROOT/y0.${i}_T_U.xy" plot_x_vs_T_U "$i" "$OF" "$EXPTDATAROOT" done echo "End" |
gnuplotに関する基本的な内容はこちらに記載しています。
代表的な$y/H=0.5$の温度と流速の実験とOpenFOAMの比較を見てみましょう。
パット見全然あっていないですね….
なぜでしょうか…
設定の見直し
こちらの文章を読むと以下のよう書いております。
実験での設定温度は明確には書いておりませんが、実験とOpenFOAMの壁面の温度からさほど差は無いように感じます。
上記の19.6℃はhotとcoldとの温度差のことです。
解析の条件設定が間違っているということもありますが、まずは最低限のことを確認してから条件の設定を見直すようにします。
最低限の確認とは以下の事です。
- 残差の確認
- 連続式の誤差の確認
- 流れ場の可視化
まずこれは確認しましょう。
ということで、初期化するために以下のコマンドを実行します。
|
1 |
./Allclean |
残差、連続式の誤差を出力するために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 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 170 171 172 173 174 |
application buoyantSimpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 1000; deltaT 1; writeControl timeStep; writeInterval 100; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; functions { 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 8000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 50; } residuals { type solverInfo; libs ( "libutilityFunctionObjects.so" ); fields ( ".*" ); } fieldMinMax1 { type fieldMinMax; libs ( fieldFunctionObjects ); // Mandatory entries (unmodifiable)// Mandatory entries (runtime modifiable) mode magnitude; fields ( U p T ); // Optional entries (runtime modifiable) location yes; writePrecision 8; writeToFile yes; useUserTime yes; region region0; enabled yes; log yes; timeStart 0; timeEnd 8000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 100; } yPlus1 { type yPlus; libs ( fieldFunctionObjects ); // Mandatory entries (unmodifiable)// Optional (inherited) entries writePrecision 8; writeToFile yes; useUserTime yes; region region0; enabled yes; log yes; timeStart 0; timeEnd 8000; executeControl timeStep; executeInterval 1; writeControl timeStep; writeInterval 100; } //https://doc.openfoam.com/2306/tools/post-processing/function-objects/field/wallHeatFlux/ wallHeatFlux1 { // Mandatory entries (unmodifiable) type wallHeatFlux; libs (fieldFunctionObjects); // Optional entries (runtime modifiable) //patches (<patch1> ... <patchN>); // (wall1 "(wall2|wall3)"); //qr qr; // Optional (inherited) entries // writePrecision 8; // writeToFile true; // useUserTime true; // region region0; // enabled true; // log true; // timeStart 0; // timeEnd 8000; executeControl timeStep; executeInterval 1; writeControl writeTime;//timeStep; writeInterval 1; } //https://doc.openfoam.com/2306/tools/post-processing/function-objects/field/heatTransferCoeff/ heatTransferCoeff1 { // Mandatory entries (unmodifiable) type heatTransferCoeff; libs (fieldFunctionObjects);//(fieldFunctionObjects); field T;//<field>; // patches ("([a-zA-Z]+-(left|right)-[0-9]+)");//(<patch1> <patch2> ... <patchN>); patches ( hot cold frontAndBack topAndBottom );//(<patch1> <patch2> ... <patchN>); htcModel fixedReferenceTemperature;//<htcModel>; TRef 297.95; //24.8deg // UInf (20 0 0); // Cp CpInf; // CpInf 4000; // rho rhoInf; // rhoInf 1.2; // Optional (inherited) entries // result <fieldResult>; region region0; // enabled true; // log true; // timeStart 0; // timeEnd 8000; executeControl timeStep; executeInterval 1; writeControl writeTime;;//timeStep; writeInterval 1; } } |
途中の結果も確認したいのでpurgeWrite 0;にしておきます。0にしておくことで、writeInterval 50;で設定されたステップごとのデータは全て出力されます。
ついでに色々と出力できるようにしてみました。
-
type continuityError;：連続式の誤差の出力 -
type solverInfo;：残差の出力 -
type fieldMinMax;：物理量の最大最小の出力 -
type yPlus;：$y^+$の出力 -
type wallHeatFlux;：熱流束の出力 -
type heatTransferCoeff;：熱伝達率の出力
設定の変更が終われば再計算を行います。
|
1 |
./Allrun |
Allrunにはメッシュ作成、計算実行、グラフ化まで一連の操作がまとめられているので、スクリプトを実行するだけで済みます。
残差の確認
残差は以下のファイルにデータが出力されているので、グラフにして確認します。
postProcessing/residuals/0/solverInfo.dat
残差はステップ数1000で収束しているような微妙な感じですね。
もしかするとステップ数が足りていない可能性があります。
連続式の誤差
連続式の誤差が大きい値を取っていると計算値が発散している場合があるため確認しておきます。
continuityError1/0/continuityError.dat
こちらも問題なさそうです。
流れの可視化
残差から1000ステップで収束しているように見えましたが、流れの方が安定になっているのかを確認しておきます。
ParaViewで確認します。
右のグラフは$y/H=0.5$での高さの温度分布になります。
なるほど。
1000ステップでは流れが安定になっていないため、実験データとOpenFOAMの結果があっていなかったのですね。
ステップ数の変更
最大ステップ数をendTime 5000;にします。
残差
連続式の誤差
可視化
では、実験とOpenFOAMの結果を見てみましょう。
代表的な高さ$y/H=0.5$での比較が以下です。
だいぶ近くなりましたね。
ただ、アニメーションを見ると5000ステップでもまだ収束していないように見えます。
その場合は、例えば緩和係数を少し大きな値にすると、収束までが速くなります。
system/fvSolution
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
relaxationFactors { fields { rho 1.0; p_rgh 0.7; } equations { U 0.7;//0.3; h 0.7;//0.3; "(k|epsilon|omega)" 0.7; } |
設定を変更して再度Allrunスクリプトにより計算を実行します。
緩和係数を大きくしすぎると、計算値が発散しやすいので注意が必要ですが、安定に計算が進んでいれば問題ないでしょう。
結果はこちらです。
波形が落ち着くのがだいぶ早くなりましたね。
壁面近傍の流速が特に実験と合っていないように見えます。
さらに改善しようとするならば、メッシュサイズ（$y^+$の確認）、乱流モデル、離散化スキームなどを見直してみると良いでしょう。
ただ、1つの実験データだけを鵜呑みにして神経質になって合わせに行くのもナンセンスです。なぜなら、実験データにも多少バラツキがあるからです。
同じ実験を行ったとしても、違う結果になることもあるので、バラツキの範囲を考えてCAE解析の結果が妥当かどうかを考える必要があります。
まとめ
本記事では、OpenFOAM を用いたキャビティ内の乱流自然対流解析の結果処理について詳しく解説しました。前回作成したモデルを用い、sample を活用して温度・速度分布のデータを出力し、実験データと比較することで、シミュレーションの妥当性を評価しました。結果が実験と一致しない場合、設定の見直しが必要となるため、残差・連続式の誤差・流れ場の可視化を確認する重要性を説明しました。
解析結果が実験と異なる原因として、計算ステップ数の不足や緩和係数の設定が影響することを確認し、計算の収束性を向上させるための調整方法を紹介しました。さらに、壁面近傍の流速の不一致を改善するための次のステップとして、メッシュサイズの最適化（$y^+$ の確認）や乱流モデルの変更が有効であることに言及しました。
計算力学技術者のための問題アプリ
計算力学技術者熱流体２級、1級対策アプリをリリースしました。
- 下記をクリックしてホームページでダウンロードできます。
- LINE公式に登録すると無料で問題の一部を閲覧できます。
※LINEの仕様で数式がずれていますが、アプリでは問題ありません。
- 計算力学技術者の熱流体2級問題アプリ作成
リリース後も試行錯誤をしながら改善に努め日々アップデートしています。
備忘録として作成の過程を綴っています。
OpenFOAMに関する技術書を販売中！
OpenFOAMを自宅で学べるシリーズを販売中です。
OpenFOAM初学者から中級者向けの技術書となっていますので、ぜひよろしくお願いいたします。
次回の技術書典18に向けて内容を考え中です。
乞うご期待！！
お勧めの参考書
本記事の内容について触れている書籍がこちらです。
CFD（流体解析）のガイドブックというタイトルだけあって、熱流体に関する内容は網羅的に書かれています。
乱流モデルの数式の展開が非常に丁寧なのはこちらの参考書です。
今まで読んだ本の中で途中式もしっかり書いてあって一番丁寧でした。
乱流モデルの話だけでなく、混相流（気液、固液）や粒子法、浅水方程式の話も乗っているので重宝しています。
乱流モデルはこちらもお勧めです。
前半は数値シミュレーションの離散化の話で、後半に乱流モデルの話が出てきます。
乱流モデルのざっくりした解説と流体全般の基礎知識にはこちらがちょうど良いでしょう。
Follow @t_kun_kamakiri
