การวิเคราะห์การนำความร้อนไม่คงตัว
แสดงการดิสครีไทซ์ตามเวลาและวิธีแก้แบบวนซ้ำสำหรับการวิเคราะห์การนำความร้อนของของแข็งด้วยวิธีไฟไนต์เอลิเมนต์ (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}\) เป็นสมการไม่เชิงเส้นและไม่คงตัว ให้ดิสครีไทซ์ตามเวลาด้วยวิธี Backward 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\) เมื่อจำนวนรอบมีน้อย (→ ดูรายละเอียดที่ การควบคุมสเต็ป)
หัวข้อที่เกี่ยวข้อง