こんにちは（@t_kun_kamakiri）
今回の記事では円管内の流れの流体解析の設定方法を示します。
2025年5月末の技術書典18に向けて、OpenFOAMを用いた熱流体解析のひとつのテーマとして円管内の流れを考えています。
- 円管壁面に熱流束$\dot{q}$一定値を与えて熱流体解析を行う
OpenFOAMでシミュレーションする設定方法を示す - 今回は乱流モデルを使用しないモデル（層流モデル）で解析結果を示す
OpenFOAM v2412（WSL2 Ubuntu 22.04）
【宣伝】円管内流れ（非圧縮流れ）
本代の円管内の流れの解析設定に入る前に宣伝になります。
本記事のメッシュ作成や解析設定について、基礎から学びたい方向けに、「円管内流れの層流解析」として、モデル作成からグラフ作成までの一連の流れをまとめたものを技術書としてまとめています。
ページ数は136ページです。
技術書では以下の内容を取り扱っています。
- 円管内流れ（層流）の速度分布の理論
- OpenFOAMで円管内流れの解析（blockMesh）
- OpenFOAMで円管内流れの解析（snappyHexMesh）
↓内容のチラ見せ
本書の解析モデルはgithubに置いているので、本書と照らし合わせながら読み進めることができます。
興味ある方はぜひ手に取ってみてください。
円管内の流れ解析設定
使用するソルバはbuoyantSimpleFoam（浮力を考慮した定常の熱流体ソルバ）です。
メッシュはblockMeshで作成しています。
さらに、流速の発達区間を設けるために助走区間を設けています。
層流における流速と温度の助走区間は、レイノルズ数$Re$、プラントル数$Pr$、および直径$d$より以下で推算できることが知られています。
- 速度助走区間：$L=0.05\,d\,Re$
- 温度助走区間：$L=0.05\,d\,Pr\,Re$
助走区間を過ぎた$z\geq 0.6$では熱流束$\dot{q}=100\text{W}/\text{m}^2$一定値を与えています。
物性値
空気を想定していますが、値はキリの良い数値に変更しています。
- モル質量：28.0[g/mol]
- 定圧比熱：1000.0 [J/kg K]
- 粘性係数：1.8E-5 [Pa・s]
- プラントル数：0.7
|
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.0; } thermodynamics { Cp 1000.0; Hf 0; } transport { mu 1.8e-5; Pr 0.70; } } |
流入の流速が$U=0.15\,\text{m}/\text{s}$、円柱の直径が$d=20\,\text{mm}$であることから、レイノルズ数$Re$が求まります。
Re=\frac{Ud }{\mu/\rho}
\end{align*}
ただし、mixture pureMixture;としていることから、 $\rho$は理想気体の状態方程式$p=\rho R_{MW}T$から決まります。
$ R_{MW}$は気体定数と呼ばれ、気体ごとに決まる定数です。
具体的には一般気体定数$R$[J/mol K]とモル質量$MW$[kg/mol]を用いて、以下のように定義されています。
R_{MW}=\frac{R}{MW}=\frac{8.3145}{28.0E-3}=296.95\, \text{J}/\text{kg K}
\end{align*}
ということなので、密度は次式で求まります。
\rho=\frac{p}{R_{MW}T}=\frac{1E5}{296.95\times 293.15}=1.148\,\text{kg}/\text{m}^3
\end{align*}
この値を用いるとレイノルズ数$Re$が求まります。
Re=\frac{Ud }{\mu/\rho}=\frac{0.3\times 20R-3}{1.8E-5/1.148}=383
\end{align*}
円管内流れではレイノルズ数$Re=2300$あたりから乱流に遷移するため、$Re=383$では層流であると判断できます。
そのため乱流モデルは使用せずに、定常解析の設定で解析を行います。
ゆえに、以下の設定で行うことにします。
- 定常解析
- 層流モデル（乱流モデル無し）
- 重力の設定は無し
また、助走区間については以下の推算ができます。
- 速度助走区間：$L=0.05\,d\,Re= 0.382$
- 温度助走区間：$L=0.05\,d\,Pr\,Re=0.267$
境界条件
乱流モデルを使用しないとなると、境界条件で使用するファイルは0/U, 0/T, 0/p_rgh, 0/pになります。
p_rghは静圧から流体質量による圧力分を差し引いた量$p_{rgh}=p-rgh$なので、0/pは0/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 |
dimensions [ 0 1 -1 0 0 0 0 ]; internalField uniform (0 0 0.3); boundaryField { #includeEtc "caseDicts/setConstraintTypes"; runwalls { type noSlip; } walls { type noSlip; } inlet { type fixedValue; value uniform (0 0 0.3); } outlet { // type pressureInletOutletVelocity; // value uniform (0 0 0.1); type zeroGradient; } } |
|
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 |
dimensions [ 1 -1 -2 0 0 0 0 ]; internalField uniform 101325; boundaryField { #includeEtc "caseDicts/setConstraintTypes"; runwalls { type fixedFluxPressure; value uniform 101325; } walls { type fixedFluxPressure; value uniform 101325; } inlet { type fixedFluxPressure; value uniform 101325; // type zeroGradient; } outlet { type fixedValue; value uniform 101325; // type totalPressure; // U U; // phi phi; // psi none; // gamma 1; // rho rho; // p0 uniform 101325; } } |
|
1 2 3 4 5 6 7 8 9 10 11 12 |
dimensions [ 1 -1 -2 0 0 0 0 ]; internalField uniform 101325; boundaryField { #includeEtc "caseDicts/setConstraintTypes"; ".*" { type calculated; value uniform 101325; } } |
|
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 |
dimensions [ 0 0 0 1 0 0 0 ]; internalField uniform 293; boundaryField { #includeEtc "caseDicts/setConstraintTypes"; runwalls { type fixedValue; value uniform 293.5; } walls { type externalWallHeatFluxTemperature; mode flux; q uniform 100; kappaMethod fluidThermo; value uniform 293; // type fixedValue; // value uniform 350; } inlet { type fixedValue; value uniform 293; } outlet { type zeroGradient; } } |
結果
こちらが温度分布がこちらです。
温度助走区間を過ぎても温度は変化し続けているのがわかります。
たまにある誤解として、「温度助走区間を過ぎた下流側では、熱的に発達した状態であるため温度変化がない」と考える人がいますが（自分がそうだった）、温度助走区間を過ぎても流体側の温度と壁面温度$T_w$に温度差がある限りは、壁面温度に漸近するまでは温度変化はし続けます。
また、本記事では扱いませんでしたが、円管内層流流れのヌセルト数は一定値であることが知られています。
- 壁面熱流束が一定の場合：$Nu=4.36$
- 壁面温度一定の場合：$Nu=3.66$
ここで、「あーなるほど、ヌセルト数がわかると$Nu=\frac{hd}{\lambda}$などから、熱伝達率$h=\frac{Nu\,\lambda}{d}$がわかるのか」….となるかというと、そう一筋縄ではいかないのです。
✅CFDで移動現象論111例題- Ansys Fluentによる計算解法 –
これどうやって混合平均温度計算しているのかと思って確認したら強引に求めてた。https://t.co/3OrHCOf3KK pic.twitter.com/GB7UJs0kq0— カマキリ🐲技術書典18 (@t_kun_kamakiri) May 5, 2025
ここにも記載しているよ円管内の流れにおいて熱流束$\dot{q}=-h(T_w-T_m(z))$という形で、混合平均温度$T_m(z)$を用いて、ヌセルト数一定値を導出しています。
混合平均温度
T_m &= \frac{\int_A \rho C_p u T \, dA }{ \int_A \rho C_p u \, dA }\\
&= \frac{ \int_0^{R} \rho C_p u(r) T \, 2\pi r\,dr}{ \int_0^R \rho C_p u(r) \, 2\pi r\,dr}
\end{align*}
- $\rho$：密度
- $u$：流速
- $T$：温度
- $A$：管断面積
- $C_p$：定圧比熱
円管内流れ（層流）の速度分布
u(r) = u_{\mathrm{max}} \left( 1 – \left( \frac{r}{R} \right)^2 \right)
\end{align*}
混合平均温度は直接OpenFOAMで出力できないため（自分が知っている限りでは）、出力方向に工夫が必要です。
$T_m(z)$のように$z$に依存するため、プログラムを実装して出力しようかなと考えています。
まとめ
今回は円管内流れの設定を行い、温度分布を可視化しました。
混合平均温度$T_m(z)$の出力をどうしようかな～と考えています。
次回の記事では、混合平均温度を出力するためにプログラムの実装を行います。
技術書典18では、本記事の「2次元円柱まわりの熱流体解析」やそのほかの熱流体のテーマを取り扱う予定です。
次回の技術書典18に向けて内容を考え中です。
乞うご期待！！
計算力学技術者のための問題アプリ
計算力学技術者熱流体２級、1級対策アプリをリリースしました。
- 下記をクリックしてホームページでダウンロードできます。
- LINE公式に登録すると無料で問題の一部を閲覧できます。
※LINEの仕様で数式がずれていますが、アプリでは問題ありません。
- 計算力学技術者の熱流体2級問題アプリ作成
リリース後も試行錯誤をしながら改善に努め日々アップデートしています。
備忘録として作成の過程を綴っています。
OpenFOAMに関する技術書を販売中！
OpenFOAMを自宅で学べるシリーズを販売中です。
OpenFOAM初学者から中級者向けの技術書となっていますので、ぜひよろしくお願いいたします。
お勧めの参考書
本記事の内容について触れている書籍がこちらです。
CFD（流体解析）のガイドブックというタイトルだけあって、熱流体に関する内容は網羅的に書かれています。
乱流モデルの数式の展開が非常に丁寧なのはこちらの参考書です。
今まで読んだ本の中で途中式もしっかり書いてあって一番丁寧でした。
乱流モデルの話だけでなく、混相流（気液、固液）や粒子法、浅水方程式の話も乗っているので重宝しています。
乱流モデルはこちらもお勧めです。
前半は数値シミュレーションの離散化の話で、後半に乱流モデルの話が出てきます。
乱流モデルのざっくりした解説と流体全般の基礎知識にはこちらがちょうど良いでしょう。
Follow @t_kun_kamakiri
