Труды Карельского научного центра РАН
№ 8. 2016. С. 24-33 DOI: 10.17076/mat392
УДК 519.6:539.2
КРАЕВАЯ ЗАДАЧА СО СВОБОДНОЙ ГРАНИЦЕЙ: ГИДРИРОВАНИЕ ЦИРКОНИЕВОГО СПЛАВА
Ю. В. Заика, Н. И. Родченкова
Институт прикладных математических исследований Карельского научного центра РАН
Одним из важных требований к изделиям из циркониевых сплавов активной зоны реакторов является низкое поглощение водорода, поскольку водородное охрупчивание может стать одной из причин разрушения циркониевой оболочки. В зависимости от уровня содержания водорода и температуры эксплуатации водород может находиться в циркониевых сплавах в виде твердого раствора или в виде гидридов. Наибольший охрупчивающий эффект оказывают гидриды, так как они могут служить участками образования и развития трещин. Проблема состоит в моделировании динамики свободной границы фазового перехода и оценке распределений концентраций в гидриде и в сплаве. В статье представлены математическая модель гидрирования циркониевого сплава с учетом фазового перехода (гидридообразования), итерационный вычислительный алгоритм решения нелинейной краевой задачи со свободной границей раздела фаз на основе неявных разностных схем и результаты вычислительных экспериментов.
Ключевые слова: гидрирование; нелинейные краевые задачи со свободной границей; разностные схемы; численное моделирование.
Yu. V. Zaika, N. I. Rodchenkova. BOUNDARY-VALUE PROBLEM WITH FREE BOUNDARY: ZIRCONIUM ALLOY HYDROGENATION
One of the most important requirements for the reactor's active zone materials (made of zirconium alloys) is low hydrogen absorptivity since hydrogen embrittlement may cause zirconium cladding damage. Depending on the hydrogen content and operation temperature, hydrogen may be present in zirconium alloys as a solid solution or as hydrides. Hydrides have the greatest embrittlement effect on alloys as they can form and enlarge cracks. The problem is to model the dynamics of the moving boundary of phase transition and to estimate the concentration distribution in hydride and in solution. This paper presents a mathematical model of zirconium alloy hydrogenation taking into account phase transition (hydride formation), the iterative computational algorithm for solving the nonlinear boundary-value problem with free phase boundary based on implicit difference schemes, and the results of computational experiments. This work is supported by the Russian Foundation for Basic Research (Project No. 15-01-00744).
Keywords: hydrogenation; nonlinear boundary-value problems with free phase boundary; difference schemes; numerical simulation.
Введение
Интерес к взаимодействию изотопов водорода с конструкционными материалами носит многоплановый характер [1—9]. Энтузиасты говорят не только об энергетике, но и о водородной экономике [6]. Некоторые математические модели дегидрирования и водородо-проницаемости, связанные с тематикой данной статьи, исследованы в [10-14]. Одним из важных требований к изделиям из циркониевых сплавов активной зоны реакторов является низкое поглощение водорода, поскольку водородное охрупчивание может стать одной из причин разрушения циркониевой оболочки. В зависимости от уровня содержания водорода и температуры эксплуатации водород может находиться в циркониевых сплавах в виде твердого раствора или в виде гидридов. Наибольший охрупчивающий эффект оказывают гидриды, так как они могут служить участками образования и развития трещин.
При разработке математической модели гидрирования авторы следовали работе [15]. Проблема состоит в моделировании динамики свободной границы фазового перехода и оценке распределений концентраций в гидриде и в сплаве. В статье представлены математическая модель гидрирования пластины из циркониевого сплава с учетом фазового перехода (гидридообразования) и итерационный вычислительный алгоритм решения нелинейной краевой задачи со свободной границей раздела фаз на основе неявных разностных схем.
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ГИДРИРОВАНИЯ
Вначале кратко опишем условия эксперимента (подробнее см. [15]). Пластина из сплава Zr — 1МЪ шлифуется с одной стороны, другая сторона практически водородонепроница-ема, торцами пренебрегаем; температура образца Т и давление газообразного водорода р поддерживаются постоянными (предпринимаются специальные меры охлаждения).
Выделим тонкий объемный слой, в котором при относительно большом давлении напуска (р ~ 2 атм) распределение Н можно считать равномерным. Только с некоторой начальной глубины ¿о начинает ощущаться диффузионное сопротивление. Когда концентрация растворенного атомарного водорода достигает определенного предела, этот слой образует начальную корку гидрида. Дальнейший перенос водорода в образец уже осуществляется сквозь растущий слой гидрида.
Обозначим: Ь — толщина пластины; ¿о — толщина слоя, в который водород абсорбиру-
ется относительно легко и еще не ощущается диффузионное сопротивление (будущая начальная корка гидрида); u(t) — концентрация H в 4)-слое (1н/см3); Q — концентрация, по достижении которой решетка перестраивается и раствор преобразуется в гидрид; c(t, x) — концентрация H в (L — .£0)-слое; ц — газокинетическая константа. Температура пластины и давление напуска водорода постоянны (Т = const, p = const). Схематически обозначения изображены на рисунке 1.
H2; ф
-фаза гидрид Zr-1Nb
Q «-фаза
раствор
u(t)
0 0 0 \c(t,x)
0 H 0 ° 0-----
Рис. 1. Этап I (начальное насыщение)
Согласно кинетической теории газов плотность Jp падающего на поверхность потока частиц (в данном случае молекул H2) связана с давлением p по формуле Герца-Кнудсена: Jp = pN2nmkT (k — постоянная Больцмана, m — масса молекулы H2). В контексте эксперимента удобно в качестве единиц измерения выбрать [x, I, L] = см, [p] = Торр. Тогда численно получаем зависимость Jp = цp, ц(Т) и 2, 474 ■ 1022/VT ([ц] = 1н2/(Торрсм2 с), [Т] = К, под корнем численное значение). Поскольку диффундирует атомарный водород, то для единообразия подсчет будем вести в атомах H: Jp = 2до. Но только малая часть H окажется в абсорбированном состоянии: Jabs = 2цsp (s ^ 1). Множитель s имеет смысл доли налетающих H, которые оказались в итоге в приповерхностном объеме. Объединяем более элементарные стадии физадсорбции, диссоциации и растворения в одну: s — эффективный коэффициент абсорбции.
Этап I: растворение H в Zr — 1Nb
Для диффузионного слоя толщины (L — i0) имеем стандартную краевую задачу:
дс д2с
« = DdX2. x € «0.L). t> (1)
c(0,x) = 0, x € [£o,L],
c(t,4)= u(t), dxc(t, L) = 0. (2)
0
L
0
Граничное условие c(t, to) = u(t) отражает непрерывность распределения H в сплаве, а dxc|L = 0 — непроницаемость стороны пластины x = L. Здесь и в дальнейшем считаем, что коэффициенты подчиняются закону Аррениуса по температуре, в частности D = D0 exp{ — Ed /[RT]}. В течение одного эксперимента поддерживается T = const.
Для концентрации u(t) запишем ОДУ, исходя из баланса потоков:
U(t)to = 2^sp — bu2 + D9xc|^, u(0) = 0. (3)
Содержательный смысл: за одну секунду через 1 см2 абсорбировалось за счет давления 2^sp атомов H, но есть встречный поток десорбции bu2 (b — эффективный коэффициент рекомбинации) и диффузионный отток. Рассогласование плотностей этих потоков идет на накопление атомов водорода в to-слое (uto). Уравнение (3) нужно рассматривать совместно с (1)-(2), поскольку u(t) определяет граничную концентрацию в (2). При небольшом p (без образования гидрида) в равновесии (когда производные равны нулю) имеем
2^sp — bu2 = 0 ^ u = Г/p, Г = л/2^s/b.
Следовательно, динамика (3) в статике согласуется с законом Сивертса u a ^p, Г — коэффициент растворимости. Подчеркнем, что речь о растворенном атомарном диффузионно подвижном водороде. В эксперименте «насыщение-дегазация» учитывается общее поглощение водорода, включая обратимый захват и гидридные фазы — коэффициент Г может иметь другой смысл и численное значение. Технически нетрудно учесть обратимый захват H в (L — to)-слое дефектами материала, но в рассматриваемой задаче считаем ловушки второстепенным фактором.
Этап II: гидридообразование и движение границы фазового перехода
По достижении u(t) порогового уровня Q = Q(T) происходит образование гидрида. Считаем, что этот переходный процесс относительно быстрый. Итак, u(t) = Q ^ to-корка превратилась в гидрид. Начинаем новый отсчет времени (t = 0) и моделируем рост гидридной фазы. К этому моменту c(0, x) = <^(x) (распределение растворенного водорода с предыдущего этапа), ^>(t0) = Q, to-слой уже гидридный и сквозь него диффундирует растворенный H. Обозначим концентрацию H в гидриде через v(t,x). Общая концентрация равна Q + v(t, x). Схематически обозначения изображены на рисунке 2.
MP
-фаза гидрид Zr-1Nb
\v(t,x) о \ 0 (X-фаза раствор
7 7/Q/ // / \ c(t,x) о N,. о о ^--- О 0 о °
/ /
/ /
0 1о КО I-
Рис. 2. Этап II (движение фазового перехода)
В пластине с растущей коркой гидрида (ж = 4(£) — граница раздела фаз, 4(0) = 40) запишем диффузионные уравнения:
| = Лg, x e (0.t(t)).
v(0, x) = 0, x e [0,to],
I = DS, x e (t(t),L),
(4)
(5)
c(0,x) = <^(x), x e [to,L].
Здесь Б* — коэффициент диффузии Н в новом материале (гидриде). Граничные условия на «входе-выходе» запишем аналогично I:
dv
2^s*p — b*v2(t, 0) = — D*—
x=o
dc dx
0.
(6)
В гидриде градиент концентрации возникает практически сразу, поскольку диффундировать сквозь ¿-фазу значительно труднее. Нет «накопительного» слоя (аналога 4о-слоя на этапе I). Граница раздела фаз ж = 4(£) уже подвижна, и для «склейки» диффузионных уравнений требуется два условия на стыке ж = 4(*).
Начнем с уравнения типа Стефана, описывающего динамику движения свободной границы раздела фаз:
[v(t,t(t)) + Q — c(t,t(t))] t(t) =
dv
= —D*—
dx
+ D^
£(t) dx
m
Появляется разрыв концентраций:
г(М) + д > Q > с(М), 4 = 4(£), í > 0.
При этом концентрация г(£,4) должна быть пренебрежимо малой (нет существенного «сопротивления» со стороны £Т-сплава). Принимаем г(£,4(£)) = 0. Поступающий поток из
26
H
2
L
фазы практически полностью уходит на формирование нового слоя гидрида (сдвиг границы х = ¿(¿)) и в раствор.
Таким образом, принимаем следующие условия на свободной границе раздела фаз:
и(М(*)) = 0,
дг>
[Я — с(М(г))] ¿(*) = ——
+ ^ е. дх
(7)
Вначале [... ] = 0 (с(0^о) = ^(¿0) = Я), потом (£ > 0) появляется разрыв концентраций ([... ] > 0) и ¿(¿) > 0. Фронт движется с концентрацией Я (общей) в сторону х = Ь. В пределе (£ ^ имеем ¿(¿) ^ Ь. Скорость движения фронта замедляется. Уровень с(£, х) выравнивается и медленно растет до значения Я в пределе, когда слоем уже можно практически пренебречь.
Вычислительный алгоритм
Этап I: растворение Н в —
Следуя технике разностных схем, введем сетку {хт = ¿0 + тЛ,х, т = 0,1,..., М} (Л,х = (Ь—¿о)/М) по пространственной переменной и сетку по времени {¿га = пЛ^, п = 0, 1,...}. Обозначим через {Ст}, {ип} приближенные значения концентраций в (Ь — ¿0)-слое (с(£га,хт)) и в ¿0-слое (и(£га)) соответственно. Рассмотрим неявную схему для уравнения диффузии (1) и неявный метод Эйлера для ОДУ (3)
Для нахождения начальных коэффициентов а1, воспользуемся следующими соображениями. Подсчитаем предварительно значения СП+1 по явной разностной схеме (в равенстве (8) справа заменяем п + 1 на п). На (п + 1)-м слое по времени аппроксимируем производную дхф0 ^ [—3СП+1 + 4СП+1 — СП+1]/2Л,ж и подставим в ОДУ (3). Принимая во внимание первое из граничных условий (2), имеем ига+1 = СП+1. В итоге получаем выражение СП+1 = /0(СП+1,С2П+1) (положительный корень квадратного уравнения):
Ь(С0га+1)2 + А1СП+1 + А 2 = 0,
А = ¿0 Ь--1 +3Д[2^х]-1,
А2 = —2^р +
+ Д[СП+1 — 4Сп+1][2^х]-1 — ¿0^"^
"V.
Зная численное значение СП+1 и выражение СП+1 = а1С"+1 + в1, получаем а1 = 0, в1 = СП+1. По «1, в1 вычисляем оставшиеся прого-ночные коэффициенты ат, вт, т = 2,..., М.
Ближайшая цель — найти значение СМ+1, необходимое для реализации метода прогонки. Рассмотрим второе краевое условие непроницаемости при х = Ь в (2). Используем разностную аппроксимацию
джс|ж=£ ~ [С^+Л — 4С^+11
с ™+1 _С™
ст_ст
и находим выражение См
/ТП+1 _
С=
сп+1 _2С"+1 + С"+1
7-, ст-1 2Ст + Ст+1 /0\
п-~-, (8) =[(4 — ам-1)вм — вм-1][(ам-1 — 4)ам + 3]
1
и™+1 — и™
¿0 и-= 2^р — Ь(ига+1)2+
%
+ д —3С0 +4С1 — с2 с™+1 = и™+1
Рассмотрим уравнения перехода с п-го слоя на
(п + 1)-й слой по времени (п ^ 0, 0 < т < М): ад =
ст-1 — [-+2]ст+1+ст+\+-ст = 0.
Значения в начальный момент времени известны: ст = и0 = 0 (0 < т < М). Следуя методу прогонки (алгоритм Томаса), ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое по времени в виде
ст+1 = ат+1Ст+\ + вт+1, т = 0,..., М — 1.
Прогоночные коэффициенты следующие (т = 1, 2,..., М — 1): ат+1 = [2 + ш — ат]-1,
вт+1 =
вт+-ст
2 + ш —
(9)
Следующий этап: с текущими приближениями значений СП+1, С^1 решаем обратным ходом прогонки трехдиагональную систему линейных алгебраических уравнений и находим новые приближения концентраций СП+1 (и остальные значения Ст+1 для т = 3,..., М — 1). Далее снова пользуемся формулой С"+1 = /0(Сп+1,Сп+1). После этого корректируем значения прогоночных коэффициентов (9), определяем С^1 = [(4 — ам-1)вм — вм-1][(«м-1 — 4)ам+3]-1 и повторяем вычисления (возвращаясь к началу абзаца) до установления граничных значений С^м (обычно 2-3 итерации). Критерием окончания вычислений на этапе I является выполнение условия и"+1 ^ Я = 0, 98 и (в программной реализации при текущих значениях параметров).
Результат окончания этапа I: образовался гидридный слой толщины ¿0 с концентрацией водорода и = Я, а в слое раствора толщины (Ь — ¿0) имеется распределение концентрации растворенного водорода <^(х).
Этап II: гидридообразование и движение границы фазового перехода
Авторы остановились на следующем варианте алгоритма: движущиеся отрезки [0,4(£)] и [4(£),Ь] преобразуются в [0,1] с последующим выбором равномерных сеток. Соответствующая замена переменных с новым отсчетом времени (¿о = 0): (¿,ж) о (¿,у), ж = 4(£)у, г(£,у) = г(£,ж(£,у)); (¿,ж) о (¿,г), ж = 4(£) + -4(£)], -(¿, г) = с(£, ж(£, г)).
Краевая задача (4)-(7) гидрирования циркониевого сплава после замены переменных:
dv D* d2v 4(t)y dv
(10)
dt 42(t) dy2 ¿(t) dy' y e (0,1), v(0,y) =0 (0 < y < 1), v(t, 1) = 0
- b*v2(t, 0) = -D ■
¿(t) dy
y=o
dc
D
д2c ¿(t)(1 - z) dc
(11) (12)
dt [L - 4(t)]2 dz2 L - 4(t) dz'
dc
c(0,z) = ^(z) (0 < z < 1), — =0, (13)
dz z=i
D dv
+
D
д c
L - ¿(t) dz
z=0
(14)
Начальное распределение ^(z) определяется первым этапом. Формально можно забыть о физическом смысле функции 4(t) как границы раздела фаз и рассматривать ее как функциональный параметр. Если формально считать распределения v(t,y), c(t, z) известными, то из условия Стефана (14) получаем нелинейное функционально-дифференциальное уравнение 4(t) = G(t, .£(•)). По решению 4(t) определяются коэффициенты модели. Поэтому вычислительный алгоритм основан на неявных разностных схемах и носит итерационный характер. Итерации будут связаны с уточнением значения 4(t) на каждом временном слое.
Введем следующие сетки: {ym = mhy, m = 0,1,...,M}, hy = 1/M - шаг по пространственной переменной y; {z^ = khz, k = 0,1,..., K}, hz = 1/K - шаг по пространственной переменной z; {tn = nht, n = 0, 1,...}, ht — шаг по времени t. После введения соответствующих нормировок оставим прежние обозначения для функций v := v/v, c := c/u, v = *p/b*, u = A/2^sp/b. Обозначим че-
рез {Vm}, {cCfc} приближенные значения концентраций в гидридном 4(^-слое (v(tn,ym)) и в (L - 4^))-слое раствора (c(tn, zk)).
1. Начальный этап: от £ = 0 к ¿1 = Д£
Будем рассматривать последовательно два слоя: сначала гидридный, потом слой раствора 1ЖЬ (это соответствует и последовательности движения растворенного водорода сквозь материал с течением времени).
1.1 Гидридный слой
Полагаем 40 = 4(0) = 0. Параметр 40 подлежит итерационному уточнению, значение 4о известно. Полагаем в краевой задаче (10)-(11)
4 = 0, 4 = 4о:
дг Б * д2г
ж = 12* ■ ду-, у е (0,1), (15)
г(0, у) = 0, у е [0,1], -(¿, 1) = 0, Б дг
2^ *р - Ь* г2(£, 0) = - —* ■ — . (16)
4о ду у=о
Для уравнения диффузии (15) в безразмерной форме (- := -/г) рассмотрим неявную разностную схему на слое 1: т = 1,..., М — 1,
c1 - co
r m r m
D vi_i - 2vm, + vm
m+1
ht
¿0
hy
(17)
Рассмотрим уравнения перехода с нулевого на первый слой по £ (п ^ 0, 0 < т < М):
ут- - к+2]т-т+кт+1+^ут=0,
■ш2 = /[Л^ Б *]. Начальные данные известны: V= 0 (0 ^ т ^ М). Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки на первом слое в виде
V1 = + вт +1, т = 0, . . . , М - 1.
Прогоночные коэффициенты (т = 1,..., М-1):
а m+1 =
ßm +1 =
1
2 + w2 - а
ßm + W2 V
2 + w2 - а
(18)
Для нахождения начальных коэффициентов а1, в1 воспользуемся следующими соображениями. Подсчитаем предварительно значения —по явной разностной схеме (в равенстве (17) справа заменяем 1 на 0). На первом слое по £ аппроксимируем производную ду-|у=о ~ [—ЭУо1 + 4—1 - —1]/2Л,у и подставляем в граничное условие (16):
О Р и - Ь2 Б * -3—1 +4—1 - -1 2^8 * - - Ь * г (Уд)2 = —----о 1 2
^0
2hy
28
В итоге получаем Уо1 = /1 (Т/^1, Т/21) (положительный корень квадратного уравнения):
М(Уо)2 + В1/1 + В2 = 0, В1 = 3Б,[2Лу ¿0]-1, В = ^[У/ — 4У/][2Лу ¿0]-1 — ргГ1.
Зная численное значение и выражение для
Т01 = а^У/ + в1, получаем а1 = 0, в1 = У)1. По а1, в1 вычисляем оставшиеся коэффициенты ат, вт, т = 2,..., М по формулам (18).
Ближайшая цель — найти значение Ум, необходимое для реализации прогонки. Из краевого условия V (¿, 1) = 0 получаем = 0.
Следующий этап: с текущими приближениями значений Том решаем обратным ходом прогонки трехдиагональную систему линейных алгебраических уравнений и находим новые приближения концентраций У^
(и остальные значения , т = 3,..., М — 1).
После этого возвращаемся к уточнению То1 = /1(111,1/21) в силу (16), корректируем прогоночные коэффициенты (18) и повторяем вычисления (возвращаясь к предыдущему абзацу) до установления распределения ■й(£1,у) [Ут, т = 0,..., М]. Шаг по времени мал, так что достаточно нескольких итераций.
1.2 Слой раствора —
Переходим к краевой задаче (12)-(13), возвращаясь к предыдущему слою (£ = 0):
дс
Ж
Б
д2с
[Ь — ¿о]2 дг2'
(19)
дс
с(0,г) = #г), г е [0,1], — =0. (20)
дг ^=1
Недостающее граничное условие определяем из уравнения Стефана (14):
¿о = 0, V = 0 ^ дгс|г=о = 0.
Для безразмерного уравнения диффузии (19) (с := с/и) рассмотрим неявную схему на первом слое по времени: к = 1,..., К — 1,
У1 У о Л*
Б
С1-1 — 2С1 + С1+1
[Ь — ¿о]2
Л2
(21)
Уравнения перехода на первый слой
(п ^ 0, 0 < к < К): Ш4 = [Ь — ¿о]2 Л2/[Л* Б],
(51-1 — [-4 + 2]С1 + ^¿+1 + -4 С£ = 0.
Начальные данные: СЦ = ) (0 ^ к ^ К). Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки на
первом слое по £ в виде С1 = 0:1+1 С1+1 + в1+1, к = 0,..., К — 1. Прогоночные коэффициенты (к = 1,..., К — 1):
«1+1
в1+1 =
1
2 + Ш4 — «1
в1 + -4 С
2 + Ш4 — «1'
(22)
Для нахождения начальных коэффициентов а1, в1 подсчитаем предварительно значения ¿12 по явной разностной схеме (в равенстве (21) справа заменяем 1 на 0). На первом слое аппроксимируем производную дг с|^=0 — [—3С0 + 4С1 — С ]/2Лг и приравниваем к нулю в силу дгс|^=о = 0 (это условие выполняется только при переходе с нулевого на первый слой по времени, в дальнейшем используется условие Стефана на подвижной границе раздела фаз). В итоге получаем Со = [4С — С2]/3.
Зная численное значение Со и выражение
С = а1 С + в1, получаем а1 = 0, в1 = Со. По «1, в1 вычисляем оставшиеся коэффициенты а1, в1, к = 2,..., К по формулам (22).
Теперь найдем значение СК, необходимое для реализации прогонки. Рассмотрим краевое условие непроницаемости (20). На первом слое по £ аппроксимируем производную д^с|г=1 - [СК-2 — 4СК-1 + ]/2Л* = 0 и подставим вместо СК_2 К_ 1 выражение = а1+1С1+1 + в1+1 при соответствующих к. В итоге получаем СК = [(4 — «к-1)вк — вк-1][(ак-1 — 4)ак + 3]-1. Следующий этап: с текущими приближениями значений С^к решаем обратным ходом прогонки трехдиаго-нальную систему линейных уравнений и находим новые приближения концентраций С^ 2
(и остальные значения С1, к = 3,..., К — 1).
После этого возвращаемся к уточнению Со = [4С1 — С ]/3, корректируем значения прогоночных коэффициентов (22) и повторяем вычисления (возвращаясь к предыдущему абзацу) до установления распределения с(£1,г), г е [0,1] [<С'1, к = 0,..., К]. Шаг по времени мал, так что достаточно нескольких итераций. Отметим, что при этом на начальном этапе понизится входная концентрация ¿(¿1, 0) < Я.
1.3 Итерационное уточнение ¿о
После этих вычислений в уравнении Стефана (14) получаем [Я — с(£1,0)] > 0 и положительную правую часть. Находим новое приближение ¿о > 0 и возвращаемся в начало излагаемого алгоритма (пункт 1.1). Заме-
тим, что теперь в пунктах 1.1, 1.2 необходимо учитывать слагаемые в правых частях диффузионных уравнений (10), (12), содержащие 4о. Формально нужно перейти к пунктам 2.1-2.3 (при п = 0) и продолжить вычисления до установления 4о. Через несколько итераций получаем установившееся значение 4о и определяем 41 = 4о + 4оЛ,£. Переходим к слою ¿1 = Д£, приняв в качестве начального приближения 4 = 4о. Далее уточняем значение 4 и распределения г(£2,у), -(¿2,г) по изложенной схеме.
2. Переход от п-го к (п + 1)-му слою по £
Опишем алгоритм нахождения распределений на (п + 1)-м слое -(£п+1, у) [—т+1, т = 0,..., М], -(¿п+1, г) [¿-"+1, к = 0,..., К], зная распределения с предыдущего слоя по времени г(£п,у), -(¿п,г). Рассмотрим последовательно два слоя материала: сначала гидридный, затем слой раствора - 1ЖЬ.
2.1 Гидридный слой
Начальное значение скорости движения свободной границы раздела фаз 4(£п) = 4(£п-1) берется с предыдущего временного слоя. Параметр 4(£п) подлежит итерационному уточнению, значение 4(£п) известно.
Для уравнения диффузии (10) в безразмерной форме рассмотрим неявную разностную схему на (п + 1)-м слое: т = 1,..., М - 1,
T/n+1 _ T/n
г m_v_m
ht
D*
T/n+1 - 2T/n+1 + T/n+1
hy
+
¿2(tn)
^^(¿ra) ym T/m+1 — T/m+1 ¿(tn) ^ 2hy '
+
(23)
W2(tn) —
42(tn)hy 4(i„)4(i„)ymhy
w3(tnj m) —
ht D*
2D*
am +1 —
Am +1 —
1 + W3
2 + W2 + (w3 - 1}am'
(1 - -Шз)вт + W2 Yn 2 + W2 + (W3 - 1)«m '
(24)
Для нахождения а1, в1 подсчитаем предварительно значения У""+1 по явной разностной схеме (в равенстве (23) справа заменяем п + 1 на п). На (п + 1)-м слое по £ аппроксимируем
ду-|у=о - [-3—П+1 + 4Т-1"+1 - -га+1]/2^у и подставим в граничное условие (11):
- b*V(/0n+1)2 —
¿(tn)
-3V01+1 + 4V1n+1 - /2n+1
2hy
В итоге получаем -"+1 = /^—"Л —1+1) (положительный корень квадратного уравнения):
Ь* г;(—о)2 + £1—о + £2 =0, £1 = 3Б*[2^у4]-1,
1
Рассмотрим уравнения перехода с n-го слоя по t на (n + 1)-й слой (n ^ 0, 0 < m < M):
[1 - W3]Tm+1 - [W2 + 2]С+1 +
+ [1+w3]Vm+1+W2C — 0,
Ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое в виде У"+1 = У—"+1 + вт+1, т = 0,..., М - 1. Прого-ночные коэффициенты (т = 1, 2,..., М - 1):
£2 = б * У - 4Щ2Лу 4] - 2^ * рг
V = У"+1. Зная значение -"+1 и выражение —о"+1 = а1—"'+1 + в1, получаем а1 = 0, в1 = УП+1. По а1, в1 вычисляем коэффициенты ат, вт, т = 2,..., М по формулам (24).
Теперь найдем значение —М+1, необходимое для реализации прогонки. Из краевого условия г(£, 1) = 0 получаем УП+1 = 0. Следующий этап: с текущими приближениями —"М1 решаем обратным ходом прогонки трехдиаго-нальную систему линейных уравнений и находим новые приближения концентраций 1—'г+1
(и остальные У"+1, т = 3,..., М - 1).
После этого возвращаемся к уточнению -"+1 = /1(—"+1, —"+1) в силу уравнения (11), корректируем значения прогоночных коэффициентов (24) и повторяем вычисления до установления распределения г(£га+1,у) [—.+1, т = 0,..., М]. Шаг по времени мал, так что достаточно нескольких итераций.
2.2 Слой раствора - 1ЖЬ
Переходим к краевой задаче (12)-(14), возвращаясь к предыдущему слою (£ = £п). Недостающее граничное условие определяем из уравнения Стефана (14).
Для уравнения диффузии (12) в безразмерной форме рассмотрим неявную разностную схему на (п + 1)-м слое: к = 1,..., К - 1,
C/n+1 _ c/n
ht
D
C<n+1 _ 2Cn+1 + C<n+1
k—1 k
[L -
h2
+
+
¿[1 - Zk] /+1 - cj—1
L-4
2hz
(25)
30
Рассмотрим уравнения перехода с п-го слоя по £ на (п + 1)-й слой (п ^ 0, 0 < к < К):
[1 — -НСП+11 — [-4 + 2]СП+1 +
Найдем значение СК+1, необходимое для прогонки. Рассмотрим краевое условие непроницаемости (13). На (п + 1)-м слое по £ аппроксимируем производную дгс|г=1 — —
1+1
? (, [ь — ¿(¿™)]2 -4(£п) =--,
^ _ [Ь — ¿(¿™)] ¿(¿™) [1 — ¿1 ] Нг
-5 (£п) к) = .
Ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое в виде СП+1 = а^С^1 + в1+1, к = 0,..., К — 1. Прогоноч-ные коэффициенты (к = 1, 2, . . . , К — 1):
«1+1 =
в1+1 =
1 + -5
2 + --4 + (-5 — 1)«1'
Я с п+1
ь и о
+
¿(¿п) и
Б
Ь — ¿(¿™)
+ [1 + ^С^! + ^С" = 0, 4СК11! + 3<СК+1]/2Л^ =0 и подставляем вместо СК+2 к- 1 соответствующие соотношения СП+1 = а^СП+Чв^ 1. В итоге имеем СК+1 =
[(4 — «к-1)вл - вк-1][(ак-1 — 4)«к + 3]-1. Следующий этап: с текущими приближениями СП+1 решаем обратным ходом прогонки трехдиагональную систему уравнений и находим новые приближения концентраций СП+1
(и остальные значения СП+1, к = 3,..., К — 1). После этого возвращаемся к уточнению
СП+1 = /2(^Г+1,Сс2П+1) в силу (14), корректируем значения прогоночных коэффициентов (26) и повторяем вычисления до установления с(£™+1, ¿) [<С'п+1, к = 0,..., К].
2.3 Итерационное уточнение ¿
После указанных вычислений из уравнения Стефана (14) находим новое приближение ¿(¿™) > 0 и возвращаемся в начало излагаемого алгоритма (этап II, пункт 2.1). Через несколько итераций получаем установившееся значение ¿(¿п) и определяем ¿(¿п+1) = ¿(¿п) + ¿(¿п)Л;. Переходим к слою ¿п+1, приняв в качестве начального приближения ¿(£п+1) = ¿(¿п). Далее уточняем значение ¿(£п+1) и распределения V(¿п+2, у), с(£п+2, ¿) по изложенной схеме.
Тестирование алгоритма
Использовались следующие входные данные: Ь = 6 х 10-2 см, ¿0 = 1,3 х 10-3, р = 1520 Торр, Т = 593 К, Я = 6,4 х
(1 — -5)в1 + -4 С"
2 + --4 + (-5 — 1)«1'
(26)
Для нахождения а1 , в1 подсчитаем предварительно значения С"'+1 по явной разностной схеме (в равенстве (25) справа заменяем п + 1 на п). На (п+1)-м слое аппроксимируем производные дус|у=1
д^с|^=о — [—3^П+1 + 4СП+1 — С-^1]^* и подставим в граничное условие (14):
¿(¿п) =
Б, V У^ — 4УЙЙ +3УП+1
2Лу
+
—3Сп+1 + 4С"+1 — С"+1 2Л '
10-20 см-3, Бо = 2,2 х 10
-3 см2 с-1, =
В итоге получаем С"+1 = /2(С"+1, С2™+1) (Ум+1, ^ = 0,1, 2 подсчитаны в пункте 2.1):
С п+1 _
Я
-6 Я + Сп+1 — 4<сп+1 + и 2 1
+ -7 [Уп- — 4Уп+11 + 3У£+1]] [-6 — 3]-1
-б(£п) = ^[Ь — ¿^¿^,
-70п) = Б Ь — ¿('га) ' ^
Б
¿(¿п)
и
Зная значение СП+1 и выражение СП+1 = а1С"+1 + в1, получаем а1 =0, в1 = СП+1. По а1 , в1 вычисляем оставшиеся коэффициенты а1, в1, к = 2,..., К по формулам (26).
35 х 103 Дж/моль, Б,0 = 1, 5 х 10-3 см2 с-1, = 59 ■ 103 Дж/моль, Ь = 5 х 10-23 см4 с-1, Ь, = 3 х 10-24 см4 с-1, 8 = 7 х 10-6, = 6 х 10-7. На рисунке 3 представлена возможность анализировать динамику распределения концентраций растворенного атомарного водорода в слоях гидрида и сплава.
Заключение
Одним из важных требований к изделиям из циркониевых сплавов активной зоны реакторов является низкое поглощение водорода, поскольку водородное охрупчивание может стать одной из причин разрушения циркониевой оболочки. Наибольший охрупчиваю-щий эффект оказывают гидриды, так как они могут служить участками образования и развития трещин.
x 102
v(t,x)
7 6 ■
5 ■ 432 ■ 1 ■
---t =15 с
-------t = 60 с
12
10-
8 ■
6 ■
4 ■
2 ■
x 101
c(t,x)
x 10-7
Л 1.2 \ 1 fT f
\ 0.8
\ 0.6
\ 0.4
\ 0.2 t [с]
\ 0 \ \\ 2 4 6 8 1 0
• \\
>v \ \ ч L = 0.06 \ \ '— |(t) -------------
0.01 0.02 0.03 0.04
0.0005 0.001 1о 0.0015
Рис. 3. Распределения концентраций растворенного водорода в слоях гидрида и сплава
8
При проведении экспериментальных исследований практически невозможно на основе доступных внешних измерений определить распределение водорода в материале в зависимости от времени, поскольку атомарный растворенный водород является чрезвычайно подвижной фазой внедрения.
В статье представлена модель гидрирования пластины из циркониевого сплава, учитывающая динамику абсорбционно-десорбционных процессов и движение свободной границы раздела фаз (гидрид-металл). На основе неявных разностных схем разработан итерационный вычислительный алгоритм решения нелинейной краевой задачи.
Работа выполнена при поддержке РФФИ (грант № 15-01-00744).
ЛИТЕРАТУРА
1. Водород в металлах / Ред. Г. Алефельд, И. Фелькль. М.: Мир, 1981. Т. 1, 506 е.; т. 2., 430 с.
2. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.
3. Писарев А. А., Цветков И. В., Марен-ков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
4. Черданцев Ю. П., Чернов И. П., Тюрин Ю. И. Методы исследования систем металл-водород. Томск: ТПУ, 2008. 286 е.
5. Изотопы водорода. Фундаментальные и прикладные исследования / Ред. А. А. Юхимчук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697 е.
6. The hydrogen economy / Eds M. Ball, M. Wietschel. Cambridge University Press, 2009. 646 p.
7. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
8. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids // Technical Physics. 1999. Vol. 44, iss. 1. P. 90-94. doi: 10.1134/1.1259257
9. Lototskyy M.V., Yartys V. A., Pollet B.G., Bowman R. C. Jr. Metal hydride hydrogen compressors: a review // International Journal of Hydrogen Energy. 2014. Vol. 39, iss. 11. P. 58185851. doi: 10.1016/j.ijhydene.2014.01.158
10. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen influence // Acta Mechanica. 2008. Vol. 195, iss. 1. P. 295-304. doi: 10.1007/s00707-007-0568-z
11. EvardE.A, Gabis I.E., Yartys V.A. Kinetics of hydrogen evolution from MgH2: Experimental studies, mechanism and modelling // International Journal of Hydrogen Energy. 2010. Vol. 35, iss. 17. P. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
12. 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. 2009. Vol. 33, iss. 10. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018
13. Zaika Yu. V., Bormatova E. P. Parametric identification of hydrogen permeability model by delay times and conjugate equations // International Journal of Hydrogen Energy. 2011. Vol. 36, iss. 1. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
®
14. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with dynamical conditions on surface // Mathematical Modelling. Nova Sci. Publishers. 2013. P. 269-302.
15. Pushilina N. S., Kudiiarov V. N., Laptev R. S. et al. Microstructure changes in Zr-1Nb alloy
References
1. Vodorod v metallakh [Hydrogen in metals]
/ Eds G. Alefel'd, J. Fel'kl'. Moscow: Mir, 1981. Vol.1 506 p.; vol. 2,403 p.
2. Vzaimodeistvie vodoroda s metallami [Interactions of hydrogen with metals]. Ed. A. P. Zakharov. Moscow: Nauka, 1987. 296 p.
3. Pisarev A. A., Tsvetkov I. V., Marenkov E. D., Yarko S. S. Pronitsaemost' vodoroda cherez metally [Hydrogen permeability through metals]. Moscow: MIFI, 2008. 144 p.
4. Cherdantsev Yu. P., Chernov I. P., Tyurin Yu. I. Metody issledovaniya sistem metall-vodorod [Methods of studying metal-hydrogen systems]. Tomsk: TPU, 2008. 286 p.
5. Izotopy vodoroda. Fundamental'nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies]. Ed. A. A. Yukhimchuk. Sarov: RFYaTs-VNIIEF, 2009. 696 p.
6. The hydrogen economy. Eds M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
7. Handbook of hydrogen storage: new materials for future energy storage. Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
8. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids. Technical Physics. 1999. Vol. 44, iss. 1. P. 90-94. doi: 10.1134/1.1259257
9. Lototskyy M.V., Yartys V.A., Pollet B.G., Bowman R. C. Jr. Metal hydride hydrogen compressors: a review. International Journal of
СВЕДЕНИЯ ОБ АВТОРАХ:
Заика Юрий Васильевич
рук. лаб. моделирования природно-технических систем, д. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: zaika@krc.karelia.ru тел.: (8142) 766312
Родченкова Наталья Ивановна
старший научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: nirodchenkova@krc.karelia.ru тел.: (8142) 766312
after pulsed electron beam surface modification and hydrogenation // Surface and Coatings Technology. 2015. Vol. 284. P. 63-68. doi: 10.1016/j.surfcoat.2015.07.082
Поступила в редакцию 25.05.2016
Hydrogen Energy. 2014. Vol. 39, iss. 11. P. 58185851. doi: 10.1016/j.ijhydene.2014.01.158
10. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen influence. Acta Mechanica. 2008. Vol. 195, iss. 1. P. 295-304. doi: 10.1007/s00707-007-0568-z
11. EvardE.A, Gabis I.E., Yartys V.A. Kinetics of hydrogen evolution from MgH2: Experimental studies, mechanism and modelling. International Journal of Hydrogen Energy. 2010. Vol. 35, iss. 17. P. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
12. 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. 2009. Vol. 33, iss. 10. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018
13. Zaika Yu. V., Bormatova E. P. Parametric identification of hydrogen permeability model by delay times and conjugate equations. International Journal of Hydrogen Energy. 2011. Vol. 36, iss. 1. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
14. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with dynamical conditions on surface. Mathematical Modelling. Nova Sci. Publishers. 2013. P. 269-302.
15. Pushilina N. S., Kudiiarov V. N., Laptev R. S., Lider A. M., Teresov A. D. Microstructure changes in Zr-1Nb alloy after pulsed electron beam surface modification and hydrogenation. Surface and Coatings Technology. 2015. Vol. 284. P. 63-68. doi: 10.1016/j.surfcoat.2015.07.082
Received May 25, 2016
CONTRIBUTORS:
Zaika, Yury
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: zaika@krc.karelia.ru tel.: (8142) 766312
Rodchenkova, Natalia
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: nirodchenkova@krc.karelia.ru tel.: (8142) 766312