Арсеньев С.А. ©
Доктор физико-математических наук Институт физики Земли им. О.Ю.Шмидта. Российская Академия Наук
НЕЛИНЕЙНАЯ ТЕОРИЯ ФРИКЦИОННЫХ ЗЕМЛЕТРЯСЕНИЙ
Аннотация
Линейная теория возникновения землетрясений и тектонического тремора в активных разломах, построенная ранее [1,2], обобщается на случай учета нелинейного трения. Получены формулы для расчетов и исследований.
Ключевые слова: землетрясения, тектонический тремор, активные разломы. Keywords: earthquakes, tectonic tremor, active faults.
ПОСТАНОВКА ЗАДАЧИ
Рассмотрим фрикционную теорию землетрясений, предложенную в работах [1,2].
Согласно этой теории колебания поверхности Земли {то есть собственно землетрясения}
возбуждаются при движении берегов активного разлома, причем генератором сейсмической
радиации является осциллирующий блок земной коры, источником энергетического питания
- движущий берег разлома, а поступление {или убыль} энергии регулируется нелинейным
элементом - силой трения, которая зависит от относительной скорости движения берегов
разлома и скорости колебаний сейсмического блока внутри разлома.
На рис.1 показана типичная зона субдукции, где океанская плита погружается под
континент со скоростью u, [1]. Сейсмическая радиация может излучаться блоками
континентальной коры, обозначаемыми цифрами 5,6,7..., поскольку они лежат на
движущейся океанской плите и могут получать энергию от этого движения. Введем массу
блока m, а смещение блока обозначим буквой x. Будем также обозначать буквами с
индексом t производные по времени, например скорость смещения dx/dt = xt , ускорение 2 2
смещения d x/dt = xtt и т.д. Ось х направлена вдоль океанской плиты 1. Она действует на блок благодаря силе трения F{V}, которая зависит от относительной скорости V = u - xt . Кроме того, блок упруго взаимодействует с соседними блоками и континентальной верхней плитой. Обозначая результирующий коэффициент жесткости буквой k, уравнение движения блока запишем в виде
mxtt + h xt + k x = F{V} . {1}
© Арсеньев С.А., 2016 г.
1 - океанская плита, 2 - аккреционная призма, 3 - континентальная плита, 4 - океан, 5,6,7 - сейсмогенные блоки литосферы, излучающие сейсмические волны, 8 - верхняя мантия. I - оффшорная береговая зона, II - зона запирания, III- переходная зона, IV-
зона тремора, V - пластическая зона
Сила трения F{V} в зоне субдукции в настоящее время неизвестна. Тем не менее, мы можем попытаться определить ее, используя механические
аналоги и лабораторные эксперименты с трибометрами {приборами для изучения трения}.
СИЛА ТРЕНИЯ
Она должна отражать скольжение океанской плиты под континент, вследствие дегидратации. Трение является скользящим, потому что вода выталкивается на поверхность погружающейся плиты, образуя силикатные растворы, играющие роль смазки. Возможно, также частичное расплавление верхней части тяжелой океанской плиты, погружающейся в мантию. В этом случае смазкой является расплавленная магма в зоне контакта океанской и континентальной плит. Таким образом, трение не может быть описано законом сухого трения Кулона-Эйлера, согласно которому функция F не зависит от площади контакта трущихся поверхностей и их относительной скорости V, [3].
Подходящей механической моделью для данного случая является скрипичная струна {аналог континента} по которой движется смазанный канифолью смычок, играющий роль океанской плиты на рис.1. Это движение возбуждает колебание струны, которое мы слышим как звук, усиливаемый резонатором - деревянным корпусом скрипки. Функция трения F{V} в данном случае хорошо известна [4]. Она имеет вид кривой в форме ямы. При V = 0 {смычок неподвижен} она максимальная Fmax .Затем, когда смычок движется с увеличивающейся скоростью V, функция спадает до минимума Fmin, а потом снова возрастает. Колебания возникают на падающем участке функции F{V}. При выбранной скорости движения смычка сила трения уравновешивается натяжением струны, но это состояние неустойчиво. Случайные колебания струны в направлении движения смычка уменьшают относительную скорость V и вызывают увеличение силы трения F. Это нарушает состояние равновесия и приводит к возникновению звуковых колебаний струны, что следует и из энергетического баланса [4].
В скрипичной струне колеблющаяся масса является распределенной, а амплитуда
звуковых колебаний малой величиной. Кайдановский Н.Л. [5] изучал в лабораторных
экспериментах возникновение и усиление колебаний отдельного блока заданной массы m,
который ставился на движущуюся, смазанную маслом поверхность стола трибометра. Блок
прижимался к этой поверхности пружиной, закрепленной на кронштейне, нависшем над
столом. Нарастающие колебания блока, достигающие больших амплитуд, измерялись
шлейфовым осциллографом. Была измерена и сила трения, которая, которая также как и в
случае скрипичной струны имела вид ямы, изображенной на рис.2. Здесь показана функция
F{V} = Fmax + Dx{s+V}-2 - Bx{s + V}-1 , {2}
аппроксимирующая графики, полученные Кайдановским Н.Л. [5]. Значения Fmax = 18000
2 2 1 1 дин, D = 47600 дин см с и B = 47000 дин см с , s = 1 см с соответствуют его
экспериментальным данным [5]. V - относительная скорость: V = u - xt , где u - скорость
движения стола трибометра.
Исследуем характерные особенности функции трения F{V}, описываемые формулой
{2}. При V = 0, F{0} = Fmax , так как
D/s2 = B/s . {3}
Здесь s = 1 см/с - масштаб скорости. Введем также масштаб массы m, масштаб частоты ю0 = 1/2
{k/m} и безразмерную скорость W = V/s. Тогда масштаб силы F* = m s ю0 и формулу {2} можно переписать в виде
ф = Fmax + X x { 1 + W}-2 - G x { 1 + W }-1 , {4}
где Ф = F/ F* , Fmax = Fmax / F* , x = D/{s2 F*}, G = B/{s F* }. Из {3} следует, что X = G, а формула {4} принимает вид
Ф = Fmax + X X {1 + W}-2 - X X {1 + W}-1 , {5}
Рис. 2. Зависимость {2} силы трения F от относительной скорости V
Из рис.1 видно, что функция F имеет минимум Fmax . Ему соответствует dF/dV = 0. Из {5} следует, что
3 2
dF/dW = - 2 X X {1 + W} + X X {1 + W} , то есть минимум достигается в точке W = 1, V = s = 1 см/c. С другой стороны, из {5} следует, что при W = 1
F{1} = Fmin = Fmax - Е/4. Отсюда находим
X =4x{Fmax - Fmin }. {6}
Существенно, что Е > 0, так как трение покоя является максимальным и Fmax > Fmin .
Функция {5}, определяющая трение, замыкает задачу. Подставляя ее в уравнение {1}, получим
^ + + y = Ф
дт2 дт max (1 + w)2 (1 + W)
где т = Шо t, 8 = h / m Шо , y = x Шо /s. Уравнение {7} c существенно нелинейной правой частью в квадратурах не интегрируется. Необходимо должным образом аппроксимировать функцию трения {5}, чтобы получить аналитическое решение.
Из данных наблюдений за спредингом океанического дна известно, что средние скорости движения океанской плиты в зоне субдукции u очень малы [6]. Максимальные скорости u достигают в некоторых регионах 6 см/год, то есть приблизительно 1 нанометр в секунду ! Геологические данные, однако, сильно осредненны по очень большим промежуткам времени. Во время землетрясений скорости подвижек плиты могут достигать заметных величин, вполне сравнимых со скоростью колебаний сейсмогенных блоков внутри разлома. Поэтому разность двух близких величин u и xt , то есть относительная скорость V, является малой величиной. При характерной скорости s = 1 см/с, является малой величиной и безразмерная относительная скорость W << 1. Таим образом, естественно разложить функцию трения в ряд Тейлора по малому параметру W
{1 + W}m = 1 + mW + ^^ W2 +... 1! 2!
m{m-1}...{m-n +1} „n
... + —---—1-- W + ...
n!
и ограничиться несколькими, например, тремя или пятью членами в этом разложении. Имеем
{ 1 + W}-2 = 1 - 2W + 3 W2 - 4 W3 + 5 W4 - 6 W5 + ... { 1 + W}-1 = 1 - W + W2 - W3 + W4 - W5 + ...
Подставляя эти формулы в выражение {5}, найдем F = Fmax - X X {W - 2 W2 + 3 W3 - 4 W4 + 5 W5 } . {8}
На рис.3 показаны три кривые, аппроксимирующие функцию трения F{V}, изображенную на рис.2. При s =1 см/c численные значения V и W совпадают. Толстой сплошной линией обозначена исходная кривая F{V} на участке V от 0 до 0,5 см/c, соответствующая рис.2. Цифрой 3 обозначено третье приближение с сохранением членов до W3 включительно в разложении {8}. Цифрой 5 показано пятое приближение, то есть кривая соответствующая полной формуле {8}. Как
19000 18000 17000 16000 15000 14000 13000 12000 11000 10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 0
0
0,1
0,2
0,3
0,4
0,5
Рис.3. Аппроксимация функции трения. Цифрой 3 - обозначено кубическое приближение, цифрой 5 - с пятой степенью. Жирная линия - точная кривая
видим, кривая 5 гораздо лучше описывает функцию трения, чем кривая 3. Тем не менее, при малой разности V во многих случаях мы можем использовать, третье более простое приближение. Поэтому представлять интерес исследование обоих случаев, и 3 и 5.
ОСНОВНОЕ УРАВНЕНИЕ Подставим разложение {8} в уравнение {7}. Получим
d2 y ~dy _
-+ У = Ф +
dr2 dz max
{9}
+ Sx{-W + 2 W2 -3W3 + 4W4 - 5W5 }
Безразмерная скорость W = w - yr зависит от безразмерной скорости плиты w = u/s и скорости колебаний блока dy/dr = yr . Чтобы раскрыть эту зависимость в уравнение {9} нужно подставить разность W = w - yr и возвести ее в степень 2,3,4,5, используя бином Ньютона. В результате получим уравнение
а
d y dr
+ + у. 2 dr
Ф
max
+ x к +
{10}
— ri d , 2 3 4 5 ->
+ XX[{ k — }Уг+ K2 Уг + k Уг + k Уг - k Уг}
где введены постоянные величины Kq
0 - w + 2w2 - 3 w3 + 4 w4 - 5 w5
K1 = 1 - 4w + 9w2 - 16w3 + 25 w4
K2 = 2 - 9w + 24w2 - 50 w3 2
k3 = 3 - 16w + 50 w ,
k4 = 4 - 25 w, k5 = 5 ,
Наконец, заменяя y на новую переменную z каноническому виду
3
y - Фп
{11} {12} {13} {14} {15}
S kq , приводим уравнение {10} к
d2 z
+ z =
{16}
dt
= хх[{ к--} 2 + к ^2 + к + к ^4 + к 251
1 1 ^ ' т 2 т 3 т 4 т 5 т 1
Существенно, что все коэффициенты к^ к2, к3, К4, к5 в уравнении {16} положительны. На рис.4 показаны зависимости коэффициентов к1 и к3 от безразмерной скорости плиты w. Как видим, коэффициент к1 не превышает единицы, а коэффициент Кз имеет небольшой минимум, но остается положительным.
0
0,2
0,4
w
i
0,6
Рис.4. Зависимость коэффициентов k-1 и k-3 уравнения {16} от безразмерной скорости плиты w
РЕШЕНИЕ ЗАДАЧИ
Величина х = 4{Fmax - Fmjn}/msroo является малой по сравнению с единицей, так как мы изучаем колебания блоков большой массы при сравнительно слабом трении скольжения
m >> 4 { Fmax-1 . {17}
s щ
Поэтому для решения уравнения {16} можно применить метод медленно меняющихся амплитуд {ММА} специально разработанный для исследования нелинейных колебаний и волн [7 - 10].
Согласно методу ММА решение уравнения {16} ищем в виде z = a{x} cos [х + 0{т}], {18}
где a{x} и 0{х} - неизвестные амплитуды и фаза колебаний, которые изменяются со временем. Далее, производная
zx = - a sin [х + 0{х}] {19}
вычисляется в предположении независимости a и 0 от времени. На самом деле они от времени зависят, поэтому необходимо требовать выполнения дополнительного условия ax cos {х + 0} - a 0х sin {х + 0} = 0 , {20}
Условие {20} нужно использовать совместно с основным уравнением {16}, которое с учетом формул {18}, {19}, {20} можно записать в виде
ax sin {х + 0} + a cos {х + 0} 0х = - х f К }, {21}
где
s . {22} f К} = {k -- К+ k х + k ^ + k К + k z5
Систему двух уравнений {21} и {20} нужно разделить на два отдельных уравнения для амплитуды a и фазы 0. Имеем
a = - X sin {х + 0} f К} , {23}
a 0х = - х cos {х + 0} f К} , {24}
В уравнениях {23}, {24}, х << 1, поэтому амплитуда и фаза являются медленными функциями времени х. Это означает, что возможно усреднение по времени правых частей {23}, {24} за период 2я, в течение которого амплитуда и фаза не меняются. Имеем
1 2p {25}
ax = -J хf{T}sint1 dV 2p 0
1 2p {26} a в =--í X f {z } cost dr ,
t 2ж t 1 1'
2Ж 0
где x1 = 0 + x. Из {19} следует, что zx = - a sin x1. Поэтому {22} имеет вид
S {27}
f {zj={k —_ }{-a} smj+к, a sin2 t1 - k a3 sin3 t1 +
+k asirf j-k a^sin51 Подставим {27} в уравнение {25}. Получим
1 2л S i
a=--f X{k— }{-a}sln2j dj-
2P 0 1 X 11 {28}
1 2p 9 ,
- — fXk.a2sin3T dr + 2p 2 1 1 2p 0
1 2p „3 • 4 ^ 1 2P 4 ■ 5 , + — f xka sin t dj--f xka sin t dj +
2p 0 3 1 1 2p 0 4 11 1 2P 5 6
+ — fx k a5 sin 6 t dj 2p 5 1 1 2p 0
Интегралы, входящие в правую часть уравнения {28} являются табличными [11] и легко вычисляются:
2p 2 3
J sin t dt^ =p; J sin = 0;
0 0
2p 4 3P 2P 5
f sin t dt =—; f sin t dt = 0;
J 1 1 4 J 11
0 4 0
2P • 6t 5p
f sin t «t =—; 0 8
В результате, уравнение {28} для амплитуды принимает вид
da , i 3 , 5 , {28}
— = са + Яа3 + aa5
dt
где
х, S> 2 3 X 5 X {29}
c=-{K-x}; 1^k3; a=—k
2 1 х 8 3 16
Аналогичные вычисления делаются и для уравнения {26}, определяющего фазу колебаний. Подставляя {27} в {26}, получим а 0т = 0, так как
2р / • т+1 х 2П
| {81пХ}т ооьХЛХ= —-Х = 0
0 \ т+1 0 при любом целом т Ф 1, [11]. Таким образом, фаза колебаний не меняется со временем и может быть произвольной.
Уравнение для амплитуды {28} решается методом разделения переменных
-г^*—V =г
2у с + 1у+ау2
и мы ввели новую неизвестную величину у = а . Отсюда, интегрируя по времени от 0 до т, найдем
iln
4 с
( 2 ^
I y 2IK,
c+Xy+ay
Здесь
1 S • i3°) —S = t 4c
{31}
S=y
y0 ax + Xl + c l - переменная интегрирования, а
c+ ly0+ay02 {32}
K1 =-2-
Уо
- постоянная интегрирования, y0 = {ao} и ao - амплитуда в начальный момент времени т = 0. Интеграл {31} вычисляется отдельно. Если 4ас > X2 , то
2 2ay +1 ^ . {33}
S= .- acrtg ,— + K2
д/4ас -X2 д/4ac -X2 Если же наоборот 4ас < X2 , то
S=—Ц М^КГ {34}
a (p - qj У-q 2
где p и q - корни уравнения a l + X l + c = 0. И наконец, в особом случае, когда 4ас = X2
S=--+К . {35}
2 ay + 1 4
В уравнениях {33}, {34},{35} константы K2, K3, K4 соответствуют начальным условиям.
АНАЛИЗ РЕШЕНИЯ Линейное приближение, когда X = 0, а = 0 рассматривалось в работах [1,2]. Уравнение для амплитуды {28} в этом случае a = da/dt = =ca легко решается a = a0 exp {ст}, где a0 -начальная, затравочная амплитуда и с - инкремент {если с>0} или декремент {если с<0} колебаний. Величина с зависит от разности двух чисел: Ti = Ski/2 {от слова Tilt - наклон, ибо это число определяет крутизну спадания характеристики трения} и Da = 8/2 {от слова Damper - демпфер, ибо оно определяется величиной h, обеспечивающей затухание колебаний за счет упругого взаимодействия блока с соседними блоками и верхней континентальной плитой - рис. 1}
с = Ti - Da . {36}
Величина c - является параметром, управляющим поведением данной динамической системы. Согласно теории [12], катастрофа - это скачкообразное изменение, возникающее в виде внезапного отклика системы на плавное изменение внешних условий, например в виде плавного изменения управляющего параметра. И действительно, в линейной теории, если Ti > Da, то с > 0 и возникающая неустойчивость типа "отрицательного трения" приводит к
землетрясению - экспоненциальному росту амплитуд колебаний блока в активном разломе [1,2].
Отметим, что число с является аналогом разности между текущим числом Рейнольдса Re (=Ti) и критическим числом Recr (=Da) в теории гидродинамической неустойчивости. Переход к турбулентности возникает, когда Re = Recr и c = 0.
Наоборот, если Ti < Da, c < 0, то амплитуда колебаний затухает по экспоненте,
причем этот процесс может идти и при отсутствии трения о плиту Fmin = Fmax , X = 0, Ti = 0.
Уравнение {1} в этом случае Ti = 0 имеет вид m xtt + h xt + k x = 0 и легко решается [1,2]: x
2 1/2
= x0 exp { - h t/2m} cos {w t + 0}, где w = [ {k/m} - {h/2m} ] . Как видим частота колебаний зависит от демпфирования h, причем линейная теория [1,2] позволяет учесть также и влияние трения о плиту
-[(Vif-^i2]1'2 {37)
m 2 m 2 ms Именно эта частота должна входить в общее решение {18}
z = a cos {w t + 0} . {38}
Смена режимов, то есть порог неустойчивости с = 0 имеет место при условии Ti = Da. В этом случае a = a0 и колебания принимают вид тектонического тремора [1,2]
x = x0 cos {W0 t + 0}, {39}
который действительно наблюдается в природе в зонах субдукции [1,2]. Частота тремора в формуле {39} совпадает с собственной частотой колебаний сейсмогенного блока w0 = {k/m}1/2, так как c = 0.
Линейная неустойчивость является, вообще говоря, слабой, так как инкремент с согласно формуле {36} определяется разностью двух близких по порядку величин Ti и Da и не может быть очень большим. Начальное, затравочное возмущение x0 здесь является любой случайной инфинитезимальной величиной, и конечный результат от нее фактически не зависит. Амплитуда возмущений растет экспоненциально при всех временах, что не соответствует наблюдениям. Слабая линейная надкритичность вблизи порога неустойчивости реально наблюдается только в моменты времени близкие к начальному [2]. Коэффициент к1 определяющий число Ti, как видно из рис.4 имеет порядок 1. Постоянные коэффициенты A и a согласно формулам {29} содержат коэффициенты к3 и к5 . Из формул {14} и {15} следует, что эти коэффициенты имеют порядок 3 и 5 соответственно {рис. 4}. Это означает, что слабая нелинейная неустойчивость с коэффициентом с в амплитудном уравнении {28} может конкурировать с сильной нелинейностью с коэффициентами A и a только при малых амплитудах в моменты времени близкие к началу землетрясения. Затем начинают доминировать нелинейные члены, которыми в уравнении {28} пренебрегать нельзя.
Если мы рассматриваем кубическое приближение {рис.2, кривая 3}, то в амплитудном
уравнении {28} надо положить a = 0. Имеем
da , 3 . {40}
—=ca+Aa dt
Это уравнение Ландау [13], которое встречается при изучении автоколебаний и турбулентности [8, 13, 14], в теории гидродинамической неустойчивости [15 - 17], в теории бифуркаций и теории катастроф [12,18]. Арнольд В.И., например, указывает: " В работе 1943 года о турбулентности, Ландау прямо выписывает уравнение бифуркации Хопфа для амплитуды теряющего неустойчивость колебания" [12, стр.96 - 97]. Арнольда В.И. , впрочем, уточнил А.С. Монин: " Отметим, что А.А.Андронов открыл эту бифуркацию на 13 лет раньше Хопфа" [14, стр.97]. Имея ввиду это уточнение, мы будем называть в дальнейшем превращение фокуса в предельный цикл - нормальной бифуркацией Андронова-Хопфа. В противном случае мы будем говорить об обратной бифуркации Андронова -Хопфа.
Уравнение Ландау {38} легко решается методом разделения переменных
аа = ат
\с\а+Ла'3
и мы полагаем здесь с > 0, Т1 > Ба. Отсюда, интегрируя по времени от 0 до 1 находим а=[_^_,1/2 ■ {41}
ехр{-2\с\т} -1К5
где
а02 {42}
К =
5 \с\ + Ла<2
- постоянная, определяемая начальной амплитудой а0 . При 1 = 0 {т =0} знаменатель в {41} максимален, а величина амплитуды а = а0 - минимальна. Затем, с течением времени знаменатель в {41} уменьшается, а амплитуда а растет. Наконец, в момент времени
t = ——ln
е 2|c|Wo
1 + Л
{43}
происходит катастрофа типа взрывной неустойчивости, при которой амплитуда а достигает бесконечно больших величин за конечный промежуток времени от 0 до te {а не бесконечный, как при линейной экспоненциальной неустойчивости [1,2]}. Происходит это потому, что скорость относительного роста амплитуд at/a = |c| + Я a2 - зависит от самих амплитуд a квадратично и рост амплитуд ускоряется. В случае линейной неустойчивости, когда l = 0 и at/a = |c| скорость роста амплитуд постоянна, ускорения этого роста отсутствуют, а бесконечно большие значения достигаются лишь при t ® ¥. Таким образом, решение {43} описывает удар землетрясения в виде взрыва амплитуд колебаний геоблока в активном разломе в определенный момент времени te, при определенных начальных амплитудах и условиях c > 0, Я >0.
Отметим, что в данном случае c > 0, Я >0 линейная и нелинейная неустойчивости действуют совместно, что и приводит к условию роста амплитуд da/dt > 0. Стационарное состояние at = 0 в этом случае невозможно. Поэтому интересно рассмотреть ситуацию, когда линейная неустойчивость отсутствует c = 0, Я >0 или даже превращается в линейную стабилизацию колебаний с < 0 при нелинейной неустойчивости Я >0.
Рассмотрим порог линейной неустойчивости, когда Ti = Da и с = 0. Уравнение at = Яа3 имеет решение
а=а (1 -1 / te f1/2, te = {2Яа2 w}-1. {44}
Это классическая взрывная неустойчивость, наблюдаемая в плазме, в пограничных слоях
гидродинамических течений, при упругих колебаниях и в других неравновесных средах [4,4
- 20]. В океане взрывной неустойчивости подвержены внутренние волны. Нарастающая
амплитуда разрушает их, образуя турбулентные пятна внутри океана, которые действительно
наблюдаются в природе [20]. В рассматриваемом случае неравновесного активного разлома,
взрывная неустойчивость проявляется в виде катастрофы - землетрясения, которое
вызывает разрушение геологических пород, самоуничтожение динамической системы и
колебаний земной поверхности. Из формулы {44} видно, что время наступления катастрофы
1/2
te зависит от начальной амплитуды а0 и частоты колебаний геоблока ш0 = {k/m} . При небольшой начальной амплитуде Оо, громоздкий и тяжелый блок с большой массой m раскачивается медленно и время наступления взрыва te затянуто. Наоборот, при больших начальных амплитудах a0 небольшие блоки обеспечивают быстрое наступление землетрясения, время te - мало.
Большой интерес также представляет докритический случай c<0, когда порог линейной неустойчивости не достигнут, Ti < Da, но нелинейный член в уравнении Ландау {40} работает, Я > 0. Интегрируя уравнение
йа А {45}
-3 = йт
-\е\а + 1а
по времени от 0 до 1 получим решение для данного случая
а= [ , К6_]1/2 • {46}
ехр (- 2\с\т) + 1Кб
где Кб - постоянная, определяемая начальной амплитудой а0
a
2
{47}
K =—^о-
6 |c| — Xa02
Важной особенностью данного случая является наличие регулярного аттрактора типа предельного цикла, на котором изменения амплитуды прекращаются и достигается стационарный режим
da Л /|С| {48} -= 0 при a =J—
dt s V l
Из решения {46} видно, что амплитуды a стремятся к стационарному значению as при возрастании времени t ® В отличие от случая линейной неустойчивости, где регулярный аттрактор имеет вид стационарной инфинитезимальной точки a = Оо , то есть тремора {39} и достигается на пороге устойчивости c = 0, Ti = Da , при нелинейной неустойчивости регулярный аттрактор имеет вид нелинейного тремора О = Os cos { w t + 0 } {49}
с конечной амплитудой as , не зависящей от начальных условий, но определяемой только параметрами c и X данной динамической системы. Из уравнения Ландау {45} видно, что
при a < as , X a2 < c , d a / d t < 0 , то есть колебания затухают со временем. Наоборот,
при a >as = {|c|/X}1/2 , da/d t > 0 амплитуды растут со временем, колебания оказываются неустойчивыми, происходит взрывная катастрофа - удар землетрясения.
При Ti ® Da, c ® 0, предельный цикл сжимается, так как as = {|c|/X}1/2 ® 0. Когда Ti достигает значения Da , величина с обращается в ноль, as = 0 и предельный цикл исчезает. Уравнение Ландау {45} превращается в уравнение at = X a3 и вновь возникает катастрофа: взрыв с амплитудой {44}. Явление исчезновения предельного цикла называется обратной бифуркацией Андронова-Хопфа, так как при прямой бифуркации происходит рождение цикла [18]. Наконец, при дальнейшем росте Ti > Da знак с меняется на с > 0, уравнение {45} превращается в уравнение Ландау {4о} с c > 0 и мы получаем стабильную взрывную неустойчивость, описываемую формулами {41}, {42}, {43}. Их удобно записать через стационарную амплитуду as
a2 a2 , {5°} a=[_as a0_]1/2
Os2 + Oq j exp (- 2|c |t)—Oq 2
1 a2 {51}
t =——1n{1+-^}. 1 1
e 2|c|WQ a2
Аналогично, для формул {46}, {47} при с< 0
2 2 . {52}
a=[7--]1/2
a^ — a22 j exp (- 2|c|t)+ a02
Рис.5. Неустойчивый аттрактор, соответствующий решению {52}, на фазовой плоскости {х, у = х^
Анализ решения {52} на фазовой плоскости показывает, что предельный цикл Андронова-Пуанкаре в данном случае неустойчив: траектории, лежащие внутри него наматываются на точку покоя, то есть возмущения с малыми амплитудами а < а затухают со временем {рис.5}. Траектории же вне предельного цикла сматываются с него и уходят в другие области фазового пространства, то есть возмущения с конечными амплитудами а > Оя растут со временем. Колебания типа нелинейного тремора {49} оказываются неустойчивыми. Монин А.С. [19] указывает на имеющиеся примеры обратной бифуркации в гидродинамике: течение Пуайзеля, возникновение стохастического странного аттрактора Лоренца и другие. После обратной бифуркации, при Т1 > Ба, с > 0 движения приобретают непериодический, случайный характер. Однако в этой области амплитуды становятся настолько большими, что уравнение {45} может быть уже не пригодным и надо использовать дальнейшие разложения функции трения, то есть более точное уравнение {28}, содержащее не только третью, но и пятую степень амплитуды.
Проанализируем уравнение {28}
а { с + Ла2 + аа4} = а!
Найдем стационарные точки, в которых а! = 0. Их две. Во - первых, точка покоя а = 0. А вторую точку надо определить из уравнения с + Лу + ау2 = 0, где у = а2 . Очевидно
Л ±
2 а
— 4 а
2
и корень со знаком минус нужно отбросить, так как у > 0. В оставшемся корне со знаком плюс, радикал
Л
Яа=
Л
4 а
2
с а
< ^^ , если с > 0 2а
то есть мы рассматриваем сверхкритический случай Т1 > Ба. Отсюда следует
у=-— + Яа < 0, 2 а
с
а
что невозможно, так как у > 0. Таким образом, при с > 0 вторая стационарная точка отсутствует, а решение {30} - {35} соответствует растущей со временем амплитуде, то есть описывает землетрясение. Фазовый портрет данной динамической системы -раскручивающийся от состояния покоя фокус.
В докритической ситуации, когда с < 0 (Т1 < Ба}, помимо стационарной точки покоя а = 0 имеется вторая стационарная точка
а1 =-— [-1 + с 2 а У
4 а \с\ + , {53}
—2 ]
Она описывает квадрат амплитуды предельного цикла - регулярного аттрактора. Согласно классификации [21] подобные предельные циклы называются полу - устойчивыми или двойными. Название полу - устойчивый возникает потому, что с внутренней стороны фазовые траектории уходят от него к состоянию равновесия а = 0, а с внешней стороны -наматываются на него {рис.6}. Название двойной предельный цикл связано с тем, что в некоторых случаях он может расщепляться на два: устойчивый внешний и внутренний неустойчивый цикл. Возбуждение колебаний в таком предельном цикле происходит жестко. Если начальная амплитуда О меньше предельного цикла ас, определяемого формулой {53}, то колебания затухают. И лишь при а > ас система может выйти на аттрактор а = ас, поскольку вне него амплитуды стремятся к предельному циклу. Величина ас не может быть большой, поскольку константы а и — одного порядка, а величина с на порядок меньше. В активных разломах колебания на этом аттракторе будут выглядеть как докритический неустойчивый тремор {49} с амплитудой определяемой формулой {53}.
Таким образом, общим бифуркационным параметром данной динамической системы, который изменяет число и характер предельных точек и циклов является параметр с = Т1 -Ба. Он регулирует поступление или отток энергии от внешнего источника в систему. При с = 0 в системе рождается регулярный аттрактор типа устойчивой точки в начале координат. При с > 0 она становиться неустойчивым фокусом, от которого амплитуда колебаний плавно нарастает при увеличении параметра с, то есть колебания возбуждаются мягко. Поскольку реальная динамическая система всегда подвержена слабым возмущениям, то колебания в виде тектонического тремора {39} возникают самопроизвольно под действием внешних возмущений: удаленные землетрясения, извержения вулканов, приливы и т.п. [1,2]. Порог
Рис. 6. Двойной, полу - устойчивый аттрактор {53} на фазовой плоскости {х, у = х1 }
самовозбуждения системы, соответствующий границе устойчивости, определяется равенством с > 0.
При с <0 в начале фазовых координат появляется устойчивая особая точка типа фокуса. Возникшие, случайные колебания быстро затухают под действием трения. Однако, если величина начального случайного возмущения превысит порог возбуждения ас определяемый формулой {53}, то родиться полу - устойчивый предельный цикл {рис.6}, в виде докритического аттрактора, на котором возможны неустойчивые автоколебания типа тектонического тремора. Небольшая флуктуация, уменьшающая амплитуду колебания до значения а < ас приводит к уходу амплитуды этого колебания в начало координат а = 0, то есть колебание затухает.
Из теории автоколебаний известно [8], что вблизи бифуркаций динамическая система очень чувствительна к малым изменениям параметров и к случайным возмущениям внешнего или внутреннего генезиса. Поэтому вблизи критических точек перехода автоколебательные системы характеризуются нестабильностью работы. Аналоги имеются и в других физических системах. Например, в лазерах или при фазовых переходах 1-го и 2-го рода, которые можно считать бифуркационными катастрофами. При этом становятся существенными случайные флуктуации. Классическая теория Ландау фазовых переходов, например, в первоначальном виде не учитывала флуктуации. Объясняя фазовые переходы возникновением одного минимума в свободой энергии Гиббса из двух минимумов, она правильно предсказывала основные особенности в поведении термодинамических величин. Однако численные значения этих величин, рассчитанные по теории Ландау отличались от наблюдаемых в эксперименте. Различия были небольшими. Например, теория предсказывала независимость теплоемкости газов при постоянном давлении от температуры. Эксперимент, однако, давал слабую зависимость, типа температура в степени 0,2 - 0,3. Эти различия удалось ликвидировать следующему поколению ученых, в частности Кеннету Вильсону, которые смогли учесть дальнодействующие флуктуации [22]. При этом был создан специальный математический аппарат - теория ренорм - группы [23]. Возможно, что эти методы можно использовать и для анализа шумов, возникающих в активном разломе в
стадии предшествующей землетрясению. Подобная работа может составить предмет будущих исследований.
Отметим, что анализ на фазовой плоскости является очень полезным для изучения общих свойств данной динамической системы. Однако заменить полное решение задачи он конечно не может. В данной работе мы получили аналитические решения, позволяющие рассчитывать поведение динамической системы типа сейсмически активного разлома при различных значениях внешних параметров. Эти решения будут использованы для интерпретации наблюдений землетрясений в наших последующих работах.
ЗАКЛЮЧЕНИЕ
Перечислим основные результаты, полученные в данной работе.
1. Линейная теория фрикционных землетрясений [1,2] обобщена на случай учета нелинейного трения. Скользящее трение в зоне разлома, имеющее минимум, представлено аналитической нелинейной функцией. Она разлагается в ряд Тейлора по малому параметру, представляющему собой разность скорости движения берега активного разлома { например скорости движения океанской плиты} и скорости колебаний сейсмогенного блока. Показано, что соответствующая аппроксимация хорошо описывает функцию трения в изучаемом диапазоне скоростей.
2. Получено новое уравнение колебаний блока в активном разломе. Оно отличается от известных в физике уравнений Рэлея и Ван-дер-Поля [24], которыми пользуются для расчетов автоколебаний малой амплитуды в динамических системах с функцией трения имеющей минимум. Главные коэффициенты при нелинейных членах в новом уравнении {16} оказываются положительными, в то время как в уравнениях Рэлея и Ван-дер-Поля они имеют разные знаки. Это обстоятельство является весьма существенным, потому что приводит к совершенно другим результатам.
3. Новое нелинейное дифференциальное уравнение колебаний второго порядка с помощью метода медленно меняющихся амплитуд {ММА} приведено к системе двух уравнений для амплитуд и фаз возникающих колебаний. Получено общее решение уравнения для амплитуд в квадратурах. Решение уравнения для фаз колебаний в рассматриваемом приближении ММА оказалось тривиальным. Фаза колебаний остается постоянной и совпадающей с начальной фазой.
4. Проведен анализ решения в кубическом приближении { в уравнении для амплитуд опускаются члены с пятой степенью}. Полученное решение в случае совместной линейной и нелинейной неустойчивости описывает удар землетрясения в виде взрыва амплитуд колебаний блока в активном разломе в определенный момент времени 1;е.
5. В обратном случае линейной устойчивости, но с кубической положительной нелинейностью, в динамической системе происходит прямая бифуркация Андронова-Хопфа: рождается предельный цикл, на котором амплитуда колебаний с течением времени выходит на стационарный режим. Однако полученный регулярный аттрактор, описывающий нелинейный тектонический тремор, оказывается неустойчивым. Траектории внутри него уходят в точку покоя, то есть колебания затухают с течением времени. А траектории вне этого предельного цикла уходят от него, амплитуды колебаний растут со временем, причем этот рост происходит очень быстро, с ускорением. Возникает взрыв, который проявляется в виде землетрясения.
6. На пороге линейной устойчивости происходит обратная бифуркация Андронова -Хопфа. Неустойчивый предельный цикл исчезает и вновь возникает катастрофа - взрыв, с амплитудой вычисляемой по полученным в работе формулам. Взрывная неустойчивость известна в физике. Она наблюдается в плазме, в гидродинамике, в электронных приборах СВЧ, в теории упругости. В случае активного разлома изучаемого в данной работе взрывная неустойчивость проявляется в виде катастрофы - землетрясения, которое вызывает разрушение геологических пород, самоуничтожение динамической системы и колебаний земной поверхности.
7. Выполнен анализ полного решения задачи с сохранением членов пятой степени в уравнении для амплитуд. В случае линейной неустойчивости вновь образуется землетрясение, описываемое полученным решением. Если же порог линейной неустойчивости не достигнут, то происходит прямая бифуркация Андронова-Хопфа, рождается полу - устойчивый двойной предельный цикл, описывающий неустойчивый тектонический тремор.
8. В результате проведенного анализа сделан вывод о важной роли введенного в работе параметра с, который представляет собой критерий неустойчивости данной динамической системы. Он регулирует поступление или отток энергии от внешнего источника в систему. Отмечается важная роль случайных колебаний и шумов в стадии предшествующей землетрясению.
Литература
1. Арсеньев С. А. Фрикционная теория землетрясений // Актуальные проблемы гуманитарных и естественных наук - 2014 - N 4 {63}. Часть II. С. 300-309.
2. Арсеньев С.А. Взрывная неустойчивость сейсмических волн // Актуальные проблемы гуманитарных и естественных наук - 2015 - N 12{83}. Часть VII. С. 141 - 150
3. Журавлев В.Ф. К истории закона сухого трения // Доклады РАН - 2010 - т.433 - N 1 - С.46 - 47
4. Басович А.Я., Новиков А.А. Неустойчивость в колебательных и волновых системах. В кн.: " Физическая энциклопедия". Под ред. А.М.Прохорова. М.: Изд-во " Большая Российская энциклопедия", 1992. С. 347 - 399.
5. Кайдановский Н.Л. Природа механических автоколебаний, возникающих при трении // Журнал технический физики - 1949 - Том XIX - Вып. 9 - С. 987 - 996
6. Кулон Ж. Разрастание океанического дна и дрейф материков. Л.: Недра, 1973. 232с.
7. Арсеньев С.А., Шелковников Н.К. Динамика морских длинных волн. М.: Изд.- во Московского университета, 1991. 188 с.
8. Ланда П.С. Автоколебания в распределенных системах. М.: "ЛИБРОКОМ", 2015. 320с.
9. Моисеев Н.Н. Асимптотические методы нелинейной механики. М.: Наука, 1981. 400 с.
10. Мигулин В.В., Медведев В.И., Мустель Е.Р., Парыгин В.Н. Основы теории колебаний. М.: Наука, 1978. 392 с.
11. Градштейн Н.С., Рыжик И.М., Таблицы интегралов, сумм, рядов и произведений. М.: Гос.изд-во " Физико - математической литература", 1962. 1100 с.
12. Арнольд В.И. Теория катастроф. М.: Наука, 1991. 128 с.
13. Ландау Л.Д. К проблеме турбулентности. Собрание трудов. Том 1. Под ред. Е.М. Лифшица. М.: Наука, 1961. С. 447 - 452 {ДАН СССР, 1944, т.44^ 8, с.339-342}.
14. Монин А.С., Яглом А.М. Статическая гидромеханика. Том 1. Санкт -1. Петербург: Гидрометеоиздат,, 1992. 695 с.
15. Churilov S.V., Shuchman I.G. Nonlinear stability pf a stratified shear flow a viscous critical layer// J. Fluid Mech. 1987. Vol.180. Pp. 1 - 20.
16. Маслоу C.A. Неустойчивость и переход в сдвиговых слоях. В кн.: "Гидродинамические неустойчивости и переход к турбулентности". Под ред. Суинни Х., Голлаба Дж. М.: Мир, 1984. С.209 - 271.
17. Churilov S.V., Shuchman I.G. Note on weakly non-linear stability theory of a free mixing layer. Proc.Roy.Soc. London A. 1987. Pp. 351 - 367.
18. Йосс Ж., Джозеф Д. Элементарная теория устойчивости и бифуркаций. М.: Мир,1983. 300 с.
19. Монин А.С., Озмидов Р.В. Океанская турбулентность. Л.: Гидрометеоиздат., 1981. 320 c.
20. Миропольский Ю З. Динамика внутренних волн в океане. Л.: Гидрометеоиздат., 1981. 303 с.
21. Андронов А.А., Витт А.А., Хайкин С.Э. Теория колебаний. М.: Гос. изд-во " Физико -математической литература", 1959. 916 c.
22. Стенли Р. Фазовые переходы и критические явления. М.: Мир,1980.
23. Ма Ш. Современная теория критических явлений. М.: Мир, 1980.
24. Ланда П.С. Нелинейные колебания и волны. М.: "Либроком", 2015. 552 с.