Файл: Численные методы и математическое моделирование.pdf

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 12.01.2024

Просмотров: 320

Скачиваний: 4

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.

104
Рис. 26. Распределение температуры в трехслойной стенке.
Рис. 27. Распределение температур для различных шагов по времени.
Тление древесины начинается выше 120 – 150°С, термическое разложение древесины начинается с 250 – 350°С, а воспламенение
(горение открытым огнем) с 450°С. Эти границы можно опреде- лить в соответствующих узлах расчетной схемы. Проще всего по- строить график изменения температуры от времени для границы между первым и вторым слоями и найти заданные температуры.

105
1   2   3   4   5   6   7   8   9   10   11

Плоская задача с внутренними источниками тепла
В данной постановке задачи остается рассмотреть вопрос о тепловыделениях внутри тела. Тепловыделения могут быть ло- кальные, привязанные к конкретной координате тела, например, горячей провод или тонкий трубопровод, или распределённые, когда тепло выделяется во всем теле и в зависимости от его тем- пературы. В первой постановке решение задачи сводится к про- стому добавлению в соответствующий узел дополнительного ис- точника тепла, как показано в основной формуле (11) последним слагаемым. Реализуем её на основе решения с граничными усло- виями второго рода.
Копируем «Лист1 (2)», вводим еще один входной параметр
«Поток тепла» со значе- нием 10 и исправляем фор- мулу для тела в узле 4, до- бавляя в него поток тепла
«+$J$6» (рис. 28):
Копи- руем фор- мулу до конца рас- четной таблицы по вре- мени и по- лучаем следую- щий гра- фик (рис.
29).
Источ- ник тепла может быть задан и как функция от температуры, которая должна браться из предыдущего слоя по времени.
Рис. 28. Формула в узле 4.
Рис. 29. Распределение температуры с локальным источником тепла в теле.

106
Чтобы построить распределение температуры с распределен- ными источниками тепла требуется найти функциональную зави- симость выделения тепла от температуры ячейки и запаса тепла в ней, который ограничен, например, для химической реакции или поглощение тепла с определенного значения температуры и в за- данном от объема количестве для плавления материала тела. Дан- ные зависимости можно решаться с использованием дополни- тельных функций, которые учитывают выделенное или погло- щённое тепло в каждой ячейке. Для этих целей надо создать еще один набор узлов, равный расчетной схеме и поместить туда име- ющийся запас тепла, который должен быть выдан или отобран из данной ячейки. В расчетные формулы для тела добавить условие, которое с определенной температуры начинает расходовать или накапливать данной тепло пока все оно не будет израсходовано.
Данная задача проще решается с использованием приемов про- граммирования.
В самом простом случае данную задачу можно решить на ос- нове предыдущего решения, в котором поток тепла вводим во все узлы тела. Задав следующие исходные данные (рис.30) получаем следующее распределение температуры по телу (рис. 31).
Рис. 30. Исходные данные для модели с распределенными источниками тепла внутри тела.
Осесимметричные модели
Рассмотрим теперь осесимметричные задачи, которые учиты- вают форму тела в виде цилиндра или сферы. Для этих целей ис- пользуем полную конечно-разностную формулу (16). Наличие оси симметрии в модели требует использования еще одного гра- ничного условия, в котором передача тепла в точке симметрии равна нулю, то есть производная от температуры по координате


107 равна нулю и, следовательно, значения температур на оси и в со- седнем узле равны друг другу.
Рис. 31. Распределение температуры с распределенными источниками тепла.
Построим решение на основе самой первой задачи с плоским телом и двумя граничными условиями первого рода. Делаем ко- пию первого листа и вносим на нем исправления в исходных дан- ных (рис. 32), исправляя заголовок и добавляя новый параметр
«Вид симметрии:».
Рис. 32. Исходные данные для решения осесимметричных задачю
Исправляем формулу в узлах тела, начиная с первого узла по координате и времени, где формула должна иметь следующий вид:
=D11+(C11-2*D11+E11)*$F$6*$C$6/$F$5^2+
+(E11-C11)*$C$6*$F$6*$J$6/(2*$F$5) где: первая строка, какая она была в первом расчете; вторая строка должна быть дописана.

108
Затем формулу копируем по всем узлам по координате тела и потом растягиваем вниз до конечного времени.
Исправляем формулу на левой границе тела, её вид для нуле- вого узла по координате и времени показан далее: =D11, где: D11
– ссылка на ячейку в соседнем узле по координате. Растягиваем формулу до последнего узла по времени и получаем графики для разных видов симметрии в задаче, когда тело с начальной темпе- ратурой 20°С опускается в горячую плотную среду (жидкость, расплавленный металл и т.п.) с температурой 100°С (рис. 33) k=0 плоскость с осью симметрии в виде плоскость; k=1 цилиндр с осью симметрии в виде линии по его оси;

109 k=2 сфера с осью сим- метрии в виде точки в её центре;
Рис. 33. Влияние форму тела на распределение температуры по нему.
Решение осесимметричных задач с различными граничными условиями ничем не отличаются от полученных ранее решений.
Следует отметить, что чем сложнее задача, тем труднее получить её решение, при рассмотренных выше решениях будет возникать проблема с подбором шага по времени, который будет стремиться к нулю, делая решение очень медленным. Для уменьшения влия- ния шага по времени на расчет надо переходить на неявный тра- фарет.
Решение задач с использованием неявного трафарета
Неявный трафарет не дает однозначного решения в каждом узле сетки, а требует на каждом шаге решать систему конечно- разностных уравнений в которой число неизвестных равно числу уравнений.
Рассмотрим данное решение на примере первой задачи (плос- кое тело с граничными условиями первого рода, которое описы- вается системой уравнений (22). Схема расчета с неявным трафа- ретом должна выглядеть следующим образом (рис. 34), предпола- гая, что все конечно-разностные преобразования строятся для но- вого слоя с неизвестными температурами. В результате для каж- дого уравнения имеем три слагаемых с температурами слева в


110 центре, и справа по трафарету, а также узел на известном слое ниже, который будет свободным членом в полученных уравне- ниях. Получаем систему конечно-разностных уравнений объемом
n-1 узлов (рис. 35). или для неявного трафарета:




































x
T
T
x
k
a
x
T
x
T
T
T
a
x
T
j
i
j
i
i
j
i
j
i
j
i
2
;
2 1
1 1
1 2
1 1
1 1
1 2
2
(14)
При использовании неявного трафарета с использование фор- мул (12, 14) получаем другое выражение:


































x
T
T
k
a
x
T
T
T
a
T
T
j
i
j
i
j
i
j
i
j
i
j
i
j
i
2 2
1 1
1 1
2 1
1 1
1 1
1

(17) в котором имеются три неизвестные величины (
1 1
1 1
1
,
,





j
i
j
i
j
i
T
T
T
) и для решения данной задачи надо использовать процедуры реше- ния систем уравнений. Для удобства решения данных систем раз- делим все неизвестные и получим коэффициенты для каждого не- известного:
j
i
j
i
j
i
j
i
T
x
k
a
x
a
T
x
a
T
x
k
a
x
a
T




















































2 2
1 2
2 1
1 2
1 2
1 1





(18)
Рис. 34. Схема неявного трафарета
Рис. 35. Система уравнений

111
В матричной форме данная система уравнений может быть по- лучена при записи уравнения (18) для всех узлов, начиная с 1.





























































x
k
a
x
a
а
k
x
a
а
x
k
a
x
a
а
T
а
T
а
T
а
T
k
k
k
k
k
k
j
i
k
k
j
i
k
k
j
i
k
k
j
i
2 9
1
при
2 1
2 2
1
,
2
,
2 1
,
,
1 1
,
1 1
,
1 1






(30)
В матричной форме данная система уравнений может быть представлена в виде матрицы коэффициентов (А) и вектора пра- вых частей уравнения (В), так как на границах температуры нам известны, эти данные переносятся в вектор В:






















































пг
n
n
j
n
j
n
j
лг
j
n
n
n
n
n
n
n
n
n
n
T
a
T
T
T
T
a
T
a
a
a
a
a
a
a
a
a
a
,
1 1
2 2
0
,
1 1
1
,
1 2
,
1 1
,
2 2
,
2 2
,
2 3
,
2 2
,
2 1
,
2 2
,
1 1
,
1 0
0


(31)
Из матрицы видно, что она сильно разряженная и имеет стро- гую трех диагональную схему расположения отличных от нуля элементов. Для данных матриц разработана своя схема решения задачи, которая называется метод прогонки. Суть решения дан- ной матрицы заключается в исключении поддиагонального члена матрицы приравниванием его к нулю путем следующего преоб- разования:
























































j
лг
j
лг
j
j
лг
j
T
a
T
a
T
a
T
a
T
a
a
a
a
a
a
a
a
T
T
a
T
a
a
a
a
a
a
a
2 1
,
1 0
,
1 1
1
,
2 0
,
1 1
1
,
1 3
,
2 1
,
1 2
,
2 2
,
1 1
,
2 2
,
1 1
,
1 2
0
,
1 1
3
,
2 2
,
2 1
,
2 2
,
1 1
,
1 1
,
1 1
,
2 0
0 0
0 0
0 0


112
Умножаем первую строку матрицы на первый элемент второй строки и вторую строку матрицы на первый элемент первой строки, потом вычитаем одну строку из другой и помещаем ре- зультат во вторую строку матрицы. Процедуру повторяем до по- следней строки матрицы.
В результате преобразования в последней строке матрицы остается один ненулевой элемент, что позволяет вычислить зна- чения температуры для последнего узла. Далее подставляем найденное значение температуры в предыдущую строку матрицы и находим следующую температуру в предыдущем узле. Повто- ряем вычисления пока не достигнем первого узла. Текст про- граммы на VBA представлена ниже с необходимыми коммента- риями:
Public Sub Progon(): Rem Имя макроса
Так как программа-макрос на первом этапе требует считыва- ния с листа начальных параметров и потом заполнения таблицы расчетными данными устанавливается заранее выбранная началь- ная ячейка
Cells(1, 1)
относительно которой организуется обраще- ние к остальным ячейкам таблицы с помощью команды
Active-
Cell.Offset(dR, dC)
, где dR
– смещение по строкам; dC
– смещение по столбца:
Cells(1, 1).Select
Tob = ActiveCell.Offset(2, 1): Rem T нач
Tsr = ActiveCell.Offset(3, 1): Rem T среды dT = ActiveCell.Offset(3, 3): Rem dTau n = ActiveCell.Offset(1, 3): Rem Число узлов nn = ActiveCell.Offset(4, 1): Rem начальный узел расчета nk = ActiveCell.Offset(4, 3): Rem конечный узел расчета
Tend = ActiveCell.Offset(5, 1): Rem Конечное время
Описываем массивы для расчета в программе:
ReDim A(1 To 3, n), B(n), T(n)
Задаемся начальной строкой для печати таблицы расчетных данных и начальным расчетным временем, далее заполняем первую строку расчетной таблицы: ik = 17: tau = 0

113
ActiveCell.Offset(ik, 0) = tau
For i = 0 To n - 1
ActiveCell.Offset(ik, i + 1) = Tob
Next i
ActiveCell.Offset(ik, n + 1) = Tsr
Основный цикл расчета в программе:
Do
В строку статуса записываем текущее расчетное время:
Application.StatusBar = "Расчетное время " & Format(tau, "#0.0")
Заполнение вектора B на рабочем листе известными значени- ями температур:
For i = 0 To n
ActiveCell.Offset(7, i + 1) = ActiveCell.Offset(ik, i + 1)
Next i
Заполнение расчетных массивов:
For i = 0 To n - 1
For j = 1 To 3
A(j, i) = ActiveCell.Offset(7 + j, i + 1)
Next j
B(i) = ActiveCell.Offset(11, i + 1)
Next i
Исключение поддиагональных элементов матрицы А и коррек- тировка значений вектора В:
For i = 0 To n - 2 a1 = A(1, i + 1) / A(2, i)
A(1, i + 1) = A(1, i + 1) - a1 * A(2, i)
A(2, i + 1) = A(2, i + 1) - a1 * A(3, i)
B(i + 1) = B(i + 1) - a1 * B(i)
Next i
Вычисление температур:
T(n) = Tsr
T(n - 1) = B(n - 1) / A(2, n - 1)
For i = n - 2 To 0 Step -1
T(i) = (B(i) - A(3, i) * T(i + 1)) / A(2, i)
Next i
Заполнение очередной строки в таблице расчетных данных:

114 ik = ik + 1: tau = tau + dT
ActiveCell.Offset(ik, 0) = tau
For i = 0 To n
ActiveCell.Offset(ik, i + 1) = T(i)
Next i
Проверка условия на завершение расчета:
Loop Until tau > Tend
Application.StatusBar = "Готово"
End Sub
Для правильной работы программы необходима строго соблю- даемая структура листа с заполненными данными (рис. 36)
Рис. 36. Вид листа с данными для применения метода прогонки.
В результате расчета можно получить решения на гораздо бо- лее длительные периоды времени. В данном примере показан нагрев пластины при наличии конвективного теплообмена с хо- лодной средой (рис. 37).