Kinetic model of the temperature-programmed desorption of ammonia to study the acidity of heterogeneous catalysts

Full Text

ВВЕДЕНИЕ

Изучение катализаторов в процессе температурно-программированной десорбции (ТПД) аммиака является стандартным методом физико-химической идентификации кислотных свойств поверхности [1–3]. Однако сам процесс чаще всего сводится к определению интегральной кислотности и разделению кислотных центров по силе на слабые, средние и сильные [4–6]. Такое описание процесса является эмпирическим и малоинформативным, поскольку не объясняет структуру центров адсорбции и десорбции аммиака и не имеет предсказательного характера [7–9]. К тому же, сравнение кислотных центров катализаторов при таком подходе возможно только для образцов, исследованных на одном и том же приборе в абсолютно одинаковых условиях. Если поменять какой-либо из параметров процесса, например, скорость нагрева, картины ТПД аммиака в значительной степени могут измениться.

В области гетерогенного катализа мониторинг специфических свойств катализатора, таких как количество кислотных, основных или металлических центров, может иметь первостепенное значение для разработки составов и способов приготовления катализаторов, что может значительно увеличить выходы и селективности целевых продуктов. Кроме того, строгое количественное описание свойств имеет основополагающее значение для развития математического аппарата материаловедения. Получение стандартизованных данных с помощью метода ТПД аммиака позволит сравнивать катализаторы, полученные различными группами исследователей, что необходимо в том числе для контроля и оптимизации промышленных катализаторов. Поэтому разработка кинетических моделей процесса десорбции аммиака, является актуальной практической задачей [10].

На сегодняшний день единственным строгим приложением данного метода является определение концентрации кислотных центров [11, 12]. Наличие двух и более перекрывающихся пиков при десорбции аммиака не позволяет корректно разделить их, поскольку не существует единой физической модели для обработки интегральных кривых: наиболее часто применяют деконволюцию кривой ТПД на Гауссовы пики [13]. Однако такой эмпирический метод обработки экспериментальных результатов в большинстве случаев лишен физического смысла [14]. Разработка аналитического метода обработки кривых ТПД аммиака позволит получать информацию не только о количестве кислотных центров, но и о их силе и типе, подвижности адсорбированных молекул аммиака, что будет являться следующим шагом эволюции теории ТПД [15, 16].

Стоит отметить, что уже опубликован ряд моделей математического описания результатов температурно-программируемой десорбции аммиака, как, например, в [17]. Предложенный авторами подход основывается на псевдоравновесных процессах реадсорбции аммиака и имеет описательный характер только для моноцентровой адсорбции, когда можно выделить максимальное значение пика десорбции, по которому принято оценивать изменение энтальпии адсорбции. Использование такого реадсорбционного приближения не оправдано в некоторых случаях некорректно.

Авторы [18] проводили аналитическую деконволюцию, основываясь на кинетике протекания десорбции, однако для исследованного оксида алюминия в низкотемпературной области также была рассмотрена реадсорбция аммиака.

Авторы [18–20] предлагают использовать для описания кинетики десорбции теорию Поляни–Вигнера:

$-\frac{d{\theta }_{{\text{NH}}_3}}{dT}=\frac{\nu }{\beta }{\theta }_{{\text{NH}}_3}^n\exp \left(-\frac{E_D}{RT}\right),$ (1)

где q – степень заполнения кислотных центров поверхности, n – наблюдаемый порядок реакции, E_D – энергия активации процесса десорбции, n – свободная частота колебания молекулы аммиака на поверхности катализатора, b = dT/dt – скорость нагрева, которая в методике проведения ТПД чаще всего является постоянной величиной. Данная модель может подойти для исследования кислотности целого ряда катализаторов, но описательный характер модели имеет ограничения, связанные с выбором наблюдаемых порядков, которые классически варьируются в диапазоне n = 0, 1, 2 [21].

Теория ТПД

Simoní Da Ros и др. в своей работе описывали десорбцию аммиака на оксиде алюминия, применяя кинетическую модель обратимой десорбции на основе простой модели Ленгмюра [18]:

$\frac{d{\theta }_{{\text{NH}}_3,i}}{dt}=k_{A,i}{P}_{{\text{NH}}_3}\left(1-{\theta }_{{\text{NH}}_3,i}\right)-k_{D,i}{\theta }_{{\text{NH}}_3,i},$ (2)

$\frac{d{C}_{{\text{NH}}_3}}{dt}=\frac{-C_{{\text{NH}}_3}}{\tau }+\frac{1-\epsilon }{\epsilon }{\rho }_p\left(-\sum _{i=1}^{NS}{N}_i\frac{d{\theta }_{{\text{NH}}_3,i}}{dt}\right),$ (3)

где k_A,i – константа скорости адсорбции, k_D,i – константа скорости десорбции аммиака с i-го типа центров, q_NH3,i и N_i – степень заполнения поверхности i-го типа кислотных центров и их интегральная концентрация, r – насыпная плотность образца, P_NH3 – текущее давление аммиака, ε – пористость, $\frac{C_{{\text{NH}}_3}}{\tau }$ – скорость массопереноса аммиака из реактора, которым в модели пренебрегают из-за малой загрузки образца (~200–700 мг).

При этом стоит отметить, что для первого типа кислотных центров (низкотемпературные пики) на оксиде алюминия авторами [18] рассматривалась кинетика с учетом процесса реадсорбции. В данной работе эксперимент проводили с постоянной подачей азота 19 мл/мин, в токе которого на выходе из реактора проводили регистрацию аммиака. После обработки данных для низкотемпературного пика были представлены значения энергий активации для адсорбции 73.05 кДж/моль и десорбции 113.1 кДж/моль, предэкспоненциальных множителей адсорбции 0.0488 атм^–1 мин^–1 и десорбции 88.23 мин^–1 с общей концентрацией адсорбционных центров 0.696 ммоль/г (при насыпной плотности образца оксида алюминия 0.96 г/см³). По этим данным зависимости концентрации выделяющегося аммиака от времени, температурных зависимостей констант скоростей адсорбции и десорбции, интегральной кислотности и насыпной плотности катализатора, можно оценить значения констант равновесия, степени заполнения поверхности и, в результате, зависимость D_гG(T) для процесса десорбции с этих центров в диапазоне температур 150–600°C.

Результат этой оценки представлен на рис. 1 (подробнее см. дополнительных материалах к статье). Известно, что, если изменение потенциала Гиббса процесса меньше нуля более чем на 3RT [22], то реакция протекает практически необратимо, т. е. в данном случае предположение о равновесности десорбции и, соответственно, реадсорбции аммиака является неверным. Стоит отметить, что полученные значения стандартных энтальпии и энтропии (40.03 кДж/моль и 62.36 Дж/(моль К) соответственно) для процесса десорбции аммиака в данной работе, согласуются с общепринятыми значениями, которые были получены в других работах [20, 23].

 

Рис. 1. Оценка изменения потенциала Гиббса для процесса обратимой десорбции аммиака с низкотемпературных центров оксида алюминия в сравнении с порогом обратимости реакции [18].

 

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

По этим причинам мы предлагаем кинетическую модель необратимой десорбции аммиака с использованием подхода, представленного в [19]. В соответствие с моделью Поляни–Вигнера по уравнению (1) для каждого i-го типа центров можно записать кинетическое уравнение десорбции аммиака как скорость уменьшения степени заполнения поверхности аммиаком от температуры при постоянной скорости нагрева b:

$-\frac{d{\theta }_{{\text{NH}}_3,i}}{dT}=\frac{{\nu }_i}{\beta }{\theta }_{{\text{NH}}_3,i}^{n_i}\exp \left(-\frac{E_{D,i}}{RT}\right),$ (4)

при этом данное кинетическое уравнение формально описывает химическую реакцию на поверхности катализатора через активированное состояние с участием n_i адсорбированных молекул аммиака в следующем виде:

$\begin{array}{c}{n}_i{\theta }_{N{H}_3,i}\stackrel{K^{\ne }}{\rightleftarrows }{\left[N{H}_3...{\theta }_i...(n_i-1) {\theta }_{N{H}_{3,}i}\right]}^{\ne }\stackrel{{\nu }_i^{\ne }}{\to }\\ (n_i-1){\theta }_{N{H}_3,i}+{\theta }_i+N{H}_3\uparrow ,\end{array}$ (5)

Подобное описание процесса встречается в ряде работ, где между реагентами и активированным состоянием принимается равновесное распределение, сильно смещенное в сторону реагентов [24]. При этом в моделировании процесса десорбции иногда принимают, что адсорбированная молекула в активированном состоянии обладает свойствами, аналогичными газообразному состоянию. Подобное допущение действительно может быть использовано при моделировании физической адсорбции по Ленгмюру, где энергия активации в точности совпадает с энтальпией десорбции. В этом случае между активированным состоянием и реагентами рассматривается равновесие как между газом и адсорбированной молекулой с некоторыми дополнительными условиями: адсорбированный газ является решеточным газом и не может обладать всеми поступательными степенями свободы (сохраняется от 0 до 2 степеней в зависимости от энергии связи молекулы с поверхностью) [25], а вращательные степени свободы в активированном состоянии полностью сохраняются как для газа и представляют собой своеобразный “волчок” [26], при этом в полном объеме учитываются колебательные степени свободы. Чаще всего в таком виде теорию активированного состояния применяют для оценки предэкспоненциального множителя, чтобы упростить расчеты и поиск наблюдаемого порядка реакции (энергию активации определяют традиционным для термического анализа методом Фридмана) [27].

В данной работе мы не используем теорию активированного комплекса, а вместо этого, опираясь на уравнение Поляни–Вигнера, определяем кинетические параметры методами математической обработки, в первую очередь, исходя из следующих предположений. Степень заполнения поверхности для i-х центров напрямую связана с их общей кислотностью, которая при адсорбции и десорбции является величиной неизменной, поэтому справедливы уравнения:

$\left\{\begin{array}{c}{h}_i=-\frac{d{H}_i}{dT}=\frac{k_0,i}{\beta }\exp \left(-\frac{E_{D,i}}{RT}\right)H_i^n,\\ h_i=-H_{\text{max},i}\frac{d{\theta }_{{\text{NH}}_3,i}}{dT},\\ H_i=H_{\text{max},i}\cdot {\theta }_{{\text{NH}}_3,i},\\ {\theta }_{{\text{NH}}_3,i}+{\theta }_i=\mathrm{1,}\\ h_{кст}=\underset{i=1}{\overset{n}{}}{h}_i,\\ H_{кст}=\underset{i=1}{\overset{n}{}}{H}_{\text{max}.i},\end{array}\right.$ (6)

где h_i – дифференциальная кислотность i-го центра, равная скорости выделения аммиака с поверхности катализатора, H_i – интегральная кислотность, равная остаточной концентрации аммиака на поверхности катализатора, H_max,i – общая кислотность i-х центров адсорбции аммиака, равная их начальной концентрации перед десорбцией, h_кат – скорость выделения аммиака со всех центров в данный момент времени и при данной температуре – именно эта величина чаще всего и представляется в виде результата ТПД NH₃, H_кат – общее значение кислотности, представляющее собой концентрацию всех центров адсорбции на исследуемом катализаторе, k_0,i – Аррениусовский предэкспоненциальный множитель, размерность которого зависит от наблюдаемого порядка реакции. Кинетическая модель, состоящая из уравнений (4)–(6) и предложенная в работе [19], описывает процесс десорбции, однако следует учитывать, что реакции десорбции должны быть простыми и иметь активационный барьер, то есть, в отличие от модели в работе [19], в которой считалось, что порядки n = 0, 1, 2. В предлагаемой нами модели порядки должны быть только целочисленными и в следующем диапазоне значений n = 1, 2, 3, что соотносится с критериями, представленными в работах С. Вязовкина [28]. При этом роль соседних адсорбированных молекул аммиака в такой модели протекания десорбции можно связать с восполнением недостатка энергии, затрачиваемой на преодоление активационного барьера процесса.

Именно из этих положений мы и будем исходить при описании процесса десорбции аммиака с поверхности катализатора. Данная работа представляет алгоритм подбора параметров уравнений кинетических кривых в соответствии с (4)-(6) и демонстрирует возможности предлагаемого подхода для получения количественных характеристик кислотности поверхности ряда материалов.

ЭКСПЕРИМЕНТАЛЬНАЯ ЧАСТЬ

Исследуемые системы

Для исследования были выбраны материалы, широко используемые в качестве носителей и катализаторов: 1) цеолиты со структурным типом MFI, для которого кривая десорбции аммиака имеет один пик, с различным химическим составом; 2) ряд образцов на основе γ-Al₂O₃ с более сложной кривой десорбции полученных из различных предшественников. К первому типу были отнесены фазовочистый силикалит^–1 (Sil^–1), ранее приготовленный и представленный в работе [29], а также титансиликалит^–1 (TS^–1) ZD-07031, производства “Zeolyst International” (Германия). Для изучения кинетики десорбции аммиака в качестве алюмооксидных систем были выбраны три образца γ-Al₂O₃: 1) промышленный гранулированный мезопористый γ-оксид алюминия A1 [30] производства АО “Ангарский завод катализаторов и органического синтеза” (Россия); 2) два образца гранулированного макро-мезопористого γ-Al₂O₃, приготовленных из псевдобемита производства ЗАО “Алтайлюминофор” (Россия) по методике экструзии пластичной массы, содержащей 30 мас. % полимерного темплата для формирования макропор, с последующим прокаливанием при температурах 700 и 800°C (образцы Yar-700 и Yar-800 соответственно) по аналогии с методиками, приведенными в работах [31, 32]. Все образцы были охарактеризованы с помощью методов рентгеновской дифракции и имеют подтвержденную структуру (рис. 2). Свойства пористой структуры были рассчитаны из изотерм низкотемпературной адсорбции-десорбции азота образцами, полученных на приборе Autosorb 6B-KR (Quantachrome, США). Предварительно образец дегазировали при 150°C в течение 3 ч. Расчет удельной поверхности проводили по уравнению Брунауэра–Эммета–Теллера (БЭТ), расчет объема пор – по величине адсорбции при максимальном достигнутом давлении. Элементный состав образцов определяли методом атомно-эмиссионный спектрометрии с возбуждением в индуктивно-связанной плазме (ИСП-АЭС) на спектрометре iCAP 6500 Duo (ThermoElemental, США). Физико-химические характеристики исследованных материалов представлены в табл. 1.

 

Таблица 1. Физико-химические свойства исследованных образцов

Образец	Sуд, м2/г	Vпор, мл/г	Элементный состав
ZD-03071	508	0.82	1.9 мас. % Ti 75 мас. % Si
Sil–1	414	0.89	47 мас. % Si
A1	240	0.59	51–53 мас. % Alа
Yar-700	186	0.60
Yar-800	157	0.63

Примечание. ^аДля метастабильных оксидов алюминия потери при прокаливании выше 700^°С меньше относительной погрешности определения элементного состава методом ИСП-АЭС; S_уд – удельная поверхность, V_пор– объем пор.

 

Рис. 2. Дифрактограммы исследуемых образов: силикалита^–1 (Sil^–1) (а) и оксида алюминия (Yar-800) (б). Данные сравнения (синий) взяты из COD для структурного типа MFI цеолита и ICDD PDF-2 для γ-Al₂O₃.

 

Перед измерением ТПД гранулированные материалы были измельчены до фракции 0.2–0.8 мм для снижения влияния диффузии.

Методика проведения эксперимента

Температурно-программируемую десорбцию аммиака проводили в кварцевом реакторе с внутренним диаметром ∅ 8 мм. Для качественного анализа использовали навеску анализируемого образца массой 200–700 мг для получения достоверных данных по десорбции, при котором вклад эффектов диффузии и теплопереноса незначителен. Образец перед измерением измельчали и просеивали до фракции 0.2–0.8 мм, затем засыпали в реактор. Регистрацию температуры проводили термопарой, помещенной в слой образца. На протяжении всего эксперимента через реактор пропускали гелий (производство “Оренбургский гелиевый завод”, филиал ООО “Газпром переработка”) с постоянной скоростью ~60 мл/мин. Эксперимент проводился в несколько стадий.

1. Тренировка образца в токе гелия. Такая подготовка проводилась для удаления с поверхности, в первую очередь, адсорбированной воды. На этой стадии реактор нагревали до температуры 600°C со скоростью 600 К/ч, затем выдерживали еще 30 минут при постоянной температуре. Затем реактор охлаждали до температуры 100°C.
2. Насыщение образца аммиаком (ООО “Чистые Газы Плюс”, Новосибирск). Для этого предварительно осушенный над реактивированном в NaOH (производство ПАО “Химпром”, Новочебоксарск) угле АГМ^–1 (производство АО “Сорбент”, Пермь) аммиак подавали со скоростью 30 мл/мин в поток гелия 60 мл/мин в течении 30 мин до полного насыщения анализируемого образца аммиаком, что проверяли с помощью масс-спектрометра RGA100 (производство “Stanford Research System”) по выходной кривой. Затем подачу аммиака прекращали и продували реактор чистым гелием при фиксированной температуре 100°C еще в течении 30 минут для удаления остатков газообразного и физически адсорбированного аммиака с поверхности анализируемого катализатора.
3. Десорбция аммиака. Этот этап проводили с приблизительно постоянной скоростью нагрева ~600 K/ч в потоке гелия ~60 мл/мин. На выходе из реактора часть потока отбирали для количественного анализа выделенного аммиака с помощью масс-спектрометра.

Регистрацию сигнала аммиака проводили одновременно с регистрацией сигнала гелия при следующих настройках ионизатора: 24 В, 10Вт. Давление в измерительной камере во время эксперимента устанавливали 10–⁵ мбар, регулируя затвор входного натекателя. Для анализа ТПД NH₃ использовали 2 регистрационных сигнала (4 а. е. м. для гелия и 17 а. е. м. для аммиака), проводя съемку в диапазоне линейности сигнала (интенсивность сигнала аммиака не превышает 10% порог от интенсивности сигнала гелия). После съемки полученный сигнал масс-спектрометра вначале синхронизовали по временной шкале с температурным, учитывая мертвый объем и задержку начала нагрева после запуска. Затем, используя заранее калиброванные значения фактора чувствительности детектора к аммиаку в гелиево-аммиачной смеси f_NH3 = 0.2055, сигналы пересчитывали в регистрируемую скорость выделения аммиака, считая газы идеальными, по формуле:

$\begin{array}{c}\overline{n_{{\text{NH}}_3}}\left[моль/\left(сг\right)\right]=\\ =\frac{f_{{\text{NH}}_3}\cdot F_{{\text{NH}}_3}}{F_{\text{He}}}{U}_{\text{He}}\cdot \frac{p_0}{RT\cdot m_{образца}},\end{array}$ (7)

где p₀ – атмосферное давление (Па), T – температура (К), R – универсальная газовая постоянная (Дж/(моль·К)), U_He – объемная скорость потока гелия (м3^/с), F_NH3, F_He – зарегистрированные сигналы аммиака и гелия в текущий момент времени соответственно, f_NH3 – фактор чувствительности детектора к аммиаку в гелии, m_{образца} – масса образца (г).

Используя полученную скорость выделения аммиака, получали удельную кислотность поверхности образца по уравнению:

$h_{кст}\left[ммоль\text{/}\right(Кмл\left)\right]=\overline{n_{{\text{NH}}_3}}\cdot \frac{\rho }{\beta },$ (8)

где r – насыпная плотность образца (г/мл), b – скорость нагрева (К/с). Полученную зависимость значения удельной кислотности h_кат(T) от температуры использовали для математического моделирования кинетики десорбции NH₃.

Для всех анализируемых образцов были получены кривые ТПД аммиака, которые представлены на рис. 3. Условия регистрации, проведения процесса десорбции и физические характеристики для каждого образца указаны в табл. 2.

 

Рис. 3. Кривые ТПД аммиака исследуемых образцов.

 

Таблица 2. Условия проведения ТПД аммиака для исследуемых образцов

Образец	UHe, см3/мин	b, К/ч	mобразца, мг	r, г/см3
ZD-03071	59.4	638	803.45	0.98
Sil–1	60.7	598	360.00	1.00
A1	59.1	753	544.26	0.96
Yar-700	58.3	639	530.86	0.94
Yar-800	61.0	646	641.17	0.96

Обозначения: U_He – объемная скорость подачи гелия, β – скорость нагрева, m– масса образца, ρ – насыпная плотность.

 

Подбор параметров

Как видно из системы (6) и уравнения (4), имеется взаимосвязь между кинетическим предэкспоненциальным множителем k_0,i (М^(1–n)с^–1) и приведенной частотой колебания вдоль координаты реакции десорбции n_i (1/K):

$k_{\mathrm{0,}i}=\beta {\nu }_i{H}_{\max ,i}^{\left(1-n\right)}.$ (9)

Используя (9), можно разрешить поставленную задачу подбора параметров уравнения (4) n_i, n, E_D,i методом нелинейной минимизации и найти q_NH3,i(T). Используя q_NH3,i(T) можно оценить и $\frac{d{\theta }_{{\text{NH}}_3,i}}{dT}\left(T\right)$, которая напрямую связана с зависимостью дифференциальной кислотности i-х центров адсорбции от температуры h_i(T) по первому уравнению системы:

$\left\{\begin{array}{c}{h}_i=-H_{\max ,i}\frac{d{\theta }_{{\text{NH}}_3,i}}{dT},\\ h_{\text{cat}}\left(T\right)={\sum }_i{h}_i.\end{array}\right.$ (10)

Сумма дифференциальных кислотностей по всем центрам даст дифференциальную кислотность h_cat для образца (второе уравнение системы (10)).

Данное моделирование проводили для конкретного типа кислотных центров, используя скрипт для Mathlab по алгоритму решения дифференциальных уравнений методом Рунге–Кутта [33, 34]:

[t, y] = ode45(@(t, y) (–v)*exp(–Ed/8.314/t)*(y).^n, [Tmin: dT: Tmax], y0);%Дифф.уравнение

dtetha= –diff(y);%Расчет производной функции для оценки скорости десорбции

m=(Tmax-Tmin)/dT;%Расчет числа выходных точек для каждого массива данных

for l=[1: m];%Цикл формирования выходных массивов

T(l)=t(l)-273.15;%Перевод температуры из Кельвинов в градусы Цельсия

Hi(l)=Hmax*y(l);%Расчет интегральной кислотности i-го центра в точке T(l)

hi(l)=Hmax*dtheta(l);%Расчет скорости десорбции в точке T(l)

end;%Логический конец скрипта

В качестве входных параметров для расчета использовали Tmin, dT, Tmax – начальная температура, температурный шаг, конечная температура в К, Hmax – максимальная концентрация i-го центра кислотности в моль/л, v – приведенная частота колебания вдоль координаты десорбции в 1/K, Ed – энергия активации десорбции в Дж/моль, n – наблюдаемый порядок десорбции аммиака, y₀ – начальное значение степени заполнения поверхности аммиаком, равное 1. В качестве выходных данных получали массивы концентрации оставшихся заполненными аммиаком центров Hi(l) (моль/л), дифференциальной кислотности h_i(l) (ммоль/(К·мл))-х кислотных центров и соответствующие этим точкам температуры в °C. Данный скрипт можно использовать в интеграции с большинством языков программирования и решать задачу описания кинетических кривых, вариативным методом подбирая для каждого типа центров параметры n_i, E_D,i, n_i, H_max,i, получая на выходе набор массивов h_i(T), сумма которых для одинаковых температур должна соответствовать экспериментальному значению кислотности.

Анализ влияния кинетических параметров на модельные кривые

Для обоснования подхода к подбору параметров кинетической системы представим анализ их влияния на вид модельной кривой.

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

 

Рис. 4. Влияние различных параметров кинетического уравнения на общий вид кинетической кривой: наблюдаемого порядка (a); предэкспоненциального множителя (приведенной частоты вдоль координаты десорбции) (б); энергии активации десорбции (в); интегральной кислотности (г).

 

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

Стоит заранее оговориться: поскольку мы приняли, что процесс десорбции – простая реакция, то порядки у нас могут быть только целочисленными, а значит будем подбирать их, исходя из общего вида кривых и соотношения между экспериментальными ветвями пиков. Такой подход должен снизить время для подбора параметров в 3^m раз, где m – количество анализируемых пиков, поскольку вариацию порядков реакции проводить мы не будем. Для других параметров ситуация обстоит проще, поскольку изменение в их δ-окрестности не меняет глобально вид кривой, а лишь предсказуемо влияет на положение или высоту пика.

Изменение частоты колебания вдоль координаты реакции приводит к смещению пика в более низкотемпературную область с увеличением высоты пика и уменьшением его ширины (рис. 4б). Для энергии активации имеем противоположный эффект: рост энергетического барьера смещает положение пика в высокотемпературную область с уменьшением высоты и уширением пика (рис. 4в). Увеличение концентрации адсорбированного аммиака влияет исключительно на площадь и высоту пика, не изменяя положения его максимума (рис. 4г).

База данных кинетических кривых ТПД

Перейдем к методу подбора кинетических параметров. После получения кривой десорбции аммиака необходимо вначале провести инициализацию параметров модели, для которой, в свою очередь, требуется база данных. Поэтому предварительно была создана база данных значений кинетических параметров и соответствующих выходных точек в координатах h_i – T с шагом 1°C и интегральной кислотностью H_max = 1 ммоль/мл. Каждый массив таких данных содержит 1000 точек для определенного набора параметров. На момент подготовки данной работы база данных содержала 260 кривых, где точки были подобраны для различных значений n = 1, 2, 3, n от 2 до 10¹⁵ 1/К (с шаговым множителем), и E_D от 20 до 170 кДж/моль (со случайным шагом). При этом большая часть кривых соответствовала выходу максимальной скорости десорбции при температурах 100–800°C, а для каждого значения одного параметра в таблице содержались более десятка наборов кинетических данных с варьируемыми другими параметрами в широком диапазоне, оставляя постоянным только H_max. База данных периодически дополняется новыми экспериментальными значениями параметров и соответствующих им точек кинетических кривых.

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

Алгоритм инициализации начальных параметров для экспериментальных кривых ТПД

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

1. Определяем максимальное значение h_max из графика.
2. В ε-окрестности от максимума пика, которая определяется набором точек пика от 0.5h_max левой ветви до 0.8h_max правой ветви, выбираются все экспериментальные точки, а затем сравниваются с ранее проанализированными кривыми десорбции аммиака с известными параметрами и стандартной кислотностью H_max = 1 ммоль/мл из базы данных. При этом не рассматриваются кривые, у которых положение максимального значения скорости десорбции h отличается от экспериментального более чем на 10°C. Остальные параметры проверяют на дополнительные условия соответствия экспериментальной кривой.
3. В диапазоне точек от 0.5h_max левой ветви до 0.8h_max правой ветви определяется диапазон температур, по которому считаются интегралы для табличной и экспериментальной кривых. По их соотношению определяется значение H_max для табличной кривой, которое используется в дальнейших вычислениях (рис. 5а).
4. Далее для каждой инициализированной модельной кривой в заданном температурном диапазоне проводится расчет коэффициента детерминации R² по следующему алгоритму:

$\left\{\begin{array}{l}\overline{h}=\frac{1}{n}\sum _j^n{h}_j,{\sigma }_h^2=\sum _j^n{(h_j-\overline{h})}^2\\ {\sigma }^2=\sum _j^n{(h_j-h(T_j\left)\right)}^2,R^2=1-\frac{{\sigma }^2}{{\sigma }_h^2},\end{array}\right.$ (11)

где h_j – экспериментально полученное значение кислотности при температуре T_j, h(T_j) – модельное значение, рассчитанное для T_j, h – среднее значение экспериментально измеренной величины для всего зарегистрированного значения скорости десорбции. Для экспериментального пика выбирается тот набор табличных параметров, который обеспечивает наибольший коэффициент детерминации, как показано на рис. 5б.

 

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

 

Если R² > 0.95 для всей кривой ТПД, то программа автоматически предлагает использовать эти параметры в дальнейшем подборе и на этом инициализация завершается. Так, например, на рис. 5, представлен вариант автоматизированной инициализации параметров для результата ТПД при анализе коммерческого образца ZD-07031 (титансиликалит^–1), для подобранных параметров которого R² = 0.991;

Если коэффициент детерминации оказался меньше 0.95, то для экспериментальных точек рассчитывается разность с выбранной моделью и стадии 1–6 повторяют для инициализации следующего набора параметров из базы данных и выбора наиболее подходящего.

Процедура инициализации заканчивается при выборе трех кинетических кривых ТПД или пока коэффициент детерминации не будет более 0.95. При этом у пользователя всегда остается возможность ручной корректировки инициализированных данных с выводом результата в виде кривой h(T) и коэффициента детерминации для экспериментальных и модельных точек. Например, на рис. 6 представлена инициализация параметров для образца оксида алюминия Yar-800. Производилась последовательная инициализация каждой модельной кривой по разности с экспериментально полученными точками.

 

Рис. 6. Инициализация кинетических параметров для кривой ТПД γ-оксида алюминия Yar-800, подбор параметров первой кривой из базы данных с проверкой сходимости модели (а) и добавлением дополнительной кинетической кривой из-за низкого значения R² после подстановки первой кривой (б).

 

Алгоритм оптимизации параметров моделирования для улучшения сходимости с экспериментальными результатами

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

$Error=\sqrt{(1-R^2)}$ (12)

Найденное значение использовалось для оценки индивидуальных отклонений параметров в цикле оптимизации по следующим формулам:

$\left\{\begin{array}{l}d{\nu }_i=Error\cdot rat{e}_{{\nu }_i}\cdot {\nu }_i,\\ d{H}_{\max ,i}=Error\cdot rat{e}_{H_{\max ,i}}\cdot H_{\max ,i},\\ d{E}_{D,i}=Error\cdot rat{e}_{E_{D,i}}\cdot \mathrm{5180,}\end{array}\right.$ (13)

где dn_i – отклонение от текущего значения n_i, dH_max,i – отклонение от текущего значения H_max,i, dE_D,i – отклонение от текущего значения энергии активации десорбции аммиака, rate_xi – индивидуальные скорости подбора параметров модели для x_i параметра, которые численно перед первой итерацией принимаются равными единице. Можно заметить, что все отклонения кроме отклонения энергии активации десорбции напрямую зависят от текущих значений самих параметров. Однако такой подход для энергии активации десорбции из-за ее значительного влияния на скорость десорбции аммиака (энергия входит в формуле в экспоненциальную часть) не позволяет задать отклонение линейно зависящим от текущего значения E_D,i. В качестве стандартного значения для расчета отклонения энергии было выбрано значение, соответствующее среднему значению RT, т. е. при температуре 350°C, равному 5180 Дж/моль.

Скорости изменения rate_xi не являются постоянными величинами по ряду причин: 1) вблизи оптимального значения параметра слишком большая скорость его изменения не позволит достичь оптимума; 2) при постоянной индивидуальной скорости подбора вдали от оптимального значения изменение параметра будет недостаточным, чтобы быстро достичь оптимума. Поэтому индивидуальные скорости rate_xi также изменялись в зависимости от того, как происходили изменения самого параметра x_i модели по следующему алгоритму:

$\begin{array}{c}rat{e}_{n+\mathrm{1,}{x}_i}=\\ =\left\{\begin{array}{l}0.9\cdot rat{e}_{n,x_i},\text{\hspace{0.33em}если\hspace{0.33em}}{x}_{n+\mathrm{1,}i}=x_{n,i}\\ rat{e}_{n,x_i},\text{\hspace{0.33em}если\hspace{0.33em}}{x}_{n+\mathrm{1,}i}=x_{n,i}\pm d{x}_{n,i}\text{\hspace{0.33em}первый раз}\\ 1.1\cdot rat{e}_{n,x_i},\text{\hspace{0.33em}если\hspace{0.33em}}{x}_{n+\mathrm{1,}i}=x_{n,i}\pm d{x}_{n,i}\\ \text{болеее раза,}\end{array}\right.\end{array}$ (14)

где индекс n – номер текущей итерации по подбору параметров. Таким образом, если параметр в моделировании активно изменяется, то скорость его изменения увеличивается, в противном случае наоборот, происходит уменьшение скорости подбора для более точной оптимизации параметра – увеличения коэффициента детерминации.

При этом в каждом цикле подбора параметров происходит перебор всех возможных вариантов. Например, для описания кинетической кривой с двумя уравнениями одновременно перебираются все параметры для двух кривых:

$\begin{array}{c}\{{\nu }_1-d{\nu }_1; {\nu }_1; {\nu }_1+d{\nu }_1\},\\ \{E_{D\mathrm{,1}}-d{E}_{D\mathrm{,1}};E_{D\mathrm{,1}};E_{D\mathrm{,1}}+d{E}_{D\mathrm{,1}}\},\\ \{H_{\text{max}\mathrm{,1}}-d{H}_{\text{max}\mathrm{,1}}; H_{\text{max}\mathrm{,1}};H_{\text{max}\mathrm{,1}}+H_{\text{max}\mathrm{,1}}\},\\ \{{\nu }_2-d{\nu }_2; {\nu }_2; {\nu }_2+d{\nu }_2\},\\ \{E_{D\mathrm{,2}}-d{E}_{D\mathrm{,2}};E_{D\mathrm{,2}};E_{D\mathrm{,2}}+d{E}_{D\mathrm{,2}}\},\\ \{H_{\text{max}\mathrm{,2}}-d{H}_{\text{max}\mathrm{,2}}; H_{\text{max}\mathrm{,2}};H_{\text{max}\mathrm{,2}}+d{H}_{\text{max}\mathrm{,2}}\}.\end{array}$

Экспериментальные точки сравниваются с каждым из полученных вариантов кинетических кривых и для каждого определяется R². Полученные значения коэффициента детерминации сравниваются между собой и определяется один самый подходящий набор параметров из 3⁶ = 729. Именно по этому набору и происходит пересчет индивидуальной скорости подбора параметров с повторением процедуры оптимизации на следующем шаге в соответствии с уравнением (14).

Подбор параметров заканчивают при одном из двух условий: 1) если получаемый для параметров коэффициент детерминации R² = 1, что означает отсутствие ошибки моделирования; 2) если все параметры перестают изменяться на протяжении более 5 итераций подряд с понижением их индивидуальных скоростей подбора, где отклонения и будут являться ошибками полученных параметров модели.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Заключительным этапом моделирования является расчет истинных кинетических предэкспонентов в соответствии с формулой (9) и вывод всех найденных параметров модели. Блок-схема описанного подбора параметров представлена на рис. 7.

 

Рис. 7. Блок-схема алгоритма подбора параметров кинетической модели ТПД NH₃.

 

Для подтверждения применимости модели необходимо продемонстрировать, что десорбция действительно может протекать по реакции третьего порядка. Более того, желательно показать, что процесс может быть моноцентровым, т. е. протекает при десорбции с центров, обладающих одним набором энергетических свойств (одинаковые предэкспоненты и энергии десорбции). Оказалось, что цеолитные системы обладают именно такими свойствами (рис. 8). В ходе подбора параметров для титансиликалита^–1 ZD-07031 и силикалита^–1 Sil^–1, обладающих одним структурным типом MFI, десорбция может быть описана одной кинетической кривой третьего порядка. Вопрос о моноцентровой модели остается дискуссионным и будет детальнее изучен в дальнейших работах с цеолитоподобными материалами.

 

Рис. 8. Кривые ТПД аммиака для образцов ZD-07031 (а) и Sil^–1 (б) с деконволюцией посредством предложенной кинетической модели.

 

Примечательно, что несмотря на одинаковый структурный тип цеолитоподобных образцов (MFI) кривые десорбции аммиака имеют разный вид. Более широкий пик десорбции для силикалита^–1 Sil^–1 означает более слабую связь аммиака с поверхностью, что в значительной степени компенсируется вариативностью адсорбции аммиака: большая часть адсорбционных мест тождественна, а, значит, и предэкспоненциальный множитель, который отвечает в большей степени именно за собственное колебание молекул на адсорбционной площадке будет также низким. Таким образом, молекулы со слабой подвижностью вдоль поверхности будут медленно десорбироваться, так как будут меньше подвержены столкновению друг с другом. Учитывая, что наблюдаемый порядок процесса десорбции превышает первый, это может означать, что и столкновения, необходимые для протекания реакций, также должны происходить реже. Столкновения необходимы молекулам адсорбата для обмена энергией и ее накопления для отрыва от поверхности – десорбции. Это согласуется с полученными кинетическими значениями параметров, приведенных в табл. 3. Несмотря на то, что изучение кислотности для цеолитов типа MFI требует более детального рассмотрения, можно считать, что десорбция аммиака по кинетическим уравнениям третьего порядка реализуется и более того, десорбирующиеся молекулы в тройных соударениях пытаются компенсировать дефицит энергии, необходимой для разрыва связи с поверхностью.

 

Таблица 3. Кинетические параметры предложенной модели для ZD-07031 и Sil^–1

Образец	Параметры модели
	ni	k0,i, M(1–n) с–1	ED,i, кДж/моль	Hmax,i, M	ni, 1/K	R2
ZD-07031 Кривая 1	3	(897 ± 2) × 109	107.6 ± 0.8	0.3090 ± 0.0013	(484.0 ± 1.1) × 109	0.9949
Sil–1 Кривая 2	3	23300 ± 270	48.9 ± 0.8	0.1450 ± 0.0008	2960 ± 40	0.9835

 

Описание кинетики десорбции аммиака с поверхности γ-Al₂O₃ также представляется довольно интересным. В предложенной серии образцов алюмооксидных носителей, несмотря на весьма различный вид зависимости скорости десорбции аммиака от температуры, разные скорости нагрева, различные значения интегральной кислотности, после обработки кривых кинетической моделью можно заметить, что основные кинетические параметры процесса десорбции имеют очень близкие значения, характеризующие сродство между анализируемыми образцами: все образцы имеют два характерных пика десорбции, у которых совпадают наблюдаемые порядки. Для первого пика (n = 3) характерные значения предэкспонента и энергии активации десорбции являются значения порядка 1 × 10⁹ M^–2с^–1 и 77 кДж/моль, для второго (n = 3) k₀ и E_D составляют соответственно 4 M^–1с^–1 и 28 кДж/моль (рис. 9, табл. 4).

 

Рис. 9. Кривые ТПД аммиака образцов A1 (а), Yar-700 (б) and Yar-800 (в) с деконволюцией посредством предложенной кинетической модели.

 

Таблица 4. Кинетические параметры предложенной модели для A1, Yar-700 и Yar-800

Образец	Параметры модели
	ni	k0,i, M(1–n) с–1	ED,i, кДж/моль	Hmax,i, M	ni, 1/K	R2
A1 Кривая 3 Кривая 4	3 2	(1.067 ± 0,009) × 109 4.22 ± 0.04	77.48 ± 0.03 28.04 ± 0.07	0.0460 ± 0.0011 0.0480 ± 0.0017	(1109 ± 9) × 104 0.967 ± 0.009	0.9846
Yar-700 Кривая 5 Кривая 6	3 2	(0,926 ± 0,003) × 109 4.37 ± 0.05	77.070 ± 0.016 28.220 ± 0.023	0.1430 ± 0.0015 0.1340 ± 0.0018	(1062 ± 4) × 105 3.275 ± 0.004	0.9995
Yar-800 Кривая 7 Кривая 8	3 2	(1,118 ± 0,003) × 109 3.89 ± 0.06	77.590 ± 0.011 27.70 ± 0.06	0.082 ± 0.006 0.138 ± 0.005	(4198 ± 4) × 104 2.991 ± 0.005	0.9830

 

На первый взгляд может показаться, что полученные данные противоречат термодинамическим принципам протекания реакций, поскольку для более высокотемпературной области по результатам моделирования наблюдается десорбция аммиака с центров с меньшим активационным барьером. В работах [35–37] можно найти объяснение подобному эффекту: для низкотемпературных пиков для оксида алюминия десорбцию, как правило, связывают с отрывом аммиака с Льюисовских кислотных центров [36], а высокотемпературные пики относят к десорбции с Бренстедовских центров [38]. При взаимодействии аммиака с Льюисовскими центрами молекула адсорбата координируется как атомом алюминия через неподеленную электронную пару азота, так и кислородами через водород (рис. 10а). Именно из-за более прочного связывания с поверхностью для десорбции аммиака с низкотемпературных центров требуется больше энергии для преодоления барьера. В работе [39] указывается, что адсорбированные молекулы на поверхности могут вести себя как двумерный газ, диффундируя вдоль поверхности с одного адсорбционного места на аналогичное другое. Для его реализации достаточным условием является меньшее значение энергии активации диффузии, которое достигается при достаточно близком расположении центров адсорбции (в этом случае считается, что энергия активации диффузии не превышает 30% от энергии активации десорбции), что для кислотных центров Льюиса легко реализуется (рис. 10б). Поскольку мы имеем кинетику десорбции второго и третьего порядка, то движение по поверхности катализатора молекул аммиака напрямую способствует их столкновению и накоплению энергии на колебательно-вращательных уровнях адсорбированных молекул, что необходимо для десорбции в газовую фазу.

 

Рис. 10. Особенности протекания поверхностной диффузии молекул аммиака и их десорбции с Льюисовских (а, б) и Бренстедовских (в, г) кислотных центров.

 

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

ЗАКЛЮЧЕНИЕ

В ходе работы был предложен универсальный подход для обработки результатов температурно-программируемой десорбции аммиака на основе кинетической модели Поляни–Вигнера. Одним из важных результатов моделирования процесса десорбции аммиака является применение кинетического уравнения для скорости реакции третьего порядка при описании экспериментальных кривых ТПД аммиака. Кроме того, предложен алгоритм, позволяющий в автоматическом режиме подбирать параметры кинетических уравнений для полученных экспериментальных зависимостей скорости десорбции аммиака h_кат(ммоль/(К мл)) от температуры T(K) Предложенный алгоритм апробирован на ряде образцов, в том числе цеолитных системах: титансиликалит^–1 и силикалит^–1. Оба образца обработаны по уравнению одной кинетической кривой третьего порядка с коэффициентом детерминации R² > 0.98.

По предложенной методике для нескольких различных образцов γ-Al₂O₃ удалось получить сходящийся результат при обработке выходных кривых ТПД NH₃ с помощью двух кинетических уравнений с коэффициентом детерминации R² > 0.98. Оказалось, что десорбцию аммиака на оксиде алюминия можно описать кинетическими уравнения третьего и второго порядков, имеющих предэкспоненциальные множители и энергии активации десорбции 1 × 10⁹ M^–2с^–1 и 77 кДж/моль и 4 M^–1с^–1 и 28 кДж/моль соответственно. Тем не менее полученные значения предэкспонентов отличаются на ±12% в ряду исследуемых алюмооксидных систем, что может быть связано с используемой моделью, где предэкспоненциальные множители при температурно-программируемых процессах не зависят от температуры. В будущих исследованиях нами будут предложены способы обработки таких кривых с уточнением математической модели и учетом влияния температуры на предэскпонент.

ФИНАНСИРОВАНИЕ

Работа частично выполнена при финансовой поддержке Министерства науки и высшего образования РФ в рамках государственного задания Института катализа СО РАН (проект FWUR-2024-0036) при выполнении синтетических и экспериментальных работ, а также при поддержке гранта Правительства Новосибирской области № гр-10 от 18.09.2023 были привлечены молодые сотрудники для проведения измерительных работ.

КОНФЛИКТЫ ИНТЕРЕСОВ

Авторы заявляют об отсутствии конфликта интересов, требующего раскрытия в данной статье.

About the authors

A. I. Lysikov

Boreskov Institute of Catalysis SB RAS (IC SB RAS); Novosibirsk State University (NSU)

Author for correspondence.
Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090; Novosibirsk, 630090

V. A. Vdovichenko

Boreskov Institute of Catalysis SB RAS (IC SB RAS); Novosibirsk State University (NSU)

Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090; Novosibirsk, 630090

E. E. Vorob’eva

Boreskov Institute of Catalysis SB RAS (IC SB RAS); Novosibirsk State University (NSU)

Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090; Novosibirsk, 630090

I. A. Shamanaeva

Boreskov Institute of Catalysis SB RAS (IC SB RAS); Novosibirsk State University (NSU)

Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090; Novosibirsk, 630090

E. V. Luzina

Boreskov Institute of Catalysis SB RAS (IC SB RAS); Novosibirsk State University (NSU)

Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090; Novosibirsk, 630090

L. V. Piryutko

Boreskov Institute of Catalysis SB RAS (IC SB RAS)

Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090

Zh. V. Veselovskaya

Boreskov Institute of Catalysis SB RAS (IC SB RAS); Novosibirsk State University (NSU)

Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090; Novosibirsk, 630090

E. V. Parkhomchuk

Boreskov Institute of Catalysis SB RAS (IC SB RAS); Novosibirsk State University (NSU)

Email: lyanig@catalysis.ru
Russian Federation, Novosibirsk, 630090; Novosibirsk, 630090

References

Supplementary files