Разностная схема нестационарного теплового поля в многослойном цилиндре
дата публикации: 2016-04-21

В соответствии с законом Фурье для нестационарного теплового поля верно следующее соотношение:

\frac{d \left| \overrightarrow{T} \right| }{dt} = A \left( t, \left| \overrightarrow{T} \right|, \overrightarrow{e} \right) \cdot \Delta \overrightarrow{T} + \Psi \left( t, \left| \overrightarrow{T} \right|, \overrightarrow{e} \right) ,

где \Delta - лапласиан, \overrightarrow{T} - вектор теплового поля (температура), A \left( t, \left| \overrightarrow{T} \right|, \overrightarrow{e} \right) = \frac{\lambda \left( t, \left| \overrightarrow{T} \right|, \overrightarrow{e} \right) }{\rho \left( t, \left| \overrightarrow{T} \right|, \overrightarrow{e} \right) \cdot C \left( t, \left| \overrightarrow{T} \right|, \overrightarrow{e} \right) } - коэффициент, в котором \lambda \left( t, \overrightarrow{e} \right) - коэффициент теплопроводности, \rho \left( t, \overrightarrow{e} \right) - плотность, C \left( t, \overrightarrow{e} \right) - теплоемкость, \Psi \left( t, \left| \overrightarrow{T} \right|, \overrightarrow{e} \right) - внутренний источник тепловой энергии, \overrightarrow{e} - вектор системы координат.

Пользуясь уже известной формулой для лапласиана в цилиндрической системе координат:

\Delta \overrightarrow{F} = \frac{1}{r} \frac{\partial (r \cdot \partial F_r / \partial r)}{\partial r} + \frac{1}{r^2} \frac{\partial^2 F_{\phi}}{\partial \phi^2} + \frac{\partial^2 F_z}{\partial z^2},

учитывая однородность цилиндра по углу и подставляя вместо произвольного вектора \overrightarrow{F} вектор теплового поля \overrightarrow{T} , получим:

\Delta \overrightarrow{T} = \frac{1}{r} \frac{\partial (r \cdot \partial T / \partial r)}{\partial r} + \frac{\partial^2 T}{\partial z^2} = \frac{\partial^2 T}{\partial r^2} + \frac{1}{r} \frac{\partial T}{\partial r} + \frac{\partial^2 T}{\partial z^2}.

В результате получаем следующее уравнение, описывающее поле:

\frac{dT}{dt} = a(t, T,r,z) \cdot \left[ \frac{\partial^2 T}{\partial r^2} + \frac{1}{r} \frac{\partial T}{\partial r} + \frac{\partial^2 T}{\partial z^2} \right] + \psi(t, T,r,z).

Таким образом, пропадает необходимость дополнять уравнение граничными условиями на стыке сред. Единственные граничные условия, которые актуальны, это температура поверхности цилиндра S , которая меняется по какому-либо закону \Phi(t,S) и начальная температура всего объема цилиндра \Theta(r,z) .

Составим конечно-разностную схему этого уравнения, для чего воспользуемся следующими формулами:

\begin{aligned} & t_i = i \cdot h_t, i \in (0,n_t), \\ & r_j = j \cdot h_r, j \in (0,n_r), \\ & z_k = k \cdot h_z, k \in (0,n_z), \\ & T(t_i, r_j, z_k) = T_{i,j,k}, \\ & a(t_i, T(t_i, r_j, z_k), r_j, z_k) = a_{i,j,k}, \\ & \psi(t_i, T(t_i, r_j, z_k), r_j, z_k) = \psi_{i,j,k}, \\ & \frac{dT}{dt} \approx \frac{T_{i,j,k} - T_{i-1,j,k}}{h_t}, \\ & \frac{dT}{dr} \approx \frac{T_{i,j,k} - T_{i,j-1,k}}{h_r}, \\ & \frac{d^2T}{dr^2} \approx \frac{T_{i,j+1,k} - 2T_{i,j,k} + T_{i,j-1,k}}{h_r^2}, \\ & \frac{d^2T}{dz^2} \approx \frac{T_{i,j,k+1} - 2T_{i,j,k} + T_{i,j,k-1}}{h_z^2}, \end{aligned}

где h_t = t_{max}/n_t - шаг интегрирования по времени, равный времени процесса отнесенному к количеству точек, h_r = R/n_r - шаг интегрирования по радиусу, равный радиусу цилиндра отнесенному к количеству точек, h_z = L/n_z - шаг интегрирования по длине, равный длине отнесенной к количеству точек.

Перепишем аналитическое уравнение к разностному виду:

\frac{T_{i,j,k} - T_{i-1,j,k}}{h_t} = \psi_{i,j,k} + a_{i,j,k} \cdot \left[ \begin{aligned} & \frac{T_{i,j+1,k} - 2T_{i,j,k} + T_{i,j-1,k}}{h_r^2} + \frac{T_{i,j,k} - T_{i,j-1,k}}{j \cdot h_r^2} + \\ & + \frac{T_{i,j,k+1} - 2T_{i,j,k} + T_{i,j,k-1}}{h_z^2} \end{aligned} \right].

Проведем преобразования:

T_{i,j,k} - T_{i-1,j,k} = h_t \cdot \psi_{i,j,k} + \frac{a_{i,j,k} \cdot h_t}{j \cdot h_r^2 \cdot h_z^2} \cdot \cdot \left[ \begin{aligned} & j \cdot h_z^2 \cdot T_{i,j+1,k} - 2 \cdot j \cdot h_z^2 \cdot T_{i,j,k} + j \cdot h_z^2 \cdot T_{i,j-1,k} + h_z^2 \cdot T_{i,j,k} - h_z^2 \cdot T_{i,j-1,k} + \\ & + j \cdot h_r^2 \cdot T_{i,j,k+1} - 2 \cdot j \cdot h_r^2 \cdot T_{i,j,k} + j \cdot h_r^2 \cdot T_{i,j,k-1} \end{aligned} \right],\begin{aligned} & c_{i,j,k}^{(1)} = 1 + \frac{a_{i,j,k} \cdot h_t}{j \cdot h_r^2 \cdot h_z^2} \cdot \left[ (2 \cdot j - 1) \cdot h_z^2 + 2 \cdot j \cdot h_r^2 \right], \\ & c_{i,j,k}^{(2)} = \frac{a_{i,j,k} \cdot h_t \cdot (1-j)}{j \cdot h_r^2}, \\ & c_{i,j,k}^{(3)} = - \frac{a_{i,j,k} \cdot h_t}{h_r^2}, \\ & c_{i,j,k}^{(4)} = c_{i,j,k}^{(5)} = - \frac{a_{i,j,k} \cdot h_t}{h_z^2}, \\ & b_{i,j,k} = T_{i-1,j,k} + h_t \cdot \psi_{i,j,k}, \\ & c_{i,j,k}^{(1)} \cdot T_{i,j,k} + c_{i,j,k}^{(2)} \cdot T_{i,j-1,k} + c_{i,j,k}^{(3)} \cdot T_{i,j+1,k} + c_{i,j,k}^{(4)} \cdot T_{i,j,k-1} + c_{i,j,k}^{(5)} \cdot T_{i,j,k+1} = b_{i,j,k}. \end{aligned}

Полученное уравнение необходимо использовать для конструирования алгоритма, который рискнем построить на python. Не уверен, что устроит скорость интерпретатора, но все же следует попробовать именно на нем.