Распространение тепла (диффузия)

В качестве примера диффузионного процесса мы рассмотрим распространение тепла в неподвижной среде. Эта задача формулируется следующим образом: заданы начальное распределение температуры, граничные условия, источники тепла и требуется найти распределение температуры в последующие моменты времени.

Уравнение теплопроводности в декартовой системе координат имеет вид.

где и (х, t) — температура среды в точке х = (х_ь х₂, ^х?>) в момент времени t; р (х, t) — плотность; с (х, t) — теплоемкость; а (х, t) — коэффициент теплопроводности и /_s(x, t) — источник тепла. Вначале мы рассмотрим основные идеи построения разностных схем для задач теплопроводности на примере простого одномерного уравнения с постоянными коэффициентами.

Разностная схема для одномерного уравнения теплопроводности с постоянными коэффициентами

В случае одномерной однородной среды предыдущее уравнение принимает вид.

где к = а/ср — коэффициент температуропроводности. Конкретная физическая задача определяется начальным условием (начальным распределением температуры) и (х, 0) = g (x) и граничными условиями, которые мы запишем в общем виде.

где /2(0 — заданные функции.

Введем разностную сетку с узлами (х", 4) и шагами /гит. Аппроксимируем производную по времени разностью впе;

ч д²гг ред относительно точки (х", 4). а производную —5- трехто;

Зх чечным разностным оператором в момент времени 4. В результате дифференциальное уравнение (11.19) заменяется разностной схемой.

+ разностная форма граничных условий.

Так как эта схема явная, то процесс вычислений организуется очень просто. Запишем предыдущее выражение в виде.

Так как распределение температуры в момент времени Ц всегда известно, то мы получили формулу для расчета распределения температуры в момент времени Рассмотренный нами пример является частным случаем более общей схемы вида.

которая включает в себя все возможные двухслойные схемы. Шаблон этой схемы показан на рис. 11.7.

A/c. 7 /.7. Шаблон для схемы (11.21).

Введение

дополнительного параметра позволяет управлять свойствами разностной схемы, но это приводит к неявным разностным уравнениям. Решение таких уравнений мы обсудим ниже.

Рассмотрим погрешность аппроксимации схемы (11.21). Подставляя точное решение в (11.21) и используя соотношения (11.14), получим.

Проведем исследование аппроксимации в точке (.r_m ?/,). Тогда величины, определенные в точке (х_п, ^₊₁), необходимо разложить в ряд Тейлора в точке, но степеням т. Окончательно мы приходим к следующему дифференциальному приближению:

из которого следует, что в общем случае схема (11.3) аппроксимирует уравнение (11.1) с точностью до членов 0(х + /г²). Выбирая параметр у определенным образом, можно повысить порядок аппроксимации схемы (11.11). Действительно, из уравнения (11.19) следует, что.

Тогда дифференциальное приближение можно записать в виде.

Отсюда следует, что при.

выражение в квадратных скобках равно нулю и разностная схема (11.18) аппроксимирует уравнение (11.19) с точностью до величин 0(т² + /г⁴).

Для анализа устойчивости схемы (11.21) подставим в нее решение вида (11.14) и получим.

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

В данном случае Х (а) принимает вещественные значения в интервале.

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

и мы приходим к условию.

В частности, для явной схемы (11.20) (у = 0) условие устойчивости приобретает вид р < 0.5. Как следует из выражения (11.22), эта схема будет иметь четвертый порядок аппроксимации при еще более строгом ограничении р = 1/6. Напротив, схема (11.21) будет устойчива при любом значении р, если у определяется выражением (11.22).

Перейдем теперь к обсуждению важного вопроса, связанного с аппроксимацией граничных условий. В случае граничных условий первого рода, когда на границе задается значение температуры, представление этих условий в разностной форме не представляет особого труда. Например, заданное при х = 0 условие и (0, t) = f (t) записывается как.

Граничные условия второго и третьего рода определяют значения теплового потока в граничных точках. В качестве примера можно привести следующие физические процессы, приводящие к условиям такого рода:

1) воздействие теплового потока на вещество.

2) конвективный теплообмен.

где и_с — температура окружающей среды; k₀ — коэффициент теплоотдачи;

3) тепловое излучение

где, а — коэффициент излучения.

Для простоты рассмотрим граничное условие (11.25). Очевидное представление этого условия в разностной форме.

непригодно, так как оно имеет первый порядок аппроксимации по h, что значительно снижает точность приближенного решения (пример 7.11). Для построения более точной аппроксимации, второго порядка по h, применим метод фиктивных областей. Для этого введем фиктивный слой с узлами (х_(, tjf). Тогда, используя центральную разность, граничное условие (11.25) можно записать в виде.

Значения температуры в граничном узле вычисляется из разностного уравнения (11.3) при п = 0:

Исключение из этих двух выражений, фиктивных значений с номером п = -1 дает разностную форму граничного условия (11.7), которую мы запишем в виде

В случае явной схемы (у = 0) это выражение непосредственно дает формулу для вычисления температуры в граничном узле.

В общем случае, когда у ^ 0, разностная схема (11.3) является неявной, так как при каждом значении п разностное уравнение включает в себя неизвестные в момент времени tk+ значения температуры в трех узлах сетки. Для того чтобы вычислить решение в момент времени ?$₊₁, необходимо рассмотреть совокупность всех разностных уравнений.

Перепишем выражение (11.21), группируя неизвестные величины в левой части, а известные в правой:

Введем обозначения.

Тогда предыдущее выражение представляет собой систему линейных уравнений вида.

с трехдиагональной матрицей.

Коэффициенты этой матрицы следуют из (11.27) и равны.

а коэффициенты первого и последнего уравнений определяются граничными условиями. Например, для граничного условия (11.24) s₀ = 1, г₀ ⁼ 0 и d₀ = /(?*+ j), а для граничного условия (11.26).

Так как элементы матрицы, А удовлетворяют условию s_n > р_п + г_п, то для решения системы уравнений можно применить метод прогонки (параграф 3.6.1). Таким образом, для выполнения одного шага по времени в схеме (11.21) требуется большее количество операций по сравнению с явной схемой (11.20). С другой стороны, в неявной схеме мы можем выбрать больший шаг по времени, поэтому при заданной длительности процесса потребуется меньшее число шагов и общее количество операций будет сравнимо или даже меньше, чем для явной схемы. Для вычисления одного временного шага по явной схеме (11.20) необходимо AN операций. Для вычислений по неявной схеме (11.21) сначала необходимо определить вектор d (AN операций), а затем решить систему уравнений (-8N операций), что в сумме дает приблизительно 12JV операций на один шаг по времени. Поэтому если в неявной схеме шаг по времени выбрать из соотношения р ~ 1.5, то общее количество операций будет сравнимо с количеством операций, требуемых для вычисления по явной схеме.

Пример 11.1 (остывание бесконечной пластины) Рассмотрим бесконечную пластину толщины /. Пусть в момент времени t = 0 она имеет температуру u_Q, а температура окружающей среды поддерживается равной нулю. Тогда распределение температуры в пластине описывается уравнением (11.1) с условиями.

Точное решение этой задачи имеет вид.

На рис. 11.8 показано распределение температуры в момент времени Т = 0.1 при к = 1, / = 1 и щ = 1. Для оценки ошибки приближенного решения введем вектор приближенного решения в момент времени ^.

и точное решение, вычисленное в узлах разностной сетки в момент времени ?*.

Тогда относительная ошибка приближенного решения определяется как.

Для решения, показанного на рис. 11.8, ошибка e (v*) 0.0347.

Рис. 11.8. Решение уравнения (11.19):

• —точное решение;——приближенное решение
• (явная схема (11.20), h = 0.1, т = 0.005)

На рис. 11.9 показано распределение температуры в момент времени Т= 0.1, полученное с помощью схемы (11.21). При этом значение т выбиралось из условия (11.22), а шаг КТ по времени выоирался из условия р = — = 1.66 666. Для это- . «/г го решения е (у») = 0.0136, что демонстрирует эффективность применения неявной схемы для решения такого типа задач.

Рис. 11.9. Решение уравнения (11.19) с условиями (11.29):

—точное решение;——приближенное решение (схема (11.21),.

h = 0.1, т= 1.66 666-Ю'², у = 0.45).