Трехмерное моделирование методом подвижных клеточных автоматов упругопластического деформирования и разрушения покрытий при контактном взаимодействии с жестким индентором

УДК 004.94:539.62

Трехмерное моделирование методом подвижных клеточных автоматов упругопластического деформирования и разрушения покрытий при контактном взаимодействии с жестким индентором

А.Ю. Смолин12, Г.М. Еремина12, В.В. Сергеев13, Е.В. Шилько12, С.Г. Псахье123

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

Three-dimensional movable cellular automaton simulation of elastoplastic deformation and fracture of coatings in contact interaction with a rigid indenter

A.Yu. Smolin12, G.M. Eremina12, V.V. Sergeev13, E.V. Shilko12, and S.G. Psakhie123

In the work, the movable cellular automaton method is applied to consider 3D models of contact interaction such as nanoindentation, sclerometry (scratching), and tribospectroscopy. The system under study is a titanium substrate with a hardening coating. The substrate and coating materials are both described in the elastoplastic approximation. It is shown that the method chosen for numerical simulation is capable of correctly describing the contact interaction of elastoplastic materials under loads of different types with explicit allowance for fracture. The possibility to detect damages in material surface layers from friction force estimates is assessed.

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

Нанесение на поверхность материала специальных покрытий является эффективным способом повышения их функциональных свойств. Например, для имплан© Смолин А.Ю., Еремина Г.М., Сергеев В.В., Шилько Е.В., Псахье С.Г., 2014

татов используются многокомпонентные биоактивные наноструктурные покрытия на основе тугоплавких соединений ТЮ, (Т^Та)(С,К) с добавлением специальных элементов (Са, Zr, Si, О, Р), которые улучшают как трибологические, так и биоактивные свойства поверхности [5-8]. Следует отметить, что эти покрытия имеют наноструктурное состояние и малую толщину.

Для изучения механических свойств покрытий и пленок в основном используются такие методы, как наноиндентирование [9-12] и склерометрия [13-15]. Наноиндентирование представляет собой контролируемое внедрение сверхтвердого наконечника (индентора) в поверхность материала. В склерометрических испытаниях (царапание) после внедрения индентора осуществляется его движение вдоль покрытия, что приводит к его разрушению. Одним из методов неразрушающего контроля качества покрытия, позволяющим определить наличие в нем дефектов, является трибоспектроскопия [16], где в качестве измеряемого и анализируемого параметра используется сила трения при движении малого контртела вдоль поверхности.

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

К настоящему времени опубликовано довольно много работ, посвященных численному моделированию процессов наноиндентирования и измерительного царапания. В зависимости от используемого метода они могут быть разделены на две группы. В работах первой группы используется метод конечных элементов и рассматривается поведение материала на макромасштабе. В этих работах изучаются особенности распределения напряжений в материале при различных значениях параметров системы, таких как форма индентора, механические свойства материала и т.д. [17-20]. Вторая группа работ посвящена применению метода молекулярной динамики (метода частиц) для изучения на микроуровне механизмов зарождения пластической деформации в непосредственной близости от вершины индентора [2124].

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

который позволил бы моделировать на различных масштабах как процесс упругопластического деформирования, так и разрушение твердых тел. Как известно, лучшими возможностями для моделирования разрушения, в том числе зарождения и роста трещин, фрагментации и перемешивания вещества, обладают методы частиц, берущие свое начало из метода молекулярной динамики. Следует подчеркнуть, что из всех таких методов только метод подвижных клеточных автоматов (метод MCA) способен корректно описывать пластическую деформацию консолидированных тел [25], что и обусловило его выбор для построения модели.

Метод подвижных клеточных автоматов [25-29] является новым эффективным численным методом, основанным на концепции дискретных элементов (частиц), которая существенно отличается от концепции численных методов классической механики сплошных сред. Следует отметить, что развитие дискретных методов началось с метода молекулярной динамики, который был разработан для изучения материалов на атомном уровне. Однако возможности атомного описания на пространственных и временных масштабах, которые представляют интерес для инженерных приложений, существенно ограничены. Это привело к развитию методов, подобных методу молекулярной динамики для мезо- и макромасштабных уровней (методов частиц), в которых структурные элементы имеют конечный размер (в отличие от атомов, которые являются точечными массами) и взаимодействуют только с ближайшими соседями. Самым известным представителем этой группы методов является метод дискретных элементов [3034], который в настоящее время широко используется для изучения механического поведения гранулированных (сыпучих) и слабосвязанных сред, в частности реологических особенностей этих систем, их разрушения и перемешивания [32, 34]. В то же время до недавнего времени применимость метода дискретных элементов для изучения механического поведения консолидированных тел была ограничена главным образом хрупкими пористыми материалами [30-32, 34] в связи с недостаточным развитием математических моделей для вычисления сил взаимодействия дискретных элементов. В частности, большинство моделей дискретных элементов используют парные (двухчастичные) потенциалы (или силы) взаимодействия. Это упрощение приводит к ряду искусственных эффектов в поведении ансамбля частиц, таких как анизотропия, навязанное упаковкой значение модуля сдвига и коэффициента Пуассона и т.д. [30, 34]. Среди проблем, возникающих вследствие этих искусственных эффектов, важно отметить неадекватность моделирования накопления необратимых деформаций (пластичности) материалов.

Последние исследования показали, что многие проблемы метода дискретных элементов, в том числе касающиеся описания консолидированных твердых тел на различных масштабах, могут быть решены с помощью использования многочастичных сил взаимодействия между элементами. В [25] описан обобщенный подход к построению многочастичных сил взаимодействия для дискретных элементов, который основан на идее, применяемой для записи межатомных потенциалов в методе погруженного атома. Так, выражение для силы, действующей на дискретный элемент i, имеющий Ni соседей, может быть записано в виде:

Ц = 22 Ца1Г + Ц°.

Эта сила представлена как суперпозиция парных компонентов Цра1Г, зависящих от пространственного расположения автомата i по отношению к соседу j, и объемнозависящей компоненты обусловленной коллективными эффектами окружения.

В рамках метода подвижных клеточных автоматов предполагается, что любой материал состоит из определенного количества элементарных объектов конечного размера (автоматов), которые взаимодействуют друг с другом и могут перемещаться в пространстве, тем самым моделируя реальные процессы деформации. Движение автоматов описывается уравнениями Ньютона-Эйлера:

г Н2К N

т = 2 ЦТ + Ц°,

I а20 N. (!)

дт з=1

где К,, 9,, т{ и 31 — радиус-вектор, вектор вращения, масса и момент инерции &-го автомата соответственно; ЦРа1Г — парная сила механического взаимодействия автоматов i и j; — объемнозависящая сила, действующая на автомат i и обусловленная взаимодействием его соседей с другими автоматами. В последнем

уравнении М у = qij (п 3 х ЦРа1Г) + К 3, где qij — расстояние от центра /&-го автомата до точки его взаимодействия (контакта) с ]&-м автоматом; п 3 = (Я3 - Я,)Гу — единичный вектор ориентации пары и г3 — расстояние между центрами автоматов (рис. 1); К 3 — крутящий момент, обусловленный относительным вращением автоматов в паре (см. ниже).

Отметим, что автоматы пары могут представлять собой части различных тел или одного консолидированного тела. Поэтому их взаимодействие не всегда является контактным в том смысле, который имеется в виду в известных задачах механики [2-4]. Так, в случае автоматов одного консолидированного тела площадка контакта автоматов будет представлять собой просто элемент плоскости, разделяющей два конечных объема материала, аналогичный граням ячеек разностной сетки в численных методах континуальной механики. По этой причине здесь и далее мы будем слово «контакт» брать в кавычки. Кроме того, как показано на рис. 1, а, размер автомата характеризуется одним параметром di, но это не означает, что форма автомата представляет собой шар. Реальная форма автомата определяется областью его «контактов» с соседями. Например, если использовать в качестве начальной ГЦК-упаковку, то автоматы будут иметь форму ромбического додекаэдра, а если кубическую, то автоматы будут представлять собой кубики.

Для локально изотропных сред объемнозависящая компонента силы Ц может быть записана через давление Р3 в объеме соседнего автоматаj следующим образом:

Ц?=-А ЪР^ П3,

где SiJ — площадь контакта &-го автомата с j-м; А — некий материальный параметр.

С другой стороны, общая сила, действующая на автомат, может быть представлена в виде суммы нормальной ЦП и касательной (сдвиговой) ЦТ компонент:

F. =£ (Fj311 - APS nj) =

= E [(FT&" (hj) - APjSj )nj + F?^ (lfar)tj ] =

= E (F" + Fj),

где ^?а1Г&и — нормальная, а рра1Г>т — касательная силы

взаимодействия в паре, зависящие соответственно от межавтоматного перекрытия Ь- (рис. 1, а) и относительного тангенциального смещения 1-1еаг (рис. 1, б), рассчитанного с учетом вращения обоих автоматов [28, 29]. Отметим что, несмотря на то, что последнее выражение в уравнении (2) формально соответствует обычной записи силы взаимодействия в методе дискретных элементов [33-34], оно принципиально отличается от последнего вследствие многочастичного характера центрального взаимодействия автоматов.

Используя процедуру осреднения для тензора напряжений в частице, изложенную в [25, 34], выражение для компонент усредненного тензора напряжений в автомате i принимает вид:

где а, X, Y, 2 лабораторной системы координат; Vi — текущий объем автомата г; п- а — а-компонента единичного вектора п -; Р- р — Р-компонента полной силы, действующей в точке «контакта» между автоматами

i и j. _

Давление Ц, или среднее напряжение а^ит, в объеме автомата может быть вычислено через компоненты тензора напряжений:

Знание компонент тензора напряжений позволяет вычислять все его инварианты в объеме автомата, в частности интенсивность напряжений:

= 1/V2 ((?„ - ^ )2 + tfyy - aZz )2 +

+ (^ -аХх)2 + 6[(аХу )2 + (а^)2 + (aixz)2])1/2. (5) Из уравнений (1)-(3) следует, что выражения для выи F?™&т определяют реологические

числения Fifair&;

свойства модельной среды.

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

и j характеризуется соответствующими напряжениями % и т\\ =ПуЯу ,

& г (6)

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

Деформацию автомата i при его центральном взаимодействии с автоматом j можно характеризовать следующей безразмерной величиной (нормальной деформацией):

q.j- dil2 dj 2 &

В общем случае каждый из автоматов представляет собой различные материалы, и перекрытие в паре некоторым образом перераспределяется между г&-м иj-м автоматами:

АЬ- = Aqlj+ Aq- = А-, /2 + А^-, dJ /2, (8)

где символ А означает приращение величины за временной шаг At численного интегрирования уравнений движения (1). Правило перераспределения деформации в паре тесно связано с выражением для вычисления силы взаимодействия автоматов. Это выражение для центрального взаимодействия похоже на соотношения Гука для диагональных компонент тензора напряжений:

А- 2С(А— + (1 - 20/К) Р,, (9)

где К — модуль объемной упругости; G — модуль сдвига материала автомата г; Ц — давление автомата г, которое может быть вычислено с помощью формул (3) и (4) на предыдущем шаге по времени или по схеме предиктор-корректор.

Для определения параметра, характеризующего деформацию сдвига в паре автоматов г^, мы воспользуемся формулой кинематики для свободного движения пары как абсолютно жесткого (недеформируемого) тела: vj - V = х г- > (10)

где г- = (R- - Rг■), V, — трансляционная скорость автомата г&; ю- — мгновенная скорость вращения пары как целого. Если мы умножим обе части уравнения (10) векторно слева на г- и пренебрежем вращением вокруг оси пары (т.е. пусть ю- • г- = 0, т.к. вращение вокруг оси пары не производит деформацию сдвига), то мы получим следующую формулу:

Г- Х (V- - Vi ) = П г- Х (V - - Vi )

Кроме такого вращения пары как целого (определяется разницей в поступательных скоростях автоматов), каждый из автоматов вращается со своей собственной угловой скоростью ю, (рис. 1, в). Разница между этими вращениями обусловливает деформацию сдвига. Так, приращение деформации сдвига автоматов г и ]

за шаг Дt определится относительным тангенциальным

смещением в точке «контакта» пары Д131еаг, деленной на расстояние между автоматами:

ДУ а +дУ ц = —

(Чц («ц - Щ)х па + q а («а - « а)х пц)д

Выражение для вычисления силы тангенциального взаимодействия подвижных клеточных автоматов похоже на соотношения Гука для недиагональных компонент тензора напряжений:

Дх„. = 2С (Ду ,3) (13)

и является чисто парным.

Разные скорости вращения каждого из автоматов также приводят к деформации относительного «изгиба» и «кручения» (последнее присутствует только в 3Э) в паре. Очевидно, что сопротивление относительному вращению в паре вызывает крутящий момент, величина которого пропорциональна разнице между поворотами автоматов:

К 3 +GJ )(9 3-9,). (14)

Формулы (1)-(4), (6)-(9), (11)—(14) описывают механическое поведение линейно-упругого тела в рамках метода подвижных клеточных автоматов. Отметим, что соотношения (8), (9), (12) и (13) записаны в приращениях, т.е. в гипоупругой форме. В работе [35] показано, что эта модель дает результаты, соответствующие численному решению методом конечных разностей обычных уравнений механики сплошной среды для изотропной линейно-упругой среды. Это обстоятельство позволило использовать метод подвижных клеточных автоматов совместно с численными методами механики сплошных сред в рамках единого дискретно-континуального подхода. В работе [29] показано, что именно учет вращения автоматов совместно с многочастичной формой центрального взаимодействия позволяет методом подвижных клеточных автоматов адекватно описывать изотропный отклик материала на приложенную нагрузку.

Для описания упругопластического поведения в рамках метода подвижных клеточных автоматов предлагается использовать теорию пластического течения, а именно модель идеальной пластичности с критерием Мизеса. Для этого к методу подвижных клеточных автоматов [25] был адаптирован известный алгоритм Уил-кинса [36, 37]. Этот алгоритм состоит в решении упругой задачи на каждом временном шаге и последующем «сбросе» компонент девиатора тензора напряжений Дар = аар -1/3 акк5ар на поверхности текучести Мизеса в случае, когда интенсивность напряжений превышает заданную предельную величину (рис. 2):

дав = Дав М,

Рис. 2. Схема работы алгоритма Уилкинса

где М = ар1/ а1п1 , аы — интенсивность напряжений; ар1 — радиус круга текучести Мизеса.

Этот алгоритм в применении к автомату & может быть записан в следующих обозначениях:

I (ааа ) = (ааа - атеап )М, + атеап &

Ыр )& = 50р М,, (15)

где а, в = X, Y, 2 и а ф в; (а^р)& — скорректированные компоненты осредненного тензора напряжений; агар — компоненты тензора напряжений, полученные в результате решения упругой задачи на текущем временном

шаге; М{ = 0^/ а[п1 — текущее значение коэффициента «сброса»; ар1 — текущее значение радиуса круга текучести Мизеса для автомата

Удельные величины нормальных и касательных сил взаимодействия корректируются с использованием текущего значения коэффициента М по аналогии с алгоритмом Уилкинса для компонент осредненного тензора напряжений (15):

К = (% - Р) М1 + Р,

{ т 3 = т 3М,,

где и т 3 представляют собой скорректированные значения удельных сил (напряжений). Легко показать, что подстановка этих соотношений в выражения (2) и (3) автоматически обеспечивает корректировку компонент осредненного тензора напряжений в автомате к кругу текучести [25].

Таким образом, реологические свойства материала автомата & определяются заданием единой кривой упрочнения а|п1 =Ф(ё^) (здесь ёП — интенсивность осредненного тензора деформаций, компоненты которого могут быть вычислены аналогично [25, 34]), эта зависимость также называется функцией отклика автомата.

Для численного интегрирования уравнений движения (1) можно использовать схему Верле в скоростной форме, модифицированную введением предиктора для оценки а,ар на текущем шаге времени.

Пару элементов можно рассматривать как виртуальный бистабильный автомат (у него существуют два состояния: связанная и несвязанная пара), что позволяет явно моделировать процессы разрушения в методе подвижных клеточных автоматов. Заданием правил перехода пары из состояния связанной в состояние несвязанной формулируется критерий разрушения моделируемого материала, который, вообще говоря, определяется физическими механизмами деформации материала. Важным преимуществом описанного выше формализма является то, что в нем возможно непосредственное применение классических критериев (Губера-Мизеса-Генки, Друкера-Прагера, Кулона-Мора, Под-горски и т.д.), которые записаны в тензорной форме. Заметим, что переключение пары автоматов в несвязанное состояние может привести к изменению сил, действующих на элементы, в частности, они не будут сопротивляться удалению друг от друга.

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

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

Наноиндентирование — это процесс контролируемого внедрения сверхтвердого наконечника определенной формы (индентора) под действием нарастающей нагрузки в плоскую поверхность образца на глубину менее 100 нм, при этом в процессе нагружения постоянно измеряется сила Р, действующая на индентор, и глубина его погружения в материал h [8]. Иногда наноин-дентирование называют измерительным индентирова-нием. На основе анализа измеряемой Р-^диаграммы можно получать такие характеристики материала, как модуль упругости, упругое восстановление и нанотвер-дость. При таком анализе в основном используется методика Оливера-Фарра [10]. Однако, как показали, например, авторы работы [11], корректно определить модуль упругости и нанотвердость тонких покрытий и пленок по данным наноиндентирования с использованием стандартной методики Оливера-Фарра возможно только при условии совпадения этих характеристик у покрытия и подложки. Очевидно, что в большинстве практически важных приложений это условие не выполняется. Одним из способов решения этой проблемы может быть использование компьютерного моделирования, в рамках которого возможно получение достаточно точных зависимостей для индентирования покрытий на любых подложках.

При разработке модели, а также ее тестирования были использованы результаты натурных экспериментов, опубликованные в работах [5-8, 38-40]. Так, геометрические характеристики модельной системы «покрытие - переходный слой - подложка» определялись по фотографии поперечного среза реального образца [38]. Полученные по результатам натурного наноинден-тирования значения модуля Юнга и предела текучести учитывались при задании параметров модельных материалов.

Механическое поведение материалов и подложки, и покрытия описывалось в рамках упругопластического приближения. Упруго-идеальнопластическое тело характеризуется следующими параметрами: плотность р, модуль сдвига G, модуль объемной упругости К и предел текучести ау. На основе данных, приведенных в [39, 40], значения указанных параметров для подложки составили р = 4500 кг/м3, G = 41 ГПа, К = 100 ГПа и ау =2 ГПа, а для покрытия р = 4700 кг/м3, G = 76 ГПа, К = 167 ГПа и ау = 15 ГПа.

В экспериментах по наноиндентированию и склерометрии обычно используются алмазные инденторы. Модуль упругости алмаза (1141 ГПа) приблизительно в пять раз превышает модуль упругости для многокомпонентных биоактивных наноструктурных покрытий (200-250 ГПа). Поэтому в данных расчетах материал индентора принимался абсолютно жестким (неде-формируемым). Форма индентора соответствовала трехгранной пирамиде Берковича с углом между гранями и осью пирамиды равным 65.3°.

Рассматривалось индентирование двух типов образцов, каждый из которых имел форму параллелепипеда. Первый тип образцов состоял полностью из нанострук-турного титана, а второй — из нескольких слоев, нижний из которых представлял собой подложку из титана, а верхний — покрытие. Толщина такого покрытия составляет 350-2000 нм [5-8], при этом общая высота образцов варьировалась от 700 до 2500 нм. Между покрытием и титановой подложкой находился промежуточный слой толщиной порядка 50 нм. В модели свойства промежуточного слоя задавались равными средним величинам между значениями для покрытия и подложки. Размер автоматов, исходя из предварительных расчетов [41, 42], брался равным 10 нм.

Общий вид построенной модели представлен на рис. 3, а в виде ГЦК-упаковки сферических частиц. Процесс нагружения имитировался путем задания одинаковой скорости движения для всех автоматов инден-тора. Для того чтобы предотвратить смещение нагружаемого материала как целого, нижний слой автоматов подложки считался жестко закрепленным.

Анализ поля скоростей в начале индентирования, когда материал деформируется упруго, показывает, что смещения в материале локализуются под индентоРис. 3. Общий вид модели (а), распределение скоростей автоматов (б, абсолютные значения в м/с), интенсивности напряжений (Па) по сечению образца при его упругой (в) и пластической (г) деформации

ром и направлены преимущественно вглубь материала (рис. 3, б). Полученные результаты согласуются с данными натурных экспериментов, которые показывают, что результаты индентирования тонкого покрытия при больших нагрузках в значительной степени определяются свойствами подложки [9, 11, 39, 40]. Распределение упругих напряжений в образце (рис. 3, в) качественно соответствует аналитическому решению для задачи Герца [4]. На рис. 3, г показано поле интенсивности напряжений при максимальной глубине проникновения, когда покрытие и подложка деформируются уже пластически. Видно, что максимальные напряжения находятся в непосредственной близости от индентора, однако из-за меньшего предела текучести у подложки она также деформируется пластически. Поля напряжений, представленные на рис. 3, г, находятся в хорошем качественном согласии с литературными данными [43, 44].

Главной характеристикой процесса измерительного индентирования является Р-^диаграмма. В случае малой глубины проникновения эта кривая в основном определяется свойствами покрытия. Чем больше глубина проникновения, тем большее влияние на кривую оказывают свойства подложки. Поэтому вначале проводилось моделирование погружения индентора в образец из титана, затем — в образец с покрытием на небольшую глубину и в заключение — индентирование на значительные глубины проникновения в образец с покрытием.

На рис. 4 представлены расчетная и экспериментальная [39] Р-^диаграммы индентирования титанового образца на глубину 650 нм. Их анализ позволяет сделать вывод о возможности нахождения свободных параметров модели для описания свойств материала

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

Для определения адекватности параметров модельного покрытия был сделан расчет с глубиной внедрения в образец с покрытием, равной 60 нм [42]. Результаты этого расчета хорошо соответствуют экспериментальным данным при максимальной нагрузке на индентор 2 мН [39, 40]. Для изучения влияния подложки на Р-^ диаграмму моделировалось погружение на большую глубину. В связи с этим использовались образцы с общей высотой 2.5 и 5.0 мкм и размером автоматов 70 и 140 нм соответственно, толщина покрытия составляла 1000 нм. Р-^диаграммы, полученные в результате моР, мЩ А

Рис. 4. Экспериментальная [39] (1) и расчетная (2) диаграммы индентирования титанового образца

Р, мН " 80 ■ 60 ■ 40 ■ 20 ■

Рис. 5. Экспериментальные [39, 40] (1, 3) и расчетные (2, 4) диаграммы индентирования титанового образца с нанострук-турным покрытием

делирования и эксперимента [39, 40], показаны на рис. 5. Согласие результатов расчетов с имеющимися экспериментальными данными свидетельствуют о том, что предложенная модель хорошо описывает особенности деформационного поведения таких систем.

Методика склерометрических испытаний (измерительного царапания) [13-15] основана на непрерывном контролируемом нагружении поверхности материала путем перемещения индентора, предварительно внедренного на определенную глубину, вдоль поверхности. При этом происходит деформирование материала вначале в упругой, а затем и в упругопластической областях до предельного состояния и последующее его разрушение. Схемы, реализующие методики индентирования и склерометрии, принципиально устроены сходным образом. Существенным отличием склерометрии является то, что при этом регистрируется перемещение инден-тора не только по нормали, но и вдоль поверхности, что дает возможность измерять тангенциальные силы, а следовательно, и коэффициент трения. Момент адгезионного или когезионного разрушения покрытия определяют при визуальном анализе микрофотографий, полученных с помощью оптического микроскопа, оборудованного цифровой камерой, а также по изменению следующих пяти параметров: акустической эмиссии, силы трения, коэффициента трения, глубины проникновения индентора и остаточной глубины царапины.

Модель процесса склерометрии (рис. 6, а) в целом подобна модели для наноиндентирования и состоит из четырех блоков: основание, подложка, переходный слой, покрытие и индентор (контртело). Образец исследуемого материала, состоящий из трех материалов, имеет геометрическую форму параллелепипеда, высота покрытия 1.5 мкм, переходного слоя — 0.2 мкм и подложки из титана — 2.5 мкм. Основание представляет собой пластину со сторонами равными сторонам образца и высотой равной размеру автомата. Высота контртела в форме конуса составляет 0.2 мкм. Процесс нагружения имитировался заданием скоростей в вертикальном (ось 2) и горизонтальном (ось Y) направлениях всем автоматам контртела. То есть вместо силы, действующей на индентор по оси 2 в эксперименте, в расчетах задавалась постоянная вертикальная скорость равная -1 м/с до погружения индентора на заданную глубину, после чего вертикальная составляющая нагрузки задавалась равной нулю. В рассматриваемом случае глубина погружения индентора составляла расстояние до промежуточного слоя. Для движения контртела вдоль модельного образца всем его автоматам задавалась постоянная скорость в горизонтальном направлении по оси Y равная -5 м/с. Отметим, что в реальных экспериментах по склерометрии скорость движения индентора относительно материала, как правило, не превышает нескольких сантиметров в секунду.

По результатам численного эксперимента были построены визуальные изображения царапины на поверхности испытанного материала. На рис. 6 показана сетка межавтоматных связей, где отрезками соединены центры автоматов, находящихся в связанном состоянии. При этом темным цветом показаны связи автоматов покрытия, более светлым — связи автоматов переходного слоя, а самым светлым — подложки. Сравнение струкРис. 6. Общий вид модели для расчетов склерометрии в виде сетки межавтоматных связей (а) и вид образца после прохождения индентора (б)

О 10 20 30 40 мкс

Рис. 7. Зависимость коэффициента трения (пунктирная линия) от времени расчета и ее сглаженная кривая (сплошная линия)

туры межавтоматных связей в исходном состоянии (рис. 6, а) и после прохождения индентора (рис. 6, б) позволяет увидеть, что вдоль пути прохождения индентора переходный слой и примыкающая к нему часть подложки деформировались пластически. При этом покрытие растрескивается и разрушается вплоть до переходного слоя, однако полного отслоения материала покрытия не происходит. Качественно полученные изображения соответствуют начальной стадии разрушения покрытия в натурном эксперименте по измерительному царапанию наноструктурных покрытий на титановой подложке, описанному в [8].

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

Полученную в результате проведения численного эксперимента зависимость (рис. 7) сравнивали с данными натурного эксперимента [8]. Средняя величина коэффициента трения, полученная в результате компьютерного моделирования (0.61), соответствует экспериментальному значению в момент разрушения покрытия (0.63, при этом в покрытии появляются большие трещины, но отслоения от подложки не происходит).

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

Механические свойства покрытий сильно зависят от наличия в них повреждений, в том числе нанораз-мерного масштаба. Поэтому актуальным является развитие методов диагностики таких повреждений. Одним из методов неразрушающего контроля качества покрытия, позволяющим оценить дефекты структуры тонких поверхностных слоев в твердом теле, является подход, связанный с использованием силы трения в качестве измеряемого и анализируемого параметра отклика системы, — трибоспектроскопия [16]. Поэтому далее в работе изучены возможности применения трибоспект-рального способа для анализа дефектности поверхностных слоев твердых тел на примере пустот нанораз-мерного масштаба различной формы и глубины их залегания в покрытии.

Методика моделирования в рамках метода МСА процессов, имеющих место при трибоспектроскопии, описана в работах [16, 45, 46]. Их сущность сводится к моделированию движения контртела по поверхности образца. Вследствие искусственной шероховатости, вызванной дискретностью метода моделирования, между контртелом и поверхностью образца возникают силы трения, а от пятна контакта генерируются упругие волны, распространяющиеся по образцу. Вследствие динамического характера такого процесса регистрируемая сила трения скольжения меняется во времени. Преобразование Фурье данной силы для бездефектного образца позволяет выделить некоторые характерные частоты [45]. К ним относятся частота, соответствующая искусственной шероховатости поверхности образца (она связана с размером элементов модели — клеточных автоматов), и собственные частоты моделируемой системы.

Следующим этапом исследования по идентификации наноскопических дефектов в материале покрытия является выявление частот, характеризующих размер дефектов и период их следования [28]. В настоящей работе исследовано влияние глубины расположения нанопор от поверхности, их формы, а также толщины образца на возможность идентификации частот, обусловленных наличием пор.

Моделировалось керамическое покрытие на титановом образце (рис. 8). Толщина покрытия была выбрана равной Н= 60 нм, размер автомата d = 3 нм, длина образца Ь = 250 нм, ширина образца М= 125 нм. В рамках используемого приближения покрытие из керамики ZrO2 считалось однородным и изотропным.

Контртело имело форму усеченного конуса с диаметром основания 60 нм, его движение моделировалось заданием для автоматов верхней поверхности контртела постоянной скорости V = 5 м/с в направлении оси У (рис. 8). При этом нижняя поверхность образца (слой автоматов, моделирующий подложку) была неподвижной, а его боковые поверхности — свободными. При

Рис. 8. Образцы с различными видами пор: общий вид моделируемого образца со сквозными трещинами, расположенными на глубине 30 нм (а), разрез образца со сквозными нанотрещинами (б), разрез образца с несквозными повреждениями, расположенными посередине (в), разрез образца с несквозными повреждениями, расположенными сбоку (г)

движении контртела проводилась регистрация силы сопротивления его движению по поверхности (силы трения Р). Затем строилось преобразование Фурье этой функции [45].

Поврежденность поверхностного слоя моделировалась заданием наноскопических нарушений сплошности покрытия. Анализировались протяженные повреждения — нанопоры. Рассматривалось периодическое расположение нанопор с расстоянием Р = 75 нм различного вида: сквозные повреждения, несквозные расположенные посередине образца и несквозные нанопоры, расположенные сбоку образца (рис. 8). При этом длина всех пор составляла w = 48 нм, высота h = 3 нм, а ширина пор, расположенных сбоку и посередине, равна их дли| 200

Частота, ГГц

Рис. 10. Спектры Фурье силы, действующей на контртело, в различное время расчета

не и составляла 48 нм. В работе изучалась возможность определения расстояния между такими дефектами.

Частота, соответствующая периоду следования нанопор, вычисляется по формуле V|P, где V — скорость движения контртела, а Р—период следования нанопор, и в нашем случае приблизительно равна 48 МГц. Данная частота хорошо видна на спектрах Фурье для составляющей силы сопротивления в направлении оси 2, независимо от глубины расположения и вида повреждений, но имеет тенденцию к уменьшению амплитуды при увеличении глубины их расположения от поверхности (рис. 9).

Частота, соответствующая размеру пор (48 нм), четко идентифицируется на спектрах горизонтальной составляющей силы сопротивления и равна 104 МГц. Следует отметить, что пик на этой частоте определяется достаточно четко на начальной стадии движения инден-тора (до прохождения второй поры, пунктирная линия на рис. 10), после чего амплитуда этого пика значительно уменьшается и на спектре появляется частота, характеризующая период расположения пор (сплошная линия на рис. 10).

Рис. 9. Спектры Фурье силы, действующей на контртело, для образцов с нанопорами, расположенными посередине образца на глубине от поверхности 6 (1), 30 (2) и 54 нм (3)

Рис. 11. Спектры Фурье вертикальной составляющей силы для образцов с различными видами нанопор: сквозные (1), посередине (2) и боковые (3); глубина пор от поверхности 30 нм

ÜJ 005 & 0Л5 & Ö25

Частота, ГГц

Рис. 12. Спектры Фурье вертикальной составляющей силы для образцов с нанопорами, расположенными сбоку образца с толщиной 60 (1), 70 (2) и 90 нм (3)

Величина амплитуды пиков также зависит от вида пор. Наибольшее значение она имеет при сквозных порах, меньшее — при расположенных посередине и минимальное — сбоку (рис. 11).

Кроме того, исследовалось влияние толщины покрытия, значения которой составляли 60, 70 и 90 нм, при условии одинаковой глубины расположения нанопор от исследуемой поверхности. Было установлено, что чем больше толщина покрытия, тем меньше амплитуда пика, соответствующего периоду следования нанопор. Для примера на рис. 12 представлены спектры Фурье для боковых трещин с различной толщиной покрытия.

Выявленный эффект влияния наноповреждений на силу трения скольжения, очевидно, связан с прогибом поверхности при прохождении контртела над трещиной: чем больше площадь трещины, тем больше значение амплитуды пика, отвечающего за период следования или размер поры [28]. В то же время полученная зависимость амплитуды пика, отвечающего за период между повреждениями, от толщины образца, по-видимому, обусловлена отражением волн от закрепленной подложки и затуханием сигнала при большой толщине образца. Следовательно, можно предположить, что для идентификации повреждений с помощью трибоспектро-скопии важен не только прогиб непосредственно в месте нахождения повреждения, но и полный путь прохождения упругих волн.

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

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

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

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

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

Данные исследования выполнены в рамках Программы фундаментальных научных исследований государственных академий наук на 2013-2020 гг., работ по государственному контракту № 16.523.12.3002, скоординированному с проектом Евросоюза ViNaT, Contract No. 295322, FP7-NMP-2011-EU-Russia, NMP.2011.1.4-5.

Литература

Landau L.D., Lifshitz E.M. Theory of Elasticity. - Oxford: Pergamon Press, 1986.

Development of the theory of contact problems in the USSR / Ed. by L.A. Galin. - Moscow: Nauka, 1976. - 493 p.

Goryacheva I.G. Mechanics of Friction Interaction. - Moscow: Nauka, 2001. - 477 p.

Popov V.L. Contact mechanics and friction. Physical Principles and Applications. - Berlin: Springer-Verlag, 2010. - 362 p.

Levashov E.A., Petrzhik M.I., Tyurina M.Ya., Kiryukhantsev-Korneev F.V., Tsygankov P.A., Rogachev A.S. Multilayer nanostructured

heat-generating coatings. Preparation and certification of mechanical and tribological properties // Metallurgist. - 2011. - V. 54. - No. 910. - P. 623-634.

Golovin I.Yu. Nanoindentation and its Possibilities. - Moscow: Mashi-nostroenie, 2009. - 316 p.

Shugurov A.R., Panin A.V., Shesterikov E.V Sclerometric study of galvanic AuNi and AuCo coatings // Tech. Phys. Lett. - 2011. - V 37. -No. 3. - P. 223-225.

Azarenkov N.A., Litovchenko S.V, Beresnev V.M., Chishkala V.A., Gravnova L.V., Ogienko S.I. Peculiarities of fracture of silicide coatings at mechanical dynamic loadings // Fiz. Inzhen. Pover. - 2012. -V. 10. - No. 4. - P. 411-417.

Psakhie S.G., Popov V.L., Shilko E.V, Smolin A.Yu., Dmitriev A.I. Spectral analysis of the behavior and properties of solid surface layers. Nanotribospectroscopy // Phys. Mesomech. - 2009. - V. 12. -No. 5-6. - P. 221-234.

Psakhie S.G., Ostermeyer G.P., Dmitriev A.I., Shilko E.V., Smolin A.Yu., Korostelev S.Yu. Method of movable cellular automata as a new trend of discrete computational mechanics. I. Theoretical description // Phys. Mesomech. - 2000. - V. 3. - No. 2. - P. 5-12.

Popov V.L., Psakhie S.G. Theoretical principles of modeling elasto-plastic media by movable cellular automata method. I. Homogeneous media // Phys. Mesomech. - 2001. - V. 4. - No. 1. - P. 15-25.

Psakhie S.G., Moiseenko D.D., Dmitriev A.I., Shilko E.V., Korostelev S.Yu., Smolin A.Yu., Deryugin E.E., Kulkov S.N. A possible method of computer-aided design of materials with a highly porous matrix structure based on the method of moving cellular automata // Tech. Phys. Lett. - 1998. - Т. 24. - № 2. - С. 154-156.

Smolin A.Yu., Roman N.V., Dobrynin S.A., Psakhie S.G. On rotation in the movable cellular automaton method // Phys. Mesomech. -2009. - V. 12. - No. 3-4. - P. 124-129.

Psakhie S.G., Smolin A.Yu., Stefanov Yu.P., Makarov P. V., Shilko E.V., Chertov M.A., Evtushenko E.P. Simulation of behavior of complex media on the basis of a discrete-continuous approach // Phys. Mesomech. - 2003. - V. 6. - No. 5-6. - P. 47-56.

Levashov E.A., Petrzhik M.I., Kiryukhantsev-Korneev F.V., Shtans-kiiD.V., Prokoshkin S.D., Gunderov D.V, Sheveiko A.N., Korotits-kiiA. V., Valiev R.Z. Structure and mechanical behavior during indentation of biocompatible nanostructured titanium alloys and coatings // Metallurgist. - 2012. - V. 56. - No. 5. - P. 395-407.

Dobrynin S.A., Kolubaev E.A., Smolin A.Yu., Dmitriev A.I., Psakhie S.G. Time-frequency analysis of acoustic signals in the audiofrequency range generated during Hadfield&s steel friction // Tech. Phys. Lett. - 2010. - V. 36. - No. 7. - P. 606-609.

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

Сведения об авторах

Cмолин Алексей Юрьевич, д.ф.-м.н., проф., внс ИФПМ CO РАН, проф. ТГУ, asmolin@ispms.tsc.ru

Еремина Галина Максимовна, инж. ИФПМ CO РАН, асп. ТГУ, anikeeva@ispms.tsc.ru

Cергеев Валерий Викторович, мнс ИФПМ CO РАН, инж. ТПУ, svalera@ispms.tsc.ru

Шилько Евгений Викторович, д.ф.-м.н., доц., внс ИФПМ CO РАН, проф. ТГУ, shilko@ispms.tsc.ru

Псахье Cергей Григорьевич, д.ф.-м.н., чл.-корр. РАН, дир. ИФПМ CO РАН, проф. ТГУ, зав. каф. ТПУ, sp@ispms.tsc.ru