-2 (34). 2005
/109
ИТЕЙНОЕК ПРОИЗВОДСТВО
The mathematical model for numerical analysis of scaling of the ingots, moving with constant speed in the furnaces with step-heating and given quantity of temperature zones, is offered. The results of computer calculations of temperature and scale layer forming at movement of ingot with constant speed in furnaces with five temperature zones are presented.
A. H. ЧИЧКО, А. С. БОРОЗДИН, БИТУ
УДК 669.27:519
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ОБРАЗОВАНИЯ ОКАЛИНЫ С ИСПОЛЬЗОВАНИЕМ СИММЕТРИЧНОЙ НЕЯВНОЙ СХЕМЫ В ДВИЖУЩЕМСЯ НАГРЕВАЕМОМ СЛИТКЕ ИЗ ХРОМИСТОЙ СТАЛИ
Цель настоящей работы — разработка математической модели и метода для численного моделирования процесса высокотемпературного окисления заготовок из хромистых сталей в условиях ступенчатого нагрева и применение разработанного метода для компьютерного моделирования процесса нагрева фрагмента слитка и роста на нем окалины в многозонных шаговых печах.
Комплексные исследования кинетики ока-линообразования для железа и сталей, выполненные в Ленинградском университете [1], позволили установить, что высокотемпературное окисление хромистых сталей в атмосфере влажного воздуха в целом подчиняется параболическому закону и может быть описано следующим уравнением:
& = Ах, (1)
где 5 — прирост массы окалины на единицу площади поверхности, кг/м2; х - время окисления, с; К - константа скорости окисления. Константа скорости окисления К при постоянном давлении кислорода в нагревающей среде выражается по закону Аррениуса следующей зависимостью [2—4]:
(2)
где L
константа; Q — энергия активации
(обычно в кал/моль); К — универсальная газовая постоянная; Т — абсолютная температура, К.
Исключение составляет начальный период окисления, когда окисление стали происходит с аномально замедленной скоростью, так называемый индукционный период. Длина индукционного периода быстро уменьшается с повышением температуры. Суммарный привес окалины на единицу поверхности во время индукционного периода для разных температур приблизительно один и тот же и составляет величину порядка 10-15 кг/м2.
Учитывая, что температура поверхности заготовки при нагревании в различных пространственных точках различна как по величине, так и по скорости изменения, для численного расчета температурного поля, входящего в уравнение (2), можно использовать дифференциальное уравнение теплопроводности [5]:
c{T)PiT)^ = f Эх дх
ду
ду
д_ ' dz (
dz
(3)
х,у,г,хе О.,
где с(Т) — функция теплоемкости, Дж/(кг-К); р(7) — функция плотности, кг/м3; ЦТ) -функция теплопроводности, Вт/(м-К); П (0=х<Аг, 0=у<У, 0=г<Д 0=х<0 — пространственно-временная область расчета.
В качестве начальных и граничных условий, описывающих взаимодействие между нагревающей средой и поверхностью объекта, равномерно перемещающегося в печи с различными тепловыми зонами, использовали следующую систему уравнений [5]:
Т0, т = 0,
T(x,y,z,T) =
71,0 < х < —,
< X < -
1)
U + L0
(4)
п — 1 п
5л 2л
т
к=1
<х<
к=1
110 /Si
xrrrf: гг гсшмггггта
(3«), 2005-
где Т0 - начальная температура заготовки, °С; п — число тепловых зон в рабочем пространстве печи; Т. (/=1, 2, п) - значение температуры в /-й тепловой зоне, °С; Ь. (/—1, 2, ..., п) — протяженность /-й тепловой зоны, м; -о — скорость перемещения заготовки вдоль рабочего пространства печи, м/с.
Уравнения (1)-(3) вместе с начальными и граничными условиями (4) представляют собой математическую модель для расчета процесса ока-линообразования в трехмерном объекте, помещенном в изменяющееся температурное поле.
Компьютерная модель построена на основе синтеза идей теории клеточных автоматов и конечно-разностных методов. Численное решение уравнения теплопроводности, которое является основным уравнением математической модели, проводили на основе симметричной неявной конечно-разностной схемы. Для этого была введена пространственно-временная сетка следующего вида:
х = пкх (п = 0, 1, ..., Л0, ут = тку (т = 0, 1, ..., М), г, = /А (/ = 0, 1, ..., ¿),
хк = кИх (к = 0, 1, ..., К), где к=ХЩ /у=У/М; к=г/Ц А=//£
Точное решение тепловой задачи Т в узле (хп, ут, I,, тк) заменяли значением сеточной функции и:
ипт1 =Т(ХП, Ут,2ьХк) •
На восьмиточечном шаблоне была составлена симметричная неявная схема, имеющая погреш-
( 2 2 2 2 I
кх + кх+ку+1гг) по экономичной схеме [7].
На каждой итерации компьютерного моделирования после расчета температурных полей слитка и рабочего пространства печи для поверхностного слоя фрагмента слитка проводили расчет слоя образовавшейся окалины. Расчет прироста окалины для моментов времени параболического закона роста окалины осуществляли подстановкой температуры поверхности в уравнение (1). Причем константу скорости окисления К рассчитывали по уравнению [3]:
1 ^ 49180 пи
\gK =--+ 9,16
2,3 RT
(5)
Для компьютерных расчетов прироста окалины для индукционного периода для /-й итерации моделирования предлагается следующее уравнение:
ind
si=s<-i + _i„d
(6)
период; — продолжительность индукционного нагрева для данной температуры, с.
Уравнение (7) построено на предположении что в индукционный период рост окалины подчиняется линейному закону. Для вычислительных экспериментов принималось равным 1 мг/см2, а продолжительности индукционных периодов представлены в табл. 1.
Таблица 1. Продолжительность индукционного периода для различных температур
г, ° с Продолжительность индукционного периода, ч
880 15,5
900 7,5
920 4,0
930 4,0
950 4,0
1000 3
1070 0,25
где 5. — количество окалины, образующееся на /-й итерации расчета, кг/м2; 51пс1 - общее количество окалины, образующейся в индукционный
Для промежуточных значений температур, отсутствующих в таблице, значения индукционных периодов рассчитывали при помощи интерполяции полиномом Лагранжа.
Для проведения расчетов на ПЭВМ данная вычислительная схема была реализована в виде набора модулей математического ядра программы для моделирования теплофизических процессов. Программное обеспечение написано в визуальной среде быстрой разработки программ Borland Delphi 6.0 на языке Object Pascal и предназначено для работы в 32-разрядных версиях ОС Windows. Разработанные модули позволяют проводить моделирование процессов нагрева и расчет слоя окислившегося металла на поверхности детали произвольной конфигурации для печей, имеющих различные по протяженности температурные зоны. При проведении вычислительных экспериментов имеется возможность варьировать как скоростью движения детали вдоль рабочего пространства печи, так и конфигурацией рабочего пространства печи. Помимо модулей расчета, программа включает в себя модули, отвечающие за 2D/3D-визуализацию рассчитанных данных и взаимодействие с пользователем.
Согласно целям поставленного исследования, была построена компьютерная модель промышленной печи, состоящей из пяти тепловых зон, с температурами 7j, Tv Tv Т4 и Т5 для каждой зоны соответственно. Протяженность каждой температурной зоны по L= Ь2=L3=L4=Ls=L=4 и. Таким образом, общая длина рабочего пространства моделируемой печи составила 20 м. Пространство печи принималось заполненным воздухом со следующими теплофизическими характеристиками: >.=0,034 Вт/(м К), с=1009 ДжДкг К), р=1,29 кг/м3 и температурой, зависящей от конкретной тепловой зоны. В качестве материала для слитка была выбрана сталь марки 40Х с тепло-физическими характеристиками, являющимися
л ГГТТгГ, ГГ гггг?
---—-2 (34). 2005
функциями от температуры слитка. В табл. 2 водности материала слитка, использованные при приведены значения теплоемкости и теплопро- моделировании тепловых полей.
Таблица 2. Теплоемкость и теплопроводность для стали 40Х
111
т,° с 0 100 200 300 400 500 600 700
с, Дж/(кг К) 496 508 529 563 592 622 634 664
Х(Т), Вт/(м К) 41,0 40,0 38,0 36,0 34,0 33,0 31,0 30,0
Т,° С 800 900 1000 1100 1200 1300 1400 1500
с, Дж/(кг К) 684 694 788 880 972 1064 1160 1120
Х(Т), Вт/(м К) 27,0 26,0 25,5 25,0 24,7 24,4 24,2 23,0
Компьютерное моделирование проводили для фрагмента слитка размерами 250x250x300 мм. Начальная температура слитка принималась равной Г0=20 °С. При численном моделировании процесса нагрева данный слиток двигался вдоль печи с равномерной скоростью г>=0,004 м/с. Оси координат выбирались так, что ось X совпадала с направлением движения слитка. Компьютерные
расчеты проводили при пространственном шаге к=к=И=к=0,02 м и временном шаге Ах=5 с. Временной шаг определялся как время, необходимое слитку, чтобы при заданной скорости движения сместиться на один шаг по оси X На рис. 1, а схематично изображены рабочее пространство моделируемой печи (вид сверху), местоположение температурных зон и первоначальное положение слитка.
Слиток ..... Зона 1 ' /гг ■ ' £ 1 Зона 2 12 Зона 3 "Тг ^ »'У •'.< . л * - . ' • % * >* ■ ; г^ П / > Ч 4 I % Зона 5
¿2 ¿4 ¿5
£
Рис. 1. Схема рабочего пространства печи и пространственная модель слитка: а — расположение температурных зон и слитка;
б - расположение контрольных точек на поверхности слитка
На поверхности слитка были выбраны три контрольные точки. В процессе моделирования для данных точек вычисленные значения температуры и окалины сохранялись в текстовый файл для последующего анализа. На рис. 1, б показано поперечное сечение слитка и расположение выбранных контрольных точек.
При моделировании проводили варьирование значениями температур в зонах рабочего пространства печи. В табл. 3 приведены различные
тепловые режимы температур, для которых проводили вычислительные эксперименты.
Моделирование выполняли для первых 3400 с процесса нагрева. За данный временной промежуток при указанной скорости движения заготовки слиток проходит вдоль всего рабочего пространства печи.
В табл. 4 приведены результаты проведенных вычислительных экспериментов по определению толщины окалины для контрольных точек слитка в различные моменты времени.
Таблица 3. Тепловые режимы для ступенчатого нагрева слитков в печи
Номер теплового режима т2,°с 7з, °С Г4, °С Г5,°С
1 1200 1200 1200 1200 1200
2 300 300 500 800 1200
3 200 600 950 1100 1200
/дгттг^ гг гсшм/^ггпт?
I 2 (34), 2005-
Таблица 4. Изменение толщины окалины (мг/см2) в контрольных точках А, В, С слитка для
различных тепловых режимов нагрева
фс Тепловой режим 1 Тепловой режим 2 Тепловой режим 3
А В С А в с А В С
0 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000
500 0,021 11,339 0,018 0,000 0,000 0,000 0,000 0,000 0,000
1000 12,299 23,011 14,131 0,000 0,000 0,000 0,000 0,000 0,000
1500 22,232 32,179 23,943 0,000 0,000 0,000 0,000 0,000 0,000
2000 30,788 39,856 32,352 0,001 0,001 0,001 0,001 0,001 0,001
2500 38,294 46,522 39,706 0,001 0,001 0,001 0,003 0,021 0,004
3000 44,963 52,457 46,240 0,001 0,001 0,001 0,029 0,056 0,033
3400 49,811 56,799 50,994 0,001 0,001 0,001 0,049 7,867 0,053
Температурно-временные диаграммы, описывающие динамику роста температуры для контрольных точек слитка, представлены на рис. 2. Анализ полученных результатов показывает, что для первого режима нагрева слитка индукционный период заканчивается достаточно быстро и рост окалины идет в основном по параболическому закону. Для второго режима нагрева максимальная температура в контрольных точка слитка не достигла 800 °С, а процесс окисления не вышел за пределы индукционного. Интересно выглядит третий режим нагрева, при котором все контрольные точки к концу нагрева имеют температуру свыше 1000 °С и только одна из трех точек перешла в период с параболическим законом роста окалины. Необходимо также отметить, что для всех режимов в точке В нагрев и, как следствие, рост окалины происходит наиболее быстро.
Представленные вычислительные эксперименты наглядно демонстрируют возможности компьютерного моделирования для анализа процессов высокотемпературного окисления, возникающих при нагреве слитков с динамическими характеристиками нагревающей среды. Предложенная математическая модель образования окалины в трехмерном слитке, движущемся с заданной скоростью в неоднородном температурном поле, позволяет учитывать конфигурацию слитка и скорость его движения. Численные расчетные зависимости по кривым нагрева слитка данного размера согласуются с известными представлениями об образовании окалины в стальных заготовках.
Литература
1. Уч. зап. Ленингр. ун-та, 1954. Вып. 14. №175.
2. Кофстад П. Высокотемпературное окисление металлов. М.: Мир, 1969.
3. Ку баше веки й О., Гопкинс Б. Окисление металлов и сплавов. М.: Металлургия, 1965.
4. Окисление металлов. Т. 1. Теоретические основы / Под ред. Ж.М. Бернара: Металлургия, 1967.
О 500 1000 1500 2000 2500 3000 3500
Рис. 2. Зависимость температуры от времени для контрольных точек слитка при его движении по температурным зонам для различных технологических вариантов: а - Т = Т2=Т = Т = Т=\200 °С; б - Т=300 °С, Г2=300 °С, Г3=500 °С, Г4=800 °С, Г5=1200 °С; в - Т=200 °С, Г2=600 °С, Г3=950 °С, Г4=1100 °С, Т5=1200 °С
5. Чичко А. Н., Бороздин А. С. Численное моделирование процесса нагрева движущегося слитка // Литье и металлургия. 2003. №4. С. 60-63.
6. БеляевН. М., Рядно А. А. Методы теории теплопроводности. М.: Наука, 1982.