3.1 Постановка задачи (постоянные теплофизические свойства)

Рис. 2.7. Зависимость удельной теплоемкости от температуры в виде точечной диаграммы

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

При добавлении линии тренда на диаграмму Microsoft Office Excel можно выбрать любой из следующих шести различных типов тренда или регрессии: прямые, логарифмические, полиномиальные, степенные и экспоненциальные линии тренда, а также линии тренда с линейной фильтрацией. Тип линии тренда, который следует выбирать, определяется типом имеющихся данных.

При аппроксимации данных с помощью линии тренда значение величины достоверности аппроксимации R2 (число от 0 до 1, которое отражает близость значений линии тренда к фактическим данным) рассчитывается приложением Excel автоматически. При необходимости полученный результат можно показать на диаграмме. Линия тренда получается наиболее точной, когда ее величина достоверности аппроксимации близка к единице.

Для того, чтобы добавить линию тренда на диаграмму необходимо щелкнуть левой кнопкой мыши по любой точке графика, а затем нажав правой кнопкой мыши выбрать Добавить линию тренда. Появится окно Формат линии тренда, в котором можно выбрать параметры линии тренда. Во вкладке Построение линии тренда отмечаем Линейная, а внизу окна ставим галочку напротив строки Показывать уравнение на диаграмме и Поместить на диаграмму величину достоверности аппроксимации и нажимаем Закрыть (рис. 2.8).

Рис. 2.8. Задание параметров линии тренда

Аналогично можнос добовить линию тренда Степенную, Полиноминальную и т.д..

Итогом этих действий станет диаграмма с различными линиями тренда и уравнениями к ним (рис. 2.9):

- линейная аппроксимация

л(Т) = -0,0239 + 57,945·Т; R2 = 0,8895;

Рис. 2.9. Аппроксимация точечной диаграммы с помощью линии тренда

- квадратичная (полиномиальная) аппроксимация

л(Т) = 69,562 - 0,056·Т + 2Е-05· Т2; R2 = 0,9438;

- экспоненциальная аппроксимация

л(Т) = 62,783·е-6Е-04·Т; R2 = 0,8754;

- степенная аппроксимация

л(Т) = 888,62·Т-0,481; R2 = 0,88.

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

Но использование этих уравнений значительно усложняет решение и приводит к возникновению ошибок.

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

Строим новую точечную диаграмму. Снова выбираем Вставка >Диаграммы > Точечная с маркерами. Выделяем область для построения диаграммы, выбираем Работа с диаграммами > Конструктор > Выбрать данные > Выбор источника данных > Добавить. В окне Изменение Ряда в строке Имя ряда ячейку можно оставить пустой, тогда ей автоматически будет присвоено имя Ряд 1), для заполнения строки Значения Х выделяем температуры от 273 К до 1123 К, для заполнения строки Значения Y выделяем значения коэффициента теплопроводности, соответствующие данному температурному диапазону. Нажимаем ОК. Снова в окне Выбор источника данных выбираем Добавить. Теперь для заполнения строки Значения Х выделяем температуры от 1123 К до 1473 К, для заполнения строки Значения Y выделяем соответствующие значения коэффициента теплопроводности. Нажимаем ОК. Далее подписываем оси и перемещаем диаграмму на отдельный лист как было описано выше.

Рис. 2.10. Сплайн функция коэффициента теплопроводности

Теперь для каждого участка можно получить простое описание в виде линейной функции. Выделяя последовательно первый и второй участок диаграммы строим линейные линии тренда с выводом уравнений для них, согласно методике описанной выше. Итогом этих действий станет диаграмма с двумя линиями тренда и уравнениями к ним (рис. 2.11).

Рис. 2.11. Диаграмма с линиями тренда и уравнениями, определяющими данные линии тренда

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

Таким образом изменение коэффициента теплопроводности для стали 20 в диапазоне температур от 273 К до 1140 К будет описываться уравнением л(Т) = 63,219-0,0325•Т Вт/(м•К), а в диапазоне от 1140 К до 1473 К уравнением л(Т) = 13,634+0,011•Т Вт/(м•К).

2.2 Определение зависимости удельной теплоемкости от температуры

2.3 Определение зависимости плотности от температуры

В итоге получаем, что в диапазоне температур от 273 К до 1473 К изменение плотности для стали 20 описывается уравнением с(Т) = 7948,2-0,3158•Т кг/м3.

Рис. 2.12. Диаграмма с линией тренда и уравнением, определяющими зависимость плотности от температуры

2.4 Определение зависимости коэффициента температуропроводности от температуры

Рис. 2.13. Диаграмма с линиями тренда и уравнениями, определяющими зависимость коэффициента температуропроводности от температуры

В итоге получаем, что в диапазоне температур от 273 К до 1073 К изменение кооэффициента температуропроводнсти для стали 20 описывается уравнением а(Т) = 1,763Е-05-1,230Е-08•Т м2/с, а в диапазоне от 1073 К до 1473 К уравнением а(Т) = -9,164Е-07+4,931Е-09•Т м2/с.

3. Решение задачи нагрева металла аналитическим методом

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

3.1 Постановка задачи (постоянные теплофизические свойства)

Решение задачи нагрева металла с постоянными теплофизическими свойствами аналитическим методом в среде MathCAD

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

Сначала надо описать постоянные величины:

ф = 1000 - время расчета или нагрева пластины, с;

r = 0,1 - толщина пластины, м;

c - удельная массовая теплоемкость стали, ;

с - плотность стали, кг/м3;

л - коэффициент теплопроводности, .

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

Встроенная функция pdesolve применяется в рамках вычислительного блока, начинающегося ключевым словом Given, и пригодна для решения различных гиперболических и параболических уравнений. Она предназначена для решения одномерного уравнения (или системы уравнений) в частных производных (того, которое определит пользователь в рамках вычислительного блока Given), зависящего от времени t и пространственной координаты х, имеет целый набор различных аргументов и работает следующим образом: Pdesolve(u,x,xrange,t,trange,[xpts], [tpts]) - возвращает скалярную (для единственного исходного уравнения) или векторную (для системы уравнений) функцию двух аргументов (x,t), являющуюся решением дифференциального уравнения (или системы уравнений) в частных производных. Результирующая функция получается интерполяцией сеточной функции, вычисляемой согласно разностной схеме:

u - явно заданный вектор имен функций (без указания имен аргументов), подлежащих вычислению. Эти функции, а также граничные условия (в форме Дирихле или Неймана) должны быть определены пользователем перед применением функции pdesolve в вычислительном блоке после ключевого слова Given. Если решается не система уравнений в частных производных, а единственное уравнение, то, соответственно, вектор и должен содержать только одно имя функции и вырождается в скаляр;

х - пространственная координата (имя аргумента неизвестной функции);

xrange - пространственный интервал, т. е. вектор значений аргумента х для граничных условий. Этот вектор должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала);

t - время (имя аргумента неизвестной функции);

t range - расчетная временная область: вектор значений аргумента t, который должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала по времени);

xpts - количество пространственных точек дискретизации (может не указываться явно, в таком случае будет подобрано программой автоматически);

tpts - количество временных слоев, т. е. интервалов дискретизации по времени (также может не указываться пользователем явно).

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

Таблица Результаты решения задачи

Рис. 3.1. График зависимости температур от времени нагрева

3.2 Постановка задачи (переменные теплофизические свойства)

4. Решение задач нагрева металла в различных многоцелевых вычислительных комплексах

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

В настоящее время компьютерные технологии стремительно развиваются. Широкое применение численных методов и наличие мощных ЭВМ позволяют почти для любой практической задачи составить математическую модель и провести ее численное исследование. Именно упрощение, связанное с использованием алгебраических, а не дифференциальных уравнений, делает численные методы широко применимыми. Лучшим способом проверки точности численного метода является сравнение с точным аналитическим решением. Численное решение дает количественное выражение закономерностей, присущих математической модели. Это позволяет предугадать поведение модели. А с помощью экспериментального исследования наблюдается сама действительность.

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

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

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

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

4.1 Решение задачи нагрева в программном комплексе COMSOL Multiphysics (Femlab)

Нажмите OK. В рабочей области на оси х появится линия, названная по умолчанию I1.

Это же окно можно открыть, если щелкнуть на кнопке Line при нажатой клавише Shift.

Нажмите OK.

Нажмите кнопку Zoom Extents, для того чтобы расположить фигуру на весь экран.

ФИЗИЧЕСКИЕ УСТАНОВКИ

Свойства подобласти

Для задания физических свойств материала пластины требуется следующее.

Откройте окно Physics>Subdomain Settings… (в последней версии программы оно вызывается клавишей F8).

В Subdomain selection выберите номер пластины. В рабочей области эта область сразу выделится красным.

Введите в поля свойства для стали в соответствии с табл. 3.1.

В остальных полях оставьте нулевые значения.

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

Введите 273 как начальное значение в поле Temperature.

Нажмите OK.

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

Откройте окно Physics > Boundary Settings… (F7).

В Boundary selection выберите номер границы. В рабочей области соответствующая граница выделяется красным.

Введите граничные условия в соответствии с табл. 4.1.1.

В нашем случае можно использовать линейный решатель, так как все физические свойства и коэффициенты приняты не зависящими от температуры.

Нажмите OK.

ГЕНЕРАЦИЯ СЕТКИ

После задания всех граничных условий определите сетку командой Mesh > Initialize Mesh или кнопкой . Получится пятнадцать узлов по длине стенки. Чтобы увеличить количество узлов, надо нажать кнопку (Mesh > Refine Mesh). При необходимости её можно нажать несколько раз.

РАСЧЕТ

В меню Solve (Решать) выберите Solver Parameters (Параметры решателя).

В списке Solver (Решатель) выберите Time dependent (Зависимый от времени).

На вкладке General в первом поле Times вместо 0:0.1:1 введите 0:20:1000. Это значит, что для вывода решения будет использоваться шаг по времени 20 секунд. Этот шаг не имеет отношения к точности решения, но для визуализации он важен. В самом решателе шаг по времени достаточно маленький и зависит от скорости изменения температуры. Расчет будет длиться до 1000 с.

Нажмите OK.

Нажмите кнопку Solve (Решать). (Solve > Solve Problem). Получите распределение температуры по толщине пластины в основном окне.

ВИЗУАЛИЗАЦИЯ

После работы решателя в главном окне будет выведено распределение температур по толщине пластины.

По результатам расчета видно, что при условиях нестационарного режима температура теплового центра пластины равна 1099,9 К.

Нагрев тела при граничных условиях I рода с учётом зависимости теплофизических свойств от температуры

Рис. Изменение температуры теплового центра пластины со временем

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

Откройте окно Postprocessing > Domain Plot Parameters.

В области Plot type включите радиокнопку Line/Extrusion Plot.

В области Solutions to use (Используемые решения) раскройте список Select via и выберите пункт Interpolated times (Интерполированные значения времени).

В поле Times введите через пробелы следующие значения времен: 20 100 200 500 900 1000.

Выберите вкладку Line/Extrusion в области Plot type, отметьте радиокнопку Line Plot.

В поле Expression должно стоять T. В списке Subdomain Selection должна быть выбрана подобласть один.

По умолчанию после нажатия ОК будут выведены разноцветные сплошные линии, без «легенды» - это неудобно при печати.

Поэтому нажмите кнопку Line settings.

В открывшемся окне в списке Line color выберите Color. В списке Line style выберите Cycle. В списке Line marker выберите Cycle.

Нажав кнопку Color, можно выбрать любой цвет.

Поставьте галочку Legend.

Нажмите OK в окне Line settings.

Нажмите OK в окне Domain Plot Parameters. Откроется график.

Рис. Распределение температур в стенке для различных моментов времени

4.2 Решение задач нагрева в программном комплексе Elcut

ELCUT - это комплекс программ для инженерного моделирования электромагнитных, тепловых и механических задач методом конечных элементов.

Редактор модели позволяет достаточно быстро описать и создать двумерную модель исследуемых объектов.

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

ELCUT позволяет решать плоские и осесимметричные задачи по следующим темам:

- электростатика;

- электрическое поле переменных токов в неидеальном диэлектрике;

- растекание токов в проводящей среде;

- линейная и нелинейная магнитостатика;

- магнитное поле переменных токов (с учетом вихревых токов);

- нестационарное магнитное поле;

- линейная и нелинейная, стационарная и нестационарная теплопередача;

- линейный анализ напряженно-деформированного состояния;

- связанные задачи.

Нагрев тела при граничных условиях I рода с постоянными теплофизическими свойствами

Рассмотрим двумерную модель однослойной пластины (приведенную к одномерной задаче) толщиной 0,1 м с граничными условиями I рода с одной стороны и II рода (условия адиабаты) c другой. Температура внешней стороны пластины равна 1000 °С. Время нагрева пластины 1000 с.

Поскольку в Elcut задать начальную температуру тела нельзя, а по умолчанию нагрев тела происходит с 0 К, то задача делиться на два этапа.

1 этап. Рассматривается двумерная модель однослойной пластины толщиной 0,1 м с граничными условиями I рода с одной стороны и II рода (условия адиабаты) c других трех сторон. Температура внешней стороны пластины равна 273 К.

2 этап. Рассматривается двумерная модель однослойной пластины толщиной 0,1 м с граничными условиями I рода с одной стороны и II рода (условия адиабаты) c других трех сторон. Температура внешней стороны пластины равна 1273 К.

1 этап

ВВОД СВОЙСТВ ЗАДАЧИ

Открыть Elcut 5.1 > Student Edition > Elcut 5.1 Student.

Чтобы создать новую задачу:

1. В меню Файл выбрать пункт Создать.

2. Отметить пункт Задача ELCUT.

3. Ввести имя задачи: нагрев пластины.

4. Указать место для задачи: y:\Tevp\2-xx\Фамилия.

5. Указать свойства задачи:

Тип задачи: Температурное поле.

Класс модели: Плоская.

Расчет: Обычный.

Файлы: Геометрия: нагрев пластины. mod.

Свойства: нагрев пластины. dht.

6. Выбрать удобные единицы измерения:

Единицы длины: Метры.

Система координат: Декартовы координаты.

ОПИСАНИЕ ГЕОМЕТРИИ

Чтобы начать работу с моделью, необходимо описать геометрию.

В меню Правка выбрать пункт Геометрическая модель. Подтвердить создание новой модели, нажатием на кнопку OK.

Чтобы добавить элемент в модель:

1. В меню Правка выбрать пункт Добавить фигуру.

2. Параметры фигуры:

Фигура: Прямоугольник; ширина: w=1, высота: h=0.1; позиция координат центра окружности: x=0.5; y=0.05 (рис. 4.2.1).

Рис. Добавление фигуры

3. Для просмотра выбрать на панели инструментов кнопку Показать все (рис. 4.2.2).

Рис. Просмотр фигуры

Чтобы дать название объектам необходимо двойным нажатием левой кнопки мыши выделить объект и заполнить появившееся окно свойств.

Выделить верхнюю поверхность пластины:

Метка: поверхность 1; шаг дискретизации: Ручной - 0.04 (рис. 4.2.3).

Рис. Задание свойств выделенной поверхности

Аналогичным образом выделить поочередно нижнюю, правую и левую поверхности пластины обозначив их соответственно поверхность 2, поверхность 3, поверхность 4. Шаг дискретизации во всех случаях задать Ручной - 0.04. Затем выделить всю пластину (рис. 4.2.4):

Метка: пластина; шаг дискретизации: Ручной - 0.04.

Рис. 4.2.4. Задание свойств пластины

ПОСТРОЕНИЕ СЕТКИ

Чтобы построить сетку в меню Правка выбрать пункт Построить сетку и нажать подпункт Во всех блоках.

ЗАДАНИЕ СВОЙСТВ МАТЕРИАЛОВ

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

Рис. Построение расчетной сетки

Рис. Свойства метки блока - пластина

ЗАДАНИЕ ГРАНИЧНЫХ УСЛОВИЙ

Двойным нажатием левой кнопки мыши выбрать в окне описания задачи название метки ребра и заполнить появившееся окно свойств метки ребра - поверхность 1: