CC BY
57
Инженерно-технические науки Engineering-technical sciences
УДК 519.711:544.431
ПРИМЕНЕНИЕ ДИСКРЕТНЫХ ЦЕПЕЙ МАРКОВА ДЛЯ МОДЕЛИРОВАНИЯ ХИМИЧЕСКОЙ КИНЕТИКИ
С.П. Бобков, Е.С. Бобкова, В.Е. Мизонов, А.П. Каткова
Ивановский государственный химико-технологический университет
В статье рассмотрены возможности применения элементов теории дискретных цепей Маркова для моделирования кинетики химических реакций. Показан физикохимический смысл математических объектов, описывающих дискретную Марковскую цепь. Рассмотрены примеры моделирования кинетики химических реакций.
Ключевые слова: химическая кинетика, цепи Маркова, моделирование.
В последнее время вырос интерес к использованию математического аппарата цепей Маркова для моделирования процессов в различных областях науки и техники. Известны успешные попытки использования данного подхода в изучении физических явлений (диффузия, теплопе-ренос) [1], технологических процессов (измельчение, смешивание, сепарация) [2], в биологии, экономике, социологии и пр. Причинами столь широкой популярности теории Марковских процессов могут являться не только желание исследователей изучать природные дискретные явления, минуя этап их формального описания с использованием абстрактных бесконечно малых и непрерывных величин, но и то обстоятельство, что среди различных видов случайных процессов Марковские процессы играют важнейшую роль, они часто встречаются на практике и имеют сравнительно простое математическое описание.
Данная статья посвящена изложению некоторых результатов использования математического аппарата дискретных цепей Маркова для моделирования химической кинетики.
Прежде всего, следует обосновать возможность использования предлагаемого подхода для описания процессов в химических системах. Рассмотрим следующие аспекты.
1. Согласно классической теории [3], любую химическую реакцию рассматривают, как последовательность столкновений между отдельными реагирующими частицами. Это обстоятельство позволяет рассматривать реакцию как дискретный процесс.
2. Чтобы столкновение частиц состоялось они, во-первых, должны одновременно оказаться в одной точке пространства, во-вторых, иметь необходимую пространственную ориентацию. Оба описанных события имеют ярко выраженный вероятностный характер. К этому следует добавить, что не все соударения приводят к реакции - частицы должны обладать достаточной энергией. Рассмотренные факторы определяют химическую реакцию как стохастический процесс.
3. Наконец, для химической реакции характерно влияние предшествующих событий (эффективных столкновений реаги-
рующих частиц) на последующие подобные события. Однако, можно утверждать, что данное влияние не распространяется по времени далее двух соседних состояний, и что химическая реакция - это процесс, для которого параметры текущего состояния в конкретный момент времени зависят от параметров предшествующего состояния в предыдущий момент времени.
Все три рассмотренных аспекта полностью характеризуют эволюцию химически реагирующей системы, как дискретную цепь Маркова и позволяют рекомендовать соответствующий математический аппарат для исследования и моделирования кинетики химических реакций.
Обычно Марковскую цепь изображают в виде ориентированного взвешенного графа, вершины которого соответствуют возможным состояниям системы А;, а дуги
- возможным переходам системы из состояния А; —> А]. Каждой дуге соответствует переходная вероятность
р||(к)=р[А|(к)/А|(1<-1)] - это условная веро-
ятность перехода системы на к-ом шаге в состояние А) при условии, что на предыдущем (к-1)-ом шаге система находилась в
СОСТОЯНИИ А;.
Полным описанием Марковской цепи служат матрицы переходных вероятностей:
Рп(к) р12(к) ... р1п(к)
Р21 (к) Р22(к) ••• Р2п(к)
Р(к)
(1)
Рш(к) Рп2(к) ••• Рт(к) Матрица ||Р(к)|| обладает важным свойством - сумма элементов каждой строки матрицы равна единице:
2>,(к) = 1 (2)
Н
Кроме того, на каждом шаге Марковская цепь характеризуется вектором-строкой вероятностей состояний |8(к)|, элемент которого Б;(к) характеризует вероятность того, что на к-ом шаге по времени система будет находиться в состоянии г.
ВД^к)^),..^)], (3)
где п - общее число состояний системы.
Очевидно, что для любого момента t сумма вероятностей всех состояний равна единице:
2>,(к)=і
(4)
Исходное состояние системы описывается вектором
|8(0)| = [81(0),82(0),...8п(0)], который обычно задается конкретными условиями процесса.
Моделирование эволюции системы производится путем определения значения вектора вероятностей состояний в каждый момент дискретного времени. Вычислять искомые значения вектора состояний удобнее всего умножением вектора-строки для предыдущего шага по времени на соответствующую матрицу вероятностей переходов:
|8(к+1)| = |8(к)|-||Р(к)||. (5)
Рассмотрим физико-химический смысл математических объектов, описывающих Марковскую цепь.
Сначала покажем, что элементы матрицы переходных вероятностей (1) характеризуют интенсивность протекания химической реакции в системе.
Эволюция любой системы, в том числе и химической, происходит вследствие воздействий на нее определенных движущих сил. В частности, химические реакции протекают под действием движущих сил термодинамического характера. Нетрудно предположить, что скорость и полнота, с которой протекает конкретная реакция, зависят, прежде всего, от интенсивности упомянутых движущих сил. Известно, что величиной, характеризующей интенсивность химической реакции, является ее скорость. Поэтому именно скорости протекающих реакций наиболее удобно принимать в качестве показателей интенсивности переходов системы между возможными состояниями. Скорости реакций принято определять, используя константы скоростей. Однако, известно, что этот показатель равен скорости реакции при единичной концентрации реагентов. Поэтому
интенсивностью переходов системы между состояниями следует считать скорость реакции, отнесенную к единичной концентрации компонентов. В свою очередь, интенсивность перехода тождественна плотности вероятности перехода, т.е. можно записать выражение для определения искомых элементов матрицы (1):
(к) = у1((к)-А1, (6)
где Уц(к) - скорость химической реакции, переводящей систему из состояния 1 в состояние у; А/ - шаг дискретного времени.
В общем случае скорость химической реакции зависит от концентраций реагентов, поэтому элементы матрицы (1) нужно рассчитывать на каждом шаге по времени. Другими словами, при моделировании кинетики химических реакций мы имеем дело с неоднородными цепями Маркова. Исключение составляют лишь химические реакции нулевого и первого порядка, где элементы матрицы переходных вероятностей будут постоянными во время всего процесса и Марковская цепь будет однородной.
Смысл элементов вектора вероятностей состояний (3) нетрудно определить из следующих соображений.
Основным параметром химического процесса, однозначно характеризующим состояние системы, следует считать концентрацию компонентов реагирующей системы. Можно предположить, что значения элементов вектора состояний должны указывать на наиболее вероятные значения концентраций реагентов в определенный момент времени. Однако, в этом случае составляющими вектора (3) должны быть не абсолютные значения концентраций реагентов, а их относительные значения, удовлетворяющие условию нормировки (4).
Алгоритм моделирования может быть следующим. Сначала задается шаг изменения модельного времени и устанавливаются начальные элементы вектора состояний 18(0)|. Далее, используя выражения (6) и (2) можно вычислить значения элементов Ру(0) матрицы переходов (1). Теперь, воспользовавшись выражением (5),
можно определить элементы вектора состояний, а, следовательно, и наиболее вероятные значения концентраций веществ, после первого шага по времени. Далее процесс вычислений повторяется. Относительно величины шага по времени, следует отметить, что выбранное значение не должно приводить к получению абсурдных результатов, т.е. элементы матрицы (1) не должны быть больше единицы или, тем более, быть отрицательными.
Рассмотрим примеры моделирования химических реакций с использованием цепей Маркова.
Пример 1. Возьмем простейший случай. Пусть необратимая химическая реакция протекает по уравнению:
2А+В -» 2У, где А и В - исходные реагенты, У - продукт реакции.
Будем считать, что вещества находятся в рассматриваемой системе в стехиометрических количествах.
Положим, что рассматриваемая реагирующая система имеет два предельных состояния. Будем считать состоянием 1 исходную смесь веществ А и В с начальными концентрациями, а состоянием 2 -систему, состоящую только из продукта реакции У. Тогда вектор вероятностей состояний (3) для данной реакции будет состоять из двух элементов:
№1 = ^(к), ЗД]. (7) Граф цепи Маркова системы приведен на рис. 1.
Рис.1. Граф цепи Маркова для примера 1
Для рассматриваемой реакции в любой момент времени должно выполняться следующее уравнение баланса:
2С 4 (к) + 2 Св(к) + ЪС¥(к) =
=2Сл(о;+2сд(о;+зсг(о;, (8)
где Сд(к), Св(к), Су(к) - текущие мольные концентрации веществ А; В, У, Сд(0), Св(0), Су(0) - начальные концентрации.
Выражение (8) позволит записать уравнения для связи элементов вектора со-
стоянии с концентрациями веществ:
Б, (к)
2СА(к) + 2Св(к) 2СА(0) + 2Св(0) + ЗСу(0)
в2(к):
ЗСу(к)
2СА (0) + 2СВ (0) + ЗС у (0)
(9)
Кроме того, можно записать:
СА(к) = 2Св(к). (10)
Выражения (9) и (10) позволят в процессе моделирования вычислить концентрации по значениям вектора состояний для каждого шага по времени.
Начальное значение вероятности перехода системы из состояния 1 в состояние 2 для рассматриваемой реакции может быть вычислено так:
Ріг(0) = уі2(0)-Аі =
=к12[СА(0)]2 [Св(0)]-А1, (11)
где Уі2(0) - начальная скорость реакции; кі2 - константа скорости реакции.
Вероятность р21 = 0 по условиям реакции, остальные элементы матрицы (1) можно вычислить с использованием выражения (2). На последующих шагах вероятности будут рассчитываться по аналогичным выражениям, но с использованием значений концентраций на соответствующем шаге.
На рис. 2 приведены результаты моделирования реакции за 15 шагов модельного времени. В качестве исходных данных были взяты следующие значения: СА(0) = 2; Св(0) = 1; Су(0) = 0; к12 = 4-103 л2/(моль2-с); А1 = 1-10"5 с.
Рис. 2. Кинетика реакции по данным примера 1
Рассмотренная химическая реакция протекает в одном направлении. Моделирующая реакцию цепь Маркова (рис.1) является поглощающей, т.к. имеет состояние, из которого выход не возможен. В соответствии с теорией цепей Маркова установившимся (стационарным) состоянием для такой цепи будет состояние 2, что полностью согласуется с основными положениями химической кинетики.
Пример 2. Рассмотрим обратимую химическую реакцию, протекающую по уравнению:
А<-> 2У,
Будем считать предельным состоянием 1 такое, при котором в системе на-
ходятся только частицы вещества А, а состоянием 2 - когда вся система состоит только из частиц вещества У.
Соответствующий граф цепи Маркова приведен на рис.З.
12
Г21
Рис.З. Граф цепи Маркова для примера 2
Должно выполняться уравнение баланса:
2СА(к) + Су(к) = 2Са(0) + Су(0) , (12)
где Сд, и Су, - мольные концентрации веществ А и У.
Будем считать предельным состоянием 1 такое, при котором в системе находятся только частицы вещества А, а состоянием 2 - когда вся система состоит только из частиц вещества У.
Запишем выражения для компонентов вектора состояний, исходя из (12): 2СА(к)
БДк) = $2(к)=
2СА(0) + СУ(0) ’
су(к)
2Са(0) + Су(0)
(13)
Вероятности переходов, определяемые скоростями соответствующих реакций, определятся таким образом:
Р12О) = к^ и
р21(0 = к21[Су(0]А1, (14)
где к!2 - константа скорости реакции, переводящей систему из состояния 1 в состояние 2; к21 - константа скорости обратной реакции.
Остальные элементы матрицы (1) вычисляются с использованием условия (2).
Результаты моделирования кинетики реакции приведены на рис. 4. Для расчета были взяты следующие исходные данные: СА(0)=0,5; Су(0)= 1; к12 = 0,5 с'1; к21 = 2 л/(моль-с); А1 = 2,5-10'2 с.
Моделирующая реакцию из примера
2 цепь Маркова является эргодической, то есть все возможные состояния достижимы. Согласно теории цепей Маркова такая цепь должна иметь установившееся (стационарное) состояние, при котором элементы вектора (3) будут иметь постоянные значения. Более того, установив-
шиеся значения элементов вектора (3) не будут зависеть от его начальных значений. Это подтверждается данными, представленными на рис.5, где приведены результаты моделирования реакции при тех же значениях кинетических констант, но с другими начальными концентрациями.
Рис. 5. Кинетика реакции при различных значениях исходных концентраций
С другой стороны, с точки зрения химической кинетики предельным состоянием системы при обратимой реакции будет равновесное состояние. Также известно, что если вещества в обратимой реакции присутствуют в стехиометрических пропорциях, то равновесные концентрации не зависят от исходных концентраций. Таким образом, результаты моделирования строго соответствуют представлениям химической кинетики и не противоречат классическому подходу.
Пример 3. В качестве следующего примера рассмотрим химическую реакцию, состоящую из двух стадий, причем первая из них является необратимой: 2А+В -» 2М 2Х+2У
Данная система будет иметь три предельных состояния. Будем считать состоянием 1 - систему, полностью состоящую из веществ А и В, состоянием 2
- систему, состоящую только из молекул промежуточного продукта М, состоянием
3 - систему, где содержатся только продукты реакции X и У.
Граф цепи Маркова для данной системы приведен на рис.6.
Рис.6. Граф цепи Маркова для примера 3
Для данной реакции справедливо следующее уравнение баланса:
4С А (к) + 4С в (к) + 6С м (к) + ЗС х (к) + ЗС у (к) = 4С А (0) + 4СВ (0) + 6СМ (0) + ЗСХ (0) + ЗС у (0) = £ С(0)
(15)
где, как и в предыдущих примерах С; -концентрации участвующих в реакции веществ.
Запишем выражения для компонентов вектора состояний, который в данном примере будет иметь три элемента:
4СА(к) + 4Св(к)
БДк) = 82(к) = $з(к)=
Ес(°)
6См(к)
Ес(°) ’
ЗСх(к) + ЗСу(к)
(16)
Ес(°)
Зависимости для определения элементов матрицы переходных вероятностей, вычисляемые через скорости реакций, запишутся так:
Р120) = к12[СА0)Г [Св0)]-А1,
ргзО) = к23Л1, (17)
РзгС) = к32[Сх(0 [СУ(0] А1.
Еще три элемента матрицы будут нулевыми, поскольку непосредственные переходы между соответствующими состояниями невозможны:
Р13(0 = Р21(0 = РЭ1(0 = 0. (18)
Остальные элементы вычисляются с использованием условия (2).
Результаты моделирования кинетики реакции приведены на рис. 7. Были взяты следующие исходные данные: Са(0)=2; Св(0)=1;
См(0)=0; Сх(0)=0;
Су(0)= 0; к^ = 0,075 л2/(моль2-с);
к2з = 0,05 с'1; к32 = 1 л/(моль-с); А1 = 0,5 с.
1,5 1
0,5 О
Рис. 7. Кинетика реакции по данным примера 3
На рис. 7. показаны результаты моделирования только за первые 16 шагов. Анализ дальнейшего протекания процесса показывает, что концентрация веществ А и В стремится к нулю, а концентрации реагентов М, X и У - стремятся к постоянным равновесным значениям.
Проведенные исследования показали возможность применения дискретных вероятностных моделей для описания эволюции химически реагирующих систем.
В начале статьи уже говорилось об универсальности и эффективности применения теории цепей Маркова для моделирования объектов различной приро-
ды. Можно еще раз констатировать поразительную универсальность этого математического инструмента, а также физическую ясность и логическую стройность получаемых моделей. Все это дает возможность рекомендовать данный подход для исследования и анализа процессов в сложных химических системах.
ЛИТЕРАТУРА
1. Мизонов В.Е. и др. // Известия вузов «Химия и хим. технол.». 2005. Т. 48, № 1, С.122.
2. Berthiaux Н., Mizonov V., Zhukov V. // Powder Technology’ 157 (2005) 128-137.
3. Эмануэль H. М. Курс химической кинетики. М.: Высш. Школа. 1984. 464 с.
Рукопись поступала в редакцию 21.03.12.
MARKOV’S DISCRETE CHAINS’ APPLICATION FOR CHEMICAL KINETICS MODELLING
S. Bobkov, E. Bobkova, V. Mizonov, A. Katkova
In the article the possibilities of Markov’s discrete chains theory elements application for chemical reactions kinetics modelling are considered. Physical and chemical sense of mathematical objects describing Markov’s discrete chain is shown. The examples of chemical reactions kinetics modelling are considered.
Key words: chemical kinetics, Markov’s discrete chains, modelling.
С, моль/л
