Математическое моделирование пластической деформации кристаллов

МИНИСТЕРСТВООБРАЗОВАНИЯ И НАУКИ УКРАИНЫ
Харьковскийнациональный университет
им.В.Н. Каразина
“К ЗАЩИТЕ”
Заведующий кафедры материалов реакторостроения
Проф*******
ВЫПУСКНАЯРАБОТА
МАГИСТРАПРИКЛАДНОЙ ФИЗИКИ
МАТЕМАТИЧЕСКОЕМОДЕЛИРОВАНИЕ ПЛАСТИЧЕСКОЙ ДЕФОРМАЦИИ И РАЗРУШЕНИЯ ГПУ КРИСТАЛЛОВ Руководитель ************ Студент **************
Харьков
2008

Содержание
  Введение 1 Применения МД для исследования пластической деформации кристаллов 1.1 Ограничения 1.2 Потенциал 1.3 Алгоритм интегрирования по времени 1.4 Процедура минимизации 1.5 Вычисление сил 1.6 Периодичность 1.7 Начальное состояние 1.8 Начальное состояние для кристалла с дефектами 1.9 Нагрузка 1.10 Уравнение для ширины ячейки моделирования 1.11 Контроль системы 1.12 Вычисление физических величин 1.13 Визуализация 2 Моделирования пластической деформации ГПУ кристаллов Заключение Список использованных источников

Анотація
В даній роботі проведено аналіз особливостей застосування методумолекулярної динаміки для моделювання пластичності кристалів. Запропонованоновий підхід до моделювання розтягування кристалів.Запропоновано динамічне рівняння для поперечного розміру комірки моделювання. Створена програма для дослідження процесупластичної деформації та руйнування кристалів. Проведено моделювання розвиткупластичної деформації ГЩУ кристалів при одноосному розтягуванні. Показанапринципова можливість імітації за допомогою цього методу кривих розтягування досконалих кристалів,зміни температури зразка, появи дислокацій, полос ковзання, поодиноких вакансійта їх скупчень, а також процесу руйнування кристалів

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

Abstract
 
Inthe given article the analysis of features of application of a moleculardynamics for simulation of a plasticity of crystals is conducted. The newapproach to simulation of strain of crystals is offered. The dynamic equationfor transverse dimensions of simulation cell  is offered. The code forinvestigation of the process of a plastic deformation and destruction ofcrystals is created. The simulation of development of a plastic deformation hcpcrystals is carried  out at monoaxial expansion. The principal possibility ofimitation with the help of this method of stress-strain curves for the perfectcrystals, temperature variation of a sample, appearance of dislocations,stripes of slide, single vacancies and their clusters, and also process ofcrystal destruction.

Введение
 
Бурный ростядерной энергетики и широкое развитие работ в области термоядерного синтезапослужили мощным стимулом интенсификации исследований в области облучаемыхматериалов. Создана новая ветвь материаловедения – атомное, радиационноематериаловедение.
Свойстваматериалов всегда были ключевым звеном, определяющим успех многих инженерныхразработок в различных областях техники. Особенно их роль возросла в последнеевремя при создании сложных конструкций, работающих в экстремальных условиях.Ядерные реакторы, устройства термоядерного синтеза – яркий пример такихконструкций. Сотни различных по составу, структуре и способам изготовленияматериалов обеспечивают их работоспособность. Но, попадая в условия высокихпотоков облучения (до 1020 нейтр/м2∙с) и большихфлюенсов (до1027-1028 нейтр/м2), онипретерпевают значительные структурные перестройки (радиационное повреждение).Следствием этих перестроек является резкое изменение всех физических свойствматериалов. Причем эти изменения носят не совсем обычный характер. Ранее ничегоподобного не встречалось в обширной практике работ с различными материалами.Так, были обнаружены абсолютно новые явления, происходящие с облученнымиматериалами и сплавами: радиационное охрупчивание, радиационное распухание,радиационное упрочнение, ускоренная диффузия, радиационно-индуцированныефазово-структурные превращения и др.
И этот перечень,судя по всему, будет продолжен, так как исследования в области термоядерногосинтеза уже выдвигают новые требования и к конструкционным материалам, которыедолжны будут работать в еще более жестких условиях (например, выдерживатьоблучение нейтронами с энергией до 14 МэВ).
Одним из основныхфакторов, определяющих свойства материалов, как известно, является структура. Вусловиях облучения она претерпевает существенные перестройки на атомарномуровне. Чтобы выявить принципиальные закономерности поведения материалов в техили иных условиях эксплуатации, создать материалы с заданными свойствами,необходимо установить связь между изменяющейся атомарной структурой материалови всей совокупностью их макроскопических свойств. Для решения этой важнойзадачи радиационного материаловедения привлекаются самые современные физическиеметоды анализа структуры материалов, начиная с уже ставшего традиционнымрентгеноструктурного анализа и кончая автоионный микроскопией,оже-спектроскопией и др.
Но в значительном числе случаев все же не удается произвести необходимыеструктурные изменения. Такая ситуация возможна при известных, но быстроменяющихся внешних условиях (например, ударное воздействие на материал припролете высокоэнергетических нейтронов). Если же точно не известны изменения вструктуре материала на всех этапах какого-либо технологического процесса, тотрудно говорить о каком-либо надежном прогнозе его поведения.
В ряде случаевсовременная аппаратура из-за своего недостаточного разрешения не позволяетнаблюдать атомные перестройки в материалах, например отдельные атомные скачкипри диффузии, растворении и рост предвыделений второй фазы в сплавах и т.д.
Чтобы преодолетьперечисленные трудности и воссоздать быстроразвивающиеся процессы в материалах,перестройки на атомарном уровне или процессы, когда доступ к материаламограничен или опасен, все чаще в атомном материаловедении привлекаетсякомпьютерный эксперимент. Обзор методов компьютерного моделированияиспользуемых в радиационном материаловедении на начало 90-х годов можно найти в[1].
Методмолекулярной динамики (МД) — один из наиболее мощных методов, используемых длякомпьютерного моделирования в радиационном материаловедении. Он позволяетпроводить детальные исследования структуры материалов исходя из первыхпринципов. В последние годы метод молекулярной динамики переживает второерождение. Связано это, во-первых, с быстрым увеличением мощности компьютеров(быстродействия, объема памяти), что позволяет проводить беспрецедентные помасштабам компьютерные эксперименты. Так, сегодня уже проведены расчетыпластической деформации нанокристаллов меди, состоящих из 100 миллионов атомов[2]. С другой стороны, были разработаны более совершенные потенциалымежатомного взаимодействия, включающие в себя также многочастичные эффекты.
Данная работаявляется продолжением защищенной в 2003 г. бакалаврской работы, в которой были проанализировано современное состояние методов математической обработкикривых растяжения реакторных материалов, а также было проведено, на примеререакторных сталей, определение основных параметров, характеризующихпластическую деформацию.
Целью работы является изучение метода молекулярной динамики иособенностей его применения к исследованию пластичности реакторных материалов,разработка алгоритма и создание программы для изучения пластичности  двумерныхГПУ-кристаллов.

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

1.1. Ограничения
 
Укажем наограничения, свойственные этому методу. Во-первых, можно сразу же спросить:почему мы используем законы Ньютона, чтобы двигать атомы, хотя известно, чтосистемы на атомном уровне подчиняются скорее квантовым законам, чемклассическим?
Простейшаяпроверка применимости классического приближения базируется на тепловой длиневолны де-Бройля, определяемой как
/>, (1)
где /> - атомная масса и /> - температура.Классическое приближение хорошо работает когда />,где /> — среднее расстояние междуатомами. Если рассматривать, например, жидкость в тройной точке, то /> порядка 0.1 для легкихэлементов, таких как Li и Ar, уменьшаясь для более тяжелых элементов.Классическое приближение плохо работает для легких систем, таких как H2, He, Ne.
Кроме того,квантовые эффекты становятся важными для любых систем, когда температура /> достаточно низка. Падениеудельной теплоемкости кристаллов ниже температуры Дебая и аномальное поведениекоэффициента теплового расширения есть хорошо известные примеры измеримыхквантовых эффектов в твердых телах. Поэтому результаты, полученные с помощью МДдолжны интерпретироваться с известной осторожностью в этих областях.
Второеограничение связано с ограниченностью используемых компьютерных ресурсов, чтоприводит к ограничению количества рассматриваемых атомов и, как следствие, кснижению точности вычисляемых физических величин. Частично эту проблему можнообойти, используя подходящие граничные условия (см. ниже). Постоянный ростмощности компьютеров также способствует смягчению этой проблемы. Так, внастоящее время, в печати имеются сообщения о расчетах с 100 млн. атомов.
Необходимоотметить, что современные мощные суперкомпьютеры являются параллельными.Поэтому для расчетов на них необходимо обеспечить эффективное распараллеливаниевычислений. Об одной новой возможности распараллеливания в рассматриваемойздесь задаче будет сказано ниже.
Ограниченноебыстродействие компьютеров накладывает ограничения на скорость деформации,используемую в вычислениях. Это связано с тем, что при решении уравненийдвижения шаг по времени должен составлять порядка 0,01 шага от периодаколебаний атомов (по порядку величины равного/>сек)для обеспечения необходимой точности вычислений. Отсюда следует, что дляобеспечения деформации ~100% за приемлемое время вычислений типичная скоростьдеформации должна составлять порядка />сек/> в то время как максимальнодостигаемая в эксперименте скорость деформации составляет 105 сек-1.Кроме очевидной возможности достижения более низкой скорости деформации,увеличивая время вычислений, есть еще одна возможность – использовать процедуруминимизации.
При нулевойтемпературе система находится в локальном минимуме внутренней энергии, а из-заотсутствия тепловых колебаний она не может покинуть этот минимум.
Процедура минимизациипозволяет деформируемой системе находиться вблизи локального минимумавнутренней энергии. При таком моделировании время не определено, так как мы нерешаем уравнений движения. Следовательно, скорость деформации не оказываетникакого влияния, если она достаточно низка чтобы исключить нагрев системы, иобеспечить сходимость процедуры минимизации. Таким образом, моделирование,основанное на процедуре минимизации, представляет собой модель идеализированногоэксперимента при нулевой температуре в пределе низкой скорости деформации,когда тепло выделяемое при деформировании, удаляется.

1.2. Потенциал
 
Для моделированияматериала необходимо задать потенциал взаимодействия составляющих его атомов.Наиболее простым является парный потенциал Леннарда-Джонса
/>, (2)
Здесь, /> - расстояние междуатомами, /> - глубина потенциальнойямы и /> связано с положениемминимума потенциала />. ПотенциалЛеннарда-Джонса качественно правильно описывает взаимодействие между атомами –сильное отталкивание на малых расстояниях, обусловленное первым слагаемым вскобках, и притяжение на больших расстояниях, за которое отвечает второеслагаемое в скобках. Он хорошо описывает ван-дер-ваальсовское взаимодействие междуатомами кристаллов благородных газов, но, вследствие своей простоты, частоиспользуется для качественного описания взаимодействия других атомов. Спотенциалом Леннарда-Джонса проведено большое количество вычислений. Онявляется стандартным в вычислениях методом МД.
Основнымиматериалами реакторостроения являются металлы – сталь, цирконий и т.д. Вметаллах природа сил взаимодействия между атомами не двухчастичная (парная) амногочастичная.  Effective Medium Theory (EMT) дает реалистическое описаниеметаллической связи с учетом её многочастичной природы [3,4].  EMT- потенциал,с вычислительной точки зрения, не намного сложнее парного потенциала, но даетнамного более реалистическое описание свойств материалов. Поскольку в даннойработе не ставится задача изучения пластических свойств конкретного материаламы будем использовать модельный потенциал Леннарда-Джонса.
Удобно при этомвыбрать в качестве единицы длины />,единицы энергии /> и единицы массы- массу атомов /> (полагаем, чтоматериал состоит из атомов одного сорта). Это приводит к ускорению вычислений.Чтобы перейти к величинам, характеризующим конкретный материал, необходимоввести соответствующие масштабные множители — /> длядлины, /> для времени, /> для скорости, /> для силы, /> (/> в двумерном случае) длянапряжения, где /> и />взяты для данногоматериала.
ПотенциалЛеннарда-Джонса простирается до бесконечности. Однако на больших расстояниях онмал. И поэтому его влияние на движение далеких атомов мало. Чтобы ускоритьвычисления эту несущественную часть потенциала отбрасывают, или, другимисловами, вводят обрезание потенциала
/>. (3)
Радиус обрезания /> традиционно выбирают /> или />. Возможные причины такоговыбора будут обсуждаться ниже. В данной работе будет использоваться радиусобрезания />. Другие необходимыеизменения в потенциале  Леннарда-Джонса будут обсуждаться в разделе,посвященном выполнению закона сохранения энергии в МД.
1.3. Алгоритм интегрирования повремени
Основнымкомпонентом программ, использующих метод молекулярной динамики, являетсяалгоритм интегрирования по времени. Он необходим, чтобы проинтегрироватьуравнения движения взаимодействующих частиц и найти их траектории.
Алгоритминтегрирования по времени основывается на методе конечных разностей, время приэтом задается на конечной сетке, шаг по времени есть расстояние междупоследовательными точками сетки. Зная положения и скорости в момент времени /> (точные детали зависят оттипа алгоритма) схема интегрирования дает те же величины в более поздний моментвремени />. Используя процедуруинтегрирования временную эволюцию системы можно прослеживать в течениидлительного времени.
Конечно, этисхемы приближенными, и, поэтому, существуют ошибки, связанные с ними. Ониклассифицируются так:
Ошибки обрывания, связанные с точностью метода конечных разностей поотношению к истинному решению. Метод конечных разностей обычно базируется наряде Тейлора, оборванном на некотором члене. Эти ошибки не зависят отпрограммной реализации метода, они присущи самому алгоритму.
Ошибкиокругления, связаны с ошибками, возникающими при программной реализацииалгоритма. Например, такие ошибки возникают из-за конечного числа цифр,используемых в компьютерной арифметике.
Оба типа ошибокможно уменьшить, уменьшая />. Длябольших /> ошибки обрываниядоминируют, но они быстро уменьшаются, когда /> уменьшается.Например, алгоритм Верле имеет ошибки обрывания пропорциональные /> для каждого временногошага интегрирования. Ошибки округления падают более медленно с уменьшением /> и доминируют в пределемалых />. Использование 64-битнойточности (соответствующую “двойной точности” в Fortrane) помогает сохранитьошибки округления минимальными.
В молекулярнойдинамике наиболее часто используемым алгоритмом интегрирования по времениявляется, вероятно, так называемый алгоритм Верле [5]. Основная идея состоит втом, чтобы записать разложение Тейлора до третьего порядка вперед и назад повремени. Пусть /> обозначаетскорость, /> - ускорение и /> - третью производную от /> по />. Тогда имеем:
/>
/>. (4)
Складывая эти 2выражения получаем
/>. (5)
Это основнаяформула алгоритма Верле. Так как мы интегрируем уравнения Ньютона, то /> есть просто сила, деленнаяна массу, и сила в свою очередь есть функция положения />:
/>. (6)
Видно, что ошибкаобрывания алгоритма, когда система эволюционирует в течении времени />, есть величина порядка />, даже если третьяпроизводная не появляется в явном виде. Этот алгоритм в тоже время являетсяпростым для программной реализации, точным и стабильным, что объясняет егобольшую популярность при МД моделировании.
Проблема с этойверсией алгоритма Верле состоит в том, что скорости явно не вычисляются. Хотяони не нужны для временной эволюции, но их знание иногда необходимо. Крометого, они нужны для вычисления кинетической энергии />,чья оценка необходима чтобы проверить сохранение полной энергии />. Это один из наиболееважных тестов, указывающих, что МД моделирование выполняется корректно. Можновычислить скорости из положений использую формулу
/>. (7)
Однако, ошибки,которые дает это выражение, порядка /> а не />.
Чтобы преодолетьэту трудность, были развиты варианты алгоритма Верле. Они дают точно ту же траекториюи отличаются переменными, которые хранятся в памяти. Leap-frog алгоритм естьодин из таких вариантов.
Лучшийреализацией того же основного алгоритма есть так называемый алгоритм Верле соскоростью, когда положение скорости и ускорения в момент времени /> получается из тех жевеличин в момент времени />следующимобразом:
/>
/>
/>
/>. (8)
Заметим, чтонеобходимо 9N ячеек памяти, чтобы сохранить 3N положений, скоростей иускорений, но нам не нужно одновременно хранить значения любой из этих величиндля двух различных времен.
1.4. Процедура минимизации
Чтобымоделировать деформацию при нулевой температуре используется процедураминимизации, которая позволяет поддерживать систему вблизи локального минимумаэнергии все время. Деформация и минимизация выполняются одновременно. Алгоритмминимизации представляет собой модифицированный алгоритм МД. После каждого шагапо времени МД для каждого атома вычисляется скалярное произведение между импульсоми силой. Для атомов, скалярное произведение для которых отрицательно, импульсзануляется, так как эти атомы движутся в направлении, в котором потенциальнаяэнергия возрастает. Таким образом, кинетическая энергия атомов удаляется, тогдакак потенциальная энергия приближается к локальному минимуму энергии вдольнаправления движения атома. Такая процедура минимизации быстро сдвигает системув окрестность локального минимума энергии, но полной сходимости не получается,так как полная сходимость требует числа шагов по времени порядка числа степенейсвободы системы. Однако, обычно увеличении числа шагов процедуры минимизацииприводит лишь к малым изменениям в эволюции системы.
1.5. Вычисление сил
Наибольшихвычислительных усилий требует вычисление сил, действующих между атомами.Поэтому оптимизации алгоритма вычисления сил необходимо уделить особоевнимание. Один из шагов в этом направлении состоит в замене сложных длявычисления выражений для сил (например, содержащих экспоненту) на легковычисляемые выражения (например, сплайны третьего порядка). Второй шаг состоитв использовании потенциалов с ограниченным радиусом действия, или, какуказывалось выше, в обрезании несущественной области потенциала, если радиусдействия потенциала бесконечен. При этом необходимо вычислить только силы,действующие со стороны ближайших атомов, т.е. находящихся внутри сферы(окружности в двумерном случае) с радиусом равным радиусу обрезания />.
Третий шагсостоит в оптимизации алгоритма поиска атомов, ближайших к данному атому. Делов том, что прямолинейный перебор всех атомов, вычисление расстояний до них иотбрасывание тех атомов, расстояние до которых превышает радиус обрезания />, требует количестваопераций пропорционального />, где /> - число атомов в системе.Следовательно, с  ростом /> числотребуемых операций быстро возрастает, и поэтому выполнение вычислений сильнозамедляется, а, для больших />,делается практически невыполнимым. Таким образом, чтобы избежать этогозамедления нужен алгоритм, для которого число требуемых операций росло бы с /> линейно, а не квадратично.В принципе такой алгоритм прост – нужно перебирать не все атомы, а толькодостаточно близкорасположенные. Такое утверждение представляет собойтавтологию, пока не конкретизировано понятие близкорасположенных атомов. Чтобысделать это, разобьем ячейку моделирования на более мелкие субячейки. Тогдаблизкорасположенные к данному атому будут атомы, которые расположены всубячейках, соседних с субячейкой, содержащей данный атом или в субячейкахсоседних с соседними.
Удобно разбитьячейку моделирования на субячейки – параллелепипеды (прямоугольники в двумерномслучае). Вследствие сильного отталкивания на малых расстояниях, атомы не могутподходить близко друг к другу. Поэтому можно выбрать такие размеры субячеек,что в каждой из них будет находится не более одного атома.
Таким образом,алгоритм поиска атомов, удаленных от данного атома на расстояние не большерадиуса обрезания />, выглядитследующим образом. По номеру атома находим координаты атома и по ним субячейку,в которой находится атом. Затем находим субячейки, удаленные от нее нарасстояние не более чем />. Атомы,расположенные в этих субячейках, и будут искомыми (см. рис.1). Чтобы найтиномер атома, хранящегося в заданной субячейке, удобно ввести массив, каждыйэлемент которого соответствует определенной субячейке. В этом элементе массивабудет хранится номер атома, расположенного в этой субячейке, или нуль, еслисубячейка пуста. Элементы этого массива обновляются на каждом шаге по времениМД. Ясно, что изложенный алгоритм обеспечивает линейный рост числа операций сростом числа атомов /> в системе.Вариации этого алгоритма используются в программах МД “Gromex”[6], “MOLDY”[7],“DL_POLY”[8] и др.
Возможна и другаяорганизация вычислений, которая будет удобна для организации параллельныхвычислений. Именно для вычисления сил, действующих на данный атом, можноперейти от суммирования по близлежащим атомам, к суммированию по близлежащимсубячейкам (см рис.1). Будем двигаться последовательно по субячейкам первогоряда. Дойдя до конца первого ряда, перейдем в начало второго ряда и т.д. 1 2 3 4 5 6 7 8 9 10 11 12
Рис.1 Схемапоиска ближайших атомов.
Если в субячейкенаходится атом, то вычисляем силу,  действующую на него, со стороны ближайшихатомов, расположенных в близлежащих субячейках. Если же субячейка пуста, топереходим к следующей. Отметим при этом, что, например, для атома находящегосяв субячейке 6 (см. рис.1) необходимо вычислить силу, действующую со стороныатомов расположенных в субячейках 1, 2, 3, 7. Силы, действующие со стороныатомов, расположенных в субячейках 5, 9, 10, 11 в силу третьего закона Ньютона,с точностью до знака уже известны. Они были вычислены, когда вычислялись силы,действующие на атомы, расположенные в этих субячейках. Таким образом, в даннойорганизации вычислений, необходимо рассматривать лишь половину близлежайшихсубячеек. Далее, при переходе к смежной субячейке 7 нет необходимостиисследовать все близлежащие субячейки для поиска находящихся в них близкорасположенных атомов. Необходимо лишь исследовать ячейки 4 и 8. И к найденным вних атомам, добавить атомы, найденные для ячейки 6, за исключением атомовнаходящихся в субячейках 1 и 6. Таким образом, информация о ближайших атомахдля данной субячейки не теряется, а используется при поиске ближайших атомовдля смежной субячейки. Это естественно приводит к ускорению вычислений.


1.6.Периодичность
 
Число атомов,помещенных в ячейку моделирования, намного меньше числа атомов входящих всостав макроскопических систем. Чтобы результаты нашего моделирования можнобыло распространить на макроскопические тела, делают допущение, чтомакроскопические системы, состоят из бесконечного числа периодическиповторяющихся ячеек моделирования. Такая периодичность может быть в одном, двухи трех направлениях в трехмерном случае и в одном и двух в двумерном случае(см. рис.2). В этой работе мы будем рассматривать только двумерные системы. Этосвязано как с повышенными требованиями к вычислительным ресурсам в случаетрехмерных систем, так и с простотой визуализации результатов расчетов в двумерномслучае. В двумерном случае ячейка моделирования представляет собойпрямоугольник. В случае периодичности в одном направлении пара противолежащихсторон отождествляется, т.е. ячейку моделирования можно представить теперь какбоковую поверхность цилиндра. В случае периодичности в двух направленияхотождествляются обе пары противоположных сторон и ячейку моделирования можнотеперь представить как боковую поверхность тора. Если атом выходит за пределыячейки моделирования, то вследствие периодичности он входит в ячейку спротивоположной стороны.
1.7. Начальное состояние
В данной работебудут исследоваться с помощью МД кристаллы. Рассмотрим размещение совершенногокристалла в прямоугольной ячейке моделирования в случае периодичности в одномнаправлении. Периодическая структура самого кристалла накладывает ограниченияна размер ячейки моделирования в направлении периодичности. Действительно, еслив вершине, находящейся на одной из сторон ячейки моделирования находится атомкристалла, то эквивалентный атом кристалла должен быть в эквивалентной вершине,находящейся на другой из тождественных сторон. Это приводит к ограничениям навозможную длину ячейки моделирования в направлении периодичности и возможныеориентации кристаллографических осей кристалла относительно сторон ячейкимоделирования. Возможны такие ориентации кристалла, при которых указанное вышетребование выполнить точно невозможно.
/>
Рис.2 Периодичность ячеекмоделирования и размещение кристалла в ячейке моделирования.
Если жеориентация кристалла выбрана удачно, то длина ячейки моделирования можетпринимать значения кратные некоторой величине. Однако, на самом деле, этиограничения не очень существенны. Для всех длин ячейки моделирования иориентаций кристалла можно найти близкие к ним значения этих величин, длякоторых условие будет выполняться точно. Рецепт состоит в том, чтобы простосовместить указанную эквивалентную вершину с ближайшим эквивалентным атомомкристалла.
Если естьпериодичность (см. рис. 2) и по второму направлению, то должно выполнятьсяаналогичное требование и для второго направления. При этом необходимо заметить,что ориентация второй стороны  для прямоугольной ячейки моделирования ужезадана, поскольку она перпендикулярна первой стороне. Поэтому её длина будеткратна некоторой величине.
Если не принятьспециальных мер при подготовке начального состояния системы, то в ней возникаютколлективные движения — колебания. Это связано с тем, что система можетоказаться в сжатом или растянутом состоянии из-за несоответствия температурысистемы с постоянной кристаллической решетки. Другими словами это тепловоерасширение (сжатие) системы. Такие колебания имеют большой период и слабозатухают. Накладываясь на исследуемый процесс (например, деформированиесистемы) они смазывают картину этого исследуемого процесса. Следовательно, отэтих колебаний необходимо избавиться. Это можно сделать несколькими способами.Во-первых, подождать пока колебания затухнут. Однако из-за большого периода ималого затухания это требует большого времени. Во-вторых, попытаться подогнатьпостоянную решетки кристалла к температуре. Опыт показывает, что, сделавнесколько попыток, можно полностью исключить колебания. В-третьих, такуюподгонку можно выполнить автоматически. О том, как это можно сделать, будетсказано ниже.
В МДмоделировании часто возникает необходимость иметь систему в состоянии,характеризуемом определенной температурой. Однако, как мы можем получитьсистему с заданной температурой? Другими словами, как мы можем контролироватьсистему?
Для изменения температуры необходимо так изменить скорости частиц, чтобыполучить желаемую температуру. В алгоритме Верле со скоростью, обсуждаемомвыше, это может быть выполнено заменой уравнения
/> (9)
на уравнение
/>, (10)
где /> желаемая температура, и /> текущая температура. Такаямодификация означает, что мы больше не следуем уравнениям Ньютона и, что полнаяэнергия больше не сохраняется.
1.8. Начальное состояние длякристаллов с дефектами
С помощью МДможно исследовать деформирование, как совершенных кристаллов, так и кристалловсодержащих дефекты, например, кристаллов подвергнутых облучению. О том, какподготовить начальное состояние для совершенного кристалла, было сказано выше.Подготовка начального состояния для облученного кристалла намного более сложнаязадача. Однако, если известны доза и спектр первично выбитых атомов, МД позволяетвыполнить моделирование каскада повреждений [9,10,11]и таким образом решить этусложную задачу. При этом описанные выше потенциалы, необходимо дополнить, чтобыучесть отталкивание на малых расстояниях, например, гладко сшивая их спотенциалом Циглера-Бирсака-Литмарка [12]. Такой подход позволяет учесть многиеявления, возникающие при облучении, но является достаточно сложным и лежит зарамками данной работы.
Можно такжеисследовать влияние определенных дефектов, возникающие при облучении ГПУ кристалловна их пластические свойства. Например, можно исследовать влияние межузельныхкластеров и дефектов Френкеля. Очевидно, что начальные состояния, содержащиетакие дефекты, легко приготовить, стартуя с начального состояния для идеальногокристалла. Для этого необходимо удалить (добавить, переместить) атомы кристаллатак, чтобы получилась конфигурация кристалла с требуемыми дефектами. Кристаллпри этом получается обычно в напряженном состоянии. Это справедливо особеннопри добавлении атомов, так как для добавленных атомов расстояния до ближайшихатомов кристалла обычно намного меньше, чем равновесные расстояния междуатомами в кристалле. Из-за сильного роста потенциала межатомного взаимодействияна малых расстояниях такие атомы обладают большой потенциальной энергией. Еслине принять специальных мер, это может вызвать разлет кристалла. Чтобы недопустить этого и обеспечить релаксацию напряжений можно использовать процедуруминимизации и последующий подогрев системы до нужной температуры.

1.9.Нагрузка
 
В данной роботерассматривалось деформирование кристаллов путем одноосного растяжения.Поскольку вдоль направления растяжения наложены периодические граничныеусловия, то отсутствуют свободные границы, к которым можно было бы приложитьнагрузку. Поэтому задается растяжение системы, и потом находится возникшеевследствие этого напряжение. МД и деформирование выполняются одновременно.После каждого шага по времени МД выполняется малое растяжение, обеспечивающеенужную скорость деформации (/> наодном шаге). Растяжение выполнялось двумя способами. В первом, традиционноиспользуемом [13], система растягивается равномерно по длине. При этомкоординаты атомов вдоль направления растяжения умножаются на масштабныймножитель />. На этот же множительумножается длина ячейки моделирования. Согласно второму способу, предложенномув данной работе, растяжение концентрируется только возле торцов системы. Этотспособ лучше соответствует экспериментальной ситуации, когда нагрузкаприкладывается к торцам системы. При этом длина ячейки моделирования умножаетсяна масштабный множитель, а координаты атомов не умножаются.

1.10.Уравнение для ширины ячейки моделирования
 
Если боковыестороны системы по отношению к растяжению свободны, то нет необходимостиследить за шириной ячейки моделирования. Если же на систему наложеныпериодические граничные условия по двум направлениям, то изменению шириныячейки моделирования необходимо уделить особое внимание. При растяжениипоявляются сжимающие в поперечном направлении напряжения и поперечный размер(ширина системы) уменьшается. Если ширину ячейки моделирования не изменять, топоявится зазор, который будет увеличиваться со временем — система разорвется впоперечном направлении.
Один из подходов[13] состоит в умножении ширины ячейки моделирования на />, при увеличении длины в /> раз. Здесь /> примерно равнокоэффициенту Пуассона. Этого, однако, может оказаться недостаточно, поэтому,вводят дополнительную оптимизацию поперечного размера системы, основанную наметоде Монте-Карло. После каждых ~20 шагов по времени МД предлагается изменениепоперечного размера системы. Если в результате этого изменения энергия системыуменьшается, изменение принимается, в противном случае отклоняется. Вследствиеэтого, точное значение, выбранное для />,становится некритичным.
В данной работепредложен и используется другой подход, основанный на динамическом уравнениидля ширины ячейки моделирования. Выше уже было отмечено, что из-запериодичности в поперечном направлении система имеет топологию цилиндра.Сжимающие в поперечном направлении напряжения приводят к уменьшению боковойповерхности цилиндра и, следовательно, к уменьшению радиуса цилиндра. Записывая2-ой закон Ньютона для движения системы как целого вдоль радиуса, имеем
/>, (11)
где /> - напряжение в поперечномотносительно растяжения направлении, /> -площадь системы. Учитывая, что ширина ячейки моделирования />, имеем для неё уравнение
/>, (12)
Чтобы, исключитьколебательные процессы, удобно ввести в правую часть уравнения слабое фиктивноезатухание  />. Решая уравнения для /> на каждом временном шагеМД, мы поддерживаем ширину ячейки моделирование вблизи равновесного положения.Очевидно также, что данный подход можно использовать для демпфированияколебаний системы рассмотренных выше.

1.11.Контроль системы
 
Правильностьработы программы МД контролировалась с помощью закона сохранения энергии:
/>, (13)
где /> - кинетическая энергияатомов системы; /> - потенциальнаяэнергия их взаимодействия; />-работа, произведенная над системой. Выполнение закона сохранения энергии оченьважно при исследовании пластичности  твердых тел. Это связано с тем, что хотятепловое равновесие устанавливается быстро, но установление механическогоравновесия требует большого времени. Поэтому при деформировании системанаходится в тепловом равновесии, но скорее не находится в механическомравновесии, т.е. является неравновесной. Следовательно, потеря или приходэнергии, вследствие невыполнения закона сохранения энергии, может существенноповлиять на характер поведения системы при деформации.
Неточноесохранение энергии связано в основном с ошибками, возникающими из-за конечногошага интегрирования по времени, а также с ошибками, возникающими из-за конечнойточности представления чисел в компьютере.
Оба типа ошибокможно уменьшить, уменьшая шаг интегрирования по времени, что, однако,увеличивает время вычислений.
Другой тип ошибоквозникает из-за использования потенциала с обрезанием. Скачок потенциала нарадиусе обрезания при пластической деформации, когда атомы могут двигаться друготносительно друга, приводит к значительному нарушению закона сохраненииэнергии. Использование потенциала без скачка (3) позволяет существенно улучшитьвыполнения закона сохранения энергии. Потенциал (3), однако, имеет скачокпроизводной (силы) на радиусе обрезания />.Это также приводит к несоблюдению закона сохранения энергии. Оно особенно яркопроявляется при уменьшении радиуса обрезания от канонических значений /> и />. Это связанно с тем, чтоканонические значения радиуса обрезания находятся в минимумах радиальногораспределения атомов гексагональной решетки. Когда же /> попадает в максимумрадиального распределения число атомов, то испытывающих действие силы (при />), то прекращающихиспытывать ее действие (при />),становиться очень большим, что и приводит к существенному несохранению энергии.Чтобы избавиться от скачка производной потенциала на радиусе обрезания /> потенциал былмодернизирован. Пусть
/>, (14)
где
/> (15)
и />,/>, />. Тогда модернизированныйпотенциал имеет вид
/> (16)
Модернизированныйпотенциал гладко сшивается (до второй производной) с потенциаломЛеннарда-Джонса на радиусе сшивки /> изануляется вместе со своей первой производной нарадиусе обрезания />. С этимпотенциалом при значениях параметров /> и /> были проведены все расчетыв данной работе.

1.12.Вычисление физических величин
 
Придеформировании системы все физические величины, такие как напряжение />, температура />,кинетическая энергия />, потенциальнаяэнергия /> характеризующиедеформируемую систему меняются. Их мгновенные значения, усредненные по малымпромежуткам времени чтобы исключить тепловые колебания, описывают состояниедеформируемой системы. В отличие от равновесных систем мы не можем теперьиспользовать усреднение по времени, а должны использовать усреднение поразличным начальным состояниям системы.
Кинетическая ипотенциальная энергия находятся как
/> (17)
/> (18)
Температура определяется как
/>, (19)
где /> - размерность системы. Вдвухмерном случае /> - среднейкинетической энергией. Выражение для тензора напряжений, основанное навириальной теореме [14,15], имеет вид
/>, (20)
где /> - />-компоненты тензоранапряжений для атома />, /> - объем, приходящийся наатом /> (/>, где /> - полный объем системы), /> - масса атома />, /> - />-компонента его импульса, /> - расстояние между атомами/> и /> (/> - компонента вектора,направленного от />-го атома к />-му атому). Это выражениедля тензора напряжений не единственное, существуют и другие его определения.Однако, когда напряжения усредняются по объему различные определения быстросходятся к макроскопическому полю напряжений. Во время моделирования  кривыенапряжение — деформация строятся после усреднения атомного напряжения по всейсистеме.
1.13. Визуализация
МД позволяетполучать огромные объемы информации, описывающие исследуемую систему во всехдеталях. Поэтому возникает задача, извлечь из этого моря информации нужнуюинформацию и предоставить ее в виде, удобном для восприятия. Например, увидетьдефекты находящиеся внутри трехмерного кристалла, невозможно, поскольку ихзакрывают атомы, находящиеся на том же луче зрения, но ближе к наблюдателю.Однако, оставив только атомы, окружающие дефекты и удалив все остальные, этоможно легко сделать. Другая возможность состоит в использовании анимации дляисследования временной эволюции деформируемой системы.
Во времядеформирования кристаллов упорядоченное расположение атомов кристалланарушается, появляются дефекты. Для исследования локального атомного порядкаобычно используется алгоритм, известный как Common Neighbor Analysis (CNA)[16,17]. Для того чтобы определить структурукристалла в этом алгоритме исследуются связи между атомами и его соседями. Дваатома считаются связанными, если расстояние между ними меньше критическогорасстояния, выбранного между первыми двумя пиками в радиальной функциираспределения. Связи классифицируются с помощью трех целых чисел />. Первое из них, />, есть число общих соседей,т.е. атомов, связанных с обоими атомами в рассматриваемой связи. Второе, />, есть число связей междуэтими общими соседями. Третье, />, естьсамая длинная цепочка, которую можно образовать из этих связей.
 Число и тип />связей, которые имеет атом,определяют локальную кристаллическую структуру. Например,  атомы в совершенномГЦК кристалле имеет 12 связей типа 421, тогда как в совершенном ГПУ кристаллеимеют 6 связей типа 421 и 6 типа 422.
 ИспользованиеCNA позволяет сделать видимыми при моделировании дислокации, границы зерен идефекты упаковки. Например, при деформировании кристаллов меди, с помощью CNAклассифицируют атомы на три класса: ГЦК, ГПУ и “другие”, где в класс “другие”попадают атомы, которые имеют число связей, отличное от 12, или тип их связейотличен от 421 и 422. Тогда внутренние дефекты упаковки видны как двесопряженные (111) плоскости ГПУ атомов, несвойственные дефекты упаковки видныкак две (111) плоскости атомов ГПУ, разделенных (111) плоскостью атомов ГЦК.Границы двойников видны как одиночные (111) плоскости ГПУ атомов. Ядрадислокации и границы зерен  состоят из атомов класса “другие”, хотя границызерен содержат также немного ГПУ атомов.
 Тепловыеколебания атомов мешают выполнять CNA. Если при моделировании температурасравнительно низкая, достаточно тщательно выбрать критическое расстояние. Длявысоких же температур необходимо предварять CNA  короткой минимизацией(достаточной чтобы уменьшить тепловые колебания атомов, но избежать движениядислокаций).
 В двухмерныхсистемах нет нужды выполнять CNA анализ — дефекты видны непосредственно. Крометого, в отличие от трехмерных систем, отсутствуют заслоняющие атомы. Поэтомудвухмерные системы удобны для анимации, т.е.  воспроизведения  временнойэволюции деформируемой системы. Для создания анимации через заданное числошагов по времени МД, используя координаты атомов, формировалося изображениесистемы, которое затем сохранялось на жестком диске в bmp-файле. Анимация достигаласьвыводом этих изображений на экран дисплея в той же последовательности, вкоторой они создавались. Одно из преимуществ анимации — это наглядность. Полядругих физических величин, например, напряжения, температуры, используяподходящую кодировку, также можно использовать для анимации.

2.Моделирования пластической деформации ГПУ-кристаллов
 
Автором быласоздана программа для изучения пластичности в двумерных кристаллах. Двумерныесистемы были выбраны, чтобы обойти проблемы, связанные с высокими требованиямик вычислительным ресурсам в случае трехмерных систем и  визуализациейрезультатов вычислений. Для решения уравнений движения использовался алгоритмВерле со скоростью. В качестве потенциала взаимодействия между атомами былвыбран потенциал Леннарда-Джонса. При моделировании вычислительная ячейкарастягивалась вдоль умножением ее продольного размера на каждом шаге по временина масштабный множитель (1+ε), где ε – малое число (0.00001),выбранное так, чтобы обеспечить требуемую скорость деформации. Координатыатомов при этом не менялись, т.е. при этом вводился зазор между атомами всмежных ячейках моделирования. При этом нагрузка прикладывалась к торцамкристалла, что лучше соответствует эксперименту. Поперечный размер системыконтролировался с помощью динамического уравнения (12). Перед деформациейсистема приводилась в равновесное состояние с заданной температурой  короткимпрогоном с помощью МД. При вычислении кривых напряжение-деформация проводилосьусреднение напряжения по атомам всей системы. С помощью закона сохраненияэнергии контролировалась правильность работы программы. Для наблюдения завременной эволюцией деформируемого кристалла использовалась анимация.Блок-схема программы приведена ниже.
Задание входных данных
·  Шаг по времени
·  Число шагов по времени
·  Число атомов по x (направление деформации)
·  Число атомов по y
·  Скорость деформации ε
·  Температура
·  Тип решетки (гексагональная)
·  Ориентация решетки
·  Граничные условия (периодические)
·  Потенциал
Задание начальных значений
·  Начальные смещения атомов
·  Начальные скорости атомов Достижение равновесного состояния с заданной температурой Цикл по атомам (Вычисление начального ускорения)
Сила Fi, действующая на i-тый атом = сумме сил, действующих со стороны соседних атомов.
Ускорение i-того атома ai (t) = Fi/mi Конец цикла по атомам Цикл по времени (Траектория + растяжение) Цикл по атомам (Вычисление нового положения и промежуточной скорости)
/>
/> Конец цикла по атомам Проверка граничных условий Растяжение Цикл по атомам (Вычисление нового ускорения и новой скорости )
Сила Fi, действующая на i-тый атом = сумме сил, действующих со стороны соседних атомов.
Ускорение i-того атома ai (t+Δt) = Fi/mi
/> Конец цикла по атомам
Контроль параметров
·  деформация
·  напряжение
·  температура Конец цикла по времени Визуализация /> /> />
Программапозволяет проводить как обычную МД, так и использовать процедуру минимизации.Программа позволяет исследовать влияние ориентации кристалла, скоростидеформирования и температуры на процесс деформирования. Она также позволяетиспользовать свободные и периодические граничные условия на боковых,относительно растяжения, сторонах кристалла.
С помощьюпрограммы проведено моделирование пластического деформирования двумерныхкристаллов с гексагональной решеткой. Такие решетки характерны для базиснойплоскости ГПУ кристаллов и для плоскости (111) ГЦК кристаллов. Начальнаятемпература  />0.025/>. На рис.3 приведена суммаполной энергия кристалла и произведенной над ним работы (на один атом кристалла)как функция времени. Видно, что с точностью до 5 знака эта величинасохраняется, что свидетельствует о правильности работы программы.
/>
Рис.3 Зависимостьсумма полной энергия кристалла и произведенной над ним работы (на один атомкристалла) как функция времени.
На рис.4приведено напряжение, поперечное относительно направления растяжения, какфункция времени. Видно, что оно почти все время поддерживается невысоким. Этосвидетельствует о правильности выбранного алгоритма регулирования ширины ячейкимоделирования.
/>
Рис.4 Зависимостьнапряжения, поперечного относительно направления растяжения, от времени.
Наконец, на рис.5приведена кривая напряжение-деформация. Эта кривая полностью соответствуетэкспериментальным кривым напряжение-деформация для совершенных кристаллов.
/>
Рис.5 Криваярастяжения совершенного кристалла
На этой кривойможно выделить несколько стадий в деформировании кристалла [18,19]. На первой,упругой стадии, происходит накопление энергии. Затем, по мере роста напряжениякристалл переходит в нелинейную область, в конце которой появляются дислокации.Начинается сброс напряжения, который заканчивается появлением полос скольжения,которые и обеспечивают дальнейший сброс напряжения. Затем происходит частичноеупрочнение кристалла с последующими сбросами напряжения. На этой стадиипроисходит накопление вакансий, их слияние в вакансионные поры. На последнейстадии поры, объединяясь, дают начало трещине, которая, развиваясь, разрываеткристалл. На рис.6-11 даны изображения кристалла, которые иллюстрируютописанную выше стадийность пластического деформирования и разрушения кристалла.

Заключение
 
1. Изученметод молекулярной динамики со всеми его важнейшими ингредиентами: потенциалвзаимодействия, граничные условия, алгоритм интегрирования по времени, заданиеначальных условий, контроль термодинамических параметров в процессемоделирования, контроль достижения термодинамического равновесия, измерениефизических величин. Изучены особенности применения метода молекулярной динамикик исследованию пластичности реакторных материалов.
2. Созданапрограмма для изучения пластичности в кристаллах. Создана программавизуализации процесса пластической деформации и разрушения кристаллов.Предложен новый подход к моделированию растяжения кристаллов, близкий киспользуемому в эксперименте. Предложено динамическое уравнение для поперечногоразмера ячейки моделирования.
3. Проведеномоделирование развития пластической деформации ГПУ кристаллов при одноосномрастяжении. Показана принципиальная возможность имитации с помощью этого методакривых растяжения совершенных кристаллов, изменения температуры образца,появления дислокаций, полос скольжения, одиночных вакансий и их скоплений, атакже процесса разрушения кристаллов.

Списокиспользованных источников
1. В.В.Кирсанов,ЭВМ-эксперимент в атомном материаловедении, Энергоатомиздат, 1990.
2. J.Schiotz,Scripta Mater., 51, 837 (2004).
3. K.W.Jacobsen,J.K.Norskov, M.J.Puska, Phys. Rev. B35, 7423 (1987).
4. K.W.Jacobsen,P.Stoltze, J.K.Norskov, Surf. Sci. 366, 394 (1996).
5. L.Verlet,Phys. Rev. 159, 98 (1967); Phys. Rev. 165, 201 (1967).
6. D. vander Spoel, E. Lindahl, B. Hess at al. Gromacs User Manual, v. 3.2, Universityof Groningen.
7. K.Refson. Moldy User’s Manual, v. 2.16, (2001).
8. W.Smith,M.Leslie, T.R.Forester. The DL_POLY_2 User Manual, (2003), DaresberyLaboratory.
9. J.B.Gibson,A.N.Goland, M.Milgram, G.H.Vineyard, Phys. Rev. 120, 1229 (1960).
10. M.P.Allen,D.J.Tildesley Computer simulation of liquids. Clarendon Press, Oxford, 1989.
11. K.Nordlund,Comput. Mater. Sci., 3, 448 (1995).
12. J.F.Ziegler,J.P.Biersack, U.Littmark, The stopping and range of ions in solids. PergamonPress, N.Y., 1987.
13. J.Schiotz,T.Vegge, F. D. Di Tolla, K.W.Jacobsen, Phys. Rev. B60, 11971(1999),cond-mat/9902165.
14. T.Egami, K.Maede, and V.Vitek, Phil. Mag. A 41,  883 (1980).
15. J.R.Ray and A.Rahman, J. Chem. Phys.  80,  4423  (1984).
16. H.Jonsson and H.C.Andersen, Phys. Rev. Lett. 60,  2295 (1988).
17. A.S. Clarke and H.Jonsson, Phys. Rev. E  47,  3975 (1993).
18. Р.ХоникомбР. Пластическая деформация металлов. М. Мир. 1972.
19. И.М.Неклюдов,Н.В.Камышанченко. Физические основы прочности и пластичности металлов. ч. 3.Пластическая деформация и разрушение кристаллических материалов. Изд.“Педагогика-Пресс” и Белгородского государственного университета, Белгород,2001