CC BY
100
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2013 Математика и механика № 4(24)
УДК 004.94 : 539.62
А.Ю. Смолин, Г.М. Аникеева, Е.В. Шилько, С.Г. Псахье
МОДЕЛИРОВАНИЕ ДЕФОРМАЦИИ НАНОСТРУКТУРНЫХ ПОКРЫТИЙ НА ТИТАНОВОЙ ПОДЛОЖКЕ ПРИ НАНОИНДЕНТИРОВАНИИ1
Представлены результаты моделирования процесса наноиндентирования упрочняющих покрытий на титановой подложке методом подвижных клеточных автоматов в трехмерной постановке. Изложены особенности метода подвижных клеточных автоматов для описания упруго-пластического поведения материалов. Сравнение результатов моделирования с экспериментальными данными подтверждает адекватность предложенной модели.
Ключевые слова: наноиндентирование, моделирование, метод подвижных клеточных автоматов, наноструктурный титан, упрочняющие биосо-вместимые покрытия.
Эффективным способом повышения функциональных свойств материалов является нанесение на их поверхность специальных покрытий. Например, для имплантатов используются многокомпонентные биоактивные наноструктурные покрытия (МБНП) на основе тугоплавких соединений ТіС (Ті,Та)(С,№) с добавлением специальных элементов (Са, 2г, 8і, О, Р), которые улучшают как трибологические, так и биоактивные свойства поверхности [1—4]. Следует отметить, что эти покрытия имеют наноструктурное состояние и малую толщину. Для изучения механических свойств покрытий и пленок в настоящее время в основном используется метод наноиндентирования [5-9]. Наноиндентирование - это процесс контролируемого внедрения сверхтвердого наконечника определенной формы (ин-дентора) под действием нарастающей нагрузки в плоскую поверхность образца на глубину менее 100 нм, при этом в процессе нагружения постоянно измеряется сила Р, действующая на индентор, и глубина его погружения в материал И [4]. На основе анализа измеряемой Р - И-диаграммы можно получать такие характеристики материала, как модуль упругости, упругое восстановление и нанотвердость. При таком анализе в основном используется методика Оливера - Фарра [6]. Однако, как показали, например, авторы работы [7], корректно определить модуль упругости и нанотвердость тонких покрытий и пленок по данным наноиндентирова-ния с использованием стандартной методики Оливера - Фарра возможно только при условии совпадения этих характеристик у покрытия и подложки. Очевидно, что в большинстве практически важных приложений это условие не выполняется. Одним из способов решения этой проблемы может быть использование компьютерного моделирования, в рамках которого возможно получение достаточно точных зависимостей для индентирования покрытий на любых подложках.
1 Работа выполнена в рамках государственного контракта № 16.523.12.3002 от 13.05.2011 по ФЦП «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007-2012 гг.» и по скоординированному с Евросоюзом проекту УіКаТ, СоПгаС N0. 295322, РР7-КМР-2011-Еи-Яи88Їа, NMP.2011.L4-5, а также проекта ІІІ.23.2.4. Программы фундаментальных исследований СО РАН на 2013-2016 гг.
Для моделирования механического поведения материала на макроскопическом уровне в рамках любого метода необходимо задавать такие параметры, как модули упругости, предел текучести, параметры деформационного упрочнения, предел прочности и т.д. Обычно эти параметры определяются из стандартных механических испытаний на сжатие-растяжение, изгиб, из измерения скорости распространения упругих волн и т.д. В этой связи следует отметить, что уникальные свойства МБНП обусловлены их сложной структурой, которая определяется технологическими параметрами их получения. Объемных материалов с такой структурой не существует. Следовательно, получить свойства материалов покрытий из обычных испытаний не представляется возможным. Существуют методики получения механических характеристик и даже диаграммы нагружения из данных по наноин-дентированию, основанные на привлечении численных расчетов методом конечных элементов [8-12]. Однако эти методики разработаны для однородных материалов, а для случая систем «покрытие - подложка» таких методик пока не существует.
К настоящему времени опубликовано довольно много работ, посвященных численному моделированию процесса наноиндентирования. В зависимости от используемого метода они могут быть разделены на две группы. В работах первой группы используется метод конечных элементов и рассматривается поведение материала на макромасштабе. Авторы изучают особенности распределения напряжений в материале при различных значениях параметров системы, таких, как форма индентора, пластические свойства материала и т.д. [9-11]. Вторая группа работ посвящена применению метода молекулярной динамики (метод частиц) для изучения на микроуровне механизмов зарождения пластической деформации в непосредственной близости от вершины индентора [13-15].
В данной работе на основе метода частиц предложена численная модель, описывающая механическое поведение МБНП на титановой подложке на мезо- и макромасштабном уровнях. Основная цель состоит в разработке модели, которая позволит изучать процессы наноиндентирования и скретч-теста (царапания) такого рода покрытий. Для этого нужен метод, который позволил бы моделировать на различных масштабах как процесс упруго-пластического деформирования, так и разрушение твердых тел. Как известно, лучшими возможностями для моделирования разрушения, в том числе зарождения и роста трещин, фрагментации и перемешивания вещества, обладают методы частиц, берущие свое начало из метода молекулярной динамики. Следует подчеркнуть, что из всех методов частиц только метод подвижных клеточных автоматов (MCA) способен корректно описывать пластическую деформацию консолидированных тел [16-20], поэтому он и был выбран для построения модели. Отметим, что данная работа представляет собой лишь первый шаг, на котором разработанная модель применяется для изучения деформации покрытия и подложки при небольшой глубине проникновения, соответствующей наноиндентированию, и процессы разрушения не учитываются.
При разработке модели наноструктурного биосовместимого покрытия на титановой подложке, а также ее тестирования были использованы результаты натурных экспериментов, опубликованные в работах [1-4, 21, 22]. Так, геометрические характеристики модельной системы «покрытие - переходный слой - подложка» определялись по фотографии поперечного среза реального образца [21]. Исходя из того, что при малых нагрузках (до 250 мН) в экспериментах не наблюдается разрушение материала покрытия [21, 22], механическое поведение материалов и подложки, и покрытия описывалось в рамках упруго-пластического при-
ближения. Полученные по результатам наноиндентирования значения модуля Юнга и предела текучести учитывались при задании параметров модельных материалов.
Упруго-идеальнопластическое тело характеризуется следующими параметрами: плотность р, модуль сдвига О, модуль объемной упругости К и предел текучести оу. На основе данных, приведенных в [21, 22], значения указанных параметров для подложки составили р = 4500 кг/м3, О = 41 ГПа, К = 100 ГПа и оу = 2 ГПа, а для покрытия - р = 4700 кг/м3, О = 76 ГПа, К = 167 ГПа и оу = 15 ГПа.
В экспериментах по наноиндентированию обычно используются алмазные пирамиды Берковича. Модуль упругости алмаза (1141 ГПа) приблизительно в пять раз превышает модуль упругости для МБНП (200-250 ГПа). Поэтому в данных расчетах материал индентора принимался абсолютно жестким (недеформируе-мым).
Описание модели
Метод МСА [16-20] является новым эффективным численным методом, основанным на концепции дискретных элементов (частиц), которая существенно отличается от концепции численных методов классической механики сплошных сред. Следует отметить, что развитие дискретных методов началось с метода молекулярной динамики (МБ), который был разработан для изучения материалов на атомном уровне. Однако возможности атомного описания на пространственных и временных масштабах, которые представляют интерес для инженерных приложений, существенно ограничены. Это побудило развитие МБ-подобных методов для мезо-и макромасштабных уровней (методов частиц), в которых структурные элементы имеют конечный размер (в отличие от атомов, которые являются точечными массами) и взаимодействуют только с ближайшими соседями. Самым известным представителем этой группы методов является метод дискретных элементов (БЕМ) [23, 24]. БЕМ в настоящее время широко используется для изучения механического поведения гранулированных (сыпучих) и слабосвязанных сред, в частности реологических особенностей этих систем, их разрушения и перемешивания [25, 26]. В то же время, до недавнего времени применимость БЕМ для изучения механического поведения консолидированных тел была ограничена главным образом хрупкими пористыми материалами [24-27] в связи с недостаточным развитием математических моделей для вычисления сил взаимодействия дискретных элементов. В частности, большинство моделей дискретных элементов используют парные (двухчастичные) потенциалы (или силы) взаимодействия. Это упрощение приводит к ряду искусственных эффектов в поведении ансамбля частиц, таких, как анизотропия, навязанное упаковкой значение модуля сдвига и коэффициента Пуассона и т. д. [24, 27]. Среди проблем, возникающих вследствие этих искусственных эффектов, важно отметить неадекватность моделирования накопления необратимых деформаций (пластичности) материалов.
Последние исследования показали, что многие проблемы БЕМ, в том числе касающиеся описания консолидированных твердых тел на различных масштабах, могут быть решены с помощью многочастичности сил взаимодействия между элементами. В [20] описан обобщенный подход к построению многочастичных сил взаимодействия для дискретных элементов, который основан на идее, применяемой для записи межатомных потенциалов в методе погруженного атома [28].
Так, выражение для силы, действующей на дискретный элемент г', имеющий N соседей, может быть записано в виде
N
Рг = > . + рп.
]=1
Эта сила представлена в виде суперпозиции парных компонентов Г^Р311, зависящих от пространственного расположения автомата г по отношению к соседу ., и объемнозависящей компоненты Ггп, обусловленной коллективными эффектами окружения.
В рамках метода МСА предполагается, что любой материал состоит из определенного количества элементарных объектов конечного размера (автоматов), которые взаимодействуют друг с другом и могут перемещаться в пространстве, тем самым моделируя реальные процессы деформации. Движение автоматов описывается уравнениями Ньютона - Эйлера:
72 О М.
тг = > Гг!™ + ,
1 Л2 £ г г
- а20г N ы
Jг-------------------= > М..
Л2 р. 4
(1)
где Я , 0г-, тг и Ji - радиус-вектор, вектор вращения, масса и момент инерции г
Г р£
автомата соответственно, - парная сила механического взаимодействия ав-
томатов г и] , Г° - объемнозависящая сила, действующая на автомат г и обусловленная взаимодействием его соседей с другими автоматами. В последнем уравнении М. = дг. (п. х .г) + К., здесь д . - расстояние от центра г-го автомата до
точки его взаимодействия («контакта») с .-м автоматом, п. = (Я. - Я)/. - единичный вектор ориентации пары и г. - расстояние между центрами автоматов (рис. 1), К. - крутящий момент, обусловленный относительным вращением автоматов в паре (см. ниже).
Отметим, что автоматы пары могут представлять собой части различных тел или одного консолидированного тела. Поэтому их взаимодействие не всегда является действительно контактным. По этой причине здесь и далее мы будем слово «контакт» брать в кавычки. Более того, как показано на рис. 1, а, размер автомата характеризуется одним параметром ф, но это не означает, что форма автомата представляет собой шар. Реальная форма автомата определяется областью его «контактов» с соседями. Например, если использовать в качестве начальной ГЦК-упаковку, то автоматы будут иметь форму ромбического додекаэдра, а если кубическую, - то автоматы будут представлять собой кубики.
Для локально изотропных сред объемнозависящая компонента может быть записана через давление Р. в объеме соседнего автомата . следующим образом [16, 18]:
N
Г = -А> РА, п.,
.=1
где 8. - площадь контакта г-го автомата с .-м, а А - некий материальный параметр.
Рис. 1. Параметры пространственного отношения пары подвижных клеточных автоматов
С другой стороны, общая сила, действующая на автомат, может быть представлена в виде суммы нормальной Ж, и касательной (сдвиговой) Ж, компонент:
М, М,
Ж = X (Г - п,) = )[( Г,) - , К + ЬГ’ (Г Ь ] =
1=1
1=1
N
= !Г + ,
,=1
(2)
где - нормальная, а ^ра1ГТ - касательная силы взаимодействия в паре, зави-
сящие соответственно от межавтоматного перекрытия hр (рис. 1, а) и относительного тангенциального смещения 1реЛ1 (рис. 1, б), рассчитанного с учетом вращения обоих автоматов [19, 20]. Отметим, что, несмотря на то, что последнее выражение в уравнении (2) формально соответствует обычной записи силы взаимодействия в методе дискретных элементов [23-27], оно принципиально отличается вследствие многочастичности центрального взаимодействия автоматов.
Используя процедуру осреднения для тензора напряжений в частице, изложенную в [20, 27], выражение для компонент усредненного тензора напряжений в центре автомате принимает вид
1
N
«Р т/ Х1і]Пі],а^у,Р
Уі 1=1
где а и р обозначают оси X, Y, Z лабораторной системы координат, Vt - текущий объем автомата i, п, а - а-компонента единичного вектора п, и Fj р - Р-компо-нента полной силы, действующей в точке «контакта» между автоматами i и j.
Давление P,, или, что то же самое, среднее напряжение ^mean, в объеме автомата может быть вычислено через компоненты тензора напряжений:
—І . —і . —і , ахх +СТ VV +СТ ZZ
P = -нг = yy zz ()
i mean З Знание компонент тензора напряжений позволяет вычислять все его инварианты в объеме автомата, в частности интенсивность напряжений:
CTi"^T2)1 (CTXx -стУу)2 +(стУу -CTZz)2 +(<4 -ctL)2 +б (y)2+(*)2 +(4)2
(5)
"72
Из уравнений (1), (2) и (3) следует, что выражения для вычисления £гр,а1Г,п и рр,т,1 определяют реологические свойства модельной среды.
В дальнейшем параметры взаимодействия подвижных клеточных автоматов для удобства рассматриваются в относительных (удельных) величинах. Так, центральное и тангенциальное взаимодействие автоматов г и р характеризуется соответствующими напряжениями Пр и Тр :
К=рр • (6)
I рт,5,.
Отметим, что в большинстве работ, посвященных описанию и использованию метода МСА, уравнения и формулы метода записаны для двумерного случая. В данной работе используется трехмерная версия метода МСА, поэтому здесь напряжение сдвига тр представляет собой вектор, лежащий в плоскости, нормальной к единичному вектору Пр .
Деформацию автомата г при его центральном взаимодействии с автоматом р можно характеризовать следующей безразмерной величиной (нормальной деформацией):
Чц - й /2
^ = , Л . (7)
В общем случае каждый автомат представляет собой различные материалы и перекрытие в паре некоторым образом перераспределяется между -м и р-м автоматами:
Щ = Ад у + Ад]г = йг /2 + 2 , (8)
где символ Д означает приращение величины за временной шаг Д/ численного интегрирования уравнений движения (1). Правило перераспределения деформации в паре тесно связано с выражением для вычисления силы взаимодействия автоматов. Это выражение для центрального взаимодействия похоже на соотношения Гука для диагональных компонент тензора напряжений:
ДПр = 20(АЪр) + (1 - 2%)Р , (9)
где К - модуль объемной упругости, О - модуль сдвига материала автомата г, Рг -давление автомата г, которое может быть вычислено с помощью формул (3) и (4) на предыдущем шаге по времени или по схеме предиктор - корректор.
Для определения параметра, характеризующего деформацию сдвига в паре автоматов г — р, воспользуемся формулой кинематики для свободного движения пары как абсолютно жесткого (недеформируемого) тела:
ир - иг = Югр Х Ггр , (10)
где Гр = (Яр - Я), и - трансляционная скорость автомата г, Юр - мгновенная
скорость вращения пары как целого. Если мы умножим обе части уравнения (10) векторно слева на Гр и пренебрежем вращением вокруг оси пары (т.е. пусть
Юр • Гр = 0, так как вращение вокруг оси пары не производит деформацию сдвига), то получим следующую формулу:
Гр Х (ир - и ) Пр Х (ир - иг )
Юр =~------2------= -I----1----- . (11)
Г 2 Г
Гу V
Кроме такого вращения пары как целого (определяется разницей в поступательных скоростях автоматов) каждый из автоматов вращается со своей собственной вращательной скоростью юг (рис. 1, в). Разница между этими вращениями
обусловливает деформацию сдвига. Так, приращение деформации сдвига автоматов г и р за шаг Д/ определится относительным тангенциальным смещением в точке «контакта» пары Д1реа1, деленной на расстояние между автоматами
Ау, + Ду= —-------= —----------:-------——---------------—— . (12)
г гу
Выражение для вычисления силы тангенциального взаимодействия подвижных клеточных автоматов похоже на соотношения Гука для недиагональных компонент тензора напряжений:
Ат, = 2С(Ау,.), (13)
и является чистым парным.
Разные скорости вращения каждого из автоматов также приводят к деформации относительного «изгиба» и «кручения» (последнее присутствует только в 3Б) в паре. Очевидно, что сопротивление относительному вращению в паре вызывает крутящий момент, величина которого пропорциональна разнице между поворотами автоматов:
К, =-(Ог + О, )(0, - 0,). (14)
Формулы (1) - (4), (6) - (9), (11) - (14) описывают механическое поведение
линейно-упругого тела в рамках метода МСА. Отметим, что соотношения (8), (9),
(12) и (13) записаны в приращениях, т. е. в гипоупругой форме. В работе [29] показано, что эта модель дает результаты, соответствующие численному решению методом конечных разностей обычных уравнений механики сплошной среды для изотропной линейно-упругой среды. Это обстоятельство позволило использовать метод МСА совместно с численными методами механики сплошных сред в рамках единого дискретно-континуального подхода. В работе [19] показано, что
именно учет вращения автоматов совместно с многочастичностью центрального взаимодействия позволяет МСА адекватно описывать изотропный отклик материала на приложенную нагрузку.
Для описания упругопластического поведения в рамках метода МСА предлагается использовать теорию пластического течения, а именно модель идеальной пластичности с критерием Мизеса. Для этого к методу МСА [20] был адаптирован известный алгоритм Уилкинса [30, 31]. Этот алгоритм состоит в решении упругой задачи на каждом временном шаге и последующем «сбросе» компонент девиатора тензора напряжений £>ар = стар-1/3 сти5ар на поверхности текучести Мизеса в
случае, когда интенсивность напряжений превышает заданную предельную величину (рис. 2):
Этот алгоритм в применении к автомату г может быть записан в следующих обозначениях:
ся в результате решения упругой задачи на текущем временном шаге,
радиуса круга текучести Мизеса для автомата г.
Удельные величины нормальных и касательных сил взаимодействия корректируются с использованием текущего значения коэффициента М по аналогии с алгоритмом Уилкинса для компонент осредненного тензора напряжений (14):
где Пу и тУ представляют собой скорректированные значения удельных сил (на-
Рис. 2. Схема работы алгоритма Уилкинса
mean
)м, +
mean
(14)
где а, в = X, У, 2 и а^р, - скорректированные компоненты осредненного
тензора напряжений, ^р - компоненты тензора напряжений, которые получают-
Мі = ст^/ст^ - текущее значение коэффициента «сброса», ст^ - текущее значение
пряжений). Легко показать, что подстановка этих соотношений в выражения (2) и (3) автоматически обеспечивает корректировку компонент осредненного тензора напряжений в автомате к кругу текучести [20].
Таким образом, реологические свойства материала автомата г определяются
заданием единой кривой упрочнения =ф(ё^) (здесь - интенсивность
осредненного тензора деформаций, компоненты которого могут быть вычислены аналогично ст^р [20, 27]), эта зависимость также называется функцией отклика автомата.
Для численного интегрирования уравнений движения (1) можно использовать схему Верле в скоростной форме, модифицированную введением предиктора для оценки ст^р на текущем шаге времени.
Пару элементов можно рассматривать как виртуальный бистабильный автомат (у него существуют два состояния: связанная и несвязанная пара), что позволяет явно моделировать процессы разрушения в методе МСА. Заданием правил перехода пары из состояния связанной в состояние несвязанной формулируется критерий разрушения моделируемого материала, который, вообще говоря, определяется физическими механизмами деформации материала. Важным преимуществом описанного выше формализма является то, что в нем возможно непосредственное применение классических критериев (Губера - Мизеса - Генки, Друккера - Пра-гера, Кулона - Мора, Подгорски и т.д.), которые записаны в тензорной форме. Заметим, что переключение пары автоматов в несвязанное состояние может привести к изменению сил, действующих на элементы, в частности, они не будут сопротивляться удалению друг от друга.
Таким образом, метод МСА позволяет моделировать механическое поведение твердого тела на различных масштабах, в том числе пластическое и вязкоупругое деформирование, разрушение, фрагментацию и дальнейшее взаимодействие фрагментов как сыпучей (гранулированной) среды.
Тестовые расчеты
Моделируемая система «индентор - покрытие - подложка» представляет собой гетерогенную структуру, составленную из материалов с существенно различными механическими свойствами: сверхпрочный алмаз (индентор), прочное покрытие и пластичная подложка из титана. Поэтому вначале целесообразно провести тестовые расчеты по одноосному сжатию каждого из материалов, из которых состоит система. Эти расчеты позволяют определить параметры функции отклика подвижных клеточных автоматов, имитирующих свойства исследуемых материалов, а также представительный объем соответствующих модельных образцов.
Для проведения тестовых расчетов генерировались трехмерные кубические образцы с ГЦК-упаковкой автоматов, которые подвергались одноосному сжатию. Исходя из среднего размера зерна наноструктурного титана и МНБП размер автомата выбирался равным 10 нм. Нагружение имитировалось заданием единой скорости движения в вертикальном направлении всем автоматам верхнего слоя автоматов и фиксацией в вертикальном направлении автоматов нижнего слоя. При этом всем автоматам нижнего и верхнего слоя разрешалось движение в горизонтальной плоскости для того, чтобы напряженное состояние в образце было более однородным.
На основе расчетных данных был проведен анализ сходимости упругих и прочностных свойств образцов с увеличением их размеров. Для этого рассчитывались модули упругости на сжатие (наклон первого линейного участка диаграммы нагружения) и пределы текучести (максимальное значение диаграммы) для каждого из образцов. Относительное отклонение расчетных модулей упругости Ег всех образцов от модуля Юнга Е в зависимости от размера основания кубических образцов Иг приведено на рис. 3, а. Относительное отклонение расчетных значений пределов текучести сг от заданного значения су приведено на рис. 3, б. На основании представленных данных можно заключить, что представительными можно считать образцы с размером основания 200 нм (20 автоматов).
(EІ - Е)/Е, % (стг - сту)1сту, %
Рис. 3. Сходимость модуля упругости (а) и предела текучести (б) тестовых образцов титана (1) и МБНП (2)
Результаты расчетов и их обсуждение
После получения параметров функции отклика автоматов для материалов моделируемой системы была построена полная модель процесса индентирования. Эта модель состояла из образца исследуемого материала и индентора. Геометрически индентор Берковича представляет собой трехгранную пирамиду, в основании которой лежит правильный треугольник. Реальные размеры инденторов производителями оборудования не сообщаются, приводятся углы наклона граней относительно оси индентора и радиус скругления вершины. Угол между гранями пирамиды Берковича и ее осью составляет 65,3°. Скругление вершины в расчетах учитывалось усечением пирамиды так, чтобы сторона получающейся треугольной грани на вершине равнялась радиусу скругления (50 нм). Высота индентора в модели определялась максимальной глубиной вдавливания, известной из экспериментальных данных, и в различных расчетах варьировалась от 90 до 950 нм.
Моделировалось два вида образцов, каждый из которых имел форму параллелепипеда, их высота варьировалась от 700 до 2500 нм. Первый вид образцов состоял полностью из наноструктурного титана, а второй - из нескольких слоев, нижний из которых представлял собой подложку из титана, а верхний - МБНП-покрытие. Толщина такого покрытия составляет 350-2000 нм [1-4, 21, 22]. Между покрытием и титановой подложкой находится промежуточный слой толщиной
порядка 50 нм. В модели свойства промежуточного слоя задавались средним между значениями для покрытия и подложки.
Общий вид построенной модели представлен на рис. 4, a в виде ГЦК-упаковки сферических частиц. Процесс нагружения имитировался путем задания одинаковой скорости движения всем автоматам индентора. Для исключения динамических эффектов эта скорость плавно менялась от 0 до -1 м/с, после чего оставалась постоянной до заданного времени погружения. Затем скорость движения индентора плавно менялась на противоположную (+1 м/с), после чего индентор с постоянной скоростью двигался вверх. Для того чтобы предотвратить смещение нагружаемого материала как целого, нижний слой автоматов подложки задавался жестко закрепленным.
в г
Рис. 4. Общий вид модели (а), распределение скоростей автоматов (б, абсолютные значения в м/с), а также интенсивности напряжений (Па) по сечению образца при его упругой (в) и пластической (г) деформации
Анализ поля скоростей в начале индентирования, когда материал деформируется упруго, показывает, что смещения в материале локализуются под инденто-ром и направлены преимущественно вертикально, а не в боковых направлениях (рис. 4, б). Упругие напряжения также локализованы под индентором и распространяются в глубь материала (рис. 4, в). Полученные данные хорошо согласуются с экспериментальными, которые показывают, что результаты индентирования тонкого покрытия в значительной степени определяются свойствами подложки. Рис. 4, г показывает поле интенсивности напряжений при максимальной глубине проникновения, когда покрытие и подложка деформируются уже пластически. Видно, что максимальные напряжения находятся в непосредственной близости от индентора, однако из-за меньшего предела текучести у подложки она также деформируется пластически. Поля напряжений, представленные на рис. 4, г, находятся в хорошем качественном согласии с данными из литературы [33, 34].
Главной характеристикой процесса наноиндентирования является зависимость силы P, действующей на индентор, от глубины вдавливания h, так называемая P - h-диаграмма. В случае малой глубины проникновения эта кривая в основном определяется свойствам покрытия. Чем больше глубина проникновения, тем большее влияние на кривую оказывают свойства подложки. Поэтому вначале проводилось моделирование погружения индентора в образец из титана, затем в образец с покрытием на небольшую глубину и в заключение - индентирование на значительные глубины проникновения в образец с покрытием.
На рис. 5 представлены расчетная и экспериментальная [21] P - h-диаграммы индентирования титанового образца на глубину 650 нм. Их анализ позволяет сделать вывод об адекватности построенной модели для описания механического поведения материала подложки.
Для определения адекватности параметров модельного покрытия был сделан расчет с глубиной внедрения в образец с покрытием равной 60 нм [32]. Результаты этого расчета сравнивались с экспериментальными данными при максимальной нагрузке на индентор 2 мН [21, 22]. Сравнение показало хорошее качественное и количественное соответствие численного расчета эксперименту, что позволило сделать вывод об адекватности предложенной модели механического поведения наноструктурных биосовместимых покрытий.
Для изучения влияния подложки на P - h-диаграмму моделировалось погружение на большую глубину. В связи с этим использовались образцы с общей высотой 2,5 и 5 мкм и размером автоматов 70 и 140 нм соответственно. Толщина покрытия составляла 1000 нм. P - h-диаграммы, полученные в результате моделирования и эксперимента [21, 22], показаны на рис. 6. Видно очень хорошее качественное и количественное соответствие между численными и экспериментальными данными.
Рис. 5. Экспериментальная [21] (1) и рас- Рис. 6. Экспериментальные [21, 22] (1, 3) и
четная (2) диаграммы индентирования ти- расчетные (2, 4) диаграммы индентирова-
танового образца ния титанового образца с наноструктурным
покрытием
Другой важной характеристикой процесса наноиндентирования является размер и форма отпечатка индентора на поверхности материала. На рис. 7 показано поперечное сечение модельного образца при проникновении индентора на глубину 950 нм (соответствующий эксперимент имеет нагрузку 100 мН). Видно, что на поверхности по контуру отпечатка индентора не образуется навалов (pile-up), что согласуется с данными, полученными в [21, 22] для больших нагрузок.
Рис. 7. Поперечное сечение модельного образца после индентирования
Заключение
Предложена модель для описания механического поведения биоактивных наноструктурных покрытий на титановой подложке методом подвижных клеточных автоматов в трехмерной постановке. На основе тестовых расчетов получены значения параметров модели. Проведено моделирование внедрения абсолютно жесткого индентора Берковича на различные глубины в образцы из титана и наноструктурного покрытия на титановой подложке. Согласие результатов расчетов с имеющимися экспериментальными данными свидетельствуют о том, что предложенная модель хорошо описывает особенности поведения таких систем. Дальнейшие развитие модели может включать в себя использование более сложных определяющих уравнений для подложки и покрытия, явный учет разрушения и внутренней структуры материалов.
ЛИТЕРАТУРА
1. Shtansky D.V., Kiryukhantsev-Korneev Ph.V., Bashkova I.A., et al. Multicomponent nanostructured films for various tribological applications // Int. J. Refractory Metals & Hard Materials. 2010. 28. P. 32-39.
2. Shtansky D.V., Gloushankova N.A., Bashkova I.A., et al. Multifunctional biocompatible nanostructured coatings for load-bearing implants // Surface and Coatings Technology. 2006. 201. P. 4111-4118.
3. Shtansky D.V., Levashov E.A., Glushankova N.A., et al. Structure and properties of CaO- and ZrO2-doped TiCxNy coatings for biomedical applications // Surface and Coatings Technology. 2004. 182. P. 101-111.
4. Левашов Е.А., Петржик М.И., Тюрина М.Я. и др. Многослойные наноструктурные тепловыделяющие покрытия. Получение и аттестация механических и трибологических свойств // Металлург. 2010. № 9. С. 66-74.
5. Головин И. Ю. Наноиндентирование и его возможности. М.: Машиностроение. 2009. 316 с.
6. Oliver W.C., Pharr GM. An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments // J. Materials Research. 1992. No. 7. P. 1564-1583.
7. Шугуров А.Р., Панин А.В., Оскомов К.В. Особенности определения механических характеристик тонких пленок методом наноиндентирования // ФТТ. 2008. Т. 5. Вып. 6. С. 1007-1012
8. Venkatesh T.A., Van Vliet K.J., Giannakopoulos A.E., Suresh S. Determination of elasto-plastic properties by instrumented sharp indentation: guidelines for property extraction // Scripta Materialia. 2000. V. 42. No. 9. P. 833-839.
9. Dao M, Chollacoop N., Van Vliet K.J, et al. Computational modeling of the forward and reverse problems in instrumented sharp indentation // Acta Materialia. 2001. 49. P. 3899-3918.
10. Bucaille J.L., Stauss S., Felder E., Michler J. Determination of plastic properties of metals by instrumented indentation using different sharp indenters // Acta Materialia. 2003. V. 51. P. 1663-1678.
11. Ogasawara N., Chiba N., Chen X. Measuring the plastic properties of bulk materials by single indentation test // Scripta Materialia. 2006. V. 54. P. 65-70.
12. Sreeranganathan A., Gokhale A., Tamirisakandala S. Determination of local constitutive properties of titanium alloy matrix in boron-modified titanium alloys using spherical indentation // Scripta Materialia. 2008. V. 58. No. 2. P. 114-117.
13. Zimmerman J.A., Kelchner C.L., Klein P.A. et al. Surface step effects on nanoindentation // Physical Review Letters. 2001. V. 87. P. 165507-165511.
14. Saraev D., Miller R.E. Atomic-scale simulations of nanoindentation-induced plasticity in copper crystals with nanometer-sized nickel coatings // Acta Materialia. 2006. V. 54. P. 33-45.
15. Mei J., Li J., Ni Y., Wang H. Multiscale simulation of indentation, retraction and fracture processes of nanocontact // Nanoscale Research Letters. 2010. V. 5. P. 692-700.
16. Псахье С.Г., Остермайер Г.П., Дмитриев А.И. и др. Метод подвижных клеточных автоматов как новое направление дискретной вычислительной механики. I. Теоретическое описание // Физическая мезомеханика. 2000. Т. 3. № 2. С. 5-13.
17. Попов В.Л., Псахье С.Г. Теоретические основы моделирования упругопластических сред методом подвижных клеточных автоматов. I. Однородные среды // Физическая мезомеханика. 2001. Т. 4. № 1.С. 15-25.
18. Psakhie S.G., Horie Y., Ostermeyer G.-P., et al. Movable cellular automata method for simulating materials with mesostructure // Theoretical and Applied Fracture Mechanics.
2001. No. 37. P. 311-334.
19. Смолин А.Ю., Роман Н.В., Добрынин С.А., Псахье С.Г. О вращательном движении в методе подвижных клеточных автоматов // Физическая мезомеханика. 2009. Т. 12. № 2. С. 17-22.
20. Psakhie S.G., Horie Y., Shilko E.V., et al. Development of discrete element approach to modeling heterogeneous elastic-plastic materials and media // Int. J. Terraspace Science and Engineering. 2011. V. 3. No. 1. P. 93-125.
21. Левашов Е.А., Петржик М.И., Кирюханцев-Корнеев Ф.В. и др. Структура и механическое поведение при индентировании биосовместимых наноструктурированных титановых сплавов и покрытий // Металлург. 2012. № 5. С.79-89.
22. Levashov E.A., Petrzhik M.I., Shtansky D.V., et al. Nanostructured titanium alloys and multicomponent bioactive films: Mechanical behavior at indentation // Materials Science and Engineering: A. 2013. V. 570. P. 51-62.
23. Cundall P.A. and Strack O.D.L. A discrete numerical model for granular assemblies // Geotechnique. 1979. V. 29. No. 1. P. 47-65.
24. Jing L. and Stephansson O. Fundamentals of Discrete Element Method for Rock Engineering: Theory and Applications. Oxford: Elsevier, 2007. 562 p.
25. Sibille L., NicotF., Donze F.V., andDarve F. Material instability in granular assemblies from fundamentally different models // Int. J. Numerical and Analytical Methods in Geomechanics 2007. V. 31. No. 3. P. 457-481.
26. Martin C.L. and Bouvard D. Study of the cold compaction of composite powders by the discrete element method // Acta Materialia. 2003. V. 51. No. 2. P. 373-386.
27. Potyondy D.O. and Cundall P.A. A bonded-particle model for rock // Int. J. Rock Mechanics and Mining Sciences. 2004. V. 41. No. 8. P. 1329-1364.
28. Daw M.S., Foiles S.M., andBaskes M.I. The embedded-atom method: A review of theory and applications // Materials Science Reports. 1993. V. 9. No. 7-8. P. 251-310.
29. Псахье С.Г., Смолин А.Ю., Стефанов Ю.П. и др. Моделирование поведения сложных сред на основе комбинированного дискретно-континуального подхода // Физическая мезомеханика. 2003. Т. 6. № 6. С. 11-21.
30. Wilkins M.L. Computer Simulation of Dynamic Phenomena. Berlin: Springer-Verlag, 1999. 246 p.
31. Уилкинс М.Л. Расчет упругопластических течений / Вычислительные методы в гидродинамике. М.: Мир, 1967. С. 212-263,
32. Psakhie S.G., Smolin A. Yu., Shilko E.V., et al. Modeling nano indentation of TiCCaPON coating on Ti substrate using movable cellular automaton method // Computational Materials Science. 2013. (в печати, http://dx.doi.org/10.1016/j.commatsci.2013.03.006)
33. Muliana A., Steward R., Haj-ali R.M., and Saxena A. Artificial neural network and finite element modeling of nanoindentation tests // Metallurgical and Materials Transactions A.
2002. 33A. P. 1939-1948.
34. Feng Z.-Q., Zei M., and Joli P. An elasto-plastic contact model applied to nanoindentation // Computational Materials Science. 2007. V. 38. P. 807-813.
Статья поступила 13.03.2013 г.
Smolin A.Yu., Anikeeva G.M., Shilko E.V., Psakhie S.G. MODELING DEFORMATION OF NANOSTRUCTURED COATINGS ON A TITANIUM SUBSTRATE UNDER NANOINDENTATION. Results of modeling nanoindentation of hardened coating on a titanium substrate by movable cellular automaton method in the 3D formulation are presented. The peculiarities of the method for describing elastic--plastic behavior of the materials are described. Comparing of the modeling results with the experimental data confirms the proposed model validation.
Keywords: nanoindentation, modeling, movable cellular automaton method, nanostructured titanium, hardening biocompatible coatings.
SMOLIN Alexey Yurievich (Tomsk State University)
E-mail: asmolin@ispms.tsc.ru
ANIKEEVA Galina Maximovna (Tomsk State University)
E-mail: anikeeva@ispms.tsc.ru
SHIL’KO Evgenii Viktorovich (Tomsk State University)
E-mail: shilko@ispms.tsc.ru
PSAKH’E Sergey Grigorievich (Tomsk State University)
E-mail: sp@ispms.tsc.ru
