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

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

ДМИТРИЕВ Максим Николаевич

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

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

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

18 АПР ¿013

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

005051968

005051968

Работа выполнена в Федеральном государственном бюджетном учреждении науки Институте нефтегазовой геологии и геофизики им. A.A. Трофимука Сибирского отделения Российской академии наук и в Федеральном государственном бюджетном образовательном учреждении высшего профессионального образования Новосибирском национальном исследовательском государственном университете.

Научный руководитель:

доктор физико-математических наук, доцент Роменский Евгений Игоревич

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

Решетова Галина Витальевна, доктор физико-математических наук, ведущий научный сотрудник Федерального государственного бюджетного учреждения науки Института вычислительной математики и математической геофизики Сибирского отделения Российской академии наук

Костин Виктор Иванович, кандидат физико-математических наук, доцент, ЗАО «Интел А/О» в г. Новосибирске, руководитель отдела по разработке программного обеспечения

Ведущая организация:

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

Защита состоится 30 апреля 2013 г. в 15 часов на заседании диссертационного совета Д 003.061.02 на базе Федерального государственного бюджетного учреждения науки Института вычислительной математики и математической геофизики Сибирского отделения Российской академии наук (ИВМиМГ СО РАН).

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

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

Автореферат разослан 25 марта 2013 г.

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

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

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

Актуальность.

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

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

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

- идеально согласованные поглощающие слои.

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

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

Научные задачи

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

2. Теоретически и экспериментально обосновать алгоритм ограничения расчетной области.

Фактический материал и методы исследования

Теоретической основой решения поставленных научных задач являются:

современная теория упругости, вязкоупругости и электроупругости;

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

теория нелинейной минимизации для аппроксимации добротности на заданном частотном диапазоне;

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

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

Защищаемые научные результаты:

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

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

3. Научно-исследовательский вариант программного обеспечения.

Научная новизна. Личный вклад.

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

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

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

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

Теоретическая и практическая значимость.

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

Апробация работы и публикации.

Основные положения и результаты, полученные автором, докладывались, обсуждались и одобрены специалистами на 9 конференциях: 10-й Международной научно-практической конференции «Геомодель-2008» (Геленджик, 2008), XLVII Международной научной студенческой конференции «Студент и научно-технический прогресс» (Новосибирск, 2009), XIII Всероссийской конференции-школе «Современные проблемы математического моделирования» (п. Абрау-Дюрсо, 2009), 4-й Международной конференции и выставке «Санкт-Петербург - 2010. К новым открытиям через интеграцию геонаук» (Санкт-Петербург, 2010), VI Международной выставке и научном конгрессе «ПЮ-Сибирь-2010» (Новосибирск, 2010), Международном

симпозиуме «Seismic Waves in Laterally Inhomogeneous Media VII» (Чешская Республика, Тепла, 2010), Молодежной международной научной школе-конференции «Теория и численные методы обратных и некорректных задач» (Новосибирск, 2010), Международной конференции «European Conference on High Order Nonlinear Numerical Methods for Evolutionary PDEs» (Италия, Тренто, 2011), XIV Всероссийской молодежной конференции-школе с международным участием «Современные проблемы математического моделирования» (п. Абрау-Дюрсо, 2011).

Результаты исследования по теме диссертации изложены в 12 опубликованных работах в том числе три статьи в ведущих научных рецензируемых журналах по перечню ВАК (Уфимский математический журнал [Дмитриев, Роменский, 2010], Сибирский журнал вычислительной математики [Дмитриев, Лисица, 2011, 2012]) и 9 материалах международных и российских конференций и симпозиумов.

Структура и объём работы.

Диссертация состоит из введения, 4 глав, заключения и списка литературы из 88 наименований. Общий объём диссертации составляет 115 страниц, включая 1 таблицу и 33 рисунка.

Благодарности.

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

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

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

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

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

Во второй главе рассматривается наиболее общая математическая модель, описывающая процессы формирования и распространения волновых полей в вязкоупругих средах. Для учета поглощения энергии используется обобщенная стандартная модель линейного твердого тела (английская аббревиатура GSLS — General Standard Linear Solid), основанная на использовании механических аналогов вязкоупругого тела [12]. Описывается конечно-разностный алгоритм решения полученной системы уравнений и приводятся результаты численных экспериментов.

В разделе 2.1 рассматривается одномерная модель GSLS [Christensen, 1982]. Показано, что исходная система уравнений может быть приведена к следующему виду:

Э и Эсг

^ dt Эх' Э<7 м

Г L \

1=1

ди . тгч _ Г„,-Г

_ _ "с/ *о7

дх jrf Г

дг, ,, ди

где введены дополнительные переменные памяти г(. Здесь и и С -скорость и напряжение соответственно, Те1 ,Та1,Мк — параметры среды. В качестве меры поглощения используется понятие добротность. Если уравнение состояния в частотной области переписать в виде & = Сё, то

добротность определяется как отношение {) =-—. В случае

1ш(С)

использования модели добротность является рациональной

1-L + У — £Г 1-

функцией частоты О(£0,Т„,Т„)=--.-——. При этом

я

возникает необходимость поиска времен релаксаций

,тє2Zel ) wTa — (Та1,Та2,...,TaL) таким образом, чтобы обеспечивалось минимальное отклонение от экспериментально измеренной функции добротности Qn(OJ):

F(те,) = JI ß"1 (to,тє, ) - Öo' (a>) I2 da^ min.

Щ

Показано, что приемлемая точность достигается при использовании трёх релаксационных механизмов. На основе исследования создана база данных времен релаксаций для широкого диапазона значений добротности. На рис. 1 показаны точная и восстановленная добротности для интервала частот от 2 до 30 кГц. Можно видеть, что на всем интервале частот относительная ошибка не превосходит 1.2%. Для проверки корректности описания амплитуды в зависимости от добротности проведена серия численных экспериментов.

Рис. 1. Точная (синим) и восстановленная (красным) функция добротности (слева). Относительная ошибка в процентах (справа).

В разделах 2.2 - 2.3 выполняется обобщение одномерной модели для трехмерного случая в декартовой и цилиндрической системах координат.

В разделе 2.4 описывается конечно-разностный алгоритм моделирования сейсмоакустических волновых полей с поглощением для трехмерных неоднородных сред с осевой симметрией в цилиндрической системе координат. Для численного решения полученной системы используется схема на сдвинутых сетках [Утеих, 1986; Костин, Решетова, Чеверда, 2008]. Основной особенностью этой схемы является

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

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

В разделах 3.1 - 3.2 приводится математическая модель, описывающая распространение волновых процессов в электроупругих средах [Kostec, Randall, 1994]. Рассмотрен один из основных классов пьезоэлектриков - PZT (цирконат-титанат свинца). Показано, что для среды, обладающей аксиальной симметрией, система уравнений, записанная в цилиндрической системе координат имеет следующий вид:

Эг Эг Эг

~ О,

г

dt Эг Э z г

+ —!-Г

Эсг . Э ш

дет, . дш

э?

где введены следующие обозначения: г Эгг £ и, Е Эи.

Д = 4 V-+4—+

с)г г ог

. Е ди, Е и. £ Эи

= 4 "эГ+4 4

. £ ЭИГ £ М £ ЭК.

Л <т = г —L 4- г — + с —-

ог г ог £ | Эм„ ди

Лст„=с44| ■

ог дг

Здесь (мг, м.) и (<7,т, <7^, <та, <Тп) - компоненты вектора скорости и

тензора напряжений соответственно, у/ - производная по времени от

£ £ £ £

электрического потенциала, р - плотность, си, с12, с13, с44 -коэффициенты тензора жесткости, £3|,е33,е15 - компоненты пьезоэлектрического тензора. Для нахождения функции \{/ решается уравнение Пуассона:

\ д( с ач-и г ат

--ге33 — ' ^

дг \ дг

г

+

\

11 Эг

где 8 +<*33А(Тгг)] + ^-[^5А(7гА >

г

й?31, ¿/15 и е^,е33 - компоненты пьезоэлектрического и диэлектрического тензора соответственно. На внутренней и внешней поверхностях пьезоэлектрического цилиндра ставятся граничные условия, задающие режим колебаний.

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

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

В разделе 3.4 представлен результат численного моделирования, полученный с использованием описанного алгоритма. В качестве модели среды была взята скважина, заполненная буровым раствором, а источник представлялся в виде пьезоэлектрического цилиндра с высокопоглощающим вязкоупругим заполнением. На рис. 2 представлена детальная схема численного эксперимента. В качестве граничных условий для потенциала выбран сигнал в виде функции Риккера с доминирующей частотой 10 кГц. На рис. 3 представлены снимки волнового поля в различные моменты времени для точечного источника (слева) и эксперимента, описанного выше (справа). Можно видеть существенное влияние конструкции источника на волновое поле.

„_1м_I

L_Q.2M_.

--!-1--Г

„OJK

—+ —

Зм

Рис. 2. Схема численного эксперимента

а) б)

1 = 0.00047512с

В) Г)

1 = 0.00095023с

0.25

0.2 0.4 0.6 0.8 Стабилизирующий параметр ['>

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

падения, = 0.01 .

15 30 45 60 75 90 Угол падения в градусах

Рис. 5. Зависимость аналитического и экспериментального коэффициентов отражения от угла падения а.

-0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 г (м) г (м) г(м) г(м)

Рис. 3. Снимки волнового поля (<7ГГ) для двух фиксированных моментов

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

Г 1.5

В четвертой главе описывается способ построения слабоотражающих граничных условий для ограничения расчетной области при моделировании волновых процессов в неограниченных областях. Рассматривается подход, известный как M-PML (от англ. Multiaxial Perfectly Matched Layer) [Kristel С. Meza-Fajardo et al., 2008] для гиперболической системы общего вида:

= (1)

at dx, dx2 M-PML в направлении х, для системы (1) записывается в виде системы:

Эм1 ,i л д и

-+ dxu - Д-= 0,

dt Эх,

(2)

ди2 ,2 . ди п — + d2u-A2—~ = 0, at ах2

где и, + и2 = и , а параметры d¡ и dn = /3d¡ есть демпфирующие функции, зависящие от х,, /3 £ [0,1] - стабилизирующий параметр. Показано, что M-PML не является идеально согласованным и обладает большей отражательной способностью по сравнению с классическим PML, при этом коэффициент отражения линейно зависит от значения стабилизирующего параметра. На основе этого исследования формулируется задача построения оптимального стабилизирующего параметра, обеспечивающего устойчивость и минимальные отражения. Получен необходимый признак устойчивости M-PML, позволяющий ограничить нижние значения стабилизирующего параметра. Выполненные численные эксперименты показывают, что этот признак не является достаточным.

В разделе 4.1 исследуется отражательная способность M-PML для системы уравнений акустики для дифференциальной и конечно-разностной постановки [2]. Предположим, что система уравнений акустики задана при х, < 0, в то время как M-PML введен в области х, > 0. Рассматривая падения плоских волн на границу раздела при х, = 0, получен теоретический коэффициент отражения в зависимости от угла падения а :

. 8ПГ(0Г) . .

&т2(а) . ч

+ СОЗ(«Г)

V о-^)2

где = с1 ■ / СО, СО - круговая частота. Аналогично получены

теоретические коэффициенты отражения для конечно-разностной постановки для схемы на сдвинутых сетках. На рис. 4 представлены абсолютные значения коэффициентов отражения в зависимости от стабилизирующего параметра для заданного угла падения и параметра . Можно видеть, что для типичных значений демпфирующих функций коэффициенты линейно зависят от стабилизирующего параметра. Для верификации полученных формул проведена серия численных экспериментов: падение плоской волны на границу раздела регулярная область — М-РМЬ для серии углов. На рис. 5 показаны аналитически и экспериментально полученные коэффициенты отражения для различных углов. Из рисунка видно, что теоретические результаты с достаточно хорошей точностью совпадают с экспериментальными.

В разделе 4.2 получено необходимое условие устойчивости М-РМЬ для произвольной гиперболической системы уравнений [3]:

Для устойчивости М-РМЬ (2) для гиперболической системы уравнений (1) необходимо выполнение следующего условия:

УКе Я2, У1(К)5,(К) + /ЗУ2(К)Б2(К)> О,

где БJ и V. - компоненты вектора медленности и вектора групповой

скорости соответственно, К - единичный волновой вектор.

В разделе 4.3 проводится численное исследование устойчивости М-РМЬ [7-10]. Рассматриваются две анизотропные упругие модели среды, которые являются общепризнанными моделями для тестирования слабоотражающих граничных условий. На рис. 6 показаны снимки волнового поля для УТ1 модели, характерной для глинистых горных пород. Можно видеть, что внутри рассматриваемого промежутка времени М-РМЬ является устойчивым, тогда как классический РМЬ является неустойчивым. В конце раздела приводится эксперимент для неоднородной изотропной среды, для которой классический РМЬ является неустойчивым.

1=0.20187с

1=0.64993с

300 600 900 1200

х1 (м)

1=0.20187с

1500 Щ

0 300 600 900 1200 х1 (м)

1=0.28065с

0 300 600 900 1200

х, (м)

1=0.64993с

0 300 Г \

600 шявш ,

900 Й тг » ш

1200 ЯЙ!'' 4 г-

1500 га У

0 300 600 900 1 200

Х1 (м)

0 300 600 900 1200 X, (м)

0 300 600 900 1 200

X, (м)

Рис. 6. Снимки волнового поля (и,) для модели глинистых пород. Сверху -классический РМЬ. Снизу - М-РМЬ со стабилизирующими параметрами Д = 0.0883 , Д = 0.0724 в направлениях Л", и х2 соответственно.

ЗАКЛЮЧЕНИЕ

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

применении современных методов вычислительной и прикладной математики.

Для ограничения расчетной области в работе детально исследован алгоритм построения слабоотражающих граничных условий М-РМЬ. Рассмотренный подход является устойчивым для произвольной неоднородной и анизотропной среды при соответствующем выборе стабилизирующего параметра. На основе исследования сформулирована задача построения оптимального стабилизирующего параметра, т.е. минимизирующего отражения от границы М-РМЬ - расчетная область и в то же время обеспечивающего устойчивость. Автором получен необходимый признак устойчивости, позволяющий оптимизировать выбор поглощающих слоев. Численными экспериментами показано, что данный признак не является достаточным. Предложенный подход обладает рядом преимуществ: простота построения, достаточно высокая эффективность, возможность производить расчет волнового поля без потери устойчивости.

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

СПИСОК ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ

Статьи в рецензируемых журналах, рекомендованных ВАК

1. Дмитриев М.Н., Роменский Е.И. \УЕ1ЧО/Рунге-Кутга метод высокой точности для моделирования упругих волн // Уфимский математический журнал. - 2010. - Т.2. - № 1. - С. 59-70.

2. Дмитриев М.Н., Лисица В.В. Применение слабоотражающих граничных условий М-РМЬ при моделировании волновых процессов в анизотропных средах. Часть I: Уровень отражений // Сиб. журн. вычисл. математики. - 2011. - Т. 14, №4. - С. 333-345.

3. Дмитриев М.Н., Лисица В.В. Применение слабоотражающих граничных условий М-РМЬ при моделировании волновых процессов в анизотропных средах. Часть И: Устойчивость // Сиб. журн. вычисл. математики. - 2012. - Т. 15, №1. - С. 45-55.

Материалы конференций

4. Дмитриев М.Н., Лисица В.В. Сравнительный анализ схем повышенной точности для системы уравнений акустики // Материалы

10-й международной научно-практической конференции «Геомодель-2008» (Геленджик, 21-26 сентября 2008 г.). - 4 с. - Электронный ресурс. - 1 электрон, опт. диск.

5. Дмитриев М.Н. Схемы повышенной точности для моделирования волновых процессов в упругих средах // XLVII Международная научная студенческая конференция «Студент и научно-технический прогресс» (Новосибирск, 9-14 апреля 2009 г.). - Новосибирск, 2009. - С. 221.

6. Дмитриев М.Н. Схемы повышенной точности для расчета упругих волн // Сборник трудов XIII всероссийской конференции-школы «Современные проблемы математического моделирования» (Абрау-Дюрсо, 14 - 19 сентября 2009 г.). - Ростов-на-Дону, 2009. -С. 205-211.

7. Dmitriev M.N., Lisitsa V.V. Stability and reflectivity of M-PML for anisotropic elastic media // «4th international Conference & Exhibition Saint Petersburg 2010» (St. Petersburg, 5-8 April 2010). - Electronic source.

8. Дмитриев M.H., Лисица B.B. Устойчивость и отражающие свойства M-PML для анизотропных упругих сред // VI Международная выставка и научный конгресс «ГЕО-Сибирь-2010» (Новосибирск, 19-29 апреля 2010 г.).- Новосибирск, 2010. - Т. 2,4.1. - С. 97-101.

9. Dmitriev M.N., Lisitsa V.V. Stability and reflectivity of M-PML for anisotropic elastic media// Workshop meeting «Seismic Waves in Laterally inhomogeneous media VII» (Chech Republic, Tepla, 21-26 June 2010). -Режим доступа свободный:

http://ccom.ucsd.edu/~dd20/downloads/dd20-abstracts.pdf.

10. Дмитриев М.Н., Лисица B.B. Ограничение расчетной области с помощью M-PML: устойчивость и уровень артефактов // II Международная школа-конференция «Теория и численные методы решения обратных и некорректных задач» (Новосибирск, 21-29 сентября 2010 г.). - Режим доступа свободный: http://math.nsc.ru/conference/onzlO/thesis/abstracts.pdf.

11. Dmitriev M.N., Romenski E.I. High accuracy Runge-Kutta WENO schemes for numerical simulation of seismic wave propagation // «European Conference on High Order Nonlinear Numerical Methods for Evolutionary PDEs: Theory and Applications» (Italy, Trento, 11-15 April 2011). -Electronic Source.

12. Дмитриев М.Н. Численное моделирование сейсмоакустических процессов // Сборник трудов XIV всероссийской конференции-школы «Современные проблемы математического моделирования» (п. Абрау-Дюрсо, 12-17 сентября 2011). - Ростов-на-Дону, 2011. - С. 118-122.

_Технический редактор Е.В.Бекренёва_

Подписано в печать 19.03.2013 Формат 60x84/16. Бумага офсет №1. Гарнитура Тайме _Печ.л. 0,9. Тираж 120. Зак. № 84_

ИНГГ СО РАН, просп. Акад. Коитюга, 3, Новосибирск, 630090

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ

УЧРЕЖДЕНИЕ НАУКИ ИНСТИТУТ НЕФТЕГАЗОВОЙ ГЕОЛОГИИ И ГЕОФИЗИКИ

им. A.A. ТРОФИМУКА СИБИРСКОГО ОТДЕЛЕНИЯ РОССИЙСКОЙ АКАДЕМИИ НАУК

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ

ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ НОВОСИБИРСКИЙ НАЦИОНАЛЬНЫЙ

ИССЛЕДОВАТЕЛЬСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

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

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

ДМИТРИЕВ Максим Николаевич

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

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель: д.ф.-м.н. Е.И. Роменский

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

ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ .............................. 4

Глава 1. ИЗУЧЕННОСТЬ РЕШЕНИЯ ЗАДАЧИ ....... 11

Глава 2. КОНЕЧНО-РАЗНОСТНОЕ МОДЕЛИРОВАНИЕ ВОЛНОВЫХ ПРОЦЕССОВ

В ВЯЗКОУПРУГИХ СРЕДАХ ............. 19

2.1 Одномерная обобщенная стандартная

линейная модель твердого тела................. 20

2.1.1 Введение переменных памяти .............. 22

2.1.2 Аппроксимация добротности

на заданном частотном диапазоне............ 23

2.2 Обобщенная стандартная линейная модель

твердого тела в трехмерном случае............... 29

2.3 Уравнения вязкоупругости

в цилиндрической системе координат.............. 36

2.4 Построение конечно-разностного алгоритма

в цилиндрической системе координат.............. 39

2.5 Численные эксперименты .................... 46

Глава 3. КОНЕЧНО-РАЗНОСТНОЕ МОДЕЛИРОВАНИЕ

ПЬЕЗОЭЛЕКТРИЧЕСКОГО ИСТОЧНИКА .... 52

3.1 Уравнения электроупругости. Общие свойства........ 53

3.2 Определяющие соотношения

в цилиндрической системе координат.............. 58

3.3 Построение конечно-разностного алгоритма

в цилиндрической системе координат............................66

3.4 Численные эксперименты ........................................69

Глава 4. ПОСТРОЕНИЕ СЛАБООТРАЖАЮЩИХ

ГРАНИЧНЫХ УСЛОВИЙ М-РМЬ......................73

4.1 Отражающие свойства М-РМЬ..................................74

4.1.1 Дифференциальная постановка..........................75

4.1.2 Конечно-разностная постановка..........................79

4.2 Устойчивость М-РМЬ..............................................87

4.3 Численные эксперименты ........................................90

ЗАКЛЮЧЕНИЕ............................103

СПИСОК ЛИТЕРАТУРЫ .....................105

ВВЕДЕНИЕ

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

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

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

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

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

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

2. Теоретически и экспериментально обосновать алгоритм ограничения расчетной области.

Теория и методы исследования

Теоретической основой решения поставленной научной задачи являются:

- современная теория упругости, вязкоупругости и электроупругости;

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

- теория нелинейной минимизации для аппроксимации добротности на заданном частотном диапазоне;

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

Основной метод исследования - математическое моделирование сейсмо-акустических волновых полей в неоднородных средах с учетом конструкции источника. При разработке численного алгоритма для электроупругой среды использовалась математическая библиотека Intel Math Kernal Library для решения системы линейных алгебраический уравнений. Для разработки научно-исследовательской версии программного обеспечения использовался язык программирования Fortran.

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

Защищаемые научные результаты:

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

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

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

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

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

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

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

Теоретическая и практическая значимость результатов.

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

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

Апробация работы и публикации. Основные положения и результаты, полученные автором, докладывались, обсуждались и были одобрены специалистами на 9 конференциях: 10-й Международной научно-практической конференции ,,Геомодель-2008"(Геленджик, 2008), XLVII Международной научной студенческой конференции "Студент и научно-технический прогресс" (Новосибирск, 2009), XIII Всероссийской конференции-школе "Современные проблемы математического моделирования "(п. Абрау-Дюрсо,

2009), 4-й Международной конференции и выставке "Санкт-Петербург -2010. К новым открытиям через интеграцию геонаук" (Санкт-Петербург,

2010), VI Международной выставке и научном конгрессе "ГЕО-Сибирь-2010"(Новосибирск, 2010), Международном симпозиуме "Seismic Waves in Laterally Inhomogeneous Media VII"(Чешская Республика, Тепла, 2010), Молодежной международной научной школе-конференции "Теория и численные методы решения обратных и некорректных задач"(Новосибирск, 2010), Международной конференции "European Conference on High Order Nonlinear Numerical Methods for Evolutionary PDEs" (Италия, Тренто, 2011), XIV Всероссийской молодежной конференции-школе с международным участием "Современные проблемы математического моделирования"(п. Абрау-Дюрсо,

2011).

Результаты исследования по теме диссертации изложены в 12 опубликованных работах в том числе 3 статьи в ведущих научных рецепзиру-

емых журналах, рекомендованных ВАК (Уфимский математический журнал [Дмитриев, Роменский, 2010], Сибирский журнал вычислительной математики [Дмитриев, Лисица, 2011, 2012]) и 9 материалах международных и российских конференций и симпозиумов.

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

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

Глава 1

ИЗУЧЕННОСТЬ РЕШЕНИЯ ЗАДАЧИ

В настоящее время численное моделирование получило весьма широкое распространение как в научных исследованиях, так и при решении практических задач, связанных с поиском и разведкой месторождений полезных ископаемых. К одной из таких задач относится акустический каротаж - метод геофизических исследований, основанный на изучении свойств горных пород, пересеченных скважиной путем анализа волновых полей, создаваемых источником, находящимся в скважине. Акустический каротаж проводят с помощью глубинного прибора, основными элементами которого являются излучатели и приёмники упругих воли, а также акустические изоляторы, предотвращающие распространение упругих волн по корпусу прибора [Ивакин, Карус, Кузнецов, 1987]. Впервые данная задача в простейшей постановке была исследована в работе Био [Вю^ 1952]. С помощью интегральных преобразований он детально исследовал волновые поля, вызванные действием точечного источника, помещенного в скважину. Развитием этих исследований является работа [Крауклис П.В., Крауклис Л.А., 1976], где проведено теоретическое исследование ряда осесимметричных задач и получены дисперсионные соотношения для основных типов волн. Однако названные исследования ограничиваются рассмотрением простейших моделей. В действительности, в большинстве случаев не представляется возможным решить поставленную задачу аналитически и единственный путь - численное решение задачи.

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

численных методах: метод конечных разностей [Лебедев, 1964; Virieux, 1986; Levancler 1988; Saenger, Gold et al, 2000; Lisitsa, Vishnevsky, 2010]; метод конечных объемов [Zhang, 1997; Zhang, 2004]; метод конечных элементов [Bathe, 1979; Moczo, Bystricky et al., 1997]; метод спектральных элементов [Chaljub, Komatitsch et al, 2007; Tromp, Komatitsch et al., 2008]; разрывный метод Га-леркина [Grote, Schneebeli et al. 2006; Käser, Dumbser, 2006; Puente, Dumbser et al., 2008; Etienne, Chaljub et al. 2010]. Методы конечных элементов, спектральных элементов и разрывный метод Галеркина основаны на дискретизации уравнений в вариационной постановке и имеют определенные преимущества для моделирования волновых процессов в областях со сложной геометрией, однако требуют существенных затрат на построение расчетной сетки. В настоящее время, пожалуй, не существует универсальных алгоритмов для построения сеток в достаточно сложных моделях. В отличие от перечисленных выше метод конечных разностей является более универсальным, так как дискретизация проводится на регулярной сетке. В последнее время широкое распространение получили так называемые сеточно-характеристические методы [Магомедов, Холодов, 1988; Кучунова, Садовский, 2008], основанные на расщеплении уравнений по пространственным переменным и решении одномерных задач методом характеристик (задач о распаде разрыва). Для повышения точности используют полиномиальную интерполяцию для нахождения потоков через границы ячеек. В работе [Дмитриев, Роменский, 2010] представлен сравнительный анализ схем данного типа с классическими конечно-разностными схемами на сдвинутых сетках для решения системы уравнений линейной теории упругости.

Среди последних следует отметить работу В.А. Чеверды с соавторами [Костин, Решетова, Чеверда, 2008], в которой на основе копечпо-разиостпого

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