Перейти до змісту

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

У цьому розділі наведено дискретизацію за часом та ітераційний метод розв’язання задачі теплопровідності твердого тіла методом скінченних елементів. Визначальні рівняння та граничні умови для континууму див. у розділі Рівняння теплопровідності.

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

Після дискретизації рівняння теплопровідності (рівняння теплопровідності (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\), коли кількість ітерацій мала (→ докладніше див. Керування кроками).

Пов’язані теми

AI-assisted translation May contain errors Official docs Status