О численном решении задач неизотермической многокомпонентной диффузии с переменными коэффициентами Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 21, № 2, 2016

О численном решении задач неизотермической многокомпонентной диффузии с переменными коэффициентами

С. А. Шанин*, А. Г. Князева

Национальный исследовательский Томский политехнический университет, Россия *Контактный e-mail: shanin_s@mail.ru

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

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

Введение

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

Так, в [1] используется неявный метод для интегрирования системы уравнений многокомпонентной диффузии, основанный на преобразовании системы к диагональному виду. Этот метод аналогичен описанному в [2]. В работе [3] рассмотрено применение

© ИВТ СО РАН, 2016

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

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

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

И, наконец, в [7] используется разработанный в настоящее время метод, основанный на явном задании выражения для свободной энергии. Соответствующие модели включают целый спектр параметров, значения которых неизвестны. В работе используются два различных численных метода: метод спектральных коллокаций и метод конечных элементов. Как правило, авторы подобных работ редко обсуждают корректность задач, однозначность выбора параметров и не обращают внимания на отрицательные концентрации и законы сохранения.

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

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

1. Описание проблемы

Систему уравнений теории многокомпонентной диффузии можно представить в виде

п— 1

^ V (Ду С)) + Д (С), к = 1, 2 ,...,П - 1,

дСк п 1

дг . 1

3 = 1

п

Е ^ = 1. к=1

В декартовой системе координат имеем одномерные уравнения

^ = Ё (С) Щ + А (С), * = 1,-1, (!)

.7 = 1 4 7

где Ск = (С1,... , Ст) — концентрации компонентов; (С) — парциальные коэффициенты диффузии и Д (С) — источники и стоки компонентов в химических реакциях (гладкие функции своих переменных).

В классической теории коэффициенты диффузии D^j — функции концентраций. Вопрос о корректности задачи теории многокомпонентной диффузии ставился в литературе неоднократно [2, 9, 10].

В неизотермических условиях к числу параметров состояния добавляется температура, а при учете напряженно-деформированного состояния — инварианты тензоров либо напряжений, либо деформаций. Для тел канонической формы напряженно-деформированное состояние иногда удается явно связать с изменением температуры и состава. Но дальнейшие преобразования приводят к появлению особенностей в системе эффективных коэффициентов переноса, которые не удается устранить такими жесткими ограничениями, как неотрицательность матрицы коэффициентов и всех ее главных миноров [10, 11]. Проблема усложняется, когда речь идет о диффузии в анизотропных средах.

Условия неотрицательности решений соответствующих задач записываются следующим образом [2, 9]:

1. Dkj(С) = 0 при Ск = 0, Сг > 0 (j = к, к = j).

2. Dkk(C) > 0 при Ск = 0, Сг > 0 (j = к). D\\ D\2 s D\m

3.

4.

Doi D

22

D

2m

Dm\ D.

m2

s Dn

> 0.

Dn D\2

D21 D22

> 0,

D22 D23 D32 D33

0,

Dm—I,m—1 Dm-D.

m— 1,m

m,m— 1

Dr,

> 0.

Эти условия позволяют предварительно оценить возможность появления отрицательного решения. Если нелинейные коэффициенты системы (1) есть следствие решения задачи о механическом равновесии, то нарушение условий 1 и 2 может быть следствием иного механизма переноса, завуалированного в коэффициентах. Пример такого явления — бародиффузия или диффузия под действием напряжений [12]. Его следствием могут быть так называемая восходящая диффузия, хорошо известная в физике [13, 14], наличие пределов растворимости, зависящих от состава и напряжений и т. д.

Если для коэффициентов, описывающих перекрестные эффекты, смена знака — дело обычное, то для диагональных коэффициентов в этом случае возникает неопределенность типа [0 • то],

дС

lim Dk ^ 0 •то, (2)

ск ^с* дх

что приводит к проблемам при численном решении.

S

2. Особенность алгоритма

Учитывая особенности процесса диффузии (наличие различных механизмов переноса, перекрестных потоков и др.), можно предложить численный алгоритм, раскрывая неопределенность вида (2) явно. Критические значения концентраций С*, при которых возникают неопределенности, можно найти из соотношения Dk(Cj) = 0, что, возможно, говорит о достижении предела растворимости [12]. Поскольку концентрация есть функция двух переменных, Ск(Ъ,х), то, определив С*, можно найти координату точки, где достигнут предел растворимости х* в момент времени

Раскроем неопределенность в (2) по правилу Лопиталя, при этом перейдем к пределу по координате

дС-

шъ (С )СЛс

Иш

дС3/дх

дх 1/Бз (С)

-И

д2С

дх2

> щ дС1

дС3 дх

ддс,

д х2

где

Ъ

И

дС3 дС3 дВ3 дх

Следовательно, уравнение диффузии в точке, где достигнут предел растворимости по элементу С к, принимает вид

дс

дг

д_ д х

( д2СЛ

V' дх2 )

+ Р(С,х, I),

(4)

с^ с*

где

Р (С,х, *) = У-Б- И

у д[вгкЩ

дх V гдх ) к=1 4 у

+ Л-( С, х, I).

к=г

Уравнение (4) совпадает по форме с уравнением Кортевега — де Вриза, которое описывает динамику уединенных волн (так называемых солитонов). Но это уравнение справедливо только в точке, где Сj ^ С*, и свидетельствует о подключении механизма переноса, отличного от диффузии.

Допустим, что новое уравнение имеет место в некоторой области изменения Сj Е [ С* — А С*; С* + АС*], и построим алгоритм численного решения задачи, сходимость которого будет зависеть от величин А С*. Для каждого элемента предел С* может быть достигнут в разных х* ив разные моменты времени. Для этого перейдем от (4) к уравнению

дС^

т

4х ) §

Ш (С,х, О,

(5)

где

а = 1, 3 = 0;

а = 0,

3 = 1;

С, Е [С* — АС*; С* + АС*],

С, Е [С* — АС*; С* + АС*].

(7)

Уравнения вида (5) решаются численно, с помощью пятиточечной прогонки [15]. Алгоритм позволяет избежать появления отрицательных концентраций, что продемонстрируем на примере одномерной задачи о перераспределении элементов в четырехком-понентной системе. В начальный момент времени ( = 0) имеем

С1 = 0, С2 = 0, Сз

{

0.2, 0,

г < До, г > По.

2

*

х—^-х

С 0.2

0.1

0.0

-0.1

-0.21—т-т-т-т--—г-т-т-т-

0.0 2.8 5.6 8.4 r-10"J 0.0 2.8 5.6 8.4 г-10"3

Рис. 1. Распределение концентрации С2 по координате в различные моменты времени, рассчитанные с учетом (б) и без учета (а) разработанного метода: при t = 300, 350, 400, 450 (кривые 1-4 соответственно)

В любой момент времени С4 = 1 — С\ — С2 — С3. На границах интервала [0,1] зададим условия вида

ж = 0 : Jx =0; С2 = 0.2; С3 = 0.2;

ж = L : Jk = 0; к = 1, 2,3.

Диагональные коэффициенты диффузии считаем зависимыми от концентрации линейно

Dkk = Dl0k [1 — Bk Ск ], Dki = const,

где D01 = 1.2 • 10-5 м2/с; Dx2 = 3 • 10-5 м2/с; Dx3 = 2 • 10-5 м2/с; D2X = 5 • 10-5 м2/с; D02 = 1 • 10-5 м2/с; D23 = 1 • 10-5 м2/с; D31 = 1 • 10-5 м2/с; D32 = 2 • 10-5 м2/с; D03 = 1 • 10-5 м2/с; В1 = 5; В2 = 5; В3 = 5; R0 = 3 • 10-3. Если Ск = 1/Вк, коэффициенты Dkk меняют знак.

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

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

3. Задача о росте покрытия

В качестве второго примера рассмотрим задачу о росте покрытия при осаждении веществ из плазмы для образца цилиндрической формы, который вращается вокруг своей оси. На поверхности полого железного цилиндра (рис. 2) (с внутренним и внешним радиусами R1 = 0.01 ми R2 = 0.02 м) формируется покрытие из хрома (Cr) и углерода (С), поступающих из плазмы. Элементы вступают в химические реакции, при этом

Рис. 2. К постановке задачи о росте покрытия

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

рс-

дТ

Ж

+ £

<Рг

г=1

дСк

-У Л + Гк, к = 1, 2,

Р-

дСк

дИ

гк, к = 4, 5,

где

31 = -р^пУС! - - рИцО^т {УТ,

Л = -РВ21УС1 - рВ22УС2 - рВ^С^т2УТ, Л = -СЛ- А2Х7С-

12 У ^2,

Т — температура; Ск — концентрации компонентов (С1 — Сг; С2 — С; С3 — Ее; С4 —

г

Сг3С2 — тугоплавкий карбид хрома; С5 — ГеС); г к = икг^к 'рг — источники и стоки

г=1

п

массы вследствие химических реакций; Q1 = - ^ тк — тепловые эффекты реак-

к=1

ций, число которых г; — стехиометрические коэффициенты компонентов к в реакциях г; тк — молярные массы компонентов; Ъ^ — парциальные энтальпии компонентов; >рг — скорости химических реакций, которые зависят от температуры по закону Ар-рениуса (например, >р1 = к1С3С2 ехр[-Е1/КГ]); с — эффективная теплоемкость; р — плотность; \т ,Игк ,А1,А2 — коэффициенты переноса; втк — коэффициенты Соре.

Из трех диффузионных потоков (подвижны только тугоплавкий металл и углерод)

п

независимы только два, и из пяти концентраций независимы только четыре ( ^ Ск = 1

\к=1

так что для концентраций железа имеем С3 = 1 - С1 - С2 - С4 - С5.

Плотность, коэффициент теплопроводности и теплоемкость определяем по правилу смеси:

Р(С) = £ Рзсз, с(С) = £ сС<, Х(С) = £ Л,С:

=1

=1

=1

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

г = : ^ = 0, = 0, = 0.

На границе раздела материалов (подложки и растущего покрытия) равны друг другу потоки тепла и температуры, а также потоки массы и химические потенциалы углерода и хрома.

г = Къ : Jq\К2- = Jq+, Т\К2- = Т+, \д2 - = \я2+, 91 \в,2 - = 91 \в,2+, ^2\Я2 - = \Я2 +, 92\Д2 - = 92\Я2 +.

Условия для химических потенциалов сводятся к условиям вида

С- = С+

С- =

2 12,

где h ,12 — коэффициенты распределения. Условия на подвижной поверхности имеют вид

7~> - т т dh frr.4 т dh dh

г = R2 + h : Jq = q0 — - ооe(T - Tw), Jx = mxyx — , J2 = m2y2 — ,

dh w rr,

где —--заданная скорость роста покрытия; q0 — суммарный источник тепла; ±w —

at

температура стенок вакуумной камеры; а0 — постоянная Стефана —Больцмана; е — степень черноты; у^ — мольная концентрация частиц у поверхности, моль/м3, которые следуют из эксперимента либо из решения внешней задачи. Для исследования процессов в покрытии и в подложке эти величины могут быть заданы.

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

Т(г, 0) = То, С1 (г, 0) = 0, С2 (г, 0) = 0, Сз (г, 0) = 1.

С 0.90.60.30.0-

г 10 ,м

1.9995 1.9998 2.0001 2.0004

1

а

Рис. 3. Распределение концентраций С, Сг, Ре (кривые 1-3 соответственно) вдоль радиуса в момент времени £ = 1200 с; сплошные кривые — традиционный алгоритм, точки — предложенный алгоритм

Более подробно модель описана в [16, 17]. В расчетах принято: Ае = 2 Вт/(м • К); ÀCr = 93.9Вт/(м • К); ÀFe = 80.1 Вт/(м • К); ЛСгзСз = 52,1 Вт/(м • К); ÀFeC = 60.1 Вт/(м • К); се = 20.8 Дж/(кг^К); = 21.1 Дж/(кг^К); CFe = 25.14 Дж/(кг-К); СсГзС2 = 35.4 Дж/(кг К); CFeC = 29.5 Дж/(кг • К); ре = 2.55 • 103 кг/м3; = 7.19 • 103 кг/м3; pFe = 7.87 • 103 кг/м3; рСгзСз = 6.68 • 103 кг/м3; pFeC = 7.6 • 103 кг/м3; тс = 12.01 • 10-3 кг/моль; mCr = 51.9 • 10-3 кг/моль; mFe = 55.84 • 10-3 кг/моль; тСтзС2 = 180 • 10-3 кг/моль; mFeC = 68 • 10-3 кг/моль; ST1 = ST2 = 10-2 K-1.

В качестве примера на рис. 3 показано распределение концентраций элементов в момент времени t = 1200 c при использовании традиционного алгоритма (уравнения параболического типа, неявная разностная схема и трехточечная прогонка) и предложенного в данной работе алгоритма (точечные кривые). На первый взгляд результаты не различаются (рис. 3, а). Но в окрестности границы раздела подложка — растущее покрытие имеется особенность (рис. 3, б). Отрицательные концентрации исчезают при учете в алгоритме физической особенности — появления потоков компонентов в направлении градиента концентрации вследствие термодиффузии, наличия градиентов других компонентов, что приводит к смене знака эффективного коэффициента диффузии.

Выводы

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

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

3. Проиллюстрирована применимость алгоритма к решению задачи формирования фазового состава растущего покрытия.

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

Благодарности. Работа выполнена в рамках Государственного задания № 11.815.2014/К. Список литературы / References

[1] Wangard, W., Dandy, D.S., Muller, J. A numerically stable method for integration of the multicomponent species diffusion equations // Journal of Computational Physics. 2001. No. 174. P. 460-472.

[2] Многокомпонентная диффузия в гетерогенных сплавах / Л.Г. Ворошнин, П.А. Витязь, А.Х. Насыбулин, Б.М. Хусид. Минск: Вышэйшая школа, 1984.

Multicomponent diffusion in heterogeneous alloys / L.G. Voroshnin, P.A. Vityaz, A.H. Nasi-bulin, B.M. Husid. Minsk: Vysheyshaya Shkola, 1984. 144 с. (In Russ.)

[3] Kozeschnik, E. Multicomponent diffusion simulation bused on finite elements // Metallurgical and Materials Transactions. 1999. Vol. 30A. P. 2576-2582.

[4] Jungel, A., Stelzer, V. Existence analysis of Maxwell — Stefan systems for multicomponent mixtures. Available at: http://arxiv.org/abs/1211.2394V1

[5] Giovangigli, V. Mass conservation and singular multicomponent diffusion algorithms // Impact of Computing in Science and Engineering. 1990. Vol. 2. P. 73-97.

[6] Thomas, H.R., Jedighi, M., Vardon, P.J. Diffusive reactive transport of multicomponent chemicals under coupled thermal, hydraulic, chemical and mechanical conditions // Geotechnical & Geological Engineering. 2012. Vol. 30. P. 841-857.

[7] Emmanuel, S., Cortis, A., Berkowitz, B. Diffusion in multicomponent systems: a free energy aproach // Chemical Physics. 2004. No. 302. P. 21-30.

[8] Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1979. 288 с.

Tihonov, A.N., Arsenin, V.Y. Methods for solution of incorrect problems. Moscow: Nauka, 1979. 288 p. (In Russ.)

[9] Вольперт А.И., Посвянский В.С. О положительности решения уравнений многокомпонентной диффузии и химической кинетики // Хим. физика. 1984. Т. 3, № 8. С. 12001205.

Volpert, A.I., Posvyanskiy, V.S. On the positive solution of multi-component diffusion and chemical kinetics // Chemical Physics. 1984. Vol. 3, No. 8. P. 1200-1205. (In Russ.)

[10] Старк Дж.П. Диффузия в твердых телах: Пер. c англ. / Под ред. Л.И. Трусова. М.: Энергия, 1980. 240 с.

Stark, J.P. Solid state diffusion. New York; London: Willey — Interscience Publ., 1976.

[11] Хаазе Р. Термодинамика необратимых процессов: Пер. с нем. / Под ред. А.В. Лыкова. М.: Мир, 1967. 544 с.

Haase, R. Thermodynamics of irreversible processes. Massachussets: Addison-Wesley, 1969.

[12] Еремеев B.C. Диффузия и напряжения. М.: Энергоатомиздат, 1984. 240 с. Eremin, V.S. Diffusion and strain. М.: Energoatomizdat, 1984. 240 p. (In Russ.)

[13] Гегузин Я.Е. Диффузионная зона. М.: Наука, 1979. 343 с. Geguzin, Y.E. The diffusion zone. Moscow: Nauka, 1979. 343 p. (In Russ.)

[14] Князева А.Г. Режимы развития из начального зародыша твердофазной реакции, лимитируемой диффузией // Физика горения и взрыва. 1996. Т. 32, № 4. С. 72-76. Knyazeva, A.G. Modes for development of the diffusion limited solid phase reaction from the initial nucleus // Combustion, Explosion and Shock Waves. 1996. Vol. 32, No. 4. P. 72-76. (In Russ.)

[15] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 591 с.

Samarskiy, A.A., Nikolayev, E.S. Methods for solving grid equations. Moscow: Nauka, 1978. 591 p. (In Russ.)

[16] Князева А.Г., Шанин С.А. Модель роста покрытия в условиях магнетронного напыления // Изв. вузов. Физика. 2010. № 1. C. 76-81.

Knyazeva, A.G., Shanin, S.A. Model for the coating growth under the conditions of magnetron sputtering deposition // Russian Physics Journal. 2010. No. 1. P. 76-81. (In Russ.)

[17] Шанин С.А., Князева А.Г., Поболь И.Л., Дениженко А.Г. Численное и экспериментальное исследование влияния технологических параметров на фазовый и химический состав карбидного покрытия, растущего в импульсной электродуговой плазме // Хим. физика и мезоскопия. 2012. Т. 14, № 4. C. 525-535.

Shanin, S.A., Knyazeva, A.G., Pobol, I.L., Denigenko, A.G. Numerical and experimental study of the influence of process parameters on the phase and chemical composition of the carbide coating, growing in a pulsed arc plasma // Chemical Physics and Mesoscopy. 2012. Vol. 14, No. 4. P. 525-535. (In Russ.)

Поступила в 'редакцию 7 июля 2014 г., с доработки —12 января 2016 г.

On the numerical solution of non-isothermal multicomponent diffusion with variable coefficients

Shanin, Sergey A., Knyazeva, Anna G.

National Research Tomsk Polytechnic University, Tomsk, 634050, Russia * Corresponding author: Shanin, Sergey A., e-mail: shanin_s@mail.ru

Purpose. The numerical implementation of the coupling models of multicomponent diffusion requires the development of special numerical algorithms which take into account the interrelation between several physical phenomena with different temporal and spatial scales. Traditional algorithms for numerical implementation for similar models are unstable and lead to the negative concentrations. The main object of this work is to develop an algorithm for solving the equation system of multicomponent diffusion with coefficients that could change sign thus leading to the ill-posed problem.

Methodology. The basic idea of the suggested algorithm consists in the determination of critical concentration coinciding which the solubility limit. Mathematically, the concentration equaling to critical one leads to the uncertainty of the type "zero to infinity". Using the L'Hospital rule we obtain the diffusion equation, which correctly describes the diffusion in the vicinity of solubility limit and eliminate the appearance of negative values for concentrations.

Findings. As a result, the developed algorithm for the numerical investigation of the problem was implemented numerically. It was shown that the separation of the feature leads to the replacement of the parabolic diffusion equation with the equation of the third order that is correct in a small vicinity of a critical concentration. The examples of numerical solution are presented for the problem of the evolution of coating composition during deposition.

Keywords : multicomponent diffusion, cross effects, associated model.

Acknowledgements. The work was supported by the Ministry of Science and Education of Russia in 2014-2016 No. 11.815.2014/K.

Received 7 July 2014 Received in revised form 12 January 2016

© ICT SB RAS, 2016