CC BY
44
Труды Карельского научного центра РАН № 3. 2010. С. 72-82
УДК 519.6: 539.2
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДЕСОРБЦИИ ВОДОРОДА С ЦИЛИНДРИЧЕСКОЙ ПОВЕРХНОСТИ
Н. И. Родченкова, Ю. В. Заика
Институт прикладных математических исследований Карельского научного центра РАН
Рассматривается дегазация цилиндрического образца, содержащего растворенный водород. Эксперимент проводится методом термодесорбции. В соответствующей краевой задаче с нелинейными динамическими граничными условиями учтены физико-химические процессы в объеме и на поверхности металла: диффузия, десорбция, растворение и захват дефектами. На основе разностных аппроксимаций разработан вычислительный алгоритм моделирования десорбционного потока, приведены результаты численных экспериментов.
К л ю ч е в ы е с л о в а : термодесорбция водорода, математическое моделирование.
N. I. Rodchenkova, Yu. V. Zaika. NUMERICAL MODELLING OF HYDROGEN DESORPTION FROM CYLINDRICAL SURFACE
Degassing of a cylindrical sample containing dissolved hydrogen is considered. The experiment is made by the thermodesorption method. In the corresponding boundary-value problem with nonlinear dynamic boundary conditions physical-chemical processes in the bulk and on the metal surface are taken into account: diffusion, desorption, capture by defects, and solution. Computational algorithm for desorption flux modelling is developed on the basis of difference approximations. The results of numerical modelling are presented. This work was supported by the Russian Foundation for Basic Research (grant 09-01-00439) and by the Russian Science Support Foundation.
K e y w o r d s : hydrogen thermodesorption, mathematical modelling.
Введение
Значительная концентрация водорода в металле приводит к водородной хрупкости [Кунин и др., 1972; Колачев, 1985]. Естественные металлургические концентрации растворенного водорода составляют от 0,1 до 100 ppm. Для измерения концентрации водорода в твердой про-
бе в условиях заводской лаборатории авторами статьи [Полянский и др., 2006] разработан анализатор водорода (АВ-1). Цилиндрический образец помещается внутрь вакуумного экстрактора из кварцевого стекла. Экстрактор помещается в печь с заданной температурой экстракции. Контакт образца и стенок экстрактора точеч-
ный, теплопроводность кварца пренебрежимо мала, поэтому теплопередача происходит за счет излучения. При нагревании образца атомарный водород диффундирует внутри и десорбируется с поверхности в молекулярной форме. С помощью масс-спектрометрического анализатора водорода фиксируется экстракционная кривая, подлежащая дальнейшей обработке (в частности, оцениваются кинетические параметры моделей). График зависимости десорбционного потока от температуры при монотонном нагреве (ТДС-спектр) обычно содержит несколько пиков. Наряду с диффузией лимитирующими факторами являются поверхностные процессы (следуем работе [Габис и др., 1987]) и захват атомов водорода различного рода дефектами (например, трещины, микрополости, включения гид-ридных фаз). Решение проблем водородного материаловедения, особенно это касается изотопов дейтерия и трития, требуют значительных затрат. Поэтому роль математического моделирования в таких задачах является достаточно весомой. Работа посвящена математическому обеспечению экспериментальных исследований. При численном моделировании десорбционного потока будем использовать параметры, характерные для алюминия и его сплавов.
Математическая модель
Уравнение нагрева. Образец имеет форму цилиндра с характерными размерами
(ГОСТ 21132.1-98): радиус основания
L = 4 10-3 м, высота H = 2 10-2 м . Концентрация в начальный момент времени с(0, r, z) = с постоянная (формируется в процессе изготовления материала). Без принципиальных изменений численного алгоритма при необходимости можно учесть снижение концентрации водорода в приповерхностном слое, например, в результате предварительной механической и термообработки. Если прогрев образца равномерный (достаточно медленный, T = T (t), [T ] = К), то динамику изменения температуры можно описать дифференциальным уравнением (ДУ) теплового баланса [Полянский и др., 2008]:
— = — 7 10-5(T + 64.3)[Te4 -T4]. dt с pV
T0 = T(0) = 293 K. Здесь Te = const - температура стенки экстрактора, S = 2nLH, V = nl}H - площадь поверхности и объем ци-
линдра, р = 2,71 103 кг ■ м-3 - объемная плотность, ~ = 1,15 ■ 103 Дж ■ кг-1К-1 - удельная теплоемкость, а = 5,67 10-8Дж ■ с-1 м-2К-4 - постоянная Стефана-Больцмана.
Поскольку не при всех Ге предположение о
равномерном прогреве справедливо с достаточной точностью, рассмотрим альтернативную распределенную модель. С учетом «трубчатой» геометрии экстрактора считаем, что нагрев идет в основном через боковую поверхность. Тем самым речь идет о нижней оценке динамики прогрева «центра» образца. Приведенное уравнение из [Полянский и др., 2008] является мажорантным сверху. Примем радиально симметричную модель:
Г(0,г) = Т0 < Г^ ге [0,L],
дГ
дг
dT (д2T 1 dT
—=к —-+------------
dt I dr r dr
= 0,
Я
dT
dr
= a[Te4 - T4(t,L)].
Здесь А - коэффициент теплопроводности (для алюминия в диапазоне Г е [300,800] А« 236 Дж ■ с-1м-1К-1); к = А(~р)-1 - температу-
*
ропроводность; время окончания счета t определяется стационаром Т(/,0) ~ Те, / > /*; а = ае,, £1 = 7 ■ 10-5(Г + 64,3) - коэффициент поглощения.
Вспомогательная задача численного моделирования нагрева: оперативно оценивать, насколько распределение Г(/, г) отличается от равномерного нагрева Г(/) при заданных Те, L, Н и теплофизических характеристиках материала. Например, при Г0 = 293К, Ге = 773К и указанных
L, Н предположение о равномерности нагрева алюминиевого образца справедливо: разность Г(/) - Г(/, 0) не превышает сотой доли градуса. И это без учета в распределенной модели поглощения тепла торцами. Максимальная температура достигается через 2,2 часа, причем в пределах часа нагрев практически линейный. В дальнейшем считаем прогрев образца равномерным.
В качестве завышенной оценки рассмотрим случай, когда торцы прогреваются столь же интенсивно, как и боковая поверхность: £ = 2жLH + 2л 13. На рис. 1 приведены графики температуры: верхний - для ДУ с учетом прогрева всей поверхности, нижний (два совпадающих графика) - для ДУ с учетом прогрева только боковой поверхности и для краевой задачи теплопроводности в бесконечном цилиндре.
L
Рис. 1. Зависимость температуры образца от времени
Диффузионная модель с учетом дефектов. Рассмотрим краевую задачу термодесорбции для цилиндра с учетом диффузии в объеме, захвата водорода дефектами, выхода из раствора на поверхность и десорбции:
Эс _ (д2 c 1 дc д2 cЛ
-=d(t )
+----*------+ '
dr2 r dr dz2
[l - W ]c(t, r, z) + a2 w(t, r, z) + a3 y(t), (1)
dw
— = a1(T)[1 - W] c - a2(T) w, W = wwi
dt
r є (0, L), z є (0, H), t є (0, t*),
dY = -a3Y(t), Y(tcrit) = Y ^ Y(t) = dt
= fexp{-a3(t-tCrlt)}, a3 = ° T < TCrlt,
c(0, r, z) = c = const, w(0, r, z) = w =
-1
max 5
(2)
ac
^ -1 , rє [0,L], zє [0,H],
dc_
dr
=0, *
dz
= 0, D(t) = D(T(t)),
b(t) = b(T(t)), £(t) = g(T(t)), d^tL(t, z)=-b(t) q^ z) - d (t) |c
c(t,L, z) = g(t)q1(t, z), z є [0,H], dq2 dc
3t
-(t, r) = -b(t) q2(t, r) - D (t) —
c(t,r,H) = g(t)q2(t,r), rє [0,L], dq3-(t, r) = -b(t) q^(t, r) + D (t) —
(3)
(4)
at
dz
c(t, r,0) = g(t) q3(t, r), t є [0, t*],
dT = 7• 10-5(T(t) + 64.3)Te4 -T4],
dt cpV
T (0) = T0 = const. Здесь c(t, r, z) - концентрация атомарного водорода (H) в металле; w(t, r, z) -концентрация H, обратимо захваченного дефек-
тами физико-химической структуры материала (например, в микрополостях); wrnax - максимальная концентрация обратимого захвата; ^(/) - концентрация Н в ловушках, которые начинают высвобождать водород только по достижении некоторой критической температуры Гсгк (характерно для включений гидридных фаз); а1 - коэффициенты поглощения и выделения Н
ловушками (а3 > 0 при Г > Гсгй); д1(/, г), д2(/, г), д3(/, г) - поверхностные концентрации
(на боковой поверхности цилиндра и на торцах); с - начальная (естественная) концентрация Н в твердой пробе; g - коэффициент соответствия концентраций атомов водорода в объеме и на поверхности (коэффициент быстрого растворения); В, Ь - коэффициенты диффузии и десорбции. В рассматриваемом температурном диапазоне (Г е [300,800]) полагаем аг > 0 константами. Изменения для случая нестационарных а1 (/) = а1 (Г (/)) непринципиальны. Для
практических целей захват учтен в простейшей интегральной форме, уточнение геометрии дефектов и их распределения в образце существенно усложнило бы модель. Для упрощения выкладок учли по одной обобщенной ловушке обратимого захвата и распада. Их может быть несколько со своими индивидуальными а) и
. Считаем, что коэффициенты диффузии и
десорбции зависят от температуры по закону Аррениуса:
В(Г) = В0ехр{-Ев /[ЯЛЬ
Ь(Г) = Ь0ехР{-ЕЬ
ЕВ Ь - энергии активации. Сокращенно обозначаем В (/) = В (Г (/)), Ь(/) = Ь(Г (/)). Отрезок времени [0, ] определяется дегазацией:
J(/) ~ 0, / > . Условия сг (/,+0, г) = 0,
сг (/, г, Н /2) = 0 следуют из симметрии.
Более точная модель растворения на поверхности (для определенности боковой) имеет форму баланса потоков:
к +(Г) с(^ L, г )[1 - д1(/, г ^ - к- (Г) ql(t, г) х х[1 - ^, L, 2 ^ = - В(Г) сг (/, L, г).
Но когда диффузия значительно медленнее растворения (температура не слишком низкая) и концентрации малы, получаем условие быстрой растворимости с ~ gq, где g = к- / к + . Если поверхность изотропна (в смысле Ек_ ~ Ек +), то
z=H/2
r=+0
r = L
z=H
z=0
параметр g слабо зависит от T . В дальнейшем обозначение Eg используем условно: это не энергия активации, а разность Ек- — Ек +, которая может оказаться и отрицательной.
Что касается ловушек, активирующихся лишь с определенной критической температуры (типа гидридных фаз): учли лишь их емкость и скорость распада. Моделирование дегидрирования - самостоятельная сложная задача, приводящая к нелинейным краевым задачам со свободными границами раздела фаз с условиями типа Стефана. Вариант с объемной десорбцией рассмотрен в [Zaika, Rodchenkova, 2008], с поверхностной - в [Zaika, Rodchenkova, 2009].
В силу симметрии начальных данных q3 = q2, и далее разностную аппроксимацию строим лишь для половины цилиндра (z е [H /2, H]) с соответствующими краевыми условиями (cz |H/2 = 0, q2 = ...). Для дефекта с
обратимым захватом (микрополость) константа W определяется в силу (2)
(д /dt = 0): «1 (T0)[1 - W/ Wmax ]с - a2(T0)W = 0.
Для ловушки типа включения гидридной фазы значения у = const, Tcrit, a3 задаются по информации о конкретном химическом составе гидрида. Наличие производных qt (накопления) соответствует представлениям о возможности миграции H по поверхности до десорбции H 2 .
Таким образом, модель нацелена на задачу анализа в целом динамики малых естественных концентраций водорода (без предварительного насыщения в лабораторных условиях). В прикладном контексте задачи основные процессы учтены минимальными математическими средствами. Дальнейшая детализация приводит к росту количества параметров. Обратная задача их оценки становится труднообозримой.
Уточнение постановки задачи. Цель работы состоит в разработке разностной схемы и вычислительного алгоритма для моделирования десорбционного потока водорода из цилиндрического образца:
J(t) = 4nb(t)^LI / q12(t,z)dz + j rq^(t,r)dr'j.
Десорбируется H2, но подсчет ведем в атомах ([J] = 1/с). Критерием правильности вычислений выбран материальный баланс:
П2Н (с + ^ + /) +
+ 4п^£| / q1(0,г)йг +| гq2(0,г)йг
pH pL
= 4п| йг I г[с(/, г, г) + м>(/, г, г) + т(?)]йг +
* Н / 2 * 0
+4п(11Н/2 ql(t, г) йг+10 г q2(t, г) йг)+10 ^т)йг.
При монотонном нагреве удобно наряду с зависимостью от / рассматривать ТДС-спектр - график J = J(Г). Обычно он содержит несколько пиков. Считают, что первый пик соответствует поверхностному водороду. Но следует соблюдать осторожность: пока десорбируется поверхностный водород, идет подкачка Н из объема. Актуальна задача оценки соответствующей поправки. При моделировании время ^ окончания поверхностного водорода определяется из
|0^(т) йт~ 4ж^Ь | 2 q1(0, г) йг+£ гq2(0, г) йг).
Численное моделирование позволяет выделить на экстракционной кривой участок, соответствующий количеству начального поверхностного водорода, оценить значения энергий активации диффузии, десорбции, параметры обратимого захвата и распада гид-ридных фаз.
Разностная аппроксимация краевой задачи
Следуя стандартной методике [Самарский, 1971], введем пространственную сетку
(r, z3): rt = ihr, i = 0,1,...,N1 = [L/hr];j z = *, j = 0,1, ...,N2 = [(H/2)/hz] і
jz
и сетку по времени
rnT=[tk = kr, к = 0,1,™,K = [t*/т]}.
Обозначим через ски приближенные значения объемной концентрации с(/к,г,г.). Анало-
^(4, г, г,), / « 7«к), ак = а, (/к),
гично
w..
Бк = Д/к), где (г, г,.) є Ц,, їк єй)т. Для уравнения
(1) рассмотрим неявную разностную схему метода переменных направлений, называемую продольно-поперечной (схемой Писмена-Рэч-форда), а для уравнения (2) - схему с весами. Переход от слоя к к слою к +1 осуществляется в два этапа. На первом этапе определяются промежуточные значения с^12 из системы уравнений
c+1/2 -c
i,_________________V _ Г)
— л-'Ь
Ґck+1/2 -2ck+1/2 +ck+1/2 1 c+1/2 ^
i+1,j________ij_________i~1J і 1 i+1j i-1,j
0,5r
k+1/2
Ж
+D
c* -2ck +c*
. hi ij1 -a+1/2[1-+ (5)
+ak+1/2wk+1/2 +0/+1
wk+1/2 - <j
0,5 t
= (1 -a)(of [1 - Wkj ] ck j - a2kwk j)+
(6)
+ a(ak+1/2[1 -r~.*+1/2]ck+1/2 -ak+1/2wf+1/2).
\ 1 L i, J i,J 2 1 ,J '
На втором этапе, пользуясь найденными
k+1/2 k+1
c,+ , находим cfj1 из системы
c*+1 - c*+1/2
iJ iJ
0,5 t
= Dk+1/2 X
f k+1/2 r, k+1/2 . k+1/2 , k+1/2 k+1/2 \
Ci+1„ - 2ci j + ci._1,j + 1 cM J - ci._1,j
-----------------1-------
i+1,j
V
h
2h
(7)
J
+ D„
ck+1 - 2ck+1 + ck+1
Ч/+1 cu j -a^p-w^ci +
h
+ a2k+1wk+1 + aY+1,
w"+1 - 1/2
i,J i,J
0,5 t
= (1 -a) x
(+1/2 [1 - ^^к+1/2]ск;1/2 - а2к+1/2 .)+ (8)
+*(+1[1 - (/с;1 - а2+4;1).
Здесь, чтобы иметь возможность использовать алгоритм прогонки на 5 -м слое по времени (5 = к +1/2; к +1), неизвестную величину Ж*.
*, 3
заменяем ее аппроксимацией из линейного по Ж". уравнения
w°j=w:1/2+4{ar1/2[1 - w:-1/2k-i/2 wma
- a2:-1/2 W :-1/2 + a: [1 - Wt:. - asws }.
(9)
A = 1 - hr (2r )-1, B, = 1 + hr (2rt )-1,
r = 2h2r1, Gkj+1/2 = 2 + ^D-1/2[1 + VjU2], V*= af[1-W*] ~+1/2
a^ ‘ +4t
-1 i,j
2a*+1/2[1 - Wk+1/2]
N, = (4T-1 - a2k )w* j
i,J 4т-1 + a2k+1/2 ’
Fk = D*
i,J D
■(hr /h)2( -2c^kj + <,-1)+ a) x x[1 + 0,5Tak+1/2Vkj ]j + h2D-+m[o.+1/2N.j + arf+1/2],
при каждом фиксированном j = 1,2,N2 -1 получаем при к > 0
A.c"+^/2 - Gk++1/2 + B.ck,+1/2 + Fi*. = 0. (10)
1 1-1,J 1, J 1, J 1 1+1,J 1, J v /
Значения в начальный момент времени (на нулевом слое) известны: c°. = c = const. Следуя
методу прогонки, ищем приближенные значения концентрации в узлах сетки на (к +1/2 )-м слое ( k > 0 ) по времени в виде
ак^2=<;ч“.2+в;^, і=ід..,N1-1, (11)
Прогоночные коэффициенты: i = 2,3,., N1,
<+1,2 = в,-У(j = - a,a;3):
e:*;11 =(А,_1Д‘;;г + fU/Gt - a,-1ak-11;3).
При r ^ +0 имеем
cr / r = (cr (t, r, z) - cr (t,0, z)) / r « crr.
Начальные коэффициенты находим из аппроксимации уравнения
ct = D(2crr + czz) - a1[1 - W ] c + a2w + a3y на (к +1/2 )-м слое, і = 1, и условия cr|+0 = 0:
<+1/2 = 1 -г (4 D*+1/2 )-1 [1 + V~k;+1/2], ви = FJ/4.
Ближайшая цель - найти значение c
В стандартных обозначениях
у = / (/, у) - это симметричная схема
У5 = У5-1/2 + {/5-1/2 + /5}Т4 для нормированного на ^тах уравнения (2) с фиксированной по времени функцией с = с(/5-1/2,г,г). Итерационную процедуру уточнения Ж*. укажем позже. Для
г, 3
определенности полагаем с = 1/2. Погрешность аппроксимации есть 0(т2 + Аг2 + Аг2) [Самарский, 1971].
Прогонка по радиусу г. Рассмотрим переход с к -го слоя на к +1/2. Выразим мгк+х12 из уравнения (6) и подставим в (5). В обозначениях
N1, j
не-
обходимое для реализации прогонки. Запишем аппроксимацию граничного условия (3) (г = L): к > 0, (я*;1/2 - )/(0,5т) =
= 0,5[-Ьк (qlк. )2 - Вксг (/,, А г.) - Ьк ~
- Вк +1/2сг (/к +1/2, А г. )]. (12)
В граничном узле с точностью до 0(Л2) имеем: га = к, к +1/2,
2Ас (/т, 4 г.) « с” -2,. - 4с” -1,. + 3с”,.. (13)
Значения концентрации на к -м слое уже известны. Для (к + 1/2)-го слоя, подставляя 4+-2. и с^+-12. из соотношения (11), имеем
+
x
+
с. «к ,і,2. і. ,,)»2- [(з ( кй'і - 4)
в;+к-12, - 4» ].
Запишем приближение компактнее:
+
л-з+<;;=(ок+1,2- - 4), в=вкн,2- 4)вк;
•к+1/2 +в.
2кгСг(К+;12,Ь, 2;) « ЛС^2 + В . (14)
Соотношения (13), (14) подставляем в (12),
N1,; = у :
обозначив ск+1/2
< = ^т1ст1,,. т = к,к +1/2 .
к +1/2 у2 +
ёЫ/2
4
Тёк+1/2 (
- +
Ак+1/2 Л 2И„
у + Г = 0.
ё*
зр
4
2К
Ч,; +
к )
Теперь найдем все значения сі
к+1/2
при
; = о і
условие
; = N2, і = 0,1,..., N1. Используя
симметрии
С, \и/2 = 0.
имеем
„к+1/2
= (4сі + - сі + ) / 3. Значения с
к+1/2
> о
(т, И малы) однозначно определяются из квадратного уравнения, аппроксимирующего граничное условие (4) при 2 = и :
ёк+И2
(с*;1;3 )2
+
2
А.
т&+1/2 к+1/2
за*
- 4с. ;
2И ' і,N2 -2 і,;
2И,
2
: ) 2 „к 1 тёк
с*+:12+
с'^ = 0.
/,N2
Зная все значения с; + , вычисляем концен-
трацию ^ :
= (4т 1^5 + а,[1 -Ж*,]с,, -аХ^,, + 1+1/2 [1 - Ж/"2] с*+1/2) (+1/2 + 4т-1 )-1.
Прогонка по переменной 2. Поскольку в цилиндрических координатах возникает особенность при г ^ +0, то переход с (к + 1/2)-го слоя на (к + 1)-й совершается в два этапа. Первый этап: і = 1, г ^ +0 (аппроксимируем
Сг / г ~ Сгг), реализуется алгоритм прогонки для уравнения
с = В(2сгг + сгг) - а1[1 - Ж] с + а2и + а3у . Второй этап: I = 2,..., Ы1 -1, г > 0, прогонка для уравнения
с( = В(сгг + сг / г + сгг) - а1[1 - Ж] с + а2и + а3у.
Уравнение для и не меняется. Уравнение (7) для 1 = 1 принимает вид
к+1 к+1/2
С1,, - С1,,
0,5 т
= 2А
к+1/2 т к+1/2 . к+1/2
2,; - 2с1,, + С0,,
к +1/2
+
1,;+1
к+1 1,1-1
(15)
- а1+1[1 -ЇЇкТС; + а^] + а3Ук+1.
Выразим ^1+ из (8) (і = 1, <7 = 1/2) и подставим в (15). В обозначениях
г = 2И2т-1, Ок; = 2 + А [1 + V] ],
+ 2И-(А+1/2В + а*(,1 -4с^-1,,))
Корни квадратного уравнения по у разных знаков (т<< 1). По физическому смыслу берем у = с^1;2 > 0 . Погрешность аппроксимации граничного условия О (И + т2), что согласуется со схемой в объеме.
Vх
а\+1/2[1 - Ж*]112] ак2 +1 + 4т-1 :
.к+1 2а,+1[1 - ж^]
Nk +1 / 2
_-1 к+Шч к+1 к+1/2
г - а2 ) а2 ^1]
'1,;
4т-1 + а2+1
2А,
'к+1/2 ^-^к+Ш
А,
(и, / н, )2(4+» - 2с-;+'" + с0+‘'2)+
;rА-:;[; + 0.5га2 +1^к,+1'2 ]с1 +И,2 Аí:;[NÍ+“2 + а>Г\
■]
к+1/2
получаем
с^^ - G1k;1c1¿;1 + c1k;+; + ЕХ]112 = 0, к > 0. (16)
Ищем приближение концентрации на ( +1 )-м слое по времени в виде: > 0 ,
<+’=«Х;+1 с?|+1+дк А.;=0,1,., N г -1, (17)
Прогоночные коэффициенты: ; = 2,3,., N 2,
1 = -1 -а*+-1 )-1.
Л*;1=(-1+Ак+11/2)(- ax;-l )-1.
Начальные коэффициенты находим из (16) при ; = 1 и условия с2 \и/2 = 0 :
< = 1 -к(2Ак+1П1 + ^у*;1],
Ок+1 77к+1/2 / о
Р1,1 = Е1,1 /2.
Разностная аппроксимация уравнения (7) для і = 2,.,N -1:
2
И
г
2
И
2
1,0
ск+1 - ск+1/2 і,І і,І
0,5 т
= Ак+1/2 Х
Ґ к+1/2 ^ к+1/2 „к+1/2 , „к+1/2 Х+1/2Л
Сі+1,; - і,; + Сі-і,; , 1 Сі +1,; - Сі-і,;
----------1------
И
2
2И
(18)
+А,,
с^++1 - 2ск+1 + с^+11 1,| +1________1,| 1./-1
Ик
- а,+1[1-Ж Х/+1]#1 + а,+У]1 + аз/+1.
Выразим ^кн/1 из (8) (а = 1/2) и подставим в
і,]
соотношение (18). В обозначениях г = 2ИУ1, + = 2 + КА-+1 [1 + V/ ],
V
к +1/2
а.*+1/2[1 - Ж*+1/2] а,+1 + 4т-1 :
„ к+1 = 2а1+;[і - ж;*;1]
і,] а2к+1 + 4т-1
N7
(4т- - а2+1/к)а2к+1 wk; 1/2
4т-1 + а*+1
А
к+1/2 ;„2
А,,
Ик х
х
/ к +1/2 т к +1/2 . к +1/2
^1,/ - 2с,,/ + Сі-1,/
ик
1 с* +1/2 - ск +1/2^ 1 і+1, / і-1, /
+ гАХ+1 [1N0.5то2 ;;Vik;*,'! ]с и, А.-11 [N1,+1/2 + а,/ +1]
г
к+1/2 ;.;
2И„
при каждом фиксированном і = 2,3,
, N1 -1
получаем
с к +1 "і, /-1
с: /.‘1 - ОХ]Іск,+ + ск++; + ек+1/2 = 0, к > 0. (19)
с к +1 ч;+1
-к +1/2 1 ,;
цель -
значение
і=1,., N1 -1. необходимое для реализации прогонки. Запишем аппроксимацию граничного условия (4)
(2 = И): к > 0, -дкі+1/2)/(0,5т) =
= 0,5[-Кк +1/2(<?2і+1/2)2 - Ак +1/2С,(^к+1/2, гі, И) - (21) -Ьк+1(^+1)2 - А*+1 с,(^к+1, г,И)].
В граничном узле с точностью до О(Ик) :
2И,с, (V, г, И) - с/
2 2^;^,^' і
т
і,N2 -2
т = к +1/2, к +1. (22)
Значения концентрации на (к +1/2 )-м слое уже известны. Для (к + 1)-го слоя, подставляя значения <&;-2 и ск^ -1 из соотношений (17)
(при = 1) и (20) (при , = 2,..., N1 -1), получаем аппроксимацию
сг (/к+1, г,, Н)« ^ [(3 + ( (-1 - 4)
+
+ Рк+1 +1 - 4)+1]
Т И 1,ы2 -і т ; -1 ^ /Ач; -1'
Запишем приближение компактнее:
а — з+о*;. -1- 4),
В -вХ+1 -1 ;(akNl -1 - 4)р£+1.
2И,с,(?к+1,г,,и) - Ас,*; + В . (2з)
Соотношения (22), (2з) подставляем в (21), обозначив с*N = у :
= ЯтЧ^,, т = Х + 1/2, Х + 1,
%- у2 +
ё*+1
Г — Кк +1/2
4 . + Ак+1А
тё*+1 2И,
у + Г = 0.
ё¿+1/2
(г )2
за,
2И,
+ — (+1В + (+1/2( 2И,
2
к +1/2 >| Лк +1/2
— 4с
/,;к -2 ^^i',;2 -
тёк +1/2 ) ;)).
Ищем концентрацию в узлах сетки на (к +1 )-м слое по ґ в виде: к > 0,
с“ "И,+; ]=0,1,., N2 -1, (20)
Прогоночные коэффициенты: / = 2,з,., N 2,
а,*; =(Оік*-; -а*+-;)-;,
вк+1 _(в*+1 + Е*+1/2 \(о*+1 а*+1 )-1
Р/,; = \Р/,;-і + Е/,;-1 ДОі,;-1 - а/,;-1 / *
Начальные коэффициенты находим из (19) при / = 1 и условия с, \и/2 = 0 :
аХі+1 = 1 - г(2А*+;)-1[1 + VV,k+1],
рк +1 т-рк +1/2 / г*
рі,і = Е,,1 12.
Ближайшая
По физическому смыслу выбираем положительный корень квадратного уравнения по у . Корни разных знаков, по крайней мере при относительно малых т (обозначаем т<< 1 ). Погрешность аппроксимации граничного условия О(И2 + т2) согласуется со схемой в объеме.
Теперь найдем все значения с,+1 при і = 0 и
і = N., / = 0,1,.,N2. Используя граничное условие на оси цилиндра (сг \+0 = 0), получаем с
ляются из уравнения
-(сІ+І/^
0+] = -ск+1)/з. Значения cN+1; > 0 опреде-
к+1
2
за,
к+1
тёк+1 2Иг
+-
кИ.
(+1 - 4Ск+1 )-
^N1 -2,; ^м.-;,;!
тёк
*+1/к = 0 с;.,;
х
+
4
аппроксимирующего условие (3) при г = Ь. брать для хорошей аппроксимации графика
(± 5 %)? Ряд сходится медленно. Напри-
По ск+1 вычисляем *, 3
1 = {4ГХ+2 + ак+1/ 2[1 - Ж™ ] с;"} -
-а2+1/2М+1/2 + а1+1[1 -Ж^/1] ск+1}((1 + 4т“1).
Коррекция Жж- (5 = к + 1/2; к + 1) состоит в том,
что можно положить ^ = м*з / мт
и повто-
рить вычисления по схеме метода переменных
направлений до установления Ж6. - Ж6, (обыч-
*,. *,3
но две-три итерации).
Классическая диффузионная задача. Для упрощенных расчетов считают диффузию единственным лимитирующим фактором (игнорируются физико-химические процессы на поверхности), и рассматривают краевую задачу I рода с нулевыми граничными концентрациями:
дс ч ( д2 с 1 дс д2 с ^
■ +-------------+ ■
дг2 г дг дг2
ге (0,Ь), 2е (0,Н), ге (0,г*),
сг (г,+0,2) = 0, с(г, Ь, 2) = 0,
с(г, г,0) = с(г, г, Н) = 0, с(0, г, 2) = с.
В качестве приближения решения используется частичная сумма ряда
c(t, г, г) = ££ Апт ехр] - Лпт } Э(Т)йт |
X
(1 • (пп г I • 81п|--------------г
Ь
Н
Апт
4(1 + (-1)п+1) с
2 -1 пп
^"пт I
Н
V Ь У
пп$п ЛСО
представляющего обобщенное решение (из-за несогласованности краевых условий). Здесь В (г) - В(Т (г)), 30, 31 - функции Бесселя первого рода нулевого и первого порядков, ¡^т -последовательные нули функции 3{)(ц). Поток Н сквозь поверхность:
8НВ(г )[1 + (-1) п+1]
1 (г) = = с ЕЕ"
п=1 т=1
пп
1+(-1)п+1
+
(-1)п+12п ([ Н^]-1 )2
• ех
р{...}.
Часто ограничиваются первым слагаемым (п = т = 1), считая остальные быстро затухающими. Но в рассматриваемой задаче асимптотика по времени не столь важна: самое интересное происходит до и в окрестности «верхушки» ТДС-пика. Сколько членов ряда необходимо
мер, при с = 1023 м 3, В0 = 2 • 10 3 м2с ’, ЕВ = 6 • 104 Дж • моль-1 (В - 2 • 10-7, Т = 773К), нужно 144 члена ряда (п = 1,3,5,7,9,11,13,15; 1 < т < 18). При четных п сумма равна нулю.
Рис. 2. Приближения диффузионного потока I (г)
На рис. 2 приведены по убыванию максимума приближения I(г) частичными суммами с
(п, т) < (NМ) = (15,18), (5,10), (3,3), (1,1). Несколькими слагаемыми не обойтись, пик только один (разве что в конкретной ситуации необходимо учесть несколько сравнимых каналов диффузии с различными В(3)).
Алгоритм и результаты моделирования
Изложим поэтапно алгоритм вычислений. Фиксируем значения Ь, Н, В0, ЕВ, Ь0, Еь, g0,
Ег , с , Мтах , 7 , ^, Te, а2 , Ор ^ , ^,
с~ . Переход от к -го слоя к к + 1 осуществляется в два этапа. Ьй этап: на (к +1/2 )-м слое по времени (к > 0) алгоритм вычислений следующий.
1. Вычисляем значения Жгк+1/2 из уравнения
(9) при 5 = к +1/ 2.
2. В соответствии с соотношениями (10), (11) прямым ходом прогонки вычисляем наборы коэффициентов а^+1/2, в+1/2.
*,3 *,3
3. Значения концентрации в граничных узлах определяем, решая квадратные уравнения относительно = 4+у > 0.
4. Обратным ходом прогонки по формуле (11) находим приближенные значения концентрации во всех внутренних узлах.
2
х
п
^+1/2 = жк+1/2
к+1/21 ',к+1/2}(а2+1/2 + 4т“1)-1.
'.У А/
тк+1/2 .
1. У
і.] ' "шах
к+1/2 - Жк+112 (2-3 итерации).
5. Используя граничные условия, доопределяем значения ск+1/2 в граничных узлах при
3 = 0 и у = Ы2, г = 0,1,...,N.
6. Из второго уравнения разностной схемы в объеме вычисляем значения концентрации в дефектах обратимого захвата:
1,2={4г-х.+а [-ж,]<,■ -
1 --ж
7. Коррекция Ж: можно положить
/ мтах и повторить п.п. 2-6 до установления Жк ..
г, 3 г, 3
Н-й этап: на (к + 1)-м слое по времени (к > 0) алгоритм следующий.
1. Вычисляем по слою значения Жгк+ из
уравнения (9) при 5 = к +1.
2. В соответствии с соотношениями (16),
(17), (19), (20) прямой прогонкой вычисляем коэффициенты дк3+1, 1, в3+1.
3. Значения концентрации в граничных узлах определяем, решая квадратное уравнение относительно у = с^ > 0.
4. Обратным ходом прогонки по формулам (17), (20) находим приближенные значения концентрации во внутренних узлах.
5. Из граничных условий: ск+1 при г = 0 и
г = N1, 3 = 0,1,., N2.
6. Из второго уравнения разностной схемы в объеме:
мк+1 = {4т~1мк+1/2 + а1+1/2[1 -Жк+1/2]ск+1/2 -^ г,3 1 *,3
—а0
,.к+1/2,, ,к+1/2
+а\+1 [1—Ж*;1 ] сгк+1 }(а|+1 + 4т“1)
+1 +1
Ж
7.
+1
Коррекция Жг
+1
г,3
можно
+1
2
-—К—1
положить
1./
новления Ж
1 / ^
1. / ' ш
+1
и повторить п.п. 2-6 до уста-
1000 3000
7000 9000 эес
Рис. 3. Влияние параметра Ь0
^)х1е-13
Рис. 4. Влияние параметра Еь
■ВДх 1е-14
Ж
к+1
Рис. 5. Влияние параметра Г0
Результаты численного моделирования
На рис. 3-10 изображены графики де-
сорбционного потока 3(г). Варьируемые коэффициенты приводятся в порядке следования максимумов 3(г) слева направо или по их убыванию. Кружком отмечается окончание начального поверхностного водорода. Фиксируем у = 1023 м-3, мтах = 1023 м-3. Параметры по умолчанию, общие для всех графиков, берутся из таблицы.
Рис. 6. Влияние параметра ЕГ
ад х 1е-14
1000 2000 3000 иес
Рис. 7. Влияние параметра g0
ад х 1е-13
Рис. 8. Влияние параметра Е£
ад, !(Ц х 1е-13
0 1000 2000 3000 1,зес 4000
Рис. 9. Сравнение 3 (£) и I (£)
Рис. 10. Влияние дефектов двух типов
Параметры моделирования, общие для всех графиков
Ъ о 1 1 0 1 ЕЪ = 12 ■ 104 Т0 = 293 а1 = 0,1 Ь = 4 ■ 10-3
40 = 10-6 II 0 -й* Тсгй = 700 а2 = 0,2 Н = 2 ■ 10-2
£0 = 103 Ея = 0 Те = 773 а3 =10-3 с = 1023
Примечание. [£0] = м1, [Ъ0, Ц0] = м2 с1,
[ЕЪ, ЕЦ , Е£ ] = Дж • МОЛЬ~1 , ^ Тсг* , Те] = К ,
[а1, а2, а3] = с “Ч [Ь, Н ] = м, [с ] = м -3.
На рис. 3 показано влияние коэффициента Ъ0 на десорбционный поток. Изменения: Ъ0 = 10-12;10-13;10-14 м2 с-1. Рис. 4
иллюстрирует влияние энергии активации Еь, изменения по убыванию максимума: ЕЪ 10“4 = 9; 10; 12 Дж • моль“1. Рис. 5: Ц0 = 5-10“5; 10-5; 10-6 м2с“1. На рис. 6 показано влияние ЕЦ : ЕЦ 10“4 = 4; 5; 6. На рис. 7 ва-
рьируется g0 = 50; 100;200 м-1. Рис. 8 отражает зависимость десорбции от значения Е£ 10-3 = 1; 0; -1. На рис. 9 представлены графики двух потоков: диффузионного - для вырожденной задачи I рода, десорбционного - для краевой задачи с динамическими граничными условиями (левый график). Значения параметров: Ц0 = 2 • 10-3, ЕЦ = 6 • 104, Ъ0 = 10-5,
Еъ = 12• 104, go = 103, / = 3 10”, = 1024 ,
а10 = 0, а3 = 0. На рис. 10 приведен график потока с учетом дефектов двух типов (обратимый захват и распад).
Работа выполнена при поддержке ОМН РАН, РФФИ (грант № 09-01-00439-а) и Фонда содействия отечественной науке.
Литература
Габис И. Е., Компаниец Т. Н., Курдюмов А. А. Поверхностные процессы и проникновение водорода сквозь металлы // Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. С. 177206.
Колачев Б. А. Водородная хрупкость металлов. М.: Металлургия, 1985. 217 с.
Кунин Л. Л., Головин А. И., Суровой Ю. И., Хох-рин В. М. Проблемы дегазации металлов. М.: Наука, 1972. 324 с.
Полянский А. М., Полянский В. А., Попов-Дюмин Д. Б., Козлов Е. А. Новый измерительный комплекс для абсолютного определения содержания водорода в материалах водородной энергетики // Альтернативная энергетика и экология. 2006. Т. 38, № 6. С. 29-31.
Полянский А. М., Полянский В. А., Яковлев Ю. А. Методы определения энергий связи водорода в твердом теле, реализованные на базе анализатора водорода АВ-1 // Взаимодействие изотопов водорода с конструкционными материалами: сборник докл. III Меж-дунар. конф. (Санкт-Петербург, 2-7 июля 2007 г.). Са-ров: РФЯЦ-ВНИИЭФ, 2008. С. 346-353.
Самарский А. А. Введение в теорию разностных схем. М.: Наука, 1971. 553 с.
Zaika Yu. V., Rodchenkova N. I. Modelling of diffusion TDS-spectrum peak of dehydriding with
СВЕДЕНИЯ ОБ АВТОРАХ:
Родченкова Наталья Ивановна
научный сотрудник, к. ф.-м. н.
Институт прикладных математических исследований КарНЦ РАН
ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910
эл. почта: NIRodchenkova@yandex.ru тел.: (8142) 766312
Заика Юрий Васильевич
зав. лаб. моделирования природно-технических систем, д. ф.-м. н.
Институт прикладных математических исследований КарНЦ РАН
ул. Пушкинская, 11, Петрозаводск, Республика Карелия,
Россия, 185910
эл. почта: zaika@krc.karelia.ru
тел.: (8142) 766312
size reduction and heat absorption effects // NATO Science for Peace and Security, Series C, Carbon Nanomaterials in Clean Energy Hydrogen Systems / Eds. B. Baranowski et al. Springer, 2008. P. 863878.
Zaika Yu. V., Rodchenkova N. I. Boundary-value problem with moving bounds and dynamic boundary conditions: diffusion peak of TDS-spectrum of
dehydriding // Applied Mathematical Modelling. Elsevier, 2009. Vol. 33, N 10. P. 3776-3791.
Rodchenkova, Natalia
Institute of Applied Mathematical Research, Karelian Research
Centre, Russian Academy of Science
11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: NIRodchenkova@yandex.ru
tel.: (8142) 766312
Zaika, Yury
Institute of Applied Mathematical Research, Karelian Research
Centre, Russian Academy of Science
11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: zaika@krc.karelia.ru
tel.: (8142) 766312
