Метод расчета упругих волн в слое с помощью контурного интегрирования Текст научной статьи по специальности «Математика»

ВЕСТНИК ЮГОРСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2011 г. Выпуск 3 (22). С. 72-80

УДК 519.632+ 550.344.094

МЕТОД РАСЧЕТА УПРУГИХ ВОЛН В СЛОЕ С ПОМОЩЬЮ КОНТУРНОГО ИНТЕГРИРОВАНИЯ

В. В. Воронин

1. Постановка задачи

В работе разрабатывается метод расчёта поля стационарных упругих волн в слоистооднородной изотропной упругой среде, состоящей из слоя со свободной границей, лежащего на подпространстве.

Пусть слой {0 < z < h} заполнен упругой средой с параметрами Ляме Я, // и плотностью р . На глубине z0 в нем (то есть в точке с координатами (0,0, z0)) располагается стационарный точечный источник, в котором сила величиной 2п действует либо по горизонтали, либо по вертикали. Будем считать, что это направление либо вдоль декартовой оси Ox, либо вдоль Oz. Причём будем считать ось Oz направленной вниз, как это часто делают в геофизике.

Нижнее полупространство {z > h} также заполнено упругой средой. Величины, относящиеся к ней, будем обозначать большими буквами: константы Ляме равны Л, М и плотность Р. На границе раздела сред действуют условия плотного контакта: равенство смещений и напряжений с обеих сторон границы.

Как всегда делают в стационарном случае, вектор смещения среды в точке R = (х, y, z) представим в виде U(R, t) = Re(w(R)e~° ) , где о - частота колебаний, а t - время.

Пусть и 0( R) - поле комплексных амплитуд смещений, которые возникали бы в однородном пространстве, целиком заполненным средой с параметрами Я, /, р (оно задаётся известными формулами, например, в [2]). Мы будем искать добавочное поле в слое с амплитудой и (R), возникающее благодаря присутствию границ при z = 0 и z = h, и этим же символом обозначать полное поле в нижнем однородном полупространстве. Таким образом, поле в слое есть и 0(R) + и (R), а в нижнем полупространстве - и (R) .

Тогда комплексная амплитуда и (R) подчиняется в слое уравнению Ляме [1, 2]:

[/А + (Я + /) graddiv + ро2]и =0. (1)

(В нижнем полупространстве, соответственно, уравнению с параметрами Л, М , Р ).

Как известно [1], поле смещений может быть представлено в виде потенциальной и со-леноидальной составляющих, которые отвечают волнам продольным и волнам поперечным:

и = (м1;и2;и3) = gradф + rotyV + rotyH. (2)

Первое слагаемое здесь отвечает продольной P-волне, а второе и третье, соответственно, поперечным волнам типа SV и SH. В вертикально-неоднородных средах задача распространения волн расщепляется: на задачу для P-SV волн и отдельную задачу для SH-волн [1]. Там же приведено рассуждение, согласно которому векторные потенциалы поперечных SV и SH-волн можно искать в виде:

Wv = Vx(0;0;m) = (my ;-Шх ;0), ун = (0;0; -у). (3)

Скалярный и векторный потенциалы должны при этом удовлетворять волновым уравнениям [1, 2]:

Аф + кр 2ф = 0; A^V + ks2ipv = 0; Афн + к2ун = 0, (4)

где волновые числа заданы формулами: k2 = ю2р/(Я + 2/); ks2 = о2р / /. Аналогичные волновые числа для нижнего полупространства будем обозначать Кр и Ks.

Поле смещений через эти функции ф, m,y будет выражаться в виде:

- = (дф д2т _ду; дф д 2m ду; дф _д2 m _д 2m) (5)

дх дzдx ду ’ ду дzдy дх ’ дz дх2 ду2 .

Такое поле смещений порождает на горизонтальной площадке с нормалью n = (0; 0; 1) напряжение, которое через введённые функции выразится в виде:

T=(t тт) 4/^+^Ui+^)+(6)

1 (T1’T2’T3) / д.- дх;,и( дz ду>,Я( дх ду ^z1 и дz \ (6)

(В нижнем полупространстве, соответственно, с параметрами Л, М). Напряжение, порожденное полем смещений и0 (R) , будем обозначать T0.

На границах слоя должны выполняться условия:

при z = 0: T + T0 = 0;

при z = h: (U + и0 )верх = Uниж ; (Т + Т0 )верх = ТниЖ . (7)

2. Сила в источнике действует по оси Ox

Перейдём в цилиндрическую систему координат (r, в, z) , где r = д/х2 + у2 , в - полярный угол в плоскости Oxy. Тогда трёхмерное расстояние от источника: R = ^r2 + (z _ z0)2 . Будем отыскивать неизвестные функции ф, m,y, описывающие добавочное поле в слое и , в виде преобразования Ханкеля (или Фурье-Бесселя [3]).

ад ад ад

ф = | Vx{k, z) kJx (kr) dk • sin в ; m = j V2 (k, z )kJ (kr )dk • sin в; у = j V3 (k, z )kJ (kr )dk • cos в. (8)

0 0 0

Аналогичные представления для поля в нижней среде пусть содержат вместо функций

Vj,V2,V3 функции Q,Q2,Q3. А функции WX,W2,W3 - пусть играют ту же роль в представлении

поля однородной среды и0 .

Подставляя выражения (8) в волновые уравнения (4), получаем:

j[V + (kp2 _k2)Vl]kJl(kr)dk = 0; j[V2 + (ks2 _k2)V2]kJl(kr)dk = 0; j[V3 + (ks2 _k2)V3]kJ1(kr)dk = 0.

0 0 0

Здесь и далее: штрих обозначает дифференцирование по переменной z.

И, следовательно, обращается в ноль содержимое квадратной скобки в каждом из интегралов. То есть функции VJ,V2,V являются решениями дифференциальных уравнений с постоянными коэффициентами.

Введем для краткости обозначения: а = Jkp1 _ k2; в = -yJks2 _ k2 . Для нижней среды аналогичные величины обозначаем А и В .

Тогда общие решения этих дифференциальных уравнений можно записать в виде:

V = b11 exp(iaz) + b12 exp(ia(h _ z)); V2 = b21 exp(iez) + b22 exp(ifi(h _ z)),

V3 = b31 exp(i fiz) + b32 exp(ie(h _ z)), (9)

где коэффициенты bj (зависящие от переменной интегрирования k ) и нужно будет отыскивать. Аналогичное представление для нижней среды запишем:

ехр(/А( 2 - к)); 0>2 = q2 ехр(/В( 2 - к)); £>3 = q3 ехр(/В( 2 - к)). (10)

Здесь в каждой формуле присутствует только одна из экспонент, так как в нижнем полупространстве должны выполняться условия излучения. То есть должны отсутствовать волны, приходящие из бесконечности или экспоненциально возрастающие по переменной г.

Заметим, что интегрирование по двойственной переменной к происходит вдоль бесконечной полупрямой. И надо считать, что под радикалами а и в понимается ветвь, где при

отрицательном значении подрадикального выражения а = г Vк2 - кг 2 ; в = ф2 - к,2 . То же

верно для соответствующих величин А и В в нижней среде.

Для получения уравнений, из которых можно получить неизвестные коэффициенты Ь и

qj, подставим интегральные представления (8) в формулы (5) и (6).

Смещения:

1 = (sin2 в - cos2 в) •I [V; + V2' + V3] Jl(kr)kdk + cos2 в • I [V; + V2' + V3 ]k2 J0 (kr )dk -

0 r 0 то

-I V3k2J0(kr)dk;

0

sin в cos в • I [ V; + V2' + V3] • (- 2 Jlk) + kJ0 (kr) 1 kdk;

u2 =

0 V r

(11)

u3 =

cos в -1 [^ + k2V2]kJ1(kr)dk.

0

Напряжения:

’ , Jl(kr) '

V 1 + v 2 + V3]

0

T; = ^{(sm2в-cos2в)•I[^ + V2 + V3 ] ——)kdk + cos2в^[^ + V2 + V3 ]k2J0(kr)dk-

0 r 0 то то л

-I V3 k2 J0 (kr)dk +1[^ + k2V2]—(cos в • J; (kr))kdk};

0 0 dx

T2 = JLI • {- sin в cos в • I [ vl + V2" + V3' ] і- 2 Jlk) + kJ0 (kr)Xdk +

0 V r

то

+ J [V1 + к 2V2]—(cos#- .^(kr ))kdk}; (12)

ад

Г3 = cos 0 - J {-A2 [V + V2 ] + (A + 2^)[V + к2V2 ]}J1 (kr)kdk.

о

Выражения, стоящие под интегралом в квадратных скобках, можно переписать, учитывая, что в силу дифференциальных уравнений

v"=-a2vi; V"=-P2V2.

Для удовлетворения условий (7) на свободной границе {z = 0} достаточно потребовать обращения в ноль скобок, участвующих под интегралами в формулах (12) (только не для добавочного поля и, а для суммарного поля и + и0):

(V + WJ-P2(V2 + W2) = 0;

-Ak2[( V + W1) + (V2 + W2) '] + (A + 2^) - [-a2 (V1 + W1) + k2 (V2 + W2)'] = 0; (13)

(V + Гзу = 0.

На границе раздела сред (2 = И) нужно потребовать, согласно (7), совпадения компонент как для напряжения, так и для смещения, для чего должны на границе совпадать также значения квадратных скобок в интегралах из формул (11).

м ■ [(VI+Ж1)' - в2 V+Ж,)] = М ■ [01 ' - В20 ];

-Лк2 ■ [(V + Ж1) + (У2 + Ж2) ' ] + (Л + 2м) ■ [-а2 (V + Ж1) + к2 (22 + Ж2 )'] =

= -Л2(Q + Q } + (Л + 2М)• (-^2QX + k2Q2 );

V + W + (V2 + W2) ' = Q + Q2; (V + Wi)' + k2 (V + W2) = Qi' + k 2Q2;

(14)

¥3 + Жз = 0з; м V + Жз) = М^ 03.

Эти соотношения (13)—(14) дают систему для определения неизвестных коэффициентов Ь и qj. Но для этого надо ещё найти правую часть системы - в эту правую часть надо перенести все слагаемые, где участвуют функции

Ж1,Ж2,Ж3, отвечающие за поле в однородной среде.

Возьмём, например, из [2], формулы для смещений в однородной среде. В случае, когда сила в источнике действует по координате х, компоненты смещения:

С — ^2 Л

пЛ

1

со2 р

k 2 exp(iksR) - д2 exp(ikpR) - exp(iksR)

R

dxz

R

R

1 д

2 f

со2р дхду

exp(iksR) exp(ikpR)

R

R

; u3

1 д

2f

со2р дхдх

exp(iksR) exp(ikpR)

R

R

Каждую компоненту можно переписать, используя интеграл Зоммерфельда [4]. Здесь следует понимать, что R = 4r2 + (z -zo)2 . Для краткости будем обозначать (z - z0) одной буквой Z.

exp(,k*vr?+Z7)=i • ? exp(/|.z > j0 (kr )kdk

sir2 + Z 2 0 4ks 2 - k2

Тогда получим выражения для компонент поля и0:

.0 ^_^{(sin2^-cos2 в) •If exp(ig|Z|) - exp(if|Zl) ] +

1 о2 n J а в r

Ы,

о р 01 a в

+ COS2 в • I f - exp(ie|Z|)^3 J0(kr)dk + I kks2 exp(ie|Z|) J0 (kr)dk}

в

в

u

2 2 со2 р

sinecosejf—J,(kr) - k 2 J 0(kr) 1-f exp(,e| Z » - exp(a| Z |)Л

в

a

kdk

u3 = „2

ор

uu

cose• sigw(z- z0)• |i[exp(ia | Z |)- exp(ie | Z |)]k2J1(kr)dk.

Сравнивая эти последние формулы с интегральными представлениями (11), нетрудно получить, что роль функций V1,V2,V3 теперь играют функции:

0

ехр(а|2|), Ж2 = . ех№\2\),

о р а о р к

1 к

К, = ---.-в ехр(в \11). (15)

сор кр

Теперь подставим выражения (9)-(10) для У1,У2,У3, Q1,Q2,Q3 в условия (13)-(14). Система при этом распадётся на систему для коэффициентов Ъп, Ъ12, Ъ21, Ъ22, я1, я2, относящихся к Р и БУ-волнам, и на систему для Ъ31, Ъ32, я3, имеющих отношение к БИ-волнам.

Первая система имеет вид:

(мп м12 Л( ЪЛ Г V л

у М21 М22 у

Я

К1 У

где Ъ =

(Ъ Л

ии

У Ъ22 у

; я =

( Я л

У Я2 у

клетки матрицы таковы:

Ы-

11

Ы

12

іа -іаешк

2цк2 - о2р (2цк2 - о2р)ешк

іцаешк -іца

(2цк2 - о2р)ешк 2цк2 - о2р

0 0 Л

0 0

-іМА МВ2

-2Мк2 + о2Р -2іМк2В

-в2 2іцк 2 в -цр2ев 2іцк 2вев

-в2ев Л -2іцк 2веівк

-цв2

-2іцк 2 в

Ы

( еіа

21

1 іве

іае‘

іак

-іа

в -івЛ

к2еів к2

(16)

Ы22 =

( -1 -іВЛ

-іА -к2

К _ " у

А правая часть, получаемая из подстановки формул (15) в соотношения (13)—(14), выглядит так:

Г П2 \

е

V

о2 р

-кеіаг0 -в— еівг0

о2р

ік еіа^ ів еівА

а к

К кеіаА - кеівА

к

(г2р-2цк2)—еаа -21цке'в а

ц(кеаА+в евА) к

(г2 р - 2 цк 2)—е1 аА - 2 I цк@е1 вА а

Здесь обозначено А = к - г0, то есть расстояние от источника до границы раздела сред. Система же для коэффициентов, относящихся к БИ-волнам, выглядит так:

(17)

(

і в _Фє:

1

iвh

-)iвh

0

_1

ІцФє

iвh ________

і ив _imb

f b л

31

b32

v q3 у

—2p

lks2 єіФz0 k e

ks єіФА

Л

iM

кв

к,.2

к

УІвА

(18)

Теперь найденные (для каждого значения переменной к) коэффициенты надо подставить в формулы (9)—(10), а результат - в формулы (11) для компонент вектора смещений и .

З. Сила действует по оси Oz

Здесь SH-волна отсутствует, то есть функция / в представлениях не участвует. Вместо представлений (8) теперь будет:

дада

ф = J V, (k, z)kJ0 (kr)dk; m = J V2(k, z)kJ0 (kr)dk. (19)

00

Вместо формул (11)-( 12):

дада

u, = cos в • J [_ V, _ V2 ]k2 J, (kr)dk; u2 = sin в • J [_ V, _ V2 ]k2 J, (kr)dk;

0

0

u

3

дада

J [V, + к2 V2 ]kJ0 (kr)dk; T, = u cos в • J [_2 V, + (в2 _ к2) V2 ]k2 J, (kr)dk;

00

да

T2 = u sin в • J [_2 V, + (в2 _ к2) V2 ]k2 J, (kr)dk;

0

да

T3 = J [_Лр2V, + 2u(_a2V, + к2V2 )]kJ0(kr)dk.

(20)

Вместо условий на границах (13)-(14), условия будут такие: при ъ = 0:

-2( V + Ж)' + (в2 - к2)(К, + Ж,) = 0;

-ЛЛр2(К1 + Щ) + 2Ц-а2(К + Ж1) + к2 (К2 + Ж2)'] = 0;

при ъ = И:

ц-2(К1 + Жху+(в2 - к2)(К + ад = щ-2оЦ+(в2 - к2)а];

-Лр2(^1 + цгг) + 2М[-а2(У1 + Щ) + к2(V + Ж2)'] = -ЛК^ + 2М[-А2£ + к^2 ];

Функции Ж1,Ж2, отвечающие за поле и0 для однородной среды, сейчас могут быть взяты в виде:

(21)

W = slgl( Z) • exp(ia | Z I); W2 = •І- ехр(ів | Z I).

~ —2P в

— 2P

Система уравнений по прежнему будет иметь вид:

f Mil M,2 ЛГ b Л V

v M21 M22 у

q

ҐїЛ

Vі у

V w у

(Ъ Л и11

где Ъ =

У Ъ22 у

; я =

( Я Л У Я2 у

клетки матрицы сейчас таковы:

Ы

11

У (

Ы-12 =

-2іа 2іаешк ц(в2 - к2) Мв2 ~-2)еівк -

2цк2 -о2р (2цк2 -ю2р)ешк 2іцк2в -2іцк2вегвк

-2іцаешк 2іца ц(в2 -к2)е,вк ц(р2-к2)

(2цк2 - о2р)еіак 2цк2 - о2р 2іцк2веівк

0 0 -

0 0

-2іцк 2 в

2іМА М(-В2 + к2)

-2М к2 +о2Р -2іМк 2В

Ы 21 =

( еіа

1 іве

івк

-ів

іае

іак

-іа

к2еівк к2

(

Ы

22

-1 -іВ

-іА -к2

А компоненты правой части таковы:

( 2аеаа +1(к2 -/2) ев

V =

о2р

в

(-о2р + 2цк2)е 0 - 2цк2е і (к2-в2) е

в0

ц(2іаешА + і

вА

в

)

У (о2р- 2цк2)еіаА + 2цк2е1 вА у

о2р

_еаА + евА

-іае

іаА

к2

в

вА

(22)

(23)

4. Контурное интегрирование

Итак, мы приходим к задаче вычисления интегралов вида:

и

ТО

| У (к) ехр( і4кр 2 - к2 г) Jv (кг )й?к.

(24)

Здесь V равно 0 или 1, а функция, обозначенная тут через У (к) содержит какой либо коэффициент из тех, что найдены при решении либо системы (16)-(17), либо системы (18), либо системы (22)-(23). Возможно с множителями в виде степени переменной к, или величины а, или в.

Заметим, что, если компоненты волнового поля предполагается рассчитывать сразу на некоторой сетке, в узлах которой расположены приёмники, то эти коэффициенты (а значит, и значение У (к) ) можно вычислить один раз - и использовать при вычислении интегралов для всей серии приёмников. То есть, можно завести массив, в компонентах которого будут накапливаться результаты численного интегрирования по каждому приёмнику - а соответствующие значения У (к) считаются один раз.

Более того, если приёмники расположены ещё и на одной глубине 2 , то и значение экспоненты, стоящей под интегралом (24), вычисляется единожды.

1

Воронин В. В. Метод расчета упругих волн в слое с помощью контурного интегрирования Теперь надо уточнить, что функция, обозначенная здесь как У (к), имеет не только слабые особенности (и точки ветвления), когда а = у1кр2 - к2 = 0 или в = ^к,,2 - к2 = 0, но и полюса. Эти полюса соответствуют волнам Релея или Стоунли для Р-БУ волн, а также волне Лява для БИ-волн [1]. Полюса образуются потому, что в некоторых точках к, лежащих на луче (0, +гс), обращаются в ноль детерминанты решаемых систем линейных уравнений.

И интеграл вида (24) надо понимать тогда в следующем смысле. Исходя из принципа предельного поглощения [1], при наличии малого поглощения полюса смещаются вверх от положительной действительной полупрямой. И значит, переходя к пределу, когда поглощение равно нулю, в интеграле надо учитывать обход соответствующих полюсов снизу (то есть включать в значение интеграла половины вычетов в соответствующих полюсах).

Находить при этом полюса и вычеты в них вовсе не обязательно. Достаточно деформировать путь интегрирования в комплексной плоскости переменной к , чтобы этот путь обходил снизу весь отрезок действительной полуоси, на котором могут находиться полюса.

Как показано в [5], полюс Лява (если он есть) находится между к, и К,, а полюса Релея

и Стоунли - не дальше, чем 1.15 к,. Следовательно, если путь интегрирования обходит снизу отрезок вещественной прямой от нуля до

к0 = 2шах{к,, К,}, (25)

то этого вполне достаточно.

Теперь представим функции Бесселя в виде:

Мкг) = Н^г) + ),

где НУ (1) и НУ (2) - функции Ханкеля 1-го и второго рода. Как известно [4], при больших значениях аргумента они имеют асимптотику:

Н <1|(кг) Чшехр(,кг - п); Н |2|(кг) Чшехр(-,кг+п).

Следовательно, на комплексной плоскости переменной к первая из них убывает экспоненциально в верхней полуплоскости, а вторая - в нижней полуплоскости.

Тогда для эффективного вычисления интеграла вида (24) его надо разбить на два:

1 гс 1 го

и = / + /2 = -1У (к )ехр(/^кр 2 - к2 г) НУ (1)(кг )йк + -1У (к )ехр(/^кр 2 - к2 г) Н У (2)(кг )йк (26)

2 0 2 0

и деформировать путь интегрирования по-разному для каждого из них.

Обозначим сейчас: к = % + щ и деформируем путь интегрирования для / в параболу, где

П = С%(% - к0) и под к0 имеется в виду величина (25). Коэффициент С возьмём таким, чтобы при заходе параболы в нижнюю полуплоскость экспонента ехр(/кг) не принимала слишком больших значений. Например, чтобы в нижней точке параболы (при к = к0 /2) значение действительной части показателя экспоненты Ке(/кг) = Ск02г /4 не превышало 8. Для этого достаточно положить

32

С

к 2г

А'0 шах

где гшах - максимальное значение величины г для той серии приёмников, где мы собирается вычислять волновое поле. Теперь сделаем замену переменной %, чтобы перевести бесконечный промежуток интегрирования (0;гс), например, в интервал (0;1) .

dk = [£'(г) + іц'(т'№т . (27)

Тем самым, первый из интегралов (26) запишется в виде:

1 1

71 = 21У (к (Г)) ехР(^кр2 - к2 (г) -г) И^1)(к (т)г)[?(т) + . (28)

0

Отметим, что хотя [%'(т) + П'Т)] ^ ГС при т ^ 1, в других сомножителях подынтегральной функции происходит экспоненциальное убывание (если одновременно не равны нулю г и (г - г0), то есть если точка наблюдения не совпадает с источником). Так что особенности

на конце промежутка интегрирования фактически нет.

Что касается пути интегрирования для второго из интегралов (26), где присутствует функция Ханкеля второго рода, то его можно деформировать так, чтобы он уходил в нижнюю полуплоскость. Например, по параболе, где ц = -С%(% + к0). И, сделав замену, аналогичную (27), а именно:

переписать его в виде, подобном (28).

Для вычисления полученных интегралов эффективным и достаточно простым является метод Ромберга [6]. Говоря кратко, он заключается в том, что вычисление интеграла ведётся простейшим методом трапеций с удвоением числе узлов при каждой последующей итерации. А в конце производится экстраполяция двух последних приближённых значений интеграла (для шага И и шага И/2) на значение шага, равное нулю.

1. Аки, К. Количественная сейсмология: теория и методы : в 2-х т. [Текст] / К. Аки, П. Ричардс. - Т.1. - [пер. с англ.]. - М. : Мир, 1987. - 520 с.

2. Купрадзе, В. Д. Методы потенциала в теории упругости [Текст] / В. Д. Купрадзе. - М. : Физматгиз, 1963. - 472 с.

3. Тихонов, А. Н. Уравнения математической физики [Текст] / А. Н. Тихонов, А. А. Самарский. - М. : Наука, 1977. - 744 с.

4. Бейтмен, Г. Высшие трансцендентные функции [Текст] / Г. Бейтмен, А. Эрдейи. - М. : Наука, 1974. - 296 с.

5. Кауфман, А. А. Введение в теорию геофизических методов. [Текст] / А. А. Кауфман, А. Л. Левшин. - Ч. 5 : Акустические и упругие волновые поля в геофизике. - Москва : Недра, 2005. - 663 с.

6. Бахвалов, Н. С. Численные методы [Текст] / Н. С. Бахвалов. - М. : Наука, 1975. - 632 с.

V. „(Т)=-_32,_т

ЛИТЕРАТУРА