МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
УДК 517.958: 539.3(3)
МОДЕЛИРОВАНИЕ ОБЛАСТИ ПОДГОТОВКИ ЯПОНСКОГО ЗЕМЛЕТРЯСЕНИЯ 11 МАРТА 2011 ГОДА Боброва М.Е., Пережогин А.С.
Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034, Камчатский край, с. Паратунка, ул. Мирная, 7 E-mail: [email protected]
Выполнены численные расчеты области влияния подготовки японского землетрясения 11.03.2011. Приведены результаты моделирования деформационного поля земной коры в окрестности произошедшего события.
Ключевые слова: напряженно-деформированное состояние, модель Миндлина, деформации земной коры
(с) Боброва М.Е., Пережогин А.С., 2012
MATHIMATICAL SIMULATION
MSC 74B05; 86-04
STRESS FIELD MODELING OF JAPANESE EARTHQUAKE 11.09.2011 Bobrova M.E., Perezhogin A.S.
Institute of Cosmophysical Researches and Radio Wave Propagation Far-Eastern Branch, Russian Academy of Sciences, 684034, Kamchatskiy Kray, Paratunka, Mirnaya st., 7 E-mail: [email protected]
Numerical calculations of the Japanese earthquake 11.03.2011 influence are executed. Deformation’s field modeling of Earth’s crust at the locality of the event are represented.
Key words: stress field, Mindlin’s problem, Earth’s crust deformation
(c) Bobrova M.E., Perezhogin A.S., 2012
Введение
Процессы деформационных изменений земной коры в пределах сейсмоактивных зон связаны с подготовкой землетрясений. Выявление размеров областей влияния готовящихся землетрясений исследовалось в работе [3]. В настоящей работе с помощью математического моделирования построена зона деформационного влияния японского землетрясения 11.03.2011 г. В основе положена статическая модель деформационного поля в рамках теории упругости [2].
Постановка задачи
Рассмотрим модель земной коры в приближении упругого однородного изотропного полупространства. Во введенной декартовой системе координат полупространство совпадает с положительным направлением оси 02. Тензоры напряжений о,- и деформаций £- и вектор смещения щ удовлетворяют системе уравнений:
dOi
д х;
+ Xi = О,
1 ( д Ui д Uj
£ij = dXj + эх;
(1)
(2)
Oij — À'Eii&ij + 2ߣij,
(3)
ij
где X, - массовые силы внутри полупространства, X, д - коэффициенты Ламэ, 8, - символ Кронекера. Граничным условием для системы (1) - (3) является свободная граница г = 0: 0^=0 = огу|г=о = 0^=0 = 0. Источник в виде комбинации девяти двойных сил помещен в точку с на оси 02 [6].
Для нахождения поля напряжений можно воспользоваться представлением Га-леркина. Компоненты тензора напряжений в упругом изотропном полупространстве могут быть выражены через частные производные вектора Галеркина [4]:
= 2(1 - V) дхЛХ + (у д - 5X2) а™ н.
,3 ( д2 \
0уу = 2(1 - V) дуА¥ + ( УЛ - ду2 )ё1У Н,
,д/ д2 .
Ozz = 2(1 - v) — АZ + í vА - ^ ) divH,
Oyz = (1 - v)( д- АY + А^ - div H,
\д z д y
д y д z д2
Ox =(1 - vH ^^|- dXdzdlv H, w д д \ д2 Oxy = (1 - v) (,^АХ + эхay) - dXdy div H,
(4)
где X, У,X - координаты вектора Галеркина Н; ахх, ауу, ахх, , аУ1, о^, оху - компонен-
ты тензора напряжений; А - оператор Лапласа, V - коэффициент Пуассона.
Проведем моделирование максимальных касательных напряжений. Для выделения не только критических, но и всех других возможных уровней напряжений, воспользуемся величиной amax = max(|o1 — 02|, |02 — 03|, |о3 — 0"i|)/2 - критерием максимальных касательных напряжений, где 01,02,03 - главные значения тензора напряжений. С помощью значения максимального касательного напряжения определим относительные деформации сдвига:
= (1 + v) 0
fcmax — ^ umax W/
E
В упругом полупространстве определим следующие уровни сдвиговых деформаций £max 10—8 — 10—7;10—7 — 10—6;10—6 — 10—5;> 10—5. Значение < 10—8 величины emax соответствует уровню приливной деформации земной коры, а значение больше чем 10—5 - образованию области разуплотнения и достижению предела прочности пород. При численном моделировании для дальнейшего анализа установим уровень сдвиговых деформаций - 10—5.
Программная реализация математической модели выполнена в системе Maxima
[9].
Используя (4) и суммы векторов Галеркина для двойных сил без момента [6] получены явные решения для тензора напряжений и сдвиговых деформаций с помощью пакета аналитических вычислений Maxima [9]:
1) Задаем коэффициенты Ламэ: lambda, mu, глубину гипоцентра, компоненты тензора сейсмического момента.
2) Вычисляем компоненты тензора напряжения oxx, oyy, oxx, ozz, oyz, ozx, oxy с помощью представления Галеркина. Представим источник как комбинацию девяти двойных сил (двойная сила в направлении оси x, двойная сила в направлении оси y, двойная сила в направлении оси z, двойная сила в направлении оси с моментом относительно оси y, двойная сила в направлении оси с моментом относительно оси z, двойная сила в направлении оси у с моментом относительно оси x, двойная сила в направлении оси y с моментом относительно оси z, двойная сила в направлении оси z с моментом относительно оси x, двойная сила в направлении оси z с моментом относительно оси y) [6].
3) Вычисляем суперпозицию решений для единичных векторов.
4) Вычислим компоненты тензора напряжения oxx, oyy, oxx, ozz, oyz, ozx, oxy в каждой точке пространства, умножив на величину сейсмического момента для соответствующей силы.
5) Приведем матрицу к диагональному виду и вычислим максимальную полу-разность главных компонент тензора. Сдвиговые компоненты можно определить по формуле (5).
6) Вычисляем радиус влияния землетрясения - максимальный радиус от эпицентра до границы области со значением деформации порядка 10—8.
Численное моделирование
Для моделирования области деформации использовались следующие параметры земной коры: V=0.25, X = 3.5 ■ 1010 Па, д=3.48■ 1010 Па, р =2900 кг/м3, £=9.8 м/с2, Б=3 ■ 106 Па, а=0.5. Задана глубина источника с = 20 км и тензор сейсмического момента для японского землетрясения 11.02.2011 г. Данные получены из Гарвардского каталога землетрясений [10].
Mpp Mpt Mpr ' -1 . 450 -0.657 4.550
M = Mtp Mtt Mtr = -0.657 -0.281 2.120
Mrp Mrt Mrr 4.550 2.120 1.730
22
10
Тензор сейсмического момента задан в Н м. Координата г направлена по радиусу к центру Земли, ї - на юг и р - на восток. Ориентация сторон света выбрана следующим образом: координата г соответствует г, координата ї соответствует у и координата р соответствует х в декартовой системе координат.
Визуализация результатов моделирования выполнена с помощью пакета построения графиков £пиріо1 [8] (рис. а). Пространственное распределение области деформационного влияния по результатам моделирования хорошо согласуется с результатами радарных наблюдений за смещением земной коры. Механизм очага сейсмического события, который заложен в математическую модель, указывает на качественное совпадение с реально наблюдаемой картиной деформационных возмущений земной поверхности. Для заданной интенсивности и глубины источника размер области
Рис. Зона относительных деформации 10 5 на поверхности полупространства (а),
данные спутника Епу1эа1 [7] (б)
дилатансии на земной поверхности (см. рис. б) составляет десятки километров.
Заключение
Результаты применения статической модели теории упругости для моделирования областей деформационного влияния землетрясений качественно согласует с экспериментальными данными о смещениях земной поверхности. В связи с этим появляется возможность оценок размеров областей подготовки землетрясений в пределах сейсмоактивных регионов.
Библиографический список
1. Аки К., Ричардс П. Количественная сейсмология: Теория и методы. М.: Мир, 1983. Т. 1. 520 с.
2. Боброва М.Е., Пережогин А.С. Моделирование поля деформаций и зон дилатансии в упругом полупространстве с комбинацией двойных сил // Вестник КРАУНЦ. Физ-мат. науки. 2011. № 1 (1). C. 29-34.
3. Добровольский И.П. Математическая теория подготовки и прогноза тектонического землетрясения. М.: Физматлит, 2009. 240 с.
4. Новацкий В. Теория упругости. М.: Мир, 1975. C. 302.
5. Пережогин А.С., Шевцов Б.М. Модели напряженно-деформированного состояния горных пород при подготовке землетрясений и их связь с геоакустическими наблюдениями // Вычислительные технологии. 2009. Т. 14. № 3. С. 48-58.
6. Mindlin R.D., Cheng D.H. Nuclei of Strain in the Semi-Infinite Solid // Journal of Applied Physics. Vol. 21. 1950. P. 926-930.
7. URL: http://www.esa.int/esaCP/SEM9PL6UPLG_index_0.html
8. URL: http://www.gnuplot.info
9. URL: http://maxima.sourceforge.net/ru
10. URL: http://www.globalcmt.org/CMTsearch.html
Поступила в редакцию / Original article submitted: 24.11.2012