熱移動

目次

はじめに

このチュートリアルでは,熱移動のモデリングを紹介する.熱移動解析を行うのに必要な支配方程式および境界条件を導出し,説明する.

熱移動はエネルギーの動きを対象とする熱工学の一分野である.熱移動の原動力となるのは温度差である.温度差はシミュレーション領域内あるいはその境界上の異なる現象から発生し,熱伝導熱対流熱放射に分類することができる.すべての効果を組み合せて,与えられた範囲で時間とともに発生する温度場の変化を熱方程式でモデル化する.

モデル化の過程はNDSolveで解くことのできる偏微分方程式(PDE)になる.さらに,このチュートリアルではさまざまなタイプの熱源の他,利用可能な熱的境界条件を使って多様な現実世界の熱相互作用をモデル化する方法の概要も紹介する.

熱移動PDEの正確さと有効性は,熱移動検証テストという別のノートブックで検証してある.

熱移動モデリングの拡張例題はモデルコレクションで見ることができる.

このノートブックで示されているシミュレーション結果のアニメーションの多くはRasterizeを呼び出して生成される.これは,このノートブックに必要なディスクスペースを減らすためである.この呼出しを使った場合の欠点は,使わなかったときほどくっきり見えないということである.

高品質のグラフィックスが得たい場合は,Rasterizeの呼出しを削除するかコメントアウトするかするとよい:

このチュートリアルで使われる記号と対応する単位は用語集のセクションにまとめてある.

このノートブックは,PDEモデルの解法におけるさまざまな面で有限要素法の機能を利用する.

有限要素パッケージをロードする:

熱方程式

熱方程式入門

エネルギー保存の法則から導出される熱方程式(1)は熱伝導媒体内の時間依存熱流束をモデル化するために使われる.

熱方程式の従属変数は温度 であり,時間 ,位置 とともに変化する.偏微分方程式(PDE)モデルは,熱エネルギーが密度 ,比熱容量 の媒体の中で,時間の経過とともにどのように移動するかを記述するものである.比熱容量とは,単位質量の物質の温度を1ケルビン上げるために必要な熱エネルギー量を指定する材料特性である.

PDEは時間微分の他にいくつかの成分でできている.まず,熱伝導率 拡散項がある.熱伝導率またはその他の数量は温度 に大きく依存することがある.これは結果として非線形熱方程式になる.

2つ目の部分は,内部の熱対流をモデル化するための流速 対流項 である.この項は,媒体で内部流れが可能な場合のみ存在する.シミュレーション媒体が固体の場合はこの項はゼロである.

は領域内の熱源を表す.これは熱源タイプのセクションで説明する.

熱方程式の導出

熱方程式を導出するために,エネルギー保存から始める.単位体積領域内で生成されたエネルギーと,領域の境界から流れてくるエネルギーの均衡を考える.

17.gif

上の図で は質量密度(単位:), は単位質量当たりの内部エネルギー(単位:)である.検査体積内の総エネルギーは積 に等しい.中央の赤い円は熱源 (単位:)を表す.力は単位時間当たりのエネルギーに等しい()ので, は領域内部で単位時間あたりに生成される熱エネルギーを示す.熱流束 (単位:) は となる境界に存在する正味エネルギー率を表す.

領域内のエネルギー均衡は次の方程式で書くことができる:

つまり,総エネルギー の時間変化率は領域内で単位時間当りに生成されたエネルギー から領域 に存在する正味エネルギー率を引いたものに等しい.

ここでエネルギー流束 は対流項 と拡散項 に分けることができる.対流項は媒体内部の可能な流束によって移動されたエネルギーを表し,流速度 に比例する:

熱移動が固形媒体で起る場合,固体は定義上内部速度場 が持てないため,対流項は に設定される.

拡散項はエネルギー勾配の結果できたエネルギー流束を表し,エネルギー拡散率 に比例する:

熱伝達のモデリングでは,拡散項(2)は温度勾配 の形式で表されることもよくある.これはフーリエ熱伝導法則として知られている:

ここでエネルギー拡散率 は熱伝導率 (単位:)で表され,マイナス符号は熱拡散が温度低下の方向であることを示している.

拡散成分は媒体のタイプに関わらず常に存在する.

拡散項の主な特性は平滑化である.これについては「拡散方程式の平滑化特性」のセクションで説明する.

(3)と(4)をエネルギー平衡方程式(5)に挿入すると,以下が得られる:

または

上の方程式はどちらも連続および不連続の密度/速度場について成り立つ.この詳細は付録の「不連続なPDE係数を持つ保存法則」でご一読いただきたい.

ここでは単位体積 の領域を考えているので,領域内の総質量は に等しい.したがって,左辺の項 は質量保存方程式として解釈することができる.特に時間微分項は時間の経過に従った単位時間当たりの領域の質量の累積(または損失)として理解できる一方,発散項の部分は質量の流入率と流出率の差を示す.

質量流入は領域の中の質量流出率に質量の累積/損失を加えたものに等しい.内部の質量生成/破壊がない場合,項 は合計するとゼロになり方程式(6)から削除され,次のようになる:

これは熱移動モデルの使用に対する制約になる.方程式(7)および結果の熱移動PDE(8)と(9)は,媒体内の質量が変化するときは使うことができないのである.

内部エネルギー は温度 に依存するため,方程式(10)は連鎖律を使って次のように書き換えることができる:

ここで項 は比熱容量 とも言われる.比熱容量は結果となる温度変化に対する,領域に加えられた/から除かれたエネルギーの比を表す.この定義を使うと,方程式(11)は熱方程式に簡約される:

流れのある媒体内の熱移動のPDE.

一般的な熱方程式は,領域内でのエネルギー保存を表し,熱移動モデルの温度場 について解くのに使うことができる.この方程式は対流項と拡散項の両方を含むので,方程式(12)は対流拡散方程式とも呼ばれる.

前の熱方程式(13)は非保存形式で書かれている.つまり, 項の質量密度 と速度 が発散演算子の外側にあるのである.対流項は発散演算子の内部にあると想像できるかもしれないが, はともに空間依存である可能性があるため,一般の熱方程式(14)は簡単に保存形式に変換できない.

しかし固体媒体では,内部速度場 はゼロに設定され,支配PDEは純粋な熱伝導方程式に簡約される.

固体媒体における熱移動のPDE.

熱方程式は円筒座標および球面座標で表すこともできる.詳細は付録の「熱方程式の特殊な場合」で説明する.

熱移動モデルの設定

熱移動モデルの空間項を直交座標で表す関数を定義すると,熱方程式の設定がより便利になる.

熱移動モデルに必要な入力は以下の通りである.

温度変数
空間依存変数
速度場
熱伝導
密度
熱容量
熱源/ヒートシンク
1Dの定常熱移動モデルを設定する:
2Dの定常熱移動モデルを設定する:
1Dの時間依存熱移動モデルを設定する:

このモデルの定義では非アクティブなPDE演算子が使われている.非アクティブな演算子については「偏微分方程式の数値解法」に記載されている.

モデルパラメータの設定

このノートブックの例題では次のモデルパラメータが使われる.これらのパラメータはシミュレーション領域,シミュレーション終了時間 ,媒体の熱特性を定義する.

領域とシミュレーション終了時間 のモデルパラメータを設定する:
方程式の特定の材料パラメータを使うために,ThermodynamicDataから関連性のあるデータを抽出する:

例によっては,熱流束 や表面温度 等の過渡パラメータに対する時間プロファイルを指示するために滑らかな階段関数 を使う場合がある.

ここで関数が到達できる最小値と最大値はそれぞれ で表す.階段の位置は で,滑らかな傾きは で制御される.

滑らかな階段関数 を定義して可視化する:

基本的な熱移動の例

次の2D定常の例[15]は熱移動モデリングの通常のワークフローを示す.

97.gif

,高さ のこのモデル領域は,高熱伝導材料に埋め込まれているセラミックストリップである.このストリップの側面は一定温度 に維持されている.ストリップの上面は の周囲環境への熱対流熱放射で熱を失っている.しかし底面は断熱されていると想定される.

ここでの目的は,セラミックストリップの定常状態温度分布を求めることである.

,高さ の矩形領域を設定する:

セラミックストリップの熱伝導率 ,熱移動係数 ,密度 ,熱容量 , 放射率 は以下で与えられる:

モデルの変数とパラメータを定義する:
左右の境界に温度境界条件を設定する:
上面に対流境界条件を設定する:
上面に熱放射境界条件を設定する:

残った底面領域には,デフォルトの断熱境界条件が適用される.

熱移動PDEモデルを解く:

流束は方程式の一部である.この理由は偏微分方程式と境界条件のセクションで説明する.

定常状態温度分布を可視化する:

定常状態では,熱移動値でモデル化された放射と対流により冷やされる上面が最低温度になる.熱は両側から媒体内に拡散するため,側面に定義された温度が最大である.

熱的境界条件の設定は以下の「熱移動の境界条件」で詳しく説明する.

焼きばめも時間非依存の応用例である.時間依存の応用例としてはレーザービーム溶接がある.

熱源タイプ

領域内の熱生成()または熱吸収()をモデル化するためには,熱方程式(16)の熱源項 を使う.熱源は,形状に基づき体積熱源層熱源点熱源に分類される.

メッシュが熱源項 の幾何学的境界条件に準拠するようにすることが大切である.このための最善の方法は,幾何学的境界条件のためのメッシュを明示的に生成することである.解法として使われる有限要素法は,それぞれの要素が材料境界を横断しないメッシュで有効である.別の方法としてMeshRefinementFunctionを使うこともできる.この場合,要素は材料境界を横断するが,メッシュが十分繊細であれば要素は小さい.両方のアプローチを組み合せることも可能である.

モデル変数とモデルパラメータを設定する:

熱源は時間 ,空間 ,温度 等の他の従属変数の関数となることができる.

体積熱源

体積熱源 は領域内の任意の形状の熱源()またはヒートシンク()をモデル化するために使うことができる.体積熱源 の単位はである.対応する熱源の強さ は単位体積に対する内部加熱または冷却の割合を表す.

体積熱源の値は常に単位で指定される.この名前は熱源の3D表現からきているが,他の次元でも使われる.

である1Dの場合を考える.体積熱源 を単位で指定すると,その値は単位の断面積 で乗算され,最終的な単位はになる.

である2D領域の場合も同様で,単位の体積熱源 は単位の厚さ で乗算され,最終的な単位は
になる.

次の2Dの例では,領域を加熱するために矩形の熱源 が導入されている.熱源の強さは に固定してある.

2D矩形領域と矩形の熱源範囲を定義する:

熱源がもっと複雑な形状の場合は,関数RegionMemberを使って熱源範囲を指定することができる.

熱源領域のRegionMember関数を設定する:

このように簡単な場合は,記号的空間座標RegionMemberFunctionを評価すると,熱源指定で直接使える簡単なブール式および効率的な時間積分を導くことができる.

RegionMemberFunctionを評価する:
強さの熱源 RegionMemberFunctionで設定する:

幾何学的により複雑な場合は,記号的空間座標でRegionMemberFunctionを評価すると未評価のまま戻る.

記号的RegionMemberFunctionが閉形式の解を戻さない場合:

直接RegionMemberFunctionを使うことは可能であるが,時間積分には効率がよくない.ブール式の方が効率的な理由は,未評価のRegionMemberFunctionはコンパイルできないがブール式は自動的にコンパイルできるからである.効率的なPDE係数を設定することに関するこちらのドキュメントも参照のこと.

熱源を扱うのに最も正確で効率的な方法は,要素マーカーを使うことである.こうするとメッシュは熱源に対する特定の部分範囲を持ち,結果的に正確な解になる.メッシュにおけるマーカーとその生成についての詳細は「要素メッシュの生成」チュートリアルに記載してある.

熱源の設定に要素マーカーを使う例は付録の「要素マーカーで熱源をモデル化する」で見ることができる.

体積熱源項 ,初期温度場 の熱移動PDEを定義する:
NDSolveValueでPDEを解く:
温度場プロットの凡例バーとContourPlotオプションを設定する:
PlotListAnimateを使って解のアニメーションを作成する:

アニメーションの質を向上させる方法についてはこちらをご覧いただきたい.

シミュレーションは の不かく乱領域から始まる.体積熱源 が領域内に置かれているので,熱エネルギーが生成され徐々に領域が加熱される.熱移動の速度は材料の熱伝導率 および熱容量 に依存する.

点熱源

点熱源 は,空間的に拡張しないと想定される内部熱源()またはヒートシンク()をモデル化するために使うことができる.点熱源は3Dおよび2Dの軸対称領域で利用することができる.

3Dモデルでは,点熱源はどこに置いてもよい.

2D軸対称モデルでは,点熱源は対称軸上で指定されている場合だけ存在する.対称軸上の点は,軸対称モデルの回転相が考慮に入れられると線熱源を意味する.

点熱源は2Dモデルでは使われない.平面外の線熱源と同等であることを意味するからである.これについての詳細は「層熱源」のセクションを参照のこと.

155.gif

赤で示された点源は茶色で示された2D軸対称領域の対称軸上で指定されている.回転する軸対称領域は点源を線源に変化させないことが分かる.

点源 の単位は常にであり,体積熱源の単位はである.一般に以下が言える:

ここで は単位の総熱量である.

単位は以下である:

これは次のように定式化できる:

ここで はディラックのデルタ関数である. は点源の位置 の各方向で単位を持つ.これで体積熱源 の式が求められる:

しかし,ディラックのデルタ関数は離散化された空間領域で解くことはできないので,数値シミュレーションで問題がある.これはディラックのデルタ関数がソースの場所 において特異だからである.2つ目の問題は有限要素法において,係数の評価はメッシュ要素内で必ず起こり,辺では起こらないという点である.従ってディラックのデルタ関数の近似が必要となる.ディラックのデルタ関数の近似処理は正則化と呼ばれる.

以下のように利用できる正規化デルタ関数 はいくつかある[17][18]:

ここで は正則化デルタ関数 のサポートを制御する正則化パラメータである.通常 はメッシュ間隔 に匹敵するサイズである. は差分 を表す.

正則化デルタ関数 を設定する:
軸対称モデルの点源

軸対称モデルの概念は前と同じであり,体積熱源の式は以下で与えられる:

ここでディラックのデルタ関数はソースの位置 で単位を与える.ここで 平面にあり,0から2 の任意の数である.対称軸の場合, を指定する必要はない.

層熱源

層熱源モデルは,モデルの形状で厚さが持てないほど薄い熱源()または熱シンク()をモデル化する.層熱源は3D,2D,2D軸対称領域で使うことができる.1D領域では層熱源はシミュレーション領域と同じ次元であり,体積源となる.

の単位は常にで指定される.

3D領域では,層熱源は形状で任意に浮いている線を含む辺で指定することができる.

2D軸対称領域では,層原は対称軸を中心とするシミュレーション領域の回転を介して長さとなる点として指定される.

2Dでもモデルの厚さDが平面外の方向にあるため,層源を表す:

192.gif

左:3D層源と平面.右:2Dで見た同じ層源.

3Dモデルの層源

3Dモデルでは,線電荷は形状のどの辺で指定することもできる:

また変形 および も可能である.

ここで はソース位置 におけるディラックのデルタ関数であり,単位を提供する. は層源の値である.

2Dモデルの線源

平面外の層源をモデル化するためには,点をソースとして指定する必要がある.

点には全方向に空間的拡張がないので,ディラックのデルタ関数 をモデリング領域のすべての次元( 等)に適用することができる.

ここで はソース位置 におけるディラックのデルタ関数(単位)であり, は層熱源の値, は単位の厚さである.

次の2Dの例では層熱源 で加えられ領域が熱せられる.層熱源の値はである.

を含むメッシュを生成する:

正則化デルタ関数 を使うために,正則化パラメータがメッシュ間隔の半分()になるようにする.

メッシュ間隔 を調べ,正則化パラメータ を設定する:
正則化デルタ関数を使って点熱源 を設定する:
層熱源 と初期温度場 で熱移動PDEを定義する:
NDSolveValueでPDEを解く:
温度場プロットについて凡例バーとContourPlotオプションを設定する:
PlotListAnimateを使って解のアニメーションを作成する:

アニメーションの品質を向上させる方法についてはこちらをご覧いただきたい.

シミュレーションは である平穏な領域から開始する.層熱源 が領域に置かれると,エネルギーが生成されすべての方向に拡散する.熱移動の速度は材料の熱伝導率 および熱容量 に依存する.

異方性材料および直交異方性材料の熱移動

前のセクションでは,熱せられた媒体は等方性がある,つまり同じ温度勾配 の場合,熱移動率は方向とは無関係であるものと想定した.しかし実際には媒体が異方性である可能性がある.異方性とは,熱が異なる割合で異なる方向に拡散することである.拡散項(19)は次のように書き換えることができる:

または

ここで は熱伝導テンソルである.およびはそれぞれ主伝導係数,非対角伝導係数と呼ばれる.

オンサガー(Onsager)[20]の不可逆過程の熱力学によると,非対角伝導率は相反定理に従わなければならない[21]:

直交異方性材料の熱移動は異方性材料の熱移動の特殊な場合だといえる.ここで材料の熱伝導率は主方向 について対称である.つまり,主方向に沿った値は非零であるが互いに等しくはないということである.非対角伝導係数はゼロである.この挙動はファイバー合成材料等で見られる.熱伝導テンソルは次のようになる:

例として層状の構造を持つ2Dの合成材料を考える:

238.gif

これは熱移動が水平方向により効率的な直交異方性の場合である.熱伝導テンソルは次のように書ける:

直交異方性材料の熱伝導テンソル を設定する:
熱源を設定する:
熱移動のPDEモデルを解く:
温度場プロットの凡例バーとContourPlotオプションを設定する:
直交異方性材料の熱移動を可視化する:

アニメーションの質を向上させる方法はここでご覧いただきたい.

前のセクションで示した例とは異なり,熱移動は水平方向の方が速く,内が高温域となる.

非線形熱移動

前のセクションでは,PDE係数,すなわち密度 ,熱容量 ,熱伝導率 は温度場 とは無関係であると想定した.しかし現実ではこれらのパラメータは,特に純金属の場合,温度で大きく変化する可能性がある.

例として初期の温度場 で温度依存の熱伝導率 の1D熱移動モデルを考える:

方程式(22)は,PDEモデルの熱伝導率 が温度 に依存しているため,非線形熱移動である.

モデル変数を設定する:
非線形伝導係数を設定する:

領域を加熱するために左側の境界に一定の熱流束 を適用する.

NeumannValueで左端に熱流束境界条件を設定する:
温度依存の熱伝導率 ,初期温度場 の非線形熱移動PDEを設定する:
NDSolveValueを使ってPDEを解く:

非線形性の効果を理解するために,線形の熱移動PDEと比較する.

一定の熱伝導率 を持つ線形熱移動PDEを設定する:
線形熱移動PDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

シミュレーションは となる不かく乱領域から始まる.一定の熱流束 を左端に適用すると,熱エネルギーが境界を越えて移動し,領域を加熱する.

非線形モデルの場合,温度 が上昇すると熱伝導率 もそれに従って上昇する.これによりさらに熱移動がスピードアップし温度場がより平坦になる.

温度依存熱容量

次は2つ目の非線形の例である.

一般に,密度 ,熱容量 ,熱伝導率 等の熱伝導材料の特性は,温度 に依存する.温度も従属変数なので,この追加の温度依存性によりPDEモデルは非線形になる.このセクションでは,熱容量 となるような温度依存,つまり非線形のモデルを提示する.

非線形比熱容量の影響を示すために,非線形比熱容量モデルの結果と線形比熱容量モデルの結果を比較する.

初期温度場20[°C],温度依存熱容量 の2D熱移動モデルを考える:

領域を定義する:
モデルの変数とパラメータを設定する:
非線形比熱容量を設定する:

比熱容量 は変数 に依存するため非線形関数である.比熱容量の単位は[],温度の単位は[]である.つまり,比熱容量の単位に合わせるために温度関数に掛け合せる暗的因子[]があるということである.

周囲温度を設定する:
領域の中で底以外の全ての辺に対流境界条件を設定する:
底に断熱境界を設定する:
領域全体は最初温度20[°C]に保たれている:

次に,システムが周囲と熱平衡となるように十分長くシミュレーションを実行する必要がある.システムが熱平衡状態になったときを検知するために,WhenEventを使う.時間積分の間に指定された事象を検知するためにNDSolveに組み込まれた特別なメカニズムであるWhenEventは大変に便利である.事象に関する詳細は「事象のある熱移動」のセクションで見ることができる.

この特別なシミュレーションのために,特定の点における現在の温度と周囲温度との差に対する許容値を選ぶことによって時間積分を中断する事象を指定する.

現在の温度と周囲温度との温度差がより小さくなったときに積分を中止する事象を設定する:

NDSolveにはシミュレーションの終了時間も必要である.平衡に達するのに必要な時間よりも大きい任意の値を選ぶ.

シミュレーションの終了時間を設定する:
パラメータの引数に依存してPDEモデルを生成する関数を設定する:
NDSolveValueを使って非線形熱移動PDEを解く:
積分が止まった時間を調べる:

次に,比熱における非線形性の影響をよりよく理解するために,初期条件として使った一定の周囲温度に比熱を設定してある線形熱移動PDEと比較する.

一定の熱容量 = の線形熱移動PDEを設定する:
NDSolveValueを使って線形熱移動PDEを解く:
積分が止まった時間を調べる:

線形モデルは非線形モデルよりも早く平衡に達する.

PlotListAnimateを使って における解のアニメーションを作成する:

アニメーションの質を向上させる方法についてはこちらをご覧いただきたい.

シミュレーションは = 20[°C]である影響を受けない領域から開始する.対流熱移動が周囲から断熱されている底辺境界以外のすべての辺に起こると,熱エネルギーは境界全体に移動し領域を熱する.領域の温度は,時間の経過とともに周囲と平衡になる.

非線形モデルでは,温度 が上昇するにつれて熱容量 も上昇する.これによって熱移動が遅くなり,線形モデルの温度プロファイルと比較して低いプロファイルとなる.このことから,線形モデルの場合はシステムと周囲との平衡がわずかに早いことが分かる.しかし領域内の温度場の勾配は,時間の経過とともに最終的に弱まり,非線形モデルのシステムでさえ周囲と平衡になる.定常状態では,線形と非線形のモデルシステムの温度プロファイルの違いはない.

事象のある熱移動

熱移動モデリングでよくあるトピックとして,動的熱負荷あるいはパルス熱負荷のシミュレーションがある.つまり,異なる条件下で熱流束がオンになったりオフになったりするのを見ることである.このような場合,WhenEventを使って熱移動モデルを構築し,効率的に解くことができる.

熱が両側から継続的に流出する1Dの屋内モデルを考える.熱は部屋を暖める領域に置かれるが,中心の温度 が閾値温度 を下回ったときにだけオンになる.過熱を避けるために,ヒーターは を超えるとオフになる.

ヒーターの条件付き加熱をモデル化するために,WhenEventを使って体積熱源 を適用する.熱源 は中心温度 または に達したときオン,オフに切り替えられる.

モデル変数を設定する:
与えられた熱源の強さWhenEventを使って条件付き熱源 を定義する:

部屋の冷却をモデル化するために,一定の冷却流束 を領域の両端に適用する.

NeumannValueを使って両端に熱流束境界条件を設定する:

次に,条件付き熱源 を持つ熱移動PDEを定義する.ヒーターは最初オフになっているため,熱源の初期値は で設定する.

条件付き熱源 を持つ熱移動PDEを設定する:
初期室温 のPDEを設定する:

この熱移動PDEを解くために,を離散変数として指定する.つまり,これはNDSolveによる一時的な積分の間,離散時間でのみ変化するということである.

条件付き熱源 を持つ熱移動PDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

シミュレーションは の均一温度から開始される.一定の冷却流束 を両側に適用したため,熱は領域から継続的に流出し,温度を下げる.熱源 は,中心温度 またはに達したときにオン/オフになるため,結果的に振動温度場となる.

動的熱源またはパルス熱源をモデル化するためにIf文を使うことも可能であるが,WhenEventではNDSolveが時間積分中に事象を検出する特別な組込みの機構を持っているということが明らかな長所となっている.この機構はIf文や同様の文でパルス熱源をモデル化する場合には利用できないことがある.

熱パルスのモデル化についての詳細は,付録のセクション熱パルスのモデル化で起こり得る問題と対処法で説明してある.

多孔質媒体における熱移動

多孔質媒体とは固体の骨組み部分と流体で満たされた多孔質領域を持つ多相オブジェクトである.多孔質材料は,その特殊な熱的機械的性質のため,制振,断熱,音吸収等多数の産業分野で広く使われている.

303.gif

多孔質媒体内の熱移動をモデル化するために,直接アプローチという一つのアプローチで,各位相に適した材料係数を持つ2つの熱方程式を含む連立PDEを構築する.一つの方程式は固体領域の温度場 を,もう一つの方程式は流体領域の温度 をモデル化する.

ここで,下付き文字 は固体相のパラメータを, は流体相のパラメータを示している.この2つの方程式は,それぞれの位相の体積分率 で結合されており,2つの位相間の熱交換は右辺にある追加の熱源/シンク項 で明示的に記述されている.

しかし,このアプローチを使うためには領域の完全な細孔構造を再生する必要がある.顕微鏡的な細孔構造は幾何学的に複雑なため,幾何学的な正確さを解決するために細かい有限要素メッシュが必要となる可能性があり,そうなるとモデルを解くのに莫大な計算コストがかかることが想定される.

別のアプローチとして,孔を顕微鏡スケールでモデル化するというものがある.この場合,熱移動モデルは平均的な温度場を使って,両方の位相の細孔構造全体を表す.このため体積平均の有効な熱特性が使われ,モデルは単独の熱方程式で表すことができる.

ここでは有効体積熱容量, は有効熱伝導率である.これらは体積分率 と各位相の特性に基づいて計算される.

多孔質媒体が完全に飽和したものと仮定すると,両位相の体積分率は で関連付けることができる.したがって方程式(23)は次のようになる.

流体の体積分率 は媒体の多孔質としても知られる.簡単にするために,次のセクションではその下付き文字 を省略する.

例として,多孔質媒体内の2D熱移動モデルを考える.流束は幅,高さで領域を通過しており,領域を加熱するために一定の熱流束が左面に適用される.熱移動における多孔質の効果を理解するために,3つの異なる多孔値 (純粋な流体媒体)を使ってモデルを解く.

固体相および流体相の密度 ,熱容量 ,熱伝導率 を設定する:
多孔質媒体の有効熱伝導率 と体積熱容量

方程式(24)は時間微分の前に有効体積熱容量を持つ.同時にこれは液相に対する材料パラメータを使う.それを満たすために,この方程式は部分的にHeatTransferPDEComponentおよび追加の時間微分項で生成される.そのため,モデル変数には時間依存変数が含まれるが,時間変数自体は含まれない.有効体積熱容量を伴う時間微分は,手動で加えられる.

モデル変数を設定する:
モデルパラメータと右への流出 を設定する:
多孔質媒体の熱移動PDEを初期温度場 で構築する:

領域を加熱するために,において左面に領域内への一定の熱流束 を適用する.

NeumannValueを使って左端に熱流束境界条件を設定する:
初期温度場 ,右の流出 の多孔質媒体の熱移動PDEを構築する:
ParametricNDSolveValueを使って,多孔率 のPDEを繰り返し解く:
多孔率 の多孔質領域の温度場を可視化する:

この温度場は 方向について対称である.

熱移動に対する多孔質の効果を調べるために, 軸に沿った温度変化を3つの異なる多孔値 (純粋な流体媒体)で比較する.

における の温度分布を調べる:

固体相の熱容量値は流体相のものより小さい()ので,固体の量が多い(多孔値 が低い)媒体はより温度変化を受けやすい.このことは,における有効体積熱容量を計算することで証明することができる.

異なる多孔値 で有効体積熱容量を比較する:

のとき,有効熱容量は最小になる.

複合次元の熱移動モデル

このセクションでは,さまざまな空間次元で定義された熱移動現象をモデル化する方法を説明する.例として,一様の初期温度を持つ2Dのセラミック板を考える.セラミックは左面に冷却流束が流れる薄いパッドを付けて冷却する.冷却流束 はパッドと熱移動係数 のセラミックとの間の温度差に比例する.

セラミックの上部には一定の熱流束 が適用される.セラミックの下面と右面は断熱されるものと仮定する.

364.gif

セラミックの温度場 について解くために最初の思いつくのは,単独の2Dシステムを構築し,冷却パッドを熱流束境界条件としてモデル化することである.しかしこの例では,冷却流束 だけでなくパッドの温度 にも依存している.パッドの温度 の値を決定するためには,薄い冷却パッド内の熱移動もモデル化する必要がある.

冷却パッドがセラミックの幅よりもずっと薄いと仮定すると, 方向におけるパッドの温度変化は無視できる.したがって,冷却パッドは1D領域としてモデル化し,セラミック板は2D領域としてモデル化することができる.つまり複合次元のモデルを使うのである.

簡単にするために,セラミック板とパッド両方の熱移動係数 と熱伝導率 ,密度 と熱容量 は1に設定する.

セラミック版の温度場 は2Dの熱方程式で記述できる.

セラミック板の温度を記述するために2Dの熱移動モデルを設定する:

パッドは 軸のに沿った1D領域として記述する.これはセラミック板の左面に一致する.1Dのパッド領域に沿ったパッドの温度 をモデル化するために,1Dの熱方程式を使うことができる.

しかし連立のPDE系を解くためにはNDSolveはすべての従属変数が同じ空間次元を持つことを必要とする.このためだけに 方向にパッド温度 のための「仮想次元」を使わなければならない.

つまり,パッド温度 を1Dのパッド領域だけでなく,セラミック板の2D領域全体で解くのである.しかし, はパッドのための仮想の次元であるため,セラミック領域内では 値に何の物理的意味もない.結果のパッドの温度場 上の1Dパッド領域でのみ有効なのである.

方向に仮想次元を持つパッドの熱移動モデルを設定する:

次に,冷却パッドとセラミック板の間の熱交換を考える必要がある.セラミックの見地からすると,熱は左のセラミック境界からパッドに失われていく.これは熱流束境界条件でモデル化することができる.

セラミックの左面に,冷却熱流束 熱流束境界条件を設定する:
熱移動係数を設定する:
熱移動値で同じ境界条件を設定する:

冷却パッドの見地からすると,パッド領域全体でセラミック板から熱を得る.これは方程式(25)の熱源項 でモデル化することができる.

エネルギー均衡の法則により,パッドの熱源 は同じ大きさであるが,セラミックの冷却流束 とは反対の符号を持つ.

熱源 を設定して,パッドの熱利得をモデル化する:

セラミックの上面に,一定の熱流束 を適用する.

セラミックの上境界に一定の熱流束を設定する:
セラミック板の初期温度 と冷却パッドの初期温度 を設定する:
連立の熱移動PDEモデルを解く:
セラミック板と冷却パッドの温度場を可視化する:

アニメーションの質を向上させるための注意事項をご覧いただきたい.

次に1Dのパッド領域ないのパッド温度 を調べ,セラミックの左面の温度 と比較する.よりよい可視化のために,プロット内の温度範囲を再スケールするカスタム関数TwoAxisPlotを定義して適用する.

に沿った におけるパッド温度 とセラミック温度 を調べる.

パッドの熱利得はセラミック温度 に依存するため,パッド温度 はセラミックの左面 と類似したパターンに従う.

横断面 のセラミックの温度場を調べる:

連立PDE系に仮想次元を導入することによって,1Dと2Dの熱方程式を含む複合次元モデルが解けるようになった.この方法は1D,2D,3Dにも同様に適用することができる.

マルチマテリアル媒体における熱移動

マルチマテリアル領域における熱移動は多数の産業アプリケーションにとって関心のある基本的問題である.考慮する形状の空間座標に沿った複合材料の材料特性を変化させると,等方性の単一材料媒体で見られるよりも複雑な熱移動メカニズムになる.以下の例は適切な境界条件を持つ複合材料内の熱移動の設定および解法を示す.

初期温度場が297[°C]である1Dの過渡熱移動マルチマテリアルモデルを考える:

1D複合材料は厚さ0.25,0.114,0.04[]の3つの異なる材料からなり,それぞれ熱伝導率が8,1.8,44[],密度が3100,2100,7800[],熱容量1050,1100,540[]等,異なる特性を持つ.このマルチマテリアル媒体の左端は1700の高温であり,対流と放射の効果によって右端で熱はなくなっていく.マルチマテリアル領域の初期温度は297に保たれる.この問題は,指定された境界条件,初期条件を持つ支配熱移動モデルの領域依存特性値を指定することで解く.

領域を定義する:
境界メッシュを設定する:

指定された領域内のマルチマテリアルの特性を指定するために,領域マーカーを使う.これらはToElementMeshの一部として指定される.メッシュ内で領域マーカーを指定することについての詳細は要素メッシュの生成モノグラフのマーカーのセクションに記載されている.

名前付きの領域マーカーを設定する:
メッシュを設定する:
モデルの材料パラメータを層で定義する:
ヘルパー関数を定義してElementMarkerを使って各層の材料パラメータを指定する:
すべての形状に対する熱移動モデルパラメータを定義する:
モデルのモデル変数varsとパラメータparsを設定する:

次のステップとして,マルチマテリアル領域の左境界でディリクレ境界条件を設定し,右境界で対流境界条件と放射境界条件を設定する.左境界は の高温で維持され,熱は右境界で周囲への対流・放射効果によってなくなっていく.

左境界で温度表面境界条件を設定する:
右境界で対流境界条件を設定する:
上面で熱放射境界条件を設定する:
領域の初期温度を = 0[ºC] = 297 [K]に設定する:

3つの層の間の熱伝導を解析するためには, から までのモデルを解く.

積分の初期時間と終了時間を定義する:
熱移動方程式を定義し,HeatTransferPDEComponentを使って解く:
熱移動PDEモデルを解き,時間/メモリ使用量を監視する:
さまざまな時間における温度場のプロットを作成する:

システムの温度は最初ずっとで維持される.左端の温度が突然に変更されると,右端の熱はその境界の時間依存温度に依存する対流・放射流束によってなくなる.マルチマテリアル領域の各部分の温度は,隣接した領域の材料特性に依存する結合における離散勾配で低くなる.時間の経過につれて,各部分の温度プロファイルは一定の勾配を持つ定常状態を実現する.

非接触式風速計または多層球の熱伝導の応用モデルで,複数材料の熱移動がさらに見られる.

相変化を伴う熱移動

熱力学における相変化とは物質がある状態(固体,液体,気体,プラズマ)から別の状態に変化する現象のことである.相変化は一定の温度と圧力で,十分なエネルギーがシステムに加わるかシステムから除去されるかするときしか起こらない.相変化に関連するエネルギーは潜熱 といわれ,物質の温度変化を引き起こすのではなく,分子構造を変化させる.相変化は相転移と呼ばれることもある.数学用語のステファン問題も有名である.

氷から水への相変化

例として,の氷の棒の周りにおける氷から水への相変化を示す以下の1Dモデルを考える.この棒は初期温度 であり,棒が溶けるように左端で一定熱流束 を適用する.右端は断熱されているものとする.

熱移動モデルでは,厳密に相変化温度 で相転移をシミュレーションするのではなく,温度区間 から において転移が起こると想定する.転移の間の物質相は,物質内部におけるもとの相の新しい相に対する割合を表す平滑化された階段関数 で説明できる.

物質相を表す平滑化された階段関数 を定義する:

温度区間 内の等価密度 および伝導率 は以下で与えられる:

氷と水の密度 および熱伝導率 を指定する:

より一般的に言うと,平坦化された階段関数は,等価密度 および熱伝導率 の引数でもある.

等価密度 および熱伝導率 を定義する:

しかし,等価比熱容量 は相転移に必要な潜熱 を説明するために余分の項 を含まなければならない.ここで は相変化の間の潜熱 の分布を表す.これは相変化の温度 付近の正則化デルタ関数で近似される:

融解潜熱 と潜熱の分布 を指定する:

の積分は,相変化に必要な潜熱 に等しい:

と潜熱 の等価性を調べる:

等価熱容量 は以下で与えられる:

等価熱容量 を定義する:

は3つのパラメータ と潜熱の分布に依存する. は上で指定した関数である.後の実験で を別の関数に置き換える.

棒の初期温度 を設定する:
左端で熱流束境界条件を一定熱流束 で設定する:

棒の右端ではデフォルトの断熱境界条件が適用される.

熱移動PDEモデルを解く:

相変化における潜熱 の効果をよりよく理解するために,上の結果と潜熱を無視した解を比較する.この場合は,すべての入力に対して0を返す関数に を設定する.

上で述べた通り,潜熱の影響を無視するために常に0となる の関数を設定している.そのために,どの入力に対しても0を返す関数を作成する.

どのような入力に対しても0を返す関数:
関数が0を返すことを確認する:
潜熱 を無視する熱移動PDEを設定して解く:
PlotListAnimateを使って解のアニメーションを作成する:

温度間隔がシミュレーションにどのような影響を与えるかという質問に答えるために,別の実験を行う.ここでは温度区間を使う.

温度区間に対する平滑化された階段関数を設定し可視化する:
新しい正則化関数を設定する:
PDEを解く:
およびでの解の絶対差をプロットする:

を使ったときの解の絶対温度差は約であることが分かる.全体の温度と比較すると,これは小さい.

管の中の液体の凍結

この2つ目の例では管の中の水の凍結をモデル化する.管の中の水は最初温度はである.管の右端はずっと同じ温度を保つ.管の左端は一定時間の間に冷やされる.その後温度は再びになる.

基本的な考え方,関数,材料パラメータは上と同じである.

領域を定義する:
メッシュを設定する:

次のステップとして,領域の左右の境界でディリクレ境界条件を設定する.先に述べたように,右境界は常に温度は5[°C]であるが左境界は温度が変化する.

の時間では,5[°C]から-5[°C]に変更するために平滑化されたステップ関数が使われる.その後の時間では温度は に保たれる.の時間では左端が に熱せられる.シミュレーションの残りの時間はずっとこの温度のままになる.

ある温度から別の温度への遷移のための平滑化関数の使用は,その間の物理的遷移を持つことである.中間値なしで1つの温度から別の温度に変化することは反物理的な動作であり,避けるべきである.

カスタム定義の関数を使って温度を変化させる時間を設定する:
温度を変化させる時間の関数を使って左側にディリクレ条件を,右側に一定条件を設定する:
領域全体の温度は最初5[°C]である:
熱移動PDEモデルを解く:
ContourPlotを使って,相境界のある温度場のプロットを作成する:
PlotとListAnimateを使って,両方の場合の相境界を持つ温度場のアニメーションを作成する:

システムの温度は最初に保たれている.左端の温度がすぐにに切り替わると,左端に氷ができ始める.この氷の層は左端の温度がに切り替わるまで,システムの中を右に向かって広がっていく.温度がに切り替わると,システム内の氷の相がすべてなくなり水の相に戻るまで,両端から氷が溶けていくのが分かる.システム全体が水の相になると,温度上昇には顕熱しか必要ないのでシステムはより早く熱せられる.最終的にシステムは定常状態温度に達する.

モデル次数低減を伴う熱移動

異なる初期データを使って同じ熱移動を再実行したい場合がある.この場合モデル次数低減が役立つことがある.モデル次数低減とはPDEの離散化を使って,より効率的に時間積分することのできる同様の離散化を求めるというものである.モデル次数低減を実行するためには,NDSolveが作成する離散化へのアクセスが必要である.これは現在NDSolveのレベルでは不可能であり,ある程度のプログラミングを要する.これについては,有限要素プログラミングチュートリアルのセクション定常係数と定常境界条件を持つ非定常偏微分方程式のモデル次数の簡約で説明する.

マルチフィジックス熱移動

熱移動は物理の別の分野と組み合されることがよくある.以下に挙げたのは,熱移動を使うマルチフィジックスの応用例である:

熱移動の境界条件

熱移動モデリングで最も一般的な境界条件はDirichletConditionNeumannValuePeriodicBoundaryConditionを使ってモデル化することができ,以下の4つのタイプに分類できる:

この4つのタイプに,以下の境界条件を分類する:

次のセクションでは,熱移動に共通するいくつかの物理的境界と,それらをDirichletConditionNeumannValuePeriodicBoundaryConditionを使ってモデル化する方法を説明する.このため,現在扱っている境界条件は常にシミュレーション領域の左側にあるものとする.例によっては左側の境界条件の挙動をよりよく示すために,右側に追加の境界条件を設定することもある.

表面温度境界条件

目的

表面温度境界条件は境界のある部分に特定の温度を設定することを目的とする.

定式化

境界上に温度 を指定すると,温度面条件は以下で与えられる:

DirichletConditionでモデル化された従属変数 の温度面境界:

導出

境界上の表面温度 を既定することは,温度境界条件と呼ばれる.表面温度 は一定でも時間依存値でもよく,熱移動PDEモデルではDirichletConditionを使って設定される.

たとえば,熱エネルギーを領域に送る加熱壁をモデル化するためには,過渡的表面温度 を左端に設定することができる.ノイマンのゼロ条件は断熱境界として右端で適用される.

ここで滑らかな階段関数を使って,表面温度 から までのプロファイルを記述する.過熱過程のシミュレーションのために,パラメータ および は任意に選ばれる.

表面温度プロファイル を滑らかな階段関数で指定する:
変数とパラメータを設定する:
DirichletConditionで温度境界条件を設定する:
不かく乱初期温度場 でPDEを設定する:
NDSolveValueでPDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

シミュレーションは である不かく乱領域から始まる.表面温度 が左端で上昇するに従って過剰な熱エネルギーが右側に移動し,領域全体の温度を上昇させる.熱移動の速度は材料の熱伝導率 と熱容量 に依存する.

熱流束境界条件

目的

熱流束境界条件は境界のある部分に流入,あるいはそこから流出する単位面積当たりの熱エネルギー率をモデル化することを目的とする.

定式化

境界 上に熱流束 が指定されると,熱流束境界条件は以下で与えられる:

NeumannValueでモデル化された従属変数 の熱流束境界:

導出

境界に垂直な熱流束 が指定され,それがゼロではない境界は熱流束境界と呼ばれる:

慣例により の前に負符号を付けて熱流束が外向きの法線 とは反対向きに指定されていることを示す.したがって, の正値は熱エネルギーが領域に入る内向きの熱流束を,負の は外向きの流束を意味する.

フーリエの熱伝導法則(26)は熱流束 と温度勾配 を関連付けるものである:

(27)を(28)に挿入すると,熱流束境界条件は次のように書ける:

熱流束 の単位は境界の次元に依存する.1D (),2D (),3D ()の領域では, の単位はそれぞれとなる.

以下の例では のとき領域を加熱するために過渡的熱流束 が左境界に適用される.

変数とパラメータを設定する:

熱流束のプロファイルは次のように定義される:

NeumannValueを使って左端に熱流束境界を設定する:
不かく乱初期温度場 でPDEを設定する:
NDSolveValueでPDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

熱流束 が左境界に適用されると,熱エネルギーはこの境界を越えて徐々に領域を加熱する.熱流束は時間 でオフになる.この時点で境界条件の特性は熱流束から断熱へと変化する.それからでこぼこの温度場が時間とともに内部熱移動により平滑化される.

アニメーションを見るとt=300において解のジャンプがあるように見える.計算されたフレーム数の影響である.

時間間隔における解の差を示す:

熱流束の値 はフーリエの法則によって温度勾配 に関連している.つまり,熱流束 は境界に垂直な温度勾配を直接制御するのである.

断熱境界条件

目的

断熱境界条件は,熱流束がない境界をモデル化することを目的とする.

定式化

断熱境界条件は以下で与えられる:

NeumannValueでモデル化された独立変数 の断熱境界:

境界のある部分で境界条件が設定されていない場合,暗示的にノイマンのゼロ境界条件が使われる.

導出

断熱境界条件は熱流束のない境界を表す:

(29)を熱流束境界条件(30)に挿入すると,断熱境界条件は次のように書くことができる:

次の例では左端の境界に断熱境界を置き,右端に熱源として一定の熱流束 を置く.

変数とパラメータを設定する:
断熱境界条件をモデル化するためには,NeumannValueを左端でに設定する:

境界のどの部分にも境界条件が指定されていない場合,デフォルトでノイマンのゼロ境界条件が使われる.これは,与えられた境界で境界条件が指定されていない場合は断熱境界がデフォルトの境界条件であるということである.

右側の境界に一定の熱流束を設定する:
不かく乱初期温度場 でPDEを設定する:
NDSolveValueでPDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

一定の熱流束を右境界に置くと,領域内の温度が徐々に上昇する.しかし左端の断熱境界では温度勾配 はずっとゼロのままである.

温度勾配 はフーリエの法則 により熱流束 と関連しているので,温度勾配がゼロであるということは境界(断熱境界)に熱流束がないことを意味する.

対称境界条件

目的

対称境界条件は,計算領域および想定される温度場がシミュレーション領域の軸について鏡面対称であるときに使われる.

定式化

対称境界条件は以下で与えられる:

NeumannValueでモデル化した従属変数 の対称境界:

境界のある部分で境界条件が設定されていない場合,ノイマンのゼロ境界条件が使われる.

導出

対称境界条件は,計算領域の範囲を完全な物理的モデル形状の対称部分領域に縮小するために使われる.これにより少ないメモリで高速の解法プロセスが可能になる.

561.gif

1Dシステムの温度場を から までで解く場合を考える.温度パターンが に対して鏡面対称であると想定されるとき,実質的にシステムの左半分だけでシミュレーション領域を構築することができる.次に で対称境界条件を適用する.

対称であるため,対称境界の温度勾配はずっとゼロのままである.これは境界の熱流束がゼロであることを意味する.したがって対称境界条件は断熱境界条件と同じである.

流出境界条件

目的

熱移動が流速 の流体媒体内で起る場合,熱が流体流れによって領域から出される出口をモデル化するために流出境界条件が使われる.

定式化

流体媒体内での熱移動をモデル化する場合,出口 における流出境界条件は以下で与えられる:

NeumannValueでモデル化された従属変数 の流出境界:

境界のある部分で境界条件が設定されていない場合,ノイマンのゼロ境界条件が使われる.

導出

熱移動を流体流れでモデル化する場合,流出口境界の拡散熱流束 はゼロに設定される.これは領域外の流れの温度場がモデル化領域の内部の流れには影響を及ぼさないことを意味する.

流出境界条件は十分に発達した流れにのみ適用される.つまり,流出口では速度プロファイル は流れの方向について変化しないのである.

乱流でよくあることであるが,流出境界で再循環がある場合,再流入する流れは領域内の流れの温度場に影響を及ぼし,拡散流束がないという仮定が破られる.このような場合は,流出境界条件はもう使えない.

流出境界条件は実質的にノイマンのゼロ条件であるため,与えられた境界で境界条件が指定されていない場合,流出境界条件が適用される.

対流境界条件

目的

対流境界条件は,境界の同じ部分に隣接する流れによって引き起こされた,境界内の熱エネルギー移動をモデル化することが目的である.

定式化

境界 上に外部温度のプロファイル と熱移動係数 があるとき,対流境界条件は以下で与えられる:

NeumannValueでモデル化された従属変数 の対流境界:

導出

境界表面に隣接して領域の外側に流れる流体がある場合,熱エネルギーの一部は流体粒子の動きを介して境界中に移動する.これは対流過熱あるいは対流冷却として知られる.

1701年,ニュートンは2つの媒体間の対流熱移動比がその温度差に比例することを発見した.したがって,対流熱流束は次のように定義される:

ここで は外部液体の温度,単位 は対流熱移動係数を表す.熱移動係数は実験的に決定され,密度 ,熱拡散係数 ,および速度 やレイリー数 等の外部流体の流れ状態等の材料特性に依存する.

対流熱移動係数 のおおよその範囲を下に挙げる:

さまざまな場合の熱移動係数 を推定するために,実験に基づく定式もいくつか構築されている[31,32].

(33)を熱流束境界条件(34)に挿入すると,対流境界条件は次のように書くことができる:

次の例では,初期温度場 の領域を加熱するために左領域に一定の外部流 を適用する.領域全体の対流熱移動は,指定された熱移動係数 によってモデル化される.

NeumannValueで対流境界条件を設定する:
不かく乱初期温度場 でPDEを設定する:
NDSolveValueでPDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

外部流が左境界から領域の外側に流出すると,過剰な熱エネルギーが境界を越えて領域内に移動してくる.対流熱流束は境界の温度差に比例するので,境界の温度 が外部温度 に近付くにつれて左境界の温度勾配 は徐々に小さくなる.

熱放射境界条件

目的

熱放射境界条件は,境界のある部分における放射による加熱あるいは冷却をモデル化することを目的とする.

定式化

境界 で周囲温度 ,表面放射率 ,シュテファン・ボルツマン定数 のとき,熱放射境界条件は次のように書くことができる:

NeumannValueでモデル化された従属変数 の熱放射境界条件:

導出

絶対零度()より高い温度を持つ物体はすべて,電磁放射を介して常に熱エネルギーを放出する.放射量は物体の温度と表面状態の両方に依存する.完全な熱放射体である黒体の場合,シュテファン・ボルツマンの法則によると放出される熱流束は物体の絶対温度の四乗に比例する:

ここで はシュテファン・ボルツマン係数である.

しかし,実際に放出される熱流束は黒体放射の熱流束より「表面放射率」という割合 だけ小さい.放射率の値はであり,放射体の物理的特性および表面条件等の要因に依存する.

シュテファン・ボルツマンの法則は次のように書き換えることができる:

(35)に基づいて,境界を超える正味の放射熱流束 を調べることにより,放射境界条件を定式化することができる:

境界から放射される外向きの放射流束 はその表面温度 と放射率 に依存する:

同時に境界は環境から入り込む放射を吸収する.この吸収される内向きの放射流束 は以下で与えられる:

ここで は周囲温度であり は環境の放射率である.追加の項 (表面吸収因子)は, 右側で増大され,境界の吸収率の主要因となる.この項は到着する総流束に対して実際に吸収された放射流束の割合を表す.

したがって,正味の境界を越える正味の放射熱流束は次で与えられる:

熱力学的平衡を満足するために,任意の物体について吸収率 はその放射率 と等しくなければならない.これはキルヒホッフの放射法則として知られる[36].

方程式(37)は次に簡約される:

周辺環境が放射率 の黒体として動作すると仮定すると,方程式はさらに次のように簡約できる:

(38)を熱流束境界条件(39)に挿入すると,放射境界条件は次で与えられる:

または

上の導出方法は絶対温度に基づいて実行された.つまり(40)の温度項の単位はケルビンである,

放射境界条件を摂氏で適用するためには(41)で単位変換をしなければならない:

ここで は絶対零度の温度を表す.

例として,左境界の周囲温度 と表面放射率 を考えてみる.境界を越える正味の放射熱流束は放射境界条件でモデル化される.

変数とパラメータを定義する:
摂氏の単位で熱放射境界条件を設定する:
不かく乱初期温度場 でPDEを設定する:
NDSolveValueでPDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

左側の方が周囲温度 が低いため,正味の放射熱流束は左側の境界からシステムから流出する.これにより領域は徐々に冷えていく.

対流境界条件と同様に,左境界の温度勾配 は境界の温度差に依存し,フーリエの法則 で計算することができる.

における左境界の温度勾配を調べる:

左端の温度勾配はとなる.

周期境界条件

目的

周期境界条件は境界の一部の温度を別の部分にマップして,領域の周期性をモデル化することを目的とする.

定式化

周期境界 の温度 をターゲットの境界にマップする関数が与えられた場合,周期境界条件は以下のように書くことができる:

PeriodicBoundaryConditionでモデル化された従属変数 の周期境界条件:

PeriodicBoundaryConditionを指定する方法についての詳細は,リファレンスページを参照されたい.

導出

周期境界条件は空間的に周期的な領域で熱移動を計算するために適用される.ターゲットとなる境界が与えられると,周期境界の温度 は,指定された関数 を使ってターゲット表面の温度 にマップすることができる.境界条件は,熱移動PDEモデルのPeriodicBoundaryConditionで設定される.

例として,熱源として薄い膜が挿入されたリングヒーターを考える.周期境界条件を使って,1D領域でシミュレーションを実行することが可能である.

653.gif

リングヒーターはリングの周長である長さ で1Dモデルに変換される.加熱過程のシミュレーションを実行するために,加熱フィルムの温度 を左境界で指定する.右境界では,境界温度 をターゲット境界にマップするために周期境界条件が適用される.

マッピング関数 を使って,右端でPeriodicBoundaryConditionを設定する:

ここで加熱フィルム上の温度プロファイル を指定するために滑らかな階段関数を使う.パラメータ は熱移動モデルのパラメータに適合するよう任意に選ばれる.

変数とパラメータを設定する:
左端で加熱フィルムの温度プロファイル を指定する:
不かく乱初期温度場 でPDEを設定する:
NDSolveValueでPDEを解く:
PlotListAnimateを使って解のアニメーションを作成する:

シミュレーションは となる不かく乱領域から始まる.加熱過程では,加熱フィルムの温度 が上昇する間周期境界上の温度 が左境界にマップされる.熱エネルギーは両端から移動し,徐々に領域を加熱する.

この他,この例ではPeriodicBoundaryConditionの代りに,右境界でも左境界のDirichletConditionを使うことができる.

周期境界条件の代りの設定:

左側が他の境界条件ならばPeriodicBoundaryConditionがそれを右に投影する.つまりPeriodicBoundaryConditionは,マップ関数 が指すところで見付けたどのような境界条件でも述語として指定されるターゲットに投射する.

付録

熱方程式の特殊な場合

定常の場合

温度場が定常状態である場合,(42)の過渡項は消失し熱方程式は下のように簡約される:

円筒座標の熱方程式

熱移動問題をモデル化する場合,モデルを直交座標で説明するのは便利ではない場合もある.そのような場合,熱方程式は円筒座標系や球面座標系を使って表すこともできる.

円筒座標を示すグラフィックス:

円筒座標系では はそれぞれ射線,方位角,垂直の方向を示す.円筒座標は直交座標ので次のように定義される:

あるいは

座標関係(43)を(44)に挿入することで,熱方程式は円筒座標で次のように表すことができる:

モデル内部の熱移動が 軸について回転対称である場合,結果の温度場 方向では不変である.よって方程式(45)は次のように簡約できる:

この場合,この対称特性を使うと,3Dの熱移動問題を2D領域でモデル化することができる.このタイプのモデルは線対称モデルといわれる.2Dの軸対称モデルで円筒座標系を利用する例は,熱移動検証テストノートブックの時間非依存の2Dの例および時間依存の2Dの例で見ることができる.

HeatTransferPDEComponentは,パラメータRegionSymmetryを指定し,それをRegionSymmetryに設定することで,熱移動方程式(46)の軸対称形式を生成することができる.

2Dの軸対称モデルで円筒座標系を利用する例は,熱移動検証テストノートブックの時間非依存の2Dの例および時間依存の2Dの例をご覧いただきたい.HeatTransferPDEComponentの参照ページにも例が記載されている.

管型反応器における放射効果の応用モデルでも軸対称モデルが使われている.

球面座標の熱方程式

球面座標を示すグラフィックス:

球面座標系では はそれぞれ射線,方位角,極の方向を示す.球面座標は直交座標で次のように定義することができる:

または

座標関係(47)を(48)に挿入することにより,熱方程式は球面座標系で表すことができる:

拡散方程式の平滑化特性

拡散方程式の基本的な動作は平滑化である.その効果を見るために,で不連続な初期条件を設定する.

不連続な初期条件を初期化してプロットする:
温度場のPDEは次で与えられる:

すべての境界の境界条件を指定するために,述語にTrueが使われる.存在する場合は,内部境界を含む可能性がある.

NDSolveValueでPDEを解く:
初期条件を可視化し,解を平滑化する:

短時間でx=1/2における不連続な初期条件の尖りが滑らかになった.この平滑化効果が拡散方程式の主な特徴であるため,熱方程式の特徴でもある.

要素マーカーを使った熱源のモデル化

熱源を扱う最も正確で効率的な方法は,要素マーカーを使うことである.メッシュが熱源についての特定の部分範囲を持ち結果として正確な解になる.このようになる大きな理由は,個々のメッシュ要素が材料境界を横断しない場合,有限要素法の方がPDEを解くのにより適しているからである.要素マーカーとメッシュでの使用法は要素メッシュの生成チュートリアルの要素マーカーに記載してある.

次の2Dの例では,矩形の熱源 を使って領域を加熱する.

熱源を残りの領域と区別する内部境界を持つ境界要素メッシュを作成し可視化する:

ToElementMesh"RegionMarker"オプションで範囲のマーカーを指定する.あるいは熱源範囲を調整するために別の最大セル測定値を指定することもできる.

内部マーカーで要素メッシュを作成し可視化する:

メッシュを生成する代りに,ブール領域関数を使って,領域の穴を挿入しないよう指定する.上と同じ範囲マーカーの過程を使う.

範囲を指定する:
メッシュを設定する:
変数とパラメータを設定する:
熱源の強さ を設定する:
熱源 を要素マーカーで設定する:
熱源項 と初期温度場 で熱移動PDEを定義する:
NDSolveValueでPDEを解く:
温度場プロットの凡例バーとContourPlotオプションを設定する:
PlotListAnimateを使って解のアニメーションを作成する:

アニメーションの質を向上させる方法はこちらでご覧いただきたい.

この結果は,「立体熱源」で提示した要素マーカーを使わなかった場合の結果と一致する.複雑な立体では,領域要素マーカーを使うと設定が簡単で,より効率的に解が計算できる.

不連続のPDE係数を持つ保存則

熱移動PDEが不連続な密度/速度場についても当てはまる理由を説明するために,質量保存から導出を始める:

密度 と流速 は接触面 で不連続であるとする:

接触面 では質量の生成も破壊もないので,質量流束は の両側で等しい:

つまり,質量流束 は領域全体で連続でなければならない.ここで方程式(49)に発散定理[50]を適用すると以下が得られる:

領域は完全に任意なので,積分を無視し質量保存方程式を微分形式で作成することができる:

すなわち,質量保存方程式(51)は不連続の密度/速度場を持つ領域でも真となる.

熱方程式は実質的にエネルギー保存方程式であり,(52)の密度 を内部エネルギー で代替することにより同じような方法で導くことができる.したがって,このチュートリアルで示した熱移動モデルは連続と不連続のどちらの密度/速度場にでも適用することができる.

熱パルスのモデル化で起こり得る問題と対処法

セクション事象のある熱移動では,時間依存の熱移動モデルで熱パルスをモデル化するためにIf文を使うことができると述べた.しかし,パルスの継続時間が非常に小さくなると,NDSolveのデフォルトの時間刻みアルゴリズムでは既定の熱パルスを検出することができない場合がある.

以下のセクションでは例を使って問題点を示し,熱パルスのモデル化での2種類の対処法を紹介する.

領域の中央で適用された周期的な熱パルスを持つ1D時間依存熱移動モデルを考える.時間区間 には継続時間の合計5つのパルスがある.

If文を使って周期的な熱パルス を設定する:
1D領域を定義する:

領域の冷却をモデル化するために,両端の温度はに固定する.

領域の両端で温度面境界条件を設定する:
領域内で熱パルスを活性化する部分を定義する:

デモンストレーションの目的で,熱伝導率 ,密度 ,熱容量 を1に設定する.

周期的な熱パルス を持つ時間依存の熱移動PDEを解く:
温度場を調べる:

オプションを使わなかったので,NDSolve における最後の2つの熱パルスを見逃した.

次で,この問題の対処法を2つ紹介する.

方法1 - MaxStepFractionを減らす

最初の方法は,NDSolveによる時間積分過程により小さい時間刻みを指定するというものである.オプション"MaxStepFraction"を使うと,合計の時間範囲の最大割合を1ステップでカバーするように設定することができる.

PDEを解き,NDSolveが時間積分に最低100ステップを使うことを確実にする:
温度場を調べる:

時間刻みが小さいと,NDSolveの間の5つの熱パルスをすべて検出する.

方法2 - WhenEventを使う

上の方法よりも優れた別の方法に,WhenEventを使って熱パルスを指定するというものがある.WhenEventではNDSolveが時間積分の間にパルスを検出するというメカニズムを備えているという明白な利点が使える.このメカニズムはIfまたは類似の文で熱パルスをモデル化する場合には利用できない可能性がある.

この場合はどちらかの式がゼロになったときにトリガされる事象を使う.トリガされると離散変数 が切り替わる.

WhenEventの詳しい使い方およびその事象検出方法はこちらを参照されたい.

WhenEventを使って周期的な熱パルス を設定する:
周期的な熱パルス を持つ時間依存の熱移動PDEを解く:

この方法では,熱パルスはNDSolveによって離散変数 として扱われる.つまりこれは時間積分中の離散時間でのみ変化するということである.解を得るためには,時間積分の刻み幅の減少は必要ではない.

温度場を調べる:

WhenEventを使うことで,NDSolveの間の5つの熱パルスをすべて検出した.

用語集

参考文献

1.  Bilbao, Stefan and Hamilton, Brian. Directional Source Modeling In Wave-Based Room Acoustics Simulation. IEEE, 2017.

2.  Peskin, Charles. The Immersed Boundary Method. Cambridge University, 2002.

3.  Churchill, Stuart W. and Chu, Humbert H.S. Correlating equations for laminar and turbulent free convection from a vertical plate. International Journal of Heat and Mass Transfer, 18 (11): 13231329, 1975.

4.  Sukhatme, S.P. A Textbook on Heat Transfer (Fourth ed.). Universities Press. pp. 257258, 2005.

5.  Riedl, M. Optical Design Fundamentals for Infrared Systems (Second ed.). SPIE Press, Bellingham, WA, 2001.

6.  Holman, J. P. Heat Transfer Tenth Edition, McGraw-Hill. pp. 111, Example 3-10 (2008).

7.  Weisstein, Eric W. Divergence Theorem, MathWorld-A Wolfram Web Resource. http://mathworld.wolfram.com/DivergenceTheorem.html.

8.  Onsager, L. Physical Review. 37, 405426, (1931); 38, 22652279, (1931).

9.  Casimir, H. B. G. Reviews of Modern Physics. 17, 343350, (1945).

10.  Hahn, D. W. and Özişik, M. N. Heat Conduction, John Wiley & Sons, Inc, ch.15 (2012).