Математическое моделирование распространения акустико-гравитационных и сейсмических волн в неоднородной модели «Земля-атмосфера» при наличии стратификации ветра в атмосфере Текст научной статьи по специальности «Физика»

Математические заметки СВФУ Апрель—июнь, 2015. Том 22, №2

УДК 517.956.3

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ АКУСТИКО-ГРАВИТАЦИОННЫХ И СЕЙСМИЧЕСКИХ ВОЛН В НЕОДНОРОДНОЙ МОДЕЛИ «ЗЕМЛЯ —

АТМОСФЕРА» ПРИ НАЛИЧИИ СТРАТИФИКАЦИИ ВЕТРА В АТМОСФЕРЕ А. А. Михайлов, В. Н. Мартынов

Аннотация. Рассматривается численный алгоритм и приводятся результаты математического моделирования распространения сейсмических и акусто-гравитацион-ных волн для совмещенной математической модели «Земля — Атмосфера» при наличии стратификации ветра в атмосфере. Распространение сейсмических волн в упругой среде записывается системой уравнений первого порядка теории упругости через взаимосвязь компонент вектора скорости смещений и компонент тензора напряжений. Система уравнений, описывающая распространение акусто-гравита-ционных волн в изотермической неионизированной атмосфере, записывается на основе уравнений Навье — Стокса через взаимосвязь компонент вектора скорости смещений, давления и изменения плотности воздуха при наличии стратификации ветра. Для численного решения поставленной задачи используется метод комплек-сирования интегральных преобразований Лагерра и Фурье с конечно-разностным методом.

Ключевые слова: уравнения Навье — Стокса, конечно-разностный метод, преобразование Лагерра, акусто-гравитационные волны, сейсмические волны.

Mikhailov A. A., Martynov V. N. Mathematical modeling of acoustic-gravity and seismic waves propagation for a heterogeneous Earth-Atmosphere model in a wind-stratified atmosphere.

Abstract: A numerical-analytical algorithm for seismic and acoustic-gravity waves propagation is applied to a heterogeneous "Earth-Atmosphere" model. Seismic wave propagation in an elastic half-space is described by a system of first-order dynamic equations of elasticity theory. The propagation of acoustic-gravity waves in the atmosphere is described by the linearized Navier-Stokes equations with the wind. The algorithm is based on the integral Laguerre transform with respect to time, the finite integral Fourier transform with respect to a spatial coordinate combined with a finite difference method for the reduced problem. The algorithm is numerically tested for the heterogeneous Earth-Atmosphere model for different source locations. Keywords: Navier—Stokes equations, Laguerre transform, finite difference method, seismic waves, acoustic-gravity waves.

Работа выполнена при финансовой поддержке РФФИ, код проекта 14-05-00867. © 2015 Михайлов А. А., Мартынов В. Н.

Введение

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

В последнее десятилетие появились теоретические и экспериментальные исследования, в которых показана высокая степень взаимосвязи между волнами в литосфере и атмосфере. В работе [1] описан эффект акусто-сейсмической индукции, при которой акустическая волна от вибратора благодаря явлению рефракции в атмосфере возбуждает интенсивные поверхностные сейсмические волны на расстоянии десятков километров. В свою очередь литосферные сейсмические волны от землетрясений и взрывов генерируют атмосферные акусто-гравитационные волны, которые особенно интенсивны в верхних слоях атмосферы с малой плотностью и ионосфере. Теоретическим исследованиям волновых процессов на границе упругого полупространства с изотермической однородной атмосферой посвящено большое количество работ, отметим только работы [2, 3]. В них установлены и исследуются свойства поверхностной волны Стоунли — Шолтэ и модифицированной волны Лэмба.

В данной статье с помощью численного моделирования мы продолжаем изучение распространения сейсмических и акусто-гравитационных волн для пространственно-неоднородной модели «Атмосфера — Земля» на основе идей, впервые предложенных Б. Г. Михайленко, который инициировал, принимал непосредственное участие и поддерживал данные исследования. В статье рассматривается численный алгоритм для решения совмещенной задачи распространения акусто-гравитационных волн в стратифицированной атмосфере и сейсмических волн в неоднородной упругой среде для декартовой системы координат. Решение подобной задачи для вертикально-неоднородной модели в цилиндрической системе координат без учета ветра в атмосфере было рассмотрено в [4]. Алгоритм решения поставленной задачи основан на применении преобразования Лагерра по времени, впервые предложенного в работе [5]. Распространение акусто-гравитационных волн в изотермической Атмосфере описывается линеаризованной системой уравнений Навье — Стокса. Предполагается что плотность Атмосферы и скорость ветра зависят от высоты. Распространение сейсмических волн в упругой среде описывается гиперболической системой первого порядка в терминах скоростей вектора смещения и компонент тензора напряжений.

Описываемый в статье алгоритм построен на комплексировании интегральных преобразований и конечно-разностного метода. При этом полагается, что параметры среды (плотность и скорости продольных и поперечных волн) имеют зависимость только по двум координатам, а по третьей координате среда однородна. Данную постановку задачи принято называть 2.5-Б задачей. Применение преобразования Лагерра по временной координате для численной реализации решения задачи можно рассматривать как аналог известного спек-

трального метода на основе преобразования Фурье, где вместо частоты ш мы имеем параметр р — степень полиномов Лагерра. Однако в отличие от преобразования Фурье применение интегрального преобразования Лагерра по времени позволяет свести исходную задачу к решению системы уравнений, в которой параметр разделения присутствует только в правой части уравнений и имеет рекуррентную зависимость. Данный метод решения динамических задач теории упругости был впервые рассмотрен в работах [5, 6], а затем развит для задач вязкоупругости [7, 8] и теории пористых сред [9]. В указанных работах рассмотрены отличительные особенности данного метода от принятых подходов и обсуждаются преимущества применения интегрального преобразования Лагерра в отличие от разностного метода и преобразования Фурье по времени.

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

Система уравнений, описывающая распространение акусто-гравитационных волн в неоднородной неионизированной изотермической атмосфере в декартовой системе координат (х, у, г) при наличии ветра, направленного вдоль горизонтальной оси х, и с вертикальной стратификацией по оси г, имеет вид

д

их дих

--\-Ух-

дг

дг

дх

+

дР ~дг

+ ъх

диг

дг

дР

+ Vx

диг

дх

др 1Й

1 дР дух

ро дх д г

1 дР Ро ду ' 1 дР рд

ро дг ро '

др , др0

дх дг

др др

ТГ7 +ух7Г = ~Р° дг дх

дих диу диг дх ду дг

дро

дг

дРо

+ Р(х, у, г, г).

(1) (2)

(3)

(4)

(5)

Здесь д — ускорение силы тяжести, ро(г) — плотность невозмущенной атмосферы, со(г) — скорость звука, г>х(г) — скорость ветра вдоль оси х, и = (их, иу, иг) — вектор скорости смещения частиц воздуха, Р и р — соответственно возмущения давления и плотности под действием распространения волны от источника массы Р(х, у, г, г) = 5(т-то)/(г), где /(г) — заданный временной сигнал в источнике. Полагаем, что ось г направлена вверх. Нулевые подиндексы для физических параметров среды означают, что их значения задаются для невозмущенного состояния атмосферы. Зависимость атмосферного давления Ро и плотности ро для невозмущенного состояния атмосферы в однородном гравитационном поле можно определить как

дРо

дг

= ~Ро9, Ро(г) = Р\ ехр(-г/Н),

где Н — высота изотермической однородной атмосферы, а рх — плотность атмосферы возле поверхности Земли, т. е. при г = 0.

д

д

и

и

у

у

д

х

2

= С

и

о

г

д

д

х

г

и

Распространение сейсмических волн в упругой среде записывается известной системой уравнений первого порядка теории упругости через взаимосвязь компонент вектора скорости смещений и компонент тензора напряжений:

дщ _ \_да1к

дг

Ро дх

(6)

дсГгк дЬ

М Гтг^ + тт1-) + А5гк (Ну и. \ дхг дхк )

(7)

Здесь А(х1, Х2,хз), р(х1, х2, хз) — упругие параметры среды, ро(х1 , х2,хз) — плотность среды, Ьц — символ Кронекера, и = (и1;и2,из) — вектор скорости смещений, агц — компоненты тензора напряжений. Е(х,у,г) = Е1вх + Е2еу + Езег описывает распределение локализованного в пространстве источника, а /(Ь) — заданный временной сигнал в источнике.

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

ди ~д1

Ро дхк

диг

рд

д х1 Ро

9Ух ^

дхз

да.

гк

дЬ

дик

Ка

+

дщ

дхк

+ АЬгк ¿1у и - Ьгк К а

дагк

+ Роди г

дх1

дР дР дРо

тг: = -Рошуи - иг——

дЬ дх дг

(8) (9) (10)

Здесь Ьц — символ Кронекера, ро(х, г) — плотность среды, А(х, г), р(х,г) — упругие параметры среды, и = (и1;и2,из) — вектор скорости смещений, ац — компоненты тензора напряжений. Е(х,у, г) = Е1ех + Е2еу + Езег описывает распределение локализованного в пространстве источника, а / (Ь) — заданный временной сигнал в источнике. Полагаем, что по оси У среда однородна.

Система (1)—(5) для атмосферы получается из системы (8)—(10), если взять а 11 = <Г22 = азз = -Р, р = 0, А = с^ро, СТ12 = а1з = агз = 0, Ка = 1. Полагая в системе (8)—(10) Ка = 0, получим систему уравнений (6), (7) для распространения сейсмических волн в упругой среде.

В нашей задаче считаем, что граница раздела сред атмосфера и упругое полупространство проходит по плоскости г = хз = 0. В этом случае условие контакта двух сред при г = 0 записывается так:

иг |г=-0 = иг |

да2

г |г=+0>

дЬ

дг

+ Родиг

г=-0

!|г=-о = ауг |г=-о =

г=+о

Задача решается при нулевых начальных данных

Щ^=о = ац |4=о = Р|4=о = р|<=о = 0, г = 1, 2, 3, ] = 1, 2, 3.

(11)

(12)

V

х

а

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

2. Численный метод решения

На первом этапе решения воспользуемся конечным косинус-синус преобразованием Фурье по пространственной координате у, в направлении которой среда считается однородной. Для каждой компоненты системы введем соответствующее косинус или синус преобразование:

а

\¥(х,г,п,1)=1 Ш^у^г^)^^) П = ° 1 (13)

о

с соответствующей формулой обращения

^ 2 N ^ У/{х, у, г, Ь) = -ЦГ{х, 0, Ь) + - У^ У/{х, п, г) сое(кпу) (14)

п п ^—/

п= 1

или

^ 2 *

Ш(х, у, = — У^ Ш(х, п, г, Ь) ъ\п(кпу), (15)

п

п= 1

гдекп = ^.

Выберем расстояние а достаточно большим и будем рассматривать волновое поле до момента времени Ь < Т, где Т — минимальное время распространения продольной волны до границы г = а. В результате данного преобразования получим N +1 независимых двумерных по пространству нестационарных задач.

На втором этапе решения к полученным таким образом N + 1 независимым задачам применим интегральное преобразование Лагерра по времени вида

Шр{х,п,г)= ! р = 0,1,2,..., (16)

о

с формулой обращения

оо !

УУ(х,п,г,Я = V, Р' ,,УУр(х,п,г)1"(Ы), (17)

р=0 (Р + а)!

где 1р(Ы) — ортогональные функции Лагерра.

Функции Лагерра 1р(М) выражаются через классические ортонормирован-ные полиномы Лагерра Ьа(М) [10]. Здесь мы выбираем а (порядок функций Лагерра) целым и положительным, тогда имеет место представление

1%(М) = (Ы)^е-^Ь^(Ы).

Для удовлетворения начальных условий (12) необходимо и достаточно положить

а > 1. Кроме того, введен параметр сдвига Н > 0, смысл и эффективность применения которого обсуждается в работах [6-8].

В результате этих преобразований решение исходной задачи (8)-(12) сводится к решению N + 1 независимых двумерных дифференциальных задач в спектральной области вида

Н „ 1 (да?, даХРХ „ \ г 1 Р 1

диХ „ д у дх д,

= ад/р -

3=0

(18)

1 /дар дар

1 I ^ уг ^ .х-

— ?/' —

2 у

+

ху

Р о дг дх

(дар даР

I ^ У XX ^ У X

Ро дг дх

- кп а*,, + Ка^

диУ

р-1

к р 1 (даР даР

_ П _ I 22 Л

2 г

Л „ л / 2°**_А 1 -

дг

+ кпарг 1 + Ка

ди£

<9ж

дх Ро даРх

Еу (п)/р - к]Т,

(19)

ц=о

р-1

дх

+ роди

дауу , р' «х^г^ + РояК

дх

дет',

дх

+ роди

р-1 ц=о

^ V _ (дЛ

2<тхУ ^ \ дх ^ п'пи,х I

р-1

Л Р ^ 2

ц=о р-1

+ кпир I =

У«'

ц=о

Ка

к Р

+

^ + Ро гм + + „Р^о = _

дх

дх

+ кпиР +

дг

+ ир

дг

Ег (п)/Р - ,

ц=о

(20)

р-1

= -^аХх, (21) ц=о

= -кЕ аЦу > (22)

ц=о р-1

= -к^~>Цг, (23)

ц=о

(24)

(25)

(26) (27)

р-1

ц=о

где /р — коэффициенты Лагерра функции источника /(Ь). Коэффициенты ир, ир, ир, арх, ару, арг, ару, арг, арг, рр в формулах (18)-(27) являются функциями от переменных (п, х, г).

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

Условие контакта двух сред при г = 0 записывается как

к

р-1

ц=о

г=-о

/ , р-1 2 ц =о

аЦг + Роди

г=+о

(28)

ир|2=-о = ир|г=+о; аХг |г=-о = |г=-о = 0.

Для решения задачи (18)-(28) воспользуемся конечным косинус-синус преобразованием Фурье по пространственной координате х и конечно-разностной аппроксимацией производных по координате г со вторым порядком точности [11]. Для этого в расчетной области введем в направлении координаты г сетки и шгг+1/2 с шагом дискретизации Дг, сдвинутые относительно друг друга на Дг/2:

= {гг = гДг ; г = 0,. .., К},

1/9 = {гг+1/2 = + г = 0, .. ., К-1}.

гг+1/2

к

ц

У

V

х

V

х

V

ц

ц

На данных сетках введем оператор дифференцирования , аппроксимирующий производную ^ со вторым порядком точности по координате 2 вида

Dzu(x, z) =

1

~Az

Az\ ( Az

u I x, z H--— U l X, z--

2 V ' 2

Определим искомые компоненты вектора решения в следующих узлах сеток: pp, up(x, z), Up(x, z), aXPx(x, z), (x, z), a1pz(x, z), ap.y (x, z) £ uzi, up(x,z),a%z(x,z),a^z(x,z) £ wz1+1/2-Выбор расположения компонент в целых и полуцелых узлах сетки осуществляется на основе разностной аппроксимации уравнений (18)—(27) и удовлетворения граничному условию (28). Для верхней и нижней границ задаются граничные условия 1-го или 2-го рода для соответствующих компонент.

По координате x воспользуемся конечным косинус-синус преобразованием Фурье аналогично сделанному ранее преобразованию по координате y с соответствующими формулами обращения:

2 M

Wp(x,n,Zi,p) = —Wo{n,Zi,p)-\— У W(m,n, zi,р) cosíkmx) (29)

7Г 7Г ¿-^

m= 1

M

W (x

,p) = — W(m, n, Zi,p) sin(kmx),

7Г ' ^

(30)

m=1

где кгп = * Следует учитывать, что среда в данном направлении неоднород-

на.

В результате получим систему линейных алгебраических уравнений, кото-

рая для г и г + 7} узлов разностной сетки записывается в виде

h

M

s=0

M

M

-J2 qi (Dz°Xpz - ks°Xp + kn<Tpy) + KaJ2 r1 (vxksuX - UpDzvx)

s=0

M

P-1

Fxfp - h^uX, (31) j=0

p-1

— u

2 у

- £ q2 {Dzapz + ksaXpy - knapm) - K^ r2Vpksupy = Fy fp - h^y, (32)

s=0 M

s=0

j=0

^ -Y, 43 (Dzapzz + ksapz + knapz) + Ka

s=0

M

—(P - y^r2vxk.

P0 ^

p s u z

s=0

p-1

Fz fp - h^ul, (33) j=0

M

M

^XX - 94 (Dzup + knupy) - 45ksupx + Ka

2 s=0

s=0

M

P0gup

- Y r2Vpks

s=0

p-1

-h

j=0

j

pp

(34)

2

h

p

к М М

2ауу ~ Е^ + ^ " Е+ Ка

в=о в=о

м

Ро дир - ^х к« а;

р

«^у у

в=о

р-1

-к

ц=о

уу

(35)

ц

м

м

в=о в=о

м

м

Родир

- Е Г2Vxk.

р

в=о

96 (ХиР + = Е ^

в=о ц=о

м р- 1

- Е №их + к«ир) = -к ^ а;

в=о ц=о

м р- 1

р-1

2

ц, ху

ц, Х2 1

оауг~12 98 + = _/г Е ^

Ка

где

м

в=о м

ц=о

7,РР Г^хкеРР + Е ® + + А»«?) + = ^ рР

в=о в=о

р-1

, (36)

ц=о

(37)

(38)

(39)

р-1

ц=о

(40)

о о

г1 = J 008 (к«х) 8ш (ктх) ¿х, г2 = J 8т (к«х)сов(ктх) ¿х, оо

о о

11

/ —--- эт (/г8ж) эт (ктх) ¿х, д2 = / —;-г соэ (/г8ж) соэ (ктх) ¿х

} ро(х,гг) } ро(х,гг)

9з =

Ро(х,гг+1/2)

008 (к«х) 008 (ктх) ¿х,

о

«4 = У А(х, гг) 008 (к«х)соз(ктх) ¿х, о

о

«5 = J [А(х, гг) + 2р(х, гг)] сов (к«х) сов (ктх) ¿х, о

о

«6 = J р(х,гг) 8т (к«х) 8т (ктх) ¿х,

к

о

1

о

q7 = j ^(ж, zi+1/2) sin (ksx)sin(kmж) dx, 0

b

qg = j ^(ж, Zj+i/2) cos (ksx) cos (kmx) dx,

qg = J Po(x,Zi)cos(ksx)cos(kmx) dx, km = ^-, ks='^-.

0

В формулах (31)—(40) приняты обозначения u,x = u£(m, n, Zj). Для других компонент аналогично. Черта над символом каждой компоненты поля означает, что рассматриваются коэффициенты преобразования Фурье по координате ж.

В результате всех сделанных преобразований получим N +1 систем линейных алгебраических уравнений, где N — количество гармоник преобразования Фурье по координате y. Представим искомый вектор решения W в следующем виде:

W (p) = (Vo(p),Vi(p),...,Vtf (p))T,

V = (px(m = 0,...,M; Zi), ^(m = 0,...,M; Zi),û£(m = 0,...,M; Zi),... )T.

Тогда для каждой n-й гармоники (n = 0,..., N) система линейных алгебраических уравнений в векторной форме может быть записана так:

A + ^E^jW(p) = F(p- 1). (41)

Последовательность компонент волнового поля в векторе решения V выбирается с учетом минимизации количества диагоналей в матрице A. При этом на главной диагонали матрицы решаемой системы уравнений специально располагаются компоненты, входящие в уравнения системы как слагаемые, имеющие в качестве сомножителя параметр h (параметр преобразования по Лагерру). За счет выбора значения параметра h имеется возможность существенно улучшать обусловленность матрицы системы. Решив систему линейных алгебраических уравнений (41) можно определить спектральные значения для всех компонент волнового поля W (m, n,p). Затем, воспользовавшись формулами обращения для преобразования Фурье (14), (15), (29), (30) и преобразования Лагерра (17), получим решение исходной задачи (8)—(12).

3. Аспекты численной реализации

В аналитических преобразованиях Фурье и Лагерра при определении значений функций по их спектру используются формулы обращения в виде сумм с бесконечным пределом. При численной реализации необходимым условием является определение требуемого количества членов суммируемого ряда для построения решения с заданной точностью. Так, например, количество гармоник в формулах обращения преобразования Фурье (14), (15), (29), (30) зависит от минимальной пространственной длины волны в моделируемой среде и размеров

b

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

Количество членов ряда в разложении по функциям Лагерра, необходимых для определения компонент поля по формуле (17), зависит от задаваемого сигнала в источнике / (£), выбора параметра ¡г и значения временного интервала моделируемого волнового поля. Как можно определить требуемое количество гармоник и выбрать оптимальное значение параметра ¡ , подробно рассмотрено в работах [6-8].

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

Для решения системы линейных алгебраических уравнений (41) наиболее эффективным оказалось использование итерационного метода сопряженных градиентов [12,13]. В этом случае для матриц систем большой размерности не требуется хранение всей матрицы в машинной памяти. Преимуществом этого метода является также быстрая сходимость к решению задачи при условии хорошей обусловленности матрицы системы. Наша матрица как раз обладает этим свойством за счет введенного параметра ¡ . Задав нужное значение ¡ , можно существенно ускорить сходимость итерационного процесса. Выбор оптимального значения ¡ в этом случае осуществляется исходя из минимизации количества гармоник Лагерра в формуле обращения (17) и уменьшения количества итераций, требуемых для нахождения решения для каждой гармоники.

Использование преобразования Фурье по пространственной координате, в направлении которой среда считается однородной, позволяет реализовать эффективное распараллеливание решения. В этом случае на каждом процессоре будет решаться независимая задача для каждой гармоники преобразования Фурье. Дополнительно при проведении расчетов на кластерных вычислительных комплексах с малым объемом оперативной памяти, доступной одному процессору, для решения больших пространственных задач (более 100 длин волн) было осуществлено распараллеливание решения двумерной пространственной задачи. На этом этапе расчетов была реализована распараллеленная версия метода сопряженных градиентов для решения системы алгебраических уравнений для каждой гармоники Фурье. На уровне входных данных, при задании модели среды, это равносильно декомпозиции исходной области на несколько подобластей двумерной задачи по координате 2. Такой подход дает возможность распределения памяти как при задании входных параметров модели, так и при

О 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 Х,км 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 Х,км

Рис. 1.

дальнейшей численной реализации алгоритма в подобластях.

4. Численные результаты

В данной статье мы рассматриваем результаты численного моделирования для двух вариантов распространения волн в среде Земля — Атмосфера при наличии ветра. Первый вариант рассматривает случай, когда скорость ветра в атмосфере постоянная и не зависит от высоты. Во втором варианте скорость ветра в атмосфере является функцией высоты. Результаты численных расчетов волнового поля даны в виде мгновенных снимков в фиксированный момент времени на рис. 1,2.

На рис. 1 представлен результат расчета волнового поля для случая, когда скорость ветра в атмосфере постоянная и равна 50 м/с. Данная скорость ветра была задана для того, чтобы получить основные физические эффекты распространения волн, не проводя расчеты на очень большие расстояния. Заданная модель среды состоит из однородного упругого слоя и слоя атмосферы, разделенных плоской границей. Физические характеристики слоев были заданы следующими:

1) атмосфера — скорость звука со = 340 м/сек, плотность в зависимости от координаты z рассчитывалась по формуле р0(z) = pi exp(-z/H), где pi = 1.225 * 10-3 г/см3, H = 6700 м;

2) упругий слой — скорость продольной волны cp = 800 м/сек, скорость поперечной волны cs = 500 м/сек, плотность ро = 1.5 г/см3.

Для расчетов использовалась ограниченная область среды размером (x,y,z) = (80 км, 80 км, 60 км). Моделировалось волновое поле от точечного источника типа центр давления, расположенного в упругой среде на глубине 1/4 длины продольной волны с координатами (x0,y0, z0) = (40км, 40км, -0.2км). Временной сигнал в источниках задавался в виде импульса Пузырева:

m = exp sin(27r/0(i - *„)), (42)

где y = 4, ¡0 = 1 Гц, ¿0 = 1.5 сек.

О 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 X (км) 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 Х(км)

Рис. 2.

На рис. 1 изображены мгновенные снимки волнового поля в момент времени t = 5 сек. для компоненты ux (x,y,z) в плоскости XZ при y = y0 = 40 км. Левый — без ветра, правый — скорость ветра в атмосфере 50 м/с. Граница раздела упругой среды и атмосферы показана сплошной линией. Из представленных на рисунке снимков видно, что в упругой среде вместе со сферической продольной волной P и конической поперечной волной S распространяется «нелучевая» сферическая волна S*, а далее следует поверхностная волна Стоунли R. В атмосфере кроме преломленных на границе конических акусто-гравитационных волн PP и SP распространяется сферическая волна P, а далее следует поверхностная волна Стоунли. Из представленных на рис. 1 снимков волнового поля заметно влияние ветра на распространение акусто-гравитационных волн в атмосфере и поверхностной волны Стоунли, а также на волновую картину в целом. Исследования результатов расчетов волнового поля и влияние на него ветра в случае, когда скорость ветра в атмосфере постоянная, приведены в работе [14]. В том числе в [14] описан открытый в результате этих исследований эффект влияния ветра на распространение поверхностной волны Стоунли. Ранее было известно лишь про влияние ветра на распространение акусто-гравитационных волн в атмосфере. Установлено, что при наличии ветра в атмосфере скорость и амплитуда сферической волны в атмосфере и поверхностной волны Стоунли зависит от направления распространения этих волн относительно направленности вектора скорости ветра.

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

1) атмосфера — скорость звука со = 340 м/сек, плотность в зависимости от координаты z рассчитывалась по формуле р0(z) = р\ exp(-z/H), где р\ = 1.225 * 10-3 г/см3, H = 6700 м;

2) упругий слой — скорость продольной волны cp = 450 м/сек, скорость поперечной волны cs = 300 м/сек, плотность р0 = 1.5 г/см3.

Для расчетов использовалась ограниченная область среды размером (x, y, z) = (40 км, 40 км, 33 км). Моделировалось волновое поле от точечного источника

типа центр давления, расположенного в упругой среде на глубине 1/4 длины продольной волны с координатами (жо,уо,^о) = (20км, 20км, -0.12км). Временной сигнал в источнике задавался по формуле (42). Скорость ветра в атмосфере была задана в виде функции

V(г) = 50 • ехр( —10 • (г - 3800)2) - 50 • ехр(-10 • (г - 7500)2) м/сек.

На рис. 2 изображены мгновенные снимки волнового поля в момент времени 4 = 40 сек для горизонтальной их компоненты скорости смещений в плоскости при у = уо = 20 км. Левый — без ветра, правый — с ветром. Граница раздела упругой среды и атмосферы показана сплошной линией.

Из представленного на рис. 2 снимка волнового поля без ветра (левый снимок) видно, что в атмосфере распространяются акусто-гравитационные коническая РР и сферическая Р волны, а далее следует поверхностная волна Сто-унли Д. На снимке волнового поля с ветром (правый снимок) видно, что в атмосфере возникают рефрагированные акусто-гравитационные волны Рг. Их возникновение обусловлено изменением скорости ветра с высотой. Падая на границу раздела атмосфера — литосфера, эти волны генерируют соответствующие продольную РгР и поперечную РгБ волны в литосфере и отраженную акусто-гравитационные волну в атмосфере. Это явление известно как эффект акустосейсмической индукции и описано, например, в работе [1].

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

Заключение

Предлагаемый подход к постановке и решению рассмотренной задачи позволяет моделировать эффекты распространения волнового поля для единой математической модели среды «Земля — Атмосфера» и исследовать процессы возникновения обменных волн на их границе. Численное моделирование таких процессов позволяет также исследовать особенности влияния ветра на распространение акусто-гравитационных волн в атмосфере и поверхностных волн Стоунли. Анализ тестовых расчетов показывает устойчивость представленного алгоритма даже для моделей сред, имеющих резкоконтрастные границы раздела слоев или содержащих тонкие слои, сравнимые с пространственной длиной волны.

ЛИТЕРАТУРА

1. Алексеев А. С., Глинский Б. М., Дряхлов С. И. и др. Эффект акустосейсмической индукции при вибросейсмическом зондировании // Докл. АН. 1996. Т. 346. № 5. С. 664—667.

2. Гасилова Л. А., Петухов Ю. В. К теории поверхностных волн, распространяющихся вдоль разных границ раздела в атмосфере // Изв. РАН. Физика атмосферы и океана. 1999. Т. 35. № 1. C. 14-23.

3. Разин А. В. Распространение сферичного акустического дельта-импульса вдоль границы газ — твердое тело // Изв. РАН. Физика Земли. 1993. № 2. C. 73-77.

4. Михайленко Б. Г., Решетова Г. В. Математическое моделирование распространения сейсмических и акустогравитационных волн для неоднородной модели Земля — Атмосфера // Геология и геофизика. 2006. Т. 47. № 5. C. 547-556.

5. Mikhailenko B. G. Spectral Laguerre method for the approximate solution of time dependent problems // Appl. Math. Let. 1999. N 12. P. 105-110.

6. Konyukh G. V., Mikhailenko B. G., Mikhailov A. A. Application of the integral Laguerre transforms for forward seismic modeling //J. Comput. Acoustics. 2001. V. 9, N 4. P. 1523-1541.

7. Mikhailenko B. G., Mikhailov A. A., Reshetova G. V. Numerical modeling of transient seismic fields in viscoelastic media based on the Laguerre spectral method // Pure Appl. Geophys. 2003. N 160. P. 1207-1224.

8. Mikhailenko B. G., Mikhailov A. A., Reshetova G.V. Numerical viscoelastic modeling by the spectral Laguerre method // Geophys. Prospecting. 2003. N 51. P. 37-48.

9. Имомназаров Х. Х., Михайлов А. А. Использование спектрального метода Лагерра для решения линейной двумерной динамической задачи для пористых сред // Сиб. журн. индустр. математики. 2008. Т. 11. № 2. С. 86-95.

10. Суетин П. К. Классические ортогональные многочлены. M.: Наука, 1974.

11. Virieux J. P-, SV-wave propagation in heterogeneous media: velocity-stress finite-difference method // Geophysics. 1986. N 51. P. 889-901.

12. Saad Y., Van der Vorst H. A. Iterative solution of linear systems in the 20th century // J. Comput. Appl. Math. 2000. N 123. P. 1-33.

13. Sonneveld P. CGS, a fast Lanczos-type solver for nonsymmetric linear system //J. Sci. Statist. Comput. 1989. N 10. P. 36-52.

14. Михайленко Б. Г., Михайлов А. А. Численное моделирование распространения сейсмических и акусто-гравитационных волн для модели «Земля — Атмосфера» при наличии ветра в атмосфере // Сиб. журн. вычисл. математики. 2014. Т. 17. № 2. C. 149-162.

Статья поступила 10 сентября 201-5 г.

Михайлов Александр Анатольевич, Мартынов Валерий Николаевич Институт вычислительной математики и математической геофизики СО РАН пр. Лаврентьева, 6, Новосибирск, 630090 alex_mikh@omzg. sscc. ru, vnm@nmsf .sscc.ru