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

Известия Тульского государственного университета Естественные науки. 2013. Вып. 1. С. 202-214

ФизикА

УДК 519.63

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

В. В. Жигунов, А. И. Лавит

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

Ключевые слова: диффузия, фазы, задача Стефана, схема Кранка-Николсон, метод Галеркина, метод конечных элементов.

Процессы диффузии часто сопровождаются фазовыми переходами. Фазы могут появляться и исчезать; их взаимное расположение, вообще говоря, неизвестно и изменяется со временем. Задача диффузии для случая нескольких фаз, притом, что положение межфазных границ неизвестно и определяется в процессе решения, называется задачей Стефана. Аналитические решения задачи Стефана [1] охватывают лишь малую часть диапазона теоретически и практически важных постановок задачи. Поэтому велика значимость численных решений и, следовательно, эффективных численных методов. Им посвящено большое число публикаций; важнейшие из них обсуждены в монографиях [1-3]. Известные методы решения разделяются на две большие группы [3]. В методах первой группы, называемых методами с выделением межфазных границ, точно удовлетворяются граничные условия на межфазных границах — условия Стефана. Эти методы сложны и становятся практически неприменимыми, когда фаз несколько, возможно их возникновение и исчезновение и (или) размерность задачи больше единицы. Вторую группу образуют методы, называемые методами сквозного счета.

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

Поэтому представляет интерес метод, сочетающий точность методов первой группы и простоту второй. Попытка разработки такого метода представлена в настоящей работе. Он основан на обобщенном методе Галеркина [4], достоинством которого является то, что в нем, как и в методе Ритца, координатные функции не обязаны удовлетворять естественным граничным условиям. Для интегрирования по времени выбрана широко применяемая конечноразностная схема Кранка-Николсон [5], координатными функциями служат обобщенные функции с компактными носителями — конечными элементами [6].

Процесс диффузии в многофазной среде описывается следующим образом. Пусть Ь — время, х^, к = 1,..., N — совокупность координат точки, в общем случае криволинейных, N — размерность задачи. Искомая функция

— концентрация диффундирующего элемента с (Ь, Xk) — определяется из решения уравнения диффузии [7]

где V — оператор Гамильтона, D (с) — коэффициент диффузии, представляющий собой некоторую известную функцию концентрации. Точкой обозначено скалярное произведение векторов.

Сформулируем условия на межфазных границах — условия Стефана. Межфазная граница — это некоторая гладкая непересекающаяся поверхность Q. Пусть M — некоторая точка на этой поверхности и n — единичная нормаль к поверхности в точке M. Направление нормали выберем так, чтобы она была направлена в сторону фазы с большим значением концентрации. Межфазная поверхность перемещается. Пусть ее скорость в точке M равна v. Обозначим область, лежащую впереди поверхности раздела (то есть там, куда указывает вектор n) знаком плюс, а область, лежащую позади этой поверхности, — знаком минус. При переходе через поверхность концентрация меняется скачком: с величины с- до величины с+. Закон сохранения массы при переходе через межфазную границу записывается в виде [7]

(i)

(с+ — c-) n • v = —n • (yDVc\+ — DVc\^j .

(2)

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

Помимо условий на межфазных границах функция с должна удовлетворять граничным условиям на поверхности тела. Пусть тело, в котором протекает изучаемый процесс, занимает область V, ограниченную непере-секающейся поверхностью 5 с единичной внешней нормалью N. На части поверхности 5, обозначаемой 5Р, заданы главные граничные условия, а на остальной ее части, обозначаемой 5п, заданы естественные граничные условия [8]

где ср (£, Xk), ,]п (Ь, Xk) — заданные функции. Последняя из них представляет собой проекцию диффузионного потока J = -В~Чс на нормаль к поверхности тела.

Должно выполняться также начальное условие

где с8 (хк) — заданная функция координат.

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

Пусть за время ДЬ межфазная граница переместится вдоль оси абсцисс на величину Дх. Тогда приращение с в точке х = а составит Дс = с- — с+ +. Очевидно, что отношение Дх/ДЬ представляет собой проекцию скорости межфазной границы на нормаль п. Таким образом, в пределе при Дх ^ 0 уравнение (4) переходит в условие на межфазной границе (2). Это значит, что искомое разрывное решение можно рассматривать как предельную точку пространства непрерывных решений уравнения (1).

5 = Бр и Бп, хк Є Бр : с = ср (і,хк), хк Є : N -Ус = — Зп ), (3)

Ь = 0 : с = с8 (хк),

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

Рис. 2. Зависимости коэффициента диффузии от концентрации: а — реальная, имеющая разрыв на межфазной границе;

Ь — аппроксимация (непрерывная зависимость); 1 — линейная аппроксимация; 2 — кусочно-постоянная аппроксимация

При построении непрерывного аппроксимирующего решения функция D (с) определяется неоднозначно. В интервале (с-, с+) ее необходимо доопределить. Каждой выбранной зависимости в упомянутом интервале будет соответствовать свое решение задачи. Напрашивается применить линейную интерполяцию (1 на рис.2). Однако более последовательно не принимать какую-либо зависимость a priori, а попытаться подобрать ее так, чтобы максимально приблизиться к выполнению условия Стефана (2). При нахождении слабого решения (см. ниже) на функцию D (с) не накладывается требование непрерывности. Поэтому допустимо считать, что в интервале (с-,с+) ее значение неизменно и равно некоторой постоянной величине D*. Чтобы понять, как влияет значение D* на решение, рассмотрим две одномерные задачи, легко решаемые аналитически.

* Граничные значения D+ и D- могут совпадать.

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

dc D д2с dt дх2

0, D

(5)

Dl, х £ (0, х^),

D-, X £ (X*, ж).

На границах интервала задаются главные граничные условия

c (t, 0) — cp0, lim c (t, х) — cpж, ^, cp^ — cons^ cp0 ^ c— < c+ ^ cpж

Условия на межфазной границе записываются как

dx*

дс1

Cl\xt-0 C-, c-\xt+0 С+,‘ с* dt Dl дх

дх

x* — 0

- D2 дс-

дх

(6)

x* +0

где с* = с+ — с- .В начальный момент концентрация постоянна с3 = срте; при этом фаза 1 отсутствует и х* = 0.

Решение задачи ищется в виде [9]

Cl — Cp0 + А1Ф

с2 — cpж B1

ш

x

Ф(х) — -2= J e-2 d£,

1 - Ф

где А\ и Б\ — постоянные. При этом уравнение (5), начальное условие и граничные условия на концах интервала удовлетворяются. Подстановка этих соотношений в первые два условия (6) дает

cp0 + А1Ф

1ф(

\2^щ;

С-,

cpж B1

1-ф( -х'^=\

с+.

Чтобы эти выражения имели смысл, необходимо, чтобы функция х* (Ь) была пропорциональна у/1, например, х* = 2А\/ОГЬ, где А — некоторая постоянная. При этом условия Стефана (6) приводят к системе трех уравнений относительно неизвестных АГ, БГ и А

~cp0

B

Ф(А)

__ Срж—c+

1 — 1—Ф(ЗА)

—nc*ß\ — A1ß exp (-X2) — B1 exp (-ß2\2)

решение которой позволяет найти значения функции с (Ь,х).

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

же, что и для первой задачи. Коэффициент диффузии — кусочно-постоянная функция

{^1, с € [Ср0, С-]

О*, с € (с-, с+) •

О, с € [с+, сРж]

Так как она имеет три участка, область х € [0, то) распадается на три интервала, решения в которых имеют различный вид. В области 1 О = О1, х € [0, х13 (£)], в области 2 О = О2, х € [х32 (¿), то), в области 3 О = О*, х € (х13 (¿), х32 (¿)). Здесь х13 (¿) — абсцисса границы между интервалами 1 и 3, а х32 (¿) — абсцисса границы между интервалами 3 и 2. В интервалах 1 и 2 решение имеет точно такой же вид, что и в первой задаче, с той лишь разницей, что в нем вместо постоянных А1 и В1 фигурируют постоянные А2 и В2. В интервале 3 оно записывается так

с3 = С*2 + Е2Ф

( х )

V2VЩ )

2у/Ш,

Постоянные А2, В2, С2, Е2 находятся из условий непрерывности функции с и диффузионного потока на границах х13 (¿) и х32 (¿). Чтобы получаемые уравнения не содержали время, необходимо, чтобы функции х13 (¿) и х32 Ю были бы пропорциональны л/1. Пусть

х13 = 2^1^07; х32 = 2^2л/в1г,

где Ц.1, 1^2 — постоянные величины. Приходим к системе шести уравнений относительно А2, В2, С2, Е2, ц1, ц2

с1 (х13) = с- ^ А2 = сф Л ^0 , с3 (х13) = с- ^ С2 + Е2Ф ^

ф(^1)

с,

с3 (х32) = с+ ^ С2 + Е2Ф ( ^I О = с+,

(/г-)

с2 (х32) = с+ ^ В2 =

О

а.

дс1

дх

дв3

дх

О,

Ж13-0

дс3

дх

ср<х> с+

[1 - Ф(р№)] , О

(8)

^ А2 = Е2\ о еХР

Ж13+0 V О1

Х32 —0

= О2 ^

дх

Ж32+0

^ В2 = Е2Рх/0ехр V О1

1 - ж) "2

(р 2 - а)

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

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

В*. Типичные зависимости с (х) приведены на рис.3. Они получены при следующих исходных данных: В1 = 1; В2 = 2; с- = 0.1; с+ = 0.2; ср0 = 0; Срж = 0.3, и значении Ь = 1. Из графиков следует, что точность аппроксимации растет с уменьшением В*. Докажем, что при В* 0 аппроксимирующее

решение совпадает с решением задачи 1. Исключим из системы (8) величины С2 и Е2. При этом два последних уравнения системы (8) примут вид

Ао =

сехр

Во =

с*ві ехр (в2 - ^2) ^

7:

Функции 1Л\ (7) и ц,2 (7) предполагаются аналитическими:

Ц-1 (7) = ^10 + ^117 + ^1272 + о (72) ; ^2 (7) = ^20 + №17 + ^2272 + о (72) •

Найдем предельные значения ^2 и В2 при 7 ^ 0. Конечные и отличные от нуля пределы получаются только при выполнении условий

^10 = ^20 = ^о; ^11 = ^21 •

Рис. 3. Зависимости с (ж) для момента времени Ь =1:1 — результат решения задачи 1; 2 — результат решения задачи 2 при В* = 0-5 (В\ + В2) = 1-5; 3 — результат решения задачи 2 при В* = 0Л

Раскрыв неопределенности, получаем

[Лсфо ехр (^°)

А2* = ІІШ А2 = - ,

7^0 1 - ехр [2^0 (^12 - ^22)]

В2 = ІІШ В2 =

пс*вМ0 ехр (в2^5)

7^0 ехр [-2^0 (М12 - М22)] - 1'

Исключив из этих равенств член exp [2^о (112 — I22)], получим уравнение л/Пв*вц,о = ^2*в exp (-^0) — В2* exp {—в2^2) •

С учетом того, что flo является предельным значением и для Ці, и для ц2, система уравнений относительно Л2*, B2*, ц0 совпадает с системой (7) относительно Ai, Bi, Л. Длина интервала 3 стремится к нулю, а соотношения для интервалов 1 и 2 оказываются такими же, как для задачи 1.

Таким образом, функция D (с) в интервале (с_,с+) должна быть равна нулю (график 2 на рис.2). Именно такая зависимость принималась при рассмотрении приведенного ниже примера.

Пусть ф (хк) — произвольная непрерывная кусочно-гладкая функция координат. Слабая формулировка задачи [6] получается из уравнения

У

дс

- -ъ- тс)

фйУ = 0, (9)

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

/ + У ЭпфдБ + У ]рф<1Б = 0. (10)

У Б„ Бр

Здесь ]р (Ь, х^) — проекция теплового потока на нормаль к поверхности Бр, то есть там, где задано главное граничное условие. В отличие от функции ,]п (Ь, х^) эта функция неизвестна.

Слабым решением уравнения (10) является функция с (Ь,хи), непрерывная в области V• В отличие от классического решения, ее производные по пространственным координатам и коэффициент диффузии V, рассматриваемый как функция пространственных координат, могут иметь разрывы первого рода. Это и позволяет получить решение задачи Стефана как слабое решение уравнения (10).

Выделим из пространства {ф} подпространство функций [ф], обращающихся в нуль при х^ € Бр. Обозначим через [в] дополнение к подпространству {ф}. Уравнение (10) при этом становится эквивалентно двум уравнениям . .

/ \%ф + тС ■ ду + /ЭпФаБ = 0, (11)

У Б„

/ (дв + ВУс -^в^дУ + I]пвйБ + I]рвйБ = 0. (12)

У Бп Бр

Так как зР входит только в уравнение (12), функция с (Ь, хк) определяется из одного уравнения (11). Затем, если это необходимо, из уравнения (12) при уже известной функции с (Ь,Хк) находится зР.

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

Для интегрирования по времени в настоящей работе применена конечноразностная схема Кранка-Николсон [5]. Пусть АЬ — заданная длина временных интервалов. Обозначим

Ьт = тАЬ, т = 0,1,2,..., ст (хк)= с (Ьт,хк),

вш = о (ст), 3т = 3 (¿т).

Для (т + 1)-го интервала можно записать приближенные равенства

д)с ст+1 _ ст 1 1

— = с-----------——, БУс = - (От+1Уст+1 + БтУст) , 1 = - (зт+1 + 1т ,

дЬ АЬ 2 v 2

где величины с индексом т — значения в начале интервала — известны.

Уравнения (11), (12) преобразуются к виду

IО+1 Уст+1 • уф + 2сАг>)ду = - / (зт+1+зт) ФЛБ -

V яп

У (' ОтУст • Уф - 2аЬф) йу, (13)

V

1 зт+1вёБ = - у зтвёБ -1 (зт+1+т ^б-

Зр Яр

2 (ст+1 - ст)

V

АЬ

-в + (От+1Уст+1 + БтУст) • У в

(IV, (14)

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

Уравнения (13), (14) решаются обобщенным методом Галеркина [4]. Пусть Ф(г) — базис подпространства {ф}. Функции ф(г) называются координатными функциями [4]. Решение задачи ищется в виде

ГО

ст+1 = £ а(г+1ф(г) (хк), (15)

г=1

где функция wm+1 (хк) задана и подчинена требованию хк € Бр : wm+1 (хк) = Ср (Ьт+1 ,Хк) •

При этом главные граничные условия (3) удовлетворяются. В отличие от классического метода Галеркина разложение (15), как и в методе Ритца, не обязано удовлетворять естественным граничным условиям.

Так как любой элемент подпространства {ф} представляет собой линейную комбинацию базисных элементов, достаточно, чтобы уравнение (13) было справедливо для всех координатных функций ф(1), I = 1, 2,.... После подстановки разложения (15) в уравнение (13) и последующего интегрирования для всех функций ф(1) получается бесконечная система уравнений относительно коэффициентов разложения ат)+1. В результате ее решения находится функция Ст+1 (хк). Отметим, что эта система уравнений нелинейна. Для ее решения можно, как это сделано в настоящей работе, ввести итерационный процесс для определения От+1, положив в первом приближении От+1 = От. При достаточно малом шаге Д£ можно обойтись без последующих итераций.

После того как определена функция ст+1 (хк), находится функция ]рт+1 (хк). Пусть 9- базис подпространства {0}. Функция ]рт+1 (хк) ищется в виде

зрт+1 = ЁЬт+%) (хк), хк € Бр.

г=1

Требование, чтобы уравнение (14) выполнялось для каждой координатной функции 9(1), I = 1, 2,..., приводит к системе линейных алгебраических уравнений относительно коэффициентов разложения Ът+1.

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

В качестве примера приложения метода рассмотрим задачу, имеющую практическое значение. Это задача о диффузии в системе никель-титан при температуре 900°С. Фазы и диапазоны изменения концентрации никеля (табл.) взяты из диаграммы состояния [12]. Значение коэффициента диффузии принималось для всех фаз одинаковым и равным О = 9.43 • 10_14 м2/с [13]. Длина пространственного интервала I = 120 мкм. На границах интервала диффузионный поток равен нулю. Начальное условие имеет вид

С(х) = I1. х € [°, 6°)

'( ) 1.0, х € [60, 120] '

Диапазоны изменения концентрации никеля

№ фазы Фаза С

1 N1 0.892 - 1

2 Т1№3 0.708 - 0.785

3 Т1№ 0.485 - 0.523

4 Т12 N1 0.292 - 0.400

5 Т1 0 - 0.085

Относительная погрешность расчета < 10_3 достигнута при разбиении пространственного интервала на 2000 конечных элементов и шаге по времени ДЬ = 10_2 с. Некоторые результаты расчетов представлены на рис.4,5.

Дт

Т

0.4 0.3 0.2 0.1 О

2 4 6 час

Рис. 4. Зависимости относительной ширины фаз от времени;

1-5 — номера фаз (таблица)

0.8

0.6

0.4

0.2

°0 0.2 0.4 0.6 0.8 х/1

Рис. 5. Изменение концентрации никеля по длине пространственного интервала для момента времени I = 8 часов

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

I, 5. График зависимости концентрации никеля от координаты x (рис.5) для момента времени t = 8 часов, то есть через 2880000 шагов по времени, свидетельствует об эффективности предложенного метода расчета: фазовые границы не расплываются и диапазоны изменения фаз соответствуют заданным интервалам.

Список литературы

1. Рубинштейн Л.И. Проблема Стефана. Рига: Звайгзне, 1967.

2. Мейрманов А.М. Задача Стефана. Новосибирск: Наука, 1986.

3. Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. М.: Едито-риал УРСС, 2003.

4. Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970.

5. Рихтмайер Р.Д., Мортон К. Разностные методы решения краевых задач. М.: Мир, 1972.

6. Зенкевич О., Морган К. Конечные элементы и аппроксимация. М.: Мир, 1986.

7. Любов Б.Я. Диффузионные процессы в неоднородных твердых средах. М.: Наука, 1981.

8. Курант Р. Уравнения с частными производными. М.: Мир, 1964.

9. Карлсроу Г., Егер Д. Теплопроводность твердых тел. М.: Наука, 1964.

10. Курант Р., Гильберт Д. Методы математической физики. Т.1. М.: Гостехиздат, 1951.

II. Диаграммы состояния двойных металлических систем // под ред. Н.П. Ляки-шева. М.: Машиностроение, 2001.

12. Bastin G.F., Rieck G.D. Diffusion in the titanium-nickel system: II. Calculations of chemical and intrinsic diffusion // Metallurgical Transactions, B.1974. V.5. №8. P.1827-1831.

Жигунов Виктор Владимирович (vzhigunov@rambler.ru), д.т.н., профессор, кафедра физики, Тульский государственный университет.

Лавит Анна Игоревна (annlavit@rambler.ru), аспирант, кафедра физики, Тульский государственный университет.

Mathematical modeling of diffusion with due regard for appearance and disappearance of phases

V. V. Zhigunov, A. I. Lavit

Abstract. A variant of shock-capturing method for Stefan problem solution that emerges with modeling of diffusion in multiphase solid body is proposed.

Specific character of the problem is that concentration of diffused element has discontinuity on interfaces. With introduced method this discontinuity is simulated without interface coordinates computing in solution process. It is achieved because of curve shape that is diffusion coefficient function where concentration is argument. Generalized Galerkin method is used because its coordinate functions not obliged to satisfy natural boundary conditions. Time integration is executed with Crank-Nicholson finite-difference scheme. Space integration is executed with finite element method. Example of calculations is confirmed efficiency of the method.

Keywords : diffusion, phases, Stefan problem, Crank-Nicholson scheme, Galerkin method, finite element method.

Zhigunov Victor (vzhigunov@rambler.ru), doctor of technical sciences, professor, department of physics, Tula State University.

Lavit Anna (annlavit@rambler.ru), postgraduate student, department of physics, Tula State University.

Поступила 21.01.2013