Модель брюсселятора является одной из самых известных математических моделей синергетики. (Название связано с тем, что она была предложена в брюссельской научной школе.) Эта модель описывает распределение по пространству и изменение со временем реагентов сравнительно узкого класса химических реакций, однако при ее исследовании были выяснены свойства диссипативных структур во многих нелинейных системах.
Из школьного курса химии известен закон действующих масс. В реакции, где два вещества, Х и Y, реагируя, дают вещество Z (Х + Y → Z), скорость изменения вещества Z пропорциональна произведению концентраций веществ X и Y. Коэффициент пропорциональности – постоянная реакции k. Обозначая через X, Y, Z концентрации соответствующих веществ, можно записать
(1) |
В самом деле, для того чтобы реакция шла, молекулы вещества X должны сталкиваться с молекулами Y. Очевидно, вероятность этого пропорциональна числу молекул X в единице объема (т.е. концентрации
). Точно так же она должна быть пропорциональна концентрации . Коэффициент пропорциональности k зависит от размеров молекул, их скоростей и т.д. Все это и отражает формула (1). Если в реакции п молекул Х взаимодействуют с одной молекулой Y, то изменение концентрации вещества Z пропорционально XnY.Обратимся теперь к самой модели. Пусть в некотором химическом реакторе превращения идут по следующей схеме:
A ↔ X, B + X ↔ Y + D, 2X + Y ↔ 3X, X ↔ E.
Концентрации веществ А и В в реакторе поддерживаются постоянными, и некоторым образом удаляются вещества D и Е, т. е. система является открытой. Будем считать, что скорости обратных реакций (k–1, k–2, k–3, k–4) гораздо меньше скоростей прямых реакций (k1, k2, k3, k4). В этих предположениях, обозначая через
концентрацию вещества X, – вещества A и т.д., получим из закона действующих масс следующую систему уравнений:Концентрации реагентов
и могут быть различными в разных точках, поэтому в уравнение входят члены D1 xx, D2 xx, учитывающие их диффузию. После несложных замен переменных, эквивалентных переходу к другой системе единицмы придем к системе уравнений в частных производных, называемых моделью брюсселятора:
Xt = A – (B+1)X + X2Y + D1Xxx, Yt=BX – X2Y + D2Yxx. | (2) |
Вещества X и Y остаются в реакторе, поэтому потребуем выполнения следующих краевых условий:
Xx(0, t) = Xx(l, t) = 0, Yx(0, t) = Yx(l, t) = 0. | (3) |
Поведение решений
Посмотрим, есть ли у уравнения (3) какие-нибудь простые решения, например не меняющиеся со временем (их называют стационарными) и однородные по пространству. При этом все производные в (3) становятся нулевыми и мы имеем систему обычных алгебраических уравнений:
А – (B + 1)X + X2Y = 0, BX – X2Y = 0
Ее единственное решение – это Х = А, Y = B/А. В наших рассуждениях оно будет играть особую роль. Будем менять концентрацию вещества B и начальные распределения концентраций X(х, 0), Y(x, 0) и смотреть, как меняется поведение решения. В этом нам опять поможет ЭВМ.
Если концентрация вещества B невелика, то независимо от начальных данных через определенное время установятся концентрации Х(x, t) = A, Y(x, t) = B/A. Оказывается, такое замечательное решение (устойчивое стационарное, на которое независимо от начальных данных выходят изучаемые распределения параметров при небольших внешних воздействиях) есть у многих нелинейных систем. Оно получило название термодинамической ветви (в случае брюсселятора это решение Х = А, Y = B/A).
На первый взгляд кажется, что такая картина будет иметь место при любых В. Однако это не так. Если зафиксировать начальные концентрации Х(х, 0), Y(х, 0) и увеличивать значение B, то мы увидим, что начиная с некоторого критического значения B происходит выход на немонотонные стационарные распределения концентраций, например такие, как показаны на рис.1 и 2.
Рис. 1. Стационарные диссипативные структуры, возникающие в модели брюсселятора.
Параметры нелинейной среды: А = 2; B = 4,6; D1 = 1,6·10–3; D2 = 8,0·10–3
Рис. 2. Распределение концентрации X.
Два различных типа структур, возможных в одной и той же нелинейной среде при задании различных начальных данных. Параметры нелинейной среды: A = 2; B = 4,6; D1 = 1,6·10–3; D2 = 8,0·10–3.
Именно для таких стационарных неоднородных по пространству устойчивых решений, возникающих вне термодинамической ветви, И.Пригожиньм и было впервые введено понятие диссипативной структуры.
Прежде чем разбираться подробнее в свойствах таких решений, подчеркнем неожиданность полученного результата. Кажется очевидным, что в реакторе распределение реагирующих веществ по горизонтали (если сила тяжести направлена по вертикали) будет однородным по пространству. Модель брюсселятора показывает, что это не так: в среде могут возникать структуры, одни реагенты могут оказаться сосредоточены в одних частях реактора, другие – в других. Здесь встает целый круг вопросов:
как меняют структуры характерные времена реакций?
какая концентрация вещества является оптимальной?
И много других. Такие вопросы возникают при решения ряда задач химической технологии.
Вернемся к модели брюсселятора. Стационарное решение Х = А, Y = B/A удовлетворяет краевой задаче при любых B. Следовательно, при B > В0 появляется несколько стационарных решений. Как говорят математики, происходит ветвление решений, или бифуркация. Аппарат теории бифуркаций, интенсивно развиваемый в настоящее время, широко используется в синергетике.
Мы зафиксировали начальные концентрации и меняли В. Поступим по-другому: зафиксируем какое-нибудь значение В > В0 и будем менять профили начальных концентраций X(х, 0), Y(x, 0). При некоторых значениях B можно наблюдать интересный эффект: при одних начальных данных имеет место выход на один стационар (стационарное решение), при других – на другой. Два стационара, возможные при одних и тех же параметрах, показаны на рис.2. Причем выход на один и тот же стационар происходит с целого класса начальных концентраций, т.е. так же, как в модели тепловых структур здесь имеет место «забывание» деталей начальных данных. А что будет, если поставить систему в положение буриданова осла – задать при тех же значениях начальные условия, приводящие к однородному решению Х(х, 0) = А, Y(x, 0) = B/A, соответствующему термодинамической ветви?
Роль флюктуаций
Если решение Х = А, Y = В/А «поставлено» идеально точно, то оно меняться не будет. Однако реально расчеты на ЭВМ дают другую картину. Даже очень малые отклонения, которые, как правило, всегда имеют место, быстро нарастают, и далее происходит выход на один из неоднородных устойчивых стационаров. Такие отклонения, называемые флюктуациями, всегда есть в физических, химических и биологических системах. Расчеты на ЭВМ показывают, что вносимые флюктуации в отличие от равновесных процессов, изучаемых классической термодинамикой, определяют всю дальнейшую судьбу нелинейной системы. Термодинамическая ветвь здесь неустойчива.
Рис. 3. Неустойчивое состояние равновесия (точка O). Флюктуация выводит шарик из равновесия; в точке M и N – устойчивое состояние равновесия.
Этот процесс можно пояснить следующим примером. Представим себе маленький шарик в желобе, форма которого показана на рис.3. Если поставить его на вершину горба, в точку О, то в соответствии с законами механики он может оставаться на вершине (это тоже стационарное решение уравнений, описывающих движение шарика), но флюктуации выведут его из равновесия и он начнет двигаться. Постепенно из-за трения энергия шарика будет уменьшаться, и в конце концов он остановится на дне желоба в точке М или N. В какой именно точке он окажется, зависит от знака флюктуации, которая вывела шарик из равновесия. Роль точки О у нас играла термодинамическая ветвь, роль равновесных положений М и N – стационарные устойчивые решения, такие, как показаны на рис.2. Можно сказать, что причиной возникновения структур являются внутренние свойства системы, а поводом – вносимые флюктуации. Такое поведение характерно для многих нелинейных неравновесных систем.
Рис. 4. Возможный вид случайной функции F(t).
Флюктуации можно учесть, добавив в правую часть уравнения (2) случайные функции. Они могут отражать процессы, в детали которых на нашем уровне описания мы не вникаем. Отвлекаясь от их конкретного вида, приведем простейший пример случайной функции. Бросаем монету с интервалом времени Δt и считаем, что если в момент времени t выпадает орел, то F(t) = α, α << 1 δо момента t + Δt, если решка – F(t) = α; β момент времени t + Δt мы опять бросаем монету. Возможный вид функции, полученной таким образом, показан на рис.4. «Возможный» потому, что точно неизвестно, когда выпадает орел, а когда решка. Функция действительно случайная. И, бросая монету, читатели могут получить функцию нисколько не хуже нарисованной здесь.