CC BY
38
Труды Карельского научного центра РАН №10. 2015. С. 42-53 DOI: 10.17076/mat147
УДК 519.6:539.2
АППРОКСИМАЦИЯ КРАЕВОЙ ЗАДАЧИ ТЕРМОДЕСОРБЦИИ ВОДОРОДА СИСТЕМОЙ ОДУ
Ю. В. Заика, Е. К. Костикова
Институт прикладных математических исследовании Карельского научного центра РАН
В рамках технологических задач водородного материаловедения (включая проект ITER) ведется интенсивный поиск различных по назначению конструкционных материалов с заранее заданными пределами водородопроницаемости. Одним из экспериментальных методов является термодесорбционная спектрометрия (ТДС). Образец, насыщенный водородом, дегазируется в условиях ваку-умирования и монотонного нагрева. С помощью масс-спектрометра регистрируется десорбционный поток, позволяющий судить о характере взаимодействия изотопов водорода с твердым телом. Интерес представляют такие параметры переноса, как коэффициенты диффузии, растворения, десорбции. В статье представлены распределенная краевая задача термодесорбции и численный метод моделирования ТДС-спектра, требующий лишь интегрирования нелинейной системы обыкновенных дифференциальных уравнений (ОДУ) невысокого порядка (по сравнению, например, с методом прямых).
Ключевые слова: водородопроницаемость, термодесорбция, нелинейные краевые задачи, динамические граничные условия, численное моделирование.
Yu.V. Zaika, E. K. Kostikova. APPROXIMATION OF THE BOUNDARY-VALUE PROBLEM OF HYDROGEN THERMAL DESORPTION BY ODE SYSTEM
One of the technological challenges for hydrogen materials science (including 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 thermodesorption 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. We are interested in such transfer parameters as the coefficients of diffusion, dissolution, desorption. The paper presents a distributed boundary value problem of thermal desorption and a numerical method for TDS-spectrum simulation, where only integration of a non-linear system of low order (compared with, e. g., the method of lines) ordinary differential equations (ODE) is required. This work is supported by the Russian Foundation for Basic Research (project N15-01-00744).
Key words: hydrogen permeability, thermal desorption, nonlinear boundary-value problems, dynamical boundary conditions, numerical simulation.
Введение
Интерес к взаимодействию водорода с различными материалами носит многоплановый характер [1-10]. Достаточно упомянуть задачи энергетики, защиты металлов от водородной коррозии, проектирования химических реакторов, ракетостроения. Гидриды позволяют удерживать большое количество водорода. С этим связаны перспективы водородных аккумуляторов и двигателей с высоким уровнем безопасности: без высоких давлений и низких температур. На обратимом легировании металлов водородом основаны пластифицирование и термоводородная обработка титановых сплавов. Некоторые частные задачи водородного материаловедения исследованы в [11-15]. Энтузиасты говорят не только о водородной энергетике, но и о водородной экономике [7].
Экспериментальный опыт показывает, что лимитирующими являются не только диффузия, но и физико-химические процессы на поверхности [1, 2]. Параметры переноса зависят и от технологических особенностей получения партии материала, поэтому вряд ли следует ориентироваться на получение «табличных данных», нужны эффективные алгоритмы обработки экспериментальных кривых. В статье остановимся на моделировании термодесорбции, учитывая лишь основные факторы и информационные возможности рассматриваемого ТДС-эксперимента.
Математическая модель переноса
Рассмотрим перенос водорода сквозь образец тестируемого металла (пластину толщины £). Физико-химическую терминологию в дальнейшем используем в минимально необходимом объеме. Кратко говорим о металлической мембране, хотя это может быть многокомпонентный сплав, интерметаллид. Считаем, что нагрев относительно медленный, практически равномерный, так что градиент температуры пренебрежимо мал и диффузионный поток можно считать пропорциональным градиенту концентрации. Концентрация растворенного водорода (в атомарном состоянии) мала. Материал достаточно однороден, чтобы пренебречь взаимодействием Н с ловушками (микродефектами кристаллической структуры, которые могут удерживать водород). Для диффузии примем стандартное уравнение:
ф, х) = В(Т)ехх(1, х), (*, х) е Яг,, (1)
где * - время, = (0, ¿*) х (0, £); е(£, х) - концентрация диффундирующего водорода (ато-
марного); В - коэффициент диффузии. Зависимость В от текущей температуры Т(*) соответствует закону Аррениуса (с предэкспонен-циальным множителем Во и энергией активации Ед; Я - универсальная газовая постоянная): В = Во ехр{-Ед /[ЯТ (£)]}.
Известны более детализированные модели переноса. В частности, можно учесть несколько каналов диффузии (транскристаллическую, по границам зерен, вдоль дефектов) с взаимообменом между ними. Однако значительное увеличение числа неизвестных априори параметров делает задачу их оценки труднообозримой. В контексте прикладных задач, выделяя только лимитирующие факторы, находят компромисс между полнотой модели и возможностями ее параметрической идентификации по экспериментальным данным.
Основные трудности численного анализа связаны не с уравнением (1), ас динамическими нелинейными граничными условиями. Пусть пластина контактирует с газообразным Н2 и поверхность является существенным потенциальным барьером (см. [2, с. 177-206; Га-бис, Компаниец, Курдюмов]). Тогда с учетом (де)сорбционных процессов краевые условия моделируются следующим образом: с(0,х) = с(х), х е [0,£], * е [0,**], (2)
ео(*) = д(ТЫ*), <*(*) = д(ТМ*), (3)
4о(£) = МТ)Ро(*) - Ъ(Т^К) + Ввх(Ь, 0) , (4)
и(*) = /лз(ТМ*) - Ь(Т)д2(г) - Вех(*,£), (5)
Ь(Т) = Ьо ехр{-Еь[ЯТ]-1}, з(Т) = ...
Здесь: ео(£) = е(£,0), е^(£) = е(*,£) - граничные объемные концентрации диффундирующего атомарного водорода; до(*), -концентрации на поверхностях (х = 0,£); д(Т) - параметр локального равновесия между концентрациями на поверхности и в приповерхностном объеме; у - кинетический коэффициент; в(Т) - параметр, отражающий тот факт, что только малая часть «налетающего» водорода окажется в форме атомов на поверхности материала (для краткости будем называть в коэффициентом прилипания, имея в виду, что он характеризует итог общего процесса физадсорбции-диссоциации-хемосорбции молекулярного газа в атомы на поверхности); ро (*), - давления молекулярного газа; Ь(Т) - коэффициент десорбции.
Кратко поясним формулы (2)-(5). Если предварительно образец обезводорожен, то с(х) = 0. Соотношения (3) означают, что объ-
емные приповерхностные концентрации Со/(£) пропорционально «отслеживают» концентрации 9о/(Ь) на поверхности, скорость растворения достаточно велика. В балансе потоков (4), (5) производные слева отражают возможность накопления Н на поверхности. Чем больше внешнее давление Н2, тем больше атомов в единицу времени попадает на единичную площадку (первые слагаемые в правых частях). Вторые слагаемые означают, что часть атомов, оказавшихся на поверхности, снова соединяются в молекулы водорода и покидают поверхность (десорбционные потоки = &(£)9о^(£)). Говоря о потоках, обычно подразумеваем их плотность. Квадратич-ность десорбции характерна для водорода, исключая экстремальные температуры. Последние слагаемые в правых частях (4), (5) соответствуют притоку или оттоку атомов Н к поверхности за счет диффузии в объеме.
В модели фигурирует как молекулярный, так и атомарный водород. Для единообразия подсчет будем вести в атомах: [с] = 1/см3, [9] = 1/см2, [^Сх] = [.] = 1/см2с (. = Ъд2). В кинетической теории газов величина ц.р определяет число частиц (в данном случае молекул Н2), соударяющихся с единичной площадкой поверхности в единицу времени. Но за счет множителя в удобно воспринимать слагаемое ц,вр как плотность потока атомов, оседающих на поверхности. Это интегральный показатель, без разделения процесса на стадии. Вместо в можно написать 2в и интерпретировать в как долю адсорбируемых атомов Н. Замечание 1. Модель (3) быстрого растворения (локального равновесия) на поверхности получается из более общих соотношений
к-(Т)до(*) - к+(Т)во(1) = -Я(Г)сх(*, 0),
к-(Т)де(1) - к+_(Т)се(1) = Б(Т)Сх(М). Коэффициенты к , к+ характеризуют интенсивность процессов растворения в объеме и выхода на поверхность. Если эти процессы в рассматриваемом диапазоне температур существенно быстрее диффузии, то, полагая в относительном масштабе Осх ~ 0, получаем соотношение (3) с д = к-/к+. Если поверхность изотропна (в смысле Ек- ~ Ек+), то д(Т) слабо зависит от температуры. Формально можно записать аррениусовскую зависимость д(Т) = до ехр{—Ед/[КГ}], Ед = Ек- - Ек+, но «энергия активации» Ед не обязательно положительна. Дополнительно можно учитывать степени заполнения поверхности и насыщенности приповерхностного объема:
к~(Т)[1 - со ,,(^ко,е(*)-
Наличие «порогового» множителя [1 — c/cmax] приводит к следующему. Если концентрация c в приповерхностном объеме близка к максимально возможной, то растворение практически прекращается. Аналогично интерпретируется другой множитель, где величина 9(t) = q(t)/qmax означает степень заполнения поверхности. В балансовых уравнениях (4), (5) можно моделировать плотность потока адсорбции атомов H (диссоциативной хе-мосорбции водорода на поверхности) выражением ßs(T)p(t)[1 — e(t)}2. Но в диапазоне малых концентраций в ^ 1, что согласовано с квадратичной десорбцией, линейностью уравнения диффузии, D = D(c). В масштабе рассматриваемого далее ТДС-эксперимента удобно выбрать [p] = Topp, откуда ß(T) = (2nmkT)-1/2 и 2.484 ■ 1022/VT. Здесь [ß] = 1h2/(Toppсм2с), [T} = K, k - постоянная Больцмана, m - масса молекулы водорода. Зависимостью кинетической константы от температуры (ß а 1/ vT) часто пренебрегают на фоне экспоненты в формуле для s(T).
Замечание 2. Не для всех материалов целесообразно придавать особую роль поверхности. В случае высокой степени пористости («рыхлости») разумно уравнения (3)-(5) заменить граничными условиями III рода:
ßs(T)po,e(t) - b(T)c2 ,e(t) = T D(T)cx
10 , £ ■
(6)
- k+(T)[1 - qo,£(t)q—ix] C0,£(t) = T D(T)c
0,£
Плотность потока абсорбции ysp (в тонкий приповерхностный слой) пропорциональна (s) давлению молекулярного водорода снаружи. Формально (6) получается из уравнений (4), (5) при малой скорости накопления на поверхности (q & 0). Коэффициент десорбции обозначается одной буквой, хотя в (6) и (4), (5) это разные величины (включая размерность). Их согласованность понимаем в смысле bsurface = g2bvolume (см. (3)). В [3] приводится более обстоятельный анализ различных стадий проникновения водорода из газа в металл и обратно. С учетом итогового «усреднения», коэффициент b в условии (6) является эффективным коэффициентом рекомбинации.
Отметим следующее. Если с обеих сторон пластины поддерживать постоянное давление насыщения р (H2) при T = const, то через некоторое время установится равновесная концентрация с растворенного атомарного диффузионно подвижного водорода. Из модели (3)—(5), приравнивая производные к нулю, получаем с = Г^р, Г = g^/ys/b (с = gq). Таким образом, модель соответствует диапазону адекватности закона Сивертса (с rc^/p).
ТДС-эксперимент и модель измерений
В камеру c лентой из исследуемого металла или сплава подается водород в газовой фазе при сравнительно большом давлении. Лента нагрета с целью увеличения скорости сорбции. После того как образец поглотит достаточное количество водорода (до состояния равновесного насыщения), он быстро охлаждается (отключается ток нагрева). При этом резко падают скорости физико-химических процессов и значительное количество водорода остается в образце. В режиме последующего постоянного вакуумирования камеры лента снова нагревается. Закон нагрева T(t) может варьироваться в широких пределах. С помощью масс-спектрометра измеряется давление молекулярного водорода в вакуумной камере, обусловленное десорбцией J(t) = b(t)q2(t):
p(t) = 0i / * J (т) exp{(r - t)0-1} dr. (7) J 0
Здесь и далее для упрощения обозначений примем сокращенную запись (знак тождества трактуется как равенство по определению):
b(t) = b(T(t)), D(t) = D(T(t)), s(t) = s(T(t)).
Коэффициент g быстрого растворения (локального равновесия «поверхность - объем») считаем константой (E&- = ) в пределах рассматриваемого ТДС-пика десорбции. Для метода ТДС выполнена симметрия:
p(t) = po/(t), q(t) = qo/(t), co(t) = c*(t), (8)
D(t)cx(t, 0) = —D(t)cx(M), c(x) = с = const.
Константа 0i зависит от площади поверхности ленты S (01 = S02), 0o и 02 определяются конкретными характеристиками экспериментальной установки, в частности, объемом камеры V и скоростью откачки вакуумной системы v (0o = V/v). В дифференциальной форме
J (t) = [p(t)/0o + p(t)]/0i.
Время t* окончания эксперимента определим условием практически полной дегазации:
p(t) и 0, t ^ t*, c(t*,x) и 0, x е [0,4
Интегро-дифференциальное уравнение типа Риккати
Принятая модель ТДС-дегазации: ct = D(t)cxx , c(0, x) = с, co/t) = gq(t), q(t) = r(t) — b(t)q2(t) + D(t)cx(t, 0),
p(t) = 01 f J(r) exp {(r — t)0-1} dr, o
r(t) = ^s(t)p(t), J(t) = b(t)q2(t).
Вакуумную систему считаем достаточно мощной, чтобы после насыщения, на этапе дегазации пренебречь ресорбцией (г(£) = 0) в динамическом граничном условии (д = ...). Разрешимость краевой задачи (и в более общем случае учета обратимого захвата в объеме) обоснована в [16]. Для обратной задачи параметрической идентификации входными данными являются функции времени р(£) и 7 (*).
В данной статье ограничимся прямой задачей численного моделирования ТДС-спектра 7 = 7 (Т). Этот этап необходим как итерационная составляющая алгоритма идентификации. Нагрев Т(*) обычно реализуют линейным (Т(*) = То + V*). Скорость нагрева V невелика (< К/с). По достижении максимальной температуры (если дегазация еще не завершилась) нагрев прекращается: Т(*) = Ттах. Поскольку включение-выключение нагрева не происходит мгновенно, то, не выделяя явно эти относительно быстрые переходные процессы, можно считать в формальных математических выкладках коэффициенты В(Т (*)), Ь(Т(*)), в(Т(*)) гладкими по * е М.
Теперь мы в состоянии математически сформулировать задачу. В [16, 17] представлены разностные схемы решения краевой задачи (в том числе и с учетом ловушек различных типов). Но чтобы сравнивать модельные и экспериментальные ТДС-спектры, нужна только поверхностная концентрация (7 = Ьд2). Естественно попытаться избежать итерационного решения краевой задачи при текущих приближениях параметров модели Во, Ед, Ьо, Еь, во, Е8, д. С этой целью проведем преобразования, придя в итоге к необходимости интегрировать лишь систему ОДУ невысокого порядка.
Сделаем замену времени * = /0íDdв
(оставив прежнее обозначение *):
ег(*,х) = Схх(*,х), е(0,х) = с, ео,^ = дд(*), (9) Сх|о = -ех|^ = <?(*) + [7(*) - г(*)]В-1(*). (10)
Считаем д(*) функциональным параметром, а (10) - дополнительным соотношением к линейной задаче (9). Сделаем замену, приводящую краевые условия в (9) к однородным:
с = е(*, х) - дд(*), с4(*, х) = Схх(*, х) + /(*), /(*) = -дд(*), с(0,х) = 0(х) = 0, с|о,£ = 0.
Запишем решение с помощью функции мгновенного источника (функции Грина) [18, гл. 2]:
с(*,х)= / с1(х,е,*)^(о de+
о
+ / / Gi- т)f(t) d^dr, Gi(x,e,t) =
.7 0 .70
£ eXP
n2n2 л . nnx . nn£
n=1
i2 V sin"T sin" £
В динамические граничные условия входит производная cx(t, 0):
Сх|о = -д(т)^'ехр |(г - ¿)} йт,
^«=1,3,5...- При т = г ряд расходится, так что подразумевается почленное интегрирование. В исходном времени г имеем
Сх(г, о) = сх(г, 0) = ех(г,£) = Сх(г,£), ех(г, 0) =
4g ^ 'Г'
Jo
С П2П2 I л
д(т)ехр|--] йв^йт.
Окончательно динамическое граничное условие запишется в форме
д(г) = г (г) - ъ(г)д2(г)- (11)
г
Г ( П2П2 1
I ¿2 у0 д(т) ехр { —°(в) йв\йт-
Полученное уравнение (11) с квадратичной нелинейностью (которая входит и в функцию г(г) при «расшифровке» р(г)) будем классифицировать как интегро-дифференциальное уравнение Риккати нейтрального типа. Уравнение эквивалентно исходной краевой задаче в том смысле, что решение д(г) однозначно определяет решение е(г,х). Аналогия с функционально-дифференциальными уравнениями нейтрального типа [19] связана с тем, что избавиться от производной д в правой части интегрированием по частям нельзя, так как появится расходящийся ряд. Нас интересует отрезок времени [г1,г2] С (0,г*), соответствующий выраженному пику термодесорбции (измерения при г близких к 0, г* малоинформативны). Для мощной вакуумной системы ресорбцией пренебрегают (г(г) = 0), что упрощает уравнение. Уравнениям Риккати посвящена обширная литература (см. [20]).
Безразмерная форма задачи
Для численного моделирования удобно перейти к безразмерным переменным, выполнив замены: г' = ¡0в(в)йв/£2, X = х/1, V = д/д (с = дд). Не меняя обозначения г, получаем:
Ь(г) ^ ^ v(г) = -Ь(Ф2(г)- (12)
- 4д£^У V(т) ехр {-п2п2[г - т]}йт,
v(0) = 1 (начальные данные). Ограничимся первыми к слагаемыми ряда. Обозначим:
хг(г) = [ V(т) ехр { - (2г - 1)2п2[г - т]} йт. ■ 'о
Продифференцируем хг по г и подставим вместо V выражение (12) (сумма до п = 2к - 1):
Хг = -IV2 - [(2г - 1)2п2 + 4д£\х-- 4д£(х1 + ... + Хг-1 + Хг+1 + ... + хк).
В итоге получаем следующую систему обыкновенных дифференциальных уравнений:
d_ dt
(v
Zi
(1
\zk)
= -b(t)v2(t)
J)
- 4gi
11
a1 1
1
/zi\
Z2
V 1 1 ... aj \zkj
ai = 1 + (2i - 1)2n2/(4gi), v(0) = 1, Zi(0) = 0 (размерность матрицы (k + 1) x k). Система быстро интегрируется численно с помощью стандартного метода Рунге-Кутты 4-го порядка (авторы пользовались пакетом Scilab 5.5). Поскольку экспериментальные погрешности измерений не оцениваются ниже 10 % и в экспонентах n2, то достаточно ограничиться небольшим порядком системы (к + 1) по сравнению, например, с методом прямых.
Процедура выбора шага интегрирования не вызывает принципиальной трудности. В эксперименте фиксируется ТДС-спектр J (T) (T о t о t'). Представление ТДС-пика на равномерной сетке из 100 узлов достаточно (и даже избыточно). Из соображений численной аппроксимации с запасом можно брать шаг порядка секунды (даже при экспериментах на очень тонких мембранах за это время никаких заметных изменений не наблюдается). Практически такая «пристрелка» шага в физическом контексте не вызывает затруднений. Другое дело, что эта равномерная сетка при пересчете на безразмерную переменную t' = J^D(s)ds/i2 дает монотонно растущий шаг и экономию объема вычислений.
Этап начального насыщения
Остановимся на начальном насыщении водородом и коэффициенте при квадратичной нелинейности 6(t). Начальное насыщение проводится при относительно высоких температуре T = T = const и давлении напуска p = const
t
1
1
(для интенсификации сорбции). После установления равновесной концентрации имеем
^(T)s(T)p = b(T)c, с = gg (g = const),
с = gb- 1/2V^Sp Texp {E6[2RTjo1}
=> с =
b(t) = с b(t)£2D-1 = с b(t)£2[gD(t)]
1
Ьо ехр {-Еь[ЕТ(*)]-1}, Еь = Еь-Ед,
Ьо = £2В0-1ТЬ^ |техр {ЕЬ[2ЯТ]-1}.
Формально Ь(*) = Ь(Т(*)) - аррениусовский параметр с энергией активации Еь - Ед. Как правило, Еь > Ед и поверхностные процессы быстрее интенсифицируются с нагревом (Ь = ... ЕьТ, В = . ..ЕдТ), так что Еь > 0.
Предэкспонента Ьо зависит от всех параметров Во, Ьо, Еь, во, Е3 (кроме Ед). Тем самым коэффициент Ь в ТДС-эксперименте меняется строго монотонно и находится в определенном ограниченном диапазоне 0 < Ь- ^ Ь(*) ^ Ь+. Первое приближение
Качественную картину отражает уже первое приближение (у второго в сумме слагаемого на порядок больше п2 в экспоненте):
= -Ь(ф2со-
- к/ г)(т) ехр { -п2[* - т]}dr, к = 4д£, Jo
vV(t) = —b(t)v2(t) — Kz(t), Z = Z1, z(t) = —b(t)v2(t) — (к + п2 )z(t),
(13)
v(0) = 1, (0) = 0. Качественный характер движения на фазовой плоскости изобра-
жен на рис. 1. Пунктиром показаны параболы
v, z = 0: z = —b(t)x
o1v2,
—b(t)[K + п2]
2 o1 v2.
v ^ —bv2, v(0) = 1 ^ v(t) ^ 1 + b(r) dr
1
показывает, что за конечное время это невозможно. Фазовая траектория движения, выхо-
дя из точки (1,0), пересекает первую параболу и в дальнейшем, не покидая «клюв» между параболами, асимптотически стремится к началу координат при t ^ (формально T = const по достижении Tmax и дегазацию можно наблюдать бесконечно).
Аналогично можно рассматривать и второе приближение уравнения типа Риккати (12):
v(t) = -b(t)v2 - KZi(t) - KZ2(t), ¿1 (t) = -b(t)v2 - [к + n2]zi - KZ2, (14) ¿2 (t) = -b(t)v2 - KZ1 - [к + 9n2] ¿2.
Результаты моделирования
Для численных экспериментов были использованы данные по никелю и нержавеющей стали марки 12Х18Н10Т. Варьирование параметров позволяет оценить степень чувствительности («производные») ТДС-спектра к выделенным лимитирующим факторам.
Данные по никелю [5]: [E] = кДж/моль; D = 7.5 х 10-3 exp{-40/RT} [см2/с]; bvol = 1.53 х 10-18exp{-43.2/RT} [см4/с]; s = 1.8 х 10-2 exp{-61.4/RT}; р = 37.4 [Торр]; £ = 0.02 [cm]; T = 770 [K], T(0) = 295 [K]; T = 0.5 [K/c]. Для определенности зафиксируем b0 surf = 1.53 х 10-14 [см2/с]; g = 100 [1/с].
Данные по стали марки 12Х18Н10Т [5, 22]: D = 3.09 х 10-4 exp{-27.78/RT}; bvol = 5.05 х 10-13 exp{-97.14/RT}; s = 7.05 х 10-2 exp{-59.51/RT}; T(0) = 300. Для определенности b0surf = 5.05 х 10-9; £ = 0.1; g = 100; p = 100; T = 850, T = 0.5 [K/c].
При моделировании потока термодесорбции численно решалась система 15 уравнений. На рис. 2, 3 параметры перечисляются в порядке убывания максимумов кривых.
Поле фазовой скорости системы (13) нестационарно, но его изменение монотонно и ограниченно в силу Ь(*) е [Ь-, Ь+]. Параболы несколько сдвигаются влево (Еь > Ед) или вправо (Еь < Ед). Качественная структура поля остается неизменной. Обоснуем неограниченную продолжимость решения v(t), ¿(¿). В соответствии с общей теорией [21, с. 110] решение либо продолжимо вперед неограниченно, либо достигнет границы треугольника (компакта). Структура поля скоростей такова, что гипотетически можно достичь только начала координат. Но оценка во время движения
Рис. 1. Движение на фазовой плоскости (v,z)
47
t
o
ЫГ(г))
800 1000
Рис. 2. Влияние энергий активации на 6(4), N1
Рис. 4. Первое приближение (13), №
ЫТ(г))
200 400 600 800 1 ООО 1200
300 400 500 600 700 800 900
Рис. 3. Влияние энергий активации на 6(4)
\ \ \
\\ ^ = о \ \ V = о
\ ^Г = 450, 500, 600
Рис. 5. Первое приближение (13), сталь
V
V
Ъ0 = 1.53 -10-1; Ъ0 = 1.53 -10-1'
------- Ъ0 = 1.53-10-1
Рис. 8. Система (14), влияние Ьо, N1
72 /Г _1
1 Г -—- \ 1 *
\ 1 1 4 /У \
\ ^_ _ /
уьт^еГЗ-''''''^~~ ~ ~~ ~
^ 12Х18Н10Т/
0 о5-„V-—_ V 06 0.8 :
-------Ь0 = 5.05 .10-1°
- Ь0 = 5.05 -10-'
Ь0 = 5.05 -10-8
Рис. 9. Система (14), влияние Ьо, сталь
1
--------------г=101
- г = 102
г = 103
Рис. 12. Система (14), влияние д, №
-------ш = 101
« = 102
Рис. 13. Система (14), влияние д, сталь
Рис. 14. ТДС-спектр, влияние Ьо, № Рис. 15. ТДС-спектр, влияние Ьо, сталь
Рис. 16. ТДС-спектр, влияние Еь, № Рис. 17. ТДС-спектр, влияние Еь, сталь
Рис. 18. ТДС-спектр, влияние д, № Рис. 19. ТДС-спектр, влияние д, сталь
Заключение
В статье представлена краевая задача с нелинейными динамическим граничными условиями для моделирования ТДС-спектра дегазации конструкционного материала, предварительно насыщенного водородом, в условиях вакуумирования и монотонного нагрева. Решение прямой задачи необходимо для качественного сопоставления с экспериментальными данными и как итерационная составляющая алгоритма параметрической идентификации. Для моделирования термодесорбцион-ного потока предложен вычислительный алгоритм, требующий (вместо приближенного решения распределенной краевой задачи с текущими приближениями параметров) лишь интегрирования нелинейной системы обыкновенных дифференциальных уравнений невысокого порядка (по сравнению, например, с методом прямых). Приведены результаты вычислительных экспериментов, использующих данные по никелю и конструкционной стали.
Работа выполнена при поддержке РФФИ (грант № 15-01-00744).
Литература
1. Водород в металлах / Ред. Г. Алефельд, И. Фелькль. М.: Мир, 1981. Т. 1,506 c. Т. 2,430 с.
2. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.
3. Писарев А. А., Цветков И. В., Марен-ков Е.Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
4. Черданцев Ю. П., Чернов И. П., Тюрин Ю. И. Методы исследования систем металл-водород. Томск: ТПУ, 2008. 286 с.
5. Изотопы водорода. Фундаментальные и прикладные исследования / Ред. А. А. Юхим-чук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697 с.
6. Varin R.A., Czujko T., Wronski Z.S. Nano-materials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi: 10.1007/9780-387-77712-2.
7. The hydrogen economy / Eds. M. Ball, M. Wietschel. Cambridge University Press, 2009. 646 p.
8. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
9. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids // Technical Physics. 1999. Vol. 44, N 1. P. 90-94. doi:10.1134/1.1259257
10. 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. P. 5818-5851. doi:10.1016/j.ijhydene.2014.01.158
11. 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
12. 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. P. 9060-9069. doi:10.1016/j.ijhydene.2010.05.092
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. Elsevier. 2009. Vol. 33, N 10. P. 37763791. doi:10.1016/j.apm.2008.12.018
14. Zaika Yu. V., Bormatova E. P. Parametric identification of a hydrogen permeability model by delay times and conjugate equations // International Journal of Hydrogen Energy. 2011. Vol. 36, N 1. P. 1295-1305. doi:10.1016/j.ijhydene.2010.07.099
15. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with dynamical conditions on surface // Mathematical Modelling. Nova Sci. Publishers. 2013. P. 269-302.
16. Заика Ю. В. Интегральные операторы прогнозирования и идентификация моделей водородопроницаемости. Петрозаводск: КарНЦ РАН, 2013. 505 c.
17. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption // Advances in Materials Science and Applications. World Acad. Publ. 2014. Vol. 3, Iss. 3. P. 120-129. doi:10.5963/AMSA0303003
18. Мартинсон Л. К, Малов Ю. И. Дифференциальные уравнения математической физики. М.: МГТУ, 2002. 368 c.
19. Хейл Дж. Теория функционально-дифференциальных уравнений. М.: Мир, 1984. 421 с.
20. Егоров А. И. Уравнения Риккати. М.: Физ-матлит, 2001. 320 с.
21. Арнольд В. И. Обыкновенные дифференциальные уравнения. Ижевск: РХД, 2000. 368 с.
22. Popov V. V., Denisov E.A. Inhibition of hydrogen permeability by TiN: evaluation of kinetic parameters / Hydrogen Materials Science and Chemistry of Carbon Nanomaterials / Eds. T. N. Veziroglu et al. Springer. 2007. P. 671-680. doi:10.1007/978-1-4020-5514-0
Поступила в редакцию 06.04.2015
51
References
1. Vodorod v metallakh [Hydrogen in metals] / Eds. G. Alefeld, J. Volkl. Moscow: Mir, 1981. Vol. 1, 506 p. Vol. 2, 430 p.
2. Vzaimodeistvie vodoroda s metallami [Interaction of hydrogen with metals] / Ed. A. P. Zakharov. Moscow: Nauka, 1987. 296 p.
3. Pisarev A. A., Tsvetkov I. V., Marenkov E. D, Yarko S. S. Pronitsaemost' vodoroda cherez metally [Hydrogen permeability through metals]. Moscow: MIFI, 2008. 144 p.
4. Cherdantsev Yu. P., Chernov I. P., Tyurin Yu. I. Metody issledovaniya sistem metall-vodorod [Methods of studying metal-hydrogen systems]. Tomsk: TPU, 2008. 286 p.
5. Izotopy vodoroda. Fundamental'nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies] / Ed. A. A. Yukhimchuk. Sarov: RFYaTs-VNIIEF, 2009. 697 p.
6. Varin R. A., Czujko T., Wronski Z. S. Nano-materials for solid state hydrogen storage. Springer, New York, 2009. 338 p. doi:10.1007/978-0-387-77712-2
7. The hydrogen economy / Eds. M. Ball, M. Wietschel. Cambridge University Press, 2009. 646 p.
8. Handbook of hydrogen storage: new materials for future energy storage / Ed. M. Hirscher. Wiley-VCH, 2010. 353 p.
9. Gabis I. E. The method of concentration pulses for studying hydrogen transport in solids. Technical Physics. 1999. Vol. 44, N 1. P. 90-94. doi:10.1134/1.1259257
10. 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. P. 5818-5851. doi:10.1016/j.ijhydene.2014.01.158
11. 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
12. 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. P. 9060-9069. doi:10.1016/j.ijhydene.2010.05.092
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. Elsevier. 2009. Vol. 33, N 10. P. 37763791. doi:10.1016/j.apm.2008.12.018
14. Zaika Yu. V., Bormatova E. P. Parametric identification of a hydrogen permeability model by delay times and conjugate equations. International Journal of Hydrogen Energy. 2011. Vol. 36, N 1. P. 1295-1305. doi:10.1016/j.ijhydene.2010.07.099
15. Zaika Yu. V., Rodchenkova N. I. Hydrogen-solid boundary-value problems with dynamical conditions on surface. Mathematical Modelling. Nova Sci. Publishers. 2013. P. 269-302.
16. Zaika Yu. V. Integral'nye operatory prognozirovaniya i identifikatsiya modelei vodorodopronitsaemosti [Integral forecasting operators and identification of hydrogen permeability models]. Petrozavodsk: KarRC of RAS, 2013. 505 p.
17. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption. Advances in Materials Science and Applications. World Acad. Publ. 2014. Vol. 3, Iss. 3. P. 120-129. doi:10.5963/AMSA0303003
18. Martinson L. K., Malov Yu. I. Differentsial'nye uravneniya matematicheskoi fiziki [Differential equations of mathematical physics]. Moscow: MGTU, 2002. 368 p.
19. Kheil Dzh. Teoriya funktsional'no-differentsial'nykh uravnenii [Theory of functional-differential equations]. Moscow: Mir, 1984. 421 p.
20. Egorov A.I. Uravneniya Rikkati [Riccati equation]. Moscow: Fizmatlit, 2001. 320 p.
21. Arnol'd V.I. Obyknovennye differentsial'nye uravneniya [Ordinary differential equations]. Izhevsk: RKhD, 2000. 368 p.
22. Popov V. V., Denisov E.A. Inhibition of hydrogen permeability by TiN: evaluation of kinetic parameters. Hydrogen Materials Science and Chemistry of Carbon Nanomaterials / Eds. T. N. Veziroglu et al. Springer. 2007. P. 671-680. doi:10.1007/978-1-4020-5514-0
Received April 06, 2015
СВЕДЕНИЯ ОБ АВТОРАХ:
Заика Юрий Васильевич
зав. лаб. моделирования природно-технических систем, д. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (8142) 766312
Костикова Екатерина Константиновна
младший научный сотрудник, к. ф.-м. н. Институт прикладных математических исследований Карельского научного центра РАН ул. Пушкинская, 11, Петрозаводск, Республика Карелия, Россия, 185910 эл. почта: [email protected] тел.: (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: [email protected] tel.: (8142) 766312
Kostikova, Ekaterina
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: [email protected] tel.: (8142) 766312
