Материал: Математическое моделирование температурных полей при теплотехнических процессах

Вышеописанные уравнения описывают алгоритм структурно-матричного метода формирования математической модели тепловой системы. При этом необходимо учитывать граничные условия.

Граничные условия I рода.

Задаются значения температур граничных узлов Т*в1 и Тв2 В этом случае на основе вышеописанных уравнений непосредственно получаем математическую модель теплового объекта в виде системы уравнений. При этом потенциалы диссипативных элементов Фдк определяются по компонентным уравнениям, а последнее уравнение позволяет определить значения тепловых потоков Ф*в1 Ф*в2 на левой в правой граничных поверхностях, необходимые для поддержания заданных значений температур граничных узлов Т*в1 и Г*в2 -Тепловые потоки Ф*в1 и Ф*в2 генерируются источниками тепловой анергии внешней среды (или представляют собой стоки - в случае отвода тепла от твердого тела). Но при граничных условиях I рода потенциалы источников внешних воздействий неизвестны, поэтому в уравнении принимается Ф*вl= 0.

Для решения полученной системы дифференциальных уравнений необходимо задать начальные условия Ti0,

Таким образом, при формировании математической модели тепловой системы с граничными условиями I рода используется динамическая модель с источниками внешних воздействий типа потока, характеризуемыми температурами граничных поверхностей Т*в1 и Т*в2.

Граничные условия II рода. Задается градиент температуры на граничных поверхностях. В одномерном случае необходимо задать значения дТ/дх при х = 0 и при х =b , где b - размер твердого тела вдоль оси х . Их можно легко свести к граничным условиям III рода. Однако при равенстве нулю градиента температуры на граничных поверхностях целесообразно непосредственно использовать граничные условия II рода. Этот случай соответствует теплоизолированной граничной поверхности.

Предположим, что пятый слой теплового объекта на рис.2 выполнен из теплоизоляционного материала, характеризуемого малым значением коэффициента теплопроводности λ5. При этом μТ5→0, а Тв2≈Тв4. Следовательно, теплообмен с окружающей средой через правую граничную поверхность будет незначительным и им можно пренебречь. В орграфе при этом граничный узел 2* может быть исключен (рис. 3, б). Составим матрицу инциденций для твердого тела с теплоизолированной правой границей и с заданными граничными условиями I рода на левой границе. Матрица Л приведена в табл.2

Таблица 2

Математическая модель объекта в этом случае будет такой же, как и в предыдущем, но только в ней ФД5=0.

Граничные условия III рода.

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

Граничные условия в первом случае имеют вид

(16) где Фвх - проекция теплового потока на ось х, АВ - площадь поверхности граничного теплообмена.

Составим граничные условия для х=0 и х=b (рис. 2), заменив при этом частные производные дT/дx отношениями конечных разностей

(17), (18)

где Т*В1, Т*В2 - температуры граничных узлов, обусловленные тепловыми потоками ФB1 и ФВ1 на граничных поверхностях.

Если тепловые потоки ФB1 и ФВ2 на левой и правой граничных поверхностях твердого тела направлены от внешней среды к телу, т. е. тепловая энергия подводится к нему с двух сторон, то проекции тепловых потоков на ось х будут равны:

Фв1х=Фв1; ФВ2Х = - ФВ2. В этом случае получаем

(19), (20)

Легко заметить, что в этом случае в системе уравнений потенциалы диссипативных элементов первого и последнего слоев твердого тела Фд1 и ФД5 должны быть заменены на потенциалы источников внешних воздействий ФВ1 и ФВ2 Эта замена на орграфе отражается исключением граничных узлов (рис.3, в). Матрица инциденций, соответствующая данному орграфу, приведена в табл.3.

Выражения используются для определения температуры граничных узлов

(21)

Таблица 3

.2 Метод конечных разностей

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

Сущность задачи заключается в покрытии области R сеткой и замене непрерывного множества независимых переменных X в области R конечным множеством точек Хк, являющихся узлами сетки. Сетка может быть прямоугольной, косоугольной, с постоянными или переменными расстояниями между узлами (шагами сетки). Наиболее часто используют сетку с постоянными расстояниями между узлами, которая называется регулярной. Алгебраизация дифференциальных уравнений в методе конечных разностей осуществляется за счет аппроксимации производных.

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

Метод прогонки является модификацией метода Гаусса для частного случая разреженных систем - системы уравнений с трехдиагональной матрицей. Он состоит из двух этапов - прямой и обратной прогонки. Прямая прогонка состоит в том, что каждое неизвестное xi выражается через xi+1 с помощью прогоночных коэффициентов Ai, Bi:

xi=Aixi+1+Bi i=1,2,3…,n-1.

Обратная прогонка состоит в последовательном вычислении неизвестных xi. Сначала находим xn , затем последовательно вычисляем все неизвестные xn-1, xn-2 , xi.

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

.3 Метод конечных элементов

При решении задач, связанных с определением напряжений и деформаций в элементах конструкций технических объектов, чаще всего используется метод конечных элементов. Дискретизация такого объекта осуществляется с использованием конечных элементов, a алгебраизация задачи заключается в выборе аппроксимирующих функций для каждого конечного элемента. Малые размеры конечных элементов позволяют использовать простые аппроксимирующие функции , причем одного и того же типа для всех конечных элементов определенной формы. Обычно в качестве  для отдельного конечного элемента применяют полиномы не выше третьей степени, например, в одномерном случае

(22)

Функции v(Х) в методе конечных элементов представляют в форме

 (23)

где коэффициенты qi имеют вполне определенный физический смысл - это значения аппроксимирующей функции в узловых точках: - функции, которые называются координатными (или функциями формы); r - число узловых точек в конечном элементе.

В общем случае аппроксимации вектора  в m-мерном пространстве выражение (23) принимает вид

где - вектор размерности ;  - вектор размерности ; N - интерполяционная матрица порядка , элементами которой являются координатные функции.

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

(24)

Минимизируя потенциальную энергию Еп, находят вектор перемещений .

Связь деформации  с перемещениями ui можно выразить в матричной форме

или более лаконично

,(25)

где S - матрица-оператор дифференцирования.

Таким образом

,где .

Обозначим  (26)

Матрица К является матрицей жесткости. Тогда. В соответствии с принципом Лагранжа дифференцируем Еп по вектору , и приравниваем к нулю. В результате получаем систему алгебраических уравнений,где  вектор правых частей, который называется вектором нагрузок.

Матрица жесткости К всей исследуемой детали составляется из матриц жесткости Кij отдельных КЭ. Матрицы Кij содержат информацию о конфигурации и упругих свойствах материала конечных элементов и подсчитываются по формуле (26), в которой при этом под R понимается подобласти, относящиеся к рассматриваемому КЭ.

. Алгоритмический анализ задачи

.1 Постановка задачи

В качестве исходной математической модели используем уравнение энергии, в котором исключим оператор конвективного переноса и учтем одномерность задачи, t = t(x, τ). В результате получим

.(27)

температурный поле моделирование теплотехнический

На левом и правом торцах объекта происходит тепловое взаимодействие со средой, и здесь необходимо задать граничные условия. Универсальным способом описать самые различные ситуации будет применение граничных условий третьего рода на левом (x = 0) и правом (x = L) торцах объекта:

;,(28), (29)

где α и tf - коэффициенты теплоотдачи и температуры окружающей среды на торцах стержня.

В этих соотношениях приравнены значения плотности теплового потока,

 поступающего из окружающей среды и вычисленного по уравнению Ньютона-Рихмана (правые части)

 и отводимого внутрь тела посредством теплопроводности и вычисленного по закону Фурье (левые части).

Следует выделить, что такое равенство справедливо при отсутствии фазовых превращений на поверхности раздела.

На правом торце задана температура, которая пульсирует около определенного значения Tf2 c амплитудой Ampl, а на левом торце - постоянная температура Tf1:

tf1=const=Tf2(1+Ampl´sin(vt))

Задача решается методом конечных разностей, о котором было сказано ранее.

 (29)

 (30)

 (31)

 (32)

 (33)

В (5) использована неявная разностная схема, которая характеризуется абсолютной устойчивостью. Выполним аппроксимацию граничных условий:

 (34)

 (35)

При i =(1..n-1):

(36)

Структура получившейся линейной системы уравнений с коэффициента-ми A, B, C, D является трехдиагональной.

Из (34) следует

Из (35) следует

При i =(1..n-1)

(37)

Система (9) позволяет найти величины t1, t2… в момент времени t=1, то есть позволяет найти распределение температуры t во внутренних узлах стержня в начальный момент времени. Далее следует многократно решать систему (8) для нахождения распределения температуры в последующие моменты времени. Для этого удобнее всего использовать метод прогонки, о котором упоминалось ранее. Для расчета в MathCAD нужно использовать программный фрагмент, при этом следует выполнить следующие этапы расчета:

Нужно задать функцию пользователя, служащую для решения системы методом прогонки

Реализовать алгоритм циклического вызова для каждого момента времени

Организовать вывод результата в матрицу, столбцы которой являются векторами решений в последний момент времени.

.2 Исходные данные

В разрабатываемом курсовом проекте используются следующие данные

Материал и размеры стержня

l = 31 (Вт/м К) Теплопроводность

с = 810 (Дж/кг К) Удельная теплоемкость

r = 7850 (кг/м ) Плотность= 0.3 (м) Длина стержня

Параметры теплового воздействия= 0.310 Амплитуда= 3 Частота изменения температуры

Шаги пространственно-временной сетки

Dx = 0.005 (м) Шаг по координате перемещения

Dt = 0.15 (с) Шаг по координате времени

Температуры внешней среды на концах стержня= 285 (К) Температура в начале стержня= 720 (К) Температура в конце стержня