Охлаждение детали после нагрева в печи

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

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

Стальной стержень охлаждается после термической обработки в печи. Начальная температура стержня составляет 700 ℃, а температура окружающей среды – 25 ℃. Длина стержня – 0.8 [м], диаметр – 0.1 [м]. Степень черноты – 0.75.

Коэффициент теплоотдачи представлен следующей зависимостью (1):

(1)\[\alpha=1.02 \cdot (T-T_{ос})^{0.33}\]

где:

T – температура детали, ℃;

\(T_{ос}\) - температура окружающей среды, ℃.

Требуется определить время снижения температуры детали до 400 ℃.

Свойства стали:

• Плотность - 7750 [кг/м³];

• Удельная теплоемкость - 483 [Дж/(кг∙℃)];

• Масса – 48.7 [кг].

1.  Уравнение для описания переходного процесса температуры и его численное решение

Уравнение для описания моделируемой системы выглядит следующим образом (2):

(2)\[m\cdot C_p \cdot \frac{dT}{dt}=\alpha \cdot F \cdot (T_{ос}-T)+\epsilon \cdot \sigma \cdot F \cdot (T_{ос}^4-T^4) \]

где:

m – масса детали, [кг];

\(C_p\) – удельная теплоемкость материала детали, [Дж/(кг∙℃)];

\(\alpha\) – коэффициент теплоотдачи, [Вт/(м\(^2\)∙℃)];

F – площадь поверхности теплообмена, [м\(^2\)];

\(\epsilon\) – степень черноты;

\(\sigma\) – постоянная Стефана – Больцмана, [Вт/(м\(^2\)∙℃\(^4\) )].

Подставим функцию для коэффициента теплоотдачи (3):

(3)\[\frac{dT}{dt}=\frac{-1.02 \cdot F \cdot(T-T_{ос})^{1.33}+\epsilon \cdot \sigma \cdot F\cdot (T_{ос}^4-T^4)}{m\cdot C_p}\]

Решим данное уравнение численно с помощью следующего кода в JupyterLite:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import csv
from scipy.integrate import odeint
#Создаем переменные
F = 0.267
cp = 483
m = 48.7
eps = 0.75
sigm  = 5.67E-8
Tamb = 25
#Определяем функцию согласно записанному уравнению
def dT_dt(T, t):
    return (-1.02*F*((T-Tamb)**1.33) + eps * sigm *  F * ((Tamb+273.15)**4 - (T+273.15)**4))/(m*cp)
#Определяем границы графика и количество точек
ts = np.linspace(0,1500,1000)
#Начальное значение температуры
T0 = 700
#Решение
Ts = odeint(dT_dt, T0, ts)
Ts = np.array(Ts).flatten()
# Строим кривую с помощью matplotlib
plt.rcParams.update({'font.size':14})
plt.xlabel('t, с')
plt.ylabel('T, °С')
plt.plot(ts, Ts, color='green')
plt.grid(True)

В результате расчета получаем следующий график изменения температуры, рис. 1:

Рис. 1 График изменения температуры

2.  Подготовка модели

После термической обработки деталь охлаждается в воздухе. Охлаждение происходит за счет конвективного теплообмена и излучения.

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

1. Метка блока - Деталь;

2. Удельная теплоёмкость - 483 [Дж/(кг∙℃)];

3. Масса - 48.7 [кг];

4. Начальная температура - 700℃.

Важно

В REPEAT VISION десятичные числа задаются через точку.

Рис. 2 Свойства блока Тепловая масса

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

Рис. 3 Свойства блока Термическое сопротивление

Моделирование теплового потока излучением выполним одноименным блоком Излучение библиотеки Теплообмен. Добавим вышеуказанные блоки на рабочую область и зададим блоку Излучение следующие значения свойств (см. рис. 4 ниже):

1. Площадь поверхности - 0.267 [м\(^2\)];

2. Коэффициент излучения - 0.00000004253 [Вт/(м\(^2\cdot ℃^4\))].

Рис. 4 Свойства блока Излучение

Блок Термическое сопротивление моделирует конвективный теплообмен с переменным коэффициентом теплоотдачи, а блок Излучение - лучистый теплообмен.

Для визуализации нам понадобятся значения тепловых потоков для каждого вида теплообмена, добавим два блока Датчик теплового потока библиотеки Теплообмен в модель. Для подключения датчиков к блоку Тепловая масса необходимо сначала добавить на рабочую область блок Узел (Теплообмен) и установить значение его свойства Количество выходных портов - 1, рис. 5:

Рис. 5 Свойства блока Узел (Теплообмен)

Добавим Источник температуры библиотеки Теплообмен, который задает граничное условие по температуре окружающей среды (в графе Метка блока напишем Тос).Для задания значения температуры воспользуемся блоком Константа библиотеки Автоматика, аналогично добавив его на рабочую область и присвоив в его свойствах Значение константы - 25, в соответствии с постановкой задачи, рис. 6:

Рис. 6 Свойства блоков

Затем для подключения блока Источник температуры к схеме добавим блок Узел (Теплообмен), задав ему следующие значения свойств (см. рис. 7 ниже):

1. Количество входных портов - 1;

2. Количество выходных портов - 3.

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

Рис. 7 Подключение блоков к схеме

На верхний порт блока Термическое сопротивление на каждом шаге расчета подается значение термического сопротивления. Согласно документации блока, для конвекции оно равно (4):

(4)\[R=\frac{1}{\alpha \cdot F}\]

где:
α – коэффициент теплоотдачи; F – площадь поверхности теплообмена.

Учитывая, что (5):

(5)\[\alpha=1.02 \cdot(T-T_{ос})^{0.33}\]

Тогда (6):

(6)\[R=\frac{1}{1.02 \cdot(T-T_{ос})^{0.33}\cdot F}\]

Реализуем данную логику с помощью блоков автоматики :

• Добавим блок Константа библиотеки Автоматика, присвоив ему следующие значения свойств (см. рис. 8 ниже):

1. Метка блока - 1;

2. Значение константы - 1.

Рис. 8 Свойства блока Константа

• Добавим блок Деление библиотеки Автоматика.

• Добавим ещё один блок Константа библиотеки Автоматика, присвоив ему следующие значения свойств (см. рис. 9 ниже):

1. Метка блока - 1.02;

2. Значение константы - 1.02.

Рис. 9 Свойства другого блока Константа

• Добавим блок Сумматор библиотеки Автоматика, присвоив ему следующее значение свойств (см. рис. 10 ниже):

1. Метка блока - Вход 1 - T, Вход 2 - Tос;

2. Коэффициенты усиления по входным портам: 1;-1.

• Добавим блок Возведение в постоянную степень библиотеки Автоматика, присвоив ему следующие значения свойств (см. рис. 11 ниже):

1. Метка блока - 0.33;

2. Показатель степени - 0.33.

Рис. 11 Свойства блока Возведение в постоянную степень

• Добавим блок Умножитель библиотеки Автоматика, присвоив ему следующее значение свойства Количество входных портов - 3, рис. 12:

Рис. 12 Свойства блока Умножитель

• Добавим последний блок Константа библиотеки Автоматика, присвоив ему следующие значения свойств (см. рис. 13 ниже):

1. Метка блока - F=0.267;

2. Значение константы - 0.267.

Рис. 13 Свойства последнего блока Константа

Для подготовки модели соединим блоки согласно схеме, рис. 14:

Рис. 14 Соединение блоков

3.  Настройка проекта

Установим следующие параметры и значения (см. рис. 15 ниже):

1. Синхронизация с реальным временем - выключить;

2. Время моделирования - 1500 [с];

3. Шаг интегрирования – 10 [мс].

Рис. 15 Настройка параметров проекта

4.  Вывод результатов расчета

Для демонстрации процесса охлаждения и тепловых потоков добавим следующие графики (см. рис. 16 ниже):

• График изменения температуры тепловой массы;

• График теплового потока конвекцией (Датчик теплового потока 5) и излучением (Датчик теплового потока 6).

Для этого необходимо:

1. Нажать на кнопку Визуализация, расположенную на нижней панели слева;

2. Два раза нажмём на кнопку Добавить, расположенную в правом нижнем углу, для добавления двух графиков;

3. В окне первого графика необходимо нажать на панель Параметры блока, ось Yи выбрать там Тепловая масса;

4. В окне второго графика необходимо нажать на аналогичную п.3 панель Параметры блока, ось Yи выбрать там оба блока Датчик теплового потока(5 и 6).

Рис. 16 Визуализация

Сохраняем проект.

5.  Запуск расчета

Запустите расчёт кнопкой «Запустить моделирование»

Расчёт должен успешно завершиться следующим сообщением:

Расчет остановлен со стороны решателя. Достигнут лимит модельного времени.

6.  Анализ результатов расчёта

1. Мы получили график изменения температуры детали, рис. 17:

Рис. 17 График изменения температуры детали

Согласно графику, время достижения ≈1241 с.

2. Мы также получили графики изменения тепловых потоков, рис. 18:

Рис. 18 Графики изменения тепловых потоков

Важно

Излучение в данном случае вносит больший вклад в процесс снятия тепла с детали.

Сравним полученные результаты с полученным ранее решением в JupyterLite, рис. 19:

Рис. 19 График сравнения полученных результатов

Согласно графику, полученные значения соответствуют аналитическому решению с помощью SciPy.

7.  Выводы

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

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