多次元非定常熱伝導
理論と物理
概要 — なぜ多次元が必要か
先生、多次元の非定常熱伝導ってFEMでしか解けないんですか？ 1次元のHeislerチャートは教科書で見たんですけど、実際の部品って3Dですよね…。
いい質問だ。実は、直交座標系で各方向の境界条件が独立なら「積の解（product solution）」が使えるんだよ。たとえば短い角柱（レンガのような形状）の温度は
のように、1次元解を各方向について掛け算するだけで求まるんだ。FEMがなくても手計算で3Dの非定常温度場が得られる。
え、掛け算だけで3Dの温度分布がわかるんですか!? それってどういう仕組みなんですか？
ポイントは変数分離法だ。多次元の熱伝導方程式は、境界条件が各軸方向で分離可能なら、各方向ごとの1D問題に分解できる。数学的には偏微分方程式が各変数の常微分方程式の積に帰着するということだ。たとえば自動車の鍛造品（短い円柱形状）の焼入れ冷却を考えてみよう。側面は円柱として、上下面は平板として扱って、それぞれの1D解を掛ければいい。
なるほど！ ということは、FEMの結果を検証するときにも使えるってことですか？
まさにその通り。積の解はFEMのベンチマーク（検証解）として非常に有用なんだ。メッシュや時間刻みの妥当性を、厳密な解析解と比較して定量的に評価できる。実務でこれを知っているかどうかで、解析結果への信頼度がまるで違うよ。
支配方程式
ではまず、多次元の非定常熱伝導の基本方程式を教えてください。
出発点はフーリエの熱伝導方程式の多次元版だ。等方性・均質な媒体で内部発熱がない場合：
ここで $\alpha = k/(\rho c_p)$ は熱拡散率 [m²/s]。この式は放物型偏微分方程式で、初期条件 $T(\mathbf{x}, 0) = T_i$ と境界条件を与えれば一意に解が定まる。
円柱座標系ではどうなるんですか？ エンジンシリンダーとか円筒形の部品が多いですよね。
円柱座標系 $(r, \phi, z)$ では：
軸対称（$\phi$ 方向に均一）で有限長の場合、$r$ 方向と $z$ 方向の2次元問題になる。これがまさに積の解の出番だ。
無次元温度 $\theta^*$、ビオ数 $\text{Bi}$、フーリエ数 $\text{Fo}$ の定義を整理しておこう：
ここで $h$ は対流熱伝達係数、$L_c$ は特性長さ（平板なら半厚、円柱なら半径）、$T_\infty$ は周囲流体温度、$T_i$ は初期温度である。
積の解（Product Solution）
積の解の数学的な根拠をもう少し詳しく教えてください。なぜ掛け算でいいんですか？
核心は変数分離の成立条件にある。解を $\theta^*(\mathbf{x},t) = X(x,t) \cdot Y(y,t) \cdot Z(z,t)$ の形で仮定して支配方程式に代入すると：
各項がそれぞれの変数のみに依存していれば、各方向の1D問題に分離される。こうして得られた $X$, $Y$, $Z$ がそれぞれ1D非定常熱伝導の解（Heisler解）になるんだ。
具体的な幾何形状ではどうなりますか？
代表的な3つのケースを挙げよう：
1. 短い角柱（Short Rectangular Bar） — 3方向の平板解の積
2. 短い円柱（Short Cylinder） — 円柱解 × 平板解
3. 半無限固体の角（Semi-infinite Corner） — 3方向の半無限固体解の積
たとえば電子部品の樹脂パッケージ（直方体）のリフロー加熱を考えると、ケース1がそのまま使えるよ。
$x$ 方向に50%冷えて、$y$ 方向にも独立に50%冷えたとする。全体としては $0.5 \times 0.5 = 0.25$、つまり75%冷えたことになる。各方向の冷え具合は独立に進行し、それぞれが掛け算で組み合わさるのが積の解の本質だ。
積の解が使えない場合ってどんなときですか？
重要なポイントだ。以下の場合は積の解が適用できない：
- 内部発熱がある場合（$\dot{q} \neq 0$）— 方程式が非斉次になり変数分離不可
- 物性値が温度依存する場合 — 方程式が非線形になる
- 境界条件が方向間で結合する場合（斜めの境界、接触抵抗の空間分布など）
- 座標系が直交しない形状（L字型、T字型など）
これらのケースでは数値解法（FEM、FDM、FVM）が必要になる。ただし、多くの実務的な検証問題は積の解で扱える形状に設定されるので、この手法を知っておく価値は大きいよ。
重ね合わせ原理（Superposition）
積の解は「掛け算」でしたけど、「足し算」で解を組み合わせることもできるんですか？
そう、重ね合わせ原理（superposition principle）だ。線形の偏微分方程式では、個別の解を足し合わせて新しい解を作ることができる。これは積の解とは別の強力な手法だよ。
例えば、ある面だけ急加熱して他面は断熱という非対称な境界条件でも、対称問題と反対称問題に分解して足し合わせれば解ける。鋼材の片面焼入れなどがまさにこのパターンだ。
積の解と重ね合わせは何が違うんですか？ 混同しそうです…。
整理するとこうなる：
| 手法 | 演算 | 使い分け |
|---|---|---|
| 積の解 | 掛け算（$\theta^* = \theta^*_1 \cdot \theta^*_2$） | 異なる空間方向の1D解を組み合わせて多次元解を構成 |
| 重ね合わせ | 足し算（$T = T_1 + T_2 - T_{\text{ref}}$） | 同じ空間で異なる境界条件の解を組み合わせる |
積の解は「空間方向の分解」、重ね合わせは「境界条件の分解」と覚えるといい。両者を組み合わせることで、かなり複雑な問題も解析的に扱えるようになるんだ。
形状係数法（Shape Factor Method）
教科書で「形状係数」っていう概念も出てきたんですけど、これは非定常とどう関係するんですか？
形状係数法は主に定常状態の多次元熱伝導で使う手法だ。導電形状係数（conduction shape factor）$S$ を定義して：
ここで $Q$ は単位時間あたりの熱流量 [W]。たとえば、地中に埋設された配管からの放熱、半導体チップから基板への熱流路など、解析解が得られている標準形状のSの値を表から引いてくるだけで熱流量を計算できる。
具体的にはどんな形状のSが表に載っているんですか？
代表的なものを挙げると：
| 形状 | 形状係数 $S$ | 条件 |
|---|---|---|
| 同心円柱（長さ $L$） | $\dfrac{2\pi L}{\ln(r_2/r_1)}$ | $L \gg r_2$ |
| 同心球殻 | $\dfrac{4\pi r_1 r_2}{r_2 - r_1}$ | — |
| 地中埋設パイプ | $\dfrac{2\pi L}{\cosh^{-1}(z/r)}$ | $L \gg z$, $z > r$ |
| 立方体の角 | $0.15 L$ | 3面が等温面 |
Incroperaの伝熱工学の教科書に数十パターンの形状係数が載っている。実務では概算段階でこれを使い、詳細設計でFEMに移行するという流れが多いよ。
形状係数は熱抵抗の概念と関係があるんですか？
まさにそうだ。多次元の導電熱抵抗は $R_{\text{cond}} = 1/(kS)$ と表せる。これを対流抵抗 $R_{\text{conv}} = 1/(hA)$ と直列・並列に組み合わせれば、回路アナロジーで複雑な系の総熱抵抗を手計算で評価できる。電子機器の放熱設計では日常的に使う手法だよ。
数値解法と実装
2D/3D Heislerチャートの使い方
Heislerチャートを多次元に拡張するには、具体的にどういう手順を踏めばいいですか？
手順を短い円柱（半径 $r_0$、半長 $L$）の例で説明しよう。全面が対流冷却（$h$, $T_\infty$）される場合：
- ステップ1： $r$ 方向について、Bi$_{\text{cyl}} = hr_0/k$、Fo$_{\text{cyl}} = \alpha t/r_0^2$ を計算し、無限長円柱のHeislerチャートから中心温度の $\theta^*_{\text{cyl}}(0,t)$ を読む
- ステップ2： $z$ 方向について、Bi$_{\text{wall}} = hL/k$、Fo$_{\text{wall}} = \alpha t/L^2$ を計算し、無限平板のHeislerチャートから中央面温度の $\theta^*_{\text{wall}}(0,t)$ を読む
- ステップ3： 短い円柱の中心温度は積で求まる
表面温度や任意の位置の温度が必要な場合は、Heislerの第2チャート（位置補正チャート）も使って場所ごとの補正をかけるんだ。
数値例で確認させてください！ たとえばステンレス鋼の短い円柱を急冷する場合は？
具体的にやってみよう。SUS304短円柱（$r_0 = 25$ mm、$L = 50$ mm = 半長）：
- $k = 14.9$ W/(m-K)、$\alpha = 3.95 \times 10^{-6}$ m²/s
- $h = 200$ W/(m²-K)、$T_i = 500$°C、$T_\infty = 25$°C
- 求めたい時刻：$t = 300$ s
$r$ 方向： Bi$_{\text{cyl}} = 200 \times 0.025 / 14.9 = 0.336$、Fo$_{\text{cyl}} = 3.95 \times 10^{-6} \times 300 / 0.025^2 = 1.896$
チャートから $\theta^*_{\text{cyl}}(0,t) \approx 0.23$
$z$ 方向： Bi$_{\text{wall}} = 200 \times 0.05 / 14.9 = 0.671$、Fo$_{\text{wall}} = 3.95 \times 10^{-6} \times 300 / 0.05^2 = 0.474$
チャートから $\theta^*_{\text{wall}}(0,t) \approx 0.61$
中心温度： $\theta^* = 0.23 \times 0.61 = 0.140$
$T_{\text{center}} = T_\infty + \theta^*(T_i - T_\infty) = 25 + 0.140 \times 475 = 91.6$°C
300秒後の中心温度は約92°Cだ。FEMで同条件を解いてこの値と比較すれば、メッシュ・時間刻みの妥当性を確認できるよ。
FEMによる多次元離散化
積の解が使えない場合はFEMに頼ることになりますよね。多次元の非定常熱伝導をFEMで解く場合のポイントは何ですか？
FEMの定式化のポイントは空間離散化と時間離散化の分離だ。弱形式（ガラーキン法）を空間に適用すると、半離散化された常微分方程式系になる：
ここで $[\mathbf{C}]$ は容量マトリクス（質量マトリクスの熱版）、$[\mathbf{K}]$ は伝導マトリクス（剛性マトリクスの熱版）、$\{F\}$ は外部熱負荷ベクトルだ。
これに対して時間方向にはオイラー法（前進/後退）やクランク-ニコルソン法などの時間積分スキームを適用する。2Dなら三角形要素や四角形要素、3Dならテトラヘドロンやヘキサヘドロン要素が使われる。
時間積分スキームの選択
時間積分のスキームはどう選べばいいんですか？ 前進オイラーとか後退オイラーとか色々ありますよね。
一般化された $\theta$ 法（ここでの $\theta$ はパラメータで、無次元温度とは別物だ）で統一的に書ける：
| $\theta$ 値 | スキーム名 | 特徴 |
|---|---|---|
| $0$ | 前進オイラー（陽解法） | 条件付き安定。$\Delta t \le \Delta x^2 / (2\alpha)$ の制約あり |
| $1/2$ | クランク-ニコルソン | 無条件安定・2次精度。ただし振動解のリスクあり |
| $2/3$ | ガラーキン法 | 無条件安定・振動抑制。商用コードで多用 |
| $1$ | 後退オイラー（陰解法） | 無条件安定・1次精度。数値拡散が大きい |
実務では$\theta = 2/3$（ガラーキン法）が最もバランスが良い。Ansysのデフォルトもこれだよ。
多次元メッシュ戦略
多次元の非定常問題では、メッシュの切り方にコツがあるんですか？
大事なポイントが3つある：
- 温度勾配が急な領域にメッシュを密にする — 対流面の近傍、角部、材料境界。特に初期のFo < 0.2 では表面近傍に急な勾配が生じるので、最低3〜5層は細かくしたい
- アスペクト比を抑える — 熱伝導でも1:10以上の扁平要素は精度低下を招く。特に角部では等方的な要素を心がける
- 対称性を活用する — 積の解が適用できる形状なら、1/4モデルや1/8モデルで計算コストを大幅削減できる。対称面には断熱条件（$\partial T/\partial n = 0$）を課す
ペクレ数 $\text{Pe} = uL/\alpha$ が大きい移流-拡散問題とは違い、純粋な熱伝導では数値的な不安定性は比較的穏やかだ。むしろ時間刻みとメッシュサイズの整合性（フーリエメッシュ数 $\text{Fo}_{\text{mesh}} = \alpha \Delta t / \Delta x^2$）が重要だよ。
実践ガイド
ベンチマーク問題 — 積の解 vs FEM
実際にFEMの検証用ベンチマークとして積の解を使うにはどうすればいいですか？
典型的なベンチマーク手順はこうだ：
- 問題設定： 短い角柱（$2a \times 2b \times 2c$）、等方性材料、全面対流冷却。初期温度 $T_i$、周囲温度 $T_\infty$
- 解析解の計算： 各方向の1D Heisler解を級数展開の最初の1項近似（Fo > 0.2 の場合）で計算し、その積を取る
ここで $\zeta_1$ は $\zeta_n \tan \zeta_n = \text{Bi}$ の第1根、$C_1 = 4\sin\zeta_1 / (2\zeta_1 + \sin 2\zeta_1)$
- FEM解との比較： 中心温度、表面温度、角部温度をそれぞれ比較。許容誤差は通常 1〜2% 以内
- メッシュ収束性の確認： メッシュを2倍に細かくして誤差が減少することを確認
角部の温度は中心より冷えるのが速いんですか？
その通り！ 角部は3方向から同時に冷却されるので最も速く温度が下がる。積の解で考えると、角部では各方向の $\theta^*$ が全て表面値（最小値）なので、その積はさらに小さくなる。逆に中心部は全方向で最も冷えにくい位置だ。焼入れで角部に割れが入りやすいのは、この温度差による熱応力が原因なんだよ。
典型的な解析ワークフロー
多次元非定常熱伝導の解析を実務でやるとき、全体の流れはどうなりますか？
基本的な流れはこうだ：
- 問題の特性把握 — まずビオ数を計算。Bi < 0.1 なら集中熱容量法で十分。Bi > 0.1 なら空間分布を考慮した解析が必要
- 形状の対称性チェック — 積の解が適用できるか判断。適用可能なら手計算で概算値を出す
- FEMモデル構築 — 対称性を活用した1/2、1/4モデル。温度勾配が急な面の近傍にメッシュを密にする
- 時間刻みの設定 — 初期の急変期（Fo < 0.2）は細かく、定常に近づいたら粗くするアダプティブ制御が理想
- 結果の検証 — 積の解やNAFEMSベンチマーク（T2, T3）との比較。エネルギー収支の確認
境界条件の設定指針
多次元の場合、境界条件の設定で気をつけることはありますか？
多次元特有の注意点がいくつかある：
- 角部・稜線の処理： 2面または3面が同時に対流冷却される角部では、熱伝達面積の二重カウントに注意。FEMソフトによっては自動で面積を適切に配分するが、手動設定の場合は面ごとに分離して設定すること
- 接触熱抵抗： 2Dや3Dでは部品同士の接触面が現れる。接触熱伝達係数 $h_c$ は圧力や表面粗さに依存し、典型的には 1,000〜50,000 W/(m²-K) の範囲
- 輻射境界条件： 高温（>300°C）では輻射の効果が無視できなくなる。輻射は $T^4$ に比例するため非線形になり、積の解は適用不可になる
- 対称面の断熱条件： $\partial T / \partial n = 0$ を忘れると計算量が無駄に増える。1/8モデルで済むところを全体モデルで計算するのは資源の無駄だ
商用ツール別の実装手順
Ansysで多次元非定常解析をやる場合、どんな設定をすればいいですか？
主要ソルバーでの設定ポイントを比較しよう：
| 項目 | Ansys Mechanical | Abaqus | COMSOL |
|---|---|---|---|
| 解析タイプ | Transient Thermal | *HEAT TRANSFER, TYPE=TRANSIENT | Heat Transfer in Solids (Time Dependent) |
| 要素 | SOLID70 (8node), SOLID90 (20node) | DC3D8, DC3D20 | 自動選択（2次がデフォルト） |
| 時間積分 | $\theta = 2/3$（ガラーキン法、デフォルト） | 後退オイラー（デフォルト） | BDF（後退差分公式） |
| 自動時間刻み | AUTOTS, ON | *CONTROLS, PARAMETERS=TIME INCREMENTATION | Time stepping: Free |
| 出力 | NSOL (節点温度), ESOL (熱流束) | .odb (Field Output) | ポストプロセスで自由に定義 |
ソフトウェア比較
多次元非定常熱伝導の対応状況
多次元の非定常熱伝導を扱えるソルバーにはどんなものがありますか？ 使い分けの基準も知りたいです。
主要なソルバーの特徴を整理しよう：
| ソルバー | 開発元 | 強み | 多次元非定常熱への適性 |
|---|---|---|---|
| Ansys Mechanical | Ansys Inc. | 構造-熱連成、APDLスクリプト | 非常に高い。熱応力連成に最適 |
| Abaqus Standard | Dassault Systemes | 非線形収束性、Subroutine拡張 | 高い。UMATNETで物性の温度依存性を自由に定義 |
| COMSOL Multiphysics | COMSOL AB | マルチフィジクス統合、GUI | 非常に高い。電気-熱-構造の3方向連成が容易 |
| OpenFOAM (laplacianFoam) | OSS | カスタマイズ性、コストゼロ | 中程度。FVM系のため基本的な熱伝導向き |
| Simcenter STAR-CCM+ | Siemens | CHT (共役熱伝達) | 高い。流体-固体の連成非定常に強い |
選定の指針
結局どれを使えばいいんですか？ 迷います…。
ユースケースで分けるのが一番だ：
- 熱応力まで求めたい（焼入れ変形、はんだ熱疲労）→ Ansys or Abaqus
- 流体との連成が主体（電子機器の自然対流冷却）→ STAR-CCM+ or COMSOL
- 電磁場-熱の連成（誘導加熱、マイクロ波加熱）→ COMSOL
- 教育・研究でコストをかけたくない → OpenFOAM + Python（積の解で検証）
- 概算で十分（設計初期段階）→ 形状係数法 + Excelスプレッドシート
ポイントは「多次元非定常の解析自体はどのソルバーでもできる」ということ。差がつくのは前後処理の効率性と他の物理場との連成のしやすさだよ。
先端技術
先端トピックと研究動向
多次元非定常熱伝導の分野で、最近の先端的な研究ってどんなものがありますか？
いくつか注目すべき動向がある：
- AIサロゲートモデル： 物理インフォームドニューラルネットワーク（PINN）で多次元非定常の温度場を学習し、リアルタイム予測するアプローチ。DeepONetやFourier Neural Operatorが注目されている。FEMの数千回のモンテカルロシミュレーションの代替として期待大
- アダプティブメッシュリファインメント（AMR）： 温度勾配が時間とともに移動する問題で、メッシュを動的に細分化・粗大化する手法。レーザー加熱やアーク溶接のシミュレーションで威力を発揮
- 非フーリエ熱伝導： 超短パルスレーザー加熱（フェムト秒スケール）では、フーリエの法則が破綻し、Cattaneo-Vernotte方程式（双曲型）が必要になる
$\tau$ は熱緩和時間（金属で $\sim 10^{-11}$ s）。通常のマクロスケール問題では無視できるが、ナノスケールやフェムト秒パルスでは本質的に重要になる。
- トポロジー最適化との統合： 多次元の温度場を目的関数として、最適な放熱経路（ヒートシンク形状）を自動設計する手法が実用化されつつある
- GPU並列計算： 多次元非定常問題は大規模行列の反復解法が支配的。CUDAベースのGPUソルバーにより、従来の10〜100倍の高速化が報告されている
PINNで温度場を予測するのはすごいですね。でも積の解みたいな古典的な解法はもう不要になるんですか？
逆だよ。AIモデルの学習データの正解値として、積の解やHeislerチャートの厳密解が不可欠なんだ。学習データの品質がAIの精度を決めるから、古典的な解析解の価値はむしろ高まっている。「新しい手法を検証するためには、古い手法が必要」という構造は、科学の本質だと思うよ。
トラブルシューティング
よくあるエラーと対策
多次元非定常の熱解析で、よくある失敗パターンを教えてください。
実務でよく遭遇するトラブルを紹介しよう：
1. 温度が振動する（オシレーション）
- 原因： 時間刻みが粗すぎる or クランク-ニコルソン法の特性
- 対策： $\text{Fo}_{\text{mesh}} = \alpha \Delta t / \Delta x^2 < 1/2$ を目安に時間刻みを細かくする。または $\theta = 2/3$ のガラーキン法に切り替え
2. 角部の温度が非物理的（$T_\infty$ を下回る）
- 原因： メッシュが粗すぎて数値拡散が過大
- 対策： 角部周辺のメッシュを2倍以上に細かくする。積の解と比較して確認
3. エネルギー保存が崩れる
- 原因： 容量マトリクスの集中化（lumped）と分散（consistent）の選択ミス
- 対策： エネルギー収支を必ずチェック。$\int_V \rho c_p (T_{\text{final}} - T_i) dV = \int_0^t Q_{\text{boundary}} dt$
4. 収束しない（非線形問題）
- 原因： 輻射境界条件の $T^4$ 項や温度依存物性値
- 対策： 初期時間刻みを十分小さく（$\Delta t_{\text{init}} \le 0.01 \times L^2/\alpha$）。緩和係数を0.5〜0.7程度に設定
「解析が合わない」と思ったら
FEMの結果が積の解と全然合わないとき、どこから疑えばいいですか？
デバッグの鉄則を教えよう：
- まず1Dに落とす — 多次元問題をいきなりデバッグするのは困難。まず1方向だけの1D問題として、Heislerチャートの値と合うか確認する
- 単位系を確認 — これが一番多い原因。特に $\alpha$ の単位 [m²/s] に対し、長さを [mm] で入力していないか。 $1 \text{mm}^2 = 10^{-6} \text{m}^2$ なのでオーダーが6桁ずれる
- 特性長さの定義を確認 — 平板のHeisler解では半厚（$L$）を使う。全厚（$2L$）を代入するとBiとFoが4倍ずれる
- 積の解の適用条件を確認 — 内部発熱項、温度依存物性、非直交形状のいずれかが入っていないか
- メッシュ収束性を確認 — 要素数を2倍にして結果が変わらなければメッシュはOK。変わるなら不十分
経験則として、解析が合わない原因の80%は入力ミスだ。物理の問題ではなく、人間のタイプミスや単位の間違い。だからこそ積の解という「答え合わせの手段」を持っていることが重要なんだよ。
その理解が大事だ。多次元非定常熱伝導は、解析解（積の解、形状係数法）と数値解（FEM）の両方を使いこなしてこそプロのCAEエンジニアと言える。次回は相変化（PCM）を伴う非定常問題を扱おう。融解・凝固のシュテファン問題だ。
関連トピック
なった
詳しく
報告