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

Статья поступила в редакцию 27.09.10. Ред. рег. № 887 The article has entered in publishing office 27.09.10. Ed. reg. No. 887

УДК 539.219.3:669.788

МОДЕЛИРОВАНИЕ ОДНОМЕРНОЙ ДИФФУЗИИ ВОДОРОДА

В МЕТАЛЛАХ

I. МЕТОД РАСЧЕТА ИЗОТЕРМИЧЕСКОЙ ДИФФУЗИИ ПРИ ПОСТОЯННОМ И ПЕРЕМЕННОМ КОЭФФИЦИЕНТЕ ДИФФУЗИИ

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

Заключение совета рецензентов: 17.10.10 Заключение совета экспертов: 27.10.10 Принято к публикации: 31.10.10

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

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

NUMERICAL MODELING OF ONE-DIMENSIONAL DIFFUSION OF HYDROGEN IN METALS

I. METHOD OF CALCULATION OF ISOTHERMAL DIFFUSION AT CONSTANT AND VARIABLE DIFFUSION COEFFICIENT

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: 17.10.10 Expertise: 27.10.10 Accepted: 31.10.10

A series of works devoted to modeling one-dimensional diffusion of hydrogen and atomic gas in the plane-parallel plate, bordering on arbitrary closed volumes with the boundary conditions of the first kind, which, however, clearly not defined. Two variants of the calculation formulas allow making a comparison and self-inspection. The method allows constructing the concentration profiles of diffusant on thickness of the plate, defining flows and concentrations of gas depending on time on both sides of the plate and the pressure in both volumes.

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

International Scientific Journal for Alternative Energy and Ecology № 10 (90) 2010

© Scientific Technical Centre «TATA», 2010

Введение

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

Настоящая серия работ посвящена отработке численных методов решения дифференциальных уравнений массопереноса с неявно заданными граничными условиями, когда невозможно использование известных аналитических и численных решений краевых задач. На практике это означает моделирование диффузии в плоскопараллельной пластине, граничащей с произвольными замкнутыми объемами или с одним общим объемом. В первой статье представлены математические основы метода для случаев постоянного коэффициента диффузии D = const и произвольной зависимости коэффициента диффузии от концентрации водорода D = f(c). Рассмотрена диффузия как одноатомного, так и двухатомного газа, диссоциирующего на входной поверхности мембраны (например водорода). Предложены два фактически разных варианта расчетных формул, что позволяет провести сравнение и самопроверку. Варьируя величины объемов и начального давления газа в них, можно одной компьютерной программой охватить случаи сорбции-десорбции и проницаемости, включая случаи ограниченных объемов приемника. В остальных работах серии приведены результаты расчетов сорбционно-десорбционного метода, метода проницаемости для D = const и D = fC) и некоторые другие практически важные ситуации. Представлен большой иллюстративный материал.

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

Рассматривается первая краевая задача для двух случаев. В первом из них пластина находится в замкнутом постоянном объеме и решается задача сорбции-десорбции. При этом пренебрегаем диффузией с торцов пластины. Во втором случае пластина контактирует с двумя постоянными замкнутыми объемами, разделяя их; рассматривается задача сорбции-десорбции и проникновения (проницаемости). Предполагается выполнение граничных условий первого рода, т.е. приповерхностная концентрация диффузан-та в пластине связана с давлением газа законом Генри (одноатомный газ) или Сивертса (двухатомный газ); соответствующие константы, КГ и КС, известны. Известны также коэффициент диффузии Б, объем резервуара V и приемника У2 (или один общий) с соответствующими начальными давлениями газа р\ и р2, толщина I и площадь пластины, начальная концентрация газа в пластине - постоянная с0 или заданная в виде функции с(х). Требуется решить диффузионную задачу, то есть построить концентрационные профили диффузанта по толщине пластины, определить в зависимости от времени потоки и концентрации исследуемого газа на обеих сторонах пластины и давления в соответствующих объемах.

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

Из определения удельного потока 3 = сЫМБ, где 5 - площадь пластины, п - число молей, t - время, и уравнения состояния идеального газар¥1Т = пЯ, где Т - температура, Я - универсальная газовая постоянная, можно получить формулу, связывающую поток газа через поверхность пластины (со стороны резервуара или приемника) с давлением в газовой фазе:

J =

V dp RTS dt '

(1)

где J - плотность потока (удельный поток), dpldt -производная давления по времени.

Эта формула справедлива для граничных условий первого, второго и третьего рода. В последних двух случаях выражение для потока входит в уравнение материального баланса. Если начало координат находится на входной стороне пластины, поток берется со знаком «минус». Введем равномерную сетку по оси х: i = 0, 1, ..., n (тогда поток на входной стороне будет обозначаться J0, а на выходной - Jn). По оси t также введем равномерную сетку j = 0, 1, ..., m. Обозначим h - шаг по оси х, т - шаг по оси t. Уравнение (1) может быть скомбинировано с первым законом Фика

J = -D £ (2)

ох

и законом Генри (А1 +1) с0,, - с,,, - А1с0,,= 0

с = кгР (3) и сп-1,, -(А, +1)ся,, + Дси, , -1 = 0; (13)

или Сивертса А2с2,, - (А, 1)с0,, - С,, = 0

с = Кс4р (4) и Сп-1,, - А <, + (А2 Сп,,-1 - 1)Сп,, = 0. (14)

с учетом соотношении Вариант 2 (интегрирование). Альтернативные

уравнения можно получить интегрированием (1) с

Лгаза Л0 и Лп Лгаза- (5)

Расчетные формулы могут быть получены в двух

разделением переменных:

альтернативных вартантах. _ fJ ш t и YMt t (15)

Вариант 1 (дифференцирование). Комбинация J RTS J rJ RTS J

(1) и (3) дает (6), а (1) и (4) - (7): °J-1 " "

fj _1 -nj_1

Левый интеграл заменяем по формуле трапеций:

-0^-1 + (0,,) = -Т;(( -р,-1)

(п=- ~гкс С0 Цт J 0 и (п=^ТкСп Цт J п • (7) и | (-1+(,,)=(- р,-1) (16)

Используя (2), получаем, соответственно, (8) и (9): и используем (2) - (4), заменив в (2) истинную производную разностной для одноатомного газа:

Jo f-1 и J _-VL_ f^ 1 ; (6)

0 RTSKr\dt n RTSKr f dt Jn , \ V

2V1 (de 1 2V2 (de

V (Эе1 _ (de

и ^f»e )_-D(de, ; (8)

Y _ C°J + e1,j_1 _ eo,j_1 )_ ^fK ( _ eo,J_1 )(17)

RTSKr{3t ,дх D z( + V

RTSKr , dt In , дх

2 hv 1,j °j 1,j°j_4 RTSKr

г v"' \^Jn d т/ \ V2 ,

---(e _e , + e . , _e , . ,|_-(e . _e . ,

2V, e f de) _ D fdel 2 ^ n,j n_1,j j n_1,j_^ RTSKr[ n,j j

RTSKf0 ,dt J° ,дх J° (18)

2 V2 fdc) nfdc) и для двухатомного газа

и RTSKC e U Jn D UJn • (9) D V

Записываем разностные аналоги уравнений (8) и (9):

-(( -c°,j + e1,j_1 _e°,j_1 ) _ -RTK ((j _e°2,j_1 ) ,(19)

2 hv 1j j 1j_ j-4 RTSKC

_ D—( - _ )_ V2 ( 2 _ 2

V e°,j _ e°,j_1 e1,j _ e°,j 2 h 'e"J e"_1,j + e"J_1 e"_1,j_ RTSKC ' e"J e"J_1 > '

_ D- '

ЯТБКГ т Н (20)

^2 сп,, сп,,-1 =-р сп,, сп-1,, ; (ю) Используя обозначения (12), окончательно полу-ЯТ;КГ т И чаем для одноатомного газа

2^ ^ e°,j _сы± _ d eu _ e°,j (Ч + l)e°,j _ q,j _ ((^ _ l)e°,j+ j) _ °; (21)

-2 LC

2V c —c ., c — c , .

^ 2 c n'] n'] = р n,J n—1.j (ii)

TSKC n T h и для двухатомного газа

Cn_1,j _ (2A + 1K,j + ((2A _ 1)en,j_1 + en_ 1,j_1 ) _ ° (22)

Обозначим константу А1 для одноатомного и А2 А2с02, + с0, -с1, - (А2с02,-1 -с0,-1 + с1,-1) = 0; (23)

для двухатомного газа (здесь V- это или У1, или У2):

V И 1 2У И 1 Сп-1,, - А2<, - Сп,, + (А2Сп2,,-1 - Сп,,-1 + Сп-1,,-1 ) = (24)

А1 =---и А2 ^-г--. (12)

1 ЛТЖГ т р ЯТБКС т р

Полученные уравнения (13) - (14) и (21) - (24) Окончательно получаем разностные уравнения могут быть встроены в известные разностные схемы для 0-го и п-го слоев, соответственно, для одноатом- для краевых задач. В нашем случае лучше всего ис-ного и для двухатомного газа: пользовать неявную схему для параболических урав-

нений [1]. Первое и последнее уравнения этой схемы

и

International Scientific Journal for Alternative Energy and Ecology № 10 (90) 2010

© Scientific Technical Centre «TATA», 2010

заменяются на соответствующие уравнения (13) - представляют собой якобиан. Такая система состав-(14) и (21) - (24) и берутся для ' = 0 и ' = п, а среднее ляется относительно поправок X,, выступающих в уравнение схемы берется при ' = 0, 1, ..., п - 1. На- роли неизвестных и имеющих вид

пример, по первому варианту для водорода разностная схема моделирования будет выглядеть так:

X = о™ - с£Л

(26)

A2C2,j (A2C0, j-1 1)с0, j C1, j = 0,

L h2 \ h2

c■ -i,j-l2+DtI c'-j + j + DTj-i =0

i = 0;

(25)

Например, из системы (25) получается

(4 (2c0kj-1) - со,j-i) +1)Xо - Xi = -F?-1), i = 0;

i = i,2,..., n -1;

Cn-i, j - A2< j + (A Cn, .-i - i)Cn, . = 0,

X'-1 "[2 + ^X.- + X,+1 =-РГ\ ' = 1,2,...,п-1; (27)

Хп-1 -(с— -— + 1)ХИ = ' = п.

В случае одноатомного газа по обоим вариантам Эта система линейных уравнений легко решается имеем системы линейных уравнений с трехдиаго- методом прогонки. Таким образом вычисляются X,

нальной матрицей, которые решаются стандартными из которых по (26) определяются искомые с« . То методами прогонки. '1

Для водорода оба варианта приводят к системе же для вт°р°го варианта моделирования: нелинейных уравнений, например, (25), которая решается методом Ньютона. Обозначим за левую часть ,-го уравнения системы (25). Тогда по методу Ньютона можно построить систему линейных уравнений, правые части которых представляют собой

функции -Fi предыдущей итерации, т.е. -Е(к-1), где Xи-1 - (А^*.4 + 1)Хп = -1)1, ' = п. к - номер итерации, а коэффициенты правой части

(2AOk-4 +1)X0 - Xi =-F(k-i), i = 0; h2

X,. ,4 2+

Dt

X. + X,+1 = -F(k-1),i = 1,2,...,n -1; (28)

Моделирование диффузии при коэффициенте диффузии Б = /(с)

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

¿c=_!( D(c)— dt dx I dx

(29)

общие аналитические решения для которого неизвестны. Среди численных решений наиболее эффективной и абсолютно устойчивой схемой является линейная схема Самарского [2, с. 281]:

-(D(c2 .-1) + 2Dfaj-1) + D(c0 .-1) + 2h2 /т)Ci,. + (D(c2 j_i) + D(cu-1))c^ = = (-2h2/т)c,,j-1 -(D(q,.-1) + D(c0,.-1))c0J, i = 1;

(D(cuj-1) + D(c-1,j-1))c-1,j -(D(ct+1,j-1) + 2D(c,.-1) + D(ci-1,.-1) + 2h2 /т)c, j + (30)

+ (D(c+1,.-1) + D(c.j-1))c+1,j =(-2h2/т)c.j-1, i = 2, ..., n-2;

(D(c„-,j-1) + D(cn_2,j-1))cn-2,j - (D(cnj-1) + 2D(cn_x + D^,j-1) + 2h2 /т)c^,j = = (-2h2 /т)cn-1,j-1 -(D(cn,j-1) + D(cn-l j-l))cn j, i = n -1.

(Там же [2, с. 282] приведена нелинейная схема, которую также можно использовать в дальнейших выкладках.) Как и в случае D = const, в эту схему могут быть встроены уравнения для входной (i = 0) и выходной (i = n) сторон пластины, полученные выше. Для D = f(c) эти уравнения примут следующий вид. Обозначим константы, соответственно, для одноатомного и для двухатомного газа:

V h

RTSKr t

A

2V h RTSKC t

(31)

(32)

Вариант 1 (дифференцирование). Формулы (13) - (14) преобразуются, одноатомный газ:

f

A

А

yD(c0J)

f

- +1

c0,j C1,j

A

D(c0,,) 0,j

co, j-i =0;

n_1, j

A

yDK j )

-+1

-o, j>

A

Cnj + D(c„j,) j

c„, j-1 = o;

двухатомный газ:

Л D(cÜJ)

-o, j

A

f А А

A2 c _ 1

I D(co,j) co,j-1 1 f

co, j _ c1,j =o;

A

VD(c», j )

2 c„, j _1 _ 1

c . = o .

n, J

" б(сп -) '

Вариант 2 (интегрирование). Формулы (21) - (24) примут вид: (одноатомный газ)

(2А, + Б(Со,- ))Со,- - б(с0 - )сх] -((2А, - б(с0 -_,)) с,,--, + б(с0 --1)с1--,) = 0 ;

б(ся, - )ся-1, - -(2 А1 + б(сп j)) ся, - +((2 А - б(сп j-1)) ся, --1 + б(сп jj-1 ) = 0 и (двухатомный газ)

А2 С2, - + Б(С, - )С, - - Б(С, - )С, - - (А2 С2, --1 + Б(С, --1 )(С, --1 - С, --1)) = 0 ; Б(ся,- )ся-1,- - А2<- - Б(сп j )сп j + (4<--1 + Б(сп j-1)(ся-1,--1 - ся,--1)) = 0 .

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

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

Приведем разностные схемы и соответствующие системы линейных уравнений метода Ньютона для разных случаев.

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

Л

А

d(c0j )

- +1

co,j c1,j

А1

D(co, j )

co, _1 = °i =o;

(D(cuj_ 1) + D(ct_1,j_l))c,_ 1,j - (D(ct+1,j_1) + 2D(c,j_1) + D(c,._1,.-1) + 2h2 /t) . + + (((c+1,j_1) + D(c,.j_1))c,.+1,j =(_2h2/t) c,.j._1, i = 1, ..., n _ 1;

(41)

n _1, j

A I A

1 +1|c„,.+—^— cn1 = o; i = n.

D(cnj)

D(cn, j )

A

D (co!.-15)

+1

Xo _ X1 =_F(k-1), i = o;

(D(c. j_1) + D(c_1,j_1)))_1 _(d(c+1,j_1) + 2D(c. j_1) + D(c_1,+ 2h2 /t)) + + (D(cj + DJ) = _ F(!-1), i = 1, 2, ..., n _ 1;

Xn_1 _

A1

VD(cnJ')

- +1

X =_F(k-1), i = n.

n n

(42)

International Scientific Journal for Alternative Energy and Ecology № 10 (90) 2010

© Scientific Technical Centre «TATA», 2010

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

А

D(c0J)

2 c2 -

S, j

í А ^

V D(Co,j)

2 Со,j-, -1

co,j- cu = 0 i =0;

(D(c. j-i) + D(c-i,j-i))c,.-i,j - (d(c,+i,j-i) + 2D(c¡ j-i) + D(c¡j-i) + 2h2 /т)c,j + + (((c,+i,j-i) + D(c, j-i))c,.+i,j = (-2h2/т)c,.,j-i, , = i, ..., n -1;

A

í A

D(c„,) V D(cn,,)

cn,j =0. = n;

(43)

A

„(*-i)^ Г - co, j-i) + i

X - Xi =-F0k-i), i = 0;

vD(c0D

(d(c,j-i) + D(c,.-i,j-i)))-i - (d(c,+i,j-i) + 2D(c,.,j+ D(c,-i,j-i) + 2h2 /т)) + + (d(c+i,j-i) + D(c,j-i))x,.+i =-F(k-i), i = i, 2, ..., n -1;

Xn-i -

A

K.D(c{nl;)

(kV2^7 -cn,j-i) +i

Xn =-F(k-i), i = n.

n n

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

(2А, + D(c0 j)) c0,. - D(c0 j - ((2( - D(c0 j)) c0,.-! + D(c0 j)с,, ..) = 0, i = 0; ((.,j-,) + D(c¡-l у-1))с,.-1,. - (((с,.+1,._,) + 2D(c,.ч) + D(c¡_l уч) + 2А2 /т),. + + ((,.-1) + D(c¡,j-l))c¡+l,j = (-2Н2 /т)с,..-,, , = 1, ..., п -1;

D(cя,. )Сп-1,. -(2( + D(Cn,.))ся,. + (( - D(cи,,-1))).-1 + D(Cn,.-1)^,.) = 0, , = п; (24 + D(c0*-1))) Х0 - D(c0^:1))X, = -, = 0;

(((с,,.-1) + D(c¡-l,.-1))х,.-1 - ((, .-1) + 2D(c¡.-1) + D(c,.-1,.-,) + 2И2 /т)),. + + ((+1,.-1) + D(c,j-l))x+l = -^-1), . = 1, 2, ..., п-1; D(cnk,^Хп-1 -(2А, + D(cj))Xn = -^к-1), . = п.

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

«. + ЩС0,. )С,. - D(Co,. )С1,. -(А2С02,.-1 + Щс0г.-1)(С,,.-1 - с,.-1)) = а . = 0; (((С,.-1) + D(c¡_l.-1))с,-1,. - (((С+1,.-1) + 2D(c¡.-1) + D(c¡-l.-,) + 2А2 /т)с,.. + + (+1..-1) + D(c¡,j-l))хl,j =(-2А2/т)е.,.-1, ¡ = 1, ..., п-1;

. )Сп-1,. - А2си2,. - D(Cn,. )Сп,. + ((2си2,.-1 + D(Cn,.-1 )(Сп-1,.-1 - Сп,.-1)) = 0. I = п;

(2А2С0+ D(c0k71)))Х0 -D(c0kj-1»)Xl =-к-1), ¡ = 0;

(,.-1) + D(c¡-l,]-1)))-1 - ((+1,.-1) + 2D(c¡.-1) + D(c¡-1,.-1) + 2Н2 /т)) +

+ ((,.-1) + D(c¡,j-l))х+l =-^(к-1), ¡ = 1, 2, ..., п -1; D«:1))Xn_l -(2А2С^.-1' + D(c;;kj-«))xи = . = п.

Системы уравнений (42), (44), (46) и (48) составляются относительно поправокX. (см. (26)).

(44)

(45)

(46)

(47)

(48)

Заключение

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

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

десорбции и проницаемости, при постоянном коэффициенте диффузии и при произвольной зависимости его от концентрации диффузанта - могут быть охвачены одной компьютерной программой. Однако для случая D = const отдельная программа ускоряет расчеты.

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

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

1. Тихонов А.Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1977.

2. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.

International Scientific Journal for Alternative Energy and Ecology № 10 (90) 2010

© Scientific Technical Centre «TATA», 2010