1次元非定常熱伝導（半無限体）

表面温度を瞬時に $T_s$ に変化させた半無限体の温度応答で、熱解析ソルバーの過渡解析機能を検証する定番問題だ。誤差補関数 erfc を含む厳密解がある。NAFEMS T1やT2ベンチマークの基盤となる理論だ。

まさにそうだ。熱伝導のCode Verificationの第一歩。空間離散化と時間離散化の両方の精度を検証できる点が構造問題にない特徴だ。

フーリエの熱伝導方程式（1次元）は

ある。表面の熱流束は

$T_i = 0$°C、$T_s = 100$°C、$k = 50$ W/(m·K)、$\rho = 7800$ kg/m³、$c_p = 500$ J/(kg·K)。$\alpha = 1.282 \times 10^{-5}$ m²/s。

AbaqusのHEAT TRANSFERステップは後退Eulerがデフォルト。TRANSIENT HEAT TRANSFERでAlpha（AMPLITUDE parameter）を設定可能。$\alpha = 0$ が後退Euler、$\alpha = 0.5$ がCrank-Nicolson。

半無限体問題では温度浸透深さ $\delta(t) = 3.6\sqrt{\alpha t}$ が時間とともに増大するから、表面近傍にメッシュを集中させる。

もちろん。空間メッシュを固定して$\Delta t$を系統的に変えれば、時間方向のRichardson外挿が可能だ。同様に、$\Delta t$を固定して空間メッシュを変えれば空間方向のGCIが得られる。両方を独立に評価するのがASME V&V 20の推奨手順だ。

Abaqus: DC1D2（1D熱伝導要素）またはDC2D8（2D）で解く。INITIAL CONDITIONS, TYPE=TEMPERATURE で初期温度、BOUNDARY で表面温度を固定。

OpenFOAMは主にCFD用だが、laplacianFoamは純粋な拡散方程式ソルバーだから熱伝導にそのまま使える。ただし固体の熱伝導だけを解くならCalculiXやCode_Asterの方が自然だ。流体と固体の連成（共役熱伝達: CHT）ではOpenFOAMのchtMultiRegionFoamが活躍する。

NAFEMS T1は「1次元非定常熱伝導」で、以下の条件だ。

時間刻みが粗すぎる場合と、メッシュが粗すぎる場合。特に後退Eulerで大きな$\Delta t$を使うと、温度の浸透速度が数値拡散で過大評価される。Crank-Nicolsonで$\Delta t$を$t_{final}/100$以下にすれば十分な精度が得られる。

Pythonでerfcを計算し（scipy.special.erfc）、FEA結果と自動比較するスクリプトを書く。

まず十分細かい時間刻みで空間収束を確認し、次に収束したメッシュで時間刻みを変えるという2段階法で5ケースに減らせる。ただし両方向が相互に影響する場合（クーラン数依存の安定化スキームなど）は独立した評価が必要になる。

Newtonの冷却法則 $q = h_{conv}(T_s - T_\infty)$ を境界条件にすると、理論解がより複雑になるが閉形式は存在する。

20要素（バイアスメッシュ）、$\Delta t = 0.32$ s（100ステップ）での結果だ。

laplacianFoamのデフォルト時間スキームはEuler（1次精度前進）だから、同じ$\Delta t$では2次精度の他ソルバーより誤差が大きい。ddtSchemes を backward（BDF2、2次精度）に変えれば精度は向上する。CFDソルバーの時間スキームのデフォルトがFEAソルバーと異なることを認識しておくべきだ。

Abaqusで同一メッシュ、$\Delta t = 1.6$ s（粗い刻み）での比較だ。

2次精度だから$\Delta t$が大きいときの精度は高い。ただしステップ変化の初期に振動（オーバーシュート）が出る欠点がある。実務では最初の数ステップを後退Eulerで計算し、その後Crank-Nicolsonに切り替える手法が使われる。Abaqusではこの切り替えが自動で行われる。

直交座標で各方向が独立なら、変数分離法で1D解の積として2D/3D解を構成できる。例えば矩形領域の4面をステップ加熱する場合、各方向のerfc解の積になる。円筒座標や球座標ではBessel関数やルジャンドル関数を含む解になる。

Stefan問題（移動境界問題）だ。固液界面が時間とともに移動する。1次元の場合、Neumann解として厳密解が存在する。界面位置は $s(t) = 2\lambda\sqrt{\alpha t}$ で進行し、$\lambda$ は超越方程式の根から決まる。

潜熱を与える温度幅（マッシーゾーン幅）の設定が結果に大きく影響する。理論解は幅ゼロ（ステップ関数）だが、数値的には有限幅が必要。幅を系統的に狭くして結果が収束するか確認する。また、時間刻みが粗いと界面が1ステップで複数要素を通過して精度が劣化する。

非線形問題になり、一般には厳密解が存在しない。しかしKirchhoff変換 $\Theta = \int_0^T k(T')dT'$ を使えば、$k(T)$が既知関数の場合に線形化できるケースがある。

MMS（製造解法）が有効だ。任意の温度場 $T(x,t)$ を仮定し、非線形熱伝導方程式に代入してソース項を逆算する。このソース項をFEAに与えて得られた数値解と仮定した解を比較する。非線形問題のCode Verificationの標準的手法だ。

非定常熱伝導で得た温度場を構造解析の熱荷重として使う。自由膨張する棒の場合、温度分布$T(x)$に対して歪みは$\varepsilon_{th} = \alpha_{th}(T-T_{ref})$、応力はゼロ（拘束がない場合）。

メッシュマッピング（熱メッシュから構造メッシュへの温度転写）の精度だ。メッシュが同一なら問題ないが、異なるメッシュ間で補間する場合、補間誤差が温度場に乗る。最悪の場合、局所的な温度ジャンプが偽の熱応力を生成する。同一メッシュでの結果を基準にして、異なるメッシュ間のマッピング誤差を定量評価すべきだ。

これはCrank-Nicolsonの既知の問題で、ステップ入力に対するGibbs現象的なオーバーシュートだ。物理的に温度が入力値を超えることはないから、非物理的な結果になる。

Abaqusは過渡熱解析で自動的に最初のいくつかのインクリメントを後退Eulerで処理する。ユーザーが意識しなくてもオーバーシュートが抑制される。ただしマニュアルの記述が不明確なので、実際の挙動をテスト問題で確認しておくべきだ。

断熱端からの反射で温度が過大評価される。温度浸透前線がモデル端に到達すると、半無限体の仮定が崩れる。

モデルを延長して再計算する。温度変化が$10^{-3}$°C以下になるまでモデル長さを増やす。実務的にはerfcの漸近挙動から必要長さを事前に推定できる。erfc(3) ≈ 0.00002 だから、$x = 3 \times 2\sqrt{\alpha t}$ で温度変化は初期変化の0.002%になる。

AbaqusとCOMSOLは時間刻みの自動適応制御を備えている。Abaqusでは*HEAT TRANSFER, DELTMX=5 で「1ステップ内の最大温度変化を5°C以内」に制限し、それに合わせて$\Delta t$が自動調整される。

adjustTimeStep yes; maxCo 0.5; と設定すればクーラン数ベースの自動調整が効く。ただしlaplacianFoamでは対流項がないからCo数の概念が直接適用されない。maxDeltaT で上限を設定し、拡散数 $Fo = \alpha\Delta t/h^2 < 0.5$ を手動で確認するのが安全だ。