熱拡散方程式の解析解 (周期解)

Introduction

小惑星の表面など, 外部から変動するエネルギーが入力されることによって温度が変化する系を考えます. まずは単純な系として物質からの熱放射については考えず, 無限に大きな物体の表面付近の熱変化を考えることにします.

物質の熱的性質を決定する物理量として以下があります.

\(\kappa\): thermal conductivity (W/m/K)
\(c_\mathrm{p}\): specific heat capacity (J/kg/K)
\(\rho\): density (kg/m³)

またこれらの積で定義される量に熱慣性 (thermal inertia) があります.

\(\Gamma = \sqrt{\kappa c_\mathrm{p} \rho}\): thermal inertia (tiu)¹

ここでは物体の表面温度と熱慣性の関係を整理します.

Thermal diffusion equation

ある物質の熱的性質が \(\kappa\), \(c_\mathrm{p}\), \(\rho\) で記述できるとします. この物質で構成された物体の表面付近の深さ方向の温度分布を考えます. もし物体が十分に大きく均質であると仮定すれば, 対称性から物体の表面に沿う方向の熱輸送については無視してよいでしょう. ここでは深さ方向の熱輸送だけを考えることにします. 物体の表面から深さ方向に \(z\) 軸をとります. また, 物体の表面を \(z = 0\) と定義します.

単位時間 (単位面積当たり) に深さ \(z\) のレイヤーから深さ \(z + \mathrm{d}z\) のレイヤーに向けて輸送される熱量は熱伝導方程式 (thermal conduction equation) によって以下のように記述できます. ここで \(\mathrm{d}z\) が十分小さいとして, 差分を微分に置き換えました. 熱勾配が深さ方向に負であれば, 深い方向に熱エネルギーが流れるので係数は負になります.

\[ \Delta Q_{z \to z + \mathrm{d}z} \left( = \kappa \frac{T(z) - T(z + \mathrm{d}z)}{\mathrm{d}z} \right) = - \kappa \frac{\mathrm{d}T}{\mathrm{d}z}(t, z). \]

また, 深さ \(z\) のレイヤーが単位時間に受け取る正味の熱量は, 上層とやりとりする熱エネルギーと下層とやりとりする熱エネルギーの和になります. 便宜的に上層からもらって下層に渡すような表記をすると, 以下のように記述できます.

\[ \Delta Q(t, z) = \Delta Q_{z - \mathrm{d}z \to z} - \Delta Q_{z \to z + \mathrm{d}z}. \]

考えているレイヤーの厚みを \(\mathrm{d}z\) とすると単位面積あたりの熱容量は \(\rho \, c_\text{p} \, \mathrm{d}z\) になります. 受け取った正味の熱量を熱容量で割ることでレイヤーの温度変化がわかります. よって深さ \(z\) のレイヤーの熱変化量は以下の式で記述できるとわかります.

\[ \frac{\mathrm{d}T}{\mathrm{d}t}(t, z) = \frac{\Delta Q(t, z)}{\rho \, c_\text{p} \, \mathrm{d}z} \frac{\kappa}{\rho c_\text{p}} \frac{1}{\mathrm{d}z}\left( \frac{\mathrm{d}T}{\mathrm{d}z}(t, z) - \frac{\mathrm{d}T}{\mathrm{d}z}(t, z-\mathrm{d}z) \right) \simeq \frac{\kappa}{\rho c_\text{p}} \frac{\mathrm{d}^2T}{\mathrm{d}z^2}. \]

熱拡散係数 \(\alpha^2 = \kappa / \rho c_\text{p}\) とおくと上記の式は以下のように記述できます. これを熱拡散方程式 (thermal diffusion equation) と呼びます.

\[ \frac{\mathrm{d}T}{\mathrm{d}t}(t, z) = \alpha^2 \frac{\mathrm{d}^2T}{\mathrm{d}z^2}. \]

この方程式は自明な解のひとつとして \(T(t, z) = T_0\) (定数解) を持ちます. 物体が等温である場合には, 熱的平衡状態に達しておりもやは熱輸送はおこなわれないことに対応します. 熱拡散方程式の解に定数を足してオフセットさせても解になります. そこで, 以下ではこの自由度については無視して, 温度は平衡温度 (平均温度) からのズレとして定義し直すことにします.

Analytical solution

この物体の温度は外部からの熱的な介入によって変動を受けます. 外部からの介入に対して物体の温度分布がどのような応答をするかを調べます. そこで表面に \(E(t)\) という (単位面積あたりの) エネルギー入力 (変動) を与えることを考えてみます. ここではエネルギー入力の大きさは周期 \(P\) (角周波数 \(\omega\)) で正弦波的に変動させることとします.

\[ E(t) = E_0 \exp (\mathrm{i}\omega t). \]

解析のしやすさのために複素数で定義しました. 外部からの介入が周期 \(P\) で変動するため, 物体の温度分布も同様な周波数で変動することが期待されます. そこで暫定的な解として, 時間方向には角周波数 \(\omega\) で, 空間方向には空間周波数 \(k\) で変動する関数系を仮定します.

\[ T(t, z) = T_0 \exp\left(\mathrm{i}\omega t + \mathrm{i}k z\right). \]

熱拡散方程式に代入して実部と虚部の係数を比較すると, 方程式を満たすための条件として以下を得ます.

\[ k_\text{r}^2 = k_\text{i}^2, \qquad \omega = - 2\alpha^2 k_\text{r}k_\text{i} \]

ここで空間周波数 \(k\) を \(k = k_\text{r} + \mathrm{i}k_\text{i} ~ (k_\text{r}, k_\text{i} \in \mathbb{R})\) と分離すると, それぞれについて以下の条件を得ます.

\[ k_\text{r}^2 = k_\text{i}^2, \qquad \omega = - 2\alpha^2 k_\text{r}k_\text{i} \]

角周波数については符号がどうであれ一般性は失わないため, \(\omega > 0\) と決めておくことにします. まだ解は一意に定まりませんが, 解が \(z\) 方向に発散しないためには \(k_\text{i} > 0\) でなければなりません. よって \(k\) は以下の式で記述できます.

\[ k = -\sqrt{\frac{\omega}{\alpha^2}} + \mathrm{i}\sqrt{\frac{\omega}{\alpha^2}} = -\sqrt{\frac{\rho c_\text{p}}{\kappa}\frac{\pi}{P}} + \mathrm{i}\sqrt{\frac{\rho c_\text{p}}{\kappa}\frac{\pi}{P}}. \]

以上より周期 \(P\) で振動する擾乱下での一般解は以下のように書けるとわかりました.

\[ T(t, z) = T_0 \exp\left( \mathrm{i}\omega t - \mathrm{i}\sqrt{\frac{\rho c_\text{p}}{\kappa}\frac{\pi}{P}} z - \sqrt{\frac{\rho c_\text{p}}{\kappa}\frac{\pi}{P}} z \right). \]

Note – Thermal Penetration Depth

一般解 \(T(t, z)\) は時間･空間的に振動する項と \(z\) 方向に減衰する項の積になっていることがわかります. ここで, 減衰項に注目すると温度の振幅は係数の逆数を特徴的なスケールとして減衰することがわかります. これを熱浸透深さ (thermal penetration depth) と呼びます.

\[ z_\text{tpz} = \sqrt{\frac{\kappa}{\rho c_\text{p}}\frac{P}{\pi}} \]

以下では, 係数 \(T_0\) を決めるために境界条件に注目します. \(z=0\) の面では外部からのエネルギー入力として \(E(t)\) が与えられています. 物体表面における輻射によるエネルギーの散逸は無視すると, エネルギーの保存則から表面層が得たエネルギーはすべて物体の内部へと流れていくと考えられます. よって, 表面の熱伝導方程式から以下の条件を得ます.

\[ E(t) = E_0 \exp(\mathrm{i}\omega t) = \Delta Q_{z \to z+\mathrm{d}z}|_{z=0} = - \kappa \frac{\mathrm{d}T}{\mathrm{d}z}(t, 0) \]

一般解がこの式を満たすように \(T_0\) を定めます. 一般解を \(z\) で微分して \(z=0\) を代入すると以下の式を得ます.

\[ \frac{\mathrm{d}T}{\mathrm{d}z}(t, 0) = T_0 \sqrt{\frac{\rho c_\text{p}}{\kappa}\frac{\pi}{P}} \left( - 1 - \mathrm{i} \right) \exp\left(\mathrm{i}\omega t \right). \]

熱伝導方程式に代入することで複素振幅としての係数 \(T_0\) が決定できます. ここで熱慣性 \(\Gamma = \sqrt{\kappa \rho c_\text{p}}\) が現れます.

\[ T_0 = \sqrt{\frac{\kappa}{\rho c_\text{p}}\frac{P}{\pi}} \frac{E_0}{\kappa(1 + \mathrm{i})} = \sqrt{\frac{P}{2\pi}} \frac{E_0}{\Gamma} \exp\left( -\mathrm{i} \frac{\pi}{4} \right) \]

また, 外部からのエネルギー入力 \(E(t)\) に応答したときの温度構造 \(T(t, z)\) は以下のとおりに求まります. \(T_0\) の偏角はどこに入れても良かったのですが, ここでは時間 \(t\) の位相として表現するのがわかりやすいと考えて整理しました. 時間にして周期の 1/8 の位相差に相当します.

\[ T(t, z) = \sqrt{\frac{P}{2\pi}} \frac{E_0}{\Gamma} \exp\left( \mathrm{i}\omega\left(t - \frac{P}{8} \right) - \mathrm{i}\sqrt{\frac{\rho c_\text{p}}{\kappa}\frac{\pi}{P}} z - \sqrt{\frac{\rho c_\text{p}}{\kappa}\frac{\pi}{P}} z \right) \]

温度変化 \(T(t, z)\) の振幅は \(E_0/\Gamma \sqrt{P}\) に比例することがわかります. 単位時間あたりの熱流入が大きければ温度変化は比例して増加します. また, 物質の性質に依らず周期に対して \(\sqrt{P}\) の依存性を持ちます. 外部エネルギーの変動がゆっくりであるほど, 長時間エネルギーを受け続けるので温度の振幅は大きくなりますが, 周期そのものではなく周期の平方根に比例して変化します.² 物質の熱的性質に関わる物理量はすべてまとめて \(\Gamma\) という形で現れます. これは物理量 \((\kappa, c_\mathrm{p}, \rho)\) がそれぞれ異なっていても, \(\Gamma\) さえ同じなら物体の温度構造は同様になることを意味しています. 逆にいうと, 温度を測定するだけでは物理量 \((\kappa, c_\mathrm{p}, \rho)\) を分離することはできず, 熱慣性 \(\Gamma\) だけがわかるとも言えます.

得られた結果を図示して確認してみます. それぞれの値は以下のように規格化しています.

\[ \kappa = 1.0, \quad c_\mathrm{p} = 1.0, \quad \rho = 1.0, \quad E_0 = 1.0, \quad P = 1.0. \]

表面温度 (赤実線) と表面へのエネルギー入力 \(E(t)\) (緑一点破線) の時間変化を図示しました. 横軸は規格化された時間を表しています. \(E(t)\) は時刻 0 でピークとなるように設定しています. 表面温度に注目すると, 温度のピークは \(t = 0.125\) にあらわれており, 周期 \(P\) の 1/8 だけ遅れていることがわかります. この遅れは外部からのエネルギー入力の振幅 \(E_0\) や, 物質の熱的性質を決めるパラメタとは無関係に決まっています.

温度構造の動的な変化を見るために横軸に深さをとって温度変化をアニメーションで図示しました. 先ほどと同様に温度を赤実線で, 熱の深さ方向への輸送量を緑一点破線で示しました. \(z=0\) にある緑色の丸は \(E(t)\) の値を表しています.³ 温度構造が熱輸送を追いかけるように変化していることがわかります.

熱慣性の単位 (\(\mathrm{J\,K^{-1}\,m^{-2}\,s^{-0.5}}\)) はあまり物理的な意味をなさないので, まとめて tiu (thermal inertia unit) と表現されることがあります. ↩
熱浸透深さ \(z_\text{tpd}\) が \(\sqrt{P}\) に依存するため, 長時間あぶられたときにエネルギーがより深くまで伝わります. 結果としてより深い領域まで温度が上昇する (実効的な熱容量が \(z_\text{tpd}\) に比例して増える) ことになり, 温度の振幅と周期の関係は \(PE_0/z_\text{tpd} \propto \sqrt{P}\) のように整理できます. ↩
図示したすべての範囲･時間で熱拡散方程式を満たしているか確認するために, 熱拡散方程式の両辺の差分を灰色の点線で示しています. 差分はゼロのまま変動していないので熱拡散方程式を満たす解であること確認できます. ↩