CC BY
54
Исследование процесса соударения сферического кластера меди с жесткой стенкой методом молекулярной динамики
А.В. Болеста, И.Ф. Головнев, В.М. Фомин
Институт теоретической и прикладной механики СО РАН, Новосибирск, 630090, Россия
В настоящей работе методом молекулярной динамики в его пропагаторной модификации изучаются процессы, происходящие при столкновении сферических кластеров меди с потенциальным барьером, имитирующим жесткую стенку в интервале начальных скоростей 100+1000 м/с. Численные исследования показали, что по характеру отражения кластера от стенки весь интервал начальных скоростей можно разбить на три области: квазиупругого удара, умеренно-необратимых деформаций и область почти полной диссипации начальной кинетической энергии. Проведен физический анализ процесса столкновения.
1. Введение
Одним из наиболее перспективных технологических способов обработки поверхностей является метод холодного газодинамического напыления [1]. При этом возникает целый ряд проблем оптимизации процесса, которые базируются на необходимости детального изучения динамики процессов соударения твердотельных кластеров с поверхностью. Особый интерес представляют явления в приповерхностных слоях на границе контакта кластеров с подложкой, так как они обусловливают появление механоактивационной связи. Так, в работах [1, 2] проведено исследование отскока сферической частицы размером 10+50 мкм от жесткой преграды на основе численного решения континуальных уравнений упругопластической среды. Основным эмпирическим параметром в задачах подобного рода является предел текучести, который может зависеть как от скорости деформаций, так и от размера частиц. В то же время, по имеющимся данным физические свойства частиц размером менее 100 нм (магнитные характеристики, теплоемкость, твердость и др.) существенно отличаются от свойств “массивных” частиц, исследуемых в континуальной механике [3]. В настоящее время в литературе отсутствуют данные об экспериментальных исследованиях высокоскоростного взаимодействия одиночных
микронных частиц с преградой, что препятствует проверке математической модели упругопластической среды и получению экспериментальных значений предела текучести в диапазоне скоростей до 1000 м/с и размеров частиц менее 40 мкм. Поэтому широкое применение для исследования неравновесных процессов в твердом теле получает метод молекулярной динамики. Так, в работе [4] исследуется процесс пластического деформирования ударной волной монокристалла с гранецент-рированной кубической решеткой. Показано, что за фронтом ударной волны кристалл необратимо деформирован, образуются дефекты упаковки — сдвиги вдоль всех возможных {111} кристаллографических плоскостей. Работы по моделированию взаимодействия сферических кластеров с поверхностью проводились в связи с приложением к задачам осаждения и имплантации пучков ионных кластеров, металлизации и очистки поверхности и т. д. [5, 6]. Диапазон скоростей в этих задачах соответствует кинетической энергии кластера от 100 эВ/атом и выше, что на несколько порядков выше, чем скорости частиц в процессе холодного газодинамического напыления. Таким образом, все это обусловило необходимость детального исследования на начальном этапе соударения кластеров с жесткой стенкой методом молекулярной динамики.
© Болеста А.В., Головнев И.Ф., Фомин В.М., 2000
2. Физическая система
Рассматривался процесс соударения трехмерных сферических кластеров меди с потенциальным барьером, имитирующим жесткую стенку. Выбор формы кластеров обусловлен экспериментальными исследованиями (см., например, [1]) и подробно обсуждался в работе [7]. Атомы расположены в узлах гранецентриро-ванной кубической решетки с постоянной решетки а = = 3.615А так, чтобы расстояние от каждого атома до центра кластера не превышало 23.5 А, при этом количество атомов в кластере N получается равным 4 921.
Для моделирования процессов в макротелах со свободной границей методом молекулярной динамики особую актуальность приобретает нахождение потенциалов межатомного взаимодействия. В случае одномерных кристаллов часто используются двухчастичные потенциалы, например потенциал Морзе [8]. Однако такой подход неприменим при моделировании процессов в 3D-кристаллах. Так, например, энергию металла невозможно представить в виде суммы парных членов, следовательно, невозможно таким способом одновременно описать энергию связи атома в бесконечном кристалле и энергетику конфигураций низкой симметрии (вакансии, поверхность). Кроме того, парные потенциалы предсказывают равенство упругих констант С12 и С44, которые значительно отличаются в реальных материалах. Это обусловливает необходимость применения многочастичных потенциалов межатомного взаимодействия. В настоящее время в задачах молекулярной динамики применяются потенциалы межатомного взаимодействия, полученные в рамках различных квантовомеханических моделей: метод внедренного атома [46], методы Финниса-Синклера, эффективной среды, метод модельного функционала электронной плотности. Все эти методы учитывают взаимодействие почти свободных электронов в металле с подсистемой остовных электронов, что выражается в многочастичности потенциалов межатомного взаимодействия. В работе выбран потенциал взаимодействия для атомов меди, предложенный Джонсоном [9], рассчитанный в рамках метода внедренного атома в приближении взаимодействия ближайших соседей:
и = £ F (р,) + 2 Х£ф(г,),
, , у
р, =Х / (гу )’
] &
/(г) = /еехР[-Р (г/ге - 1)], Ф(Г) = феехр[-У (г/ге - 1)], F(Р) = -Ес (1 - 1п х)х - 6фе7,
где х = (Р/Ре)ар; 7 = (Р/Ре); а, в, У, /е, Фе, Ес — константы; ге = а/42 — равновесное расстояние меж-
ду атомами; ре = 12 /е. Учет взаимодействия только ближайших соседей (в первой координационной сфере гранецентрированной кубической решетки находится 12 атомов) позволил значительно сократить время компьютерного счета.
В качестве потенциального барьера, имитирующего стенку, использовалась отталкивающая ветвь потенциала Леннарда-Джонса с параметрами е = 3 -10-14 эрг и
d= 1 А:
Ж(х,) = е У/х,)12.
Стенка расположена в плоскости х = 0, кластер налетает на стенку со стороны положительных значений х, начальная скорость кластера У0 отрицательна.
3. Модель и численный метод
Функция Г амильтона описанной выше системы имеет вид:
н = Е
С 2 Р 2 2 ^
Рх, + £у1 + Ря_
2т 2т 2т
V у
Численный расчет траекторий осуществлялся с помощью пропагаторной модификации метода молекулярной динамики [10, 11]. Была реализована схема второго порядка точности по временному шагу с использованием “списков Верлета”. Использовался целый ряд обез-размеривающих множителей: координаты измерялись
в 10
гия-
м
в
время — в 10
-13 с, масса — в 10-27 кг, энер-10-21 Дж, скорость — в 103 м/с, сила — в 10-11Н, давление и напряжение — в 109 Па. Шаг по времени менялся от 0.01 до 0.001 (от 10-15 до 10-16 с), так чтобы численная ошибка в конце счета, которая составляла 10-2 — 1 %, позволяла получить достоверные значения рассчитываемых макропараметров, число временных шагов при расчете одного процесса столкновения составляло от 104 до 105.
Моделирование начальных координат и импульсов кристалла при температуре Т> 0 К осуществлялось следующим способом. Вначале атомы располагались в положениях объемного равновесия с постоянной решетки, одинаковой для всех атомов и с нулевыми импульсами. Затем проводилось численное интегрирование динамики этой системы с периодическим обнулением всех импульсов частиц. Эта процедура позволяет получить кристалл с температурой Т=0. Далее проводился разогрев системы хаотически направленными импульсами силы постоянной амплитуды до любого заданного значения температуры. Параметры действующей случайной силы подбирались в соответствии с требованием равновесного распределения импульсов по Максвеллу в каждый момент процесса разогрева. Другим опробованным способом получения кристалла с заданной температурой является разброс начальных импульсов и координат методом Монте-Карло в соответствии с рас-
Рис. 1. Зависимость от времени: А — силы взаимодействия со стенкой В — скорости центра масс кластера Ус. Начальная скорость кластера У0 = 200 м/с
пределением Гиббса. Однако для получения распределения координат этим методом, ввиду значительной трудности диагонализации матрицы вторых производных потенциала, зависящего от 3N координат, необходимо прибегать к различным приближениям. Например, можно пренебречь производными потенциала вида д 2и/дх, дху, где номера частиц / и j соответствуют разным частицам. Полученные таким образом начальные координаты и импульсы отвечают некоторой неравновесной системе, которой требуется некоторое время для релаксации к равновесному состоянию. С учетом вышеизложенного, авторами был выбран первый способ моделирования начальных координат и импульсов атомов кристалла при заданной температуре.
4. Результаты
Расчеты проводились в интервале начальных скоростей ¥0 = 100+1000 м/с и начальных температур кристалла Т = 0+300 К. Для анализа динамики процесса столкновения рассчитывались следующие макропараметры системы: скорость центра масс кристалла сила взаи-
модействия со стенкой полная внутренняя энергия кристалла Е{п, ее кинетическая Ek и потенциальная и составляющие. Приведенные ниже результаты для большей наглядности относятся к начальной температуре Т=0К.
Рис. З. Изображение кластера (проекция на плоскость XY), охлажденного после столкновения со стенкой с начальной скоростью V0 = = 2GG м/с
По характеру отражения кластера от стенки весь интервал скоростей можно разбить на три области.
1. V0 < 300 м/с — область квазиупругого удара. B качестве примера ниже приведены результаты для V0 = =2GG м/с. На рис. l дана зависимость скорости центра масс кластера и силы взаимодействия со стенкой от времени. Хорошо видны моменты начала и конца взаимодействия (Fw = G и Vc = const), что позволяет определить и длительность соударения Tc. Как видно, конечная скорость несколько ниже начального значения 2GG м/с, то есть часть кинетической энергии центра масс дисси-
Рис. 2. Зависимость от времени: А — внутренней энергии Е1п ; В — ее кинетической составляющей Ek; С — ее потенциальной составляющей и. Начальная скорость кластера У0 = 200 м/с
Рис. 4. Изображение атомов кластера, охлажденного после столкновения со стенкой с начальной скоростью У0 = 200 м/с, в экваториальном шаровом слое высотой 2 А, лежащем в плоскости ХУ. Прозрачные кружки соответствуют атомам с потенциальной энергией больше -550
Рис. 5. Зависимость от времени: А — силы взаимодействия со стенкой В — скорости центра масс кластера V . Начальная скорость кластера У0 = 400 м/с
пирует во внутреннюю энергию кристалла (рис. 2). На этом же рисунке приведены тепловая энергия Ек (кинетическая энергия атомов в системе центра масс кристалла) и изменение потенциальной энергии взаимодействия атомов Ли = и - и(ґ = 0). Как видно, после столкновения эти энергии сравнялись и каждая составляет половину изменения внутренней энергии ЛЕіп = Еіп --Еіп (ґ = 0), что является необходимым условием равновесности конечного состояния.
На рис. 3 изображен кластер после столкновения в момент времени £ = 100 (проекция в плоскости ХУ), охлажденный до нулевой температуры тем же способом (см. выше), что и при подготовке начальных координат и импульсов. Видна правильная ГЦК-структура решетки, то есть отскок произошел без нарушения идеального строения кристалла. На рис. 4 изображены атомы в экваториальном шаровом слое высотой 2 А, лежащем в плоскости ХУ. Прозрачные кружки соответствуют атомам с потенциальной энергией больше определенного значения (-550). Видно, что все они сосредоточены на поверхности кластера, что соответствует наличию определенной поверхностной энергии.
2. 300 м/с < У0 < 500 м/с — область умеренно-необратимых деформаций. Ниже приведен анализ столкно-
Рис. 7. Изображение кластера (проекция на плоскость ХУ), охлажденного после столкновения со стенкой с начальной скоростью У0 = = 400 м/с
вения кластера со скоростью ¥0 = 400 м/с. На рис. 5 изображена зависимость скорости центра масс кластера и силы взаимодействия со стенкой от времени. В отличие от первого интервала скоростей здесь наблюдается значительная диссипация кинетической энергии центра масс, а максимальное значение силы возросло вдвое. Доля тепловой энергии Ek после столкновения (рис. 6) значительно снизилась по сравнению с запасенной потенциальной энергией А и и, как видно, обмен энергией между этими долями прекратился. Конечная структура кластера (рис. 7), хотя его форма и близка к сферической, претерпела значительные изменения. Видны харак-
Рис. 6. Зависимость от времени: А — внутренней энергии Еіп; В — ее кинетической составляющей Ek; С — ее потенциальной составляющей и. Начальная скорость кластера У0 = 400 м/с
Рис. 8. Изображение атомов кластера, охлажденного после столкновения со стенкой с начальной скоростью У0 = 400 м/с, в экваториальном шаровом слое высотой 2 А, лежащем в плоскости ХУ. Прозрачные кружки соответствуют атомам с потенциальной энергией больше -550
Рис. 9. Изображение атомов кластера, охлажденного после столкновения со стенкой с начальной скоростью ¥0 = 400 м/с, в экваториальном шаровом слое высотой 2 А, лежащем в плоскости, полученной вращением плоскости XX вокруг оси X на 45 °. Прозрачные кружки соответствуют атомам с потенциальной энергией больше -550
терные границы конуса, составляющие примерно угол 45° к оси X, отделяющие область с нарушениями идеальной ГЦК-решетки от “нормальной” структуры в “конце” кластера. Наибольшую качественную информацию дает энергетический анализ атомов в шаровом слое (рис. 8 и 9). Видно формирование конических границ с повышенной потенциальной энергией атомов, обусловленной нарушением идеальной структуры. Сравнение шарового слоя, лежащего в плоскости XY (рис. 8), и слоя, лежащего в плоскости, полученной вращением плоскости XY вокруг оси X на 45° (рис. 9), позволяет заключить, что в кристаллической решетке произошли сдвиги вдоль плоскостей {111}.
Результаты для начальной скорости V0 = 500 м/с показывают характер изменения процессов во втором интервале скоростей. Если для зависимости скорости кластера и силы взаимодействия с барьером от времени (рис. 10) мы имеем результаты аналогичные столкно-
Рис. 11. Зависимость от времени: А — внутренней энергии Еіп; В — ее кинетической составляющей Ек; С — ее потенциальной составляющей и. Начальная скорость кластера V0 = 500 м/с
вению со скоростью 400 м/с, то в поведении энергетических характеристик (рис. 11) заметны новые аспекты. Прежде всего необходимо отметить, что абсолютное значение тепловой энергии хаотического движения слабо изменилось по сравнению с приростом потенциальной составляющей внутренней энергии, что говорит о повышении роли формирования дефектов кристаллической структуры в процессе диссипации энергии упорядоченного движения центра масс кластера. Кроме того, в момент времени t = 25, после резкого падения силы взаимодействия с барьером, формируется плато на графике полной потенциальной энергии кластера, то есть не видно столь сильного влияния упругой составляющей. Форма кластера после отскока так же близка к сферической (рис. 12). Энергетический анализ в шаровом слое (рис. 13) показывает, что произошло формирование множества микрокристаллитов сложной формы с границами, образованными пересечением ряда кристаллографических плоскостей {111}.
Рис. 10. Зависимость от времени: А — силы взаим°действия со стен- Рис. 12. Изображение кластера (проекция на плоскость ХУ), охлаж-
кой ^; В — скорости центра масс кластера Ус. Начальная скорость денного после столкновения со стенкой с начальной скоростью У0 =
кластера У0 = 500 м/с =500 м/с
Рис. 13. Изображение атомов кластера, охлажденного после столкновения со стенкой с начальной скоростью У0 = 500 м/с, в экваториальном шаровом слое высотой 2 А, лежащем в плоскости ХУ. Прозрачные кружки соответствуют атомам с потенциальной энергией больше -550
3. У0 > 500 м/с. Эта область скоростей характеризуется, прежде всего, почти полной диссипацией начальной кинетической энергии центра масс. Далее в качестве примера приведена динамика столкновения с начальной скоростью У0 = 800 м/с. Из рис. 14 видно, что скорость центра масс кластера в конце процесса столкновения близка к нулю. Несмотря на резкое падение силы взаимодействия с барьером в момент времени t = = 15, наблюдается дальнейший значительный рост потенциальной энергии кластера с формированием плато (рис. 15). Время отскока резко возрастает, так как конечная скорость кластера близка к нулю и он находится длительное время в поле потенциального барьера. Геометрическая форма кластера после столкновения сильно отличается от сферической и близка к цилиндрической (рис. 16). Анализ структуры и энергетики кристалла в шаровом слое (рис. 17) показывает, что вся центральная часть находится в квазиаморфном состоянии с
Рис. 15. Зависимость от времени: А — внутренней энергии Е 1п; В— ее кинетической составляющей Ек; С — ее потенциальной составляющей и. Начальная скорость кластера У0 = 800 м/с
повышенной потенциальной энергией атомов, а начальная структура ГЦК-решетки сохранилась в форме тонкостенного цилиндра на периферии кластера.
Большой интерес представляет анализ поведения важнейших макрохарактеристик в зависимости от начальной скорости. Прежде всего, это относится к выделению области необратимых деформаций. Физически это явление соответствует нахождению исследуемой системы в конечном состоянии в локальном экстремуме потенциальной энергии и, который лежит выше начального значения и0 при температуре Т = 0 К, соответствующего идеальному кристаллу:
Дийт = и^ Т=0К - ио| Т=0К-
Для расчета конечной энергии и^ Т=0к проводилось охлаждение кристалла, отскочившего от стенки, тем же способом, что и при получении начальных данных (см. выше). На рис. 18 представлена зависимость необратимого приращения потенциальной энергии крис-
Рис. 14. Зависимость от времени: А — силы взаимодействия со стенкой Fw; В — скорости центра масс кластера V . Начальная скорость кластера У0 = 800 м/с
Рис. 16. Изображение кластера (проекция на плоскость ХУ), охлажденного после столкновения со стенкой с начальной скоростью У0 = = 800 м/с
У
X
40
Рис. 17. Изображение атомов кластера, охлажденного после столкновения со стенкой с начальной скоростью У0 = 800 м/с, в экваториальном шаровом слое высотой 2 А, лежащем в плоскости ХУ. Прозрачные кружки соответствуют атомам с потенциальной энергией больше -550
талла Ли1гг от начальной скорости. Как видно, при У0 < < 300 м/с эта величина равна нулю, то есть деформации носят структурно обратимый характер, а само столкновение является квазиупругим. Термин “квазиупругий” здесь применяется ввиду наличия диссипации части начальной кинетической энергии центра масс в тепловую (рис. 19). Коэффициент диссипации энергии равен
Лес = Ес0 -Ее? -100%,
С Ес0
где Ес0 и Ес£ — кинетическая энергия центра масс кристалла соответственно в начале и в конце процесса столкновения, и достигает 66 % при У0 = 300 м/с. Второй интервал скоростей от 300 до 500 м/с характеризуется ростом коэффициента диссипации почти до 100 % при уменьшении доли конечной тепловой энергии в полной внутренней энергии (рис. 20). Следовательно, увеличение диссипации происходит за счет нарушения идеальной структуры кристаллической решетки. Свыше
ДЦгг-10'21, Дж 300000
150000
о
100 300 500 700 900 У0, м/с
Рис. 18. Зависимость необратимого приращения потенциальной энергии кристалла Ли1п. от начальной скорости
ЛЕС, %
100
50
0
100 300 500 700 900 У0, м/с
Рис. 19. Зависимость коэффициента диссипации энергии кристалла ЛЕс от начальной скорости
500 м/с начинается близкий к линейному рост конечной тепловой энергии при сохранении почти квадратичной зависимости потенциальной энергии (рис. 20), накапливаемой при необратимых деформациях кристаллической решетки (рис. 18).
Представляет интерес качественное сравнение времени столкновения кластера тс, полученного в рамках дискретно-атомной механики (рис. 21), с результатами континуальной механики. В работах [1, 2] приведены графики зависимости времени контакта сферических частиц из меди и алюминия с твердой преградой при ударе, рассчитанные для упругих и упругопластических частиц. В случае упругих частиц результаты отлично согласуются с приведенными на рис. 21 в первой области скоростей столкновения 100 м/с < V < 300 м/с (с учетом линейной зависимости тс от радиуса частиц). Эти значения тс близки к аналитическим результатам, полученным в квазистатическом приближении в рамках контактной теории Герца. Однако в случае упругопластических частиц имеются некоторые различия в результатах. Прежде всего, авторы [1, 2] указывают, что уже при скорости удара 100 м/с частицы испытывают пластическую деформацию и, в результате, на графике зависимости времени столкновения от скорости виден минимум при скорости, равной приблизительно 150 м/с. Далее время контакта начинает возрастать с увеличением скорости, что отличается от полученной зависи-
ЛЕ|П'10'21, Дж 500000
250000
0
100 300 500 700 900 У0, м/с
Рис. 20. Зависимость от начальной скорости: А — конечной внутренней энергии Еп; В — ее кинетической составляющей Ек; С — ее потенциальной составляющей и{
тс-Ю'13, с
100 300 500 700 900 V0, м/с
Рис. 21. Зависимость времени столкновения кластера тс со стенкой от начальной скорости
мости на рис. 21 (время контакта испытывает минимум при У0 = 500 м/с). Для комментария к этим отличиям следует обратиться к упомянутым во введении существенным различиям в физических свойствах нанокрис-таллических кластеров и “массивных” частиц. Авторы [1, 2] указывают, что при увеличении предела текучести минимум на графике зависимости тс от У0 достигается при больших начальных скоростях.
5. Заключение
Динамическое исследование процесса соударения кластера с жесткой бесструктурной стенкой позволило выявить три интервала начальных скоростей У0, принципиально отличающихся по механизму диссипации кинетической энергии центра масс Ес0. В интервале до 300 м/с наблюдается перекачка энергии Ес в хаотическую энергию теплового движения атомов кристалла Е к и их потенциальную энергию взаимодействия и, значения которой, как показали расчеты, соответствуют упругому приближению. Действительно, при отражении от стенки не наблюдается перенос энергии Е к в энергию центра масс Ес в соответствии со вторым началом термодинамики, а вся конечная энергия Ес набирается за счет запасенной упругой составляющей потенциальной энергии в момент, когда Ес была равна нулю. Кластер сохраняет только такую долю и, которая равна Ек, также в соответствии с основными положениями статистической теории равновесных состояний. С увеличением скорости У0 (интервал 300+500 м/с) канал переноса энергии в хаотическое тепловое движение практически исчерпан и вся энергия Ес запасается в потенциальной энергии взаимодействия атомов и. Ее величина уже такова, что после прекращения взаимодействия
кластера со стенкой образуются выделенные плоскости, в которых атомы имеют повышенную остаточную потенциальную энергию. Система находится в локальных экстремумах, значения которых значительно выше начальной энергии U0. Конечная энергия центра масс, как и время взаимодействия со стенкой, монотонно падают с увеличением V0 в этом интервале, но, в то же время, сохраняются конечная форма кластера и ГЦК-структура решетки. Это область умеренных упругопластических деформаций. И наконец, при Vc > 500 м/с значение начальной кинетической энергии Ec достигает такой величины, что ее полная диссипация в момент остановки сопровождается нарушением идеальной кристаллической структуры практически во всем объеме с одновременной сильной деформацией начальной формы. Таким образом, достигается попадание системы в очень высокие локальные экстремумы потенциальной энергии кристалла. Это область сильных пластических деформаций.
Литература
1. Алхимов А.П., Гулидов А.И., Косарев В.Ф., Нестерович НИ. Осо-
бенности деформирования микрочастиц при ударе о твердую преграду // ПМТФ. - 2000. - Т. 41. - № 1. - С. 204-209.
2. Гулидов А.И., Фомин В.М., Фурсенко Р.В. Время контакта сферической частицы с твердой преградой при ударе // Численные методы решения задач теории упругости и пластичности: Труды 16 межресп. конф. - Новосибирск, 1999. - С. 66-70.
3. Гусев А.И. Эффекты нанокристаллического состояния в компактных металлах и соединениях // УФН. - 1998. - Т. 168. - № 1. -С. 55-83.
4. Holian B.L., Lomdahl PS. Plasticity induced by shock waves in nonequilibrium molecular-dynamics simulations // Science. - 1998. -V. 280. - P. 2085-2088.
5. Yorizane K., Muramoto T, Yamamura Y. Computer studies on reflection and sputtering due to low energy cluster impacts // Nucl. Instr. and Meth. in Physics Research B. - 1999. - V. 153. - P. 292-297.
6. Hsieh H., Averback R.S., Sellers H., Flynn C.P. Molecular-dynamics simulations of collisions between energetic clusters of atoms and metal substrates // Phys. Rev. B. - 1992. - V. 45. - No. 8. - P. 4417-4430.
7. Зенгуил Э. Физика поверхности. - М.: Мир, 1990. - 536 с.
8. Царегородцев А.И., Горлов Н.В., Демьянов Б.Ф., Старостен-ков М.Д. Атомная структура антифазной границы и ее влияние на состояние решетки вблизи дислокации в упорядоченных сплавах со сверхструктурой L12 // Физика металлов и металловедение. - 1984. - Т. 58. - Вып. 2. - С. 336-343.
9. Johnson R.A. Alloy models with the embedded-atom method // Phys. Rev. B. - 1989. - V. 39. - P. 12554-12559.
10. Головнев И.Ф., Головнева Е.И., Конев А.А., Фомин В.М. Физическая мезомеханика и молекулярно-динамическое моделирование // Физ. мезомех. - 1998. - Т. 1 - № 2. - С. 21-33.
11. Головнев И. Ф., Головнева Е.И., Фомин В.М. Моделирование процессов соударения твердых тел методом молекулярной динамики // Докл. РАН. - 1997. - Т. 356. - № 4. - С. 466-469.
