CC BY
29
Раздел V. Математическое моделирование экосистем
А.С. Черепанцев
КРИТИЧЕСКОЕ СОСТОЯНИЕ В OFC-МОДЕЛИ ПРИ РАЗЛИЧНЫХ ДОПОЛНИТЕЛЬНЫХ УСЛОВИЯХ
Интерес к предложенной Олами, Федер, Кристиансен (OFC) [1] модели само-организованного критического состояния, возникающего в простейшей системе связанных ячеек, во многом определяется возможностью описания с ее помощью скольжения тектонических блоков при рассмотрении их как отдельных объектов, обладающих локальными напряжениями [2]. Данная модель носит неконсервативный характер, т.е. общая энергия системы уменьшается со временем, и для поддержания системы в динамически активном состоянии требуется наличие внешнего источника энергии. Поведение системы [3] определяется, наряду с законом передачи энергии между соседними ячейками, рядом дополнительных параметров и условий. В данной работе рассмотрены два таких условия - условия на границе рассматриваемой области и условия перераспределения энергии при передаче энергии соседним ячейкам. Последнее можно задать в двумерной модели с помощью задания различных регулярных сеток на плоскости.
Физическая модель. Моделирование организации сейсмического процесса в пространственно-временной области послужило одной из причин появления рассматриваемой OFC-модели. Механическая модель [2] (stick-slip model) представлена на рис. 1.
Она представляет собой двумерную квадратную решетку, состоящую из одинаковых блоков. Каждый из блоков связан с соседними четырьмя блоками упругой пружиной. Блоки расположены на неподвижной платформе и соединены упругой связью с движущейся выше Рис. 1. Модель прерывистого скольжения плитой. Сила, действующая на отдельный блок с координатой (i,j), определяется как
F,j = fj + f+w + f-i,j + fuj+i + fuj-i,
где fl - сила, управляющая движением и образованная упругой связью с верхней плитой, f.±1 j±1 - силы, образованные упругими связями блока (i,j) с соседними блоками. В случае отсутствия силы трения между основанием блока и неподвижной платформой движение блоков было бы непрерывно. В случае же наличия
трения движение носит скачкообразный характер. Если сила Ftj не превышает силы трения покоя, то блок остается неподвижным. Когда же сила F.} становится больше силы трения покоя, блок перескакивает в новое положение, такое же, как
суммарная сила р / = 0. Для записи уравнений движения определим основные параметры системы. Пусть Кх, К2 - коэффициенты упругости пружин по х и у направлениям соответственно. Пусть ¡¡, 12 - нулевая длина пружин при отсутствии растягивающих сил. Будем предполагать, что блоки под действием приложенных сил могут деформироваться, при этом среднее расстояние между блоками ак>1к , к = 1,2. Пусть (х. /, у. /) определяет смещение блока (', ]) из положения равновесной конфигурации. Предполагается, что пружины подчиняются линейному закону Гука упругой деформации и движение верхней плиты происходит вдоль оси Ох.
Тогда координаты силы, действующей на блок (, ]), определяются как
= -Кьх1,1 + + ./¿-и + /х,/+\ + -С-,
ру = /++и+/-и+/У+1+/¿-1. (1)
Для проекций сил взаимодействия блоков справедливо равенство
/1+1, / = I /1+\, /1 с°® +1, ] = /+1, ] |
С учетом того, что \/1+1,] | = К1 (1-А) = К1 (^Л/(а1 + х‘+1,/ - х.,/)2 + (у.+1,/ - у.,])2 -110 =
получаем: /х+х, = К1
а1 + х.+1 - х -
11 (а1 + х+1 - х.)
л/(а1 + х+1 - х)2 + (у<+1- У)2
Аналогично можно получить:
/-1,/ =- К1
/ х+1 = К 2
/х-1=-К 2
а1 + х,] - х.-и -
11 (а1 + х,] - х,--1,])
л/(а1 + хм - х-1, ■ У + (у, ]- У.-1, ■)2 12 (х, /+1 - х, / )____________________
д/к]+1 - х.,] )2 +(а2 + У,]+1 - У,] )2
12 (х,] - х,]-1 )
-1 )2
\2 + У,/ - У.,/- -
Подставляя полученные соотношения в (1), получим выражение для составляющей силы по координате х. Аналогично рассуждая для у-координаты силы справедливо
р У/ =- К 2
12 (
а 2 + У,/+1- У,/
(х,]+1 - х,])2 +(а2 + У,/+1 - У,/)2
12 (а2 + У, / - У, /-1 ) (х,/ - х,/-1)2 +(а2 + У, / - У, /-1)2
_____________11 (У+1,/ - У, /)______________
(а1 + х,-+и - х,/)2 + (У.+1,/ - У,/)2 л/(а1 + х,/ - х-1,/)2 + (У,/ - У-1,/)
--К1 -2У,/ -У+1,/ -У-1,/ +
V
11 (У,/ - У-1,/)
При соскальзывании блока (і, 1), новое значение F¡j ® 0. Тогда изменение компонент сил, действующих на соседние блоки, определяется новыми значениями (хі: 1, у¡, 1). Так, для блока (і +1, 1) имеем
Р +1,1 - Р +і,; =-К (- (У,; - X, ]) -
°1 + хі*1,1 - У. 1
)
¡1 (а1
+ У+1,1 - Хг,1
)
)2 + (Уг +1,1 - У, 1 )2 л/(а1^Хі+1^/™Хі,/)^^^^/™-Уі,/)
(а1 + У+1,1 - Уг, 1,
где (у ., У; .) - новые координаты блока (і, 1).
Соответствующее изменение силы, действующей на блок (і, 1), равно
Ут т т І \ ( ¡1(а1 + X+1 ,■ - У,-)
Ри - рх = 0 - ру = К (у. 1 - X, 1) + к ^у -- '+ 1и г+1,1 ¡-11
¡1(а1 + У+1,1 - У, і)
(2)
(а1 + У+1,1 - У, 1} + (у.+1,і - .Уі, 1} ¡1 (а1 + У, 1 - У-1,1) ,
(а1 + У+1,1 - У, 1У + (уі+1,і - Уі, 1У д/(а1 + У,і - У-1,і ^ +(уі,і - Уі-1,1У
¡1 (а1 + У, 1 - У-1,1)
-к
2(У, і - У, і) +
¡2 (У,і+1 - У,і )
(а1 + У,і - У-1,і )2 +(у,і - Уі-1,і У ^ [ , 1 ,і л/(У,і+1 - У,і У +(а2 + Уі,і+1 - У,і )2
____________¡2 (У, 1+1 - У,і )___________________¡2 (У, 1 - У, 1 -1)________________________+
(У,і+1 - У,і У +(а2 + Уі,і+1 - Уі,і У т/(У,і - У, 1 -1У +(а2 + Уі,і - Уі,і-1У
____________¡2 (У,і - У,і-1У_______________
(У,і - У, 1 -1У + (а2 + Уі,і - Уі,і-1)2
(3)
Линеаризуем данные выражения с учетом ак >> Ах. 1; ак >> Ау. 1; к=1,2. Тогда (2), (3) можно представить в виде
Р+1, і - Р+1, і » к 1(У, і - у , і);
- К (У,1 - У,1)+ 2К1(У,1 - У,1) + 2К 2(У,1 - У,1) - К2-^(х/ - У,!).
Определим константы внутренней деформации 5к = (ак - ¡к)/ак , коэффициенты анизотропии системы а = К2/К1 , к = К1/К1. Тогда новое значение действующей х-компоненты силы можно представить через ее предшествующие значения
рх = рх +а рх Л+1,1 Г 1+1,1 + “И., 1 .
Аналогично можно получить соотношения и для изменения остальных компонент силы:
Р-1,1 = р,-и + а , ^±1 = р‘х,1 ±1 +52^«1 , где а1 = - 1
РУ = РУ +-ґі±1,1 Гі±1,і ^ а - і,і
Р" , РУ±1 = РУ±1 +«2, где «2 =
2(1 + 5 2ст) + к 1
2(1 + 51ст)
Полученная модель переходит в модель ОБО в случае, если ¥У = 0 и 5 2ст = 1.
Первое условие означает, что можно пренебречь движением по у-координате. Второе условие определяет существование зависимости меду деформацией по оси Оу и анизотропией упругих свойств системы. Параметр а1, определяющий поведение модели ОБО, в данном случае определяется физическими параметрами анизотропии и внутренних напряжений.
Граничные условия в OFC-модели. Принято выделять три основных типа граничных условий в рассматриваемой модели.
Свободные условия на границе, в случае рассмотрения квадратной решетки, определяются как условия, при которых энергия сброса у граничных элементов распределяется между соседними элементами - двумя для угловых элементов в
4
количестве 2аЕтах и тремя - для граничных элементов в количестве — аЕтах . При
этом общая энергия, переданная соседним ячейкам при сбросе, остается постоянной величиной, равной 4аЕтах < Етах для любой области, будь то внутренняя или граничная области. Для двух противоположных границ условия могут быть записаны как
где] = 0,...Ь - 1; к= 2 при] = 0, Ь - 1; к = 4/3 при] = 1,..., Ь - 2.
Для двух остальных противоположных границ условия формулируются аналогично.
Открытые условия на границе определяются тем, что все элементы рассматриваемой системы имеют четыре соседних элемента, т. е. рассматривается система элементов, являющаяся частью большей системы с теми же параметрами. Вне зависимости от положения, соседние элементы при сбросе элемента (,, ]) получают добавку величиной аЕтах. Физическая реализуемость такой системы определяется условиями на границе окружающей системы
В данном случае общая потеря энергии системы (ее диссипативность) определяется не только потерей энергии при сбросе (а<1/4), но и оттоком энергии через границы.
Периодические условия на границе формулируются для систем, в которых отсутствует характерный размер системы Ь. В модели с такими условиями решетка сворачивается в тор и условия на границе решетки размером ЬхЬ выглядят следующим образом:
Е1,з ® Е1,з + каЕ0,з
Е1,з ® Е1,з + аЕ0,з
Е,
Е
Ь-1
Ь-1,3
> Е,
Етах ® Е,,0 ® Е,,0 + аЕ,,Ь-1 .
'тах
Как показали модельные расчеты, при значениях а<0,25 в такой системе не происходит развития самоорганизованного состояния и соответственно степенного характера функции плотности распределения размеров сброса. Сброс отдельной ячейки практически никогда не сопровождается нарушением устойчивости сосед -них ячеек и не происходит лавинообразного нарастания пространственного масштаба сброса. В этом случае в модели повторяется периодически одна и та же последовательность сбросов.
В качестве критерия достижения критического состояния рассмотрим характер распределения сбросов по размерам в течение заданного числа итераций. Устойчивое критическое состояние определяется существованием степенной зависимости распределения с постоянным значением показателя степени / в диапазоне масштабов \,2,...Smax , где Smax<L2. Под масштабом сброса понимается количество связанных ячеек, участвующих в сбросе в течение отдельной итерации.
Результаты моделирования представлены на рис.2, 3.
Анализ полученных зависимостей позволяет сделать ряд выводов.
1. Выбор граничных условий является существенным параметром модели, определяющим вид полученного предельного состояния. Показатель степени / распределения сбросов по масштабам значимо отличен для моделей со свободными и открытыми границами. При периодических граничных условиях критическое состояние системы не достигается.
2. Сходимость к критическому предельному состоянию при свободных граничных условиях выше, чем при открытых границах.
3. Предельное распределение сбросов при открытых граничных условиях в двойном логарифмическом масштабе является идеально линейным, в отличие от зависимости в случае свободных границ. В последнем случае распределение носит достаточно изрезанный характер.
С целью выявления причин такой разницы поведения предельного распределения рассмотрим временное изменение такого параметра системы, находящейся в критическом состоянии, как средняя энергия <Е> по ансамблю частиц.
На рис. 4 представлены вариации <Е> .
Рис. 2. Выборочная плотность распределения f сбросов по масштабу S при заданных граничных условиях: а - свободные, б - открытые. Параметры модели: Ь=64; а=0,15; Етах=30; ЛЕ=0,1. Кривая 1- число итераций п = 250; 2 - п = 103; п = 3103; п = 104; п = 3104; п = 105; п = 3105
Рис. 3. Предельное распределение сбросов при различных граничных условиях: а - выборочная плотность распределения сбросов для ОЕС-модели L = 64, а = 0,2, Етах = 30; ЛЕ = 0,1 : 1 - свободные границы при п = 107; 2 - открытые границы при п = 107; 3 - периодические границы при п = 106; б, в, г - распределение энергии по ячейкам после указанного числа итераций соответственно при свободных, открытых, периодических граничных условиях
б
г
<Е: 18 4 16 14
200
400
600
800
Рис. 4. Вариации среднего значения энергии системы и оценки корреляционной размерности системы со свободными граничными условиями (а), (в) и с открытыми граничными условиями - (б), (г)
Как видно из рис. 4, а, в случае свободных граничных условий устойчивое кри -тическое состояние определяет близкое к стационарному поведение системы в
г
целом, при котором энергия распределена по дискретному набору масштабов. Поэтому система характеризуется изрезанным характером кривой распределения по масштабам (рис. 3,а). Более сложное поведение во времени показывает система при заданных открытых граничных условиях (рис. 4,б). В этом случае режим является апериодическим, определяя возможность сбросов любых масштабов и соответственно гладкость прямой распределения событий (рис. 3,а). Сложный характер поведения системы в этом случае можно показать и путем оценки размерности динамической системы по виду изменения корреляционного интеграла (рис. 4,в,г). На различный характер поведения системы указывает и анализ распределения энергии в предельном случае (рис. 3,б,в). В случае свободной границы система демонстрирует однородный характер поведения на всей решетке, т.е. пространственные структуры больших масштабов, имеющие одинаковое значение энергии, могут включать и граничные точки, что, по всей видимости, определяется свойством постоянной величины диссипации энергии в любой ячейке решетки. При этом в случае открытых граничных условий поведение системы на границе и во внутренних точках различно. На границе поведение системы остается случайным и пространственные структуры данную область не затрагивают. Это, в частности, определяет то, что возможный максимальный масштаб сброса в системе с открытыми границами меньше максимального масштаба в сисеме со свободными границами.
Регулярные решетки в OFC-модели. Рассмотрим, насколько зависят параметры степенных распределений от выбора типа регулярной решетки, плотно покрывающей ограниченную область поверхности. В качестве таких решеток на плоскости могут быть рассмотрены треугольная решетка, состоящая из равносторонних треугольников, квадратная решетка и шестиугольную решетка, состоящая из правильных шестиугольников и также плотно покрывающая поверхность.
Для треугольной сетки правила эволюции системы определяются как
Еі,2 і > Етах —
Еі Лі — 0
Е,,2і±1 — Еі,2і±1 + аЕі,2і Еі-1,2і — Еі-1,2і + аЕі,2і
Е > Е —
і,2 і+1 _ тах
Еі,2і+1 — 0 Еі,2і+1±1 — Еі,2і+1±1 + аЕі,2і+1 Еі+1,2і+1 — Еі+1,2і+1 + аЕі,2і+1
1=0.1 - 1; ]=0...(Ь - \)/2.
Модель становится консервативной при а = —. Для шестиугольной сетки правила эволюции определяются как
Еи — 0
Е > Е —
^і,і — тах ^
Е —— Е + аЕ
Еі ,і ±1 — Еі, і ±1 + аЕі і Еі+1,і+1 — Еі+1, і+1 + аЕі і Еі-1,і+1 — Еі-1, і+1 + аЕі і
,і,]=0...Ь - 1.
Рис.4. Зависимость числа итераций п, необходимых для появления сброса размером 5. а - зависимости для: 1-треугольной сетки; 2 - квадратной сетки; 3 - шестиугольной сетки. б -начальный участок зависимости в полулогарифмическом масштабе
Модель становится консерватив-1
ной при а = — .
6
В качестве характеристики перехода системы к критическому состоянию рассмотрим число итераций, необходимых для возникновения сброса заданного масштаба. Соответствующие зависимости п(я) для различных сеток представлены на рис. 5,а. Можно выделить два участка зависимости - начальный переходной участок, с экспоненциальной зависимостью роста (рис. 5,б) в диапазоне малых масштабов 5 = 1 - 20 и степенной рост зависимости для 5>20. Экспоненциальная зависимость п ~ вр характерна для систем со случайным несвязанным поведением отдельных ячеек сетки [4]. Степенная же зависимость
С ^
п ~ п0 5 свойственна системе при
наличии взаимодействия отдельных элементов. Расчетные значения параметра диссипации а:
ак
= 3 =3
'“квад — ^треугольная — ^шестиугольная
4 2
= 0,15 определяют одинаковое значение диссипации энергии при отдельном сбросе: 0,4£тах. Этим определяется общий вид кривых и соответственно скорость сходимости к предельному состоянию. В случае же задания одинаковых значений параметра а для всех решеток, наивысшей скоростью будет обладать шестиугольная решетка. Тип решетки при фиксированном значении диссипации влияет также в незначительной степени
на степенной показатель: Хтреугольная = 1,5 ± °1; Хквадратная = 1,6 ± °1;
Хшестиугольная 1,8 ± °1
Таким образом, основным параметром, определяющим скорость сходимости к предельному состоянию, является не тип решетки, а величина параметра диссипации энергии в модели ОБС.
БИБЛИОГРАФИЧЕСКИМ СПИСОК
1. Olami Z., Feder H., Christensen K. Self-organized criticality in a continuous, nonconservative cellular automaton modeling earthquakes, Phys.Rev. Lett. 68, P. 1244 -1247. 1992.
2. Burridge R., Knopoff L. Model and Theoretical Seismicity, Bull. Seism Soc. Am. 57, P.341 - 371, 1967.
3. Черепанцев А.С. Связь пространственно-временных параметров в OFC-модели тектонических процессов // Изв. ТРТУ. Тематический выпуск «Экология -2006». 2006. С. 5.
4. J.J.Binney, N.J.Dowrick, A.J.Fisher, M.E.Newman. The Theory of Critical Phenomena. Oxford University Press. 1992.
СВЯЗЬ ПРОСТРАНСТВЕННО ВРЕМЕННЫХ ПАРАМЕТРОВ В OFC-МОДЕЛИ ТЕКТОНИЧЕСКИХ ПРОЦЕССОВ
Фрактальные структуры встречаются в различных физических системах: от моделей образования снежных хлопьев до распределения галактик. Турбулентность является примером системы, проявляющей фрактальные свойства. Диссипация энергии в такой системе происходит не пространственно инвариантно, а каскадом на определенных пространственных масштабах.
Другой известной системой, демонстрирующей фрактальные свойства, является сейсмичность. Система тектонических плит имеет масштабно-инвариантные закономерности: измерение распределения плит по размерам. Распределение энергии землетрясений по частоте повторений, описываемое законом Гуттенберга -Рихтера, также дается степенной функцией. Наряду с этим и геометрическое распределение эпицентров (проекций центров землетрясений на поверхность) и их временная последовательность показывают фрактальность структуры.
Фактически сейсмичность имеет много общего с турбулентностью. Обе являются открытыми динамическими системами с большим количеством независимых элементов, взаимодействующих между собой. Обе системы управляются внешними процессами, с постоянным притоком внешней энергии и диссипацией через каскад масштабных размеров. Диссипация энергии в обеих системах показывает дискретный характер во времени и пространственно-временную фрактальную организацию.
Модель скольжения тектонических блоков при рассмотрении их как отдельных объектов, обладающих локальными напряжениями, предложили Олами, Фе-дер, Кристиансен (ОБС) [1]. Данная модель носит неконсервативный характер, т.е. общая энергия системы уменьшается со временем, и для поддержания системы в динамически активном состоянии требуется наличие внешнего источника энергии.
Пусть дана кубическая решетка размерностью ё и размером Ьа. Поставим в соответствие каждой г-й ячейке некоторый динамический параметр Ег. В простейшем случае под Ег будем понимать внутреннюю энергию, запасенную в г-й ячейке. Предположим, что в единицу времени все ячейки получают одну и ту же добавочную величину приращения энергии:
Такое изменение во времени энергии любой отдельной ячейки, в случае отсутствия влияния соседних ячеек происходит до тех пор, пока Е1 < Етах, где Етах -некоторое пороговое максимальное значение упругой энергии, при превышении которого ячейка сбрасывает накопленную энергию, часть из которой передается соседним ячейкам:
А.С. Черепанцев, С.Ф. Черепанцев
Et ® E{ +DE, i=1,... ,N.
(1)
(2)
где индекс кк определяет соседние ячейки.
