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

2013 Механика № 3

УДК 539.3

И.Л. Исупова, П.В. Трусов

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

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

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

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

I.L. Isupova, P.V. Trusov

Perm National Research Polytechnic University, Perm, Russian Federation

MATHEMATICAL MODELING OF PHASE TRANSFORMATIONS IN STEEL UNDER THERMOMECHANICAL LOADING

The article is devoted to the mathematical modeling of phase transformations in steels under thermomechanical loading. To construct the model, a multilevel approach based on the use of internal variables in its structure - the parameters characterizing the state and evolution of meso- and micro-

structure of the material - has been applied. The proposed model can be used to describe diffusionless (martensitic) and diffusion phase transformations. When modeling the diffusion phase transformations take into account the fact that, along with the restructuring of the crystal lattice can be redistribution of carbon and alloying elements. Statement of the general problem of individual facilitated the release of subtasks - namely, determining the stress-strain state, temperature, and describe the redistribution of alloying elements, which can be formulated in relation to independent productions. For the considered sub-tasks include different types of models. So much for the problem of determining the stress-strain state uses a direct model of the second type, and for the problems of heat conduction and diffusion - a direct model of the first type. In this case, the problem of describing the redistribution of atoms of carbon and alloying elements is considered only at the meso level, since it is at this scale diffusion processes are significant. The paper presents the general structure of the two-level model. For the problem of determining the stress-strain state, problems of heat conduction and diffusion are all considered statement on the scales. Special attention is given to obtaining a kinetic equation for the change of the volume fractions of the coexisting phases.

Keywords: phase transformation, steel, multi-level model, the internal variables.

Введение

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

На сегодняшний день в литературе представлено большое количество работ, посвященных математическому моделированию как без-диффузионных (мартенситных), так и диффузионных фазовых превращений в сталях [2]. В большинстве работ, посвященных описанию термомеханических процессов, используются так называемые прямые модели первого типа [6, 7], когда каждому зерну ставится в соответствие совокупность нескольких конечных элементов и для каждого из элементов напрямую используется та или иная модель для описания фазовых переходов. Использование таких моделей для моделирования реальных процессов в трехмерной постановке требует существенных вычислительных ресурсов. Еще один из недостатков большинства таких моделей - это отсутствие прямого включения механизмов деформирования и фазовых превращений и их носителей, присущих в реальном материале различным структурным уровням. В моделях, где процессы фазовых превращений рассматриваются на различных структурных уровнях, моделирование обычно ведется в области, состоящей из нескольких зерен.

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

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

1. Структура двухуровневой модели

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

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

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

На макроуровне I ставится и решается краевая задача в обычном смысле (включающая в свою постановку уравнения равновесия, уравнения теплопроводности, необходимые краевые условия). Для задачи определения напряженно-деформированного состояния используется прямая модель второго типа (рис. 1), т. е. каждой макроскопической точке интегрирования ставится в соответствие элемент макроуровня II [6, 7]. Элемент макроуровня II рассматривается как статистическая выборка соответствующих элементов мезоуровня; параметры нагружения для модели нижнего уровня назначаются с верхнего уровня, а параметры эволюционирующей структуры, текущие физико-механические свойства, неупругие деформации на макроуровне II, которые далее используются при решении краевой задачи на макроуровне I, определяются тем или иным способом осреднения соответствующих параметров мезоуровня.

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

Для температурной задачи используется прямая модель первого типа [6, 7] (рис. 2).

КЭ анализ

Рис. 2. Схематическое представление масштабных уровней для температурной задачи

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

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

Процесс диффузии на больших масштабах происходит очень медленно, поэтому на макроуровне I задача о перераспределении атомов углерода и легирующих элементов не ставится и не решается, а принимается гипотеза об однородном распределении концентрации во всем теле. При решении задачи на макроуровне II предлагается ис-

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

Отметим, что величины макроуровня обозначаются заглавными буквами, а мезоуровня - соответствующими строчными буквами.

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

При построении несимметричной теории встает вопрос о выборе несимметричной меры скорости деформации, которая должна быть независимой от выбора системы отсчета. В качестве меры скорости изменения деформационного состояния предлагается использовать градиент скорости относительного движения частиц [11], что согласуется с используемым определяющим соотношением - законом Гука в скоростной релаксационной форме, записанным с позиции наблюдателя в системе отсчета, связанной с решеткой аустенита (компоненты тензора упругих свойств «привязаны» к решетке аустенита, в качестве меры скорости изменения напряженного состояния принята коротационная («решеточная») производная тензора напряжений Коши).

2. Постановка задачи на макроуровне I

На макроуровне I необходимо определить перемещения, напряжения, деформации и распределение температуры, удовлетворяющие следующей системе уравнений [4]:

уравнению равновесия в скоростях

У-Ё-У-

определяющему соотношению

2сЯ = 2 + 2 • а + аТ • 2 = П :

кинематическим соотношениям

Zr = Z - а, Z = УУ;

уравнению теплопроводности

К0-У-

граничным условиям

V = УГ на Гу ;

N • 2-

0 = 0г на Г<;

N •

N •

/О

при начальных условиях

и

У • 2о + Во = 0,

N • 2о = То,

0

В соотношениях (1)—(3) 2 - тензор напряжений Коши макроуровня; В — вектор распределенных объемных сил; V — скорость перемещения; П — тензор упругих свойств; А — тензор термического расширения; Zг, Zгr” — транспонированный градиент относительной

скорости перемещений, его неупругая составляющая; а — тензор спина макроуровня; верхний индекс «СЯ» обозначает коротационную производную. Вся поверхность тела Г разделена на две части — Г2 и Гу . На Г2 задаются силовые граничные условия, а на Гу — кинема-

тические. В (5)—(6) и (10)—(12) УГ — значение скорости перемещения на части границы Гу ; N — вектор внешней нормали к поверхности; Т —

вектор поверхностных сил; и — вектор перемещений; нижний индекс «0» указывает на значения величин в начальный момент времени. В уравнениях, описывающих температурные поля ((4), (7)—(9), (13)), 0,00,0 Г — температура, ее значение в начальный момент времени и на границе; К — теплоемкость; У — мощность внутреннего теплового источника; Л — тензор теплопроводности; Г0 — часть границы, на которой задается температура; ГО — часть границы, на которой задается проекция О вектора теплового потока Є = Л •УО на внешнюю нормаль; Г0/ — часть границы, на которой задается условие теплообмена /о

с окружающей средой по закону Ньютона, к — коэффициент теплообмена, 0Ш — температура окружающей среды.

Тензор упругих свойств П, тензор термического расширения А, тензор спина а, неупругая составляющая транспонированного градиента относительной скорости перемещений Zгr”, теплоемкость К, мощность внутреннего теплового источника У и тензор теплопроводности Л в каждой точке интегрирования определяются из моделей для представительного объема макроуровня, о которых будет говориться ниже.

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

3. Определение напряженно-деформированного состояния.

Постановка задачи на макроуровне II

На макроуровне II рассматривается представительный объем материала, состоящего из совокупности кристаллитов — элементов мезо-уровня. Конститутивная модель материала на рассматриваемом уровне принимается в виде

2сЯ = 2 + 2 • а + ат • 2 = П :

где М — число элементов мезоуровня, необходимых для описания представительного объема макроуровня; г — индекс элемента мезоуровня (далее опускается). На макроуровне явными внутренними переменными [9] являются неупругая составляющая градиента относительной скорости перемещений Zгr”, тензор спина П, тензор упругих

свойств П и тензор термического расширения А. В качестве неявных внутренних переменных, входящих в замыкающие уравнения (15)—(18), используются неупругая составляющая градиента относительной скорости перемещений £гП, тензор спина ю, тензор упругих свойств п и тензор термического расширения а отдельных кристаллитов. Роль эволюционных уравнений, в которые вышеперечисленные характеристики кристаллитов входят в качестве переменных, играют соотношения мезоуровня.

4. Определение напряженно-деформированного состояния.

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

где о, осг — тензор напряжений Коши и его коротационная производная; £г = £ - ю = уУ - ю = угV — транспонированный градиент относительной скорости перемещений (индексом г обозначены величины, характеризующие относительное движение, фиксируемое подвижным наблюдателем в жесткой подвижной системе отсчета мезоуровня); у — скорость перемещений.

Постановка задачи на мезоуровне

о °г = о - ш • о + о • ш = п :

Неупругая составляющая градиента относительной скорости перемещений определяется двумя составляющими: пластической и трансформационной:

N п

?Г = е 4, Е р,“ М' ’ + !: 4, С',

1 '=1 1

И ,

где 4, = V - объемная доля соответствующей фазы (|Гг| - объем, занимаемый ,-й фазой, |К| - объем области, занимаемый отдельным зерном); С* - градиент трансформационной деформации, описывающий превращение начальной фазы в ,-ю; N - количество фаз, испытывающих пластическую деформацию; Ы2 - количество фаз, которые могут

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

и термодинамические потенциалы, отличные друг от друга); у(' ) - скорость сдвига по '-й системе скольжения ,-й фазы (принято удвоенное количество систем скольжения, т.е. системы скольжения с положительным и отрицательным направлениями сдвига считаются различными); р') = Ь('^п(') - ориентационный тензор '-й системы скольже-

ния; п, - единичный вектор нормали к плоскости скольжения; Ь, -

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

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

N

п = Е4,-п,-,

, =1

N

а = 14, а,,

,=1

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

что тензоры пг. являются симметричными по парам индексов uijki = пkliJ, но в общем случае не являются симметричными внутри пар

индексов: пт * пjm, пу1к * пуы

На мезоуровне явными внутренними переменными являются неупругая составляющая градиента относительной скорости перемещений ZГ, тензор спина ю, тензор упругих свойств п и тензор термического расширения а  В качестве неявных внутренних переменных вы-

ступают скорости сдвига по системам скольжения у; ', сопротивления сдвигу тк, ориентационные тензоры o, объемные доли сосуществующих в зерне фаз и скорости их изменения 4г-, для которых записываются эволюционные уравнения. Уравнения для пластической и трансформационной составляющих градиента относительной скорости перемещений являются замыкающими для скоростей сдвигов и скоростей изменения объемных долей фаз в рамках элемента мезо-уровня. К группе замыкающих уравнений относятся также соотношения для определения тензора упругих свойств и тензора термического расширения в многофазной системе, спина решетки.

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

/

(к) / т

у

Т

где т(к^ - критическое напряжение сдвига по k-й системе скольжения;

c г

•

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

Сопротивление сдвигу т1' в каждый момент процесса определяется следующими соотношениями:

-

& (к)_

=уД

а

У

1=1

1 г

V 1=1

(Л

у

- Д 0

Гт

\т

—1

Ёу‘Л, ¥,>1, У(1)^0,

1

С! 0’

т

С! 0

где Тс, о - начальный предел текучести в к-й системе скольжения ,-й

фазы; Е1 - параметр, который играет роль нормирующего множителя, в качестве которого можно выбрать любой из диагональных компонент тензора упругих свойств; у,, Б1, ш1, тс. - материальные константы; а^ - модули упрочнения, которые обычно разделяют на два вида:

1 г

первый (1=к) определяет влияние на критическое напряжение в данной системе скольжения, предшествующей деформации, по этой же системе скольжения (коэффициенты деформационного упрочнения), а второй Цфк) - на движение дислокаций в данной системе скольжения, которое оказывают дислокационные процессы в других системах (коэффициенты латентного упрочнения); у( 1 - накопленный сдвиг по

системе скольжения 1, у(1) = | у(1 т. Физические причины, приводящие

о

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

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

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

і=1

где И1 характеризует скорость упрочнения аустенита из-за увеличения

объемной доли новой фазы; N - количество фаз, которые могут образоваться при фазовом превращении (например, различные варианты мартенсита, которые будут оказывать влияние на сопротивление сдвигу в аустените).

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

т

с , с , 1 с А 4 у

т^)’1 - сопротивление сдвигу в новой 1-й фазе (для мартенсита - это

различные варианты) с учетом влияния предварительной пластической

деформации аустенита; т^ ) - критическое напряжение сдвига аусте-

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

Изменение ориентации зерна определяется по его спину в результате интегрирования соотношения

ш = о•от .

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

1 N1 п . .

ш = »’ - ^Е-е’Сп'1’ - п!к)ь!к>),

2 2=1 к=1

где w* = -2-! уУ-Уу ) = 1! ? - ^т).

Экспериментально установлено, что для металлов существенную роль в поворотах решетки играет несовместность скольжения дислокаций в соседних зернах, поэтому модель для описания поворотов решетки можно усложнить за счет введения дополнительной составляющей, характеризующей ротацию собственно решетки кристаллита, обусловленную взаимодействием с окружением, для описания которого вводится внутренняя переменная - действующее на зерно момент-ное напряжение [10].

Изменение объемной доли описывается следующим кинетическим соотношением:

0 84, 0 54;,

где ¥ - свободная энергия Гиббса, - кинетические коэффициенты,

0 - температура. Вопрос получения кинетического уравнения (29) в рамках термодинамики необратимых процессов рассматривается ниже.

5. Согласование определяющих соотношений для задачи

определения напряженно-деформированного состояния

Согласование определяющих соотношений для задачи определения напряженно-деформированного состояние основано на подходе, предложенном в работе [10]. Основная идея данного подхода заключается в том, что при использовании на мезоуровне и макроуровне II однотипных определяющих соотношений (в рассматриваемой работе это закон Гука в скоростной релаксационной форме) определяющее соотношение на макроуровне II получается осреднением определяющего соотношения мезоуровня при наложении априорных связей между однотипными параметрами и установлении связей по остальным параметрам из сопоставления определяющего соотношения макроуровня II и осредненного определяющего соотношения мезоуровня. На основе использования данного подхода получаются соотношения для определения явных внутренних переменных в определяющем соотношении макроуровня II через неявные внутренние переменные и параметры состояния системы, определяемые из соотношений мезоуровня:

а и

г=<?>>

2 = ю ;< о >,

а =

А =

0 =

Ъгп = < £ ш > ^а/д^ + п-1 ;< § . п/ . ш/ + ю/ + а/& ) > _

-П-1:( <ю/ • о/ >-<о/ • ю/ >),

где символ означает осреднение соответствующих величин при пе-

реходе на макроуровень II; символ «'» обозначает флуктуации величин по отношению к своим средним значениям по представительному объ-

рый на произвольный тензор второго ранга Т действует следующим

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

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

к0-У-( X-У0) = у, е( і = О)=0о,

к - теплоемкость; X - тензор теплопроводности; у - мощность внутреннего источника тепла; 90 - значение температуры в начальный момент времени.

Значения температуры и ее градиента определяются суммами средних по объему макроуровня величинами и некоторыми флуктуациями [15]:

9 =

У9 =

где оператор осреднения ^ обладает следующим свойством:

Принимается, что средние значения равны значениям температуры и ее градиента на макроуровне I:

Если воспользоваться простейшим способом осреднения - осреднением по объему, то условие (41) можно записать в следующем виде:

|у9'ёК =[9' пёГ = 0.

V Г

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

9 + =9' ^9 +-9 = У©

Г+ Г- Г+ Г- \Г+ Г4 '

знаки «+» и «-» соответствуют противоположным границам.

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

N

X = 14, X, ,

э, г i=1

N

К = 14гК- ,

ог г г=1

где X ,, к? - значения соответствующих величин в ,-й фазе.

Мощность внутреннего теплового источника в зерне определяется как

N1 п N2

У = & )У

г=1 к=1 г=1

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

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

К = <к>.

Л =

г=< у) .

7. Постановка задачи диффузии для точки интегрирования

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

15' = У {I5, К • К - « ?е)),

ср

ср = сгр на Гс,

где ср, ср, с^ - концентрация р-го легирующего элемента, ее значение

в начальный момент времени и на границе; Ар, d р - тензоры диффузии

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

8. Кинетическое уравнение для изменения объемной доли в рамках термодинамики необратимых процессов

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

1-#

Цу - положительный кинетический параметр, связанный со скоростью

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

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

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

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

где е - внутренняя энергия.

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

Принцип Онзагера постулирует линейную релаксацию средних значений экстенсивных термодинамических переменных к их равновесным значениям. Если обозначить отклонения этих величин от их равновесных значений через а,, то поведение а, во времени описывается линейным феноменологическим уравнением [5]

где = —г— средний термодинамический поток; X, - термодинами-

&

ческая сила, а кинетические коэффициенты Ц удовлетворяют соотношениям взаимности Онзагера (Ц = ).

2,..., с1, с2,...

ёа

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

^ = -2g..ага. ,

2 , ]

gj - компоненты положительно определенной матрицы. При этом соотношения взаимности Онзагера имеют место, если термодинамические силы определены соотношением [5]

X, =§^ = -£g,аJ .

да, т

Заметим, что производная по времени от энтропии (57) (производство энтропии) может быть представлена в виде

^=-!^л £=!«.

т. е. описывается билинейным выражением относительно потоков

и термодинамических сил, входящих в феноменологические уравнения, для которых справедливы соотношения Онзагера.

Таким образом, вычисление производства энтропии дает рецепт выбора потоков и сопряженных им термодинамических сил. С учетом (56) соотношение (59) можно переписать в виде

ёМ

ді

=ЕЕЬх.х, г О,

Ь - положительно определенная матрица.

В случае записи кинетического уравнения для изменения объемной доли в качестве обобщенной силы можно принять = хт, а в

дЕ .

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

д,

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

N

ЕЕ,= 1,

,=1

необходимо применить технику множителей Лагранжа.

Выражение (55) запишется следующим образом:

д,

где

lN

Ее . -1

т=1

Л

й¥,

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

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

5£

5Е

- + Х = 0, т = 1..N.

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

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

дЕ

д,

5Е, 5Е

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

ёе = 0с1^ +

С де ^ Чдс1 у

ёс +... +

4^1 У

дяе

:ёяе.

Используя связь между свободной и внутренней энергиями, е = ¥ + 0£, получаем следующее соотношение:

ё¥ = —ё0 +

С де ^ Чдс1 У

ёс +... +

ёЕ1 +... +

ч^1 У

дяе

: ёяе.

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

^1 У

и т.д.

^1 У

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

д$

(

дґ

1 1 5¥

0 5$, 0 5$

Л

І У

где X. = -^¥ -~“ - движущая сила фазового превращения у-й фазы

в ,-ю.

Если Xу > 0 и 0 < ^ <1, 0 <Е,у < 1, то происходит фазовое превращение у-й фазы в /-ю, объемную долю которой можно определить из решения уравнения (70). Если же X. > 0 и Е1/ = 1, Е,. = 0, то превращение не происходит.

8. Определение свободной энергии многофазной системы

Предположим, что удельную свободную энергию можно представить в следующей форме:

¥

где - упругая составляющая удельной свободной энергии, представляющая собой квадратичную форму от упругой составляющей меры деформационного состояния qe ; ¥сн - химическая составляющая,

N

определяемая как ¥ссн = Х^г¥ (с1,с2,...,0),¥1 - удельная свободная

/ =1

энергии для отдельной фазы; ¥5 - поверхностная энергия.

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

¥г (,с2,..., 0) = ¥гчистый + ¥гсмеш + ¥гюаим,

где ¥/ ЧИСТЫй - свободная энергия в отдельной фазе, ¥/ смеш - изменение свободной энергии за счет смешивания компонент, ¥вза™ - изменение

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

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

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

Х1р, X2гр, ХзР - константы.

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

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

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

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

К (0)=х;г+х2 ; 0+хзр 01п (0),

Я = кЫа, Ыа - число Авогадро. Для всей системы получим

п

;=1

ррсмеш = Ясмеш _ 0^смеш = ср 1п

V Р=1

где Ыт - количество молей вещества.

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

химическая связь между атомами при описании изменения их расположения не важна.

Третье слагаемое (72) описывает изменение свободной энергии за счет химического взаимодействия отдельных компонент. Будем рассматривать только парное взаимодействие и применим модель [16] для описания взаимодействия компонент друг с другом. Тогда изменение свободной энергии за счет химического взаимодействия компонент определится следующим образом:

1 п _1 п    / ч

рвзаим _нвзаим 0*взаим _ ^ ^ С*

г г i ^ ^ р ] р] \ / 5 ^ '

V р=1 ]=р+1

£,взаим = о, так как энтропия системы полностью определяется ее конфигурационной составляющей, Тр] (0) - параметр, описывающий

взаимодействие между р и] компонентами в соответствующих фазах.

Третье слагаемое в выражении (71) для свободной энергии (поверхностная энергия) в самом общем виде может быть записано следующим образом [17]:

Р* = 1хХ ,

к

где хк - энергия границы между фазами г и], занимающей область Бк.

В работе [18] предложено простое соотношение для определения поверхностной энергии, записанное в терминах объемной доли,

Р* =

г =1 1 г=1

х

где Ш = 1; I - толщина границы. Равенство (78) выведено в предположении о том, что кристаллы новой фазы представляют собой пластинки.

Заключение

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

Работа выполнена при финансовой поддержке гранта РФФИ (проект №13-01-96006-урал_а), а также гранта ФЦП «Научные и научно-педагогические кадры инновационной России на 2009-2013 годы» (мероприятие 1.2.2, Соглашение 14.В37.21.0382).

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

1. Волегов П.С., Никитюк А.С., Янц А.Ю. Геометрия поверхности текучести и законы упрочнения в физических теориях пластичности // Вестник Перм. гос. техн. ун-та. Сер. Математическое моделирование систем и процессов. - 2009. - № 17. - С. 25-33.

2. Исупова И.Л., Трусов П.В. Обзор математических моделей для описания фазовых превращений в сталях // Вестник Перм. нац. исслед. политехн. ун-та. Механика. - 2013. - № 3. - С. 157-192.

3. Лурье А.И. Нелинейная теория упругости. - М.: Наука, 1980. -

512 с.

4. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритмы, приложения. - М.: Наука, 1986. - 232 с.

5. Пригожин И., Кондепуди Д. Современная термодинамика. От тепловых двигателей до диссипативных структур / пер. с англ. Ю.А. Данилова и В.В. Белого. - М.: Мир, 2002. - 461 с.

6. Трусов П.В., Швейкин А.И. Многоуровневые физические модели моно- и поликристаллов. Статистические модели // Физическая мезомеханика. - 2011. - № 4. - С. 17-28.

7. Трусов П.В., Швейкин А.И. Многоуровневые физические модели моно- и поликристаллов. Прямые модели // Физическая мезоме-ханика. - 2011. - Т. 14. - № 4. - С. 5-30.

8. Конститутивные соотношения и их применение для описания эволюции микроструктуры / П.В. Трусов, В.Н. Ашихмин, П.С. Воле-гов, А.И. Швейкин // Физическая мезомеханика. - 2009. - Т. 12. -№ 3. - С. 61-71.

9. Трусов П.В., Ашихмин В.Н., Швейкин А.И. Двухуровневая модель упругопластического деформирования поликристаллических материалов // Механика композиционных материалов и конструкций. -2009. - Т. 15. - № 3. - С. 327-344.

10. Многоуровневые модели неупругого деформирования материалов и их применение для описания эволюции внутренней структуры / П.В. Трусов, А.И. Швейкин, Е.С. Нечаева, П.С. Волегов // Физическая мезомеханика. - 2012. - Т. 15. - № 1. - С. 33-56.

11. Трусов П.В., Нечаева Е.С., Швейкин А.И. Применение несимметричных мер напряженного и деформированного состояния при построении многоуровневых конститутивных моделей материалов // Физическая мезомеханика. - 2013. - Т. 16. - № 2. - С. 15-31.

12. Beese A.M., Mohr D. Anisotropic plasticity model coupled with lode angle dependent strain-induced transformation kinetics law. Submitted for publication, 2011.

13. Cherkaoui M., Berveiller M., Sabar H. Micromechanical modeling of martensitic transformation induced plasticity (TRIP) in austenitic single crystals // International Journal of Plasticity. - 1998. - Vol. 14. - No. 7. -P. 597-626.

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

15. Ozdemir I., Brekelmans W.A.M., Geers M.G.D. Computational homogenization for heat conduction in heterogeneous solids // International Journal for Numerical Methods in Engineering. - 2008. - Vol. 73. - No. 2. -P. 185-204.

16. Redlich O., Kister A.T. Algebraic representation of thermodynamic properties and the classification solutions // Ind. Eng. Chem. - 1948. -Vol. 40. - No. 2. - P. 345-348.

17. Wagemaker M., Mulder F.M., Van Der Ven A. The role of surface and interface energy on phase stability of nanosized insertion compounds // Adv. Mater. - 2009. - Vol. 21. - P. 2703-2709.

18. Wang J.J., Van Der Zwaag S. Stabilization mechanisms of retained austenite in transformation-induced plasticity steel // Metall. Mater. Trans. A. - 2001. - Vol. 32. - No. 6. - P. 1527-1539.

References

1. Volegov P.S., Nikityuk A.S., Yants A. Yu. Geometriya poverkh-nosti tekuchesti i zakony uprochneniya v fizicheskikh teoriyakh plastich-nosti [The geometry of the yield surface and hardening laws in physical theories of plasticity]. Vestnik Permskogo gosudarstvemmogo polytech-nicheskogo universiteta. Mekhanika, 2009, no. 17, pp. 25-33.

2. Isupova I.L., Trusov P.V. Obzor matematicheskikh modeley dlya opisaniya fazovykh prevrashcheniy v stalyakh [Review of mathematical models on phase transformations in steels]. Vestnik Permskogo natsional-nogo issledovatelskogo polytekhnicheskogo universiteta. Mekhanika, 2013, no. 3, pp.

3. Lurie A.I. Nelineynaya teoriya uprugosti [Nonlinear theory of elasticity]. Moscow: Nauka, 1980, 512 p.

4. Pozdeev A.A., Trusov P.V., Nyashin Yu.I. Bolshie uprugo-plasti-cheskie deformatsii: teoriya, algoritmy, prilozheniya [Large elastic-plastic deformation: theory, algorithms and applications]. Moscow: Nauka, 1986, 232 p.

5. Kondepudi D., Prigogine I. Modern thermodynamics: from heat engines to dissipative structures. John Wiley & Sons, 1998, 508 p.

6. Trusov P.V., Shveykin A.I. Mnogourovnevye fizicheskie modeli mono- i polikristallov. Statisticheskie modeli [Multilevel physical models of single- and polycrystals. Statistical models]. Fizicheskaya mezomekhanika, 2011, no. 4, pp. 17-28.

7. Trusov P.V., Shveykin A.I. Mnogourovnevye fizicheskie modeli mono- i polikristallov. Pryamye modeli [Multilevel physical models of single- and polycrystals. Direct models]. Fizicheskaya mezomekhanika, 2011, vol. 14, no. 4, pp. 5-30.

8. Trusov P.V., Ashikhmin V.N., Volegov P.S., Shveykin A.I. Konsti-tutivnye sootnosheniya i ikh primenenie dlya opisaniya evolyutsii mik-rostructury [Constitutive relations and their application to the description of microstructure evolution]. Fizicheskaya mezomekhanika, 2009, vol. 12, no. 3, pp. 61-71.

9. Trusov P.V., Ashikhmin V.N., Shveykin A.I. Dvukhurovnevaya model uprugoplasticheskogo deformirovaniya polikristallicheskikh materia-lov [Two-level model for polycrystalline materials elastoplastic deformation]. Mekhanika kompozitsionnykh materialov i konstruktsiy, 2009, vol. 15, no. 3, pp. 327-344.

10. Trusov P.V., Shveykin A.I., Nechaeva E.S., Volegov P.S. Mnogourovnevye modeli neuprugogo deformirovaniya materialov i ikh primenenie dlya opisaniya vnutrenney struktury [Multilevel models of inelastic deformation of material and their application for description of internal structure]. Fizicheskaya mezomekhanika, 2012, vol. 15, no. 1, pp. 33-56.

11. Trusov P.V., Nechaeva E.S., Shveykin A.I. Primenenie nesimmet-richnykh mer napryazhennogo i deformirovannogo sostoyaniya pri postroe-nii mnogourovnevykh konstitutivnykh modeley materialov [Non-symmetric stress-strain measures using when construct multilevel constitutive material models]. Fizicheskaya mezomekhanika, 2013, vol. 16, no. 2, pp. 15-31.

12. Beese A.M., Mohr D. Anisotropic plasticity model coupled with lode angle dependent strain-induced transformation kinetics law. Submitted for publication, 2011.

13. Cherkaoui M., Berveiller M., Sabar H. Micromechanical modeling of martensitic transformation induced plasticity (TRIP) in austenitic single crystals. International Journal of Plasticity, 1998, vol. 14, no. 7, pp. 597-626.

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

15. Ozdemir I., Brekelmans W.A.M., Geers M.G.D. Computational homogenization for heat conduction in heterogeneous solids. International Journal for Numerical Methods in Engineering, 2008, vol. 73, no. 2, pp.185-204.

16. Redlich O., Kister A.T. Algebraic representation of thermodynamic properties and the classification solutions. Ind. Eng. Chem., 1948, vol. 40, no. 2, pp. 345-348.

17. Wagemaker M., Mulder F.M., Van der Ven A. The role of surface and interface energy on phase stability of nanosized insertion compounds. Adv. Mater., 2009, vol. 21, pp. 2703-2709.

18. Wang J.J., Van Der Zwaag S. Stabilization mechanisms of retained austenite in transformation-induced plasticity steel. Metall. Mater. Trans. A, 2001, vol. 32, no. 6, pp. 1527-1539.

Об авторах

Исупова Ирина Леонидовна (Пермь, Россия) - аспирант кафедры математического моделирования систем и процессов Пермского национального исследовательского политехнического универ-

ситета (614990, г. Пермь, Комсомольский пр., 29, e-mail: enoty-forever@yandex. ru).

Трусов Петр Валентинович (Пермь, Россия) - доктор физикоматематических наук, профессор, заведующей кафедрой математического моделирования систем и процессов Пермского национального исследовательского политехнического университета (614990, г. Пермь, Комсомольский пр., 29, e-mail: tpv@matmod.pstu.ac.ru).

About the authors

Isupova Irina Leonidovna (Perm, Russian Federation) - postgraduate student of Department of Mathematical Modeling of Systems and Processes, Perm National Research Polytechnic University (29, Komsomolsky av., 614990, Perm, Russian Federation, e-mail: enotyforever@yandex.ru).

Trusov Peter Valentinovich (Perm, Russian Federation) - Doctor of Physical and Mathematical Sciences, Professor, Head of Department of Mathematical Modeling of Systems and Processes, Perm National Research Polytechnic University (29, Komsomolsky av., 614990, Perm, Russian Federation, e-mail: tpv@matmod.pstu.ac.ru).

Получено 21.07.2013