CC BY
13
Труды Карельского научного центра РАН
№ 7. 2019. С. 42-52 DOI: 10.17076/mat1094
УДК 519.6:539.2
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ
ДВИЖЕНИЯ ГРАНИЦЫ ГИДРИДООБРАЗОВАНИЯ
В ОБОЛОЧКЕ ИЗ ЦИРКОНИЕВОГО СПЛАВА
Н. И. Родченкова1, К. В. Грудова2
1 Институт прикладных математических исследований КарНЦ РАН, ФИЦ «Карельский научный центр РАН», Петрозаводск, Россия
2 Отдел комплексных научных исследований КарНЦ РАН, ФИЦ «Карельский научный центр РАН», Петрозаводск, Россия
Одним из важных требований к изделиям из циркониевых сплавов активной зоны реакторов является низкое поглощение водорода, поскольку водородное охрупчивание может стать одной из причин разрушения циркониевой оболочки. В зависимости от уровня содержания водорода и температуры водород может находиться в циркониевых сплавах в виде твердого раствора или в виде гидридов. Наибольший охрупчивающий эффект оказывают гидриды, так как они могут служить участками образования и развития трещин. Проблема состоит в моделировании динамики свободной границы фазового перехода и оценке распределений концентраций в гидриде и в сплаве. В статье представлены математическая модель гидрирования циркониевого сплава с учетом фазового перехода (гидридообразования) и итерационный вычислительный алгоритм решения нелинейной краевой задачи со свободной границей раздела фаз на основе неявных разностных схем.
Ключевые слова: гидрирование; нелинейные краевые задачи со свободной границей; разностные схемы; численное моделирование.
N. I. Rodchenkova, K. V. Grudova. MODELING OF THE MOVEMENT OF THE FREE BOUNDARY OF HYDRIDE FORMATION IN A CYLINDRICAL SHELL MADE OF A ZIRCONIUM ALLOY
One of the most important requirements for the materials (zirconium alloys) used in the reactor active zone is low hydrogen absorptivity, since hydrogen-induced 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 initiate and enlarge cracks. The problem is to model the dynamics of the moving boundary of phase transition and to estimate the concentration distribution in the hydride and the solution. This paper presents a mathematical model of zirconium alloy hydrogenation taking into account the phase transition (hydride formation) and the iterative computational algorithm for solving the nonlinear boundary-value problem with free phase boundary based on implicit difference schemes.
Key words: hydrogenation; nonlinear boundary-value problems with free phase boundary; difference schemes; numerical simulation.
Введение
Общие проблемы водородной энергетики и водородного материаловедения обстоятельно изложены в монографиях [1, 2, 4-6, 8, 9, 12]. Некоторые модели дегидрирования и водоро-допроницаемости, связанные с тематикой данной статьи, исследованы в [11, 13-15]. Одним из важных требований к изделиям из циркониевых сплавов активной зоны реакторов является низкое поглощение водорода, поскольку водородное охрупчивание может стать одной из причин разрушения циркониевой оболочки. В зависимости от уровня содержания водорода и температуры водород может находиться в циркониевых сплавах в виде твердого раствора или в виде гидридов. Наибольший охрупчивающий эффект оказывают гидриды, так как они могут служить участками образования трещин.
При разработке математической модели гидрирования авторы следовали работам [7, 10]. Проблема состоит в моделировании динамики свободной границы фазового перехода и оценке распределений концентраций в гидриде и в сплаве. В статье представлены корректировка математической модели гидрирования, поставленной в [3], для полого цилиндрического образца циркониевого сплава с учетом фазового перехода (гидридообразования) и итерационный вычислительный алгоритм решения нелинейной краевой задачи со свободной границей раздела фаз на основе неявных разностных схем.
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ГИДРИРОВАНИЯ
Вначале кратко опишем условия эксперимента (подробнее см. [10]). Внешняя поверхность полой цилиндрической оболочки из сплава 2т — шлифуется, внутренняя поверхность практически водородонепроницае-ма, основаниями пренебрегаем; температура образца Т и давление газообразного водорода р поддерживаются постоянными (предпринимаются специальные меры охлаждения).
Выделим тонкий объемный слой у внешней поверхности, в котором при относительно большом давлении напуска (р ~ 2 атм) распределение Н можно считать равномерным. Только с некоторой начальной глубины (радиуса) ро начинает ощущаться диффузионное сопротивление. Когда концентрация растворенного атомарного водорода достигает определенного предела, начинается образование зародышей гидридной фазы, и этот слой отно-
сительно быстро преобразуется в начальную корку гидрида. Этот переходный процесс «от зародышей к корке» в приповерхностном объеме считаем практически мгновенным в масштабе времени последующего медленного движения фронта гидрирования на заметную глубину. Дальнейший перенос водорода в образец уже осуществляется сквозь растущий слой гидрида со значительно меньшей скоростью.
Обозначим: R - внешний радиус полого цилиндра; r0 - внутренний радиус; (R — р0) -толщина слоя, в который водород абсорбируется относительно легко и еще не ощущается диффузионное сопротивление (будущая начальная корка гидрида); u(t) - концентрация H в (R — ро)-слое (1н/см3); Q - концентрация, по достижении которой локально решетка перестраивается и возникают зародыши гидридной фазы; um = umax - принципиально возможная концентрация (пресыщенный раствор), Q < umax; Qh - концентрация атомов водорода в гидридной фазе (химически связанный водород, образующий гидрид как вещество); c(t,r) - концентрация растворенного H в (ро — г0)-слое; ц - газокинетическая константа. Температура цилиндра и давление напуска постоянны (Т = const, p = const).
Согласно кинетической теории газов плотность Jp падающего на поверхность потока частиц (в данном случае молекул H2) связана с давлением p по формуле Герца-Кнудсена:
Jp = p/yj 2nmkT
(k - постоянная Больцмана, m - масса молекулы H2). В контексте эксперимента удобно в качестве единиц измерения выбрать [r, р, R] = см, [p] = торр. Тогда численно получаем зависимость
Jp = fip, ц(Т) и 2.474 • 1022/vT
(М = 1н2/(торрсм2 с), [Т] = К, под корнем безразмерное численное значение). Поскольку диффундирует атомарный водород, то для единообразия подсчет будем вести в атомах H: Jp = 2ц'[). Только малая часть H окажется в абсорбированном состоянии: Jabs = 2^sp (s < 1). Множитель s имеет смысл доли налетающих H, которые оказались в приповерхностном объеме. Объединяем более элементарные стадии физадсорбции, диссоциации и растворения в одну: s - эффективный коэффициент абсорбции.
Этап I: растворение H в Zr-lNb
Для диффузионного слоя толщиной (ро — Го) имеем стандартную краевую задачу:
1 dc д2е 1 дс . . ч
d ' д = дГ2 + Г' дГ'Г е (Г0,р0) (1)
с(0, г) = 0, r е [го,ро],
с|ро = u(t), дгc\r0 =0. (2)
Граничное условие c(t, р0) = u(t) отражает непрерывность распределения H в сплаве, а дгc(t, Го) = 0 - непроницаемость внутренней поверхности цилиндра r = го. Здесь и в дальнейшем считаем, что коэффициенты подчиняются закону Аррениуса по температуре, в частности D = Do exp{ —ed/[RT]}. В течение одного эксперимента T = const.
Для концентрации u(t) запишем ОДУ, исходя из баланса потоков:
2^sp
1 — W
um -
+ DP
дг
1 — u-
um ?2
-bu2
2nRH
2проН = [nR2H — np20H] ^u,
Ро dt
2 h 2 , про д)c — bu + DR • дг
Ро
R2 — ро du
2R
dt
(3)
Содержательный смысл: за 1 с через 1 см2 абсорбировалось за счет давления 2^вр атомов Н (множитель [1 — ии^,1 ]2 - для торможения роста и(Ь) - вероятность диссоциативно абсорбировать уменьшается с ростом и), но есть встречный поток десорбции Ьи2 (Ь -эффективный коэффициент рекомбинации) и диффузионный отток. Рассогласование плотностей этих потоков идет на накопление атомов водорода в (Я—р0)-слое ([Я2—р^][2Я}-1 и). Уравнение (3) нужно рассматривать совместно с (1)-(2), поскольку и(Ь) определяет граничную концентрацию в (2). При небольшом р (без образования гидрида) в равновесии (когда производные равны нулю) имеем
2^вр [1 — иит1]2 — Ьи2 = 0 ^ и = Гур(1 +Тит1^Р)-1, г = V2ц.зЬ-1.
Следовательно, когда ит велико, динамика (3) в статике согласуется с законом Сиверт-са и а ур, Г - коэффициент растворимости. Подчеркнем, что речь — о растворенном атомарном диффузионно подвижном водороде. В эксперименте «насыщение-дегазация» учитывается общее поглощение водорода, включая
обратимый захват и гидридные фазы - коэффициент Г может иметь другой смысл и численное значение. Технически нетрудно учесть обратимый захват Н в (р0 — г0)-слое дефектами материала, но в рассматриваемой задаче считаем ловушки второстепенным фактором.
Этап II: гидридообразование и движение границы фазового перехода
По достижении концентрации и(Ь) порогового уровня Q = Q(T) происходит образование зародышей гидрида (¿-фазы). Дальнейшее поступление абсорбированного атомарного водорода расходуется на рост зародышей до образования сплошной корки гидрида (с концентрацией Qh > Q химически связанного водорода) и поддержание уровня с(Ь, р0) = Q в а-фазе раствора. При этом с входной стороны (внешней поверхности) зародыши растут и смыкаются быстрее, общая концентрация водорода уже начинает превышать уровень Qh за счет дополнительного растворения в гидриде. Когда приток к границе гидридного микрослоя г = ро начинает превышать отток в сплав Zr—1NЬ, появляется «движущая сила», и граница фазового перехода начинает смещаться, слой гидрида растет.
Считаем, что этот переходный процесс относительно быстрый. Разумеется, речь не о строгой последовательности описанных процессов, а о преимущественных стадиях. Детальное уточнение представляет самостоятельный интерес. Но для задачи моделирования динамики границы фазового перехода описание начальных данных имеет скорее смысл технической процедуры организации начальных вычислений, поскольку на длительных временах проявляется эффект сглаживания, характерный для диффузионных дифференциальных уравнений.
За новый отсчет времени (Ь = 0) удобно принять момент, когда в (Я — ро)-слое достигается концентрация и(Ь) = Q. Быстро образуется корка гидрида с концентрацией Qh. К этому моменту с(0, г) = г) (распределение растворенного водорода в а-фазе с предыдущего этапа), ф(р0) = Q, (Я — р0)-слой уже гид-ридный и сквозь него диффундирует растворенный Н. Речь о диффузионно подвижном атомарном водороде в отличие от химически связанного (гидрид является новым «самостоятельным» материалом). Обозначим концентрацию диффундирующих атомов Н в гидриде через у(Ь,г). Общая концентрация равна Qh + у(Ь,г).
В плане второстепенного уточнения можно для диффузионного уравнения (1) незначительное время продолжить вычисления с гра-
44
2
ничными условиями с(Ь, ро) = Q, дгс(Ь, го) = 0 (пока зерна срастаются в корку), чтобы «подправить» распределение ^>(г). Далее подобные переходные процессы не учитываем.
В цилиндре с растущей коркой гидрида (г = р(Ь) - граница раздела фаз, р(0) = р0) запишем диффузионные уравнения:
1 ду д 2у 1 ду , , , *
Б • ду = +1 • дТ' г е (4)
у(0,г) = 0, г е [ро'Щ, дс д 2с 1 дс
в • дь = дТ2 + Т • дТ' г е {го'р(г))' (5)
с(0, г) = р(г), г е [го,ро].
Здесь Б* - коэффициент диффузии Н в новом материале (гидриде). Граничные условия на «входе-выходе» запишем аналогично I:
2^s*p
дс дг
1 -
Го
v(t, R)
2 - b*v2(t, R) = -D* ^ дг
R
0.
[v(t,p(t)) + Qh - c(t,p(t))]dV
dv
dr
дс
2npH At + D—
p(t) дг
p(t)
2npHAt,
dV = H[np2 (t) - np2(t + At)] = -Hn2pp At + o(At), [v{t, p(t)) + Qh - c(t, p(t))] p(t)
D —
D* M
дг
- dдг
p(t) дг
p(t)
¿-фазы практически полностью уходит на формирование нового слоя гидрида (сдвиг границы г = р(Ь)) и в раствор. При с(Ь, р) <Q тонкий слой гидрида на границе раздела фаз становится неустойчивым и частично распадается, поэтому полагаем с(Ь, р) = Q.
Таким образом, принимаем следующие условия на свободной границе раздела фаз:
у(Ь, р(Ь)) = 0, с(Ь, р(Ь)) = Q, Ь > 0, (6)
дv
[Qh - сР] Р = -D* gr
+ D^
р дг
t ^ ts > 0.
(7)
В гидриде градиент концентрации возникает практически сразу, поскольку диффундировать сквозь ¿-фазу значительно труднее. Нет «накопительного» слоя (аналога (Я — ро)-слоя на этапе I). Когда баланс «приток-отток» станет положительным, граница раздела фаз г = р(Ь) станет подвижной и для «склейки» диффузионных уравнений потребуются дополнительные условия на стыке г = р(Ь).
Начнем с уравнения типа уравнения Стефана, описывающего динамику движения свободной границы раздела фаз:
Подчеркнем, что уравнение Стефана (7) «подключается» к модели, когда правая часть (баланс потоков) станет положительной. Так, условно выбран Ь = 0. До этого момента краевые задачи в слоях решаются численно независимо при р = ро. Фронт движется (Ь ^ Ь3) с общей концентрацией Qh в сторону г = го.В пределе имеем р(Ь) ^ го (Ь » 1). Уровень с(Ь,г) выравнивается, не превышая значение Q. Тонким «остаточным» слоем 2г — 1ЫЪ в расчетах уже можно пренебречь.
Преобразование модели: переход к безразмерному виду
Этап I. Сделаем замену времени t = (R -ro)2 D_1t', где (R-ro)2D_1 - характерное время диффузии, и замену независимой переменной r = ro + фо - ro], c(t',z) = c(t(t'),r(z)). Введем соответствующую нормировку концентраций и оставим прежнее обозначение для функции с: с := с/и, и := и/и, и = \j2^sp/b.
Краевая задача (1)-(3) после замены переменных примет следующую форму:
дс 7 д2 с 7 . . дс .
etc = di + ^ aZ' z e (0-1)' (8)
di s (R - , d2(z)= di(po - ro)
(po - ro)2 с(0, z) = 0, z e [0,1], дс
ro + z(po - ro):
c)z
z=o
= 0, a(t', 1) = U(t'),
(9)
dU 2/j.sp 2R(R - ro)2
Появляется разрыв концентраций:
у(Ь,р) + Qh > Qh >с(Ь,р), р = р(Ь).
При этом концентрация у(Ь, р) должна быть пренебрежимо малой (нет существенного «сопротивления» со стороны ^г-сплава). Принимаем у(Ь, р(Ь)) = 0. Поступающий поток из
dt'
— Ьи
+
и D(R2 - po)
[i - и]2
ис
2R(R - ro)2 D(R2 - po)
2po (R - ro)2 дс
po - ro R2 - P2 c)z
z=1
U(0) = 0. (10)
45
р
v
max -
Заметим, что вследствие равенства сомножителей при [1 — и]2 и и2 уравнение (10) можно переписать в виде
Вычислительный Алгоритм
^ = d3(1 - 2и) + d4 ^ , и(0) = 0, (11 dt' dz z=i
_ 2bu R(R - tq)2 D
d3 — —ft • —^-2—, d4 — d3 —
D
R2 - PO
Этап II. Сделаем замену времени t = ^ t'
Vdd
и пространственной переменной в слоях гидрида и раствора:
Этап I: растворение Н в Zг — 1NЬ
Следуя технике разностных схем, введем сетку [гт = тНх, т = 0,1,..., М} (Нг = 1/М) по пространственной переменной и сетку по времени {Ь'п = иН?,п = 0, 1,...}. Обозна-Ьи Я(р0 — г0)' чим через {Ст}, {ип} приближенные значения концентраций в (р0 — г0)-слое (с(Ь'п^т)) ив (Я — р0)-слое (и(Ь'п)) соответственно. Рассмотрим неявную схему для уравнения (8) и неявный метод Эйлера для ОДУ (11):
Ро
(4): t = p(t') + y[R - p(t')], v ^ V(t', y);
(5): t = tq + z[p(t') - To], c ^ c(t',z).
Введем соответствующие нормировки и оставим прежние обозначения для функций V := V/V, c := c/и, V = y/2^s*p/b*, и = y/2ßsp/b. Дополнительно введем функцию
= fy = d ы «<)•
где точка означает дифференцирование по t'.
Краевая задача (4)-(7) после замены переменных принимает следующий вид:
' R - To Y d2v dy2
cn+l cn cn+l — 2Cn+l + cvn+l
Cm - = d Cm-l 2Cm + cm+l
ht
+ d2(Zm)
h2
cvn+l _ cvn+l
cm+l cm-l
2hz
(17)
Uc n+l _ Uc n
ht
= d3 (1 - 2IJn+l)
+ d4
cJn+l _ лст+l + 3(Jn+l
CM—2 4CM-1 + 3CM
M-1
M
2h z
(18)
dV dt
R - p
+
X(t')(1 - y) + /D* (R - tq)2(R - p)-1
1 - Rp
l
D p + y(R - p)
dV dy, (12)
y e (0,l),J(0,y) = 0,y e [0,l],J(t',0) = 0,
r)v
d* [1 - 2 J(t', 1)= - -
b* v[R - p(t')]
, d* —
D
(13)
dc= [D iR - tq\2 d2J
dt' V D* V p - tq) ^ dZ2
Рассмотрим переход с п-го на (п + 1)-й слой по времени (п ^ 0, 0 < т < М): = Н2Н- 1й-1, ■^2(гт) = 0.5НХ (12(ХтУ1-1,
[1 — Ш2Ы]Ст+-\ — Ш1 +2]Ст+1 + [1+Ш2Ы]Ст++1+шСт = 0.
Значения в начальный момент времени известны: Ст = и0 = 0 (0 ^ т < М). Следуя методу прогонки (алгоритм Томаса), ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое по времени в виде
Сш1 = ат+1Ст++1 + вт+1, т = 0,...,М — 1.
Прогоночные коэффициенты следующие: (т = 1,...,М — 1),
+
\(t')z 1 - ToP-1
+
D (R - ro)2(p - tq)
l
To + z(p - To)
dc
dz'
am+l =
1 + W2(zm)
2 + Wl + (W2(z m) - 1) am
(14)
О dc
c(0,z)= ф(z),z e [0,1],c(t', 1) = О, dz
0,
(15)
0h_ c(+' i)i \(t) = _ DV (R - tq)2 dV
LU c(t, \(t )= V Du p(R - p) • dy
+
D (R - To)2 dc D* p(p - tq) dz
z=l
(t' > t's;
(16)
(1 — т т
Рт+1 = 7Т----ч-ТТ-. (19)
2 + + (Ш2^т) — 1)ат
Для нахождения начальных коэффициентов а1, в1 воспользуемся следующими соображениями. Подсчитаем предварительно значения С'п+1 по явной разностной схеме (в равенстве (17) справа заменяем п +1 на п). Рассмотрим первое краевое условие непроницаемости при z = 0 в (9). Используя разностную
46
o
аппроксимацию дхс\х=о & [—ЭС^1 + 4Сп+1 — С'п+1]/2Нг = 0, получаем выражение С'П+1 = /1(СП+1, С21+1) = (4СП+1 — С2(+1)/Э. Зная численное значение С^1 и выражение (7,П+1 =
тС'П+1 + вь получаем а1 =0, в1 = С'П+1. По а1, в1 вычисляем оставшиеся коэффициенты ат, вт, т = 2,...,М по формулам (19).
Ближайшая цель - найти значение СМ+1, необходимое для реализации обратного метода прогонки. Рассмотрим разностную аппроксимацию уравнения (11):
Ц/п+1 — иг'
ы>
= дз(1 — 2£/"+1)
Сп+1 _ 4Сп+1 + ЭСп+1 См-2 4СМ-1 + ЭСМ
+ d4
2НХ
Используя соотношение Ст+1 = ат+1Ст+11 + вт+1 при М — 2, М — 1 и второе условие (9) (ип+1 = СМ+1), получаем СМ+1 = АА
1
А = 2НХ [ип + Ы, дз]
+ К? [(ам-1 — 4)вм + вм-1 ], А-2 = 2НХ [1 + 2Нг> дз] — Ы> д4[(ам -1 — 4) ам + 3].
Следующий этап: с текущими приближениями значений СП+1, СП+1 решаем обратным ходом прогонки трехдиагональную систему линейных уравнений и находим новые приближения концентраций СП+1 (и остальные значения СП+1 для т = Э,...,М — 1). Далее снова пользуемся формулой СП+1 = /1(Сп+1, С%+1). После этого корректируем значения прогоноч-ных коэффициентов (19), определяем С^1 по приведенной выше формуле и повторяем вычисления (возвращаясь к началу абзаца) до установления граничных значений С'п+м (обычно 2-3 итерации). Критерием окончания вычислений на этапе I выбрано условие ^п+1 ^ Q/u = 0, 98.
Результат окончания I этапа: образовался гидридный слой толщины (Я — ро) с концентрацией химически связанного водорода Qh, а в слое раствора толщины (ро — го) имеется распределение концентрации растворенного водорода ^(г), р(ро) = Q.
Этап II: гидридообразование и движение границы фазового перехода
Начальное распределение ф(г) (р(г)) определяется первым этапом. Формально можно забыть о физическом смысле функции р(Ь') как границы раздела фаз и рассматривать ее
как функциональный параметр. По решению р(Ь') (зная Л(Ь7)) определяются коэффициенты модели. Поэтому вычислительный алгоритм основан на неявных разностных схемах и носит итерационный характер. Итерации будут связаны с уточнением значения Л(Ь') на каждом временном слое.
Введем следующие сетки: {ут = тНу, т = 0,1,...,М}, Ну = 1/М - шаг по пространственной переменной у; {г^ = кНг, к =
0.1...., К}, Нх = 1/К - шаг по пространственной переменной г; {Ь'п = иН?, п = 0, 1,...},
H, - шаг по времени Ь'. Обозначим через {1^}, {Сп} приближенные значения концентраций в гидридном (Я — р(Ь'))-слое (У(Ь'п,ут)) и в (р(Ь') — го)-слое раствора (с(Ь'п,гк)).
1. Начальный этап: переход с нулевого на первый слой и нахождение Л > 0
При р = ро будем рассматривать последовательно и независимо две краевые задачи: сначала для гидридного слоя, потом для слоя раствора в сплаве Zг — 1ЖЬ. Граница раздела фаз на данном этапе неподвижна, полагаем в (12), (14) Л = 0 (р = 0). Опишем алгоритм перехода с нулевого на первый слой по Ь'.
I.1 Гидридный слой
Рассмотрим краевую задачу (12)-(13) при Л = 0, р = ро:
ду д'
д2у
ш = д1 • + д2(у)
ду
^ У е (0,1), (20)
д1 — \ ——
д2(У) =
д1(Я — ро) ро + у(Я — ро),
У(0,у) = 0, у е [0,1], у(Ь', 0) = 0,
д*[1 — 2 у^ 1)] = — ду
= Ь*у [Я — ро]
, д* — г-. •
1 В*
(21)
Для уравнения диффузии (20) рассмотрим неявную разностную схему на первом слое: т = 1,...,М — 1,
и1 - V0
у т у т
ы>
= д
1 •
и 1 - 2и1 + V1
т— 1 т 1 г т
т+1
Н2у
V 1 ,, - V1
+ д2(ут) • т+1 т-1. (22)
Рассмотрим уравнения перехода с нулевого на первый слой по времени (0 < т < М): Wl — КуН-, 1д-1, т2(ут) — 0.5Нуд2(ут)д-1,
[1 — W2(Уm)]Um -1 — [2+ ^1]ит + [1 + ■и2(ут)\и1+1 + WlUm = 0.
Начальные данные: ит — 0 (0 ^ т ^ М). Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки в виде
1т = ат+1Уг1+1 + вт+1, т = 0,...,М — 1. Прогоночные коэффициенты (т = 1,..., М—1):
am+1 —
ßm+1 —
1 + w2(ym)
2 + Wi + (W2(y т ) - 1)am
(1 - W2(y ))ß m
2 + Wi + (W2(y m ) - 1)am
(23)
Для нахождения начальных коэффициентов а1, в1 воспользуемся краевым условием у(Ь',0) = 0. Зная значение и1 = 0 и выражение
V)1 = а1 V1 + в1, получаем а1 = 0, в1 = V,1 = 0. По а1, в1 вычисляем оставшиеся прогоночные коэффициенты ат, вт, т = 2,...,М, по формулам (23).
Теперь найдем значение Ум, необходимое для реализации прогонки. На первом слое по Ь' аппроксимируем производную
dv dy y=i
Vi-2 - 4Vl_i + 3 V
1
AI
2hy
и подставляем вместо 1 выражения
УЦ = аи+1У1+1 + вк+1 при соответствующих к. Рассмотрим краевое условие (21):
d* [1-2VM] — -
aM-1aMVM + aM-1ßM
2hy
+
ßM-1 - 4(амV1 + ßM) + 3V1
м
2hy
В итоге получаем:
V1 — (4 - aM-1 )ßM - ßM-1 - 2hy d* M (aM-1 - 4)aM + 3 - 2hy d* '
Используя VM, решаем обратным ходом прогонки трехдиагональную систему линейных алгебраических уравнений и находим значения концентрации V^, m — 1,...,M - 1.
1.2 Слой раствора Zr-1Nb
Переходим к краевой задаче (14)—(15), возвращаясь к предыдущему слою(t'0 — 0):Ло — 0,
д с д2 с д с
- — d1W + <Ы*)-Тг, г e (0,i) (24)
d1 — \ тг"
d2(z) —
d1(po - го) ro + г(ро - ro):
О дс
с(0,г) — ^(г), с(f, 1) — О, дг
z=0
— 0. (25)
Для уравнения (24) рассмотрим неявную схему на первом слое: к = 1,...,К — 1,
1 о С - С — d1
Ch - 2Ck + С1
fc+1
ht
hz2
1 1 + d2(zfc)• Cfc+12h C k-1. (26)
Уравнения перехода с нулевого на первый слой по времени (0 < к <К): w1 — Н^Н- д-1, W2(zk) — 0.5Нг д2(гк )д-1,
[1 — W2(Zk )]Ск_1 — [2 + Wl] С к + [1 + W2 (гк )]Ск+1 + WlCk = 0.
Начальные данные: Ск = ф(гк) (0 ^ к ^ К). Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки в виде С1к = ак+1Ск+1 + вк+1, к = 0,...,К — 1. Про-гоночные коэффициенты (к = 1,...,К — 1):
afc+1
ßk+1 —
1 + W2 (г к )
2 + W1 + (W2(zk) - 1)afc:
(1 - W2 (гк))ßk + W1 СО 2 + W1 + (w2(zo) - 1)ak'
(27)
Для нахождения начальных коэффициентов а1, в1 подсчитаем предварительно значения С^2 по явной разностной схеме (в равенстве (26) справа заменяем 1 на 0). Рассмотрим краевое условие непроницаемости при г = 0 в (25). Используя разностную аппроксимацию
дси=о & [—ЭСо1 + 4С-1 — С2,\/2НХ = 0, получаем выражение Со = /1(С1,С1) = (4С1 — С21)/3. Зная численное значение Со и выражение Со = а1С1 + в1, получаем а1 = 0, в1 = С^. По а1, в1 вычисляем оставшиеся прогоночные коэффициенты ак, в к, к = 2,...,К по формулам (27).
Значение СК, необходимое для реализации прогонки, находим из краевого условия при г = 1 в (25): С1К = Q/u.
С текущими приближениями значений С, СК1 решаем обратным ходом прогонки трех-диагональную систему линейных уравнений и находим новые приближения концентраций С11 2 (и остальные значения СК1 для к = Э,...,К — 1). После этого возвращаемся к уточнению Со = /1(С1,С21), корректируем значения прогоночных коэффициентов (27) и повторяем вычисления (возвращаясь к началу
48
абзаца) до установления граничных значений С1 к (обычно 2-3 итерации).
1.3 Нахождение Л0 > 0
После указанных вычислений аппроксимируем производные на первом слое
дуУ\у=о и [—Щ1 + Щ1 — Щ\/2НУ,
дхс\2=1 и [СК_2 — 4СК-1 + 3СК]/2Нг
и подставляем их в граничное условие (16):
Qh v 1
L U
K
\о
V (R - ro)2 -Щ1 + Щ1 - V21
D u po(R - po)
2hb
+ iD. (R - ro)2
D* po(po - ro)
cK_2- 4CK~i + 3<VK
2h*
Отсюда находим приближение Ло. На данном начальном этапе возможны два варианта.
1. Л0 ^ 0 (на начальном этапе распределение V практически нулевое по сравнению с распределением с и линейная комбинация градиентов в (16) неположительна). По физическому смыслу задачи полагаем Ло = 0 и возвращаемся в начало излагаемого алгоритма (пункт1.1), используя в качестве начальных данных текущие распределения концентраций УШ = Ут, т = 0,...М, С0 = С1, к = 0,...,К.
2. Л0 > 0. Принимаем это число в качестве уточненного значения Л0 и завершаем итерации. Далее вычисляем р(0) = Л0р0 ив силу
(
Л(Ь0 ) = (1п р(Ь%,
^ \п(р(г'0 + аь')) и 1п р(г'0) + Л(г'0) • аь'
определяем 1пр(Ь'1), = Ь'0 + АЬ'. По значению 1п р(Ь'1) находим р(Ь'1) и переходим к пункту 2, используя в качестве начальных данных распределения УШ, т = 0,...М, С1, к = 0,...,К.
2. Переход от п-го к (п + 1)-му слою по Ь'
Опишем алгоритм нахождения распределений на (п + 1)-м слое (п > 1) у(Ь'п+1,у) [УШ+1,
т = 0,...,М], с (Ь'п+1, z) [С П+1, к = 0,...,К], зная распределения с предыдущего слоя по времени V(Ь'п,у), с(Ь'п^). Рассмотрим последовательно два слоя материала: сначала гидрид-ный, затем слой сплава Zг — 1NЬ.
Начальное приближение Л(Ь'п) и Л(Ь'п_ берется с предыдущего временного слоя. Параметр Л(Ь'п) и распределения у{Ь'п+1,у),
с(Ь'пподлежат итерационному уточнению, значение р(Ь'п) известно.
2.1 Гидридный слой
Для уравнения диффузии (12) в безразмерной форме
§ = (1(Ь') • дУу + (2(И,У) • ду, у € ^ 1),
di(t) = */ D
R - ro R - p(t/)_
, d2(f,y)
Kt')(l - y) , D(R - ro)2(R - p(t'))
1
+
1 - Rp-1(t/)' V D p(t/) + y(R - p(t'))
рассмотрим неявную разностную схему на (n+ 1)-м слое: m = 1,...,M - 1,
VU+1 _ Vn Vn+1 - 2Vn+1 + Vn+1
'■m. vm , /./ n r m-1 rm 1 r m+1
ht
= d1(t'n)-
h2y
Vn+1 _ Vn+1 + d2(t'n,ym) ■ m+\- m-1. (28
2hy
Рассмотрим уравнения перехода с n-го на (n + 1)-й слой по времени (n ^ 1, 0 < m < M):
wl(t'n) = h2yh- '^Ю, W2(t'n,ym) — 0.5hy (l2(t'n,ym)d-l(t'n),
[1 - W2(t'n,ym)]Vm+l - [2+ Wl(t'n)]Vm+l + [1 + W2(t'n, ym)]VmX\ + Wl Vm = 0.
Начальные данные Vm берутся из пункта 1.3. Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки в виде
Vm+l = am+lVm+l + ßm+l, m = 0,...,M - 1.
Прогоночные коэффициенты (m = 1,..., M-1):
a i = _1 + W2(t'n,ym)_
m+l 2+ Wl(t'n) + (W2(t'n,ym) - 1)am '
в = (1 - W2(t'n,ym))ßm + Wl(t'n) Vm (2g) m+l 2 + Wl(t'n) + (W2(t'n,ym) - 1)am ' Для нахождения начальных коэффициентов al, ßl воспользуемся краевым условием V(t', 0) = 0. Зная значение Vn+l = 0 и выражение Vn+l = alVln+1 + ßl, получаем al = 0, ßl = Vn+l = 0. По al, ßl вычисляем оставшиеся прогоночные коэффициенты am, ßm, m = 2,..., M, по формулам (29).
Теперь найдем значение V'M+l, необходимое для реализации прогонки. На (n + 1)-м слое по t' аппроксимируем производную
dV dy y=1
Vn+1 4V
n+1 M-1
+ 3VM
n+1
M
2hy
49
2
и подставляем вместо им+2 м-1 выражения
Vkn+1 = ак+1икп+11 + вк+1 при соответствующих к. Рассмотрим краевое условие (13):
А [1 от>п+11 ам-1амимм+1 + ам-1вм д* [1—2УM \ =--
2 Ну
+
вм-1 — 4(ам !м+1 + вм) + ЭУм+1
' м_
2 Ну
Л(Ь')г
1 — гор-1(Ь')
+
Б (Я — го)2(р(Ь') — го)
1
го + г(р(Ь') — го)
рассмотрим неявную разностную схему на (п+ 1)-м слое: к = 1,...,К — 1,
сс п+1 — Сп Сп+1 — 2 С п+1 + Сп+1
ы
= д1(0 •
Н2
ак+1Сп+11 + вк+1, к = 0,...,К — 1. Прогоноч-ные коэффициенты (к = 1, 2,...,К — 1):
ак+1
1 + W2(t'n,Zk)
2 + Wl(tn) + ^2(Ь'п ,гк) — 1)ак'
В итоге получаем: ип+1 = (4 — ам-1)вм — вм-1 — 2Ну д*Ю
м = (ам-1 — 4)ам + 3 — 2Ну д*(Ь'п) '
Используя У^1, решаем обратным ходом прогонки трехдиагональную систему линейных алгебраических уравнений и находим значения концентрации ит+1, т = 1,...,М — 1. .
2.2 Слой раствора Zг — 1NЬ
Переходим к краевой задаче (14)-(15), возвращаясь к предыдущему слою (Ь' = Ь'п).
Для уравнения диффузии (14) в безразмерной форме
д д2 д Ж = д1 С) а? + д2ы г е (0,1),
"'«^цБ ( Щ—оо, )2
С п+1 _ С п+1
+ д2(Ь'п,гк)• к+12^ к-1 . (30)
Рассмотрим уравнения перехода с п-го слоя по Ь' на (п + 1)-й слой (п ^ 1, 0 < к < К):
Wl(tn) — Н2К- 1д-1(Ьп),
W2(t'n,Zk) — 0.5НХ д2(Ь'п,гк )д- 1(Ь'п),
[1 — W2(t'n, гк)]Сп+11 — [2 + Wl(tn)]Cn+1 + [1 + W2 (С гк )]Сп+11 + Wl(tn)Cn = 0.
Ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое в виде Сп+1 =
о = (1 — W2(t'n,Zk))вк + Wl(tn) Сп вк+1 2 + Wl(tn) + (w2(tn г) — 1)ак. 1 ;
Для нахождения начальных коэффициентов а1, в1 подсчитаем предварительно значения Спп+1 по явной разностной схеме (в равенстве (30) справа заменяем п + 1 на п). Рассмотрим краевое условие непроницаемости при г = 0 в (15). Используя разностную аппроксимацию
дгси=о & [—ЭС^+1 +4Сп+1 —Сп+1]/2Нг = 0, получаем выражение С5+ = /)(Сп+1, С21+1) = (4Сп+1 — Сп+1)/Э. Зная численное значение Сп+1 и выражение С^1 = аС'^1 + в1, получаем а1 =0, о, = Сп+1. По а1, о, вычисляем оставшиеся прогоночные коэффициенты ак, вк, к = 2,...,К по формулам (31).
Значение СК+1, необходимое для реализации прогонки, находим из краевого условия
при г = 1 в (15): СК+1 = Q/u■
С текущими приближениями значений С(п+1, СК+1 решаем обратным ходом прогонки трехдиагональную систему линейных уравнений и находим новые приближения концентраций Сп+1 (и остальные значения СК+1 для к = Э,...,К — 1). После этого возвращаемся к уточнению С^1 = ЛС^1, Сп+1), корректируем значения прогоночных коэффициентов (31) и повторяем вычисления (возвращаясь к началу абзаца) до установления граничных значений итерации).
2.3 Итерационное уточнение Л(Ь'п)
После указанных вычислений аппроксимируем производные на (п + 1)-м слое
дуУ \у=о & [—Э^^1 + 4ип+1 — У2n+1]/2Нy,
дгС\г=1 & [СК+-2 — 4СК+-1 + ЭСК+1]/2Нг
и подставляем их в граничное условие (16):
Qh Сп+1
и К
л(ь'п) = —/ Б* и
(Я — го)2
Б и р(Ь'п)(Я — р)
—ЗУ^1 + 4У1п+1 — У^1
2 Ну
+
Б (Я — го)2 СК+-\ — 4СК+-\ + ЭСК+1 Б* р(р — го) 2Нг
®
Отсюда находим новое приближение Л(Ь'п) > 0 и возвращаемся в начало излагаемого алгоритма (пункт 2.1). Через несколько итераций получаем установившееся значение Л(Ь'п).
Далее вычисляем значение производной
рЮ = лЮрЮ и в силу (
Л(Ь'п) = ( 1п р(Ь%п
^ 1п(р(Ь'п + АЬ')) и 1п р(Ь'п) + Л(Ь'п) • АЬ'
определяем 1п р(Ь'п+1), Ь'п+1 = Ь'п + АЬ'. По значению 1п р(Ь'п+1) находим р(Ь'п+1) и переходим к слою Ь'п+2, приняв в качестве начального приближения Л(Ь'п+1) и Л(Ь'п). Далее уточняем значение Л(Ь'п+1) и распределения у(Ь'п+2,у), с(Ь'п+2^) по изложенной схеме (пп. 2.1-2.3).
Заключение
Одним из важных требований к изделиям из циркониевых сплавов активной зоны реакторов является низкое поглощение водорода, поскольку водородное (гидридное) охрупчива-ние может стать одной из причин разрушения циркониевой оболочки.
При проведении экспериментальных исследований практически невозможно определить распределение водорода в материале в зависимости от времени, поскольку атомарный растворенный водород является чрезвычайно подвижной фазой внедрения.
В статье представлена модель гидрирования цилиндрической оболочки из циркониевого сплава, учитывающая динамику абсорбционно-десорбционных процессов и движение свободной границы раздела фаз (гидрид-металл). На основе неявных разностных схем разработан итерационный вычислительный алгоритм решения краевой задачи.
Финансирование исследований осуществлялось из средств федерального бюджета на выполнение государственного задания КарНЦ РАН (Институт прикладных математических исследований, Отдел комплексных научных исследований КарНЦ РАН).
Литература
1. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.
2. Водород в металлах / Ред. Г. Алефельд, И. Фелькль. М.: Мир, 1981. Т. 1, 2. 506, 430 с.
3. Заика Ю. В., Родченкова Н. И., Гру-дова К. В. Численное моделирование динамики свободной границы гидридообразования // Труды КарНЦ РАН. 2018. № 7. С. 14-24. ао1: 10.17076/та1843
4. Изотопы водорода. Фундаментальные и прикладные исследования / Ред. А. А. Юхим-чук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697 с.
5. Писарев А. А., Цветков И. В., Марен-ков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
6. Черданцев Ю. П., Чернов И. П., Тюрин Ю. И. Методы исследования систем металл-водород. Томск: ТПУ, 2008. 286 с.
7. Denisov E. A., Kompaniets M. V., Kompaniets T. N, Bobkova I. S. Peculiarities of hydrogen permeation through Zr-1%Nb alloy and evaluation of terminal solid solubility // Journal of Nuclear Materials. 2016. Vol. 472. P. 13-19. doi: 10.1016/j.jnucmat.2016.01.022
8. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
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. 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
11. Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface // International Journal of Hydrogen Energy. 2011. Vol. 36, iss. 1. P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
12. The hydrogen economy / Eds. M. Ball, M. Wietschel. Cambridge: Cambridge University Press, 2009. 646 p.
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. 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
15. 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.
Поступила в редакцию 03.06.2019
51
References
1. Vzaimodeistvie vodoroda s metallami [Interactions of hydrogen with metals]. Moscow: Nauka, 1987. 296 p.
2. Vodorod v metallakh [Hydrogen in metals]. Moscow: Mir, 1981. Vol. 1, 2. 506, 403 p.
3. Zaika Yu. V., Rodchenkova N. I., Grudova K. V. Chislennoe modelirovanie dinamiki svobodnoi granitsy gidridoobrazovaniya [Numerical simulation of the dynamics of the free boundary of hydride formation]. Trudy KarNTs RAN [Trans. KarRC RAS]. 2018. No. 7. P. 14-24. doi: 10.17076/mat843
4. Izotopy vodoroda. Fundamental'nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies]. Sarov: RFYaTs-VNIIEF, 2009. 697 p.
5. 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.
6. 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.
7. Denisov E. A., Kompaniets M. V., Kompaniets T. N., Bobkova I. S. Peculiarities of hydrogen permeation through Zr-1%Nb alloy and evaluation of terminal solid solubility. Journal of Nuclear Materials. 2016. Vol. 472. P. 13-19. doi: 10.1016/j.jnucmat.2016.01.022
8. Handbook of hydrogen storage: new materials for future energy storage. Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
СВЕДЕНИЯ ОБ АВТОРАХ:
Родченкова Наталья Ивановна
старший научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований КарНЦ РАН, Федеральный исследовательский центр «Карельский научный центр РАН» ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (8142) 766312
Грудова Ксения Васильевна
младший научный сотрудник Отдел комплексных научных исследований КарНЦ РАН, Федеральный исследовательский центр «Карельский научный центр РАН» ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (8142) 766312
9. Lototskyy M. VYartys 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. 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
11. Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface. International Journal of Hydrogen Energy. 2011. Vol. 36, iss. 1. P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
12. The hydrogen economy. Eds. M. Ball, M. Wietschel. Cambridge: Cambridge Univ. Press, 2009. 646 p.
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. 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
15. 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.
Received June 06, 2019
CONTRIBUTORS:
Rodchenkova, Natalia
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: [email protected] tel.: (8142) 766312
Grudova, Kseniia
Department of Multidisciplinary Scientific Research,
Karelian Research Centre,
Russian Academy of Sciences
11 Pushkinskaya St., 185910 Petrozavodsk,
Karelia, Russia
e-mail: [email protected]
tel.: (8142) 766312
