אנליזת הולכת חום חולפת
מוצגות דיסקרטיזציית הזמן ושיטת הפתרון האיטרטיבית לאנליזת הולכת חום במוצקים באמצעות שיטת האלמנטים הסופיים (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}\]
נבחן שיפור של וקטור הטמפרטורה \(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\) מוגדל (לפרטים ראו בקרת צעדים).
נושאים קשורים