Монте-Карло моделирование активности углерода в ОЦК-железе

Сплавы на основе железа, несмотря на рост производства пластмасс, керамики и композитов, продолжают оставаться основными конструкционными материалами. Это связано с тем, что с помощью относительно несложной термомеханической обработки могут быть получены метастабильная ОЦК-фаза, перенасыщенная атомами С, обладающая высокой прочностью, пластичностью и рядом важных для практики физических свойств (электро- и теплопроводность, жаростойкость, магнитная проницаемость и т. п.). Одним из основных факторов, определяющих эти свойства, является параметр взаимодействия между атомами углерода. К сожалению, энергии взаимодействия атомов растворенного вещества не могут быть непосредственно измерены в эксперименте. Основным источником информации о взаимодействии были термодинамические данные о растворах углерода в железе [1].

Оценки энергии взаимодействия атомов углерода в ОЦК-железе были получены в [2–6].

Подробный обзор взаимодействий между атомами углерода в ОЦК-железе приведен в [6]. Из-за низкой растворимости углерода в α-железе (менее 0,02 ат. %), на практике у нас будет почти бесконечно разбавленный раствор, где взаимодействия C–C не появляются. Прямой расчет зависимости энергии взаимодействия С–С от активности углерода не является продуктивным. Авторы первых работ [7, 8], где были сделаны такие расчеты, получили отрицательную энергию взаимодействия С–С. Однако они использовали неправильное координационное число для междоузельных подрешеток [9]. Позднее Бхадешиа сделал более точные расчеты [10]. Он определил энергию взаимодействия ближайших соседей углерода в железе (ECC = 1,54 эВ), используя квантовохимическую теорию активности. Это значение согласуется с энергией ECC = 1,78 эВ, полученной в [6]. В работах [2–5] получены очень различные значения энергии С–С взаимодействия в ОЦК-железе с использованием методов первопринципного компьютерного моделирования и псевдопотенциального прибли- жения для решетки железа. К сожалению, до сих пор неясно, какие данные верны.

В настоящей работе мы проанализировали текущую ситуацию и решили проблему, вычислив активность углерода в ОЦК-железе. Для расчётов использовался метод Монте-Карло с доступным набором энергий С–С взаимодействия. Для этого нам необходимо было ответить на следующие вопросы:

• 1.    Почему данные по активности углерода воспроизводятся с различными наборами параметров С–С взаимодействия?

• 2.    Какое из предложенных значений параметров С–С взаимодействия следует рассматривать как наиболее подходящее для расчета термодинамических свойств ОЦК-железа?

Монте-Карло метод

Для моделирования активности углерода в железе мы использовали метод, предложенный в [11] и примененный для расчета углеродной активности в железе [12]. Расчет углеродной активности в железе проводили по следующей формуле:

a С = a conf exp ( A GjkT ) , (1) где a С – активность углерода в железе; a conf – конфигурационная часть активности и неконфигурационная часть. Для вычисления неконфигурационной части использовались данные работы [13].

Идея метода заключается в следующем. Химический потенциал системы представляется через отношение статистических сумм для двух значений числа частиц – N и N + 1. Удается свести данное отношение к усредненному по ансамблю ^e вАЛ' ^, где в = 1/kT, k – константа Больцмана, T – абсолютная температура, ΔEi – изменение энергии системы из N частиц при введении виртуального атома углерода в узел i при фиксированной позиции других атомов. Усреднение выполняется по обычной процедуре Монте-Карло, в результате чего приходим к следующему выра- жению для конфигурационной части активности [14, 15]:

a conf = M ( N C + 1 )

/ Г M Г m Vi z Z exp ( / I j = 1 L i = 1

A EilkT ) ^ . (2)

Суммирование по i выполняется по всем пустым октаэдрическим междоузлиям m =3 N Fe – N C и по числу шагов j ( j =1, 2, …, M ), используемых в расчете методом Монте-Карло, N _Fe и N _C – число атомов железа и углерода в решетке. Для моделирования использовали систему 20×20×20 элементарных ячеек ОЦК-решётки железа. Сходимость результата усреднения достигали при числе шагов Монте-Карло порядка 500.

Вычисление активности углеродав ОЦК-железе

Как упоминалось выше, энергии C–C взаимодействия были получены с использованием первопринципного моделирования [3–6] (см. таблицу). Авторы [3–5] использовали су-перячейку из 128 атомов ОЦК-железа. В наших расчетах мы рассматривали суперячейку из 54 атомов (рис. 1). Расчеты проводили для широкого диапазона температур ( Т = 955, 975, 1000, 1026, 1056, 1070, 1086, 1121 К). Однако ниже мы представили только данные для температуры 1000 К, поскольку не наблюдалось каких-либо качественных различий в результатах. Растворимость углерода в ОЦК-железе составляет 0,1 ат. %. Мы вычислили зависимость активности углерода в ОЦК-железе от концентрации, используя данные из таблицы, и сравнили полученные кривые (рис. 2).

Все рассмотренные исследования отлично согласуются с экспериментом [16] (рис. 2а), поскольку растворимость углерода в ОЦК-железе мала. Так как атомы С расположены на достаточном расстоянии друг от друга, энергия взаимодействия С–С мала и конфигурационная часть активности в уравнении (1) близка к 1. Чтобы сравнить различные параметры взаимодействия C–C, мы экстраполи-

Энергии взаимодействия атомов углерода в ОЦК-железе, эВ

Координационная сфера	1	2	3	4a	4b	5	6	7a	7b	8a	11
Домэйн [3]	–	0,65	0,09	0,09	1,67	–0,13	–	–	–0,14	–	–0,16
Бароух [4]	1,91	0,82	0,23	0,02	1,52	0,06	0,07	0,02	0,11	0,02	–
Рубан [5]	1,95	0,96	0,34	0,25	1,9	–0,07	–0,02	–0,03	0,1	0,05	0,02
Данная работа [6]	1,78	0,66	0,12	–0,13	2,89	–0,14	–0,1	–0,08	0	0,03	–0,19

Физическая химия и физика металлургических систем

Рис. 1. Возможные положения атомов углерода в ОЦК-железе. Черные круги обозначают атомы железа; квадраты – атомы углерода. Начальный атом находится в положении 0

Рис. 2. Зависимость активности углерода в ОЦК-железе от концентрации. Кривые были получены методом Монте-Карло с использованием данных, приведенных в таблице

ровали кривую активности углерода в более широкий диапазон концентраций (не достижимых в реальном ОЦК-железе). После этого мы видим, что рассматриваемые исследования приводят к различным кривым активности (см. рис. 1). Кривые, рассчитанные из [3] и наших данных [6], почти совпадают, а другие две кривые (из [4] и [5]) лежат выше. Согласие наших результатов с Рубаном [3] можно объяснить тем, что он не учитывал все возможные расположения атомов углерода. Мы моделировали суперячейку из 54 атомов, в то время как в [4, 5] число атомов в супер-ячейке было в два раза больше. Это могло бы означать, что на энергии взаимодействия атомов углерода, расположенных в 5–11 координационных оболочках, влияют периодические граничные условия. Мы хотели проверить эту гипотезу и выяснить необходимость рассмотрения таких далеких пар C–C для расчета углеродной активности. Мы провели численные эксперименты с учетом взаимодействия атомов углерода, расположенных только в первых четырех координатных оболочках. Результаты подтвердили, что, не учитывая неточных значений энергий взаимодействия С–С, мы можем получить кривую активности, которая совпадает с одной, рассчитанной из более точных исследований [5] (рис. 2б). Но, что более интересно, мы обнаружили, что взаимодействия атомов углерода, расположенные в 5–11 координационных оболочках, вообще не важны. Действительно, кривые активности, рассчитанные с учетом C–C взаимодействий как в 11, так и только в четырех координатных оболочках, почти совпадают. Более того, можно отметить, что две различные параметризации [4, 5] энергий взаимодействия С–С дали нам одни и те же кривые активности, когда атомы углерода были расположены в первых четырех координатных оболочках. Поэтому анализ кривых активно- сти не позволяет получить однозначную информацию о взаимодействии С–С в ОЦК-железе.

Заключение

Мы провели обзор литературы по взаимодействию С–С и методам расчета углеродной активности в ОЦК-железе. Мы сравнили рассматриваемые данные с использованием Монте-Карло-моделирования углеродной активности в ОЦК-железе. Для ОЦК-железа было установлено, что для вычисления кривой активности углерода достаточно установить энергии взаимодействия для атомов углерода в первых четырех координационных сферах. Мы показали, что экспериментальная кривая активности может быть воспроизведена с использованием нескольких наборов параметров С–С взаимодействия. Поэтому проблема описания взаимодействия между атомами углерода математически некорректна.

Исследование выполнено за счет гранта Российского научного фонда №16-19-10252.