Математическое моделирование фазовых превращений в сталях в рамках подхода диффузионной границы Текст научной статьи по специальности «Физика»

2012 Механика № 3

УДК 536.425

И.Л. Исупова

Пермский национальный исследовательский политехнический университет, Пермь, Россия

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФАЗОВЫХ ПРЕВРАЩЕНИЙ В СТАЛЯХ В РАМКАХ ПОДХОДА ДИФФУЗИОННОЙ ГРАНИЦЫ

В сталях наблюдаются все известные для твердого состояния фазовые превращения: полиморфное с широким спектром морфологических и кинетических особенностей; эвтектоидный распад (перлитное превращение); распад пересыщенных твердых растворов внедрения и замещения; упорядочение с изменением ближнего и дальнего порядка в аустените и мартенсите. Важная особенность данных систем заключается также в резко различающейся диффузионной подвижности металлических атомов и углерода, поэтому при превращениях перестройка кристаллической решетки может происходить наряду с диффузионным перераспределением углерода и легирующих элементов.

В представленной работе в рамках подхода диффузионной границы производится моделирование изменения структурного состояния сплавов, происходящего при термомеханической нагрузке. В подходе «диффузионной границы» форма и взаимное расположение областей, которые составляют микроструктуру, описываются непрерывными по пространству и времени функциями, переменными фазового поля, поэтому положение границы с течением времени неявно определяется изменением данных параметров. Изменение во времени переменных фазового поля описывается кинетическими уравнениями, которые были получены в рамках термодинамики необратимых процессов. Фаза рассматривается как область, которая в равновесном состоянии имеет вполне определенный термодинамический потенциал, отличный от потенциала других фаз. С целью выделения областей системы с различными значениями термодинамического потенциала вводится совокупность параметров, определяющих доли различных фаз. Учитываются также особенности твердофазного состояния, когда сильное межатомное взаимодействие вызывает при превращениях возникновение полей упругих напряжений, а стремление системы к снижению энергии упругой деформации обусловливает действие различных релаксационных механизмов, которые влияют на форму, ориентировку, взаимное расположение и внутреннюю структуру кристаллов новых фаз. В работе представлен алгоритм решения рассмотренной задачи и результаты численного эксперимента по моделированию изотермических мартенситных переходов в сталях.

Ключевые слова: стали, фазовые переходы, неоднородные системы, структура, неупругое деформирование.

I.L. Isupova

Perm National Research Polytechnic University, Perm, Russian Federation

MATHEMATICAL MODELING OF PHASE TRANSFORMATIONS IN STEELS IN FRAMEWORKS OF DIFFUSE INTERFACE DESCRIPTION

In steels all known solid state phase transformations are observed: polymorphic transformation with wide spectrum of morphological and kinetics features, eutectoid decomposition (pearlite transformation), decomposition of supersaturated solid solutions, short- and long - range ordering in austenite and martensite. The important feature of these systems is different diffusion mobility of metal and carbon atoms. Therefore, reorientation of the crystal lattice may occur simultaneously with the diffusion redistribution of carbon and alloying elements.

In this paper in the context of diffuse interface description is investigated microstructure evolution in steels during thermo-mechanical loading. In the frameworks of diffuse interface model the shape and mutual distribution of the grains is represented by functions that are continuous in space and time, the phase-field variables. Therefore, the evolution of the shape of the grains, or in other words the position of the interfaces, as a function of time, is implicitly given by the evolution of the phase-field variables. The temporal evolution of the phase- field variables is described by a set of kinetic equations, which are received in terms of thermodynamics of irreversible processes. Phase is thought as a labeling device for identifying regions with different state functions, and in this work the definition of a phase is taken to be a region of microstructure characterized by a homogeneous state function that differs from the functions of the other phases. Phase fraction variables are introduced for the purpose of associating regions of the system with different state functions. Also features of solid state are taken into account. Strong inter-atomic transformation is the causes of elastic stress and various relaxation mechanisms that influence on shape, orientation and structure of new phase. In the article the algorithm of the considered problem decision and results of the isothermal martensite transformation numerical experiments are submitted.

Keywords: steels, phase transformations, non-uniform systems, structure, inelastic deformation.

Введение

В системах Fe-C и Fe-Me-C наблюдаются все известные для твердого состояния фазовые превращения: перлитное, промежуточное (бейнитное) и мартенситное. Перлитное превращение при постоянной температуре начинается после некоторого инкубационного периода и при достаточной выдержке завершается полным распадом аустенита. Мартенситное превращение начинается без инкубационного периода с большой скоростью и при постоянной температуре прекращается по достижении определенной степени. Промежуточное (бейнитное) превращение, как и перлитное, характеризуется наличием инкубационного периода и температурной зависимостью. Промежуточное превращение, как и мартенситное, затухает при сохранении некоторой доли остаточного аустенита; степень превращения тем больше, чем ниже тем-

пература превращения. Изменение состава стали вызывают значительные изменения кинетики перлитного и промежуточного превращений, а также температурных интервалов превращения во всех трех областях.

Важная особенность системы Бе-Ме-С заключается в резко различающейся диффузионной подвижности металлических атомов и углерода. Значение коэффициента диффузии углерода в аустените на 4-5 порядков больше, чем коэффициенты диффузии (самодиффузии) атомов легирующих элементов. При превращениях переохлажденного аустенита переход гранецентрированной кубической (ГЦК) кристаллической решетки в объемно-центрированную тетрагональную (ОЦТ) решетку может происходить наряду с диффузионным перераспределением углерода и легирующих элементов. Переход у ^ а определяется полиморфизмом железа и сводится к перестройке решетки на границе раздела фаз (перемещению межфазной границы). Этот переход может осуществляться как по так называемому нормальному (если межфазная граница некогерентная), так и по мартенситному (если эта граница когерентная) механизмам.

Скорость диффузии легирующих элементов в аустените резко уменьшается с понижением температуры, и ниже 500-450 °С перераспределение их в процессе превращения практически исключается; скорость же диффузии углерода при этих температурах еще значительна.

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

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

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

1. Равновесная термодинамика концентрационно-неоднородных систем. Подход диффузионной границы

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

Отмеченные проблемы применения классической термодинамики для описания систем, имеющих сложную микроструктуру, привели к необходимости введения некоторых модификаций в классическую термодинамику. Одним из таких модифицированных подходов является равновесная термодинамика концентрационно-неоднородных систем [3-5, 9, 15]. Основная идея данного подхода заключается в том, что для адекватного описания гетерогенных систем необходимо в формализм классической термодинамики вводить характерные размерные масштабы. Эта идея реализована в подходе диффузионной границы, где для описания системы с границами, имеющими некоторую конечную толщину, в качестве независимых переменных вводятся градиенты параметров состояния [4].

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

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

2. Идентификация фаз. Доля фазы

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

С целью выделения областей системы с различными значениями термодинамического потенциала вводится совокупность параметров Ф, (I = 1,..., N, N - количество фаз), определяющих доли различных фаз. Совокупность характеризует распределение фаз и определяет фазовый состав в каждой точке рассматриваемой области. Значение параметра может изменяться от 0 до 1; ф, = 0 соответствует области, где

нет , фазы, ф, = 1 соответствует однофазной области. Таким образом,

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

многофазных системах в каждой материальной точке должно выполняться следующее равенство [13, 14]:

IФ, = 1. (1)

г-1

При описании многофазных систем делается предположение о том, что каждой фазе соответствует свой термодинамический потенциал (например, энергия Гиббса), который является функцией независимых термодинамических параметров ({х} ,{ф1, Ф2,..., фж}) [8]. При этом

термодинамический потенциал многофазной системы (фазы находятся в одинаковых внешних условиях) есть величина аддитивная; она зависит от доли фаз, составляющих систему. Таким образом, состояние термодинамического равновесия системы, состоящей из нескольких фаз, каждая из которых имеет соответственно доли ф(., характеризуется минимумом термодинамического потенциала:

/ ({Х} , {ф^ Ф2 — Фж }) = ! Ф/ ({х}). (2)

г-1

3. Определение термодинамического потенциала многофазной системы (неизотермический случай)

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

Предположим, что термодинамический потенциал бесконечно малого объема в неоднородной системе зависит как от значений параметров состояния в рассматриваемой области, так и от их значений в соседних областях. В связи с этим в выражение для потенциала вводятся не только параметры состояния, но и их градиенты (в работе [1] подобным образом было получено выражение для свободной энергии негомогенной системы).

Для анализа необходимо определить, какие параметры в конкретных условиях исследуемой проблемы влияют на поведение системы. Речь идет о независимых параметрах, от которых зависит термодинамический потенциал. Как было сказано выше, термодинамический потенциал гетерогенной системы в значительной степени зависит от ее фазового состава. В рассматриваемой задаче параметрами, характеризующими пространственное распределение фаз, являются фазовые доли ф1, ф2,.... Предполагается, что потенциал зависит не только от локальных значений фазовых долей, но и от их градиентов определенного порядка. Для того чтобы рассматривать системы с границами, имеющими некоторую конечную толщину, необходимо в формализм классической термодинамики ввести характерные размерные масштабы. Это делается путем введения градиентов фазовых долей Уф1,Уф2,...,УУф1,УУф2,.... Один из важнейших

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

Возьмем энтропию и разложим ее в ряд Тейлора вблизи точки

50 - 5(ф1,ф2,...,и1,и2,...,е0 (ге,0,ф1,ф2,...,и1,и2,...),0,0,...,0,0,...), т.е. точки, в которой все параметры постоянны, а градиенты равны нулю:

5(ф1 ,ф2 ,... ,и1 ,и2 ,... ,е0 (^ ,0,ф1 ,ф2 ,... ,и1 ,и2 ,...),

Уф1 ,Уф2 ,.. .,УУф1 ,УУф2 ,...) =

— 50 (ф1 ,ф2 ,... ,и1 ,и2 ,... ,е0 (^ ,0,ф1 ,ф 2 ,... ,и1 ,и2 ,...)) -55 V? 55 V? 55 Х7Х7

5Уф1 ’ ф1 5Уф2 ’ ф2 ... 5УУф1 : ф1 (3)

55 : УУф2 -... - 1у* • •Уф - ^Уф2 • ^ -Уф2 -

5УУф2 2 5Уф15Уф1 2 5Уф25Уф2

52 5

-...-Уф,------------Уф2 -...

5Уф15Уф2

Здесь е0 - внутренняя энергия.

дэ

В общем случае коэффициенты -------------- являются векторами,

дУфі

а ———,-----------------------------------------—- - тензорами второго ранга (і, 1=\..Ы). Обозначим

дУУфі дУфідУф,

дэ дэ д 2э

дУфі дУУфг дУфідУф] шем (3) в виде

через Р(і), Н (і), К (ч) соответственно и перепи-

5 (ф^ ф2 , •••, ul, и2 — е0 ( £ Є , ^ Фl, Ф2,•••, U1, «2 , •••) , УФ1, Уф 2 ,^ ••, УУф^ УУФ2,•••) =

50 (фl, Ф2,•••, u2,•••, е0 (, ^ Фl, ф2 , •••, ul, «2 , •••))

N N

(4)

N Л N N

-I Р(,) - Уф, - т I (Уф, )Т • К(»! • (Уф 1) -1 И"1: УУф1'<.

,=1 2 г,у =1 г=1

Последнее слагаемое в выражении (4) можно записать следующим образом:

Н

(і)

: УУф(і) = У.(Н(і) • Уфі )-(ух • Н(і)) • Уфг

N Т 5и(') (5)

-Е(уфу) -5ф--(уф,)■

Запишем выражение для энтропии с учетом равенства (5), два последних слагаемых которого приведены в терминах Р(1) и К(ч):

5 = 1( 50 (фl, Ф2 , •••,«1, и2 — е0 (, ^ Фl, Ф2,•••, ul, «2 ,•••))

V

N 1 N N ч (6)

1 Р(,) - Уф, -Т 1 (у<ф )т • К(") (Уф,) - 1у (Н(і) - Уф,)) ЛV,

і=1 2 і,,=1 і=1

Воспользовавшись теоремой о дивергенции, последнее слагаемое равенства (6) мы можем записать следующим образом:

]Ху-( И(,)^ф, ) = п •( И(,)-Уф,) <15. (7)

V г=1 5 I =1

Внешнюю поверхность 5 можно выбрать таким образом, чтобы поток Уф, через нее был бы равен нулю, поэтому последнее слагаемое из выражения для свободной энергии можно исключить. Так как плотность энтропии должна быть инвариантна к преобразованию Уф, —— —Уф, (симметрия отражения), то Р(1) = 0 для всех ,. С учетом всех рассуждений, приведенных выше, для выражения энтропии окончательно имеем

5 = Д ^0 (фр ф2,-’и1 ’ и2 ’■■■, е0 (£Й , ^ Фl, ф2,•••, Щ,~))-

V

1 N

-Х(уф,)т К1’)^(уф,)

2 =1 у

(8)

dV.

4. Кинетические уравнения

Изменение фазовой доли описывается следующим кинетическим уравнением [2, 12]:

^ = I*5-, (9)

дг 5ф,

где I - положительный кинетический параметр, связанный с подвижностью границ ц.

Это соотношение можно получить в рамках термодинамики необратимых процессов, которые могут возникать по различным причинам. В термодинамике все величины, которые могут вызвать необратимое явление, носят названия сил и обозначаются X (температурный градиент, градиент концентрации и т. п.). Силы являются причиной возникновения некоторых потоков, которые обозначаются через 3 (поток тепла, диффузионный поток и т. п.). В общем случае любая сила может вызвать любой поток. Это обстоятельство позволяет выразить любое необратимое явление в общем виде следующим феноменологическим соотношением:

3 = IX.

В нашем случае в качестве обобщенной силы можно принять

-55 = X, а в качестве обобщенной скорости - = 3.

5ф, дг

Для того чтобы в многофазной системе в каждой точке выполнялось равенство (1), необходимо применить технику множителей Лагранжа. Тогда выражение (9) запишется следующим образом:

дг

5ф,

(10)

где

5' = 5 + Ы£ф, — 1

dV,

(11)

где X - множители Лагранжа.

Для определения X нужно воспользоваться следующей системой уравнений:

55

5ф,

- + Х = 0,, = 1.^.

(12)

Из (12) находим

(13)

Подставляя последнее выражение в (12), получим следующее выражение для определения изменения фазовой доли в многофазной системе:

55 55

5ф, 5ф’ у

I

=-£ N’1

дsn

1

дф, 2 к г=1

Куф. )т

дК(к ,г) дф,

(Уфг ) +

+£у^( К(,к )^Уфк )■

к=1

+1 £(Уфк )т дК

дф ’ 2 ^к)

(к Г)

дф’

•(уф, ) —

(14)

£у( к ’ ^)

к=1

Изменение концентрации к-го компонента (углерода или легирующих элементов) определяется следующим образом:

ик =-У-1Ф; ]

/=1

к

(15)

где ^ - диффузионный поток к-го компонента в г фазе, определяемый согласно термодинамике необратимых процессов:

]к=мк -V

г г

у

(16)

где Мк - диффузионная мобильность соответствующего компонента в г фазе.

Подставив выражение для диффузионного потока (16) в уравнение диффузии (15), получим следующее соотношение для изменения концентрации углерода (легирующих) элементов:

N

V

!фгМ* -V

г=1

(17)

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

йео = 0ё^о +

где0л

К дщ у

( део Л ЧдФ1 У

ёф1 +... +

део Чд£< У

ё£е.

(18)

Используя связь между свободной и внутренней энергиями е0 = /0 + 05о, можно получить

У> = -^0 +

\ди1 у

"део Л

ЧдФ1 у

ёф1 +... +

део чд^ у

ё£е.

(19)

Из этих соотношений можно установить следующие зависимо-

сти:

д^о Л 1 Г део Л = 1 'ж \ У

дщ у = -0 чдщ1 у = 0 Удщ1 у

до '' 1 ( део " 1 (д

дФ1 у = -0 д 1 у = 0 1дФ1 у

и т.д.,

и т.д.

(2о)

Теперь в терминах свободной энергии кинетические уравнения

(14), (17) запишутся следующим образом:

дг N7=1

}-1 V

0 дф,

£^( К (ік >.Уф,

к-1

.1 / 0 дф

+ £у.( К »к ».Уф,

к-1

(21)

Кк -V-

N

і=1

Еф, м к -V

1 д/с

\\

с і

Ч0д»ік уу

(22)

При записи выражения (21) предполагалось, что тензор К(іі) не зависит от фазового состава.

Запишем плотность свободной энергии в виде суммы плотностей упругой /Е и химической /сн составляющих свободной энергии:

/с (фр ф2>-,К1> «2,-,£Є , 0) = Ґ1 + /СН =

N / \ N

= Еф, (в; :С,: в;) + Еф/ (щ,..,0),

(23)

=1 ' г=1

где / - удельная свободная энергия для отдельной фазы.

Определим удельную свободную энергию для отдельной фазы / . Для этого воспользуемся следующей моделью [6, 1 о]:

г ( п\ /’чист . /*смеш . /’взаим /ллч

у(щ,u1,..., 0) = / + / + / , (24)

где /чист - это свободная энергия в отдельной фазе; /смеш - изменение свободной энергии за счет смешивания компонент; /юаим - изменение свободной энергии за счет химического взаимодействия отдельных компонент.

Первое слагаемое в (24) может быть определено следующим образом:

п • / ч

/= Е К/ (0),

к=1

(25)

где ик - концентрация к-го компонента стали (железа, углерода, легирующих элементов, примесей), измеряемая в молярных долях; /I (Т) -

свободная энергия к-го компонента в I фазе, которая может быть достаточно точно аппроксимирована следующей функцией [7]:

/к (е) = Ак + ВкТ + СкТ 1п (е). (26)

Второе слагаемое в выражении (24) отвечает за изменение свободной энергии за счет смешивания компонент. Случайное изменение расположения атомов является причиной возникновения конфигурационной энтропии, изменение которой можно вычислить через логарифм числа возможных перестановок (для атомов разновидности 1 и 2):

С . Л

конф

= к 1п () = к 1п

П!

V п1!п2! У

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

5конф = _^ (^ 1п (и1) + и2 1п (и2 )) ,

где Я = ША, ЫА - число Авогадро.

Для всей системы мы получим

5 конф

= _ ик 1п (ик )•

к=1

Тогда изменение свободной энергии за счет смешивания компонент определится следующим образом:

/смеш = ксыеш _ е5смеш = яе£ик 1п (ик). (27)

к=1

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

Третье слагаемое (24) описывает изменение свободной энергии за счет химического взаимодействия отдельных компонент. Будем рассматривать только парное взаимодействие и применим модель Редли-ска-Кистера для описания взаимодействия компонент друг с другом.

Т огда изменение свободной энергии за счет химического взаимодействия компонент определится следующим образом:

п-1 п .

/взаим / взаим т? взаим х-* V-* тг /Л\ /ло\

= к - Т -1 I ики]Ьк] (0), (28)

к-1 ] -к+1

где 5взаим - 0, так как энтропия системы полностью определяется ее конфигурационной составляющей; Ик] (0) - параметр, описывающий взаимодействие между к и у компонентами в соответствующих фазах.

5. Температурная задача

Тепловой поток при термической обработке может быть описан законом теплопроводности Фурье:

д0

к —+ У- ч - 2, ч --ХУ0, (29)

дг

где к - рег (р - плотность, ег - удельная теплоемкость); 0 - температура; ч - тепловой поток; г - мощность внутреннего теплового источника; X - теплопроводность. В данном случае мощность внутреннего

теплового источника определяется скрытой теплотой фазового перехода и находится следующим образом:

2 -Iн, Ф,, (30)

г -1

где Н, - скрытая теплота фазового перехода.

Теплоемкость и теплопроводность в многофазной системе определяются по правилу смеси:

к-1 Ф,к, (31)

, -1

Х-1 Ф, X, (32)

, -1

где к,,X, - теплоемкость и теплопроводность в отдельной фазе.

6. Описание неупругого деформирования поликристаллического образца

В качестве базового определяющего соотношения используется обобщенный закон Гука в скоростной релаксационной форме:

аг = а - ю • а + а • ю = c

: (d - din),

d =1 (Уу + Уут), (33)

I УУ + V’Т

2

V-а = 0

где а - тензор напряжений Коши; ^ $п - тензор деформации скорости и его неупругая составляющая; с - тензор упругих свойств, ю - тензор спина решетки.

Тензор скорости неупругих деформаций многофазной системы определяется следующим образом:

= I ф, I у(*) р1") + I ф,< + «0, (34)

к=1

1

•(*)

У Г = У 0

х( *)

х( к)

(х(к)), * = 1...П,і = 1..^1, (35)

где I* - градиент трансформационной деформации, описывающий превращение одной фазы в другую по г-му варианту; N - количество фаз, испытывающих пластическую деформацию; И2 - количество фаз,

• (*) (*)

испытывающих превращение; у г ;, х^' - накопленный сдвиг и критическое напряжение сдвига по к-й системе скольжения;

Р<‘> = 2(.«*'Ь(‘) + Ь(к'п(к') - ориентационный к-й системы сколь-

жения; п(к) - единичный вектор нормали к плоскости скольжения; Ь(к) -единичный вектор по направлению вектора Бюргерса, характеризующий направление скольжения; у0 . - скорость деформирования, являющаяся параметром материала и характеризующая скорость сдвига по системе скольжения при касательном напряжении, равном критическому напряжению сдвига; х(к) = р/к): о - сдвиговое напряжение, действующее в системе скольжения к; хс ^к) > 0 - сопротивление сдвигу в

системе скольжения к, представляющее собой материальную функцию параметров нагружения; т. - чувствительность кристаллов к скорости приложения нагрузки; а - тензор термического расширения.

т

с і

Для определения тензора упругих свойств и тензора термического расширения в многофазной системе используется правило смеси:

с = 1 ф;с;, (36)

і=1

а = 1 ф;а, (37)

і=1

С, аі - тензор упругих свойств и тензор термического расширения в отдельной фазе.

Для вычисления спина решетки используется модель Тейлора. При использовании модели Тейлора тензор спина решетки определяется как разность тензора вихря и антисимметричной части тензора скоростей неупругих сдвигов по системам скольжения элемента:

1 N п / ч

ю І ^ ік} Ик )ьік} + )п(к})’ (3 8)

2 і=1 1 ’

к=1

где w = 2 (уу -(Уу )г ) - тензор вихря.

Тензор поворота, который переводит лабораторную систему координат в систему координат, связанную с решеткой аустенита, получаем в результате решения следующего уравнения:

ю = о • сг. (39)

7. Алгоритм решения

Для решения системы уравнений (21), (22), (29), (33) используется следующий конечно-элементный алгоритм. Общее время наблюдения ґ разбивается на М малых шагов At. При решении этих уравнений делаются следующие предположения: а) для уравнения (21) доли фаз, напряжения, температура и концентрации (углерода и легирующих элементов) фиксированы и известны с предыдущего временного шага (уравнение решается для (N-1) сосуществующих фаз, а фазовая доля Ы-й фазы определяется из условия (1)); б) для уравнений (22), (29), (33) все фазовые доли фиксированы и известны из решения N уравнений (21) на данном временном шаге. Благодаря этим предположениям на каждом шаге интегрирования каждое из уравнений (21) может быть решено отдельно от уравнений (22), (29), (33).

Ниже приведен подробный конечно-элементный алгоритм реше-

ния.

Разбиваем исследуемую область на конечные элементы и определяем шаг по времени. При этом разбиение по пространству и времени зависит от специфики структуры исследуемого материала и характерного времени процесса превращения. Пусть на некотором временном шаге t нам известны напряжения, доли фаз, температура и концентрации углерода (легирующих элементов). Необходимо определить их значения на шаге t + At.

1) Решаем отдельно каждое из N уравнений (21) и определяем доли фаз на временном шаге t + At:

[с - ]{«ф}

Здесь [с* ] = А[с*е ];

Е([ к

к=1

(* у

к

{фк } + {б/.

к

'(*)'

А [к*(1к)е]; {о } = А{е;};

С ] = Л N Г [ N ] й¥е;

к

(к )е

|[В]т [к(;)][В]№,;

о.}=1Л ■М' ]т-(АТС"+ АЕ) «V

0

где [ N ] - матрица функций формы, [ В] - матрица производных функций формы, знак А означает ансамблирование.

2) При моделировании диффузионных фазовых превращений, в процессе которых происходит перераспределение атомов углерода (легирующих элементов), необходимо решить уравнение (22)

[С ]{£ }4 Е [ в]т ф *А [м; ][ в]{

3) При рассмотрении неизотермических фазовых превращений решаем уравнение (29) для определения распределения температуры

с ]{¥ Ик ]{®}-{2}=0

где

[С] = А[с*]; [К] = А[к-], {2} = А [г-]; [с*] = {[^ к[N]ЛУ*;

У*

[К*]=/[ВГ И[В]■

У*

где к, [X] - теплоемкость и матрица коэффициентов теплопроводности материала.

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

{дЕ'}' *"=2 ф, 2 ду с д р<11.

г к=1

Для определения скорости сдвиговой деформации по каждой системе скольжения применяем уравнение

+д?

(t+Дt \

х(к )г ) ’

У (+ ?=у

0 г

Х(к )г

х(к)

с ,

х(ГД' = р(к): оt, р(к) =1 (п(к)ъ(к) + ъ(к)п(к)).

(к)г г 7 г \ I I II)

где , 2'

Теперь скорость неупругих деформаций в многофазной системе определяется из уравнения

К г+д=2ф+дt {д<}t+д+£ дф, к}+«-+ддв,

г г

N

Л+д^ л+Дt„

где « =2ф, «г.

г =1

Далее решаем систему уравнений для определения изменения деформаций:

[К]{дц}+д = {Щ‘+д.

Здесь [К ] = А [к* ] - матрица жесткости; [К* ] = А в. Г [С][ в. ] Лу,,

У*

[ В, ] - стандартная конечно-элементная В-матрица; [С] - матрица упругих свойств (тензор упругих свойств в многофазной системе опреде-

ляется как c'+A' = ^ф'+A'c;); {Au} - вектор изменения деформаций;

'=1

R} = А{я}е - вектор нагрузок;

{Re}'+" = !{[Б,]T [C]';A' {Ae'"}'+A' +[Б,]T [T]'

V,

-[Б, ]T {а}'} dV, + J{AF}d5,;

-,'+а' {а}'

[T ] =

о о о -2Аюі2 -2Аю13 о

о о о 2Аюі2 о -2Аю23

о о о о 2Аю13 2Аю23

A®12 -A®12 о о -Аю23 -Аю13

Аю13 о -Аю13 Аю23 о -А®і2

о Аю23 -Аю23 □ю13 А®і2 о

A»=2 ((y(A“ ))T -v(Au))-2 £ ф' *A]T Ay;;£(n(k )b(k)- b(kin

2 V ' 2 i=1 k=1

Теперь можно определить напряжения на шаге ' + A':

A' = [c]'+A ([Б,]{Au}-{Ae"}'+A')-[T] ;А{а}' + {а}'.

(k)

Для интегрирования соотношения (39) применяется следующий алгоритм. Сначала определяется вектор угловой скорости / как ассо-

г 1 г

циированный вектор для тензора спина ю, / = ^ є: ю , по которому

определяется единичный вектор оси скорости (приращения) поворота ДГ+Д

e=

Ar

'+A'

(Ar = /A ') и угол P , на который необходимо повернуть

решетку из текущего положения. Ортогональный тензор ротации AR восстанавливается по найденным значениям компонент вектора e , относительно которого тензор ротации AR поворачивает все пространство, и угла такого поворота P :

AR' = cos P'I + (l - cos P') ee + sin p'e x I.

S

F

С использованием найденного тензора ротации переопределяется ориентационный тензор:

Ог+д = дя. ог.

8. Применение модели к описанию мартенситных переходов

Рассматривается изотермическое мартенситное превращение в стали, состав которой приведен в табл. 1. Предполагается, что в процессе превращения может образовываться мартенсит трех различных типов (рис. 1). Пластически деформироваться может только аустенит. Все параметры, необходимые для проведения моделирования, представлены в табл. 2.

Таблица 1

Химический состав стали 204М

С Сг Мп Б Р N1 Бі N Бе

Массовый процент 0,038 15,4 8,92 0,001 0,017 2,75 0,44 0,236 Ва1.

Молярная доля 0,00172 0,16145 0,08851 0,00002 0,00029 0,02554 0,00854 0,00918 0,70473

Таблица 2

Параметры для моделирования мартенситного перехода [11]

Аустенит

Параметр решетки, А ау= 3,591

А Упругие модули С , Гпа кА = 242,0; к А = 156,5; < = 136,0

Скорость деформирования, с-1 г 0 = 10-3

Чувствительность т = 0,005

Сопротивление сдвигу, МПа тс = 150

Мартенсит

Параметры решетки, А аа = 2,875

Упругие модули <, сМ , сМ , ГПа (в базисе, связанном с собственной решеткой) < = 497,0; кМ = 405,0; кМ = 265,0; кМ = 617,0; кМ = 263,0; кМ = 287,0

Рассматриваются следующие варианты начальных условий для задачи определения фазовых долей трех разновидностей мартенсита:

1. Во всей расчетной области все фазовые доли различных разновидностей мартенсита равны нулю за исключением нескольких узлов, в которых заданы равные доли для всех типов (рис. 2, а).

Вид решетки аустенита и мартенсита

Матрицы упругих свойств аустенита и мартенсита (упругие свойства мартенсита определены в базисе, связанном с его собственной решеткой)

Трансформационные матрицы (определены в базисе, связанном с решеткой аустенита)

"12аа аа

Л1 =-, Лз = —

кА кА КА 0 0 0

к2" кА КА 0 0 0

кА К2А кА 0 0 0

0 0 0 кА 0 0

0 0 0 0 кА 0

0 0 0 0 0 кА

кГ кГ кГ 0 0 0

кМ кГ кГ 0 0 0

кГ кГ кГ 0 0 0

0 0 0 кГ 0 0

0 0 0 0 кГ к6 0

0 0 0 0 0 кГ

Л3 0 0

0 Л1 о

0 0 л1

к кГ к2Г 0 0 0

к кГ кГ 0 0 0

кг КзГ кГ 0 0 0

0 0 0 кГ 0 0

0 0 0 0 К 0

0 0 0 0 0 кГ к6

к кГ кГ 0 0 0

кГ кГ кГ 0 0 0

кГ кГ кГ К4 0 0 0

0 0 0 кГ к6 0 0

0 0 0 0 кГ 0

0 0 0 0 0 кГ

Л1 о о

0 Л3 0

0 0 л1

/з =

Л1 0 0

0 Л1 0

0 0 Лз

а

а

А

с

с, =

с =

2

С =

Рис. 1. Превращение кубической решетки аустенита в тетрагональную решетку мартенсита

2. Во всей расчетной области все фазовые доли различных разновидностей мартенсита равны нулю за исключением нескольких узлов, в которых находится единственный вариант мартенсита (третий) (рис. 2, б).

ъ дф,-

Все уравнения решаются при условии равенства потока —— нулю

дп

вдоль всей границы образца (п - внешняя нормаль к границе образца).

Граничные условия для задачи определения напряженно-деформированного состояния следующие:

1. Равенство нулю нормальной составляющей изменения перемещений Ли | „ = 0.

п '^Ли

2. Равенство нулю касательной составляющей поверхностных нагрузок ЛРХ |5лр = 0 .

а б

Рис. 2. Начальные условия для задачи определения фазовой доли (справа приведена градиентная шкала для доли фаз)

На рис. 3-6 приведены результаты расчета.

Шт

а

Ф'Ш*

в

Рис. 3. Распределение мартенсита в монокристалле: а - первого варианта; б -второго варианта; в - третьего варианта (начальные условия приведены на рис. 2, а)

Рис. 4. Распределение мартенсита в поликристалле: а - первого варианта; б -второго варианта; в - третьего варианта (начальные условия приведены на рис. 2, а)

б

Рис. 5. Распределение мартенсита в монокристалле: а - первого варианта; б -второго варианта; в - третьего варианта (начальные условия приведены на рис. 2, б)

Рис. 6. Распределение мартенсита в поликристалле: а - первого варианта; б -второго варианта; в - третьего варианта (начальные условия приведены на рис. 2, б)

При наличии в образце нескольких областей, которые потенциально могут стать местами образования любой из разновидностей мартенсита, в монокристалле наблюдаются все типы мартенсита, причем их доля примерно одинакова. При рассмотрении поликристалла с теми же начальными условиями и с учетом того, что в процессе деформирования могут изменяться ориентации зерен, мы видим, что зерна могут сориентироваться так, что становится невозможным образование некоторых разновидностей мартенсита. В рассмотренном случае образовались только 2 варианта мартенсита.

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

смотренном выше примере помимо третьей разновидности мартенсита образовалась и первая.

Данные для вычисления химической энергии были взяты из работ [7, 10].

Заключение

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

Работа выполнена при финансовой поддержке ФЦП «Научные и научно-педагогические кадры инновационной России на 20092013 годы» (мероприятие 1.2.2, Соглашение 14.B37.21.0382).

Библиографический список

1. Исупова И.Л., Трусов П.В. Математическое моделирование фазовых превращений в сталях при термомеханической нагрузке // Вестник ПНИПУ. Механика. - Пермь: Изд-во Перм. нац. исслед. политехн. ун-та, 2011. - № 4. - С. 62-78.

2. Allen S.M., Cahn J.W. Microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening // Acta Metallur-gia. - 1979. - Vol. 27(6). - P. 1085-1095.

3. Cahn J.W. Adopting thermodynamics to material science problems // J. Phase Equilibria. - 1994. - Vol. 15. - P. 373-379.

4. Cahn J.W., Hilliard J.E. Free energy of a non-uniform systems.

I. Interfacial free energy // J. Chem. Phys. - 1958. - Vol. 28. - P. 258-266.

5. Cahn J.W., Hilliard J.E. Free energy of a non-uniform systems. III. Nucleation in a two-component incompressible fluid // J. Chem. Phys. -1959. - Vol. 31. - P. 688-699.

6. Chuang Y.-Y., Chang Y.A. A Thermodynamic Analysis and Calculation of the Fe-Ni-Cr Phase Diagram // Metall. Mater. Trans. - 1987. -Vol. 18A. - P. 733.

7. Dinsdale A. T. SGTE Data for Pure Elements // Calphad. - 1991. -Vol. 15. - P. 317-425.

8. Guggenheim E.A. Thermodynamics: an advanced treatment for chemists and physicists. - North-Holland Publishing Company. 1967.

9. Larche F.C., Cahn J.W. The interactions of composition and stress in crystalline solids // Acta Metall. - 1985. - Vol. 33. - P. 331-367.

10. Miettinen J. Calculation of solidification-related thermophysical properties for steels // Metall. Mater. Trans B. -1997. - Vol. 28A. -P. 281-297.

11. Myoung-Gyu Lee, Sung-Joon Kim, Heung Nam Ham Crystal plasticity finite element modeling of mechanically induced martensitic transformation (MIMT) in metastable austenite//International Journal of Plasticity. - 2010. - Vol. 26. - P. 688-710.

12. Steinbach I., Apel M. Multi-phase field model for solid state transformation with elastic strain // Physica D. - 2006. - Vol. 217. -P. 153-160.

13. A phase field concept for multi-phase systems / I. Steinbach [et al.] // Physica D. - 1996. - Vol. 94. - P. 135-147.

14. The multi-phase field model with an integrated concept for modeling solute diffusion / J. Tiaden [et al.]// Physica D. - 1998. - Vol. 115. -P. 73-86.

15. Wang J., Chen L.-Q., Khachaturyan A.G. Kinetics of the Strain-Induced Morphological Transformation in Cubic Alloys with Miscibility Gap // Acta Metall. Mater. - 1993. - Vol. 41. - P. 279-296.

References

1. Isupova I.L., Trusov P.V. Matematicheskoe modelirovanie fasovyh perehodov v ctalyah pri termomechanicheskoi nagruske [Mathematical modeling of phase transformations in steels during thermomechanical loading]. Vestnik Permskogo Natsinalnogo Issledovatelskogo Politehnicheskogo Universiteta. Mehanika, 2011, no. 4, pp. 62-78.

2. Allen S.M., Cahn J.W. Microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metallur-gia, 1979, vol. 27(6), pp. 1085-1095.

3. Cahn J.W. Adopting thermodynamics to material science problems. J. Phase Equilibria, 1994, vol. 15, pp. 373-379.

4. Cahn J.W., Hilliard J.E. Free energy of a non-uniform systems. I. Interfacial free energy. J. Chem. Phys, 1958, vol. 28, p. 258.

5. Cahn J.W., Hilliard J.E. Free energy of a non-uniform systems.

III. Nucleation in a two-component incompressible fluid. J. Chem. Phys, 1959, vol. 31, p. 688.

6. Chuang Y.-Y., Chang Y.A. A Thermodynamic Analysis and Calculation of the Fe-Ni-Cr Phase Diagram. Metall. Mater. Trans, 1987, vol. 18A, p. 733.

7. Dinsdale A.T. SGTE Data for Pure Elements. Calphad. 1991, vol. 15, pp. 317-425.

8. Guggenheim E.A. Thermodynamics: an advanced treatment for chemists and physicists, North-Holland Publishing Company, 1967.

9. Larche F.C., Cahn J.W. The interactions of composition and stress in crystalline solids. Acta Metall, 1985, vol. 33, p. 331-367.

10. Miettinen J. Calculation of solidification-related thermophysical properties for steels. Metall. Mater. Trans B, 1997, vol. 28A, p. 281.

11. Myoung-Gyu Lee, Sung-Joon Kim, Heung Nam Ham Crystal plasticity finite element modeling of mechanically induced martensitic transformation (MIMT) in metastable austenite. International Journal of Plasticity, 2010, vol. 26, pp. 688-710.

12. Steinbach I., Apel M. Multi-phase field model for solid state transformation with elastic strain. Physica D, 2006, vol. 217, pp. 153-160.

13. Steinbach I., Pezzolla F., Nestler B., Seebelber M., Prieler R., Schmitz G. J., Rezende J. L. A phase field concept for multi-phase systems. Physica D, 1996, vol. 94, pp. 135-147.

14. Tiaden J., Nestler B., Diepers H., Steinbach I. The multi-phase field model with an integrated concept for modeling solute diffusion. Physica D, 1998, vol. 115, pp. 73-86.

15. Wang J., Chen L.-Q., Khachaturyan A.G. Kinetics of the Strain-Induced Morphological Transformation in Cubic Alloys with Miscibility Gap. Acta Metall. Mater, 1993, vol. 41, pp. 279-296.

Об авторах

Исупова Ирина Леонидовна (Пермь, Россия) - аспирант кафедры математического моделирования систем и процессов Пермского национального исследовательского политехнического университета (614990, г. Пермь, Комсомольский пр., 29, e-mail: [email protected]).

About the authors

Isupova Irina Leonidovna (Perm, Russian Fderation) - postgraduate student of Department of Mathematical Modeling of Systems and Processes, Perm National Research Polytechnic University (29 Komsomolsky av., Perm, 614990, Russian Fеderation, e-mail: [email protected]).

Получено 9.09.2012