ПЗ - Численное модлирование распределения температурного поля
Введение…………………………………………………………………. 3
1. ОБЗОР ЧИСЛЕННЫХ МЕТОДОВ МОДЕЛИРОВАНИЯ В МЕХАНИКЕ……………………………………………………………...
4
1.1 Общая характеристика численных методов…………………….. 4
1.2 Погрешности в численных методах……………………………… 6
1.3 Численные методы в механике……………………………………. 7
1.4 Методы решения численных методов в термодинамике………. 8
2. АЛГОРИТМИЧЕСКИЙ АНАЛИЗ ЗАДАЧИ………………………. 12
3. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ ЗАДАЧИ……………………….. 17
3.1 Решение задачи с помощью пк MathCAD………………………… 17
3.2 Решение задачи с помощью пк ФОРТРАН………………………. 20
3.3 Решение задачи с помощью пк ANSYS…………………………… 24
Заключение………………………………………………………………. 31
Список использованных источников…………………………………. 32
Приложение А Код программы в среде MathCAD………………….. 33
Приложение Б Код программы в среде ФОРТРАН…………………. 35
////
2. АЛГОРИТМИЧЕСКИЙ АНАЛИЗ ЗАДАЧИ
Если требуется найти распределение температуры в твердом теле в случае, когда температура зависит от двух или трех пространственных координат, самый очевидный подход - попытаться получить точное решение основного уравнения. Уравнение стационарной теплопроводности в твердом теле с постоянным коэффициентом теплопроводности при отсутствии внутреннего тепловыделения имеет вид:
Это уравнение Лапласа. Уравнение Лапласа является линейным дифференциальным уравнением в частных производных. Известно несколько стандартных методов его решения. Один из них, например, метод разделения переменных, особенно полезен для решения задач теплообмена.
После того как тем или иным методом распределение температуры найдено, тепловой поток определяется с помощью закона Фурье. В двух- и трехмерных системах этот закон удобнее всего выразить в векторной форме:
где - градиент температуры (скалярной величины).
Градиент скалярной величины, такой, как температура, является вектором, т. е., согласно векторной записи закона Фурье, плотность теплового потока - это вектор. Обычно не рассматривается плотность теплового потока как вектор, поскольку она имеет размерность мощности на единицу площади, а ни одна из этих величии не является вектором. Однако удобно вообразить, что тепло «течет» в некотором направлении; поэтому величину часто называют вектором плотности теплового потока.
Вектор плотности теплового потока обладает важным геометрическим свойством, присущим градиентам: он направлен по нормали к изотерме, линии постоянной температуры, во всех точках твердого тела. Для иллюстрации этого свойства на рисунке 2 показаны несколько изотерм и типичных векторов плотности теплового потока в точках А, В и С твердого тела. Длина каждого из трех векторов плотности теплового потока пропорциональна местному градиенту температуры. Это значит, что в области тесного расположения изотерм градиент велик и плотность теплового потока также велика. В области, где расстояние между изотермами больше, плотность теплового потока соответственно меньше. На рисунке 1 плотность теплового потока в точке А больше, чем в точке В, где градиент температуры меньше.
/////
Решение: сначала нанесем на тело сетку с квадратными ячейками, как показано на рисунке 6. Пронумеруем узлы с неизвестными температурами от 1 до 16. Основная ячейка сетки представляет собой квадрат со стороной Δх = Δу = 250 мм. Температуры в узлах на 1-ой поверхности известны.
Рисунок 6 – Создание сетки
Система 16 разностных уравнений баланса энергии записывается следующим образом:
Узел 1:
(q∙∆x)/∆x+T_2-T_1=0;
Узел 2:
T_2=〖0,25T〗_1+〖0,25T〗_3+〖0,25T〗_7+250;
Узел 3:
T_3=〖0,25T〗_2+0,25T_4+〖0,25T〗_8+250;
Узел 4:
T_4=0,25T_3+〖0,25T〗_5+〖0,25T〗_9+250;
Узел 5:
(1000∙T_10)/2+T_4+(20∙h∙∆x)/k-T_5∙(2+(h∙∆x)/k)=0;
Узел 6:
(q∙∆x)/∆x+T_7-T_6=0;
Узел 7:
T_7=0,25T_2+〖0,25T〗_6+〖0,25T〗_8+〖0,25T〗_11;
Узел 8:
(T_9∙T_13)/2+T_3+T_7+(20∙h∙∆x)/k-T_8∙(3+(h∙∆x)/k)=0;
Узел 9:
(T_8∙T_10)/2+T_4+(20∙h∙∆x)/k-T_9∙(2+(h∙∆x)/k)=0;
Узел 10:
(T_5∙T_9)/2+T_4+(20∙h∙∆x)/k-T_10∙(2+(h∙∆x)/k)=0;
Узел 11:
(q∙∆x)/∆x+T_12-T_11=0;
Узел 12:
T_12=0,25T_7+〖0,25T〗_11+〖0,25T〗_13+〖0,25T〗_15;
Узел 13:
(T_8∙T_16)/2+T_12+(20∙h∙∆x)/k-T_13∙(2+(h∙∆x)/k)=0;
Узел 14:
T_14=0,0008T_11+0,0008T_15+27,955;
Узел 15:
(T_14∙T_16)/2+T_12+(20∙h∙∆x)/k-T_15∙(2+(h∙∆x)/k)=0;
Узел 16:
(T_13∙T_15)/2+(20∙h∙∆x)/k-T_16∙(1+(h∙∆x)/k)=0.
Полученную систему из 16 уравнений используем для определения неизвестный температур в программе MathCAD с помощью функции Find() (см. приложение А):
////
3.2 Решение задачи с помощью пк ФОРТРАН
Для численного решения задач на ЭВМ очень удобен итерационный метод, основанный на непосредственном определении температуры в каждом узле из разностного уравнения баланса энергии для этого узла. Например, если мы рассматриваем уравнение баланса энергии для внутреннего узла двумерного твердого тела, то получаем уравнение:
Разрешая это уравнение относительно температуры в узле 0, получаем:
Это соотношение типично для внутреннего узла в твердом теле с постоянными теплофизическими свойствами при отсутствии внутреннего тепловыделения, если применяется сетка с квадратными ячейками. Аналогичное соотношение получается для температуры в узле, расположенном на границе тела. Например, если узел 0 находится на границе, где происходит конвективный теплообмен с окружающей средой. температуру Т0 можно найти из уравнения:
Итак, температуру в каждом узле можно выразить через температуры в соседних узлах. Число полученных соотношений равно числу узлов с неизвестными температурами.
При использовании итерационного метода последовательно выполняются следующие четыре операции.
Операция 1. Выводит разностные уравнения, записав баланс энергии для каждого узла с неизвестной температурой. Из каждого уравнения выражают в явном виде температуру узла, для которого составлялся баланс энергии. Уравнения для всех внутренних узлов одинаковы по форме. Уравнения для граничных узлов будут различными в зависимости от типа граничных условий в конкретной задаче.
Операция 2. Задают ряд значений температур во всех узлах. Если задача будет решаться вручную, разумная начальная оценка всех температур позволит снизить затраты времени на вычисление истинных значений температуры в каждом узле. Если проводится численный расчет на ЭВМ, удобно принять все начальные температуры равными нулю.
Операция 3. Вычисляют новые значения температур, используя уравнения, полученные при операции 1. Как только получено новое значение какой-либо температуры, немедленно заменяют ее старое значение новым, так что новые значения температур в узлах все время вычисляют с использованием самого последнего приближения для остальных температур. Это позволяет уменьшить время сходимости решения к конечным стационарным значениям температур. Этот частный вид итерационного метода часто называют методом Гаусса—Зайдем.
/////
3.3 Решение задачи с помощью пк ANSYS
1. Создание плоскости
Preprocessor – Modeling – Create – Areas – Rectangle – By Dimensions
2. Создание материалов
Preprocessor – Material Props – Material Models/ Thermal –
Conductivity – Isotropic
Material – New Model… (Создание нового материала)
3. Создание конечного элемента
Preprocessor – Element Type – Add/ Add…/ Thermal Mass – Solid/
Quad 4 node 55
4. Разбитие плоскости на конечные элементы
Preprocessor – Meshing – MeshTool/ Set… (Выбор материала) / Mesh
(Создание сетки)
5. Приложение параметров
5.1 Задание конвекции
Thermal – Convection – On Lines
(К боковой поверхности стержня)
Вторая строка – коэффициент теплообмена
Четвертая строка – Температура внешней среды
////
Приложение Б
Код программы в среде Фортран
!Программа iter
program iter
Implicit None
real, dimension (16) :: T
real, dimension (16) :: TT
real summ, Dx, h, k
real TOLER
integer I
!Допустимое отклонение
TOLER = 1
!Назначаем на первом шаге расчёта все температуры равные 0
do I = 1, 16
T(I) = 0
enddo
!Пока не достигли допустимого отклонения делаем цикл
!summ = 20
Dx = 0.25
h = 500
k = 0.2
do while (summ > 16)
!Узел 1:
TT(1) = 800*Dx/k+T(2)
!Узел 2:
TT(2) = 0.25*T(1)+0.25*T(3)+0.25*T(7)+250
!Узел 3:
TT(3) = 0.25*T(2)+0.25*T(4)+0.25*T(8)+250
////...
Студент получает работу
После доработок преподаватель принимает работу и студент доволен
И только после этого эксперт получит свою оплату
