5. Нестационарная теплопроводность (неявный метод)

Основной недостаток явного численного метода состоит в том, что разностное уравнение баланса энергии для каждого узла должно удовлетворять критерию стойкости. Неявный метод основан на том, что у уравнения баланса энергии температуры определяются в моменты времени τ +Δτ.

Пример 5.1. Решить пример 4.1. неявным числовым методом. Параметры плиты и окружающей среды такие же.

Решение. Разбивка плиты на элементарные объемы отвечает предшествующей задаче. Для граничного с окружающей средой елементарного объема 1 неявная форма разностного уравнения баланса энергии имеет вид

Т₁^τ+Δτ [1+2Fo(1+Bi)] - 2Fo [ Т₂^τ+Δτ + Bi Т_l^τ+Δτ ] - Т₁^τ =0.

Для других внутренних узлов, обозначенных индексом i, применимое следующее соотношение

[1 + 2Fo] Т_i^τ+Δτ - Fo [Тt_i-1^{τ +Δτ} + Тt_i+1^τ+Δτ ] - Т_i^τ = 0.

Будем решать систему алгебраических уравнений методом обращенной матрицы, поэтому запишем все разностные уравнения баланса энергии в матричной форме: A1 ^. T = B.

Уравнение баланса энергии для граничного элементарного объема в матричной форме имеет вид

Т₁^τ+Δτ [1+2Fo(1+Bi)] - 2Fo Т₂^τ+Δτ = 2Fo Bi Т_l^τ+Δτ + Т₁^τ ,

а для всех внутренних элементарных объемов уравнения баланса записывается следующим образом

- Fo Т_i-1^{τ +Δτ} + [1 + 2Fo] Т_i^τ+Δτ - Fo Т_i+1^τ+Δτ = Т_i ^τ .

Если принять значение шагов по координатам и времени как и в предшествующей задаче, то значение чисел Фурье и Био будут Fo = 0,25 и Bі = 0,5. Ограничим область решения 10 элементарными объемами. В таком случае используем 10 разностных уравнений баланса энергии:

узел 1: 1,75 Т₁^τ+Δτ - 0,5 Т₂^τ+Δτ = 25 + Т₁^τ ;

узел 2: - 0,25 Т₁^τ+Δτ + 1,5 Т₂^τ+Δτ - 0,25 Т₃^{τ +Δτ} = Т₂^τ ;

узел 3: - 0,25 Т₂^τ+Δτ + 1,5 Т₃^τ+Δτ - 0,25 Т₄^{τ +Δτ} = Т₃^τ ;

узел 4: - 0,25 Т₃^τ+Δτ + 1,5 Т₄^τ+Δτ - 0,25 Т₅^{τ +Δτ} = Т₄^τ ;

узел 5: - 0,25 Т₄^τ+Δτ + 1,5 Т₅^τ+Δτ - 0,25 Т₆^{τ +Δτ} = Т₅^τ;

узел 6: - 0,25 Т₅^τ+Δτ + 1,5 Т₆^τ+Δτ - 0,25 Т₇^{τ +Δτ} = Т₆^τ;

узел 7: - 0,25 Т₆^τ+Δτ + 1,5 Т₇^τ+Δτ - 0,25 Т₈^{τ +Δτ} = Т₇^τ ;

узел 8: - 0,25 Т₇^τ+Δτ + 1,5 Т₈^τ+Δτ - 0,25 Т₉^{τ +Δτ} = Т₈^τ;

узел 9: - 0,25 Т₈^τ+Δτ + 1,5 Т₉^τ+Δτ - 0,25 Т₁₀^{τ +Δτ} = Т₉^τ;

узел 10: - 0,25 Т₉^τ+Δτ + 1,5 Т₁₀^τ+Δτ - 0,25 Т₁₁^{τ +Δτ} = Т₁₀^τ.

В эту систему уравнений входит значение температуры в узле 11. Температуру в узле 11 можно найти из граничного условия, которое требует, чтобы она оставалась равной начальной температуре 473 К, если время, на протяжении которого отыскивается решение, не слишком велико. Отметим, что на основании полученного решения предшествующей задачи температура в узле 5 практически не меняется в течении первой минуты. Поэтому, приняв температуру в узле 11 равной Т₀, мы не получим серьозных ошибок. Тогда, уравнение энергии для узла 10 будет иметь вид:

- 0,25 Т₉^τ+Δτ + 1,5 Т₁₀^τ+Δτ = 50 + Т₁₀^τ.

Дальнейшее решение задачи будем выполнять в пакете программ Mathcad, используя матричный метод. Для этого все коэффициенты у систем и уравнений запишем в виде двух матриц:

.

Обращенная матрица имеет вид

Для всех следующих шагов вычислений значения величин в матрице М^-1 остаются постоянными. Ниже приведенные результаты расчеты вычислений для значений времени 13,3; 26,6; 39,9; 43,2 и 56,5 с.

На основании проведеннях расчетов строим графики изменения температуры внутр. Рассматриваемого твердого тела ( рис. 5.1).

Рис. 5.1. Зависимость изменения температуры внутри твердого

тела по времени.

ЗАДАЧИ

Задача 5.1. Решить задачу 4.4. неявным методом.

Задача 5.2. Нагревается брус прямоугольного пересечения (коэффициент температуропроводносто а, м²/с; коэффициент теплопроводности λ, Вт/(м^. К)). Начальная температура бруса Т₀= 300 К. Нагревание осуществляется с одной стороны и заканчиваются в тот момент, если на противоположной стороне бруса температура поверхности достигнет значение Т, К . На торцах бруса имеет место конвективний теплообмен ( коэффициент теплоотдач α = 100 Вт/(м^{2 .} К)) с окружающей средой, которая имеет температуру Т_ж = 300 К. Нагреватель, который имеет температуру поверхности Т_н = 673 К, плотно прилегает к поверхности бруса (термическим сопротивлением воздушной прослойки между брусом и нагревателем можно пренебречь). Какое время должный работать нагреватель, чтобы температура бруса на противно-положной стороне стала равной Т, К. Исходные данные для решения задачи приведены в табл. 5.1.

Таблица 5.1.

5 10 4 9 3 8 2 7 1 6 l1 q2 l2 q3 q1 q4 α, Тж	№	λ,Вт/(м*К)	a, м2/с	T,K	l1,м	l2,м
	1	25	2,5.10-5	373	0,10	0,20
	2	30	2,0.10-5	363	0,12	0,22
	3	28	3,0.10-5	403	0,08	0,24
	4	40	2,2.10-5	383	0,10	0,18
	5	45	2,8.10-5	393	0,12	0,20
	6	42	4,0.10-5	380	0,14	0,26
	7	35	3,5.10-5	410	0,10	0,24
	8	38	2,7.10-5	405	0,12	0,18
	9	29	2,1.10-5	396	0,14	0,20