Математическое моделирование процесса сейсморазведки с учётом различия реологических свойств отдельных частей геологического массива сеточнохарактеристическими методами тема диссертации и автореферата по ВАК РФ 00.00.00, доктор наук Голубев Василий Иванович
- Специальность ВАК РФ00.00.00
- Количество страниц 233
Оглавление диссертации доктор наук Голубев Василий Иванович
ВВЕДЕНИЕ
ГЛАВА 1. ОПРЕДЕЛЯЮЩИЕ УРАВНЕНИЯ И РЕОЛОГИЧЕСКИЕ СООТНОШЕНИЯ
1.1. Введение
1.2. Акустическая среда
1.3. Линейно упругие среды
1.3.1. Изотропная среда
1.3.2. Трансверсально-изотропная среда
1.4. Флюидонасыщенная среда
1.5. Континуальная модель слоистой среды
1.6. Выводы
ГЛАВА 2. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ДИНАМИЧЕСКИХ ЗАДАЧ
2.1. Введение
2.2. Сеточно-характеристический численный метод
2.3. Одномерное линейное уравнение переноса
2.4. Пористые флюидонасыщенные среды
2.5. Акустические и линейно упругие среды
2.5.1. Акустическая среда
2.5.2. Изотропная линейно упругая среда
2.5.3. Трансверсально-изотропная среда
2.6. Периодические среды
2.7. Пространственное расщепление
2.8. Компактные схемы для акустических сред
2.8.1. Одномерная задача
2.8.2. Двумерная задача
2.9. Выводы
ГЛАВА 3. ГРАНИЧНЫЕ И КОНТАКТНЫЕ УСЛОВИЯ
3.1. Введение
3.2. Граничные условия
3.2.1. Акустические среды
3.2.2. Изотропные линейно упругие среды
3.2.3. Анизотропные (УТГ) линейно упругие среды
3.2.4. Пористые флюидонасыщенные среды
3.3. Контактные условия
3.3.1. Контакт изотропной и акустической сред
3.3.2. Контакт анизотропной и изотропной сред
3.3.3. Контакт флюидонасыщенной и изотропной сред
3.3.4. Контакт флюидонасыщенной и акустической сред
3.4. Исследовательское программное обеспечение
3.5. Выводы
ГЛАВА 4. СЕЙСМИЧЕСКИЕ ВОЛНЫ В ПОРИСТЫХ ФЛЮИДОНАСЫЩЕННЫХ СРЕДАХ
4.1. Введение
4.2. Сравнение с акустической и упругой моделями
4.3. Контакт со средами с различной реологией
4.4. Выводы
ГЛАВА 5. СЕЙСМИЧЕСКИЕ ВОЛНЫ В АНИЗОТРОПНЫХ СРЕДАХ
5.1. Введение
5.2. Волны в однородных поперечно-изотропных средах
5.3. Контакт поперечно-изотропных и изотропных сред
5.4. Реалистичные геологические модели
5.5. Осреднённые анизотропные модели
5.6. Выводы
ГЛАВА 6. ВОЛНЫ В ТРЕЩИНОВАТЫХ СРЕДХ. КОНТИНУАЛЬНАЯ МОДЕЛЬ
6.1. Введение
6.2. Сравнение различных моделей трещиноватых сред
6.3. Анизотропия отклика от континуальной модели
6.4. Учёт повреждаемости среды
6.5. Учёт влияния гравитационной силы
6.6. Трещиноватый объект в модели Marmousi II
6.7. Трёхмерная модель баженовской свиты
6.8. Выводы
ГЛАВА 7. ВОЛНЫ В ТРЕЩИНОВАТЫХ СРЕДХ. ЯВНОЕ ВЫДЕЛЕНИЕ ГРАНИЦ
7.1. Введение
7.2. Двумерные постановки задач
7.2.1. Изотропные линейно упругие среды
7.2.2. Поперечно-изотропные среды
7.2.3. Пористые флюидонасыщенные среды
7.3. Трёхмерные постановки задач
7.3.1. Волны в тонких пластинах
7.3.2. Явное выделение криволинейных границ
7.3.3. Оценка анизотропии отклика
7.4. Выводы
ГЛАВА 8. ОБРАТНЫЕ СЕЙСМИЧЕСКИЕ ЗАДАЧИ
8.1. Введение
8.2. Задача миграции в приближении Борна
8.2.1. Акустические среды
8.2.2. Изотропные линейно упругие среды
8.2.3. Результаты численных расчётов
8.3. Полноволновая задача миграции
8.4. Применение нейронных сетей
8.5. Выводы
ГЛАВА 9. ЗАДАЧИ ГЛОБАЛЬНОЙ СЕЙСМИКИ
9.1. Введение
9.2. Моделирование очага землетрясения
9.3. Сейсмические воздействия на объекты
9.4. Пример использования компактной схемы
9.5. Нагрузка ледового образования
9.6. Выводы
ЗАКЛЮЧЕНИЕ
СПИСОК ЛИТЕРАТУРЫ
Введение
Рекомендованный список диссертаций по специальности «Другие cпециальности», 00.00.00 шифр ВАК
Разработка численных методов для моделирования распространения упругих волн в неоднородных средах2015 год, кандидат наук Фаворская, Алена Владимировна
Метод исследования пространственных волновых явлений в средах со сложной структурой с помощью вычислительных экспериментов2019 год, доктор наук Фаворская Алена Владимировна
Численное решение пространственных динамических задач механики неоднородных деформируемых сред2014 год, кандидат наук Голубев, Василий Иванович
Моделирование распространения динамических волновых возмущений в гетерогенных средах на высокопроизводительных вычислительных системах2022 год, доктор наук Хохлов Николай Игоревич
Развитие схем на основе квазиодномерного подхода для решения задач аэроакустики на неструктурированных сетках2013 год, кандидат наук Бахвалов, Павел Алексеевич
Введение диссертации (часть автореферата) на тему «Математическое моделирование процесса сейсморазведки с учётом различия реологических свойств отдельных частей геологического массива сеточнохарактеристическими методами»
Актуальность темы исследования
Развитие численных методов решения гиперболических систем уравнений находит непосредственное применение в таких актуальных практических задачах как поиск и разведка месторождений полезных ископаемых, неразрушающий контроль современных сложно построенных материалов (например, композитных), оценка последствий и прогнозирование глобальной сейсмической активности планеты.
В настоящее время интерес к выявлению признаков трещиноватых сред в волновом поле возрастает и связан, во многом, с поиском месторождений нефти и газа, приуроченных к карбонатным резервуарам. Развитие компьютерных технологий и алгоритмов анализа данных стимулирует их использование для детального изучения волнового поля на месторождениях, где трещиноватость предполагается или зафиксирована по результатам бурения скважин. Среди технологий, которые активно внедряются в граф интерпретации, можно отметить азимутальную инверсию и полноазимутальную миграцию. Данные технологии апробированы на месторождениях Оренбургской области (месторождения Царичанское и Кувайское), Сербии (месторождение Турия). Благодаря этим технологиям стало возможным получить качественно новые характеристики волнового поля - получить информацию о том, как волновое поле меняется в разных направлениях. Однако открытым остается вопрос о том насколько сильно трещиноватость может влиять на волновое поле и какие особенности волнового поля могут быть использованы для изучения резервуара. Естественным выбором, для ответа на эти вопросы, является использование методов моделирования волнового поля. Исследование волнового поля на моделях анизотропных, трещиноватых и пористых средах выполняются длительный период времени. Эти работы посвящены как численным расчётам волнового поля, так и исследованиям на физических моделях. Отметим, что, как правило, единичные трещины не представляют практического интереса в свете поиска и описания свойств резервуара. Важным представляется случай, когда трещин множество и их направление и плотность может существенно изменять фильтрационно-емкостные свойства резервуара. В связи с этим важной задачей является разработка континуальных механико-математических моделей таких сред и методов их расчёта для проведения компьютерных экспериментов.
Важнейшей прикладной задачей разведочной геофизики является определение структуры подповерхностного пространства и локализация контрастных границ геологических слоёв. Математически это соответствует алгоритмам построения миграционного изображения, в которых совместно численно решаются прямая и сопряжённая волновые задачи. С увеличением производительности современных многоядерных и многопроцессорных вычислительных систем появляется возможность использовать всё более подробные математические модели для повышения точности восстановления структуры среды. Развитие данного направления может значительно снизить затраты на разработку нефтяных и газовых залежей, обеспечить рентабельность для нетрадиционных углеводородов.
На протяжении длительного времени учёные исследуют процессы, протекающие в земной коре и приводящие к инициации сейсмической активности. Последствиями крупных землетрясений являются значительное разрушение жилых и промышленных сооружений, а также человеческие жертвы. Отметим, что можно выявить два уровня протекающих процессов - неупругие процессы непосредственно в очаге землетрясения и в его окрестности, а также процесс распространения сейсмических волн в дальней зоне. В связи с невозможностью проведения контролируемого полномасштабного экспериментального исследования данных процессов, развитие методов их математического моделирования является актуальной задачей.
Степень разработанности темы исследования
Компьютерное моделирование широко применяется в области сейсмической разведки месторождений нефти и газа. Распространение сейсмических волн в геологическом массиве описывается гиперболической системой уравнений в частных производных. Множество используемых численных методов можно разделить по следующему признаку: в частотной или во временной области рассматривается определяющая система уравнений. В настоящем обзоре отметим основные численные методы, используемые для описания процесса распространения сейсмических волн во временной области.
Хорошие обзоры численных методов, применяющихся для решения уравнений гиперболического типа представлены в работах [Рождественский и др. 1968; Магомедов и др. 1988; Куликовский и др. 2001]. В работе [Riemann 1860] были исследованы характеристические свойства системы, описывающей поведение баротропного газа. В дальнейшем множество работ было посвящено исследованию свойств уравнений газовой динамики и общих систем уравнений гиперболического типа [Русанов 1953; Ferry 1954; Петровский 1961; Курант 1964;
Тихонов и др. 1966; Рождественский и др. 1968; Годунов 1971; Lax 1972; Годунов и др. 1998; Куликовский и др. 1998].
Наиболее известным является семейство конечно-разностных методов, в основе которых лежит аппроксимация производных искомых функций по времени и по пространству конечными разностями на заданном шаблоне схемы. Наиболее распространены на практике схемы с аппроксимацией повремени второго порядка, и по пространству - четвёртого порядка. Это связано с тем, что для уменьшения численной дисперсии при вычислительно осуществимых расчётах на больших пространственных сетках необходимо использование аппроксимации пространственных производных с порядком не ниже четвёртого. В частности, характерными параметрами детальности сетки можно считать требование наличия не менее 10 узлов на минимальную длину волны при дискретизации со вторым порядком, и не мене 3 узлов на минимальную длину волны при дискретизации с четвёртым порядком. В работе [Lines et al. 1999] представлена общая формула для оценки области устойчивости конечно-разностной схемы произвольного порядка аппроксимации по времени и пространству для акустического уравнения и случая постоянной плотности среды. Для изотропной линейно-упругой среды схемы на сдвинутых сетках полностью соответствуют структуре уравнений в переменных скорости-напряжения. При этом возможно ставить задачи по обеспечению сохранения повышенного порядка выбранной схемы для всех пространственных направлений и длин волн. Отметим, что не последнее значение играет и сохранение достаточной скорости расчёта, чтобы решать практически значимые задачи. Одним из последних достижений стало предложение конечно-разностных схем с оптимальными свойствами для акустических волн [Chen et al. 2020b] и для упругих волн [Zhou et al. 2022]. Сохранение достигнутых свойств для сильно неоднородных сред остаётся вопросом, требующим систематического исследования. Кроме того, значительной проблемой является корректный учёт рельефа дневной поверхности. В этом направлении успехи были достигнуты с использованием схемы Лебедева [Puente et al. 2014] и совмещённых сеток (работы научных групп Wei Zhang и Xiaofei Chen). Наиболее распространёнными являются три подхода: 1) сохранение шаблона схемы в узлах вблизи свободной границы и введение дополнительных фиктивных узлов; 2) изменение шаблона схемы вблизи границы без введения фиктивных узлов; 3) использование других, не конечно-разностных, алгоритмов вблизи границы и сшивка с конечно-разностной схемой для основного объёма области интегрирования. В работах [Zhang et al. 2012; Sun et al. 2016, 2018, 2019] предложен подход, основанный на построении фиктивных узлов и совмещённой криволинейной сетки вблизи границы. Ввиду того, что сетка может деформироваться во всех трёх направлениях, возможно точное покрытие достаточно сложного рельефа. Схема устойчива и достаточно точна даже для больших градиентов высот. Метод наложенных сеток был
успешно применён в работах [Zang et al. 2021], который основан на наложении деформированных сеток на основную прямоугольную и переинтерполяции значений искомых функций. Метод погруженных границ также был успешно обобщён в работе [Gao et al. 2015] для конечно-разностных методов на сдвинутых сетках. С изменением шаблона вблизи свободной поверхности [Shragge et al. 2020] была построена новая система уравнений для вертикально деформированных сеток и соответствующие ей граничные условия. В направлении сшивки конечно-разностной схемы с другими численными методами положительные результаты для разрывного метода Галёркина были получены в работе [Lisitsa et al. 2016]. Второй проблемой для конечно-разностного метода является корректный учёт наличия контрастных границ материала внутри одной расчётной ячейки, так называемая проблема "sub-sell resolution". Для решения данной проблемы расчётный метод должен явно зависеть от положения и ориентации границы внутри ячейки. Отметим два различных подхода к этой проблеме. В работах [Mittet et al. 2017, 2021] отмечается, что необходимо согласовывать расчётную сетку материала и рассчитываемое волновое поле. Например, возможно создать мелкую сетку с использованием фильтра низких частот, чтобы удалить волновые числа, не соответствующие основной грубой сетке. Другим подходом является предложенное в работах [Moczo et al. 2014, 2019; Kristek et al. 2017, 2019; Gregor et al. 2021] орторомбическое представление сильно неоднородных упругих, вязкоупругих и пороупругих сред. Его основа -использование схемы на сдвинутых сетках четвёртого порядка аппроксимации по пространству и второго порядка по времени. Данный метод остаётся вычислительно эффективным как для гладкого изменения параметров материала, так и при наличии резких контрастов. Отметим проведённую детальную верификацию относительно результатов, получаемых с использованием метода спектральных элементов, разрывного метода Галёркина, аналитических и полуаналитических решений. Разные авторы проводили детальный анализ устойчивости конечно-разностных схем, например, детально был исследован класс конечно-разностных схем для системы уравнений Био [Alkhimenkov et al. 2020]. Открытыми на сегодня остаются вопросы по созданию корректных и устойчивых алгоритмов, позволяющих проводить моделирование при наличии контакта грубых и мелких пространственных расчётных сеток, а также с локальным измельчением (по пространству и времени) [Kostin et al. 2015; Li et al. 2015; Fan et al. 2015; Nie et al. 2017].
Развитию метода конечных элементов и его применению для численного решения задач деформируемого твердого тела посвящено множество работ [Бате 2010; Зенкевич 1975; Оден 1976; Oden 1983; Zienkiewicz 2000]. Использование базисных функций высокого порядка в методе конечных элементов приводит к построению метода спектральных элементов. Теоретическому исследованию данного класса методов посвящён ряд исследовательских работ,
например, [Warburton et al. 1999; Komatitsch et al. 2008; Capdeville et al. 2003; Bernardi et al. 1997; Kanevsky et al. 2006]. Его основы применительно к рассматриваемым в работе задачам изложены в работе [Moczo et al. 2014]. Он примечателен тем, что объединяет геометрическую мощь метода конечных элементов и точность спектральных методов. Однако, данный метод не лишён и недостатков. Например, если границы слоев не могут быть явно выделены на этапе построения расчётной сетки, возможно значительное падение точности расчётов [Chaljub et al. 2015]. С практической точки зрения основные проблемы возникают при попытке описать распространение поверхностных волн, рассеивающихся на сильных горизонтальных неоднородностях (края долин и бассейнов). Для корректного представления мелких неоднородностей (мельче размера одного элемента сетки), необходимо построение осреднённых (гомогенных) моделей эффективных сред. Для этого может быть использован подход Two-Scale Homogenization (TSH), который приводит к полностью анизотропной эффективной среде, определяющие уравнения которой могут быть успешно решены методом спектральных элементов [Capdeville et al. 2020]. Для таких сред могут использоваться полиномы высокого порядка - вплоть до 40-го. В том случае, если геометрия расчётной области не может быть с достаточной точностью покрыта структурной расчётной сеткой, либо эффективная гомогенная модель не может быть построена, метод спектральных элементов допускает эффективное комбинирование с методом конечных элементов, с разрывным методом Галёркина [Brun et al. 2021]. Известны примеры объединения различных методов для описания водных включений со сложной геометрией и нелинейным поведением [Terrana et al. 2018; Brissaud et al. 2017]. Преимуществом этого является возможность явно выделить границу между твердым телом и флюидом, что не получается сделать при применении методики гомогенизации. Отметим, что широкую популярность приобрёл метод спектральных элементов в динамических задачах моделирования процесса землетрясения.
Широко применяется в сейсмических задачах и разрывный метод Галёркина, изначально предложенный для численного решения задачи транспорта нейтронов с использованием методов интегрирования по времени Рунге-Кутты высоких порядков. В дальнейшем он был успешно обобщён на гиперболические системы общего вида [Hesthaven et al. 2008]. В работе [Галеркин 1915] данный метод был подробно описан, тогда как ранее уже применялся для решения задач упругости в работе [Бубнов 1913]. В виду этого он иногда называется методом Бубнова-Галёркина. Его подробное теоретическое обоснование было представлено в работе [Келдыш 1942]. Благодаря локальности используемого пространственного шаблона, он позволяет использовать криволинейные и неструктурные сетки, состоящие из треугольников и тетраэдров, что облегчает включение в расчёт сложных геологических структур и топографии дневной поверхности, такой как вулканы, осадочные бассейны, разломы и трещиноватые зоны,
границы раздела геологических слоёв и контрастные внутренние границы [Mercerat et al. 2015; Gabriel et al. 2020]. Отличительной чертой данного метода является использование численных потоков через границы элементов, что не требует какой-либо непрерывности искомых полей и позволяет естественно работать с нелинейными контактными условиями, например, в задачах глобальной сейсмики [Pelties et al. 2012]. Для волновых задач преимуществом является низкая численная дисперсия и величина диссипации, зависящая от размера ячеек и степени используемых полиномов. Объединение разрывного метода Галёркина с подходом Arbitrary high-order DERivative (ADER) для интегрирования по времени позволило построить высокоточные одношаговые схемы. Развитием данного подхода является обобщение на нелинейные задачи с использованием апостериорных лимитеров и конечно-объёмных схем [Reinarz et al. 2019]. Несмотря на высокую вычислительную сложность разрывного метода Галёркина (по сравнению с конечно-разностными методами), активное использование высокопроизводительных вычислительных систем сделало его применимым для решения реальных задач [Wilcox et al. 2010; Breuer et al. 2014; Heinecke et al. 2014]. Этого позволили добиться следующие особенности метода: простота и эффективность написания локальных схем интегрирования по времени [Uphoff et al. 2017]; простое внедрение метода адаптивного измельчения расчетной сетки [Burstedde et al. 2011]; низкоуровневые оптимизации для увеличения скорости расчётов [Uphoff et al. 2020]. Методы ADER-DG были успешно обобщены на различные реологии сред: вязкоупругая и анизотропная [Uphoff et al. 2016; Wolf et al. 2020].
Отдельно необходимо упомянуть методы, основывающиеся на теореме Сомильяна (C. Somigliana), связывающей поле смещений внутри рассматриваемой области с интегралами по её границе от смещений и напряжений. К ним относятся метод граничного элемента и непрямой метод граничного элемента. Первый из них следует из теоремы Сомильяна и основывается на взаимности. Второй метод является следствием линейности задачи, поскольку поле состоит из суперпозиции распределенных нагрузок, которые являются промежуточными неизвестными. Его можно рассматривать как математическую запись принципа Гюйгенса. Несмотря на то, что данные методы могут быть применены во временной и в частотной области, в сейсмологии превалирует использование последней. Дискретизация границ, применение граничных условий и использование преобразования Фурье приводит к расчётным формулам метода [Kawase et al. 1988; Sanchez-Sesma et al. 1991]. В работе [Perton et al. 2016] представлена реализация непрямого метода граничного элемента, которая позволяет моделировать распространение упругих волн в сложных областях, состоящих из однородных подобластей с нерегулярными границами или с плоскими слоями. Несмотря на элегантность методов граничного элемента, область их применения ограничена, поскольку функция Грина может быть эффективно вычислена только для однородных сред или для сред с постоянным градиентом. Таким образом,
моделирование реалистичных геометрий и реологий выходит за рамки возможностей методов. Тем не менее, универсальность и внутренняя экономичность методов граничных элементов делают их полезными для тестирования и проверки различных численных методов.
Для задач, решения которых могут содержать разрывы, успешно применяются схемы ENO (essentially non-oscillatory), изначально предложенные в классической работе [Harten et al. 1987]. Конечно-разностные и конечно-объёмные схемы основываются на интерполяции дискретных данных полиномиальными или другими простыми функциями. С возрастанием длины используемого пространственного шаблона возрастает и порядок аппроксимации построенной схемы. Стандартные конечно-разностные схемы используют фиксированный пространственный шаблон. Данный метод хорошо работает на гладких решениях, однако, при наличии разрывов решения линейные схемы выше первого порядка аппроксимации являются немонотонными и порождают осцилляции на разрыве. Для борьбы с данным эффектом возможно, например, вводить введение в схему искусственной вязкости. Недостатком данного подхода является необходимость настройки параметра вязкости под конкретную расчётную задачу, чтобы достаточно погасить осцилляции вблизи разрыва и незначительно снизить точность расчёта в гладких областях. Другим подходом является использование в схеме ограничителей (лимитеров). Например, в областях с большими градиентами возможен переход на монотонную схему первого порядка аппроксимации. Аккуратное построение лимитера может обеспечить построение TVD схемы, в том числе для нелинейных одномерных задач. Недостатком данных схем является падение точности расчёта вблизи гладких экстремумов решения. Достаточно подробно два указанных подхода рассмотрены в работах [Sod 1985; LeVeque 1990]. Для случая, когда решения задачи является кусочно-гладкой функцией, что характерно для гиперболических задач, был предложен другой подход, получивший название ENO [Harten et al. 1987]. Он основывается на выборе локального пространственного шаблона схемы исходя из локального анализа поведения численного решения. Данный подход получил развитие в работах многих авторов. В работах [Shu et al. 1988, 1989] построены робастные TVD схемы с использованием методов Рунге-Кутты. Оптимизация точности и устойчивости данных схем путём выбора процедуры изменения шаблонов исследовалась в работах [Fatemi et al. 1991; Shu 1990]. Идея использования линейных комбинаций схем с различными шаблонами была представлена в работах [Jiang et al. 1996; Liu et al. 1994] и получила название WENO (weighted essentially non-oscillatory). Схемы ENO успешно комбинируются с другими вычислительными методами: multiresolution [Bihari et al. 1995], спектральный метод [Cai et al. 1993]. Альтернативными подходмами к решению проблемы монотонизации расчётных схем являются MUSL-схемы, CABARET-схемы. Однако, данные схемы имеют приблизительно первый порядока локальной и интегральной сходимости в областях влияния ударных волн [Остапенко
1997, 2000; Casper et al. 1998; Engquist et al. 1998]. При этом классические немонотонные схемы, использующие аналитические функции численных потоков, сохраняют повышенный порядок сходимости [Ковыркина и др. 2010]. В работе [Зюзина и др. 2018] предложено комбинировать два данных подхода. Таким образом, изначально во всей области использется немонотонная схема, которая имеет повышенный порядок сходимости в областях влияния ударных волн. В окрестностях больших градиентов построенное решение корректируется путем численного решения внутренних начально-краевых задач по одной из NFC-схем (nonlinear flux correction). При этом выбор областей, в которые требуется проведение корректировки, возможен как простейшим градиентным методом, так и более продвинутыми методами с использованием вейвлет разложения [Arandiga et al. 2008] или Фурье разложения [Gelb et al. 2002].
Компактные (продолженные) численные схемы основаны на использовании дифференциальных [Grudnitskii 1977; Aoki 1991] или интегральных [Rogov 2012; Аристова и др. 2012] следствий исходной системы уравнений и позволяют получить повышенный порядок аппроксимации исходной дифференциальной задачи без расширения пространственного шаблона. Особенно актуально данное их свойство при решении линейных гиперболических задач с разрывными коэффициентами. Использование двухточечного по пространству шаблона обеспечивает возможность для произвольно заданной геологической модели, расположив узлы сетки в местах разрыва коэффициентов, решать задачу интерполяции гарантированно для непрерывных искомых функций. Схемы данного типа активно развивались в работах отечественных учёных А.И. Толстых [Tolstykh 1990], Б.В. Рогова [Rogov 2012], Е.Н. Аристовой [Аристова и др. 2012], М.Н. Михайловской [Михайловская и др. 2012] и зарубежными специалистами Yabe T. [Yabe et al. 1991; Fukumitsu et al. 2015; Yabe et al. 2004], Ito K. [Ito et al. 2014] Компактные сеточно-характеристические продолженные схемы для одномерного и многомерного линейного уравнения переноса с постоянными коэффициентами были предложены в работах Хохлова Н.И., Голубева В.И. и Петрова И.Б. [Голубев и др. 2016a]. Из последних достижений в данной области отметим построение бикомпактной схемы второго-четвёртого порядка для нестационарных уравнений Навье-Стокса [Брагин и др. 2021b], обобщение схем на линейное многомерное уравнение конвекции-диффузии [Брагин и др. 2021a].
По-видимому, прямой метод характеристик был изначально предложен в работе [Massau 1899]. В дальнейшем он был обобщен на случай двумерных [Жуков 1960] и трехмерных [Butler 1960; Магомедов 1966] задач. Основным его недостатком является сгущение расчетных узлов в отдельных областях расчетной области, что приводит к снижению точности расчета вдали от них. В дальнейшем в работах в [Магомедов и др. 1969; Кукуджанов 1976] был предложен метод обратных характеристик (сеточно-характеристический метод) на фиксированных расчетных
сетках. Вышеперечисленные методы успешно применялись для получения гладких решений задач газовой динамики. Для обеспечения возможности корректного расчета разрывных решений были предложены сеточно-характеристические методы с явным выделением разрывов [Шевелев 1986; Холодов 2008]). Одним из важных свойств численных методов является монотонность [Fridrichs 1954]. К монотонным схемам относится, например, схема первого порядка аппроксимации [Холодов 1978]. В [Федоренко 1962, Петров и др. 1984] предложены монотонные гибридные схемы повышенного порядка аппроксимации. Общий подход к построению таких схем с использованием пространства неопределенных коэффициентов был предложен в [Холодов и др. 1980]. Также для повышения порядка сходимости без расширения расчетного шаблона в [Холодов и др. 2006] предложена монотонная схема с использованием продолженных систем. В работе [Холодов Я.А. и др. 2018] на основе сеточно-характеристического критерия монотоности предложен универсальный алгоритм построения монотонных сеточно-характеристических схем. Отметим, что существуют и другие квазимонотонные схемы [Вязников и др. 1989; Головизнин и др. 1998]. Изначально сеточно-характеристический метод был применен для задач газовой динамики, а к задаче описания динамических процессов в сплошных средах он адаптирован сравнительно недавно. Среди первых работ, посвященных его применению для численного решения многомерных уравнений динамики упругих, а также упругопластических и вязкоупругих, сред можно отметить работы [Петров и др. 1984Ь; Кондауров и др. 1978]. В работе [Петров и др. 1990] с использованием сеточно-характеристического метода решались трехмерные волновые задачи теории упругости. Работа [Петров и др. 1990], по-видимому, является одной из первых отечественных работ, в которой рассматривались волновые задачи геофизики. В [Иванов и др. 1990] сеточно-характеристический метод обобщен на случай численного решения вязкоупругих волновых задач, в [Петров и др. 1989] - для исследования сложных волновых процессов в средах многослойной структуры. В дальнейшем разработаны методы до пятого порядка аппроксимации, как на структурных [Голубев и др. 2013], так и на тетраэдральных расчетных сетках для трехмерных задач [Петров и др. 2013]. В [Левянт и др. 2005] сеточно-характеристический метод на треугольных сетках расширен на случай моделирования геологической среды с явным выделением неоднородностей (в том числе трещин) и измельчения расчетной сетки вблизи них для повышения точности расчета. В работе [Петров и др. 2014] сеточно-характеристический метод на структурных сетках был успешно применён для численного решения задач железнодорожной безопасности. В работе [Астанин и др. 2016] с его помощью был произведён трёхмерный расчёт сейсмических волн, инициированных падением метеорита. В работе [Петров, 2017] был представлен сеточно-характеристический метод на системах вложенных иерархических сеток, который позволил проводить расчёт областей,
Похожие диссертационные работы по специальности «Другие cпециальности», 00.00.00 шифр ВАК
Численные методы и алгоритмы расчета волновых сейсмических полей в средах с локальными осложняющими факторами2018 год, доктор наук Лисица Вадим Викторович
Численное моделирование волновых процессов в гетерогенных твердых деформируемых средах2011 год, кандидат физико-математических наук Квасов, Игорь Евгеньевич
Численное моделирование динамических процессов в твердых телах на основе схем повышенной точности1998 год, доктор физико-математических наук Богульский, Игорь Олегович
Моделирование деформационных динамических процессов в задачах сейсморазведки2012 год, кандидат физико-математических наук Панкратов, Сергей Александрович
Численное моделирование распространения упругих волн в неоднородных анизотропных и пористых средах2000 год, кандидат физико-математических наук Крючкова, Виктория Валерьевна
Список литературы диссертационного исследования доктор наук Голубев Василий Иванович, 2022 год
ь -
Рис. 8.15. Сейсмограммы нулевых удалений для модели с одной наклонной границей: а) горизонтальная компонента, Ь) вертикальная компонента
Рассмотрим следы на миграционных изображениях от описанных волн. На рисунке 8.16 цвет границы соответствует её источнику, отдельные части изображения чередуются для предотвращения наложения, а ложные границы отмечены цифрами на изображениях по р и б компонентам.
г ггап5№опа1 г б
Рис. 8.16. Положение восстановленных границ на миграционных изображениях: а) акустическое миграционное изображение, Ь) упругое изображение по р-компоненте, с) упругое изображение по Б-компоненте
Поскольку алгоритм акустической миграции содержит лишь одну характерную скорость распространения волн, задержка отражённого сигнала однозначно определяет глубину залегания восстановленной границы: прошедшая и S-волны возвращаются обратно со скоростью с+. То же самое относится и к границам 4, 5 и 6, 7 в р-компоненте изображения упругой миграции. Границы 1 и 2 обусловлены р-волнами, а граница 3 обусловлена s-волнами, распространяющимися в обратном направлении со скоростью прошедшей волны.
На миграционном изображении по s-компоненте р-волны переходят в s-волны, что обусловливает границы 1 и 3; р-волны становятся прошедшими, что обусловливает границу 2; проходящие волны становятся s-волнами, что вызывает границу 4, а s-волны становятся проходящими, что вызывает границы 5 и 6.
8.3. Полноволновая задача миграции
Численные сеточно-характеристические методы на структурных сетках, используемые в настоящей работе, позволяют с высокой точностью моделировать распространение волн в акустических, изотропных и анизотропных линейно упругих, трещиноватых и насыщенных пористых средах. В настоящем разделе изложен подход, позволяющий построить миграционное изображение среды, используя данные, регистрируемые на дневной поверхности, на основе результатов полноволнового моделирования.
Близость приближенного решения т (вектор распределения параметров модели по пространству) обратной задачи к истинной модели среды будем оценивать с использованием функционала невязки в следующем виде:
Х(т) = )/||5(жг, Ь; т) - й(хг(8.49)
Здесь хг - положение г-го приёмника на поверхности, 5 и й - синтетическая и полевая сейсмограммы скорости смещения, суммирование ведется по всем приёмникам в регистрационной сети. В случае, когда полевые данные представляют собой набор сейсмограмм для разных конфигураций источников, можно выбрать й так, что показания одного сейсмометра из й будут суммой показаний этого сейсмометра по всему набору. Заметим, что соответствующие синтетические сейсмограммы могут быть получены в результате моделирования, в котором источник задан как одновременное действие источников, расположенных в местах расположения сейсмоприёмников при полевом исследовании. Это возможно в силу линейности по внешнему возмущению уравнений линейной динамической теории упругости.
Для изотропной упругой среды вариация функционала невязки может быть вычислена с помощью одного из следующих интегралов по всему пространству [Ьио et а1. 2013]:
8Х = ¡(Кр81п р + Кк 81пк + Кр 51п р)аНх, (8.50)
8Х = | {кг51п 7 + КСр81п с+ + КСб81п с5) й3х, (8.51)
где р - плотность среды, к и д - модуль всестороннего сжатия и модуль сдвига, с+ и с5 -скорости продольной и поперечной волн, 2 = рс+ - импеданс, Ка - ядро параметра а.
При этом миграционное изображение может быть сформировано, используя ядро из любых указанных параметров. В настоящей главе для этого использовалось ядро импеданса которое определяется из выражений
К, = Кр +Кк + Кр, (8.52)
Кр (ж) = р(х) I -)у(х, О (8.53)
Кк(ж) = -к(х) ¡[V ■ з}(х,-Ь)][У ■ 5(ж,0] dt, (8.54)
(ж) = -2^(х) I -0: Б(ж, О dt, (8.55)
О = ^ - ^ I, ^: О = ^О,-. (8.56)
Здесь 5 и V - смещение и скорость, и V* - сопряженные величины, I - единичный тензор, и в выражениях (8.9) - (8.12) используются параметры фоновой модели среды.
Для получения сопряженного поля скоростей используется сопряженный источник - в правую уравнения движения добавляется величина
^(х, Ь) = ^^(Хг, -) - й(х,-ь)]6(х - хг). (8.57)
При этом решение данной сопряжённой гиперболической системы также проводится численно сеточно-характеристическим методом. В данной работе для численного интегрирования по времени для получения полей смещений из полей скоростей была использована формула Симпсона:
I Г(№ - |0' Г(т - Н (Г(0) + Г (Т) + 2Щ_1 Г(2кт) + 4Щ=1 Г([2к - 1]т)), (8.58) х(2пт) = ¡02пт у(в)ав - НТ^2{у([к - 1]т) + 4у(кт) + у([к + 1]т)}, (8.59) в([2п + 1]т) = ¡0[2п+1]°у(в)ав - НТк«0,2{у([к - 1]т) + 4у(кт) + у([к + 1]т)},(8.60)
5}(-Т + 2пт) = ¡('+2птV*(в)йв - 012«(;2{^+ (-Т + [к- 1]т) + 4у*(-Т + кт) + V*(-Т + [к + 1]т)}, (8.61)
в*(-Т + [2п + 1]т) = ¡_'+[2п+1]ту}(в)ав - Н 1к«0,2{^+ (-Т + [к - 1]т) + + кт) +
V* (-Т + [к + 1]т)}, (8.62)
где Т - длина записи всей сейсмограммы, а т - шаг интегрирования по времени. Для численного расчёта пространственных производных использовалась формула второго порядка точности
д$Ях) „ (8.63)
Построение миграционного изображения производится в несколько этапов. Сначала для фоновой модели среды решается прямая задача с конфигурацией источников и приемников, соответствующей произведённым наблюдениям. В результате этого получаются синтетические поля и сейсмограммы. Последние вместе с полевыми сейсмограммами используются на втором этапе для моделирования сопряженного поля. Наконец синтетическое и сопряженное поля используются для расчета ядер и ядра импеданса, в частности [Golubev 2019].
Исследовался процесс распространения сейсмических волн в трещиноватой среде в двумерной постановке. Использовались три различные модели, каждая из которых представляла собой однородную прямоугольную область размерами 8000 х 2900 м со следующими упругими характеристиками: скорость распространения продольных волн Ср -4051 м/с, скорость распространения поперечных волн С? - 2272 м/с, плотность р - 2272 кг/м3. Расчетная сетка имела шаг 5 м и содержала миллион узлов. Отличия между моделями заключались в различной геометрии трещиноватых включений. В первой модели содержались газонасыщенные вертикальные трещины: три из них протяженностью 200 м и одна протяженностью 300 м. Во второй модели содержались две флюидонасыщенные горизонтальные трещины протяженностью 300 м, расположенные рядом друг с другом. В третьей модели содержалась одна субвертикальная (отклонение 10 градусов) флюидонасыщенная трещина протяженностью 300 м.
В качестве начального возмущения использовалась группа точечных источников, расположенных на дневной поверхности через каждые 10 м и создающих вертикальную нагрузку вида импульс Рикера с пиковой частотой 30 Гц. Сейсмометры располагались на дневной поверхности через каждые 10 м и фиксировали значения вертикальной и горизонтальной компонент скорости смещения с шагом 0.6 мс на протяжении 3.8 с. По зарегистрированным на них компонентам скорости строились синтетические сейсмограммы.
В дальнейшем для всех трёх моделей были построены миграционные изображения в соответствии с изложенным в данной главе алгоритмом. В качестве фоновой модели использовалась однородная среда с теми же упругими характеристиками [Golubev et а1. 2017с; Go1ubev et а1. 2018d]. Объёмные волновые поля, полученные решением прямой и сопряжённой задачи, сохранялись с шагом 2.5 мс. Полученные миграционные изображения приведены на Рисунках 8.18-8.20.
Рис. 8.18. Миграционное изображение модели с четырьмя вертикальными газонасыщенными трещинами. А - с наложением истинного положения трещин, Б - без наложения
Рис. 8.19. Миграционное изображение модели с двумя горизонтальными флюидонасыщенными трещинами. А - с наложением истинного положения трещин, Б - без наложения
Рис. 8.20. Миграционное изображение модели с наклонной флюидонасыщенной трещиной. А - с наложением истинного положения трещин, Б - без наложения
Анализ миграционных изображений показывает, что края как газонасыщенных, так и флюидонасыщенных трещин восстанавливаются в виде крестообразных артефактов. Эта особенность позволяет восстановить пространственное положение центра трещины и определить ориентацию её плоскости. Необходимо отметить также наличие множества ложных границ на восстановленном изображении, что в целом является свойством всех алгоритмов миграции. Поиск возможностей их устранения является предметом дальнейшего научного исследования.
В работе исследовалось распространение сейсмических волн в трещиноватой слоистой среде в двумерной постановке. За основу модели был принят шестислойный вмещающий массив, мощностные и упругие характеристики которого взяты из работы [Голубев и др. 2016Ь]. Дополнительно на глубине 1650 м располагались две флюидонасыщенные субвертикальные макротрещины протяжённостью 100 м, наклонённые на угол 10 градусов к вертикали. Расчетная сетка имела шаг 5 м и содержала примерно 1 миллион узлов. В качестве начального возмущения использовалась группа точечных источников, расположенных на дневной поверхности через каждые 10 м и создающих вертикальную нагрузку вида импульс Рикера с пиковой частотой 30 Гц. Сейсмометры располагались на дневной поверхности через каждые 10 м и фиксировали значения вертикальной и горизонтальной компонент скорости
смещения с шагом 0.45 мс на протяжении 3.3 с. Для подавления нефизичных отражений от боковых и нижней границ модели использовались расширения фоновой модели.
На Рисунке 8.21 приведены миграционные изображения, полученные с использованием описанного в работе подхода. Границы раздела слоёв на них отсутствуют, поскольку фоновая модель была изначально слоистой. На изображении при этом отчётливо выделяются места расположения двух трещин (рисунок 8.21а). При ближайшем рассмотрении миграционного изображения можно убедиться в том, что ориентация трещин в точности соответствует истинной. В дальнейшем правая трещина была добавлена в фоновую модель, и также рассчитано миграционное изображение (рисунок 8.21б). В данном случае оно не содержит особенностей в месте расположения правой трещины, что подтверждает корректность работы и возможность использования трещиноватых фоновых моделей.
Рис. 8.21. Миграционное изображение слоистого геологического массива, содержащего две флюидонасыщенные субвертикальные трещины. В качестве фоновой среды использовалась слоистая модель без трещин (а) и с одной трещиной (б).
Как видно, представленный алгоритм позволяет оценивать положение линейных неоднородностей в геологическом массиве. Однако, при проведении полевой сейсмической разведки и последующей обработки экспериментальных данных, положения отражающих горизонтов не известны точно. Как правило, результатом является сейсмический куб параметров, что вносит ошибки в положения границ слоёв. Было изучено их влияние на рассмотренный в работе алгоритм.
Фоновая модель представляла собой слоистый массив, характеристики которого представлены в таблице 8.1. Нумерация слоёв в модели идёт с возрастанием глубины их залегания [Golubev et а1. 2017Ь, Golubev et а1. 2020е]. Значение плотности принималось постоянным и равным 2500 кг/м3. На Рисунке 8.22 представлен общий вид расчётной области. На каждый геологический слой строилась отдельная структурная сетка, между которыми использовались условия полного слипания. В частности, разбиение последнего слоя на два с одинаковыми параметрами были использовано для иллюстрации корректности данного подхода. На глубине 1650 м располагались две флюидонасыщенные трещины протяжённостью 100 м и отклоненные на угол 10 градусов от вертикали. Для этого в соответствующем слое строилась криволинейная расчётная сетка, имеющая шаг порядка 5 м. В качестве источника сигнала использовался набор точечных источников на дневной поверхности, имеющий временную зависимость в виде импульса Рикера с частотой 30 Гц. Источники и приёмники в модели располагались с шагом 10 м. При этом регистрировались все три компоненты вектора скорости. Физическое время численного эксперимента составлял 3.3 с, шаг по времени составлял 0.45 мс.
Табл. 8.1. Упругие свойства геологических слоёв
Номер слоя Толщина, м Скорость Р- Скорость S-
волны, м/с волны, м/с
1 500 3500 1750
2 1000 4500 2250
3 300 5000 2500
4 100 4000 2000
5 300 5500 2750
6 800 5500 2750
Рис. 8.22. Рассматриваемая модель геологической слоистой среды, дополненная двумя субвертикальными трещинами, заполненными жидкостью (внутри синего слоя). Границы слоёв
изображены черными горизонтальными линиями.
На Рисунке 8.23 представлено миграционное изображение, получаемое при задании в качестве фоновой модели среды исходной слоистой модели. Как видно, в местах расположения трещин проявляются артефакты, позволяющие оценить их пространственное положение. При этом не возникает горизонтальных артефактов на границах раздела геологических слоёв.
'А
Юлпр
8.000е-18
4.25е-18
15е-19
-3,25е-18
Рис. 8.23. Миграционное изображение, построенное с использованием точной слоистой модели среды. В местах расположения трещин наблюдаются артефакты.
Для оценки влияния неточности фоновой модели на миграционное изображение был проведен следующий численный эксперимент. Исходная многослойная модель была размыта с помощью фильтра Гаусса вдоль вертикальной оси с окном 20 м. Полученная модель и использовалась в алгоритме полноволновой миграции. Из-за несоответствия положения отражающий горизонтов полученное миграционное изображение, представленное на Рисунке 8.24, содержит артефакты не только в правильных местах расположения трещин, но и вблизи границы слоев. При этом их амплитуды значительно превышают полезный сигнал от трещин.
Рис. 8.24. Миграционное изображение, построенное с использованием размытой слоистой модели среды. В местах расположения трещин наблюдаются артефакты, однако, они значительно меньше по интенсивности артефактов вблизи контактных границ.
8.4. Применение нейронных сетей
В данном разделе представлены результаты решения задачи локализации трещиноватой среды в сложно построенной геологической модели. Был использован подход на основе глубоких свёрточных нейронных сетей.
Используемая геологическая модель была основана на стандартной упругой двумерной модели Магто^ II и покрывала область 17 х 3.5 км. Она содержит почти 200 тонких геологических слоёв, отличающихся упругими характеристиками. Самый верхний слой в модели соответствует водному объёму и описывается в рамках упругой модели с сильно сниженной скоростью распространения продольных волн. Диапазоны параметров модели были следующими: плотность - от 1010 до 2627 кг/м3, скорость распространения продольных волн -от 1028 до 4700 м/с, скорость распространения поперечных волн - от 1 до 2802 м/с (см. рис.
8.25). Внутрь упругой модели помещался трещиноватый объём размером 65 тысяч м2, динамическое поведение которого описывалось в рамках континуальной модели с предельным сдвиговым напряжением на уровне 10 % от используемого сейсмического сигнала (см. рис.
8.26). Исходные упругие данные имели шаг 1.25 м и были использованы для построения модели на сетке 5 м путём процедуры переинтерполяции. Исходя из условия Куранта, шаг по времени фиксировался на уровне 1 мс. В качестве источника сигнала использовалась приповерхностная продольная волна (импульс Рикера, частота 30 Гц). На рисунке 8.27 представлена синтетическая сейсмограмма, построенная по вертикальной компоненте зарегистрированного на дневной поверхности сигнала. Чётко идентифицируются субгоризонтальные отражённые Р-волны, отклик от дна, а также кратные волны. Высокая сложность модели обуславливает сложность регистрируемой волновой картины и трудность локализации трещиноватого объекта. При этом полное волновое поле слабо отличается от поля, посчитанного для исходной упругой модели Магто^ II, что свидетельствует о низкой амплитуде полезного отклика.
1028 3000 4700
P-wave velocity, m/s
Рис. 8.25. Пространственное распределение скорости продольных волн в модели Marmousi
Г) 5 Г) Г) Г) 1ПППП 15 Г) Г) Г)
Рис. 8.26. Место расположения трещиноватого объекта в модели (чёрный цвет)
Рис. 8.27. Синтетическая сейсмограмма на дневной поверхности модели Marmousi II с трещиноватым объектом
В качестве входных данных для нейросети была сформирована база из 740 образцов (модель среды - синтетическая сейсмограмма). При от образца к образцу изменялось только местоположение трещиноватого объекта. Использовался 341 равноудалённый источник на верхней поверхности расчётной области, на каждом из которых было зарегистрировано по 301 временному отсчёту. Данные синтетические сейсмограммы были распределены случайным образом на обучающую и верификационную выборки в соотношении 70/30. Использовалась архитектура UNet [КоппеЬе^ег et а1. 2015], ранее демонстрировавшая хорошие результаты на подобных задачах упругих/акустических сред. При обучении использовалось до 100 эпох в рамках алгоритма оптимизации Адама [К^та et а1. 2014] с постоянной скоростью обучения 10-4 и мини-батчами размером 10. При обучении использовалась функция потерь MSE.
Наличие верификационной выборки позволило оценить качество предсказаний. В частности, использованная нейросеть способна достаточно хорошо локализовывать часть трещиноватых объектов (см. рис. 8.28). Практически отсутствует шум на восстановленном изображении, размеры и положение объекта выявляется чётко. Однако не все прогнозы
одинаково высокого качества. При нахождении трещиноватой области вблизи контрастных границ, на восстановленном изображении присутствуют артефакты. Результаты данного текста представлены на рисунке 8.29. Размер и расположение объекта восстанавливаются правильно. Изображение также содержит значительный шум в большой области вокруг реального местоположения трещин. Этот шум может быть устранён на этапе постробработки, используя простое отсечение по амплитуде или более сложные методы. Качество прогнозов существенно снижается при увеличении глубины залегания объекта. Полезный отклик значительно уменьшается по амплитуде, а сигнал от надлежащих контактных границ существенно его экранируют. Тем не менее, иногда используемая нейросеть всё же смогла достаточно точно локализовать объект на значительной глубине (порядка 2 км), тогда как ручная интерпретация сейсмограмм в этих случаях вряд ли возможна. Рисунок 8.30 демонстрирует пример восстановленного изображения. В ходе дальнейших расчётов не было получено стабильных результатов по восстановлению зон трещиноватости, расположенных на глубине более 2 км. Типичный результат прогноза для таких случаев содержал только шум без локализованной яркой области.
Рис. 8.28. Пример хорошего предсказания. Верхнее изображение показывает записанный сейсмический отклик. Среднее изображение - реальное расположение трещиноватого объекта. Нижнее изображение - предсказанное расположение трещиноватого объекта.
Sample: seismogram_7000_1000
Data: 301 records, 341 receivers per sample
0 50 100 150 200 250 300
Real media (white depicts area of interest)
o.o 1.0 2.0 3.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.012.0 13.0 14.0 15.0 16.0 17.0
Predicted media
o.o 1.0 2.0 3.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.012.0 13.0 14.0 15.0 16.0 17.0
Рис. 8.29. Пример предсказания по зашумленному сигналу. Изображение шумное, но хорошо видно реальное местоположение зоны трещиноватости.
Sample: seismogram_5800_1900
Data: 301 records, 341 receivers per sample
Real media (white depicts area of interest)
o.o 1.0 2.0 3.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.012.0 13.0 14.0 15.0 16.0 17.0
Predicted media
o.o 1.0 2.0 3.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.012.0 13.0 14.0 15.0 16.0 17.0
Рис. 8.30. Пример предсказания глубоко залегающего трещиноватого объекта. Изображение сильно размыто, но местоположение по-прежнему определяется.
8.5. Выводы
В данной главе изложен подход, позволяющий решать обратную задачу сейсмической разведки - построение миграционного изображения. Он основан на использовании приближения Борна и явном построении оператора решения прямой динамической задачи. Представлены вычислительные формулы для двух различных постановок задач: акустическая среда и изотропная линейно упругая среда. На основе ряда моделей, обладающих различными структурными особенностями, показаны преимущества использования алгоритмов упругой миграции перед алгоритмами акустической миграции. Также представлен подход, позволяющий успешно провести локализацию отдельных геологических трещин на основе прямого сеточно-характеристического солвера, обеспечивающего учёт сейсмического отклика от трещин. Дополнительно продемонстрировано успешное применение современных нейронных сетей к задаче поиска трещиноватых включений в геологическом массиве. Для описания процесса формирования сейсмического отклика от трещиноватого объекта использовалась обсуждаемая в работе континуальная модель среды с плоскостями скольжения и отслоения. Продемонстрированы примеры успешной локализации целевого объекта, обсуждены ограничения разработанной методики.
Глава 9. Задачи глобальной сейсмики
9.1. Введение
В главе представлены численные решения задач об инициации сейсмической активности для случая двумерной и трёхмерной постановки. Для описания упругих волн, распространяющихся из гипоцентра землетрясения в дальней зоне, использованы две различные механико-математические модели: «модель подвижки по разлому» и «модель двойной пары сил». Проведено сопоставление результатов расчётов с сейсмотрассами, полученными с использованием других численных методов. В полной трёхмерной постановке решена задача формирования сейсмических волн в очаге землетрясения и их воздействия на композитный трубопровод, расположенный на дне моря. Получено численное решение задачи о нагружении ледового острова. Представлены результаты применения компактной сеточно-характеристической схемы повышенного порядка точности для расчёта сейсмических волн в неоднородных моделях золоторудного месторождения Сухой Лог и Магшош1. Данный алгоритм позволят подробнее воспроизвести структуру сигналов-откликов от геологических границ.
9.2. моделирование очага землетрясения
Для описания процесса инициации сейсмической активности в геологическом массиве, рассматривая распространение сформировавшихся сейсмических волн на достаточном удалении от него, необходимо использовать модели очага землетрясения - сейсмического источника [Вгеш et а1. 2019]. Наиболее простой и интуитивно понятной из существующих является «модель подвижка по разлому». Она основывается на представлении о массиве как двухслойной среде, содержащей сформировавшийся ранее плоский разлом. Благодаря возрастающим региональным напряжения в некоторый момент времени происходит небольшой проскальзывание блоков вдоль него [Favorskaya et а1. 2017а]. Таким образом, существует некоторый объём упругой среды, в котором одна её часть движется с заданной начальной
скоростью V, а вторая - обратно направленной скоростью —V [Go1ubev et а1. 2018Ь]. На Рисунке 9.1 представлена иллюстрация, поясняющая имеющиеся в модели параметры. В частности, необходимо задание трёх углов для описания ориентации разлома. Отметим, что абсолютная
величина и направление начальной подвижки V могут быть оценены на основе карт сейсмических смещений, построенных на дневной поверхности.
Рис. 9.1. Очаг землетрясения - «модель подвижки по разлому» [Golubev et а1. 2013].
Она была успешно использована в двумерной задаче о распространении сейсмических волн при землетрясении на шельфе. Расчётная модель состояла из четырёх геологических слоёв с горизонтальными границами, упругие параметры и толщины которых задавались согласно данным из работы [Заславский и др. 2008]. Их покрывал водный слой, описываемый в рамках акустической системы уравнений. Общий размер расчётной области составлял 4 х 1.5 км, а пространственный шаг расчётной сетки был равен 5 м. Очаг землетрясения располагался на глубине 1.15 км, а его объём принимался равным 50 х 150 м. Все три угла, определяющие ориентацию плоскости разлома, полагались равными 45 градусов. На Рисунке 9.2 представлено распределение модуля скорости в рассматриваемой модели в фиксированный момент времени, достаточный для формирования всех типов сейсмических волн. Исходные продольные и поперечные волны, рождённые вблизи очага, распространяются к дневной поверхности моря, отражаются от неё и от морского дна. Все они несут в себе не только информацию о гипоцентре землетрясения, но и о структуре геологического массива. Развивающиеся сегодня высокоточные методы прямого моделирования позволяют оценить их раздельное и совместное влияние на регистрируемый на сейсмостанциях сигнал [Favorskaya et а1. 2018а,Ь].
Рис. 9.2. Распределение модуля скорости в двумерной слоистой модели.
Рассмотренная двумерная модель была трансформирована в трёхмерную путём увеличения её размеров до 4 х 4 х 1.5 км. При этом вдоль добавленной оси параметры слоёв оставались постоянными, а размер очага землетрясения увеличился до 50 х 150 х 150 м. В связи с тем, что давление на дневной поверхности акустического слоя в явном виде поддерживалось равным нулю, для отображения волновых процессов на поверхности использовалась карта максимальных смещений. Она представлена на Рисунке 9.3 для всех трёх компонент вектора смещений. Явно видна азимутальная анизотропия сигнала, что связано с наклоном плоскости разлома. Кроме того, вертикальная компонента примерно в десять раз превышает горизонтальные компоненты сигнала.
-2000 -1500 -1000 -500 0 500 1000 1500 2000 X-axis, m
Рис. 9.3. Пространственное распределение максимальных смещений на дневной поверхности: х-компонента (а), у-компонента (b) и z-компонента (с).
Второй распространённой моделью очага землетрясения является «модель двойной пары сил» [Aki et al. 1980; Frankel 1993]. В рассмотрение вводится тензор сейсмического момента, определяемый как
Mtj(0 = M0(t) * (п>dj + rijd>), (9.1)
где n - вектор нормали к плоскости разлома, d вектор смещения (проскальзывания), M0(t) - зависящая от времени величина сейсмического момента. Для задания данного
источника в сеточной модели в фиксированном наборе узлов расчётной сетки задаются зависящие от времени внешние силы. Если индексы точки расположения гипоцентра равны О.,], к), то
(Мхх\
Р±х = ± \Мху ) т (1±1,],к), (9.2)
\МХ2/
(МУх\
р±У = ± 1МУУ ) Ш (и±1,к), (9.3)
\МУ&/
(М&х\
Р±& = ±\М&У ) Ш (1,],к±1). (9.4)
Положение плоскости разлома и расположение векторов сил в «модели двойной пары сил» представлено на Рисунке 9.4.
Рис. 9.4. Очаг землетрясения - «модель двойной пары сил».
По построению модели известно, что сумма всех внешних сил, приложенных к среде, равна нуля, так же, как и суммарный момент сил. Таким образом, данный источник не приводит к движению среды как целого, либо её вращению. При этом использование различных временных зависимостей M0(t) определяет интенсивность генерируемого сейсмического сигнала и его частотный спектр.
Реализация данной модели в программном комплексе была протестирована на результатах расчётов, опубликованных Arthur Frankel для трёхмерной модели San Bernardino Valley [Frankel 1993]. В указанной работе использовалась конечно-разностная схема в смещениях-напряжениях, имеющая четвёртый порядок точности по пространству и второй - по времени. Рассматривалась полностью упругая постановка, с заданием условия свободной границы на верхней плоскости и неотражающих условий на всех остальных границах [Bagaev et al. 2019].
Размеры расчётной области были 37 x 16 x 7 км. Источник сигнала располагался на глубине 6 км, в точке с координатами (20 км, 14 км, 1 км). Временная зависимость тензора сейсмического
(t-1.5)2 ^
момента имела вид M0(t) = е 025 . Нормаль к поверхности разлома имела координаты n = (0,-1,0), вектор скольжения был равен d = (1,0,0). Все параметры соответствовали случаю M5 землетрясения из работы [Frankel 1993]. Использовалась кубическая расчётная сетка с шагом 100 м. Шаг по времени выбирался из условия Куранта и составлял 12 мс. При этом моделировалось 48 с физического времени. В Таблице 9.1 представлены упругие свойства трёхслойной модели. Для проведения сопоставления горизонтальная компонента скорости Vx была записана в точке (20 km, 6.3 km, 0) на дневной поверхности. Схема расстановки источника и приёмника изображена на Рисунке 9.5.
Рис. 9.5. Горизонтальное (справа) и вертикальное (слева) сечения модели. Отмечены положения источника (*) и приёмника (х).
Таблица 9.1. Упругие свойства слоёв
Диапазон глубин, км Скорость P-волны, км/с Скорость S-волны, км/с Плотность, кг/м3
0 - 1 1.1 0.6 2000
1 - 3 3.5 2.0 2600
3 - 7 5.0 2.9 2600
Сопоставление результатов расчётов с данными, представленными в работе [Frankel 1993], приведено на Рисунке 9.6. Хорошо видно, что динамические и кинематические характеристики первых вступлений воспроизведены хорошо. Различия в дальнейшей части сигнала могут быть связаны с различным порядком точности используемых численных схем, и некоторое расхождение, безусловно, было внесено оцифровкой положений источника / приемника по рисункам из работы [Frankel 1993].
Рис. 9.6. Величина Vx, рассчитанная сеточно-характеристическим методом (сплошная линия), и данные из работы [Frankel 1993] (точки). Расстояние от приёмника до гипоцентра составляла 7.7 км.
Было также проведено сопоставление результатов моделирования с данными, полученными Robert W. Graves в работе [Graves 1996] с использованием солвера на сдвинутых сетках. Модель состояла из одного геологического слоя с параметрами: скорость P-волн 4000 м/с, скорость S-волн 2300 м/с и плотность 1800 кг/м3. Размеры расчётной области были 100 x 100 x 25 км. Параметры «модели двойной пары сил» были следующие: амплитуда 5e12, нормаль
плоскости разлома n = (0, — -у-, и вектор d = (0, -у-, -у). Источник располагался в точке (50
км, 45 км, 2.5 км). Временная зависимость представляла собой треугольный импульс с углом наклона 45 градусов и протяжённостью 2 с. Пространственный шаг сетки составлял 250 м, шаг по времени был равен 25 мс. Всего моделировалось 16 с физического времени. Для сопоставления выбиралась вертикальная компонента скорости смещения в точке (55 км, 50 км, 0) на дневной поверхности. Схема расстановки источника и приёмника изображена на Рисунке 9.7.
Рис. 9.7. Горизонтальное (справа) и вертикальное (слева) сечения модели. Отмечены положения источника (*) и приёмника (х).
На рисунке 9.8 показано сравнение результатов, полученных сеточно-характеристическим методом, с данными, приведенными в работе [Graves 1996]. Основной пик хорошо воспроизводится и по времени прихода, и по амплитуде сигнала. Различия в регистрируемых сигналах могут быть связаны с различным порядком точности используемых схем. Также следует отметить, что в работе [Graves 1996] использовались сдвинутые сетки, что приводит к увеличение количества узлов в источнике с 6 до 27.
Рис. 9.8. Величина Vz, рассчитанная сеточно-характеристическим методом (сплошная линия), и данные из работы [Graves 1996] (точки).
9.3. Сейсмические воздействия на объекты
В последнее время были открыты нетрадиционные месторождения, разведка и разработка которых традиционными методами не представляется возможной. К ним относятся: сланцевый газ, сверхглубокие залежи, месторождения Арктического шельфа. При этом, при разработке шельфовых месторождений оптимальным является организация нефтепровода (его прокладка по морскому дну) для трансфера нефти от места добычи к местам хранения и переработки. Ввиду того, что доступ к элементам трубопровода ограничен в период его эксплуатации (наличие навигационного периода, сложность и дороговизна глубоководных работ) к его надёжности предъявляются повышенные требования. Одной из причин аварий трубопровода является естественная и техногенная сейсмическая активность - землетрясения. При этом из очага землетрясения, в близи которого происходят процессы как хрупкого, так и вязкого разрушения, к дневной поверхности распространяется набор сейсмических упругих волн, разрушающий объекты наземной и подземной инфраструктуры месторождения. Одной из возможностей предотвращения природной катастрофы, связанной с разливом нефти, является использование современных полимерных композиционных материалов для изготовления элементов трубопровода повышенной стойкости.
Изучением процесса разрушения композиционных материалов занимается множество научных коллективов по всему миру. Подробные обзоры экспериментального, аналитического и численного исследования полимерных композитов приведены в [Abrate 1991, 2016]. Основные результаты серии международных проектов WWFE (World Wide Failure Exercise), целиком посвященных разработке и сравнению критериев объемного разрушения полимерных композиционных материалов (ПКМ) с армированием длинными волокнами, приведены, например, в [Hinton et al. 2004, 2013]. К настоящему моменту было разработано множество различных критериев разрушения, с той или иной степенью достоверности описывающих натурные эксперименты. Наиболее часто применяемыми являются критерии Хашина, Пака, Цая-Хилла, Цая-Ву и Друкера-Прагера. Подробное описание этих критериев и их сравнение друг с другом и экспериментом приведено в [Beklemysheva et al. 2016]. Эти критерии используются во многих коммерческих расчетных пакетах.
Наиболее опасным типом разрушения является BVID - Barely Visible Impact Damage, слабо видимые повреждения при ударе [Abrate 1991; Lopresto et al. 2013]. Обычно они представляются собой трещины в матрице, в том числе протяженные трещины между слоями композита. Такие повреждения могут появиться даже после слабого удара (с энергией порядка 10-100 Дж) и не заметны невооруженным глазом, но существенно снижают остаточную прочность композита, особенно при испытаниях на сжатие [Caprino 1984; Husman et al. 1975].
В настоящей главе рассматривается задача распространения сейсмических волн из очага землетрясения в шельфовом регионе и разрушения композитного нефтепровода под их воздействием [Beklemhysheva et а1. 2018, 2019, 2021]. Математическая постановка задачи включает в себя систему уравнений теории акустики и линейной упругости с учётом анизотропии материала [Рейт et а1. 2015]. Применён расчётный алгоритм оценки объёма разрушений нефтепровода, основанный на последовательном решении двух задач (распространение сейсмической волны от очага землетрясения к морскому дну и её динамическое воздействие на элемент трубопровода), позволяющий оценить возможность разрушения композитной трубы при землетрясении заданной интенсивности.
Задача о разрушении композитного нефтепровода под действием сейсмического импульса является вычислительно сложной из-за большой разницы в характерных пространственных размерах интересующих нас объектов. Общий размер расчетной области, в которой распространяется сейсмический импульс, составляет несколько километров, а толщина композитной трубы находится в пределах от нескольких миллиметров до нескольких сантиметров. Объединение в одном расчете очага землетрясения и элемента трубопровода потребовало бы существенного локального измельчения сетки, что привело бы к увеличению расчетного времени. При этом трубопровод в силу своих размеров оказывает достаточно слабое воздействие на общую волновую картину, которым можно пренебречь. В этом случае задача разбивается на два отдельных этапа (Рисунок 9.9).
На первом этапе моделируется распространение сейсмического импульса от очага землетрясения к дневной поверхности на шельфе. Модель среды была плоскослоистой, включающей верхний слой жидкости, в полной трёхмерной постановке [Заславский и др. 2008]. Она состояла из пяти слоёв с различными параметрами, соответствующими слагающим их породам. Первый слой описывался в рамках теории акустики и имел плотность 1000 кг/м3, скорость распространения продольных волн 1500 м/с, толщину 250 м. Все остальные слои считались упругими и описывались уравнениями линейной теории упругости. Второй слой (слой осадков) имел плотность 1500 кг/м3, скорость распространения продольных волн 1600 м/с, скорость распространения поперечных волн 60 м/с, толщину 50 м. Третий слой имел плотность 2100 кг/м3, скорость распространения продольных волн 2500 м/с, скорость распространения поперечных волн 1000 м/с, толщину 300 м. Четвёртый слой имел плотность 2500 кг/м3, скорость распространения про-дольных волн 3500 м/с, скорость распространения поперечных волн 1300 м/с, толщину 400 м. Пятый слой (кристаллический фундамент) имел плотность 2500 кг/м3, скорость распространения продольных волн 4000 м/с, скорость распространения поперечных волн 2500 м/с, толщину 500 м. В качестве источника сейсмического сигнала использовалась «модель подвижки по разлому» [Go1ubev et а1. 2012] с
геометрическими размерами источника - 50 х 150 х 150 м, углами наклона плоскости разрыва rake=dip=slip=45 градусов и глубиной залегания - 1150 м. Шаг расчётной сетки составлял 5 м, физическое время моделирования - 1,25 с.
1000
2000
-1000
2000
-2000 -2000
Рис. 9.9. Постановка задачи, общий вид
-100С
2000
По результатам расчетов на первом этапе снимались данные в точке, где предполагалось залегание элемента трубопровода. Эти данные использовались на втором этапе для моделирования разрушения фрагмента композитной трубы.
Для моделирования разрушения композитной трубы был взят объем грунта вокруг короткого прямого отрезка трубы. Радиус трубы - 25 см, толщина - 6 см. грань куба - 1 м. Ось трубы направлена вдоль оси Х, вертикаль - ось Y. Общий вид расчетной области приведен на Рисунке 9.9 (справа). Для моделирования изогнутого анизотропного материала трубы, необходимо в каждом расчетном узле хранить либо рассчитывать на каждом шаге реологические матрицы и ряд внутренних параметров метода, что крайне неэффективно либо по времени, либо по используемой оперативной памяти. В связи с этим, труба была разбита на 16 сегментов, материал каждого из которых (углепластик, свойства материала приведены в [Beklemysheva et а1. 2016]) повернут относительно соседних. Рассматривались два случая -направление волокон вдоль оси Х и перпендикулярно ей. Сейсмическое воздействие моделировалось как граничное условие на нижней стороне расчетной области - для его постановки использовались компоненты тензора напряжений, полученные на первом этапе. Остальные границы считались поглощающими, моделируя эффективно вмещающий массив.
На первом этапе использовался метод на структурных сетках с третьим порядком точности по пространству [Go1ubev et а1. 2015с]. Использовался метод явного выделения контактных границ геологических слоёв, что обеспечивает корректный расчёт не только кинематических характеристик сейсмического сигнала, но и динамических характеристик (амплитуд) [Голубев и др. 2016Ь]. Необходимо отметить, что использованное контактное условие между первым слоем (жидкостью) и вторым слоем (слоем осадков) позволяет с
высокой точностью воспроизвести наблюдения, осуществляемые при шельфовой сейсморазведке.
На втором этапе использовался гибридный метод 1-2 порядка на нерегулярных тетраэдральных сетках, которые позволяют быстро и с хорошей точностью приближать сложную геометрию, для линейно упругого анизотропного материала [Bek1emysheva et а1. 2016]. Применение тетраэдральных сеток позволяет моделировать геометрию любой сложности, что дает возможность быстро и эффективно использовать данный метод для расчета прочности различных участков трубопровода (соединений труб, креплений и т. д.). Используются критерии разрушения композитов на основе определенного набора параметров материала, доступных для непосредственного измерения: критерии Хашина, Пака, Цая-Хилла, Цая-Ву и Другера-Прагера.
Для решения первой части задачи была проведена серия расчётов с возрастающей интенсивностью землетрясения. На Рисунке 9.10 приведены волновые картины, полученные в последовательные моменты времени. Продольные и поперечные волны от источника начинают распространяться к дневной поверхности и проходя дно порождают набор акустических волн разных частот и амплитуд.
Рис. 9.10. Распределение модуля скорости в шельфовой модели. Слева - 0,375 с (волны в третьем слое); справа - 0,875 с (волны в воде).
В работе была проведена калибровка параметров модели источника на основании серии прямых расчёт с оценкой максимальных смещений поверхности дна и использования шкалы Рихтера [Richter 1935]. На Рисунке 9.11 представлено распределение максимального за всё время расчёта вертикального смещения дна для двух расчётов. В первом случае скорость среды
в области разлома задавалась равной 1 м/с, во втором случае - 100 м/с. На основе шкалы Рихтера данные землетрясения соответствуют 1 баллу и 3 баллам.
Для дальнейшей оценки сейсмостойкости нефтепроводов проводилась запись напряжённо-деформированного состояния среды на удалении 2 км от очага землетрясения на глубине 350 м. Рассчитывался полный тензор напряжений и вектор скорости точек среды. В дальнейшем эти данные выступали в качестве граничных значений в процессе моделирования динамического воздействия на элемент нефтепровода.
1500 1000 500
4 °
о
-500 -1000 -1500
_20-2000-1500-1000-500 0 500 1000 1500 2000 -5000-1500-1000-500 0 500 1000 1500 2000
X-axis, m X-axis, m
Рис. 9.11. Распределение maxDz в слое донных осадков шельфовой модели. Слева - при начальной скорости источника 1 м/с (1 балл), справа - 100 м/с (3 балла).
Для решения второй части задачи был проведен ряд расчетов на основе данных, записанных в указанной точке геологического массива, для различной интенсивности землетрясения, различного направления волокон и различных критериев разрушения. Наиболее интересные результаты демонстрирует критерий Хашина, так как он позволяет анализировать механизмы, по которым произошло разрушения, и не требует дополнительных внутренних параметров модели, значения которых отсутствуют в открытой печати.
На Рисунке 9.12 приведены области разрушения для различной силы землетрясения и различных механизмов разрушения в критерии Хашина, волокна направлены параллельно оси трубы. При силе землетрясения в 1 балл разрушений не наблюдалось. При силе землетрясения в 2 балла происходят небольшие расслоения в нижней части трубы. При силе землетрясения в 2.5 балла расслоения происходят по всей стенке трубы, часть матрицы композита растрескивается. При силе землетрясения в 3 балла расслоения и растрескивание матрицы происходят по всему объему, часть волокон в нижней части трубы рвется.
На Рисунке 9.13 приведены области разрушения для различной силы землетрясения и различных механизмов разрушения в критерии Хашина, волокна направлены перпендикулярно оси трубы. При силе землетрясения в 2 балла происходят существенные расслоения в нижней и верхней частях трубы. При силе землетрясения в 2.5 балла расслоения происходят по всей
стенке трубы, часть матрицы композита в нижней части трубы растрескивается. При силе землетрясения в 3 балла расслоения и растрескивание матрицы происходят по всему объему, разрыва волокон не происходит.
Г'
Рис. 9.12. Области разрушения по критерию Хашина. Волокна направлены вдоль оси трубы. Сила землетрясения: слева - 3 балла, в середине 2.5 балла, справа - 2 балла. Механизмы разрушения: снизу - расслоение, в середине - растрескивание матрицы, сверху - разрывы волокон.
Рис. 9.13. Области разрушения по критерию Хашина. Волокна направлены вдоль оси трубы. Сила землетрясения: слева - 3 балла, в середине 2.5 балла, справа - 2 балла. Механизмы разрушения: снизу - расслоение, в середине - растрескивание матрицы, сверху - разрывы волокон.
Анализ результатов численного моделирования свидетельствует об изменении картины разрушения при изменении направления укладки волокон в композите. При укладке волокон параллельно оси трубы расслоение при землетрясении магнитудой 2 балла меньше, чем при укладке волокон перпендикулярно оси, но при землетрясении магнитудой 3 балла происходит разрушение волокон, которого нет при укладке волокон перпендикулярно оси. При магнитуде землетрясения менее 2 баллов разрушений в трубе не наблюдалось.
9.4. Пример использования компактной схемы
Представленная в работе компактная сеточно-характеристическая схема повышенного порядка точности была успешно применена для численного решения задачи распространения сейсмических волн в неоднородном геологическом массиве [Golubev et al. 2022a].
Первой из рассмотренных моделей была модель крупнейшего по запасам золоторудного месторождения в мире Сухой Лог. Месторождение Сухой Лог находится в Восточной Сибири, в 900 км от г. Иркутск. В настоящее время оно считается самым крупным золоторудным месторождением в России. Наиболее подробное геологическое описание можно найти в работах [Буряк и др. 1997; Distler et al. 2004; Wood et al. 2006; Yudovskaya et al. 2016].
На рисунке 9.14 изображены его основные геологические тела. В соответствии с [Wood et al. 2006], тело Хомолхо состоит из углеродистых сланцев и алевролитов (Нижнее Хомолхо, hm1); известняков с прослоями углистых известняков (Среднее Хомолхо, hm2); углеродистых сланцев, филлитов и алевролитов (Верхнее Хомолхо, hm3). Тело Имнях разделено на Нижний Имнях (im1), который состоит преимущественно из известняковых песчаников, сланцев и алевролитов, и Верхний Имнях, который состоит из известняка с прослоями сланцев и доломитов [Wood et al. 2006]. Месторождение расположено в ядре узкой антиклинали южной оконечности, расположенной в висячем крыле надвига южной оконечности. Антиклиналь полого погружается на запад-северо-запад, а ее осевая поверхность падает под углом 20-30° к северо-северо-востоку. Минерализованная зона образует наклонное пластообразное тело, параллельное осевой плоскости складки. Зона основного оруденения образует плиту размерами 2000 х 700 м, падающую по ее короткой оси. Её толщины варьируется в диапазон от 70 до 100 м. Основная рудная зона состоит из вкрапленных пирита и параллельных слоистости пирит-кварцевых прожилков, залегающих в Неопротерозойских черных углеродистых сланцах и алевролитах Хомолхинской свиты (hm).
При построении геологической модели использовалась цифровая модель Сухого Лога, представленная в работе [Malovichko et al. 2019]. Её трёхмерный вид и двумерное сечение представлены на рис. 9.15. Верхняя дневная поверхность считалась плоской. Вся модель представлена кубами с размером 25 м. Исходная модель содержала данные по электропроводности. На основе открытой геофизической и геологической информации они были пересчитаны в данные о плотности и скорости распространения P-волн.
North
V
East
Рис. 9.14. Основные геологические тела, использованные в численном эксперименте. Здесь «im» — формация Имнях, «hm2» — формация Среднего Хомолхо. Нижний Хомолхо (hm1), который здесь не показан, расположен между im и hm2. Верхний Хомолхо (hm3) включен в фоновую модель.
Рис. 9.15. Цифровая модель Сухого Лога, использованная в численном эксперименте. Представлены данные (а): Исходная 3D-модель. Цифры соответствуют номерам в табл. 9.1; (Ь) и (с) - 2D модели плотности и скорости Р-волн, соответственно.
Плотность и электрические свойства различных пород в пределах Сухого Лога были достоверно установлены в ходе разведочной кампании, предпринятой в 1970-х годах, которая включала в себя 847 скважин, 11 км подземных проходок, 100 км траншей и многие десятки тысяч образцов. Напротив, данные о механических свойствах пород немногочисленны. Было проведено ограниченное количество лабораторных измерений скоростей продольных и поперечных волн. Установлено, что для ряда фаций существует корреляция скорости и плотности продольных волн с коэффициентом корреляции, равным 0,72. К сожалению, данные о скорости волн в пределах рудного тела отсутствуют. Однако скорость в сланцах основного рудного тела, вероятно, будет ниже, чем во вмещающих породах. В таблице 9.1 приведены свойства материалов, присвоенные каждому телу, использованному при создании модели.
Таблица 9.1. Параметры материалов всех описываемых геологических тел
# Геологическое тело Плотность, кг/м3 Скорость P-волн, м/с
1 Четвертичные отложения 1600 800
2 Имнях (т), нормальная часть 2750 5500
3 Имнях (т), перевернутая часть 2750 5500
4 Основное рудное тело 2860 4500
5 Нижнее Хомолхо ^т1) 2770 4800
6 Среднее Хомолхо ^т2) 2710 5500
7 Верхнее Хомолхо ^т3) 2740 4900
На основе описанной геологической модели была создана численная модель для расчёта сейсмических волн. Сетка содержала порядка 1 300 000 узлов с пространственным шагом 5 м. Шаг по времени составлял 0,5 мс и ограничивался условием устойчивости Куранта. Источник сейсмического сигнала использовался в виде приповерхностной Р-волны с импульсом Рикера частотой 30 Гц. Всего было выполнено 1500 временных шагов.
Технология OpenMP использовалась для ускорения процедуры моделирования. Подпрограммы StepX и StepY были параллелизованы, так как их выполнение занимает почти все время работы программы. В зависимости от количества доступных вычислительных ядер сетка разбивается на набор непрерывных непересекающихся частей. Каждое ядро вычисляет фиксированное количество строк в функции StepX и столбцов в функции StepY. В базовой версии параллельного алгоритма все выделение памяти выполнялось главным потоком. Ускорение в 10 раз было достигнуто на 64-ядерном процессоре AMD. Следует отметить, что по стандарту в системах с произвольным доступом к памяти (NUMA) время доступа к лично выделенной памяти для потока меньше, чем к выделенной другими потоками. Эта возможность позволила оптимизировать алгоритм за счет выделения памяти по потокам и добавления
дополнительного шага предварительного копирования. Это несколько усложнило алгоритм программирования, однако на той же компьютерной системе было достигнуто ускорение почти в 40 раз.
Преимущества предложенной компактной сеточно-характеристической схемы повышенного порядка точности перед сеточно-характеристической схемой первого порядка отражены в сравнении рассчитанных волновых полей (рис. 9.16 а, Ь). Например, кривизна Р-волны, прошедшей через Нижнее и Среднее Хомолхо, реконструировалась лучше (см. стрелки 1 на рис. 9.16). Отраженная волна малой амплитуды в пределах основного рудного тела (рис. 9.16а, стрелка 2) была практически полностью подавлена численной диссипацией метода первого порядка. Напротив, она легко прослеживается в отклике, рассчитанном по компактной схеме (рис. 9.16Ь, стрелка 2). То же наблюдение справедливо и для объемных волн, распространяющихся в низкоскоростных четвертичных отложениях (стрелки 3). В целом рассеянные волны, отраженные от контактных границ, образуют сложную картину, значительно размытую на рис. 9.16а, но сохранившуюся на рис. 9.16Ь. Поскольку современные алгоритмы полноволновой инверсии используют динамические параметры сейсмических сигналов, применение предложенной компактной схемы в составе прямого решателя позволит повысить точность решения обратных задач.
(a) (b)
Рис. 9.16. Волновые поля, полученные при 2D моделировании по схеме первого порядка (а) и по схеме повышенного порядка (b).
Второй рассмотренной моделью является стандартная модель Marmousi [Versteeg 1994], представленная на рисунке 9.17а. Она служит часто используемым бенчмарком для оценки работы алгоритмов решения прямых и обратных задач сейсмики. Геометрия этой модели основана на профиле через прогиб North Quenguela в бассейне Cuanza. Геометрическая и скоростная модели были созданы для получения комплексных сейсмических данных. Распространение акустических волн в этой модели было успешно смоделировано предложенным численным алгоритмом. Сетка содержала примерно 2 384 000 узлов с
пространственным шагом 5 м. Шаг по времени был установлен равным 0,8 мс для выполнения условия устойчивости Куранта. Использовался точечный источник, заглубленный на 5 м ниже дневной поверхности с вейвлетом Рикера 30 Гц. Всего было выполнено 2500 временных шагов. Полученное волновое поле, наложенное на модель плотности, представлено на рис. 9.17Ь. Результаты расчётов подтверждают возможности численной схемы по разрешению подробной структуры сигналов-откликов в неоднородных моделях.
10000
15000
15000
1028 3000 4700
Р-\л/ауе уе1осНу, т/Б
(а)
(Ь)
Рис. 9.17. Пространственное распределение скорости продольных волн в модели Магтош1 (а) и волновое поле в ней в фиксированный момент времени (Ь).
9.5. Нагрузка ледового образования
Освоение Арктики - одна из актуальных современных проблем, стоящих перед научно-техническим сообществом. Полярные районы России включают в себя несколько архипелагов, но в основном это воды Северного Ледовитого океана. Природные ресурсы этого региона огромны, но их добыча обычными методами, применяемыми в других местах, невозможна из-за климатических ограничений. Специфические для Арктики конструкции и подходы, включая
специальные ледостойкие платформы, являются дорогостоящими по сравнению с их аналогами, используемыми в более умеренном климате. Чтобы сделать разработку большинства месторождений рентабельной, затраты должны быть значительно снижены. Наиболее перспективными направлениями развития, направленными на решение этой проблемы, являются искусственные ледовые сооружения, в том числе ледяные острова. Их можно использовать для добычи нефти и природного газа. Сравнение искусственных ледяных островов с обычными конструкциями можно найти в работе [Горгуца и др. 2017], методы строительства и обслуживания обсуждаются в работе [Morgan et al. 2015]. Технология создания ледообразования была реализована в канадском море Бофорта [Riley 1975], ее преимущества и недостатки исследованы в [Karelov et al. 2015].
Доказано, что конструкции из искусственного льда экономически выгодны. С другой стороны, их строительство и эксплуатация намного сложнее, чем использование традиционных методов, из-за новизны и отсутствия соответствующего опыта. Устойчивость ледяных конструкций следует рассматривать с точки зрения как теплового, так и упругого воздействия. Термические взаимодействия включают плавление; следовательно, необходимо решать задачу Стефана. Задача Стефана — это особый вид параболической системы дифференциальных уравнений в частных производных, используемый для описания эволюции систем, содержащих различные фазы, границы между которыми могут перемещаться во времени [Олейник 1960; Meirmanov 2011]. Механические нагрузки, вызванные ударными воздействиями при бурении, а также давление со стороны техногенных сооружений описываются уравнениями механики сплошной среды.
В настоящем разделе исследуется распределение напряжений и распространение упругих волн в ледовом образовании [Petrov et al. 2021b]. На первом этапе вычисляется распределение температуры после его построения и эволюции в течение года. Затем эти результаты, в том числе форма острова и поле температур, используются при расчете упругих процессов. Эти данные могут быть использованы в качестве верхней оценки возможных нагрузок на остров за каждый месяц и ограничивают период времени безопасной эксплуатации острова.
Коллегами (Муратов М.В., Конов Д.С.) с использованием неявной схемы Письмена-Рекфорда была численно решена задача Стефана для ледового образования заданной формы. Результатами расчёта являлись поля температур и геометрия острова (граница раздела фаз лёд-вода). Для описания динамического поведения льда под нагрузкой с учётом изменения температуры по пространству использовалась вязкоупругая модель Максвелла [Argatov 2012] и формула Берденникова [Берденников 1948]:
Е(Т) = 108(87.6 - 0.21Г - 0.0017Г#), (9.5)
а коэффициент Пуассона считался неизменным и равным величине 0,295. Зоны повреждения во льду определялись на основе критерия Мизеса. В модели использовалось значение 2.2 * 106 Па для предельного напряжения сдвига.
На рисунке 9.18 представлена геометрия расчётной области в задаче Стефана. Цифрами обозначены различные среды: 1 - лёд, 2 - подстилающий массив, 3 - вода, 4 - воздух. Для численного моделирования в качестве граничных условий использовались экологические данные Баренцева моря, в том числе средние температуры воды и воздуха по месяцам. На рисунке 9.19 представлено распределение температуры в последовательные моменты времени (по месяцам). Переход через ноль градусов соответствует границе острова.
® 300 m
ф Ют
<3> 8 т
(2) : 5 т
Рис. 9.18. Геометрия расчётной области в задаче Стефана.
О 20 40 60 80 100 120 140 160 0 20 40 60 80 100 120 140 160
Temperature distribution in the end of july Temperature distribution in the end of august
Рис. 9.19. Поля температур в течение года по месяцам.
На рисунке 9.20 представлен вид итоговой расчётной области в задаче о динамическом нагружении ледового острова. Параметры, использованные в расчёте, представлены в таблице 9.2. На верхней границе использовалось условие свободной поверхности. Дополнительно на нижней и боковой границах задавались неотражающие граничные условия. К верхней границе ледяного образования прикладывалась постоянная вертикально направленная сила с величиной контактной площадки 5 м. Расчёт вёлся на установление [Golubev et а1. 2022Ь].
Таблица 9.2. Параметры для задачи о нагружении ледового острова
# Материал Плотность, кг/м3 Скорость Р-волн, м/с Скорость S-волн, м/с
2 Вода 1000 1500 1
3 Дно 2500 1806 316
4 Осадочные породы 2500 2250 1000
Рис. 9.20. Геометрия расчётной области в задаче о нагружении ледового острова.
Для точного описания геометрии ледового образования использовался набор из структурных криволинейных сеток. Пространственный шаг по вертикальной оси находился в пределах от 0,04984 м до 0,04985 м. В средней части пространственный шаг по горизонтальной оси был равен 1 м. В левой и правой частях - в пределах от 0,35 м до 0,8 м вверху и от 1 м до 1,2 м внизу и посередине (см. рисунок 9.21). Водная область состояла из двух частей: криволинейной сетки, ближайшей к кромке льда, и прямоугольных сеток. Эти сетки имели такой же пространственный шаг по вертикальной оси, как и во льду. Пространственный шаг по горизонтальной оси был примерно равен 1 м. Грунт состоял из двух слоев: водного дна и осадочных пород. Его сетка имела прямоугольную форму и была разделена на несколько частей, чтобы обеспечить такой же пространственный шаг по горизонтальной оси, как и в областях выше. Пространственный шаг по вертикальной оси был равен 1 м. Поле температуры интерполировалось из исходной неструктурной сетки. Динамическая упругая задача решалась с фиксированным шагом по времени: т = 0,01 мс.
August
10 8 6 4 2 0
О 50 100 150 200 250 300
Рис. 9.21. Ледяное образование, покрытое расчётной сеткой. Её криволинейная часть темно-серого цвета, прямоугольная часть светло-серого цвета.
Было рассчитано влияние оборудования для геофизических исследований на ледообразование весом 1000 кг. На рисунке 9.22 представлено волновое поле через короткое время после начала расчета. На рисунке 9.23 представлено волновое поле спустя долгое время после начала расчета. Эти эксперименты показывают, что процесс плавления, заключающийся в распределении температуры и изменении геометрической формы острова, не оказывает существенного влияния на наблюдаемые волновые поля. Было установлено, что ледовый остров не разрушается под действующей нагрузкой.
vnoim
00е+00 0.001 0.002 0.003 4 5е-03
Рис. 9.22. Представлен модуль скорости через короткое время после начала расчета. Для января - вверху, для августа - внизу.
у_по»т
0.0е+00 0.001 0.002 О.ООЗ 4.5еОЗ
Рис. 9.23. Представлен модуль скорости через длительное время после начала расчета. Для января - вверху, для августа - внизу.
9.6. Выводы
В главе представлены численные решения задач об инициации сейсмической активности для случая двумерной и трёхмерной постановки задач. Для описания упругих волн, распространяющихся из гипоцентра землетрясения в дальней зоне, использованы две различные механико-математические модели: «модель подвижки по разлому» и «модель двойной пары сил». Проведено сопоставление результатов расчётов с сейсмотрассами, полученными с использованием других численных методов. Выявлено, что динамические и кинематические характеристики сигналов, как правило, воспроизводятся хорошо. Различия в части сигналов могут быть связаны с различным порядком точности используемых численных схем, и неточностями оцифровки положений источника / приемника.
В главе в полной трёхмерной постановке решена задача формирования сейсмических волн в очаге землетрясения и их воздействия на композитный трубопровод, расположенный на дне моря. Анализ результатов численного моделирования свидетельствует об изменении картины разрушения при изменении направления укладки волокон в композите. Получено численное решение задачи о нагружении ледового острова. Представлены результаты успешного решения прямой задачи распространения сейсмических волн в неоднородных моделях золоторудного месторождения Сухой Лог и Магтошь Использование компактной сеточно-характеристической схемы повышенного порядка точности позволят подробнее воспроизвести структуру сигналов-откликов от геологических границ.
Заключение
Основные результаты работы:
1. Разработаны сеточно-характеристические методы на структурных расчётных сетках повышенного порядка аппроксимации с использованием расширенного пространственного шаблона и компактного пространственного шаблона совместно с дифференциальными следствиями исходной системы уравнений. Обеспечено сохранение повышенного порядка аппроксимации многомерной схемы при использовании полиномиальной интерполяции заданного порядка на этапе решения линейных одномерных уравнений переноса путём применения многостадийных методов расщепления по пространственным направлениям. Выполнена адаптация сеточно-характеристических методов для моделирования процесса динамического нагружения пористых флюидонасыщенных сред.
2. Изучена континуальная модель линейно-упругой среды, содержащей плоскости возможного скольжения. Построены явно-неявные расчётные схемы для численного решения полулинейных гиперболических систем уравнений, описывающих динамическое поведение сред данного типа. Обеспечено сохранение повышенного порядка расчётной схемы с опорой на использование сеточно-характеристических методов для построения численного решения линейно-упругой задачи.
3. Получены выражения для построения миграционных изображений в рамках приближения Борна для акустических и линейно-упругих сред, содержащих слабо контрастные неоднородности. Применён сеточно-характеристический метод на структурных сетках для расчёта миграционного изображения в линейно-упругих трещиноватых средах. Обоснована возможность включения трещиноватых объектов в фоновую модель геологической среды.
4. Изучены волновые процессы, протекающие в пористых флюидонасыщенных средах под влиянием динамической нагрузки. Исследованы сигналы-отклики, генерируемые на границе: акустическая среда - пористая среда и линейно-упругая - пористая среда. Обоснована возможность скважинной диагностики наличия зон повышенной пористости и проницаемости с использованием высокочастотных источников. Выполнен расчёт процесса распространения сейсмических волн в анизотропных неоднородных средах. Исследован процесс взаимодействия продольных и поперечных волн с заглублёнными трещиноватыми структурами, описываемыми как в рамках континуальных моделей, так и с помощью методов явного выделения границ отдельных трещин. Обоснованы преимущества применения методов упругой миграции перед
методами акустической миграции. Проведено моделирование процесса распространения сейсмических волн из очага землетрясения и их взаимодействия с инфраструктурными объектами.
Автор выражает благодарность коллективам кафедры информатики и вычислительной математики, кафедры вычислительной физики, лаборатории прикладной вычислительной геофизики МФТИ и в особенности научному консультанту Петрову И. Б. за поддержку научных исследований и плодотворное обсуждение научных результатов.
Список литературы
Алексеев, А.Д. Природные резервуары нефти в отложениях баженовской свиты на западе широтного приобья: дис. канд. геолого-минералогических наук: 25.00.12. - ЗАО «Моделирование и мониторинг геологических объектов им. В.А. Двуреченского», Москва, 2009 - 185 с.
Алхименков, Ю.А., Баюк, И.О. (2013) Границы применимости параметров Томсена для трещиноватого карбонатного коллектора // Технологии сейсморазведки. - 2013. - № 4. -С. 36 - 48.
Аристова, Е.Н., Рогов, Б.В. (2012) О реализации граничных условий в бикомпактных схемах для линейного уравнения переноса // Матем. моделирование. - 2012. - Т. 24, № 10. - С. 3
- 14.
Астанин, А. В., Дашкевич, А. Д., Петров, И. Б., Петров, М. Н., Утюжников, С. В., Хохлов, Н. И. (2016) Моделирование влияния головной ударной волны Челябинского метеорита на поверхность Земли // Матем. моделирование. - 2016. - Т. 28, № 8. - С. 33 - 45.
Бате, К. Ю. (2010) Методы конечных элементов. - М.: ФИЗМАТЛИТ, 2010. - 1022 с.
Берденников, В.П. (1948) Изучение модуля упругости льда // Труды ГГИ. - 1948. - Вып. 7 (61).
- С. 13 - 23.
Брагин, М.Д., Рогов, Б.В. (2021a) Бикомпактные схемы для многомерного уравнения конвекции-диффузии // Ж. вычисл. матем. и матем. физ. - 2021. - Т. 61, № 4. - С. 625 -643.
Брагин, М.Д., Рогов, Б.В. (2021b) О точности бикомпактных схем в задаче о распаде вихря Тейлора-Грина // Ж. вычисл. матем. и матем. физ. - 2021. - Т. 61, № 11. - С. 1759 - 1778.
Брадучан, Ю.В., Гурафи, Ф.Г., Захаров, В.А. и др. (1986) Баженовский горизонт Западной Сибири (стратиграфия, палеография, экосистема, нефтеносность). - М.: Наука, 1986. -216 с.
Бреховских, Л.М. (1973) Волны в слоистых средах. - М.: Наука, 1973. - 502с.
Бубнов, И.Г. (1913) Отзыв о сочинениях проф. Тимошенко, удостоенных премии Д.И.Журавского // Сб. Института инженеров путей сообщения. - 1913 - Вып. 81, Петербург.
Бураго, Н.Г., Глушко, А.И., Ковшов, А.Н. (2000) Термодинамический метод получения определяющих уравнений для моделей сплошных сред // Изв. РАН. МТТ. - 2000. - № 6. - С. 4 - 15.
Бураго, Н.Г., Ковшов, А.Н. (2001) Модель дилатирующей разрушающейся среды // Изв. РАН. МТТ. - 2001. - № 5. - С. 112 - 117.
Бураго, Н.Г., Кукуджанов, В.Н. (2004) Численное решение задач континуального разрушения. -М.: Изд-е ИПМ РАН, 2004. - 40с.
Буряк, В.А., Хмелевская, Н.М. (1997) Сухой Лог - одно из крупнейших золоторудных месторождений мира (генезис, закономерности размещения оруденения, критерии прогнозирования). - Владивосток: Дальнаука, 1997. - 156 с.
Войнов, О.Я., Голубев, В.И., Жданов, М.С., Петров, И.Б. (2016) Построение метода упругой миграции сейсмических данных в приближении Борна // Труды МФТИ. - 2016. - Т. 8, № 2 (30). - С. 60-66.
Вязников К.В., Тишкин В.Ф., Фаворский А.П. (1989) Построение монотонных разностных схем повышенного порядка аппроксимации для систем уравнений гиперболического типа // Математ. моделирование. - 1989. - Т. 1, № 5. - С.95 - 120.
Галеркин, Б.Г. (1915) Стержни и пластинки. Ряды в некоторых вопросах упругого равновесия стержней и пластинок // Вестник инженеров. - 1915. - Т. 1. - С. 897 - 908.
Годунов, С.К. (1971) Уравнения математической физики. - М.: Наука, 1971. - 416 с.
Годунов, С.К., Роменский, Е.И. (1998) Элементы механики сплошных сред и законы сохранения. - Новосибирск: Научная книга, 1998. - 280 с.
Головизнин, В.М., Карабасов, С.А. (1998) Нелинейная коррекция схемы Кабаре // Матем. моделирование. - 1998. - Т. 10, № 12. - С. 107 - 123.
Голубев, В. И., Петров, И. Б., Хохлов, Н. И. (2016a) Компактные сеточно-характеристические схемы повышенного порядка точности для трёхмерного линейного уравнения переноса // Матем. моделирование. - 2016. - Т. 28, № 2. - С. 123 - 132.
Голубев, В.И., Петров, И.Б. (2016Ь) Опыт расчета сейсмических откликов от криволинейных геологических границ на основе их явного выделения // Технологии сейсморазведки. -2016. - №4. - С.45 - 51.
Голубев, В.И., Петров, И.Б., Хохлов, Н.И. (2013) Численное моделирование сейсмической активности сеточно-характеристическим методом // Ж. вычисл. матем. и матем. физ. -2013. - Т. 53, № 10. - С. 1709 - 1720.
Голубев, В.И. (2014) Методика отображения и интерпретации результатов полноволновых сейсмических расчётов // Труды МФТИ. - 2014. - Т. 6, № 1 (21). - С. 154 - 161.
Горгуца, Р.Ю., Курило, Е.Ю. (2017) Строительство искусственных ледовых островов в условиях Арктики // Гидротехника. XXI ВЕК. - 2017. - № 4. - С. 56 - 59.
Жуков, А.И. (1960) Применение метода характеристик к численному решению одномерных задач газовой динамики // Тр. МИАН СССР. - 1960. - № 58. - С. 4 - 150.
Заславский, Ю.М., Кержаков, Б.В., Кулинич, В.В. (2008) Вертикальное сейсмическое профилирование на морском шельфе // Акустический журнал. - 2008. - Т.54, №3. -С.483 - 490.
Зволинский Н.В., Шхинек К.Н. (1984) Континуальная модель слоистой среды // Изв. АН СССР. МТТ. - 1984. - № 1. - С.5 - 14.
Зенкевич, О. (1975) Метод конечных элементов в технике. - М.: Мир, 1975. - 540 с.
Зюзина, Н.А., Ковыркина, О.А., Остапенко, В.В. (2018) Монотонная разностная схема, сохраняющая повышенную точность в областях влияния ударных волн // Доклады Академии наук. - 2018. - Т. 482, № 6. - С. 639 - 643.
Иванов, А. М., Хохлов, Н. И. (2018) Параллельная реализация сеточно-характеристического метода в случае явного выделения контактных границ // Компьютерные исследования и моделирование. - 2018. - Т. 10, № 5. - С. 667 - 678.
Иванов, В.Д., Петров, И.Б., Суворова, Ю.В. (1990) Расчет волновых процессов в наследственных вязкоупругих средах // Механ. композитных материалов. - 1990. - № 3. - С. 447 - 450.
Келдыш, М.В. (1942) О методе Галеркина решения краевых задач. - Известия АН СССР, 1942.
Ковыркина, О. А., Остапенко, В. В. (2010) О сходимости разностных схем сквозного счета // Доклады Академии наук. - 2010. - Т. 433, № 5. - С. 599 - 603.
Кондауров, В.И. (2001) Тензорная модель континуального разрушения и длительной прочности упругих тел // Изв. РАН. МТТ. - 2001. - № 5. - С.134 - 151.
Кондауров, В.И., Кукуджанов, В.Н. (1978) Об определяющих уравнениях и численном решении некоторых задач динамики упругопластических сред с конечными деформациями // Сб. по числ. методам в механ. деформируемого твердого тела. - 1978. - С. 84 - 122.
Кондауров, В.И., Никитин, Л.В. (1989) Модель континуального разрушения вязкоупругих тел // Изв. АН СССР. МТТ. - 1989. - № 3. - С. 131 - 139.
Кондауров, В.И., Фортов, В.Е. (2002) Основы термодинамики конденсированной среды. - М.: МФТИ, 2002. - 336 с.
Конторович, А.Э., Нестеров, И.И., Салманов, Ф.К., Сурков, В.С., Трофимук, А.А., Эрвье, Ю.Г. (1975) Геология нефти и газа Западной Сибири. - М.: Недра, 1975. - 680 с.
Кристенсен, Р. (1972) Введение в механику композитов. - М.: Мир, 1972. - 334 с.
Кукуджанов, В.Н. (1999) Микроскопическая модель разрушения неупругого материала и ее применение к исследованию локализации деформаций // Изв. РАН. МТТ. - 1999. - №5. -С. 72 - 87.
Кукуджанов, В.Н. (1985) Численное моделирование динамических процессов деформирования и разрушения упругопластических сред // Успехи механики. - 1985. - Т.8, Вып.4. - С. 21 - 65.
Кукуджанов, В.Н. (1976) Численное решение неодномерных задач распространения волн напряжений в твердых телах // Сообщения по прикл. матем. ВЦ АН СССР. - 1976. -Вып. 6. - С. 11 - 37.
Куликовский, А.Г., Погорелов, Н.В., Семенов, А.Ю. (2001) Математические вопросы численного решения гиперболических систем уравнений. - М.: Наука, 2001. - 608 с.
Куликовский А.Г., Свешникова Е.И. (1998) Нелинейные волны в упругих средах. - М: Моск. лицей, 1998.
Курант Р. (1964) Уравнения с частными производными. - М.: Мир, 1964. - 823 с.
Ландау, Л.Д., Лифшиц, Е.М. (1987) Теоретическая физика, т. 7. - М.: Наука, 1987. - 248 с.
Левянт, В.Б., Петров, И.Б., Челноков, В.Б. (2005) О природе отклика рассеянной сейсмической энергии от зоны диффузной кавернозности и трещиноватости в массивных породах // Геофизика. - 2005. - № 6. - С. 5 - 19.
Магомедов, К.М. (1966) Метод характеристик для численного расчета пространственных течений газа // Ж. вычисл. матем. и матем. физ. - 1966. - Т. 6, № 2. - С. 313 - 325.
Магомедов, К.М., Холодов, А.С. (1988) Сеточно-характеристические численные методы. М.: Наука, 1988. - 287 с.
Магомедов, К.М., Холодов, А.С. (1969) О построении разностных схем для уравнений гиперболического типа на основе характеристических соотношений // Ж. вычисл. матем. и матем. физ. - 1969. - Т. 9, № 2. - С. 383 - 396.
Михайловская, М. Н., Рогов Б. В. (2012) Монотонные компактные схемы бегущего счета для систем уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. - 2012. - Т. 52, № 4. - С. 672 - 695.
Молотков, Л.А., Хило, А.Е. (1986) Эффективные модели слоистых сред с линейными контактами общего вида // Зап. научн. семинаров ЛОМИ. - 1986. - Т. 156. - С. 148 - 157.
Новацкий, В. (1975) Теория упругости. - М.: Изд. "Мир", 1975. - 872 с.
Оден, Дж. Конечные элементы в нелинейной механике сплошных сред. - М.: Мир, 1976. - 464 с.
Олейник, О. А. (1960) Об одном методе решения общей задачи Стефана // Докл. АН СССР. -1960. - Т. 135, № 5. - С. 1054 - 1057.
Остапенко, В. В. (2000) О построении разностных схем повышенной точности для сквозного расчета нестационарных ударных волн // Ж. вычисл. матем. и матем. физ. - 2000. - Т. 40, № 12. - С. 1857 - 1874.
Остапенко, В. В. (1997) О сходимости разностных схем за фронтом нестационарной ударной волны // Ж. вычисл. матем. и матем. физ. - 1997. - Т. 37, № 10. - С. 1201 - 1212.
Петров, И. Б., Фаворская, А. В., Хохлов, Н. И. (2017) Сеточно-характеристический метод на системах вложенных иерархических сеток и его применение для исследования сейсмических волн // Ж. вычисл. матем. и матем. физ. - 2017. - Т. 57, № 11. - С. 1804 -1811.
Петров, И. Б., Фаворская, А. В., Хохлов, Н. И., Миряха, В. А., Санников, А. В., Голубев, В. И. (2014) Мониторинг состояния подвижного состава с помощью высокопроизводительных вычислительных систем и высокоточных вычислительных методов // Матем. моделирование. - 2014. - Т. 26, № 7. - С. 19 - 32.
Петров, И. Б., Стогний, П. В., Хохлов, Н. И. (2021) Математическое моделирование 3D-динамических процессов вблизи трещины с использованием модели Шонберга // Докл. РАН. Матем., информ., проц. упр. - 2021. - Т. 500. - С. 40 - 44.
Петров, И.Б., Тормасов, А.Г. (1990) О численном исследовании трехмерных задач обтекания волнами сжатия препятствия или полости в упругопластическом полупространстве // Докл. АН СССР. - 1990. - Т. 314, № 4. - С. 817 - 820.
Петров, И.Б., Тормасов, А.Г., Холодов, А.С. (1989) О численном изучении нестационарных процессов в деформируемых средах многослойной структуры // Известия АН СССР. Серия МТТ. - 1989. - № 4. - С. 89 - 95.
Петров, И.Б., Тормасов, А.Г., Холодов, А.С. (1990) Об использовании гибридизированных сеточно-характеристических схем для численного решения трехмерных задач динамики деформируемого твердого тела // Ж. вычисл. матем. и матем. физ. - 1990. - Т. 30, № 8. -С.1237 - 1244.
Петров, И.Б., Фаворская, А.В., Санников, А.В., Квасов, И.Е. (2013) Сеточно-характеристический метод с использованием интерполяции высоких порядков на тетраэдральных иерархических сетках с кратным шагом по времени // Матем. моделирование. - 2013. - Т. 25, № 2. - С. 42 - 52.
Петров, И.Б., Холодов, А.С. (1984а) О регуляризации разрывных численных решений уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. - 1984. - Т. 24, № 8. - С. 1172 - 1188.
Петров, И.Б., Холодов, А.С. (1984Ь) Численное исследование некоторых динамических задач механики деформируемого твердого тела сеточно-характеристическим методом // Ж. вычисл. матем. и матем. физ. - 1984. - Т. 24, Вып. 5. - С. 722 - 739.
Петровский И.Г. (1961) Лекции об уравнениях с частными производными. - Физматгиз, 1961. -400 с.
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 1) Программа для расчета сейсмических процессов в гетерогенных геологических средах (SEISMIC_MODELING) /
И. Б. Петров, Н. И. Хохлов, В. И. Голубев; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)». - № 2013619847, 2013617990 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 2) Программа для построения геологических разрезов путём миграции сейсмических данных (SeisMig) / В. И. Голубев, О. Я. Войнов, И. Б. Петров; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)». - № 2015663470, 2015660394 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 3) Программа для трёхмерного морского сейсмического и электроразведочного моделирования и совместной инверсии ^ЕМЛ) / И. Б. Петров, Н. И. Хохлов, Н. Б. Явич, М. С. Маловичко, В. И. Голубев, А. В. Фаворская; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)». - № 2016614469, 2016617024 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 4) Программа для моделирования процесса сейсмической разведки месторождений углеводородов на основе реологии насыщенной пористой среды (STESSHA) / В. И. Голубев; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)». - № 2017617167, 2017614222 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 5) Программа по решению задач сейсмики с учётом динамики разрушенного материала (SeismicDestroy) / В. И. Голубев, Н. И. Хохлов, Д. П. Григорьевых; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)». - № 2018616552, 2018613291 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 6) Программа для гибридных высокопроизводительных вычислительных систем по моделированию процесса сейсмической разведки (SeismicSimulationHybrid) / Н. И. Хохлов, В. И. Голубев, А. М. Иванов, И. Б. Петров; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-
технический институт (государственный университет)». - № 2018617401, 2018614757 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 7) Программа для построения геологических разрезов на основе упругой миграции сейсмических данных (ElasticMig) / В. И. Голубев, А. В. Фаворская; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)». - № 2018616261, 2018613526 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 8) Модуль расчёта волновых полей 'таСШШ" / А. М. Иванов, Н. И. Хохлов, В. И. Голубев, И. Б. Петров, А. В. Екименко; О. с ограниченной ответственностью «Газпромнефть Научно-Технический Центр» (ООО «Газпромнефть НТЦ»). - № 2019664986, 2019663550 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 9) Программа для построения миграционных изображений гетерогенных геологических сред (SeismicElasticMigration) / В. И. Голубев, О. Я. Войнов; В. И. Голубев. - № 2019665498, 2019661830 (Российская Федерация).
Свидетельство о гос. регистрации программы для ЭВМ. (ЭВМ 10) Программа для моделирования процесса сейсмической разведки в горизонтально слоистых анизотропных средах ^ЕКМО_УТГ) / В. И. Голубев; федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт (государственный университет)». - № 2021618231, 2021616901 (Российская Федерация).
Работнов Ю.Н. (1979) Механика деформируемого твердого тела. - М.: Наука, 1979. - 744 с.
Рождественский Б.Л., Яненко Н.Н. (1968) Системы квазилинейных уравнений и их приложения к газовой динамике. - М.: Наука, 1968. - 546 с.
Русанов В.В. (1953) Метод характеристик для пространственных течений газа // Теоретическая гидромеханика. - 1953. № 11, Вып. 3. - С. 71-77.
Салганик, Р.Л. (1987) Приближение сплошной среды для описания деформирования слоистого массива // Изв. АН СССР. МТТ. - 1987. - №3. - С.48 - 56.
Салганик, Р.Л. (2005) Термо- и пороупругое деформирование многослойной структуры, слои которой работают на изгиб (континуальное приближение) // Изв. РАН. МТТ. - 2005. -№3. - С.60 - 65.
Седов Л.И. (1976) Механика сплошной среды. - М.: Наука, 1976. - 576 с.
Стогний, П. В., Петров, Д. И., Хохлов, Н. И., Петров, И. Б. (2018) Численное моделирование сеточно-характеристическим методом влияния ледовых образований на сейсмические отклики // Матем. моделирование. - 2018. - Т. 30, № 8. - С. 107 - 115.
Стогний, П. В., Хохлов, Н. И., Петров, И. Б. (2020) Численное моделирование распространения упругих волн в геологических средах с газовыми полостями с использованием сеточно-характеристического метода // Сиб. журн. вычисл. матем. - 2020. - Т. 23, № 3. - С. 325 -338.
Тихонов А.Н., Самарский А.А. (1966) Уравнения математической физики. - М.: Наука, 1966. -724 с.
Федоренко, Р.П. (1962) Применение разностных схем высокой точности для численного решения гиперболических уравнений // Ж. вычисл. матем. и матем. физ. - 1962. - Т. 2, № 6. - С. 1122 - 1128.
Холодов, А.С. (1980) О построении разностных схем повышенного порядка точности для уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. - 1980. - Т. 20, № 6. - С. 1601 - 1620.
Холодов, А.С. (1978) О построении разностных схем с положительной аппроксимацией для уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. - 1978. - Т. 18, № 6. - С. 1476 - 1492.
Холодов, А.С. (2008) Численные методы решения уравнений и систем гиперболического типа. Энциклопедия низкотемпературной плазмы (серия "В"). Т. VII 1, Ч. 2. Математическое моделирование низкотемпературной плазмы. - М.: ЯНУС К, 2008. - С. 141-174.
Холодов, А.С., Холодов, Я.А. (2006) О критериях монотонности разностных схем для уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. - 2006. - Т. 46, № 9. - С. 1638 - 1667.
Холодов, Я. А., Холодов, А. С., Цыбулин, И. В. (2018) Построение монотонных разностных схем для систем уравнений гиперболического типа //Ж. вычисл. матем. и матем. физ. -
2018. - Т. 58, № 8. - С. 30 - 49.
Хохлов, Н. И., Петров, И. Б. (2019) Применение сеточно-характеристического метода для решения задач распространения динамических волновых возмущений на высокопроизводительных вычислительных системах // Труды ИСП РАН. - 2019. - Т. 31, № 6. - С. 237 - 252.
Шевелев Ю.Д. (1986) Пространственные задачи вычислительной аэрогидродинамики. - М.: Наука, 1986. - 368 с.
Abrate, S. (2016) Damage in laminates from low-velocity impacts // Dynamic Deformation, Damage and Fracture in Composite Materials and Structures. - 2016. - P. 35 - 69.
Abrate, S. (1991) Impact on Laminated Composite Materials // Applied Mechanics Reviews. - 1991. -V. 44, № 4. - P.155 - 190.
Aki, K., Richards, P. G. (1980) Quantitative Seismology: Theory and Methods. — San Francisco: Freeman, 1980.
Alkhimenkov, Y., Khakimova, L., Podladchikov, Y. (2020) Stability of discrete schemes of Biot's poroelastic equations // Geophysical Journal International. - 2020.
Aoki, T. (1991). A universal solver for hyperbolic equations by cubic-polynomial interpolation // Computer Physics Communications. - 1991. - V. 66. - P. 233 - 242.
Arandiga, F., Baeza, A., Donat, R. (2008) Vector cell-average multiresolution based on Hermite interpolation // Adv. Comput. Math. - 2008. - V. 28. - P. 1 - 22.
Argatov, I. (2012) Mathematical modeling of linear viscoelastic impact: Application to drop impact testing of articular cartilage // Tribology International. - 2012. - V. 63. - P. 213 - 225.
Bagaev, R.A., Golubev, V.I., Golubeva, Yu.A. (2019) Full-wave 3D earthquake simulation using the double-couple model and the grid-characteristic method // Computer Research and Modeling. -
2019. - V. 11 (6). - P. 1061 - 1067.
Bakhvalov, N.S., Panasenko, G. (1989) Homogenisation: Averaging Processes in Periodic Media: Mathematical Problems in the Mechanics of Composite Materials. - Springer, 1989. - 402 p.
Bakulin, А. (2003) Intrinsic and layer-induced vertical transverse isotropy // Geophysics. - 2003. - V. 68, № 5. - P. 1708 - 1713.
Batdorf, S. B., Budiansky, B. (1954) Polyaxial stress-strain relations of strain-hardening metal // J. of Applied Mech. - 1954. - V. 21.
Beklemysheva, K.A., Ermakov, A.S., Petrov, I.B., Vasyukov, A.V. (2016) Numerical simulation of the failure of composite materials by using the grid-characteristic method // Mathematical Models and Computer Simulations. - 2016. - V. 8, №5. - P. 557 - 567.
Beklemysheva, K., Golubev, V., Petrov, I., Vasyukov, A. (2021) Determining effects of impact loading on residual strength of fiber-metal laminates with grid-characteristic numerical method // Chinese Journal of Aeronautics. - 2021. - V. 34 (7). - P. 1 - 12.
Beklemysheva, K.A., Golubev, V.I., Vasyukov, A.V., Petrov, I.B. (2019) Numerical Modeling of the Seismic Influence on an Underwater Composite Oil Pipeline // Mathematical Models and Computer Simulations. - 2019. - V. 11 (5). - P. 715 - 721.
Beklemysheva, K.A., Vasyukov, A.V., Golubev, V.I., Zhuravlev, Y.I. (2018) On the Estimation of Seismic Resistance of Modern Composite Oil Pipeline Elements // Doklady Mathematics. -2018. - V. 97 (2). - P. 184 - 187.
Bernardi, C., Maday, Y. (1997) Spectral Methods // in the Handbook of Numerical Analysis V, P.G. Ciarlet and J.-L. Lions Eds., North-Holland, 1997. - P. 209 - 485.
Berryman, J.G. (2008) Exact seismic velocities for VTI and HTI media and extended Thomsen Formulas for stronger anisotropies. - Lawrence Berkeley National Laboratory, 2008. - 73 p.
Bihari, B., Harten, A. (1995) Application of generalized wavelets: an adaptive multiresolution scheme // Journal of Computational and Applied Mathematics. - 1995. - V. 61. - P. 275 - 321.
Biot, M. A. (1956) Theory of propagation of elastic waves in a fluid-saturated porous solid, I: Low-frequency range // The Journal of the Acoustical Society of America. - 1956. - V. 28 (2). - P. 168 - 178.
Blokhin, A. M., Dorovsky, V. N. (1995) Mathematical Modeling in the Theory of Multivelocity Continuum. - New York: Nova Science Publishers Inc, 1995.
Обратите внимание, представленные выше научные тексты размещены для ознакомления и получены посредством распознавания оригинальных текстов диссертаций (OCR). В связи с чем, в них могут содержаться ошибки, связанные с несовершенством алгоритмов распознавания. В PDF файлах диссертаций и авторефератов, которые мы доставляем, подобных ошибок нет.