CC BY
143
Численное моделирование трехмерных динамических
задач сейсмологии
М.М. Немирович-Данченко
Институт геофизики СО РАН, Новосибирск, 630090, Россия
Предложен численный метод трехмерного расчета деформирования и разрушения на основе модели гипоупругохрупкой среды. Решена задача об излучении сейсмических волн от подошвы активного гравитационного оползня. Проанализированы колебания частиц среды в поле акустической эмиссии. Показано, что колебания точек поверхности оползня имеют преимущественно горизонтальную поляризацию. Проведено сравнение с экспериментом.
1. Введение
Поведение геологической среды, как в естественных условиях, так и при искусственных воздействиях, — это поведение сложно построенной существенно неоднородной среды. Постоянный интерес к воспроизведению (моделированию) динамических процессов, протекающих в таких средах, приводит исследователей к численным методам. Значительно выросшие за последнее десятилетие вычислительные мощности компьютеров позволяют от численного решения одномерных и двумерных задач перейти к моделированию пространственных процессов. В сейсмологии настоятельная потребность в решении трехмерных задач вызвана, в первую очередь, необходимостью адекватно описывать динамические процессы в условиях критического состояния геосреды. Это процессы, протекающие в приразломных зонах, вблизи подошвы гравитационных оползней и т. п.
Несмотря на развитие различных классов численных методов, конечно-разностные методы продолжают активно использоваться при решении прямых задач сейс-мики [1-4]. В настоящей работе простая конечно-разностная схема предлагается для решения динамических задач для произвольно-неоднородных сложно построенных сред. В качестве примера использования такой схемы моделируется излучение сейсмических волн активным гравитационным оползнем.
2. Основные уравнения и построение численной схемы
Для описания динамических процессов, происходящих при деформировании горных пород, нужно иметь в виду, что неоднородное строение и обязательное наличие дефектов (микротрещины, поры) приводит к разрушению среды. Как показано в [5], хрупкое разрушение, происходящее во многих локальных частицах среды, описывается такой же кривой ст-е, как для пластического течения. Поэтому для построения численной схемы использована система уравнений, описывающая поведение гипоупругохрупкой среды. Ранее гипоупру-гохрупкая модель использована нами для двумерного моделирования процессов деформирования и разрушения геоматериалов [6].
Для расчетов использованы первый закон Коши (уравнения движения в напряжениях)
Эстг7
= р(Х, Y, Z)щ (1)
дх]
и определяющие соотношения гипоупругой среды
(X, Y, Z)ёи. (2)
V ^стг7
Здесь Оу =—+- - 0.л ак) - ^ д — производная Яу-
1 Ощ Ои :
манна; — тензор напряжений; =— —- +-------------------—
2 дxj дхг
\ /
© Немирович-Данченко М.М., 2002
есть тензор деформаций; — спин-тензор; р — плот-
ность; щ — вектор смещения; ащ — тензор модулей упругости; t — время; XI — декартовы координаты X, Г, Z.
Хрупкость описывается как образование берегов трещин (отрыва и сдвига), свободных от соответствующих напряжений.
Вся среда, для которой проводится моделирование (расчетная область), разбивается координатными плоскостями на расчетные ячейки. Будем ниже рассматривать простейший случай разбиения с равными шагами по координатным осям — разбиения на кубики.
Элементарная ячейка приводится на рис. 1. Цифрами обозначены узлы ячейки (точки пересечения координатных плоскостей). В узлах определяются смещения и их производные, а деформации и напряжения — в центрах ячеек.
По скоростям смещений, рассчитанным по формуле
(1), можно вычислить скорости деформаций
1
£ ч
/ ди{ дич Л
—+ —-дхч
дх,-
При этом пространственные производные определяются с помощью осреднения величин по соответствующим граням кубической ячейки. Например, формула для определения компоненты ехх выглядит так:
-((их )2 + (их )з + (их )б +
(3)
-(их)5 - (их)8).
4Лх
+ (их )7 - (их )і - (их )4
По рассчитанным таким образом компонентам тензора скоростей деформаций определяются компоненты тензора напряжений (инкрементального тензора напряжений) по формуле (2).
И, наконец, имея уравнения движения и зная значения напряжений в восьми соседних ячейках можно определить новые значения скоростей смещений. Скорости смещений, как уже было сказано, вычисляются в узлах с использованием значений напряжений в восьми
элементарных ячейках, непосредственно окружающих узел.
Полученная таким образом численная схема является условно устойчивой. Для проверки устойчивости использовано число Куранта ЛЛх, где Vmax — максимальная продольная скорость во всей расчетной области.
Дискретизация среды по пространству порождает, кроме хорошо известных проблем аппроксимации и устойчивости, еще одну проблему. Так, Л.И. Слепяном [7] доказаны две теоремы: одна утверждает, что взяв исходный импульс достаточно плавным, можно обеспечить его распространение без изменения формы на наперед заданное расстояние. Вторая — в некотором смысле обратная ей — доказывает, что какую бы плавную нагрузку мы ни взяли, найдется такое время (расстояние), свыше которого проявится дискретность.
Покажем, как за счет правильного выбора исходного сигнала можно существенно снизить влияние дискретизации сплошной среды.
Примем для простоты, что выполнено равномерное пространственное разбиение среды на расчетные кубические ячейки с ребром Ах. Такая дискретная система имеет граничную (наивысшую) частоту собственных пространственных колебаний 1
К.
2Лх
(4)
Это так называемая частота Найквиста. Значение этой частоты и ее наличие в спектре начальных импульсов в значительной степени, как будет показано ниже, определяет точность численного решения.
Под спектром функции g(t) будем понимать функцию (7(ш), связанную с g( {) парой преобразований Фурье [8]
* (*) = 2^1 ° (ш)е]Ю‘
П"“ (5)
G (ю) = | g (ї )еЖ.
Рис. 1. Элементарная расчетная ячейка
Здесь используются общепринятые обозначения ш = = 2= т 5 Т — основной период колебаний. Формулы (5) содержат переменную t. Это может быть как время, что обычно и подразумевается (тогдаf— временная частота с размерностью Т-1), но может быть и пространственная переменная (тогда вместо f пишут к, это пространственная частота с размерностью Ь~1).
Обычно начальные импульсы задаются либо контактным способом, либо в виде исходного импульса как функции времени. Будем считать, что исходный сигнал — это импульс скорости смещения частиц (а не давления), тогда оба способа нагружения становятся однотипными.
Рис. 2. Прямоугольный импульс (а), его спектр Фурье (б) и изменение его вида в численном решении (в)
Пусть дана прямоугольная призма из модельного однородного материала (Ур = 3 000 м/с, У8 = 1500 м/с, р = = 2 000 кг/м3). Длина призмы вдоль осиXравна Ь. Грани с нормалями, параллельными осям OZ и ОУ, закреплены, торец X = Ь свободен от напряжений. На торце X= = 0 задается одинаковая во всех точках функция Ух (‘). Рассмотрим два варианта задания этой функции.
а) Выберем в качестве исходной функции Ух (г) прямоугольный импульс с амплитудой У0 и длительностью Т0 = 0.007 с, направленный вдоль оси X.
Длительность этого импульса в пространстве (то есть длина волны) равна
V
Л = -р = VpT0 = 21 м.
/
(6)
Функцию Ух (ї) будем представлять в виде
1
Ух (X).
Найдем пространственный спектр такой функции:
Л
Л/ 2
G(юх) = h І е-чЮхХах = У0Л-
-Л/ 2
БІП Ю х —
Л
Ю х —
(7)
Здесь юх = 2пk, к — пространственная частота.
На рис. 2, а приведен график функции Ух (х), а на рис. 2, б — амплитудный спектр, т. е. график функции \G(юх). Рассмотрим рис. 2, б подробнее. При шаге пространственной дискретизации ЛХ = 1 м ширина импульса Ух (х) составляет 21 пространственный диск-
рет: X = 21 ЛХ Условие равенства нулю функции (7) записывается в виде
Ю хЛ А 1 О
х = пи, где п = 0, 1, 2,
(8)
Частота Найквиста определяется формулой (4), на
ш х X 2п 21ЛХ
этой частоте -----=-= 10.5п. Эта точка пока-
2 4ЛХ
зана на рис. 2, б стрелкой.
Правее этой точки график спектра приводить нет смысла, так как при использовании формул (2) для дискретных функций частота Найквиста является граничной (теорема Котельникова). Заметим еще, что по формуле
(7) ордината точки графика рис. 2, б при ш = 0 численно равна площади исходного импульса У0 X.
Что получится, если описываемый импульс будет использован в численном расчете? Представим себе, что в произвольной точке плоскости X = 0 установлен датчик, постоянно записывающий скорость колебаний частицы, к которой он прикреплен. Время, требуемое упругой продольной волне, чтобы добежать до границы X = Ь, отразиться от нее и зарегистрироваться датчиком,
равно Т21 = #. На рис. 2, в приведен график функции Ур
Ух (£) для времени 0 < г < 9Т21. Хорошо видно, как сильно искажается форма и частотное содержание исходного прямоугольного импульса. Кроме того, возникают проблемы с оценкой скорости распространения такого импульса в среде — она задана и не должна меняться, но при появлении нескольких экстремумов вместо одного можно определить только групповую скорость.
Рис. 3. Амплитудные спектры для сигнала, отраженного первый раз (сплошная линия) и девятый раз (штриховая линия); 1 — максимумы в области низких частот
На рис. 3 приведены временные амплитудные спектры для сигнала, отраженного первый раз и девятый раз. Для волны, пришедшей после девятикратного отражения, хорошо видно появление новых максимумов в области низких частот, они помечены на рисунке цифрой 1.
Рассмотрим теперь пример другого исходного импульса.
Возьмем в качестве исходного импульса так называемый колоколообразный импульс (рис. 4, а):
Ух (х) = У„є(х-х-> (9)
На рис. 4, б приводится его пространственный спектр, а на рис. 4, в — расчет распространения в модельном материале (рис. 4, в построен аналогично рис. 2, в). Спектр функции (9), как и спектр всякой ограниченной функции, бесконечен, но достаточно быстро и гладко убывает, что для значений частоты, близких к частоте Найквиста, дает числа порядка 10-7 от начальной амплитуды. Рис. 4, в демонстрирует преимущества
Рис. 4. Колоколообразный импульс (а), его пространственный спектр (б) и его вид в численном решении (в)
такого импульса при численном моделировании — стабильно сохраняются амплитуда и частотная характеристика его (т.е. и форма сигнала, и его содержание, и гармонический состав).
Итак, если задавать исходный импульс, как в первом примере, амплитудный спектр его в окрестностях пространственной частоты Найквиста не только не равен нулю, но даже одного порядка с величиной У0Х, что генерирует паразитные колебания на этой частоте и приводит к значительной дисперсии упругих волн. Импульс с быстро и плавно убывающим спектром позволяет без использования специальных приемов проводить расчеты для большого числа длин волн.
Можно, конечно, и прямоугольный импульс выбрать так, чтобы частота Найквиста находилась дальше от основной части спектра. Например, если ширина исходного импульса будет 100 дискретов, то граничная частота станет 50л. Но, с одной стороны, для этого понадобится, соответственно, увеличивать число расчетных шагов, что в пространственном случае подчас технически невозможно; с другой стороны, амплитуда в окрестности частоты Найквиста будет порядка 2 % от основной амплитуды. Этих процентов хватает, чтобы возбудить колебания на резонансной частоте.
Справедливо утверждение — чем короче импульс, тем шире его спектр, и наоборот; предельный случай — дельта-функция имеет бесконечно протяженный спектр с равномерной плотностью. Технически же и при физическом моделировании, и в прикладных вопросах необходимо создавать очень короткие по времени и в то же время очень мощные импульсы. Но нужно, чтобы эти импульсы имели как можно более узкий спектр. Такая же проблема встает и в компьютерном моделировании. И решать ее можно здесь так же, как и в технике, а именно: необходимо выбирать импульс, для которого произведение ширины спектра на временную длительность минимально. Подробно этот вопрос изложен в
[8] и здесь далее не обсуждается.
Ниже приведены формулы сигналов, использование которых при конечно-разностном моделировании может быть рекомендовано. Все эти функции непрерывны вместе со своими производными.
Импульс из работы [9]:
Y = e(c(t-0»2cos(2nf (t -t0)), с = fk/2.15.
(10)
Вид сигнала сильно зависит от небольшого изменения любых параметров.
Один из многочисленных импульсов Рикера (пример взят из работы [10], он назван там «a two-loop Ricker wavelet»):
Y = 2nfVe(-(t - T ))e-2(nf (t-T ))2. (11)
Для этого импульса характерно исключительно плавное начало и плавное затухание. Обычно полагают время задержки Т = 1/f, тогда вид сигнала не зависит от изменения частоты /
Импульс в виде затухающей синусоиды имеет вид:
Y = e"P(t-T) sin(2nft),
e=4, t=1/f.
qT
(12)
Варьируя в показателе экспоненты величину q (а это не что иное, как добротность), можно менять число видимых периодов сигнала.
Выше мы полагали, что воздействие задано одинаковым во всех точках плоскости X = 0, т.е. рассмотрели пространственную задачу с одним актуальным измерением. Это было сделано исключительно для простоты демонстрации изменений в спектре исходного сигнала при многократном прохождении через расчетную призму. В случае же необходимости задавать начальные импульсы только в части поверхности расчетной области нужно иметь в виду, что все, сказанное выше, будет в полной мере относиться и к пространственному распределению начальных импульсов.
Идеальным часто считается точечный источник, то есть в одной точке плоскости X = 0, в точке ^0, 20) задается сигнал У^): V^, 2, г) = 8^0, 20)У(г), 8 — дельта-функция Дирака. В этом случае пространственный спектр воздействия есть прямая линия, и на всех частотах возбуждаются колебания одной амплитуды, в том числе и на частоте Найквиста. Чтобы сгладить эти численные артефакты, нужно сделать амплитудный спектр быстро убывающим. Хороший результат дает уже исследованная колоколообразная функция:
f Z) = exp
(Y - Y0)2 + (Z - Z0)2 2d 2
. (13)
При такой записи точка ^0, 2 0) — центр импульсного источника, а d имеет смысл полуширины («радиуса») источника.
В заключение раздела заметим, что первую из упомянутых теорем Слепяна можно сформулировать следующим образом: «Чтобы уменьшить влияние дискретизации сплошной среды на распространение упругой волны, необходимо выбрать такой исходный сигнал, пространственный спектр которого плавно убывает так, что на частоте Найквиста амплитуда его близка к нулю».
3. Расчет волнового поля в однородной изотропной среде от действия сосредоточенного источника
В качестве тестовой задачи была решена задача Лэмба [11]. На небольшую область верхней грани призмати-
200 х 200 х 100 м. Шаг дискретизации — 1м. В качестве исходного импульса был взят импульс Рикера (формула
(11)) с частотой 150 Гц и сглаживанием по верхней грани призмы по формуле (13).
На рис. 6 приводится численный снимок векторного поля скоростей смещений для части расчетной среды, вырезанной из первоначальной призмы по плоскостям, показанным на рис. 5 тонким пунктиром. Момент времени, для которого показано поле скоростей, выбран так, чтобы продольные и поперечные колебания хорошо разделились.
На рис. 6 хорошо видно, как поляризованы колебания частиц в тех или иных волнах. Заметим, что если источник находился бы не в центре верхней грани призмы, а, например, на краю, то вдоль ближней к источнику грани побежала бы еще поперечная волна с колебаниями, поляризованными в горизонтальной плоскости ^Н-волна).
Коническая волна—это, собственно, результат взаимодействия продольной волны со свободной поверхностью, из-за чего происходит переизлучение попереч-
Рис. 6. Численный снимок векторного поля скоростей смещений в фиксированный момент времени: 1 — фронт продольной волны; 2 — фронт поперечной волны с колебаниями, поляризованными в вертикальной плоскости ^У-волна); 3 — коническая волна; 4 — волна Рэлея
Рис. 5. Геометрия тестовой задачи. Стрелкой сверху показано действие источника
ческого тела действует источник возмущения (рис. 5). Все остальные грани тела свободны от напряжений. Среда однородна и изотропна, скорость продольной волны У,^ = 3 000 м/с, скорость поперечной волны У8 = = 1500 м/с, плотность 2000 кг/м3, размеры призмы
г\
Рис. 7. Геометрия задачи об активном гравитационном оползне
ной волны. Последняя бежит, очевидно, с меньшей скоростью и постепенно отстает, в результате чего ее фронт в пространстве выглядит как часть конуса, касательного к фронту поперечной волны.
Таким образом, можно утверждать, что предложенный метод адекватно описывает возбуждение и распространение упругих волн в пространственных моделях.
4. Численное моделирование сейсмоакустической эмиссии от формирующейся поверхности скольжения активного гравитационного оползня
Как показывает опыт, гравитационные оползни образуются в результате некоторых интенсивных воздействий естественного или техногенного характера. Как правило, для таких оползней характерно так называемое «зеркало скольжения» — некоторый слой, часто очень тонкий, обладающий пониженным сопротивлением сдвигу. Эти очевидные соображения и легли в основу численного моделирования оползневого процесса в той его части, которая, на наш взгляд, достаточна для получения предварительных представлений о характере акустической эмиссии, сопутствующей активному оползню.
Для численного моделирования тело оползня было представлено в виде однородного изотропного прямоугольного параллелепипеда с квадратным основанием (рис. 7). К верхней грани параллелепипеда, то есть к плоскости Z = 0, приложено постоянное по времени и одинаковое во всех точках грани усилие, направленное вдоль оси Y. Боковые грани модели свободны от напряжений. При расчетах вся область разбивается на кубические ячейки, каждая из которых является элементом расчетной сетки.
Нижняя грань моделирует формирующуюся поверхность скольжения. Для имитации ослабленной зоны предполагается, что в начальный момент времени половина ячеек на нижней грани может свободно без трения скользить в направлении оси Y, остальные ячейки жестко закреплены. Координаты закрепленных ячеек задаются случайным образом.
При расчетах принимается, что по достижении в жестко закрепленной ячейке некоторого порогового (критического) значения напряжения (компоненты аУ2 тензора напряжений), она получает возможность свободно двигаться в направлении Y, то есть возникает трещина сдвига. Конкретная величина критического напряжения в нашей постановке задачи принципиальной роли не играет, так как влияет лишь на амплитуду излучаемых импульсов. Таким образом, задача решается в динамической постановке с переменными граничными условиями на грани Z = Н.
При расчетах были выбраны следующие параметры модели: толщина параллелепипеда Н = 20 м, сторона его основания L = 200 м, скорость продольной волны Ур = 1 000 м/с, скорость поперечной волны У8 = 500 м/с, плотность р = 2200 кг/м3. Шаг расчетной сетки (длина ребра кубических ячеек) равнялся 0.5 м.
Численное моделирование показало, что нагружение верхней грани приводит к концентрации напряжений на нижней грани. По достижении в некоторой расчетной
Рис. 8. Численный снимок векторного поля скоростей смещений для некоторого момента времени. Стрелками справа показано направление деформирования
Рис. 9. Колебания центральной точки верхней грани в двух плоскостях (поляризационные кривые) (X, У) (а) и (У, Z) (б)
ячейке критического напряжения происходит образование трещины сдвига. В свою очередь, образование трещин сдвига влечет за собой высвобождение упругой энергии, возникают своего рода источники сейсмических волн, случайным образом распределенные по подошве модели. В качестве иллюстрации на рис. 8 приведен пример расчета пространственной волновой картины в векторном виде для некоторого фиксированного момента времени — численный снимок векторного поля скоростей смещений. Этот момент соответствует удвоенному времени прихода прямой волны от первого по времени источника. Таким образом исключается влияние границ расчетной области.
Для изучения особенностей поляризации импульсов сейсмоакустической эмиссии на верхней грани модели были проанализированы результаты расчетов для трех компонент скорости смещения в точке этой поверхности с координатами (0, £/ 2,0), отмеченной на рис. 7 белым кружком. Примеры поляризационных кривых, построенных по результатам этих расчетов для плоскостей (X, У) и (У, Z) и отображающих траектории движения частиц модельной среды на поверхности Z = 0, приведены на рис. 9. Кривые наглядно иллюстрируют, что направление преимущественной поляризации колебаний примерно совпадает с направлением оси У.
Характер поляризационных кривых, приведенных на рис. 9, можно сравнить с полевыми наблюдениями. Так, Ю.И. Колесниковым с коллегами [12] проводились наблюдения акустической эмиссии на оползне в долине реки Суусамыр (Северный Тянь-Шань, Киргизия). Ими показано, что в точках наблюдения, находящихся на поверхности средней части оползня (что исключает краевые эффекты), в основном, характерными являются малая вертикальная составляющая колебаний и их суб-горизонтальная поляризация. Это очень хорошо подтверждает характер кривой на рис. 9, б.
Итак, расчеты подобного рода могут лежать в основе экспериментальных методик, направленных на обнаружение геологических тел, находящихся в состоянии активного деформирования и разрушения с излучением сейсмических волн. Подобное поведение геологических тел часто предшествует перерастанию медленного разрушения в катастрофически быстрое.
5. Выводы
В статье предложен численный метод трехмерного расчета деформирования сред с возможным разрушением и излучением сейсмических волн. В основе метода лежит система уравнений и численная реализация, описывающая поведение упругохрупких тел.
Рассмотрены спектральные свойства исходных сигналов и их связь с возможностью использования метода в расчетах для большого числа длин волн.
Решена задача об излучении сейсмических волн от подошвы активного гравитационного оползня. Показано, что колебания точек поверхности оползня должны иметь преимущественно горизонтальную поляризацию.
Работа поддержана Российским фондом фундаментальных исследований (грант № 00-05-65337-а).
Литература
1. Inoue T, Miyatake T. 3-D simulation of near-field strong ground motion
based on dynamic modeling // Bulletin of the Seismological Society of America. - 1995. - V. 88. - No. 6. - P. 1445-1456.
2. Krebes E.S., Quiroga-Goode G. A standard finite-difference scheme for the time-domain computation of anelastic wavefields // Geophysics. - 1994. - V. 59. - No. 2. - P. 290-296.
3. Magistrate H., Day S. 3-D simulations of multi-segment thrust fault rupture // Geophysical Research Letters. - 1999. - V. 26. - No. 14. -P. 2093-2096.
4. Vai R., Castillo-Covarrubias J.M. Elastic wave propagation in an irregularly layered medium // Soil Dynamics and Earthquake Engineering. - 1999. - V. 18. - P. 11-18.
5. Zapperi S., Vespignani A., Stanley H. Eugene. Plasticity and avalanche behaviour in microfracturing phenomena // Nature. - 1997. - V. 388. -P. 658-660.
6. Немирович-Данченко М.М. Модель гипоупругой хрупкой среды: применение к расчету деформирования и разрушения горных пород // Физ. мезомех. - 1998. - Т. 1. - № 2. - С. 107-114.
7. СлепянЛ.И. Нестационарные упругие волны. - Л.: Судостроение, 1972. - 374 с.
8.Харкевич А.А. Спектры и анализ. - М.: Гос. изд-во технико-теоретической литературы, 1957. - 236 с.
9. Seriani G., Priolo E., Carcione J., Padovani E. High-order spectral element method for elastic wave modeling // 62-th Annual International Meeting and Exposition, Expanded Abstracts. - Tulsa: Society of Exploration Geophysicists, 1992. - P. 1285-1288.
10. Nielsen P. Numerical modelling of seismic waves: on the elimination of grid artifact / Norsk Hydro Research Center, N-5020, Bergen, Norway, 1994. - 47 p.
11. Lamb H. On the propagation of tremors over the surfase of an elastic solids // Phil. Trans. Roy. Soc. of London. S. A. - 1904. - V. 203. -P. 1-42.
12. Колесников Ю.И., Бабушкин С.М., Дучков А.Д. и др. Изучение геофизическими методами структурных и геодинамических особенностей оползневого склона в долине р. Суусамыр // Геология и геофизика. - 2001. - Т. 42. - № 10. - С. 1574-1584.
