【徹底解説】chtMultiRegionFoamの設定ファイルを一つずつ読む（輻射・発熱面つき）
こんにちは（@t_kun_kamakiri）
金属のかたまりの側面にヒーター（ヒートマット）を貼り付けて温めると、ブロックの中はどんな温度分布になり、どのくらい変形（熱膨張）するでしょうか。しかも、ヒーターを途中で切ったら変形はどう戻っていくでしょうか。このシリーズでは、その一部始終をOpenFOAM（熱・流れ・輻射）とFrontISTR（構造）を組み合わせて計算しています。計算対象は、1000mm角の空気の中に置いた200×200×400mmの金属ブロックで、側面の一部に貼ったヒートマットを0〜5分ON（20W）、5〜10分OFFにします。
↑計算対象のモデルです。左がOpenFOAMの流体領域（空気）、右がFrontISTRの構造モデルで、いずれも計算開始時（0秒、一様20℃）の状態です。
実際に計算すると、次のように動きます。
↑左がOpenFOAMで解いた流体（空気）の温度場、右がその温度を受け取ったFrontISTRの構造モデル（変形）です。ヒートマット側から温度が上がって空気が対流し、ブロックがじわじわ反っていく様子が見えます。
本記事はそのシリーズの中でも、OpenFOAM側の設定ファイルを一つずつ読み解く「徹底解説編」です。chtMultiRegionFoam（共役熱伝達＋輻射＋発熱面）のケースについて、各設定ファイルの中身をブロック単位で「このキーワードは何をしているのか」まで潰していきます。自分のケースへ改造するときの辞書として使ってください（計算の全体像と結果は概要編、温度の受け渡しプログラムは連成の詳細編にまとめています）。
- caseSettings／blockMeshDict／topoSetDictの各キーワードの意味
- thermophysicalProperties（thermoType）の読み方と、流体（空気）と固体（鋼材）の違い
- radiationProperties（fvDOM）とturbulenceProperties、controlDictの主要設定
- fvSchemes／fvSolutionで何を指定しているか（輻射項の離散化・線形ソルバー・PIMPLE）
- changeDictionaryDictでの結合境界条件（turbulentTemperatureRadCoupledMixed）と発熱面の設定
caseSettings：寸法と条件を一元管理する
まずsystem/include/caseSettingsです。ここはOpenFOAMの標準ファイルではなく、寸法や発熱条件をまとめておくために自分で用意したファイルです。#includeで他の辞書から読み込み、#eval #{ ... #}でその場で数式評価します。
|
1 2 3 4 5 6 7 8 9 10 11 12 |
domainLx 1.0; // 流体ドメイン X [m] solidLx 0.2; // 固体 X [m] solidLz 0.4; // 固体 Z（高さ）[m] cellSize 0.05; // メッシュ1辺の目安 [m] heaterWattage 20.0; // ヒーター出力 [W] heaterOnTime 300; // ON継続 [s] heaterMatWidth 0.10; // ヒートマット幅 [m] heaterMatHeight 0.20; // ヒートマット高さ [m] nx #eval #{ round($domainLx / $cellSize) #}; // 格子数を自動計算 solidX0 #eval #{ 0.5 * ($domainLx - $solidLx) #}; // 固体を中央に配置 |
#eval #{ round(...) #}：辞書の中で四則演算や関数を使える。ここでは寸法÷セルサイズで格子数を出している$domainLx：同じファイル内や#include元で定義した変数を参照する記法
これで、cellSizeを変えるだけで格子数が、solidLxを変えるだけで固体の位置が自動追従します。
blockMeshDict：背景メッシュ
最初は固体も流体も区別しない、1つの直方体メッシュだけを作ります。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
#include "include/caseSettings" convertToMeters 1; // 座標の単位はm vertices // 直方体の8頂点 ( (0 0 0) ($domainLx 0 0) ($domainLx $domainLy 0) (0 $domainLy 0) (0 0 $domainLz) ($domainLx 0 $domainLz) ($domainLx $domainLy $domainLz) (0 $domainLy $domainLz) ); blocks ( hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1) ); |
#include "include/caseSettings"：先ほどのcaseSettingsファイルをここで読み込む行。これにより、$domainLxや$nxといった変数がこのblockMeshDict内でも使えるようになる（寸法をcaseSettings側で一元管理できるのはこのおかげ）convertToMeters 1：頂点座標をそのままメートルとして扱う（mmで書くなら0.001にする）hex (0..7)：8頂点で1つの六面体ブロックを定義($nx $ny $nz)：各方向の分割数（caseSettingsから）simpleGrading (1 1 1)：格子の粗密比。1 1 1で等間隔
blockMeshを実行しParaViewで確認すると以下となります。
topoSetDict：3つのcellZoneへ分ける
1つのメッシュを、位置で指定した箱の中のセルごとにcellZoneへ登録します。ここが今回のケースの肝です。
|
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 |
// 背景メッシュ(region0)を、固体直方体・ヒートマット・流体の3つの // cellZoneに分割する(互いに排他的である必要がある。splitMeshRegions // -cellZones はメッシュ上の全cellZoneを領域分割に使うため、重複があると // FATAL ERRORになる)。 // heaterMatは固体の+X面の一部を占める薄い直方体で、solidの範囲からは // 差し引く(action subtract)ことで、solid/heaterMatが重ならないようにする。 // box座標は system/include/caseSettings と共有しており、blockMeshDict の // セル境界と厳密に一致する(cellSizeがdomain/solid/heaterMatの寸法の // 整数分の1のため)。 #include "include/caseSettings" actions ( // heaterMat: 固体+X面の一部、厚さ1セル分の薄い直方体 { name heaterMatCellSet; type cellSet; action new; source boxToCell; box ($heaterX0 $heaterY0 $heaterZ0) ($heaterX1 $heaterY1 $heaterZ1); } { name heaterMat; type cellZoneSet; action new; source setToCellZone; set heaterMatCellSet; } // solid: 固体直方体全体からheaterMatを除いたセル { name solidCellSet; type cellSet; action new; source boxToCell; box ($solidX0 $solidY0 $solidZ0) ($solidX1 $solidY1 $solidZ1); } { name solidCellSet; type cellSet; action subtract; source cellToCell; set heaterMatCellSet; } { name solid; type cellZoneSet; action new; source setToCellZone; set solidCellSet; } // fluid: 固体直方体(heaterMat含む)の外側の全セル { name fluidCellSet; type cellSet; action new; source boxToCell; box ($solidX0 $solidY0 $solidZ0) ($solidX1 $solidY1 $solidZ1); } { name fluidCellSet; type cellSet; action invert; } { name fluid; type cellZoneSet; action new; source setToCellZone; set fluidCellSet; } ); |
boxToCell：指定した直方体box (最小点) (最大点)の中に中心があるセルを選ぶcellZoneSet+action new：選んだセル集合を新しいcellZoneにするaction subtract：既存ゾーンから引く。3ゾーンが重ならないようにするための必須処理
このあとsplitMeshRegions -cellZones -overwriteで本当に別リージョンへ分割すると、solid_to_fluidやheaterMat_to_solidといった結合パッチが自動生成されます。
変数はすべて#include “include/caseSettings”でセットされているので、寸法を変更したい場合はこちらの数値を変更します。
topoSetコマンドを実行すると以下のようにsetsフォルダができており、こちらにセルゾーンなどが作られています。
|
1 2 3 4 5 6 7 8 |
tree -L 1 constant/polyMesh/sets/ constant/polyMesh/sets/ ├── fluid ├── fluidCellSet ├── heaterMat ├── heaterMatCellSet ├── solid └── solidCellSet |
うまくセルゾーンができているかをParaViewで確認することができます。
このtopoSetを実行することで各流体領域と固体領域の下地となるセルゾーンができました。
ただ、このままでは領域は分割されておらず（blockMeshのまま）次のsplitMeshにより領域分割する必要があります。
regionProperties：どのゾーンを固体／流体として解くか
|
1 2 3 4 5 6 |
// constant/regionProperties regions ( fluid (fluid) solid (solid heaterMat) // ヒートマットも固体扱い ); |
splitMeshRegionsで分けた各cellZoneを、流体ソルバーで解くか固体ソルバーで解くかを振り分けるファイルです。fluidグループにはfluidを、solidグループにはsolidとheaterMatを入れています。ヒートマットは発熱する部品ですが、扱いとしては熱伝導だけの固体なので、固体側に含めています。
境界条件の設定と領域分割
マルチリージョンソルバの境界条件の設定手順はとても面倒ですが、こちらの方法をお勧めします。
まずは、0フォルダに以下のファイルを用意し
|
1 2 3 4 5 6 7 8 9 10 11 12 |
tree -L 1 0.orig 0 ├── G ├── IDefault ├── T ├── U ├── cellToRegion ├── fluid ├── heaterMat ├── p ├── p_rgh ├── qr |
中身をすべて以下とします。
|
1 2 3 4 5 6 7 8 |
boundaryField { ".*" { type calculated; value $internalField; } } |
それをコピーして0フォルダを作ります。
|
1 |
cp -r 0.orig 0 |
splitMeshRegions ：領域分割
領域分割には以下のコマンドを実行します。
|
1 |
splitMeshRegions -cellZones -overwrite |
これによりセルゾーンから完全にメッシュの領域分割ができました。
さらに0フォルダの境界条件ファイルも自動で作ってくれます。
|
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 |
tree -L 2 0 0 ├── G ├── IDefault ├── T ├── U ├── cellToRegion ├── fluid │ ├── G │ ├── IDefault │ ├── T │ ├── U │ ├── cellToRegion │ ├── p │ ├── p_rgh │ └── qr ├── heaterMat │ ├── G │ ├── IDefault │ ├── T │ ├── U │ ├── cellToRegion │ ├── p │ ├── p_rgh │ └── qr ├── p ├── p_rgh ├── qr └── solid ├── G ├── IDefault ├── T ├── U ├── cellToRegion ├── p ├── p_rgh └── qr |
境界条件の設定
次にsystem/<領域>/changeDictionaryDictにそれぞれの境界条件をしておきます。
例えば
|
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 |
U { internalField uniform (0 0 0); boundaryField { ".*" { type noSlip; } } } T { internalField uniform 293.15; boundaryField { ".*" { type fixedValue; value uniform 293.15; } "fluid_to_.*" { type compressible::turbulentTemperatureRadCoupledMixed; Tnbr T; kappaMethod fluidThermo; kappa none; qrNbr none; qr qr; value uniform 293.15; } fluid_to_heaterMat { // heaterMatの外側面はmode=powerの熱量境界条件(ケーシングで // 断熱されている想定)としたため、fluid側はここでは // 単純な断熱壁として扱う(solid_to_fluid等とは異なりTnbrで // 結合しない。文字列完全一致でfluid_to_.*正規表現より優先)。 type zeroGradient; } } } p_rgh { internalField uniform 1e5; boundaryField { ".*" { type fixedFluxPressure; value uniform 1e5; } } } p { internalField uniform 1e5; boundaryField { ".*" { type calculated; value uniform 1e5; } } } G { internalField uniform 0; boundaryField { ".*" { type calculated; value uniform 0; } } } qr { internalField uniform 0; boundaryField { ".*" { type calculated; value uniform 0; } } } IDefault { internalField uniform 0; boundaryField { ".*" { type wideBandDiffusiveRadiation; value uniform 0; } } } |
※各領域で適切な境界条件設定しておいてください。
このようにしつつ、以下のコマンドを実行します。
ひとつずつ実行してもいいですし、
|
1 2 3 |
changeDictionary -region fluid changeDictionary -region solid changeDictionary -region heaterMat |
for分にしても良いです。
|
1 2 |
方法B：forループで1行（Allrun.preと同じ方式） for r in fluid solid heaterMat; do changeDictionary -region $r; done |
設定のまととめ
以下のコマンドをスクリプトにしておくと便利でしょう。
Allrun.pre
|
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 |
#!/bin/sh cd "${0%/*}" || exit # Run from this directory . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions #------------------------------------------------------------------------------ runApplication blockMesh runApplication topoSet # Restore initial fields (0.orig -> 0) restore0Dir runApplication splitMeshRegions -cellZones -overwrite # Remove fluid-only / radiation-only fields from the solid-family regions # (solid + heaterMat, per constant/regionProperties) # (avoid brace expansion: /bin/sh may be dash, which does not support it) for region in $(foamListRegions solid) do for f in nut alphat epsilon k U p_rgh G IDefault qr do rm -f 0/$region/$f rm -f processor*/0/$region/$f done done for region in $(foamListRegions) do runApplication -s $region changeDictionary -region $region done echo echo "End" #------------------------------------------------------------------------------ |
次回からは以下のコマンドでメッシュ作成から解析設定まで（境界条件の設定まで）を行ってくれます。
|
1 |
./Allrun.pre |
g：重力の向き（自然対流を駆動する）
|
1 2 3 |
// constant/g dimensions [0 1 -2 0 0 0 0]; // 加速度の次元 [m/s^2] value (0 0 -9.81); // -Z方向（高さ方向の下向き）に重力 |
今回は空気の温度差による浮力（自然対流）を計算するため、重力が必要です。ヒートマットで温められた空気は軽くなって上昇し、周囲の冷たい空気が入れ替わりに降りてくる、という流れがこのgによって生まれます。高さ方向を$z$にとっているので、下向きの(0 0 -9.81)を与えています。もし重力を与えないと浮力が働かず、空気の対流が起きない（熱が伝導と輻射だけで伝わる）計算になってしまいます。
thermophysicalProperties：熱物性の指定方法
熱物性はthermoTypeで「どんなモデルの組み合わせで熱を計算するか」を指定し、mixtureで具体的な数値を与えます。まず流体（空気）です。
|
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 |
// constant/fluid/thermophysicalProperties thermoType { type heRhoThermo; // 密度ベースの熱物性 mixture pureMixture; // 単一成分 transport const; // 粘性・熱伝導は一定 thermo hConst; // 比熱Cp一定 equationOfState perfectGas; // 理想気体（浮力に必要） specie specie; energy sensibleEnthalpy; } mixture { specie { molWeight 28.9; // 空気の分子量 } thermodynamics { Cp 1000; // 比熱 [J/(kg K)] Hf 0; } transport { mu 1.8e-05; // 粘性係数 Pr 0.7; // プラントル数 } } |
ポイントはequationOfState perfectGasです。理想気体にすることで温度によって密度が変わり、自然対流（浮力）が生まれます。次に固体（鋼材）です。
|
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 |
// constant/solid/thermophysicalProperties thermoType { type heSolidThermo; // 固体用 transport constIso; // 等方な熱伝導 thermo hConst; equationOfState rhoConst; // 密度一定 ... } mixture { transport { kappa 50; // 熱伝導率 [W/(m K)] } thermodynamics { Cp 460; // 比熱 [J/(kg K)] Hf 0; } equationOfState { rho 7850; // 密度 [kg/m3] } } |
流体との違いは、heSolidThermo（固体用）、constIso（等方熱伝導率kappaで指定）、rhoConst（密度一定）です。ヒートマットも同じ鋼材相当にしています（constant/heaterMat/thermophysicalProperties）。
turbulenceProperties：乱流モデル
乱流は今回simulationType laminar（層流）としています。
|
1 |
simulationType laminar; |
乱流モデルの設定はすべての流量領域に設定する必要があります。
固体領域には設定は不要です。
radiationProperties：輻射（fvDOM）
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
// constant/fluid/radiationProperties radiation on; radiationModel fvDOM; // Discrete Ordinates Method fvDOMCoeffs { nPhi 3; // 方位角（周方向）の分割数 nTheta 5; // 極角（天頂方向）の分割数 tolerance 1e-3; // 輻射反復の収束判定 maxIter 10; } solverFreq 10; // 10回の流れ反復ごとに輻射を1回解く absorptionEmissionModel constantAbsorptionEmission; constantAbsorptionEmissionCoeffs { absorptivity absorptivity [ m^-1 ] 0.01; // 空気はほぼ透明 emissivity emissivity [ m^-1 ] 0.01; } |
nPhi × nTheta：輻射の方向を何本の離散方向で近似するか。増やすほど正確だが重くなるsolverFreq：毎ステップ輻射を解くと重いので、何回かに1回にする間引きabsorptivity/emissivityを0.01：空気は輻射をほぼ吸収・放射しない（面同士の輻射交換が主役）
各壁面の輻射率（グレー拡散面としての扱い）はconstant/fluid/boundaryRadiationPropertiesで個別に、固体表面0.85、外周壁0.9などと設定しています。固体側はradiation off（内部では輻射を解かない）です。
fvDOMは何を解いているのか（数式）
fvDOMは放射伝達方程式（RTE: Radiative Transfer Equation）を解いています。ある位置 $\mathbf{r}$ で、ある方向 $\mathbf{s}$ に進む放射強度を $I(\mathbf{r},\mathbf{s})$ とすると、散乱の無いグレー媒体では次の式が成り立ちます。
\mathbf{s}\cdot\nabla I(\mathbf{r},\mathbf{s})
=
\kappa_a
\left(
\frac{\sigma T^4}{\pi}
–
I(\mathbf{r},\mathbf{s})
\right)
\end{align*}
左辺は「その方向へ進むと強度がどれだけ変化するか」、右辺は「その場の温度で放射される分（$\sigma T^4/\pi$、$\sigma$はステファン・ボルツマン定数）から、吸収される分（$I$）を引いたもの」です。$\kappa_a$ は吸収係数で、radiationPropertiesのabsorptivityに対応します。
この式は「方向 $\mathbf{s}$」という連続的な向きについての方程式なので、そのままでは解けません。そこで向きを有限個の代表方向 $\mathbf{s}_i$ に離散化するのがDOM（離散座標法）です。各方向について、強度 $I_i$ の移流方程式を1本ずつ解きます。
\mathbf{s}_i\cdot\nabla I_i
=
\kappa_a
\left(
\frac{\sigma T^4}{\pi}
–
I_i
\right),
\qquad
i=1,\dots,N
\end{align*}
この方向の本数 $N$ が、まさにfvDOMCoeffsのnPhiとnThetaで決まります。OpenFOAMでは
N
=
4\,n_\phi\,n_\theta
\end{align*}
で、今回は $n_\phi=3$、$n_\theta=5$ なので $N=4\times3\times5=60$ 方向です（実際、実行ログにも「fvDOM : Allocated 60 rays」と出ます）。1方向あたり1本の移流方程式なので、この左辺の離散化がfvSchemesのdiv(Ji,Ii_h)、各 $I_i$ を解く線形ソルバーがfvSolutionのIiに対応します。
全方向の強度を足し合わせた（立体角 $\omega_i$ で重み付き和を取った）ものが入射放射 $G$ です。これがfvSolutionに出てきたGフィールドの正体です。
G
=
\int_{4\pi} I\, d\Omega
\;\approx\;
\sum_{i=1}^{N} \omega_i I_i
\end{align*}
最終的に、エネルギー方程式（温度の式）へ加わる輻射による発熱・吸熱は、放射で出ていく分 $4\sigma T^4$ と、入ってくる分 $G$ の差で決まります。
\nabla\cdot\mathbf{q}_r
=
\kappa_a
\left(
4\sigma T^4
–
G
\right)
\end{align*}
壁面では、面から出ていく強度が「自分が温度で放射する分」と「入ってきた放射を反射する分」の和になります。グレー拡散不透明面（opaqueDiffusive）なら、放射率 $\epsilon$ と反射率 $1-\epsilon$ を使って次のように書けます（$q_{in}$ は入射熱流束）。
I_w
=
\frac{1}{\pi}
\left[
\epsilon\,\sigma T_w^4
+
(1-\epsilon)\,q_{in}
\right]
\end{align*}
この $\epsilon$ がboundaryRadiationPropertiesで面ごとに設定した放射率（固体表面0.85、外周壁0.9）です。なお、これらの輻射方程式を毎ステップ解くと重いので、solverFreq 10で「10回の流れ反復ごとに1回だけ輻射を解く」と間引いています。まとめると、radiationProperties・fvSchemesのdiv(Ji,Ii_h)・fvSolutionのIi／G・boundaryRadiationPropertiesは、すべてこの一連のRTE計算のための設定というわけです。
controlDict
時間や保存の設定はcontrolDictです。
|
1 2 3 4 5 6 7 8 9 |
// system/controlDict application chtMultiRegionFoam; endTime 600; // 10分 deltaT 0.01; // 初期の時間刻み writeControl adjustableRunTime; // 実時間で区切って保存 writeInterval 5; // 5秒ごとに保存 adjustTimeStep yes; // 時間刻みを自動調整 maxCo 0.6; // 流体のクーラン数上限 maxDi 10.0; // 拡散数の上限 |
adjustTimeStep yesにすると、maxCo（クーラン数）とmaxDi（拡散数）を超えないようdeltaTが自動で調整されます。安定性を気にせず回せる一方、保存時刻を揃えたいのでwriteControlはadjustableRunTime（5秒ちょうどで保存）にしています。この5秒刻みの保存時刻が、そのままFrontISTRへ渡す時刻になります。
fvSchemes：離散化スキーム
各項をどう離散化するかの指定です。流体側の要点だけ抜き出します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
// system/fluid/fvSchemes ddtSchemes { default Euler; // 時間：1次陰的 } divSchemes { div(phi,U) Gauss upwind; // 移流：風上 div(phi,h) Gauss upwind; // エンタルピー移流 div(Ji,Ii_h) Gauss linearUpwind grad(U); // ★輻射強度の移流 } laplacianSchemes { default Gauss linear corrected; // 拡散 } snGradSchemes { default corrected; } |
div(Ji,Ii_h)が輻射（fvDOM）を使うときに必要な項です。これが無いと輻射計算でエラーになります。固体側のfvSchemesはもっと単純で、熱伝導に相当するlaplacian(alpha,h)だけを指定します（移流が無いため）。
fvSolution：線形ソルバーとPIMPLE
各変数を解く線形ソルバーと、反復制御（PIMPLE）を指定します。
|
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 |
// system/fluid/fvSolution solvers { p_rgh { solver GAMG; smoother GaussSeidel; tolerance 1e-7; relTol 0.01; } "(U|h)" // Uとhをまとめて同じ設定で解く { solver PBiCGStab; preconditioner DILU; tolerance 1e-7; relTol 0.1; } Ii // 輻射強度（fvDOM） { solver GAMG; smoother symGaussSeidel; tolerance 1e-4; maxIter 10; } G // 入射放射（Iiと同じ設定を流用） { $Ii; } } PIMPLE { momentumPredictor yes; nCorrectors 2; // 圧力補正の回数 pRefCell 0; // 圧力の基準（閉じた領域なので必要） pRefValue 1e5; } |
p_rgh（浮力を分離した圧力）はGAMG（マルチグリッド）で解くIi／Gは輻射（fvDOM）の変数。輻射を使うときに必要pRefCell／pRefValue：外に開いていない閉じた領域では圧力の基準点が必要
固体側は温度（エンタルピーh）だけを解くので、fvSolutionもhのソルバー1つとPIMPLEの非直交補正だけの、とても短い内容です。
【補足】changeDictionaryDict：境界条件を確定させる
splitMeshRegions直後の境界条件は仮の状態なので、リージョンごとにchangeDictionaryで本来の条件へ書き換えます。まず流体側の温度Tです。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
// system/fluid/changeDictionaryDict の T ".*" // 外周壁は20℃固定 { type fixedValue; value uniform 293.15; } "fluid_to_.*" // 固体との結合面 { type compressible::turbulentTemperatureRadCoupledMixed; Tnbr T; // 相手側の温度と結合 kappaMethod fluidThermo; qr qr; // 輻射熱流束を渡す } fluid_to_heaterMat // ヒートマット外側＝断熱 { type zeroGradient; } |
結合面に使うcompressible::turbulentTemperatureRadCoupledMixedが、共役熱伝達の中核です。両側の温度・熱流束が釣り合うように解かれます。主なキーワードは次の通りです。
Tnbr T：相手リージョンの温度フィールド名（隣と温度をやり取り）kappaMethod：熱伝導率をどこから取るか（fluidThermo／solidThermo）qr／qrNbr：輻射熱流束を相手と受け渡すためのフィールド名
固体側は、底面などはzeroGradient（断熱）、結合面は同じくturbulentTemperatureRadCoupledMixedです。そしてヒートマットの外側面にだけ、時間変化する熱量境界条件を与えます。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
// system/heaterMat/changeDictionaryDict の T heaterMat_to_fluid { type externalWallHeatFluxTemperature; mode power; // 熱量[W]で与える Q table ( (0 20) (300 20) // 0〜300秒: 20W (302 0) (600 0) // 302秒以降: 0W ); kappaMethod solidThermo; value uniform 293.15; } |
externalWallHeatFluxTemperatureをmode powerで使うと、面に総熱量[W]を与えられます。Qを(時刻 値)のテーブルにすることで、矩形波状のON/OFFが表現できます。
温度測定
温度測定用にprobeを設定しておきます。
probeの座標もinclude/caseSettingsで設定しています。
座標位置はParaViewで確認すると良いでしょう。
solid/probes
|
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 |
#include "../include/caseSettings" solidSurfaceProbes { type patchProbes; libs (sampling); region solid; patches (solid_to_fluid floor); fields (T); writeControl timeStep; writeInterval 1; probeLocations ( ($solidX1 $heaterYc $heaterZc) // P1: ヒートマット近傍(solid側、+X面のノッチ縁) ($solidX0 $heaterYc $heaterZc) // P2: 反対側面(-X面) ($solidXc $heaterYc $solidZ1) // P3: 上面中央 ); } heaterMatSurfaceProbes { type patchProbes; libs (sampling); region heaterMat; patches (heaterMat_to_fluid); fields (T); writeControl timeStep; writeInterval 1; probeLocations ( ($solidX1 $heaterYc $heaterZc) // P0: ヒートマット外側面(発熱面そのもの) ); } |
これらの温度時刻歴をグラフにすると以下になります。
グラフの読み方
- 上段は温度そのもの、下段は全領域共通の初期条件293.15 K（20℃）からの温度変化です。
- P0は発熱面なので加熱中に最も高温になり、300秒のヒーター停止後に低下します。
- P1は固体側なのでP0より遅れて上昇し、停止直後もしばらく上昇します。
- P2・P3の変化が小さいことから、今回の時間範囲では熱が反対側や上面まで十分に広がっていないと読めます。
- 300秒直後のP0の急な温度低下は、発熱量を20 Wから0 Wへ段階的に切り替えた条件に対応します。
P0だけ約0.5 ℃高い位置から始まる理由
P0はヒートマットのセル内部ではなく、externalWallHeatFluxTemperatureで20 Wを与える外側境界面の温度です。初期内部温度は全領域で293.15 Kですが、境界条件は最初の計算時刻から熱流束を満たす壁面温度を作ります。ヒーター面積は$0.10\times0.20=0.020\,\mathrm{m^2}$なので、熱流束は$q”=20/0.020=1000\,\mathrm{W/m^2}$です。境界面から隣接セル中心までの距離$d=0.025\,\mathrm{m}$、熱伝導率$k=50\,\mathrm{W/(m\,K)}$より、
$$T_{wall}-T_{cell}\approx\frac{q”d}{k}=\frac{1000\times0.025}{50}=0.5\,\mathrm{K}$$
となります。このため最初の記録はP0が約293.65 K、その他が約293.15 Kです。これは初期条件をP0だけ高く設定したためではなく、20 Wの熱流束境界条件により壁面とセル中心の間に必要な温度勾配によりセル間の保管温度として境界面温度を出力しているためです。P0をヒートマットの代表温度として使う場合は、境界面温度であることに注意してください
まとめ
chtMultiRegionFoamのケースを、設定ファイル単位で読み解いてきました。改造の勘所を整理すると次の通りです。
- 寸法・条件は
caseSettingsに集約し、blockMeshDict／topoSetDictから#includeで共有する - 領域は
topoSet（boxToCell＋subtractで非重複）→splitMeshRegionsで分割する - 浮力を出すため流体は
perfectGas、輻射はfvDOM＋div(Ji,Ii_h)＋Ii/Gソルバーがセット - 結合面は
turbulentTemperatureRadCoupledMixed、発熱面はexternalWallHeatFluxTemperature（mode=power＋テーブル）
OpenFOAMの温度分布をFrontISTRへ渡す部分（座標合わせ・逆距離加重補間・入力生成）は、連成の詳細編にまとめています。あわせて読んでいただければと思います。
Follow @t_kun_kamakiri
