Компьютерное моделирование пространственных динамических процессов в сложных областях интегрирования с использованием наложенных сеток тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Стецюк Владислав

• Стецюк Владислав
• кандидат наук
• 2023

Оглавление диссертации кандидат наук Стецюк Владислав

Введение

Глава 1. Определяющие уравнения

1.1 Введение

1.2 Уравнения теории упругости

1.3 Теория упругости в изотропной среде

1.4 Параметры среды

1.5 Теория акустики

Глава 2. Математическая модель

2.1 Введение

2.2 Сеточно-характеристический метод

2.2.1 Матричные уравнения линейной теории упругости в декартовых координатах

2.2.2 Матричные уравнения линейной теории упругости в криволинейных координатах

2.2.3 Уравнения теории акустики в матричном виде

2.2.4 Граничные условия

2.3 Методы описания топографии

2.4 Методы интерполяции

2.4.1 Метод ближайшего соседа

2.4.2 Метод обратных расстояний

2.4.3 Метод Сибсона

2.4.4 Барицентрические координаты

2.4.5 Метод параметрического представления

2.4.6 Метод локального разложения

2.4.7 Метод глобального полинома

2.4.8 Кригинг

Глава 3. Программный комплекс

3.1 Введение

3.2 Интерполяция как предобработка

3.3 Интерполяция наложенных сеток

Стр.

3.3.1 Процесс интерполяции наложенных сеток

3.3.2 Проверка вхождения точки в сетку

3.3.3 Разбиение пространства на ячейки

3.3.4 Коллинеарность и компланарность

3.4 Реализации методов интерполяции

3.4.1 Методы ближайших соседей и обратных расстояний

3.4.2 Барицентрические координаты и параметрическое представление

3.4.3 Метод локального разложения

3.4.4 Метод естественной окрестности

3.4.5 Зависимость погрешности от количества точек

3.5 Многопоточная реализация интерполяции

3.6 Генерирование реалистичных синтетических карт высот

3.7 Построение сеток для описания поверхности

Глава 4. Численное моделирование распространения волн в

геологических средах

4.1 Введение

4.2 Валидация метода и программного комплекса на двумерных сетка

4.3 Валидация метода и программного комплекса на трехмерных сетках

4.4 Проверка сеточной сходимости

4.5 Моделирование сложной двумерной топографии

4.5.1 Моделирование сложной синтетической топографии

4.5.2 Моделирование реалистичной синтетической топографии

4.6 Моделирование сложной трехмерной топографии

4.6.1 Моделирование синтетической топографии

4.6.2 Моделирование реальной топографии

4.7 Моделирование границы раздела сред на шельфе

Заключение

Список литературы

Список рисунков

Список таблиц

Введение диссертации (часть автореферата) на тему «Компьютерное моделирование пространственных динамических процессов в сложных областях интегрирования с использованием наложенных сеток»

Введение

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

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

Общей характеристикой всех этих методов является то, что они инициируют короткий (по времени) волновой импульс. При численном моделировании

вместо записей сигнала от настоящего источника возмущений, для удобства и простоты обычно подставляются синтетические импульсы. Чаще всего при этом используются Рикера или Ормсби [7]. Использование коротких импульсов не является обязательным условием, существуют также технологии сейсморазведки с использованием постоянно работающих источников с известной частотой [5; 8].

В любом случае, после инициации волны сейсмоприемники фиксируют как саму волну, так и отклик на нее. Причиной возникновения этого отклика служит неоднородность среды. Наличие под поверхностью слоев различных пород приводит к тому, что на границах между ними происходит резкое изменение физических свойств среды, где возникают вторичные волны, их которых и состоит измеряемый отклик [9; 10]. Далее этот отклик используется для решения задачи сейсмической инверсии, которая состоит в построении приблизительной модели среды. Для ее решения сперва берется некоторое приближение по известным геологическим свойствам среды в рассматриваемом географическом регионе, после чего это приближение итеративно уточняется так, чтобы минимизировать отличие прогнозируемого отклика модельной среды от экспериментально измеренного [11—15].

Другим примером задачи томографии, который схож с сейсмической томографией, является неразрушающее исследование внутренней структуры тел на предмет наличия внутренних дефектов, трещин и полостей [16; 17]. Дополнительную схожесть этим задачам придает тот факт, что на практике сейсморазведка чаще всего используется для обнаружения залежей полезных ископаемых, в частности нефти и газа. Геологические структуры, насыщенные нефтью и газом называются коллекторами и обычно имеют пористую или трещинноватую структуру [18], которая также создают отраженные сейсмические волны [19]. Существенное отличие этих задач состоит в том, что при дефектоскопии приемники могут быть закреплены по всей поверхности тела и использоваться для регистрации проходящих волн, в то время как в сейсмологии это потребовало бы погружения приемников в толщу породы на значительную глубину, что значительно увеличило бы стоимость исследования, поэтому в основном исследуются отраженные волны. В терминах дефектоскопии [20] можно сказать, что сейсморазведка полагается в основном на эхо-метод, а использование теневого метода затруднительно.

Другой важной задачей сейсмики является моделирование землетрясений и анализ сейсмостойкости строений. При проведении таких исследований обыч-

но можно полагаться на ранее построенную модель среды под поверхностью, а интерес в основном вызывает построение картин волнового поля и анализ характеристик волн, приходящих на поверхность. В отличии от сейсморазведки, для которой характерно неглубокое залегание источника волн, гипоцентр землетрясения обычно лежит на значительной глубине (землетрясение считается «нормальным» если эта глубина не превышает 70 километров, а для глубокофокусных она превышает 300 километров [21]). Различные модели землетрясений и анализ их воздействия на строения приведены в [22--27].

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

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

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

Поверхность земли при моделировании обычно описывается при помощи граничных условий свободной поверхности. При достижении поверхности, волна отражается от нее и продолжает распространение вглубь среды. В случае плоской поверхности особых трудностей моделирование этого процесса не вызывает. Однако в случаях, когда поверхность не плоская а имеет рельеф, при отражении от нее происходит изменение формы волны, интерференция и прочие эффекты, пренебрежение которыми существенно снизит точность моделирования [28; 29].

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

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

К наиболее известным методом решения задачи сейсмического моделирования относятся метод конечных элементов, конечно-разностные методы и сеточно-характеристический метод.

Метод конечных элементов [30; 31] основан на разбиении моделируемой области на некоторое количество элементов (непересекающихся подобластей) и задании на каждом из них аппроксимирующей функции простого вида. Дополнительно на эти функции накладывается требование равенства значений на границах, чтобы результирующая кусочная аппроксимационная функция на всей области была непрерывной. При использовании метода конечных элементов важно удачно выбрать вид аппроксимирующих функций. Например, при моделировании волновых процессов удачный выбор базиса может привести к тому, что разложение функции значений в моделируемой области окажется точным,

то есть кусочная аппроксимирующая функция не будет отличаться от истинных значений. Важными подклассами метода конечных элементов являются метод спектральных элементов и метод Галеркина, использующие в качестве базиса тригонометрические многочлены, полиномы Эрмита, Чебышева или Лежанд-ра [32—35]. Одним из преимуществ метода конечных элементов является то, что при его использовании значение функции всегда известно во всей области, а не только в дискретном наборе точек (узлов сетки). Также в [36—39] показана возможность очень точного моделирования идеальных поглощающих границ методом конечных элементов при помощи метода PML (perfectly matched layer).

Метод конечных разностей — это сеточный метод, который основан на представлении в решаемых дифференциальных уравнениях производных в виде определяемых используемой разностной схемой выражений, не содержащих операции дифференцирования. При использовании этого метода важно выбрать или построить подходящую разностную схему с учетом свойств заменяемого ей дифференциального оператора. При этом метод достаточно прост в реализации и не очень требователен к вычислительным ресурсам [40].

Сеточно-характеристический метод [41—43], используемый в данной работе, также основан на использовании сеток, однако вместо замены дифференциальных операторов разностными схемами он основан на замене переменной, при помощи которой решаемое уравнение сводится к уравнениям переноса, и нахождении характеристик, вдоль которых этот перенос происходит.

Несмотря на большое количество уже проведенных исследований, дальнейшее развитие численного моделирования в сейсмологии остается актуальной задачей. Большое количество ее практических применений порождает спрос со стороны промышленности. При этом постоянно возрастают требования к точности моделирования, для удовлетворения которых требуется затрачивать все больше вычислительных ресурсов. При этом прирост производительности процессоров на ядро (поток) замедлился. Частично это компенсируется развитием векторизации и SIMD-инструкций, которое приводит к появлению более быстрых алгоритмов для решения часто встречающихся задач, например генерации случайных чисел [44], сортировки [45; 46] или умножения матриц [47], а также используется и в самом моделировании [48]. Тем не менее, необходимо больше внимания уделять выбору менее ресурсоемких методов моделирования, оптимальной их реализации и использованию массово-параллельных архитектур. Также стоит отметить стремительное развитие и удешевление гетерогенных вычислительных

систем. Для ускорения вычислений начали активнее использоваться графические процессоры и специализированные аппаратные модули. Повысилось и удобство их использования, помимо развития CUDA [49—51], ОрепАСС [52; 53] и Ореп^ [54; 55] поддержка ускорителей появилась в ОрепМР [56; 57]. Представляет интерес также технология SYCL [58—60].

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

Для достижения поставленной цели необходимо было решить следующие задачи:

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

2. Разработка вычислительного метода, позволяющего использовать метод наложенных сеток в сочетании с сеточно-характеристическим методом на структурированных сетках

3. Реализация предложенного вычислительного метода в виде программного комплекса, предназначенного для численного моделирования волновых процессов в неоднородных средах

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

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

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

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

Научная новизна:

1. Предложена модификация сеточно-характеристического метода с использованием наложенных сеток для численного решения пространственных динамических задач в сложных областях интегрирования

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

3. Исследована возможность использования предложенного метода для моделирования распространения динамических процессов в областях с реалистичным рельефом поверхности

4. Исследовано влияние выбора методов интерполяции на точность результатов моделирования с использованием предложенного метода

Методология и методы исследования.

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

Основные положения, выносимые на защиту:

1. Предложена модификация сеточно-характеристического метода с использованием наложенных сеток для численного решения пространственных динамических задач в сложных областях интегрирования

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

3. Исследована возможность использования предложенного метода для моделирования распространения динамических процессов в областях с реалистичным рельефом поверхности

4. Исследовано влияние выбора методов интерполяции на точность результатов моделирования с использованием предложенного метода

5. Разработка программного комплекса для численного решения пространственных задач о распространении волновых процессов в гетерогенных средах в сложных областях интегрирования

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

Апробация работы. Основные результаты работы докладывались на:

1. 62-й Всероссийской научной конференции МФТИ (г. Долгопрудный,

2019)

2. 63-й Всероссийской научной конференции МФТИ (г. Долгопрудный,

2020)

3. 64-й Всероссийской научной конференции МФТИ (г. Долгопрудный, 2021)

4. 65-й Всероссийской научной конференции МФТИ (г. Долгопрудный, 2023)

5. Конференции «The 3rd BRICS Mathematics Conference July 21-26, 2019 Innopolis University, Russia» (г. Казань, 2019)

6. Конференции «Quasilinear Equations, Inverse Problems and Their Applications 2019» (QIPA-2019)

7. Конференции «Quasilinear Equations, Inverse Problems and Their Applications 2020» (QIPA-2020)

8. Международной конференции "Математические идеи П.Л. Чебышёва и их приложения к современным проблемам естествознания"(г. Обнинск, 2021)

9. Конференции «EAGE 2020 Annual Conference Exhibition Online» (Онлайн, 2020)

10. Конференции «10th International Conference "Distributed Computing and Grid Technologies in Science and Education"» (GRID'2023) (г. Дубна, 2023)

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

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

Публикации. Основные результаты по теме диссертации изложены в 4 печатных изданиях [29; 61—63], из них 3 — входящих в Scopus, WoS или перечень рекомендованных ВАК [61—63].

Объем и структура работы. Диссертация состоит из введения, 4 глав, заключения. Полный объём диссертации составляет 127 страниц, включая 64 рисунка и 7 таблиц. Список литературы содержит 161 наименование.

Глава 1. Определяющие уравнения

1.1 Введение

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

Вектор смещений в точке пространства, заданной радиус-вектором х обозначим как и(х, £). Можно ввести тензор деформации Коши-Грина [64]:

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

где ^ — компонента плотности поля сил по этой оси, — тензор напряжений, р — плотность.

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

1.2 Уравнения теории упругости

(1.1)

(1.2)

3 3

(1.3)

к=1 1=1

Тензор С^ы называется тензором жесткости [65]. Поскольку тензоры £ и а симметричны, то в каждом из них есть не более 6 независимых компонент. Кроме того, мы предполагаем, что существует дифференцируемая по £ функция плотности энергии деформации W, и что а^ = , поэтому тензор С также симметричен. Тогда можно представить этот тензор как матрицу Сар в нотации Фойгта [66]:

Сав =

С11

С12 С13 С14 С15 С16

С12 С22 С23 С24 С25 С26

С13 С23 С33 С34 С35 С36

С14 С24 С34 С44 С45 С46

С15 С25 С35 С45 С55 С56

С16 С26 С36 С46 С56 С16

(1.4)

Этот тензор содержит 21 независимую переменную, но если известен тип анизотропии среды, то их количество можно уменьшить.

1.3 Теория упругости в изотропной среде

В данной работе рассматриваются только изотропные среды. В них вид тензора С упрощается, для его описания достаточно всего двух параметров — Л и ц (параметры Ламе) [65]:

Сав =

Л + 2ц Л Л 0 0 0

Л Л + 2ц Л 0 0 0

Л Л Л + 2ц 0 0 0

0 0 0 ц 0 0

0 0 0 0 ц 0

0 0 0 0 0 ц

Закон Гука принимает вид:

агг = (Л + 2ц)£ц + Л £ык = Л £кк + 2Ц£

к=1

агз = Ц^ ^ (&гк+ Ь]к) £к1 = 2Ц£^

(1.5)

(1.6) (1.7)

Можно объединить эти случаи в одно уравнение:

3

= Л£кк + (1.8)

к=1

Или, возвращаясь от тензора деформации к вектору смещений:

3

дик (диг ди. \

= дх-к + 4 щ + ах;) (1-9)

к=1 4 7

Введем вектор скорости V как производную по времени вектора смещений и. Тогда можно записать систему определяющих уравнений так:

р I = у. «от + Г

д сг

— = Л(У . V)! + р(У ® V + (У ® у)т)

1.4 Параметры среды

(1.10)

Кроме рассмотренного выше способа описания свойств упругой среды с помощью задания параметров Ламе и плотности, на практике встречаются и другие варианты [67]. Вместо параметров Ламе могут использоваться объемный модуль упругости К и модуль сдвига О:

К=Лр

(1.11)

О = р

или модуль Юнга Е и коэффициент Пуассона V:

Е = р(3Л + 2р)

Л + р (1.12) Л

V =

2Л + 2р

Также возможна параметризация через задание скоростей продольных (Р-wave) и поперечных (S-wave) волн в среде:

Ср

'Л + 2ц Р

(1.13)

С* = ,/ р

Иногда используется модуль Р-волны М:

М = Л + 2ц (1.14)

1.5 Теория акустики

В точке с радиус-вектором х в момент времени Ь обозначим давление, плотность и скорость как рА(х,Ь), рА( х,Ь) и "А(х,Ь) соответственно. Запишем два основных уравнения гидродинамики в этой точке [68]. Уравнение непрерывности:

дрА + V- (рау~л) = 0 (1.15)

Уравнение Эйлера:

+ ("А V)vЛ = - — Vp (1.16)

дЬ рл

Мы предполагаем, что до начала распространения волн среда находится в состоянии покоя. Обозначим давление, плотность и скорость в состоянии покоя как ро( х, Ь), ро( х, Ь) и "о( х, Ь), причем по определению щ(х, Ь) = 0. Тогда можно выделить изменения в параметрах, происходящие из-за распространения волновых возмущений:

р(х, Ь) = рл(х, Ь) — ро(х, Ь)

р( х, Ь) = рл( х, Ь) — ро( х, Ь) (1.17)

у(х, Ь) = тГА(х, Ь) — щ(х, Ь) = уа(х, Ь) Подставим эти соотношения в (1.15) и (1.16):

д д р 0 = — (ро + р) + V - ((ро + р)") = дЬ + PоV - V + V - (р") (1.18)

дУ ду

0 = (ро+р)—+(Ро+р)( уУ)у+У(ро +р) = (ро+р)д^+(Ро+Р)( УУ) У+УР (1.19)

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

др + роУ- V =0 (1.20)

дЬ

ду 1

— + - Ур = 0 (1.21)

дЬ ро

Процесс распространения волн в идеальной жидкости является адиабатическим. Можно записать связь плотности с давлением и скорость звука в среде:

р=(Ю). р

(1-23)

Если подставить эти соотношения в уравнение непрерывности, то оно примет следующий вид:

др + рос2У - V =0 (1.24)

дЬ

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

ср = с (1.25)

с, = 0 (1.26)

р = 0 (1.27)

ср =к\Л (1.28)

К = Л (1.29)

Тензор напряжений в этом случае имеет следующий вид:

Оу- = -рЬг3 (1.30)

Также можно рассмотреть случай, когда акустическая среда до начала распространения волн находится не в состоянии покоя, а в ней присутствуют установившиеся течения, то есть v0 = 0, однако задачи в такой постановке в данной работе не исследовались.

Возможна ситуация, когда изменения плотности акустической среды не могут считаться малыми величинами. В данной работе из акустических сред рассматривается только вода. На глубине 500 метров давление примерно равно 5 МПа. Сжимаемость воды составляет около 5 - 10-1оПа-1 [69], то есть изменение плотности на такой глубине можно оценить как 0.25%. Избыточное давление, создаваемое распространяющимися волнами от источника также явно не превышает 1 МПа, то есть изменение плотности действительно мало.

В некоторых случаях акустическое приближение можно использовать и в упругих средах. Применимость такого приближения исследована в работах [70; 71].

Глава 2. Математическая модель

2.1 Введение

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

2.2 Сеточно-характеристический метод

2.2.1 Матричные уравнения линейной теории упругости в декартовых

координатах

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

Соберем все неизвестные в (1.10) в один вектор и запишем ее с учетом введенного обозначения [72]:

9 =

VI У2 У3 Ол 022 033 023 013 012

д

98 97 9б

т

д_ дь

94 =

д д д

р дЪ91 = 94 + дх1 0 99 + " дх2

д д д

р дЪ92 = - 0 99 + дх1 я 95 + " дх2

д д д

р дЪ93 = -о 98 + дх1 я 97 + " дх2

Список литературы диссертационного исследования кандидат наук Стецюк Владислав, 2023 год

источников

3.4.5 Зависимость погрешности от количества точек

На Рис. 3.18 представлен график зависимости погрешности интерполяции от количества точек с известными значениями. Погрешность вычислялась по норме L2, количество точек с известными значениями варьировалось от 240 до 2160. Для тестирования использовалась функция ftest(r) = sin7 (0.7r) в области [0; 10] х [0; 10] — аналогичная использовавшейся в примерах, но с большим показателем степени.

Метод kNN-IDW использовался с параметрами k = 4, p = 1.5. Можно видеть, что он вносит наибольшую погрешность. При этом если количество известных точек невелико, то его погрешность сопоставима с погрешностью метода барицентрических координат, поэтому его использование может быть осмысленно в случаях, когда интерполяция проводится из сетки с большим шагом в сетку с существенно меньшим.

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

Дискретный метод Сибсона при любом количестве точек-источников уступает в точности обычному методу Сибсона, что ожидаемо. Его точность можно улучшить, уменьшив шаг дискретизации, но это существенно уменьшит скорость его работы.

Метод локального разложения по базису радиальных функций воспроизводит значения функции точнее всего, однако это достигается за счет большего времени работы. Также следует отметить, что значения £ подстраивались вручную для снижения погрешности (для базиса гауссиан использовалось £ = 0.67, для обратных мультиквадриков — £ = 1 и для ядра Матерна О2 — £ = 1.67).

3.5 Многопоточная реализация интерполяции

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

Процесс интерполяции в реализованном программном комплексе можно разбить на три этапа:

1. Инициализация сеток и структур данных

2. Независимое рассмотрение точек целевой сетки и вычисление интерполяционных коэффициентов для них

3. Объединение и сохранение полученных результатов

При этом в большей части используемых методов, используемые при интерполяции структуры данных никак не модифицируются, что значительно облегчает параллельную реализацию. Это неверно в случае метода естественной окрестности, где происходит добавление и удаление точек в диаграмму Вороного, однако параллельная реализация этого алгоритма также возможна, как показано в [138—140].

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

12345678 12345678

Число ПОТОКОВ Число ПОТОКОВ

Рисунок 3.19 — Графики зависимости ускорения и эффективности использования

многопоточности от количества потоков

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

Графики зависимости времени интерполяции от количества используемых потоков приведены на Рис. 3.19.

3.6 Генерирование реалистичных синтетических карт высот

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

>Г №7»

О 20 40 60 80 100 120 140 0 20 40 6 0 80 100 120 140 0 2 0 40 60 80 100 120 140 0 20 40 60 80 100 120 140

а) 2 октавы

б) 4 октавы

в) 8 октав

г) Взвешенная

сумма

Рисунок 3.20 — Двумерный шум Перлина с разным количеством октав

За основу будущего рельефа берется шум. В реализованном генераторе используется шум Перлина [141; 142] — градиентный шум с плавными переходами и хорошим выделением областей. Примеры шума Перлина с разным количеством октав представлены на Рис. 3.20. Поскольку в реальном рельефе присутствуют особенности разных масштабов, мы генерируем шум Перлина с разным количеством октав, после чего складываем полученные значения с разными весами. Мы используем шумы с 2, 4, 8 и 16 октавами, и коэффициенты 1, 1, 8 и 16 соответ-

ственно.

150

150

150

а) Степень 0.9

б) Степень 1.25

в) Степень 1.45

Рисунок 3.21 — Влияние показателя степени на шум

Далее полученные значения возводятся в степень. Это нужно для более явного разделения областей с большими значениями. На Рис. 3.21 приведены примеры карт высот с разными показателями степени. Видно, что выбор большего показателя приводит к появлению явно выраженных "пиков"и превращению остальной карты высот в "равнину". Значение показателя 1.25 обычно дает достаточно хорошие результаты.

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

H

gauss

(х,У) =

L + exp

(x X center ) + (y У center )

1 + L

И (x,y)

(3.6)

2

U

0.125 0.100 0.075 0.050 0 0.025 25 0.000 -0.025 -0,050

125 150

150

125 150

150

а) Ь = 0.05, а = 0.2 б) Ь = 0.05, а = 0.35 в) Ь = 0.05, а = 0.65 Рисунок 3.22 — Влияние параметра а на синтетическую карту высот

125 150

150

125 150

150

а) Ь = 0.2, а = 0.35 б) Ь = 0.35, а = 0.35 в) Ь = 0.5, а = 0.35

Рисунок 3.23 — Влияние параметра Ь на синтетическую карту высот

Параметр Ь задает вклад "гауссианизации"в значения, а а определяет ширину пиков. Примеры карт высот для разных значений этих параметров приведены на Рис. 3.22 и 3.23.

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

Нте^атге(х,уУ)

Н {х,у),

Н(х,у) ^ М

М + (Н(х, у) - М) • С, Н(х, у) <М

(3.7)

Здесь М — медианное значение функции в области, а С — параметр, задающий вклад преобразования. Иллюстрация вклада этого параметра в результат представлена на Рис. 3.24. Универсальные значения параметров а, Ь и С найти не

150 150 150

а) С = 0.35 б) С = 0.6 в) С = 0.9

Рисунок 3.24 — Влияние параметра Ь на синтетическую карту высот

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

Последним шагом является нормировка по высоте, после чего карта высот построена и может быть использована для построения сетки и проведения вычислительных экспериментов.

Рисунок 3.25 — Вариант сетки, описывающей синтетическую топографию

3.7 Построение сеток для описания поверхности

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

Рисунок 3.26 — Построение сетки при помощи нормалей

Для построения сетки необходимо вычислить координаты узлов, входящих в ее следующие слои. Нулевым слоем мы считаем саму границу. Первый вариант вычисления координат узлов следующего слоя — провести нормали к поверхности (кривой) предыдущего слоя, сместиться вдоль этих нормалей на какое-то расстояние и зафиксировать новые узлы. На Рис. 3.26 приведена иллюстрация построения сетки таким способом. При этом может иметь смысл выбрать в нулевом слое узлы так, чтобы расстояние между ними было постоянным или вдоль кривой, или по одной из координат. У такого варианта построения есть преимущество, связанное с тем, что ячейки сетки по форме остаются близки к прямоугольникам. Его основной недостаток состоит в том, что соотношение сторон этих прямоугольников может быть далеко от 1.

Рисунок 3.27 — Построение сетки при помощи вертикального переноса

Другой метод проиллюстрирован на Рис. 3.27 и заключается в том, что следующий слой получается смещением по вертикальной оси относительно

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

Рисунок 3.28 — Построение сетки при помощи нормалей с ограничением смещения

Для того, чтобы частично компенсировать недостатки этих методов, можно использовать некий промежуточный вариант: при построении следующего слоя вычислять координаты новых узлов при помощи построения нормалей, но ограничивать смещение узлов относительно предыдущего слоя по осям, перпендикулярным вертикальной. В этом случае ячейки получаются более близкими к прямоугольным, чем при использовании вертикального сдвига, а их соотношение сторон ближе к 1 чем при использовании только нормалей. Пример построения сетки таким способом представлен на Рис. 3.28

а) Без сглаживания б) Со сглаживанием

Рисунок 3.29 — Сглаживание слоев сетки при помощи фильтра Савицкого-Голея

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

В некоторых случаях имеет смысл использовать сглаживание поверхностей слоев. Чаще всего достаточно применить его только к верхнему слою, но не всегда. Для сглаживания в реализованном в рамках данной работы программном средстве используется фильтр Савицкого-Голея [143—145]. Его использование проиллюстрировано на Рис. 3.29.

Глава 4. Численное моделирование распространения волн в геологических

средах

4.1 Введение

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

4.2 Валидация метода и программного комплекса на двумерных сетка

Для проверки корректности моделирования с использованием предложенного метода в двумерном случае, было проведено сравнение результатов с квазианалитическим решением и с другим известным программным комплексом.

Моделируемое пространство разделено плоскостью на два полупространства. Одно из них заполнено однородной изотропной упругой средой с известными характеристиками, другое — вакуумом. Под плоской поверхностью среды расположены источник возмущений и приемник, их положения известны. Рассматривалось две постановки задачи: с источником, направленным к поверхности (задача Лэмба [146]) и со всенаправленным источником сжатия (задача Гарвина [147]).

Для большей наглядности, в расчете с использованием наложенных сеток граница полагалась наклоненной на 27°. Параметры вычислительного эксперимента приведены в таблице 1

Таблица 1 — Параметры валидационного численного эксперимента по численному моделированию задач Гарвина и Лэмба в двумерном случае

Параметр Значение

Размер основной сетки 601 х 601 узел

Шаг основной сетки 2 м

Размер наложенной сетки 1001 х 11 узлов

Шаг наложенной сетки 2м

Скорость продольных волн ср 2500 м/с

Скорость поперечных волн с8 1558 м/с

Плотность среды р 1500 кг/м3

Шаг по времени 5•10"4с

Расстояние от источника до границы 250 м

Расстояние от приемника до границы 400 м

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

Рисунок 4.1 — Иллюстрация к постановке задач Гарвина и Лэмба в двумерном пространстве и картина волнового поля через 400 мс в задаче Лэмба

Компонента vx

_1 ---- ex2dvael - rect -

/

1

1 УЛ

л 1 V/

\]

0.0 0.1 0.2 0.3 0.4 0.5 0.6

Компонента чу

_ _ fael -

л - rect

/ \ - error

/ \

J \

т

1

1 *

0.0 0.1 0.2 0.3 0.4 0.5 0.6

Рисунок 4.2 — Сравнение результатов моделирования задачи Гарвина

На Рис. 4.1 приведены иллюстрация используемой постановки и картина волнового поля в эксперименте по моделированию задачи Лэмба через 400 мс после точки отсчета. Зеленый квадрат соответствует положению источника, желтый квадрат — положению приемника. В качестве источника использовался вейвлет Рикера с пиковой частотой 14.5 Гц. Во всем, кроме типа источника, постановка задач была идентична. Для описания границы использовалась наложенная сетка. Перенос значений осуществлялся с использованием метода локального разложения и билинейного базиса.

Результаты для сравнения были получены с использованием программных пакетов ex2ddir [148], ex2dvael [149] и specfem2d [150—152].

Для построения графиков сравнения было выполнено приведение компонент скорости, зарегистрированных приемником, к осям до поворота. На Рис. 4.2 представлены показания приемника в задаче Гарвина, а на Рис. 4.3. Также там представлен график погрешности результатов моделирования предложенным методом относительно квазианалитического решения, полученного при помощи пакета ex2dvael в задаче Гарвина и ex2ddir в задаче Лэмба. Можно видеть, что и в задаче Гарвина, и в задаче Лэмба погрешность мала.

Рисунок 4.3 — Сравнение результатов моделирования задачи Лэмба 4.3 Валидация метода и программного комплекса на трехмерных сетках

Валидация предложенного метода в трехмерном случае заключалась в сравнении результатов численного решения задачи Гарвина без использования наложенных сеток (отражение волны от границы кубической сетки) и с использованием повернутой сетки для описания границы. В обоих случаях использовался один и тот же вычислительный комплекс. Параметры эксперимента приведены в таблице 2.

Угол поворота наложенной сетки составлял 17°. В качестве источника использовался вейвлет Рикера с пиковой частотой 14.5 Гц. Для построения графиков сравнения было выполнено приведение компонент скорости, зарегистрированных приемником, к осям до поворота. Как видно на Рис. 4.4, различия в результатах моделирования незначительны.

Таблица 2 — Параметры валидационного численного эксперимента по численному моделированию задачи Гарвина в трехмерном случае

Параметр Значение

Размер основной сетки 301 х 301 х 301 узел

Шаг основной сетки 2м

Размер наложенной сетки 421 х 421 х 11 узлов

Шаг наложенной сетки 2м

Скорость продольных волн ср 2500 м/с

Скорость поперечных волн с8 1558 м/с

Плотность среды р 1500 кг/м3

Шаг по времени 5 • 10"4с

Расстояние от источника до границы 60 м

Расстояние от приемника до границы 100 м

Расстояние между источником и приемником вдоль оси Ох 50 м

Расстояние между источником и приемником вдоль оси Оу 100 м

4.4 Проверка сеточной сходимости

Для измерения порядка сеточной сходимости использовалась постановка эксперимента как на Рис. 4.5. В области с изотропной однородной средой задана начальным условием продольная волна, распространяющаяся вдоль оси Ох. Часть области дополнительно покрыта наложенной сеткой. По обеим осям на основную сетку наложены периодические граничные условия.

Параметры двух постановок, в которых был проведен эксперимент, приведены в таблице 3. Порядок сеточной сходимости определялся исходя из погрешности, то есть разности между начальным и конечным состоянием, вычисленной по норме Ь\. На Рис. 4.6 представлены графики погрешностей для различных методов интерполяции при угле наклона наложенной сетки равном 30°. Можно видеть, что использование наложенной сетки увеличивает погрешность, однако сходимость все еще наблюдается.

На 4.7 представлены графики зависимости погрешности от угла поворота наложенной сетки. Можно видеть, что погрешность максимальна при угле поворота наложенной сетки близком к 90°, что связано с тем, что в этом случае

Компонента V„

0.00 le-12

0.5 0.0 -0.5 -1.0 -1.5 -2.0

псе

---- refere - overse - error

/ Л

Компонента V»

0.10 0.15 0.20

Компонента v7

/Л ---- reference - overset -

/ Л - error

\J

---- refere _ псе it -

- error

Л

\

^—'

vy v

Рисунок 4.4 — Сравнение результатов моделирования трехмерной задачи Гарвина

Таблица 3 — Параметры численного эксперимента по определению порядка

сеточной сходимости

Параметр

Постановка 1 Постановка 2

Размер основной сетки 360 х 360 узлов 720 х 720 узлов

Шаг основной сетки 0.025 м 0.0125 м

Размер наложенной сетки 80 х 80 узлов 160 х 160 узлов

Шаг наложенной сетки 0.025 м 0.0125 м

Скорость волны с 10 м/с

Плотность среды р 20 кг/м3

Шаг по времени 5 • 10"4с 25 • 10"5с

Рисунок 4.5 — Иллюстрация к постановке эксперимента по измерению порядка

сеточной сходимости

Errors Li, 30'

_1_1 Без наложенной сетки Билинейный базис - Барицентр, координаты - RBF Inverse Multiquadric - Метод Сибсона

200 400 600 800 1000 1200

Рисунок 4.6 — Зависимость погрешности от количества узлов основной сетки

Error i_i

_1_ - idle bilinear

- baryce - rbfjm - sibson ntric q

0 25 50 75 100 125 150 175

Рисунок 4.7 — Зависимость погрешности от угла поворота наложенной сетки

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

Orders Li

2.9

О 25 50 75 100 125 150 175

Рисунок 4.8 — Зависимость порядка сходимости от угла поворота наложенной

сетки

Итоговый график зависимости порядка сеточной сходимости от угла поворота наложенной сетки представлен на Рис. 4.8. Можно видеть, что при небольших углах поворота порядок сходимости мало отличается от случая без наложенной сетки, а при углах порядка 90° порядок сходимости наименьший, то есть этот график качественно аналогичен графику для погрешностей. Также интересно, что интерполяция методом барицентрических координат практически не отличается от билинейной ни по вносимой погрешности, ни по порядку сходимости.

Рисунок 4.9 — Иллюстрация к постановке модельной задачи

4.5 Моделирование сложной двумерной топографии 4.5.1 Моделирование сложной синтетической топографии

Т| " "

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

Иллюстрация к постановке задачи приведена на Рис. 4.9. Верхняя граница наложенной сетки повторяет положение свободной поверхности. Зеленой точкой отмечено положение точечного источника, а белыми точками справа — положения приемников. Всего показания снимались со 120 приемников, расположенных с шагом 20 м по оси Оу (глубине).

В таблице 4 приведены используемые в расчете параметры среды, размеры сеток и шаг по времени. На Рис. 4.10 приведены картины волнового поля в различные моменты времени. На них видно, как первичная волна от источника постепенно достигает свободной поверхности и отражается от нее. Сначала отраженная волна в общих чертах повторяет изогнутый контур поверхности, это хорошо заметно по волновым полям через 360 и 440 мс после точки отсчета. Из-за формы поверхности под ней возникает интерференция, в результате которой

Таблица 4 — Параметры численного эксперимента по моделированию сложной

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

Параметр Значение

Размер основной сетки 1000 х 500 х узлов

Шаг основной сетки 4м

Размер наложенной сетки 1000 х 25 узлов

Шаг наложенной сетки « 4 м

Скорость продольных волн ср 2500 м/с

Скорость поперечных волн с3 1558 м/с

Плотность среды р 1500 кг/м3

Шаг по времени 2 • 10"4с

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

В этой же постановке было проведено моделирование с другим типом источника: вместо источника напряжения использовался источник силы, направленный вертикально вниз. По показаниям, регистрируемым приемниками, были построены синтетические сейсмограммы, приведенные на Рис. 4.11. По первой приходящей волне на сейсмограмме для компоненты скорости уу виден приемник, находящийся на одной глубине с источником, так как на нем компонента скорости уу первой волны, приходящей от источника до отражения равна нулю, а графики показаний остальных приемников для этой волны симметричны относительно него.

Показания этого приемника были использованы для сравнения полученных результатов с результатами моделирования задачи в такой же постановке при помощи программного пакета specfem2d [150—152]. Результаты сравнения приведены на Рис. 4.12. Видно, что результаты моделирования предложенным методом с хорошей точностью совпадают с результатами моделирования при помощи specfem2d.

а) 120 мс

в) 360 мс

д) 520 мс

б) 240 мс

г) 440 мс

е) 600 мс

ж) 680 мс з) 760 мс

Рисунок 4.10 — Картины волнового поля в эксперименте по моделированию син тетической двумерной топографии в различные моменты времени

а) Компонента скорости ух б) Компонента скорости уу

Рисунок 4.11 — Сейсмограммы компонент скорости в эксперименте по моделированию синтетической двумерной топографии

Компонента vx

---- specfem2d - rect

- différé rice

л \л/ 1 !К1л(\

ЛГИ W [ У ^ \J Va-

V 11 V г

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

Компонента vy

л ---- specfem2d - rect

- différé псе

j А

V \

Рисунок 4.12 — Сравнение показаний приемника

4.5.2 Моделирование реалистичной синтетической топографии

При помощи предложенного метода было проведено численное моделирование в области с реалистичной топографией. В качестве источника данных использовалась синтетическая модель Canadian Foothills, часто используемая для тестирования методов сейсмического моделирования и решения задач сейсмической инверсии [153—156].

На Рис. 4.13 представлена картина слоев с различными свойствами среды. Сверху темно-синим цветом выделена наложенная сетка, описывающая свободную границу. Размер основной расчетной сетки составлял 1200 х 447 узлов с шагом 20.83 м, размер наложенной сетки — 1200 х 7 узлов с шагом примерно 21 м в среднем. Шаг по времени был равен 1 мс

В качестве источника возмущения использовался импульс Рикера пиковой частотой 15 Гц. Источник находился на глубине 1470 м под поверхностью, на расстоянии 12550 м от левой границы области и 12450 м от правой. Приемники находились на глубине 100 метров под поверхностью, расстояние между ними по оси Ox составляло 100 м. Всего моделировались показания 247 приемников, самый левый из которых находился на расстоянии 150 метров от левой границы области моделирования, а самый правый на расстоянии 150 метров от правой границы области.

На Рис. 4.14 представлены картины волнового поля в различные моменты времени. По картинам волнового поля через 625 мс и 750 мс после точки отсчета наглядно видно, как формируется отраженная от поверхности волна и происходит интерференция, в результате которой формируются отдельные волновые фронты,

Рисунок 4.13 — Физические параметры различных слоев среды в модели Canadian

Foothills

а) 625 мс

в) 875 мс

д) 1250 мс

б) 750 мс

г) 1000 мс

е) 1500 мс

ж) 1750 мс з) 2000 мс

Рисунок 4.14 — Картины волнового поля в эксперименте с двумерной моделью Canadian Foothills в различные моменты времени

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

На 4.15 представлены синтетические сейсмограммы компонент скорости, регистрируемой приемниками. Так же, как и в картинах волнового поля, в них видна основная волна большой амплитуды и большое количество ряби. Из-за того, что приемники располагались на небольшой глубине, разделение основной волны на первичную (от источника) и отраженную от свободной поверхности затруднительно.

а) Компонента скорости vx мс б) Компонента скорости vy

Рисунок 4.15 — Синтетические сейсмограммы компонент скорости, полученные в эксперименте с двумерной моделью Canadian Foothills

4.6 Моделирование сложной трехмерной топографии

4.6.1 Моделирование синтетической топографии

Для проверки корректности работы предложенного метода в областях со сложной синтетической топографией был проведен тестовый расчет на сетке, сгенерированной и использовавшейся в качестве иллюстрации в разделе 3.6. В качестве источника возмущения использовался точечный всенаправленный источник, в качестве сигнала был выбран вейвлет Рикера с пиковой частотой 15 Гц. Источник располагался вблизи поверхности на вершине горы. Параметры расчета приведены в таблице 5.

В ходе вычислительного эксперимента были получены картины волнового поля. На Рис. 4.16 представлены картины волнового поля на поверхности. Видно, как волна постепенно распространяется от вершины, и ее интенсивность убывает по мере распространения. Также видно, что форма фронта волны не является сферической, а интенсивность убывает неравномерно, что обусловлено формой рельефа поверхности.

Таблица 5 — Параметры численного эксперимента по моделированию реальной

трехмерной топографии

Параметр Значение

Размер основной сетки 136 х 136 х 151 узел

Шаг основной сетки 30 м

Размер наложенной сетки 150 х 150 х 5 узлов

Шаг наложенной сетки « 30 м

Скорость продольных волн ср 2850 м/с

Скорость поперечных волн с8 1650 м/с

Плотность среды р 2400 кг/м3

Шаг по времени 2•10"3с

На Рис. 4.17 представлены картины волнового поля в двух срезах основной сетки, пересечение которых лежит недалеко от вершины горы. Поскольку амплитуда волны достаточно быстро убывает, для наглядности две последних картины построены с другой цветовой карты. На картинах волнового поля видна первичная

в) 400 мс г) 560 мс

Рисунок 4.16 — Картины волнового поля на поверхности в эксперименте по моделированию синтетической трехмерной топографии в различные моменты

времени

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

а) 80 мс

в) 400 мс

д)720 мс

б) 240 мс

г) 560 мс

е) 880 мс

Рисунок 4.17 — Картины волнового поля в основной сетке в эксперименте по моделированию синтетической трехмерной топографии в различные моменты

времени

4.6.2 Моделирование реальной топографии

а) Местоположение на карте. Скриншот Google Maps [157]

б) Фрагмент топографической карты. Скриншот QGIS [158]

Рисунок 4.18 — Гора Нефин на физико-географической и топографической карте

С использованием описанного в данной работе метода было проведено моделирование распространения в области горы Нефин (см. Рис. 4.18). Наложенная сетка использовалась для описания рельефа поверхности. В качестве источника возмущений рассматривалась плоская волна, распространяющаяся вертикально к поверхности. Среда под поверхностью полагалась однородной. Карта высот была взята из данных проекта Copernicus DEM [159]. Параметры и подробности расчета приведены в таблице 6.

На Рис. 4.19 представлены картины волнового поля в разные моменты времени. Наложенная сетка представлена на них в виде поверхности, а основная сетка — двумя срезами по разным осям, пересечение которых находится под вершиной горы. Плоская волна начинает достигать поверхности примерно через 1.1 -1.2 с после точки отсчета, в волновой картине на поверхности начинают появляться ненулевые значения. Пик волны достигает поверхности примерно через 1.35 с, а примерно через 1.65 с он достигает вершины горы. По волновым картинам через 1.6 и 1.8 с видно, как каждый из склонов горы формируют отдельную отраженную волну, формой повторяющую поверхность. Далее по волновым картинам через 2.0 и 2.2 с можно видеть, как происходит интерференция волн, в результате которой возмущение под вершиной имеею большую амплитуду, чем амплитуда

Таблица 6 — Параметры численного эксперимента по моделированию реальной

трехмерной топографии

Параметр Значение

Размер основной сетки 215 х 220 х 160 узлов

Шаг основной сетки 25 м

Размер наложенной сетки 227 х 234 х 10 узлов

Шаг наложенной сетки « 25 м

Скорость продольных волн ср 2850 м/с

Скорость поперечных волн с3 1650 м/с

Плотность среды р 2400 кг/м3

Шаг по времени 2•10"3с

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

а) 0.4 с

в) 1.4 с

д) 1.8 с

б) 1.2 с

г) 1.6 с

е) 2.0 с

ж) 2.2 с з) 2.4 с

Рисунок 4.19 — Картины волнового поля в эксперименте по моделированию реальной трехмерной топографии в различные моменты времени

4.7 Моделирование границы раздела сред на шельфе

а) Местоположение на карте. Скриншот Google Maps [160]

б) Фрагмент топографической карты. Скриншот QGIS [158]

Рисунок 4.20 — Моделируемая область шельфа на физико-географической и топографической карте

Наложенные сетки могут использоваться не только для описания границы раздела между средой и вакуумом, но и на границах между средами. В этом расчете наложенная сетка использовалась для описания шельфовой топографии, то есть границы раздела упругой среды шельфа и акустической среды (воды). На Рис. 4.20 приведены положение моделируемой области на карте, а также фрагмент топографической карты. Карта глубин с разрешением 12 м была взята из

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

На Рис. 4.21 представлены визуализации областей с различными свойствами среды. Можно видеть, что без использования наложенной сетки граница раздела областей имеет характерную лестничную структуру.

В данном случае можно было бы применить другой подход к описанию границы раздела сред при помощи наложенных сеток. Вместо одной наложенной

[161].

Таблица 7 — Параметры численного эксперимента по моделированию рельефа

океанического шельфа

Параметр

Значение

Размер основной сетки Шаг основной сетки Размер наложенной сетки Шаг наложенной сетки

Скорость продольных волн в среде шельфа ср

Скорость поперечных волн в среде шельфа с8

Плотность среды шельфа р1

Скорость волн в воде с

Плотность волы р2

Шаг по времени

180 х 220 х 120 узлов 6 м

192 х 234 х 11 узлов

« 6 м

2850 м/с

1650 м/с

2400 кг/м3

1500 м/с

1050 кг/м3

10"3с

а) Только основная сетка б) С наложенной сеткой

Рисунок 4.21 — Области с различными свойствами среды в постановках с одной сеткой и с наложенной сеткой в эксперименте по моделированию шельфового рельефа

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

а) Только основная сетка б) С наложенной сеткой

Рисунок 4.22 — Срезы волнового поля через 380 мс в эксперименте по моделированию шельфового рельефа

На Рис. 4.22 представлены срезы волнового поля через 380 мс после точки отсчета моделирования. Видно, что картины качественно схожи, но при использовании наложенной сетки амплитуда шума, следующего за отраженной от границы раздела сред волной, меньше. Причиной этого служит то, что при использовании лестничной сетки волна от источника, отражаясь от лестничной границы, испытывает дифракцию и интерференцию. Аналогичные явления возникают при использовании лестничной сетки для описания свободной поверхности.

В данном вычислительном эксперименте приемники были расположены вблизи поверхности. Общее количество приемников — 3150, они были сгруппированы в 45 рядов по 70 приемников. Ряды были расположены вдоль оси Ох, расстояние между приемниками составляло 12 метров, и смещены друг относительно друга на 24 метра по оси Оу. Для наглядности сейсмограммы строились отдельно для каждого ряда.

На Рис. 4.23 представлены сейсмограммы компоненты скорости уг. Так же, как и по картинам волнового поля, можно видеть, что использование наложенной сетки позволяет значительно уменьшить амплитуду паразитных волн.

На Рис. 4.24 представлены сравнения значений графиков значений компонент скорости ух и уг на отдельном приемнике. Это сравнение еще раз подтверждает сделанный ранее вывод о том, что использование наложенной сетки позволяет заметно уменьшить количество шума и погрешность, вносимую паразитными волнами, делая тем самым более различимой общую картину.

а) Только основная сетка

Рисунок 4.23

б) С наложенной сеткой — Синтетические сейсмограммы компоненты скорости в эксперименте по моделированию шельфового рельефа

- No overset

I А Single overset

J 4 Шш, ^/ЧА/

J 1/ V ^ ттл/ч/ч

V

1

0 о о 2 0 4 0 .6 0 а 1. 0 1 2 1 4 1.6

а) Компонента скорости ух

| _ No overset Single overset

A

11

0.0 0.2 0.4 0.6

0,8 1.0

1.4 1.6

б) Компонента скорости Рисунок 4.24 — Сравнение графиков компонент скорости на отдельном приемнике в эксперименте по моделированию шельфового рельефа

Заключение