1次元非定常熱伝導の変数分離法

理論と物理

概要 — 変数分離法とは

• クエンチ（急冷）処理の内部温度予測 — 鋼材を水中に投入した後、中心温度がいつマルテンサイト変態点を通過するか。これを変数分離法で秒単位で予測できる
• FEM解の検証（V&V） — NAFEMS T2ベンチマーク問題は変数分離法の厳密解を参照値にしている
• 設計初期の概算 — コンクリート壁の日射応答、電子基板の突然通電時の温度上昇を、FEMを回す前に手計算で見積もる

支配方程式と無次元化

ここで $\alpha = k / (\rho c_p)$ は熱拡散率 [m²/s]、$k$ は熱伝導率、$\rho$ は密度、$c_p$ は比熱である。

解析を統一的に扱うため、無次元温度 $\theta$ と 無次元数 を導入する:

ここで $T_i$ は初期温度、$T_\infty$ は周囲流体温度、$L$ は代表長さ（平板の半厚）、$h$ は熱伝達係数である。$\theta=1$ が初期状態、$\theta=0$ が定常状態（周囲温度と一致）に対応する。

変数分離の手順

これを熱伝導方程式に代入すると:

これにより2つの常微分方程式が得られる:

境界条件と固有値方程式

厚さ $2L$ の平板で、中心 $x=0$ に対称条件、両面 $x=\pm L$ に各種境界条件を適用した場合:

固有関数の直交性を利用すると、展開係数 $C_n$ が初期条件から一意に決まる。平板で初期温度が一様 $T_i$ の場合:

最終的な解は以下の無限級数で表される:

フーリエ数と1項近似

• 第1項: $\exp(-\zeta_1^2 \times 0.2) = \exp(-0.493) \approx 0.611$ — まだ大きい
• 第2項: $\exp(-\zeta_2^2 \times 0.2) = \exp(-4.44) \approx 0.012$ — もう1%程度
• 第3項: $\exp(-\zeta_3^2 \times 0.2) = \exp(-12.3) \approx 4.5 \times 10^{-6}$ — 事実上ゼロ

だから「1項近似（one-term approximation）」で十分なんだ。ハイスラーチャートはまさにこの1項近似をグラフ化したものだよ。

1項近似を適用した中心温度 $\theta_0$（$x=0$）の式:

円柱と球への拡張

数値解法と実装

解析解の数値計算手順

1. 固有値の計算: $\zeta_n \tan(\zeta_n) = \text{Bi}$ をNewton法で解く。$n$ 番目の根は区間 $((n-1)\pi, \; (n-0.5)\pi)$ に必ず存在するので、Brent法で挟み込むのが確実
2. 展開係数の計算: $C_n = 4\sin(\zeta_n) / (2\zeta_n + \sin(2\zeta_n))$ を各 $n$ について計算
3. 級数の合算: 所望の $x$ と $t$（Fo）に対して $\theta = \sum C_n \cos(\zeta_n x/L) \exp(-\zeta_n^2 \text{Fo})$ を計算
4. 打ち切り判定: $|C_n \exp(-\zeta_n^2 \text{Fo})| < \epsilon$（例えば $10^{-10}$）になったら級数を打ち切る

from scipy.optimize import brentq
import numpy as np

def eigenvalues_plate(Bi, N=50):
    """平板・第3種BCの固有値をN個求める"""
    zeta = []
    for n in range(1, N+1):
        a = (n - 1) * np.pi + 1e-12
        b = (n - 0.5) * np.pi - 1e-12
        f = lambda z: z * np.tan(z) - Bi
        zeta.append(brentq(f, a, b))
    return np.array(zeta)

def theta_plate(x_L, Fo, Bi, N=50):
    """無次元温度 theta(x/L, Fo, Bi)"""
    zn = eigenvalues_plate(Bi, N)
    Cn = 4 * np.sin(zn) / (2*zn + np.sin(2*zn))
    return np.sum(Cn * np.cos(zn * x_L)
                  * np.exp(-zn**2 * Fo))

FEMとの比較検証

1. 粗いメッシュ×大きい時間刻みで解く
2. メッシュを2倍に細分化して再解析 → 温度変化を確認
3. 時間刻みを半分にして再解析 → 温度変化を確認
4. 両方が1%以内に収まるまで繰り返す
5. 最終結果を変数分離法の厳密解と比較

時間積分スキームの影響

実践ガイド

クエンチ（急冷）解析への応用

• $\text{Bi} = h \cdot r_0 / k = 5000 \times 0.025 / 50 = 2.5$
• 固有値方程式 $\zeta_1 J_1(\zeta_1)/J_0(\zeta_1) = 2.5$ を解くと $\zeta_1 \approx 1.814$
• 中心温度が300°Cに達する時刻: $\theta_0 = (300 - 20)/(850 - 20) = 0.337$ から Fo を逆算
• $0.337 = C_1 \exp(-\zeta_1^2 \text{Fo})$ → $\text{Fo} \approx 0.152$ → $t = \text{Fo} \cdot r_0^2 / \alpha \approx 4.8$ 秒

わずか5秒弱で中心が300°Cまで冷えることが分かる。冷却速度から焼入れ硬さを予測する連続冷却変態（CCT）図との組み合わせが実務のポイントだ。

断熱材・多層壁の設計検討

• 等価物性法: 多層壁を単一層に換算（等価 $\alpha$ を計算）して変数分離法を適用
• 各層で個別に変数分離: 界面で温度と熱流束の連続条件を課す（かなり数式が複雑になる）
• 実務ではFEMが無難: 層数が3以上なら1次元FEMの方が圧倒的に楽

単層の場合——例えば「厚さ50mmのグラスウール断熱材の片面を急に加熱したら、反対面は何分後に温まるか？」——にはそのまま使える。

商用ツールでの実装手順

ソフトウェア比較

ツール別の非定常熱伝導対応

オープンソースツール

• Python + SciPy: 前述のスクリプトで固有値計算から温度分布プロットまで完結。実質これが最も柔軟
• Wolfram Mathematica: DSolve コマンドで記号的に導出可能。学術論文用に美しい数式出力が得られる
• Octave / MATLAB: fzero で固有値を求め、besselj で円柱問題にも対応
• hteqn（Pythonライブラリ）: 平板・円柱・球の変数分離解を関数呼び出し一発で計算。ハイスラーチャート値との自動照合機能付き

先端技術

PINNと変数分離法の融合

• 空間方向のサブネットワーク: 固有関数に近い基底を学習
• 時間方向のサブネットワーク: 指数減衰を学習
• 結果: 従来のPINNより桁違いに少ないデータで収束し、物理的に正しい解を保証できる

変数分離法の数学的構造が、200年の時を経てニューラルネットの「帰納バイアス」として復活したわけだ。

非フーリエ伝導への拡張

トラブルシューティング

級数が収束しない場合

FEMと解析解が合わない場合

1. 変数分離法側のチェック:
   • 固有値を十分な個数（少なくとも20個）計算したか？
   • $L$ の定義は半厚か全厚か？ 教科書によって異なる
   • 展開係数 $C_n$ の式は境界条件に対応したものを使っているか？
2. FEM側のチェック:
   • 物性値を温度非依存に設定したか？ 温度依存テーブルが紛れ込んでいないか？
   • メッシュ収束は確認したか？ 特に表面近く（大きな温度勾配）は要細分化
   • 時間刻みは十分小さいか？ 自動時間刻みの最大値制限は適切か？
   • 対称面の境界条件は正しいか？ 断熱≠自由面
3. 両者が一致しない場合: まず第1種BC（表面温度固定）で検証。これは最も単純で間違えにくい。第3種BCの検証は後回しにする