Сер. 4. 2009. Вып. 2
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
УДК 550.34.01+550.348.098 С. В. Лебедев, В. А. Павлов
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ВОЗДЕЙСТВИЯ ЗЕМЛЕТРЯСЕНИЙ И МОЩНЫХ ВЗРЫВОВ НА АТМОСФЕРУ
Введение. Изучение распространения звука в атмосфере началось с зарождения акустики. В конце XVIII - начале XIX вв. У. Дарем (Англия) изучал зависимость скорости звука от скорости ветра, Бьянкони (Италия) и Ш.-М. Кондамин (Франция) изучали влияние температуры на скорость звука. На тот факт, что сила тяжести и стратификация геофизических сред сильно модифицируют распространяющиеся в них звуковые волны, учёные впервые обратили внимание в конце XIX века. С этого момента были начаты теоретические исследования свойств распространения акустических волн в атмосфере.
Задача об эволюции акустической волны в неоднородной атмосфере является одной из интересных и ёмких задач нелинейной акустики. В последнее время большой интерес вызывает численное моделирование такого рода задач. Это связано с появлением технической возможности решения такого класса задач с помощью всевозможных математических пакетов и пакетов собственно моделирования. Интерес, проявляющийся к такому роду исследованиям, в первую очередь связан с возможностью применения получаемых результатов на практике для анализа и оценки изменений в атмосфере и ионосфере, а также возможностью сравнения этих результатов с данными, полученными на практике: данные сейсмографов при мощном взрыве и землетрясении, данные со спутников или зондов. Не менее интересно решение и обратных задач: когда по косвенным признакам можно определить не только место проведения взрыва, но и охарактеризовать источник или причину, породившую такого рода изменения.
Постановка задачи. В работах [1, 2] на основании разработанной методики строится приближённое представление нестационарной акустической волны над эпицентром землетрясения с учётом сферической расходимости фронта в неоднородной атмосфере. Задача решается поэтапно (см. приложение). Для рассматриваемых импульсов на начальной стадии числа Рейнольдса велики по сравнению с единицей. Это позволяет на начальном этапе учесть влияние нелинейности при пренебрежении диссипацией. На этой стадии решение описывает эволюцию акустического импульса от уровня, на котором задаётся профиль возмущения, вплоть до высот, когда профиль приобретает асимптотическую треугольную форму (10). При подъёме вверх число Рейнольдса уменьшается. Следующим этапом решения является построение представления для величин скачка положительной и отрицательной фазы импульса для ситуации, когда мы ещё не достигли высот, где число Рейнольдса равно единице, т. е. уровня с существенными потерями. Последним этапом решения является описание эволюции положительной и отрицательной фаз треугольного импульса на высотах, где число Рейнольд-са меньше единицы и учёт диссипации производится в линейном приближении (13).
© С. В. Лебедев, В. А. Павлов, 2009
В настоящей работе представлено численное моделирование на базе пакета А^ SYS 10.0 процесса распространения акустического импульса в слое неоднородной атмосферы высотой вплоть до 200 км. Полученные результаты сравниваются с работами [1, 2], посвящённым этой проблеме и дающими аналитическое описание процессов происходящих, в атмосфере над эпицентром землетрясения с учётом слабой нелинейности и диссипации.
На основе результатов моделирования появляется возможность не только сравнить отдельные этапы эволюции, но и ответить на вопросы о том, что происходит в зоне перехода от одного этапа к другому. Тем самым удаётся описать более полно картину происходящего в слоях атмосферы. Волновой процесс описывается следующей системой уравнений газовой динамики:
Р
др
+ УР - рдеу - цДг
0;
+ ё1у(рг) = 0;
дЬ
д(Р СРТ) сЪ
= А\у(кУТ) + Ж у + Ек + Ф +
дР
дГ
(1) (2) (3)
Здесь г - скорость среды, Ь - время, Р - давление, р - плотность, ц - коэффициент динамической вязкости, д - ускорение силы тяжести, Т = 288 К, Ср - удельная теплоёмкость при постоянном давлении, к - коэффициент теплопроводности;
- слагаемое, описывающее вязкую работу, причём по по-
Ж
V
«з Ц
_а_1±_ , _й__нь.
дх^ дх^ 'Л/ ,
вторяющимся индексам подразумевается суммирование по трём ортогональным на-
правлениям («1 = гх, «2 = гу, «з
Ек = -ш
(к г,1) С' дх V 2 ( /
ду
(I г,2)
Ср ду V 2 ( /
дг
к д С„ дг
(Ь V*)
кинетическая энергия; Ф = ц ризующая диссипацию, связанную с вязкостью. Невозмущённые поля имеют следующий вид:
Эх,
+
дх-
дх.
величина, характе-
ро = ро(0) ехр(-у/Н), Ро = Ро(0) ехр(-у/Н), г = 0, Н = const,
(4)
где у - текущее значение высоты, Н - высота однородной атмосферы.
Моделируемая область задачи представляется в виде прямоугольного участка атмосферы в плоскости (х, у) лежащей в пределах от у = X = 5 • 103 м, где расположен источник возмущения, вплоть до у = 2 • 105 м по высоте. По протяжённости в горизонтальном направлении рассматриваемая зона представляет участок в 80 км, в центре которого расположен источник протяжённостью в 400 м.
Граничные условия задаются следующим образом: на вертикальных стенках задаются нулевые значения горизонтальной составляющей скорости гх =0. На нижней стенке моделируемой области, за исключением области расположения источника, задаются нулевые значения вертикальной составляющей скорости гу =0.
и(0-
¡//¡т
0
Рис. 1. Профиль возмущения скоро-
сти среды на уровне у = X
На верхней стенке задаётся давление, рассчитанное по барометрической формуле для высоты у = 2 • 105 м. На участке, соответствующему источнику, на уровне у = X задаётся возмущение скорости, изменяющееся во времени гу (Ь) (см. (6)) и имеющее профиль, изображённый на рис. 1, где г(Ь1) = тахгу(Ь, X).
т
2
Моделирование ведётся на элементе ANSYS FLOTRAN-FLUID141, позволяющем описать распределение температуры, плотности, скорости и давления в двумерных и трёхмерных потоках для однофазной среды в соответствии с законами сохранения массы, импульса и энергии (1)—(3). Использование квадратичного элемента с четырьмя узлами совместно с упорядоченной сеткой, приводит к тому, что рабочая область представляет собой набор из 390 000 квадратов со стороной 200 м. Форма и размер элемента подбирались, исходя из имеющихся характеристик компьютера, на котором производились расчёты. Увеличение количества элементов приводит к увеличению используемой памяти (RAM) и приводит к необходимости использования машин с памятью, превышающей 2046 Мб. С другой стороны, необходимо было учесть тот факт, что высотное распределение давления задавалось в узлах, т. е. увеличивая размер элемента мы автоматически загрубляем используемую модель среды. Также была учтена связь размера элемента, шага по времени и характерной скорости распространения, которая превышает скорость звука в воздухе. Рассматривалась ситуация, соответствующая «слабой нелинейности». С учётом этой связи, волна не должна проходить более одного элемента за единицу времени. В свою очередь, единица времени - шаг по времени, проходила процедуру подбора на основании критерия - точности, необходимой и достаточной для сравнения получаемых в ходе решения построенных алгоритмов результатов с реальными физическими процессами. Немаловажным фактором при выборе размера элемента явилось условие наименьшего влияния стенок моделируемой области. При достаточно близком расположении стенок в процессе распространения акустической волны наступает ситуация, когда отражённые от боковых стенок части фронта волны складываясь, порождают вторую волну.
Источники землетрясений отличаются большим разнообразием: различная мощность, глубина залегания, форма очага землетрясения, различная структура среды и т. д. По этой причине представляет интерес введение в рассмотрение некоторого «стандартного», «эталонного» землетрясения. Одним из существенных параметров представляется характеристика землетрясения по «эквивалентной» мощности, приписываемой землетрясению. Предлагается приписывать землетрясению мощность такого наземного точечного взрыва, который производит в атмосфере возмущения близкие по своему характеру к возмущениям от землетрясения. Предлагается воздействие землетрясения заменить точечным взрывом на уровне поверхности земли. При этом, на основе работ [2, 3], удаётся связать «площадь» S+ положительной фазы начального импульса и выделяемую мощность ю в килотоннах (кт) тротилового эквивалента.
где ¿1, ¿2 - время начала и конца импульса, размерность функции - длина. Для можно дать следующее описание (в метрах):
где множитель 2 обусловлен влиянием границы раздела земля-воздух. Для X = 5^10 км параметр В порядка единицы. Так, при X =10 км В « 1, при X = 4 км В « 3. При этом остаётся открытым вопрос об определении той доли энергии землетрясения, которая «выделилась» в атмосферу. На основании полученных связей можно классифицировать моделируемый процесс как землетрясение, чья выделяемая мощность лежит в интервале 103 + 104 кт. Такой способ описания землетрясения использовался ранее [3]. Если
S+(A,, ю) « 2В(^)ю1/3,
иу(1, у), 24 м/с
22 20 18 16 14 12 10 8 6 4 2 0
у -- 8 I Л0
20
Ж
25
30
шз
60
о(г1) = 40 м/с
8103 ]
и(г1) = 20 м/с
ю = 103 кт
и(г1) = 15 м/с ю = 422 кт
20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 г, с
Рис. 2. Зависимость профиля скорости уу (Ь,у) от времени для разных высот у =
^90 км над эпицентром землетрясения для трёх значений выделяемой мощности ю = 422, 103, 8- 103 кт
ю
же сравнивать с более общей классификацией, то такого рода землетрясения относятся к «обычным», наиболее часто происходящим в природе землетрясениям. Землетрясения, чья выделяемая мощность лежит в интервале 104 + 106 кт - «сильные» землетрясения. В силу того, что «сильные» землетрясения довольно редки, интерес именно к «обычным» землетрясениям наиболее оправдан с точки зрения необходимости прогнозирования и исследований.
Результаты и их обсуждение. Одним из результатов численного моделирования процесса распространения акустической волны в неоднородной изотермической атмосфере (4) является зависимость профиля скорости от высоты над источником (рис. 2).
На данном рисунке в виде импульсов представлены профили вертикальной составляющей скорости «у (¿, X), снимаемые на различных высотах, по мере распространения возмущения вверх, непосредственно над эпицентром возмущения для случая выделяемой мощности 8 • 103 кт. Необходимо отметить, что на рис. 2 приведены только положительные части профилей скорости, где через точки максимума проведена огибающая. Также приведены огибающие для случая ю = 103 кт и ю = 422 кт, полученные при выполнении аналогичных построений.
Специфика подаваемых в качестве источника возмущений заключалась в следующем. При фиксированном значении длительности положительной фазы ¿2 начального импульса за величину отражающую изменения в значениях выделяемой мощности можно принять значения амплитуды этого импульса. Таким образом, на рис. 2 приведены результаты моделирования поставленной задачи для случаев когда профили скорости, подаваемые в области источника, имеют амплитудные значения, соответственно, «(¿1) = 15 м/с, «(¿1) = 20 м/с, «(¿1) = 40 м/с и дают значения выделяемой мощности в 422, 103, 8 • 103 кт. Представляемая модель позволяет строить зависимости профиля
скорости от высоты с шагом, равным 200 м в горизонтальном и вертикальном направлениях, т. е. в каждом узле упорядоченной сетки. В данном случае кривые построены с шагом по высоте, равным 1000 м.
Характерными особенностями эволюции являются следующие.
1) Амплитуда акустического импульса начинает убывать согласно закону у-1 ехр[у/2 • Н], где у - расстояние, отсчитываемое вдоль вертикального направления, вплоть до достижения минимума на расстоянии, равном удвоенной высоте однородной атмосферы [1]. Для данной модели на участке высот от 8 • 103 км до 14103^15103 м, в зависимости от значения начальной амплитуды, затухание амплитуды описывается законом 1/у с хорошей степенью точности - это соответствует ситуации отсутствия неоднородности среды. На высотах свыше 15 • 103 м амплитуда акустического импульса продолжает монотонно убывать, подчиняясь закону у-1 ехр[у/2 • Н] пока не достигает своего минимума. Высота, на которой реализуется минимум утт, не зависит от выделяемой мощности и достигается через 50,25 с после начала распространения импульса.
2) После достижения минимума, по мере распространения возмущения вверх начинается постепенное возрастание амплитуды акустического импульса с постепенным искривлением профиля. При этом происходит укручение фронта волны и образование ударной волны. Отметим, что для интервала высот от момента достижения ут-ш до момента образования ударного фронта существенно конкурирующее влияние вязкости. Поэтому, на данном этапе был выполнен поиск высоты, на которой начинает проявлять себя вязкость. В частности, при увеличении последней в 1,5 раза для задачи с начальной амплитудой г(Ь1) = 20 м/с, Ю = 103 кт даваемое построенным алгоритмом решение не отличалось от задачи с нормальными условиями по вязкости вплоть до высот равных 6 • 104 м. Начиная уже с 6,5 • 104 м, отличия в профиле скорости, вносимые изменённым параметром вязкости, становятся все более ощутимыми. Влияние вязкости можно разделить на две составляющие: «размытие» фронта ДЬ и уменьшение амплитуды. На рис. 3а, -б приведены зависимости ДЬ от текущего значения высоты у и Дгтах также от у. Здесь ДЬ = Ь1 — ¿2, ^ - время прихода фронта в заданную точку на высоте у над эпицентром землетрясения по уровню гу (Ь,у) = 0,1 м/с - сплошная линия и гу (Ь,у) = 1 м/с - пунктирная линия для задачи с увеличенной вязкостью; Ь2 - время прихода фронта для задачи с нормальными условиями на вязкость; Дгтах = = г(Ь\) — г(Ь2) и г(Ь\) - максимальное значение профиля скорости на высоте у для задачи с нормальной вязкостью, г(Ь2) - для задачи с увеличенной в 1,5 раза вязкостью.
По мере того, как волна проходит выше уровня у = 65 км, наблюдается увеличение ДЬ, значения которой выходят на некоторый уровень вплоть до момента достижения высоты первоначального опрокидывания волны. При этом точки фронта лежащие ниже чувствуют на себе влияние вязкости раньше нежели точки фронта, лежащих выше. Из рис. 3а видно, что для точек фронта, лежащих на уровне гу(Ь,у) = 1 м/с, размытие начинается только после достижения высоты более 70 км, что на 5 км выше, чем для уровня гу (Ь,у) = 0,1 м/с. Сложившаяся ситуация такова, что чем большее расстояние пройдено волной, тем большая часть фронта чувствует на себе влияние вязкости. При этом для Дгтах наблюдаются малые колебания относительно нулевого уровня Дгтах - когда амплитудные значения профилей скоростей для задачи с увеличенной вязкостью совпадают с такими же параметрами для задачи с нормальными условиями. Проявляется некоторая конкуренция между «размытием» фронта и изменением амплитуды. В области высоты первоначального опрокидывания волны наблюдается резкое увеличение ДЬ и уменьшение Дгтах. Таким образом, в какой-то момент
Дг, с 2
1,5
0,5
0
« \ » V * \ 1 \
t 1 * 1 * \
\ * 1 % Л
60 70 80 90 100 110 120 130 140
До , м/с
тах
4
60 70 80 90 100 110 120 130 140
у, км
Рис. 3. Зависимость временного масштаба «размытия» фронта ДЬ (а) и величины Дгтах (б) от высоты:
высота у =60 + 140 км, сплошная линия соответствует уровню Ьу(Ь,у) = 0,1 м/с, пунктирная линия — уровню Ьу(Ь,у) = 1 м/с, мощность ш = 103 кт
3
1
мы получаем ситуацию, когда профиль вертикальной составляющей скорости для задачи с повышенной вязкостью начинает распространяться быстрее, имея большую, чем при нормальных условиях, амплитуду.
После образования ударного фронта различия в решениях, полученных при различных значениях вязкости, становятся более существенными, что отражается на амплитудных значениях вертикальной компоненты скорости и на величине размытия фронта. Для задачи с увеличенной вязкостью амплитуда получаемых профилей скорости становится меньше, а ДЬ больше по сравнению с получаемыми решениями для задачи с нормальным условием для вязкости.
При учёте нелинейности и пренебрежении диссипацией постепенное укручение профиля волны приводит к образованию ударного фронта на высоте у0 (9), (10). Для реализуемого алгоритма при одновременном учёте нелинейности и диссипации была получена зависимость величины уо от мощности ш. Для «(¿1) = 15 м/с, ш = 422 кт, уо = 150 — 152 км; для у(Ь1) = 20 м/с, ш = 103 кт, у0 = 125 — 127 км; для = 40 м/с, ш = 8 • 103 кт, уо = 104 — 106 км. Из полученных результатов следует, что высота начала образования ударной волны обратно пропорциональна выделяемой мощности. Чем меньше ш, тем выше необходимо распространиться акустической волне, прежде, чем она начнёт приобретать асимптотически треугольную форму.
Сравнивая получаемые решения в моделировании с аналитическими представлениями для данной области высот, необходимо отдельно обсудить зависимость длительности положительной фазы импульса от высоты. По мере распространения импульса вверх, длительность его положительной фазы монотонно увеличивается - на протяжении высот до момента образования треугольного профиля прирост составляет порядка 1 с на 10 км пройденного возмущением пути. После образования ударной волны прирост становится более значительным. В данном случае имеется два механизма, отвечающих за увеличение длительности импульса. Первый из них связан с тем, что при высотах,
сравнимых с высотами уо, где задаётся поле источника, нелинейность проявляется слабо и приводит она к фазовой модуляции. При дальнейшем подъёме от уровня источника происходит изменение длительности положительной фазы импульса. Второй механизм связан с наличием влияния вязкости на фронт акустического импульса, которое появляется на высотах гораздо меньших, чем высота образования ударной волны. Подобного рода совместное влияние двух механизмов приводит не только к большему увеличению длительности импульса, но и к увеличению высоты опрокидывания волны.
Амплитуда продолжает расти с момента образования ударного фронта вплоть до момента, когда начинает проявляться диссипация. На основании (11), (12) можно выяснить, на какой высоте над источником проявления нелинейности и диссипации становятся равноправными. В области высот, лежащих выше этого значения, можно приближённо провести учёт диссипации в линейном приближении (13). Это приведёт к тому, что под влиянием диссипации, начиная с некоторой высоты, [1] амплитуда акустического импульса начинает спадать с одновременным размытием переднего фронта, обусловленным влиянием вязкости. В получаемом при численном моделировании задачи для мощностей ю = 422 ^ 8 • 103 кт и для конкретной формы начального импульса (рис. 1) решении такой точки спада амплитуды выявить не удаётся. Ситуация такова, что вплоть до 200 + 250 км амплитуда акустического импульса, также как и длительность, продолжают расти.
На высотах свыше 100 + 150 км для стандартной атмосферы (см. CIRA-61 [5]) величина H(у) перестаёт быть постоянной. При этом на высотах, где выполняется приближённое равенство H(у) « const для скорости акустического импульса как функции высоты, наблюдается наличие только одного минимума, что совпадаёт как с теоретическими результатами, так и с результатами, полученными при моделировании. В итоге, не наблюдается спад амплитуды акустического импульса на высотах свыше 150 км. Происходит это по причине того, что для моделирования выполнено условие (4), в частности H(у) = const для любой высоты. Также свой вклад в происходящее вносит совместный учёт слабой нелинейности и потерь. В связи с этим возникает необходимость использования более приближённых к естественным условиям моделей атмосферы, где такие параметры как высота однородной атмосферы, теплопроводность, температура - перестают быть константами и становятся функциями от высоты.
Приложение. Для рассматриваемых импульсов на начальной стадии числа Рей-нольдса велики по сравнению с единицей. Это позволяет учесть нелинейность при пренебрежении диссипацией. На «малых» высотах (Хо ^ у ^ уо; Уо - начало образования ударной волны, Хо - масштаб возмущённой области) вертикальную скорость движения атмосферы опишем линейным приближением
где Н - высота однородной атмосферы, а о = \ZjgH = -у/уРо(0)/ро(0) - линейная скорость звука, у - отношение теплоёмкостей при постоянных давлении и объёме, д - ускорение свободного падения, Ро, ро - невозмущённые значения давления и плотности атмосферы, X - уровень на котором задаётся профиль возмущения
(5)
(6)
Функция exp(y/2H) учитывает неоднородность атмосферы в линейном приближении, а множитель у-1 описывает сферическую расходимость волнового фронта при у ^ уо над эпицентром землетрясения.
Влияние нелинейности на одномерную волну опишем уравнениями газовой динамики
dp | 9{pv) = п.
dt + ду dv dv 1 dP dt ду p ду ^
dP dP
dv
-m+ll—y+ar9d-y
где а - скорость звука. Используя линейное представление для возмущений полей РЛ = Рл — Р0(у), рЛ = рл — ро(у) и линейное представление (5), получим нелинейное представление полей
X
У
■ ехр
У ~ X 2 Н
д№.
дх '
Р'
X
Ро(у) У
exp
У - X
P'
X
р0(у) У
exp
2Н
у -X 2 Н
(/№ | 1 df(
V 2Н ' a0 дт J ' /Дт)(2-у) у df(T)\
V 2Н а0 di J '
где Т = T(y,t), и т(у, t) выбирается так, чтобы на кривых T(y,t) = const выполнялось приближённое характеристическое уравнение
dt dy
1
a0
1
P'
Р'
2ао VPo(y) Ро(у)
(7)
Впервые идея такой нелинейной модификации была сформулирована Л. Д. Ландау [4]. Интегрируя (7) при T(y, t) = const, получим связь между y, t и т:
t - — = (у)
ао
(у-1 )/(т) у-10/ 4Яо0 2оо <9т
+ T(y,t),
(8)
где ф(у) = // 1/у'exp(y7(2Я))в.у'.
Для высот X ^ у ^ у0 имеет место следующее представление: при у « X нелинейность проявляется слабо
ао
А = Хе нф(у)
'^-fft-^ 2а,о dt \ ао
4Нао V ао
и приводит к фазовой модуляции. При дальнейшем подъёме от уровня у = X будет происходить изменение длительностей положительной и отрицательной фаз импульса при постепенном укручении фронта импульса.
v
v
2
а
о
На высоте у = уо начинается образование ударной волны при появлении огибающей семейства линий г = ¿(у, т): дг(у, т)/дт = 0, что даёт соотношение
У + 1д2/ у - 1 /
20д СЙ2
4На0 дт
-Хе знф(у)
Условие первоначального опрокидывания волны имеет вид:
+ 1 = 0.
ф(г/о) а2(У0)
= 2
тах
V
дг2
(у+1)Хе 2н
(9)
(10)
Для высот у0 ^ у ^ ух (ух - уровень, где число Рейнольдса равно единице, т. е. уровень, где существенны потери) загрубляем, соответственно, (8) и (10) к виду:
— -г + х:
ао
У +1
2а0 j
Хе 2нф(у)СМ-ат
ЩУ> 4од сН \сН
31 ат
(11) (12)
Соотношение (12) позволяет найти представление для у(г,у):
, ч X и->- с¥ X У~>• г>и, у) = - е — и ±- е 2Н , У ' У ск у
4а0 / у(г', Х)А'
т
Х(у+ 1)(1?(у)е-ш
Соотношение (11) даёт выражение для длительности положительной и отрицательной фазы импульса вдали от точки уо.
Для области высот у > ух произведём учёт диссипации в линейном приближении
УЛ У-У1
У
дюю(ух, ю)е
V ю*(у)/
(13)
где множитель ух/у ехр[(у — ух)/(2Н)] учитывает изменение поля при увеличении высоты; у(ух, ю) - спектр двуполярной функции у(ух,г) с треугольным профилем огибающей, множитель ехр(—ю2/ю2) приближённо описывает диссипацию монохроматической волны; г' = г — г+(у) — Ух/ах — (у — ух)/ао, а,\ = (dг/dy)-1, а,\ - сверхзвуковая скорость перемещения фронта ударной волны.
Заключение. Подводя итоги рассмотрения результатов численного моделирования эволюции акустического импульса, возбуждаемого землетрясением или мощным взрывом, отметим, что благодаря реализованному алгоритму, удалось исследовать не только области преимущественного влияния нелинейности и потерь, но и область высот равноправного влияния этих факторов. Удалось описать численно решение для областей не только непосредственно над эпицентром землетрясения, но и удалённых от него. Были исследованы основные этапы эволюции акустической волны по мере её распространения над эпицентром возмущения. При этом использовалась одна из простейших моделей атмосферы. Интерес для будущего исследования состоит в том, чтобы, используя проверенную методику расчётов, моделировать процессы распространения акустического импульса в более сложных по своим свойствам средах. Данная работа может служить начальным этапом в такого рода исследованиях.
!
2
2
Литература
1. Павлов В. А. Акустический импульс над эпицентром землетрясения // Геомагнетизм и аэрономия. 1986. Т. 26. № 5. С. 807-815.
2. Безрученко Л. И., Павлов В. А. О возмущении ионосферы над эпицентром землетрясения // Вестн. С.-Петерб. ун-та. Сер. 4: Физика, химия. 1997. Вып. 4. С. 138-141.
3. Броуд Г. Расчёты взрывов на ЭВМ / пер. с англ. М., 1976. 272 с.
4. Ландау Л. Д. Об ударных волнах на далёких расстояниях от места возникновения // Прикл. математика и механика. 1945. Т. 9. Вып. 4. С. 18-22.
5. Справочник по геофизике / пер. с англ. М., 1965. 356 с.
Принято к публикации 8 декабря 2008 г.