Химическая кинетика реакций¶

Изучение динамических процессов в химии на макроскопическом масштабе часто относится к химической кинетике.

Эта дисциплина затрагивает закономерности протекания химических реакций во времени, их зависимость от внешних условий и механизмы химических превращений.

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

Pkg.add("DifferentialEquations")

В рамках небольшого проекта по изучению простых моделей химической кинетики мы пройдем четыре примера:

Абстрактная химическая реакция для трех веществ как модель DifferentialEquations,
Переход на библиотеку Catalyst для упрощения синтаксиса модели,
Переход к реалистичным названиям веществ (окисление сероводорода),
Добавление уравнений в модель: учет влияния температуры на скорость протекания реакции (модель Аррениуса).

Реакции как дифференциальные уравнения¶

В первом примере мы рассчитываем последовательную химическую реакцию, где вещество A превращается в B, а затем B превращается в C. Процесс описывается системой обыкновенных дифференциальных уравнений (ОДУ), где скорость изменения концентрации каждого вещества зависит от констант реакции $k_1$ и $k_2$. Уравнения выглядят так:

$\frac{dA}{dt} = -k_1·A$ (вещество A расходуется),
$\frac{dB}{dt} = k_1 \cdot A - k_2 \cdot B$ (вещество B образуется из A, но распадается в C),
$\frac{dC}{dt} = k_2 \cdot B$ (вещество C накапливается).

Такая модель позволяет предсказать, как со временем изменяются концентрации всех участников реакции.

using DifferentialEquations

# Определяем систему ОДУ: dA/dt = -k1*A, dB/dt = k1*A - k2*B, dC/dt = k2*B
function kinetics!(du, u, p, t)
    A, B, C = u
    k1, k2 = p
    du[1] = -k1 * A          # dA/dt
    du[2] = k1 * A - k2 * B  # dB/dt
    du[3] = k2 * B           # dC/dt
end

# Начальные условия и параметры
u0 = [1.0, 0.0, 0.0]  # [A0, B0, C0]
p = [0.5, 0.2]        # константы скорости k1, k2
tspan = (0.0, 10.0)   # интервал времени

# Решаем систему
prob = ODEProblem(kinetics!, u0, tspan, p)
sol = DifferentialEquations.solve(prob, Tsit5())

# Визуализация
plot(sol, label=["A" "B" "C"], xlabel="Время", ylabel="Концентрация", lw=2)

Такие графики наблюдаются при расчете системы с параметрами $k_1$ = 0.5, $k_2$ = 0.2 и при начальной концентрации ( $[A]_0$ = 1, $[B]_0$ = 0, $[C]_0$ = 0 ). Для ее решения задачи использовался стандартный метод численного интегрирования Tsit5().

Мы показали, как решить типичный пример для лекций по химической кинетике средней сложности. Теперь изучим, как упростить его описание и сделать всю методологию более пригодной для решения реалистичных задач. Перейдем к описанию этой реакции при помощи команд для предметной области (Domain Specific Language, DSL), представленном в библиотеке Catalyst.

Применение библиотеки Catalyst¶

Зададим ту же самую последовательную реакцию $A \to B \to C$ при помощи более короткого описания, заданного командами из библиотеки Catalyst.jl (сперва установив эту библиотеку).

Pkg.add("Catalyst")

Библиотека содержит много средств, упрощающих решение задач химической кинетики. Например, символьный API для создания цепочек реакций (@variables, @species, Reaction, ReactionSystem), инструменты анализа реакционных сетей (классы связности, дефицит, обратимость, обнаружение законов сохранения) и визуализация сетей через GraphMakie и CairoMakie, конвертация реакционных систем в ОДУ, стохастические ДУ или прыжковые процессы, интеграция с пакетами для расчетов на GPU и т.д.

Но основной прием, которым мы воспользуемся – это команда Catalyst DSL (языка для описания реакций), базовым блоком является конструкция @reaction_network begin .. end.

rn = @reaction_network begin
    k1, A --> B
    k2, B --> C
end

Каждая реакция задается отдельной строчкой внутри блока begin..end. Первым элементом описания является скорость протекания реакции, затем отношение входящих в реакцию символов, обозначающих химические вещества. Из технических особенностей нам нужно будет задать следующие дополнительные параметры:

начальные условия,
количество шагов (ограниченный вектор времени),
значения параметров.

Полученную систему можно оформить как ОДУ (объект ODEProblem), рассчитать при помощи команды solve и вывести на график при помощи команды plot.

using Catalyst, DifferentialEquations
using Catalyst: solve # Неважно, из какой библиотеки брать эту функцию

# Определяем реакционную сеть
rn = @reaction_network begin
    k1, A --> B
    k2, B --> C
end

# Задаем параметры и начальные условия
p = [:k1 => 0.5, :k2 => 0.2]           # константы скорости
u0 = [:A => 1.0, :B => 0.0, :C => 0.0] # начальные концентрации
tspan = (0.0, 10.0)                    # интервал времени

# Создаем и решаем задачу ОДУ
prob = ODEProblem(rn, u0, tspan, p)
sol = solve(prob, Tsit5())

# Визуализация
plot(sol, label=["A" "B" "C"], xlabel="Время", ylabel="Концентрация", lw=2)

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

Реалистичные названия переменных¶

Сделаем полшага в сторону более удобного описания реакций. Чтобы изобразить реалистичную реакцию нужно всего лишь поменять названия переменных.

Рассмотрим реакцию окисления сероводорода ( $H_2S$ ) сначала до серы ($S$), затем до диоксида серы ($SO_2$) в присутствии кислорода. Упрощенно: ( $H_2S \to S \to SO_2$). Это учебный пример, где все стадии имеют кинетику первого порядка.

using Catalyst, DifferentialEquations
using Catalyst: solve

# Определяем реакционную сеть
rn = @reaction_network begin
    k1, H2S --> S      # H₂S → S
    k2, S --> SO2      # S → SO₂
end

# Параметры и начальные условия
p = [:k1 => 0.3, :k2 => 0.1]               # константы скорости (условные, 1/с)
u0 = [:H2S => 1.0, :S => 0.0, :SO2 => 0.0] # начальные концентрации (моль/л)
tspan = (0.0, 20.0)                        # интервал времени (с)

# Решаем задачу ОДУ
prob = ODEProblem(rn, u0, tspan, p)
sol = solve(prob, Tsit5())

# Визуализация
plot(sol, label=["H₂S" "S" "SO₂"], xlabel="Время (с)", ylabel="Концентрация (моль/л)", lw=2 )

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

Дополняем описание реакции моделю Аррениуса¶

Напоследок изучим пример, в котором показано, как дополнить пример с окислением $H_2S$, добавив в неё учёт дополнительных зависимостей, а именно зависимость скорости протекания реакций от температуры.

В первом приближении, химические реакции между веществами происходят только в результате столкновения молекул этих веществ (простая модель столкновений). Чтобы преодолеть энергетический барьер, молекулы должны обладать некой минимальной энергией активации, отсюда возникает модель $k(T) = A \cdot e^{-E_a/RT}$, где $A$ характеризует частоту столкновений, $T$ – температура, $R$ – универсальная газовая постоянная, $E_a$ – энергия активации. Количество молекул, входящих в химическую реакцию, определяется из распределения Больцмана как пропорциональное $exp(-\frac{E_a}{RT})$.

Теперь у нашей реакции реалистичная скорость протекания, значительно больше скорости в прошлом примере. Новый интервал наблюдений отражается на параметреtspan.

using Catalyst, DifferentialEquations
using Catalyst: solve
using Plots.PlotMeasures

# Функции для констант скорости (например, уравнение Аррениуса)
k_1(T) = 1e5 * exp(-5000 / (8.314 * T))  # A * exp(-Ea/RT), условные параметры
k_2(T) = 2e4 * exp(-3000 / (8.314 * T))  # разные A и Ea

# Реакционная сеть с функциональными константами
rn = @reaction_network begin
    k_1(T), H2S --> S    # k1 зависит от температуры T
    k_2(T), S --> SO2    # k2 тоже зависит от T
end

# Параметры и начальные условия
tspan = (0.0, 0.0002)
p = [:T => 273.0]
u0 = [:H2S => 1.0, :S => 0.0, :SO2 => 0.0]

# Решаем ОДУ
prob1 = ODEProblem(rn, u0, tspan, p)
sol1 = solve(prob1, Tsit5())

# График
p1 = plot(sol1, label=["H₂S" "S" "SO₂"], xlabel="Время (с)", ylabel="Концентрация (моль/л)", lw=2,
          title="Реакция окисления H₂S при 273 K", titlefont=font(11), guidefont=font(8),
          bottommargin=30px)

# Параметры и начальные условия
p = [:T => 400.0]
u0 = [:H2S => 1.0, :S => 0.0, :SO2 => 0.0]
prob2 = ODEProblem(rn, u0, tspan, p)
sol2 = solve(prob2, Tsit5())

p2 = plot(sol2, label=[false false false], xlabel="Время (с)", ylabel="Концентрация (моль/л)", lw=2,
          title="Реакция окисления H₂S при 400 K", titlefont=font(11), guidefont=font(8))

plot(p1, p2, layout=(2,1))

Заключение¶

Эти примеры демонстрируют два подхода к моделированию химической кинетики в Engee с применением языка Julia. Первый вариант — ручное задание системы дифференциальных уравнений — дает полный контроль над расчетами и полезен для понимания основ. Второй способ, с использованием Catalyst.jl, позволяет описывать реакции в лаконичном формате, автоматически генерируя уравнения. Это особенно удобно для сложных систем, где важна скорость разработки и минимизация ошибок.

Кроме этого мы изучили, как в модель добавить дополнительные уравнения и сделать всю методологию более пригодной для решения реальных задач химической кинетики.