CC BY
23
MS С 80А30
ИССЛЕДОВАНИЕ КРИТИЧЕСКОЙ ПОВЕРХНОСТИ СТОХАСТИЧЕСКОЙ МОДЕЛИ ХИМИЧЕСКОЙ КИНЕТИКИ БИНАРНОЙ АВТОКАТАЛИТИЧЕСКОЙ РЕАКЦИИ. СИЛЬНО АСИММЕТРИЧНЫЙ СЛУЧАЙ
Аннотация. Изучается стохастическая модель бинарной автокаталитической химической реакции на основе уравнений химической кинетики со стохастически возмущенными параметрами, в которой возмущение описывается обобщенным случайным процессом белохх) шума. Исследуется стационарная плотность распределения первохх) порядка елучайжнх) процесса, описывающих) эволюцию концентрации оджнх) из компонентов реакции. С качественной точки зрения исследуется критическая поверхность в трехпараметричееком пространстве термодинамических параметров состояния смеси, при пересечении которой происходит индуцированный шумом фазовый переход в стохастической системе, то есть перестройка унимодальной плотности распределения в бимодальную. Развита теория возмущений в случае сильнохх) превосходства концентрации оджнх) катализатора на концентрацией другого. Она позволяет вычислить приближенно фазовые диаграммы при малых значениях параметра, характеризующих) такую асимметрию.
Ключевые слова: химическая реакция, стохастическая модель, фазовый переход, бифуркация, плотность распределения, критическая поверхность, теория возмущений.
1. Введение. В работе |1| была рассмотрена стохастическая модель химической кинетики, исследовавшаяся ранее в |2, 3| и описывающая протекание автокаталитических химических реакций определенного тина между двумя реагентами с учетом их термодинамических флуктуаций, В рамках этой модели эволюция во времени относительной доли Хг концентраций реагентов описывается посредством стохастического дифференциального уравнения
со свободными параметрами модели А € К, а > 0 а Е [0,1] и с мультипликативным белым шумом Стохастический дифференциал dwt в этом уравнении понимается
но Стратоновичу. Таким образом, эволюция системы представляет собой неоднородный диффузионный марковский сну чайный процесс в пространстве относительных концентраций. Этот процесс обладает финальной плотностью распределения р(х) для случайной величины Хг, которая, в зависимости от значений трех параметров (А, а2, а) модели является либо унимодальной, либо бимодальной. Финальная плотность определяется явной формулой |1|
Фам Минь Туан, Ю.П. Вирченко
Белгородский государственный университет, ул. Победы, 85, Белгород, 308015, Россия, e-mail: [email protected]
dxt = [а — xt + Axt (1 — xt)] dt + axt(1 — xt)dwt
(1)
где К-в(•) - модифицированная функция Бесселя второго рода с показателем (—в) и в = 2(2а + Л — 1)/а2, Переход от унимодальной плотности к бимодальной происходит на некотором двумерном многообразии Е в пространстве параметров (Л, а2, а), которое называется критической поверхностью. Такой переход физически трактуется как динамический фазовый переход в смеси химических реагентов. Уравнение дня критической поверхности находится из условия равенства нулю производных р'(х) = 0, что эквивалентно уравнению
а2
Б(х) = а — х + А;г(1 — х) ——#(1 — — 2х) = 0 , (3)
и р''(х) = 0. При этом совместное решение такой системы уравнений для концентрации х обязано находиться на интервале [0,1].
В работе [1] было показано, что критическая поверхность Е в пространстве наборов (Л, а2, а)
фазовый переход, определяется уравнением
Р(Л, а2, е) = Л4 + Л2( 1 — 5а2 — а4/2^ — Ле(9а4 + 18а2 — 4Л2) — 4а2 (1 — а2/4) 3 — 27а4е2 = 0 ,
(4)
где е = а — 1/2 е [—1/2,1/2]. Кроме того, в работе [1] было показано, что параметры критической поверхности, удовлетворяющие этому уравнению, должны, дополнительно, обладать свойством
\Г(\ «2 < 1 га сЛ - + 2а2) + 36 еа2
Если ввести функции С±(Л,а2,е) = (1 ± С(Л,а2,е))А
А = 4Л2 + 3а4 — 12а2 ,
то это условие перепишется в более простом виде
С+(Л,а2,е)С_(Л,а2,е) > 0,
так, что допустимая область дня расположения точек бифуркации представляет собой объединение областей (5±, границами которых являются соответственно гиперболы
С±(Л, а2, е) = (2а2 + 1 ± 2Л)2 — (а2 + 2(4 ^ 9е))2 + 4(4 ^ 9е)2 — 1 = 0 .
На это обстоятельство пе было обращено должного внимания в предшествующих рабо-
Е
одпой связной частью поверхности, определяемой уравнением (4).
В настоящем сообщении мы разрабатываем метод вычисления фазовых диаграмм термодинамически равновесного состояния системы химических реагентов при значе-а
то есть к 0 и 1. Ввиду симметрии поверхности относительно одновременной замены а ^ 1 — аиЛ ^ — Л, нам достаточно изучить только окрестность точки а = 0.
2. Вырожденный случай. Рассмотрим сечение критической поверхности в вырожденных случаях при е = ^1/2,
Теорема 1. Сечение поверхности Е с полуплоскостью {(Л,а2,е = —1/2) : Л е М,а2 > 0} ({(Л,а2,е = 1/2) : Л е М,а2 > 0}) состоит из (двукратной) полупрямой Л0_} = 1 + а2/2 (Л0+) = —(1 + а2/2)), а2 > 1 и кривых Л*-) = (—а2 ± 4а)/2 (Л*+) = (а2 ± 4а)/2), а2 > 1, для которых общими с цилиндром точками являются Л = ±3/2 а2 = 1 □ Уравнение (4) при значениях е = ±1/2
Л4 + Л2( 1 — 5а2 — а4/^ Т Л(9а4 + 18а2 — 4Л2)/2 — 4а2 ( 1 — а2/4)3 — 27а4/4 = 0
допускает явное разложение па множители
(Л ± (1 + а2/2))2(Л2 ^ а2Л + а2(а2 — 16)/^ =0
и поэтому - точное алгебраическое выражение всех решений, которые являются вещественными, представляется следующим образом.
При е = —1/2 список решений состоит из двукратного решения Л0_) = 1 + а2/2, то есть представляет собой «удвоенную» прямую, и решений
Л*-) = (—а2 ± 4а)/2 ,
которые представляет собой кривую, составленную из двух но.иубесконечпых ветвей а2 > 0 Л = 3/2 а2 = 1
Л
изменяется от —то до 2 (в точке Л = 2 функция Л* достигается максимума).
Точно также, при е = 1/2 список решений состоит из «удвоенной» прямой Л0+) =
—(1 + а2/2)
Л*+) = (+а2 ± 4а)/2 .
Л
[—2, то) ((—3/2,1) - точка склейки двух полубесконечных ветвей). ■
Естественно, что в обоих случаях кривые, определяемые теоремой, находятся вне эллипса А = 4Л2 + 3а4 — 12а2 = 0, Это проверяется подстановкой выражения для Л*±) в уравнение для эллипса, что приводит к неравенству А = 4а2(а ± 1)2 > 0. Таким
а = 0, Л = 0
а = 1,Л = ^3/2.
3. Теория возмущений. Нашей задачей является изучение критической поверхности вблизи предельных значений е = ±1/2. Принимая во внимание симметрию критической поверхности, мы изучим только случай е = -1/2 (а = 0).
Покажем сначала, что изучение сечений поверхности, то есть кривых А(а2), при по-
а
а
Рис. 1. В плоскости (Л, а2) при £ = -1/2 изображены: 1 — эллипс, ограничивающий расположение кривой раздела фаз; 2 — параболы л! \ «сшитые» в точке А = (3/2,1) и расположенные полностью
вне эллипса; 3 — сдвоенная прямая Л0
Произведем подстановку е = -1/2 + а в уравнение (4) Р(А, а2, -1/2 + а) = 0. Будем искать решение в виде разложения А(а2) = А0(а2) + а1/2А1 (а2) + аА2(а2) + 0(а3/2), где А0(а2) = 1 + а2/2 - одна из ветвей решения этого уравнения (Ао-^) при а = 0. В результате, получаем
Р(Л0(а2),(72,-1/2) = 0, = 0,
\дА/ Ао,а2,-1/2
(так как прямая А0 (а2) двойная) и
р = а(т)
\9е У Ао,-1/2
+ [аЛ? + 2а3/2Л1Л2] +
2 V дА У Ао,а2,-1/2
Ао,а2 ,-1/2
\д\д£ )
Ао ,а2 ,-1/2
или, после подстановки явных выражений дня частных производных, имеем Р = —а(Л0(9а4 + 18а2 — 4Л0) — 27а4) + (а2 — 1)2[аЛ? + 2а3/2Л1Л2] + 2(а2 + 1)а3/2Л3 —
— 6а3/2(а2 — 1)(а2 + 2)Л1.
Слагаемые порядка а дают выражение Л2(а2) = 4(а2 — 1) при а2 = 1. Баланс слагаемых
а3/2
(а2 — 1)2Л2 + (а2 + 1)Л? — 3(а2 — 1)(а2 + 2) = 0,
что с учетом явного выражения для Л1(а2), дает
2- а2
Л2(а2) =
а2 1
Мы видим, что в нервом приближении появляется расщепление А^сг2) = ±2-\Дт2 — 1, однако, уже в следующем приближении появляется малый знаменатель вблизи точ-
а2 = 1
выбрав в качестве нулевого приближения другую ветвь (Л*-)) из числа тех, которые слипаются в точке (3/2,1), Л0(а2) = (—а2 + 4а)/2, то для Л1(а2) получим уравнение 4а3(а — 1)3 + (а2 — 1)2Л2( а2) = 0, у которого нет решения при а2 > 1 (интересующая нас часть кривой, которая является границей раздана фаз. расположена в полуплоскости с
а2 > 1)
В связи с описанной выше невозможностью построения теории возмущений посредством прямого разложения в ряд но нолуцелым степеням, мы применим другой подход дня изучения поведения границы раздана фаз. Совершим подстановку
Л =(3 + и + г>)/2 , а2 = 1 + V — и (5)
так, что уравнение (4) принимает вид
и2^2 + 4и) + а{36и2 + 6^(и + V — 3) — 4(и3 + V3)} — 27а2(1 + V — и)2 = 0 ,
а уравнение эллипса в этих переменных, и2 + V2 — uv + 3и = 0. Теперь, чтобы произвести расщепление, подберем правильный масштаб изменения переменных и, V в окрестности точки (0, 0) пр и а ^ 0. Для этого произведем замену и ^ а"и, V ^ а6^, оде а > 0 Ь > 0 постоянные показатели
а2(а+6)и^2 + 4а3аи3 + 36а2а+1и2 + 6а2а+6+^и2 + 6а26+а+1^2 — 18аа+6+1^—
—4а3а+1и3 — 4а36+^3 — 27а2(1 + V — и)2 = 0 .
1 Заметим, что на наличие сложности при построении теории возмущений указывает уже то. что при
£ = —1/2 имеется пересечен не осп а2 = 0, а пр и £ = —1/2 такого пересечения нет, так как уравнение для точек пересечения Л4 + Л2 + 4£А3 = 0 имеет единственное решение Л = 0. Следовательно части прямых
±(а2/2 + 1), которые расположены при а2 < 1 должны бифуркационным образом деформироваться так, чтобы их точки пересечения с а2 =0 исчезли.
Дня нахождения правильного вида асимптотики, нужно разбить все имеющиеся наименьшие показатели степеней 2а + 26, 3а 2а + 1 а + 6 +1 36 +1, 2 у перемен ной а на группы так, чтобы имелась группа, по крайней мере, из двух слагаемых, с одинаковыми минимальными значениями степеней. Это возможно сделать единственным способом а = 2/3, 6 = 1/3. В результате, приравнивая коэффициент при минимальной степени а
4(и3 - V3 - 27) + («V - 9)2 = 0 . (6)
При V — то неявная функция м(у) не может быть ограниченной, так как, в противном случае, получаем противоречие. Поэтому имеются асимптотики кривой «(V) — то при V — то Возможны следующие случаи: и ~ V2, и ~ V1/2. Они следуют из сравнения слагаемых со старшими степенями и^2, и3 и V3 и из соображений компенсации слагаемых с наибольшими показателями степени. Рассмотрим оба случая. Положим и = av2 + /. Подстановка с удержанием главных степеней V6, которые должны исчезать, дает условие а = -1/4. После этого находим / = -2^ + о^-1), как главное
1/2
жений симметрии V = и2/4 + 2/и + о(и-1).
Общий вид кривой, определяемой уравнением (6), на плоскости (и^) показан на рис. 2, где расположены сверху вниз вдоль оси V обе ее связные компоненты и кубическая парабола, ограничивающая кривую снизу. Такой ее вид обусловлен следующим качественным анализом.
Кривая переходит сама в себя при растяжении а-2/3 вдоль о си и н а-1/3 вдоль о си V. Кроме того, она расположена выше кубической параболы V3 - и3 + 27 = 0 на рте. 2,2) Заметим, что кривая симметрична относительно замены и ^ -V.
Проверим кривую (6) па самопересечение. В случае наличия самопересечения доли
ванном V,
ф(и) = 4и3 + и^2 - 4v3 - 18^ - 27 = 4(и3 - V3 - 27) + (^ - 9)2 = 0 .
Используя алгоритм Евклида для этого полинома и его производной по и, ^'(и)/2 = 6и2 + uv2 - 9^ ^^^^^^^ ^^^^^^^ ^^^^^^^^^^ применения алгоритма) в виде (V3 - 27)3/81. Следовательно, единственная точка, в которой самопересечение возможно, определяется точкой при V = 3.
Далее, кривая имеет точки пересечения с диагональю и = -V. Наличие таких точек позволяет сделать заключение о ее двухсвязности. Совершая подстановку и = -V в уравнение, получаем (и - 1)(и + 3)3 = 0, то есть имеется трехкратный корень и = -3, V = 3, н поэтому точка (-3, 3) является особой, а также корень и = 1, V = -1.
"Заметим, что ио = 9 ± 2\/27+" о3 - и3. Тогда при обращении в нуль подкоренного выражения происходит смена знака в формуле, определяющей кривую. При этом та ветвь кривой, для которой такое положение реализуется, касается в точках указанного изменения знака кубической параболы.
Рис. 2. В плоскости (и, у) изображены: 1 - ветвь двухсвязной кривой (6), содержащая «касп» и описывающая участок границы раздела фаз; 2 - нефизическая ветвь кривой (6); 3 - кубическая парабола, касающаяся нефизической ветви в точках и = 3[(1 ± \/5)/2]1//3, которые находятся как совместные решения системы иу = 9, и3 — у3 = 27.
В точке (1, — 1) кривая пересекает трансверсально диагональ V = —и, так как в этой точке невозможно ее самопересечение. Это указывает на то, что кривая двухсвязна. Если бы она была односвязна, то дуги кривой, связанные с этим пересечением диагонали, при их продолжении, ввиду симметрии кривой относительно диагонали, должны соединиться с продолжением дуги кривой, проходящей через точку (—3, 3), то есть иметь, по крайней мере, две точки таких соединений. Если соединение состоит в том, что указанные продолжения кривой переходят друг в друга, то кривая должна быть заключена в ограниченной области плоскости, а это, как указано выше, невозможно. Если же соединение происходит посредством самопересечения — пересечения двух ветвей кривой, то это невозможно, так как такое самопересечения может реализоваться только в точке (—3, 3) на диагонали. Таким образом, отрезок дуги, пересекающий диагональ в (1, — 1),
(—3, 3)
(—3, 3)
кривой па этой ее связной части в окрестности этой особой точки и тин возникающей при этом особенности. С этой долыо перейдем к полярным координатам с центром в точке (-3, 3), то есть произведем в уравнен ни замену u ^ u —3 и v ^ v+3. В результате, получаем следующее
u2v2 + (u — v)(4(u2 + v2) + 10uv) = 27(u + v)2 ,
а затем положим u = pcos p, v = p sin p. Тогда, после деления уравнения на р2, имеем
^р2 sin2 2<р + р(cos р - sin р)(4 + 5 sin 2р) - 27(1 + sin 2р) = 0 . (7)
Квадратное уравнение дает два решения
2 г з/21
р = —j— (sin р - cos p)(4 + 5 sin 2p) ± (2(2 + sin 2p)) ' , (8)
sin2 2p L -I
(+)
ввиду отрицательности последнего слагаемого в уравнении, выражение в правой части (8) положительно. Так как решение с p(p) > 0 единственно, то в исследуемой точке реализуется особенность тина «касп», так как при наличии самопересечения в этой точке, должно быть, по меньшей мере, два решения p(p). Заметим, что характер особенности при р — то удается выяснить только при учете слагаемого, пропорционального р2, так как свободный член ~ (1 + sin 2p) уравнения обращается в ноль при p = —п/4, 3п/4.
Найденные значения углов определяют возможные направления подхода кривой к особой точке p — —п/4, 3п/4 (отвечающей кратному решению) при р — 0. Для исследования того, какое из них реализуется на самом деле, введем отклонения ф = p + п/4 и ф = p — 3п/4 и изучим поведение решений уравнения в ок рестности угла ф = 0 в обоих случаях. Так как
sin 2р = -1 + 2ф'2 + о(ф2), cos р - sin р = =рл/2 ± ф2 + о(ф2),
где верхние знаки в этой формуле соответствуют p = 3п/4, а нижние — p = — п/4, то
ф - 0
ф - 0
р2 ± 4л/2р - 216V'2 + о(ф2) = 0 .
Отсюда следует, что при знаке (+) получается р(р) = 21\f2ф2 + о(ф2). При знаке (—) имеется только одна возможность получить решение p(p), обращающееся в нуль при ф = 0, Но оно, ввиду требования неотрицательности p(p), дает только лишь изоли-
ф=0
кривая p(p) может подходить к кратной точке толь ко под углом p = 3п/4, то есть кратная точка не является точкой самопересечения кривой. Следовательно, в кратной точке реализуется особенность кривой в виде «касна».
Вернемся теперь к исходным переменным Л, а2 на изученной нами связной части кривой. Учитывая все проделанные замены переменных в процессе перехода к кривой, описываемой формулой (8), находим
Л = | [3(1 + 0'1/3 - 0'2/3) + Q'1/3p sin <р - а2/3р cos у] ,
а2 = 1 + 3(а1/3 + а2/3) + а1/3р sin < - а2/3 р cos <
с функцией р(<), определяемой (7). Эти формулы параметрически задают приближенно, с точностью до а2/3, кривую раздела фаз при малых значениях а, Заметим, что эти выражения для Л, а2 при < = 3п/4 удовлетворяют уравнению эллипса с той точностью, с которой они получены. Таким образом, проблемы, возникшие при построении пря-
а
асимптотического разложения вблизи а = 0.
Если перейти к исходным переменным па другой связной части в точке пересечения
диагонали (1, —1), то получим, что а2 = 1 — а1/3 — а2/3(1 + o(1)). Таким образом, на
этой связной части а2 < 1, и поэтому она не описывает кривую раздела фаз в нашей
а2
Литература
1. Фам Минь Туан, Вирчснко Ю.П. Анализ стохастической модели химической кинетики бинарной автокаталитической реакции /7 Belgorod State University Scientific Bulletin Mathematics & Physics. 2013. 12(155);31.' C.183-185.
2. Хоретхемке В., Лсфсвр P. Индуцированные шумом переходы: Теория и применение в физике, химии и биологии / Пер. с англ. /М.: Мир, 1987. 400 с.
3. Arnold L., Horsthemke W., Lcfcvcr R. White and coloured external noise and transition phenomena in nonlinear systems /7 Zs. Phvs. 1978. B29. P.367-373.
4. Фам Минь Туан, Вирчснко Ю.П. Анализ фазовой диаграммы в стохастической модели химической кинетики бинарной циклической реакции /7 Belgorod State University Scientific Bulletin Mathematics & Physics. 2013. 26 (169);33. C.57-63.
CRITICAL SURFACE INVESTIGATION OF CHEMICAL STOCHASTIC MODEL OF BINARY AUTOCATALYTIC REACTION. LARGE ASYMMETRIC CASE Pham Minh Tuan, Yu.P. Virchenko
Belgorod State University, Pobedy Str., 85, Belgorod, 308015, Russia, e-mail: [email protected],
Abstract. The stochastic model of binary chemical reaction is studied on the basis of chemical kinetics equations with stochastically perturbed parameters. The perturbation is described by generalized random process named "white noise". The stationary probability distribution density depended on relative concentration is investigated. The critical surface in the three thermodynamical parameter space of the mixture state is observed qualitatively. The noise induced phase transition in the stochastic system is under consideration. At this the reconstruction of unimodal distribution density into the bimodal one is occurcd at intersection of the critical surface. In the case of the large overcoming of the one component concentration over the other concentration component the perturbation theory is developed. It permits to calculate approximately critical surface at small values of parameter that characterizes such an asymmetry.
Key words: chemical reaction, stochastic model, distribution density, phase transition, bifurcation, critical surface, perturbation theory.
