Перейти к содержанию

Нестационарный анализ теплопроводности

В этом разделе рассматриваются дискретизация по времени и итерационный метод решения задачи теплопроводности твёрдых тел методом конечных элементов (FEM). Основные уравнения и граничные условия для сплошной среды см. в разделе Уравнение теплопроводности.

Дискретизированное уравнение (исходная точка)

Дискретизация уравнения теплопроводности (уравнение теплопроводности (gov_he_main)) методом Галеркина даёт

\[\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}\]

Рассмотрим уточнение вектора температуры \(T_{t=t_0+\Delta t}^{(i)}\), который приближённо удовлетворяет уравнению \(\eqref{eq:2.4.13}\), для получения более точного решения \(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\), когда число итераций мало (→ подробности см. в разделе Управление шагами).

Связанные разделы