Аннотация результатов, полученных в 2019 году
Были разработаны компактные продолженные схемы для системы динамических уравнений линейной теории упругости в двумерной постановке. Они основаны на расщеплении исходной системы по пространственным направлениям и дополнении исходных уравнений их дифференциальными следствиями. На примере одномерного нестационарного уравнения исследована сходимость разностной схемы, а также предложен механизм, позволяющий сохранить её порядок аппроксимации для случая разрывных коэффициентов. Были построены математические модели и явно-неявные схемы расчета динамики слоистых и блочных сред при нелинейных контактных условиях с малым параметром вязкости с учётом возможности проскальзывания на межслойных (межблочных) границах. Для перехода к континуальной модели среды введена дополнительная неизвестная функция, имеющая смысл распределенной скорости проскальзывания, и использованы соотношения теории скольжения типа Батдорфа-Будянского, которые применялись многими авторами для построения моделей неупругих сред с непрерывным распределением плоскостей скольжения. Полученные системы уравнений континуальных моделей с дискретными наборами плоскостей скольжения являются полулинейными гиперболическими. Дифференциальная часть этих уравнений совпадает с системой уравнений упругой среды, а нелинейные свободные члены в определяющих уравнениях отражают нелинейные контактные условия на плоскостях скольжения. Для построения численных схем решения предложена явная аппроксимация уравнений движения и неявная аппроксимация определяющих уравнений, содержащих малый параметр вязкости в знаменателе свободных членов. Аналитическое решение алгебраической системы разностных уравнений с использованием разложения по малому параметру естественным образом приводит к разностным формулам, которые представляют собой корректировку упругого шага расчета по времени, учитывающую нелинейное скольжение на межслойных (межблочных) границах. Была проведена последовательная программная реализация разработанных моделей и численных схем. Ввиду того, что реалистичные геологические модели могут содержать десятки и сотни миллионов расчётных узлов, изначально в архитектуру программы была заложена возможность реализации параллельных расчётов с использованием современных технологий (OpenMP, MPI, CUDA, OpenCL), что планируется осуществить на следующих этапах проекта. В качестве языка разработки был выбран C++, позволяющий существенно оптимизировать скорость расчётов. Был разработан модуль, осуществляющий корректировку после упругого шага по времени. В расчётных формулах участвуют также и значения тензора напряжений с предыдущего временного слоя, поэтому в корректоре создаётся копия расчётной сетки, хранящая их. В каждом узле сетки дополнительно к вектору скорости и компонентам тензора напряжений хранится и величина скорости проскальзывания границ слоев. Разработанная авторами модель и явно-неявные схемы для описания поведения слоистых сред для нелинейных контактных условий была применена к геолого-геофизической модели месторождения Западной Сибири. В этом нефтедобывающем регионе наибольший интерес представляет изучение Баженовской свиты и её аналогов. По некоторым оценкам, содержащиеся в ней ресурсы нефти могут составлять до 500 млрд тонн. Коллекторы Баженовской свиты могут быть приурочены к микрослоистым, листоватым аргиллитам с межслоевой пустотностью с высоким содержанием органического вещества. Пустотное пространство коллекторов сформировано в основном в результате автофлюидоразрыва пород по слоистости, возникающего под действием процессов преобразования органического вещества из твердого состояния в жидкое. Для описания данного пласта было исследовано применение стандартной упругой и новой континуальной модели. Для точечного источника получены типичные сейсмограммы общего пункта возбуждения (ОПВ). На сейсмограммах выделяются основные типы волн: прямая продольная волна, прямая поперечная волна от источника, последовательность отражённых волн от геологических границ. Интерес представляет отрицательное отражение — это отражение от кровли Баженовской свиты. Был выполнен анализ волнового поля, заключавшийся в построении графиков изменения амплитуд отраженной волны для каждой трассы. Такая методика анализа широко распространена при интерпретации сейсмических данных и позволяет судить о том, какая часть энергии продольной волны претерпевает обмен на границе раздела геологических слоёв. В свою очередь коэффициент обмена волн зависит от соотношения скоростей продольных и поперечных волн, которое является важным индикатором литологии и фильтрационно емкостных свойств. Было выявлено, что энергия отражения на малых удалениях больше, а её уменьшение происходит существенно быстрее, при использовании континуальной модели, чем для модели упругой среды.

Аннотация результатов, полученных в 2020 году
Были проведены эмпирические оценки порядка сходимости построенных компактных продолженных схем для двумерных динамических задач. Подтверждено достижение третьего порядка аппроксимации при распространении возмущения вдоль каждой из координатных осей. По норме L1 и Linf достигнут второй порядок аппроксимации для произвольного направления распространения волны. Таким образом, на компактном шаблоне был достигнут повышенный порядок сходимости. Были построены неявные схемы расчета определяющих уравнений слоистых и блочных сред для нелинейных контактных условий с малым параметром вязкости с учетом возможного отслоения. Условия контактного взаимодействия, зависящие от физико-механических свойств прослоек, были приняты в виде условий скольжения трения Кулона с малой вязкой добавкой при сжатии и слабого отслоения с малыми коэффициентами упругой связи при растяжении. Были введены дополнительные непрерывные функции координат и времени, имеющие смысл распределенных скоростей скольжений и отслоений. С использованием соотношений теории скольжения был учтен вклад скоростей скольжений и отслоений в тензоры скоростей неупругой деформации соответственно. Система замыкалась уравнениями движения. Для её решения использовался явно-неявный метод. Неявная аппроксимация применялась для определяющих уравнений, которые содержат малый параметр в знаменателе свободного члена, уравнения движения аппроксимировались явно. В результате аналитического решения неявной аппроксимации уравнений для касательных напряжений и аналогичной аппроксимации уравнений для нормальных напряжений были получены корректировочные формулы после упругого шага с учётом слабого отслоения. Была проведена верификация расчётного алгоритма на известных аналитических и полуаналитических решениях: сейсмический отклик от свободной границы в упругом полупространстве, вертикально направленная внешняя сила и два контактирующих полупространства с различными упругими свойствами. При измельчении расчётной сетки подтверждено совпадение результатов. Было проведено сопоставление результатов, получаемых при моделировании трещиноватой геологической среды с явным выделением границ неоднородностей и в рамках континуального подхода. Для случая вертикально ориентированных трещин с увеличением количества трещин (пространственной плотности их расположения) сейсмические отклики становятся всё более схожими. Однако в случае наклонных трещин наблюдаются отличия, которые объясняются тем, что в данной постановке реализуется режим растяжения в континуальной модели, который не рассматривается при явном выделении контактных границ. Это нелинейное поведение материала существенно влияет на отраженные и прошедшие волны. Был проведён расчёт сейсмического поля в модели Marmousi2 с добавленным трещиноватым объектов (в рамках континуальной модели слоистой среды). Выявлено, что для случая одного точечного источника на сейсмограммах гиперболические оси отраженных волн практически не видны; также регистрируется большое количество дифрагированных волн. Сейсмограммы для моделей с трещинами и без трещин очень похожи. На разностной сейсмограмме хорошо виден конус рассеянных волн. Его амплитуда на порядок меньше амплитуды волн, отраженных от геологических границ. Аналогичная картина наблюдается при анализе данных с источником в виде продольной волны. На внешних частях сейсмического разреза наблюдаются вытянутые синфазные оси отраженных волн, отраженная волна от морского дна и кратные волны хорошо видны. Центральная часть практически лишена регулярных отражений. Как и в случае с точечным источником, рассеянные волны регистрируются на разностной сейсмограмме. Была построена геологическая модель, основанная на свойствах и геометрии слоёв, характерных для Западной Сибири. Скорость распространения продольных волн определялась по результатам акустического каротажа, а скорость S-волн и плотность рассчитывались по эмпирическим соотношениям. Область интереса - низкоскоростной пласт, соответствующий баженовской свите. Он описывался с использованием нелинейной континуальной модели ввиду мелкослоистости. Результаты расчётов показали, что полученное волновое поле не зашумленное, паразитные отражения от краев расчетной области имеют низкую интенсивность и не усложняют интерпретацию сигнала. Целевое отражение легко распознается как отрицательная фаза.

Аннотация результатов, полученных в 2021 году
В ходе реализации проекта были построены компактные сеточно-характеристические схемы для численного решения линейных гиперболических систем уравнений. При построении интерполяционного полинома третьей степени для одномерного линейного уравнения переноса используются хранимые в расчётных узлах значения искомых функций и их пространственных производных. Кроме значительного возрастания времени вычислений показано, что возрастает и объём потребляемой оперативной памяти. Всё это позволяет сделать вывод о целесообразности применения схем на расширенных шаблонах совместно с согласованной схемой пространственного расщепления для численного решения реальных практических задач сейсморазведки. Построены явно-неявные схемы, объединяющие схемы для явного расчета упругого шага с неявным шагом аппроксимации определяющих уравнений нелинейной системы. Неявная аппроксимация применяется только для тех уравнений, которые содержат малый параметр в знаменателе свободного члена, остальные уравнения аппроксимируются явно. Отметим, что построенная процедура корректировки упругого решения применима и для тех случаев, когда параметр вязкости у прослоек, обеспечивающих режим скольжения контактных границ, не мал. Разработан вычислительный алгоритм, позволяющий проводить расчёт сейсмических процессов в VTI средах, содержащих трещиноватые включения. При этом может использоваться как однобереговая модель трещины, так и двухбереговая. На возникающих внутренних границах (два борта трещины) применяются физически корректные условия скольжения. Разработаны параллельные алгоритмы для систем с общей и распределённой памятью, выполнена их программная реализация на языке C++. Отметим, что использование структурных сеток позволило фактически снять ограничение на количество узлов в общей расчётной сетке задачи. Для этого декомпозиция производится в пространстве индексов. Фактическая же генерация сеток производится раздельно на каждом вычислительном узле. Проведена серия верификационных расчётов на простейших трёхмерных моделях и на моделях, построенных с привязкой к реальным геологическим объектам, в том числе при наличии трещиноватых включений. Результаты проведённых компьютерных расчётов трещиноватых сред позволяют сделать вывод, что эффект межтрещинного взаимодействия, видимый при явном выделении отдельных трещин, не воспроизводится при использовании континуальной модели, а амплитуда сигнала-отклика увеличивается с числом моделируемых неоднородностей. Построены многостадийные алгоритмы расщепления по пространственным направлениям, сохраняющие повышенный порядок сходимости многомерной схемы. Они обеспечивают достижение третьего порядка сходимости итоговой схемы. Отметим, что сохраняется возможность использования пространственной локальности задачи, то есть в программе хранятся лишь значения на одном временном слое и не увеличивается потребление оперативной памяти.