CC BY
73
УДК 539.3
Моделирование напряженно-деформированного состояния и потери устойчивости термобарьерного покрытия при тепловом ударе
П.А. Люкшин1, Б.А. Люкшин1,2, Н.Ю. Матолыгина1, С.В. Панин1
1 Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия 2 Томский университет систем управления и радиоэлектроники, Томск, 634050, Россия
В работе проведено компьютерное моделирование деформационного поведения термобарьернык покрытий. Показана возможность возникновения неустойчивостей, имеющих периодический характер. На примере термического нагружения медного образца с керамическим покрытием проведены исследования их характерного периода от свойств сопрягаемык материалов.
Ключевые слова: термобарьерное покрытие, компьютерное моделирование, потеря устойчивости, теплопроводность, напряженно-деформированное состояние
Simulation of stress-strain state and stability loss of a thermal coating under thermal shock
P.A. Lyukshin1, B.A. Lyukshin12, N.Yu. Matolygina1 and S.V. Panin1
1 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia
2 Tomsk University of Control Systems and Radioelectronics, Tomsk, 634050, Russia
The paper presents the results of computer simulation of the strain-induced behavior of thermal barrier coatings. The possibility of periodic instabilities to arise is demonstrated. The characteristics period of the instabilities in relation to the properties of conjugate materials is studied by the example of thermal loading of a Cu specimen with a ceramic coating.
Keywords: thermal barrier coating, computer simulation, stability loss, thermal conductivity, stress-strain state
1. Введение
Проблема разработки термобарьерных покрытий [1-5] на сегодняшний день не является полностью решенной несмотря на многолетнюю историю исследований в области газовой динамики, механики деформируемого твердого тела, механики разрушения, физического материаловедения и неразрушающего контроля (контроля целостности/состояния). Актуальность исследований в данном направлении связана не только с использованием термобарьерных покрытий в деталях оборудования ответственного назначения в энергетике и двигателестроении, но и с непрерывным поиском новых материалов, применение которых позволяет повышать рабочие температуры газовых потоков, что повышает коэффициент полезного действия энергетических установок, а также увеличивать ресурс работы. С точки зрения фундаментальных исследований термобарьерные покрытия представляют особый класс объек-
тов, которые не только должны обеспечивать защиту подложки от воздействия высокой температуры, обладать высокой адгезией, жаропрочностью и жаростойкостью, но и сохранять аэродинамические характеристики изделия в условиях многократного циклического приложения термомеханических нагрузок. Таким образом, исследования, направленные на повышение свойств термобарьерных покрытий, всегда будут актуальными, а компьютерное моделирование деформационного поведения термобарьерных покрытий позволяет глубже понять физику процессов разрушения покрытий и оптимизировать их проектирование.
Основная проблема функционирования термобарьерных покрытий связана с различием коэффициентов линейного температурного расширения хрупкого керамического покрытия, нанесенного на поверхность пластической подложки, и самой подложки. При температурном воздействии деформация подложки оказы-
© Люкшин П.А., Люкшин Б.А., Матолыгина Н.Ю., Панин С.В., 2011
вается существенно выше, что приводит к возникновению мощных концентраторов напряжений на границе раздела «покрытие-подложка», релаксация которых протекает путем формирования трещин. В более общей постановке указанная проблема может быть сведена к различию модулей упругости двух сопряженных сред, а расчет деформации может быть проведен для случая слоистого композита, что достаточно хорошо отработано в рамках аппарата теории упругости [6].
Расчет распределения напряжений на границе двух зерен в упруго нагруженном кристалле был проведен в [7], где было показано, что распределение имеет периодический осциллирующий характер с максимумом в середине двойного стыка и постепенным затуханием по направлению к краю зерна. Другим важным результатом, касающимся получения периодических решений при термическом нагружении тонкой пленки на массивной подложке, следует считать работу [8], в которой показано, что распределение напряжений в такой системе подчиняется синусоидальному распределению и зависит от толщины пленки.
Развитием этих результатов на случай двумерного распределения следует считать результаты работы [9], полученные при моделировании поведения тонких оболочек. Было показано, что вследствие потери устойчивости может быть получено регулярное двумерное распределение деформаций. Наглядное экспериментальное подтверждение «шахматного» распределения напряжений и деформации при термическом циклировании медного образца с наноструктурным теплозащитным покрытием получено в работе [10], что явилось еще одним подтверждением концепции о подобии процессов мас-сопереноса, развивающихся в объектах живой и неживой природы и способных иметь периодический («шахматный») характер распределения [11].
В настоящей работе была поставлена задача с использованием «классических» подходов механики деформируемого твердого тела на примере термического нагружения медного образца с керамическим покрытием показать возможность возникновения неустойчи-
востей, имеющих периодический характер, и провести исследования их характерного периода от свойств сопрягаемых материалов.
2. Описание параметров модели
Рассматривается пластинка AEDB с термопокрытием CDEF (рис. 1), которая подвергается интенсивному нагреву со стороны ED. Толщина пластинки значительно меньше двух других линейных размеров. За материал термопокрытия взят оксид алюминия А1203, материал подложки — медь. Такая постановка позволяет анализировать в плоском приближении характер распределения температуры, напряжений и далее возможность потери устойчивости термобарьерных покрытий.
Сетка конечных элементов приведена на рис. 2 и содержит 7200 элементов и 3 721 узел. При решении задачи о распределении температуры в расчетной области система линейных алгебраических уравнений, соответствующая этой сетке, содержит 3 721 уравнение. Нестационарная задача теплопроводности решалась по неявной схеме [12].
Для расчета были использованы следующие физикомеханические характеристики материалов неоднородной пластинки. Для оксида алюминия удельная теплоемкость С = 1106 Дж/(кг-°С), плотность р = 3970 кг/м3, коэффициенты теплопроводности в направлении осей х и у равны Кхх = Куу = 40 Вт/(м • °С). Для меди удельная теплоемкость С = 380 Дж/(кг-°С), плотность р = = 8 900 кг/м3, коэффициенты теплопроводности Кхх = = Куу = 385 Вт/(м-°С).
3. Решение задачи теплопроводности
Для определения поля температуры в двухслойной пластинке методом конечных элементов [13] решается нестационарная задача теплопроводности для области ABCDEF. На границе материалов с различными физикомеханическими характеристиками (линия FC) ставятся условия идеального теплового контакта: равенства тем-
Рис. 1. Поперечное сечение прямоугольной пластинки
Рис. 2. Сетка конечных элементов для решения задачи теплопроводности и определения напряженно-деформированного состояния
= 0.
(1)
пературы и тепловых потоков [14]. На кромках АЕ и BD задаются условия симметрии:
дТ
дХ „л
ил ae wx BD
Постановка таких условий (иногда их называют условиями теплоизоляции) означает, что вводимые границы искусственно сужают размеры области, подвергаемой анализу, а в реальной ситуации примыкающая область имеет точно такие же характеристики, как расчетная.
На кромках АВ, где поддерживается температура 0 °С, и ED, где температура равна 700 °С, задаются условия Дирихле:
Т\лв = ° TL = 70°. (2)
В момент времени t = 0 температура во всей области ABCDEF была равна нулю:
Т(х, у) = 0. (3)
При решении задачи методом конечных элементов уравнение теплопроводности с граничными условиями (1) и (2) сводятся к минимизации функционала [12]:
Х = J °.5
V
А, = ср,
к,.\ IX I +к„.
( дТ)2 дТ„
— +2К—Т dy dt
dV, (4)
Рис. 3. Поверхность, иллюстрирующая распределение температуры (а), изолинии температуры в двухслойной пластинке в момент времени t = 400 мкс (6)
где Кхх, Куу — коэффициенты теплопроводности; с — удельная теплоемкость; р — плотность; V — объем области.
Поле температуры и изолинии в неоднородной пластинке в момент t = 400 мкс приведены на рис. 3.
4. Анализ напряженно-деформированного состояния
Интенсивный нагрев неоднородной пластинки вызывает возникновение температурных напряжений. Для их вычисления решается задача напряженно-деформированного состояния. Предполагается, что на границе контакта слоев (линия FC) существует идеальная адгезия: равны перемещения и напряжения. Задача напряженно-деформированного состояния решается с использованием метода конечных элементов [12].
Потенциальная энергия упругого тела при воздействии на него поверхностных сил Рх, Ру, Р2 при наличии начальных деформаций е0 равна [12]
П = X | 1 {и}т [В(е) ]т р(е) ][В(е) ]{и}&У -
е У( е) 2
-X I {и}т[в(е)]т -
е У( е)
-ZJ {U}T[^(e)]T
es( e)
P(e)
X P(e)
Py
P(e)
dS,
где {Ц} — вектор перемещений; [ В(е)] — матрица, связывающая деформации и перемещения, содержащая производные от функций формы; [ D(e)] — матрица упругих характеристик; {е0} — вектор начальной деформации; [ N(е)] — матрица, содержащая функции формы; Рх, Ру, Рг — интенсивности поверхностной нагрузки вдоль осей х, у, г; S(е) — площади элементов, на которые действует поверхностная нагрузка; V — объем упругого тела; 51 — площадь упругого тела.
Разрешающие уравнения метода конечных элементов получаются из условия минимизации функционала полной потенциальной энергии механической системы П:
-^ = 0.
д{и }
В результате минимизации функционала П получаем систему линейных алгебраических уравнений:
[К]{Ц} = {^,
где [К] — глобальная матрица жесткости механической системы; ^} — вектор сил.
Глобальная матрица жесткости сетки конечных элементов равна сумме матриц жесткости отдельных элементов. Глобальный вектор сил^} равен сумме векто-
ров сил отдельных элементов {f(е)}:
{^} = £ {f(е)}.
е
Вектор сил одного конечного элемента при наличии поверхностных сил Рх, Ру, Р2 и начальных деформаций запишется следующим образом:
{/(е)} = I [В(е)]г[D(е)]{ео^У +
У (е)
+ | [N(e)]
5 (e)
P Є) p(e) py
p( e)
d5.
Вектор начальной деформации в случае плоского напряженного состояния равен
{e0} = aAT
где а — коэффициент температурного расширения; АТ— перепад температуры (нагрев) одного конечного элемента.
В случае плоского напряженного состояния интеграл, связанный с тепловым расширением (начальной деформацией), равен
I [ B(-'f [ С«](ЄоМГ = EiOr
v 2(1 - v)
bi
(5)
ьк ск
ь = Уj - Ук , С = хк - х] , где Е—модуль Юнга; V—коэффициент Пуассона; а( — толщина элемента. Коэффициенты bj, Сj, Ьк, ск в (5) получаются круговой подстановкой.
В случае плоского деформированного состояния вектор начальной деформации {е0} и матрица упругих коэффициентов [D(e)] изменяются [13]. Интеграл, связанный с тепловым расширением, в случае плоской деформации равен
V
С
I [ в(е) ]т [&е>]>,чцу=Е’аТ
У 2(1 2у)
ьк ск
Подобно проведенному выше расчету число конечных элементов в сетке составляло 7 200, число узлов было равно 3 721, число алгебраических уравнений, в
результате решения которых получается вектор перемещения, равно 7 442. В результате решения задачи были получены поля перемещений, деформаций, напряжений в поперечном сечении пластинки ABCDEF. При решении плоской задачи теории упругости в области ABCDEF задавались условия симметрии. На кромке АВ перемещение Vвдоль оси у равно нулю, т.е. У\АВ = 0, на линии 0102 (рис. 4, б) перемещение и вдоль оси х равно нулю, т.е. и | = 0.
Для расчета напряженно-деформированного состояния были использованы следующие физико-механические характеристики материалов: для оксида алюминия А1203 модуль упругости Е = 382 ГПа, коэффициент температурного расширения а = 7 • 10-6, коэффициент Пуассона V = 0.35; для меди модуль упругости Е = = 123 ГПа, коэффициент температурного расширения а = 16.5 • 10-6, коэффициент Пуассона V = 0.35.
На рис. 4, а приведена поверхность, иллюстрирующая распределение напряжений стп, вызванных неравномерным распределением температуры по толщине пластинки, на рис. 4, б — соответствующие изолинии напряжений стп, на рис. 4, в — кривая, характеризую-
0.0030
у, м 0.0010-
0.0005
Ш
^ -У А|2°3
CU
0.0015
0.0030 X, м
Рис. 4. Поверхность напряжений аи (а), изолинии температурных напряжений (6), распределение напряжений аи вдоль линии О1О2 (в)
щая распределение напряжений стп вдоль линии 0102 (рис. 4, 6). Из приведенных на рис. 4 результатов следует, что в слое FCDE (оксид алюминия) возникают сжимающие напряжения порядка 450...470 МПа, причем сжимающие напряжения действуют вдоль осей х и г. Таким образом, слой оксида алюминия, нанесенный на медную подложку и подверженный неравномерному нагреву, находится в условиях двустороннего сжатия.
Поскольку на границе раздела возникают значительные сжимающие напряжения, возникает вопрос, могут ли температурные напряжения быть причиной выпучивания тонкого слоя оксида алюминия на медной подложке. Для исследования этого вопроса были проведены расчеты на устойчивость термопокрытия.
5. Устойчивость термопокрытия
В данном разделе слой оксида алюминия на медной подложке моделируется как пластинка на упругом винк-леровском основании [15]. Принимается, что полученные выше распределения напряжений в тонком слое термобарьерных покрытий одинаковы для всех направлений в плоскости, параллельной поверхности термобарьерных покрытий.
Уравнение устойчивости пластинки на упругом основании имеет вид [15, 16]:
д w д w
^ 4 + ;л 2^ 2
дх дх ду
w
Л
ду4
+ ах
w
дх
+ 2т
где D = -
д 2-
д 4-
ду4
(6)
Eh
—-----цилиндрическая жесткость; h —
12(1-v2)
толщина пластинки; k — коэффициент постели упругого основания; v — коэффициент Пуассона; w — нормальный прогиб; ах, а и т — нормальные и касательное напряжения, действующие в срединной поверхности пластинки в момент потери устойчивости.
Однородное уравнение (6) решается совместно с однородными граничными условиями. На кромке x = = const возможны два варианта:
условия шарнирного опирания д2 w
w = О’ ЭхГ = ()
условия жесткого защемления дw
w = 0, — = 0. (8)
дх
Граничные условия на кромке y = const записываются следующим образом:
условия шарнирного опирания
w = 0’ дт = 0 ду2
(9)
условия жесткого защемления
w = 0, ^ = о.
’ ду
(10)
Уравнения устойчивости (6) совместно с граничными условиями (7)—( 10) представляют собой систему однородных дифференциальных уравнений в частных производных с однородными граничными условиями. Условие нетривиальности решения системы однородных уравнений определяет критическую нагрузку пластинки на упругом основании под действием сжимающих температурных напряжений а х = ау.
Таким образом, вопрос нахождения критической нагрузки на пластинку на упругом основании сводится к вопросу о нахождении собственных чисел дифференциального уравнения в частных производных. В работе эта проблема решалась с применением метода конечных разностей. Заменяя в уравнении (6) дифференциальные операторы их конечно-разностными аналогами во всех узлах конечно-разностной сетки, получим однородную систему линейных алгебраических уравнений, в которой число уравнений меньше числа неизвестных [15, 17]. Заменяя в граничных условиях (7)-( 10) дифференциальные операторы их конечно-разностными аналогами, получим систему линейных алгебраических уравнений, которая добавляется к уже существующей. В итоге после дискретизации уравнений устойчивости и граничных условий получается однородная система линейных алгебраических уравнений, в которой число уравнений равно числу неизвестных. Для существования нетривиального решения этой системы необходимо, чтобы ее определитель был равен нулю.
В настоящей работе для вычисления определителя системы линейных алгебраических уравнений использовался метод Гаусса [18]. Следует отметить, что для вычисления значения определителя в процедуре Гаусса достаточно выполнить прямой ход, обратный ход делать нет необходимости. Значение напряжения, при котором определитель равен нулю (точнее говоря, интервал напряжений, в котором определитель меняет знак), принималось критическим. На практике определялось не численное значение определителя, а его знак для различных значений напряжений ах, ау. Определив интервал напряжений, в котором определитель меняет свой знак, легко определить критическую нагрузку с любой разумной степенью точности. После того как определено собственное число (критическое напряжение), находится собственный вектор системы алгебраических уравнений. Для этого полагалось, что одна компонента (с номером п) собственного вектора равна константе, другие компоненты собственного вектора определялись в результате решения системы (п - 1) уравнений. Таким образом, находился собственный вектор системы линейных алгебраических уравнений (форма потери устойчивости) с точностью до постоянного множителя.
Для проведения расчета использовали конечно-разностную сетку, содержащую 3600 ячеек и 3721 узел. При этом напряжения ах и ау, найденные для узлов сетки конечных элементов (рис. 2) могут быть использованы в конечно-разностной сетке.
На рис. 5 приведены критические нагрузки и формы потери устойчивости пластинки (покрытия) при одновременном сжатии в двух направлениях х и у. Предполагалось, что в докритическом состоянии напряжения ах и ау равны друг другу. Алгоритм нахождения критических напряжений изменяется незначительно, если соотношение напряжений ах и ау будет произвольным.
При решении задачи устойчивости ставились два вида граничных условий: шарнирного опирания и жест-
кого защемления. Коэффициент постели упругого основания в случаях, рассмотренных на рис. 5, равен нулю. Аналитическое решение нахождения критических напряжений пластинки при шарнирном опирании и жестком защемлении по краям приведено в [17].
Критические нагрузки и формы потери устойчивости пластинки на упругом основании приведены на рис. 6. В качестве граничных условий использовано шарнирное опирание по всем кромкам. Отношение толщины пластинки к длине стороны вдоль оси х и у равно к!а = к/Ъ = У200. Модуль упругости Е = 382 ГПа, коэффициент Пуассона V = 0.35.
Из рис. 6 нетрудно видеть, что с увеличением жесткости упругого основания растет число волн на одной
ах-ау- 45.66 МПа (метод конечных разностей) ах = иу = 45.99 МПа (аналитическое решение)
0 00
0.00
Рис. 5. Критические напряжения и формы потери устойчивости квадратной пластинки под действием равномерных сжимающих напряжений ст^, ст^: шарнирное опирание (а), жесткое защемление по всем четырем кромкам (б). Кривые нормальных прогибов вдоль линии ОО2 (еУ-1 — шарнирное опирание, 2 — жесткое защемление по всем четырем кромкам
0.00 ' I '----------1--------1---------
0.00 0.03
0.00 0.03
Рис. 6. Критические напряжения пластинки на упругом основании при двустороннем сжатии. Отношение к/а = к/Ь = 1/200
Рис. 7. Критические напряжения пластинки на упругом основании при двустороннем сжатии. Отношение к а = к/Ь = 1/100
и той же площадке. Так, если коэффициент постели k = = 0.3 ГН/м3, то число полуволн при потере устойчивости равно 2 вдоль каждой из осей, если коэффициент постели к = 3 ГН/м3, то число полуволн равно 3 вдоль каждой из осей, если коэффициент постели к = 30 ГН/м3, то число полуволн возрастает до 5 вдоль каждой из осей, причем размеры площадки не изменяются.
Если отношение толщины пластинки к длинам сторон вдоль осей х и у увеличить до ка = к/Ъ = 1/100, то число полуволн при потере устойчивости уменьшается при той же жесткости основания и размерах пластинки. На рис. 7 приведены критические напряжения и формы
потери устойчивости пластинки с отношением толщины к длине и ширине пластинки к а = к/Ъ = 1/100.
6. Выводы
При равномерном нагревании двухслойной пластинки в термопокрытии действуют растягивающие напряжения, в подложке — сжимающие. Данные напряжения являются следствием различных коэффициентов температурного расширения термопокрытия и подложки и могут достигать значений, опасных с точки зрения нарушения прочности и/или устойчивости покры-
тия. Так, при тепловом ударе в термопокрытии возникают значительные сжимающие напряжения порядка 500 МПа.
При потере устойчивости термобарьерного покрытия оно ведет себя как сжатая пластинка на упругом основании. Число образующихся при потере устойчивости волн, располагающихся в шахматном порядке, растет с увеличением жесткости упругого основания и уменьшением отношения толщины пластинки к ее ширине и длине. В этом отношении термобарьерное покрытие ведет себя как сжатая тонкая упругая пластинка на упругом основании Винклера.
Литература
1. Evans A.G., Mumm D.R., Hutchinson J.W., Meier G.H., Pettit F.S. Mechanisms controlling the durability of thermal barrier coatings // Progr. Mater. Sci. - 2001. - V. 46. - P. 505-553.
2. Koolloos M.F.J. Behaviour of Low Porosity Microcracked Thermal Barrier Coatings under Thermal Loading. Ph. D. Thesis. - Eindhoven: Technische Universiteit Eindhoven, 2001. - 168 p.
3. Nissley D.M. Thermal barrier coating life modeling in aircraft gas turbine engines // J. Therm. Spray Tech. - 1997. - V. 6. - No. 1. - P.91-98.
4. Xie L., Ma X., Jordan E.H., Padture N.P., Xiao D.T, GellM. Deposition mechanisms of thermal barrier coatings in the solution precursor plasma spray process // Surf. Coat. Tech. - 2004. - V. 177-178. -P. 103-107.
5. Jadhav A.D., Padture N.P., Jordan E.H., Gell M., Miranzo P., Ful-lerE.R. Jr. Low-thermal-conductivity plasma-sprayed thermal barrier coatings with engineered microstructures // Acta Mater. - 2006. -V. 54. - P. 3343-3349.
6. Немировский Ю.В., Мищенко А.В. Влияние выбора материалов и структуры конструкции на пластическое деформирование и разрушение слоистык стержневых систем // Физ. мезомех. - 2004. -Т. 7. - Спец. вып. - Ч. 1. - С. 180-183.
7. Гриняев Ю.В., Панин В.Е. Расчет напряженного состояния в упру-гонагруженном поликристалле // Изв. вузов. Физика. - 1978. -№ 12. - С. 95-101.
8. Cherepanov G.P. On the theory of thermal stresses in a thin bonding layer // J. Appl. Phys. -1995. - V. 78. - P. 6826-6832.
9. Черепанов О.И. Численное решение некоторых квазистатичес-ких задач мезомеханики. - Новосибирск: Изд-во СО РАН, 2003. -180 с.
10. Панин В.Е., Сергеев В.П., Панин А.В., Почивалов Ю.И. Наноструктурирование поверхностных слоев и нанесение нанострук-турнык покрытий — эффективный способ упрочнения современные конструкционные и инструменталынык материалов // ФММ. -2007. - Т. 104. - № 6. - С. 650-660.
11. Панин Л.Е., Панин В.Е. Эффект «шахматной доски» и процессы массопереноса в интерфейсных средах живой и неживой природы // Физ. мезомех. - 2007. - Т. 10. - № 6. - С. 5-20.
12. Зенкевич О., Морган К. Конечные элементы и аппроксимация. -М.: Мир, 1986. - 318 с.
13. СегерлиндЛ. Применение метода конечнык элементов. - М.: Мир, 1979. - 392 с.
14. Плятт Ш.Н. Расчеты температурных полей бетонных сооружений. - М.: Энергия, 1974. - 407 с.
15. Тимошенко С.П., Войновский-Кригер С. Пластинки и оболочки. -М.: Наука, 1966. - 636 с.
16. Прочность, устойчивосты, колебания: Справочник в 3-х т. Т. 3. -М.: Машиностроение, 1968. - 567 с.
17. КармишинА.В., ЛясковецВ.А., МяченковВ.И., ФроловА.Н. Статика и динамика тонкостенных оболочечных конструкций. - М.: Машиностроение, 1975. - 376 с.
18. Зенкевич О. Метод конечнык элементов в технике. - М.: Мир, 1975. - 541 с.
Поступила в редакцию 17.09.2010 г.
Сведения об авторах
Люкшин Петр Александрович, к.ф.-м.н., снс ИФПМ СО РАН, petrljuk@ispms.tsc.ru Люкшин Борис Александрович, д.т.н., внс ИФПМ СО РАН, lba2008@yandex.ru Матолыгина Наталья Юрьевна, к.ф.-м.н., нс ИФПМ СО РАН, ksa@ispms.tsc.ru Панин Сергей Викторович, д.т.н., доц., зав. лаб. ИФПМ СО РАН, svp@ispms.tsc.ru
