科学

熱伝導方程式の解き方は?変数分離法と数値解法の基本を解説!(解析解・差分法・陽解法・陰解法)

熱伝導方程式の解法の全体像
当サイトでは記事内に広告を含みます
いつも記事を読んでいただきありがとうございます!!! これからもお役に立てる各情報を発信していきますので、今後ともよろしくお願いします(^^)/

熱伝導方程式は、物体の内部で温度がどのように変化するかを表す代表的な偏微分方程式です。

棒を加熱したときの温度分布、金属部品の冷却、建築物の断熱設計、電子機器の放熱解析など、多くの場面で利用されます。

解き方には、条件を満たす式を導く解析解と、計算機を使って温度を段階的に求める数値解法があります。

本記事では、変数分離法による解析解と差分法による数値計算を中心に、陽解法と陰解法の違い、計算時に注意したい安定条件まで順に解説します。

熱伝導方程式の解法の全体像

熱伝導方程式の解法の全体像

それではまず、熱伝導方程式の解き方と選び方について解説していきます。

熱伝導方程式の基本式

一次元の熱伝導方程式は、棒の長さ方向だけに温度が変化すると仮定した場合に用いられます。

温度を u、位置を x、時間を t、熱拡散率を αとすると、基本式は次の形です。

∂u/∂t = α∂²u/∂x²

u は位置と時間によって変化する温度です。

α は熱が材料内へ広がる速さを表す熱拡散率です。

左辺は時間に対する温度変化を表し、右辺は温度分布の曲がり具合を表します。

高温部と低温部の差が大きいほど熱は移動しやすく、時間が経つにつれて温度差は小さくなります。

この現象を数式化したものが熱伝導方程式です。

熱拡散率は、熱伝導率、密度、比熱から求められます。

材料によって値が異なるため、銅、アルミニウム、鉄、樹脂、コンクリートでは温度変化の速さも変わります。

解析解と数値解法の違い

解析解とは、方程式を数学的に変形し、温度を式として直接表す方法です。

条件が比較的単純であれば、変数分離法やフーリエ級数を使って温度分布の一般的な形を得られます。

一方で、複雑な形状、材料の違い、局所的な発熱、時間ごとに変わる境界条件を扱う場合には、解析解だけで対応するのは難しくなります。

そのような場合に使われるのが数値解法です。

数値解法では、対象を細かい格子に区切り、各時刻の温度を繰り返し計算します。

代表例が有限差分法であり、計算式の作り方によって陽解法と陰解法に分かれます。

分類 特徴 向いている場面
解析解 式の形で温度分布を表せる 単純な棒、平板、円柱など
変数分離法 位置と時間を分けて解く 境界条件が整った基本問題
差分法 格子点ごとに温度を計算する 実務的な温度解析
陽解法 計算が比較的簡単 時間刻みを小さくできる場合
陰解法 安定性に優れる 長時間の計算や粗い時間刻み

初期条件と境界条件の重要性

熱伝導方程式を解くには、方程式だけでは足りません。

開始時点の温度を示す初期条件と、端部や表面の状態を示す境界条件が必要です。

初期条件の例としては、棒全体が一様に 20℃である、中央部だけが高温である、といった設定があります。

境界条件には、端の温度を一定に保つ条件、熱の出入りを止める断熱条件、外気と熱を交換する条件などがあります。

熱伝導問題では、式の解き方よりも先に初期条件と境界条件を明確にすることが重要です。

条件が曖昧なままでは、計算できても現実の挙動と一致しにくくなります。

端部を常に 0℃に保つ場合は、温度固定の境界条件として扱います。

端部から熱が逃げない場合は、温度の傾きが 0になる断熱条件を設定します。

実際の設計では、周囲の空気への放熱も考える必要があるため、熱伝達係数を含む境界条件が用いられることもあります。

変数分離法による解析解

続いては、変数分離法を使った解析解の考え方を確認していきます。

変数分離の考え方

変数分離法では、温度 u を位置だけの関数 X と時間だけの関数 T の積として考えます。

つまり、u を X と T の積の形に置き換えます。

u(x,t)= X(x)T(t)

位置による変化を X、時間による変化を T として分けて扱います。

この形を熱伝導方程式へ代入すると、位置に関する式と時間に関する式を分けることができます。

偏微分方程式を二つの常微分方程式に分けられる点が、変数分離法の大きな利点です。

ただし、どの問題にもそのまま使えるわけではありません。

境界条件が整理されており、形状や物性値が単純な場合に特に有効な方法です。

両端温度固定の棒の問題

長さ L の棒があり、両端の温度を 0℃に固定する場合を考えます。

このとき、位置に関する関数は正弦関数で表されます。

温度分布は、複数の正弦波を重ね合わせたフーリエ級数の形になります。

u(x,t)= Σ An sin(nπx/L)exp(-αn²π²t/L²)

An は初期温度分布から決まる係数です。

n は 1、2、3 と続く整数です。

exp の部分は時間の経過に伴う温度変化を表します。

時間が大きくなるほど指数関数の値は小さくなり、棒の内部にあった温度差は徐々に消えていきます。

高次の成分ほど早く減衰するため、細かな温度の凸凹は先に滑らかになります。

この性質は、熱が温度分布をならしていく現象とよく対応しています。

フーリエ係数の求め方

解析解に含まれる An は、初期温度分布に応じて決まります。

初期温度を f(x)とすると、フーリエ係数は正弦関数との積分によって求めます。

初期温度分布が単純な直線や一定値であれば、計算は比較的進めやすいでしょう。

一方で、温度分布が区分的であったり、局所的な加熱が含まれたりすると、係数の計算は複雑になります。

解析解は温度変化の仕組みを理解するうえで非常に有用ですが、係数を多く求める必要がある問題では数値計算のほうが実務的なこともあります。

変数分離法では、初期条件を満たすために複数の固有関数を重ね合わせます。

一つの正弦関数だけで表せない温度分布でも、級数として近似できる点が重要です。

解析解を使う際は、係数の項数を増やすほど初期状態の再現性が上がる一方、計算量も増えることを意識してください。

また、境界条件と固有関数の組み合わせが合っているかを確認することも欠かせません。

有限差分法による数値解法

続いては、差分法によって熱伝導方程式を計算する基本を確認していきます。

格子分割と差分近似

有限差分法では、連続的な棒や平板を小さな区間に分け、格子点ごとの温度を扱います。

位置方向の刻み幅を Δx、時間方向の刻み幅を Δt とします。

温度の時間微分は、ある時刻と次の時刻の温度差で近似できます。

空間に関する二階微分は、中央の格子点と左右の格子点の温度差から近似します。

この置き換えにより、偏微分方程式は加減乗除を中心とした計算式になります。

手計算では少数の格子点に限られますが、表計算ソフトやプログラムを使えば多数の格子点でも処理できます。

記号 意味 設定時の注意点
Δx 位置方向の格子幅 小さいほど空間分布を細かく表せる
Δt 時間方向の刻み幅 大きすぎると計算が不安定になりやすい
u i n 位置 i、時刻 n の温度 添字の意味を統一して管理する
α 熱拡散率 材料と温度条件に応じた値を使う
r αΔt/Δx² 陽解法の安定性を左右する値

差分式の導出

一次元問題で時間を前進差分、空間を中心差分で近似すると、次の形になります。

u i n+1 = u i n+r(u i+1 n-2u i n+u i-1 n)

r = αΔt/Δx²

次の時刻の温度を、現在の自分と左右の温度から求める式です。

右側の括弧内は、中央点と隣接点の温度差をまとめたものです。

中央点が周囲より高温なら値は負になり、次の時刻では温度が下がる方向に働きます。

反対に、中央点が周囲より低温なら値は正になり、周囲から熱を受け取って温度が上がります。

この繰り返しによって、熱が高温側から低温側へ広がる様子を表現できます。

差分式は物理現象を細かな時間ステップへ分解したものと考えると理解しやすくなります。

計算手順の整理

差分法を使う場合は、計算の順番を固定しておくとミスを減らせます。

まず、解析する長さ、格子幅、時間刻み、計算終了時間を決めます。

次に、すべての格子点へ初期温度を設定し、端部の境界条件を与えます。

その後、内部の格子点について差分式を使い、一つ先の時刻の温度を計算します。

新しい温度を全格子点分そろえたら、次の時刻へ進みます。

同じ時刻の温度だけを使って次時刻の温度を求めることが、陽解法の基本です。

計算途中で更新済みの温度を混ぜると、意図しない計算結果になるおそれがあります。

表計算ソフトで試す場合は、時刻を行、位置を列に割り当てると温度変化を追いやすくなります。

グラフを作成し、各時刻の温度分布を折れ線で確認すると、異常な振動や不自然な発散にも気付きやすいでしょう。

陽解法の安定条件

続いては、陽解法で特に重要となる安定条件を確認していきます。

陽解法の計算構造

陽解法は、現在時刻の温度が分かれば、次の時刻の温度を直接計算できる方法です。

連立方程式を解く必要がないため、考え方が分かりやすく、初学者が差分法を学ぶ際にもよく使われます。

格子点ごとの計算式も単純であり、プログラム化しやすい点が魅力です。

ただし、時間刻みを自由に大きくできるわけではありません。

条件を守らずに計算すると、実際には起こらない温度の振動や極端な増減が現れます。

これを数値的不安定と呼びます。

クーラン条件に近い安定性の目安

一次元の熱伝導方程式を陽解法で計算する場合、一般的には r が 2分の1以下となる必要があります。

つまり、熱拡散率、時間刻み、格子幅の関係を確認しなければなりません。

αΔt/Δx² ≦ 1/2

Δt ≦ Δx²/2α

格子幅を半分にすると、許容される時間刻みはおよそ4分の1になります。

この関係は、格子を細かくするほど計算回数が大幅に増えることを意味します。

精度を高めるために Δx を小さくした結果、計算時間が増える点は、数値解析における重要なトレードオフです。

二次元や三次元ではさらに条件が厳しくなるため、設定前に安定条件を必ず確認してください。

発散を防ぐ設定方法

計算結果が途中から大きく正負へ振れる場合、まず時間刻みが大きすぎないかを確認します。

熱拡散率の単位も見直す必要があります。

熱伝導率をそのまま熱拡散率として使ってしまうと、まったく異なる結果になるため注意が必要です。

熱拡散率の単位は、通常 m²/s で扱います。

長さを mm で入力し、熱拡散率を m²/s のまま使うと単位が一致しません。

このような単位の不統一も、計算結果がおかしくなる原因です。

陽解法では、安定条件、単位系、境界条件の三つを計算開始前に確認してください。

温度が発散した場合は、まず Δt を小さくすることが基本の対処になります。

安定条件を満たしていても、格子が粗すぎれば急な温度変化を十分に表現できません。

安定性と精度は別の観点であり、両方を確認する姿勢が大切です。

陰解法と陽解法の比較

続いては、陰解法の特徴と陽解法との違いを確認していきます。

陰解法の基本構造

陰解法では、次の時刻の温度を含む式を立て、それらをまとめて解きます。

陽解法のように一つずつ直接計算するのではなく、全格子点の未知温度を連立方程式として扱う方法です。

代表的な手法には、後退差分を使う完全陰解法があります。

一回の時間ステップごとに行列計算が必要になるため、計算式は陽解法より複雑になります。

しかし、時間刻みに対する安定性が高く、長い時間を扱う計算では有利です。

計算の簡単さを優先するなら陽解法、安定性を重視するなら陰解法という考え方が基本になります。

計算負荷と精度の違い

陽解法は一回の計算が軽いものの、安定条件のために非常に細かな時間刻みが必要になる場合があります。

陰解法は一回ごとに連立方程式を解く必要がありますが、より大きな時間刻みを設定できることがあります。

どちらが速いかは、格子数、解析時間、必要な精度、使用する計算環境によって変わります。

比較項目 陽解法 陰解法
次時刻の計算 直接求めやすい 連立方程式を解く必要がある
実装の分かりやすさ 比較的高い やや複雑
安定条件 時間刻みに制約がある 一般に安定性が高い
一回の計算負荷 小さい 大きい
長時間計算 条件によって負担が増える 有利になる場合がある

精度については、時間刻みと空間刻みの両方が影響します。

陰解法だから常に高精度というわけではなく、問題に合った格子設定が必要です。

クランクニコルソン法の位置付け

陽解法と陰解法の中間的な考え方として、クランクニコルソン法があります。

現在時刻と次時刻の情報を平均的に使うため、時間方向の精度を高めやすい手法です。

ただし、陰解法と同様に連立方程式を解く必要があります。

急激な条件変化がある場合には、時間刻みの選び方によって振動が見られることもあります。

実務では、単純な試算なら陽解法、安定した長時間解析なら陰解法、精度との両立を狙う場合はクランクニコルソン法というように使い分けることがあります。

手法名だけで決めるのではなく、材料特性、温度範囲、計算時間、許容誤差を整理して選ぶことが重要です。

熱伝導解析における実践上の注意点

続いては、熱伝導方程式を実際の解析へ使う際の注意点を確認していきます。

物性値と単位系の確認

熱伝導解析では、熱伝導率、密度、比熱、熱拡散率などの物性値を使用します。

これらの値は温度によって変化することがあり、常温の値を高温環境へそのまま使うと誤差が大きくなる場合があります。

特に金属は高温になると物性値が変化し、樹脂や断熱材では材料内部の状態によって値が変わることもあります。

入力値の信頼性は計算結果の信頼性を左右する要素です。

また、mm、cm、m の混在や、秒と分の取り違えにも注意してください。

単位をあらかじめ SI単位系へ統一すると、計算式の確認がしやすくなります。

境界条件の現実性

実物の部品は、完全に断熱されていることも、表面温度が完全に一定であることも多くありません。

空気への放熱、接触面を通じた熱移動、放射による熱移動などが同時に起こります。

周囲との熱交換を考慮する場合は、対流熱伝達の境界条件を設定します。

高温域では熱放射の影響も無視できなくなるため、温度の四乗に関係する放射熱を考慮する必要があるでしょう。

解析モデルは、現実をすべて再現するためではなく、目的に必要な現象を適切に取り入れるためのものです。

過度に単純化すると誤差が増え、過度に複雑化すると入力条件の不確かさが増える場合があります。

たとえば冷却時間を見積もりたいだけなら、主要な放熱経路を押さえた一次元モデルでも有効なことがあります。

局所的な温度上昇や熱応力まで評価するなら、二次元や三次元の有限要素解析が必要になる場合もあります。

結果の検証と可視化

数値計算の結果は、数値が出た時点で正しいと判断してはいけません。

温度が物理的にあり得る範囲か、時間とともに不自然な振動をしていないかを確認します。

単純な条件については、変数分離法による解析解と差分法の結果を比較する方法も有効です。

格子幅や時間刻みを半分にして再計算し、結果が大きく変わらないかを調べる作業は格子収束の確認につながります。

計算条件を細かくしても結果がほとんど変わらない状態を確認できれば、設定の妥当性を判断しやすくなります。

温度分布をグラフ、等温線、動画などで可視化すると、熱の移動経路や高温部の残り方を直感的に把握できます。

熱伝導方程式の解法のまとめ

熱伝導方程式は、材料内部の温度分布が時間とともにどのように変化するかを表す基本的な方程式です。

単純な形状と境界条件であれば、変数分離法とフーリエ級数を使い、解析解として温度変化を表せます。

複雑な条件や実務的なモデルでは、有限差分法などの数値解法が役立ちます。

陽解法は計算式が分かりやすい反面、αΔt/Δx² が 2分の1以下となる安定条件を守る必要があります。

陰解法は連立方程式を解く手間があるものの、時間刻みを大きく取りやすく、長時間の解析で有力な選択肢になります。

どの方法を選ぶ場合でも、初期条件、境界条件、物性値、単位系を正しく設定することが欠かせません。

まずは一次元の棒を対象に、解析解と差分法の両方を試してみると、熱が拡散する仕組みと計算手順を深く理解しやすくなるでしょう。