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

УДК 539.3, 539.32, 53.096

Исследование влияния температуры и размеров наноструктур на процесс разрушения под внешними механическими воздействиями

И.Ф. Головнев, Е.И. Головнева, А.В. Уткин

Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, 630090, Россия

В работе исследовано влияние температуры и размеров нанообъектов на механические характеристики и разрушение. Использован метод молекулярной динамики. Показано, что наименьший исследуемый поперечный размер nyxnz = 5x5 кристаллических ячеек является граничным для наноструктур с идеальной кристаллической решеткой и при больших размерах эти характеристики стремятся к аналогичным параметрам макрообъектов. Наблюдается линейное уменьшение модуля Юнга с ростом температуры для всех рассматриваемых систем. Значение приложенного напряжения, при котором происходит разрушение (критическая нагрузка), уменьшается с ростом температуры для всех изученных систем по линейному закону.

Ключевые слова: наноструктура, влияние размера, влияние температуры, разрушение, молекулярно-динамическое моделирование

A study into the effect of temperature and size of nanostructures on fracture

under external mechanical loads

I.F. Golovnev, E.I. Golovneva, and A.V. Utkin

Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia

The paper studies the effect of temperature and size of nano-objects on mechanical characteristics and fracture using molecular dynamics simulation. It is shown that the smallest studied transverse size nyxnz = 5x5 of lattice cells is boundary for nanostructures with a perfect crystal lattice, and at large sizes these characteristics approach the same parameters of macroobjects. Young's modulus and the value of applied stress at fracture (critical load) decrease linearly with temperature rise for all the systems considered.

Keywords: nanostructure, size effect, temperature effect, fracture, molecular dynamics simulation

1. Введение

Фундаментальная задача исследования влияния температуры на механические характеристики твердых тел поднимается уже давно. В монографии [1] подробно рассмотрены вопросы определения долговечности под нагрузкой и отдельная глава посвящена исследованию зависимости долговечности твердых тел под нагрузкой от температуры. В данной монографии, как и в других работах тех лет, поднятый вопрос исследовался в рамках континуального рассмотрения. Атомная природа разрушения вещества изучалась уже в те годы, но в силу ограниченности вычислительных возможностей численные расчеты не проводились. По мере развития вычислительной техники появилась возможность под-

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

Исследованию влияния температуры на механические характеристики с помощью метода молекулярной динамики в мировой науке уделяется достаточно большое внимание. Например, в работе [2] исследуется влияние температуры на предел сдвиговой прочности А1 и Си. В работах [3, 4] исследуется влияние температуры на разрушение углеродных нанотрубок и графена соответственно. Авторы статьи [5] исследовали влияние тепловых флуктуаций на критическое напряжение и разрушение микроструктуры при постоянной внешней нагрузке твердых тел с парным межатомным потенциалом Леннарда-Джонса, который не описывает ме-

© Головнев И.Ф., Головнева Е.И., Уткин А.В., 2017

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

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

2. Приготовление начального состояния системы

Для исследования процесса разрушения строился идеальный кристалл меди в форме прямоугольного параллелепипеда с числом кристаллических ячеек пх = = 50, пу = п2 = 5 (пу = п2 = 7) вдоль соответствующих осей. Для моделирования взаимодействия атомов выбран потенциал, построенный на основе метода внедренного атома [9].

В работе использовалась широко известная скоростная модификация Верле второго порядка точности с шагом по времени 10-16 с. В случае изолированной системы погрешность по энергии не превышала 10-5 % на интервале времени 50 пс.

Первоначально атомы находились в узлах идеальной кристаллической ГЦК-решетки. Затем с помощью метода искусственной вязкости [10] определялись координаты и импульсы системы в состоянии минимума потенциальной энергии. Найденные значения координат и импульсов использовались в качестве начальных данных для разогрева системы до нужной температуры.

3. Моделирование разогрева системы и закрепления структуры в зажимах

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

чайными импульсами (аналог поглощения теплового излучения). Атомы системы через временной интервал

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

тта =^п. (1)

С помощью датчика случайных чисел определяли компоненты приращения импульса атомов по схеме

Др§г- =ДРо(0.5 -1), 1 = х, у, г, (2)

где / — случайные числа в интервале (0, 1). Поскольку число атомов наноструктуры N конечно, то для того чтобы полный импульс системы не изменялся, использовался прием, позволяющий перейти к приращениям компонент импульса Др^, сумма которых по всей системе равна нулю:

N

р§=£Др^, (3)

,=1

ДРЦ =Д?Ц - . (4)

Кроме того, для того чтобы система не приобретала дополнительный момент импульса Мсь, после пересчета всех импульсов со случайными добавками рассчитывали полный момент импульса

N

Мсь = Е[г,р,- ] (5)

i=l

и находили (с учетом, что система — кластер и начало координат совпадает с центром масс) компоненты угловой скорости ^, которые появились за счет конечного и не всегда очень большого числа атомов системы:

мсь = . (6)

После этого пересчитывали скорости всех атомов по формуле

V,- = у,- -[«Г] (7)

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

Приращение полной кинетической энергии системы за временной интервал тта можно представить в виде

Щ = Е

,=1

N (р,- .Др,-)2 N РГ_ N р,- -Др,- N Др,

-Ер-= Е

-Е-

(8)

2т г=12т г=1 т г=1 2т Поскольку импульсы атомов и случайные приращения импульсов не коррелированны, то первое слагаемое в правой части (8) равно нулю и приращение кинетической энергии имеет вид

Щ =Е

Ар^

2т

(9)

Если учесть, что приращение пропорционально амплитуде Др0 (т.е. Др 1 = Др0£,, где £, — случайный вектор с длиной меньшей единицы), то производная кинетической энергии по времени может быть записана в виде

= Др1 ЕЕ

Дt N^1 ,=12т

Таким образом, скорость разогрева определяется двумя параметрами: амплитудой случайного импульса Др0 и числом шагов по времени, за которые происходит разогрев системы.

Далее будет более удобно применять систему единиц, используемую в программе. Единица длины Ь — 10-8 см, единица времени Т — 10-13 с, единица массы т — 10-24 г. Следовательно, единица скорости:

V = (Ь • 10-8 )/(Т • 10-13) см/с = V • 105 см/с = V км/с. Единица силы:

^=тф 2 =т •10-27 • Ь210-10 • -К^-М=т?-ю- н.

Т2 •Ю-26 с2 Отсюда можно получить единицу импульса:

р = р •Ю-24Н• с.

В программе использовался шаг по времени 10-16 с. Следовательно, в программных единицах шаг по времени: т = 0.001.

Кроме того, используется единица энергии:

^ т• Ь2 т•Ю-24 • Ь2 •Ю-16 г• см2

Е =_=___=

Т2 Т2 •Ю-26 с2

= Е •Ю-14 эрг = Е •Ю-21 Дж.

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

1. В первом блоке происходит разогрев системы на величину ДТС по методу, описанному выше. Численные исследования показали, что наиболее оптимальным является выбор ДТС от 5 до 25 К.

2. Во втором блоке система эволюционирует свободно и релаксирует к термодинамическому равновесию в течение времени тЛге1.

3. В третьем блоке система также свободна и проводится усреднение по тепловым флуктуациям в течение времени тЛау. После вывода информации расчет возвращается к первому блоку.

В работе система «разогревалась» со следующими значениями параметров: амплитуда случайного импуль-

са Др0 = 4, число шагов по времени между стохастическим прогревом Лп = 20, время частичной релаксации (число шагов по времени) Лге1 = 500, время усреднения Лау = 500.

В работе использовалась следующая модель зажимов для моделирования одноосного растяжения с постоянным напряжением.

Атомы неподвижной правой грани помещены во внешний потенциал вида

V(г,) = к/4(х, - х°)4 + к/4 (у, - у0)4 + к/4(2, - 2г0)4.

Атомы левой грани помещены во внешний потенциал вида

V (г,) = к/ 4( у - у0)4 + к/ 4( 2,. - 2 )4,

, 0 0 0Ч

где (х,, у,, 21) — координаты г-го атома после охлаждения, к = 5000.

На каждый атом левой грани действует сила вдоль оси X:

Л = ст0

где ст0 — приложенное растягивающее напряжение; 5 — площадь торцевой грани (левой), к которой приложено напряжение; щ — количество атомов в крайней левой атомной плоскости.

На всех этапах расчетов (разогрев, релаксация, квазистатическое растяжение с постоянным напряжением) проводился мезоанализ для нахождения распределения необходимых характеристик по кристаллу. Подробно описание методики проведения мезоанализа приведено в работах авторов [18-20].

4. Моделирование квазистатического нагружения

Моделирование квазистатического нагружения в данной работе базировалось на использовании медленного роста внешней нагрузки до заданного значения, в сочетании с искусственной вязкостью [10] на этапе такого роста. В результате было получено стационарное состояние системы под воздействием внешней статической нагрузки. Однако если система в конечном состоянии должна иметь некоторую температуру, то использование вязкости при этом исключается.

В результате сравнительного анализа ряда вариантов была выбрана следующая модель расчета построения ст-е-диаграмм.

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

2. С помощью метода искусственной вязкости система приводится в минимум потенциальной энергии с учетом наличия свободной поверхности и с температурой близкой к 0 К.

3. Производится разогрев структуры до нужной температуры с использованием координат и импульсов атомов после охлаждения.

0.016"

0.008-

0

0 40 80 120 и пс

0.016-

0.008-,-,-

О 40 80 120 и пс

Рис. 1. Зависимость относительного удлинения от времени. Напряжение 1.0 (а, б), 12.2 ГПа (в); температура 0 (а, в), 100 К (б)

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

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

6. К левому подвижному зажиму прикладывается напряжение, линейно нарастающее в течение 50 пс до заданного значения а0.

7. В течение 100 пс система эволюционирует при заданном напряжении и при этом строится зависимость одноосной деформации всего образца 8^) от времени.

8. Если на этапе релаксации напряжения не происходит разрушения системы, то проводится усреднение 8 ^) по времени.

Все расчеты были проведены для медной наноструктуры размером пх х пу х пг = 50x5x5 кристаллических ячеек вдоль соответствующих осей координат.

Е, ГПа 90-г-

50-1-1-1-

0 200 400 Г, К

Рис. 3. Зависимость модуля Юнга от температуры. Размер 50x5x5 кристаллических ячеек

На рис. 1, а приведена зависимость относительного удлинения от времени для напряжения 1 ГПа и температуры 0 К, а на рис. 1, б — для температуры 100 К.

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

Расчеты предельного значения напряжения при квазистатическом растяжении проводились для внешних напряжений а0 с интервалом 0.1 ГПа, что позволило определить с этой точностью предельную нагрузку для данной структуры и заданной температуры.

5. Исследование влияния размера наноструктуры на механические характеристики в исследуемом интервале температур

Размер нанокристалла, для которого получены все вышеприведенные результаты, составлял пх х пу х п2 = = 50x5x5 кристаллических ячеек вдоль соответствующих осей координат. Были проведены расчеты квазистатического растяжения вдоль оси X этой базовой системы и нанокристалла, поперечный размер которого составляет пу x пг = 7x7, его площадь поперечного сечения примерно в 2 раза больше этой величины для базовой системы.

На рис. 2 приведена а-8-диаграмма для исходного кристалла с сечением 5x5 кристаллических ячеек для трех различных температур. Как видно, зависимость а от 8 линейна для всего интервала исследуемых темпе-

Рис. 2. Зависимость напряжения от величины одноосной деформации. Температура структуры: 0 (1), 300 (2), 500 К (3). Размер 50x5x5 кристаллических ячеек

Рис. 4. Зависимость критического напряжения от температуры. Размер 50x5x5 кристаллических ячеек

Рис. 5. Зависимость напряжения от величины одноосной деформации. Температура структуры: 0 (1), 300 (2), 500 К (3). Размер 50x7x7 кристаллических ячеек

а, ГПа 8

0

а, ГПа 4 2 0

Рис. 6. ст-е-диаграммы для идеального кристалла с поперечным размером 5x5 (сплошная линия) и 7x7 кристаллических ячеек (пунктирная линия) при температуре 0 (а), 500 (б)

ратур, а угол наклона уменьшается с увеличением температуры. Это позволило построить зависимость модуля Юнга от температуры (рис. 3), который линейно зависит от температуры: Е = -0.05021608392Т + + 82.90025641.

Подробный анализ полученных результатов позволил определить зависимость напряжения, при котором наступает разрушение (критическое напряжение), от температуры (рис. 4). Эта характеристика также линейно уменьшается с температурой, и аппроксимация, построенная с помощью метода наименьших квадратов, дает аналитическую зависимость: стсг = -0.01Т + 11.58.

Для случая увеличения площади поперечного сечения примерно в два раза (параметры системы 50x7x7 кристаллических ячеек вдоль осей координат X, У, Т) ст-е-диаграмма приведена на рис. 5. Видно, что зависимости линейны. Влияние поперечного размера приведено на ст-е-диаграммах (рис. 6) для температур 0 и 500 К. Аналогичное совпадение ст-е-диаграмм для наноструктур разного поперечного размера наблюдается для всех температур в исследуемом интервале. Очевидно, то же самое можно ожидать и для модуля Юнга (рис. 7). Кроме того, зависимости критического напряжения от температуры также совпали с большой точностью для температур выше 100 К для двух структур с разными поперечными сечениями (рис. 8).

Е, ГПа

Рис. 7. Зависимость модуля Юнга от температуры. Поперечный размер 5x5 (сплошная линия) и 7x7 кристаллических ячеек (пунктирная линия)

1284-1-,-,-

0 200 400 Г, К

Рис. 8. Зависимость критического напряжения от температуры для структуры с идеальной кристаллической решеткой. Поперечный размер 5x5 (сплошная линия) и 7x7 кристаллических ячеек (пунктирная линия)

Увеличение площади поперечного размера наноструктуры в два раза не повлияло на рассчитываемые квазистатические механические характеристики.

6. Выводы

Сравнение полученных механических характеристик бездефектных систем разного поперечного размера показало следующее.

Увеличение площади поперечного размера наноструктуры в два раза не повлияло на рассчитываемые квазистатические механические характеристики. Следовательно, наименьший исследуемый поперечный размер пу x п2 = 5x5 кристаллических ячеек является граничным для наноструктур с идеальной кристаллической решеткой и при больших размерах эти характеристики стремятся к аналогичным параметрам макрообъектов. При Т = 0 К модуль Юнга бездефектных кристаллов совпадает со справочными экспериментальными значениями для литой меди. Наблюдается линейное уменьшение модуля Юнга с ростом температуры для всех рассматриваемых систем.

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

Работа выполнена при поддержке грантов РФФИ №№ 14-01-00465, 17-01-00068.

Литература

1. Регелъ В.Р., Слуцкер А.И., Томашевский Э.Е. Кинетическая природа прочности твердых тел. - М.: Наука, 1974. - 240 с.

2. Iskandarov A.M., Dmitriev S.V., Umeno Y. Temperature effect on ideal

shear strength of Al and Cu // Phys. Rev. B. Condens. Matter. - 2011. — V. 84. — P. 224118(7). — doi 10.1103/PhysRevB.84.224118.

3. Tang Ch., Guo W., Chen Ch. Structural and mechanical properties of partially unzipped carbon nanotubes // Phys. Rev. B. — 2011. — V. 83. — P. 075410(6). — doi 10.1103/PhysRevB.83.075410.

4. Zhao H., Aluru N.R. Temperature and strain-rate dependent fracture strength of graphene // J. Appl. Phys. — 2010. — V. 108. — P. 064321. — doi 10.1063/1.3488620.

5. Yamamoto A., Kun F., Yukawa S. Microstructure of damage in thermally activated fracture of Lennard-Jones systems // Phys. Rev. E. — 2011. — V. 83. — P. 066108. — doi 10.1103/PhysRevE.83.066108.

6. Subramanian K.R.S. Sankaranarayanan, Venkat R. Bhethanabotla, Babu Josrph. Molecular dynamics simulation of temperature and strain rate effects on elastic properties of bimetallic Pd-Pt nanowires // Phys. Rev. B. — 2007. — V 76. — P. 134117 (13). — doi 10.1103/PhysRevB. 76.134117.

7. Liu Y., Wang F., Zhao J., Jiang L., Kiguchi M., Murakoshi K. Theoretical investigation on the influence of temperature and crystallo-graphic orientation on the breaking behavior of copper nanowire // Phys. Chem. Chem. Phys. — 2009. — V. 11. — P. 6514—6519. — doi 10.1039/b902795e.

8. Wang F., Sun W., Gao Y., Zhao J., Sun Ch. Investigation on the most probable breaking behaviors of cupper nanowires with the dependence of temperature // Comp. Mater. Sci. — 2013. — V 67. — P. 182—187. — doi 10.1016/j.commatsci.2012.07.048.

9. Voter A.F. Embedded Atom Method Potentials for Seven FCC Metals: Ni, Pd, Pt, Cu, Ag, Au, and Al // Los Alamos Unclassified Technical Report # LA-UR 93-3901. — 1993.

10. Головнева Е.И., Головнев И.Ф., Фомин В.М. Моделирование квазистатических процессов в кристаллах методом молекулярной динамики // Физ. мезомех. — 2003. — Т. 6. — № 6. — С. 5—10.

11. Болеста А.В., Головнев И.Ф., Фомин В.М. Плавление на контакте при соударении кластера никеля с жесткой стенкой // Физ. мезомех. — 2001. — Т. 4. — №. 1. — С. 5—10.

12. Schneider T., Stoll E. Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions // Phys. Rev. B. — 1997. — V. 17. — No. 3. — P. 1302—1322.

13. Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F., DiNola A., Haak J.R. Molecular dynamics with coupling to an external bath // J. Chem. Phys. — 1984. — V. 81. — No. 8. — P. 3684—3690.

14. Andersen H.C. Molecular dynamics simulations at constant pressure and/or temperature // J. Chem. Phys. — 1980. — V. 72. — P. 2384—2393.

15. Nose S. A unified formulation of the constant temperature molecular dynamics methods // J. Chem. Phys. — 1984. — V. 81. — No. 1. — P. 511— 519.

16. Hoover W.G. Canonical dynamics: Equilibrium phase-space distributions // Phys. Rev. A. — 1985. — V. 31. — No. 3. — P. 1695—1697.

17. Martyna G.J., Klein M.L., Tuckerman M. Nose—Hoover chains: The canonical ensemble via continuous dynamics // J. Chem. Phys. — 1992. — V. 97. — No. 4. — P. 2635—2643.

18. Golovnev I.F., Golovneva E.I., Fomin V.M. The influence of the surface on the fracture process of nanostructures under dynamic loads // Comp. Mater. Sci. — 2015. — V. 97. — P. 109—115. — doi 10.1016/ j.commatsci.2014.10.022.

19. GolovnevI.F., GolovnevaE.I., Fomin V.M. Molecular-dynamics analysis of the rotary field formation in the nanostructure during stretching at a constant deformation velocity // Comput. Mater. Sci. — 2015. — V. 110. — P. 302—307. — doi 10.1016/j.commatsci.2015.08.012.

20. Головнев И.Ф., Головнева Е.И., Фомин В.М. Молекулярно-динами-ческое исследование роли поверхности в процессе разрушения наноструктур // Физ. мезомех. — 2014. — T. 17. — № 6. — C. 45—51.

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

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

Головнев Игорь Федорович, к.ф.-м.н., снс, снс ИГПМ СО РАН, golovnev@itam.nsc.ru Головнева Елена Игоревна, к.ф.-м.н., снс ИТПМ СО РАН, elena@itam.nsc.ru Уткин Андрей Вячеславович, к.ф.-м.н., снс ИТПМ СО РАН, utkin@itam.nsc.ru