半無限固体の非定常熱伝導

理論と物理

概要と適用場面

• 溶接の入熱直後：アーク溶接やレーザー溶接で母材に入った熱が拡がっていく初期段階。HAZ（熱影響部）の温度分布はまさに半無限体の解で見積もれる
• レーザー加工の表面温度：レーザーが照射された瞬間から数ミリ秒間の表面温度上昇。加工深さの制御に直結する
• 地中の凍結深さ：冬季に地表面温度が急降下したとき、凍結がどこまで進むかを推定する。水道管の埋設深さの設計根拠になる
• 鋳造・鍛造の急冷：高温の金属を水や油で焼入れしたときの表面近傍の冷却速度。マルテンサイト変態が起きるかどうかの判定に使う

支配方程式と相似変数

ここで $\alpha = k/(\rho c_p)$ は熱拡散率 [m²/s] だ。鋼で約 $1.2 \times 10^{-5}$、アルミで約 $9.7 \times 10^{-5}$ m²/s。この値が大きいほど温度変化が速く伝わる。

• 初期条件：$T(x, 0) = T_i$（全体が均一温度）
• 表面境界条件：$T(0, t) = T_s$（$t > 0$ で表面温度が一定に変化）
• 遠方条件：$T(\infty, t) = T_i$（十分遠方は初期温度のまま）

そして決定的なのが相似変数の導入だ：

この変数変換をすると、偏微分方程式が常微分方程式に帰着される。$x$ と $t$ という2つの独立変数が $\eta$ 1つにまとまるんだ。これをボルツマン変換と呼ぶ。

解析解の導出

境界条件は $\theta(0) = 0$, $\theta(\infty) = 1$。これは誤差関数 erf で解けて：

元の温度に戻すと、半無限固体の温度分布の最も重要な公式が得られる：

ここで $\operatorname{erfc}(\eta) = 1 - \operatorname{erf}(\eta)$ は相補誤差関数だ。

つまり $\eta = 2$（$x = 4\sqrt{\alpha t}$）より奥ではほぼ温度変化がない。これが浸透深さの根拠になる。

浸透深さと表面熱流束

例えば鋼（$\alpha \approx 1.2 \times 10^{-5}$ m²/s）に溶接で瞬間的に入熱した場合：

• $t = 1$ 秒後：$\delta \approx 4\sqrt{1.2 \times 10^{-5}} \approx 14$ mm
• $t = 10$ 秒後：$\delta \approx 44$ mm
• $t = 100$ 秒後：$\delta \approx 139$ mm

板厚がこの $\delta$ より十分大きければ、半無限体として扱ってOK。板厚 20 mm の鋼板なら、溶接後 $t \lesssim 2$ 秒程度は半無限体近似が有効だ。

$1/\sqrt{t}$ に比例するから、初期に猛烈な熱流束が必要で、時間とともに減衰する。$t \to 0$ で $q \to \infty$ になるのは、初期の温度勾配が表面で無限大になることに対応している。実際にはどんな熱源にも有限の出力があるから、ごく初期は一定温度条件より一定熱流束条件のほうが現実的だね。

境界条件の3パターン

Case 2 の表面温度は $T(0,t) = T_i + \frac{2q_0''}{k}\sqrt{\frac{\alpha t}{\pi}}$ で、$\sqrt{t}$ に比例して上昇する。レーザー加工では出力が一定なので Case 2 のほうが実態に近いことが多い。

半無限体近似の妥当性判定

• $\mathrm{Fo} < 0.05$：半無限体近似の誤差は 1% 未満。安心して使える
• $0.05 < \mathrm{Fo} < 0.2$：裏面の影響が出始める。精度に注意
• $\mathrm{Fo} > 0.2$：もはや半無限体ではない。有限厚さの解（フーリエ級数解）を使うべき

例えば厚さ 50 mm の鋼板（$\alpha = 1.2 \times 10^{-5}$）なら $\mathrm{Fo} = 0.05$ に対応する時間は $t = 0.05 \times 0.05^2 / (1.2 \times 10^{-5}) \approx 10.4$ 秒。つまり溶接後 10 秒以内なら半無限体として計算して問題ない。

数値解法と実装

解析解が使えない場面

• 温度依存物性：鋼の熱伝導率は 0°C で約 60 W/(m·K)、800°C で約 30 W/(m·K) まで低下する。非線形になるので解析解はない
• 2次元・3次元効果：溶接ビードのように加熱領域が有限幅なら、横方向への熱拡散を考慮する必要がある
• 時間変化する境界条件：溶接トーチの移動、レーザーのパルス照射など
• 相変化：溶融池の形成や凝固を伴う場合（Stefan問題）
• 複合材料・異種材接合：界面での熱抵抗が存在する場合

FEMによる半無限固体の離散化

1. 十分大きな解析領域を取る：浸透深さ $\delta = 4\sqrt{\alpha t_{\max}}$ の2〜3倍の領域を確保し、遠方境界に $T = T_i$（Dirichlet条件）を設定する。最もシンプルで確実
2. 無限要素を使う：AnsysのINFIN111やAbaqusのCINPE4など、遠方境界の減衰を内蔵した特殊要素を境界に配置する。解析領域を大幅に小さくできる

熱伝導のFEM定式化では、Galerkin法により以下の離散系が得られる：

ここで $[C]$ は熱容量マトリクス、$[K]$ は熱伝導マトリクス、$\{F\}$ は熱負荷ベクトルだ。構造解析の $[M]\{\ddot{u}\} + [K]\{u\} = \{F\}$ と形が似ているだろう？ 1階の時間微分か2階かの違いだ。

メッシュ戦略

Ansys なら Bias Factor を 10〜20 程度に設定すると、自動的にいい感じの不等分割になる。Abaqus では edge seed の single bias を使う。

時間積分スキーム

• 初期の時間刻み：$\Delta t_1 = \Delta x_{\min}^2 / (6\alpha)$ 程度（表面最小要素の熱拡散時間の1/6）
• 時間刻みの成長：各ステップで 1.2〜1.5 倍ずつ増やしていく（幾何級数的増大）
• 自動時間刻み制御：商用ソルバーの自動制御に任せるのが最も安全。Abaqus の *HEAT TRANSFER ステップや Ansys の AUTO TIME STEPPING を有効にする

陰解法（backward Euler, Crank-Nicolson）なら安定性の制約はないが、精度のために初期の細かい刻みは依然として必要だ。

無限要素の活用

無限要素を使えば、解析領域を浸透深さの 1.5 倍程度まで縮小でき、要素数を大幅に削減できる。ただし、無限要素の配置方向が誤っていると意味のない結果になるから、必ず理論解と比較して検証すること。

実践ガイド

溶接HAZの温度予測

入熱量 $Q$ [J/m] の線熱源が厚板表面を移動する場合、ローゼンタールの移動点熱源解から導かれる冷却速度の近似式：

ただし簡略版としては、鋼板厚 $d$ > 浸透深さのとき：

例えば、入熱 $Q = 1.5$ kJ/mm、予熱なし（$T_0 = 20$°C）、板厚 25 mm の軟鋼なら $t_{8/5} \approx 12$ 秒。これが 6 秒以下だと水素割れリスクが高まり、30 秒以上だと靱性が低下する。溶接施工条件の設計に直結する数値だ。

レーザー加工の表面温度

例えば、SUS304（$k = 16$ W/(m·K)、$\alpha = 4.2 \times 10^{-6}$ m²/s）にパワー密度 $I_0 = 10^8$ W/m² のレーザーを 1 ms 照射すると：

$T(0, 1\text{ms}) = 20 + \frac{2 \times 10^8}{16}\sqrt{\frac{4.2 \times 10^{-6} \times 10^{-3}}{\pi}} \approx 20 + 1,450 = 1,470$°C

融点 1,400°C を超えるから溶融が始まる。照射時間を変えることで加工深さを制御できるわけだ。ただし融解が始まると潜熱の効果で実際の温度上昇は鈍化するから、精密な予測には相変化モデル（Stefan問題）が必要になる。

地中の凍結深さ推定

$\operatorname{erfc}(\eta_f) = 5/20 = 0.25$ → $\eta_f \approx 0.81$

土壌の熱拡散率を $\alpha \approx 5 \times 10^{-7}$ m²/s とすると、1週間（$t = 6 \times 10^5$ 秒）で：

$x_f = 2 \times 0.81 \times \sqrt{5 \times 10^{-7} \times 6 \times 10^5} \approx 0.89$ m

つまり約 90 cm まで凍結が進む。水道管の凍結防止には、この深さ以上に埋設する必要がある。日本の寒冷地の水道管埋設深さが 1.0〜1.2 m なのは、まさにこの計算に基づいているんだ。ただし実際には土壌の水分凍結に伴う潜熱放出があるため、凍結進行は erfc 解より遅くなる。より精密にはノイマン解（Stefan問題の解析解）を使う。

接触温度の計算

ここで $e = \sqrt{k \rho c_p}$ は熱浸透率（thermal effusivity）。単位は J/(m²·K·s^{1/2})。

例えば 200°C の鋼（$e \approx 14,000$）と 20°C の銅（$e \approx 37,000$）を接触させると：

$T_c = (14000 \times 200 + 37000 \times 20)/(14000 + 37000) \approx 69$°C

接触面温度は銅寄りになる——熱浸透率が大きい方が「温度を押し付ける力」が強いんだ。人がステンレスの手すりを冷たく感じ、木の手すりを温かく感じるのも、同じ原理。手の熱が手すりに奪われるかどうかは、手すり材料の熱浸透率で決まる。

ソフトウェア比較

商用ツールでの設定方法

1. 一定温度境界の Case 1 で FEM 結果と erfc 解を比較
2. 表面温度の時間履歴、深さ方向の温度分布の両方で一致を確認
3. 表面熱流束 $q_s'' = k(T_s - T_i)/\sqrt{\pi\alpha t}$ との比較（数値微分の精度確認）
4. メッシュ収束性：3水準以上のメッシュ密度で解が収束することを確認

NAFEMS の T2 ベンチマーク（1次元非定常熱伝導）がまさにこの問題だ。

オープンソースでの実装

解析解の例（Python）：

import numpy as np
from scipy.special import erfc

alpha = 1.2e-5   # 鋼の熱拡散率 [m²/s]
Ti, Ts = 20, 800  # 初期温度, 表面温度 [°C]
t = 5.0           # 時間 [s]
x = np.linspace(0, 0.05, 100)  # 深さ [m]

T = Ti + (Ts - Ti) * erfc(x / (2 * np.sqrt(alpha * t)))
delta = 4 * np.sqrt(alpha * t)  # 浸透深さ
print(f"浸透深さ: {delta*1000:.1f} mm")

FEniCS で同じ問題を解くなら、1次元メッシュ + Dirichlet BC + backward Euler で数十行のコードで済む。商用ソルバーの結果検証に活用できるよ。

先端技術

非フーリエ熱伝導

このスケールではCattaneo-Vernotte 方程式（双曲型熱伝導方程式）が必要：

$\tau_q$ は熱緩和時間（金属で $10^{-11}$〜$10^{-13}$ 秒程度）。この方程式は波動方程式と拡散方程式のハイブリッドで、有限速度の「熱波」が伝播する。フェムト秒レーザー加工、半導体のホットスポット解析、極低温物理で重要。

逆問題への応用

• 鋳造の界面熱伝達係数の同定：鋳型内の温度履歴から溶湯-鋳型間の熱伝達係数を逆算
• ロケットエンジンの熱負荷推定：壁面温度の計測から表面熱流束を逆算
• 急冷焼入れの冷却曲線解析：内部温度の計測から表面の冷却速度を推定

半無限固体の解析解があるおかげで、逆問題のカーネル関数（グリーン関数）が既知となり、正則化法（Tikhonov 正則化など）との組み合わせで安定した逆解析が可能になる。機械学習（PINN: Physics-Informed Neural Network）と組み合わせた最新の手法も研究が進んでいる。

トラブルシューティング

よくある失敗と対策

検証用ベンチマーク

1. NAFEMS T2：1次元非定常熱伝導。一定温度境界条件で 32 秒後の深さ 0.01 m の温度を検証。正解は erfc 解で計算可能
2. Case 1 自作ベンチマーク：鋼（$k=50$, $\rho=7850$, $c_p=500$）で $T_i=20$, $T_s=500$°C。$t=10$ s 後の温度分布を erfc 解と比較。$x=5$ mm での正解は $T \approx 434$°C
3. Case 2（一定熱流束）：同じ物性で $q_0''=10^6$ W/m²。表面温度の時間履歴を解析解と比較
4. 表面熱流束の検証：Case 1 で FEM の表面熱流束出力と $k(T_s-T_i)/\sqrt{\pi\alpha t}$ を比較。数値微分の精度チェック

いずれも2%以内の誤差で一致すれば、メッシュと時間刻みは十分だと判断できる。