Переходная теплопроводность: теплоёмкость, начальное решение и шаг времени

Почему тепловые процессы медленные

Электромагнитные поля в устройствах устанавливаются за доли секунды, тепловые — за минуты и часы. Обмотка двигателя после включения греется десятки минут, бетонная стена «чувствует» суточный ход температуры с задержкой в несколько часов, а при коротком замыкании провод успевает перегреться быстрее, чем сработает защита. Все эти задачи — о том, как температурное поле меняется во времени: за это отвечает левая часть уравнения теплопроводности, которой мы до сих пор пренебрегали.

Нестационарные задачи важны и по другой причине: в них проявляется разница между скоростью нагрева (как быстро) и степенью неравномерности температуры в теле (насколько поверхность опережает центр). Эти два вопроса управляются двумя разными безразмерными числами, и путать их — самая частая ошибка.

Нестационарное уравнение и температуропроводность

Полное уравнение теплопроводности, выведенное в уроке о законе Фурье: $$\rho_mc\,\frac{\partial T}{\partial t}=\nabla\cdot (\lambda\,\nabla T)+q_V .$$ Здесь $\rho_m$ — плотность массы, кг/м³, $c$ — удельная теплоёмкость, Дж/(кг·К). Произведение $\rho_mc$ — объёмная теплоёмкость — показывает, сколько энергии надо вложить в единицу объёма, чтобы поднять его температуру на один кельвин; оно характеризует инерционность тела. Для однородной среды без источников уравнение делится на $\rho_mc$: $$\frac{\partial T}{\partial t}=a\,\Delta T,\qquad a=\frac{\lambda}{\rho_mc}.$$

Любопытно, что у воздуха температуропроводность выше, чем у стали: тепло он проводит плохо, но и запасать его почти не способен.

Для полной постановки задачи к уравнению и граничным условиям добавляется начальное условие — распределение температуры в момент $t=0$: $$T(\vec r,0)=T_0(\vec r).$$ Чаще всего начальное распределение однородно: $T_0=\mathrm{const}$. Как и в стационарной задаче, удобно работать с превышением $\vartheta=T-T_\infty$ над температурой среды — уравнение для него то же.

Характерное время

Уравнение $\partial T/\partial t=a\,\Delta T$ содержит единственный размерный параметр $a$ с единицей м²/с. Если характерный размер тела $L$, то из $a$ и $L$ можно составить ровно одну величину с размерностью времени: $$\tau\sim\frac{L^2}{a}.$$ Это оценка времени, за которое температурное возмущение проходит расстояние $L$ внутри тела. Квадратичная зависимость от размера — главный качественный факт нестационарной теплопроводности: вдвое более толстая пластина прогревается в четыре раза дольше. Для стальной пластины $L=10$ мм: $L^2/a=10^{-4}/(1{,}3\cdot10^{-5})\approx7{,}7$ с; для бетонной стены $L=0{,}1$ м: $10^{-2}/(7{,}4\cdot10^{-7})\approx1{,}35\cdot10^4$ с, почти четыре часа. Аналогия с электромагнитной диффузией поля в проводник здесь полная: уравнение то же, роль $a$ играет $1/(\mu\sigma)$.

Два безразмерных числа

Рассмотрим тело с характерным размером $L$ и теплопроводностью $\lambda$, омываемое средой с коэффициентом теплоотдачи $\alpha$. У процесса два «сопротивления»: внутреннее, кондуктивное $\sim L/\lambda$, и внешнее, конвективное $\sim1/\alpha$.

Число Био характеризует форму температурного профиля (насколько тело неизотермично), число Фурье — стадию процесса. Решение любой задачи о прогреве тела простой формы с однородным начальным условием и конвекцией на границе имеет вид $\vartheta/\vartheta_0=f(x/L,\ \mathrm{Bi},\ \mathrm{Fo})$; для пластин это и есть то, что показывает лаборатория.

Модель тела с сосредоточенной ёмкостью

При $\mathrm{Bi}<0{,}1$ перепад температуры внутри тела не превышает нескольких процентов перепада между телом и средой, и тело можно считать точкой с одной температурой $T(t)$. Такой подход называют моделью с сосредоточенной теплоёмкостью.

Условие применимости $\mathrm{Bi}<0{,}1$ записывают через характерный размер $L_c=V/S$: для пластины толщиной $2L$ это $L$, для цилиндра радиуса $a$ — $a/2$, для шара — $a/3$. Заметим, что $\tau$ в этой модели не зависит от $\lambda$ — вся динамика определяется ёмкостью тела и теплоотдачей, а не теплопроводностью. Показатель экспоненты можно записать и через безразмерные числа: $t/\tau=\alpha S t/(\rho_mcV)=\mathrm{Bi}\cdot\mathrm{Fo}$ при $L=V/S$.

Прогрев пластины: метод разделения переменных

Пластина толщиной $2L$ ($-L\le x\le L$) с начальной температурой $T_0$ с момента $t=0$ омывается с обеих сторон средой $T_\infty$ с теплоотдачей $\alpha$. Задача симметрична, достаточно половины $0\le x\le L$. Для $\vartheta=T-T_\infty$: $$\frac{\partial\vartheta}{\partial t}=a\frac{\partial^2\vartheta}{\partial x^2};\qquad \frac{\partial\vartheta}{\partial x}\Big|_{x=0}=0;\qquad -\lambda\frac{\partial\vartheta}{\partial x}\Big|_{x=L}=\alpha\,\vartheta(L,t);\qquad \vartheta(x,0)=\vartheta_0 .$$

Ряд сходится тем быстрее, чем больше $\mathrm{Fo}$: члены убывают как $e^{-\mu_n^2\mathrm{Fo}}$, а $\mu_2>\pi$, так что уже при $\mathrm{Fo}>0{,}2$ второй член меньше первого примерно в $e^{-(\mu_2^2-\mu_1^2)\cdot0{,}2}$ раз — на практике это менее $1$–$2$ %. Отсюда одночленное приближение (стадия регулярного режима): $$\frac{\vartheta(x,t)}{\vartheta_0}\approx C_1\cos\frac{\mu_1x}{L}\,e^{-\mu_1^2\mathrm{Fo}},\qquad \mathrm{Fo}>0{,}2 .$$ В регулярном режиме температура в каждой точке убывает по одной экспоненте с темпом $\mu_1^2a/L^2$, а отношение превышений температуры поверхности и центра постоянно и равно $\cos\mu_1$. Предельные случаи поучительны: при $\mathrm{Bi}\to0$ уравнение $\mu\operatorname{tg}\mu=\mathrm{Bi}$ даёт $\mu_1^2\approx\mathrm{Bi}$, $C_1\to1$, и одночленное приближение переходит в модель сосредоточенной ёмкости $e^{-\mathrm{Bi}\cdot\mathrm{Fo}}$; при $\mathrm{Bi}\to\infty$ (граничное условие первого рода) $\mu_1\to\pi/2$, $C_1\to4/\pi$.

Полуограниченное тело

Если тело толстое, а время мало ($\mathrm{Fo}\ll1$), возмущение с поверхности не успевает дойти до противоположной стороны, и тело можно считать полуограниченным ($x\ge0$). При скачке температуры поверхности от $T_0$ к $T_s$ в момент $t=0$ решение выражается через дополнительную функцию ошибок [Тихонов]: $$\frac{T(x,t)-T_0}{T_s-T_0}=\operatorname{erfc}\frac{x}{2\sqrt{at}},\qquad \operatorname{erfc}z=\frac{2}{\sqrt\pi}\int_z^\infty e^{-u^2}\dd u .$$ Профиль зависит только от комбинации $x/\sqrt{at}$: глубина прогрева растёт как $\sqrt{at}$. Так как $\operatorname{erfc}(1{,}8)\approx0{,}011$, на глубине $x\approx3{,}6\sqrt{at}$ изменение температуры составляет около $1$ % — это условная глубина проникновения. Для бетона за сутки $\sqrt{at}=\sqrt{7{,}4\cdot10^{-7}\cdot86\,400}=0{,}25$ м, глубина проникновения около $0{,}9$ м: суточные колебания температуры воздуха почти не ощущаются на глубине метра. Плотность потока на поверхности $q(0,t)=\lambda(T_s-T_0)/\sqrt{\pi at}$ убывает как $1/\sqrt t$ и в начальный момент неограниченна — скачок температуры поверхности физически требует бесконечно большой теплоотдачи.

Свойства, зависящие от температуры

Если $\lambda$, $c$ или $\rho_m$ зависят от $T$, уравнение $$\rho_m(T)\,c(T)\,\frac{\partial T}{\partial t}=\nabla\cdot \big(\lambda(T)\,\nabla T\big)+q_V$$ нелинейно. Разделение переменных и суперпозиция теряют силу: решение при удвоенном перепаде не равно удвоенному решению, собственные функции не существуют. Преобразование Кирхгофа, спасавшее стационарную задачу, здесь не помогает: коэффициент при производной по времени, $\rho_mc\,\lambda_0/\lambda(T)$, остаётся зависящим от температуры. Остаётся численное решение, в котором на каждом шаге по времени свойства берутся при текущей температуре. Для стали в диапазоне до $500$ °C изменение $\lambda$ и $c$ достигает десятков процентов, так что при больших перегревах пренебрегать им нельзя.

Численное решение: шаги по времени

Разобьём пластину на узлы $x_i$ с шагом $\Delta x$, время — на слои $t^n$ с шагом $\Delta t$; $T_i^n$ — температура узла $i$ в момент $n$. Вторую производную заменяем разностью $(T_{i+1}-2T_i+T_{i-1})/\Delta x^2$, первую по времени — $(T_i^{n+1}-T_i^n)/\Delta t$. Вопрос только в том, в какой момент брать пространственную разность.

• Явная схема: берём её на известном слое $n$, $$T_i^{n+1}=T_i^n+\frac{a\Delta t}{\Delta x^2}\left(T_{i+1}^n-2T_i^n+T_{i-1}^n\right).$$ Каждое новое значение вычисляется напрямую. Но схема устойчива лишь при $$\frac{a\,\Delta t}{\Delta x^2}\le\frac12 :$$ иначе ошибки округления растут от шага к шагу и решение «разваливается» в пилу. Смысл условия: за шаг по времени тепло не должно диффундировать дальше соседнего узла. Измельчение сетки вдвое требует уменьшения $\Delta t$ вчетверо — для тонких сеток явная схема становится непомерно дорогой.
• Неявная схема: разность берём на новом слое $n+1$. На каждом шаге решается система линейных уравнений с трёхдиагональной матрицей (это недорого), зато схема устойчива при любом $\Delta t$; точность по времени — первого порядка.
• Схема Кранка–Николсон: полусумма явной и неявной. Безусловно устойчива и имеет второй порядок точности по времени, поэтому её применяют чаще всего [Самарский].

Устойчивость не означает точности: при слишком большом $\Delta t$ неявные схемы дают гладкое, но неверное решение. Правило — $\Delta t$ должен быть малым по сравнению с характерным временем $L^2/a$ процесса, который мы хотим разрешить; проверка — повторный расчёт с уменьшенным шагом (сходимость).

Лаборатория

Лаборатория моделирует пластину толщиной $2L$ с начальной температурой $T_0$, которая с момента $t=0$ охлаждается или нагревается конвекцией ($\alpha$, $T_\infty$) с обеих сторон (используется симметрия) либо получает заданную температуру поверхности. Численное решение строится схемой Кранка–Николсон (60 узлов по половине толщины, предвычисленная временная шкала) и сравнивается с аналитическим рядом из 30 членов. Слайдер времени и воспроизведение позволяют проследить эволюцию профиля. Выводятся $a$, $\mathrm{Bi}$, $\mathrm{Fo}$, температуры центра и поверхности (численно и по ряду), накопленное на единицу площади тепло и его отношение к предельному значению $\rho_mc\cdot2L\,(T_\infty-T_0)$, время достижения центром $90$ % перепада; графики $T(x)$ в текущий момент и истории $T_{\text{ц}}(t)$, $T_{\text{п}}(t)$. Опция зависимости $\lambda(T)$ доступна только в численном решении — аналитический ряд для неё не существует.

Что стоит проделать:

• Задайте стальную пластину из примера ($\mathrm{Bi}=0{,}02$). Убедитесь, что кривые центра и поверхности сливаются, а время достижения $90$ % совпадает с расчётом.
• Увеличьте $\alpha$ до $5000$ (вода) при той же пластине: $\mathrm{Bi}=1$. Поверхность теперь заметно опережает центр; вычислите по одночленному приближению $\vartheta_{\text{ц}}/\vartheta_0$ при нескольких $\mathrm{Fo}$, сравните с показаниями лаборатории и найдите, с какого $\mathrm{Fo}$ они совпадают.
• Удвойте толщину при прочих равных и измерьте, во сколько раз выросло время достижения $90$ %. Повторите при малом и при большом $\mathrm{Bi}$ и объясните разницу.
• Переключитесь на заданную температуру поверхности ($\mathrm{Bi}\to\infty$) и проверьте, что $\vartheta_{\text{ц}}/\vartheta_0$ в регулярном режиме убывает как $(4/\pi)e^{-\pi^2\mathrm{Fo}/4}$.
• Включите $\lambda(T)$ и посмотрите, как аналитический ряд расходится с численным решением на стадии большого перегрева и сходится с ним к концу процесса.

Итоги

• Нестационарное уравнение $\rho_mc\,\partial T/\partial t=\nabla\cdot (\lambda\nabla T)+q_V$ требует начального условия; скорость выравнивания температуры задаёт температуропроводность $a=\lambda/(\rho_mc)$, а характерное время $\tau\sim L^2/a$ растёт квадратично с размером.
• Число Био $\mathrm{Bi}=\alpha L/\lambda$ определяет форму профиля (неизотермичность тела), число Фурье $\mathrm{Fo}=at/L^2$ — стадию процесса.
• При $\mathrm{Bi}<0{,}1$ работает модель сосредоточенной ёмкости: $\vartheta=\vartheta_0e^{-t/\tau}$, $\tau=\rho_mcV/(\alpha S)$ — аналог $RC$-цепи.
• Прогрев пластины описывается рядом по собственным функциям $\cos(\mu_nx/L)$ с корнями уравнения $\mu\operatorname{tg}\mu=\mathrm{Bi}$; при $\mathrm{Fo}>0{,}2$ достаточно первого члена.
• Зависимость свойств от температуры делает задачу нелинейной и доступной только численно; явная схема устойчива при $a\Delta t/\Delta x^2\le1/2$, неявные схемы (в том числе Кранка–Николсон) устойчивы всегда, но точность всё равно требует проверки сходимости.

Далее — связанные электромагнитно-тепловые задачи: потери как источники тепла и обратное влияние температуры на проводимость.