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

УДК 517.958: 544.034.54

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

© В. С. Зверев

Уральский государственный университет им. А. М. Горького Россия, 620083 г. Екатеринбург, пр. Ленина, 51.

Тел: (343) 350 74 29.

E-mail: v-s-zverev@yandex. ru

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

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

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

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

Рис 1. Распределение диффузионных потоков в радиальной геометрии.

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

В обоих случаях u(t,r,z) - концентрацию диффузанта на поверхности через, а внутри - w(t,r,z). Приведем сначала модель для радиальной геометрии и на ее примере опишем основные этапы исследования. Диффузионный поток по поверхности подложки в системе координат рис. 1 будет описываться уравнением (1), в котором учитывается цилиндрическая форма подложки и отток диффузанта с поверхности.

du ^ 1 Э ( du ^ D2 Эw

— = D-------1 r— I+——

Эt r Эг ^ Эг J S Эz

u(t,rd) = u0, u(0, z) = 0.

t > 0, r > r,, z > 0

(1)

Процесс диффузии внутри подложки может быть описан следующим образом:

Эw (t, r, z) Э2 w(t, r, z)

D2

Эt 2 Эz2

w(0, r, z) = 0, w(t, r,0) = u(t,r)

t > 0, r > rd , 0 < z <S(t,r)

D — 2 Эz

z =S

= —hw|

(2)

lz = S

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

ЭЕ(^ r) hw(t, r, z)

Эt

j

S(0, r) = 0.

z = S(t, r)

(3)

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

границы, то есть можно считать X малой величиной. Тогда можно считать характерное время диффузионного процесса внутри подложки тв << тх времени движения границы реакции, а значит в подзадаче (2) можно не искать точное решение уравнения параболического типа с подвижной границей, которое является весьма громоздким, и полагать, что w»0. Это является основанием использования метода дифференциальных рядов, описанного например в [5], и того, что в нем можно ограничиться первым приближением. Строгое обоснование этих шагов для разных геометрии расположения реагентов выходят за рамки данной статьи.

Для дальнейшего изучения введем безразмерные переменные и функции:

к

Л =---г.

А

_ и _ w и =

и и,

_ , к , =

, = —, 5 (т,р) = —£(<, г), где <0 = ,2

и А к ип

Тогда из того, что и уравнения (2) полу-Л

1 --

а значит, подставляя это выра-

жение в (1) и (3), придем к системе

1 Э ( ди | р и р Эр [Р Эр) 1 + 5

(4)

= 5 (1+5)

Перепишем первое уравнение системы (4), проинтегрировав по г и учтя начальные условия:

т

и = Др| Шт- Р5, где, как и в предыдущем выраже-

0

нии Р = Б2ф1 кди0 - безразмерный параметр модели.

Тогда, использовав факт 5 << 1 и второе соотношение в (4), сведем исходную задачу (1-3) к следующей нелинейной краевой задаче на луче:

1А (рА Г 5+52' и=р5, рЭрГ Эр1 2 11 ’

(5)

5(р ® ~,т) = 0, 5(р = рл ,т) = £„(т) = -1 + V1 + 2т, где рл = , а 80(т) является решением диффе-

ренциального ур авнения и(т,ра) = 5 (т,рл) (1 + 5 (т,рл))

при условиях п(т,рл) = 1 и 5 (0,р) = 0.

Задача (5) ввиду её нелинейности решалась численно. Для применения численных методов луч [ра,¥) заменой р = ра/р был переведен в отрезок [0, 1]:

рЗ р4

і^8'+р^8" = Р(71+2^-1), я = 5 + 572'

(6)

ра ра

Функция / (V) = р (V1 + 2у -1) непрерывна и ограничена вместе со своими производными на любом конечном отрезке изменения параметров, допус-

каемых физическим смыслом задачи, то есть неотрицательности функции, что обеспечивает сходимость конечно-разностного метода для краевых задач. Введем на отрезке [0,1] равномерную сетку с шагом к количеством точек N = 1/к и заменим производные разностными выражениями

8(р,) » 8м~,8‘-1, 8(р,) » ^ - 28‘ + 8‘-1. (7)

2к

Подставляя их в дифференциальное уравнение (6), получим систему алгебраических уравнений:

«,й,+1 + Рй, + Ші-1 = ї(й) = р (уі 1 + 2й, -1)

(8)

Уравнения (8) решается методом Ньютона:

й™ = +дм

(.+1) + ^(.+1) = г(£«) + /й(й«)дм 1 <, < N -1 (9)

«й (+1+1) + д'г) = Д1? = 0

Линеаризованная система в свою очередь решается методом прогонки.

В качестве критерия окончания итерации для решения (8) выбрано условие||8® -8,*+ц|<е, что

согласно [6] является хорошей оценкой для метода Ньютона решения нелинейных систем, ввиду его квадратичной скорости сходимости.

Длина поверхностного слоя 1(т) нового вещества определяется из алгебраического соотношения

5 (т, 1з (т)) = 0. Корень этого уравнения находился методом деления отрезка пополам, а также с учетом монотонности функции 5, которая следует из неотрицательности величин и, 5 и второго уравнения в системе (4). Алгоритм был реализован в системе МаЛСаё. Результаты численных расчетов представлены на рис. 2.

Перейдем теперь к рассмотрению продольной геометрии расположения реагентов. Отличительные особенности этого случая хорошо видны на рис. 3.

Рис. 2. Длина поверхностного слоя при различных параметрах в радиальной геометрии.

г

ит =

и

Рис. 3. Распределение диффузионных потоков в продольной геометрии.

Приведем сразу модель, описывающую этот случай:

Эи = ^ 32и А2 3щ дt 1 З22 8 Эг г=

u(t,0) = и0, и(0,2) = 0.

t > 0, 0 < г < г , 2 > 0,

(70)

Эм (т, г, 2) 1 Э Г Эм (Т, г, 2) 1

---------------- = А,— — | г—Ц------------------------- I, Т > 0, Х(Т, 2) < г < г,, 2 > 0,

Эt

ЭХ(Т, 2) = кщ(Т, г, 2)

Эг

Эм

2 Эг

(77)

= кщ

г = 2

г = 2

Эт

ф

г = Х(Т, 2)

, 0<Х< г,,Х(0,2) = г,. (72)

Уравнение (70) вместе с начальными условиями описывает диффузионный поток по поверхности подложки. В описании процесса диффузии в подложки учитываем цилиндрическую форму, поэтому в уравнении (77) в отличие от предыдущего случая оператор Лапласа в цилиндрических координатах. Соотношение (72) определяет закон движения границы.

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

Т к 2

т = —, р=— г,ц =

А

А1Т0

- и _ щ „ к „ а ф

и = —, щ = —, 5 = — X, где Т0 = А2ф

А

к и„

К модели для продольной геометрии применимы такие же рассуждения о характерных временах

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

(73)

1 - 51п (5/р,) | р

Подставив (73) в (70) и (72), придем к системе:

' Р иБ______

/ р, )

(74)

5 = —

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

* Г Г 1п р -11 - 51 =-р з 2-рр),

Эч\ 2 | рр 2) J 2р, ’ (75)

5 (т, л ® ~) = р,, 5 (т,ц = 0) = 50 (т).

Аналогично предыдущему случаю 50(т) , которое определяется из решения дифференциального уравнения и(т,0) = (5(т,0) 1п(5/рр)-1)5(т,0) при условиях и (т,0) = 1 и 5 (0,^) = р,, а именно уравнения (76).

(76)

Уравнение не зависит от свободной переменной, поэтому можно понизить его порядок: сделаем замену и перейдем к системе уравнений первого порядка.

'35

8

Эл 51п (5/р,)-1

*=-2р:(5 2-р; >

(77)

Запишем

её

38 _ Р (рр -52)(51п(5/р,)-1)

симметрической форме Это уже уравнение с

35 2р, §

разделяющимися переменными. Тогда учитывая здесь условие 5(т,Л ® ¥) = р, и решив квадратное уравнение относительно д, можно найти:

0

в

и

и

0

54 Г 5 1 1 2 52 Г 5 1 1 53 ^ з Г 3 2

Р\ 4 [“р, 4] + рр 2 [“р, 21 +3 р + А[ 16рр + 3

51п — -1

Учитывая, что функция 5 в данной геометрии возрастает и на интервале 0 < 5 < р, 51п (5 /р,)-1 < 0, необходимо выбрать знак минус. Теперь можно записать формальное решение исходной задачи. Знак минус учтем, поменяв пределы интегрирования.

На рис. 4 представлен график 1,(т) для продольной геометрии, а на рис. 5 объединены случаи разных геометрии.

Длина поверхностного слоя определяется из условия 5(т,I, (т)) = р, или, принимая во внимание рассходимость интеграла в точке 5=р,, можно выбрать 5(т,^(т)) = р, -£. Тогда:

I

V 1п

-1 | + 12р>21 21п| — | -1 | + -48р> + 48С

.-3v4| 41п|

■''/4^ V I I р,

=— л

(79)

1, (т) = 4Цт

5с(т)

V 1п

-1

р,-е -3v | 41п

-1| + 12р,2v2

21п

-1| + 16v3 - 48р> + 9р4 + 32р,3

(20)

V

1

5

V

V

V

Рис.4. Длина поверхностного слоя при различных параметрах в продольной геометрии.

р = 1

--------- р, = 10

......... ра = 4

--------- ра = 1

________ра = 0.2

Рис.5. Сравнение длин поверхностного слоя для продольной (сплошная линия) и радиальной геометрии.

Заключение

Функции заданные неявным образом, а также величины связанные с ними исследовались с помощью математического пакета МаЛСАБ. И в радиальной, и продольной геометрии в результате моделирования поверхностной реакционной диффузии была получена система, состоящая из двух уравнений в частных производных параболического типа. Одно из уравнений имеет подвижную границу, движение которой описывалось обыкновенным дифференциальным уравнением.

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

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

Работа выполнена при финансовой поддержке РФФИ, грант №07-01-96091-р_урал_а.

ЛИТЕРАТУРА.

1. Нейман А. Я. // Ионы в твердофазном состоянии. 1996. Т. 83.

С 263-273.

2. Нейман А. Я., Утюмов В. Ю., Карпов С. Г., Костиков Ю. П., Шиятова М. В. // Поверхность. 2000. № 3. С. 52-61.

3. Нейман А. Я. // Журнал физической химии. 2001. Т. 75. №12. С. 2119-2134.

4. Нейман А. Я., Шиятова М. В., Карпова С. Г., Костиков Ю. П. // Поверхность. 1996. №11. С. 20-28.

5. Карташов Э. М. Аналитические методы в теории теплопроводности твердых тел: Учеб. Пособие. 3-е изд., перераб. и доп. М.: Высш. шк., 2001. -550с.

6. Калиткин Н. Н. Численные методы. М. Наука, 1978. -512 с.

Поступила в редакцию 06.09.2008 г.