Geum.ru
   Авторефераты по всем темам  >>  Авторефераты по техническим специальностям

На правах рукописи

РЕШЕТОВА ГАЛИНА ВИТАЛЬЕВНА

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СЕЙСМИЧЕСКИХ И СЕЙСМОАКУСТИЧЕСКИХ ВОЛНОВЫХ ПОЛЕЙ В РАЗНОМАСШТАБНЫХ И РЕЗКОКОНТРАСТНЫХ СРЕДАХ

05.13.18 - математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени доктора физико-математических наук

НОВОСИБИРСК - 2010

Работа выполнена в Учреждении Российской академии наук Институте вычислительной математики и математической геофизики Сибирского отделения РАН (ИВМиМГ СО РАН)

Научный консультант: доктор физико-математических наук, академик РАН Михайленко Борис Григорьевич

Официальные оппоненты:

доктор физико-математических наук, профессор Лаевский Юрий Миронович доктор физико-математических наук, профессор Роменский Евгений Игоревич доктор физико-математических наук, профессор Крукиер Лев Абрамович

Ведущая организация: Учреждение Российской академии наук Институт вычислительного моделирования Сибирского отделения РАН (ИВМ СО РАН, г. Красноярск)

Защита состоится 15 февраля 2011 г. в 15 час. на заседании диссертационного совета Д 003.061.02 при Учреждении Российской академии наук Институте вычислительной математики и математической геофизики СО РАН.

Адрес: пр-т Ак. Лаврентьева, 6, Новосибирск, 630090.

Факс: (383) 3308783, e-mail: kosova@rav.sscc.ru

С диссертацией можно ознакомиться в библиотеке ИВМиМГ СО РАН.

Автореферат разослан ____ ноября 2010 г.

Ученый секретарь диссертационного совета д.ф.-м.н. С.Б. Сорокин

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

К ним относится, например, изучение сейсмоакустических волновых полей при акустическом каротаже скважин. Действительно, даже в простейшем случае открытой (необсаженной) скважины приходится иметь дело с двумя характерными размерами - диаметром скважины 0.1-0.2 м и максимальным расстоянием источник - приемник 10-12 м. При наличии обсадки приходится также учитывать влияние стальной трубы, отделяющей заполненную жидкостью скважину от вмещающей среды. Толщина стенок трубы составляет первые сантиметры и является еще одним масштабом задачи.

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

За последние десятилетия накоплен огромный экспериментальный материал об ионосферных возмущениях, вызванных сильными землетрясениями и техногенными воздействиями, такими как промышленные взрывы, подземные ядерные испытания. Возбужденные в литосфере сейсмические колебания генерируют в атмосфере акустогравитационные волны, которые регистрируются дистанционным радиозондированием. Этот эффект может быть положен в основу методов локации землетрясений и подземных ядерных взрывов. Однако лишь в небольшом количестве работ сделаны попытки с теоретических позиций объяснить механизмы этих связей в системе литосфераатмосфера-ионосфера. Численное моделирование способствует выявлению особенностей формирования волнового процесса и созданию на этой основе надежных систем контроля. Но его выполнение сопряжено с существенными вычислительными затратами из-за больших перепадов скоростей сейсмических волн в литосфере и акустогравитационных волн в атмосфере (в 15-20 раз), а также огромных перепадов плотностей в упругом полупространстве и в верхних слоях атмосферы (значение плотности экспоненциально падает с высотой).

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

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

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

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

3. На основе интегрального преобразования Лагерра разработать спектрально-разностный метод численного моделирования сейсмических волновых полей для трхмерных неоднородных сред с резкоконтрастными границами раздела. Создать на этой основе научно-исследовательский вариант параллельного программного обеспечения.

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

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

- теория конечно-разностных схем для аппроксимации начальнокраевых задач теорий упругости и вязкоупругости;

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

- методы функционального анализа для построения алгоритмов интерполяции в задаче пространственного измельчения шагов конечноразностной сетки с сохранением порядка аппроксимации конечноразностной схемы;

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

- теория параллельных вычислений, в том числе методы пространственной декомпозиции областей для разработки параллельных версий алгоритмов, а также использование последних эффективных разработок программного интерфейса обмена сообщениями в стандарте MPI2, в частности, новые возможности параллельных коллективных обменов и ввода/вывода в среде MPI-2 I/O.

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

Окончательная верификация программного обеспечения проводилась в нефтяных и сервисных компаниях (Schlumberger, Total, ООО РН-КрасНИПИНефть) и Югорском НИИ Информационных технологий путм сопоставления результатов численного моделирования и лабораторных и натурных наблюдений.

Защищаемые научные результаты 1. Теоретически и экспериментально обоснованный конечноразностный метод численного моделирования сейсмоакустических волновых полей при акустическом каротаже в трхмерных неоднородных средах с поглощением, разработанный с использованием адаптивных пространственных сеток в цилиндрической системе координат. Научно-исследовательский вариант параллельного программного обеспечения (с графическим интерфейсом пользователя).

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

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

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

1.1. Интегро-интерполяционным методом (методом баланса) построена консервативная конечно-разностная схема для численного моделирования обобщенной модели стандартного линейного тела GSLS для нескольких релаксационных механизмов, проведен анализ устойчивости схемы и исследованы ее дисперсионные свойства для минимизации численной дисперсии скорости.

1.2. В ходе реализации метода разработаны оригинальные подходы для эффективного решения ряда ключевых задач:

- использование экономичного -метода, обеспечивающего заданное поведение добротности в модели вязкоупругой среды;

- ограничение расчетной области слабо отражающими слоями в цилиндрической системе координат с помощью модифицированного метода PML (аббревиатура английского Perfectly Matched Layer) без расщепления переменных по азимуту;

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

- радиальное и азимутальное измельчение сеток, адаптирующихся к неоднородностям численной модели среды;

- организация параллельных вычислений методом декомпозиции расчетной области с использованием библиотеки MPI-2.

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

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

Создан научно-исследовательский вариант параллельного программного обеспечения для численного расчта сейсмических волновых полей в трхмерных неоднородных разномасштабных упругих средах, основывающийся на новых возможностях параллельных коллективных обменов и параллельного ввода/вывода больших массивов данных с использованием библиотек MPI-2 I/O (Input/Output):

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

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

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

3. Разработан, теоретически и экспериментально обоснован и протестирован новый спектрально-разностный метод численного моделирования сейсмических волновых полей для трхмерных неоднородных сред с резкоконтрастными границами раздела, основывающийся на отделении времени с использованием интегрального преобразования Лагерра (совместно с Б.Г.Михайленко).

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

Создан научно-исследовательский вариант параллельного программного обеспечения для численного моделирования волновых полей в трхмерных неоднородных средах с резкоконтрастными границами раздела. Численными экспериментами установлено существование предсказанной ранее теоретически поверхностной волны Стоунли-Шолтэ и модифицированной волны Лэмба.

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

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

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

- для исследования проявлений неоднородности зоны проникновения в сейсмоакустических волновых полях и для оценки на этой основе фильтрационных свойств коллектора;

- для оценки качества выполненного гидроразрыва пласта;

- для исследования скользящих волн (creeping waves) на стенке скважины для определения анизотропии прискважинной зоны.

2. Разработанный автором конечно-разностный метод с локальным пространственно-временным измельчением сеток позволяет устойчиво выделять и корректно описывать тонкие эффекты взаимодействия сейсмических волн с кавернозно-трещиноватыми коллекторами.

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

- для определения фильтрационно-емкостных характеристик кавернозно-трещиноватых коллекторов по полю рассеянных волн;

- при моделировании волновых полей в средах с детальным описанием неоднородных включений (конструкций).

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

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

- раннее обнаружение цунамигенного землетрясения;

- контроль подземных ядерных взрывов.

Реализация результатов Результаты диссертации составили основу проекта Численное моделирование взаимодействия сейсмических волн с кавернознотрещиноватыми коллекторами в трхмерно-неоднородных средах реалистичного строения, занявшего первое место в Конкурсе проектов в сфере высокопроизводительных вычислений под девизом Невозможное стало возможным, организованного корпорацией Intel совместно с Российской корпорацией нанотехнологий (РОСНАНО) в 2010 году.

Представленный в диссертации спектрально-разностный метод на основе интегрального преобразования Лагерра получил развитие в серии работ зарубежных и отечественных исследователей (Lk, Chung, Sarkar, Jung, А.А. Михайлов, А.Ф. Мастрюков, В.Н.Мартынов).

Эффективность этого метода для численного моделирования электромагнитных волновых полей продемонстрирована в работе Chung et al.

(2003), в которой убедительно показано, что применение преобразования Лагерра дает ускорение времени вычислений в 100 раз по сравнению с применением обычных конечно-разностных схем при численном моделировании двумерной системы уравнений Максвелла.

Результаты, представленные в диссертации, были включены в Основные результаты исследований Сибирского отделения РАН по приоритетным направлениям развития науки и техники и в Список достижений СО РАН в 2003, 2005, 2006 годах.

Основным направлением научно-исследовательских работ, выполняемых в ИВМиМГ СО РАН с участием диссертанта, являлось развитие данного подхода к численному моделированию подземных и наземных ядерных взрывов в рамках проекта по программе СО РАН 1.4.1.1. Математическое моделирование природных и техногенных геофизических полей в средах сложной геометрии и реологии.

Все исследования, проводимые по теме диссертации, являются составной частью планов НИР Института, а их выполнение постоянно поддерживалось Российским фондом фундаментальных исследований в рамках проектов № 00-05-65323, № 04-05-64177, № 06-0564149, № 07-05-00538 № 10-05-00233.

Внедрение научных результатов Созданный диссертантом программный продукт по моделированию сейсмических волн в двумерно-неоднородных сложнопостроенных упругих средах был внедрен в мобильный программноаппаратный комплекс, созданный в рамках комплексного проекта Разработка комплексной технологии поиска и разведки углеводородов в сложнопостроенных, глубокозалегающих месторождениях по Госзаказу 2005-РИ-00.0/009/202, шифр - RU.IBMMG.00103-01 34 (Акт о внедрении от 25.11.2006).

Разработанное автором программное обеспечение для моделирования трхмерных сейсмоакустических волновых полей в средах с затуханием (акустический каротаж) в настоящее время используется инженерными технологическими центрами компании Schlumberger (в рамках контракта с ИНГГ СО РАН, Акт о внедрении от 01.08.2009).

Разработанное диссертантом программное обеспечение для моделирования сейсмических волновых полей в трехмерных неоднородных средах с субсейсмическими неоднородностями (кавернознотрещиноватые коллекторы) используется ООО Геола с целью изучения проявлений ориентации коридоров трещиноватости по полю рассеянных волн (Акт о внедрении от 14.09.2010).

Апробация, публикации, объем и структура диссертации Результаты диссертационной работы известны научной общественности. Всего по теме диссертации автором лично и в соавторстве опубликовано более 100 работ, в том числе 38 статей, из которых 15 - в ведущих рецензируемых научных журналах из перечня ВАК. Результаты диссертационной работы докладывались и обсуждались на Международных научных конференциях в России и за рубежом.

Основные:

- Генеральные Ассамблеи Европейского геофизического общества (EGU) - 2000, 2002, 2004;

- Коференции Европейской ассоциации геофизиков (EAGE) - 2000, 2003, 2007, 2008, 2010;

- Конференции сообщества геофизиков-исследователей (SEG) - 2006, 2008, 2010;

- Международные конференции "Математические и численные аспекты теории распространения волн" - 2007, 2009;

- 7-я Международная европейская конференция по вычислительной математике, Австрия, Грац, 2007;

- 8-я Международная конференция "Математические и численные аспекты теории распространения волн", Великобритания, Рединг, 2007.

Диссертация состоит из введения, трех глав и заключения. Содержит 298 страниц, в том числе 98 рисунков. Библиография содержит 171 наименование.

Благодарности Успешному проведению исследования способствовала поддержка академиков РАН А.Н. Коновалова и Б.Г. Михайленко, оказавших большое влияние на формирование научных взглядов соискателя.

Автор глубоко благодарна своему коллеге и соавтору д.ф.-м.н.

В.А. Чеверде за содержательные и плодотворные обсуждения, моральную помощь при выполнении работы.

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

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

СОДЕРЖАНИЕ РАБОТЫ

Глава 1. Численное моделирование сейсмоакустических волновых полей при акустическом каротаже в трхмерных неоднородных средах с поглощением Первый раздел посвящн истории вопроса и описанию современных методов численного моделирования сейсмоакустических волновых полей при акустическом каротаже. Сформулированы особенности постановки, связанные с наличием в упругом пространстве скважины, заполненной жидкостью (Biot, 1952; Крауклис П.В. и Крауклис Л.А., 1976). Введено понятие поглощения сейсмической энергии и прокомментированы математические модели, позволяющие описывать волновые процессы с затуханием (Chistensen, 1982; Robertsson et al., 1994).

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

Во втором разделе формулируется математическая постановка задачи и вводятся основные модели среды с поглощением: модель Максвелла, модель Фойгта, стандартная линейная модель тврдого тела (Curtin and Sternberg, 1962; Christensen, 1982). Наиболее подробно исследуется последняя модель, для которой получены частотные зависимости фазовой скорости и добротности. Далее вводится модель обобщнного стандартного линейного тврдого тела и доказывается, что для разумных значений добротности (Q>10) для е описания достаточно двух релаксационных механизмов.

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

Уравнения движения:

ur 1 r rr rr rz ;

t r r z r z zz rz uz rz ; (1) t r r z r u z 2 r r ;

t r r z r Закон Гука:

L rr P S S r M (1 L ) 2(1 L ) 21 L ur rrrl ;

t l L zz P S S z M (1 L ) 2(1 L ) 21 L uz rzzl ;

t l ur L P S S M (1 L ) 2(1 L ) 21 L 1 u r ;

l t r r l (2) ur L rz S 1 L uz rrzl ;

t r z l z u L S 1 L 1 uz rzl ;

t r z l 1 ur u L r S 1 L u rrl ;

t r r r l r r0, 0, z z0 где введены обозначения div u, F(t).

2 r Уравнения на переменные памяти:

rrrl 1 ur P S S (M 2 ) 2 rrrl ;

t l r rzzl 1 uz P S S (M 2 ) 2 rzzl ;

t l z rl 1 P S u 1 ur S (M 2 ) 2 rl ;

t l r r (3) rrzl 1 ur uz S rrzl ;

t l z r rzl 1 S 1 ur u rzl ;

t l r z rrl 1 S u 1 ur u rrl ;

t l r r r Система уравнений (1)-(3) решается с нулевыми начальными данными ur t0 u t0 uz t0 0; t0 0 (4) rr t0 zz t0 r t0 rz t0 z tи условиями контакта на границе скважины r R0, заключающимися в непрерывности нормальных компонент вектора скорости смещений и тензора напряжений:

[ur ] 0;[ ] [ ] ;[ ] ;[ ] 0. (5) rR0 rr rR0 rr rR0 rz rR0 r rRДля ограничения расчтной области автором реализована оригинальная версия PML, разработанная при е непосредственном участии (Костин, Решетова, Чеверда, 2006).

Конечно-разностная аппроксимация задачи (1)-(5) строится на сдвинутых сетках (Virieux, 1986) методом баланса (Самарский, 1977).

Полученная при этом явная конечно-разностная схема является условно устойчивой (доказан спектральный признак устойчивости) и имеет второй порядок аппроксимации всюду вне границ разрыва упругих параметров среды. На границах разрыва она аппроксимирует условия согласования типа (5) с первым порядком аппроксимации, однако е решение сходится к решению исходной задачи со вторым порядком (Самарский, 1977; Lisitsa, Podgornova and Tcheverda, 2010).

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

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

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

Особенностью конструкции скважины является наличие весьма сложной структуры в окрестности е границы - обсадной колонны, цементного кольца, а также трхмерной неоднородной зоны проникновения. Для детального описания таких особенностей при расчетах вводятся специальные транзитные зоны, в которых происходит уменьшение/увеличение шагов по радиусу по экспоненциальному закону: hi q1hi1 (рис.2). Конечно, при этом аппроксимация становится уже первого порядка. Однако за счт выбора оптимального для заданного частотного диапазона коэффициента изменения шагов удатся добиться приемлемого уровня погрешности.

б) а) Рис.1. Согласование сеток при азимутальном измельчении. Красным отмечены узлы, в которых должны быть вычислены дополнительные значения переменных на грубой и мелкой сетках, по известным узлам (зеленым).

Рис.2. Сетка с периодическим измельчением по азимуту и квазиравномерным изменение шагов по радиусу в окрестности границы скважины.

r Устранение математической особенности при. Используемая сетка начинается с окружности r 0.5hr и, следовательно, не содержит ось r 0. Однако для вычисления некоторых компонент вектора скорости смещений и тензора напряжений необходимо знание компонент ur, , при r 0. Так как скважина заполнена жидкоr rz стью, внедиагональные компоненты тензора напряжений внутри не равны нулю, поэтому в нуле надо определить только радиальную компоненту вектора скорости смещений. Для этого используется интерполяция полиномами Лагранжа по шести соседним узлам.

Организация параллельных вычислений производится декомпозицией расчтной области на непересекающиеся дискообразные элементы, каждый из которых приписывается своему процессорному элементу. На каждом временном шаге соседние процессорные элементы обмениваются текущими значениями векторов скорости смещений и тензора напряжений, рассчитанными на соприкасающихся гранях элементарных дисков (рис.3). Высокая эффективность и масштабируемость параллельного программного обеспечения была достигнута благодаря организации асинхронных вычислений. Для этого использовалось сочетание неблокирующих MPI-процедур межпроцессорного обмена Isend/Irecv и специальной организации вычислений внутри элементарной области - из е самой внутренней точки по направлению к соседним процессорным элементам. Оценка точности вычислений как для идеально-упругой среды, так и среды с поглощением проводилась путм сравнения численного решения с точным, а также выделением таких характеристик волнового поля, как времена вступления отдельных волн (продольных и поперечных) и уменьшение их амплитуды.

В заключение представлены результаты расчтов для нескольких реалистичных сейсмоакустических моделей, одна из которых показана на рис. 4. Расчтная область - цилиндр радиусом 1 м и высотой 6 м, радиус скважины (1) - 0.1 м, толщина охватывающей е стальной обсадной колонны (2) - 0.01 м, цементного кольца (3) - 0.04 м, толщина вертикальной трещины (7) - 0.02 м. Параметры модели:

1. буровой раствор в скважине: 1000кг / м3,Vp 1500 м / с,Q 65;

2. обсадная колонна: 7830 кг / м3,Vp 5600 м / с, Vs 3270 м / с,Q 100;

3. цементное кольцо: 2400 кг / м3,Vp 4200 м / с, Vs 2425 м / с,Q 80;

4. упругая среда №4: 2400 кг / м3,Vp 4989 м / с, Vs 2605 м / с,Q 100;

5. упругая среда №5: 2400 кг / м3,Vp 3208 м / с, Vs 1604 м / с,Q 60;

6. упругая среда №6: 2400 кг / м3,Vp 2650 м / с, Vs 1219 м / с,Q 15;

7. вертикальная трещина заполнена тем же самым буровым раствором, что и скважина.

Рис.3. Декомпозиция расчётной Рис.4. Модель строения среды и констобласти. рукции скважины для моделирования проявления вертикальной трещины в сейсмоакустических полях: а) вертикальное сечение; б) горизонтальное сечение; в) увеличение в окрестности скважины.

Рис.5. Серия моментальных снимков разности компонент тензора наrr пряжения, рассчитанных для модели с вертикальной трещиной.

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

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

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

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

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

Согласование сеток. Интересующие нас рассеянные волны, возникающие на микроструктуре пласта-коллектора, имеют амплитуду порядка 1% от интенсивности падающей волны. Поэтому разрабатываемый конечно-разностный метод должен обладать искусственными отражениями от границы смены шагов на уровне 0.0001 0.001% от интенсивности падающей волны. Для обеспечения этого уровня выполнено измельчение по времени и по пространству на разных поверхностях, охватывающих целевую область (рис.6а).

б) в) а) Рис.6. а) Взаимное расположение поверхности с измельчением шага по времени (красный прямоугольник) и шага по пространству. б) Согласование шагов по времени. в) Согласование шагов по пространству.

Измельчение шага по времени. Для одномерного случая измельчение шага по времени при фиксированном шаге по пространству представлено на рис. 6б и достаточно очевидно. Его модификации в двумерной и трхмерной постановках отличаются только использованием соответствующего пространственного шаблона, оставляя неизменной структуру вычислений (Reshetova et al., 2009).

Измельчение шага по пространству. Для согласования сеток с различными шагами по пространству используется интерполяция на основе БПФ. Чтобы упростить объяснение, рассмотрим двумерную реализацию этого подхода. На рис. 6в изображено взаимное расположение крупной и мелкой сеток. Треугольниками обозначены узлы, в которых задаются компоненты векторов скорости смещения (ориентация треугольника совпадает с направлением компоненты вектора), а прямоугольниками и косыми крестиками - узлы, в которых задаются соответственно диагональные ( и ) и внедиагональные ( ) xx zz xz компоненты тензора напряжений. Рассмотрим шаг, на котором по заданным компонентам вектора скоростей смещения вычисляются компоненты тензора напряжений. Как видно, для того чтобы вычислить компоненту на границах мелкой сетки, необходимо знать векторы xz скоростей смещения в узлах, отмеченных красным треугольником. Но эти узлы не принадлежат крупной сетке, поэтому искомые значения в них должны быть получены интерполированием, для чего используется преобразование Фурье. Определяющими факторами этого выбора явились экспоненциальная точность такой интерполяции и высокое быстродействие БПФ. Единственное отличие в трхмерном случае - необходимость выполнения двумерной интерполяции.

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

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

б) а) Рис.7. а) Сверху вниз: одномерная, двумерная и трёхмерная декомпозиции области. б) Ускорение при использовании одномерной (чёрная), двумерной (синяя) и трёхмерной (красная) декомпозиции расчётной области.

Взаимодействие процессорных элементов в группе и между группами. Для организации параллельных вычислений используется пространственная декомпозиция области на соприкасающиеся параллелепипеды, каждый из которых приписывается своему процессорному элементу (рис.7а). Используемый пространственный шаблон требует для расчта скоростей и напряжений в точке M в последующий момент времени их знания в некоторой окрестности этой точки в предыдущий момент. Итак, возникает необходимость организации пошаговых обменов между соприкасающимися процессорами. Время, затрачиваемое на такие обмены, относится к непроизводительным потерям, и чем больше объм пересылаемых данных, тем значительнее эти потери и, следовательно, ниже эффективность разработанного программного обеспечения. В свою очередь, объм пересылаемых данных, а значит, и эффективность программного обеспечения зависят от геометрии декомпозиции области. Как видно из рис.7б, наилучшей эффективностью обладает трхмерная декомпозиция, что объясняется тем фактом, что среди всех параллелепипедов фиксированного объма куб имеет наименьшую площадь поверхности и, следовательно, при этом пересылается наименьший объм данных.

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

Взаимодействие между группами процессорных элементов.

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

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

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

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

Максимально возможная степень асинхронности вычислений при организации межгрупповых обменов и обменов внутри групп обеспечивается за счт использования неблокирующих процедур MPI Isend/Irecv и специальной организации вычислений - на каждом шаге по времени они начинаются из самой внутренней точки области, приписанной конкретному процессорному элементу.

Рис.8. а) Две группы процессоров. б) Взаимодействие между ними Пример моделирования рассеянных волн на коридорах трещиноватости. В заключение приводятся результаты численного моделирования для реалистичной модели карбонатного резервуара с коридорами трещиноватости. Такого типа коллекторы характерны для месторождений нефти и газа в Восточной Сибири, в частности в Юрубчено-Тохомской зоне. В основу построенной численной модели резервуара положена схематическая модель, представляющая погруженный в однородную среду однородный слой, содержащий в себе два слоя с коридорами трещиноватости. Параметры вмещающей среды выбраны соответствующими типичной модели карбонатных пород: Vp 4500 м / c, Vs 2500 м / с, 2500 кг / м3.

Резервуар содержит область трещиноватости коридорного типа с изменчивостью концентрации трещин от 0 в цельных фациях до максимального значения 0.3 (рис.9). Для дискретизации резервуара использовался пространственный шаг 0.5 м. В итоге получается массив из 4000х4000х400 точек. Декомпозиция области была выполнена таким образом, чтобы на каждый из процессорных элементов приходился куб 200х200х200 точек (примерно 1Gb RAM), что потребовало 8процессорных элементов. Вмещающая среда заполняет куб 2000х2000х2000 м и дискретизуется с шагом 2.5 м, что соответствует массиву 800х800х800 точек. Выбирая такую же загруженность процессорных элементов, получаем необходимость использования процессорных элементов. Таким образом, для моделирования взаимодействия сейсмических волн с тонкой структурой данного резервуара необходимо 864 процессорных элемента.

а) б) Рис.9. Распределение коридоров трещиноватости в резервуаре. а) вид сбоку, б) вид сверху.

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

а) б) Рис.10. Трёхкомпонентные сейсмограммы, вычисленные на профилях, параллельных (а) и ортогональных (б) коридорам трещиноватости.

Глава 3. Спектрально-разностный метод численного моделирования сейсмических волновых полей в резкоконтрастных средах на основе преобразования Лагерра по времени Описанные в предыдущих главах численные методы моделирования сейсмических волновых процессов основаны на явных конечноразностных схемах. Эти схемы являются условно устойчивыми, то есть величина шага по времени связана с шагом по пространству так называемым условием Куранта Ch /Vmax, где через Vmax обозначена максимальная скорость в среде. В то же время, все явные схемы обладают численной дисперсией, что ведт к необходимости выбора шага по пространству, подчиняющегося неравенству h Vmin /(k ), где есть доминирующая частота сигнала, а k - некоторый множитель, зависящий от выбранной схемы (так, для схемы на сдвинутых сетках второго порядка k=20). Из приведнных соотношений видно, что в средах с высокими контрастами скоростей применение явных конечно-разностных схем, особенно при необходимости проведения моделирования волновых полей на большие расстояния и, следовательно, для больших времн, связано со значительными затратами вычислительных ресурсов.

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

Указанного недостатка лишн подход с применением преобразования Лагерра по времени. В этом случае получается строго отрицательно определнный эллиптический оператор. Более того, этот оператор не зависит от параметра разделения, что существенно уменьшает затраты на вычисление и хранение соответствующей матрицы. Этот подход предложен и обоснован при непосредственном участии автора (Конюх (Решетова), Михайленко, Михайлов; 2001; Mikhailenko, Mikhailov, Reshetova, 2003). Одно из приложений подхода, опирающегося на преобразование Лагерра, подробно рассматриваемое в диссертации, связано с моделированием взаимодействия сейсмических и акустогравитационных волновых полей для модели неоднородной среды литосфера - атмосфера.

История вопроса. При математическом моделировании распространения сейсмических волн при землетрясениях и взрывах обычно исходят из того, что поверхность Земли граничит с вакуумом, что позволяет рассматривать е как свободную от напряжений. При этом не учитываются эффекты взаимодействия волн на дневной поверхности, приводящие к возбуждению акусто-гравитационных волн в неоднородной атмосфере. Однако появились теоретические и экспериментальные работы, в которых приводятся экспериментальные данные о существовании взаимозависимости волновых процессов в литосфере и атмосфере. В частности, А.С. Алексеевым с соавторами описан эффект акусто-сейсмической индукции (Alekseev et al., 1996), заключающийся в возбуждении акустической волной от мощного вибратора интенсивных сейсмических поверхностных волн, распространяющихся на расстояние в десятки километров. Теоретическим исследованиям волновых процессов на границе упругого однородного полупространства с изотермической моделью однородной атмосферы посвящены работы А.В. Разина (1993), Л.А. Гасиловой и Ю.В. Петухова (1999). В них с использованием метода контурного интегрирования и результатов исследования дисперсионного уравнения для собственных решений построенной модели установлено существование поверхностной волны Стоунли-Шолтэ и модифицированной волны Лэмба.

Постановка задачи. Система линеаризованных уравнений, описывающая распространение акусто-гравитационных волн для модели неоднородной атмосферы в цилиндрической системе координат (r, z), имеет вид:

uz P ur P P 0 P 0 g; 0 ; c0 uz u ;

z t z t r t t z z (6) ur u 1 z 0 ur u F(r, z, t).

z t r r z z Здесь ur, u - компоненты скорости движения частиц воздуха в волне, z P, - соответственно, возмущения давления и плотности под действием источника массы F(z, r,t) (z z0) f (t) c0 (z) - скорость звука в ат, мосфере, константа g - ускорение силы тяжести. Нулевыми индексами обозначены величины, относящиеся к среднему невозмущенному состоянию атмосферы. Зависимость значений гидростатического давления P0 и плотности атмосферы от высоты определяется соотz ношением P0, 0 ~ exp(z / H ), где H - высота изотермической однородной атмосферы. Для покоящейся атмосферы в однородном поле силы тяжести имеют место формула: P0 g и (z) exp(z / H), где 1 - 0 0 плотность при z 0 (граница литосфера - атмосфера).

Распространение сейсмических волн для модели неоднородной литосферы в цилиндрической системе координат записывается в виде системы уравнений теории упругости первого порядка:

uz rz zz rz 0 0 Fz (r, z,t);

t r z r ur rr rz rr 0 0 Fr (r, z, t);

t r z r uz ur ur zz 2 ;

t z r r (7) ur ur ur rr 2 ;

t r z r ur uz ur 2 ;

t r z r ur uz rz .

t z r В этих уравнениях ur, u - компоненты скорости смещения упругой z среды, ,,, - компоненты тензора напряжений, z - плотzz rr rz ность упругой среды, z и z - коэффициенты Ламе. При этом Fz, Fr - составляющие силы F(z, r,t) Fz (z, r,t)ez Fr (z, r,t)er, описывающие распределение локализованного в пространстве источника типа лцентра давления, вертикальной силы или диполя без момента (если источник находится в упругой среде). Например, в случае источника ти па лцентра давления компоненты силы F(z, r,t) имеют вид:

(r) d d (r) Fz (z, r,t) (z z0 ) f (t); Fr (z, r, t) (z z0 ) f (t).

2r dz dr 2r Функция f (t) задат изменчивость источника по времени. Для формального описания источников используются выражения, содержащие обобщенные функции и производные от них. При численном решении для обобщенных функций вводятся соответствующие гладкие аналоги. Задачи (6) и (7) решаются при нулевых начальных данных в атмосфере и упругом полупространстве, а также и при условии согласования на границе литосфера - атмосфера:

P P0 zz uz лит uz атм; uz, 0. (8) rz лит t t z лит атм Алгоритм решения. Последовательное применение преобразований Фурье-Бесселя по радиусу и Лагерра по времени сводит задачу к решению системы линейных обыкновенных дифференциальных уравнений для набора правых частей. Для отыскания численного решения использовалась конечно-разностная схема на сдвинутых сетках, обеспечивающая второй порядок аппроксимации. В полученной системе алгебраических уравнений параметр разделения р (степень полинома Лагерра) присутствует только в правой части, а матрица системы от него не зависит. Это дат возможность использовать прямые методы решения, например, на основе разложения Холецкого, выполняемого только один раз. При этом решение системы строится одновременно для всех правых частей или, другими словами, сразу для всех значений параметра разделения. Этот процесс повторяется для всех корней уравнения Бесселя, что позволяет вернуться в пространственную область по радиальной переменной.

Спецификой поставленной задачи является необходимость моделирования совместных волновых процессов в атмосфере и упругом неоднородном полупространстве, скорости распространения волн в которых отличаются в 10-15 раз. Поэтому для изучения волновых явлений необходимо рассматривать волновые поля для больших промежутков времени, а следовательно, и для больших пространственных областей. Для ограничения расчетной области в атмосфере и упругом полупространстве вводятся специальным образом сконструированные достаточно узкие (меньше длины волны) поглощающие слои, учитывающие специфику задачи, связанную с использованием преобразования Лагерра (Решетова, Чеверда, 2006). С их помощью удатся обеспечить устойчивость алгоритма вне зависимости от длительности расчетов по времени.

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

Для понимания особенностей волновых процессов рассматривается однородное упругое полупространство, граничащее с атмосферой, в которой скорость звука постоянна, а плотность экспоненциально падает с высотой. Скорость распространения продольных волн в упругом полупространстве Vp 3000м / с, скорость поперечных волн Vs 1760м / с, плотность 0 2.3г / см3. Скорость распространения акустических волн в атмосфере c0 340м / с, плотность 0 0001225 г / смz (при ). Источник типа лцентра давления расположен в упругом полупространстве.

На рис. 11а приведен численный снимок в момент времени t 20с для горизонтальной компоненты скорости смещения в упругом полупространстве и давления в атмосфере. Источник типа лцентра давления расположен в полупространстве на расстоянии 3000м (, где - доминирующая длина волны P) от границы. В первых z вступлениях в упругом полупространстве регистрируется продольная P - волна от источника и отраженная PP - волна от границы с атмосферой и далее обменная - волна. В атмосфере регистрируется коPS ническая волна.

Волновая картина усложняется, если источник типа лцентра давления находится вблизи свободной поверхности на расстоянии 750м (1/ 4 ). На рис. 11б приведен численный снимок в тот же момент для той же самой горизонтальной компоненты скорости смещения в упругом полупространстве и давления в атмосфере. В упругом полупространстве, кроме классических волн, возникает нелучевая S*волна. Как видно из рисунка, она возникает и в атмосфере. Кроме того, по границе литосфера-атмосфера распространяется поверхностная волна Стоунли-Шолтэ, а в атмосфере - серия конических волн.

б) а) Рис.11. Численный снимок горизонтальной компоненты волнового поля в момент времени 20 с для двух положений источника: 3000 м (а) и 750 м (б).

ЗАКЛЮЧЕНИЕ Результатом исследования является разработка новых эффективных методов, алгоритмов и программного обеспечения для численного моделирования сейсмических и сейсмоакустических волновых полей в трхмерных неоднородных разномасштабных резкоконтрастных средах сложного строения.

1. На основе конечно-разностного подхода разработан и теоретически и экспериментально обоснован новый метод, создан оригинальный алгоритм и реализовано программное обеспечение для численного моделирования волновых полей акустического каротажа. Для максимально точного описания наиболее контрастной границы (скважина/вмещающая среда) используется цилиндрическая система координат с квазинеравномерным радиальным и адаптивным азимутальным измельчением разностной сетки. Их использование стало возможным, благодаря оригинальной процедуре согласования адаптивных сеток по азимуту на основе интерполяции с помощью преобразования Фурье.

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

- выполнение измельчения по пространству и по времени на различных поверхностях, охватывающих целевую область, что обеспечило устойчивость метода;

- использование для согласования сеток пространственной интерполяции на основе преобразования Фурье, что обеспечило е экспоненциальную точность и низкий уровень возникающих при этом артефактов (менее 0.1% от амплитуды падающей волны);

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

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

- свести решение исходной задачи к решению системы уравнений с отрицательно определнной матрицей, что обеспечивает уверенную сходимость итерационных методов;

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

- обеспечить возможность ограничения расчетной области устойчивыми PML без расщепления.

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

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

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

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

Основные публикации по теме диссертации 1. Конюх (Решетова) Г.В., Михайленко Б.Г., Михайлов А.А. Численное моделирование сейсмических полей в вязкоупругих средах на основе спектрального метода Лагерра // Математическое моделирование. -2001. -Т.13 (2). -С.61-70.

2. Ельцов И.Н., Кашеваров А.А., Решетова Г.В., Чеверда В.А.

Проявление неоднородностей зоны проникновения в геофизических полях вдоль ствола скважины // Геофизика. -2004. -Т.6. -С.17-21.

3. Ельцов И.Н., Решетова Г.В., Чеверда В.А. Численное моделирование процессов распространения сейсмоакустических полей с учетом неоднородности зоны проникновения // Физическая мезомеханика. -2005. -Т.8. -С.99-105.

4. Михайленко Б.Г., Решетова Г.В. Численно-аналитический метод решения задачи о распространении сейсмических и акустогравитационных волн для неоднородной модели Земля-Атмосфера // Сибирский журнал вычислительной математики. -2006. -Т.9. -С.37-46.

5. Михайленко Б.Г., Решетова Г.В. Математическое моделирование распространения сейсмических и акусто-гравитационных волн для неоднородной модели Земля-Атмосфера // Геология и геофизика, 2006, т.47(5), с. 547-556.

6. Решетова Г.В., Чеверда В.А. Использование преобразования Лагерра для построения идеально подходящих поглощающих слоев без расщепления // Математическое моделирование. -2006. -Т.18 (10). С.91-101.

7. Костин В.И., Решетова Г.В., Чеверда В.А. Численное моделирование трехмерного акустического каротажа с использованием многопроцессорных вычислительных систем // Математическое моделирование. -2008. -Т.20(9). -С.51-66.

8. Гилбо Ж., Ланда Е., Решетова Г.В., Хайдуков В.Г., Чеверда В.А.

Численное моделирование сейсмических волновых полей в двумернонеоднородных упругих разномасштабных средах (карстовые включения) // Технологии сейсморазведки. -2008. -Т.3. -С.19-28.

9. Konyukh (Решетова) G.V., Krivtsov Y.I., Mikhailenko B.G. Numerical-Analytical Modeling of Seismic Wave Propagation in Vertically Inhomogeneous Media//Applied Mathematics Letters. -1998. -V.11(1). P.99-104.

10. Konyukh (Решетова) G.V., Mikhailenko B.G., Mikhailov A.A.

Application of the integral Laguerre transforms for forward seismic modeling // Journal of Computational Acoustics. -2001. -V.9(4). -P.1523-1541.

11. Mikhailenko B.G., Mikhailov A.A., Reshetova G.V.Numerical viscoelastic modeling by the spectral Laguerre method // Geophysical Prospecting. -2003. -V.51. -P.37-48.

12. Mikhailenko B.G., Mikhailov A.A., Reshetova G.V.Numerical modeling of transient seismic fields in viscoelastic media based on the Laguerre spectral method // Pure and Applied Geophysics. -2003. -V.160. -P.

1207-1224.

13. Pissarenko D.V., Reshetova G.V., Tcheverda V.A. 3D finitedifference synthetic acoustic logging in cylindrical coordinates // Geophysical Prospecting. -2009. -V.57(3). -P. 367-377.

14. Mikhailenko B.G., Reshetova G.V. Mathematical modeling of seismic and acousto-gravitational waves in a Heterogeneous earthatmosphere model // Journal of Computational and Applied Mathematics. 2010. -V.234. -P. 1678-1684.

15. Pissarenko D.V., Reshetova G.V., Tcheverda V.A. 3D finitedifference synthetic acoustic log in cylindrical coordinates // Journal of Computational and Applied Mathematics. -2010. -V. 234. -P. 1766-1772.

Подписано к печати 03.11.20Формат 60х84/16. Объем 2 печ.л.

Тираж 100 экз. Зак. 172.

ОOO Омега Принт, пр-т. Ак. Лаврентьева, 6, Новосибирск 6300    Авторефераты по всем темам  >>  Авторефераты по техническим специальностям