CC BY
24
Труды Карельского научного центра РАН
№8. 2017. С. 31-47 DOI: 10.17076/mat646
УДК 519.6:539.2
ОБРАТНАЯ ЗАДАЧА ИДЕНТИФИКАЦИИ СПЕКТРОВ ТЕРМОДЕСОРБЦИИ ВОДОРОДА
Ю. В. Заика, Е. К. Костикова
Институт прикладных математических исследований Карельского научного центра РАН, Петрозаводск
В рамках технологических задач водородного материаловедения (включая проект ITER) ведется интенсивный поиск различных по назначению конструкционных материалов с заданными пределами водородопроницаемости. Одним из экспериментальных методов является термодесорбционная спектрометрия (ТДС). Образец, насыщенный водородом, дегазируется в условиях вакуумиро-вания и монотонного нагрева. С помощью масс-спектрометра регистрируется десорбционный поток, позволяющий судить о характере взаимодействия изотопов водорода с твердым телом. В статье представлены распределенная краевая задача термодесорбции с динамическими граничными условиями и численный метод моделирования ТДС-спектра, требующий лишь интегрирования нелинейной системы обыкновенных дифференциальных уравнений (ОДУ) невысокого порядка (по сравнению, например, с методом прямых). Предложен метод оценки коэффициентов растворения и десорбции в условиях реального динамического взаимодействия поверхностных и диффузионных процессов, без искусственного разделения исследований на режимы лимитирования диффузией (DLR) и поверхностных реакций (SLR).
Ключевые слова: водородопроницаемость; поверхностные процессы; термодесорбция; динамические граничные условия; параметрическая идентификация; численное моделирование.
Yu. V. Zaika, E. K. Kostikova. INVERSE PROBLEM OF IDENTIFICATION OF HYDROGEN THERMODESORPTION SPECTRA
One of the technological challenges for hydrogen materials science (including the ITER project) is the currently active search for structural materials with various potential applications that will have predetermined limits of hydrogen permeability. One of the experimental methods is thermal desorption spectrometry (TDS). A hydrogen-saturated sample is degassed under vacuum and monotone heating. The desorption flux is measured by mass spectrometer to determine the character of interactions of hydrogen isotopes with the solid. The paper presents a distributed boundary-value problem of thermal desorption with dynamical boundary conditions and a numerical method for TDS spectrum simulation, where only integration of a nonlinear system of low order (compared with, e. g., the method of lines) ordinary differential equations (ODE) is required. A method for estimating dissolution and desorption coefficients in a real dynamic interaction of surface and diffusion processes without the artificial division of studies into diffusion limited regime (DLR) and surface limited regime (SLR) is suggested. This work was supported by the Russian Foundation for Basic Research (project № 15-01-00744).
Keywords: hydrogen permeability; surface processes; thermal desorption; dynamical boundary conditions; parameter identification; numerical simulation.
Введение
Интерес к взаимодействию изотопов водорода с различными материалами носит многоплановый характер [1-4, 11, 12, 15-17, 23, 24]. Достаточно упомянуть перспективы развития энергетики, проблему защиты конструкционных материалов от водородной коррозии, задачи проектирования химических реакторов, ракетостроения. Некоторые аспекты водородного материаловедения, связанные с тематикой данной статьи, исследованы в [6-9, 13, 21, 26, 30]. Работа является непосредственным продолжением [19, 27-29].
Экспериментальный опыт показывает, что лимитирующими являются не только диффузия, но и физико-химические процессы на поверхности [1, 2]. Параметры переноса зависят и от технологических особенностей получения партии материала, поэтому вряд ли следует ориентироваться на получение «табличных значений», нужны эффективные алгоритмы обработки экспериментальных данных. В статье остановимся на моделировании термодесорбции, учитывая лишь основные лимитирующие факторы и информативность рассматриваемого варианта ТДС-эксперимента.
Работы [6, 13, 21, 27] посвящены проблеме интерпретации ТДС-пиков. Анализ причин различных всплесков термодесорбции принципиально важен при выборе реакторных конструкционных материалов, контактирующих с изотопами водорода. В [6, 13] представлен достаточно подробный обзор. При численном моделировании без учета динамики поверхностных процессов ТДС-пики вынужденно интерпретируются как следствие захвата диффундирующего атомарного водорода дефектами структуры материала (ловушками) с различными энергиями связи и(или) как проявление многоканальной диффузии [5]. Разнообразие причин этим не исчерпывается. В [27] показано, что к появлению пика десорбции могут приводить и определенные сочетания скоростей поверхностных и диффузионных процессов. Таким образом, проблема интерпретации ТДС-спектра усложняется.
Данная работа, являясь продолжением [27], использует ту же математическую модель. Поэтому авторы будут стремиться к более сжатому изложению вопросов, подробнее освещенных в [27]. Основным результатом настоящей статьи является метод параметрической идентификации по экспериментальным данным. Трудности решения обратных задач математической физики хорошо известны, име-
ется обширная математическая литература и ряд специализированных журналов (inverse problems, ill-posed problems ...). В экспериментальной практике, чтобы избежать вычислительной некорректности многопараметрических обратных задач, стараются свести оценку параметров к раздельному исследованию DLR (diffusion limited regime) и SLR (surface limited regime). Но реальный материал эксплуатируется в условиях динамического взаимодействия «поверхность-объем». Для оценки коэффициента диффузии имеются подробно разработанные методики. Значительно сложнее определить параметры десорбции и растворения (если искусственно не занижать температуру, чтобы «отключить» объемные процессы). В статье представлен алгоритм, который позволяет оценить коэффициенты десорбции и растворения, когда диффузия активно взаимодействует с поверхностными процессами. Авторы стремились к использованию минимального математического аппарата, чтобы алгоритм был доступен для реализации на уровне использования свободно распространяемых математических пакетов (без необходимости разработки специализированного программного обеспечения).
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПЕРЕНОСА
Уравнение диффузии
Рассмотрим диффузию водорода в образце тестируемого материала (пластине толщиной Для определенности подразумеваем металл, сплав. Материал однороден, чтобы пренебречь взаимодействием атомов H с ловушками. Нагрев достаточно медленный. Для диффузии примем стандартное уравнение:
ct(t, x) = D(T)cxx(t, x), (t, x) e Qt., (1)
где Qt„ = (0,t*) x (0,^); c(t,x) — концентрация диффундирующего H. Зависимость коэффициента диффузии D от температуры T = T(t) соответствует закону Аррениуса: D(T) = Do exp{-ed/[RT ]}.
Можно учесть несколько каналов диффузии (транскристаллическую, по границам зерен, вдоль дефектов) с взаимообменом между ними. Однако значительное увеличение числа неизвестных априори параметров делает задачу их оценки труднообозримой. Информативность ТДС-эксперимента ограничена, так что можно считать D в уравнении (1) некоторым интегральным эффективным показателем.
«В погоне» за описанием множества пиков термодесорбции удобна модель
дс д2с
dt = дХ2 -
(2)
V = 1
- У^ U- [1 - Zv] c(t, x) - a+Zv(t, x)
dz,
= а-(Т)[1 - с(*,ж) - а+(Т^ (*,ж),
где zv(*,х) — концентрации атомов водорода, захваченного дефектами различных типов; а^ — коэффициенты поглощения и выделения Н ловушками; Zv = zv(*, ж)^тах — степень насыщения ^тах = тах zv). Для практических целей захват учтен в простейшей «интегральной» форме, уточнение геометрии дефектов и их распределения существенно усложнило бы модель. Если дефект, например, не микрополости, а включения гидридной фазы, то на этапе ТДС-дегазации коэффициент а-(Т) тождественно равен нулю, а значение а+(Т) положительно лишь после достижения критической температуры: Т(*) ^ ТСГй.
За счет различных энергий связи (соответствующих коэффициентов Еа) можно добиться заданного количества ТДС-пиков. При этом такая модель требует специализированного программного обеспечения. Алгоритмы на основе разностных схем и результаты численного моделирования представлены в работах [19, 28]. В данной статье ограничимся уравнением (1) (Б = ).
Динамические граничные условия
Особенность модели связана с динамическими нелинейными граничными условиями. Пусть пластина контактирует с газообразным Н2 и поверхность является существенным потенциальным барьером (см. [1, с. 177-206; Га-бис, Компаниец, Курдюмов]). Тогда с учетом (де)сорбционных процессов краевые условия моделируются следующим образом:
с(0,х) = с(ж), х е [0,4* е [0,**], (3)
сс(*) = д(Т)до(), ф) = д(ТЫ*), (4)
<?с(*) = ^(Т)в(Т)рс(*)-Ь(ТШг)+Всх(1, 0), (5) де(*) = ^(Т)в(Т)р(*)-Ь(Т)д2(*)- Бсх(*,¿), (6) Ь(Т )= Ьс ехр{-Еь[КТ ]-1}, з(Т) = вс ехр{-Е3[КТ]-1}.
Здесь: сс(*) = с(*, 0), с^(*) = с(*, £) — граничные объемные концентрации Н (знак
тождества часто используем в смысле равенства по определению); qo/(t) — концентрации на поверхностях (x = 0,^, q = dq/dt); g(T) — параметр локального равновесия между концентрациями на поверхности и в приповерхностном объеме (коэффициент быстрого растворения); ß — газокинетический коэффициент; po/(t) — давления молекулярного газа (H2); b(T) — коэффициент десорбции. Математическая корректность таких динамических граничных условий доказана в [31].
В модели фигурирует как молекулярный, так и атомарный водород. Для единообразия подсчет будем вести в атомах: [с] = 1/см3, [q] = 1/см2, [Dcx] = [bq2] = 1/см2с. В кинетической теории газов величина ßp определяет число частиц (в данном случае H2), соударяющихся с единичной площадкой поверхности в единицу времени. Но за счет множителя s удобно воспринимать слагаемое ßsp как плотность потока атомов, оседающих на поверхности (без разделения процесса на стадии). Вместо s можно написать 2s и интерпретировать s как долю адсорбируемых атомов H.
В масштабе рассматриваемого далее ТДС-эксперимента удобно выбрать [p] = mopp. Тогда получаем выражение ß(T) = (2nmkT)-1/2 и 2.474 х 1022/уГ. Здесь [ß] = 1h2/(mopp см2с), [T] = K (под корнем численное значение температуры), k — постоянная Больцмана, m — масса молекулы водорода. Зависимостью кинетической константы от температуры (ß а 1/VT) часто пренебрегают на фоне экспоненты в зависимости s(T).
Баланс потоков «поверхность—объем»
Модель (4) быстрого растворения (локального равновесия) на поверхности получается из более общих соотношений
k~(T)[1 - co^cmL] qo,e(t)--k+(T)[1 - qoAt)q-L]coAt) = TD(T)c
(7)
\o,£'
Коэффициенты к-, к+ характеризуют интенсивность растворения в объеме и выхода на поверхность. Если эти процессы в рассматриваемом диапазоне температур существенно быстрее диффузии и концентрации далеки от максимальных, то, полагая в относительном масштабе Бсх ~ 0, получаем соотношение (4) с д = к-/к+. Если поверхность изотропна (в смысле Ек- ~ Ек+), то д(Т) слабо зависит от температуры. Формально можно записать д(Т) = дс ехр{-Ед/[КТ}], Ед = Ек- - Ек+, но «энергия активации» Ед не обязательно положительна. Наличие «порогового» множителя [1 - с/стах] приводит к замедлению рас-
творения при во,е ^ cmax. Аналогично интерпретируется другой множитель, где величина в{Ь) = q(t)/qmax означает степень заполнения поверхности. В балансовых уравнениях (5), (6) можно моделировать плотность потока адсорбции атомов H выражением ^sp[1 — в]2. Но при ТДС-дегазации в диапазоне достаточно высоких «рабочих» температур в ^ 1.
Эти упрощения согласуются с ограниченной информативностью ТДС-эксперимента. При большом количестве параметров легче аппроксимировать экспериментальные данные. Но теряется единственность оценок. Это может привести к существенным ошибкам при экстраполяции результатов «с тонкой пластины на стенку», что неявно предполагается.
ТДС-эксперимент
Уточним рассматриваемый вариант ТДС-эксперимента. В камеру c пластиной из исследуемого металла или сплава подается водород в газовой фазе при сравнительно большом давлении. Образец нагрет с целью увеличения скорости сорбции. После равновесного насыщения он быстро охлаждается (отключается ток нагрева). В режиме последующего постоянного вакуумирования камеры образец медленно нагревается. Нагрев обычно линейный: T(t) = To + vt. Скорость нагрева v невелика. По достижении максимальной температуры (если дегазация еще не завершилась) нагрев прекращается: T(t) = Tmax. С помощью масс-спектрометра измеряется давление p(t) молекулярного водорода в вакуумной камере, обусловленное десорбцией J(t) = Jdes(t) = b(t)q2(t). Здесь и далее для упрощения обозначений примем сокращенную запись
b(t) = b(T(t)), D(t) = D(T(t)), s(t) = s(T(t)).
Обобщенный коэффициент g быстрого растворения (локального равновесия «поверхность-объем») считаем константой (Ek- = Ek+) в пределах диапазона нагрева. Для метода ТДС выполнена симметрия:
p(t) = po/(t), q(t) = qo/(t), co(t) = ce(t), (8)
D(t)cx(t, 0) = —D(t)cx(t,i), c(x) = с = const.
Время t* окончания эксперимента определим дегазацией: p(t) & 0, t ^ t*, c(t*,x) & 0, x £ [0,4 Информативность начальной и конечной стадии эксперимента невелика, так что достаточно ограничиться t*, когда поток из образца упадет на порядок по сравнению с максимумом. Экспериментальными данными считаются графики плотности десорбции J(t)
или ТДС-спектры J(T) (T(t) о t) при различных условиях насыщения и скоростях нагрева. «Пересчет» р(-) ^ J(■) конкретизируется характеристиками экспериментальной установки. Современное оборудование позволяет достигать достаточно глубокого вакуума (10-9 - 10-7 mopp). Поэтому на этапе насыщения слагаемое P = ßsp (P ^ 1) является определяющим, а на этапе дегазации ресорб-цией можно пренебречь (P ^ 1).
Подчеркнем, что описан один из вариантов насыщения. Используются также криогенные технологии и электрохимический способ. Целесообразно процесс насыщения согласовывать с условиями эксплуатации конструкционного материала. Соответствующие технические корректировки далее непринципиальны.
Анализ равновесия
Кратко проанализируем равновесные соотношения в рамках принятого уровня детализации моделирования. Пусть давление и температура постоянны. Формально равновесие характеризуется равенством нулю всех производных. Ориентируясь на широко применяемый закон Сивертса, будем интересоваться пропорциональностью величине л/р.
1. На поверхности выполнено соотношение
q
^вр[1 — в{\2 — bq2 = 0, 6>1 =
5шах
Здесь в1 — степень заполнения поверхности (далее у в^ аналогичный смысл). Параметр в можно заменить (переобозначить) на 2в, если интерпретировать в как долю адсорбируемых атомов из газовой фазы. Выразим в1:
2 _ Vs
V =
bq
2
max
^ pv2[i - öi]2 = 02 ^
=
v Vp
1 + V -p
< 1, V —p < 1 ^ 01 x -p.
2. В равновесии «поверхность-объем»
(02 = c/cmax, k1 = k qmax, k2 = k+cmax,
Y = k1/k2) из соотношений (7) получаем
kl0l[1 - 02] = k202[1 - 01] ^ 02 =
Y0i
1 + [y - 1]01'
3. Учтем емкость ловушек (ограничимся одним типом, т = 1). Воспользуемся уравнением (2), добавив для симметрии множитель насыщения [1 — в2\ при г = . Тогда (вз = г/гшах, а = а+гшах/[а-Сшах\) получим
02[1 - 03] = 0з[1 - 02]а ^ 03 =
02
а + [1 - а] 02'
Эксперимент «насыщение-дегазация» дает информацию о некоторой средней концентрации с в объеме V = Б£ (торцами пренебрегаем):
Ус = 2б>1<5тах£ + 02 стах V + 0зZmaxV. Нормируем с на стах и рассмотрим функцию
в(р) = = 2^ + 02 + 03 —. (9)
стах ^стах стах
Нас интересуют графики в осях (^р, в). Численные результаты показаны на рисунках 1, 2 (а = а+/а-, стах = 1018Н/см3, I = 0.1 см). Параметров влияния достаточно много, но обычно можно оценить порядки величин и понять, каково «расстояние до прямолинейности» (до области адекватности закона Сиверт-са).
Кривые в = /(-^р) имеют характерный вид «нарастающей волны». Это заметно лишь в широком диапазоне давлений. Для указанных параметров наиболее ярко точка перегиба
видна на кривых, маркированных звездочками. Анализ слагаемых в формуле (9) для общей концентрации показывает, что точка перегиба определяется выходом первого и третьего слагаемого на режим насыщения с дальнейшим преобладанием роста 02. Подчеркнем, что в более узких диапазонах давлений графики практически прямолинейны, что соответствует участкам классического закона Сивертса.
Выскажем некоторое предостережение. В экспериментальной практике при оценке во-дородопроницаемости материалов часто для плотности проникающего потока Н используется приближение Ричардсона
= РоиКРт. (10)
Мембраны обычно очень тонкие, с относительно высокой проницаемостью. Для градиента концентрации можно принять равенство сх = [соИ, -с;п]/^. Тогда в рамках закона Фика «точной» является формула 3 = -Б[со^ - сп]/£ Здесь 3 (*) — плотность проникающего потока атомарного водорода. Но непосредственно граничные объемные концентрации не измеряются, и приходится пользоваться квазиравновесным приближением, подставляя вместо с равновесную концентрацию с = Г^р в соответствии с законом Сивертса (Г — коэффициент растворимости, кратко - растворимость). В силу неравенств с;п > сп, со^ < со^ такая подстановка завышает второй множитель в выражении 3 = Б[сп - сои^/£: с;п - со^ > сп - coиt. Чтобы сохранить равенство, нужно формально занизить значение Б.
Итак, если по экспериментальным данным {р(*),3(*)} найдено подгонкой значение Ф (проницаемость) в соответствии с формулой (10), то Ф < БГ. По заданным числам Ф и Г из формулы Ф = БГ (которую обычно используют) получаем нижнюю оценку коэффициента диффузии Б. Но это еще не все. В режиме стационарной (квазистационарной) проницаемости участвует в основном растворенный диффузионно подвижный атомарный водород. В рамках принятой модели (4)-(6) имеем с = Г^р, Г = д л/^в/Ь. Эксперимент «насыщение-дегазация» дает значение общей концентрации с > с и завышенное (для задачи проницаемости) значение Гтах: с = Гтах^р.
Таким образом, можно подгонкой оценить число Ф (проницаемость) по формуле (10). Эта информация ценна на практике как удобный коэффициент пересчета давлений в поток. Если же взять значение Б из одного эксперимента (статьи), Г — из другого источника (Гтах), то получаем иерархию, строго говоря, трех разных чисел: Ф < БГ < БГтах. Для
материала с высоким уровнем захвата в ловушки можно на порядок завысить вычисленную проницаемость по сравнению с реальной. Проницаемость Ф (как параметр в соотношении (10)) имеет S-образную (аррениусовскую) форму кривой насыщения по порядку давлений. Лишь при относительно высоких давлениях (когда граничные концентрации приближаются к сивертсовским) имеем Ф ~ Dr. Необходимо контролировать, что подразумевается под одними и теми же терминами в разных публикациях и по какой модели каким методом тот или иной коэффициент вычислялся (коэффициенты не измеряются).
Функционально дифференциальное ТДС-уравнение
Идентификация ТДС-спектров необходима не только для выявления причин пиков термодесорбции, но и для численной экстраполяции и обобщения результатов, полученных на лабораторных образцах (обычно - — доли мм).
В модели (2) за счет ловушек с параметрами а± можно получить любое количество пиков. Могут ли различные пики возникать при дегазации практически однородного материала? Для выяснения этого вопроса ограничимся простейшим уравнением диффузии (1), но с динамическими симметричными граничными условиями (4)—(6), (8). Поверхность считаем изотропной в смысле д = const в рамках диапазона нагрева. Ресорбцией при ваку-умировании пренебрегаем. Итак, мы ограничились минимальным числом параметров для модели, которая учитывает динамическое взаимодействие поверхностных процессов с диффузией. Далее работаем в этом приближении.
Чтобы сравнивать модельные и экспериментальные ТДС-спектры с целью оценки параметров, нужна только поверхностная концентрация (далее J = Jdes = bq2). Естественно попытаться избежать итерационного решения краевой задачи при текущих приближениях D0, Ed, b0, Eb, s0, Es, д. Проведем преобразования, придя в итоге к интегрированию лишь системы ОДУ невысокого порядка.
Вывод уравнения типа Риккати
Принятая модель ТДС-дегазации:
ct = D(t)cxx , c(0, x) = с, co/{t) = gq(t),
q(t) = —b(t)q2(t) + D(t)cx(t, 0), J(t) = b(t)q2(t).
Сделаем замену времени t' = j^Ddr: ct(t,x) = cxx(t,x), c(0,x) = с, co,^ = gq(t), (11)
(оставляем прежнее обозначение Ь). Считаем q(Ь) функциональным параметром, а (12) — дополнительным соотношением к линейной задаче (11). Сделаем замену, приводящую краевые условия в (11) к однородным:
с = с(Ь, х) - gq(t), с*(Ь, х) = схх(г, х) + /(Ь),
/(Ь) = —gq(Ь), с(0,х) = ф(х) = 0, с|о/ = 0.
Запишем решение линейной краевой задачи с помощью функции мгновенного источника (функции Грина):
c(t,x)= / / Gi(x,i,t — т)f (т) didr, Jo Jo
n=1
, . . 2 ^ r n2n2 л nnx nn£ Gl(x, £ t) = exP{---T^ sm— sm—
В граничные условия входит cx(t, 0):
-I 4g f4 . . w r n2 n2 . Cx |0 = —у у0 Ч(т^ exP ( -p - (
n2 n2
т — tn dT,
Е' = Е„=1,3,5...' При т = Ь ряд расходится, так что подразумевается почленное интегрирование. В исходном времени Ь имеем Сх(Ь, 0) = Сх(Ь, 0) = вх(Ь,£)= ЙхМ),
4д '[t ( n2n2 f t *]
cx(t, 0) = — -f J <?(т) exp | — J D(s) dsj dr.
Динамическое граничное условие запишется в интегро-дифференциальной форме
q(t) = —b(t)q2(t) —
(13)
cx|o = — cxli = q(t) + J (t)D (t)
(12)
Г-^ Г п2п2 Г* . , ч ,
- — уо q(т) ехР\ —¿Г у °(з)
Полученное уравнение (13) с квадратичной нелинейностью будем классифицировать как функциональное уравнение Риккати нейтрального типа. Уравнение эквивалентно исходной краевой задаче в том смысле, что решение q(Ь) однозначно определяет решение с(Ь,х). Аналогия с ФДУ вида Х(Ь) = Т[Ь, х(Ь), Х(Ь), Х[о,*], Х[0,*]] нейтрального типа [10] связана с тем, что избавиться от производной q в правой части интегрированием по частям нельзя, так как появится расходящийся ряд. Нас интересует отрезок времени [Ь1,Ь2] С (0, Ь*), соответствующий пику термодесорбции (измерения при Ь ~ 0, Ь* малоинформативны). Уравнениям типа Риккати (включая матричные в теории оптимального управления) посвящена обширная литература.
t
Безразмерная форма задачи
Для моделирования удобно перейти к безразмерным переменным, выполнив замены:
¿' = )йтД2, х' = ж/^, V = с/Т (с = дТ).
Не меняя обозначения ¿, имеем:
v(t) = —b(i)v2(i) —
(14)
-
)(т) exp { — п2п2 [t — тj} dr
где v(0) = 1, b(t) = #(i)^2/D(i) - безразмерный параметр квадратичной «десорбции».
Этап начального насыщения
Остановимся на выбранном варианте (для определенности) начального насыщения пластины растворенным водородом и соответствующем уточнении коэффициента при квадратичной нелинейности b(t). Насыщение проводится при относительно высоких температуре T = const и давлении напуска p = const (для интенсификации сорбции). После установления равновесной концентрации имеем:
^(T)s(T)p = b(T)c2, с = #<?(# = const), ^ с = Texp {Eb[2RTj-1},
b(t) = с b(t)^2D-1 = с b(i)^2[gD(i)j-1, ^ b(t) = bo exp {—Eb[RT(t)j-1},
bo = Texp {Eb[2RTj-1}.
Формально Ь(£) = Ь(Т(¿)) — аррениусовский параметр с энергией активации Еь — Ед. Как правило, Еь > Ед и поверхностные процессы быстрее интенсифицируются с нагревом (Ь = ...ЕьТ, I) = ...ЕдТ), так что Еь > 0. Предэкспонента Ь0 зависит от всех парамет- ©(в) = 4^^ ехр { ров Ь0, Еь, в0, Е5 (кроме Ед). Тем самым коэффициент Ь в ТДС-эксперименте меняется строго монотонно и находится в ограниченном диапазоне 0 < Ь- ^ Ь(£) ^ Ь+.
В дальнейшем будем ориентироваться на предел поверхностной концентрации порядка Стах ~ 1015 — 1016 (монослой в рамках геометрической статики). Если на этапе насыщения
5, то при
модели динамики порог концентрации стах, по-видимому, может быть и выше, если уточнять, что понимать под «кипящим» поверхностным слоем (при достаточно высокой температуре). Фактическое значение Т несколько условно, поскольку на него сильное влияние оказывает предварительная подготовка эксперимента (вакуумирование перед началом нагрева). Но основная часть водорода находится в объеме, уравнение диффузии (как и сам процесс) обладает сглаживающим эффектом, и нас интересуют выраженные пики термодесорбции (начальный и конечный отрезки времени эксперимента малоинформативны). По этим причинам в начальных вычислениях достаточно правильно оценить порядок «эффективной» концентрации Т, требования к высокой точности здесь некритичны (и их практически не реализовать).
Подчеркнем, что обсуждался один из возможных способов начального равномерного насыщения. Можно придерживаться независимой методики. По окончании эксперимента (по зарегистрированной динамике десорбци-онного потока) подсчитывается общее содержание водорода в образце. В модели (по Б, I и текущей оценке д) определяется начальное (равномерное) насыщение и соответствующее значение С Наконец, можно использовать данные по растворимости (с = с(р, Т)), поскольку обычно материал исследуется комплексно.
Выделение интегрируемой особенности
Функциональное уравнение термодесорбции (14) имеет особенность, которая затрудняет его численное решение. Функция
22 —n п sj,
Е' - Е
n=1,3,5...
для выбранных параметров q > 1015
численном моделировании учитываем степень заполнения поверхности и вместо ^вр = ЬТ2 используем соотношение ^вр [1 — ТС—ах]2 = ЬТ2 для вычисления начального значения Т.
Вместе с тем такое априорное ограничение (из соображений геометрического расположения неподвижных «шариков» на плоской поверхности) представляется дискуссионным. В
имеет конечные значения при в > 0. Ряд быстро сходится при больших в, а при формальной подстановке в = 0 (когда переменная интегрирования т достигает верхнего предела ¿) получаем расходящийся ряд. «Спасает» почленное интегрирование. При этом слагаемые ряда имеют порядок 0(п-2) (п = 2г — 1, г ^ 1), что приводит к медленной сходимости. Поставим следующую задачу: аппроксимировать уравнение системой ОДУ невысокого порядка, чтобы можно было эффективно использовать стандартные пакеты прикладных программ (например, свободно распространяемый БсИаЪ). Для явного выделения интегрируемой особенности проведем преобразования, используя аппарат тэта-функций Якоби. Точ-
t
o
нее, нас интересует функция
те
93(t, x) = 1+2^exp{—n2n2t} cos(2nnx), t > 0.
n=1
При x = 0 имеем другое представление [25]: 93(t, 0) = 1 + 2V~ exp{—n2n2t} =
f ^n=1
Vnt
exp
n
T
Ряд слева быстро сходится при больших значениях Ь, а ряд справа — при малых Ь (которые нас сейчас и интересуют). Если определить 6(Ь) = ^ехр{—пп2Ь} (п е Ъ, Ь > 0), то получим соотношение 6(1/Ь) = лД,6(Ь), или лД^ехр{-пп2Ь} = ^ехр{—пп2/Ь} (п е Ъ), известное как функциональное уравнение для тэта-функции [14].
Проведем вспомогательные преобразования, разбивая промежуточную сумму на нечетные и четные слагаемые:
1 2 ^^
X>xp{- у} = 1+2Е exp{-n2n2t}:
n=1
Vnt i t
= 1 + 2E + 2 e exp{—k2n24t}
fc=i
/ 1 ^ 2V/ + ^=Vexp
n 4t
Отсюда (после вычитания из первого ряда последнего и удвоения) получаем искомое разложение при в > 0:
e(s)=4^e xp{
-n2n2s
1 - П >}=1-Q, Q(s)
n=1
дробь (слабая особенность под интегралом) быстро убывает от бесконечности (т = Ь) до практически нуля (т = Ь — 1). При Ь > 1 нижний предел интегрирования можно заменить на Ь — 1. Таким образом, предысторию, которую можно не учитывать, нетрудно установить в исходном физическом времени по связи 1 = ¿2 — = Б(Т(т))йт/£2. Компактное функционально дифференциальное (ю = йю/йЬ) ТДС-уравнение (15) с начальными данными )(0) = 1 заменяет нелинейную краевую задачу с динамическими граничными условиями в том смысле, что для построения ТДС-спектра нужна формально только динамика поверхностной концентрации (десорбции).
Численный метод и результаты моделирования
Для определенности ориентируемся на опубликованные данные для никеля и стали (марки 12Сг18№10Т1) [3], вольфрама [4] и бериллия [20, 22]. Оценки существенно зависят от условий эксперимента и подготовки образцов, так что значения воспринимаем как модельные для численных иллюстраций.
Общие величины: £ = 0.1 см, То = 300 К. Принятые значения параметров модели для стали: Ь0 = 1.3 х 10-9 см2/с, Еь = 80 х 103 Дж/моль, Б0 = 3.1 х 10-4 см2/с, Ев = 31 х 103, д = 50см-1, с = 5 х 1017 1/см3 (1 = ат.Н), Т = 1 К/с; для никеля: Ь0 = 1.53х10-14, Еь = 43.2х 103, Бо = 7.5-х10-3, Ев = 40х 103, д = 102, в0 = 1.8 х 10-2, Е3 = 61.4 х 103, с = 9.9 х 1017 1/см3, р = 37.4торр, Т = 770, Т = 0.5; для вольфрама: Ь0 = 6 х 10-4, Еь = 69.56 х 103, Б0 = 4.1 х 10-3, Ев = 37.63 х 103,
с =
q = — exp{- 1/(4s)}. Ряд Q быстро сходится при малых значениях s. При s ^ +0 имеем П ^ 0 и интегрируемую особенность в « 1/у^. График функции Q(s) имеет S-образный («аррениусовский») вид кривой насыщения и уже П(1) ~ 0.9996. Функция Q(s)/vns монотонно возрастает до максимума (~ 0.828 при s ~ 0.334) с последующим монотонным затуханием. В представлении Q(s) с помощью ряда достаточно ограничиться небольшим числом слагаемых (5-8). В получаемом вместо (14) с помощью описанного представления e(s) уравнении
v(t) = —b(t)v2(t) — xf1 —Vv(t) dr (15) Jo у n[t — т]
38
д = 104, в0 = 9 х 10-3, Еа = 15 х 10 7.8 х 1014, р = 670 торр, Т = 1300, Т = 5; для бериллия: Ь0 = 3.1 х 10-9, Еь = 57.43 х 103, Б0 = 3 х 10-3, Ев = 28 х 103, д = 2 х 102, в0 = 1.44 х 10-4, Е3 = 1.82 х 103, с = 5.2 х 1017, р = 760 торр, Т = 1150, Т = 5.
Основную роль в динамике дегазации играет квадратичная десорбция. Поэтому аппроксимировать будем интегральное слагаемое в уравнении (15). При этом горизонт последействия равен Н < 1 (т е [Ь, Ь — Н] в безразмерном времени, ф(1) ~ 0.999). Ориентируясь на формулу трапеций численного интегрирования и график гладкой функции Q(в)/^/■кв (рис.3), фиксируем Н ~ 0.3—0.4. Замена Ь = ЦБйт/£2 ориентирована на характерное время диффузии £2/Б, так что такое Н соответствует существенному отрезку физического времени в эксперименте. Будем последовательно на отрез-
ках безразмерного времени длины Н аппроксимировать ТДС-уравнение (15) системой ОДУ.
1 0.8 0.6 0.4
0.2 0
- ^^^^ Q(s) f'-Q( ' )d
4 Q(s) N - q= cxpi l/(4.v)!.
n=1 7V> 10 s ~ 0.334
0 0.2 0.4 0.6 0.8
Рис. 3. Функции Q(s) и Q(s)/y/ns
Из-за неограниченной (но интегрируемой) особенности основной вклад в интеграл вносят значения г)(т), т « t. Поэтому с учетом вогнутости функции vv используем квадратичное приближение vv(t) « vv(t) + A[t — т] + B[t — т]2. Рассмотрим текущий отрезок t G [kh, (k + 1)h], k ^ 0. Условия
vv(t)|fch = vv(kh), / vv(t) dT = v(t) — v(kh) Jkh
(полагаем v(t) « v(kh) + vv(t)[t — kh]) определяют искомые значения A(t), B(t) (const по т):
vv(t) — vv(kh) r vV(т) « vv(t) + 2 ( ; ) [t — т] —
- 3
t—kh v(t) — v(kh)
[t — т]2, (16)
[£ - кН]2
т е [кН, £] (£ е [кН, (к + 1)Н]). Представляем интеграл в (15) суммой: т е [0, кН], т е [кН, £]. Во втором слагаемом интегрируем в аналитическом виде особенность ¿(т)/\/£ — т с подстановкой (16), а для интеграла без особенности (^(+0)^^+0 = 0) применяем теорему о среднем и формулу трапеций:
г г
' : ¿(т) ат = ¿(4)/ —, , , ат^
kh
V^t—т]
kh
v^ît—^I
оставим справа несколько слагаемых с учетом показателя п2 (для определенности п = 1, 3). При этом отбросили отрицательные слагаемые (V < 0). Компенсируем это заменой кН на £ и введем дополнительно переменные г1,2 (£):
= [ ехр{—(2г — 1)2п2[£ — т]}-(т) ат,
Jo
щ(г) = — (2г — 1)2п2+ -(£), ^1,2(0) = 0. В итоге вместо (15) получаем систему ОДУ:
-(£) = — Ь(£)-2(£) + А [15 + ^гК-(кН) —
уП 15 4 — 4К [^1(£)+ М2(*)], (17)
го 1(£) = ... — [4х+п2]ад1(£)— 4Кг2(£), гг2(£) = ... — 4К ш1(£) — [4К + 9п2]г2(£),
Ь =-... Ь ,]-, и(0) = 1, М1,2(0)=0,
-, + к Г32 Q(ifc)l /Г-
1+v^L 15 r~ JVtk
к
к =
1 + к Г 32 Q(tfc) ] /г-1 + vH 15 ^Jv k
, tk = t — kh.
vv(kh) + vv(t) 1 ^ , ,4 ,-—
« 1 'V • - • Q(t — kh)V t — kh.
2 y n 2
Для аппроксимации первого интеграла kh kh
Троеточие обозначает первую строку первого уравнения (—Ь... ¿(кН)). Ряд для ф(з) сходится очень быстро, в программной реализации достаточно 5-8 членов.
Для вычисления искомой функции -(£) = 5 (£)/<? на текущем отрезке безразмерного времени (в системе (17) £ = е [кН, (к + 1)Н], к ^ 0) используем стандартный метод Рунге -Кутты 4-го порядка для интегрирования ОДУ (авторы пользовались пакетом БсИаЪ). В случае двухпиковой структуры ТДС-спектра для повышения точности моделирования целесообразно использовать дополнительно переменные гз,4. После возвращения к физическому времени получаем модельный ТДС-спектр 3(Т) = Ь(Т)^2(£) (£' о £ О Т(£)).
Качественный вид ТДС-спектра металлических конструкционных материалов типичен. На рисунках 4-7 представлены численно смоделированные ТДС-спектры для указанных выше значений параметров водородопроница-емости конструкционных материалов при варьировании скоростей нагрева. Принятая модель и без учета дефектов (ловушек) позволяет получать различный вид двухпиковых графиков (когда интенсивнее низкотемпературный всплеск либо же пики сравнимы). В верхней части графиков для сравнения отмечено изменение так называемого транспортного параметра Ш = ^Ь^с/^ = ^Ьс/(^д2), который является определяющим при исследовании водородопроницаемости мембран [4]. Ори-
s
ентировочно Surface Limited Regime (SLR) соответствует W < 10-2, а лимитирование диффузией (DLR) соответствует W > 104. Для бериллия и стали низкотемпературный ТДС-пик происходит в диапазоне, когда сопоставимую роль играют диффузия и процессы на поверхности. Высокотемпературный пик наблюдается уже в зоне DLR. Для никеля поверхностные процессы и диффузия оказывают сравнимое воздействие во всем диапазоне температур, а на графике при низких температурах видим лишь ступень, напоминающую пик.
J(T)-IQ'14 "3 -4
Beryllium T = 1Q, 5, 2-5
300 400 500 600 700 800 900 1 000
Рис. 4- ТДС-спектры, бериллий
J(T) 'Ю-14 40
16
50
60
12
_ Nickel X \
T = 1, 0.5, 0.25 / \
TNi
W
0
300 400 500 600 700 800 900 1 000
Рис. 5. ТДС-спектры, никель
Подробнее остановимся на рисунке 7, поскольку такой ТДС-спектр качественно отличается от предыдущих. Разумеется, приведенные выше значения параметров для вольфрама являются условными. Помимо микропримесей существенное влияние имеют особенности обработки поверхности образца. Тем более что в экспериментальной практике пластины тонкие, объем мал и ярче влияние поверхност-
ных процессов. Это одна из причин разброса количественных оценок параметров водородо-проницаемости. Другая причина (и не последняя): при обработке экспериментальных данных используются разные модели (хотя формально коэффициенты носят одно и то же название). На принятых модельных данных обнаружен узкий всплеск с последующим затяжным выходом на второй пик, менее выраженный. Модель такова, что в образце нет достаточно емких ловушек (стандартное уравнение диффузии), но более детальное внимание уделяется поверхности (см.(4)—(6), пластина тонкая). Сначала активно десорбируется приповерхностный водород. По мере нагрева медленно активируется диффузия (градиент концентрации значителен, подкачка к поверхности растет). Но обычно такие явления объясняют наличием ловушек, а не динамикой взаимодействия «поверхность-объем».
J(T)-10 12 10 8 6 4 2
,-13
10 10
10
10
12Cr18Ni10Ti DLR ^^
T = 2, 1, 0.5
J Ts) TSt
W
0
300 400 500 600 700 800 900 1 000
Рис. 6. ТДС-спектры, сталь
J(T) -10 4
-12
810
910
300 400 500 600 700 800 900 1 000
Рис. 7. ТДС-спектры, вольфрам
В работе [18] на рисунке 1 имеется удивительно похожий экспериментальный график (помеченный черными кружками). Хотя там речь идет о примесях в вольфраме и дейтерии. Смысл нашего модельного результата в том, что помимо стандартного упоминания ловушек следует иметь в виду и другие возможные причины различных ТДС-пиков.
На рисунках 8-11 отслеживается раздельно динамика поверхностных процессов и диффузии. Высокотемпературным пикам способствует существенно возросший с нагревом коэффициент десорбции и активизация диффузионного притока из объема. Пики при относительно низких температурах происходят при сравнительно небольшом диффузионном притоке к поверхности, но больших концентрациях в приповерхностном объеме.
10
10"' -2
10
И
10
10:
-6
10
10
эе—эе Л?-У- Beryllium Nickel
ЬЛу* i i ^ i
/1 /'Be /Ь/ь0 TNi i i i т2 TBe
0~ -*--
10 Vq)2 ¡
\(q/q)2
ю"3 V i Л 1 \ i \ i \ i
-6 10 \ i i Pv 1
-9 10 yS 1 yf / 1 / /Ъ/ъ0 -12Cr18M10T/
/ # 1 » —»- Tungsten
1012 тА т1 Tw >st< / Ts2 T 2 1 W
стины) наименее выражен для никеля - достаточно быстро устанавливается концентрация, близкая к равномерной. В бериллии и вольфраме, напротив, уже в начале эксперимента происходит значительное снижение концентрации у поверхности и прогиб профиля концентрации в дальнейшем при нагреве сохраняется. Для стали заметим, что лишь после достижения достаточно высокой температуры начинается активная дегазация образца.
4 JDiff ■ 10 /c
6 5 J - —D — J°iff dx x=t / Tungsten ч ¡r I \ /\ 1 X s \ 1 N
4 Beryllium/ у s \ 1 ■f \ 1 r \ i \ i \ i
3 \ i \i Pk
2 i \ i \ i \ i \
1 i \ i \ i \
0 J/Tbs , ТВе T 2 1 W
Рис. 10. Диффузионный приток к поверхности
4
JDiff ■ 10 4/c
300 400 500 600 700 800 900 1 000
Рис. 8. Поверхностные параметры (Бе, N1)
300 400 500 600 700 800 900 1 000
300 400 500 600 700 800 900 1 000
Рис. 11. Диффузионный приток к поверхности
Представим в компактной форме алгоритм численного моделирования ТДС-спектра.
Рис. 9. Поверхностные параметры С^ сталь)
Рисунок 12 иллюстрирует динамику концентрации (недоступную для измерений) внутри исследуемого материала. Прогиб профиля концентрации (от середины к краю пла-
• Задаем значения параметров модели «о, Е8, Ь0, Еь, £>о Ед, д. При условиях насыщения рс, Т определяем величины с, &(£'), к = согласно изложенному на стр.37. На задании значения в(Т) можно «сэкономить», если известна равновесная
Beryllium
с / с 1
0.8
0.6 0.4
0.4
0.3 0-2 X/С
1100 о
12Cr18Ni10Ti
S / с
\
1
0.8
I lu^, ^^Ш/г Л
900 ^r o.l 1100 0
0-2 x/f
1100 0
Рис. 12. Динамика распределения объемной концентрации внутри исследуемого материала
объемная концентрация с растворенного диффузионно подвижного водорода. Возможны другие сценарии начального насыщения (см. стр. 37). Соответствующая корректировка начальных данных определяется особенностями реализации конкретного ТДС-эксперимента.
• В безразмерном времени поэтапно (на отрезках длины h) численно интегрируем ODE-систему невысокого порядка вида (17). В [27] представлены и другие ODE-аппроксимации. Это стандартная операция для любого современного математического пакета (включая freeware).
• По функции v(t') = q(t')/q в исходном физическом времени строим график плотности термодесорбции J(t) = b(T(t))q2(t) или ТДС-спектр J(T).
Обратная задача параметрической идентификации
Представленный алгоритм численного моделирования позволяет быстро сканировать различные сценарии и условия эксплуатации материала (включая закон нагрева и экстраполирование результатов с ростом £). Это полезная статистическая информация при выработке стратегии экспериментальных исследований. Для относительно новых материалов (различных сплавов) сначала нужно оценить параметры водородопроницаемости. Здесь мы сталкиваемся с обратными задачами математической физики, трудности решения которых хорошо известны. Будем считать, что коэффициент диффузии Б(Т) известен (обычно используется метод Дайнеса - Бэррера в БЬИ-режиме проницаемости). Поставим задачу: по информации .](¿) (потоку десорбции) оценить
поверхностные параметры Ьо, Еь, д. При этом мы находимся в реальных условиях, когда динамика десорбции тесно связана с диффузией в объеме. Задача оценки коэффициента объемной десорбции (Ьуо! = Ь/д2), когда накоплением на поверхности можно пренебречь (д ~ 0), рассмотрена в работе [29].
Поскольку функция = 0(Т(¿)) из-
вестна, то целесообразно сразу перейти к безразмерному времени ориентированному на характеристическое время диффузионного переноса Чтобы не усложнять выкладки, оставим прежнее обозначение ¿. Представим нелинейное слагаемое в системе (17) в следующем виде: Ь(ф2(£) = а(к) 3(¿)/[£(£)с], 3(*) = (¿). Здесь а — функция параметра к, получаемая преобразованиями в силу обозначений в (14) и (17). Элементарные, но несколько громоздкие формулы опускаем. Функция 3(¿) (£ = ¿') известна по измерениям. После подстановки в (17) получаем систему трех линейных ОБЕ. Для повышения точности вычислений целесообразно добавить переменные -ш3,4. Измерения обычно зашумлены, но функция 3(¿) входит в правую часть ОБЕ и при интегрировании происходит сглаживание.
Обратим внимание на то, что правые части уравнений теперь зависят только от одного оцениваемого коэффициента растворения к = д^. Это позволяет провести декомпозицию алгоритма параметрической идентификации.
Задание значения д позволяет в итоге рассчитать зависимость = д) в исходном физическом времени ¿. Теперь вспомним, что
J = bq2 ^ b = Jq 2 ^ ln bo —
R • 103 T
103
— = ln (18)
где ^(t) = J/q2, T = T(t). В силу монотонности нагрева t о T(t), что позволяет ввести координаты X = 103/T, Y = ln На плоскости (X, Y) получаем параметрически заданную кривую X (t), Y(t). Судя по соотношению (18), это должна быть прямая линия с отрицательным угловым коэффициентом. Отсюда следует критерий выбора «правильного» значения g: кривая X(t), Y(t) = ln{J(t)/q2(t;g)} должна оказаться отрезком прямой на плоскости (X, Y). Эта подзадача скалярная, варьируется только g.
Формально продолжаем отрезок прямой до пересечения с осями координат. Пересечение с осью Y (X = 0) дает значение ln b0. Точка пересечения с осью абсцисс X = 103 х R ln b0/Eb определяет энергию активации Еь.
Трудоемкость алгоритма определяется итерационным обращением к подпрограмме численного решения начальной задачи для систе-
мы ODE. Это требует определенных усилий и навыков работы с математическими пакетами, но использование стандартной встроенной операции несравнимо проще и быстрее итерационного решения исходной нелинейной краевой задачи с тремя варьируемыми параметрами (требующей специализированного программного обеспечения).
-1 с >12 2 1 1 И03
x = 103R y e;1 Tungsten T
y = ln bo 1.2 x
-10 K- L x=*
-20 0.8;
-30
ln{J(T(f)) q-2(T(t),g)}
Рис. 13. Оценка параметров десорбции, вольфрам
0
0.5
1.5
2.5
103
-25 -30 -35 -40 -45 -50
y - In 60 12Cr18№10T/ T
1.2 к
/ 0.8 -
- 103Ry - Eb ~ - 2.128
ln{J(T(t)) q2(T(t),g)}
Рис. 14. Оценка параметров десорбции, сталь
Предложенный метод выпрямления проиллюстрирован на рисунках 13-16. Для тестирования алгоритма были построены модельные спектры (при g = const). Для получения более точных оценок на двухпиковых спектрах численно решалась система из пяти уравнений вида (17) для переменных v(t), wi-4(i). Дополнительно проводилась серия экспериментов на спектрах, в которые добавлена случайная погрешность до 20% (использовалась стандартная функция Scilab генерирования случайных чисел, распределенных по равномерному зако-
ну). Алгоритм идентификации на основе интегрирования системы ОДУ (с подстановкой 7(¿) в правые части) показал помехоустойчивость обработки входных данных.
Рис. 16. Оценка 6, д (N1, зашумленный спектр) Начальное значение д дает баланс
+ 2<? = + 2д-2] = 2 7(т) йт.
./о
Эксперименты показали, что точность такой оценки снижается при д ^ 1. На графиках представлены кривые с использованием «правильного» численного значения к, при котором полученная кривая близка к отрезку прямой, и кривые, получаемые при отклонении от к на 20%. Относительная погрешность тестирования алгоритма идентификации не превышает нескольких процентов (сказываются вычислительные погрешности решения как прямой, так и обратной задачи). Заметна высокая чувствительность (особенность типа «вер-
тикальный клюв») при превышении «истинного» значения к. Математическая особенность возникает из-за того, что формально функция г(£) меняет знак, принимая и отрицательные значения. Авторы не заботились о том, чтобы избежать такого «нефизического» перехода, поскольку это ярко свидетельствует о неправильном выборе значения к. Такие же численные эксперименты проведены для за-шумленных графиков. С более высокой точностью восстанавливаются энергии активации, поскольку энергетические характеристики активнее влияют на десорбцию по мере нагрева. С математической точки зрения совпадения модельных и экспериментальных кривых добиваться не обязательно (с учетом экспериментальных погрешностей в десятки процентов). ОДУ-аппроксимация вполне адекватна.
Заключение
В статье анализируется проблема идентификации спектров термодесорбции водорода из конструкционных материалов, применяемых в ядерной энергетике. На качественном уровне идентификация состоит в выявлении причин всплеска десорбции. Обычно ТДС-пики связывают с высвобождением водорода из ловушек с различными энергиями связи. Как показано на простейшей модели диффузии (однородный материал), учет динамики поверхностных процессов может дать двухпи-ковую структуру даже на очень тонких экспериментальных образцах. Тенденция списывать все на «теорию различных ловушек» понятна, но объем образца практически нулевой для проявления емкости захвата.
Другая часть проблемы — параметрическая идентификация. Имеется большое разнообразие моделей водородопроницаемости, включая различные граничные условия (что влечет вариативность оценок коэффициента диффузии даже по порядку). ТДС-эксперимент имеет естественные ограничения информативности. Результаты должны корректно (единственным образом) численно экстраполироваться, например, с ростом толщины материала (экспериментальные образцы очень тонкие из-за ограничений времени измерений). В статье нелинейная краевая задача (стандартное уравнение диффузии, но с динамическими граничными условиями) сведена к функционально дифференциальному уравнению для поверхностной концентрации, поскольку для построения ТДС-спектра требуется лишь динамика десорбции. Предложен эффективный вычислительный алгоритм, ориентированный на использование математи-
Рис. 15. Оценка 6, д (Бе, зашумленный спектр) 0 0.5 1 1.5 2 2.5 3 1рз
-28 Ы'юкв! 7
-Ъ2Л У = 1п Ьо 1.2-
-36 -40 V _ 103ЯУ~ чХ х = ^ ~ - 6.122
// 0.8 '
-44
-48
ВД7-(/)) ^(ДО.д)}
ческих пакетов (в том числе freeware). Основной итоговый результат статьи — геометрически прозрачный метод решения обратной задачи идентификации поверхностных параметров, когда десорбция динамически взаимосвязана с диффузией в объеме.
Работа выполнена при поддержке РФФИ (грант №15-01-00744).
Литература
1. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.
2. Водород в металлах / Ред. Г. Алефельд и И. Фелькль. М.: Мир, 1981. Т. 1, 506 с.; т. 2, 430 с.
3. Изотопы водорода. Фундаментальные и прикладные исследования / Ред. А. А. Юхим-чук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697с.
4. Писарев А. А., Цветков И. В., Марен-ков Е.Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
5. Andronov D. Yu., Arseniev D. G., Polyanskiy A.M., Polyanskiy V.A., Yakovlev Yu. A. Application of multichannel diffusion model to analysis of hydrogen measurements in solid //Int. J. Hydrogen Energy. 2017. Vol. 42. P. 699-710. doi: 10.1016/j.ijhydene.2016.10.126
6. Belyaev A.K., Polyanskiy A.M., Polyanskiy V.A., Sommitsch Ch., Yakovlev Yu. A. Multichannel diffusion vs TDS model on example of energy spectra of bound hydrogen in 34CrNiMo6 steel after a typical heat treatment // Int. J. Hydrogen Energy. 2016. Vol. 41. P. 86278634. doi: 10.1016/j.ijhydene.2016.03.198
7. Castro F. J., Meyer G. Thermal desorption spectroscopy (TDS) method for hydrogen desorption characterization (I): theoretical aspects// J. Alloys and Compd. 2002. Vol. 330332. P. 59-63.
8. Denisov E. A., Kompaniets M. V., Kompa-niets T. N., Bobkova I. S. Peculiarities of hydrogen permeation through Zr-1%Nb alloy and evaluation of terminal solid solubility // J. Nucl. Mater. 2016. Vol. 472. P. 13-19. doi: 10.1016/j.jnucmat.2016.01.022
9. Evard E. A., Gabis I. E., Yartys V. A. Kinetics of hydrogen evolution from MgH2: experimental studies, mechanism and modelling // Int. J. Hydrogen Energy. 2010. Vol. 35. P. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
10. Hale J. Theory of Functional Differential Equations. Springer-Verlag, 1977. 366 p.
11. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
12. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen
influence// Acta Mechanica. 2008. Vol. 195. P. 295-304. doi: 10.1007/s00707-007-0568-z
13. Legrand E., Oudriss A., Savall C., Bouhat-tate J., Feaugas X. Towards a better understanding of hydrogen measurements obtained by thermal desorption spectroscopy using FEM modeling// Int. J. Hydrogen Energy. 2015. Vol. 40. P. 2871-2881. doi: 10.1016/j.ijhydene.2014.12.069.
14. Lang S. Elliptic functions. Addison-Wesley publishing. 1973. 326 p.
15. Lototskyy M.V., Yartys V.A., Pollet B.G., Bowman Jr. R. C. Metal hydride hydrogen compressors: a review // Int. J. Hydrogen Energy. 2014. Vol. 39. P. 5818-5851. doi: 10.1016/j.ijhydene.2014.01.158
16. McRae G.A., Coleman C.E., Leitch B.W. The first step for delayed hydride cracking in zirconium alloys//J. Nucl. Mater. 2010. Vol. 396. P. 130-133. doi: 10.1016/j.jnucmat.2009.08.019
17. Mieza J. I., Vigna G. L., Domizzi G. Evaluation of variables affecting crack propagation by delayed hydride cracking in Zr-2.5Nb with different heat treatments // J. Nucl. Mater. 2011. Vol. 411. P. 150-159. doi: 10.1016/j.jnucmat.2011.01.101
18. Oya Y., Inagaki Y., Suzuki S. et al. Behavior of hydrogen isotope retention in carbon implanted tungsten// J. Nucl. Mater. 2009. Vol. 390--391. P. 622-625. doi: 10.1016/j.jnucmat.2009.01.175
19. Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface// Int. J. Hydrogen Energy. 2011. Vol. 36. P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
20. Samsonov A. V., Korenkov A. Yu., Gabis I. E., Kurdyumov A. A. Limiting role of desorption in hydrogen transport across a deposited beryllium film// Tech Phys. The Russian J. Appl. Phys. 1998. Vol. 43, no. 1. P. 114-116.
21. Schmid K., Rieger V., Manhard A. Comparison of hydrogen retention in W and W/Ta alloys // J. Nucl. Mater. 2012. Vol. 426. P. 247-253. doi: 10.1016/j.jnucmat.2012.04.003
22. Tazhibaeva I. L., Shestakov V. P., Klepi-kov A. Kh., Romanenko O. G., Chikhray Y. V., Kenzhin E.A., Zverev V. V. Modeling of Hydrogen release from irradiated beryllium // National Nucl. Center of the Republic of Kazakhstan Bull. 2000. Vol. 1. P. 37-41.
23. The hydrogen economy. Eds. M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
24. Varin R.A., Czujko T., Wronski Z.S. Nano-materials for solid state hydrogen storage. NY: Springer, 2009. 338 p.
25. Whittaker E. T., Watson G. N. A Course of Modern Analysis. Cambridge University Press, 1996. 612 p.
26. Zaika Yu. V., Bormatova E. P. Parametric identification of hydrogen permeability model by delay times and conjugate equations // Int. J. Hydrogen Energy. 2011. Vol. 36. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
27. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermal desorption by ODE-approximation// Int. J. Hydrogen Energy. 2017. Vol. 42. P. 405-415. doi: 10.1016/j.ijhydene.2016.10.104
28. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption // Adv. in Mater. Sci. and Appl. 2014. Vol. 3, iss. 3. P. 120129. doi: 10.5963/AMSA0303003
29. Zaika Yu. V., Kostikova E. K. Determination of effective recombination coefficient by
thermodesorption method // Int. J. Hydrogen Energy. 2014. Vol. 39. P. 15818-15826. doi: 10.1016/j.ijhydene.2014.07.117
30. Zaika Yu. V., Rodchenkova N. I. Boundary-value problem with moving bounds and dynamic boundary conditions: diffusion peak of TDS-spectrum of dehydriding // Appl. Math. Model. 2009. Vol. 33, iss. 10. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018
31. Zaika Yu. V. The solvability of the equations for a model of gas transfer through membranes with dynamic boundary conditions // Comp. Maths. Math. Phys. 1996. Vol. 36, iss. 12. P. 17311741.
Поступила в редакцию 23.05.2017
References
1. Vzaimodeistvie vodoroda s metallami [Interaction of hydrogen with metals]. Ed. A. P. Zakharov. Moscow: Nauka, 1987. 296 p.
2. Hydrogen in metals. Eds. G. Alefeld, J. Volkl. Berlin: Springer-Verlag, 1978.
3. Izotopy vodoroda. Fundamental'nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies]. Ed. A. A. Yuk-himchuk. 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. Andronov D. Yu., Arseniev D. G., Polyanskiy A. M., Polyanskiy V. A., Yakovlev Yu. A Application of multichannel diffusion model to analysis of hydrogen measurements in solid. Int. J. Hydrogen Energy. 2017. Vol. 42. P. 699-710. doi: 10.1016/j.ijhydene.2016.10.126
6. Belyaev A.K., Polyanskiy A.M., Polyanskiy V.A., Sommitsch Ch., Yakovlev Yu. A. Multichannel diffusion vs TDS model on example of energy spectra of bound hydrogen in 34CrNiMo6 steel after a typical heat treatment. Int. J. Hydrogen Energy. 2016. Vol. 41. P. 8627-8634. doi:10.1016/j.ijhydene.2016.03.198
7. Castro F. J., Meyer G. Thermal desorption spectroscopy (TDS) method for hydrogen desorption characterization (I): theoretical aspects. J. Alloys and Compd. 2002. Vol. 330-332. P. 59-63.
8. Denisov E. A., Kompaniets M. V., Kompa-niets T. N., Bobkova I. S. Peculiarities of hydrogen permeation through Zr-1%Nb alloy and evaluation of terminal solid solubility. J. Nucl. Mater. 2016. Vol. 472. P. 13-19. doi: 10.1016/j.jnucmat.2016.01.022
9. Evard E. A., Gabis I. E, Yartys V. A. Kinetics of hydrogen evolution from MgH2: experimental studies, mechanism and modeling. Int. J. Hydrogen Energy. 2010. Vol. 35. P. 9060-9069. doi: 10.1016/j.ijhydene.2010.05.092
10. Hale J. Theory of Functional Differential Equations. Springer-Verlag, 1977. 366 p.
11. Handbook of hydrogen storage: new materials for future energy storage. Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
12. Indeitsev D. A., Semenov B. N. About a model of structure-phase transfomations under hydrogen influence. Acta Mechanica. 2008. Vol. 195. P. 295304. doi: 10.1007/s00707-007-0568-z
13. Legrand E., Oudriss A., Savall C., Bouhattate J., Feaugas X. Towards a better understanding of hydrogen measurements obtained by thermal desorption spectroscopy using FEM modeling. Int. J. Hydrogen Energy. 2015. Vol. 40. P. 2871-2881. doi: 10.1016/j. ijhydene. 2014.12.069
14. Lang S. Elliptic functions. Addison-Wesley publishing, 1973. 326 p.
15. Lototskyy M. V., Yartys V. A., Pollet B. G, Bowman Jr. R. C. Metal hydride hydrogen compressors: a review. Int. J. Hydrogen Energy. 2014. Vol. 39. P. 5818-5851. doi: 10.1016/j.ijhydene.2014.01.158
16. McRae G.A., Coleman C.E., Leitch B.W. The first step for delayed hydride cracking in zirconium alloys. J. Nucl. Mater. 2010. Vol. 396. P. 130-133. doi: 10.1016/j.jnucmat.2009.08.019
17. Mieza J. I., Vigna G. L., Domizzi G. Evaluation of variables affecting crack propagation by delayed hydride cracking in Zr-2.5Nb with different heat treatments. J. Nucl. Mater. 2011. Vol. 411. P. 150-159. doi: 10.1016/j.jnucmat.2011.01.101
18. Oya Y., Inagaki Y., Suzuki S. et al. Behavior of hydrogen isotope retention in carbon implanted tungsten. J. Nucl. Mater. 2009. Vol. 390--391. P. 622-625. doi: 10.1016/j.jnucmat.2009.01.175
19. Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface. Int. J. Hydrogen Energy. 2011. Vol. 36. P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
20. Samsonov A. V., Korenkov A. Yu., Gabis I. E., Kurdyumov A. A. Limiting role of desorption in hydrogen transport across a deposited beryllium film. Tech Phys. The Russian J. Appl. Phys. 1998. Vol. 43, no. 1. P. 114-116.
21. Schmid K., Rieger V., Manhard A. Comparison of hydrogen retention in W and W/Ta alloys. J. Nucl. Mater. 2012. Vol. 426. P. 247-253. doi: 10.1016/j.jnucmat.2012.04.003
22. Tazhibaeva I. L., Shestakov V. P., Klepi-kov A.Kh., Romanenko O.G., Chikhray Y.V., Kenzhin E.A., Zverev V. V. Modeling of Hydrogen release from irradiated beryllium. National Nucl. Center of the Republic of Kazakhstan Bull. 2000. Vol. 1. P. 37-41.
23. The hydrogen economy. Eds. M. Ball, M. Wietschel. Cambridge Univ. Press, 2009. 646 p.
24. Varin R.A., Czujko T., Wronski Z. S. Nano-materials for solid state hydrogen storage. NY: Springer, 2009. 338 p.
25. Whittaker E. T, Watson G. N. A Course of Modern Analysis. Cambridge Univ. Press, 1996. 612 p.
СВЕДЕНИЯ ОБ АВТОРАХ:
Заика Юрий Васильевич
рук. лаб. моделирования природно-технических систем, д. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: zaika@krc.karelia.ru тел.: (8142) 780059
Костикова Екатерина Константиновна
научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: kostikova@krc.karelia.ru тел.: (8142) 766312
26. Zaika Yu. V., Bormatova E. P. Parametric identification of hydrogen permeability model by delay times and conjugate equations. Int. J. Hydrogen Energy. 2011. Vol. 36. P. 1295-1305. doi: 10.1016/j.ijhydene.2010.07.099
27. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermal desorption by ODE-approximation. Int. J. Hydrogen Energy. 2017. Vol. 42. P. 405-415. doi: 10.1016/j.ijhydene.2016.10.104
28. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption. Adv. in Mater. Sci. and Appl. 2014. Vol. 3, iss. 3. P. 120-129. doi: 10.5963/AMSA0303003
29. Zaika Yu. V., Kostikova E. K. Determination of effective recombination coefficient by thermodesorption method. Int. J. Hydrogen Energy. 2014. Vol. 39. P. 15818-15826. doi: 10.1016/j.ijhydene.2014.07.117
30. Zaika Yu. V., Rodchenkova N. I. Boundary-value problem with moving bounds and dynamic boundary conditions: diffusion peak of TDS-spectrum of dehydriding. Appl. Math. Model. 2009. Vol. 33, iss. 10. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018
31. Zaika Yu. V. The solvability of the equations for a model of gas transfer through membranes with dynamic boundary conditions. Comp. Maths. Math. Phys. 1996. Vol. 36, iss. 12. P. 1731-1741.
Received May 23, 2017
CONTRIBUTORS:
Zaika, Yury
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Russia e-mail: zaika@krc.karelia.ru tel.: (8142) 780059
Kostikova, Ekaterina
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Russia e-mail: kostikova@krc.karelia.ru tel.: (8142) 766312
