Расчет динамики термоупругих напряжений в керамическом клапане методом конечных элементов Текст научной статьи по специальности «Физика»

УДК 519.6

РАСЧЕТ ДИНАМИКИ ТЕРМОУПРУГИХ НАПРЯЖЕНИЙ В КЕРАМИЧЕСКОМ КЛАПАНЕ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ

© В. И. Ткачев1*, В. В. Чудинов1, Н. Д. Морозкин2

1 Башкирскогий Государственный Университета, Бирский филиал Россия, Республика Башкортостан, 452453 г. Бирск, ул. Интернациональная, 10.

Тел./факс: +7 (34784) 4 04 09.

2Башкирский Государственный Университет Россия, Республика Башкортостан, 450076 г. Уфа, ул. Заки Валиди, 32.

Тел./факс: +7 (347) 229 63 70.

E-mail: tvivlad@gmail.com

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

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

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

Пусть Ое Я2 - занимаемая клапаном область с границей Ге дО. Распределение температуры области описывается уравнением теплопроводности

дТ д2Т д2ТЧ

^ эТ = *(д*2 +д7},

(1)

где Т(х,у,/), х,уеО, /е [0,Т],р - плотность

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

Начальные и граничные условия имеют вид

Т (X, y)|t=0=T1 m^yl а(т - Тж), дп

(2)

где а - коэффициент теплообмена, п - вектор внешней нормали к границе области, Тх - температура расплавленного металла.

Перемещения, возникающие в области О при нагреве, описываются уравнениями

Е д ,ди дуч

--(— + —) +

2(1 -у) дх дх ду

E

2(1 + Л)

E д

. аТЕ дТ

Дм + —--= 0

1 - у дх

(3)

,ди Эу. (— + —) + 2(1 -у) ду дх ду

Е . аТЕ дТ „

+-Ду + —--= 0

2(1 + у) 1 -у ду

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

ми напряжениями, возникающие в следствии изменения температуры, поэтому граничные условия для уравнения термоупругости (3) примут вид

\<ххПх + <хуПу = 0 \&xynx + <yyny = 0

где < - тензор напряжении.

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

Решение

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

Ш = с • вгпРгт ,

где вг - число Грасгофа, Рг - число Прандтля. Параметр т определяется выражением

т = 0.3 + (0.02/Рг1/3) .

Значения постоянных с и п зависят от величины числа вг. При Ог = 102 -109 (ламинарный режим) значение с = 0.52 и п = 0.25 если же вг > 109 с = 0.105 и п = 1/3 (турбулентный режим теплообмена) [7]. Само число Грасгофа определяется по

формуле Gr =

в' (Т -Т^) • g • l3

здесь в - коэффи-

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

+

2

V

* автор, ответственный за переписку

Используя значение числа Нуссельта, коэффициент теплообмена можно вычислить по формуле

к ■ Ш

а =

l

Задачу (1)-(2) будем решать методом конечных элементов. Для этого от задачи теплопроводности перейдем к вариационной задаче, воспользовавшись принципом стационарности дополнительной энергии [1]. Вариационная постановка задачи эквивалентная уравнению теплопроводности (1) записывается в виде

( дТ ^

dq

da) =

1 у ,r

п дТ

^ (| рсТ + | к ■ Т,

i=1 п п

п дТ

= |(а-(Т-Тх)—)€Г

i=l г

через запятую обозначены производные по переменным х и у, г = 1,2 . Интеграл в правой части запишем в виде суммы интегралов

n дТ

^Sqi (JpcTddrdü+ Jk - T,

1=1 a dqi a

f дТ ^

dq

da) =

1 / ,r

n дТ дТ

= -Z%i (JaT-dF- Ja- T -—dT)

i=1 г dqi г dqi

. (4)

Искомое температурное поле

Т = Т(д1,д2,...,дп,х,у,/) рассматривается как функция, зависящая от п обобщенных координат qi, которые являются неизвестными функциями времени. Систему (4) из п уравнений с п неизвестными функциями от обобщенных координат qi можно рассматривать как дополнительную форму уравнений Лагранжа, описывающих закон сохранения энергии.

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

— = Z qi (t) N (x, У) i=1

(5)

где qi (/)- неизвестные коэффициенты, Ni (х, у)-базисные функции. В теории МКЭ обобщенные координаты отождествляются с неизвестными весовыми коэффициентами на каждом узле сетки.

Подставляя решение (5) в уравнение (4) получим

( \ |рс ■ NN^€¡1 +

i=1 j=1

+ Jk - NhrN],rq1da + a

+ J(oNiNJqi )dr

. г

= ZJo-T^dr (6)

j=1r

Решение вариационной задачи (6) будем искать в линейном приближении. Для каждой вершины элемента строим базисные функции вида

N(х, у) = ах + Ьу + с . (7)

Для первого узла элемента коэффициенты получим из системы уравнений

а1 х1 + Ь1 у1 + с1 = 1

• а1 х2 + Ь1 у 2 + с1 = 0,

а1 х3 + Ь1 у3 + с1 = 0 где (х1, у1), (х2, у2), (х3, у3) координаты вершин элемента. Аналогичным образом находятся базисные функции для двух других вершин элемента.

Запишем уравнение (6) используя базисные функции вида (7) для одного элемента Пе с вершинами i ] к. Для удобства интегралы запишем в матричной форме, тогда

рс | N]Т [N]€Пе =

= Pc J

ae

= Pc J

ae

Ni Nj

Nk

N,2

NiNj

[Ni Nj Nk We =

NiNj NiNk

N 2

NjNk

NkNi NkNj N2

dae

матрица производных для одного треугольного элемента имеет вид

[B] =

следовательно

dNi dN; dNk

dx dx dx

dNi dNj dNk

dy dy dy

ai aj ak Ь bj bk

k J [ B]T [B]dae = kae

ae

ai aj ak bi bj bk

ai bi aj bj

ak bk

= kae

a2 + bj

aiaj + bibj aiak + bibk

ajai + bjbi aj + b j ajak + bjbk

akai + bkbi akaj + bkbj

aj + bl

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

aJ [ Ni ]T [ Nj ]dre =aJ

r r

Nj

NNj 0

n1nj N2

0

dre.

Таким образом, уравнение (6) для одного элемента запишем в виде

a

a

n n

рс | [N]т [N] ЙПе + к | [В]т [В]ЙПе +

пе пе

+ а$ [ Ni ]т [Ы} ]ЙГе =а\Тх йГе

птгл^ ^-Л^е

Г г

(8)

Просуммируем по всем элементам области уравнение (8)

т

Е (рс |[N]т [N]ЙПе • [д] + к |[В]т [В]ЙПе [д] +

е=1 пе пе

т

+ а\ [ N ]т [ N] ]ЙГе [д]) = Еа\ т„ ЙГе

Ге е=1 Г

Обозначим матрицы коэффициентов

т

С = Ерс |[N]т [N]йПе ,

е=1 пе

т

В = Е(к \[В]т [В]ЙПе +а\[N]т [N] ]ЙГе),

е=1 пе Г

т

Р = Еа\т^ЙГе,

е=1 Ге

тогда (8) запишется в виде системы С[д] + В[д] = Р

Получили систему из т линейных дифференциальных уравнений. Заменив производную по времени, конечной разностью получаем

с с

(С + В)д] = р + Т т

(9)

где т - шаг по времени.

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

Задачу термоупругости (3) также заменим эквивалентной вариационной задачей, основанной на принципе виртуальных (возможных) перемещений. Запишем вариационное уравнение для случая, когда на тело действуют поверхностные силы / и термические нагрузки, возникающие на неравномерном тепловом поле.

Вариационная задача равносильная уравнениям (3) имеет вид:

| (Леде + 20ец8ец )йП-п

- 12(Я + 0)ат(т-т0)деёП. =, (10) п

= \^,]3и1п]аг

Г

где е = £хх + £уу - деформация, Я = —уЕ _ , хх ^ (1 -у)2

О-

Е

- постоянные Ламе для плоского на-

2(1 + у)

пряженного состояния, щ = и, и2 = V, i, ] = 1,2. Уравнение (10) описывает изменение полной тер-

В =

Е

1 -У

£ = £

£

X су

£у ]т =

ди ^

— + —

ду дх

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

Запишем (10) в матричной форме. Постоянные Ламе Я и О запишутся в виде матрицы констант упругости

у 0 " у 1 0 0 0 (1 -у)/2

Соотношения Коши, выражающие деформации через перемещения примут вид

ди дv дх ду

Перемещения запишем в виде вектора {и} = [и Vт. Решение задачи будем искать на неструктурированной расчетной сетке с треугольными ячейками, использованной для решения задачи теплопроводности. Перемещения аппроксимируем кусочно-линейными функциями вида

п п

и = Е (х у)и,, v = Е (x, y)vi , i=\ ,=1

где щ и vi - неизвестные коэффициенты.

Соотношения Коши для одного элемента с узлами i ] к имеет форму

[ В] =

N 0 0 дМк 0

дх дх дх

0 д.N 0 0 дЩ

ду ду ду

дЫг дЫг дыз дЫ} дык дЫк

(11)

дх ду дх ду дх ду Таким образом, уравнение (11) для одного элемента можно записать

I [В]т [В][В]ЙПе =

пе

= | [ В]т [В]{£т }ЙПе + ¡{Р}(и }ЙГе

пе ге

где £т = ат Дт{1 1 0}т - вектор температурных деформаций. Суммируя (11) по всем элементам получим систему линейных алгебраических уравнений

К{и} = ^,

где {и} - вектор неизвестных коэффициентов для функции перемещений, К - матрица жесткости, Г -вектор узловых нагрузок. Элементами матрицы жесткости и вектора узловых нагрузок являются

ке = IБTВБJdпе, /е = IЫ'рЙГе. пе ге

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

Предел прочности на растяжение Максимальные растягивающие напряжения

-1-1-'-

0 5 t, с

---Предел прочности на сжатие

-Максимаьлные сжимающие напряжения

0

tc 10

6)

Рис. 1. Динамика максимальных растягивающих

и сжимающих термоупругих напряжений.

Описанный алгоритм реализован в среде про-

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

Результаты вычислительного эксперимента

Керамический клапан имеет следующие теп-лофизические свойства: ц = 0.25, E = 2.76 1011 Па, р = 3600 кг/м3, с = 920 Дж/кг-K, Ь = 25 Вт/м-K,

а = 8 10-6 1/K . Предел прочности на растяжение и сжатие равны соответственно ор = 2.21 108 Па, Oc = 2.484 109 Па. Расчеты проведены при температуре металла T^= 1400 °C , в течение времени T = 10 с.

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

T(х,y)|t=0 = 20 °C рис. 1а и T(x,y) |t=0 = 350 oC рис. 1 б.

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

Для проведения расчетов была построена сетка на 622 узлах, состоящая из 1098 элементов рис. 2. На рис. 3 показано распределение температуры в клапане через 10 сек. после погружения при

T(х,y)|t=0 = 350 °C . На рис. 4 и 5 приведено распределение термоупругих напряжений по главным осям тензора напряжений, которые вычисляются по формуле:

5 = (0 + 0 )/2 ±

0

-°yy )/2)2 +<

Рис. 2. Расчетная сетка.

«Оа2Я*5,18БНМ -535} 136,94543601 5J7i«MJ)?751M 1W217E91.209595 153S SOS23.34)673 21Э0М1 B5.J73752 252517407.GEJ533 301S5M«.737308 351343ЕЭ1.883387

^m 4ЯВ л'л. оегобб Рис. 4. Распределение максимальных главных напряжений (Па) при t = 10 с.

ноз 1295 11» 1Kb

9Н0 Е75 776 665 56Q

Рис. 3. Распределение температуры (°C) при t = 10 с.

23726323,628SB37 -23325624,5163669 гй3775йда Л 76 -1369295?! ,106363 ■1SU481163,351219 -244033*17,536169 -237565365,9411 Z -351137314,006071 -««89262.331021 -45^1210.575072

Рис. 5. Распределение минимальных главных напряжений (Па) при t = 10 с.

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

ЛИТЕРАТУРА

1. Био М. Вариационные принципы в теории теплообмена М.: Энергия, 1975.

2. Коваленко А. Д. Основы термоупругости. Киев: Наукова думка, 1975, 301 с.

3. Новацкий В. Динамические задачи термоупругости. М.: Мир, 1970.

4. Боли Б., Уэйнер Дж. Теория температурных напряжений. М.: Мир, 1964.

5. Сегерлинд Л. Применения метода конечных элементов. М.: Мир, 1979.

6. Кислицын А. А. Основы теплофизики. Тюмень, 2002, 152 с.

7. Исаченко В. П., Осипова В. А., Сукомел А. С. Теплопередача. М.: Энергия, 1975. 488 с.

Поступила в редакцию 21.02.2014 г.

ISSN 1998-4812

Вестннк EamKHpcKoro yHHBepcHTeTa. 2014. T. 19. №1

13

CALCULATION OF THERMAL STRESSES DYNAMICS IN THE CERAMIC VALVE WITH FINITE ELEMENT METHOD

© V. I. Tkachev1*, V. V. Chudinov1, N. D. Morozkin2

1Bashkir State University, Birsk Branch 10 Internatsionalnaya St., 452453 Birsk, Republic of Bashkortostan, Russia. 2Bashkir State University 32 Zaki Validi St., 450074 Ufa, Republic of Bashkortostan, Russia.

Phone: +7 (34784) 4 04 09. E-mail: tvivlad@gmail.com

The dynamics of thermal stresses in the ceramic valve during operation of melting unit is studied. Thermal stresses calculation is carried out in a two-dimensional approximation with the heat transfer at the boundary between the ceramic valve and the liquid metal. The necessity of such calculations is connected with the fact that the ceramic valve gets destroyed from the appearing thermal stresses while working. The thermal stresses are calculated with the finite element method, which allows to consider valve geometrical peculiarities more properly.

Keywords: thermal field, heat exchange, thermoelastic stress, finite element method, mathematical modeling.

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 teploobmen [Variational Principles in the Theory of Heat Transfer]. Moscow

2. Kovalenko A. D. Osnovy termouprugosti [Fundamentals of Thermoelasticity]. Kiev: Naukova dumka, 1975,

3. Novatskii V. Dinamicheskie zadachi termouprugosti [Dynamic Problems of Thermoelasticity]. Moscow: Mir, 1970.

4. Boli B., Ueiner Dzh. Teoriya temperaturnykh napryazhenii [Theory of Thermal Stresses]. Moscow: Mir, 1964.

5. Segerlind L. Primeneniya metoda konechnykh elementov [Application of the Finite Element Method]. Moscow: Mir

6. Kislitsyn A. A. Osnovy teplofiziki [Fundamentals of Thermal Physics]. Tyumen', 2002,

7. Isachenko V. P., Osipova V. A., Sukomel A. S. Teploperedacha [Heat Transfer]. Moscow: Energiya, 1975.

: Energiya, 1975.

, 1979.

Received 21.02.2014.