CC BY
64
Моделирование квазистатических процессов в кристаллах методом молекулярной динамики
Е.И. Головнева, И.Ф. Головнев, В.М. Фомин
Институт теоретической и прикладной механики СО РАН, Новосибирск, 630090, Россия
В рамках метода молекулярной динамики введены искусственные силы вязкости, действующие на атомы, движущиеся в пространстве. Дополнительное использование модели плавного увеличения внешней нагрузки со временем, позволило провести моделирование конечных состояний нанообъектов, соответствующих квазистатическому нагружению. Приведены примеры расчета механических квазистатических характеристик бездефектных нанокристаллов меди.
1. Введение
Как известно, уравнения механики деформируемого твердого тела выводятся традиционно из законов сохранения массы, импульса и энергии. При этом приходится феноменологически вводить не только такие общие законы, как, например, закон Фурье для теплопроводности, зависимость ст(е), но и экспериментально определять ряд механических параметров, например модуль Юнга, коэффициент Пуассона и т.д. Для экспериментального нахождения таких параметров в случае макроскопических тел существует три основных типа экспериментов по механическому нагружению образцов: а) растяжение образца с постоянной скоростью движения границы; б) растяжение (сжатие) с постоянной внешней силой, описываемой функцией Хевисайда; в) приложенное напряжение возрастает с постоянной скоростью или, в предельном случае, бесконечно медленно (квазистатически). Однако провести такие эксперименты для наночастиц не представляется возможным, по крайней мере в настоящее время. Поэтому единственно возможным способом получения механических характеристик таких объектов является их численный расчет методом молекулярной динамики с использованием потенциалов, найденных в рамках квантовой механики. В предыдущей работе [1] авторы провели численный расчет для случая растяжения бездефектного кристалла при постоянной скорости движения свободной границы. Однако при моделировании процессов с медленно нарастающим
внешним воздействием методом молекулярной динамики встает проблема расчетного времени на компьютере. Поскольку шаг по времени т при интегрировании уравнений движения порядка 10-16 с, то реальная длительность моделируемых процессов на современных компьютерах составляет максимально 10-9-10-7 с. В связи с этим особую актуальность приобретает разработка способа моделирования квазистатического процесса в рамках метода молекулярной динамики.
2. Моделирование явлений в нанокристаллах, находящихся под воздействием внешней силы
Рассматривался трехмерный кристалл меди в форме прямоугольного параллелепипеда, атомы левой грани закреплены на неподвижной стенке, к атомам правой (незакрепленной) грани прикладывается растягивающая сила, действующая вдоль оси X. Рассмотрены два случая: а) F = F0Н(^ (И(1) — функция Хевисайда); б) F = (F0 /^ X При t > t0 полагаем F = F0. С помощью величины t0 можно регулировать скорость роста внешней силы. Подробное описание физической системы приведено в работе [2]. Отличие состоит в том, что зажим на правой грани моделируется потенциалом Морзе
Ж (г) = D
5Г V I / 5
+ Ds
е~2а‘(у -у>) - 2е
«„( Уi - у")
е-2«в(2< -Ч ) _ 2е~а‘(2' -2‘ )
+
© Головнева Е.И., Головнев И.Ф., Фомин В.М., 2003
закрепляющим атомы грани только в плоскости YZ. Это позволяет атомам грани двигаться вдоль оси X под действием внешней силы. После этого система охлаждается, а полученные координаты и импульсы атомов используются в качестве начальных данных.
Воздействие внешней силы учитывается внешним потенциалом:
Wf =-
]=N - N 8 +1
где N. — количество атомов на правой грани, а ^а — сила, действующая на отдельный атом. Используемая математическая модель и численный метод приведены в работе [3].
Для расчета локальных компонент тензора напряжения использовалось выражение, полученное в рамках кинетической теории [4]:
1
а В
1 ^+1 ра
і * а
(1)
Здесь ра — а-я компонента импульса /-го атома; ^ — а-я компонента силы, действующей между атомами с номерами / и j'; гЦ — в-я компонента радиус-вектора Гу; Ус — объем мезоячейки, в которой находится подсистема атомов. Форма и размеры мезоячейки подбирались исходя из физической и геометрической структуры рассматриваемой задачи. Так, в нашем случае, изучаемым объектом является однородный кристалл в форме прямоугольного параллелепипеда, вытянутого вдоль оси X, и внешнее воздействие также направлено вдоль оси X. Кристалл разбивался плоскостями, перпендикулярными оси X, между которыми находились две атомные плоскости. Разделяющие плоскости параллельны атомным плоскостям [100] ГЦК-структуры, и на протяжении всего процесса нагружения их расположение относительно атомных плоскостей не меняется, т.е. мезо-структуры рассматриваются в лагранжевых координатах. Для определения размеров использовалось условие совпадения плотности меди в мезоячейке с экспериментальными данными. В этом случае по оси X длина мезоячейки совпала с размером кристаллической ячейки а, а по осям Y и Z длина ребра параллелепипеда рассчитывалась по формуле: I = па + а/ 2 (п — число ячеек по оси Y или Z).
Наибольший интерес представляет продольная компонента тензора напряжения а хх, т.к. в рассматриваемой задаче моделируется одноосное растяжение вдоль оси X. Для контроля правильности выбора мезообъема расчет компоненты ахх дополнительно проводился по следующему выражению, имеющему простой физический смысл:
Здесь £—площадь грани мезоячейки, перпендикулярной направлению действия силы, а именно, оси X;
— х-я компонента силы, действующей на /-ю частицу со стороны j-й, при условии, что частицы находятся по разные стороны от грани. (Эту формулу правомерно использовать, когда тепловая часть мала.)
Компоненты тензора деформации в процессе нагружения определялись следующим образом. Для расчета длины мезоячейки Lcx проводилось усреднение по координатам х атомов левой и правой плоскостей, находилась их разность и добавлялась величина а/2. По осям Y и Z находились средние значения координат граничных атомов, бралась их разность и для расчета поперечных размеров Lcy и Lcz также добавлялась величина а/ 2. Далее расчет диагональных компонент тензора деформации проводился традиционным способом:
с = Lci - Lci 0
Lci 0
где Lci0, Lci — начальный размер мезоячейки по соответствующей оси.
Численные расчеты проводились для внешнего напряжения а0 в интервале значений от 0.5 до 20 ГПа. Прежде всего, было проведено сравнение результатов расчета компоненты ахх при использовании выражений (1), (2). На рис. 1 приведены графики зависимости этих величин от количества шагов по времени для мезоячейки, в которой находились 20-я и 21-я атомные плоскости. Как видно, в случае малых вкладов потокового члена в выражение (1) результаты идентичны. Для иллюстрации влияния линейного роста на локальные параметры на рис. 2 приведено сравнение зависимости а хх в той же ячейке от числа шагов по времени для импульсного нагружения внешней силой и в случае ее линейного роста. (Время возрастания нагружения t0 = 2 • 10-11 с.) В качестве тестового примера показано поведение внешнего приложенного напряжения на свободной границе. При импульсном возмущении отчетливо виден волновой характер процессов в кристалле, на который накладывается осцилляция атомных плоскостей. При линейном росте напряжения атомные осцилляции практически исчезли, но волновой характер движения сохра-
Рис. 1. Компонента тензора напряжения: сплошная линия — расчет по формуле Зубарева (1); штриховая линия — расчет по формуле (2)
ахх, ГПа
О 100000 200000 Ыт, 10“16с
Рис. 2. Сплошная черная линия — внешнее нагружающее напряжение; штриховая линия — напряжение внутри выделенной мезоячей-ки; серая линия — напряжение внутри выделенной мезоячейки, когда график внешнего напряжения имеет вид функции Хевисайда
нился, хотя и с гораздо меньшей амплитудой. На рис. 3 приведена зависимость скорости свободной границы стержня от числа шагов по времени. Как видно, осцил-ляционно-волновой характер процессов в кристалле влияет и на макрохарактеристики. Очевидно, выбор линейно возрастающего внешнего напряжения обеспечивает улучшение адекватности моделирования квазиста-тического процесса. Влияние длины интервала возрастания внешней силы F проиллюстрировано на рис. 4. Видно, что с ростом волновой характер процессов в кристалле меньше сказывается на получаемых характеристиках — увеличение t0 приводит к уменьшению амплитуды волновых процессов. Однако требование ограничения времени счета на компьютере обусловливает создание другой модели.
3. Введение искусственной вязкости в бездиссипативную механическую модель
В работе [3] проведено подробное обоснование правомерности использования бездиссипативных сил в пропагаторной механике. В реальном физическом эксперименте по квазистатическому растяжению стержня волновой характер процессов также присутствует, но на характерных временах проведения эксперимента энергия волнового движения диссипирует в тепловую энергию стержня, зажимов и т.д. При моделировании методом молекулярной динамики учесть все эти процессы невозможно и возникает необходимость ввести искусственную вязкость («вязкий эфир») для моделирования влияния внешней системы — термостата.
Предполагаем, что на каждый ^й атом действует сила со стороны «вязкого эфира»
^ =-ЧР1 ■
Очевидно, что введение такой силы не скажется на конечном состоянии, в котором должно выполняться условие статического равновесия. В то же время, силы вязкости должны погасить как атомные осцилляции, так и волновые процессы.
М\_, км/с 0.4-
Рис. 3. Скорость свободного конца стержня. Черная линия — линейное возрастание внешнего нагружения; серая линия — внешнее напряжение имеет вид функции Хевисайда
ахх, ГПа и
ю-
0 100000 200000 N.. 10-16с
Рис. 4. Черная штриховая линия — внешнее напряжение растет линейно, при этом t о = 2 • 10-11 с; черная сплошная линия — напряжение внутри выделенной мезоячейки при росте внешнего напряжения в течение временного интервала 10 = 2 • 10-11 с; серая штриховая линия — внешнее напряжение растет линейно, при этом 10 = 0.5 • 10-11 с; серая сплошная линия—напряжение внутри выделенной мезоячейки при росте внешнего напряжения в течение временного интервала 10 = 0.5 • 10-11 с
стхх, ГПа 11
0 40000 80000 Ыт, 10-16с
Рис. 5. Штриховая линия — внешнее напряжение в виде функции Хевисайда; черная линия — локальное напряжение при учете вязкости; серая линия — локальное напряжение без учета вязкости
Действительно, даже для импульсного нагружения внешней силой свободного конца стержня наблюдается сильное гашение как атомных осцилляций в ячейке (рис. 5), так и волнового движения свободной границы (рис. 6). Если же к вязким силам добавить линейный
Рис. 6. Черная линия — скорость свободной границы при учете вязкости; серая линия — скорость свободной границы без учета вязкости, при этом внешняя нагрузка подается в виде функции Хевисайда
Рис. 7. Черная линия — зависимость внешней приложенной силы от числа шагов по времени; штриховая линия — сила в сечении двадцатой ячейки с учетом вязкости; серая линия — сила в сечении двадцатой ячейки без учета вязкости
рост внешней силы, то процесс отвечает всем условиям квазистатики. На рис. 7 приведены зависимости внешней приложенной силы (черная линия) и силы в сечении двадцатой мезоячейки без вязкости (серая линия) и с вязкостью (штриховая линия). Как видно, введение вязкости позволило убрать атомные осцилляции и влияние волновых процессов. Сила в этой ячейке (как и в остальных) возрастает также линейно с задержкой на время прихода возмущения в мезоячейки. На рис. 8 приведена зависимость скорости свободного конца стержня от времени для случая линейного нагружения с учетом вязкости и без нее. Введение вязкости позволило полностью устранить волновые явления и получить в конце счета покоящийся нанокристалл. В качестве примера на рис. 9 приведена зависимость коэффициента Пуассона от числа шагов по времени для случая квазистатического нагружения с вязкостью. Видно, что в предельном случае установления статического равновесия его значение выходит на постоянную величину, совпадающую с экспериментальными данными для меди с точностью до 2%.
Разработанная методика введения искусственной вязкости для моделирования квазистатических процес-
сов позволила рассчитать предельную статическую а(е) диаграмму бездефектного нанокристалла меди. Численный эксперимент проводился следующим образом. На внешней границе задавалось напряжение, возрастающее линейно со временем. При этом были выбраны значения t0 = 2 • 10-11 с, а коэффициент вязкости V = = 0.1. После выхода на стационарное состояние статического равновесия бралось окончательное значение относительного удлинения стержня (рис. 10). Окончательная диаграмма а(е) приведена на рис. 11. Для сравнения там же приведена зависимость а(е), полученная при динамической нагрузке (движение свободной границы с постоянной скоростью) [1]. Видно, что в пределах точности расчетов эти функции совпадают. Это подтверждает вывод работы [1]: для бездефектных нанокристаллов функция а(е) не зависит от скорости деформации.
Для иллюстрации преимущества введения искусственной вязкости перед реальным учетом диссипации волновой энергии в зажимы за счет теплопередачи был проведен следующий численный эксперимент. К описанному выше нанокристаллу, было приложено внешнее напряжение, линейно возрастающее в течение ин-
Рис. 8. Зависимость скорости свободного конца стержня от времени для случая линейного нагружения с учетом вязкости (черная линия) и без нее (серая линия)
Рис. 9. Зависимость коэффициента Пуассона от времени для случая квазистатического нагружения с учетом вязкости в мезоячейке в середине кристалла
Рис. 10. Относительное удлинение стержня для различных значений внешнего растягивающего напряжения
Рис. 11. Сплошная линия — результат квазистатического моделирования с учетом вязкости; пунктирная линия — аппроксимирующий полином a = -106.735є2 + 108.5є
тервала t0 = 2 • 10-11 с до величины а 0 = 5 ГПа. Затем в течение времени t0 = 10-10 с (106 шагов по времени) это напряжение оставалось постоянным. Процесс теплопередачи через левый зажим моделировался следующим образом. Через каждые 200 шагов по времени импульсы атомов монослоя около неподвижного зажима толщиной в две кристаллические ячейки приравнивались нулю. В качестве примера на рис. 12 приведено поведение компоненты тензора напряжений аи относительного удлинения со временем во внутренней мезо-
ахх, ГП; 8
4
0
0 40000 80000 NT, 1016 с
8
0.08
0.04
0.00
0 40000 80000 Мт, 1016 с
Рис. 12. Зависимость компоненты тензора ахх (а) и относительного удлинения внутренней мезоячейки кристалла (б) от числа шагов по времени
ячейке кристалла. Отчетливо видно затухание волновых процессов. Компонента axx стремится к значению внешнего напряжения, а колебания относительного удлинения мезоячейки происходят около значения є = = 0.0538. При моделировании конечного квазистатичес-кого состояния с помощью введения искусственной вязкости предельное значение для случая a0 = 5 ГПа было равно 0.0528. Отличие составляет менее 2 %, однако время расчета с вязкостью в четыре раза меньше.
4. Заключение
Таким образом, искусственное введение сил вязкости, действующих на атомы при их движении в пространстве, и использование модели плавного увеличения внешней нагрузки со временем позволило провести моделирование конечных состояний нанообъектов, соответствующих квазистатическому нагружению. Проведенные расчеты проиллюстрировали возможность получения необходимых параметров для моделирования квазистатических процессов в наноструктурах методами механики деформируемого твердого тела. Предложенный подход позволил закрыть «белое пятно» — исследование состояний системы для квазистатических процессов методом молекулярной динамики. Отчетливо видна аналогия между взаимосвязью молекулярной динамики и континуальной механики и статистической физикой и термодинамикой. В рамках термодинамики возможно решение широкого класса проблем, если с помощью методов статистической физики будут получены уравнение состояния и выражение для внутренней энергии системы. В континуальной механике базовые уравнения удастся замкнуть, если в рамках молекулярной динамики рассчитать необходимые механические динамические и квазистатические характеристики изучаемых материалов.
Благодарности
Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант № 0201-00371), гранта Е02-4.0-222 по фундаментальным исследованиям в области естественных и точных наук Министерства образования Российской Федерации, гранта Президента Российской Федерации для поддержки ведущих научных школ № НШ-2282.2003.1.
Литература
1. Головнева Е.И., Головнев И.Ф., Фомин В.М. Молекулярно-динамичес-
кий анализ динамического разрушения наноструктур // Физ. мезо-мех. - 2003. - Т. 6. - № 2. - С. 37-46.
2. Головнев И.Ф., Конева Е.И., Фомин В.М. Численное моделирование разрушения бездефектных кристаллов при динамических нагрузках // Физ. мезомех. - 2001. - Т. 4. - № 5. - С. 5-11.
3. Головнев И.Ф., Головнева Е.И., Конев А.А., Фомин В.М. Физическая мезомеханика и молекулярно-динамическое моделирование // Физ. ме-зомех. - 1998. - Т. 1. - № 2. - С. 21-33.
4. Зубарев Д.Н. Неравновесная статистическая термодинамика. - М.: Нау-
ка. - 1971. - 415 с.
Simulation of quasistatic processes in crystals by a molecular dynamics method
E.I. Golovneva, I.F. Golovnev, and V.M. Fomin
Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, 630090, Russia
In the framework of a molecular dynamics method we introduce the artificial viscosity forces acting on atoms that move in space. The model of a gradual increase of the external load with time is additionally used, which allows us to simulate the finite states of nanoobjects corresponding to quasistatic loading. Examples are given of how mechanical quasistatic characteristics of defect-free copper nanocrystals are calculated.
