CC BY
29
УДК 536.24
УПРАВЛЕНИЕ ТЕПЛОВЫМ АГРЕГАТОМ ПРИ ОХЛАЖДЕНИИ КЕРАМИЧЕСКИХ ИЗДЕЛИЙ С ОГРАНИЧЕНИЕМ НА ТЕРМОУПРУГИЕ
НАПРЯЖЕНИЯ
© В. И. Ткачев
Башкирский государственный университет, Бирский филиал Россия, Республика Башкортостан, 452453 г. Бирск, ул. Интернациональная, 10.
Тел./факс: +7 (34784) 4 04 09.
Email: tvi-vlad@mail. com
Исследуется процесс управления охлаждением теплового агрегата в процессе обжига керамических изделий с учетом ограничений на термоупругие напряжения. Предложен подход, позволяющий вычислить температурный режим охлаждения изделия, не превышая при этом пределов прочности материала. Для вычисления температурных полей и термических напряжений в трехмерных моделях изделий использован метод конечных элементов. На примере держателя спирали проведен вычислительный эксперимент и определен температурный режим, при котором охлаждение агрегата осуществляется в течение заданного времени без превышения допустимых значений термоупругих напряжений.
Ключевые слова: Температурное поле, термические напряжения, метод конечных элементов, управление температурой.
Постановка задачи
Пусть ОеЯ3 - занимаемая изделием область с границей ГеШ рис. 1 Распределение температуры области описывается уравнением теплопроводности
дТ
РЛ — = 4 (V 2T),
dt
(1)
где Т(х,у^,0, х,у^еО, ^[0, Т ], ри - плотность материала изделия, Си- теплоемкость изделия, А« - коэффициент теплопроводности изделия.
Рис. 1. Керамический держатель спирали.
Начальные и граничные условия имеют вид
Т(х.ул01= о=То, дТ(X'У'2) 1г = а(Т - Т) , (2) дп
где а — коэффициент теплообмена, п - вектор внешней нормали к границе области, Тп — температура печи.
Поле напряжений, возникающее в области О при охлаждении, описывается уравнениями
XX + ХУ + XZ + X = = 0,
дх dY dz
дс„, do\„,
+ УУ + + Y - 0,
дх д dz
XX + yz + zz + 7 - 0,
дх dY dz
(3)
здесь ау — компоненты тензора напряжений. Границу поверхности держателя спирали считаем свободной от механических нагрузок, следовательно, граничные условия для уравнения термоупругости (3) примут вид
а XX «х + а ху п у + а Х2 п2 = 0 аух «х + ауу «у + аyz «2 = 0 ,
^ zx «x + ^ zy « y + ° zz «z
= 0
где ау- тензор напряжений, п, Пу, п - направляющие косинусы внешней нормали к граничной поверхности.
Необходимо определить режим изменения температуры печи, при котором охлаждаемое изделие достигнет конечной температуры за заданное время и не разрушится от возникающих термических напряжений.
Применение метода конечных элементов к вычислению термических напряжений в произвольной трехмерной области
Для решения поставленной задачи (1)-(2) применим метод конечных элементов в сочетании с методом Галеркина. Приближенное решение ищем в виде
_ п
Т = X Я, С)N (х, у, г) , (4)
1=1
где qi(t) - неизвестные коэффициенты, N (х,у^) -базисные функции.
Применяя метод Галеркина к уравнению (1) и подставляя в него приближенное решение (4) получим
„ „ / . л
II J[Njf
i=i j=in it
РиСи
dT ~dt"
•Ли (v 2T)
dQ. = 0, (5)
J
где [А£|Т - вектор базисных функций.
Наличие в (5) вторых производных исключает возможность применения линейных базисных
функций, поэтому представим второе слагаемое в подынтегральном выражении уравнения (5) в виде
| [ N ] ЛиУ2 ТОП = | у([ N ] ЛиУТ -а а
-1 V[N ] Л^ТОП
а
Преобразуем первый интеграл в правой части по теореме Остроградского-Гаусса получим
/ N]Т ЛЦЧТ )йП = | 4 -д—dГ , п г дп
тогда подынтегральное выражение в (5) запишется
I [ N ]т
Р и си
дг
дТ Яи (V 2 Т) | = | Ри Си +
+ I V[ N ] Яи VTdn -| Яи — ¿Г
п г дп
В результате, с учетом граничных условий (2) и решения (4) при постоянных ри, Си, ^и уравнение (5) примет вид
£ £ Ри си I]т Щ ]dП■
¿=1 )=1 V п
Яи I ) ]т У[ N V п
дд1 (()
31
(6)
+ а/[И) ]т [И1 ]dГ 9, (О - а/[И) ]т Гп<^ = 0.
г ) г )
Дискретизацию расчетной области О проводим с помощью треугольных элементов, для чего строим конечно-элементную модель из т непересекающихся элементов с п узлами. На полученной модели ищем приближенное решение задачи теплопроводности (6).
Так как, базисные функции определены на элементе, то вклад отдельного элемента с номерами узлов г, ], к, I в (6) можно записать следующим образом
Ре Се Я*Е ]Т [Ме М пе д
+ (Ль \[Бе ]т [ВЕ
пе
+ а ¡[МЕ ]т [МЕ цге )[де 1Г
= а |[МЕ]Т Т ,
(7)
где [Ы]=[Ы N N Ы],
[Ве ] =
Просуммируем (7) по всем элементам области и запишем в виде
С[д ] + В[д] = Р,
т
здесь С = X Р и Си / [Ые ]Т [Ые №е,
В = X (Ли I[Ве ]Т [Ве + а |[Ые ]т [Ыеуге),
е=1 гЕ
дМ дМк и
дх дх дх дх
дМ, дМу дМк т, . V
ду ду ду ду
дМ, дМ дМк дМ, ж
_ дг д дг дг _
р = ха |[Nе]тТп.
е=1 ге
Получили систему из п линейных дифференциальных уравнений. Заменив производную по времени, конечной разностью получаем
(- + В)д' = Р + - д'-1, (8)
т т
где т — шаг по времени.
Таким образом, исходная задача сводится к системе линейных алгебраических уравнений. Для линейного уравнения теплопроводности матрицы С и В всегда являются симметричными, с преобладающими элементами на главной диагонали, поэтому для решения этой системы уравнений используем метод итераций.
Задачу термоупругости (3) решим, пользуясь принципом виртуальных (возможных) перемещений [3-5]. Уравнение (3) для случая, когда на тело действуют поверхностные силы и термические нагрузки, возникающие на неравномерном тепловом поле в тензорной форме примет вид
| (Ле8е + 2^£у8£у )dО -п
- |(3Л + 2^)ат(Т - Т0)8edQ. = , (9) п
= 1 ^/А'«/dГ
г
где е=Ехх+Еуу+Е22 - объемная деформация,
Л =-^-, ц = —Е--постоянные Ламе,
(1 + V )(1 - 2у ) 2(1 + V)
5м; — виртуальные перемещения, % — компоненты тензора деформаций г,}=1,2,3.
Левая часть уравнения (9) описывает возможную работу внутренних термомеханических напряжений, а правая возможную работу граничных нагрузок.
Обозначим ы=Н1, У=У2, w=wз и запишем в виде вектора
(Ц)=[и V w]T Приближенные функции перемещений находим в виде
п
= Е ^ У' 2)мг'
I=1
п
* = Е (х' У'2 К-'
г =1 п
= Е (х' У'2 №'
г =1
где и, V, Wi, — неизвестные коэффициенты.
т
п
п
+
+
Г
Решение задачи находим на конечно-элементной модели, которая использовалась для решения задачи теплопроводности.
Запишем (9) в матричной форме для одного элемента. Постоянные Ламе запишутся в виде матрицы констант упругости (см. приложение). Соотношения Коши для одного элемента с узлами i} k I [Ве] имеет форму приведено в приложении.
Таким образом, уравнение (9) для одного элемента можно записать
|[Ве]т[Б][Ве] йО,е{ие} = пе
= I [Bе f [D]{s^ }dQe + |{N e}{P}dr'
, (10)
где {е^ }=атАТ{1 1 1 0 0 0}т - вектор температурных деформаций, {Ре} - вектор поверхностных нагрузок. Суммируя (10) по всем элементам получим систему линейных алгебраических уравнений К{Ц}=Р,
где {Ц}- вектор неизвестных коэффициентов для функции перемещений, К- матрица жесткости, F -вектор узловых нагрузок. Матрицей жесткости и вектором узловых нагрузок являются соответственно
к = £ |[Бе]т[Б][Бе ] ,
е=1пе
м
F = £ |[Бе ]т}<Юе + | [Ые}{Ре}<Ге.
Симметричность матрицы и преобладание главной диагонали позволяет применить для решения системы линейных алгебраических уравнений метод итераций.
Оценка коэффициента теплообмена
В высокотемпературных печах передача тепла происходит одновременно благодаря излучению и конвекции. При учете теплопередачи обоих видов для удобства вводится понятие о суммарном коэффициенте теплопередачи
а2=а„+ а*.
Здесь аи - коэффициент теплопередачи излучением, ак— коэффициент теплопередачи конвекцией.
При расчете печей с лучепрозрачной средой (например, электрических печей) обычно используется формула
1
1
1
F
•х,
T
4 Г T 44
и. i i ст 100 ) {100
= сп
T 100
4
T
ст
100
где Со - коэффициент излучения абсолютно черного тела, Еи, Ест, - интегральная степень черноты изделия и стенки печи соответственно, Ри
— площадь поверхности изделия, FCT — площадь внутренней поверхности печи, Ти — температура изделия, Тст— температура стенки печи [6].
Коэффициент теплопередачи излучением, можно вычислить по формуле
T Y (T Y
-*И I _ ст I
100 ) у 100)
Т - Т
и ст
В рассматриваемом случае конвективная составляющая теплообмена, определяется оценками для свободной конвекции в ограниченном пространстве [6, 7]
ак=(Е]ДсрУ1,
Ек- коэффициент конвекции, Аср — теплопроводность среды, I — расстояние между поверхностями. Коэффициент конвекции для случая От'Рт>1000 определяется выражением
ек=0.18(ОгРт)025 для случая СгРг<1000, Ек=1, т.е. в этом случае конвекция отсутствует.
Число Грасгофа, характеризующее интенсивность конвекции, вычисляется по формуле
Ог = (Ги - Тст)813
T -v2
ср г ср
где g - ускорение свободного падения, Тср- температура среды, Vср— кинематическая вязкость среды.
Алгоритм управления охлаждением печи
1. Задаем начальное распределение температуры изделия То, скорость нагрева Ун и охлаждения печи Уо, конечную температуру Тз;
2. Вычисляем распределение температуры в изделии для времени =+Лt;
3. Если минимальная температура изделия Тшв>Тз, то переходим к пункту 7;
4. Для полученного распределения температуры вычисляем максимальные термические напряжения аШах;
5. Если аШах< адоп, то с учетом величины У о снижаем температуру печи и переходим к пункту 2; иначе к пункту 6 (адоп - предельное допустимое напряжение);
6. Повышаем температуру печи с учетом Ун и переходим к пункту 2;
7. Вывод результатов.
Результаты вычислительного эксперимента
В соответствии с приведенным выше алгоритмом, проведен вычислительный эксперимент, в котором рассчитывается режим управления охлаждением керамического держателя с учетом ограничений на термические напряжения.
Теплофизические характеристики материала полагаем постоянными: ц=0.22, Е=1.46-10и Па, р=3800 кг/м3, с=1046 Дж/кгК, А=9.2Вт/мК, ат=1.6340-5 1/К. Предел прочности на растяжение и сжатие согласно работам [8, 9] равны ар=5.12-107 Па, ас=5108 Па соответственно.
Q
Г
+
V Е ст
и
X
4
Максимальная скорость охлаждения печи Vo=60 °C/min, максимальная скорость нагрева также равна VH=60 °C/min. Расчет охлаждения производится до максимальной температуры изделия Гз=50 °C.
Критические термонапряжения оцениваются, согласно первой теории прочности, применяемой для хрупких материалов. Опасным состоянием материала считается достижение одним из главных напряжений предела прочности [10].
Из приведенных рис. 2-3 видно, что керамический держатель спирали при заданной скорости охлаждения можно охладить за время t<70 мин., не превышая при этом предела прочности.
п—1—I—1—I—1—г
О 20 40 60 мин Рис. 2. Динамика максимальной температуры изделия (1), динамика минимальной температуры изделия (2), динамика температуры печи (3).
Рис. 3. Динамика максимальных растягивающих напряжений.
Полученный температурный режим позволяет избежать разрушения изделия вследствие, возникающих термических напряжений и сократить время охлаждения изделия в печи.
ЛИТЕРАТУРА
1. Био М. Вариационные принципы в теории теплообмена. М.: Энергия, 1975. 209 с.
2. Сегерлинд Л. Применения метода конечных элементов. М.: Мир, 1979. 393 с.
3. Коваленко А. Д. Основы термоупругости. Киев: Наукова думка, 1975. 301 с.
4. Новацкий В. Динамические задачи термоупругости. М.: Мир, 1970. 256 с.
5. Боли Б., Уэйнер Дж. Теория температурных напряжений. М.: Мир, 1964. 512 с.
6. Мастрюков Б. С. Теплотехнические расчеты промышленных печей М.: Металлургия, 1972. 368 с.
7. Кислицын А. А. Основы теплофизики. Тюмень: изд-во ТюмГУ, 2002. 152 с.
8. Гавриш Д. И. Огнеупорное производство. Справочник. М.: Металлургия, 1965. Т. I. 573 с.
9. Стрелов К. К., Мамыкин П. С. Технология огнеупоров. М.: Металлургия, 1978. 370 с.
10. Саргсян А. Е. Сопротивление материалов, теории упругости и пластичности. М.: Высшая школа, 2000. 286 с.
Общий вид матрицы констант упругости
Приложение.
D =
E(1 -v ) (1 + v )(1 - 2v )
1 V /(1 -v ) V /(1 -v ) 0 0 0
V /(1 -v ) 1 V /(1 -v ) 0 0 0
V /(1 -v ) V /(1 -v ) 1 0 0 0
1 - 2v 2(1 -v )
1 - 2v 2(1 -v )
1 - 2v 2(1 -v )
0
0
0
0
0
0
0
0
0
Соотношения Коши для одного элемента с узлами i j k l.
[Be ] =
0 0 0 0 dNk 0 0 ^NL 0 0
дх дх дх дх
0 Щ 0 0 дNj 0 0 дNk 0 0 щ 0
ду IT ду ду
0 0 Щ 0 0 дNj 0 0 дЩ 0 0 Щ
ду ду ду
Щ д^ 0 дNj дNj 0 дNk д^ 0 щ дЩ 0
ду дх дх ду дх ду дх
дМ, 0 дN UNj 0 дNj дNk 0 д^ щ 0 Щ
dz дх дх дz дх дz дх
0 щ Щ 0 дNj öNj 0 дNk дNk 0 дЫ1 дN
дz ду дz ду дz ду
Поступила в редакцию 03.09.2014 г.
THE THERMAL UNIT CONTROL IN COOLING PROCESS OF CERAMIC PRODUCTS WITH RESTRICTIONS ON THERMAL STRESSES
© V. I. Tkachev
Birsk branch of Bashkir State University 10 Internatsionalnaya St., 452453 Birsk, Republic of Bashkortostan, Russia.
Phone: +7 (34784) 4 04 09. Email: tvi-vlad@mail. com
The control process of thermal unit in the process of firing ceramic products taking into account the constrains on thermal stresses is studied. The approach allowing calculating the temperature regime of cooling products, without exceeding the material ultimate strength is proposed. The temperature fields and thermal stresses in three-dimensional models are calculated using the finite element method. On the example of a spiral holder, a computational experiment is conducted. The temperature regime, wherein the cooling unit is carried out within a specified time without exceeding the allowable values of thermo elastic stresses is defined.
Keywords: temperature field, thermal stresses, finite element method, temperature control. Published in Russian. Do not hesitate to contact us at bulletin_bsu@mail.ru if you need translation of the article.
REFERENCES
1. Bio M. Variatsionnye printsipy v teorii teploobmena [Variational Principles in the Theory of Heat Transfer]. Moscow: Energiya, 1975.
2. Segerlind L. Primeneniya metoda konechnykh elementov [Application of Finite Element Method]. Moscow: Mir, 1979.
3. Kovalenko A. D. Osnovy termouprugosti [Basics of Thermoelasticity]. Kiev: Naukova dumka, 1975.
4. Novatskii V Dinamicheskie zadachi termouprugosti [Dynamic Problems of Thermoelasticity]. Moscow: Mir, 1970.
5. Boli B., Ueiner Dzh. Teoriya temperaturnykh napryazhenii [Theory of Thermal Stresses]. Moscow: Mir, 1964.
6. Mastryukov B. S. Teplotekhnicheskie raschety promyshlennykh pechei [Heat Engineering of Industrial Furnaces]. Moscow: Metal-lurgiya, 1972.
7. Kislitsyn A. A. Osnovy teplofiziki [Basics of Thermal Physics]. Tyumen': izd-vo TyumGU, 2002.
8. Gavrish D. I. Ogneupornoe proizvodstvo. Spravochnik [Refractory Production. Handbook]. Moscow: Metallurgiya, 1965. T. I.
9. Strelov K. K., Mamykin P. S. Tekhnologiya ogneuporov [Technology of Refractories]. Moscow: Metallurgiya, 1978.
10. Sargsyan A. E. Soprotivlenie materialov, teorii uprugosti i plastichnosti [Strength of Materials, Theory of Elasticity and Plasticity]. Moscow: Vysshaya shkola, 2000.
Received 03.09.2014.
