寝室が早く冷えるエアコン吹出角度を調べてみた
OpenFOAMのbuoyantPimpleFoamで遊んでみました。
このメモはなに？
趣味でやってるシミュレーションで遊んだ内容をまとめたチラシの裏です。
2024年4月末頃から入門書籍を購入したり公式wikiの"3 weeks"チュートリアル(※下記リンク参照)を頼りにopenfoamを覚え始めて約2ヶ月が経ち、ようやく自分で3Dモデルの作成から計算・可視化までできるようになりました。
備忘のため、およびアウトプットすることによる自身への定着のためにこのメモにまとめてみました。
目次
シミュレーションの概要
夏場の寝室をエアコンで早く冷やすには吹き出し口の上下方向の角度はどのくらいが良いのか？大まかな傾向が知りたいなーとふと興味が湧き調べてみたくなりました。
openfoamの練習題材にちょうど良いなと思いシミュレーションしてみました。
いつも寝室のエアコンを「強」でONにして部屋を出て少し待ってから入るのでそれを想定した条件にします。
条件表と実際の部屋の写真です。
実際の自宅エアコンの仕様書には強運転時の風量の記載が無かったため流量値は類似機種の仕様書から算出しました。
上下方向の影響を見たかったのと計算時間の都合で、左右方向の角度は45°で決め打ちしました。
他の条件は固定して上下方向のみを0°, 30°, 45°, 60°の4通りに変化させて30秒後の部屋の冷え方を比較します。
計算環境はwindows11 + WSL2(ubuntu22.04LTS)にopenfoam2312をインストールしたものを使用しました。
手順
言わずと知れたPENGUINITIS様の中の記事、OpenFOAMの概要を参考に部屋の3Dモデルを作るところからスタートしました。
もちろん一発でうまくいくわけもなく、計算が発散してメッシュ作成～計算の実行までを行ったり来たりしてました。
部屋の3Dモデル作成
メッシュ作成
計算条件の設定
計算の実行
結果の可視化と評価
実測結果との比較(未完)
部屋の3Dモデル作成
自室の3DモデルをFreeCADを使って作成していきます。
作成に際し、Xsim様のFreeCAD使い方メモが大変参考になりました。
FreeCADでのモデル作成からopenfoamで読み込めるSTL形式への変換までをこの項で解説します。
これが今回シミュレーション対象となる部屋の写真です。寝室の奥に細長い書斎スペースがある間取りになってます。
入って左手側の壁は一面クローゼットになっており、基本的に引き戸は閉めたままになっているのでこの面は壁として扱うことにしました。
また窓枠やカーテンなど細かな箇所は私のモデリングスキルの無さもあり省略することにしました。
寸法は自宅購入時にもらった間取り図を参考に、またベッドなどの家具の寸法はメジャーで実測しました。
かなり単純なモデルで、ある程度慣れた人なら1時間以内で作れそうです。・・・がFreeCAD初心者の私はだいぶ苦戦しました。
まずは部屋の壁(床+天井)部分をPartDesignベンチのスケッチ→押し出しで作成し、面に分解します。
さらに分解した面にスケッチ→押し出しでエアコン本体やテーブル等を作成していきます。エアコンの吹き出し口と吸い込み口は境界条件設定用に別々の面として保存します。
境界面の設定を細かく行えるように面をまとめておき、面ごとにまずast形式でエクスポートします。
※このときast形式でエクスポートしないとバイナリ形式で保存される様で、エディタで編集できなくなるみたいなので要注意
ast形式で出力したファイルをエディタで直接開き、1行目と末尾行の名前の部分がデフォルトでは"Mesh"となっているので任意の面の名前に変えます。
名前を変えたら拡張子を.ast→.stlに変えておきます。これをすべての面について実施します。
実施出来たらFreeCADですべてのSTLファイルをインポートしてみて正しく表示されているかどうか確認します。
ここまで出来たらメッシュの作成に移ります。
メッシュの作成
作成したSTLファイルからcfMeshを使ってメッシュを作成します。
snappyHexMeshも試しましたが、checkMeshの結果は問題無いにも関わらず計算がcfMeshよりも安定しませんでした。
とはいえ、この辺は私の設定スキル不足によるものが大きいと思います。
このメッシュの作成が今回最も苦戦した部分かもしれません。
作り始めた時はメッシュ品質がそこまで重要だと考えておらず、checkMeshの結果が悪くても気にせずに計算しようとしてました。
もちろんうまく計算できないので計算条件側で何とかならないか模索してましたが当然うまくいくはずもなく、一時はこのシミュレーション自体を諦めかけてました。
趣味とはいえここまでやって諦めるのも何となく嫌だったので色々調べているうちにメッシュから作り直した方が良いことに気付き、checkMesh --AllFields⇒paraviewで可視化 で良くないメッシュを調べてそれが無くなるように作り直しを重ねました。。。
前項のstlファイルの取り扱い方など、主な手順は公式チュートリアル "3 weeks" series Day_6 - Meshing with cfMeshを参考にしました。
メッシュ作成前にこのタイミングでケースディレクトリを作ってしまいます。チュートリアルのheatTransfer\buoyantPimpleFoam\hotRoomをベースにしました。stlディレクトリを作成しその中に作成したstlファイルを全て突っ込んでおきます。
また、paraview用のダミーファイル(.foamファイル)も作成しておきます。
ケースディレクトリは下記のようになります。
myroom_ac_0deg/
├── 0/
│ ├── U
│ ├── T
│ └── ...
├── constant/
│ ├── g
│ ├── thermophysicalProperties
│ └── turbulenceProperties
├── stl/
│ ├── aircon_body.stl
│ ├── blowout_face.stl
│ └── ...
├── system/
│ ├── controlDict
│ ├── fvSchemes
│ ├── fvSolution
│ └── meshDict
└──myroom_ac_0deg.foamまず作成した面ごとのstlファイルを一つのファイルにまとめます。
stlディレクトリ下でcatコマンドで全てのstlファイルをcombined.stlにまとめてしまいます。
$ cat aircon_body.stl blowout_face.stl ・・・ ceiling.stl >> combined.stl openfoamでは[mm]などの単位の情報が失われ、全て[m]として扱われるため必要に応じスケールの変更を行います。
※下記はmm→mに変更する場合
$ surfaceConvert -scale 0.001 combined.stl combined_m.stl今回のケースでは特徴線情報を得るためにsurfaceFeatureEdgesでfmsファイルに変換してからメッシュ生成するよりも、stlから直接メッシュ生成したほうがメッシュ品質が良かったです。
fmsに変換する場合は下記コマンドで変換して、変換後のfmsファイルの先頭部分にある面のtypeをwallやpatchに変更します。（デフォルトは全てempty）
$ surfacefeatureEdges combined_m.stl combined_m.fmsSTLから生成する場合はメッシュ生成後にpolymesh/boundaryファイルを直接編集して吹き出し口/吸い込み口をwall->patchに変更しておきます。
ようやくmeshDictの作成に取り掛かります。
前述の通りトライアンドエラーを繰り返してたどり着いたのでなぜこの数字になったのか？について明確に解説できないのですが、ポイントとしては下記の通りです。
流速の早いエアコン周りの領域のメッシュが細かくなるようにobjectRefinementで設定
形状の変化の大きな境界面も個別設定
それ以外の領域は計算時間削減のため可能な限りメッシュを粗くする
meshDict作成に際してはチュートリアルに加えてcfMeshによるメッシュ作成入門(秋山善克 氏)やcfMeshの公式UserGuideを参考にさせて頂きました。
// meshDict
surfaceFile "./stl/combined_m.stl";
maxCellSize 0.5;
boundaryCellSize 0.05;
boundaryLayers
{
nLayers 1;
thicknessRatio 1.2;
maxFirstLayerThickness 0.01;
}
localRefinement
{
"aircon_body"
{
cellSize 0.01;
}
"blowout_face"
{
cellSize 0.01;
}
"suction_face"
{
cellSize 0.01;
}
"pc_top"
{
cellSize 0.01;
}
"beds"
{
cellSize 0.05;
}
"table"
{
cellSize 0.05;
}
}
objectRefinements
{
// エアコンに最も近い部屋の角を中心とした半径2mの球
blowarea
{
type sphere;
additionalRefinementLevels 3;
centre (6.25 0.0 2.35);
radius 2.0;
}
}cfMeshのcartesianMeshでメッシュを作成。
$ cartesianMesh生成できたらcheckMesh -writeAllFieldsでメッシュ品質を確認しその結果を0ディレクトリに出力。これをparaviewで可視化して非直行性やひずみの大きい部分のメッシュが少なくなるように調整を繰り返しました。
$ checkMesh -writeAllFieldscheckMesh結果の抜粋です。
Checking geometry...
Overall domain bounding box (-6.06539e-06 -3.36573e-05 -0.000146672) (6.25023 3.55001 2.3)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (-1.94175e-16 5.33067e-16 -1.7697e-15) OK.
Max cell openness = 3.31634e-16 OK.
Max aspect ratio = 8.21867 OK.
Minimum face area = 2.54506e-06. Maximum face area = 0.228904. Face area magnitudes OK.
Min volume = 3.28301e-09. Max volume = 0.134939. Total volume = 34.0866. Cell volumes OK.
Mesh non-orthogonality Max: 78.4497 average: 4.76809
*Number of severely non-orthogonal (> 70 degrees) faces: 7.
Non-orthogonality check OK.
<<Writing 7 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 2.42207 OK.
Coupled point location match (average 0) OK.
Mesh OK.非直交性が良くない面はベッドとエアコン下テーブルが隣接している部分に集中して存在していました。
objectRefinementsでこの部分だけを細かくしたり、いろいろ試してみましたがうまくいかず・・・モデルの作り方も悪かったのでしょう。
とりあえず作り直しの目安と言われるMax80°未満に抑えられ、その他項目でも特にエラーが発生していないのでこれで行くことにしました。
なお、総セル数は約31万セルでした。
作成したメッシュをparaviewで確認した結果です。とりあえず意図した通りのメッシュになっているようです。
計算条件の設定
heatTransfer/buoyantPimpleFoam/hotroomのチュートリアルケースをベースに境界条件の設定を行いました。
設定に際し、OpenFOAMのGoogleGroupの事例が参考になりました。
ソルバ： buoyantPimpleFoam
乱流モデル： RANS 標準k-ε
輻射： 設定なし
重力： あり (0 0 -9.81)
Uの設定
自宅エアコンの仕様書には「強」運転時風量の記載が無かったため、類似機種の仕様書をダイキン工業HPの仕様書検索ページより調べました。
「強」運転時に14.5[$${m^3}$$/min]≒0.241667[$${m^3}$$/s]より、この値とモデリングしたエアコンの吹き出し口・吸い込み口の面積から流入・流出風量が同じになるように風速を計算。
（実測時のミスで吹き出し口を実際よりちょっと大きくしてしまったことにだいぶ後で気づきました＼(＾o＾)／）
吹き出し風速
= 風量0.241667[$${m^3}$$/s] / (吹出口横幅 0.64[m] × 吹出口縦幅0.14485[m])
= 2.60686[m/s]
吸い込み風速
= 風量0.241667[$${m^3}$$/s] / (吹出口横幅 0.64[m] × 吹出口縦幅0.22[m])
= 1.71638[m/s]
吸い込み口は-z方向のみの風速(0 0 -1.71638)になりますが、吹き出し口はxyz各方向に分解した値（球面座標）に直しておきます。
0° ・・・ (-1.84333 1.84333 0)
30°・・・ (-1.59637 1.59637 -1.30432)
45°・・・ (-1.30432 1.30432 -1.84333)
60°・・・ (-0.921665 0.921665 -2.25761)
// U(0deg)
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
".*"
{
type noSlip;
}
blowout_face
{
type fixedValue;
value uniform (-1.84333 1.84333 0);
}
suction_face
{
type fixedValue;
value uniform (0 0 -1.71638);
}
}当初はflowRateInletVelocity, flowRateOutletVelocityを用いて設定する予定でしたが、flowRateInletVelocityでは吹き出し角度を変えることができない（方向は面の法線方向で固定？）ようなので断念。また、吸い込み側だけflowRateOutletVelocityで設定してシミュレーションしてみると風速が大きくなりすぎる部分が発生したり発散したりでこれまた断念。
上記の設定でなんとか発散することなく計算できましたが（圧力の残差がやや大きかった・・・）、いま振り返るとやはり自由流入流出面を設けておいた方が良かったと思います。
実際の部屋でもドアの隙間や換気口から空気が出入りしてるはずですし。
なお、シミュレーション中の連続の式の残差sum_localはおおよそ1e-6くらいの値で推移してました。
k, epsilonの設定
公式チュートリアル "3 weeks" series Day_7の資料で紹介されていたk,epsilonの推定式を用いて値の見積もりを行いました。
この辺りの乱流の理論に関わる部分は難しくてあまり理解できなかったです。要復習です。
k：乱流運動エネルギー(∝乱れの大きさ)
$$
k = \frac{3}{2}(流入流速×乱流強度)
$$
ε：乱流散逸率(∝乱れが消えていく速さ)
$$
ε = \frac{0.09^\frac{3}{4}k^\frac{3}{2}}{0.07×代表長さ}
$$
(Cμを0.09, 混合長を0.07×代表長さとした場合)
乱流強度を5%、吹き出し口の短辺長さ0.14885[m]を代表長さとして上式にあてはめ計算すると
k = 0.0254840
ε = 0.0659276
となり、これをもとにk, epsilonファイルを作成しました。
// k
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0.0254840;
boundaryField
{
".*"
{
type kqRWallFunction;
value uniform 0.0254840;
}
blowout_face
{
type turbulentIntensityKineticEnergyInlet;
intensity 0.05;
value uniform 0.0254840;
}
suction_face
{
type zeroGradient;
}
}// epsilon
dimensions [0 2 -3 0 0 0 0];
internalField uniform 0.0659276;
boundaryField
{
".*"
{
type epsilonWallFunction;
value uniform 0.0659276;
}
blowout_face
{
type turbulentMixingLengthDissipationRateInlet;
mixingLength 0.0101395;
value uniform 0.0659276;
}
suction_face
{
type zeroGradient;
}
}Tの設定
部屋の温度を308[K]（≒35[℃]）、エアコン吹き出し口の温度は295[K]（≒22[℃]）、ノートPC上面の温度を313[K]（≒40[℃]）として下記のように設定。
（ノートPCは結局ほとんど結果に影響しませんでした。。。）
// T
dimensions [0 0 0 1 0 0 0];
internalField uniform 308;
boundaryField
{
".*"
{
type zeroGradient;
}
blowout_face
{
type fixedValue;
value uniform 295;
}
pc_top
{
type fixedValue;
value uniform 313;
}
}その他のパラメータの設定
その他のパラメータは基本的にhotroomのチュートリアル通りに設定しました。
constant/thermophysicalProperties, turbulencePropertiesはそのままでOK、constant/gは重力が-z方向にかかるように(0 0 -9.81)に修正しておきます。
修正を忘れてY方向に重力を働かせたまま計算してしまったことがありました＼(＾o＾)／
計算走らせる前に落ち着いて条件の見直しをしておきましょう。
// p_rgh
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 100000;
boundaryField
{
".*"
{
type fixedFluxPressure;
value uniform 100000;
}
}// p
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 100000;
boundaryField
{
".*"
{
type calculated;
value uniform 100000;
}
}// nut
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
".*"
{
type nutkWallFunction;
value uniform 0;
}
blowout_face
{
type calculated;
value uniform 0;
}
suction_face
{
type calculated;
value uniform 0;
}
}// alphat
dimensions [1 -1 -1 0 0 0 0];
internalField uniform 0;
boundaryField
{
".*"
{
type compressible::alphatWallFunction;
value uniform 0;
}
blowout_face
{
type calculated;
value uniform 0;
}
suction_face
{
type calculated;
value uniform 0;
}
}計算の実行
自分のノートPCで6CPU並列計算で1水準あたりおよそ6時間程度かかりました。
contolDict内の設定でdeltaT自動調整機能をonにし、maxCoはやや大きいですが10に設定しました。
$ decomposePar
$ mpirun -np 6 buoyantPimpleFoam -parallel 夜中寝てる間に走らせていたのですが、ノートPCのファンの音がうるさいと妻に怒られたのとCPU温度が90℃近くになってPCの寿命的にもよろしくないなと思い始めたので試しにGoogleCloudのVMインスタンスe2-highcpu-16を使って計算させてみました。自分のPCの半分以下の時間で計算が終わり軽く感動しました。でも若干のお金がかかるので多用は控えたいところ。。。
なお必要に応じ、バックグラウンドで処理を継続させるためにnohupコマンドを付記して実行させた方がよいかもしれません。
$ nohup -np 6 buoyantPimpleFoam -parallel &実行ディレクトリ下のnohup.outファイルに結果が出力されるようになるので、定期的にtail等でファイルの中をのぞいて進捗を確認できます。
コードの中身と支配方程式
こちらにかなり詳しい解説が載っていました。UEqnの方はなんとかイメージが持てるのですがEEqn, pEqnの方はさっぱりです。
ちゃんと理解できるように勉強しておきたいです。
Fumiya Nozaki. CFD WITH A MISSION - buoyantPimpleFoam and buoyantSimpleFoam in OpenFOAM
結果の可視化と評価
いよいよお楽しみの可視化タイムです。
Paraview(5.12.1)がすぐ強制終了して落ちちゃうのが難点。
メモリに十分余裕がある状態でも落ちてました。特にParticleTracerやStreamTracerを使ったときに落ちやすかったです。少し古いバージョン(4.4)を使用することでやや改善しましたがそれでも落ちるときは落ちました。
色々試してるうちに並列処理を行わない1CPUの仮想マシン下だと処理は重いけど落ちずに処理ができることに気づき、仕方なくvirtualBoxで1cpuのubuntu22.04desktop環境を上記の落ちやすいフィルタ使用時のみ使ってます。
setting>general内の各種項目の変更、グラフィックドライバ更新、pvserver版の使用、・・・あたりを試してもだめでした。このparaviewすぐ落ちる問題は良い解決法があれば是非教えてほしいです。
t=0sでの状態の確認
t=0[s]での状態です。部屋全体は308[K]となっています。
全体をwireframe(90%透過)で表示させてます。
t=30sでの比較
t=30[s]でのエアコン吹き出し口からの流線と等温度面を比較しました。
等温度面は304[K]です。
0°では天井部分ばかり冷えて部屋に風が下りてきてないようです。また、一部の冷風が吸い込み口に吸われてしまっている様子も見えます。
30°だとちょうどベッドを超えたくらいの位置に冷風が下りてきてます。部屋の入口側に向かって風が固まっている感じがします。
45°, 60°はベッドに冷風が当たりいい具合に(？)部屋およびベッド上に拡散しています。今回のシミュ結果の中ではこのへんの角度が一番よさそうです。
冷房は水平方向が良いとよく聞くので、0°か30°がいい結果が出ると思ってましたが意外でした。
今回は30秒だけでしたがもう少し長い時間、もしくはbuoyantSimpleFoamで定常状態を見るなどすればまた別の考察ができると思います。
自分が寝る場所の温度の時間変化
paraviewでの時系列プロットの練習としてやってみました。
自分はいつもエアコンから遠い側のベッドで寝ているので、そのベッド中央部の5cm上空（下記赤十字マーク：x=5.25 y=2.65 z=0.50）の温度の時間変化を調べてみました。
もちろん実際はONにしてすぐ冷たい風は出ないわけですが。。。
（そのへんの現実との差異や解析の妥当性はさておき）今回の結果からは、大まかな傾向として自分の寝る場所は0°, 30°より45°, 60°の方が早めに冷えそうだよーと言えそうです。
当たり前と言えばそれまでですが、ベッドに風がぶつかる・ぶつからないでここまで顕著に冷え方が変わってくるのが見えて面白いです。
ParticleTracer動画
ParticleTracerの動画を試しに作ってみました。
流れが動画で見えると楽しいです＾ｑ＾
※左上：0° , 右上：30°,
左下：45°, 右下：60°
実測結果との比較(未完)
前項の自分が寝てる場所の温度の時間変化を実測しようと思ってましたがまだやれてません＼(＾o＾)／
今回無視した要因が多いのと私の未熟さもあり実測と一致することはおそらくないと思いますが、各角度での大まかな傾向くらいは確認したいなーと思ってます。
たしかに体感では水平方向よりも斜め下向きのほうがすぐに涼しくなってる気がしますが（風が当たるからそりゃそうか）、所詮体感なのでちゃんと定量的に見ておきたいです。
感想
たのしかったです。
もう本当にこの一言に尽きます。
久しぶりに夢中で取り組めた気がします。
後から振り返ると「あーこれまずいなー」と思うところだらけですが(震え声)。
今後ものんびり楽しみながら勉強していこうと思います。