Содержание к диссертации
Введение
1 Первый раздел. Математические модели сред, рассматри ваемыхвдиссертации 24
1.1 Двумерные уравнения акустики 25
1.2 Двумерные уравнения линейной теории упругости 27
1.3 Трехмерные уравнения линейной теории упругости 29
1.4 Трехмерные уравнения Эйлера 33
2 Второй раздел. Используемые численные методы 36
2.1 Численные методы решения уравнений акустики и двумерной линейной теории упругости 36
2.2 Численные методы решения трехмерных уравнений линейной теории упругости 41
2.3 Численные методы решения трехмерных уравнений Эйлера 43
2.4 Граничные условия 46
2.4.1 Граничные условия, используемые в двумерном сейсмическом решателе 46
2.4.2 Граничные условия,используемые в трехмерном сейсмическом решателе 47
2.4.3 Граничные условия, используемые в газодинамическом решателе 48
3 Третий раздел. Результаты численного моделирования Челябинского события 50
3.1 Поверхностные волны 50
3.1.1 Волны Рэлея 50
3.1.2 Волны Лява 52
3.2 Моделирование сейсмического отклика на основания возмущения, вызванного инфразвуковыми волнами 54
3.2.1 Постановка задачи 54
3.2.2 Анализ результатов и сравнение с реальными данными 56
3.3 Моделирование сейсмического отклика, вызванного головной ударной волной 59
3.3.1 Постановка газодинамической части 59
3.3.2 Постановка сейсмической части 59
3.3.3 Возможные сценарии удара Челябинского метеорита 61
3.3.4 Анализ результатов и сравнение с реальными данными 63
Заключение 71
Список литературы
- Двумерные уравнения линейной теории упругости
- Трехмерные уравнения линейной теории упругости
- Граничные условия, используемые в двумерном сейсмическом решателе
- Моделирование сейсмического отклика на основания возмущения, вызванного инфразвуковыми волнами
Введение к работе
Актуальность и степень разработанности темы. Проблема астероидно-кометной угрозы становится все более актуальной. Так, с увеличением численности населения Земли возрастает вероятность фатальных последствий в результате столкновения с достаточно крупным метеоритом. Ученые-астрономы сегодня могут только относительно достоверно определить, когда и какие космические тела войдут в атмосферу Земли. Столкновение с метеоритом в районе г.Челябинск наглядно продемонстрировало уровень опасности такого явления природы. Следует заметить, что несмотря на относительно крупный размер (около 15-17 м в поперечнике), он не был обнаружен заблаговременно. В целом, класс объектов размером 10-300 м представляет значительную угрозу, поскольку на данный момент регистрация таких объектов производится хуже всего, а разрушительный эффект от таких объектов может быть весьма значительным. В результате взаимодействия с атмосферой Земли такие объекты обычно лавинообразно разрушаются, а затем испаряются с образованием ниспадающего к поверхности Земли высокоскоростного расширяющегося потока газа. В результате столкновения с поверхностью Земли образуется ударная волна, а точнее – серия ударных волн, распространяющихся на большие расстояния. Как следствие, характерный размер области взаимодействия превышает характерный размер метеорита на порядки. При этом основной эффект от столкновения обусловлен воздействием головной ударной волны от взрыва метеорита в атмосфере. По последним оценкам, частота прилета космических объектов, которые могут нести существенную угрозу, превышает раз в сто лет. Когда подобные объекты обнаруживаются слишком поздно, не возможно существующими средствами воспрепятствовать их столкновению с Землей. В этом случае необходимо иметь достаточно полную картину вероятных последствий разрушений от возможного столкновения космического тела с Землей с учетом как атмосферного, так и сейсмического эффектов. Эффективное численное моделирование возможных последствий позволит подготовить наземные и подземные объекты с целью минимизации разрушений (например, метро, АЭС, плотины гидроэлектростанций). Данная проблема особенно актуальна для России, занимающей обширную
территорию в 1/7 суши всей Земли. Оперативное численное моделирование позволит избежать применения дорогостоящих мер борьбы с космическим телом в случае, когда зона поражения по результатам моделирования придется на ненаселенный район. Еще одним важным аспектом применения оперативного численного моделирования снижение вероятности ядерного конфликта в случае столкновения с метеоритом. Так, известно, что Южная Корея изначально интерпретировала сейсмические данные, полученные от взрыва Челябинского метеорита, как вероятный ядерный удар со стороны Северной Кореи. Сегодня вопросы защиты от астероидно-кометной угрозы ставятся на самом высоком уровне. В подтверждение этому в конце 2013 года в журнале Scientifc American вышла статья United Nations to Adopt Asteroid Defense Plan (ООН принимает план по защите от астероидов). В ней подчеркивается, что несмотря на высокую степень опасности «ни одно государство в мире сегодня не поставило своими ведомствам конкретную задачу по защите планеты от астероидов». Толчком в изучении этой проблемы послужило падение Челябинского метеорита 15 февраля 2013 года (1613 пострадавших). Оценка размера Челябинского метеорита дала величину 15-17 м в поперечнике. По последним оценкам размер Тунгусского метеорита мог достигать 25 метров и более. В свою очередь, частота прилета на Землю космических объектов диаметром 30 метров оценивается как один раз в сто лет, для меньших объектов частота прилета значительно выше. Научная значимость проблемы предупреждения астероидно-кометной угрозы подтверждается проведением многочисленных конференций по данной тематике в последнее десятилетие. В частности, следует отметить конференции “Planetary Defense Conference” 2004, 2007, 2009, 2011 и 2013 годов, проходившие в США, Испании и Румынии; конференция, посвященная Тунгусскому метеориту в Москве в 2008 году; две конференции, посвященные астероидно-кометной угрозе в Санкт-Петербурге в 2005 и 2009 годах.
На данный момент, рассматривая вопрос входа космического тела в атмосферу Земли, большинство авторов выделяют и исследуют какое-то одно характерное явлений, касающееся этой задачи.Такой подход позволяет более тщательно рассмотреть выбранное явление, но не дает полной картины о всем явлении входа космического тела в атмосферу с последующим вза-
имодействием с поверхностью Земли. Существуют математические методы моделирования, позволяющие достаточно полно описывать самые различные процессы, связанные с вхождением метеорита и его последующим влиянием. Но стоит отметить, что предложенные подходы носят оценочный характер и построены, например, на обобщении наблюдений. Поэтому разработка комплексного метода моделирования поверхностных волн, вызванных влиянием движущегося из газовой среды объекта на упругую среду в рамках трехмерной постановки через сопряженное решение задачи газ-твердое тело позволит избежать оценочных предположений и моделировать явления, попадающие в рассматриваемый класс, с высокой точностью. В частности, явления входа космического тела в атмосферу Земли.
Цели и задачи работы Целью работы является разработка комплексного метода моделирования поверхностных волн, вызванных влиянием движущегося из газовой среды объекта на упругую среду в рамках двухмерной и трехмерной постановки через сопряженное решение задачи газ-твердое тело, и алгоритма анализа реальных данных наблюдений на основании их математической модели. Применительно к явлению столкновения метеорита с Землей, разработанный метод моделирования позволит лучше верифицировать характеристики входа метеорита в атмосферу Земли, а разработанный алгоритм интерпретации даст механизм распознавания природы удара одновременно как на основании атмосферных, так и сейсмических данных. Задачи работы были следующие:
-
Разработка метода моделирования поверхностных сейсмических волн, вызванных акустическим инфразвуковым возмущением через решение сопряженной задачи газ-твердое тело в рамках двухмерной постановки. Исследование возникновения и эволюции низкочастотной сейсмической волны в приповерхностном слое земной коры, инициированной инфразвуковой волной от взрыва;
-
Разработка метода моделирования поверхностных волн, вызванных ударной волной через решение сопряженной задачи газ-твердое тело в рамках трехмерной постановки. Создание трехмерной модели, имитирующей столкновение метеорита с Землей, учитывающей как
газодинамическую, так и сейсмическую часть на основе современных численных методов;
-
Комплексное исследование влияния головной ударной волны метеорита на поверхность Земли;
-
Разработка алгоритма интерпретации явления входа тела в газовую среду на основании анализа реальных сейсмических данных. Верификация и уточнение параметров входа Челябинского метеорита путем сравнения с реальными данными;
-
Разработка двухмерного программного комплекса для моделирования акустических и сейсмических явлений.
Методы исследования. Процесс прогрессивного разрушения твердого тела при входе в газовую среду на большой скорости оценивается с помощью подхода, разработанного в работах Тирского Г. А. и Ханукаевой Д. Ю. Для моделирования газодинамической части задачи используется численный метод с применением технологии динамических подвижных сеток, разработанный в работах Утюжникова С. В. и Руденко Д. В. Сейсмическая задача рассчитывается с помощью сеточно-характеристического метода, разработанного в работах Холодова А. С. и Петрова И. Б. Программный комплекс для решения задачи сейсмики написан для реализации счета на вычислительных системах с параллельной архитектурой. Разработанные модели и результаты вычислений апробированы на решении тестовых задач и путем сравнения с данными наблюдений.
Научная новизна. Немногочисленные существующие модели исследуемого в диссертации явления в основном либо слишком приближенные (носят оценочные характер, основаны на обобщении наблюдений), либо строятся при упрощающих предположениях (источник моделируется точечным взрывом, не учитывается сейсмический эффект). В работе разработан метод моделирования и 3D-модель (с проведением предварительных оценок на основании 2D-модели) сразу для всей задачи (трехмерное моделирование атмосферных и сейсмических эффектов). Применяется модель взрыва, достаточно подробно учитывающая все происходящие при этом явления.
Сопряжение задач газ-твердое тело позволяет задать источник сейсмических волн физически правдоподобно, не прибегая к каким-либо приблизительным предположениям и оценкам. Численное решение задачи выполнено на динамически подвижных сетках с использованием параллельных вычислений. Стоит отметить, что трехмерное моделирование осложняется именно большим объемом вычислений. В частности, мировой лидер в области моделирования удара космических тел и его последствий - Sandia National Laboratory (USA) – только приступает к трехмерному моделированию, поскольку они проводят расчеты на стационарной сетке, охватывающей всю атмосферу Земли. Использование динамически подвижных сеток значительно ускоряет процесс вычисления, поскольку допускает локализацию моделируемого явления.
Теоретическая и практическая значимость. В рамках работы создана модель столкновения метеорита с Землей и прогноза его сейсмического эффекта на поверхности планеты. Задача рассмотрена на примере падения Челябинского метеорита. Челябинский метеорит относится к классу космических тел кометной природы, которые разрушаются и полностью или почти полностью испаряются в атмосфере Земли. К таким телам также относится хорошо известный Тунгусский метеорит. Основное воздействие от таких космических тел осуществляется посредством образованной после взрыва головной ударной волны. Выбор именно этого класса объектов связан со следующими причинами: такие объекты прилетают на Землю наиболее часто среди всех космических тел, они хуже всего регистрируются (могут быть обнаружены слишком поздно), разрушительный эффект от таких объектов может быть весьма значительным. Значимость работы, в том числе и теоретическая, состоит в создании инструмента для распознавания природы удара. Сравнение сейсмических результатов, полученных на выходе численного эксперимента, с данными, записанными сейсмографами, позволит интерпретаторам однозначно выделять их в отдельный класс сейсмических событий, вызванных именно влиянием головной ударной волны метеорита (а, например, не взрывом атомной бомбы). Также предложенный подход позволяет уточнить параметры входа в атмосферу космических тел на основании анализа сейсмических данных.
Положения, выносимые на защиту Положения, выносимые на защиту, совпадают с основными результатами диссертации, приведенными в конце автореферата.
Степень достоверности и апробация работы. Высокая степень достоверности результатов работы подтверждена корректностью использованного математического аппарата и сравнениями с данными реальных наблюдений. Результаты диссертации доложены, обсуждены и получили одобрение специалистов на следующих конференциях и научных семинарах:
54-я, 55-я, 56-я, 57-я, 58-я научные конференции МФТИ (Москва, 2011-2015);
Международная конференция "Потоки и структуры в жидкостях"(Санкт-Петербург, 2013);
Семинары кафедры информатики и вычислительной математики МФТИ (Москва, 2013-2016);
Семинары лаборатории математического моделирования нелинейных процессов в газовых средах МФТИ (Москва, 2011-2015).
Публикации по теме диссертации. Результаты диссертации опубликованы в девяти статьях, две из которых [1 2], изданы в журналах, рекомендованных ВАК, одна [2] в издании, входящем в Scopus.
Личный вклад соискателя в работы с соавторами. Все результаты диссертации и соответствующие им части публикаций получены лично соискателем. Реализован программный код для решения двумерных задач акустики и сейсмики [3, 4]. В [1, 6] на основании предложенного в диссертации метода моделирования, разработана модель и проведены исследования возникновения и эволюции низкочастотной сейсмической волны в приповерхностном слое земной коры, инициированной инфразвуковой волной от взрыва (двумерные расчеты). В [2, 5, 7, 8, 9] разработан метод моделирования поверхностных волн, вызванных ударной волной и предложена
трехмерная модель, имитирующая столкновение метеорита с Землей, учитывающая как газодинамическую, так и сейсмическую часть. Проведена серия газодинамических расчетов. Проведена серия сейсмических расчетов. Так же в [2] проведено исследование влияния головной ударной волны метеорита на приповерхностный слой Земли, выполнена верификация и уточнение параметров входа Челябинского метеорита в атмосферу путем сравнения с реальными данными.
Структура и объем работы. Диссертация состоит из введения, трёх разделов, заключения, списка литературы. Полный объем диссертации - 82 страницы текста с 19 рисунками и 1 таблицей. Список литературы содержит 99 наименований.
Двумерные уравнения линейной теории упругости
Поскольку основное исследование диссертации строится на моделировании Челябинского события, отдельное внимание уделим работам, посвященным изучению этого явления.
К настоящему времени помимо первых и наиболее полных работ по Челябинскому событию [18, 17], вышел ряд работ отечественных и зарубежных авторов, исследующих отдельные части этого сложного с точки зрения физики явления [44, 45, 46, 47, 48, 49, 50, 51]. В работе [44] производится реконструкция орбиты Челябинского метеорита до столкновения с Землей, что помогает классифицировать такого рода объекты, что в свою очередь, способствует проведению оценки частоты входа объектов подобного класса в атмосферу Земли. В [45] авторы проводят достаточно подробное исследование состава метеорита. В статье [46] моделируется волна Рэлея на расстоянии 4000 км от источника. На основании информации из анализа сейсмических диаграмм, реконструкции траектории входа из видеонаблюдений, производится оценка магнитуды сейсмического источника, которая по заключению авторов составляет 3.7. На основании теоретической модели в [47] проводится исследование возмущения геомагнитных полей, вызванное прохождением метеорита через ионосферу Земли. Показано, что полученные результаты моделирования хорошо согласуются с наблюдаемыми данными. В работе [48] проводится исследование инфра-звуковых возмущений, зафиксированных широкой сетью станций, и оценивается их энергия. Оценка параметров движения Челябинского метеорита на основе простых моделей проводится в [49]. В используемых моделях учитывается абляция и механическое дробление метеорита при входе. Несмотря на достаточно подробный анализ, полученные в работе результаты следует считать весьма приближенными по причине недостаточно подробного учета механизма дробления. В работе [50] представлено численное моделирование теплового разрушения Челябинского метеорита при входе в атмосферу Земли в рамках подхода, включающего расчет траектории движения с учетом сопутствующих движению метеорита физических процессов (лучисто-конвективный теплообмен, вдув продуктов разрушения во встречный поток газа, прогрев, разрушение). Определялось поле течения. Решение задачи позволяет точнее определить траекторию движения метеорита от точки входа в атмосферу Земли до полного разрушения. В работе [51] проводится моделирование начальной стадии разрушения метеорита в плотных слоя атмосферы. Для решения задачи используется упругопла-стическое приближение. В серии расчетов варьируются характеристики материала метеорита и угол входа. Моделирование показывает, что падение скоростного напора является следствием разрушения метеорита и что динамика скоростного напора зависит не только от аэродинамических, но и прочностных свойств материала, составляющего метеорит.
Основным подходом решения системы уравнений акустики и линейной теории упругости в моделях, рассматриваемых в диссертации, является метод расщепления многомерной системы уравнений в частных производных гиперболического типа на несколько одномерных. Для решения одномерных систем на основании их гиперболичности используется переход к инвариантам Римана, благодаря которому одномерная система распадается на несколько (пять для двухмерной задачи и девять для трехмерной) независимых скалярных уравнений переноса, решаемых конечно-разностными схемами. Рассмотрим иные подходы, которые используются для решения систем уравнений акустики и линейной теории упругости. Прежде чем об суждать сами численные методы, рассмотрим некоторые типы сеток, на которых эти численные методы определяются. По способу заданий сеточных функций в узлах сетки можно выделить три типа сеток (рис. 1): 1. обычная сетка (conventional) - сетка, на которой все сеточные функции определены в каждом сеточном узле; 2. частично разделенная сетка (partly-staggered) - сетка, на которой смещения среды или скорости смещения среды определены в одних узлах, а компоненты тензора напряжения - в других; 3. разделенная сетка (staggered) - сетка, на которой смещения среды или скорости смещения среды, сдвиговые компоненты тензора напряжения и нормальные компоненты тензора напряжения определены на разных сеточных узлах;
Также различают системы уравнений упругости в зависимости от входящих в них физических переменных. Уравнения могут включать смешения среды, скорости смещения среды (или то и другое вместе) и компоненты тензора напряжения. В зависимости от типов используемых сеток и решаемых на них уравнениях могут быть использованы различные численные методы. В самом начале развития численных методов для задач сейсмики широкое распространение получили численные методы на обычных сетка в силу более естественной аппроксимации уравнений [52, 53, 54, 55, 56, 57]. Но такие численные методы могут быть неустойчивы для моделей сред с сильным различием в скоростях распространения волн. С этой проблемой могут эффективно бороться методы на разделенных сетках [58, 59]. К тому же их использование требует меньшего использования памяти при численной реализации, нежели схем на обычных сетках [60] и может быть получен более высокий порядок аппроксимации, чем при использовании метода расщепления [61, 62, 63, 64, 65]. Недостатком схем на разделенных сетках является их неудовлетворительное использование для задач с анизотропными средами [66]. В этом случае большую применимость показывают схемы на частично-разделенных сетках [67, 68]. В целом, сравнивая численные методы для решения задач сейсмики на обычных сетках и на разделенных сетках (или частично разделенных), стоит отметить, что несмотря на все описанных преимущества последних и высокий порядок предлагаемых численных методов,часто оказывается, что оператор метода сложно имплементировать на границу, и тем самым высокий порядок аппроксимации снижается за счет ошибок более низкого порядка с границ [69, 70]. Подробное исследование вопроса постановки и имплементации граничных условий для двумерных задач производится в работе [71]. Поэтому при выборе численного метода стоит также учитывать простоту его имплементации в программный код.
Трехмерные уравнения линейной теории упругости
Для решения системы уравнений акустики и двумерной линейной теории упругости (1) применяется метод расщепления. Вместо системы уравнений (1) решаются одномерные задачи вида: d\J d\J о—1 Aj— = 0, (7) at axi Пусть L - двумерный оператор перехода от решения на временном слое п к решению на временном слое п + 1: Un+ = L(Ai,A2)\Jn.
Пусть Li - одномерный оператор перехода от решения на временном слое п к решению на временном слое п + 1, соответствующий одномерной системе с матрицей АІ. Фактически, оператором Li определяется разностная схема, применяемая для решению конкретной одномерной системы уравнений (7): Un+1 = Li(Ai)\Jn. В двумерном случае выражение оператора L через операторы L\ и L2 имеет простой вид: L{A\)A2) = 2( 2) 1( 1). Такой подход обеспечивает второй порядок точности решению (при выборе численного метода второго порядка точности для решения одномерной системы).
Для решения одномерных систем (7) используется переход к инвариантам Римана. Поскольку система уравнений (7) является системой уравнений в частных производных гиперболического типа, то матрица системы Aj обладает полным набором собственных векторов с действительными собственными числами. Тогда эту систему можно представить в виде: После перехода к инвариантам Римана полученные независимые одномерные уравнения переноса решаются сеточно-характеристический методом.
Для решения одномерных уравнений переноса удобно использовать сеточно-характеристические схемы [16, 76]. Схемы такого типа хорошо учитывают природу уравнений гиперболического типа. После перехода к инвариантам Римана каждому независимому уравнению переноса соответствует своя характеристическая прямая, вдоль которой сохраняется решение. В этом случае решение уравнения переноса оказывается достаточно простым. Из узла на + 1 временном слое, на котором требуется получить решение, опускается характеристика до пересечения с временным слоем , на котором решение уже известно. Если характеристика пересекает слой между узлами расчетной сетки, то проводится интерполяция значения там по соседним узлам сетки. Поскольку решение на характеристике сохраняется, то решение в узле на + 1 слое будет определяться по полученному путем интерполяции на слое значению. Качество проведенной интерполяции будет определять точность численного метода. Схематически подход, используемый в сеточно-характеристическом методе, изображен на рисунке 2.
Выпишем для уравнения (8) гибридную монотонную сеточно-характеристическую схему второго порядка, основанную на комбинации центрально-разностной схемы и схемы с ориентированными разностями. Представим ее в потоковой форме п+1 и = ипт - г(/ ,1 - /m-i), і где потоки имеют вид: J m-\-i ам г-1 + (1 - а - /3)г + /Зг +1 , а 0 агС+2 + (1 - « - /3)г +1 + /Зг , а J m-h &и т-2 + (1 - а - /3)м г-1 + /Зг , а 0 aum+i + (1 - « - (3)uvm + /Зм г-1 , а 0 r - число Куранта. В такой записи получается целое семейство схем, зависящих от параметров а и /3. Требование минимума аппроксимационной вязкости приводит к условию для класса схем с ориентированными разностями (/3 = 0):
Из этих соотношений для каждого класса схем (с заданными таким образом и ) можно составить условия монотонности схем. Для класса схем с ориентированными разностями условие монотонности примет следующий m+i = 0.5(1 + )С, + 0.5(1 — ) +1, m_i = 0.5(1 + ) г_1 + 0.5(1 — }ат. Выписанная явная монотонная гибридная схема будет обладать вторым порядком аппроксимации как по пространству, так и по времени (для гладких решений). Условие устойчивости схемы 0 Р 1. 2.2 Численные методы решения трехмерных уравнений линейной теории упругости
Также как и для уравнений двумерной линейной теории упругости, чтобы решить систему (2), используется метод расщепления. Вместо системы уравнений (2) решаются одномерные задачи вида (7).
При определенном выборе коэффициентов использование метода расщепления обеспечивает второй порядок точности решению. Классическим в этом смысле можно назвать следующий подход: 1 —v L(Ai,А2,A3) = - у Li{Ai)Lj{Aj)Lk{Ak). Но поскольку такой подход требует 18 применений операторов Li и хранения данных промежуточных шагов для них, он является весьма затратным с точки зрения программной реализацией.
Альтернативным здесь является применение подхода, предложенного в [77]. Чтобы сохранить второй порядок точности и при этом обеспечить простоту применения самого метода, используется предложенная в работе комбинация:
Для решения одномерных систем (7) используется переход к инвариантам Римана. После перехода к инвариантам Римана, система уравнений (7) распадается на девять независимых скалярных уравнений переноса. Полученные независимые одномерные уравнения переноса решаются с помощью TVD-схемы с ограничителем МС (monotonized central) [78]. Такая схема обладает вторым порядком точности. Выпишем TVD-схему для уравнения переноса (8) в виде: и7+ = uvm — r(/m+i — /m_i), 2 2 где потоки / имеют вид: Jm+l = U \ фшт){1 — Г) [U + і — U ), 2 2 JTO_i = мто Н—фшт){1 — г){ит — um_i)} 2 2 r - число Куранта, ф - ограничитель потока. Параметр qm определяется как Qm = В расчетах для TVD - схемы использовался ограничитель MC: 1 фш) = тах\0, min(2q} -(1 + q), 2)1, lim 0(g) = 2. Анализ такой TVD-схемы для гладкого, треугольного и прямоугольного сигналов показывает его хорошую применим как по качеству аппроксимации (точность, монотонность) так и по скорости счета (в сравнении с другими ограничителями).
Граничные условия, используемые в двумерном сейсмическом решателе
В результате расчета была получена серия синтетических сейсмограмм, записанных на поверхности. Проведем сравнительный анализ сигналов и спектров сигналов, полученных по сейсмической трассе на расстояние 1600 км с реальными данными, полученными на сейсмической станции. На рисунке 9 представлено сравнение реального и модельного сигнала. Стоит заметить, что ноль по оси абсцисс - абсолютное время от момент инициации источника.
Сравнивая сигналы чисто визуально, можно видеть их достаточно хорошее совпадение, что подтверждает, что выбранное моделирование воспроизводит реальную природу сигнала, импульс Берлаге хорошо подходит для моделирования инфразвуковых волн, вызванных взрывом метеорита. Это означает, что такой импульс может быть использован и при моделировании схожих явлений. Теперь сравним спектры сигналов. На рисунке 10 представлены спектры реального и модельного сигнала.
На этом рисунке видно, что формы спектров и доминантная частота (0.03 Гц) обоих сигналов фактически совпадают. Учитывая, что спектр ин-56 фразвуковой волны задавался по фактически зарегистрированному диапазону [18], а спектр синтетической сейсмической волны соответствует спектру реальной волны, это также подтверждает пригодность выбранной модели для изучения как непосредственно Челябинского события, так и допустимость к исследованию явления подобного класса.
В результате исследования сейсмического эффекта в приближении модели точечного взрыва метеорита, вычислена динамика распространения волн как в воздухе, так и в приповерхностном слое Земли. Сравнительный анализ результатов показывает, что даже такая грубая, по современным меркам, модель (двумерное моделирование, точечный источник модели взрыва) способна дать физичный результат, что подтверждается близостью форм и спектров модельного сигнала и реального сигнала, зарегистрированного на сейсмостанции в г. Обнинск.
Безусловно, подробные трехмерные расчеты с более физичной моделью взрыва метеорита способны обеспечить более точные результаты (что и будет показано далее). Но обычно такие расчеты требуют значительных
Спектр модельного (1) и реального (2) сигналов. числительных ресурсов, что увеличивает время счета (В подобных расчетах, проводимых автором работы, это время идет на дни для одного расчета на кластере). И может оказаться, что проведенный расчет дает неверный результат в силу недостаточно точного задания входных данных. В свою очередь модель, рассматриваемая в работе требует значительно меньших мощностей (несколько часов расчетов на ПК). Поэтому, поскольку она достоверно воспроизводит физическую природу, такая модель может быть эффективно использована для предварительного конфигурирования входных данных более точных моделей исследования сейсмических эффектов, вызванных взрывом крупного космического тела в атмосфере Земли. 3.3 Моделирование сейсмического отклика, вызванного головной ударной волной
В рамках газодинамической модели используется постановка, предложенная в [97, 98]. Задача решается в стратифицированной атмосфере, в которой плотность и давление определяются из следующих соотношений:
На высоте взрыва помещается объем, содержащий однородную газовую смесь со следующими параметрами: масса газа полагается равной массе космического тела; скорость газа скорости входа в атмосферу; плотность газа соответствует плотности вещества космического тела; статическое давление - давлению торможения на данной высоте. Отношение внутренней энергии вещества космического тела к его кинетической энергии в рамках данной постановки однозначно определяется отношением плотности атмосферы на высоте взрыва к плотности вещества тела. Приемлемость такой модели для трехмерной задачи распространения ударных волн от взрыва метеорита вплоть до взаимодействия с поверхностью Земли подтверждается расчетами, произведенными в [14].
Физическая область модели земной поверхности представляет собой параллелепипед. Физика распространения сейсмических волн в этой области определяется уравнениями линейной динамической теории упругости. На верхней границе физической области, соответствующей поверхности Земли, было определено граничное условие “свободная поверхность”, на остальных же границах заданы поглощающие граничные условия. Поскольку сей смическая станция, на которой был записан реальный сигнал и с данными которой производится сравнение, находится на расстоянии от Челябинска почти 1600 км, то размеры физической области были выбраны следующие: 2000 км в длину, 800 км в ширину и 150 км в глубину. Такой существенный запас в размерах обоснован тем, что численная реализация поглощающих граничных условий не идеальна, и некоторые отраженные волны все равно возвращаются в расчетную область, а данный выбор не позволяет отраженным волнам влиять на основное возмущение (рис. 12). Сейсмический источник (возмущение от головной ударной волны) был помещен на поверхности в равном удалении на 400 км от трех боковых границ области и на 1600 км от четвертой боковой границы. Для задания геологических свойств среды были приняты следующие предположения. На основании реальных данных и из оценки средней скорости поверхностной волны Рэ-лея [18],инициированной головной ударной волной Челябинского метеорита, среднее значение длины волны составило примерно 100 км. Для такой достаточно большой длины волны, учитывая структуру земной коры Русской равнины и Уральских гор (рис. 11), основное влияние оказывают более глубокие, базальтовые и гранитные слои, несмотря на различие в структуре приповерхностного слоя [99]. На основании этого на всей протяженности распространения волны от источника до момента регистрации среда рассматривается как однородная с наиболее приемлемыми средневзвешенными значениями продольной и поперечной скорости распространения волн = 7.18 км/с и = 3.4 км/c [88]. Рассмотрим вопрос передачи данных из газодинамической части в сейсмическую. Передача данных осуществлялась на основании прямого механизма, то есть избыточное давление из распределения на поверхности Земли в газодинамическом расчете передавалось на прямую как функция, определяющая упругие деформации в сейсмической модели, и формировало сейсмический источник [91]. При передачи данных на области 40 км 40 км производился переход с более подробной газодинамической сетки (равномерная квадратная на поверхности с шагом 400 м) на менее подробную сейсмическую (равномерная квадратная на поверхности с шагом 4 км).
Моделирование сейсмического отклика на основания возмущения, вызванного инфразвуковыми волнами
Перед основным исследованием было проведено предварительное изучение зависимости параметров сейсмического сигнала (амплитуды и спектра), получаемого в модели от параметров входа. Основными параметрами входа, определяющими конфигурацию явления, надо считать высоту взрыва, скорость входа, угол входа, плотность метеоритного вещества и характерный диаметр. Стоит отметить, что такие параметры входа, как скорость и угол известны достаточно точно из данных наблюдений, и небольшие их изменения не приводят к чувствительному изменению параметров сигна Рис. 12: Распространение волны Рэлея в трехмерном вычислительном домене. Провизуализирована нормальная компонента скорости смещения среды.. Поэтому для варьирования в исследовании были выбраны параметры входа с большей неопределенностью в значении: высота взрыва, плотность и диаметр метеорита. Было получено, что при фиксированных значениях скорости и угла входа, к уменьшению амплитуды сигнала приводит по отдельности уменьшение плотности, уменьшение диаметра и увеличение высоты взрыва. Так же стоит отметить, что варьирование этих трех параметров не приводит к изменению доминантной частоты сигнала. На основание оценок плотности, диаметра и высоты взрыва была проведена серия расчетов с различными значениями этих параметров. В качестве основного результата здесь будут представлены и проанализированы три варианта (см. Таблицу 1). Поясним, почему представлены именно эти три варианта. Вариант 1 наиболее точно воспроизводит реальный сигнал среди всей серии тестовых расчетов. Вариант 2 дает сейсмический отклик, который получается при выборе значений параметров, взятых в качестве наиболее вероятных в согласии с [18, 17]. Вариант 3 дает минимальный сейсмический отклик при выборе высоты взрыва, плотности и скорости из допустимых диапазонов (оценка снизу). Оценка сверху здесь не приводится, поскольку Номер вариан- Высота взрыва Плотность кос- Характерный та h, км мического вещества, г/см3 диаметр, м даже Вариант 1 немного завышает амплитуду сигнала по сравнению с реальным сигналом. Для всех вариантов: скорость входа в атмосферу V = 18.6 км/c, угол входа 18 град. Скорость и угол входа соответствуют значениям, представленным в [18, 17]. Так же там дается оценка высоты взрыва 23 - 30 км, в качестве оценки плотности (для хондрита) можно выбрать 3-3.7 г/см3 (в зависимости от состава), в качестве характерного диаметра -16-19 м. На основании исследований, проведенных в [18, 17] (анализ свечения, состава фрагментов метеорного вещества) Вариант 2 можно считать основным (с наиболее вероятными параметрами).
Сравним результаты моделирования для каждого варианта с реальными данными, полученными на сейсмостанции в г. Обнинск, Московская обл.. На рис. 13 слева представлен отфильтрованный сигнал, зарегистрированный сейсмодатчиками. Стоит обратить внимание на структуру сигнала при t 490 cиt 560 с. Для этих моментов времени возмущение определяется уже не воздействием головной ударной волны метеорита, а маломощным землетрясением, поэтому при сравнении сигналов и спектральном анализе не стоит учитывать эти части структуры сигнала. Справа на рисунке представлен спектр реального сигнала. Пик спектральной картины приходится на частоту примерно равную 0.03 Гц.
Для начала сравним сигналы, полученные по трем модельным вариантам, с реальным по форме и амплитуде. Видно (рис. 14), что форма сигнала во всех случаях одна и та же, то есть можно полагать, что принципы, заложенные в реализацию взрыва метеорита и последующего распространения ударной волны и ее влияния на поверхность Земли отражают Рис. 13: Реальный сигнал, записанный на сейсмостанцие (слева) и его спектр (справа) реальную природу явления. Вариант 3 отражает минимально возможное воздействие, которое могла бы оказать головная ударная волна. Заметим, что для него выбраны минимальные значения из диапазона оценок для плотности и диаметра метеорита, а также максимальная оценка высоты взрыва метеорита. Наиболее близки к реальному сигналу по амплитуде варианты 1 и 2, причем вариант 2 – это те параметры входа, которые принято считать наиболее приемлемыми значениями оценки [17]. В варианте 1 значение высоты взрыва увеличено до верхней границы оценки высоты (30 км).
Для наглядности сравним Варианты 1 и 2 с реальным сигналом по отдельности. Из рис. 15 видно, что Вариант 2 наиболее хорошо воспроизводит первый пик сигнала, но при этом имеет ощутимо выделяющийся второй пик (минимум). Вариант 1 (рис. 16)лучше воспроизводит второй и третий пики. Для выбора более предпочтительного варианта проведем спектральный анализ модельных сигналов на основании быстрого преобразования Фурье и сравним их с реальным сигналом (рис. 17). На основании графика можно заключить, что все три варианта имеют одну и ту же доминантную частоту, которая практически совпадает с доминантной частотой реального сигнала. Анализируя амплитуды спектральных картин, можно заключить, что Вариант 1 более точно воспроизводит реальный сигнал, поскольку спектральная амплитуда Варианта 1 оказывается Рис. 14: Реальный сигнал и модельные сигналы трех вариантов достаточно близкой к спектральной амплитуде реального сигнала. Спектральная амплитуда Варианта 2 оказывается существенно превышающей спектральную амплитуду как реального сигнала, так и Варианта 1. Тем не менее, поскольку спектральная амплитуда Варианта 1 все же немного превышает спектральную амплитуду реального сигнала, то для коррекции параметров входа необходимо уточнить некоторые из трех параметров входа (плотность, диаметр, высота взрыва). Поскольку для уменьшения амплитуды необходимо увеличение высоты взрыва, а высота в Варианте 1 выбрана как верхняя граница соответствующей оценки, этот параметр оставлен без изменений. Также большинство исследователей (в том числе и авторы статей [18, 17]) склонны считать, что метеорит состоял из однородного вещества с плотностью, точно определенной на основании исследования фрагментов метеорита [17]. Поэтому единственный параметр за счет которого можно провести коррекцию в рамках Варианта 1 – уменьшение диаметра метеорита (выбрать диаметр меньше 18 м. Тем не менее такое уменьшение должно быть несущественным, менее метра. Этот вывод можно сделать на основании сравнения со спектром Варианта 3).
Существуют некоторые исследователи, считающих что значительная часть метеорита состояла изо льда (по причине того, что было найдено мало фрагментов метеорита). Как видно из Варианта 3, уменьшение средней плотности при незначительном уменьшении диаметра приводят к значительному уменьшению амплитуды сейсмического сигнала. В случае же структуры лед-хондрит, значение плотности станет существенно меньше плотности Варианта 3 (в зависимости от того, конечно, какую часть метеорита считать ледяной), что приведет к амплитуде сигнала значительно меньшей, чем у реального. Поэтому, на основании сейсмического анализа можно прийти к заключению, что такая структура является невозможной.