Przejdź do treści

Nieustalona analiza przewodzenia ciepła

W tej sekcji przedstawiono dyskretyzację czasu oraz iteracyjną metodę rozwiązania analizy przewodzenia ciepła w ciałach stałych metodą elementów skończonych (MES). Równania rządzące i warunki brzegowe kontinuum opisano w Równaniu przewodzenia ciepła.

Równanie zdyskretyzowane (punkt wyjścia)

Dyskretyzacja równania przewodzenia ciepła (równanie przewodzenia ciepła (gov_he_main)) metodą Galerkina prowadzi do

\[\begin{equation} K T + M \frac{\partial T}{\partial t} = F \label{eq:2.4.8} \end{equation}\]

gdzie

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

Tutaj \(K\), \(M\), \(F\) i \(N\) oznaczają odpowiednio macierz przewodzenia ciepła (wraz z wkładami brzegowymi od konwekcji i promieniowania), macierz masy, wektor obciążenia cieplnego i macierz funkcji kształtu. Definicje symboli właściwości materiałowych (\(\rho\), \(c\), \(k_x, k_y, k_z\), \(Q\), \(hc\), \(hr\) itd.) są zgodne z Równaniem przewodzenia ciepła.

Dyskretyzacja czasu i rozwiązanie iteracyjne

Równanie \(\eqref{eq:2.4.8}\) jest nieliniowym równaniem nieustalonym. Przy zastosowaniu wstecznej metody Eulera do dyskretyzacji czasu, gdy znana jest temperatura w chwili \(t=t_0\), temperaturę w chwili \(t=t_0+\Delta t\) oblicza się z następującego równania.

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

Rozważmy poprawienie wektora temperatury \(T_{t=t_0+\Delta t}^{(i)}\), który w przybliżeniu spełnia równanie \(\eqref{eq:2.4.13}\), tak aby otrzymać dokładniejsze rozwiązanie \(T_{t=t_0+\Delta t}^{(i)+1}\).

W tym celu najpierw przedstawia się wektor temperatury w następującej postaci.

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

Iloczyn macierzy przewodzenia ciepła i wektora temperatury, macierz masy oraz powiązane człony przybliża się w następujący sposób.

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

Po podstawieniu równania \(\eqref{eq:2.4.14}\), równania \(\eqref{eq:2.4.15}\) i równania \(\eqref{eq:2.4.16}\) do równania \(\eqref{eq:2.4.13}\) oraz pominięciu członów drugiego i wyższych rzędów otrzymuje się następujące równanie.

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

Macierz współczynników po lewej stronie jest dalej przybliżana za pomocą następującego równania.

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

Tutaj \(K_{T_{t=t_0+\Delta t}}^{(i)}\) jest macierzą sztywności stycznej.

Ostatecznie temperaturę w chwili \(t=t_0+\Delta t\) można obliczyć, wykonując iteracje z użyciem następującego równania.

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

W szczególności w analizie ustalonej iteracje wykonuje się z użyciem następującego równania.

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

W analizie nieustalonej, ponieważ do dyskretyzacji czasu stosowana jest metoda niejawna, wybór przyrostu czasu \(\Delta t\) zasadniczo nie podlega ograniczeniu co do jego wielkości. Jeśli jednak przyrost czasu \(\Delta t\) jest zbyt duży, wzrasta liczba iteracji potrzebnych do osiągnięcia zbieżności. Ogólnie zbyt duży przyrost czasu \(\Delta t\) zwiększa liczbę iteracji. W implementacji automatyczne sterowanie przyrostem monitoruje wielkość wektora residuum, zmniejszając \(\Delta t\), gdy zbieżność jest wolna, i zwiększając \(\Delta t\), gdy liczba iteracji jest mała (→ szczegóły: Sterowanie krokiem).

Tematy powiązane