A multiscale model of creep in steels with account for the microstructure

Full Text

ВВЕДЕНИЕ

Цель настоящей работы — разработка физической модели для расчетов скорости радиационной ползучести в сталях. Рассматривается движение дислокаций в неоднородном профиле локальных деформаций и напряжений и взаимодействие дислокаций с образовавшимися под действием облучения точечными дефектами и элементами микроструктуры. Меняющаяся со временем микроструктура и напряжения, создаваемые системой протяженных радиационных дефектов, могут существенно менять приложенное извне напряжение и вместе с этим динамику перемещения дислокаций. Поэтому учет взаимодействия дислокаций с кластерами дефектов и дислокаций необходимо проводить в рамках одного последовательного атомистического подхода.

Обычные феноменологические оценки механизма ползучести в литературе практически не учитывают сами механизмы преодоления барьеров скольжения, создаваемые различными препятствиями, за счет потока точечных дефектов. Предполагается, что скорость переползания дислокаций однозначно связана с эффективностью стока на них и никак не зависит от типа препятствия. Между тем в процессе перестройки микроструктуры, связанной с зарождением и ростом пор и петель и изменением локальных упругих полей, вклад препятствий каждого типа в процессе переползания меняется. А если учесть, что радиационная ползучесть в сталях может зависеть и от других факторов (движение зерен, реакции в частичных дислокациях и т.п.), становится ясно, что улучшение предсказательной силы расчетов ползучести невозможно без понимания физики взаимодействия дислокаций с препятствиями различного типа.

Расчет микроструктуры и ее влияние на деформацию начинается с расчетов концентрации дефектов различного типа, образовавшихся под облучением. Далее учитывается их диффузия, образование кластеров, массивов дислокаций и дислокационных петель. В гранецентрированных кубических металлах (ГЦК) это в основном тетраэдры дефектов упаковки (ТДУ) (stacking fault tetrahedra (SFT)), поры и петли. А в объемно-центрированных кубических металлах (ОЦК) — это также петли и нанопоры. В процессе облучения эти дефекты создают локальные одноосные, гидростатические и сдвиговые напряжения γ ^р, препятствуя движению дислокаций, приводя к коллективному механическому отклику материала, в частности, упрочнению при облучении. Изменения микроструктуры влияют на следующий по масштабу уровень локализации деформации и напряжений, обычно выявляемый расчетами в рамках метода конечных элементов (МКЭ).

Многомасштабная процедура определения скорости ползучести начинается с молекулярно-динамического (МД) моделирования предельных напряжений преодоления дислокацией кластеров и протяженных дефектов. Эти и другие параметры взаимодействия отдельной дислокации с дефектами различного типа использовали далее для определения скорости скольжения дислокации. Затем решали уравнение кинетики системы взаимодействующих дислокаций в профиле дефектов каждого типа с заданной концентрацией с помощью программы дискретной динамики дислокаций (ДДД, Discrete Dislocation Dynamics, DDD). Решение данного уравнения определяло скорость ползучести по механизму скольжения, вызванного переползанием (Climb-Induced Glide (CIG)).

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

Результаты предыдущих исследований, описанные, например, в [2–4], не привели к прояснению в полной мере затронутых в настоящей работе вопросов.

МЕТОДИКА МОДЕЛИРОВАНИЯ

Алгоритм построения дислокаций в кристаллах с ОЦК- и ГЦК-структурой был подробно описан в работе [5]. Сначала проводили минимизацию полной энергии системы, состоящей из атомов одного сорта, поскольку достижение локального минимума, соответствующего наличию квадруполя дислокаций, происходит, как правило, быстрее, чем для сплава. Визуализацию дефектов и анализ дислокационных реакций проводили каждые 100 МД временных шагов. При этом расположение атомов ядра дислокации и включений определяли из анализа параметров отклонения центра симметрии (c_sym), как описано в [5]. Атомы с координацией укладки, отличной от основной структуры, также визуализировали параметром c_sym, доступным в LAMMPS.

Для моделирования взаимодействия краевой дислокации с небольшими кластерами точечных дефектов использовали точные потенциалы взаимодействия атомов легирующих элементов, имеющиеся в базе “Interatomic Potentials Repository” [6]. Для сплава FeCr применяли потенциал, описанный в [7], а для FeCrNi — в [8]. Многочисленные расчеты такого типа, известные из литературы, как правило, относятся либо к чистым металлам, либо к грубым аппроксимациям, использующим методы континуальной механики. LAMMPS обеспечивает последовательную процедуру релаксации МД ячейки с дислокациями, которая возможна только с применением квадруполя дислокаций, где можно наложить периодические граничные условия.

Объем моделирования для ОЦК решетки в сплаве FeCr составлял 0.7×22.7×24.7 нм и содержал около 33 000 атомов, а для ГЦК-решетки в сплаве FeCrNi — 0.88×22.2×22.4 нм и около 37 000 атомов. Использованы периодические граничные условия.

Другое использованное нами преимущество LAMMPS — применение различных типов термодинамических ансамблей, без которых какая-либо последовательная интерпретация экспериментальных данных по эволюции и устойчивости микроструктур попросту невозможна.

Релаксацию размеров кристаллита с дислокациями при заданном внешнем тензоре напряжений рассчитывали с помощью изострессового-изоэнтальпического ансамбля минимизацией энергии взаимодействия атомов в ОЦК/ГЦК-решетке, содержащей в узлах необходимое количество атомов Fe и случайно распределенных атомов Cr и Ni заданной концентрации. В большинстве наших расчетов применяли гибридный потенциал eam/alloy eam/fs; alloy FeCr_d.eam.alloy [7]. Для чистого железа применяли потенциал Мишина (Fe_Mishin2006.eam.alloy) [9], а для системы FeCrNi — потенциал Fe–Ni–Cr_fcc.eam.fs [8].

На рис. 1 с помощью параметра c_sym (расчет фактора отклонения от центра симметрии) проведено сравнение упругих полей и формы ядра дислокации (ЯД) для чистого железа и сплава Fe12%Сr. Это сравнение показывает, что наличие второй компоненты сплава делает ЯД менее прямолинейным, более размытым, даже без температуры.

 

Рис. 1. Поле деформаций вокруг ядер дислокаций в чистом Fe (а) и в сплаве Fe–12%Cr (б).

 

Понятно, что такой эффект существенно сказывается на скорости скольжения дислокации в кристалле. И речь идет не просто об изменении скорости дислокации в бездефектном кристалле, которая практически не сказывается на скорости ползучести. Речь идет о сильном различии значений напряжений отрыва из-за искажений линии ядра по сравнению с чистыми металлами, что приводит к качественно иной картине движения дислокаций. Рассчитанные напряжения являются базовыми компонентами в расчетах ДДД, результаты которых свидетельствует о принципиальном отличии движения дислокаций в сталях от того, что наблюдается в чистом железе.

В расчете скорости переползания дислокаций ключевым препятствием для всех работ по ползучести является эффект химического окружения на подвижность ТД в сталях. Нами был развит новый метод расчета сил стока ТД в поле дислокации с приложенным внешним напряжением, использующий ОК МК в интерактивном комплексе программ LAMMPS. Эта схема моделирования точно учитывала локальные поля при определении энергии седловой точки (ЭСТ) без привлечения приближений теории упругости, таких как метод тензора поляризации, что позволило последовательно усложнять модель на случай раздельной диффузии компонентов на стоки.

Для включения LAMMPS в схему симуляции нами был разработан специальный скрипт, являющийся суперпозиций программ, написанных на языке FORTRAN и C++, который позволял отслеживать энергии миграции в процессе моделирования отдельной траектории. Другая независимая программы отслеживала траектории отдельных ТД. Для ускорения расчета эффективности стока была предложена модель “серого” стока с радиусом поглощения/отражения $R_{\text{g}}<r_{\text{s}}$, в которой существенно сокращена диффузионная длина в области, представляющей наибольший интерес для нашей задачи. Модель серого стока проиллюстрирована на рис. 2.

Влияние упругих полей дислокации и внешнего напряжения учитывалось только внутри окружности этого радиуса, что существенно снижало время моделирования.

Симуляцию траекторий производили в этом небольшом объеме для расчета коэффициента поглощения/отражения и числа прыжков до поглощения или выхода из серой зоны. Затем эти результаты объединяли с независимыми универсальными расчетами числа прыжков в нейтральной области (независимых от концентрации компонент сплава и температуры).

Модель серого стока, показанная на рис. 2, включает в себя две области, приходящиеся на одну дислокацию и отличающиеся типом прыжков. Первая область имеет форму квадрата размером $L$, а внутри ограничена окружностью радиусом серого стока $R_{\text{g}}$. Вторая область внутри круга радиуса $R_{\text{g}}$ содержала черный сток — ядро дислокации (ЯД). Полный объем всей ячейки длиной $L_{\text{z}}$ в направлении дислокации составлял $L^2{L}_{\text{z}}$.

 

Рис. 2. Модель серого стока.

 

Процедура подсчета прыжков на траектории ТД типа 1 начинаются из распределенной равномерно случайной точки в нейтральной области и заканчивается после пересечения окружности радиусом $R_{\text{g}}$. Таким образом, процесс поглощения ТД на радиусе $R_{\text{g}}$ рассматривали как нейтральный сток. Траектории 2 и 4 включали поглощение на ядре дислокации. Первые начинаются на окружности $R_{\text{g}}$, вторые — из распределенной равномерно случайной точки внутри круга серой области. Траектории 3а и 3в являются циклическими, без поглощения, они начинающиеся и заканчивающая на окружности $R_{\text{g}}$. И наконец, траектории типа 4а и 4в начинаются внутри круга с поглощением на дислокации или выходом из него. Их доля определяется отношением площадей $q=\pi R_{\text{g}}^2/L^2$. Траектории, достигающие внешней границы ячейки $L\times L$, испытывали отражение. Траектории каждого типа $i\in 1÷4a\mathrm{,4}b$ характеризуется средним числом прыжков $k_i$ до момента захвата, входа или выхода из “серой” зоны. Таким образом, вероятность отражения серым стоком определяется из значений среднего числа траекторий каждого типа.

Моделирование траекторий в наиболее затратной по времени серой зоне давало значения коэффициента поглощения и среднего числа прыжков для каждого типа траекторий. Начальные положения ТД задавали с равномерным распределением в объеме моделирования и прослеживали до момента достижения границы серого стока. Траектории другого типа, где ТД перемещались вне серой зоны, начинались с окружности $R_{\text{g}}$ и заканчивались там же.

Во всех расчетах использовали значение радиуса серой области $R_{\text{g}}=17.5d\approx 50\text{\hspace{0.17em}нм.}$

Методом МК рассчитывалось количество прыжков из случайной точки нейтральной области до отражения серым стоком. Полученное количество прыжков давало отношение эффективностей стока $\xi \left(\sigma \right)/\xi \left(0\right)$. Таким образом, скорость ползучести становилась универсальной функцией напряжения $\sigma$ с коэффициентом, зависящим от концентрации дислокаций (для данного значения $R_{\text{g}}$). Пример расчета эффективностей стока при Т=800 К при внешнем напряжении σ=200 МПа и плотности дислокаций r=3·10¹⁴ м^-2 приведен в табл. 1.

 

Таблица 1. Эффективности стока для чистого Fe и сплавов FeCr и FeCrNi

Композиция	Тип ТД	ξ(0)	ξ(σ)	1 − ξ(σ) / ξ(0)
Fe	Вак	0.0102	0.0102	0
	СМА	0.0102	0.01	0.006
Fe88Cr12	Вак	0.012	0.0114	0.014
	СМА	0.0078	0.0072	0.030
Fe60Cr30Ni10	Вак	0.0234	0.0216	0.012
	СМА	0.009	0.008	0.032

 

Наиболее простая оценка скорости ползучести с участием скольжения и переползания дислокаций в среде с препятствиями основана на модели непрозрачных препятствий. Такая оценка следует из комбинации набора высоты и скольжения между двумя препятствиями, которая приводит к известному уравнению Орована:

$\dot{\epsilon }=\rho V_{\text{c}}\frac{L}{h}b$. (1)

Здесь V_c — скорость переползания, L — расстояние между препятствиями, h — высота препятствия,  — плотность дислокаций,  — модуль вектора Бюргерса. В более детальной картине взаимодействия дислокаций с СТД различного типа определяющим параметром является напряжение преодоления дислокациями препятствий.

В приближении линейного натяжения теории дислокаций напряжения отрыва и прохождения дислокации через препятствие зависят от критического напряжения отрыва τ_c = α Gb/L, где G — модуль упругого сдвига, α — коэффициент усиления. В случае системы препятствий для дислокации данное выражение переходит в выражение:

${\tau }_{\text{c}}=\frac{\alpha }{\lambda }=\alpha \sqrt{Nd}$, (2)

где N — плотность препятствующих объектов, d – их диаметр,  — средняя длина до столкновения дислокации с препятствием.

Напряжение сдвига, вызванное системой пор с плотностью N_il и средним диаметром D_il, определяется в модели дисперсного барьерного упрочнения:

${\tau }_{il}=\alpha \mu b\sqrt{N_{il}{D}_{il}}$. (3)

Коэффициент усиления α, отражающий силу барьера, устанавливали из МД-моделирования, используя формулу (2). Для этого определяли напряжение срыва дислокации как функцию расстояния до центра задерживающего препятствия. В большинстве расчетов он был близок к значению α~0.43. Сравнение этой формулы для пор с формулой Бэкона–Кокса–Скаттергуда (BKS) [10]:

${{\tau }_{vo}}_0=\frac{\mu b}{2\pi L_{vo}}{\left[\ln \left(\frac{L_{vo}}{b}\right)\right]}^{-\frac{1}{2}}{\left[\ln \left(\frac{D_{vo}}{b}\right)+C\right]}^{\frac{3}{2}},$ (4)

где $L_{vo}=1/\sqrt{N_{vo}{D}_{vo}}$ есть среднее расстояние между порами, $D_{vo}^{\text{'}}=D_{vo}{L}_{vo}/(D_{vo}+L_{vo})$, C — постоянная, давало значение α~0.7.

МД-моделирование показало, что, как и в чистом металле, критическое напряжение сдвига для прохождения через пору уменьшается с повышением температуры. Если оценивать температурный эффект путем подгонки данных МД согласно модели [9]:

${\tau }_{vo}={\tau }_{vo0}\left(1-\exp \left(\frac{E_t}{kT}\right)\right)$, (5)

то для Fe12%Cr энергия термической активации получается равной 0.027 эВ. При этом средний размер нанопор в FeCrNi полагали равным D_vo=10 нм, а среднее расстояние между ними L_vo =1 мкм.

Для взаимодействия петли с краевой дислокацией среднее значения параметра взаимодействия по ориентациям петли соответствовало величине α = 0.07. Для взаимодействия петли с винтовыми дислокациями выбирали α = 0.6 независимо от того, являются ли векторы Бюргерса b изначально параллельными или наклонными, т.е. все петли считали сильными препятствиями для винтовой дислокации. По значению α определяли эффективную скорость скольжения дислокации, применяемую в качестве входного параметра в расчетах ДДД.

Модельный кристаллит в МД-расчетах релаксировали в рамках изоэнтальпического-изотермического ансамбля при температуре T =300 К с шагом MD по времени 1 фс в течение 100 пс. Затем прикладывали напряжение сдвига. Соответствующее приложенное разрешенное напряжение сдвига было рассчитано из соотношения τ = F_x/(L_xL_y), где F_x — сила сдвига, оказываемая атомами внутренней области на верхней границе. Краевая дислокация скользила по плоскостям (110) в сплаве Fe12%Cr со стационарной скоростью 29.2 нм/н. Моделирование без наличия петель показало, что при T = 300 К для этого требуется напряжение τ = 2.1 МПа. Моделирование считали завершенным, когда скользящая дислокация либо преодолевала барьер, создаваемый включением, либо достигала границы ячейки.

Для верификации алгоритма динамики дислокаций сначала было выполнено моделирование движения дислокаций в кристалле без преципитатов и других дефектов.

Напряжение Пайерлса дислокации в отсутствие кластера составляло ~5 МПа. Первоначально наблюдали притяжение дислокации к кластеру с падением напряжения сдвига почти до нуля, но дальнейшая деформация приводила к линейному увеличению напряжения до критического значения 15 МПа. Далее дислокация становилась изогнутой, освобождаясь от петли/встречной дислокации, так что и приложенное напряжение падало до напряжения Пайерлса. Напряжение отрыва (critical resolving shear stress, СRSS) t_с ~ 22 МПа фиксировали далее для определения скорости скольжения в ДДД.

Чтобы вызвать скольжение дислокации, к МД-ячейке прикладывали сдвиговое напряжение 0, 30, 50 и 80 МПа. МД-ячейка предварительно релаксировала с помощью NτT ансамбля при температуре T = 300 или 600 К. Установившаяся скорость v_D изолированной дислокации при этой скорости деформации оценивали из соотношения Орована.

Важный вопрос во всей схеме влияния микроструктуры на механические свойства — какие локальные напряжения действуют на дислокацию, заставляющие ее двигаться и преодолевать препятствия. Общее упрочнение, вызванное различными видами дефектов облучения, рассматривается на основе метода суперпозиции. Он также реализуется в методе ДДД. Движение дислокаций приводит к деформации протяженной области кристалла. В рамках линейной теории упругости на каждую дислокацию действует сила Пича–Кехлера. При облучении поле напряжений, испытываемое дислокацией, выражается как σ = σ₀ + σ_irr, где σ₀ — приложенное поле. Поле σ_irr, индуцированное неподвижными дефектами облучения, равно

${\sigma }_{irr}=\left(\begin{array}{c}{\tau }_{irr}(\mathbf{b}\otimes \mathbf{n}\mathbf{+}\mathbf{b}\otimes \mathbf{n})\mathrm{sign}\left({\mathrm{\tau }}_{\mathrm{g}}\right), \left|{\mathrm{\tau }}_{\mathrm{g}}\right|>{\mathrm{\tau }}_{\mathrm{irr}}\\ {\tau }_g(\mathbf{b}\otimes \mathbf{n}\mathbf{+}\mathbf{b}\otimes \mathbf{n})\mathrm{sign}\left({\mathrm{\tau }}_{\mathrm{g}}\right), \left|{\mathrm{\tau }}_{\mathrm{g}}\right|<{\tau }_{irr}\end{array}\right).$ (6)

Здесь b и n — вектор Бюргерса и вектор нормали к плоскости скольжения; ⊗ обозначает тензорное произведение; t_g пропорционально силе Пича–Кехлера при отсутствии дефектов.

В классической формулировке трехмерного (3D) ДДД криволинейные линии дислокации разбиваются на небольшие сегменты, и их эволюция отслеживается в плоскостях скольжения и в объеме 3D-моделирования. Для конкретных расчетов нами использован код RheoMan [11], находящийся в открытом доступе. Использована 2D-версия программы, в которой учитывается скольжение лишь в одной плоскости, а взаимодействие дислокаций друг с другом и с неподвижными размерными дефектами учитывается предельными напряжениями срыва. Предполагается, что дислокации представляют собой идеальные прямые отрезки, параллельные друг другу и перпендикулярные плоскости огранки. Двигаясь в своих плоскостях скольжения, они меняли локальный профиль напряжений. А по их смещениям вычисляли пластическую деформацию, произведенную за определенный временной шаг. Условные препятствия размером h были распределены равномерно с заданной плотностью. Алгоритм движения дислокаций учитывал преодоление препятствия в случае превышения локального напряжение в точке соприкосновения над предельным значением t_с.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

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

 

Рис. 3. Накопленная деформация как функция времени для механизмов скольжения и переползания скольжением.

 

На рис. 4 показаны начальная и окончательная сетки дислокаций в Fe12%Cr при наложении растягивающей нагрузки σ_yy. На рис. 5 изображены начальное (равномерное) и конечное распределение напряжений. Видно, как в результате эволюции перераспределение плотности дислокаций приводит к значительному перераспределению напряжений. Нижняя кривая на рис. 3 показывает составляющую поля деформаций, полученную при рассмотрении только скольжения, а верхняя кривая — и переползания и скольжения. В обоих случаях можно идентифицировать стенки дислокаций, перпендикулярные направлению скольжения, поскольку дислокации с одинаковым знаком имеют тенденцию к самоорганизации вдоль этого направления. Более того, наблюдается также образование диполей, пар дислокаций, характеризующихся одним и тем же вектором Бюргерса, но с противоположным знаком.

 

Рис. 4. Распределение дислокаций до облучения (а) и после облучения (б) в ячейке моделирования ДДД. Крестиками показано положение дислокаций.

 

Рис. 5. Распределение напряжений до облучения (а) и после облучения (б) в ячейке моделирования ДДД, соответствующих распределению дислокаций, показанному на рис. 4.

 

Пример расчета скорости радиационной ползучести для сталей ЭП823 и ЭК164 приведены в таблице 2. Экспериментальные данные, приведенные в табл. 2, получены авторами настоящей работы на экспериментальном оборудовании АО “ИРМ”. Для расчетов использованы средние значения концентрации СТД, близкие к известным из экспериментов.

 

Таблица 2. Данные по скорости радиационной ползучести в сталях

Давление, МПа	Скорость деформации, %/с, эксперимент	Скорость деформации, %/с, переползание	Скорость деформации, %/с, переползание и скольжение
ЭП823
0	0.196 10−8	0	0.1 10−9
50	0.588 10−8	0.4 10−8	0.45 10−8
80	1.176 10−8	0.6 10−8	0.7 10−8
ЭК164
55	0.500 10−8	0.2 10−8	0.23 10−8

 

Размер ДДД ячейки принимали равным экспериментальному значению размера зерна.

ЗАКЛЮЧЕНИЕ

1. Развитый выше многомасштабный метод моделирования эволюции микроструктуры выявил принципиальное отличие моделей движения дислокаций в сталях и чистых металлах, вызванное нерегулярной формой ядра дислокаций в концентрированных растворах.
2. Не имеющая аналогов схема моделирования движения дислокаций в оболочке LAMMPS дает для скорости ползучести в модельных сплавах FeCr и FeCrNi значения, близкие к экспериментальным.

Данная работа финансировалась за счет средств бюджета института. Никаких дополнительных грантов на проведение или руководство данным конкретным исследованием получено не было.

Авторы данной работы заявляют, что у них нет конфликта интересов.

About the authors

K. Y. Khromov

National Research Center "Kurchatov Institute"

Author for correspondence.
Email: khromov_ky@nrcki.ru
Russian Federation, Moscow, 123098

V. A. Ryabov

National Research Center "Kurchatov Institute"

Email: khromov_ky@nrcki.ru
Russian Federation, Moscow, 123098

A. V. Kozlov

JSC “Institute of Nuclear Materials”

Email: khromov_ky@nrcki.ru
Russian Federation, Zarechny, Sverdlovsk region, 624250

V. L. Panchenko

JSC “Institute of Nuclear Materials”

Email: khromov_ky@nrcki.ru
Russian Federation, Zarechny, Sverdlovsk region, 624250

References

Supplementary files