8-3. 熱伝導方程式

熱伝導方程式は、物体内の温度分布の時間変化を記述する偏微分方程式である。熱の伝わり方を理解し、予測するために不可欠なツールであり、工学、物理学、材料科学など、さまざまな分野で応用される。制御工学では、温度制御や熱管理が必要なシステムにおいて、熱伝導方程式を使ってプロセスの動作をモデル化し、制御設計を行う。

熱伝導は、高温の物体から低温の物体へと熱エネルギーが移動する現象である。この熱の移動は、原子や分子の振動、自由電子の運動などによって起こる。熱伝導の基本的な法則は、「熱流束（単位時間あたりに単位面積を通過する熱量）は、温度勾配に比例し、熱伝導率に依存する。」というフーリエの法則として知られている。

１次元の熱伝導方程式

熱伝導方程式の最も基本的な形は式(1)となる。$$\frac{\partial u(x,t)}{\partial t} = \alpha \frac{\partial^2 u(x,t)}{\partial x^2} \;\;\; \cdots (1)$$ここで、\(u(x,t)\)は位置\(x\)と時間\(t\)における温度、\(\alpha = \frac{k}{\rho c}\)は熱拡散係数（\(k\)は熱伝導率、\(\rho\)は密度、\(c\)は比熱容量）である。この式は、温度変化の時間微分が温度の空間２階微分に比例することを表している。つまり、温度勾配が大きいほど、温度は早く平衡に向かうということを表している。

３次元の熱伝導方程式

3次元空間での熱伝導方程式は式(2)のように拡張される。$$\frac{\partial u(x,y,z,t)}{\partial t} = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \right)$$また、熱源や外部入力がある場合は、熱源項 \(f(x,y,z,t)\)を加える。$$\frac{\partial u}{\partial t} = \alpha \nabla^2 u + f(x,y,z,t)$$ここで、\(\nabla^2\)はラプラシアン、\(f(x,y,z,t)\)は外部熱源である。
この方程式は、熱が高温部から低温部へと拡散していく過程を示している。温度分布は時間とともに平滑化され、最終的には空間的に一様な平衡温度に達することになる。

熱伝導方程式の解

熱伝導方程式を解くには、以下の条件が必要である。
１）初期条件：時刻 \(t = 0\)の温度分布を指定する。$$u(x, 0) = u_0(x)$$これは、物体の初期状態を表し、時間発展の出発点となる。\(u(x, 0)\)は 時刻\(t = 0\)における位置 \(x\) の温度、 \(u_0(x)\)は、位置\(x\)の初期温度分布を表す関数である。
初期条件の例としては、
・物体全体が均一な温度の場合は、\(u(x, 0) = T_0\) (\(T_0\) は定数)
・物体の一部が加熱または冷却されている場合は、\(u(x, 0) = \text{特定の関数}\)
などである。
２）境界条件： 領域の端での温度や熱流束を指定する。
・ディリクレ境界条件（温度指定）：境界における温度を指定する。$$u(a, t) = T_a$$​ \(a\)は境界の位置、\(T_a\)は時間 \(t\) における境界の温度である。
・ノイマン境界条件（熱流束指定）：境界における熱流束 (温度勾配) を指定する。$$-k\frac{\partial u}{\partial n}(b, t) = q_b$$ \(k\)は熱伝導率、\(n\) は境界\(b\)の法線方向、\(q_b\)は時間\(t\)における境界の熱流束である。
・ロビン境界条件（熱伝達境界）：境界における温度と熱流束の関係を指定する。$$-k\frac{\partial u}{\partial n}(b,t) = \alpha\left(u(b,t)- u(\infty)\right)$$ \(\alpha\)は熱伝達率、\(u(\infty)\)は周囲の流体の温度である。

熱伝導方程式を解くためには、初期条件と境界条件を適切に組み合わせる必要がある。つまり、
・初期条件：物体の初期状態を指定
・境界条件：物体と外部環境との相互作用を指定
なので、これらの条件を適切に設定することで、物理的に意味のある解を得ることができる。

熱伝導方程式の解法は、問題の条件（形状、境界条件、初期条件など）によって大きく異なる。代表的な解法は以下である。
１．解析的解法
解析的解法は、数学的な手法を用いて厳密解を求める方法で、比較的単純な形状や境界条件の問題に適用できる。というか、単純な問題でしか解析的には解が求まらない。
・変数分離法：熱伝導方程式を時間と空間の関数に分離し、それぞれの関数に対する微分方程式を解く方法。主に、直方体や円柱などの単純な形状の問題に適用される。
・ラプラス変換：時間に関する微分方程式を代数方程式に変換し、解を求めた後、逆変換によって元の解を得る方法。時間的に変化する境界条件や熱源を持つ問題に有効である。
・グリーン関数法：点熱源に対する解（グリーン関数）を用いて、任意の熱源や境界条件を持つ問題の解を求める方法。複雑な形状や境界条件の問題にも適用できるが、グリーン関数を求めるのが難しい場合がある。

２．数値解法
数値解法は、コンピュータを用いて近似解を求める方法で、複雑な形状や境界条件の問題に適用できる。実用的な方法であるが、求めた解の適正を慎重に判断する必要がある。特に、計算精度や安定性に注意が必要である。また、計算結果の妥当性を検証するために、解析解や実験結果と比較することが肝要である。
・差分法：微分方程式を差分方程式で近似し、格子状の点における温度を時間発展的に計算する方法。比較的単純なアルゴリズムで実装できるが、計算精度や安定性に注意が必要である。
・有限要素法：対象を小さな要素に分割し、各要素における温度分布を近似的に計算する方法。複雑な形状や境界条件の問題に柔軟に対応できるが、計算コストが高くなる。
・有限体積法：対象を小さな体積要素に分割し、各要素における熱エネルギーの収支を計算する方法。流体解析と相性が良く、流体と熱伝導が連成する問題に有効である。

熱伝導方程式の解析例

「無限長の棒に対する１次元熱伝導方程式」を考える。
無限に長い棒（十分に長い棒）があり、時間とともに温度が変化するとする。熱伝導方程式は式(2)で表される。$$\frac{\partial u(x,t)}{\partial t} = \alpha \frac{\partial^2 u(x,t)}{\partial x^2}$$ここで、\(\alpha\)は熱拡散係数 (\(\alpha = \frac{k}{\rho c}\)​)、\(u(x,t)\)は位置\(x\)と時間\(t\)における温度、である。
初期条件として、棒の初期温度分布を与える。$$u(x, 0) = f(x)$$
境界条件として、無限長の棒なので境界条件は自然に消える。（物理的には \(|x| \to \infty\)で温度は有限である。）
以上より、温度分布\(u(x,t)\)を時間の関数として求めることを考える。
分離変数法を使う。解を次の形に仮定する。$$u(x,t) = X(x)T(t)$$これを元の偏微分方程式に代入すると$$X(x)\frac{dT(t)}{dt} = \alpha T(t)\frac{d^2 X(x)}{dx^2}$$両辺を\(X(x)T(t)\)で割ると$$\frac{1}{T(t)}\frac{dT(t)}{dt} = \alpha \frac{1}{X(x)}\frac{d^2 X(x)}{dx^2} = - \lambda$$ここで、分離定数 \(\lambda\)を導入している。これで次の2つの常微分方程式に分解できる。
・空間方向の方程式$$\frac{d^2 X(x)}{dx^2} + \frac{\lambda}{\alpha} X(x) = 0$$
・時間方向の方程式$$\frac{dT(t)}{dt} + \lambda T(t) = 0$$
空間方程式は2階の線形常微分方程式で、一般解は$$X(x) = A \cos(kx) + B \sin(kx) \quad (\text{ただし} \, k = \sqrt{\frac{\lambda}{\alpha}})$$時間方程式は指数関数解を持つので$$T(t) = C e^{-\lambda t}$$従って、解は次のように表せる。$$u(x,t) = \left[A \cos(kx) + B \sin(kx)\right] e^{-\alpha k^2 t}$$
初期条件\(u(x, 0) = f(x)\)を満たすように、解をフーリエ級数展開する。$$u(x,t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \hat{f}(k) e^{jkx} e^{-\alpha k^2 t} \, dk$$ここで、\(\hat{f}(k)\)は初期温度分布のフーリエ変換$$\hat{f}(k) = \int_{-\infty}^{\infty} f(x) e^{-jkx} \, dx$$である。

具体例として、初期温度分布が局所的なガウス関数だとする。$$f(x) = e^{-x^2}$$そのフーリエ変換は、$$\hat{f}(k) = \sqrt{\pi} e^{-k^2/4}$$これを温度分布の一般解に代入すると$$u(x,t) = \frac{1}{\sqrt{4\alpha t + 1}} e^{-\frac{x^2}{4\alpha t + 1}}$$つまり、時間が経つと温度分布はガウス関数の形を保ったまま広がり、ピーク温度が低下する。これは熱拡散の典型的な挙動である。

Scilabによる数値シミュレーション

無限長の棒に対する１次元熱伝導をScilabで数値シミュレーションする。

有限差分法 を使い、空間2階微分を中心差分で離散化する。$$\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \quad \rightarrow \quad u_i^{n+1} = u_i^n + s \left( u_{i+1}^n - 2u_i^n + u_{i-1}^n \right)$$