МОДЕЛИРОВАНИЕ ОДНОМЕРНОЙ ДИФФУЗИИ ВОДОРОДА В МЕТАЛЛАХ. V. ПРОНИЦАЕМОСТЬ ПРИ НАЛИЧИИ ОБРАТИМОЙ ХИМИЧЕСКОЙ РЕАКЦИИ ВТОРОГО ПОРЯДКА Текст научной статьи по специальности «Физика»

Статья поступила в редакцию 17.04.11. Ред. рег. № 977

The article has entered in publishing office 17.04.11. Ed. reg. No. 977

УДК 539.219.3:669.788

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

1 2 В.Н. Лобко , И.Н. Бекман

владимирский государственный университет 600000 Владимир, ул. Горького, д. 87 Тел. (4922)47-98-67, e-mail: lobko_vn@laser-2.vpti.vladimir.ru 2МГУ им. М.В. Ломоносова 119991 Москва, ГСП-1, Ленинские горы, д. 1, стр. 3 Тел. (495)939-32-12

Заключение совета рецензентов: 27.04.11 Заключение совета экспертов: 28.04.11 Принято к публикации: 30.04.11

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

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

NUMERICAL MODELING OF ONE-DIMENSIONAL DIFFUSION OF HYDROGEN IN METALS. V. PERMEABILITY IN THE PRESENCE OF REVERSIBLE CHEMICAL REACTION

OF THE SECOND ORDER

V.N. Lobko1, I.N. Beckman2

'Vladimir State University 87 Gorky st., Vladimir, 600000, Russia Tel. (4922)47-98-67, e-mail: lobko_vn@laser-2.vpti.vladimir.ru 2Moscow State University 1 Leninskiye Gori, Moscow, GSP-1, 119991, Russia Tel. (495)939-32-12

Referred: 27.04.11 Expertise: 28.04.11 Accepted: 30.04.11

In the first article the mathematical foundations of the method for cases of constant diffusion coefficient and its dependence on the arbitrary concentration of hydrogen were presented. In this work a similar problem was solved in the presence of reversible chemical reaction of second order between the diffusing gas and the material of the plate in the diffusion volume with a constant diffusion coefficient. The influence of the equilibrium constant of chemical reaction and the ultimate degree of conversion were considered.

Keywords: diffusion, hydrogen, metals, modeling, coefficient of diffusion, permeability, finite difference method, concentration dependent of diffusion coefficient, profile of concentration of hydrogen.

Введение

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

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

и коэффициента диффузии, зависящего от концентрации водорода. Приведено сравнение расчетов по двум вариантам алгоритмов.

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

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

Теория метода

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

дс dt

D (1 + K рс )2

д2с

KvC2m +(1 + Kpc)

(1)

где Кр - константа равновесия химической реакции; С2т - предельная степень превращения (в дефектной среде - максимальная емкость ловушек частиц диффундирующего газа).

Введем равномерную сетку по оси х: / = 0, 1, ..., п и по оси /:] = 0, 1, ..., т. Обозначим И - шаг по оси х, т - шаг по оси /. Применим метод конечных разностей в варианте, приводящем к неявной схеме:

с. . — с. . ,

D (1 + KP c,,j )2 с,—u — 2си + с,

i+1, j

K C m +(1 + Kp с,, j )

h

(2)

Полученная схема представляет собой систему нелинейных уравнений, которую можно решить итеративным методом Ньютона. Перепишем уравнение (2) в более удобном виде и введем

F — с,, j —! ) —

Dt

(1 + K Р с, j )2

-(с—1,j — 2с,j + с+1,j ) = 0 •

h KpC2m +(1 + Kpс,,j )

Вычислим соответствующие производные для реализации метода Ньютона:

dFt dFt =— Dt (1 + Kp с, j )2 ; дс,—1, j дс, + 1, j h 2 KpC2m +(1 + Kp с,, j )2'

(3)

(4)

dF 2Dt д

(

dc.j h dcu

Kp2 с3 j + 2Kp < j + с,, j

Л

+1 =

= 1 +

vKp2 < j + 2Kp^. j + (1 + KpC2 m),

2Dt (3KXj + 4Kp с. j + 1)(Kp2 j + 2Kp с, j + (1 + K£2т)) — 2K p (Kf . j +1)( K\ < j + 2Kp с,2 j + с, y) h2

(( < j + 2K p^. j + (1 + K pC2m ))2

Уравнением (3) можно аппроксимировать внутренние узлы сетки; на входной ( = 0) и выходной ( = п) сторонах пластины используем выражения, полученные ранее [1] по 1-му и 2-му вариантам моделирования:

Л =—_[-г-] и Лп = _ (первый вариант);

RTSKr\dt) 0 п

RTSKr [ dt

т 2F' (дс Л г 2^2 (Эс Л ( й )

J = —____'_„ с I — I и Jn =____:_„ сп | — I (второй вариант)

RTSKC 01 dt

RTSKC п [dt

в комбинации с первым законом Фика. Обозначим константы: для одноатомного газа

International Scientific Journal for Alternative Energy and Ecology № 4 (96) 2011

© Scientific Technical Centre «TATA», 2011

T

A =■

V_h J_

RTSKT t D

(5)

и для двухатомного газа

A =■

V h 1

RTSK2 t D

(6)

где V- это или входной У1, или выходной У2 объем).

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

Первый вариант моделирования; одноатомный газ:

(А +1)со,, - с,,, - Лс,,,-1 = 0; г = 0;

От (1 + Кр с,,,)

((- )-~ (- 2е-+)=0;'=12 n -1; -i,j -(A+j+Aicn,j-i =0;' =n

(Ai +1) X0 - Xi =-F(k-1); i = 0; dF (k-1) dF(k-1) dF(k k-1)

dFi .x„ , -X,. +dF-X..+, = -Fi(k-1); i = 1, 2, ..., n-1;

de.

dc

de,.,

"i'-1,j i,j ' + 1, i Xn-1 -((1 + 1)Xn =-Fn(k-1); i = n

Первый вариант моделирования; двухатомный газ

A2eo2,j- (A2eo,j-1-1) eo,j- e1,j =0;' =0;

Dt (1 + Kp C, j )2

(.j - C,j-1 )-

h KpC2m +(1 + KP C, j)

en-1,j - A2e„2,j +((2en,j-1 - 1)en,j = 0; ' =

( - 2e,,j + e,+u) = 0; i = 1, 2, ..., n -1;

(7)

(8)

(9)

(A2 (2e0k^1) - e0,j-1) +1) X0 - X1 = -F0k-1); i = 0;

dF.(k-1) dF<k-1) 3F.(k-1) (k 1)

■ Xt-1 —r-X,' + -X,' +, = -F(k-1); i = 1, 2, ..., n -1;

' i+1

dC+1, j

X„ = - F„( k-1); i = n

de<-ij " de¡,J

Xn-1 -((2 (j - en,j-1)+1).....

Второй вариант моделирования; одноатомный газ:

(2A +1) e0,j - el,j - ((2A1 -1) e0,j-1 + e,^) = 0; i = 0;

j- e, j-1)-

Dt

(1+Kp e, j )2

(-1,j - 2c,,j + e+1,j) = 0; i = 1, 2, ..., n -1;

h KpC2m +(1 + Kp c,, j)"

en-1,j - (2A1 + 1)en,j + ((2A - 1)en,j-1 + Cn-1,j-1) = 0; i = n (2A +1) X0 - X1 =-F0k-1); i = 0;

dF (k-1) 3F.(k-1) 3F.(k-1) (k n ^-Xi-1 -Xt + —.-Xi+1 = -F.(-1); i = 1, 2, ..., n-1;

X„-1 -(2A1 +1) Xn = -F(k-1); i = n

(10)

(11)

(12)

Второй вариант моделирования; двухатомный газ:

Л-2,1 +—о, 1- -1,1- (Л-2,1 -1- —о, 1 -1+-1,1 -1) =0;' =0; , , пт (1 + Ко..) , .

( - о.,1-1)-И--У / ' \2 (о.1 - 2С,1 + о,+1,1 ) = 0; г = 1, 2, ..., п -1; (13)

И КрС2т +(1 + Кр о, ,1)

оп-1,1- Л2<1- оп, 1+((2<1 -1- оп, 1 -1 + оп-1,1 -1) =0;' = п

^Л^ +1) Х0 - X! =--1); , = 0;

ЭР(к-1) э^(к-1) э^(к-1)

^-X. 1 -X. -X.+1 = -^(к-1); , = 1, 2, ..., п -1; (14)

"V ,-1 "V , "V , + 1 , 7 1 1 1 7 V /

д0,-1,1 до,,1 до,+1,1

Xn-l -(2Л2-п^ + 1) =-к-1); , = п

Системы метода Ньютона (8), (10), (12) и (14) дополняются выражениями для производных (4). Они решаются обычным методом прогонки относительно поправок:

X, = - о?/, (15)

где о(к - искомые приближения концентраций последующей итерации.

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

Ниже приведены примеры моделирования диффузии водорода в металлах при наличии обратимой реакции второго порядка по первому и второму вариантам моделирования в сравнении со случаем D = const = 5,66-10-10 м2/с без реакции. В программе рассчитываются концентрационные профили диффу-занта, временные зависимости концентраций c(t) и потоков J(t) на входной и выходной поверхностях пластины, а также - давлений p(t) в резервуаре и приемнике. Поскольку наиболее показательным параметром являются выходные потоки J(t) метода проницаемости, приведена именно эта зависимость. Во всех последующих примерах имели место параметры: толщина пластины, 1,5-10-3 м; температура, 680,7 К; площадь пластины, 3,14-10-4 м2; коэффициент диффузии, 5,66-10-10 м2/с; константа Сивертса,

о ^ п s о ^

2,45-10- моль/м/Па , ; объем резервуара, 5-10- м; объем приемника, 5-10-5 м3; давление в резервуаре, 8,95 -104 Па; начальное давление в приемнике, = 0 Па; начальная концентрация в пластине, 1-10-52 моль/м3; точность метода Ньютона, 10-6; число разбиений по оси х, 100; конечное время 10000 с; число разбиений по оси t, 500.

На рис. 1 представлены результаты моделирования для C2m = 5 моль/м3 при константах равновесия Кр = 0,01 м3/моль; Кр = 0,5 м3/моль; Кр = 1 м3/моль. Наличие реакции второго порядка приводит (как и в последующих примерах) к замедлению выходного потока («вертикальная» часть кривой), причем затем все зависимости потока от времени, J(t), достигают

одинакового значения. Заметный вклад реакционной диффузии наблюдается при Кр порядка нескольких десятых м3/моль. Дальнейшее увеличение Кр дает меньший эффект замедления.

Рис. 1. Моделирование диффузии с реакцией второго порядка. Выходной поток J(t) при C2m = 5 моль/м3; влияние Кр.

1 (тонкая линия) - без реакции 2-го порядка; 2 - Кр = 0,01 м3/моль; 3 - Кр = 0,5 м3/моль; 4 - Кр = 1 м3/моль; 1-й и 2-й варианты моделирования (расхождений нет) Fig. 1. Numerical modeling of diffusion with reaction of the second order. Output flux J(t) at C2m = 5 mol/m3; influence of Кр. 1 (thin line) - without reaction of 2nd order; 2 - Кр = 0,01 m3/mol; 3 - Кр = 0,5 m3/mol; 4 - Кр = 1 m3/mol; 1st and 2nd versions of numerical modeling (there aren't divergences)

На рис. 2 и 3 показано влияние емкости среды -параметра С2т. При выбранных значениях С2т происходит сильное замедление выходного потока газа, причем этот эффект больше при большем значении константы равновесия реакции Кр.

International Scientific Journal for Alternative Energy and Ecology № 4 (96) 2011

© Scientific Technical Centre «TATA», 2011

27

2 6 f,c, х10э 10

Рис. 2. Моделирование диффузии с реакцией второго порядка. Выходной поток J(t) при Кр = 0,5 м3/моль; влияние C2m. 1 (тонкая линия) - без реакции 2-го порядка; 2 - C2m = 5 моль/м3; 3 - C2m = 15 моль/м3; 1-й и 2-й варианты моделирования (расхождений нет) Fig. 2. Numerical modeling of diffusion with reaction of the second order. Output flux J(t) at Кр = 0,5 m3/mol; influence of C2m. 1 (thin line) - without reaction of 2nd order; 2 - C2m = 5 mol/m3; 3 - C2m = 15 mol/m3; 1st and 2nd versions of numerical modeling (there aren't divergences)

Заключение

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

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

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

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

27

Рис. 3. Моделирование диффузии с реакцией второго порядка. Выходной поток J(t) при Кр = 1 м3/моль; влияние C2m. 1 (тонкая линия) - без реакции 2-го порядка; 2 - C2m = 5 моль/м3; 3 - C2m = 15 моль/м3; 1-й и 2-й варианты моделирования (расхождений нет) Fig. 2. Numerical modeling of diffusion with reaction of the second order. Output flux J(t) at Кр = 1 m3/mol; influence C2m. 1 (thin line) - without reaction of 2nd order; 2 - C2m = 5 mol/m3; 3 - C2m = 15 mol/m3; 1st and 2nd versions of numerical modeling (there aren't divergences)

1. Лобко В.Н., Бекман И.Н. Моделирование одномерной диффузии водорода в металлах. I. Метод расчета изотермической диффузии при постоянном и переменном коэффициенте диффузии // Альтернативная энергетика и экология - ISJAEE. 2010. № 10. С. 36-42.

2. Лобко В.Н., Бекман И.Н. Моделирование одномерной диффузии водорода в металлах. II. Сорбция и десорбция в плоскопараллельной пластине, находящейся в замкнутом объеме // Там же. 2010. № 11. С. 38-42.

3. Лобко В.Н., Бекман И.Н. Моделирование одномерной диффузии водорода в металлах. III. Проницаемость через плоскопараллельную пластину, контактирующую с замкнутыми объемами // Там же. 2011. № 4. С. 15-19.

4. Лобко В.Н., Бекман И.Н. Моделирование одномерной диффузии водорода в металлах. IV. Проницаемость при коэффициенте диффузии, зависящем от концентрации диффузанта // Там же. № 4. С. 20-24.

5. Бельнов В.К., Бекман И.Н., Сафонов М.С. Режимы диффузионной кинетики газа в микрогетерогенных материалах // ЖФХ. 1995. Т. 69, № 9. С. 1566-1571.

6. McNabb A., Foster P.K. A new analysis of the diffusion of hydrogen in iron and ferritic steel // Trans. Metall. Soc. AIME. 1963. Vol. 227, No. 3. P. 618.

7. Norgett M.J., Lidiard A.B. Radiation Damage in Reactor Materials, Diffusion of Inert Gases in Ionic Crystals. IAEA: Vienna. 1969. 61.