Численное моделирование динамики свободной Гранины гидридообразования Текст научной статьи по специальности «Физика»

Труды Карельского научного центра РАН № 7. 2018. С. 14-24 DOI: 10.17076/mat843

УДК 519.6:539.2

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ СВОБОДНОЙ ГРАНИЦЫ ГИДРИДООБРАЗОВАНИЯ

Ю. В. Заика, Н. И. Родченкова, К. В. Грудова

Институт прикладных математических исследований КарНЦ РАН, ФИЦ «Карельский научный центр РАН», Петрозаводск, Россия

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

Ключевые слова: гидрирование; нелинейные краевые задачи со свободной границей; разностные схемы; численное моделирование.

Yu.V. Zaika, N.I. Rodchenkova, K. V. Grudova. NUMERICAL MODELLING OF DYNAMICS OF FREE BOUNDARY OF HYDRIDE FORMATION

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. The study was carried out under state order to the Karelian Research Centre of the Russian Academy of Sciences (Institute of Applied Mathematical Research KarRC RAS).

K e y w o r d s: hydrogenation; nonlinear boundary-value problems with free phase boundary; difference schemes; numerical simulation.

Введение

Интерес к взаимодействию изотопов водорода с конструкционными материалами носит многоплановый характер [1-5, 8, 9, 11]. Энтузиасты говорят не только об энергетике, но и о водородной экономике [11]. Некоторые модели дегидрирования и водородопроницаемости, связанные с тематикой данной статьи, исследованы в [7, 12-15]. Одним из важных требований к изделиям из циркониевых сплавов активной зоны реакторов является низкое поглощение водорода, поскольку водородное охруп-чивание может стать одной из причин разрушения циркониевой оболочки. В зависимости от уровня содержания водорода и температуры водород может находиться в циркониевых сплавах в виде твердого раствора или в виде гидридов. Наибольший охрупчивающий эффект оказывают гидриды, так как они могут служить участками образования трещин.

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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ГИДРИРОВАНИЯ

Вначале кратко опишем условия эксперимента (подробнее см. [10]). Пластина из сплава Zr — 1МЪ шлифуется с одной стороны, другая сторона практически водородонепроница-ема, торцами пренебрегаем; температура образца Т и давление газообразного водорода р поддерживаются постоянными (предпринимаются специальные меры охлаждения).

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

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

Обозначим: L — толщина пластины; ¿о — толщина слоя, в который водород абсорбируется относительно легко и еще не ощущается диффузионное сопротивление (будущая начальная корка гидрида); u(t) — концентрация H в 4)-слое (1н/см3); Q — концентрация, по достижении которой локально решетка перестраивается и возникают зародыши гидрид-ной фазы; Qh — концентрация атомов водорода в гидридной фазе (химически связанный водород, образующий гидрид как вещество); c(t, x) — концентрация растворенного H в (L — ^о)-слое; у — газокинетическая константа. Температура пластины и давление напуска водорода постоянны (Т = const, р = const).

Согласно кинетической теории газов плотность Jp падающего на поверхность потока частиц (в данном случае молекул H2) связана с давлением p по формуле Герца - Кнудсе-на: Jp = рД/2nmkT (k — постоянная Больц-мана, m — масса молекулы H2). В контексте эксперимента удобно в качестве единиц измерения выбрать [x, £, L] = см, [р] = торр. Тогда численно получаем зависимость Jp = ур, у(Т) и 2, 474 ■ 1022/VT ([у] = 1н2/(торрсм2 с), [Т] = К, под корнем безразмерное численное значение). Поскольку диффундирует атомарный водород, то для единообразия подсчет будем вести в атомах H: Jp = 2ур. Только малая часть H окажется в абсорбированном состоянии: Jabs = 2ysp (s ^ 1). Множитель s имеет смысл доли налетающих H, которые оказались в приповерхностном объеме. Объединяем более элементарные стадии физадсорбции, диссоциации и растворения в одну: s — эффективный коэффициент абсорбции.

Этап I: растворение H в Zr — 1Nb

Для диффузионного слоя толщины (L — i0) имеем стандартную краевую задачу:

dc d2c

- = DdX2. x ^ «Ь-Ч. t> 0 (1)

c(0,x) = 0, x e [£o,L],

c(t,4)= u(t), 3xc(t,L) = 0. (2)

Граничное условие c(t,i0) = u(t) отражает непрерывность распределения H в сплаве, а dxc|L = 0 — непроницаемость стороны пластины x = L. Здесь и в дальнейшем считаем, что коэффициенты подчиняются закону Аррениуса по температуре, в частности D =

D0 exp{ — Ed/[RT]}. В течение одного эксперимента поддерживается T = const.

Для концентрации u(t) запишем ОДУ, исходя из баланса потоков:

U(t)4 = 2^sp — bu2 + , u(0) = 0. (3)

Содержательный смысл: за 1 с через 1 см2 абсорбировалось за счет давления 2^sp атомов H, но есть встречный поток десорбции bu2 (b — эффективный коэффициент рекомбинации) и диффузионный отток. Рассогласование плотностей этих потоков идет на накопление атомов водорода в ^о-слое (й£0). Уравнение (3) нужно рассматривать совместно с (1)— (2), поскольку u(t) определяет граничную концентрацию в (2). При небольшом p (без образования гидрида) в равновесии (когда производные равны нулю) имеем

2^sp — bu2 = 0 ^ u = Г/p, Г = л/2^sb-1.

Следовательно, динамика (3) в статике согласуется с законом Сивертса u a /р, Г — коэффициент растворимости. Подчеркнем, что речь — о растворенном атомарном диффузионно подвижном водороде. В эксперименте «насыщение-дегазация» учитывается общее поглощение водорода, включая обратимый захват и гидридные фазы — коэффициент Г может иметь другой смысл и численное значение. Технически нетрудно учесть обратимый захват H в (L —^0)-слое дефектами материала, но в рассматриваемой задаче считаем ловушки второстепенным фактором.

Этап II: гидридообразование и движение границы фазового перехода

По достижении концентрации u(t) порогового уровня Q = Q(T) происходит образование зародышей гидрида (¿-фазы). Дальнейшее поступление абсорбированного атомарного водорода расходуется на рост зародышей до образования сплошной корки гидрида (с концентрацией Qh > Q химически связанного водорода) и поддержание уровня c(t,i0) = Q в а-фазе раствора. При этом с входной стороны зародыши растут и смыкаются быстрее, общая концентрация водорода уже начинает превышать уровень Qh за счет дополнительного растворения в гидриде. Когда приток к границе гидридного микрослоя x = ^о начинает превышать отток в сплав Zr — 1Nb, появляется «движущая сила», и граница фазового перехода начинает смещаться, слой гидрида растет.

Считаем, что этот переходный процесс относительно быстрый. Разумеется, речь не о

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

За новый отсчет времени (Ь = 0) удобно принять момент, когда в ¿о-слое достигается концентрация и(Ь) = Быстро образуется корка гидрида с концентрацией Qh. К этому моменту с(0, х) = ^>(х) (распределение растворенного водорода в а-фазе с предыдущего этапа), ^(¿0) = Q, ¿о-слой уже гид-ридный и сквозь него диффундирует растворенный Н. Речь о диффузионно подвижном атомарном водороде в отличие от химически связанного (гидрид является новым «самостоятельным» материалом). Обозначим концентрацию диффундирующих атомов Н в гидриде через -и(Ь,х). Общая концентрация равна Qh + у(Ь,Х).

В плане второстепенного уточнения можно для диффузионного уравнения (1) незначительное время продолжить вычисления с граничными условиями с(Ь,^о) = Q, дхс(Ь, Ь) = 0 (пока зерна срастаются в корку), чтобы «подправить» распределение <^(х). Далее подобные переходные процессы не учитываем.

В пластине с растущей коркой гидрида (х = ¿(Ь) — граница раздела фаз, ¿(0) = ¿0) запишем диффузионные уравнения:

I = D.g, x е (°,ед),

v(0,x)=0, x е [0,4],

Ж = D£, x е №),L),

(4)

(5)

c(0,x) = <^(x), x е [l0,L].

Здесь Б* — коэффициент диффузии Н в новом материале (гидриде). Граничные условия на «входе-выходе» запишем аналогично I:

dv

2^s*p — b*v2(t, 0) = — D* —

ж=0

dc dx

0.

В гидриде градиент концентрации возникает практически сразу, поскольку диффундировать сквозь ¿-фазу значительно труднее. Нет «накопительного» слоя (аналога ¿о-слоя на этапе I). Когда баланс «приток-отток» станет положительным, граница раздела фаз х = ¿(Ь)

16

L

станет подвижной и для «склейки» диффузионных уравнений потребуются дополнительные условия на стыке х = ¿(*).

Начнем с уравнения типа уравнения Стефана, описывающего динамику движения свободной границы раздела фаз:

KM(í)) + Qh — c(t,¿(t))]¿(t) =

dv

dx

+ ^^

£(t) dx

m

Появляется разрыв концентраций:

и(м) + д >д >с(м), ¿ = ¿(*).

При этом концентрация должна быть

пренебрежимо малой (нет существенного «сопротивления» со стороны ^г-сплава). Принимаем ¿(¿)) = 0. Поступающий поток из ¿-фазы практически полностью уходит на формирование нового слоя гидрида (сдвиг границы х = ¿(*)) и в раствор. При с(^) < д тонкий слой гидрида на границе раздела фаз становится неустойчивым и частично распадается, поэтому полагаем = д.

Таким образом, принимаем следующие условия на свободной границе раздела фаз:

•и(м(*)) = 0, с(м(*)) = д, * ^ 0, (6)

dv

[Qh - Q]¿(í) = —D ^

+ ^^ ¿ dx

(7)

* ^ > 0. Подчеркнем, что уравнение Стефана (7) «подключается» к модели, когда правая часть (баланс потоков) станет положительной. Так, условно выбран £ = 0. До этого момента краевые задачи в слоях решаются численно независимо при ¿ = ¿о. Фронт движется (* ^ с общей концентрацией О, в сторону х = Ь. В пределе имеем ¿(*) ^ Ь (* » 1). Уровень с(£, х) выравнивается, не превышая значение д. Тонким «остаточным» слоем в расчетах уже можно пренебречь.

Преобразование модели: переход к безразмерному виду

Этап I. Сделаем замену времени t = L2 D-11', где L2D-1 - характерное время диффузии, и замену независимой переменной x = ¿o + z[L — ¿o], c(t', z) = c(t(t'),x(z)). Введем соответствующую нормировку концентраций и оставим прежнее обозначение для функции с: с := с/u, и := u/u, U = A/2ysp/b.

Краевая задача (1)-(3) после замены переменных примет следующую форму:

дс dt7

L

LL — ¿o

2 д 2 с

^ С(0, Z) = 0, (8

z е [0- 1Ь 1

z=1

= 0, c(t', 0) = U(t'), (9)

L2

L2u

du

dt7 = D¿ou D¿o

bu2 +

+

L2

дс

¿o[L — ¿o] dz

z=0

, U(0) = 0.

(10)

Заметим, что вследствие

2yspL2[D¿ou]-1 = bL2u[D¿o]

2

1

уравнение (10) можно переписать в виде

^ = d1 (1 — u2) + d2 Iе , u(0) = 0, (11)

dt' dz z=o

= ^[£¿0]-1, ¿2 = Ь2[¿о(Ь — ¿о)]-1. Этап II. Сделаем замену времени

Г2 .

* = ^дв* и пространственной переменной:

(4): х = ¿(*')У, V ^ ,у);

(5): х = ¿(¿')+ — ¿(¿О], с ^ с(^).

Введем соответствующие нормировки и оставим прежние обозначения для функций V := с/г;, с := с/и, V = л/2^*р/Ь*, и = л/2^р/Ь. Дополнительно введем функцию

л<(/> * щ = £¿^

где точка означает дифференцирование по

Краевая задача (4)-(6) после замены переменных принимает следующий вид:

д-с /D* / L\2 д2-с ' д-с

dt7 v "DVTJ • ду2 +A(t)y• ду-

t' > 0, у е (0,1), -(0, у) = 0, v(t', 1) = 0,

(12)

д-с

d* f1 — г;2(t', 0)] = — -

M(t')v , d* — n ,

y=o D*

(13)

дс

IÍ7

L i2 д2с+ A(t')(1 — z) дс

U — ¿J Iz2 L¿-1(t') — 1 Iz'

(14)

Q дс

e(0,z) = ^(z), é(t', 0) = Q, ^ =0, (15)

и az z=1

Qh — a(t', 0)

u

D

+ W-ÍT- •

A(t')= WD*-ГLf д-с A(t) = —V"DTUliJ ду

+

L2

дс

Dt ¿(í')[L — ¿(t')] 9z

z=o

(t' > ts). (16)

1

Вычислительный Алгоритм

Этап I: растворение Н в —

Следуя технике разностных схем, введем сетку {гт = тЛ, т = 0,1,..., М} (Л^ = 1/М) по пространственной переменной и сетку по времени {Ь'п = пЛ^, п = 0, 1,...}. Обозначим через {Ст}, |С/п} приближенные значения концентраций в (Ь — ¿0)-слое (с(ЬП, £т)) и в ¿0-слое (и(ЬП)) соответственно. Рассмотрим неявную схему для уравнения диффузии (8) и неявный метод Эйлера для ОДУ (11)

<n+1 _ <r

hf

U П+1 _ U r.

L

LL—¿0

2 Ara+1 oAra+1 i Ara+1 2 <т- 1 — 2<т + Cm+1

h2

ht

= d1[1 — (f>r+1)2] +

т+1

(17)

(18)

+ d2

—3<C(?r+1 + 4Cr+1 — <r+1 2h '

Рассмотрим уравнения перехода с п-го слоя на (п + 1)-й слой по времени (п ^ 0, 0 < т < М): = Л-- ^Ь — ¿0]2Ь"2,

ст- — к+2]Ст+1+¿тй+^¿т = 0.

Значения в начальный момент времени известны: (Ст = и0 = 0 (0 < т < М). Следуя методу прогонки (алгоритм Томаса), ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое по времени в виде

СП1 = ат+1^т+1 + вт+1, т = 0,..., М — 1.

Прогоночные коэффициенты следующие: (m = 1,..., M — 1): ат+1 = [2 + W1 — ат]

-1

вт+1

вт + W1 <

2 + w1 — ап

(19)

Для нахождения начальных коэффициентов а1, в1 воспользуемся следующими соображениями. Подсчитаем предварительно значения СП+1 по явной разностной схеме (в равенстве (17) справа заменяем п + 1 на п). Подставляя эти значения в (18) и используя второе условие (9) (£/п+1 = СП+1), в итоге получаем выражение Сп+1 = МС^1, Сп+1) (положительный корень квадратного уравнения):

¿1(СП+1)2 + + А 2 = 0,

А = Л-1 + 3 ]-1, А2 = —+

+ d2[<<r+1 — 4<<r+1][2hz ]-1 — U r[hf ]

1

Зная численное значение (7Г+1 и выражение ((Г+1 = а1<('Г+1 + в1, получаем а1 = 0, в1 = ¿Г+1. По а1, в1 вычисляем оставшиеся прого-ночные коэффициенты ат, вт, m = 2,..., M.

Ближайшая цель — найти значение С1^1, необходимое для реализации обратного метода прогонки. Рассмотрим первое краевое условие непроницаемости при z = 1 в (9). Используя разностную аппроксимацию

dzc|z=1 « [<<М+-12 — 4<^М+-11 + 3<Cjr+1]/2hz = 0

и соотношение <<т+1 = ат+1<<т+:1 + вт+1 при M — 2, M — 1, находим выражение <<Г+1 =

[(4 — ам-1)вм — вм-1][(ам-1 — 4)ам + 3]-1.

Следующий этап: с текущими приближениями значений <<Г+1, <5]Г+1 решаем обратным ходом прогонки трехдиагональную систему линейных уравнений и находим новые приближения концентраций <1+1 (и остальные значения 0ПП+1 для m = 3,..., M — 1). Далее снова пользуемся формулой <<Г+1 = /1(<<г+1, <<2!'+1). После этого корректируем значения прогоноч-ных коэффициентов (19), определяем <5Г+1 =

[(4 — ам-1)вм — вм-1][(ам-1 — 4)ам + 3]-1 и повторяем вычисления (возвращаясь к началу абзаца) до установления граничных значений <<ом (обычно 2-3 итерации). Критерием окончания вычислений на этапе I выбрано условие Ur+1 ^ Q/u = 0,98.

Результат окончания I этапа: образовался гидридный слой толщины ¿0 с концентрацией химически связанного водорода Qh, а в слое раствора толщины (L — ¿0) имеется распределение концентрации растворенного водорода <p(x), ^>(¿0) = Q.

Этап II: гидридообразование и движение границы фазового перехода

Начальное распределение ^(z) (<^(x)) определяется первым этапом. Формально можно забыть о физическом смысле функции ¿(t') как границы раздела фаз и рассматривать ее как функциональный параметр. По решению ¿(t') (зная A(t')) определяются коэффициенты модели. Поэтому вычислительный алгоритм основан на неявных разностных схемах и носит итерационный характер. Итерации будут связаны с уточнением значения A(t') на каждом временном слое.

Введем следующие сетки: {ут = mhy, m = 0,1,...,M}, hy = 1/M - шаг по пространственной переменной y; {z^ = khz, k = 0,1,..., K}, hz = 1/K - шаг по пространственной переменной z; {t1 = nhf, n = 0, 1,...},

18

Н^ — шаг по времени Обозначим через (Кп), (СП) приближенные значения концентраций в гидридном ¿(*/)-слое ("С(*П,Ут)) и в (Ь — ¿(*/))-слое раствора (с(*П, Zfc)).

1. Начальный этап: нахождение Л > 0

При ¿ = ¿о будем рассматривать последовательно и независимо две краевые задачи: сначала для гидридного слоя, потом для слоя раствора в сплаве 1Ж&. Граница раздела фаз на данном этапе неподвижна, поэтому полагаем в (12), (14) Л = 0^ = 0). Опишем алгоритм перехода с п-го на (п + 1)-й слой по

1.1 Гидридный слой

Рассмотрим краевую задачу (12)-(13) при

Л = 0, ¿ = ¿о:

dv

'D* / L л2 д2с

dt7 = V D" ■ UJ ■ , уе (0,^

v(0, у) = 0, у е [0,1], V(t7, 1) = 0,

(20)

dv

d* f1 - г,2^ 0)] = --

_ M0V

' d* = -n-. (21) y=0 D*

V^1 _ Vn v m v m

(22)

l л vm+1 - 2im+i + vm+1

m+1

am+1 =

1

2 + W2 - a

вт+1

вт + W2 Vm

2 + W2 - a„

(23)

(п + 1)-м слое по £ аппроксимируем производную дуV |у=о - [—3)га+1 + 4У1га+1 — )га+1]/2Ну и подставляем в граничное условие (21):

d* [1 - (VT+1 )2] = -

-3V0n+1 + 4V1ra+1 - V2n+1

2hy '

В итоге получаем )га+1 = /2()1га+1, )га+1) (положительный корень квадратного уравнения):

(V0n+1)2 + B1V0n+1 + B2 = 0,

n+1

1

Для уравнения диффузии (20) рассмотрим неявную разностную схему на (п + 1)-м слое: т = 1,...,М — 1,

Рассмотрим уравнения перехода с n-го на (n + 1)-й слой по t (0 < m < M):

vm+1 - [w2+2] vm+1+vm+1+w2Vn = 0,

w2 = [,^/L]2^D/D* h^/ht. Начальные данные

известны: 1/° = 0 (0 ^ m ^ M). Следуя методу прогонки, ищем приближенные значения концентрации в узлах сетки в виде Vm+1 = am+1Vm+1 + вт+1, m = 0,...,M-1. Прого-ночные коэффициенты (m = 1,.. .,M-1):

В = 3[2Ну й*] 1, й* = МосД В = [)га+1 — 4)га+1][2Ну й*]-1 — 1.

Зная численное значение )га+1 и выражение для Т)о"'+1 = а1У1га+1 + в1, получаем а = 0, в1 = !)ога+1. По «1, в1 вычисляем оставшиеся прогоночные коэффициенты ат, вт при т = 2,..., М по формулам (23).

Ближайшая цель — найти значение )П+1, необходимое для реализации прогонки. Из краевого условия г)(*/, 1) = 0 ^ )М+1 = 0.

Следующий этап: с текущими приближениями значений КП+т1 решаем обратным ходом прогонки трехдиагональную систему линейных алгебраических уравнений и находим новые приближения концентраций )г+1

(и остальные К^1, т = 3,..., М — 1).

После этого возвращаемся к уточнению

)га+1 = /2(КГ+1,)Г+1) в силу (21), корректируем прогоночные коэффициенты (23) и повторяем вычисления (возвращаясь к предыдущему абзацу) до установления распределения с(*га+1,У) [Кст+1, т = 0,..., М]. Шаг по времени мал, достаточно нескольких итераций.

1.2 Слой раствора —

Переходим к краевой задаче (14)-(15), возвращаясь к предыдущему слою Л = 0,

дс dt7

L

L - 4

2.

dz2'

(24)

О дс

с(0^) = с^, 0) = ^ =0. (25)

и дz ^=1

Для уравнения диффузии (24) рассмотрим неявную схему на (п+1)-м слое: к = 1,.. .,К—1,

Для нахождения начальных коэффициентов «1, в1 воспользуемся следующими соображениями. Подсчитаем предварительно значения )г+1 по явной разностной схеме (в равенстве (22) справа заменяем п + 1 на п). На

Cn+1 _ Cn

(26)

L

L - 4

Cn+1_2C n+1 + Cra+1

k—1 k

h2

2

Уравнения перехода на (п + 1)-й слой

(0< к <К): = [(Ь — ¿0)/Ь]2^Б*/БЛ2/Л^,

С—1 — И + 2]Сп+1 + сп+1 + ^4 сп = 0.

Начальные данные: С0 = ) (0 ^ к ^ К). Следуя методу прогонки, ищем приближенные значения концентрации в узлах на (п + 1)-м

слое по в виде СП+1 = а^Сп^1 + вк+1, к = 0,..., К — 1. Прогоночные коэффициенты (к = 1,...,К — 1):

ак+1

А+1

1

2 + -ш4 — ак

вк: + ^4 Сп

2 + -ш4 — ак'

Зная численное значение С^1 = Q/u (первое граничное условие в (25)) и выражение С701+1 = а1Сп+1 + в1, находим начальные коэффициенты а1 =0, в1 = Сп+1. По а1, в1 вычисляем оставшиеся коэффициенты , вк, к = 2,..., К, по формулам (27).

Теперь найдем значение СК+1, необходимое для реализации прогонки. Рассмотрим краевое условие непроницаемости в (25). На (п + 1)-м слое по ¿' аппроксимируем производную

0,..., М], с(4+1, г) [Сп+1, к = 0,..., К], зная распределения с предыдущего слоя по времени г)(Ьп,у), ¿(¿п,г). Рассмотрим последовательно два слоя материала: сначала гидридный, затем слой сплава — с раствором Н (это соответствует и последовательности движения растворенного диффузионно подвижного водорода сквозь материал с течением времени).

Начальное приближение Л(^п) ~ Л(ьп_ 1) берется с предыдущего временного слоя. Параметр Л(4) подлежит итерационному уточнению, значение ¿(¿п) известно.

2.1 Гидридный слой

Для уравнения диффузии (12) в безразмерной форме рассмотрим неявную разностную (27) схему на (п + 1)-м слое: т = 1,..., М — 1,

тсп+1 _ тсп Т>п+1 _ сп+1

Тт Тт = Л(Ьп) Ут ■ т+\— т-1 + (28)

дс дг

2=1

СК+2 — 4СК+\ +3СК+1

2Лг

0

и подставляем вместо СК+2 к- 1 выражения Сп+1 = ак+1Сп+11 + вк+1 при соответствующих к. В итоге получаем СК+1 =

= [(4 — ак-1)вк — вк-1][(ак-1 — 4)ак + 3] 1.

Используя СсКп+1, решаем обратным ходом прогонки трехдиагональную систему линейных алгебраических уравнений и находим значения концентрации Сп+1, к = 1,..., К — 1.

1.3 Итерационное уточнение Л

После проведенных вычислений из уравнения (16) находим новое приближение Л. Если Л ^ 0, то возвращаемся в начало излагаемого алгоритма (пункт 1.1), в противном случае переходим к пункту 2 с текущими распределениями концентраций в качестве начальных данных и начальным значением Л > 0.

2. Переход от п-го к (п + 1)-му слою по ¿'

Опишем алгоритм нахождения распределений на (п + 1)-м слое ^(^п+1, У) [ст+1, т =

Л

2Л

у

+ /Б

+ V Б

Ь

¿(4)

2 — 2сп+1 + сп+1

т+1

Рассмотрим уравнения перехода с п-го слоя на (п + 1)-й слой по Ь' (п ^ 0, 0 < т < М):

[1 — — К +2]сп+1 +

+ [1 + ^з]С+1 + = 0,

, / ч ¡Б

^2(Ьп) = \! Б

¿(4) I2

Ь

Л/'

^(С т) ^ ) ■ ^ ■ Л(Ьп) Ут

2

Ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое в виде Сп+1 = ат+1^п+1 + вт+1, т = 0,..., М—1. Прогоночные коэффициенты (т = 1, 2,..., М — 1):

1 + ^з

а +1 =

2 + ^2 + (^з — 1)ат'

(1 — ™з)в + ^2 сп

Рт+1 = —-—7-ТТ-. (29)

2 + "Ш2 + (^з — 1)ат

Для нахождения а1, в1 подсчитаем предварительно значения с1г+1 по явной разностной схеме (в равенстве (28) справа заменяем п + 1 на п). На (п + 1)-м слое по ¿' аппроксимируем

ду£|у=0 - [—3СТ+1 + 4Сп+1 — с2п+1]/2Лу и подставляем в граничное условие (13):

¿*[1 — (СГ+1 )2] = —

—3С0п+1 + 4С1п+1 — с2п+1

2ЛУ '

20

В итоге получаем У0га+1 = /2 (У1га+1, У2га+1) (положительный корень квадратного уравнения):

СЮ)2 + £1^0 + В2 = 0, £ = 3[2ДУ 4]

1

1

£2 = [Ю - 4У1][2Лу 4]-1 - 1, 4 = МгШ

V = у^1, £ = £(^). Зная значение У0га+1 и выражение = а^^1 + въ получа-

ем а1 =0, в1 = ^0га+1. По «1, в1 вычисляем оставшиеся прогоночные коэффициенты ат, вт, т = 2,..., М, по формулам (29).

Теперь найдем значение УМ+1, необходимое для реализации прогонки. Из краевого условия £(£', 1) = 0 получаем ТМ+1 = 0. Следующий этап: с текущими приближениями УП+г1 решаем обратным ходом прогонки трехдиаго-нальную систему линейных уравнений и находим новые приближения концентраций Т/^1

(и остальные т = 3,..., М — 1).

После этого возвращаемся к уточнению У0га+1 = /2(УГ+1, УТ+1) в силу уравнения (13), корректируем значения прогоночных коэффициентов (29) и повторяем вычисления до установления распределения ^(¿П+1,у) [ТП+1, т = 0,..., М]. Шаг по времени мал, так что достаточно нескольких итераций.

2.2 Слой раствора —

Переходим к краевой задаче (14)-(15), возвращаясь к предыдущему слою (¿' = ¿П).

Для уравнения диффузии (14) в безразмерной форме рассмотрим неявную разностную схему на (п + 1)-м слое: к = 1,..., К — 1,

СП+1 — Сп А(4) [1 — ^] — С71П+11

Л*

В

+ V В

— 1

2Л

+

Ь

1_Ь — £(4Л

2 _ 2Сп+1 + С™+1

к+1

Л2

(30)

Рассмотрим уравнения перехода с п-го слоя по на (п + 1)-й слой (п ^ 0, 0 < к < К):

[1 — ^С—1 — и + 2]Сп+1+

+ [1 + + ^С? = 0,

в

™4(4) = \/ "ТТ

ь—£(4)

Ь

Л2

л*/

^(4; к) = ^4 (4)

£(4) Л* Ь — £(4)

Ищем приближенные значения концентрации в узлах сетки на (п + 1)-м слое в виде

СП+1 = а^СП^1 + вк+1, к = 0,

К 1.

Прогоночные коэффициенты (к = 1,..., К—1):

ак+1

А+1 =

1 + ^5

2 + -Ш4 + (^5 — 1)ак;

(1 — + ^4 Сп

2 + -Ш4 + (^5 — 1)ак'

(31)

Зная численное значение С0 = ф/м (первое из граничных условий (15)) и выражение Сц = + в1, находим начальные коэффи-

циенты а1 = 0, в1 = С. По а1, в1 вычисляем оставшиеся прогоночные коэффициенты ак, , к = 2,..., К, по формулам (31).

Найдем значение СП+1, необходимое для прогонки. Рассмотрим краевое условие непроницаемости с|г=1 = 0 в (15). На (п + 1)-м слое по аппроксимируем производную

^с|г=1 » [СП— — 4СП+11 + 3СП+1]/2Лг = 0

и подставляем вместо С^— к- 1 выражения

С-1 = «к СП+1 + вк

при соответствующих к. В итоге СП+1 = = [(4 — ак-1)вк — вк-1][(ак-1 — 4)ак + 3]-1.

Используя СП+1, решаем обратным ходом прогонки трехдиагональную систему линейных алгебраических уравнений и находим значения концентрации СП+1, к = 1,..., К — 1.

2.3 Итерационное уточнение А(4)

После указанных вычислений аппроксимируем производные на (п + 1)-м слое

ду%=1 « [ТМ- — 4уП+11 + 3УМ+1]/2Лу, ^с|г=0 « [—3С01+1 + 4СП+1 — СП+1]/2Лг

и подставляем их в граничное условие (16):

А(4) =

ь м 0

'В* V

~о М

Ь

й

2 Т/га+1 4Т/п+1 + ЧТ/га+1 2 ТМ- 2 — 4ТМ- 1 + 3ТМ

2Л

+

у

В Ь2 —3Сп+1 +4Сп+1 — Сп+1

В* £[Ь — £]

2Л2

2

2

Отсюда находим новое приближение Л(4) > 0 и возвращаемся в начало излагаемого алгоритма (пункт 2.1). Через несколько итераций получаем установившееся значение Л (¿п). В итоге считаем, что число Л(^) фиксировано.

Далее вычисляем ¿(4) = Л(4Ж4) и в силу в

Л(4) = ^1п ¿(¿')1*п ^ ^ ln(¿(tn + ль')) - 1п ¿(4) + Л(4) ■ Д*

определяем 1п¿(*п+1), 4+1 = 4 + Д*'. По значению 1п ¿(¿п+1) находим ¿(¿п+1) и переходим к слою 4+2, приняв в качестве начального приближения Л(*п+1) — Л(*п). Далее уточняем значение Л(*п+1) и распределения й(4+2,у), с(*п+2,г) по изложенной схеме (пп. 2.1-2.3).

Тестирование алгоритма

При тестировании алгоритма использовались следующие входные данные: Ь = 6 х 10-2 см, ¿0 = 1,3 х 10-:з см, р = 1520 торр, Т = 593 К, Q = 6,14 х 1020 см-3, Qh = 6, 5 х 1021 см-:з, Б0 = 2, 2 х 10-3 см2 с-1, = 35 х 10з Дж/моль, = 1, 5 х 10-:з см2 с-1, Ед = 59 х 10з Дж/моль, Ь = 5, 5 х 10-24 см4 с-1, Ь* = 3 х 10-27см4 с-1, 8 = 7 х 10-7, =6 х 10-8. На рисунке представлена возможность анализировать динамику распределения концентраций растворенного водорода в слоях гидрида и сплава с помощью программного обеспечения, разработанного в среде ВеПаЪ.

Заключение

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

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

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

Финансовое обеспечение исследований осуществлялось из средств федерального бюджета на выполнение государственного задания КарНЦ РАН (Институт прикладных математических исследований КарНЦ РАН).

Распределения концентраций растворенного водорода в слоях гидрида и сплава Distribution of dissolved hydrogen concentration across the hydride and the alloy layers

Литература

1. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.

2. Водород в металлах / Ред. Г. Алефельд, И. Фелькль. М.: Мир, 1981. Т. 1, 506 с.; т. 2., 403 с.

3. Изотопы водорода. Фундаментальные и прикладные исследования / Ред. А. А. Юхим-чук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697 с.

4. Писарев А. А., Цветков И. В., Марен-ков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.

5. Черданцев Ю. П., Чернов И. П., Тюрин Ю. И. Методы исследования систем металл-водород. Томск: ТПУ, 2008. 286 с.

6. 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

7. Evard E.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

8. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010.

9. Lototskyy M. V., Yartys V. A., Pollet B. G, Bowman R. C. Jr. Metal hydride hydrogen

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. Izotopy vodoroda. Fundamental'nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies]. Sarov: RFYaTs-VNIIEF, 2009. 697 p.

4. 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.

5. 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.

6. Denisov E. A., Kompaniets M.V., Kompaniets T. N., Bobkova I. S. Peculiarities of hydrogen permeation through Zr-1%Nb alloy and

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. The hydrogen economy/ Eds. M. Ball, M. Wietschel. Cambridge University Press, 2009.

12. 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

13. 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

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. 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

Поступила в редакцию 17.04.2018

evaluation of terminal solid solubility. Journal of Nuclear Materials. 2016. Vol. 472. P. 13-19. doi: 10.1016/j.jnucmat.2016.01.022

7. Evard E. 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

8. Handbook of hydrogen storage: new materials for future energy storage. 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. The hydrogen economy. Cambridge Univ. Press, 2009. 646 p.

12. 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

13. 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

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. 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

Received April 17, 2018

СВЕДЕНИЯ ОБ АВТОРАХ:

Заика Юрий Васильевич

рук. лаб. моделирования природно-технических систем, д. ф.-м. н. Институт прикладных математических исследований КарНЦ РАН, Федеральный исследовательский центр «Карельский научный центр РАН» ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: zaika@krc.karelia.ru тел.: (8142) 780059

Родченкова Наталья Ивановна

старший научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований КарНЦ РАН, Федеральный исследовательский центр «Карельский научный центр РАН» ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: nirodchenkova@krc.karelia.ru тел.: (8142) 766312

Грудова Ксения Васильевна

стажер-исследователь Институт прикладных математических исследований КарНЦ РАН, Федеральный исследовательский центр «Карельский научный центр РАН» ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: grudova@krc.karelia.ru тел.: (8142) 766312

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) 780059

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

Grudova, Kseniia

Institute of Applied Mathematical Research,

Karelian Research Centre,

Russian Academy of Sciences

11 Pushkinskaya St., 185910 Petrozavodsk,

Karelia, Russia

e-mail: grudova@krc.karelia.ru

tel.: (8142) 766312