Расчет сейсмического волнового поля в слое методом контурного интегрирования

Расчет сейсмического волнового поля в слое методом контурного интегрирования ...

ЛИТЕРАТУРА
1. Купрадзе В.Д. Методы потенциала в теории упругости. М.: Физматгиз, 1963, 472 с.
2. Аки К., Ричардс П. Количественная сейсмология: Теория и методы. В 2-х т. Т.1 Пер. с англ. М.: Мир, 1987, 520 с.
3. Аки К., Ричардс П. Количественная сейсмология: Теория и методы. В 2-х т. Т.2 Пер. с англ. М.: Мир, 1987, 360 с.
4. Шерифф Р., Гелдарт Л. Сейсморазведка: В 2-х т. Т.1 Пер. с англ. М.: Мир, 1987, 448 с.
5. Саваренский Е.Ф. Сейсмические волны. М.: Недра, 1972

ВВЕДЕНИЕ
Сейсмология - это наука, в которой фактические данные представляют собой записи механических колебаний Земли, называемые сейсмограммами. Колебания могут быть вызваны искусственными взрывами или естественными причинами – землетрясениями и извержениями вулканов. Оба этих природных явления привлекали внимание человечества в течение многих веков, и даже сегодня вызывают наравне с научным интересом чувство мистического страха.
Огромный прогресс, достигнутый сейсмологией за последние сто лет, был стимулирован главным образом получением постоянно улучшающихся данных. Наиболее важные шаги в этом направлении были сделаны учеными, хорошо владеющими методами математической физики. Каждое поколение сейсмологов продвигалось к получению количественных данных, преодолевая барьеры на пути расчето в сначала с помощью механических калькуляторов, а в последнее время благодаря достижениям в области цифровых микропроцессоров.
Результаты современной сейсмологии используются в инженерном деле для проектирования сейсмостойких сооружений, в разведке полезных ископаемых и при поиске месторождений нефти и природного газа. Другая область использования сейсмологических данных связана с важными политическими, экономическими и социальными проблемами, касающимися обнаружения ядерных взрывов, и с сокращением возможного сейсмического риска путем выявления сейсмически опасных участков при строительстве крупных промышленных сооружений и плотин.

Используя формулы (4.22) запишем для составляющих векторa следующие выражения:
(4.22)
Преобразуем (4.19):
, где
Используя интеграл Зоммерфельда, распишем:
где ,, k - некая комплексная функция, а выбирается таким образом, чтобы:
Таким образом,
т.к. и по формуле () , то принимает вид:
(4.23)
Используя (4.5), (4.23) и обратное преобразование Ханкеля, находим :
(4.24)
Для удобства распишем:
, (4.23)
Таким образом,
(4.24)
Найдем :
Выражая , получим:
(4.25)
Найдем выражение :
Находим :
(4.26)
Повторив рассуждения для случая, когда сила в точечном источнике действует параллельно оси z (отмечается индексом (3), также поставленным сверху), получим:
(4.17)
Исходя из граничных условий и условий непрерывности, при z=0, должны быть выполнены следующие условия равновесия:
()
Используя () подставим в () получим:
Исходя из граничных условий (3.4) на границе при переходе из слоя в полупространство должны быть непрерывны смещения:
()
()
Подставим в ():
1.
2.
3.
4.
5.
6.
где уравнения в нижней среде, равные
Смещения непрерывны:
Напряжения:
Итак, непрерывны:
Таким образом:
Условия непрерывности:
z=h
z=h
z=h
z=h
Правые части:
1.
2.
1
2
3
4
Для нахождения неизвестных коэффициентов составим системы уравнений
- для источника, ориентированного по оси у P и SV волн
=
= (24)
где, =h-z0 , , ;
- для SH волн
=
=
Получаем три компоненты вектора смещения:
(25)
(26)
(27)
В случае колебаний источника по оси z имеем:
для P и SV-волн (SH-волны отсутствуют)
=
= (28)
(29)
(30)
(31)
Таким образом, нахождение компонент смещения сводится к расчету интегралов. Данные интегралы будем находить, используя метод контурного интегрирования.
5. МЕТОД КОНТУРНОГО ИНТЕГРИРОВАНИЯ
Метод контурного интегрирования основан на преобразовании интеграла на промежутке в интегрирование по некоторому контуру, где медленно убывающая функция на концах заменяется экспоненциально убывающей функцией. В качестве контура выбирается кривая следующего вида:
(5.1)
Рис 5.1. Контур интегрирования
Запишем уравнение гиперболы:
, где (5.2)
Выразим коэффициенты B и P:
, , где (5.3)
Используя равенство и подставляя его в (5.2) получим:
где
Рис.5.2. Разбиение
Выражаем и :
Таким образом, кривая (5.1) принимает вид:
(5.4)
где при , а P и B находятся по формулам (5.3).
Таким образом, используя метод контурного интегрирования, интегралы (25)-(27) и (29)-(31) принимают следующий вид (например, для ):
6. ДИСКРЕТНАЯ ПОСТАНОВКА
6.1. Алгоритм решения
1) Строим прямоугольную сетку z на r расчетной области. Узлы данной сетки являются приемниками.
Pис. 6.1. Прямоугольная сетка расчетной области
2) Задаем начальное приближение U=0; N=16; h=1/8
3) Строим равномерное разбиение промежутка [-1,1] с шагом h
4) Находим значение функции k в узлах этого разбиения
5) Определяем необходимые параметры
6) Находим неизвестные коэффициенты и в системах (18), (20).
7) Используя алгоритм Ромберга, находим интегралы (25)-(27) и (29)-(31) по следующей схеме:
где (iz, ir) – координаты приемника, Nt – тип источника, ju-тип компоненты, m-итерация - подынтегральная функция
7. Если выполняется условие , то решение найдено, в противном случае берем h=h/2 и повторяем вычисления начиная с пункта 3.
6.2. Алгоритм Ромберга
Прямая задача рассчитывалась методом Ромберга. Методы приближенного интегрирования основаны на использовании геометрической интерпретации значения определенного интеграла, как площади криволинейной трапеции, ограниченной осью абсцисс, прямыми x = a, x = b и кривой f(x)
Алгоритм:
Сначала строится несколько интерполирующих полиномов, которые заменяют на интервале интегрирования подынтегральную функцию f(x). В качестве первой итерации полиномы вычисляются по 1, 2 и 4 интервалам.
Интеграл от каждого полинома с известными коэффициентами легко вычисляется аналитически. Таким образом, определяется последовательность интегралов от интерполирующих полиномов: Например, по правилу трапеций:
Из-за интерполяции по разному числу точек вычисленные интегралы несколько отличаются друг от друга. Причем, чем больше точек используется для интерполяции, тем интеграл от интерполяционного полинома ближе к искомому интегралу, стремясь к нему в пределе бесконечного числа точек. Поэтому определенным образом осуществляется экстраполяция последовательности до нулевой ширины элементарного интервала. Результат этой экстраполяции J0 принимается за нулевое приближение к вычисляемому интегралу.
Осуществляется переход к новой итерации с помощью еще более частого разбиения интервала интегрирования, добавления нового члена последовательности интерполирующих полиномов и вычисления нового (N-го) приближения Ромберга JN.
Чем больше количество точек интерполяции, тем ближе очередное приближение Ромберга к вычисляемому интегралу и, соответственно, тем меньше оно отличается от приближения предыдущей итерации. Как только разница между двумя последними итерациями |JN - JN-1| становится меньше погрешности (абсолютной, или в единицах |JN|), итерации прерываются, и JN принимается в качестве результата интегрирования.
ЗАКЛЮЧЕНИЕ
ЛИТЕРАТУРА
1. Купрадзе В.Д. Методы потенциала в теории упругости. М.: Физматгиз, 1963, 472 с.
2. Аки К., Ричардс П. Количественная сейсмология: Теория и методы. В 2-х т. Т.1 Пер. с англ. М.: Мир, 1987, 520 с.
3. Аки К., Ричардс П. Количественная сейсмология: Теория и методы. В 2-х т. Т.2 Пер. с англ. М.: Мир, 1987, 360 с.
4. Шерифф Р., Гелдарт Л. Сейсморазведка: В 2-х т. Т.1 Пер. с англ. М.: Мир, 1987, 448 с.
5. Саваренский Е.Ф. Сейсмические волны. М.: Недра, 1972
Электронные информационные ресурсы:
6.
ПРИЛОЖЕНИЕ 1
ПРОГРАММНАЯ РЕАЛИЗАЦИЯ
П.1. Общее описание программы
Для расчета сейсмического волнового поля методом контурного интегрирования был написан комплекс функций, позволяющий проводить численные эксперименты решения задачи для различных входных параметров.
Для разработки данного комплекса была выбрана математическая среда MatLab. Все функции написаны на m-языке.
Описание функций:
1) Main.m – главная функция, которая осуществляет ввод входных параметров и расчет волнового поля.
2) Ugol.m – функция, которая просчитывает угол между источником и приемниками.
3) PodI2.m – функция, которая рассчитывает подынтегральные выражения для источника, ориентированного по оси y.
4) PodI3.m – функция, которая рассчитывает подынтегральные выражения для источника, ориентированного по оси z.
5) IstPoY.m – функция, которая решает систему уравнений для источника по оси y и находит неизвестные коэффициенты.
6) IstPoZ.m – функция, которая решает систему уравнений для источника по оси z и находит неизвестные коэффициенты.
П.2. Текст программы
1) Main.m
%вводим параметры
global m l M l p kp ks M L P Kp Ks h0 d w z0
% disp('Введите входные параметры:')
% h=input('h-толщина слоя')
% z0=input('z0-координата источника')
% w=input('w-хрень какая-то')
% p=input('ro')
% vp=input('vp')
% vs=input('vs')
% P=input('P')
% Vp=input('Vp')
% Vs=input('Vs')
%вводим параметры
h0=0.5;
z0=0.2;
w=2*pi*0.02;
p=2;
vp=2;
vs=1.3;
P=2.5;
Vp=2.5;
Vs=1.6;