CC BY
10
Труды Карельского научного центра РАН
№7. 2020. С. 46-56
DOI: 10.17076/mat1246
УДК 519.87:539.2
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПИКОВ ТЕРМОДЕСОРБЦИИ ВОДОРОДА
Ю. В. Заика, Е. К. Костикова
Институт прикладных математических исследований КарНЦ РАН, ФИЦ «Карельский научный центр РАН», Петрозаводск, Россия
При исследовании спектров термодесорбции (ТДС-спектров) изотопов водорода используются различные математические модели, в частности, в форме реакции первого порядка для усредненной по объему концентрации (в случае лимитирования диффузией), второго порядка (с учетом рекомбинации атомов водорода на поверхности в молекулы) и более детализированные модели, отражающие многостадийность переноса (диффузия в объеме, выход на поверхность, рекомбинация, десорбция). Интерес представляет анализ структуры спектра с целью выявления соответствующих каждому пику причин и «движущих сил» физико-химического характера. В статье приведены анализ симметричности унимодальных ТДС-спектров и результаты численного моделирования.
Ключевые слова: термодесорбция водорода; численное моделирование.
Yu. V. Zaika, E. K. Kostikova. NUMERICAL MODELING OF HYDROGEN THERMAL DESORPTION PEAKS
Various mathematical models are used for investigate the thermal desorption spectra (TDS spectra) of hydrogen isotopes, in particular, in the form of a firstorder reaction for the volume-averaged concentration (in the case of limitation by diffusion), second-order models (taking into account the recombination of hydrogen atoms on the surface into molecules) and more detailed models reflecting multistage transport (diffusion in the buk, mifration to the surface, recombination, desorption). It is interesting to analyze the structure of the spectrum in order to identify the causes and «driving forces» of a physical and chemical nature corresponding to each peak. The article presents an analysis of the symmetry of unimodal TDS spectra and the results of simulation.
Keywords: hydrogen thermodesorption; numerical modeling.
Введение
Интерес к взаимодействию изотопов водорода с различными материалами носит многоплановый характер [1—5]: защита конструкционных материалов от водородной коррозии, транспортировка углеводородного сырья, ракетостроение, водородная энергетика, перспективы термоядерного синтеза.
Одним из наиболее информативных методов исследования кинетики взаимодействия новых материалов с водородом является тер-модесорбционная спектрометрия (ТДС), позволяющая «сканировать» материал в нестационарных условиях в широком диапазоне рабочих температур. Математические модели термодесорбции и водородопроницаемости с уче-
том различных стадий переноса и численные методы решения краевых задач, разрабатываемые авторами статьи, изложены в [11, 15-18].
Трудности анализа спектров типичны для некорректных обратных задач математической физики: неединственность решения, неоднозначность интерпретации и высокая чувствительность к погрешностям входных данных. На экспериментальном спектре наблюдается результат наложения множества процессов, подлежащих идентификации. Сначала необходима идентификация по физико-химическим процессам, а затем уже параметрическая идентификация в соответствии с принятой физически обоснованной моделью.
Статья содержит некоторые результаты численного моделирования, ориентированные на экспериментальные данные и теоретические положения кинетики взаимодействия изотопов водорода с различными материалами, изложенные в работах [4, 7-10, 12-14].
МОДЕЛЬ РЕАКЦИИ ПОРЯДКА а Е [1, 2]
В качестве отправной точки рассуждений рассмотрим уравнение реакции первого порядка X(t) = -K(T)X(t), где X(t) - текущая усредненная концентрация водорода в образце, K(T) = K0 exp{—Q/RT} (аррениусовская зависимость кинетического коэффициента от температуры), T(t) = T0 + ßt (равномерный линейный нагрев с относительно невысокой скоростью ß, [ß] = K/s).
Для полноты изложения, чтобы пояснить смысл обобщения, приведем кратко схему вывода уравнения, когда десорбция лимитирована диффузией. Рассмотрим краевую задачу ct = D(T(t))cxx [t> 0, x Е (0, h)], c = co > 0 [t = 0], c = 0 [t > 0, x = 0,h]. Для растворенного водорода речь идет об атомах, [c] = 1/см3. В начальный момент диффузант распределен равномерно. Затем вследствие ваку-умирования мгновенно (в относительном масштабе времени) устанавливаются нулевые концентрации в подповерхностном объеме пластины. После перехода к безразмерным переменным t' = /qD(T(r)) dr/h2, x' = x/h, c' = c/co получаем формально D = 1, h = 1, c0 = 1 (в текущих промежуточных выкладках оставляем прежние обозначения, не усложняя их штрихами). Тогда для усредненной переменной c(t) = f0 c(t, x) dx (см., например, [6, с. 27]) справедливо представление
Е
«.= 1,3,5...
8
-Г~2 exP i--
п2п2 У r
n2t
1
r
п
2
При анализе ТДС-спектров нас интересует производная dc/dt, которая при t = 0 не определяется почленным дифференцированием ряда. Это естественное следствие несогласованности при t = 0 начального и граничного условий. Строго говоря, ряд представляет так называемое обобщенное решение краевой задачи. Спустя непродолжительное время можно ограничиться первым слагаемым: c(t) & 8exp{-t/r}/п2 (t > t0 > 0), где параметр r приобретает смысл «времени релаксации» (уменьшения в e раз). Чтобы компенсировать отброшенные слагаемые и ретроспективно «вписаться» в начальные данные c(0) = c0 = 1, увеличим коэффициент с 8/п2 до 1. Теперь для аппроксимации X(t) & c(t) (t ^ 0) примем модель X(t) = X(t)/r, X(0) = X0 = 1, и возвратимся к исходным переменным (t,x):
X(t) = -K(T)X(t), X0 = 1, K = n2h-2D(T). Здесь мы оставили нормировку c(t) = j'0h c(t,x) dx/h ^ c(t)/c0, т. е. теперь X(t) Е (0,1) — усредненная доля (от исходной равномерной концентрации водорода c0), оставшаяся в образце к моменту времени t > 0. Модель работоспособна, когда в приложениях нас интересует интегральный поток десорбции без детализации последовательности физико-химических процессов. Интерпретация K(T) позволяет для пористых и порошкообразных материалов говорить обобщенно об эффективных характеристиках переноса: коэффициенте диффузии D и характерном диффузионном пробеге h. В принципе можно не связывать K(T) с диффузией, ограничившись выводом уравнения на уровне «скорость десорбции пропорциональна текущей усредненной концентрации». В любом случае постулируется аррениусовская зависимость K (T) = K0 exp{-Q/RT} с предэкспонентой (частотным множителем K0) и энергией активации Q.
С другой стороны, в теории адсорбции-десорбции используется уравнение Поляни -Вигнера 0(t) = -M(T)&n(t), где в — степень заполнения поверхности. В частности, если речь идет о термодесорбции диссоциативно хемосорбированного двухатомного газа, то n = 2. Модель можно адаптировать и для усредненной по объему концентрации, используя понятие эффективного коэффициента рекомбинации [5]. Для этого уравнение dc/dt = -b(T)c2 разделим на c0 и получим X = —K (T)X2, где K = bc0. Концентрация c0 определяется начальным насыщением образца. При этом с ростом давления насыщения из газовой фазы температура максимума модельного потока десорбции будет уменьшаться.
47
Преимущественно адсорбция водорода характерна при относительно низких температурах, а описанное выше лимитирование диффузией - при высоких. В реальности диффузия и десорбция с поверхности взаимосвязаны (динамика процессов на поверхности диктует граничные условия для уравнения диффузии). Такие более детализированные модели уже требуют разработки специализированного математического обеспечения [16, 18].
Синтезируя приведенные рассуждения, рассмотрим усредненную модель (для равномерного нагрева):
г!Х
_ = -в-1К(Т)Ха(Т), X(То) = 1, а е [1, 2],
К(Т)= Ко ехр{^[ЯТ]-1}, Т е [То,Т],
т (г) = То + вг, йТ = вйг.
Здесь Т0 — начальная температура (обычно комнатная), когда десорбция водорода практически отсутствует; Т* — конечная температура (когда десорбция уже пренебрежимо мала на фоне максимума потока); безразмерная переменная X(Т) имеет смысл доли усредненной по объему образца концентрации с от концентрации равномерного начального насыщения со. В силу г о Т можно записать X(Т) = с(Т)/со. Кинетический коэффициент К(Т) ([К] = 1/в) заранее не связываем явной формулой с коэффициентом диффузии Б(Т) или коэффициентом рекомбинации Ь(Т). Если предполагается обрабатывать данные с различными условиями равномерного насыщения, то целесообразно явно выделить зависимость от со в форме К = КС(Т)Са-1.
Параметр а е [1,2} позволяет учитывать степень участия в процессе насыщения ассоциативной хемосорбции, растворения в объеме с последующей рекомбинацией атомов водорода в молекулы при термодесорбции. Иными словами, применяем усреднение не только по концентрации, но также по лимитирующим (во взаимодействии) процессам диффузии и рекомбинации. Что касается самой «дробности», то в математической физике теория уравнений диффузии с дробными производными разработана. Усредненные по концентрации модели (если нас интересует общий поток без детализации) адекватно описывают основную часть изолированного всплеска на ТДС-спектре (когда росту К(Т) активно «противодействует» уменьшение X(Т), что и формирует пик). Начальные и конечные участки графика термодесорбции (когда поток относительно мал и измеряется с меньшей точностью) малоинформативны.
СВОЙСТВА МОДЕЛЬНОГО ТДС-ПИКА
Унимодальность и метод Киссинджера
Нас интересует зависимость ш(Т) = -йХ/йТ = К (Т )Ха (Т )/в — нормированный поток термодесорбции (в долях X = с/со на градус температуры). Знак тождества (=) часто используем по контексту в смысле равенства по определению. Если оперируем временем, то рассматриваем у(г) = -^С(г) = К(Т(г)^а(г). При этом ш(Т) = у(Т)/в с учетом г о Т = То + рг. График ш(Т) представляет ТДС-спектр. Поскольку используется (безразмерная) нормировка X(То) = 1, то
должно быть 5 = /Т* ш(Т) йТ = 1. Строго говоря, нужно формально интегрировать на (0, +го), но мы рассматриваем отрезок [То,Т*}, вне которого поток пренебрежимо мал. Итак, модельный ТДС-пик нормирован по площади. Регистрируемый поток (например, в единицах шо1/еш3в) нормируется интегралом (который равен со). Говоря о потоках, слово «плотность» обычно опускаем по контексту.
Интегрируя уравнение йX/йT = —KXa/в в квадратурах, получаем: ш(Т) = -сШ/йТ,
ш(Т )=<
в-1К(Т)ехр\ - в-11Т К(Т) йТ к
в-1К(Т) 1 + (а -1) в-1!ТоКйТ
где первая формула соответствует а = 1, а вторая - а е (1,2]. Представления согласованы: правосторонний предельный переход а ^ +1 во второй формуле дает первую. Начальная температура То относительно низкая (обычно комнатная, когда К (То) ^ 1 и поток еще пренебрежимо мал). С ростом температуры аррениусовский коэффициент К (Т) монотонно выходит на горизонтальную асимптоту, а множитель с интегралом монотонно убывает. Это дает колоколообразный график ш(Т).
Далее производную по Т (в отличие от точки сверху для производной по времени г) будем обозначать штрихом. Приведем результаты вычислений (с учетом К1 = KQ/RT2):
ш'(Т) = ш(Т)р(Т), а е [1,2], ш"(Т ) = ш(Т )[^2 (Т) + <р'(Т)],
ш"'(Т) = ш(Т) &3(Т) + 3<р(Т)<р'(Т) + <р"(Т)]. Здесь приняты обозначения
№) = ^ - ав-1 К(Т)К-1(Т),
К(Т; а) =
гТ
1 + (а - 1) в-1 К(Т) йТ1
^ То
48
1 -а
В контексте обратной задачи параметрической идентификации отметим, что
К(Т; а) = X 1-а(Т), К(Т;1) = 1.
Для поиска температуры Тт = Ттах экстремума (который является единственным, глобальным максимумом) запишем уравнение ™'(Тт) = 0, откуда ^ = 0 ^ < = 0,
K(Tm) = О?- K(Tm; а)
aim.
z =
Q
RTm
(1)
Выделение безразмерной переменной г в выкладках целесообразно, поскольку дробь вида —Е/ЕТ фигурирует в экспоненте закона Ар-рениуса. В дальнейшем рассмотрении потребуются следующие производные:
■ы"{Тт) = <ш(Тт)<р'(Тт)
= ■ю(Тт)Тт-2г [ — а-1 г — 2] < 0, (2)
w'"(Tm) = w(Tm)^"(Tm)
(3)
Подчеркнем, что зависимость Тт(Коа, в) определяется алгоритмически численным решением уравнения (4), т. е. у переменной г = Q/ЕTm числитель и знаменатель не являются независимыми при фиксированных Ко.
На рисунках 1-4 показаны возможные значения параметра Ко в зависимости от скорости нагрева в для широких диапазонов пиковых температур Треак-Треак (согласно данным [7]) и энергий активации основных диффузионных и хемосорбционных процессов в графите QI-QIV (в соответствии с табл. 1 [4]).
q, kJ/mol
In (V)
= w(Tm)T-iг [а-2(а — 2) г2 + 2г + 6.
Рассмотрим сначала (1) как уравнение для Т = Тт. Перепишем его в форме
f(Т) = аТ2К — (а — 1)Е\ТтК<1Т = Е. (4)
Геометрически требуется найти пересечение графика функции f (Т) с горизонтальной прямой, определяемой значением вQ/Е. Поскольку из соотношений К(Т) = К0 exp{—Q/ЕT}, К'(Т) = К(Т)Q/[ЕT2] следует неравенство f'(Т) = К(Т)[2аТ + Q/Е] > 0 (положительна и вторая производная), то решение единственно. Существование обеспечивается условием К (То) ^ 1 (формально можно рассматривать Т0 = 0 и Т > 0). Увеличение в (правой части уравнения (4)) смещает температуру максимума спектра (Тт) вправо.
Уточним априорное ограничение на выбор предэкспоненты (частотного множителя) Ко, диктуемое уравнением (4). Значение То определено как температура (обычно комнатная), при которой поток еще незначителен на фоне предстоящего максимума при Т = Тт > То. Поэтому в рамках модели ТДС-пика должно быть выполнено неравенство f (То) < вQ/Е, т. е. Ко < вQ ехр{Q/ЕTо}/[аЕТ2].
На практике (при а = 1) для оценки параметров Ко, Q используют зависимость регистрируемого значения Тт от скорости нагрева в. После чего полученные оценки сопоставляются с экспертными данными о параметрах возможных физико-химических процессов. Скорости нагрева известны, а Тт определяются по экспериментальному графику.
600 г'тоо
soor* 1000
Qi
i юг
Рис. 1. Значения \i\[K0/ß]: — Tm, Q & Qj Fig. 1. Values \n[K0/ß\: - Tm, Q & Qj
q, -1=2 kl
IiW)
160-
500 tek700
900 г
1 3001- 3k 1 450
t,J/o
Рис. 2. Значения \n[K0/ß]: - Tm, Q & Qjj Fig. 2. Values \n[K0/ß]: - Tm, Q & Qjj
Классического метода Киссинджера и приближения ТДС-1 (а = 1) не всегда достаточно для удовлетворительной аппроксимации. В статье [7] на рис. 3 представлены пять ТДС-спектров дейтерия, полученных при различных скоростях нагрева графита 180-880.
49
к
Ii
Q, kJ/mol
In W)
800 Гр?,к 1 000
1 200 tpl, 1 400
лтт, k
1 600TPik 1700
Рис. 3. Значения \n[K0/в]: — Tm, Q « Qm Fig. 3. Values \п[К0/в]: - Tm, Q « Qm
Q, kJ/mol
liWfl
Tm, K
1 100 1 200 1 зоотр^ 1 400 1 500 1 600TPik 1 700
Рис. 4. Значения \п[К0/в]: — Tm, Q « QIV Fig. 4. Values \n[K0/p]: — Tm, Q « Qiv
Спектры имеют по 4 вершины (пики пронумерованы по возрастанию температуры). Авторы, используя метод Киссинджера, оценивают энергию активации высокотемпературных пиков 3 и 4 по сдвигу пиковой температуры при варьировании скорости нагрева. По мнению авторов, пик 3, вероятно, лимитируется диффузией, а пик 4 — процессами высвобождения из дефектов. Но при наложении на экспериментальный график пиков с энергиями активации, полученными методом Киссинджера, видно существенное расхождение: теоретические пики оказываются заметно уже экспериментальных. Кроме того, отмечается, что полученная методом Киссинджера энергия активации пика 3 (3.48 эВ) довольно высока, чтобы быть энергией активации диффузии. Поэтому на этот пик, по мнению авторов, оказывают сильное влияние процессы захвата
и высвобождения ловушками. Для пика 4 та же ситуация: теоретическая кривая (по Киссинджеру) не воспроизводит экспериментальные спектры. Авторы предполагают, что широкий высокотемпературный пик 4 обусловлен кластерами дефектов сложной структуры. А потому пик 4 может разлагаться на сумму (до 10!?) пиков десорбции из ловушек с различными энергиями активации.
Выполнен расчет энергий активации высокотемпературных пиков 3-4 (в обозначениях [7]) с использованием аналога метода Киссинджера, но в предположении, что пики обусловлены реакциями второго порядка (а = 2). Результаты представлены на рис. 5. С учетом К-1 (Тт; а) = Xа-1(Тт) «соотношение Киссинджера» для различных а записывается как
\n
в1/о T 2
1
aRT„
Q+\n
Q ■ (Kowa-\Tm)f/a
Это соотношение следует из подстановки в формулу (1) усредненной текущей концентрации, определяемой выражением X = (ш/[К/в})1/а, где нормированный поток ш известен по результатам эксперимента. В координатах { — 1/[аRTm}, чение энергии активации есть угловой коэффициент линейной регрессии. Как показали численные эксперименты, при расстановке пиковых температур для различных скоростей нагрева (данные взяты с рис.3 [7]) в координатах «типа Киссинджера» точки группируются близко к линии регрессии и при а = 1, и при а = 2 (промежуточные значения дают такую же картину). Полученные значения Q несколько меньше, чем в [7] (при а = 1).
-0.' -15.0
055
-0.05
-0.045
-0.04
-0.035
-0.03
-15.5-
-16.0-
-17.0-
-17.5.
Peak 3 / x = l-1/(2RTm)
(~1300KV
t Peak 4 /
300.08 kJ/mol/ (~1600KV
(3.11 eV) P 4
J /532.98 kJ/mol
y = 300.08x-1.818 f (5.53 eV)
/ y = o32.98x+3.6656
У = in
Рис. 5. Определение энергии активации, а = 2 Fig. 5. Determination of activation energy, а = 2
50
ш
Q
IV
Анализ симметрии ТДС-пика
Одной из общепринятых методик является разложение спектра на сумму гауссовских кривых (например, в пакете прикладных программ Origin). Более адекватным является разложение на сумму ТДС-пиков, генерируемых реакциями 1-2 порядков. Но они, строго говоря, несимметричны. Возникает задача поиска критериев симметричности и диапазонов параметров, где отклонения от гауссианов практически несущественны.
Качественное отличие ТДС-пика (а Е [1, 2]) от гауссовской кривой вида
G(T) =
A
V2n
ехр i -
а
(T - Tm)2 2а2
лишь в возможной существенной (в широком диапазоне параметров) несимметричности относительно вертикали Т = Тт. В этом случае унимодальный ТДС-спектр хорошо аппроксимируется склейкой двух гауссианов
Gi,2(T) = w(Tm)exp{ -
(T - Tm)2 2а2
T ^ T„
Следуя «правилу двух сигм», такую аппроксимацию достаточно рассматривать на отрезке [Т-,Т+] = [Тт—2аъТт+2а2]П[То,Т]. Практический интерес представляет «яркая» часть спектра, начальный и конечный этапы (когда поток незначителен на фоне максимума) малоинформативны.
На рис.6 представлено отношение а2/а1 для модельных спектров при различных а Е [1, 2] и энергиях активации QI-QIV, характерных для процессов хемосорбции и диффузии в углеродных материалах (см. табл. 1 [4]).
Численные эксперименты в широком диапазоне {Ко[4] показывают, что, за редким исключением, модельный спектр несимметричен. Это соответствует экспериментальным данным. Аппроксимация спектров гаусси-анами позволяет получать приближения параметров с последующим уточнением в рамках рассматриваемых моделей. Пик, соответствующий реакции первого порядка (а = 1), может быть практически симметричным. Возможен вариант «относительно крутой восходящий фронт и пологий нисходящий», но обычно
igM
ауау, Ql - 20 kJ/mol
ig [КсМ 16.0
14.0 -
12.0 -
10.0-
oyby, Qu - 120 kJ/mol
12.0-
а ю.о
Рис. 6. Отношение ширины нисходящего фронта к ширине восходящего Fig. 6. The ratio of the width of the descending front to the width of the ascending
пологий восходящий фронт сменяется крутым нисходящим. Для спектра реакции второго порядка (а = 2) характерен относительно крутой восходящий фронт и пологий нисходящий.
Перейдем к критерию симметричности графика и(Т) относительно вертикали Т = Тт. Для этого воспользуемся разложением
и(Тт + А) = ы(Тт) + и'(Тт)А + 2 и",А2
+ 1 и'"(Тт)А3 + 1 и(4)(Тт)А4 + о(А4). Ь 24
Численные эксперименты подтверждают, что в физически реальном диапазоне Т е [Тт-,ТГ+] пятую производную уже можно не учитывать. В экстремуме и'(Тт) = 0, а слагаемые с А2 и А4 дадут те же значения при замене А на —А. Следовательно, для оценки симметричности нужно оценить и''' (Тт) в сравнении с и''(Тт). Критерий симметричности (с учетом (1), (2))
А
3
\w"(Tm)\ » V\w'"(Tm)\, А е [-2а!, 2а2],
Tm = Tm( Ко, Q; a,ß), ai = ai (Tm,Q; a,ß),
преобразуется к виду (z = Q/RTm)
\А\ \6 + 2z - a-2(2 - a) z
3 Tm
2 + z
« 1. (5)
Строго говоря, спектр всегда несимметричен. Поскольку третья производная и'''(Тт) может менять знак (в широком диапазоне параметров), то более пологой может оказаться как восходящая, так и нисходящая ветвь.
При фиксированном Q спектр второго порядка шире, значения и''(Тт) и и'''(Тт) больше при а = 2, чем при а = 1. Варьирование порядка реакции а е [1,2] сопровождается сдвигом пиковой температуры. Изолированный пик, в котором пологий восходящий фронт сменяется резким нисходящим, не следует связывать с реакцией второго порядка.
Проиллюстрируем предложенный критерий на рис. 7. Среди представленных допустимых значений а и г только малая часть точек оказывается на «высоком хребте», где соотношение параметров дает практически симметричный спектр. Отметим, что при а = 2 заведомо и'''(Тт) = и(Тт)(2г + Ь)г/Тт > 0, что объясняет выводы численных экспериментов о несимметричности ТДС-спектра реакции второго порядка (восходящий фронт ниже). Вопрос лишь в том, насколько это отклонение существенно в конкретной задаче. Для а < 2 спектры, близкие к симметричным, получаются при и'''(Тт) к 0, откуда г(а) к а2[1 + л/1 + Ьа~2(2 — а)]/(2 — а). Формально г(а) ^ при а ^ 2—, но гипотетическое превышение порога г > 20 — 30 требует дополнительного физического обоснования.
В частности, для реакции первого порядка г к Q/[RTm] = 1 + \[Ч (а = 1), что соответствует относительно низким значениям энергии активации Q. Даже если выполняется равенство и"'(Тт) = 0 (г2 — 2г — Ь = 0), то при соответствующем значении г выполняется и(5\Тт) = и(Тт)гТ-5р4(г), р4(г) = (г2 — 2г — Ь)(9г2 + 20г — 42) — Ь0г — 132 < 0.
Рис. 7. Критерий (5) симметричности пиков Fig. 7. Criterion (5) of peaks symmetry
МОДЕЛЬ С ДИНАМИЧЕСКИМИ ГРАНИЧНЫМИ условиями
Выше рассмотрена модель динамики десорбции в терминах усредненной по объему образца концентрации. Рассмотрим более детализированную модель, явно разделяя объемные и поверхностные процессы (следуя работе [1, с. 177-206]). Ресорбцией пренебрегаем. Для тонкой пластины толщины £ краевая задача ТДС-дегазации примет следующий вид:
ct(t, x) = D(T)cxx, t e (0, t), x e (0, £), c(0,x) = co, x e [0,£], co/(t) = g(T)q(t), q(t) = -b(T)q2(t)+ D(T)cx(t, 0), J(T) = b(T)q2(t), T(t) = To + ßt, ß> 0.
Здесь c(t, x) — концентрация растворенного атомарного водорода; q(t) — поверхностная концентрация; D, b, g — (аррениусовские) коэффициенты диффузии, десорбции, быстрого растворения (квазиравновесия поверхностной и приповерхностной объемной концентрации); J(t) = J(T(t)) — плотность десорбции (атомов, рекомбинировавших в молекулы).
52
2
Модель позволяет на качественном уровне воспроизвести кривые, подобные экспериментальным пикам 1-4 (см. Fig. 5 в [8]). Результаты численного моделирования представлены на рис. 8 и 9, в том числе в логарифмической шкале. Фиксированы следующие значения параметров: t = 0.1 см, b0 = 103 см2/с, Eb = 250 кДж/моль, Do = 0.1 см2/с, ED = 125, go = 104 см"1, Eg = 0, Co = 4.8 х 1018 см"3, T0 = 300 К. Скорость нагрева варьируется: T = ß = 0.0083, 0.0167, 0.0333, 0.0667, 0.1 K/c. При меньшей скорости нагрева максимумы пиков достигаются при меньшей температуре (но за большее время).
./(7) /?"'•! 0
V14
/ AVeak #2
Peak
Рис. 8. Модельные спектры, влияние в Fig. 8. Model spectrum, effect of в
в [16, 18]. «Раздвоение» (бифуркация) спектра появляется при определенном соотношении энергий активации диффузии и десорбции, когда ни один из процессов не является строго лимитирующим и проявляется взаимообусловленность их динамики. В «синтезированной» модели для усредненной концентрации использование параметра а & [1,2} позволяет учитывать интегрально вклад (доли) диффузионных и поверхностных процессов.
Представление в линейных координатах на рис. 8 показывает несопоставимость масштабов пиков, обозначенных номерами 2 и 3 на Fig. 1,6 [8], с основным ТДС-всплеском. Выраженными эти «пики» становятся лишь при использовании логарифмической шкалы.
Более того, причина модельного пика 2-3 известна. После десорбции с поверхности и из подповерхностного слоя поток локально падает (пик 1). Затем поток начинает снова расти. Причина: дальнейший нагрев и образовавшийся большой градиент концентрации у поверхности. Падение градиента замедляет скорость диффузионной подкачки, и образуется естественное плечо 2-3, которое не является следствием каких-либо реакций 1-2 порядков в объеме (ловушек с «энергиями связи» нет). Это в модели, когда ответ на вопрос о причинах изгибов спектра известен после численного решения краевой задачи. Приписывание каждому экспериментальному локальному пику энергий связи в объеме требует дополнительного физического обоснования.
Модель динамического взаимодействия «поверхность-объем» позволяет воспроизвести экспериментальный спектр с пиками 1-4 (Fig. 5 [8]). Один из способов добавления пятого широкого высокотемпературного пика в модельном спектре — учет дефектов:
dtc = Ddze -
[1 - Zv] c -
a+ zv
Рис. 9. Спектры в логарифмической шкале Fig. 9. Spectrum in logarithmic scale
Алгоритм решения краевых задач термодесорбции на основе разностных схем (включая учет обратимого захвата диффузанта различного рода ловушками) подробно описан в [15]. Аппроксимации в классе ODE представлены
dtzv = a-(T)[1 - Zv] c(t,x) - a+(T)z„ (t,x).
Здесь zv(t,x) — концентрации атомов водорода, захваченного дефектами различных типов; aJ — коэффициенты поглощения и выделения H ловушками; Zv = zv(t,x)/zmax — степень насыщения (zmax = max zv). Для практических целей захват учтен в простейшей «интегральной» форме, уточнение геометрии дефектов и их распределения существенно усложнило бы модель. Если дефект, например, не микрополости, а включения гид-ридной фазы, то на этапе дегазации соответствующий коэффициент a- (T) = 0, а значение
a+(T) > 0 лишь после достижения критической температуры: T(t) ^ Tu-it.
53
о
600
800
1 000
1 200
1 400
1 600
m
a
Например, для модельного воспроизведения высокотемпературного пика 5 (Fig. 5 [8]) в краевую задачу достаточно добавить одну ловушку типа «разложения гидридной фазы» (в процессе дегазации а- =0), динамика концентрации в которой определяется соотношением dtz = -a+(T)z. Параметры: z(0) = 1018 см-3;
= 50 c ; Ea = 150; дефект начинает выделять водород при достижении критической температуры Tcrit = 1200 K.
Воспримем теперь ТДС-спектр, генерируемый динамической моделью, как экспериментальный и применим к пикам ##1-3 на рис. 9 (по аналогии с пиками 1, 4, 5 в нумерации [8]) методику разложения на «составляющие» элементарные реакции и методику Киссинджера. Соответствующие этим локальным пикам температуры: T\-3 & 650,1300,1600 K. Плечо 2-3, заметное лишь в логарифмической шкале, относим к естественному замедлению диффузии-десорбции и не воспринимаем как проявление реакций в объеме. Обработка «экспериментальных данных» в априорном предположении наложения реакций 1-го порядка (для определенности) дает следующие результаты (см. рис. 10). Энергия активации, соответствующая Tm(Peak #1), имеет значение 264.5 кДж/моль. Это соответствует Eb = 250, так что первый локальный пик спектра можно считать именно десорбционным. На его долю приходится малое количество выделившегося водорода (менее 0.2%). Второй выраженный пик Tm(Peak #2) дает «диффузионное» значение 125.3 & 125 = Ed. Высокотемпературному пику соответствует (естественно в силу модели) 150 = Ea. При этом в координатах Киссинджера {1/Tm, \и[в/Тт]} имеем практически прямые во всех трех случаях.
ln
ßT„;]
[E"] = kJ mol-1 Peak #Ц
Peak #2 \
Peak #3 = 125241 \
\ \ 1 E'= 264.491
£3*=149.926 -1 T m
Рис. 10. Оценка энергий активации Fig. 10. Estimation of activation energies
Таким образом, в приведенном модельном примере методика разложения на элементарные реакции (в терминах усредненной объемной концентрации) дает хорошие приближения, но не для «энергий связи» захваченного в объеме водорода, а для энергий активации физико-химических процессов десорбции, диффузии, разложения «гидридной фазы» (или оксида на поверхности при соответствующем увеличении эффективного коэффициента b) в их динамическом взаимодействии.
Финансирование исследований осуществлялось из средств федерального бюджета на выполнение государственного задания КарНЦ РАН (Институт прикладных математических исследований КарНЦ РАН).
Литература
1. Взаимодействие водорода с металлами / Ред. А. П. Захаров. М.: Наука, 1987. 296 с.
2. Водород в металлах / Ред. Г. Алефельд, И. Фелькль. М.: Мир, 1981. Т. 1, 2. 506, 430 с.
3. Изотопы водорода. Фундаментальные и прикладные исследования / Ред. А. А. Юхим-чук. Саров: РФЯЦ-ВНИИЭФ, 2009. 697 с.
4. Нечаев Ю. С. О природе, кинетике и предельных значениях сорбции водорода углеродными наноструктурами // Успехи физических наук. 2006. Т. 176(6). С. 581-610.
5. Писарев А. А., Цветков И. В., Марен-ков Е. Д., Ярко С. С. Проницаемость водорода через металлы. М.: МИФИ, 2008. 144 с.
6. Шьюмон П. Диффузия в твердых телах. М.: Металлургия, 1966.
7. Atsumi H., Takemura Y., Miyabe T, Konishi T, Tanabe T, Shikama T. Desorption of hydrogen trapped in carbon and graphite // J. Nucl. Mater. 2013. Vol. 442(1-3). P. S746-S750. doi: 10.1016/j.jnucmat.2013.03.041
8. Atsumi Н., Kondo Y. Retention and release of hydrogen isotopes in carbon materials priorly charged in gas phase // Fusion Eng. Des. 2018. Vol. 131. P. 49-53. doi: 10.1016/j.fusengdes.2018.04.039
9. 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. Hydrog. Energy. 2015. Vol. 40. P. 2871-2881. doi: 10.1016/j.ijhydene.2014.12.069
10. Nechaev Yu. S., Alexandrova N. M., Shurygina N. A., Cheretaeva A. O., Pisarev A. A. On the kinetic analysis of the hydrogen thermal desorption spectra for graphite and advanced carbon nanomaterials // Fullerenes, Nanotubes
54
-16
-17
-18
-19
and Carbon Nanostructures. 2019. Vol. 28, iss. 2. P. 147-149. doi: 10.1080/1536383X.2019.1680982
11. Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface // Int. J. Hydrog. Energy. 2011. Vol. 36(1). P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
12. Vyazovkin S., Burnham A. K., Criado J. M., Perez-Maqueda L. A. ICTAC Kinetics Committee recommendations for performing kinetic computations on thermal analysis data // Thermochim. Acta. 2011. Vol. 520. P. 1-19. doi: 10.1016/j.tca.2011.03.034
13. Vyazovkin S., Chrissafis K., Di Lorenzo M. L., Koga N. ICTAC Kinetics Committee recommendations for collecting experimental thermal analysis data for kinetic computations // Thermochim. Acta. 2014. Vol. 590. P. 1-23. doi: 10.1016/j.tca.2014.05.036
14. Wei F.-G., Enomoto M, Tsuzaki K. Applicability of the Kissinger's formula and comparison with the McNabb-Foster model in simulation of thermal desorption spectrum // Comp. Mater. Sci. 2012. Vol. 51. P. 322-330. doi: 10.1016/j.commatsci.2011.07.009
15. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption // Adv. Mater. Sci. Appl. 2014. Vol. 3, iss. 3. P. 120-129. doi: 10.5963/AMSA0303003
16. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermal desorption by ODE-approximation // Int. J. Hydrog. Energy. 2017. Vol. 42(1). P. 405-415. doi: 10.1016/j.ijhydene.2016.10.104
17. 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. Modelling. 2009. Vol. 33. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018
18. Zaika Yu. V., Sidorov N. I., Fomkina O. V. Identification of hydrogen permeability parameters of membrane materials in an aggregated experiment // Int. J. Hydrog. Energy. 2020. Vol. 45. P. 7433-7443. doi: 10.1016/j.ijhydene.2019.04.098
Поступила в редакцию 18.05.2020
References
1. Vzaimodeistvie vodoroda s metallami [Interaction of hydrogen with metals]. Moscow: Nauka, 1987. 296 p.
2. Vodorod v metallah [Hydrogen in metals]. Moscow: Mir, 1981. Vol. 1, 2. 506, 430 p.
3. Izotopy vodoroda. Fundamental'nye i prikladnye issledovaniya [Hydrogen isotopes. Fundamental and applied studies]. Sarov: RFYaTs-VNIIEF, 2009. 697 p.
4. Nechaev Yu. S. O prirode, kinetike i predel'nykh znacheniyakh sorbtsii vodoroda uglerodnymi nanostrukturami [On the nature, kinetics and limit values of hydrogen sorption by carbon nanostructures]. Uspekhi fizicheskikh nauk [Advances in Physical Sci.]. 2006. Vol. 176(6). P. 581-610.
5. Pisarev А. А., Tsvetkov I. V., Мarenkov Е. D., Yarko S. S. Pronitsaemost' vodoroda cherez metally [Hydrogen permeability through metalls]. Moscow: MIFI, 2008. 144 p.
6. Sh'yumon P. Diffuziya v tvyordykh telakh [Diffusion in solids]. Moscow: Metallurgiya, 1966.
7. Atsumi H., Takemura Y., Miyabe T., Konishi T., Tanabe T., Shikama T. Desorption of hydrogen trapped in carbon and graphite. J. Nucl. Mater. 2013. Vol. 442(1-3). P. S746-S750. doi: 10.1016/j.jnucmat.2013.03.041
8. Atsumi Н., Kondo Y. Retention and release of hydrogen isotopes in carbon materials priorly charged in gas phase.
Fusion Eng. Des. 2018. Vol. 131. P. 49-53. doi: 10.1016/j.fusengdes.2018.04.039
9. 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. Hydrog. Energy. 2015. Vol. 40. P. 2871-2881. doi: 10.1016/j.ijhydene.2014.12.069
10. Nechaev Yu. S., Alexandrova N. M., Shurygina N. A., Cheretaeva A. O., Pisarev A. A. On the kinetic analysis of the hydrogen thermal desorption spectra for graphite and advanced carbon nanomaterials. Fullerenes, Nanotubes and Carbon Nanostructures. 2019. Vol. 28, iss. 2. P. 147-149. doi: 10.1080/1536383X.2019.1680982
11. Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface. Int. J. Hydrog. Energy. 2011. Vol. 36(1). P. 1239-1247. doi: 10.1016/j.ijhydene.2010.06.121
12. Vyazovkin S., Burnham A. K., Criado J. M., Perez-Maqueda L. A. ICTAC Kinetics Committee recommendations for performing kinetic computations on thermal analysis data. Thermochim. Acta. 2011. Vol. 520. P. 1-19. doi: 10.1016/j.tca.2011.03.034
13. Vyazovkin S., Chrissafis K., Di Lorenzo M. L., Koga N. ICTAC Kinetics Committee recommendations for collecting experimental thermal analysis data for kinetic computations. Thermochim. Acta. 2014. Vol. 590. P. 1-23. doi: 10.1016/j.tca.2014.05.036
©
14. Wei F.-G., Enomoto M, Tsuzaki K. Applicability of the Kissinger's formula and comparison with the McNabb-Foster model in simulation of thermal desorption spectrum. Comp. Mater. Sci. 2012. Vol. 51. P. 322-330. doi: 10.1016/j.commatsci.2011.07.009
15. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermodesorption. Adv. Mater. Sci. Appl. 2014. Vol. 3, iss. 3. P. 120-129. doi: 10.5963/AMSA0303003
16. Zaika Yu. V., Kostikova E. K. Computer simulation of hydrogen thermal desorption by ODE-approximation. Int. J. Hydrog. Energy. 2017. Vol. 42(1). P. 405-415. doi: 10.1016/j.ijhydene.2016.10.104
СВЕДЕНИЯ ОБ АВТОРАХ:
Заика Юрий Васильевич
руководитель лаб. моделирования природно-технических систем, д. ф.-м. н.
Институт прикладных математических исследований
КарНЦ РАН, Федеральный исследовательский центр
«Карельский научный центр РАН»
ул. Пушкинская 11, Петрозаводск,
Республика Карелия, Россия, 185910
эл. почта: zaika@krc.karelia.ru
тел.: (8142) 780059
Костикова Екатерина Константиновна
научный сотрудник, к. ф.-м. н.
Институт прикладных математических исследований
КарНЦ РАН, Федеральный исследовательский центр
«Карельский научный центр РАН»
ул. Пушкинская, 11, Петрозаводск,
Республика Карелия, Россия, 185910
эл. почта: kostikova@krc.karelia.ru
тел.: (8142) 766312
17. 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. Modelling. 2009. Vol. 33. P. 3776-3791. doi: 10.1016/j.apm.2008.12.018
18. Zaika Yu. V., Sidorov N. I., Fomkina O. V. Identification of hydrogen permeability parameters of membrane materials in an aggregated experiment. Int. J. Hydrog. Energy. 2020. Vol. 45. P. 7433-7443. doi: 10.1016/j.ijhydene.2019.04.098
Received May 18, 2020
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
Kostikova, Ekaterina
Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Sciences 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia
e-mail: kostikova@krc.karelia.ru tel.: (8142) 766312
