暫態熱傳導分析
本節說明使用有限元素法（Finite Element Method, FEM）進行固體熱傳導分析時的時間離散化與迭代求解法。連續體的控制方程式與邊界條件請參閱熱傳導方程式。
離散化方程式（起點）
將熱傳導方程式（熱傳導方程式 (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}\) 是非線性且暫態的方程式。
現在以後向 Euler 法對時間進行離散化，當時刻 \(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\)，以進行自動增量控制（→ 詳情請參閱步驟控制）。
相關項目