CC BY
48
УДК 538.913
DOI: 10.14529/met160206
МОЛЕКУЛЯРНО-ДИНАМИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ ДВУХОСНЫХ ДЕФОРМАЦИЙ НА РАСТВОРИМОСТЬ ВОДОРОДА В ОЦК-ЖЕЛЕЗЕ С ИСПОЛЬЗОВАНИЕМ ЕАМ-ПОТЕНЦИАЛОВ
Д.А. Емелин, А.А. Мирзоев
Южно-Уральский государственный университет, г. Челябинск
Комплекс негативных воздействий водорода на металл называют водородной деградацией. Процессы водородной деградации существенно зависят от особенностей диффузии и растворимости водорода в конкретных материалах. Многие промышленные изделия при изготовлении сохраняют существенные остаточные напряжения (сварные трубы, сварные швы, несущие балки). В этом случае особую ценность для прогнозирования процессов деградации имеют зависимости растворимости и коэффициента диффузии водорода от температуры образца и напряжений, приложенных к образцу. В этой работе мы приводим результаты тестирования потенциалов Картер на воспроизведение основных энергетических характеристик, хорошо изученных методами первопринципного моделирования, а именно: энергии растворения водорода и величины диффузионных барьеров. После этого приводятся результаты исследования зависимости энергии растворения атомов водорода от величины двухосной деформации. Особый интерес представляет вопрос о воспроизведении потенциалом Картер перескока водорода из тетраэдриче-ских пор в октаэдрические поры ОЦК-железа под действием двухосных деформаций. Результаты моделирования сопоставляются с аналогичными результатами, представленными в соответствующих статьях и хорошо согласуются с результатами расчета из первых принципов. Потенциал B воспроизводит переход водорода из тетрапор в октапоры под действием двухосных напряжений.
Ключевые слова: молекулярная динамика; энергия растворения водорода; ОЦК-железо.
Введение
Комплекс негативных воздействий водорода на металл называют водородной деградацией [1-2]. Процессы водородной деградации существенно зависят от особенностей диффузии и растворимости водорода в конкретных материалах. Современный обзор исследований влияния таких особенностей можно найти в [3-5].
Многие промышленные изделия при изготовлении сохраняют существенные остаточные напряжения (сварные трубы, сварные швы, несущие балки). В этом случае особую ценность для прогнозирования процессов деградации имеют зависимости растворимости и коэффициента диффузии водорода от температуры образца и напряжений, приложенных к образцу. Однако долгое время подобные исследования проводились отрывочно, и построение теории влияния напряжений на поведение водорода в металлах далеко от завершения.
Экспериментально подтверждено [5], что в разбавленных твердых растворах водорода в ОЦК-железе атомы водорода при комнатной температуре располагаются в тетраэдрических порах решетки железа. Однако в экспериментальных работах с принудительным насыщением образца водородом с помощью электрохимических методов или водородной атмосферы с повышенным давлением (см. обзор [6]) было показано, что часть атомов водорода оказывается в октаэдрических порах, причем их количество увеличивается с ростом температуры. Результаты ряда экспериментов (см. [7] и ли-
тературу, приведенную в источнике) показали, что в присутствии растягивающих напряжений эта доля может быть увеличена посредством локальной кластеризации водорода. Кластеры имеют форму плоскостей с периодом, совпадающим с параметром ОЦК-железа [7]. Этот эффект был обнаружен и в отсутствие внешних растягивающих напряжений. Его связывают с увеличением концентрации водорода. При повышении температуры происходит рост кластеров, а при температурах выше 500 К кластеры распадаются (см. [8, 9]). Недавние результаты моделирования зависимости энергии растворения водорода от одноосной, двухосной и трехосной деформаций с использованием первопринципных методов [10] показали возможность увеличения доли растворенного в октапорах водорода при приложении двухосных напряжений к суперячейке. Было показано, что при относительном двухосном сжатии образца более чем на 4 % и растяжении более чем на 6 % водород начинает переходить из тетрапор в окта-поры.
В связи с этим возникает вопрос о влиянии данного явления на механизмы диффузии водорода в ОЦК-железе, подвергнутом деформации. Поскольку экспериментальное изучение степени относительного заполнения окта- и тетрапор ОЦК-железа является достаточно сложным, то представляется разумным привлечь для рассмотрения данного вопроса методы прямого компьютерного моделирования. В литературе имеется множество данных [2-6] о коэффициенте диффузии водорода,
полученных как из эксперимента, так и методами компьютерного моделирования, которые отличаются на порядки. Для внесения дополнительной ясности в исследуемый вопрос необходимо использовать методики моделирования диффузионных процессов. Они требуют рассмотрения системы с большим числом атомов (порядка ~103). Поэтому подходящей методикой является метод молекулярной динамики. При его использовании встает вопрос о модели межчастичного взаимодействия. В настоящее время наиболее точными потенциалами для моделирования считаются потенциалы погруженного атома (ЕАМ-потенциалы), позволяющие качественно учитывать влияние окружения на межчастичные взаимодействия. Для молекулярно-динамического моделирования системы Fe-H в настоящее временя разработаны ЕАМ-потенциалы (Руда, Вен, Картер) и МЕАМ-потенциал Ли (см. ссылки в [11]). В работах Пак-стона и Чу-Чун Фу [12, 13] показано, что ЕАМ-потенциалы Картер [11] приводят к результатам, хорошо согласующимся с рядом данных натурного и компьютерного эксперимента. В частности, ЕАМ-потенциалы Картер хорошо воспроизводят диффузионные барьеры и энергию растворения водорода в чистом, бездефектном ОЦК-железе. В своей работе Картер представила четыре типа потенциалов. Однако из-за того, что тестирование потенциалов проводилось для ограниченного набора свойств системы Fe-H (зависимость диффузионных барьеров и энергии растворения водорода от упругих деформаций, энергия связи водорода и вакансий, энергия связи водорода и ядра винтовой дислокации), выделить наиболее точный из построенных потенциалов не удалось. Поэтому в данной работе мы использовали все 4 указанных потенциала (А, В, А', В') с целью определения оптимального для моделирования процесса диффузии водорода в ОЦК-железе.
В этой работе мы проводим тестирование потенциалов Картер на воспроизведение основных энергетических характеристик, хорошо изученных методами первопринципного моделирования, а именно: энергии растворения водорода и величины диффузионных барьеров. После этого проведем исследование зависимости энергии растворения атомов водорода от величины двухосной деформации. Результаты моделирования сопоставляются с аналогичными результатами, полученными в статьях [10] и [12, 13].
Методика моделирования
Для получения результатов использовалась реализация метода молекулярной динамики в программном пакете LAMMPS [14]. В работе тестируются 4 межатомных потенциала Картер. Воспроизведение энергетических параметров проводилось на системе с суперячейкой 10x10x10 атомов железа и одним атомом водорода. Минимиза-
ция потенциальной энергии и объема структур проводилась методом сопряженных градиентов с точностью 10-25 по энергиям и по силам. Энергия и сила измерялись в электронвольтах (эВ) и элек-тронвольтах/ангстрем (эВ/А). При расчете энергий растворения водорода в заданной конфигурации при фиксированном атоме водорода на нем обнулялась сила.
Перед проведением численных экспериментов один атом водорода помещался в суперячейку при 0 К, после чего проводилась минимизация полной энергии системы и действующих на атомы сил. Данная процедура показала, что предпочтительной позицией для водорода являются тетраэд-рические поры, что согласуется с наблюдаемой ранее картиной в экспериментах и при моделировании системы ОЦК-Fe-H [5].
Ключевую роль при моделировании динамических характеристик системы Fe-H играет воспроизведение диффузионных барьеров и энергий активации водорода. Для него возможны три направления перескока: октапора - тетрапора - окта-пора (O-T-O), тетрапора - октапора - тетрапора (T-O-T) и тетрапора - локальный минимум - тетрапора (T-S-T) (рис. 1). В работе [6] было показано, что влияние на коэффициент диффузии имеют только переходы T-S-T, T-O-T. Соответствующие им энергии активации равны изменению энтальпии процесса перескока водорода
AE = AH =AU + pAV,
где AU - изменение потенциальной энергии; p - давление; AV - изменение объема системы. В работе [15] при помощи ab initio методов были рассчитаны значения AU в выделенных направлениях [100]O (T-O-T) и [101]t (T-S-T) (рис. 1) как функция смещения вдоль этих направлений. Аналогичный, более свежий расчет был проведен в работе [13]. В настоящей статье с целью тестирования ЕАМ-потенциалов мы провели аналогичные расчеты методом МД-моделирования в пакете LAMMPS.
Для воспроизведения энергий активации и потенциального барьера при моделировании процесса диффузии достаточно рассчитать потенциальную энергию системы в точках S, T, O после релаксации. Моделирование системы проводилось при фиксированном положении водорода. Энергетические барьеры рассчитывались по формулам:
AUt_s_T = Erelax (FewH(S)) - Erelax (FewH(T));
AUt-o-t = Erelax (Fe„H(O)) - Erelax (Fe„H(T)) , где AUT-s-T и AUT-O-T - потенциальные барьеры перехода (T-S-T) и (T-O-T); Erelax (FewH(S)),
Erelax (FeWH(T)) , Erelax (FeWH(O)) - соответствующие энергии системы с водородом, находящимся в точках S, T, O, после минимизации энергии и объема. Вместе с тем была рассчитана энер-
Физическая химия и физика металлургических систем
Рис. 1. Особые точки потенциальной энергии взаимодействия Fe-H. Квадраты - октапоры, треугольники -тетрапоры, черные незакрашенные точки - локальный минимум между тетрапорами
Результаты воспроизведения энергетических характеристик системы Fe-H
Характеристика DFT Эксп. [16] DFT+MD [13] Потенциалы A, B Потенциалы A', B'
ESoi, эВ 0,23 [13] 0,30 - 0,2365 0,2472 0,2746 0,2514
AUT-S_T, эВ 0,042 [15] 0,044 [13] - 0,040 0,0443 0,0539 0,0432 0,0643
AUT_O_T, эВ 0,035 [13] - 0,049 0,0561 0,0491 0,0625 0,0741
AUT-S_T -AUT_O_T , эВ 0,009 [13] - 0,009 -0,0118 0,0048 -0,0193 -0,0098
гия растворения водорода (водород располагался в тетраэдрической поре)
Esol = E^е„Н) -) -0,5Е(^), где для энергии связи молекулы водорода E(H2) принято значение -4,651 эВ. Результаты расчета представлены в таблице.
Проведенное тестирование показывает, что только потенциал B хорошо воспроизводит разницу ДU между барьерами и величину самих барьеров. Остальные результаты моделирования ДUT-0-т и ДUT—8-т неплохо согласуются с результатами [13, 15], но приводят к неверным результатам для величины ДUT—8-т -ДUT-0-T . В связи с этим дальнейшее моделирование проводилось с использованием потенциала B как наиболее точного.
Влияние деформаций на энергию растворения водорода
Особую важность в воспроизведении перехода из тетрапоры в октапору имеет зависимость энергии растворения водорода от двухосной деформации. Энергия растворенного водорода в де-
формированном образце рассчитывалась по
формуле
Edissol = E — E
е _ -^е^е^ Н -^е^ед, '
где ЕеР-елгН - энергия системы в присутствии напряжений с водородом; Ее р. - энергия системы
в присутствии напряжений без водорода.
Для рассмотрения различающихся конфигураций водорода в матрице ОЦК-железа требуется провести классификацию октапор и тетрапор. Любая занятая водородом октапора вызывает максимальное искажение решетки в направлении малой диагонали октаэдра. Эти искажения могут быть направлены по осям х, у, z или [100], [010], [001]. По отношению к ним разделим октапоры на три типа: октапоры x-, у-, z-типа. В то же время заполненная водородом тетрапора вызывает искажения решетки по двум ребрам тетраэдра. Среди всех возможных искажений, вызываемых водородом, находящимся в тетрапоре, неэквивалентных всего три: ху, xz, yz. Соответственно, будем различать тетрапоры ху-, xz-, ук--типа.
[001]
-[100]
, [010]
[0011
[100]
• *[010]
[001]
[100]
♦•[010]
1Г (хг)
2 Г (ху)
10 (г)
2 О (х)
Рис. 2. Используемые конфигурации (черные точки - атомы железа, белые - атом водорода)
Проведем деформацию образца по направлениям [100] и [010] одновременно. Нагрузки были выбраны так, чтобы результирующее относительное удлинение (сжатие) образца находилось в области 0,9^1,1 параметра суперячейки. Из-за наличия симметрии в структуре чистого ОЦК-железа нам достаточно рассмотреть среди октапор окта-поры x- и г-типа, среди тетрапор - ху- и хг-типа (см. рис. 2)
Результаты расчета представлены на рис. 3 и 4. Значения получены после проведения минимизации полной энергии системы для атома водорода во всех направлениях и для атомов железа по направлению [001].
При отсутствии напряжения диффузия водорода в бездефектном ОЦК-железе при комнатной температуре протекает по тетрапорам [5]. При двухосном сжатии по направлениям [100] и [010] менее чем на 7 % энергия растворения водорода в тетрапорах растет быстрее, чем энергия растворения водорода в октапорах типа 10. В этой области деформаций при моделировании обнаруживается тенденция к переходу водорода в октапору типа 1О. Она проявляется в виде смещения равновесного положения водорода из тетрапоры типа 1Т в сторону октапоры 10 по прямой линии. При сжатии более 7 % разница в энергии растворения водорода между октапорами типа 10 и тетрапорами 1Т становится существенной. В этой области происходит переход водорода из тетрапоры 1Т в октапору 10. В области растяжений менее 8 % у водорода в тет-рапоре типа 2Т появляется смещение атома водорода в сторону октапоры типа 20. Растяжения более 8 % приводят к переходу водорода в тетра-поре 2Т в октапору 20.
Полученные нами результаты позволяют высказать следующее положение. Предположим, что при комнатной температуре в ОЦК-железе возможно локальное повышение концентрации водорода в порах хг-типа. Индуцированные водородом в тетрапорах хг-типа двухосные растяжения складываются в конструктивные локальные двухосные
деформации решетки по направлениям х, г. В этом случае положение хг-тетрапор смещается в сторону г-октапор. При критических локальных концентрациях водорода, вызывающего деформации хг-типа более 8 %, начнет проявляться механизм перехода водорода из хг-тетрапор в г-октапоры. После перехода из равновесного положения вблизи хг-тетрапоры в г-октапору водород вызывает максимальные искажения вдоль оси г. Это искажение конструктивно складывается с другими искажениями в направлении г, вызванными водородом в хг-тетрапорах. Такие конструктивные искажения активируют механизм перехода для ближайших атомов водорода, находящихся в хг-тетра-порах, в г-октапоры. Для каждого атома водорода, находящегося в октапоре г-типа, энергия в ближайших тетрапорах хг-, уг-типа выше, чем в октапоре. Следовательно, для водорода в г-октапорах увеличивается длина диффузионного перескока. Она активируется при определенных температурах. При температурах ниже этого барьера часть локального скопления атомов водорода запирается в г-октапорах. Такой переход водорода из хг-тетра-поры в г-октапоры аналогичен попаданию его в высокоэнергичные ловушки, создающие вокруг себя такие же ловушки для водорода. Этот процесс со временем может привести к локальному скоплению водорода в г-октапорах - его кластеризации. Кластеры будут иметь периодическую структуру с периодом, соответствующим периоду решетки ОЦК-железа.
Аналогичную характеристику кластеру дают в работе [7]. Образованные кластеры с конструктивным искажением решетки вдоль оси г могут привести к понижению пороговой деформации, необходимой для фазового перехода железа из ОЦК- в ГЦК-фазу. Влияние водорода на этот переход экспериментально подтверждено [3]. Особенно важно то, что кластеризация водорода при двухосных деформациях более 10 % может привести к образованию микротрещин и дислокаций в окрестности кластера, что понизит предел текучести.
Физическая химия и физика металлургических систем
Обсуждение результатов и выводы
Результаты проведенных численных экспериментов подтверждают, что разработанные EAM-потенциалы Картер хорошо воспроизводят основные энергетические характеристики: энергию растворения водорода, величину потенциальных барьеров, разницу между ними. В присутствии двухосных деформаций наблюдается хорошее согласие результатов моделирования зависимости энергии растворения с использованием потенциала B в области двухосного сжатия с результатами [10]. Также результаты проведенных численных экспериментов показали, что потенциал B пригоден для воспроизведения механизма перехода водорода из тетрапоры в октапору. Полученные результаты согласуются с результатами расчета из первых принципов [10], а также в некоторых экспериментальных исследованиях [5, 7-9]. Проявление данного эффекта особенно важно учитывать при моделировании процесса диффузии водорода, так как он существенно повлияет на зависимость коэффициента диффузии от двухосных деформаций. Возможно также, что указанный эффект может служить причиной водородной кластеризации и влиять на фазовый переход ОЦК-ГЦК в железе с растворенным в нем водородом.
Работа поддержана грантом Российского научного фонда (проект № 16-19-10252).
Литература/References
1. Арчаков Ю. И.Водородная коррозия стали. М.: Металлургия, 1985. 192 с. [Archakov Yu.I. Vodorodnaya korroziya stali (Hydrogen Corrosion of Steel). Moscow, Metallurgiya Publ., 1985. 192 p.]
2. Колачев Б.А. Водородная хрупкость металлов. М.: Металлургия, 1985. 218 с. [Kolachev B.A. Vodorodnaya khrupkost' metallov [Hydrogen Embrit-tlement of Metals]. Moscow, Metallurgiya Publ., 1985. 218 p.]
3. Gaseous Hydrogen Embrittlement of Materials in Energy Technologies. Vol. 1. The Problem, Its Characterisation and Effects on Particular Alloy Classes. Somerday B.P., Gangloff R.P. (Eds.). Philadelphia, Woodhead Publishing Ltd., 2012. 500 p. DOI: 10.1533/9780857095374
4. Gaseous Hydrogen Embrittlement of Materials in Energy Technologies. Vol. 2. Mechanisms, Modelling and Future Developments. Somerday B.P., Gangloff R.P. (Eds.). Philadelphia, Woodhead Publishing Ltd., 2012. 840 p. DOI: 10.1533/9780857093899
5. Fukai Y. The Metal-Hydrogen System. Berlin et al., Springer Verlag, 2005. 500 p.
6. Kiuchi K., McLellan R.B. The Solubility and Diffusivity of Hydrogen in Well-Annealed and Deformed Iron. Acta Metallurgica, 1983, vol. 31, no. 7, pp. 961-984. DOI: 10.1016/0001-6160(83)90192-X
7. Fujita F.E. The Role of Hydrogen in the Fracture of Iron and Steel. Trans. Japan Inst. Metals, 1976, vol. 17, pp. 232-238. DOI: 10.2320/matertrans1960.17.232
8. Панасюк В.В., Андрейкив А.Е., Харин В.С. Модель роста трещины в деформированных металлов под действием водорода. Физико-химическая механика материалов. 1987. Т. 23, № 2. С. 111-124. [Panasyuk V.V., Andreikiv A.E., Kha-rin V.S. A Model of Crack Growth in Deformed Metals under the Action of Hydrogen. Materials Science, 1987, vol. 23, no. 2, pp. 111-124. DOI: 10.1007/BF00718130]
9. Ткачев В.И. Механизм обратного влияния водорода на механические свойства сталей. Физико-химическая механика материалов. 1999. Т. 35, № 4. С. 29. [Tkachov V.I. Mechanism of Reversible Effect of Hydrogen on Mechanical Properties of Steel. Materials Science, 1999, vol. 35, no. 4, pp. 477-484. DOI: 10.1007/BF02365744]
10. Zhou H.B., Jin S., Zhang Y., Lu G.H., Liu F. Anisotropic Strain Enhanced Hydrogen Solubility in bcc Metals: The Independence on the Sign of Strain. Physical Review Letters, 2012, vol. 109, no. 13, 135502. DOI: 10.1103/PhysRevLett.109.135502
11. Ramasubramaniam A., Itakura M., Carter E.A. Interatomic Potentials for Hydrogen in a-Iron Based on Density Functional Theory. Physical Review B, 2009, vol. 79, no. 17, 174101. DOI: 10.1103/PhysRevB.79.174101
12. Paxton A.T., Elsasser C. Electronic Structure and Total Energy of Interstitial Hydrogen in Iron: Tight-Binding Models. Physical Review B, 2010, vol. 82, no. 23, 235125. DOI: 10.1103/PhysRevB.82.235125
13. Hayward E., Fu C.C. Interplay Between Hydrogen and Vacancies in a-Fe. Physical Review B, 2013, vol. 87, no. 17, 174103. DOI: 10.1103/PhysRevB.87.174103
14. Plimpton S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of Computational Physics, 1995, vol. 117, pp. 1-19. DOI: 10.1006/jcph.1995.1039
15. Jiang D.E., Carter E.A. Diffusion of Interstitial Hydrogen into and Through bcc Fe from First Principles. Physical Review B, 2004, vol. 70, no. 6, 064102. DOI: 10.1103/PhysRevB.70.064102
16. Hirth J.P. Effects of Hydrogen on the Properties of Iron and Steel. Metallurgical Transactions A, 1980, vol. 11, no. 6, pp. 861-890. DOI: 10.1007/BF02654700
Емелин Дмитрий Анатольевич, аспирант кафедры общей и теоретической физики, Южно-Уральский государственный университет, г. Челябинск; emelin_d_a@mail.ru.
Мирзоев Александр Аминулаевич, д-р физ.-мат. наук, профессор кафедры общей и теоретической физики, Южно-Уральский государственный университет, г. Челябинск; mirzoevaa@susu.ru.
Поступила в редакцию 4 апреля 2016 г.
DOI: 10.14529/met160206
MOLECULAR DYNAMIC MODELING OF THE BIAXIAL STRAIN EFFECT ON HYDROGEN SOLUBILITY IN BCC Fe USING EAM POTENTIALS
D.A. Emelin, emelin_d_a@mail.ru, A.A. Mirzoev, mirzoevaa@susu.ru South Ural State University, Chelyabinsk, Russian Federation
The complex of negative effects of hydrogen on metal is termed hydrogen degradation. The processes of hydrogen degradation considerably depend on the characteristics of hydrogen diffusion and solubility in a specific material. Many manufactured articles keep their residual stresses during production (e.g. welded pipes and seams, supporting girders). In this case dependences of solubility and diffusivity of hydrogen on temperature and applied stress are of particular value for the prediction of hydrogen degradation. This paper presents the results of MD tests performed with Carter EAM potentials for the reproduction of general energy characteristics, namely hydrogen solution energy and diffusion barrier values, that are well examined with ab initio methods. The research on the relation between the hydrogen solution energy and the biaxial strain was also carried out. Of great interest is the capability of Carter potential to reproduce the hydrogen jump from tetrahedral to octahedral interstices of bcc Fe under the influence of the biaxial strains. Simulation results are compared with previous ab initio calculations and are in good agreement. Carter's potential B is capable of reproducing hydrogen transition from tetrahedral to octahedral interstices under biaxial stress.
Keywords: molecular dynamics; hydrogen solution energy; bcc Fe.
Received 4 April 2016
ОБРАЗЕЦ ЦИТИРОВАНИЯ
Емелин, Д.А. Молекулярно-динамическое моделирование влияния двухосных деформаций на растворимость водорода в ОЦК-железе с использованием ЕАМ-потенциалов / Д.А. Емелин, А.А. Мирзоев // Вестник ЮУрГУ. Серия «Металлургия». - 2016. -Т. 16, № 2. - С. 40-45. DOI: 10.14529/met160206
FOR CITATION
Emelin D.A., Mirzoev A.A. Molecular Dynamic Modeling of the Biaxial Strain Effect on Hydrogen Solubility in BCC Fe Using EAM Potentials. Bulletin of the South Ural State University. Ser. Metallurgy, 2016, vol. 16, no. 2, pp. 40-45. (in Russ.) DOI: 10.14529/met160206
