Phân tích dẫn nhiệt không ổn định
Trình bày rời rạc hóa theo thời gian và phương pháp giải lặp cho phân tích dẫn nhiệt trong vật rắn bằng phương pháp phần tử hữu hạn (Finite Element Method). Về phương trình chi phối và điều kiện biên trong môi trường liên tục, xem Phương trình dẫn nhiệt.
Phương trình rời rạc hóa (điểm xuất phát)
Khi rời rạc hóa phương trình dẫn nhiệt (phương trình dẫn nhiệt (gov_he_main)) bằng phương pháp Galerkin, ta được
\[\begin{equation} K T + M \frac{\partial T}{\partial t} = F \label{eq:2.4.8} \end{equation}\]
trong đó
\[\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}\]
Trong đó, \(K\), \(M\), \(F\), \(N\) lần lượt là ma trận dẫn nhiệt (bao gồm các số hạng đối lưu và bức xạ từ đóng góp biên), ma trận khối lượng, vectơ tải nhiệt và ma trận hàm dạng. Định nghĩa các ký hiệu tính chất vật lý (\(\rho\), \(c\), \(k_x, k_y, k_z\), \(Q\), \(hc\), \(hr\), v.v.) theo Phương trình dẫn nhiệt.
Rời rạc hóa theo thời gian và phương pháp giải lặp
Phương trình \(\eqref{eq:2.4.8}\) là phương trình phi tuyến và không ổn định. Rời rạc hóa theo thời gian bằng phương pháp Euler lùi; khi nhiệt độ tại thời điểm \(t=t_0\) đã biết, nhiệt độ tại \(t=t_0+\Delta t\) được tính bằng phương trình sau.
\[\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}\]
Xét việc cải thiện vectơ nhiệt độ \(T_{t=t_0+\Delta t}^{(i)}\) thỏa xấp xỉ phương trình \(\eqref{eq:2.4.13}\) để thu được nghiệm chính xác hơn \(T_{t=t_0+\Delta t}^{(i)+1}\).
Trước hết, biểu diễn vectơ nhiệt độ như sau.
\[\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}\]
Tích của ma trận dẫn nhiệt và vectơ nhiệt độ, ma trận khối lượng, v.v. được biểu diễn xấp xỉ như phương trình sau.
\[\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}\]
Thay các phương trình \(\eqref{eq:2.4.14}\), \(\eqref{eq:2.4.15}\), \(\eqref{eq:2.4.16}\) vào phương trình \(\eqref{eq:2.4.13}\) và bỏ các số hạng bậc hai trở lên, ta được phương trình sau.
\[\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}\]
Tiếp theo, ma trận hệ số ở vế trái được đánh giá xấp xỉ bằng phương trình sau.
\[\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}\]
Trong đó, \(K_{T_{t=t_0+\Delta t}}^{(i)}\) là ma trận độ cứng tiếp tuyến.
Cuối cùng, có thể tính nhiệt độ tại thời điểm \(t=t_0+\Delta t\) bằng cách lặp theo phương trình sau.
\[\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}\]
Đặc biệt, trong phân tích ổn định, phép lặp được thực hiện bằng phương trình sau.
\[ 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}\]
Trong phân tích không ổn định, vì rời rạc hóa theo thời gian dùng phương pháp ẩn nên nói chung việc chọn bước thời gian \(\Delta t\) không bị ràng buộc về kích thước. Tuy nhiên, nếu bước thời gian \(\Delta t\) quá lớn, số lần lặp cần để hội tụ sẽ tăng. Nói chung, khi bước thời gian \(\Delta t\) quá lớn, số lần lặp tăng. Trong triển khai, độ lớn của vectơ phần dư được theo dõi và dùng điều khiển bước tăng tự động: giảm \(\Delta t\) khi hội tụ chậm và tăng \(\Delta t\) khi số lần lặp ít (→ xem Điều khiển bước để biết chi tiết).
Chủ đề liên quan