非定常熱伝導解析
有限要素法（Finite Element Method）による固体の熱伝導解析の時間離散化と反復解法を示す。連続体での支配方程式・境界条件は熱伝導方程式を参照。
離散化方程式（出発点）
熱伝導方程式（熱伝導方程式の式 (gov_he_main)）を Galerkin 法によって離散化すると、
\[\begin{equation}
K T + M \frac{\partial T}{\partial t} = F
\label{eq:2.4.8}
\end{equation}\]
ただし、
\[\begin{equation}
K = \int\left( k_x \frac{\partial N^T}{\partial x}\frac{\partial N}{\partial x}
+ k_y \frac{\partial N^T}{\partial y}\frac{\partial N}{\partial y}
+ k_z \frac{\partial N^T}{\partial z}\frac{\partial N}{\partial z} \right) dV
+ \int hc N^T N ds + \int hr N^T N ds
\label{eq:2.4.9}
\end{equation}\]
\[\begin{equation}
M = \int \rho c N^T N dV
\label{eq:2.4.10}
\end{equation}\]
\[\begin{equation}
F = \int Q N^T dV - \int q_s N^T dS
+ \int{hc} T c N^T dS
+ \int{hcTr} ({T+Tr}) ({T^2 + T r^2}) N^T dS
\label{eq:2.4.11}
\end{equation}\]
\[\begin{equation}
N = (N^1, N^2, \ldots, Ni)
\label{eq:2.4.12}
\end{equation}\]
ここで、\(K\), \(M\), \(F\), \(N\) はそれぞれ熱伝導マトリクス（境界寄与の対流・輻射項を含む）、質量マトリクス、熱荷重ベクトル、形状関数行列である。物性記号（\(\rho\), \(c\), \(k_x, k_y, k_z\), \(Q\), \(hc\), \(hr\) 等）の定義は熱伝導方程式に従う。
時間離散化と反復解法
方程式\(\eqref{eq:2.4.8}\)は非線形かつ非定常の方程式である。
いま、時間に関して後退オイラー法により離散化して、時刻\(t=t_0\)における温度が既知のとき時刻\(t=t_0+\Delta t\)での温度を次式を用いて計算することにする。
\[\begin{equation}
K_{t=t_0+\Delta t} T_{t=t_0+\Delta t} + M_{t=t_0+\Delta t} \frac{T_{t=t_0+\Delta t} - T_{t=t_0}}{\Delta t} = F_{t=t_0+\Delta t}
\label{eq:2.4.13}
\end{equation}\]
ここでの式\(\eqref{eq:2.4.13}\)を近似的にみたす温度ベクトル\(T_{t=t_0+\Delta t}^{(i)}\)
を改善して、精度の良い解\(T_{t=t_0+\Delta t}^{(i)+1}\)を求めることを考える。
そのために、まず、温度ベクトルを次のようにあらわす。
\[\begin{equation}
T_{t=t_0+\Delta t}=
T_{t=t_0+\Delta t}^{(i)} + \Delta T_{t=t_0+\Delta t}^{(i)}
\label{eq:2.4.14}
\end{equation}\]
熱伝導マトリクスと温度ベクトルとの積、質量マトリクスなどを次式のように近似的にあらわす。
\[\begin{equation}
K_{t=t_0+\Delta t} T_{t=t_0+\Delta t} =
K_{t=t_0+\Delta t}^{(i)} T_{t=t_0+\Delta t}^{(i)} +
\frac{\partial \big(K_{t=t_0+\Delta t}^{(i)} T_{t=t_0+\Delta t}^{(i)}\big) }
{\partial T_{t=t_0+\Delta t}^{(i)} } \Delta T_{t=t_0+\Delta t}^{(i)}
\label{eq:2.4.15}
\end{equation}\]
\[\begin{equation}
M_{t=t_0+\Delta t} =
M_{t=t_0+\Delta t}^{(i)}
+ \frac{\partial M_{t=t_0+\Delta t}^{(i)}}{\partial T_{t=t_0+\Delta t}^{(i)}}
\Delta T_{t=t_0+\Delta t}^{(i)}
\label{eq:2.4.16}
\end{equation}\]
式\(\eqref{eq:2.4.14}\)、式\(\eqref{eq:2.4.15}\)、式\(\eqref{eq:2.4.16}\)を式\(\eqref{eq:2.4.13}\)に代入して二次以上の項を省略すると次式を得る。
\[\begin{equation}
\bigg(\frac{M_{t=t_0+\Delta t}^{(i)}}{\Delta t}
+ \frac {\partial M_{t=t_0+\Delta t}^{(i)} }
{ \partial T_{t=t_0+\Delta t}^{(i)} }
\frac{T_{t=t_0+\Delta t}^{(i)} - T_{t=t_0}}{\Delta t}
+ \frac{\partial \big(K_{t=t_0+\Delta t}^{(i)} T_{t=t_0+\Delta t}^{(i)}\big)}
{\partial T_{t=t_0+\Delta t}^{(i)}} \bigg)
\Delta T_{t=t_0+\Delta t}^{(i)} \\\
= F_{t=t_0+\Delta t} - M_{t=t_0+\Delta t}^{(i)}
\frac{T_{t=t_0+\Delta t}^{(i)} - T_{t=t_0}}{\Delta t}
- K_{t=t_0+\Delta t}^{(i)} T_{t=t_0+\Delta t}^{(i)}
\label{eq:2.4.17}
\end{equation}\]
さらに左辺の係数マトリクスを次式を用いて近似評価する。
\[\begin{equation}
K^{(i)} = \frac{M_{t=t_0+\Delta t}^{(i)}}{\Delta t} +
\frac{\partial \big( K_{t=t_0+\Delta t}^{(i)} T_{t=t_0+\Delta t}^{(i)} \big)}{\partial T^{(i)}_{t=t_0+\Delta t}}
= \frac{M_{t=t_0+\Delta t}^{(i)}}{\Delta t} + K_{T_{t=t_0+\Delta t}}^{(i)}
\label{eq:2.4.18}
\end{equation}\]
ここで\(K_{T_{t=t_0+\Delta t}}^{(i)}\)は接線剛性マトリクスである。
結局次式を用いて反復計算を行うことによって時刻 \(t=t_0+\Delta t\)での温度を計算することができる。
\[\begin{equation}
K^{(i)} \Delta T_{t=t_0+\Delta t}^{(i)} = F_{t=t_0+\Delta t} - M_{t=t_0+\Delta t}^{(i)} \frac{T_{t=t_0+\Delta t}^{(i)} - T_{t=t_0}}{\Delta t} - K^{(i)} T_{t=t_0+\Delta t}^{(i)}
\label{eq:2.4.19}
\end{equation}\]
特に定常解析においては次式を用いて反復計算を行う。
\[
K_T^{(i)} \Delta T_{t=\infty}^{(i)} = F_{t=\infty} - K_T^{(i)} \Delta T_{t=\infty}^{(i)}
\]
\[\begin{equation}
T_{t=\infty}^{(i+1)} = T_{t=\infty}^{(i)} + \Delta{T}_{t=\infty}^{(i)}
\label{eq:2.4.20}
\end{equation}\]
非定常解析において時間増分\(\Delta t\)の選び方は、時間に関する離散化に陰解法を採用しているので、一般にその大きさの制約を受けない。ただし時間増分\(\Delta t\)が大きすぎると、反復計算における収束回数は増加する。
一般に、時間増分\(\Delta t\)が大きすぎると反復回数が増加する。実装上は残差ベクトルの大きさを監視し、収束が遅ければ\(\Delta t\)を縮小、反復回数が少なければ\(\Delta t\)を拡大する自動増分制御が用いられる（→ 詳細はステップ制御を参照）。
関連項目