瞬态热传导分析
本节给出采用有限元法（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}\) 是非线性瞬态方程。
现在，对时间采用后退 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\)（→ 详情请参见分析步控制）。
相关项目