CC BY
58
УДК 539.32
Влияние анизотропии на напряженно-деформированное состояние и потерю устойчивости керамического защитного покрытия
при тепловом ударе
П.А. Люкшин1, Б.А. Люкшин1, 2, Н.Ю. Матолыгина1, С.В. Панин1, 3
1 Институт физики прочности и материаловедения СО РАН, Томск, 634055, Россия 2 Томский университет систем управления и радиоэлектроники, Томск, 634050, Россия 3 Национальный исследовательский Томский политехнический университет, Томск, 634050, Россия
В работе с использованием компьютерного моделирования исследовано деформационное поведение термобарьерных покрытий. Обсуждается механизм возникновения неустойчивостей в таких покрытиях, основанный на их представлении в виде пластины, находящейся на упругом основании. Потеря устойчивости проявляется в виде образования двоякопериодической системы зон экструзии и интрузии, что качественно согласуется с известными результатами экспериментальных исследований. На примере моделирования термического нагружения медного образца с защитным керамическим покрытием исследованы характерные особенности возникающей картины потери устойчивости и ее зависимости от свойств сопрягаемых материалов. Получена оценка влияния анизотропии деформационно-прочностных свойств материала покрытия на характер возникающей неустойчивости.
Ключевые слова: термобарьерные покрытия, потеря устойчивости, тепловой удар, анизотропия
Anisotropy effect on stress-strain state and loss of stability in protective ceramic coating under thermal shock loading
P.A. Lyukshin1, B.A. Lyukshin1 2, N.Yu. Matolygina1, and S.V. Panin1, 3
1 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634055, Russia 2 Tomsk State University of Control Systems and Radioelectronics, Tomsk, 634050, Russia 3 National Research Tomsk Polytechnic University, Tomsk, 634050, Russia
The behavior of thermal barrier coatings was investigated by computer modeling. The mechanism of instability generation in such coatings modeled as a plate lying on an elastic foundation was studied. The loss of stability is manifested as the formation of a doubly periodic system of extrusion and intrusion zones, which qualitatively agrees with known experimental results. The thermal loading of a copper specimen with a protective ceramic coating was modeled to investigate the features of the observed pattern of stability loss and its dependence on the properties of the conjugate materials. The effect of the anisotropy of deformation and strength properties of the coating material on the observed instability pattern was estimated.
Keywords: thermal barrier coatings, stability loss, thermal shock, anisotropy
1. Введение
В классическом понимании термобарьерные покрытия предназначены для тепловой защиты деталей, ра-
ботающих при высоких температурах (до 1300 °С). Их толщина составляет от 100 до 500 мкм; в них применяется внутреннее охлаждение, в результате чего на поверхности металлической подложки (как правило, суперсплав на никелевой основе) температура снижается до 100-300 °С. Помимо подобных экстремальных условий внешнего воздействия, к которым также сле-
дует отнести термоциклирование, циклическое механическое нагружение, время их безотказной работы в промышленных турбинах должно составлять не менее 30000 ч [1, 2].
Несмотря на то что термобарьерные покрытия являются предметом глубокого исследования, включая аспекты механики, физики, материаловедения, химии и пр., на сегодняшний день проблема разрушения окончательно не решена. Это связано с тем, что в процессе эксплуатации в термобарьерных покрытиях развивают-
© Люкшин П.А., Люкшин Б.А., Матолыгина Н.Ю., Панин C.B., 2015
ся сложные многоуровневые процессы, вызванные различием коэффициентов термического расширения слоев; высокотемпературным окислением металлов, входящих в состав композиции; изменениями состава, структуры, морфологии и других структурных составляющих термобарьерных покрытий [1].
В научно-технической литературе, посвященной изучению поведения термобарьерных покрытий, отмечается, что такое большое количество процессов, какие свойственны для термобарьерных покрытий, сложно обнаружить в других технических системах. В результате для правильного понимания, описания и проектирования термобарьерных покрытий следует рассматривать не менее дюжины различных процессов, среди которых выделяют: диффузию, фазовые превращения, окисление, спекание, пластическую деформацию и ползучесть, тепловое расширение, излучение, разрушение, усталость ит.д. [1]. Безусловно, при разработке моделей поведения композиций с термобарьерными покрытиями важно учитывать все эти процессы, однако большинство исследований ограничивается лишь некоторыми из них.
Подход, развиваемый авторами настоящей работы, не стал исключением из этого правила. В работе развивается модель деформационного поведения теплозащитного покрытия на основе аппарата механики деформируемого твердого тела, с анализом возможных вариантов распределения напряжений и деформаций, возникающих вследствие приложения термической нагрузки. Подход основан на представлении термобарьерного покрытия пластиной, размещенной на упругом основании, деформация которой может развиваться по механизму потери устойчивости (подобной потере устойчивости стержня при сжатии [2]). Подобные подходы описаны в литературе, например [3, 4]. В процессе отслоения покрытие может принимать форму гофра, а период формирующейся структуры является функцией величины адгезии покрытия и основы [5]. Количественные оценки величин сдвиговых напряжений, возникающих на границе «пленка - подложка» при растяжении и приводящих к отслаиванию, получены в работе [6].
С позиции развития модельных представлений о развитии деформационных процессов в термобарьерных покрытиях одной из наиболее исчерпывающих работ следует считать обзорную статью [7], в которой даны исчерпывающая классификация и описание механизмов разрушения термобарьерных покрытий, включая оценку напряжений и деформаций.
Двумя наиболее распространенными методами формирования термобарьерных покрытий являются воздушно-плазменное напыление, когда покрытие имеет характерную чешуйчатую структуру, и осаждение из паровой фазы при плавлении мишени электронным
лучом (ЕВ PVD), когда покрытие имеет характерное столбчатое строение. Для каждого из этих методов характерны свои закономерности разрушения покрытий. Для метода воздушно-плазменного напыления, когда нанесению покрытия предшествует пескоструйная подготовка поверхности, искусственные шероховатости являются местом возникновения концентраторов напряжений, вызывающих локальное гофрирование покрытия, что завершается ростом трещин в этой области
[8]. В модельных условиях статического нагружения (растяжение) поведение покрытий, нанесенных методом воздушно-плазменного напыления, рассмотрено в
[9]. Более подробное и систематическое изложение закономерностей разрушения покрытий, полученных воздушно-плазменным напылением, приведено в работе
[10]. С другой стороны, для покрытий, с формированных методом осаждения из паровой фазы при плавлении мишени электронным лучом, характерно отслаивание значительных по размеру фрагментов как целого. Результаты подобных исследований изложены в [11], а также в работе [12], где стойкость термобарьерных покрытий пытались повысить за счет формирования подслоя, содержащего лантан.
Одним из ключевых аспектов, определяющих разрушение термобарьерных покрытий независимо от метода их формирования, является постепенное формирование термически выращенного оксидного слоя, появление которого обусловлено проникновением кислорода с поверхности покрытия к границе между керамическим слоем на основе диоксида циркония, стабилизированного иттрием, и адгезионным подслоем, например из оксида алюминия. Численная модель разрушения термобарьерных покрытий при возникновении термически выращенного оксидного слоя и последующей осцилляции напряжений на таком интерфейсе предложена в работе [13]. В статье [14] предложен механизм формирования термически выращенного оксидного слоя в покрытии, сформированном методом воздушно-плазменного напыления.
Поскольку характер растрескивания термобарьерных покрытий связан с распределениями напряжений и деформаций, анализ картины разрушения позволяет понять развитие деформационных процессов и обсуждать возможные способы снижения риска формирования трещин. В работах [15, 16] показано, что характерное расстояние между трещинами (пространственный период растрескивания) имеет определенное значение. При этом указанный период растрескивания является функцией ряда параметров, включая сдвиговую прочность на границе раздела «покрытие - основа» [17, 18]. Поскольку термобарьерные покрытия являются многослойными структурами, для адекватного описания процессов их разрушения используются модели, учитывающие несоответствие упругой и пластической
деформации в отдельных слоях [19]. Для получения количественной оценки напряженно-деформированного состояния в термобарьерных покрытиях применяют классические подходы механики разрушения [20]. В работе [21] было проведено исследование периодического распределения адгезионных трещин на распределение напряжений на границе раздела и характер последующего разрушения. Развитием этих исследований следует считать работу [22], где характер разрушения термобарьерных покрытий анализировали исходя из периодического растрескивания поверхности покрытия.
После того как в термобарьерном покрытии возникает одна трещина или система трещин, их наличие в значительной степени влияет на поведение всей системы. Для получения количественной оценки такого влияния рассчитывают плотность поверхностных трещин [23] либо влияние уже существующей трещины и ее морфологии на характер ее распространения и разрушения покрытия в целом [24, 25]. Поскольку характер распределения напряжений и деформаций вблизи мест возникновения трещин либо внутренних границ раздела (включая границы сплэтов в покрытиях, полученных воздушно-плазменным напылением, а также термически выращенный оксидный слой в них [26]), необходимо использовать модифицированные способы задания расчетных сеток, включая способы ее измельчения вблизи мест концентрации напряжений [27].
Еще одним ключевым фактором, определяющим характер распределения концентраторов напряжений и последующего разрушения термобарьерных покрытий, является профиль границы раздела, что несколько более актуально для метода воздушно-плазменного напыления формирования покрытий. При этом в моделях предлагается учитывать толщины термически выращенного оксидного слоя, а также морфологию интерфейса [28]. Развитием таких исследований следует считать работу, в которой для анализа процесса растрескивания термобарьерных покрытий в определенных пределах меняли шероховатость на интерфейсе [29]. В некотором предельном рассмотрении шероховатость на границе раздела «покрытие - основа» может быть представлена как набор вершин и впадин. В работе [30] анализировали распределение напряжений и деформаций в областях таких вершин и впадин и последующего роста трещин от этих структурных неоднородностей. Также было проведено исследование влияния распределения напряжений на структурных неоднородностях на границе раздела на зарождение и распространение трещин, когда при повышении температуры система оказывается под действием растягивающих напряжений [31].
В завершение краткого обзора отметим, что для случая нетермобарьерных покрытий периодические решения в характере распределения напряжений и деформаций были получены в работах [32, 33]. В работе
[34] ««шахматное» распределение характера растрескивания и отслоения фрагментов наноструктурного теплозащитного покрытия было выявлено при термическом циклировании образца с подложкой из меди.
В настоящей работе с использованием классических подходов механики деформируемого твердого тела на примере термического нагружения медного образца с керамическим покрытием показана возможность возникновения неустойчивостей, имеющих периодический характер, и проведено их исследование в зависимости от свойств сопрягаемых материалов, в том числе от анизотропии деформационно-прочностных свойств покрытия, а также от толщины покрытия.
Представлены следующие этапы анализа: решение задачи о распределении температур в покрытии и подложке в переходном режиме под действием теплового удара; оценка параметров напряженно-деформированного состояния покрытия под действием температурных нагрузок; моделирование процесса потери устойчивости покрытия с примерами расчетов при различных геометрических и жесткостных параметрах покрытия и подложки; оценка напряжений в покрытии после потери устойчивости.
2. Описание параметров модели
В плоской постановке анализируется расчетная область, включающая в себя подложку АВСГ и термопокрытие CDEF (рис. 1), которая подвергается интенсивному нагреву (тепловому удару) со стороны поверхности ED. Принимается, что размеры АЕ = BD значительно меньше двух других линейных размеров, что и позволяет проводить анализ в плоской постановке. При проведении расчетов для материала термопокрытия использованы характеристики оксида алюминия А1203, для материала подложки — меди.
Для оксида алюминия удельная теплоемкость С = = 1106 Дж/(кг-°С), плотность р = 3970 кг/м3, коэффициенты теплопроводности в направлении осей X и У
Рис. 1. Конфигурация расчетной области
= 0.
(1)
равны Кхх = Куу = 40 Вт/(м-°С). Для меди коэффициенты теплопроводности Кхх = Куу = 385 Вт/(м-°С), удельная теплоемкость С = 380 Дж/(кг-°С), плотность р = 8900 кг/м3.
3. Решение задачи теплопроводности
Для определения поля температуры методом конечных элементов [35-37] решается нестационарная задача теплопроводности для области АБВЕ. На границе материалов с различными физико-механическими характеристиками (линия РС) задаются условия идеального теплового контакта: равенства температуры и тепловых потоков [35]. На кромках АЕ и ББ задаются условия симметрии:
дТ_
дх АЕ дх ВБ
Постановка таких условий отражает тот факт, что вводимые границы АЕ и ББ искусственно ограничивают размеры области, подвергаемой анализу, когда в реальной ситуации примыкающая область имеет точно такие же характеристики, как расчетная.
На кромках АВ и ЕБ задаются условия Дирихле:
т\ав = т° аео = т* (2)
Принимается, что в момент времени t = 0 температура во всей области АБВЕ постоянна и равна начальному значению:
Т (х, у) = Т° (3)
При решении задачи методом конечных элементов уравнение теплопроводности с граничными условиями (1) и (2) сводятся к минимизации функционала [35]:
=^т
дх
Х = 1 °-5
V
К
дТ_ дх
+ К
2
уу
ду
+ 2Хд-Гт
Эt
dV +
+ | qTdS + \ | (Т - Тсш )2d
(4)
где — площадь поверхности, где задан поток тепла; S2 — площадь поверхности, где происходит конвективный обмен тепла; Кхх, К — коэффициенты теплопроводности; Т — температура; 1 =с р, с — удельная теплоемкость, р — плотность; V — объем области; q — поток тепла; h — коэффициент теплообмена; Теш — температура окружающей среды.
Сетка конечных элементов содержит 7200 элементов и 3721 узел. При решении задачи о распределении температуры в расчетной области система линейных алгебраических уравнений, соответствующая этой сетке, содержит 3721 уравнение. Нестационарная задача теплопроводности решалась по неявной схеме [35]. Поле температуры и изолинии в расчетной области в момент t = 0.05 с приведены на рис. 2. Толщина покрытия составляла 0.5 мм, подложки — 11 мм.
4. Расчет температурных деформаций и напряжений в анизотропной среде
Для расчета напряженно-деформированного состояния анизотропной пластинки используется метод конечных элементов. Для этого записывается потенциальная энергия упругого тела при воздействии на него поверхностных сил Рх, Ру, р при наличии начальных деформаций 8° в виде функционала энергии П, из условия минимума которого ЭП
д{и }
где {и} — вектор перемещений. Получаем систему линейных алгебраических уравнений [ К ]{и} =
где [К] — глобальная матрица жесткости механической системы; {Р} — вектор нагрузок.
а
т, ос( 600-
0.03
К м
0.01
0.00
о2
-50
-50—
01
0.03
Рис. 2. Поверхность, иллюстрирующая распределение температуры (а), изолинии температуры (б), изменение температуры по центру пластинки (линия 0201) (в) в момент времени t = 0.05 с, координата у отсчитывается от точки 02к точке 01
Глобальная матрица жесткости сетки конечных элементов равна сумме матриц жесткости отдельных элементов:
[ к ] -£ [ k(е)],
е=1
матрица жесткости отдельного элемента записывается в виде
[к(е)] = | [В(е)]т[^е)]•[В(е)] ¿V.
v (е)
Глобальный вектор сил равен сумме векторов сил отдельных элементов {/(е)}:
{р}-£ {/(е)},
е=1
вектор сил одного конечного элемента при наличии поверхностных сил Рх, Ру, Р2 и начальных деформаций запишется следующим образом:
г Р (е) ^
{f(e)} — J [B(e)]T[D(e)]• {80}dV + J [Me)]T
y(e ) S(e )
P (e) y
p ( e)
d S.
Вектор начальной (температурной) деформации в случае плоского напряженного состояния в локальной системе координат Х1, У1, совпадающей с осями орто-тропии, равен:
{80 } =
À1 AT À 2 AT 0
где , X2 — коэффициенты температурного расширения вдоль осей ортотропии 01Х1, 0171; АТ — изменение температуры (нагрев, охлаждение) одного конечного элемента (рис. 3).
Вектор температурной деформации {е0} в глобальной системе координат ОХУ связан с вектором деформации {в0} в локальной системе координат 01Х171 соотношениями [38]:
^х 0
y0 Y ху 0
cos2 а
sin2 а
2sin а cos а
8Xl0
л
sin2 а
- sin а cos а sin а cos а -2 sin а cos а cos2 а - sin2 а
cos2 а
у 10 Yxl у 10
(5)
Равенство (5) в сокращенной форме можно записать в виде
Рис. 3. Взаимная ориентация глобальной ХОУ и локальной Х1ОУ1 систем координат
{80} = [T*Г{80},
где
[T ] =
cos а sin2 а
sin2 а
cos2 а
sin а cos а - sin а cos а
-2 sin а cos а 2 sin а cos а cos2 а - sin2 а
Матрицу [Т ] назовем матрицей трансформации для деформаций [38-40].
Закон Гука для анизотропного материала запишется в виде
{а} = [ D ]({е} - {ео}).
Матрица упругих характеристик [D] в случае, когда локальные оси (ортотропии) ОХ171 совпадают с глобальными осями ОХУ, равна
" dli 4 0 *
[ D] — [ D1] — d12 d 22 0 . (6)
0 0 d\3
Здесь
dl =-
d — "
-, di2 — ■
E1V 2
-,^22- 9 ^2 - , - G,
1 -У1У2 1 -У1У2 1 -У1У2
где Е1, Е2 — модули упругости вдоль осей Х15 71; У1, V2 — коэффициенты Пуассона; G — модуль сдвига.
Матрица упругих характеристик [D] в случае, когда локальные оси (ортотропии) ОХ171 не совпадают с глобальными осями ОХУ, равна [40]:
[D] - [Т] • [• [Т*]-1,
где [ D1 ] — матрица упругих характеристик в локальных осях ОХ171; [Г] — матрица трансформации для напряжений [40]:
[T ] =
cos а sin2 а
sin2 а
cos а
2 sin а cos а -2 sin а cos а
- sin а cos а sin а cos а cos2 а - sin2 а
Следует отметить, что [Т*]-1 - [Т]т, т.е. обратная матрица трансформации деформаций равна транспонированной матрице трансформации для напряжений [40].
С учетом вышеприведенного замечания матрица упругих характеристик может быть представлена в следующем виде:
[D] - [Т] • [д] • [Т]т
Следовательно, матрица упругих характеристик анизотропного тела будет полностью заполненной и запишется в виде
dn di2
*13
[D] =
В матрице [D] ее элементы равны:
di2 d 22 d 23
d13 d23 d33
" k 0 ci
0 ci bi
1 bj 0 cj
2A 0 cj bj
bk 0 ck
0 ck bk
- У( к ), ci = x( k )
d11 = d111 cos4 a + d22 sin4 a + + 2(d/2 + 2d33 )sin2 a cos2 a,
: d111 sin4 a + d22 cos4 a +
2_____2 _
d12 '
+ d\2 (cos4 a + sin4 a) - 4d]3 sin2 a cos2 a, d13 = [d22 sin2 a - d111 cos2 a +
+ (d\2 + 2d33 )(cos2 a - sin2 a)] sin a cos a, d23 = [d22 cos2 a - d111 sin2 a -- (d/2 + 2d33 )(cos2 a - sin2 a)] sin a cos a,
(7)
d33 = (d111 - 2d/2 + d22) sin2 a cos2 a +
+ d]3 (cos2 a - sin2 a)2. В случае плоского напряженного состояния анизотропной среды интеграл, связанный с тепловым расширением (начальной деформацией), входящий в вектор сил, равен [36]:
I [ B( e)]T[ D e )]{80}dr.
V
В формуле (8)
(8)
[ B
.(ebT _
2 A
'k
0
bj
0
bk 0
0
ci
0
cj
0
ci
b
bj ck
(9)
где Ь - у(/) - у(k), с - х(k) - х(Коэффициенты Ь^, с^, Ьк, ск в формуле (9) получаются круговой подстановкой.
На рис. 4, а, б приведены поверхность напряжений а11 и изолинии напряжений, соответствующие температуре, приведенной на рис. 2. На рис. 4, в приведена кривая изменения напряжения ап вдоль прямой О О л.
Граничные условия при решении задачи о напряженно-деформированном состоянии задавались следующим образом. Перемещение вдоль оси у на линии АВ равно 0, перемещение вдоль оси х в точке О1 равно 0:
1AB
= 0, и | a = 0.
ап, МПа
0.000
0.005 0.010 Y м
Рис. 4. Поверхность, иллюстрирующая распределение напряжений а11 в неравномерно нагретом теле (а), изолинии напряжений а11 (б), изменение напряжений а11 по центру пластинки (линия О2О1) (в) в момент времени t = 0.05 с, координата у отсчитывается от точки О2 к точке О1
5. Устойчивость термопокрытия
Слой оксида алюминия на медной подложке моделируется как пластинка на упругом винклеровском основании [39]. Уравнение устойчивости анизотропной пластинки на упругом основании имеет вид [38, 39, 41]:
Э4 г д4 г д4 г
А1 тг + 4 Аб т-Т + 2(Оп + 2 Б66) т-2Г"2 +
ох ох оу ох оу
+ 4 Л
Э 4
Э 4-
г ^ д г 7 д г + Ао —т- + ^ к—- +
Э2-
2б ^ 3 Л22 ^
дхду ду дх
д2-
э2-
+ 2тхук-+ аук—- - кг — 0.
ху у 2
ЭхЭу у Эу2
Здесь
к/ 2
а 1 к
а22к
г 2 т а11 к А11 = ] а11= "Т^' Л22 =-рТ -к 2 12 12
Л = а12^ П = ^3 Л12 = , „ ' Лбб
12
а1Лк3
Л - а1б
Л1б
А =
12
а2б к3
12 ' Л26 12
(10)
h — толщина пластинки; коэффициенты аи, а22,..., абб входят в физические соотношения для анизотропной пластинки, связывающие напряжения и деформации следующим образом [41]:
Рис. 5. Оси ортротропии материалаХ1, Y1 повернуты относительно глобальных осей X, Y на угол а
а х = а118 х + а128 у + а1б У ху' а у = а 218 х + а228у + а 2б У ху'
(11)
т ху = а618 х + а62 8у +
ху
Для ортотропного материала коэффициенты а11, а22,абб в случае, когда оси ортотропии материала Х1, У1 совпадают с глобальными осями X, Y, равны
а11 = а11 = л ' а22 = а22 л 1 -V1V2 1
= 0 = Е2 V1 = 0 = г а12 = а12 = л ' абб = абб = 1 -У1У2
= 0 = 0 = 0 = 0
а1б = а1б = 0' а2б = а2б = 0'
где Е1, Е2 — модули упругости вдоль осей Х1 и 71; У1, V2 — коэффициенты Пуассона; G—модуль сдвига.
а^о — а^о — ■
ж/т
0.030
0.015
У, м
0.000
ж/т
0.030
0.015
Y, м
0.030
0.015
X. м
0.030
0.015
X. м
Рис. 6. Критические напряжения и формы потери устойчивости ортотропной пластинки при жестком защемлении кромок (а) и при шарнирном опирании кромок (б). к а — к/Ъ —1/100; ап = а22 = 0.126 (а), 0.0524 ГПа (б); а = 0
W/W
0.030
0.015
W/W
0.030
0.015
Y, м
W/W
0.030
0.015
Y, м
W/W
0.030
0.015
Y, м
0.000
0.000
0.000
0.000
0.030
0.030
0.015
X, м
0.030
0.015 X, м
0.030
0.015
X, м
Рис. 7. Критические напряжения и формы потери устойчивости ортотропной пластинки на упругом основании. Шарнирное опирание кромок (а, г), жесткое защемление сторон (б, в). а11 = 0.407 (а), 0.428 (б), 0.321 (в), 0.306 ГПа (г); а22 = 0.306 (а), 0.322(б), 0.428 (в), 0.407 ГПа (г); а = 0 (а, б), П2 (в, г)
В случае если оси ортотропии материала X1, Y1 составляют с глобальными осями X, Y угол a (рис. 5), то коэффициенты a11, a 22,..., a66 в (11) преобразуются аналогично коэффициентам d11?..., d33 в формуле (7).
Компоненты деформации ех, еу, у^ связаны с прогибом срединной поверхности ^ соотношениями
d 2 w dx2
d2 w _ d2 w
8 У =- Z ' УХУ =-2 Z
ЭхЭу
(12)
0.030
0.015
К, м
ж/жг
0.030
0.015
Y, м
ж/ж
0.030
0.015
Y, м
ж/ж
0.030
0.015
К м
0.000
0.000
0.000
0.000
0.030
0.030
0.015
X, м
0.030
0.015
X, м
0.030
0.015
X, м
Рис. 8. Критические напряжения и формы потери устойчивости ортотропной пластинки на упругом основании при жестком защемлении кромок. Оси ортотропии Х101Y1 повернуты относительно глобальных осей 0ХУ на угол а, ап = а22 = 0.434 (а, г), 0.518 ГПа (б, в); к! а — к/Ъ —1/200 (а, г), 1/100 (б, в); а—л/ 4 (а, б), -л/4 (в, г)
Входящая в (12) величина z есть расстояние от срединной поверхности пластинки до текущей точки.
Однородное уравнение (10) решается совместно с однородными граничными условиями. Если кромки пластины жестко заделаны, то граничные условия имеют вид
дг дг
х = 0, м = 0, — — 0; х = а, м = 0, — — 0; дх дх
(13)
(14)
дг дг
у = 0, м = 0, — — 0; у = Ь, м = 0, — — 0.
Эу Эу
Уравнения устойчивости (10) совместно с граничными условиями (13), (14) представляют собой систему
однородных дифференциальных уравнений в частных производных с однородными граничными условиями. Условие нетривиальности решения системы однородных уравнений определяет критическую нагрузку пластинки на упругом основании под действием сжимающих температурных напряжений ах, ау.
Таким образом, вопрос нахождения критической нагрузки для пластинки, лежащей на упругом основании, сводится к вопросу о нахождении собственных чисел дифференциального уравнения в частных производных. В работе эта проблема решалась с применением метода конечных разностей. Заменяя в уравнении (10) дифференциальные операторы их конечно-разностными аналогами во всех узлах конечно-разностной сетки, получим однородную систему линейных алгебраических уравнений, в которой число уравнений меньше числа неизвестных [39, 40]. Заменяя в граничных условиях (13), (14) дифференциальные операторы их конечно-разностными аналогами, получим систему линейных алгебраических уравнений, которая добавляется к уже существующей. В итоге после дискретизации уравнений устойчивости и граничных условий получается однородная система линейных алгебраических уравнений, в которой число уравнений равно числу неизвестных. Для существования нетривиального решения этой системы необходимо, чтобы ее определитель был равен нулю.
В настоящей работе для вычисления определителя системы линейных алгебраических уравнений использовался метод Гаусса [42]. Следует отметить, что для вычисления значения определителя в процедуре Гаусса достаточно выполнить прямой ход, обратный ход делать нет необходимости. Элементы матрицы, стоящие на главной диагонали после прямого хода, перемножали, и получали значение определителя системы, соответствующее определенным значениям сжимающих напряжений. Значение напряжений, при котором определитель равен нулю, принимали критическим. На практике вычисляли не величину определителя, а его знак для различных значений напряжений ах, ау. Интервал напряжений, в котором определитель меняет знак, определяет критическую нагрузку с любой разумной степенью точности. После этого находили собственный вектор системы алгебраических уравнений. Полагалось, что одна компонента (с номером п) собственного вектора равна константе (например 1), другие компоненты собственного вектора находили в результате решения системы п - 1 уравнений. Таким образом, находили собственный вектор системы линейных алгебраических уравнений (форма потери устойчивости) с точностью до постоянного множителя.
6. Примеры расчета устойчивости ортотропной пластинки на упругом основании
Для расчета использовалась конечно-разностная сетка, содержащая 3600 ячеек и 3721 узел. При этом могут
быть использованы напряжения ах и а у , найденные для узлов сетки конечных элементов. Следует отметить, что при дискретизации уравнений устойчивости анизотропной пластинки (10) методом конечных разностей получается система алгебраических уравнений с несимметричной матрицей. При использовании метода конечных элементов матрица системы алгебраических уравнений симметрична всегда. Технология решения системы с несимметричной матрицей изложена в [43].
На рис. 6 приведены критические нагрузки и формы потери устойчивости пластинки (покрытия) при одновременном сжатии в двух направлениях X и У. Предполагалось, что в докритическом состоянии напряжения ах и ау равны друг другу. Алгоритм нахождения критических напряжений изменяется незначительно, если соотношение напряжений ах и ау будет произвольным.
При решении задачи устойчивости ставились два вида граничных условий: шарнирного опирания и жесткого защемления. Коэффициент постели упругого основания в случаях, рассмотренных на рис. 6, равен нулю. Расхождение между критическими напряжениями, полученными численно в данной работе и аналитическими, составляет 5 %. Аналитическое решение нахождения критических напряжений ортотропной пластинки при шарнирном опирании по краям приведено в [41, 44].
Критические нагрузки и формы потери устойчивости ортотропной пластинки на упругом основании приведены на рис. 7. В качестве граничных условий использовано шарнирное опирание или жесткое защемление по всем кромкам. Отношение толщины пластинки к длине стороны вдоль оси X или У равно И/а - И/Ь --1/200. G = 80 ГПа, модули упругости Е1 - 380 ГПа, Е2 - 190 ГПа, коэффициенты Пуассона v1 - 0.30, V2 -= 0.15. Коэффициент постели К = 1010 Н/м3, коэффициенты линейного температурного расширения а1 -= 6-10-6К-1, а2 - 9•Ю-6К-1.
Из рис. 7 нетрудно видеть, что тепловой удар вызывает потерю устойчивости термопокрытия. Форма потери устойчивости напоминает шахматную доску, причем
Рис. 9. Нормальный прогиб Ж в средине пластинки при потере устойчивости от воздействия температуры. Граничные условия —жесткое защемление по всем кромкам. И а - И/Ь -1/100
0
Рис. 10. Прогибы W ортотропной пластинки при нагреве Отношение к а — к/Ъ —1/100; а11 — напряжения вдоль оси Упругое основание отсутствует
вдоль оси ортотропии с большим модулем упругости наблюдаются три полуволны, вдоль оси с меньшим модулем укладывается пять полуволн. При повороте осей ортотропии на 90° соответственно поворачивается и картина волнообразования при потере устойчивости.
0.000
0.000
0.000
0.000
[ потере устойчивости. Кромки пластинки жестко защемлены. а99 — напряжения вдоль оси у, а19 — сдвиговые напряжения.
На рис. 8 приведены критические напряжения и формы потери устойчивости в случае, когда оси ортотропии 01X1У1 повернуты относительно глобальных осей 0ХУ на 45°, коэффициент постели К = 1010 Н/м3. В этом случае наблюдается анизотропия, матрица упру-
0.000
Рис. 11. Прогибы пластинки без упругого основания (а) и 1010 Н/м3 (б)
гих характеристик [D] полностью заполнена. Из рис. 8 нетрудно видеть, что картина волнообразования симметрична относительно осей ортотропии XlOlYl. При повороте осей ортотропии на 90° картина волнообразования поворачивается также на 90°. Если отношение толщины пластинки к длине ее стороны равно Н/a = = 0.01, то вдоль одной оси ортотропии число полуволн равно 1, вдоль другой оси — 3. Если отношение толщины пластинки к длине стороны равно На = 1/200, то вдоль одной оси ортотропии число полуволн равно 4, вдоль другой оси — 3.
7. Напряженно-деформированное состояние анизотропной пластинки после потери устойчивости в результате теплового удара
В результате нагрева ортотропная пластинка увеличивает свои размеры вдоль осей X и Y на AZ1 и AL2, причем AL =a1 ATL10, AL2 = а2 ATL20. На рис. 9 представлен нормальный прогиб посредине пластинки W (Х=0.015м), вычисленный из условия, что длина кривой ABC равна L0 + a1 ATL0:
60 I-2-2
Lo +aiATL,o = EV(+i - )2 + (м+i - w,.)2. (15)
i=1
Параметры напряженно-деформированного состояния ортотропной пластинки при прогибе W, вычисленном исходя из уравнения (12), приведены на рис. 10.
На рис. 11 приведены прогибы одной клетки «шахматной доски», возникающей при потере устойчивости
0-1- ' ——^
0 0.005 0.010
на упругом основании (б). Коэффициент постели K = 0 (а) и
пластинки. Нетрудно заметить, что упругое основание уменьшает прогибы пластинки примерно в 2 раза. Параметры напряженно-деформированного состояния пластинки на упругом основании после потери устойчивости приведены на рис. 12. Можно отметить, что уменьшение толщины покрытия (пластины) и увеличение жесткости (модуля упругости) подложки приводит к росту числа полуволн при фиксированных размерах расчетной области и уменьшению их амплитуд (рис. 12). Величины напряжений, возникающих за счет потери устойчивости, зависят не столько от абсолютных значений величины прогиба, сколько от кривизны поверхности покрытия, возникающей после потери устойчивости, т.е. от вторых производных величины прогиба по осям вдоль поверхности. Кроме того, напряжения зависят от толщины пластины. В этом смысле толщина влияет на максимальные значения напряжений двояким образом: с одной стороны, с уменьшением толщины покрытия увеличение числа волн при потере устойчивости ведет к росту кривизны, с другой — мо-ментные напряжения уменьшаются, поскольку уменьшается так называемая цилиндрическая жесткость пластины.
Это обстоятельство можно использовать для проведения параметрических исследований и определения рациональных значений толщины покрытия в зависимости от деформационно-прочностных и теплофизи-ческих характеристик материалов покрытия и подложки.
аш МПа 600
а22, МПа 600
а12, МПа 200
ж, мкм 100
0.03
X, м
0.03
X, м
0.03
X, м
0.03
X, м
Рис. 12. Напряжения вдоль осей X, Y, сдвиговые напряжения в пластинке на упругом основании при потере устойчивости при нагревании. Прогиб пластинки после потери устойчивости. к а — к/Ъ —1/200. Коэффициент постели К = 1010 Н/м3
8. Выводы
При тепловом ударе в начальные моменты времени, когда поцложка находится еще в исходном непрогретом состоянии, в термопокрытии возникают значительные сжимающие напряжения, которые могут служить причиной потери его устойчивости как тонкостенной пластинки, моделирующей это покрытие, находящейся на упругом основании (подложке).
При потере устойчивости покрытия экстремальные значения напряжений расположены в порядке, напоминающем шахматную доску. Если физико-механические характеристики покрытия анизотропны, меняется общая картина потери устойчивости по сравнению с изотропным покрытием. «Клетки шахматной доски» вытягиваются вдоль оси, которой соответствует больший модуль упругости.
При потере устойчивости пластинки на упругом основании число волн на поверхности покрытия увеличивается с увеличением жесткости упругого основания и уменьшением отношения толщины пластинки к ее ширине и длине. При проведении расчетов для конкретного ортотропного покрытия показано, что для квадратной пластины вдоль оси с большим модулем упругости располагается меньшее число полуволн, нежели вдоль оси ортотропии с меньшим модулем упругости.
При потере устойчивости анизотропной пластинки, оси ортотропии которой повернуты относительно глобальных осей на 45°, выпуклости и вмятины располагаются симметрично относительно осей ортотропии.
Работа выполнена в рамках Программы фундаментальных научных исследований государственных академий наук на 2013-2020 гг., п. III.20.1.3, а также гранта Президента РФ по поддержке ведущих научных школ № НШ-2817.2014.1.
Литература
1. Padture N.P., Gell M, Jordan E.H. Thermal barrier coatings for gasturbine engine applications // Science. - 2002. - V. 296. - P. 280284.
2. Clarke D.R., Levi C.G. Materials design for the next generation thermal barrier coatings // Annu. Rev. Mater. - 2003. - V. 33. - P. 383417.
3. AhmedF., Schmid C., Schaufler J., Durst K. Study on the deformation mechanics of hard brittle coatings on ductile substrates using in-situ tensile testing and cohesive zone FEM modeling // Surf. Coat. Tech. -2012. - V. 207. - P. 163-169.
4. Fukumasu N.K., Angelo C.M., Ignat M, Souza R.M. Numerical study of tensile tests conducted on systems with elastic-plastic films deposited onto elastic-plastic substrates // Surf. Coat. Tech. - 2010. -V. 205.- P. 1415-1419.
5. Faou J.Y., Parry G., Grachev S., Barthel E. How does adhesion induce the formation of telephone cord buckles // Phys. Rev. Lett. -2012. - V. 108. - P. 116102.
6. Chen B.F., Hwang J., Chen I.F., Yu G.P., Huang J.H. A tensile-film-cracking model for evaluating interfacial shear strength of elastic film on ductile substrate // Surf. Coat. Tech. - 2000. - V. 126. - P. 91-95.
7. Evans A.G., Mumm D.R., Hutchinson J.W., Meier G.H., Pettit F.S. Mechanisms controlling the durability of thermal barrier coatings // Prog. Mater. Sci. - 2001. - V. 46. - P. 505-553.
8. Mao W.G., Dai C.Y., Zhou Y.C., Liu Q.X. An experimental investigation on thermo-mechanical buckling delamination failure characteristic of air plasma sprayed thermal barrier coatings // Surf. Coat. Tech. -2007. - V. 201. - P. 6217-6227.
9. Qian L.H., Zhu S.J., Kagawa Y., Kubo T. Tensile damage evolution behavior in plasma-sprayed thermal barrier coating system // Surf. Coat. Tech. - 2003. - V. 173. - P. 178-184.
10. Schlichting K.W., Padture N.P., Jordan E.H., Gell M. Failure modes in plasma sprayed thermal barrier coatings // Mater. Sci. Eng. A. -2003. - V. 342. - P. 120-130.
11. Zhao X., Wang X., Xiao P. Sintering and failure behaviour of EB-PVD thermal barrier coating after isothermal treatment // Surf. Coat. Tech. - 2006. - V. 200. - P. 5946-5955.
12. Xie X.Y., Guo H.B., Gong S.K., Xu H.B. Thermal cycling behavior and failure mechanism of LaTi2Al9O19/YSZ thermal barrier coatings exposed to gas flame // Surf. Coat. Tech. - 2011. - V. 205. -P. 4291-4298.
13. Karlsson A.M., Evans A.G. A numerical model for the cyclic instability of thermally grown oxides in thermal barrier systems // Acta Mater. -2001. - V. 49. - P. 1793-1804.
14. Rabiei A., Evans A.G. Failure mechanisms associated with the thermally grown oxide in plasma-sprayed thermal barrier coatings // Acta Mater. - 2000. - V. 48. - P. 3963-3976.
15. Erdogan F., Ozturk M. Periodic cracking of functionally graded coatings // Int. J. Eng. Sci. - 1995. - V. 33. - P. 2179-2195.
16. Schulze G.W., Erdogan F. Periodic cracking of elastic coatings // Int. J. Solids Struct. - 1998. - V. 35. - P. 3615-3634.
17. Erdogan F. The crack problem for bonded nonhomogeneous materials under antiplane shear loading // J. Appl. Mech. - 1985. - V. 52. -P. 823-828.
18. Agrawal D.C., Raj R. Measurement of the ultimate shear strength of a metal-ceramic interface // Acta Metall. - 1989. - V. 37. - P. 12651270.
19. McGuigan A.P., Briggs G.A.D., Burlakov V.M., Yanaka M., Tsukaha-ra Y An elastic-plastic shear lag model for fracture of layered coatings // Thin Solid Films. - 2003. - V. 424. - P. 219-223.
20. Aktaa J., SfarK., MunzD. Assessment ofTBC systems failure mechanisms using a fracture mechanics approach // Acta Mater. - 2005. -V. 53. - P. 4399-4413.
21. Wu C.W., Chen G.N., Zhang K., Luo G.X., Liang N.G. The effect of periodic segmentation cracks on the interfacial debonding: study on interfacial stresses // Surf. Coat. Tech. - 2006. - V. 201. - P. 287-291.
22. Fan X.L., Xu R., Zhang W.X., Wang T.J. Effect of periodic surface cracks on the interfacial fracture of thermal barrier coating system // Appl. Surf. Sci. - 2012. - V. 258. - P. 9816-9823.
23. Yang L., Zhong Z.C., Zhou Y.C., Lu C.S. Quantitative assessment of the surface crack density in thermal barrier coatings // Acta Mech. Sinica. - 2014. - V. 30. - P. 167-174.
24. Zhou B., Kokini K. Effect of pre-existing surface crack morphology on the interfacial thermal fracture of thermal barrier coatings: a numerical study // Mater. Sci. Eng. A. - 2003. - V 348. - P. 271-279.
25. Zhou B., Kokini K. Effect of preexisting surface cracks on the interfacial thermal fracture of thermal barrier coatings: an experimental study // Surf. Coat. Tech. - 2004. - V. 187. - P. 17-25.
26. Fan X.L., Zhang W.X., Wang T.J., Sun Q. The effect of thermally grown oxide on multiple surface cracking in air plasma sprayed thermal barrier coating system // Surf. Coat. Tech. - 2012. - V. 208. - P. 7-13.
27. Fan X.L., Zhang W.X., Wang T.J., Liu G.W., Zhang J.H. Investigation on periodic cracking of elastic film/substrate system by the extended finite element method // Appl. Surf. Sci. - 2011. - V 257. - P. 67186724.
28. Ranjbar-Far M., Absi J., Mariaux G., Smith D.S. Crack propagation modeling on the interfaces of thermal barrier coating system with different thickness of the oxide layer and different interface morphologies // Mater. Design. - 2011. - V. 32. - P. 4961-4969.
29. Zhang W.X., Fan X.L., Wang T.J. The surface cracking behavior in air plasma sprayed thermal barrier coating system incorporating interface roughness effect // Appl. Surf. Sci. - 2011. - V. 258. - P. 811— 817.
30. Xu R., Fan X.L., Zhang W.X., Song Y., Wang T.J. Effects of geometrical and material parameters of top and bond coats on the interfacial fracture in thermal barrier coating system // Mater. Design. - 2013. -V. 47. - P. 566-574.
31. Mao W.G., Dai C.Y., Yang L., Zhou Y.C. Interfacial fracture characteristic and crack propagation of thermal barrier coatings under tensile conditions at elevated temperatures // Int. J. Fract. - 2008. -V. 151.- P. 107-120.
32. Черепанов О.И. Численное решение некоторых квазистатических задач мезомеханики. - Новосибирск: Изд-во СО РАН, 2003. -
180 c.
33. Cherepanov G.P. On the theory of thermal stresses in a thin bonding layer // J. Appl. Phys. - 1995. - V. 78. - P. 6826-6832.
34. Панин B.E., Сергеев В.П., Панин A.B., Почивалов Ю.И. Нано-структурирование поверхностных слоев и нанесение нанострук-турных покрытий — эффективный способ упрочнения современных конструкционных и инструментальных материалов // ФММ. -2007. - Т. 104. - № 6. - С. 650-660.
35. Zienkiewicz O.C., Morgan K. Finite Elements and Approximation. -Dover Publications, 2006. - 328 p.
36. SegerlindL.J. Applied Finite Element Analysis. - New York: John Wiley and Sons, 1984. - 427 p.
37. Бате К.-Ю. Методы конечных элементов / Под ред. Л.И. Турча-ка. - М.: Физматлит, 2010. - 1024 с.
38. Прочность, устойчивость, колебания. Справочник в 3-х т. Т. 3. -М.: Машиностроение, 1968. - 567 с.
39. Тимошенко С.П., Войновский-Кригер С. Пластинки и оболочки. -М.: Наука, 1966. - 636 с.
40. Тетерс Г.А., Рикардс Р.Б., Нарусберг В.Л. Оптимизация оболочек из слоистых композитов. - Рига: Зинатне, 1978. - 240 с.
41. Лехницкий С.Г. Анизотропные пластинки. - М.: ОГИЗ, Гостех-издат, 1947. - 355 с.
42. Zienkiewicz O.C., Taylor R.L., Zhu J.Z. The Finite Element Method: Its Basis and Fundamentals. - Butterworth-Heinemann, 2005.
43. Сосис П.М. Алгоритмический язык АЛГ0Л-60 и применение его в строительной механике. - Киев: Бущвельник, 1965. - 172 с.
44. Абовский Н.П., Андреев Н.П., Деруга А.П. и др. Численные методы в теории упругости и теории оболочек. - Красноярск, 1986. - 384 с.
Поступила в редакцию 24.03.2015 г.
Сведения об авторах
Люкшин Петр Александрович, к.ф.-м.н., снс ИФПМ СО РАН, petrljuk@ispms.tsc.ru
Люкшин Борис Александрович, д.т.н., внс ИФПМ СО РАН, зав. каф. ТУСУР, lba2008@yandex.ru
Матолыгина Наталья Юрьевна, к.ф.-м.н., нс ИФПМ СО РАН, ksa@ispms.tsc.ru
Панин Сергей Викторович, д.т.н., зам. дир. ИФПМ СО РАН, проф. ТПУ, svp@ispms.tsc.ru
