多層球の熱伝導

はじめにPDE熱移動モデルを解く
領域とメッシュ生成後処理と可視化
熱移動モデル用語集
初期条件と境界条件参考文献

はじめに

化学工学，機械工学，航空工学等，多数の産業分野では，複数材料の層で構築された装置が使われている．複数層材料の重要な特徴の一つに，いくつかの熱特性が組み合せられるというものがある．数的なシミュレーションを実行し，そのような合成物の熱分布を得ることは，実際のシステムでどのように動作するかを知るヒントになるため重要である．

このモデルは，多層球を介した3D熱移動が解析された[Singh et al., 2016]に基づいている．問題は，熱接触における熱分析を行い，時間 t における中が空洞の同心球の n 個の層の温度場分布 T(t,x,y,z)を求めることである．立体の2D横断面は n=3の場合について示している．

1.gif

球の各層には固有の材料特性がある．例えば，i 番目の層の熱伝導率は Ki，質量密度は ρi，比熱容量はCpi である．このシミュレーションでは，3層まで，つまり中が空洞の3層の同心球までに制限するので，i =1,2,3となる．[Singh et al., 2016]で述べてあるように，異なる2種類の境界条件を課す．一つは多層球の外接面，もう一つは内接面に設定する．このモデルでは，中空の球の外接面は常に非線形熱源 Qexternal(θ,ϕ)の影響を受け，最初の層の内接面は一定温度 Tinternal=0[ºC]に設定される．シミュレーションは，球のすべての層について初期温度 Tinitial=0[ºC]で開始される．層と層の間の熱移動（各層の温度変化）は，想定通りそれぞれの熱特性に依存する．

熱移動解析の理論的情報はについては「熱移動」に記載されている．

有限要素パッケージとOpenCascadeLinkをロードする：

領域とメッシュ生成

以下のプロットは同心球の3つの層の二次元横断面である．半径 ri とそれぞれの層のモデルパラメータのラベル κi，ρi，Cpi が示されている．

2.gif

各層の半径の値を指定する：
OpenCascadeLinkの関数を使って，4つの同心球を生成する：
球を面に変換して，立体の内部境界を定義する：

次に面を結合する．これにより4つの領域0<r<r0，r0<r<r1,...を持つオブジェクトになる．この多層球は領域0<r<r0が空洞である．この条件はToElementMeshのオプションRegionHolesを使って，メッシュ生成の過程で指定することができる．

4面と4領域を持つオブジェクトを定義する：
境界メッシュを生成する：

この3つの層はマーカー1，2，3でラベル付けする．メッシュでのマーカーおよびその生成については，要素メッシュの生成チュートリアルをご参照いただきたい．

層の名前とマーカーを関連付ける：
完全なElementMeshを作成する：

熱移動モデル

温度場分布 T(t,x,y,z)を求めるために過渡熱移動方程式(1)を使う：

このPDEモデルでは，多層球の内面(r=r0)および外面(r=r3)に2つの境界条件が適用される．内面では，一定温度がTinternal=0[ºC]に等しいという表面境界条件が課せられる．外面では値 Q(θ,ϕ)=q0θ2(π-θ)2ϕ2(π-ϕ)2（θ は球座標の極角度，ϕ は方位角）の非線形熱流源が適用される．

層ごとのモデルの材料パラメータを定義する：
ヘルパー関数を定義し，ElementMarkerを使って各層の材料パラメータを指定する：

効率的なPDE係数の設定についてのこちらのヒントも参照のこと．

立体全体について熱移動モデルパラメータを定義する：
モデルのモデル変数 vars とパラメータ pars を設定する：

初期条件と境界条件

このPDEモデルでは，層1の内面は Tinternal=0[ºC]の温度面境界条件に設定される．層3の外面は，極角度 θ と方位角 ϕ に依存する，不均一な時間非依存の熱流束境界条件の影響を受ける．[Singh et al., 2016]では問題の記述は球座標で定式化されている．もとは球座標で定義された境界条件を適用するために，球座標から直交座標への変換を行う．

球座標から直交座標への変換規則を定義する：
熱源の大きさ q0 を指定し，もとの境界条件を直交座標に変換する．

まず，HeatTemperatureCondition関数を使ってDirichletConditionを設定する．

r=r0における温度を常に0[ºC]に設定する：

球の外側のNeumannValueはHeatFluxValue関数を使って生成する．熱流源を層3の外面に適用する．

r=r3の0<=θ<π，0<=ϕ<=πの範囲である外側に熱流束を生成する：
多層球の初期温度を Tinitial=0[ºC]に設定する：

以下のプロットは，問題を定式化するために使われる熱流境界条件と座標系を示している．熱流は y=r3 でピーク値の q0となり領域 y<0で0に減少する．

4.gif

PDE熱移動モデルを解く

3つの層の間の熱伝導を分析するために，モデルを tinitial=0sから tinitial=0s まで解く．

積分の初期および最終時間を定義する：
HeatTransferPDEComponentを使って解くための熱移動方程式を定義する：
熱移動PDEモデルを解き，時間／メモリ使用量を監視する：

後処理と可視化

多層球の温度場分布 T(t,x,y,z)を一次元，二次元，三次元切断で可視化する．

一次元切断：半径方向の温度分布

θ=π/2に沿った半径方向および方位角(ϕ)方向a) ϕ=0，b) ϕ=π/2，c) ϕ=π における過渡温度分布解をプロットする．

6.gif

温度をプロットする時間を定義する：

線形部分 ϕ=0& θ=π/2は x>0& y=0に，ϕ=π/2 & θ=π/2 は x=0&y>0に，ϕ=π & θ=π/2は x=0& y<0に対応する．

{x,y,-y}方向のプロットについて，多層球の半径切断を定義する：
時間の凡例値を生成する：
3つの半径切断における温度場をプロットする：
解を表示する：

二次元切断：温度等高線プロットの分布

球の z=0横断面における温度場等高線を可視化する．システム全体の温度分布は時間とともに増加する．t∈{5,10,20,30,45,50}のそれぞれの時間についてプロットに対する9個の温度値のリストを等高線として定義する．

7.gif

各時間 t に対する9個の温度等高線値のリストを6つ作成し，それを時間のリストに関連付ける：
多層球の横断面を構築する：
等高線の色スタイル：
ContourPlotを使う関数を定義し，時間 t における温度場をプロットする：
プロットを配列で表示する：
さまざまな時間における等高線プロットをすべて表示する：

三次元切断：温度密度プロットのアニメーション

球の三次元断面における温度場分布の時間アニメーションを実行する．熱源が領域 y>0から適用されると，以下のアニメーションから分かるように温度は y 方向で上昇し始める．各層の熱伝導率は κ3>κ2>κ1という関係なので，密度プロットのアニメーションで見られるように，層3（外側の層）は熱を伝達し，内側の層よりも早く温度が上昇する．

温度範囲の値，およびBarLegendによる凡例バーの形式を定義する：

SliceContourPlot3Dのプロット範囲の制限を使う．スライス面の中心の球"CenterCutSphere"がプロットされるメッシュ領域と全く同じである場合，面の一部が表示されない可能性があるからである．したがって，プロットの範囲はメッシュの範囲よりも小さくなければならない．

プロット範囲の制限と密度プロットのオプションを定義する：

以下の3Dアニメーションには上記のディスク空き容量の平均が必要な場合がある．このノートブックに必要なディスクスペースを減らすために，このノートブックでは可視化で関数Rasterizeを呼び出している．この呼出しを使った場合の欠点は，使わなかったときほどくっきり見えないということである．高品質のグラフィックスが得たい場合は，Rasterizeの呼出しを削除するかコメントアウトするかするとよい．

アニメーションを表示する：

用語集


参考文献

1.  Suneet, S., Jain, P. K., & Uddin, R. Analytical Solution for Three-Dimensional, Unsteady Heat Conduction in a Multilayer Sphere. ASME. Journal of Heat Transfer, 2016, 138(10), 101301. doi: 10.1115/1.4033536.