Моделирование распространения сейсмических волн в флюидонасыщенных средах

Алексей Шевченко

Моделирование распространения сейсмических волн в флюидонасыщенных средах

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

Информатика

Дипломы

Вуз: Московский физико-технический институт (государственный университет) (МФТИ)

ID: 5f2d017ecd3d3e0001b48949

UUID: d8938ff0-baac-0138-1361-0242ac180006

Язык: Русский

Опубликовано: больше 3 лет назад

Просмотры: 5

10.1

Алексей Шевченко

Московский физико-технический институт (государственный университет) (МФТИ)

Комментировать 11

Рецензировать 0

Скачать - 1,4 МБ

Поделиться работой

Федеральное государственное автономное образовательное учреждение высшего образования «Московский физико-технический институт (национальный исследовательский университет)» Физтех-школа Аэрокосмических Технологий Кафедра логистических систем и технологий Направление подготовки / специальность: 27.03.03 Системный анализ и управление (бакалавриат) Направленность (профиль) подготовки: Теория и математические методы системного анализа и управления в технических, экономических и социальных системах МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ СЕЙСМИЧЕСКИХ ВОЛН В ФЛЮИДОНАСЫЩЕННЫХ СРЕДАХ (бакалаврская работа) Студент: Шевченко Алексей Владимирович _________________________ (подпись студента) Научный руководитель: Голубев Василий Иванович, канд. физ.-мат. наук, доц. _________________________ (подпись научного руководителя) Консультант (при наличии): _________________________ _________________________ (подпись консультанта) Москва 2020

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

Содержание 1 Введение 4 2 Описание модели 6 3 Численный метод 8 3.1 Одномерный случай . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Двумерный случай на прямоугольных сетках . . . . . . . . . . . . . . . 10 3.3 Двумерный случай на криволинейных сетках . . . . . . . . . . . . . . . 11 3.4 Граница . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.5 Контакт . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4 Результаты расчётов 8 21 4.1 Задача нагружения пористой среды . . . . . . . . . . . . . . . . . . . . . 21 4.2 Задача морской сейсмической разведки . . . . . . . . . . . . . . . . . . . 22 4.3 Задача прискважинной диагностики . . . . . . . . . . . . . . . . . . . . 25 5 Заключение 27 Литература 31 3

1 Введение На данный момент мировое потребление нефти и газа высоко. Истощение существующих резервуаров делает задачу поиска и разведки их месторождений актуальной. Основным способом решения этой задачи является сейсморазведка, использующая упругие волны и особенности их распространения в геологических массивах. Однако интерпретация результатов полевых работ представляет собой сложную и трудоёмкую задачу, при решении которой могут быть использованы результаты численного решения прямых задач. Для решения обратных задач сейсмической разведки необходимо иметь возможность быстро и точно решать прямые задачи распространения сейсмических волн. Развитие вычислительной техники и, в частности, высокопроизводительных вычислительных систем позволяет эффективно решать всё более и более сложные задачи. Одним из направлений «усложнения», помимо увеличения расчётной области и уменьшения пространственного шага сетки, является выбор более сложных физикоматематических моделей, точнее описывающих геологическую среду. Так, от акустических уравнений можно перейти к теории линейной упругости, хорошо зарекомендовавшей себя при описании многих геологических сред. Тем не менее, возникает необходимость в описании сред с несколько иными характеристиками: анизотропные [1], трещиноватые [2], насыщенные среды [3]. Последние естественным образом возникают в задачах, связанных с непосредственным моделированием нефтенасыщенных пластов, пористость которых может быть, например, в пределах 10–30 % для Баженовской свиты; также может быть полезен учёт водонасыщенности донных осадков в задачах моделирования шельфа, в том числе в Арктике. Именно моделирование процесса распространения упругих волн в насыщенных геологических средах является темой данной работы. Ещё в середине 20-го века советским учёным Я. Френкелем [4] была предложена модель пористой насыщенной среды, в которой распространяются две продольных волны. Но модель была узкоспециальной и не получила широкого распространения. В 1951 году была предложена модель, описывающая пористую насыщенную среду путём пересчёта осреднённых упругих параметров [5]. Благодаря незначительной дополнительной вычислительной сложности по сравнению с моделью линейной 4

упругости, она завоевала достаточно большую популярность, хоть она и не способна представить качественно новые результаты. В 1956 году М. Био, ссылаясь на вышеупомянутую работу Френкеля, предложил модель двухфазной пористой флюидонасыщенной среды, впоследствии получившую его имя. В связи с ограничениями вычислительной техники модель Био долгое время оставалась теоретической, хотя исследователи подтверждали корректность его уравнений, выводя их другими способами, например, основываясь на микроструктуре пор [6] или более строгой технике осреднения по объёму [7]. С развитием вычислительной техники модель стала весьма широко использоваться в различных задачах, например, моделировании костей человека [8], моделировании акустических волн в снежном покрове [9], моделировании поверхностных волн в геологической среде [10]. В 1989 году В. Доровским была предложена нелинейная континуальная теория фильтрация, развитая в дальнейшем в работе [11]. Её линеаризация стала известна как модель Доровского. В работе [12] было проведено сравнение этой модели с моделью Био-Джонсона (последняя является одной из модификация исходной модели Био) на основе экспериментальных данных. Было показано, что обе модели демонстрируют схожие результаты и достаточно точно моделируют реальную насыщенную среду. Модель Доровского успешно применялась для прямых и обратных задач сейсмики [13], для моделирования волн Стоунли в пористой насыщенной формации вокруг скважины [14]. Целью работы было численное исследование процесса распространения динамических возмущений во флюидонасыщенных средах. Были поставлены следующие задачи: • исследовать аналитически определяющую систему уравнений; • адаптировать сеточно-характеристический метод для решения уравнений модели Доровского; • разработать граничные условия и алгоритм их постановки; • разработать контактные условия и алгоритм их постановки; • реализовать программный код для моделирования динамики флюидонасыщенных сред; 5

• получить численные решения нескольких типов задач, таких как: задача о нагружении пористой среды, задача о морской сейсмической разведке, задача о сейсмическом мониторинге прискважинной зоны. Результаты работы были представлены на конференции, проводимой European Association of Geoscientists and Engineers (EAGE) — Геомодель-2019 [15], конференции МФТИ-62, а также опубликованы в статьях [16] и [17], входящих в список ВАК, Scopus и WoS. 2 Описание модели Рассматривается изотропная упругая двухфазная среда, состоящая из твёрдого скелета и насыщенных флюидом взаимосвязанных пор. Предполагается, что в каждом достаточно малом объёме есть и упругий скелет, и насыщающая жидкость. Для моделирования такой среды используется линеаризованная система уравнений, известная как континуальная теория фильтрации или модель Доровского [11]: ~ut + 1 1 (∇ · h)T + ∇p = F~ , ρs ρ0 1 ∇p = F~ , ρ0 ρs ρf T ht + µ ∇ ⊗ ~u + (∇ ⊗ ~u) + λ − K (∇ · ~u) − K(∇ · ~v ) I = 0, ρ0 ρ0 ~vt + pt − (K − αρ0 ρs )(∇ · ~u) + αρ0 ρf (∇ · ~v ) = 0. (1) (2) (3) (4) Здесь неизвестными функциями координат и времени являются: вектор ~u — скорость твёрдой фазы (скелета), вектор ~v — скорость флюида, h — тензор, противоположный по знаку тензору напряжений в твёрдом каркасе σ, p — давление флюида. Правая часть образована вектором приложенной внешней силы F~ . В систему входят параметры: ρs = (1 − β)ρs0 , ρf = βρf 0 , ρ0 = ρs + ρf , где ρs0 , ρf 0 — физические плотности скелета и флюида соответственно, а β — пористость среды, определяемая как отношение объёма пор ко всему объёму. Также в систему входят K, µ, α — упругие параметры насыщенной пористой двухфазной среды, которые однозначно могут быть выражены через 3 скорости рас6

пространения волн в пористой среде (о них будет сказано подробнее далее): 2 скорости продольных волн cp1 , cp2 (cp1 > cp2 ) и скорость поперечных волн cs : µ = ρs c2s , ρ0 ρs 2 8 ρf 2 K= cp1 + c2p2 − c − 2ρf 3 ρ0 s s 8 ρs 2 1 α3 = 2 c2p1 + c2p2 − c + 2ρ0 3 ρ0 s s (5) ! 64 ρ ρ f s 4 (c2p1 − c2p2 )2 − c , 9 ρ20 s (6) ! 64 ρ ρ f s (c2p1 − c2p2 )2 − c4 , 9 ρ20 s (7) α = ρ0 α3 + K/ρ30 . (8) Модель Доровского, так же, как и модели Роменского [3], Френкеля [4], Био [5] учитывает три скорости распространения волн (две скорости продольных волн, одна скорость поперечных). Вторая, медленная, продольная волна (добавочная по сравнению, например, с линейной теорией упругости) появляется за счёт двухфазности среды. Система линейных уравнений в частных производных, определяющая модель Доровского, является гиперболической. В случае любой размерности пространства (1D, 2D, 3D) матрица при производной по любому направлению (подразумевается каноническая форма записи) имеет полный набор собственных векторов, при этом абсолютные значения собственных чисел являются скоростями распространения волн. Модель Доровского построена на базовых физических принципах: законах сохранения, квазилинейности уравнения движения жидкости, инвариантности относительно преобразований Галилея, согласованности с условиями термодинамического равновесия. Законы сохранения энергии и массы упругого тела являются следствием приведённых уравнений. Подобное построение модели позволяет «связывать» её с другими физическими теориями, например, электромагнетизмом, как было сделано в работе [18]. Интересным отличием модели Доровского от модели Био является наличие всего трёх упругих параметров среды (K, µ, α), которые взаимно однозначно связаны со скоростями волн. Таким образом, для определения этих параметров достаточно экспериментального нахождения вышеуказанных скоростей. 7

3 Численный метод Поскольку система уравнений Доровского является гиперболической, для её численного решения может быть использован сеточно-характеристический метод, хорошо зарекомендовавший себя в задачах, связанных с моделированием геологических сред. Использование этого метода позволяет моделировать области, которые включают в себя среды, описываемые различными уравнениями, в рамках единого подхода и численного метода. Для нахождения отдельных формул аналитически использовалось программное обеспечение Wolfram Mathematica. Все описанные ниже методы для модели Доровского были реализованы в программном комплексе RECT Лаборатории Прикладной Вычислительной Геофизики МФТИ, написанном на С++. Эти методы были распараллелены на основе функционала RECT c использованием технологий MPI и OpenMP. 3.1 Одномерный случай В одномерном случае систему уравнений Доровского можно записать следующим образом: ~qt + A~qx = f~. (9) T Здесь ~q = u v hxx p — вектор неизвестных функций, компонентами которого являются скорость твёрдой фазы u, скорость жидкой фазы v, единственный в одномерном случае элемент минус тензора напряжений hxx и давление флюида p. f — правая часть, определяемая приложенной внешней силой. Матрица A определяется параметрами среды:  0 0    0 0 A=  4/3µ + Kρf /ρ0 −Kρf /ρ0  −K + αρ0 ρs αρ0 ρf 1/ρs 1/ρ0 0 0 0    1/ρ0  .  0   0 (10) Рассмотрим случай однородной постановки задачи, то есть при нулевой правой части. Гиперболичность системы уравнений означает, что матрица A имеет полный набор собственных векторов. За счёт перехода к базису из собственных векторов 8

можно получить спектральное разложение матрицы: A = Ω−1 ΛΩ, (11) где столбцами матрицы Ω−1 являются собственные векторы матрицы A, Ω·Ω−1 = I — единичная матрица, Λ = diag(λ1 , . . . , λn ) — диагональная матрица из собственных чисел, расположенных в том же порядке, что и собственные векторы в столбцах матрицы Ω−1 . Для одномерного случая рассматриваемой системы уравнений Доровского Λ = diag(−cp1 , cp1 , −cp2 , cp2 ). (12) Подставляя (11) в (9) и домножая обе части на матрицу Ω слева, получим Ω~qt + ΛΩ~qx = ~0. (13) Поскольку матрица Ω определяется параметрами среды, которые мы можем считать постоянными во всей области, можно вносить матрицу под знак дифференциала. Тогда введём замену ω ~ = Ω~q, где ω ~ — вектор инвариантов Римана. В новых переменных ω ~ система примет вид ω ~ t + Λ~ωx = 0. (14) Таким образом, поскольку Λ — диагональная матрица, в терминах инвариантов система дифференциальных уравнений в частных производных вырождается в набор независимых уравнений переноса вида wt + λwx = 0. (15) Для нахождения значения каждого инварианта на следующем шаге по времени используется постоянство искомой функции на характеристиках с помощью формулы ωi (t + dt, x) = ωi (t, x − λi · dt), (16) где значение ωi (t, x − λi · dt) определяется путём интерполяции требуемого порядка. Таким образом, общий алгоритм решения следующий: зная начальные значения, выполняем итеративно нужное число шагов по времени dt, при этом на каждом шаге: 1. переходим к инвариантам Римана; 9

2. находим новые значения инвариантов с помощью интерполяции; 3. выполняем обратную замену — переход к исходным переменным. Шаг по времени dt выбирается из условия Куранта: dt < h λmax , (17) где h — шаг сетки, λmax — максимальное собственное число. Для быстроты и точности расчёта матрицы Ω, Ω−1 были найдены аналитически. 3.2 Двумерный случай на прямоугольных сетках В двумерном случае система может быть записана в виде ~qt + Ax ~qx + Ay ~qy = f~. (18) T — вектор неизвестных, состоящий из Здесь ~q = ux uy vx vy hxx hyy hxy p компонент скоростей скелета и флюида, компонент тензора h, давления p; f~ — правая часть, определяемая внешней силой; матрицы Ax , Ay определяются параметрами среды согласно выражениям ниже:  0    0    0    0 Ax =   4  3 µ + ρρf K 0   2  − 3 µ + ρρf K 0    0  −K + αρ0 ρs 0 0 0 1 ρs 0 0 0 0 0 0 0 1 ρs 0 0 0 0 0 1 ρs 0 0 0 0 0 0 0 − ρf0 K 0 ρ 0 0 0 0 − ρf0 K 0 ρ 0 0 0 µ 0 0 0 0 0 0 0 0 0 0 αρ0 ρf 10  1 ρ0   0   0   0 ,  0   0   0  0 (19)

 0 0   0 0   0 0   0 0 Ay =    0 − 23 µ + ρρf K 0   ρf 4 0 µ + ρ0 K 3   µ 0  0 −K + αρ0 ρs 0 0 0 0 0 0 1   ρ0  0 0 0 1 ρs 0 0 0 0 0 0 0 0 0 0 0 − ρf0 K 0 ρ 0 0 0 − ρf0 K 0 ρ 0 0 0 0 0 0 0 0 0 0 0 αρ0 ρf  1 ρs   0  1   ρ0  . 0   0   0  0 (20) Заметим, что система уравнений является гиперболической, cобственные числа для Ax , Ay совпадают и образуют диагональную матрицу Λ = diag(0, 0, −cs , cs , −cp2 , cp2 , −cp1 , cp1 ). (21) Для решения данной системы используется метод расщепления по направлениям: на каждом шаге по времени 1. делаем шаг по X: решаем ~qt + Ax ~qx = ~0; 2. делаем шаг по Y : решаем ~qt + Ay ~qy = ~0 с начальными значениями после выполненного шага по X; 3. если задача не однородная, решаем ~qt = f~ с начальными значениями после выполненного шага по Y . Матрицы преобразования к инвариантам Римана и обратные к ним, одинаковые во всех точках сетки, для шага по X и шага по Y были найдены аналитически, что экономит потребляемую программой память и не требует дополнительного времени вычислений на каждом шаге. 3.3 Двумерный случай на криволинейных сетках Рассмотрим систему в виде (18). Опишем метод решения такой задачи на криволинейных структурных сетках. Для решения этой системы в двумерном случае сеточнохарактеристическим методом необходимо провести расщепление по направлениям. 11

Для прямоугольных сеток это были направления осей Ox, Oy; но для криволинейных сеток требуется расщепление вдоль произвольного направления Oξ. После поворота исходной системы координат Oxy на угол ϕ вокруг точки (x, y) (эта операция выполняется для каждого узла сетки), уравнения перейдут в ~q̃t + Aξ ~q̃ξ + Aη ~q̃η = f~, (22) ~q̃ = ~q̃(t, ξ, η) = ~q(t, x(ξ, η), y(ξ, η)), (23) Aξ = Ax cos ϕ + Ay sin ϕ, (24) Aη = −Ax sin ϕ + Ay cos ϕ. (25) где Рассмотрим однородную систему вдоль направления Oξ: ~q̃t + Aξ ~q̃ξ = 0. (26) Направление Oξ будем задавать единичным направляющим вектором ~n0 = (cos ϕ, sin ϕ)T . Перпендикулярный ему вектор ~n1 определяется формулой ~n1 = (− sin ϕ, cos ϕ)T . Определим симметричные тензоры 2 × 2 1 Nij = (~ni ⊗ ~nj + ~nj ⊗ ~ni ), i, j ∈ {0, 1}. 2 (27) Гиперболичность уравнений означает, что матрица Aξ (вдоль любого направления Oξ, не только матрицы Ax , Ay ) имеет полный набор собственных векторов. Таким образом, решение однородной задачи вдоль произвольного направления полностью аналогично решению задачи в одномерном случае. Матрица Λ одинакова для любого направления и имеет вид Λ = diag(0, 0, −cs , cs , −cp2 , cp2 , −cp1 , cp1 ). (28) Важным на практике отличием от случая прямоугольных сеток является то, что матрица Aξ зависит от угла ϕ, определяемого в задаче направлением осей криволинейной сетки и, следовательно, отличного в каждом узле сетки. С одной стороны, это делает выражения для матриц Ω, Ω−1 более громоздкими. С другой стороны, эти матрицы также отличны в каждом узле сетки, однако их хранение для каждого узла приводит к большим затратам памяти, а повторное вычисление на каждом шаге — к 12

неприемлемо большому времени расчёта. Эту проблему возможно решить, записав преобразования Ω~q и Ω−1 ω ~ в тензорной записи, разделяя части, связанные с параметрами модели и с углом ϕ. Ещё одним улучшением алгоритма является вычисление не самих новых значений инвариантов Римана ω ~ , а только прибавки ω ~ (n+1) − ω ~ (n) . Это позволяет вычислять только те инварианты, для которых собственные числа отличны от нуля. Таким образом, для рассмотренной модели необходимы только 6 инвариантов Римана. Их вычисление из исходных переменных выполняется по следующим формулам: √ ω1,2 = ∓ µρs (~n1 , ~u) + N01 : h, 1 N00 : h + d4 · p, ρs 1 = ±cp1 · (~n0 , ~u) ± d1 · (~n0 , ~v ) − N00 : h + d5 · p, ρs (29) ω3,4 = ∓cp2 · (~n0 , ~u) ∓ d0 · (~n0 , ~v ) + (30) ω5,6 (31) где в каждом выражении используются или верхние знаки, или нижние. Введённые коэффициенты di являются комбинациями параметров среды, определяемые следующим образом: h r1 = − 6Kρ0 ρf 6αρ2s − 4µ + 3αρ2f ρs + 3αρ3s + 4µρs + ρ20 9α2 ρ2f ρ2s + 6αρf ρs · i1/2 2 2 2 2 2 3αρs − 4µ + 4µ + 3αρs + 9K ρ0 , (32) r1 )+αρ20 d0 = ρf αρ0 (3Kρf + (3 (K − αρ0 ρs ) − 3K + 4µ)+3K (K − 3 (K − αρ0 ρs )) / q √ 6 (K − αρ0 ρs ) / ρ0 ρs (3Kρf + ρ0 (−3 (K − αρ0 ρs ) + 3K + 4µ) − 3Kρs − r1 ), (33) 3Kρf )+αρ20 d1 = ρf −αρ0 (r1 − (3 (K − αρ0 ρs ) − 3K + 4µ)+3K (K − 3 (K − αρ0 ρs )) √ q 6 (K − αρ0 ρs ) / ρ0 ρs (3Kρf + ρ0 (−3 (K − αρ0 ρs ) + 3K + 4µ) − 3Kρs + r1 ); (34) 3K 2 − αρ20 (3K + µ) d2 = , 6Kµ αρ20 , 2K ρ0 (3 (K − αρ0 ρs ) + 4µ) + r1 d4 = , 6ρ0 ρs (K − αρ0 ρs ) d3 = d5 = r1 − ρ0 (3 (K − αρ0 ρs ) + 4µ) 6ρ0 ρs (K − αρ0 ρs ) 13 (35) (36) (37) (38) /

Обратное преобразование от добавочной части инвариантов Римана к исходным неизвестным производится по следующим формулам: 1 1 ~u += √ (ω2 − ω1 )~n1 + [d1 (ω3 − ω4 ) + d0 (ω5 − ω6 )] ~n0 , 2 µρs g0 ~v += 1 [cp1 (ω4 − ω3 ) + cp2 (ω6 − ω5 )]~n0 , g0 p += g1 · (ω3 + ω4 + ω5 + ω6 ), h += (ω1 + ω2 ) · N01 + g3 [(ω3 + ω4 )(g4 N11 + g6 I) + (ω5 + ω6 )(g5 N11 − g7 I)], (39) (40) (41) (42) где gi — комбинации параметров среды, определяемые следующим образом: g0 = 2(cp1 d0 − cp2 d1 ), g1 = 1 , 2d4 + 2d5 g2 = d2 − d3 , g3 = 1 , (d2 − d3 )(2d4 + 2d5 ) (43) (44) (45) (46) g4 = 2d3 d5 ρs + 1, (47) g5 = 1 − 2d3 d4 ρs . (48) Общая схема решения двумерного уравнения на каждом шаге следующая: 1. шаг по Oξ: решаем однородное уравнение вдоль направления Oξ, заданного углом ϕ; 2. шаг по Oη: решаем однородное уравнение вдоль направления Oη, заданного углом ϕ + π/2 с начальными значениями после шага 1; 3. решить уравнение ~qt = f~ (если исходное уравнение неоднородно) с начальными значениями после шага 2. Таким образом, имея заданные начальные условия, получаем прямой алгоритм с итеративными шагами по времени. 14

3.4 Граница Приведённые выше формулы верны, строго говоря, только во внутренних точках расчётной области. Согласно формулам (12), (21) в матрице Λ (для случая любой размерности) есть и положительные, и отрицательные собственные числа. Следовательно, для нахождения значений инвариантов на характеристиках необходимо интерполировать значения функции и «слева», и «справа» от рассматриваемой точки сетки по формуле (16). Таким образом, в граничных точках некоторые характеристики выходят за границу расчётной области, и для корректного определения значений соответствующих инвариантов Римана требуются дополнительные граничные условия. В силу линейности модели рассматриваются линейные граничные условия, которые могут быть записаны следующим образом: ~ D·ω ~ = d, (49) где D — прямоугольная матрица, ω ~ — вектор инвариантов Римана; само равенство должно выполняться в граничных точках. При этом условия (49) должны задавать все те и только те инварианты, для которых характеристики выходят из расчётной области. Обычно граничные условия формулируются не в терминах инвариантов Римана, а в терминах исходных физических переменных. Тогда граничные условия записываются в таком виде: B · ~q = ~b, (50) где B — прямоугольная матрица с линейно независимыми строками, число которых равно числу выходящих за границу характеристик. Обозначим с помощью ~qin вектор неизвестных на (n + 1)-ом шаге по времени, посчитанный в граничной точке только по внутренним характеристикам, а Ωout — прямоугольная матрица, в столбцах которой находятся собственные векторы, соответствующие выходящим характеристикам. Тогда можно вывести формулу граничной корректировки, позволяющую правильно посчитать значение на границе на (n + 1)-ом шаге по времени: ~q(n+1) = ~qin + Ωout (B · Ωout )−1 · (~b − B~qin ). 15 (51)

Далее рассмотрим несколько примеров различных использующихся граничных условий и формул корректировки. Формулы гран. условий приведены в векторной записи, что позволяет их использовать для случая любой размерности. Отметим, что для описанного сеточно-характеристического метода на структурных сетках граничная корректировка всегда происходит во время шага по координате, направление которой перпендикулярно самой границе. Обозначим ~n — единичный вектор внешней нормали к границе, ~τ — единичный вектор касательной. Для многих задач необходима постановка условий так называемой «свободной границы». В работе были рассмотрены случаи открытых и закрытых пор. В случае открытых пор должны выполняться следующие условия на границе: отсутствие давления жидкости в порах, отсутствие нормального напряжения, отсутствие касательного напряжения. Эти условия могут быть записаны следующим образом:   p = 0, (52)  ||h|| · ~n = ~0. Если поры закрыты, т.е. насыщающая жидкость не взаимодействует непосредственно с внешней средой, то по-прежнему должно отсутствовать касательное напряжение, но давление не обязано занулиться. Вместо этого полное нормальное напряжение будет равно нулю. В силу закрытости пор, нормальные компоненты скоростей скелета и флюида должны совпадать. Таким образом, получаем следующие условия:     −||h|| · ~n − p = 0,    (53) ||h|| · ~τ = 0,      (~u, ~n) − (~v , ~n) = 0. Ещё одно используемой граничное условие: заданная на границе скорость скелета V~ совместно с условием равенства нормальных к границе скоростей скелета и флюида (непроницаемость). Эти условия задаются следующим образом:   ~u = V~ , (54)  (~u − ~v , ~n) = 0. Также используется ещё одно граничное условие, задающее на границе нормальное напряжение скелета с обратным знаком H и давление флюида P , а касательное 16

напряжение скелета равно нулю:     (h · ~n, ~n) = H,    (h · ~n, ~τ ) = 0,      p = P. (55) Для случая прямоугольных сеток выражения для граничной корректировки могут быть достаточно легко найдены по формуле (51), для этого необходимо записать условия в матричном виде (50), что несложно сделать, поэтому в данной работе мы их не приводим. Однако для случая криволинейных сеток угол ϕ, определяющий направление ~n, вообще говоря, отличен в каждой точке, и многие компоненты матрицы B не зануляются, в связи с чем выражения получаются значительно более громоздкими, а сами формулы корректировки различными в различных точках сетки. В двумерном случае они были записаны аналитически с помощью тензорных обозначений, что позволяет при численном моделировании избежать как хранения матриц для каждого узла, так и их вычисления на каждом шаге расчёта. Для свободной границы в случае открытых пор (52) должна выполняться следующая корректировка: ~u = ~uin + c1 · N00 : hin − c2 pin · ~n + c0 · N00 : hin · ~n1 , (56) ~v = ~v in + c3 · N00 : hin + c4 pin · ~n, h = (I + c5 · N00 ) : hin − c6 pin · N11 , (57) p = 0. (59) (58) Использованные коэффициенты ci определяются следующим образом: √ c0 = 1/ µρs , c1 = d0 − d1 , (d4 + d5 )g0 g1 ρs d1 d4 + d0 d5 , (d4 + d5 )g0 g7 cp1 − cp2 c3 = , (d4 + d5 )g0 g1 ρs c2 = c4 = cp1 d4 + cp2 d5 , (d4 + d5 )g0 g7 17 (60) (61) (62) (63) (64)

c5 = g5 − g4 − 2, (d4 + d5 )g2 ρs (65) d4 g4 + d5 g5 g1 . (d4 + d5 )g0 g7 (66) c6 = Для граничных условий (54) формулы корректировки следующие: ~u = V~ , (67) ~v = ~n1 · (~n1 , ~v in ) + ~n0 · (~n0 , V~ ), (68) p = pin + (~n0 , e5 (~uin − V~ ) + e6 (~v in − V~ )), (69) h = hin +(~n0 , ~uin −V~ )·(e2 I+e3 N11 )+e4 ·(~n1 , ~uin −V~ )N11 +(~n0 , ~v in −V~ )·(e0 I+e1 N11 ). (70) Использованные выше коэффициенты ei определяются следующими выражениями: e0 = 2g1 ρs (d1 d4 + d0 d5 ), (71) e1 = 2g1 g2 (d0 g4 − d1 g5 ), (72) e2 = 2g1 ρs (cp1 d4 + cp2 d5 ), (73) e3 = 2g1 g2 (cp2 ∗ g4 − cp1 ∗ g5 ), (74) √ e4 = 2 µρs , (75) e5 = 2g1 (cp2 − cp1 ), (76) e6 = 2g1 (d0 − d1 ). (77) Для граничных условий (55) формулы корректировки следующие: ~u = ~uin + c1 (N00 : hin − H) − c2 (pin − P ) · ~n + c0 N01 : hin · ~n1 , (78) ~v = ~v in + c1 (N00 : hin − H) − c2 (pin − P ) · ~n + c0 N01 : hin · ~n1 , (79) h = I : hin + c5 (N00 : hin − H) − 2H − c6 (pin − P ) · N11 + S · I, (80) p = P. (81) 18

3.5 Контакт Для многих практических задач необходимо моделировать в одном расчёте слои с различными реологическими характеристиками. Представляется целесообразным делать это путём создания нескольких расчётных сеток, каждая из которых описывает соответствующий слой своими уравнениями. Такой подход приводит к необходимости решать задачу контакта нескольких сред, описываемых различными физическими моделями. В данной работе были построены контактные условия между моделью Доровского и акустической моделью, моделью Доровского и линейно упругой моделью. Рассмотрим построение контактных условий. Пусть контактируют две сетки, обозначаемые нижними индексами a и b. В работе система уравнений в каждой из сеток решалась сеточно-характеристическим методом. Для каждой из сеток в граничных точках часть характеристик выходят из рассматриваемой сетки, следовательно, для их задания нужны граничные условия. Мы можем выбрать некоторым образом (каким именно, уточним далее) условия в форме (50), по которым строится корректор (51). Линейные контактные условия в общем виде записываются следующим образом:   ~qa C   = ~c, ~qb (82) где C — прямоугольная матрица размерами R × (na + nb ), na — размер вектора ~qa , nb — размер вектора ~qb , R — суммарное число выходящих характеристик для обеих сеток. Корректировка на границе (51) может быть записана в следующем виде: (n+1) ~ui ~ = Φi~uin i + fi , i ∈ {a, b}. (83) Матрицы Φi известны, векторы f~i неизвестны. Подставляя (83) в (82), получим систему линейных алгебраических уравений относительно f~i . Выбор граничных условий для нужных (выходящих) характеристик гарантирует однозначную разрешимость этой системы. Таким образом, получаем следующий алгоритм обработки контактных узлов на каждом шаге по времени (рассмотрен шаг по координате, направляющий вектор которой перпендикулярен контактной границе): 19

1. выполнить все вычисления для шага (для граничных узлов могут использоваться дополнительные ghost-узлы); 2. решить СЛАУ (82) относительно [f~a f~b ] в каждом граничном узле; 3. выполнить граничную корректировку (83) с fi с шага 2 для обеих сред. Отметим, что для сред a и b граничная корректировка должна выполняться с противоположными нормалями, поскольку формулы обоих граничных корректоров выведены в предположении внешней нормали. Ниже приведены условия (4 уравнения в двумерном случае, 5 в трёхмерном), которые должны выполняться в каждой точке на контакте между пористой средой, описываемой моделью Доровского, и средой, описываемой уравнениями акустики:    (−||h|| · ~n, ~n) = −Pa ,       (−||h|| · ~n, ~τ ) = 0, (84)    (1 − β)(~u, ~n) + β(~v , ~n) = (V~a , ~n),       ρf · p = P . a βρ0 Здесь Pa — давление в акустической среде, V~a — скорость частиц в акустической среде. Аналогичные условия на контакте с линейно упругой средой (5 уравнений в двумерном случае, 7 в трёхмерном):    ~u = ~velastic ,       (−h · ~n0 , ~n0 ) − p = (T · ~n0 , ~n0 ), (85)    (−h · ~n0 , ~n1 ) = (T · ~n0 , ~n1 ),      (~u, ~n ) − (~v , ~n ) = 0. 0 0 Здесь ~velastic — скорость частиц в упругой среде, T — симметричный тензор напряжений в упругой среде. 20

4 Результаты расчётов 4.1 Задача нагружения пористой среды Для демонстрации свойств уравнений Доровского и выбранного метода решения на прямоугольных сетках, а также сравнения с другими использующимися физическими моделями были проведены расчёты для задачи динамического нагружения среды. При этом среда моделировалась различными моделями: акустической, линейно упругой, пористостой в рамках модели Доровского. Параметры расчёта были следующими: • Область — прямоугольник 10 км х 3 км; • Равномерная прямоугольная сетка с шагом 5 м; • Шаг по времени dt = 2 мс, было выполнено 500 шагов; • Использовался сеточно-характеристический метод второго порядка; • Источник — начальные условия в сферической области (глубина 200 м, радиус 50 м) — шаровой тензор напряжений 104 Па; • Параметры акустической среды: cp = 2000 м/с, ρ = 1500 кг/м3 ; • Параметры упругой среды: cp = 2000 м/с, cs = 1300 м/с, ρ = 1500 кг/м3 ; • Параметры пористой среды: cp1 = 2000 м/с, cp2 = 450 м/с, cs = 1300 м/с, ρfs = 1500 кг/м3 , ρff = 1000 кг/м3 , пористость β = 10% (через них пересчитываются K, µ, α). На представленных рисунках 1, 2 приведены волновые поля (отображается модуль скорости) для части расчётной области через 1 с после начала расчёта. Во всех моделях хорошо видно отражение от свободной границы. Видны различные типы волн: • в акустической модели — только продольные волны (исходная и отражённая от свободной границы); 21

(a) Акустическая среда (b) Линейно упругая среда Рис. 1: Волновое поле через 1 с после начала расчёта • в линейно упругой модели — продольные и поперечные волны (как исходные, так и отражённые от свободной границы), а также поверхностные (волна Рэлея); • в модели Доровского — два типа продольных волн, поперечные волны, поверхностные волны (включая исходные и отражённые). 4.2 Задача морской сейсмической разведки Была рассмотрена упрощённая задача морской сейсморазведки: распространение сейсмической волны от пневмопушки в морской воде и её отражение от подстилающего дна — донных осадков. Рассматривался двумерный случай. Задача решалась на двух прямоугольных сетках: верхняя описывала слой воды уравнениями акустики, нижняя — слой водонасыщенных осадков с помощью модели Доровского. На контакте сеток использовались условия (84). Параметры модели следующие: • верхняя сетка: прямоугольник 1000 м х 250 м; • нижняя сетка: прямоугольник 1000 м х 50 м; • параметры акустической среды: c = 1500 м/с, ρ = 1000 кг/м3 ; 22

Рис. 2: Волновое поле для модели Доровского через 1 с после начала расчёта • параметры насыщенной среды: cp1 = 2000 м/с, cp2 = 450 м/с, cs = 1300 м/с, ρfs = 1500 кг/м3 , ρff = 1000 кг/м3 , пористость β = 10% (через них пересчитываются K, µ, α); • шаг сетки h = 1 м; • шаг по времени dt = 0.4 мс из условия Куранта (17); Действие пневмопушки моделировалось точечным источником, заглублённым на 3 м в акустической среде, со временной зависимостью в виде импульса Рикера (с основной частотой 30 Гц). Для решения использовался сеточно-характеристический метод с интерполяцией схемой Русанова с монотонизатором minimax третьего порядка. На рисунке 3 представлено распределение давления в воде, моделируемой уравнениями акустики, и во флюиде, насыщающем поры донных осадков, моделируемых 23

моделью Доровского, в последовательные моменты времени. Видно, что вначале от точечного источника распространяется сферическая волна. Она взаимодействует с донными осадками, за счёт чего появляется отражённая волна, распространяющаяся обратно к дневной поверхности. При этом видно, что в насыщенной среде инициируются две продольные волны с различными скоростями распространения согласно свойствам модели. (a) 160 мс после начала расчёта (b) 220 мс после начала расчёта Рис. 3: Распределение давления 24

4.3 Задача прискважинной диагностики Была рассмотрена прямая задача исследования прискважинной зоны с помощью акустического зондирования в двумерной постановке. Схема задачи представлена на рисунке 4. Были использованы три расчётные сетки, каждая из которых моделировала среду соответствующими уравнениями. Внутренняя подобласть — круг радиусом 0.2 м. Предполагается, что скважина заполнена водой, потому эта подобласть моделировалась уравнениями акустики с параметрами ρ = 1000 кг/м3 , c = 1500 м/с. Вокруг этой подобласти находится кольцо прискважинной зоны шириной 0.3 м. В силу насыщения жидкостью, поступающей из скважины, для моделирования этого слоя использовалась модель Доровского со следующими параметрами (в большей части кольца): cp1 = 2000 м/с, cp2 = 400 м/с, cs = 1300 м/с, ρs = 1800 кг/м3 , ρf = 100 кг/м3 , пористость β = 10 %. Третий слой представляет собой внешнее кольцо из породы, моделируемой уравнениями линейной упругости с параметрами ρ = 2000 кг/м3 , cp = 2000 м/с, cs = 1300 м/с. Был рассмотрен случай наличия неод- Рис. 4: Схема задачи нородности — наличие трещины в пористой среде. Для нефтедобывающих скважин подобные трещины могут служить проводящими путями, по которым углеводороды поступают в скважину из окружающей породы. Были рассмотрены трещины, ориентированные вдоль радиуса, с отношением раскрытости к протяжённости 1:10. В 25

одном случае трещина направлена под углом 90◦ (на рисунке 4), в другом — 180◦ . В расчёте трещины моделировались заданием подобласти внутри второй сетки (внутреннее кольцо) в форме эллипса с другими значениями параметров: µ = 1.1e9 Н/м2 , α = 2017.983999 м5 ·кг·с2 , K = 1.524554690e9 Н/м2 , ρs = 1100 кг/м3 , ρf = 450 кг/м3 , пористость β = 45 %. Между скважиной и прискважинной средой были поставлены контактные условия (84), между насыщенной прискважинной средой и окружающей её извне породой — контактные условия (85). На внешней границе использовались неотражающие граничные условия. Задача решалась на структурных криволинейных сетках с максимальным шагом ячейки 5 мм, шаг по времени dt = 0.435 мкс был выбран из условия Куранта. Для решения отдельных уравнений переноса использовалась схема Русанова третьего порядка точности. Для моделирования импульса использовался точечный источник в центре акустической среды со временной зависимостью в форме импульса Пузырёва с несущей частотой 40 кГц. Сохранялись значения давления на границе скважина–насыщенная среда. На рисунке 5 приведена нормированная разность давлений по периметру скважины для случаев наличия одной из трещин и случая однородной насыщенной среды в момент времени, соответствующий приходу волны от источника. Видны изменения амплитуды пришедшей волны в случае наличия трещины на азимуте, соответствующем этой трещине. Также проведённые эксперименты показали, что при наличии неоднородности после взаимодействия приходящей волны с ней инициируется поверхностная волна на контакте между скважиной и насыщенной окружающей средой. Эта волна распространяется со скоростью около 1620 м/с. На рисунке 6a представлен их вид. Сплошная линия обозначает случай ориентации трещины на 90◦ , штрихпунктирная — на 180◦ . Из-за отличия в характеристиках насыщенной среды и окружающей её извне породы при падении прямой волны на границу между ними возникает отражённая волна, движущаяся обратно к центру скважины. Поскольку подобласть трещины моделировалась параметрами, отличными от параметров в остальной среде, то волна, 26

Рис. 5: Нормированное отклонение давления через 152 мкс после начала расчёта пришедшая с азимута трещины, приходит позже (задержка на 30 мкс из-за разных значений cp1 ) и с изменённой амплитудой. Распределение нормированного отклонения давления в момент прихода отражённой волны приведено на рисунке 6b (обозначения такие же, как и на рисунке 6a). 5 Заключение В работе был исследован процесс распространения сейсмических волн в пористых флюидонасыщенных средах. Для описания их динамического поведения использовалась модель Доровского. Система определяющих уравнений была аналитически исследована, были найдены собственные числа и собственные векторы. Для её численного решения был адаптирован сеточно-характеристический метод на прямоугольных сетках в двумерном и трёхмерном случае, а также на структурных криволинейных сетках в двумерном случае. Были построены корректные граничные и контактные условия, задающиеся в явном виде. Поддержка вышеописанных методов была добавлена в программный комплекс RECT лаборатории прикладной вычислительной геофизики, написанный на С++, с поддержкой параллельных вычислений на 27

(a) 200 мкс после начала расчёта (b) 462 мкс после начала расчёта Рис. 6: Распределение нормированного давления по периметру скважины основе технологий MPI и OpenMP. Было проведено численное моделирование задачи динамического нагружения пористой среды, результаты которого продемонстрировали наличие различных типов упругих волн, ожидаемых в данной модели. Было показано, что использование 28

модели Доровского качественно учитывает все волны, распространяющиеся в линейно упругой среде (с отличающимися амплитудами), а также предусматривает наличие нескольких дополнительных волн, обусловленных двухфазностью среды. Таким образом, данная модель пористой среды является в некотором смысле обобщением теории линейной упругости. Были корректо поставлены и явным образом заданы граничные условия «дневной поверхности». Был проведён расчёт, демонстрирующий возможность использования модели для решения задач морской сейсморазведки. В этом случае учёт пористости среды может иметь большое значение, поскольку слой донных осадков может быть достаточно тонким, а границы с другими слоями — достаточно контрастными, чтобы медленная продольная волна, описываемая моделью Доровского, могла быть замечена на записываемых сейсмограммах. Также было проведено моделирование задачи в двумерной радиальной постановке об акустической диагностике прискважинной среды. С точки зрения вычислительной математики, данная задача демонстрирует успешную реализацию численного метода для решения задачи сеточно-характеристическим методом на существенно криволинейных структурных сетках, а также реализацию контактных условий между акустической и пористой средами, линейно упругой и пористой средами. С точки зрения практической значимости, пласт породы вокруг скважины вследствие процедуры бурения во многих задачах может быть моделирован корректно только моделями насыщенных сред (например, моделью Био или использованной в данной работе моделью Доровского), а само акустическое зондирование важно для исследования трещиноватости окружающей среды, что требуется в некоторых технологиях нефтедобычи. Результаты работы были представлены на конференции Геомодель-2019 [15], конференции МФТИ-62, а также опубликованы в двух статьях [16] и [17], входящих в перечень ВАК, Scopus и WoS. В дальнейшем планируется применить разработанный подход для расчёта распространения сейсмических волн в реалистичной модели нефтяного месторождения. Использование уравнений Доровского для описания продуктивного пласта обосновано наличием в нём флюида (нефти). Постановка контактных условий в явном 29

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

Список литературы [1] Carcione J.M. Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media. — 3rd edition. — Elsevier, 2014. — 690 p. [2] Квасов И.Е., Левянт В.Б., Петров И.Б. Решение прямых задач сейсморазведки в трещиноватых средах методом сеточно-характеристического моделирования (с целью исследования возможности прямого обнаружения трещиноватых зон). — М.: ООО «ЕАГЕ Геомодель», 2018. — 332 с. [3] Romenski E., Reshetova G., Peshkov I., Dumbser M. Modeling wavefields in saturated elastic porous media based on thermodynamically compatible system theory for multiphase mixtures. // arXiv.org: physics. — 2019. [Электронный ресурс]. arXiv:1910.04207. URL: https://arxiv.org/abs/1910.04207 (дата обращения: 12.06.2020). [4] Frenkel J. On the Theory of Seismic and Seismoelectric Phenomena in a Moist Soil // Journal of Physics. 1944. — №3(5). — P. 230–241. [5] Biot M.A. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range. // The journal of the acoustical society of America. 1956. — №28(2). — P. 168–178. [6] Burridge R., Keller J.B. Poroelasticity equations derived from microstructure // The Journal of the Acoustical Society of America. 1981. — №70(4). P. 1140–1146. [7] Pride S.R., Gangi A.F., Morgan F.D. Deriving the equations of motion for porous isotropic media // The Journal of the Acoustical Society of America. — 1992. — №92. — P. 32–78. [8] Chen H., Gilbert R.P., Guyenne P. A Biot model for the determination of material parameters of cancellous bone from acoustic measurements // Inverse Problems. — 2018. — №34(8). [9] Sidler R. A porosity-based Biot model for acoustic waves in snow // Journal of Glaciology. — 2015. — №61(228). — P. 789–798. 31

[10] Sidler R., Carcione J.M., Holliger K. Simulation of surface waves in porous media // Geophysical Journal International. — 2010. — №183(2). — С. 820–832. [11] Blokhin A.M., Dorovsky V.N. Mathematical Modeling in the Theory of Multivelocity Continuum. — New York: Nova Science Publishers Inc., 1995. — 182 p. [12] Dorovsky V.N., Perepechko Yu.V., Fedorov A.I. Stoneley waves in the Biot-Johnson and continuum filtration theories // Geologiya i geofizika. — 2012. — №53(5). — P. 621–632. [13] Imomnazarov Kh.Kh., Kholmurodov A.E. Direct and inverse dynamic problems for SH-waves in porous media // Mathematical and Computer Modelling; 2007. — №45(3-4): 270–280. [14] Sinev A.V., Romensky E.I., Dorovsky V.N. Effects of a mudcake on Stoneley waves in a fluid-filled porous formation around a borehole // Russian Geology and Geophysics. — 2012. — №53(8). — P. 823–828. [15] Голубев В.И., Шевченко А.В., Петров И.Б. Моделирование распространения сейсмических волн в флюидонасыщенных геологических породах сеточнохарактеристическим методом. Геомодель-2019 (21-я конференция по вопросам геологоразведки и разработки месторождений нефти и газа “Геомодель-2019”, 9-13.09.2019), Геленджик. [16] Голубев В.И., Шевченко А.В., Петров И.Б. Об учёте водонасыщенности донных осадков в задаче морской сейсмической разведки. // Доклады Академии наук. — 2019. — №488(3). — C. 248–252. [17] Петров И.Б., Голубев В.И., Шевченко А.В. О задаче акустической диагностики прискважинной зоны. // Доклады Академии наук. — 2020. — №492. — C. 92–96. [18] Dorovsky V.N., Dorovsky S.V. An electromagnetoacoustic method of measuring electric conductivity and ζ-potential. // Russian Geology and Geophysics. — 2009. — №50(6). — P. 572–578. 32 Powered by TCPDF (www.tcpdf.org)