ОПРЕДЕЛЕНИЕ ВРЕМЕННОГО ШАГА ИНТЕГРИРОВАНИЯ ПРИ МОДЕЛИРОВАНИИ ЛАЗЕРНОГО
_____________ВОЗДЕЙСТВИЯ НА МЕТАЛЛЫ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ______________________
УДК 519.6, 519.85
ОПРЕДЕЛЕНИЕ ВРЕМЕННОГО ШАГА ИНТЕГРИРОВАНИЯ ПРИ МОДЕЛИРОВАНИИ ЛАЗЕРНОГО ВОЗДЕЙСТВИЯ НА МЕТАЛЛЫ МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
БЕСОГОНОВ ВВ., АНДРЕЕВ В. Вен., АНДРЕЕВ В.Вяч.
Институт прикладной механики УрО РАН,426067, г. Ижевск, ул.Т.Барамзиной, 34
АННОТАЦИЯ. Рассмотрены способы оценки и автоматического расчета оптимального размера временного шага интегрирования при решении задач методом молекулярной динамики для ансамблей атомов с большим диапазоном скоростей, характерным для процессов лазерной обработки материалов.
КЛЮЧЕВЫЕ СЛОВА: - молекулярная динамика, шаг интегрирования, лазерная обработка.
ВВЕДЕНИЕ
Нанокристаллические металлические материалы в настоящее время широко исследуются экспериментально, теоретически и методами математического моделирования. Рассмотрим воздействие на нанопленки и наночастицы монохроматического импульсного лазерного излучения. Особенностью воздействия лазерного излучения на металлы является передача скин-слою атомов достаточно большой энергии за малый промежуток времени. Физические процессы лазерного воздействия на металлы на макро-уровне экспериментально достаточно хорошо исследованы [1 - 3]. Однако экспериментально трудно, а зачастую невозможно исследовать процессы, происходящие в металле в период времени ( 10-13 ^ 10-9) с после воздействия. Наиболее удобным инструментом исследования процессов в этом интервале времени является эффективно используемый в настоящее время метод молекулярной динамики (ММД).
ОСНОВНЫЕ УРАВНЕНИЯ МОЛЕКУЛЯРНОЙ ДИНАМИКИ
Основой метода молекулярной динамики является расчет траекторий движения отдельных частиц, взаимодействие между которыми задается потенциалом взаимодействия. Движения частиц описываются уравнениями классической механики Ньютона (1)
d2 г. (г) ^ ч
). (1)
где . = 1,2,3,...,N; N - число частиц; г. - радиус-вектор . -ой частицы; Mi - масса . -ой частицы; ¥{ - суммарная сила, действующая на ьую частицу со стороны остальных частиц.
Сила ¥ определяется как производная от потенциала по расстоянию между частицами (2)
¥ (Г) = ~щр, (2)
дг
где Г = {г1, г2, г3,..., гп}; Г. - радиус-вектор . -ой частицы; и (г) - потенциальная энергия,
зависящая от взаимного расположения всех частиц.
Начальные условия задаются набором координат и скоростей частиц для момента
г = о
{X.,У,,2,,УХ1 V х . = 1,N},
где N - число частиц; X,У,Z- координаты частиц; ¥х,¥У,¥г - проекции скоростей частиц на оси координат.
ВРЕМЕННОЙ ШАГ ИНТЕГРИРОВАНИЯ
При расчетах методом молекулярной динамики временной шаг интегрирования выбирается как компромисс между двумя противоречивыми требованиями:
- с одной стороны шаг интегрирования должен быть достаточно мал, чтобы обеспечить устойчивость модели и необходимую точность расчетов;
- с другой стороны шаг интегрирования должен быть достаточно большим, чтобы время, необходимое для расчета задачи, было минимальным и позволило рассчитать задачу за физически приемлемое время.
Задача моделирования лазерного воздействием на металлы характерна тем, что некоторые атомы в результате воздействия мощного лазерного импульса получают настолько большую энергию, что временной шаг интегрирования для ансамбля атомов, рассчитываемого ММД и необходимый, чтобы обеспечить устойчивость модели, становится очень мал. Поэтому для расчета достаточно небольшой модели требуется значительное машинное время.
Определение оптимального шага интегрирования для кристаллита при 300 К.
Рассмотрим ансамбль атомов при температуре 300 К. Расчетную область представим в виде куба размером из 864 атомов с периодическими граничными условиями по всем координатным осям. Потенциал взаимодействия для алюминия рассчитан метом погруженного атома [4]. Расчеты выполняются с использованием программы LAMMPS [5].
На рис. 1 и 2 представлены графики проекции траектории движения контрольного атома по оси Х для шагов интегрирования 2^10"14 с и 1 •Ю"14 с.
Коор. х
в А 2,4-
1.8 --------------■-------------1-------------<-------------■-------------»_►
0 20 40 60 80 100 Шагов интегр.
Рис. 1. Проекция траектории движения контрольного атома по оси Х для At=2•10-14 с
Коор. X в А
2,4 -
1,7 ---------------1--------------1---------------1---------------!_►
О 50 100 150 200 Шагов интегр.
Рис. 2. Проекция траектории движения контрольного атома по оси Х для At=1•10-14 с
Ошибка в траектории движения контрольного атома за время интегрирования 2 пс для шагов интегрирования 2^10"14 с и Г10-14 с составляет 0,03437 А, что достаточно много для скоростей движения атомов при температуре 300 К. При такой ошибке интегрирования модель становится не устойчивой.
На рис. 3 и 4 представлены графики проекции траектории движения контрольного атома по оси Х для шагов интегрирования 1 •Ш-14 с и 1 •Ш-15 с.
Коор. х
в А
Коор. X в А
О 200 400 600 800 1000 Шагов интегр.
Рис. 4. Проекция траектории движения контрольного атома оси Х для At=1•10-15 с
Ошибка в траектории движения контрольного атома за время интегрирования 1 пс для шагов интегрирования Г10-14 с и Г10-15 с составляет 0,03067 А, что также много для скоростей движения атомов при температуре 300 К. При такой ошибке интегрирования модель становится не устойчивой.
Все расчеты для кристаллита с температурой 300 К приведены в табл.1.
Таблица 1
Погрешность траектории движения контрольного атома в проекции по оси Х для различных шагов интегрирования при температуре 300 К
Размер шага интегрирования At, с Количество шагов интегрирования Координаты атома по оси Х, А Погрешность положения атома, А
2-10-14 100 1,93985 1,97422 - 1,93985 = 0,03437
110-14 200 1,97422
1-10-14 100 2,05091 2,02024 - 2,05091 = 0,03067
110-15 1000 2,02024
110-15 100 2,06981 2,06978 - 2,06981 = 0,00003
1-10-16 1000 2,06978
Из табл. 1 видно, что для ансамбля атомов с температурой 300 К оптимальный размер временного шага составляет около Г10-15 с, так как дальнейшее уменьшение шага интегрирования, даже в 10 раз, не вызывает появления погрешности в четвертом знаке после запятой. Практические расчеты показывают, что при шаге интегрирования Г10-15 с, для ансамбля атомов при температуре 300 К, модель устойчива и расчетные параметры модели хорошо совпадают с экспериментальными данными [6 - 8].
Определение оптимального шага интегрирования для кристаллита при 3000 К, 10000 К, 20000 К и 40000 К.
Рассмотрим определение оптимального размера шага интегрирования для такого же, что и в предыдущих случаях, кристаллита. Расчеты по определению оптимального шага интегрирования представлены в табл. 2 - 5.
Таблица 2
Погрешность траектории движения контрольного атома в проекции по оси Х для различных шагов интегрирования при температуре 3000 К
Размер шага интегрирования At, с Количество шагов интегрирования Координаты атома по оси Х, А Погрешность положения атома, А
1-10-14 100 8,76825 3,21855
1-10-15 1000 11,9868
1-10-15 100 -11,41193 0,00003
1-10-16 1000 -11,41196
1-10-16 100 10,76782 0,00002
1-10-17 1000 10,76784
Из табл. 2 видно, что для ансамбля атомов с температурой 3000 К оптимальный размер временного шага составляет ~10-15 с, так как дальнейшее уменьшение шага интегрирования в 10 раз не вызывает появления погрешности в четвертом знаке после запятой.
Таблица 3
Погрешность траектории движения контрольного атома в проекции по оси Х для различных шагов интегрирования при температуре 10000 К
Размер шага интегрирования At, с Количество шагов интегрирования Координаты атома по оси Х, А Погрешность положения атома, А
5-10-15 100 7,49896 1,22628
1-10-15 500 6,27268
1-10-15 100 6,85705 0,00068
1-10-16 1000 6,85637
1-10-16 100 5,935011 0,000003
1-10-17 1000 5,935014
Таблица 4
Погрешность траектории движения контрольного атома в проекции по оси Х для различных шагов интегрирования при температуре 20000 К
Размер шага интегрирования с Количество шагов интегрирования Координаты атома по оси Х, А Погрешность положения атома, А
5-10-15 100 0,17131 7,587171
1-10-15 500 -6,87004
1-10-15 100 -7,96812 0,00177
1-10-16 1000 -7,96989
1-10-16 100 -6,919271 0,000002
1-10-17 1000 -6,919273
Таблица 5
Погрешность траектории движения контрольного атома в проекции по оси Х для различных шагов интегрирования при температуре 40000 К
Размер шага интегрирования с Количество шагов интегрирования Координаты атома по оси Х, А Погрешность положения атома, А
510-15 100 1,04015 6,89865
110-15 500 -5,8585
110-15 100 0,788149 0,036163
110-16 1000 0,751986
110-16 100 -0,715613 0,000002
110-17 1000 -0,715611
Из табл. 3-5 видно, что для ансамбля атомов с температурой 10000 К, 20000 К и 40000 К оптимальный размер временного шага составляет ~10-16 с, так как дальнейшее уменьшение шага интегрирования в 10 раз не вызывает появления погрешности в четвертом знаке после запятой.
Практические расчеты показывают, что при выбранных значениях шагов интегрирования, когда ошибка при уменьшении шага интегрирования не вызывает появления погрешности в четвертом знаке после запятой, модель устойчива.
В табл. 6 сведены данные, которые показывают, что шаг интегрирования связан с максимальной скоростью движения атома в ансамбле (или температурой кристаллита) таким образом, что расстояние, проходимое атомами, имеющими максимальную скорость, не должно превышать вполне определенного значения. В нашем случае это расстояние составляет порядка (0,01 - 0,05) А.
Таблица 6
Максимальная скорость атома в ансамбле в зависимости от температуры кристаллита
Температура кристаллита, К Максимальная скорость атомов, А /пс Шаг интегрирования, фс Расстояние за 1 шаг интегрирования, А
300 14,28 1 0,01428
3000 49,467 1 0,049467
10000 97,71 0,1 0,009771
20000 136,1469 0,1 0,013614
40000 188,33 0,1 0,018833
ЗАКЛЮЧЕНИЕ
На основании данных, сведенных в табл. 6, можно сделать вывод, что модель ансамбля атомов является устойчивой, если расстояние, которое проходит самый быстрый атом за 1 шаг интегрирования, не превышает 0,01 А.
Так как при моделировании лазерного воздействия на материалы скорость атомов в ансамбле будет изменяться достаточно быстро и в большом диапазоне, то целесообразно использовать механизм автоматического подбора временного шага основанный на том, чтобы самый быстрый атом в ансамбле за один шаг интегрирования не мог пройти расстояние большее установленной величины 0, 01 А.
СПИСОК ЛИТЕРАТУРЫ
1. Рыкалин Н.Н, Углов А.А., Кокора А.Н. Лазерная обработка материалов. М. : Машиностроение, 1975. 296 с.
2. Анисимов С.И., Имас Я.А., Романов Г.С., Ходыко Ю.В. Действие излучения большой мощности на металлы. М. : Наука, 1970. 367 с.
3. Рэди Дж. Действие мощного лазерного излучения. М. : Мир, 1974. 470 с.
4. Cai J. and Ye Y.Y. Simple analytical embedded-atom-potential model including a long-range force for fcc metals and their alloys // Phys. Rev. B 54. 1996. P. 8398-8410.
5. Pieter in 't Veld LAMMPS Users Manual. URL: http://lammps.sandia.gov (дата обращения 15.05.1909).
6. Chudinov V.G., Andreev V.V. Atomic Mechanism of primary defect heterogeneous nucleation under radiation in FCC metals with a large stacking energy // Nucl. Mater. 1990. V.172. P. 106-113.
7. Chudinov V.G., Andreev V.V. Atomic Mechanism of Solid-Liquid Transition in F.C.C. Metals with a Large Stacking-Fault Energy // Phys. stat. sol. A. 118. 1990. C. 415.
DEFINITION OF THE TIME STEP OF INTEGRATION AT MODEL OPERATION OF LASER ACTION ON METALS THE METHOD OF THE MOLECULAR DYNAMICS
Besogonov V.V., Andreev V.Ven., Andreev V.Vjach.
Institute of Applied Mechanics Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. Expedients of an estimate and self-acting calculation of the optimum size of a time step of integration surveyed at the solution of problems by a method of the molecular dynamics for ensembles of atoms with very major gamut of velocities the reference for processes of laser processing materials.
KEYWORDS: the molecular dynamics, step of integration, laser processing.
Бесогонов Валерий Валентинович, кандидат технических наук, старший научный сотрудник ИПМ УрО РАН, тел. (3412) 20-74-33, e-mail: besog@udman.ru
Андреев Вячеслав Вениаминович, кандидат технических наук, старший научный сотрудник ИПМ УрО РАН, Андреев Вячеслав Вячеславович, аспирант ИПМ УрО РАН