CC BY
8
УДК 550.834
DOI 10.33764/2618-981X-2022-2-2-246-253
Численное решение уравнений Био в квазистатической постановке для оценки сейсмического затухания в анизотропных средах
С. А. Соловьев1 *, М. А. Новиков2, В. В. Лисица1
1 Институт математики им. С. Л. Соболева СО РАН, г. Новосибирск, Российская Федерация
2 Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, г. Новосибирск,
Российская Федерация
* e-mail: novikovma@ipgg.sbras.ru
Аннотация. Представлен алгоритм численного апскейлинга упругих свойств трещинно-пори-стых сред в низкочастотном диапазоне. Для этого разработан алгоритм решения уравнений Био в квазистатическом состоянии. В статье описываются следующие аспекты: формулируется задача апскейлинга пороупругих сред для получения частотно-зависимого тензора жесткости, соответствующего вязкоупругим средам; конечно-разностная аппроксимация рассматриваемой краевой задачи; особенности численного решения СЛАУ с численным экспериментом и анализом производительности. Проведены численные эксперементы для анализа влияния связности трещин и свойств их наполнителя на частотно-зависимое затухание в трещиновато-пористой флюидонасыщенной среде.
Ключевые слова: пороупругость, модель Био, конечно-разностная аппроксимация, сейсмическое затухание, прямые методы решения СЛАУ
Numerical solution of Biot equations in quasistatic state for seismic attenuation estimation in anisotropic media
S. A. Solovyev1 *, M. A. Novikov2, V. V. Lisitsa1
1 Institute of Mathematics SB RAS, Novosibirsk, Russian Federation 1Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk,
Russian Federation
* e-mail: novikovma@ipgg.sbras.ru
Abstract. An algorithm for numerical upscaling elastic properties of fractured-porous media in low-frequency range is presented. To do so, the algorithm to solve Biot equations in quasi-static state is developed. Paper describes the following aspects: stating the problem of the poroelastic media upscaling to obtain frequency-dependent stiffness tensor, corresponding to the visco-elastic media; the finite-difference approximation of the considered boundary-value problem; peculiarities of the numerical solution of the SLAE with numerical experiments and performance analysis. Numerical experiments are performed to observe the influence of fracture connectivity and fracture-filling material properties on frequency-dependent seismic attenuation in fractured porous fluid-saturated media.
Keywords: poroelasticity, Biot model, finite-difference approximation, seismic attenuation, direct methods for SLAE
Введение
Современное повышение актуальности развития технологий захоронения CO2 [1, 2], разведки геотермальной энергии [3, 4] ставит сложные задачи перед методами сейсмического мониторинга. Перспективным сейсмическим атрибу-
том для оценки транспортных свойств трещиноватых коллекторов является частотно-зависимое затухание. Известно, что изменение свойств флюида или геометрии порового пространства существенно влияют на индуцированные волной потоки флюидов (wave-induced fluid flow, WIFF), возникающие при распространении сейсмической волны в трещиновато-пористых средах [5,6]. Обычно эти потоки разделяют на два типа - между вмещающей породой и флюидом (fracture-to-background WIFF, FB-WIFF) и между трещинами (fracture-to fracture WIFF, FF-WIFF), причем у обоих типов свои характерные частоты и факторы, влияющие на их интенсивность [6,7,8] В теоретических и большинстве численных исследований WIFF рассматриваются простые модели среды с попарным пересечением трещин [6,7], кроме работы [9] с применением более реалистичных дискретных систем трещин. Одна из причин рассмотрения глобальной связности трещин -отсутствие эффективного численного алгоритма для моделирования экспериментов по нагружению образца. Требование репрезентативности размера модели [10,11] и достаточной детальности для отражения стрыктуры неоднородности. В случае трещинноватых сред репрезентативный размер связан с характерной длиной перколяции [12,13], а дискретизация зависит от характерного раскрытия трещины и, как правило, превышает 1000 узлов в направлении каждой координатной оси. Поскольку восстановление эффективных упругих свойств среды требует серии экспериментов для каждой реализации модели среды, то необходимо разработать эффективный алгоритм численного апскейлинга, применимый для решения ряда крупномасштабных задач.
В данной работе представлен алгоритм численного апскейлинга упругих свойств трещинноватых-пористых флюидонасыщенных сред. В работе представлена постановка задачи апскейлинга трещиновато-пористой флюидонасыщен-ной среды методом восстановления частотно-зависимого тензора жесткости, соответствующего эквивалентной вязкоупругой однородной среде. Для решения задачи разработан алгоритм численного решения системы уравнений Био в квазистатической постановке на основе их конечно-разностной аппроксимации. Алгоритм применен на примере трещиноватых сред разной длины перколяции для оценки сейсмического затухания в них и анализа влияния связности трещин (длины перколяции) на затухание.
Постановка задачи апскейлинга
Задача апскейлинга заключается в нахождении эффективных упругих свойств однородной среды, в каком-то смысле эквивалентной рассматриваемой неоднородной среде. В частности, в работе рассматривается неоднородная поро-упругая флюидонасыщенная среда, а под задачей апскейлинга понимается восстановление эффективных свойств вязкоупругой среды, эквивалентной рассматриваемой среде.
Анизотропная пороупругая флюидонасыщенная среда описывается уране-ниями Био в квазистатической постановке в плоскости Oxz [14, 15]:
д_
дх д_
дх
С дих
'11
С55 (
дх дих
С
13
д_Щ дг
+ а
м (
ди.
дг д_ дх д_ дг
дх
)
+ ■
д
дг
С
дх дих
д,
) С55(
дг
дих ди„
13
дх
С
33
д_Щ
дг
■ + а>
М (
дг
дх дм?„
дх дг
) )
= 0, = 0,
М (а1 М (а1
дих
ди
—— + а 2 —
дх дг
г ) + М ("''х
дих
дх
а
диг
дг
■)+м (
дх дг дх дг
■ Л
= т—w
к
х
• Л
№-W
г •
Здесь их, иг - компоненты вектора смещения твердых частиц упругого скелета, wx, wz - компоненты вектора относительного смещения флюида, Су -компоненты тензора жесткости флюидонасыщенного скелета, Л - вязкость флюида, к-1, ^2 - проницаемость в направлениях осей Ох и Ог соответственно, М -модуль взаимодействия между фазами, а1 и а1 - анизотропные параметры Био-Уиллиса [15,16], Ф - частота.
Методом апскейлинга восстанавливается тензор жесткости
С л Г Си /V СС13 /V СС15
С(ю) = /V С31 /V С33 /V С35
/V V С51 /V С53 /V С55
а
хх
Vе7 х,
а
хг
а
хг у
хх
V8 хг
'хг
'хг у
С £
в прямоугольной области О = [ ,¿2] х [Ь\,¿2] однородной вязкоупругой
/V А
среды, для которого результирующие средние деформации £ и напряжения о в рассматриваемой области и в такой же области с неоднородным пороупругим материалом (определяемые через смещения, полученные решением уравнений
/V
(1) [15,16] совпадают. Сам тензор С(ю) должен удовлетворять уравнениям
д_ дх
д_
дх
С„(ю) ^ + Аз(ш)
дх
дг
д + —
дг
С55(Ф)
дух
дг дх
С55 (®)
дг дх
+ ■
д_
дг
ду ду
С13(Ю) + Сзз(ш) ду
дх дг
= 0, = 0.
Для единственности решения задачи нахождения С(ф) необходимо рассмотреть три краевые задачи для уравнений Био. Предлагается решение трех задач базовых нагружений образца, определяемых граничными условиями:
Таблица 1
Граничные условия трех краевых задач по нагружению
Одноосное нагружение вдоль оси Ох Одноосное нагружение вдоль оси Ог Сдвиговое нагружение
ахх = 1, х е {¿х,Ц}, а^ = 0, г е {¿1,¿2}, ахг = 0, (х, г) е дБ а хх = 0, х е {Ц, а^ = 1, г е {¿1,¿2}, ахг = 0, (х, г) е дБ а хх = 0, х е {Ц, а^ = 1, г е{2?,¿2}, ахг = 1 (X, г) едБ
Непротекание: wx = 0, х е {,¿2}; = 0, г е , ¿2}
В результате решения трех краевых задач (уравнения (1) + краевые условия
/V
из Таблицы 1) девять искомых компонент С(ю) находятся решением СЛАУ, где
коэффициенты уравнений - осредненные компоненты тензоров деформаций £ и напряжений о.
Численное решение уравнений Био
Для численного решения уравнений (1) они аппроксимируются со вторым порядком на сдвинутых сетках с Ых х N г ячейками, где уравнения и пространственные производные аппроксимируются, а коэффициенты уравнений и компоненты смещений определяются в разных группах узлов (Рисунок 1) [16,17].
Пороупругость
п-п\до = О £
С(х,у),ах(х.у),а2(х.у), М(х,у),кх(х.у),кг(х.у)
О
Вязкоупругость
а ■ п\дв = ст0 д
Я
С{ ш)
о
Рис. 1. Схематичное изображение эквивалентных пороупругой и вязкоупругой областей (слева); расположение узлов аппроксимации уравнений и производных в узлах одной ячейки расчетной сетки (справа)
Вместе аппроксимацией граничных условий (Таблица 1) конечно-разностные уравнения определяют СЛАУ из 2 • х - 1)Nг +2 • Nх (Nг - 1) уравнений:
Ах = Ь, А е С *х*, Ь е С * (2)
Для численного решения системы (2) используется решатель MKL PARDISO на основе прямых методов решения СЛАУ.
Численные эксперименты
Верификация - слоистая среда. Первый ряд численных экспериментов -восстановление эффективного тензора жесткости для слоистой среды - проводится для верификации алгоритма. Для сравнения с аналитическими оценками численные оценки затухания и фазовой скорости (вдоль оси Ox ) вычисляются через компоненты тензора C( ю ) и средняя плотность материала в среде р
Vp =WQi/р, Q = 3Qi
Рассматривается слоистая среда со слоями двух материалов. Свойства материалов взяты из [19], сама модель и аналитические оценки - из работы [18]. Сравнение оценок представлено на Рис 2. Наблюдается достаточно хорошее соответствие оценок как для затухания, так и для скорости.
2600 2580 2560 2540 2520 2500 2480 2460
10° 102 10° 102 Частота, Гц Частота, Гц
Рис. 2. Сравнение оценок фазовой скорости (слева) и затухания (справа) в слоистой среде поперек слоев, полученных численно (красная сплошная линия) и аналитически (черная пунктирная линия) в зависимости от частоты
Влияние связности трещин на затухание. Второй набор численных экспериментов проводится для наблюдения влияния степени связности трещин на сейсмическое затухание вследствие "ЖР. Рассматриваются 6 наборов (по 8 реализаций) трещиноватых моделей с горизонтальными и вертикальными трещинами с разной длиной перколяции (примеры трех моделей представлены на Рис. 3). Детали генерации моделей представлены в [12]. Трещины заполнены анизотропным материалом, в зависимости от их ориентации (Таблица 2).
Полученные оценки затухания представлены на Рис. 3. Два пика затухания на относительно низкой частоте и высокой частоте соответствуют пикам затухания вследствие ЕВ-"ТРЕ и ЕЕ-"ТРЕ, соответственно. С увеличением связности
трещин флюидопотоки между трещин становятся интенсивнее, о чем свидетельствует рост высокочастотного пика. Также с ростом связности заметно падение пика затухания от флюидопотоков типа БВ-ОД^Т и его смещение по частоте, вероятно вызванные изменением характерного размера неоднородности. Результаты для разной связности трещин демонстрирует перспективность применения разработанного алгоритма для выявления зависимостей между структурой тре-щиноватости и затуханием.
Таблица 2
Физические свойства трещиноватой среды
Параметр Вмещ. порода Горизонт. тр-ны Вертик. т-ны Пересечения
С11, Па 6.910-1010 3.896 -1010 4.643 •Ю10 3.896 -1010
С^ Па 7.159 -1010 1.9811010 1.981 •Ю10 1.981 • 1010
Па 6.910 -1010 4.643 •Ю10 3.896 •Ю10 3.896 -1010
^ Па 3.097 -1010 1.131 -1010 1.1311010 1.1311010
М , Па 2.010-1010 9.429 •Ю9 9.429 •Ю9 9.402 •Ю9
а 1 м , Па 5.953 •Ю9 6.502 -109 6.051 •Ю9 6.502 -109
а 2м , Па 5.953-109 6.051 •Ю9 6.502 •Ю9 6.502 •Ю9
г| , Па • с к1 м2 11012 7.072 •Ю6 1.887-109 7.072 •Ю6
^ , Па • с к1 м2 11012 1.887-109 7.072 •Ю6 7.072 -106
х10"5
10° ю5 ю10
Частота, Гц
Рис. 3. Примеры трещиноватых моделей разной связности (слева). Осредненные по реализациям численные оценки частотно-зависимого затухания для трещиноватых моделей с ростом связности, соответствующем порядку цветов: черный, розовый, голубой, синий, зеленый, красный
Заключение
Разработан алгоритм численного апскейлинга, основанный на конечно-разностной аппроксимации системы уравнений Био анизотропной неоднородной
пористой флюидонасыенной среды в квазистатической постановке. Численным решением трех задач по нагружению образца пороупругой среды восстанавливается эффективный тензор жесткости эквивалентной вязкоупругой среды, с помощью которого вычисляются частотно-зависимые затухание и фазовая скорость. Результаты численных экспериментов по нагружению моделей трещиноватой среды и последующей оценки затухания показывают значительное влияние флюидопотоков обоих типов (FB-WIFF и FF-WIFF), а рост связности трещин приводит к существенному качественному изменению затухания, вкладов флюидопотоков разных типов на общее затухание. Однако, на высоких частотах сильное влияние на затухание оказывает рассеяние на системе трещин, поэтому в будущем для полноты оценки нужно проводить как моделирование нагруже-ния образцом, так и основанного на решении уравнений пороупругости в динамической постановке распространения сейсмической волны.
С. А. Соловьев разработал алгоритм численного апскейлинга, В. В. Лисица сформулировал постановки краевых задач нагружения образца при поддержке Российского научного фонда № 19-77-20004, М. А. Новиков провел анализ результатов численных экспериментов. Численное моделирование проводилось с использованием вычислительных ресурсов Сибирского суперкомпьютерного центра СО РАН (кластер NKS-30T).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Menke H.P., Reynolds C.A., Andrew M.G., Pereira Nunes J.P., Bijeljic B., Blunt M.J. 4d multi-scale imaging of reactive flow in carbonates: Assessing the impact of heterogeneity on dissolution regimes using streamlines at multiple length scales // Chemical Geology. - 2018. - Vol. 481. - P. 23-37.
2. Huang F., Bergmann P., Juhlin C., Ivandic M., Luth S., Ivanova A., Kempka, T., Henninges, J., Sopher, D., Zhang, F. The First post-injection seismic monitor survey at the ketzin pilot co2 storage site: results from time-lapse analysis // Geophysical Prospecting. - 2018. - Vol. 66. - № 1. - P. 6284.
3. Salaun N., Toubiana H., Mitschler J.B., Gigou G., Carriere X., Maurer V., Richard A. Highresolution 3d seismic imaging and refined velocity model building improve the image of a deep geothermal reservoir in the upper rhine graben // The Leading Edge. - 2020. - Vol. 39. - № 12. -P. 857-863.
4. Marty N.C.M., Hamm V., Castillo C., Thiery D., Kervevan C. Modelling water-rock interactions due to long-term cooled-brine reinjection in the dogger carbonate aquifer (paris basin) based on in-situ geothermal well data // Geothermics. - 2020. - Vol. 88. - P. 101899.
5. Muller T.M., Gurevich B., Lebedev M. Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks - a review // Geophysics. - 2010. - Vol. 75. - № 5. -P. 75A147-75A164.
6. Rubino, J.G., Muller, T.M., Guarracino, L., Milani, M., Holliger, K. Seismoacoustic signatures of fracture connectivity // Journal of Geophysical Research: Solid Earth. - 2014. -Vol. 119. - № 3. - P. 2252-2271.
7. Kong L., Gurevich B., Zhang Y., Wang Y. Effect of fracture fill on frequency-dependent anisotropy of fractured porous rocks // Geophysical Prospecting. - 2017. - Vol. 65. - № 6. - P. 16491661.
8. Guo, J., Rubino, J.G., Glubokovskikh, S., Gurevich, B. Effects of fracture intersections on seismic dispersion: theoretical predictions versus numerical simulations // Geophysical Prospecting. - 2017. - Vol. 65. - № 5. - P. 1264-1276.
9. Hunziker J., Favino M., Caspari E., Quintal B., Rubino J. G., Krause R., Holliger K. Seismic attenuation and stiffness modulus dispersion in porous rocks containing stochastic fracture networks // Journal of Geophysical Research: Solid Earth. - 2018. - Vol. 123. - № 1. - P. 125-143.
10. Ovaysi, S., Wheeler, M., Balhoff, M. Quantifying the representative size in porous media // Transport in Porous Media. - 2014. - Vol. 104. - № 2. - P. 349-362.
11. Bazaikin, Y., Gurevich, B., Iglauer, S., Khachkova, T., Kolyukhin, D., Lebedev, M., Lisitsa, V., Reshetova, G. Effect of ct image size and resolution on the accuracy of rock property estimates // Journal of Geophysical Research: Solid Earth. - 2017. - Vol. 122. - № 5. - P. 3635-3647.
12. Novikov M.A., Lisitsa V.V., Bazaikin Y.V. Wave propagation in fractured porous media with different percolation length of fracture systems // Lobachevskii Journal of Mathematics. - 2020. - Vol. 41. - № 8. - P. 1533-1544.
13. Xu, C., Dowd, P.A., Mardia, K.V., Fowell, R.J. A connectivity index for discrete fracture networks. Mathematical Geology. - 2006. - Vol. 38. - № 5. - P. 611-634
14. Biot, M.A. Theory of propagation of elastic waves in fluid-saturated porous solid. i. low-frequency range // Journal of the Acoustical Society of America. - 1956. - Vol. 28. - P. 168-178.
15. Cheng, A.H.D. Material coefficients of anisotropic poroelasticity // International Journal of Rock Mechanics and Mining Sciences. - 1997. - Vol. 34. - № 2. - P. 199-205.
16. Masson Y. J., Pride S. R. Finite-difference modeling of biot's poroelastic equations across all frequencies // Geophysics. - 2010. - Vol. 75. - № 2. - P. N33-N41.
17. Virieux J. P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method // Geophysics. - 1986. - Vol. 51. - P. 889-901.
18. Rubino, J.G., Caspari, E., M'uller, T.M., Milani, M., Barbosa, N.D., Holliger, K. Numerical upscaling in 2-d heterogeneous poroelastic rocks: Anisotropic attenuation and dispersion of seismic waves // Journal of Geophysical Research: Solid Earth. - 2016. - Vol. 121. - № 9. - P. 6698-6721.
19. Carcione, J., Picotti, S. P-wave seismic attenuation by slow-wave diffusion: Effects of inhomogeneous rock properties // Geophysics. - 2006. - Vol. 71. - № 3. - P. O1-O8.
© С. А. Соловьев, М. А. Новиков, В. В. Лисица, 2022
