CC BY
15
УДК 519.63
Д. Н. Ворощук1, В. А. Миряха1'2, И. Б. Петров1'2, А. В. Санников1'2
1 Московский физико-технический институт (государственный университет) 2Институт автоматизации проектирования РАН
Исследование сейсмического отклика от кластера субвертикальных макротрещин разрывным методом
Галеркина
Разрывный метод Галеркина на неструктурированных сетках адаптирован и реализован для моделирования волновых откликов от систем субвертикальных макротрещин в карбонатных породах для численного решения прямых задач сейсморазведки. В работе проводится сравнение сейсмических откликов для нескольких механико-математических моделей трещиноватого коллектора. Модели различаются способом задания коллектора: явное выделение макротрещин c параметрами среды в области коллектора, совпадающими с вмещающей средой, либо отличными от неё. Показана возможность учёта межтрещинных волновых взаимодействий с помощью используемой в работе модели трещиноватого слоя, исследованы волновые явления, образующиеся в результате взаимодействия сейсмических импульсов с трещиноватыми коллекторами.
Ключевые слова: разрывный метод Галеркина, контактное условие скольжения, флюидонасыщенная трещина, отраженные волны.
D.N. Voroshchuk1, V. A. Miryaha1'2, I. B. Petrov1'2, A. V. Sannikov1'2
Moscow Institute of Physics and Technology (State University) 2Institute of Computer Aided Design of the Russian Academy of Science
Seismic response investigation of a cluster of subvertical cracks by the discontinuous Galerkin numerical method
The discontinuous Galerkin method for an unstructured mesh is developed and implemented for modelling wave responses from subvertical macrocracks in corbonate media for numerical simulation of a direct problem in seismic search. In this paper different mathematical models of cluster are simulated and compared in terms of seismic responses. This model differs in the approach to cluster representation: explicit allocation of cracks in the cluster region with parameters which either coincide or differ from the media. It is demonstrated that intercracks interaction can be seen in the presented in the work model cracks layer. Also, the work investigates wave phenomena created by seismic wave interaction with the cracks cluster.
Key words: DGM, exploration seismology, glide contact condition, fluid-filled crack, converted waves.
1. Введение
В последние годы методы численного решения прямых задач в геофизике приобретают все более важную роль из-за высокой стоимости поискового бурения. Вместе с тем растет необходимость в повышении детализации геометрии геологических срезов и расширения диапозона исследуемых частот. Следовательно, для решения актуальных задач сейсмической разведки необходимо использовать методы с высоким порядком точности. Исторически сложилоть так, что при моделировании сейсмических волн наиболее часто используется формулировка уравнений линейной теории упругости в виде дифференциальных уравнений второго порядка для перемещений среды, т.к. позволяет учесть основые
протекающие в задаче процессы. Большое количество исследований используют данную формулировку. В последнее время все чаще можно встретить альтернативную формулировку в виде гиперболической системы законов сохранения первого порядка. Основные трудности при численном моделировании волновых процессов в сейсморазведке представляет высокочастотная составляющая волн на больших временах. Для таких задач традиционные конечно-разностные методы показывают себя не лучшим образом Для методов низких порядков свойственна численная вязкость, вызывающая размазывание волновых фронтов. Для гиперболических систем уравнений разрабатываются численные методы дискретизации, позволяющие получать высокий порядок точности как по пространству, так и по времени, без расширения шаблона метода. Одним из таких численных методов является разрывный метод Галеркина.
2. Постановка численных экспериментов
При проведении численного эксперимента рассматривалась геометрия расчётной области, приведенная на рис. 1. Параметры области: ширина 7000 м, глубина 3000 м Размеры коллектора составляли 2000 м на 100 м. Глубина залегания коллектора равнялась 1500 м (см. рис. 1).
Область коллектора задавалась двумя различными способами: последовательностью параллельных макротрещин с совпадающими и отличающимися параметрами среды в области коллектора и за его пределами. Далее для краткости будем называть их первой и второй моделями. Для простоты интерпретации сейсмических откликов в предложенных моделях также будет рассмотрен отклик от модели однородного изотропного ограниченного пласта. Использовалась линейно-упругая модель среды со следующими параметрами: скорость распространения продольных волн для вмещающей среды - ср = 3200 м, поперечных - с8 = 1780 М, плотность породы - р = 2450 м?. Для случая однородного изотропного слоя его параметры: ср = 2400 м, с8 = 1330 м, Р = 2450 М3. Для второй модели коллектора параметры среды в области трещиноватости следующие: ср = 2880 м, с8 = 1600 м, р = 2450 ^.
г м3
Д {} п Трещиноватый М 1500 м. коллектор
3000 м. | | $100 М.
Г| <^2000 м.^> V
7000 м.
Рис. 1. Геометрические параметры модели
Логичным кажется начать более детальное рассмотрение всех способов задания трещиноватой области с самого простого - модели однородной изотропной среды. Конечно, следует упомянуть, что более точное соответствие эксперименту будет достигнуто в случае использования пространственно-независимых параметров среды. Примером такого подхода может служить [14]. Однако в рамках текущей статьи мы останемся рассматривать случай пространственно-независимых параметров среды коллектора, так как результаты будут использоваться исключительно для качественного сравнения моделей. Одним из способов задания коллектора является явное задание трещин в породе [1-4]. При проведении численного эксперимента с дискретной моделью рассматриваем случай параллельных трещин. Похожие исследования волнового отклика от системы распределенных случайным
образом и пересекающихся трещин представлены в [14,15]. Предположительно более корректной моделью коллектора является вторая, являющаяся, по сути, комбинацией подходов по выделению макротрещин и заданию однородного изотропного слоя. Для наглядности, при сравнении результатов будем в будущем использовать именно первую модель и однородный изотропный слой.
Для моделирования трещин и разломов в коллекторах использовалась модель бесконечно-тонкой флюидонасыщенной трещины [1]. В указанной модели не требуется сильно измельчать расчетную сетку, что в свою очередь позволяет избежать связанных с этим вычислительных затрат. Геометрические параметры дискретной и комбинированной моделей коллектора приведены на рис. 2.
Рис. 2. Расположение трещин в кластере
При построении расчётной сетки использовались специализированные библиотеки [5], [6]. Мелкость сетки (см. рис. 3) задавалась неравномерной для увеличения точности расчетов в области кластера трещин и экономии машинного времени.
Рис. 3. Пример построения расчётной неструктурированной сетки
В качестве исходной задавалась плоская волна с амплитудой в форме импульса Берлаге [2] длиной 200 м.
В сформулированной постановке численного эксперимента использовались геометрические параметры, соответствующие проведенным расчетам с помощью сеточно-характеристического метода [3], где исследовался волновой отклик от системы субвертикальных макротрещин. Детальное изучение отклика от одной еубвертикальной макротрещины рассмотрено в [1].
3. Контактное условие скольжения для разрывного метода Галеркина
В основе используемого в работе программно-вычислительного комплекса [7, 8] лежит разрывный метод Галеркина для неструктурированных сеток [9]. Метод имеет ряд положительных особенностей [10,11], основной из которых для авторов является возможность реализации высокого порядка аппроксимации по пространственным координатам и времени, что необходимо в силу волновой специфики решаемой задачи.
В данной работе в качестве системы базисных полиномов использовались полиномы Лагранжа 4-го порядка. Интегрирование по времени проводилось с помощью интегратора Дорманда-Принца с адаптивным шагом по времени.
Получим определяющие соотношения для контактного условия скольжения, лежащего в основе модели бесконечно-тонкой флюидонасыщенной трещины. Соотношения для любых контактных условий (полное слипание, скольжение) для разрывного метода Галёркина сводятся в конечном итоге к решению одномерной задачи Римана [12] при учёте дополнительных соотношений. В частности, для скольжения это непрерывность нормальной составляющей скорости и отсутствие сдвиговых напряжений на контакте. В [13] подробно расписано, откуда возникает необходимость решать задачу Римана при решении линейной системы уравнений упругости с помощью разрывного метода Галеркина и как её приближённо можно свести к одномерной. Там приведено решение только для случая полного слипания. Опуская подробности, приведём краткое описание численного метода и получим выражение для численного потока для условия скольжения в двумерном случае, чего ранее авторами не встречалось.
Система уравнений упругости в матричном виде в двумерном случае для изотропного пространства в переменных напряжения и скорости имеет вид [12]:
дир . диа „ диа , .
_Е + Ли ^ + ВРа ^ = 0, (1)
где и - вектор из 5 неизвестных переменных и = (ахх, ауу, аху, ух,уу)т. Здесь и далее предполагается суммирование по повторяющимся индексам. Собственные значения матриц Ара и Вра таковы:
«1 = -Ср, 82 = -С3, =0, = С3, ,в5 = Ср.
Область интегрирования разбита на треугольники т М и пронумерована. Предполагаем, что матрицы Ара и Вра постоянны внутри Т м . Решение в треугольнике Т М аппроксимируется линейной комбинацией из 1 (Ж+ 1)(Ж+2) не зависящих от времени полиномиальных функций Ф,(х,у) степени не выше N, образующих базис с носителем Т М , и зависящих от времени коэффициентов:
. („,(тА („п. л - Лт)
: (*,»,*) = йррт^у). (2)
После умножения 1 на базисную функцию Ф,, проинтегрировав по треугольнику
Т (т),
получаем
/ ф, ^ + / ф, (Ап ^ + в„ ^) ^ = 0. (3)
]Т(т) т УТ(т) V дх дУ )
Далее, применив формулу интегрирования по частям, получаем
' Ф, ^^ + ^ [ Ф^'^Б -
(4)
йУ = 0.
/Т(т) ^ ~1 У (дТ(т)) Р
' з
/т (Л 1ф,АРа Ыа + 9фвРа (^а)
Второе слагаемое появилось ввиду разрывности решения иь и матриц Ара, Вра на границе треугольника Т(т) в общем случае. Здесь Рр''' - численный поток через ]-е ребро треугольника в глобальной системе координат, а через (дТ(т)) . обозначены стороны треугольника
Т(т), з = 1, 2, 3 (см. рис. 4).
Пусть Т.т - матрица перехода в систему координат X'У (преобразование ортогонально), связанную с з-м ребром треугольника Т(т) (см. рис. 4). Тогда поток Рр можно найти приближенно:
Рр* =ТтзАа1 иьи (5)
где и^ - решение одномерной задачи Римана (см. рис. 5). На рис. 5 характеристики среды, отмеченные символом (-), соответствуют треугольнику Т(т), символом - соседнему по ,7-му ребру, - неизвестные состояния в задаче Римана в системе координат X'У'.
Запишем пределы решения слева и справа от рассматриваемого ребра ] (см. рис. 4):
м*) = №') чгфГ, (6)
Рис. 4. Элемент расчётной сетки
Рис. 5. Решение одномерной задачи Римана
Применяя для каждой характеристики условие скачка Ранкина—Гюгонио (Rankine-Hugoniot jump condition) и учитывая непрерывность нормальной составляющей скорости и отсутствие сдвиговых напряжений на контакте, получаем следующую систему уравнений:
u - ua = air^
ua - ub = Q^r^
uc - ud = «4Г+
ud - u+ = «5rs
h "xx _ _c "xx,
.b xy _ _c _ "xy = 0,
Vl = VCx = 0,
где см - неизвестные коэффициенты, Г1 - столбцы матрицы Кря, составленной из собствен'
л(т) Лтл)
ных векторов матриц Ард и Аря :
Rpq —
Решая систему уравнений 7, получаем выражение для ub:
ь
i х- + 2ß- 0 0 0 А+ + 2ц,+ \
X- 0 1 0 Х+
0 Р- 0 Р+ 0
с- 0 0 0 -с+
V 0 с- 0 —С- 0 /
u
a1
a.2
u — airi
СХ2Г2,
X
1 {Vx Vx P + (^xx ®xx)
c+pp+ + Cvpx
(9)
1
c-
(°xy \
c7p- )
Для случая полного слипания, когда выполняется условие непрерывности всех компонент скорости и тензора напряжения на контакте, коэффициент а\ совпадает со случаем скольжения, т.е. продольные волны распространяются в этих случаях одинаково. Коэффициент же а2 для полного слипания имеет вид
СУ.2
1 К — vy.)cxp+ + (^xy — <7+)
(10)
С- С+Р+ + С- р-
Видно, что а2 для скольжения (см. формулу 9) получается из формулы 10 путём зануле-ния продольной скорости и сдвигового напряжения в граничащей ячейке, т.е. с+ = а+у = 0, что равносильно контакту с акустической средой, в которой нет поперечных волн. Отсюда становится ясным физический смысл модели бесконечно-тонкой флюидонасыщенной трещины.
Объединие выражений 4, 5, 6, 9, 10 завершает построение численной схемы, более подробный вывод которой приведён в [7,9,13].
с
4. Результаты численных экспериментов
Результаты представлены в виде серии волновых картин для фиксированных моментов времени и сейсмограмм. В статье приводятся результаты проведенных расчетов для случая однородной изотропной среды и первой модели трещиноватого коллектора. На рис. 6 представлены расчётные сейсмограммы горизонтальной (слева) и вертикальной (справа) составляющей скоростей (см. рис. 1). Далее для краткости будем называть их vx- и vy-сейсмограммами соответственно.
На рис. 8 представлена волновая картина для фиксированного момента времени, соответствующая сейсмограмме 6. На рис. 9 представлена волновая картина для схемы при задании кластера первым методом, соответствующая сейсмограмме 7.
5. Анализ расчётных волновых картин и сейсмограмм
Рассмотрим подробнее приведенную на рис. 6 сейсмограмму. На vy-сейсмограмме отчетливо видны продольные р-волны, отраженные от верхней (II) и нижней (III) границы кластера. За ними располагается последовательность двукратных волн, отраженных от дневной поверхности (VII). В рамках рассматриваемой модели представляется возможным, исходя из значения скорости упругих волн в среде и интервала времени между фиксацией волн, полученного с сейсмограммы, вычислить толщину кластера, которая с высокой точностью совпадает исходной.
Рис. 6. Модель однородной изотропной среды
Рис. 7. Первая модель трещиноватого коллектора
На ух-сейсмограмме отчетливо видны два симметричных (в силу симметрии геометрических параметров задачи) волновых фронта гиперболической формы (I и VI), которые связаны с продольными и обменными волнами, дифрагированными от краев кластера, соответственно. Аналогичный эффект можно заметить и на уу-сейсмограмме, однако в сравнении с компонентами II и III его вклад невелик. Менее отчетливо указанные эффекты (I, VI) видны на сейсмограмме на рис. 7, где вследствие отсутствия границы раздела сред между кластером и вмещающей породой отражение происходит на краевых субвертикальных макротрещинах. Следует отметить не только количественное, но и качественное отличие указанных фронтов - при сравнении (I, VI) на рис. 6 и 7 наблюдаются перераспределения фазовых характеристик сигналов. В обоих приведенных случаях (рис. 6, 7) по расстоянию между вершинами фронтов I можно судить о расположении и горизонтальных размерах кластера.
На ух-сейсмограмме, изображенной на рис. 7, отмечены последовательности горизонтальных плоских волновых фронтов (IV), образованных в результате многократных отражений между трещинами кластера (обменные волны), по которым можно судить о трещиноватой структуре вмещающей породы. В реальности трудно детектировать обменные волны из-за наличия диссипаций в среде. Зависимость параметров периодической структуры обменных волн от длины падающей волны и геометрических характеристик трещиноватой области является предметом отдельного исследования.
Все приведенные в этой части статьи рассуждения подкрепляются соответствующими картинами волновых полей (рис. 8, 9), которые в сравнении с сейсмограммами несут не только дополнительную информацию, но и позволяют более наглядно представить картину происходящего. Авторы намеренно не опираются в своих рассуждениях на волновые картины, так как в реальной жизни нет возможности их получить при проведении геолого-разведывательных работ.
Рис. 8. Пример волновой картины для модели однородной изотропной среды
Рис. 9. Пример волновой картины для первой модели трещиноватого коллектора 6. Результаты
• Разрывный метод Галёркина на неструктурированных расчётных сетках адаптирован и реализован для моделирования отклика от субвертикальных макротрещин в карбонатных породах.
• Получены формулы для численного потока в случае контактного условия скольжения, для решения задачи контактного разрыва в модели флюидонасыщенной трещины.
• С помощью данного метода получены волновые картины в фиксированные моменты времени и сейсмограммы от системы субвертикальных макротрещин.
• Получены краевые волны, обусловленные взаимодействием падающего сейсмического импульса с краями трещиноватого пласта.
• Показана возможность детального исследования волновых процессов, инициированных падением сейсмического импульса на систему субвертикальных макротрещин, и решения прямых задач сейсморазведки с помощью разрывного метода Галёркина.
• Проведено сопоставление предложенных моделей. Показана возможность учёта межтрещинных волновых взаимодействий при моделировании трещиноватого слоя.
Следует отметить, что путём сравнения расчётных сейсмограмм с натурными можно делать выводы о геометрических и механических свойствах трещиноватых коллекторов.
Работа выполнена при поддержке гранта РФФИ 15-37-20673 мол_а_вед.
Литература
1. Kvasov I.E., Petrov, I.B. Numerical study of the anisotropy of wave responses from a fractured reservoir using the grid-characteristic method // Mathematical Models and Computer Simulations. 2012. V. 4, N 3. P. 336-343.
2. Левянт В.Б., Петров И.Б., Муратов М.В. Численное моделирование волновых откликов от системы (кластера) субвертикальных макротрещин // Технологии сейсморазведки. 2012. № 1. C. 5-21.
3. Муратов М.В., Петров И.Б. Расчет волновых откликов от систем субвертикальных макротрещин с использованием сеточно-характеристического метода // Математическое моделирование. 2013. V. 25, № 3. C. 89-104.
4. Левянт В.Б., Петров И.Б., Голубев В.И., Муратов М.В. Численное 3D моделирование объемного волнового отклика от систем вертикальных макротрещин // Технологии сейсморазведки. 2014. № 2. C. 5-21.
5. Shewchuk J.R. A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator. https://www.cs.cmu.edu/ quake/triangle.html
6. Karypis G. METIS, Serial graph partitioninge. 1998. http://www.cs.umn.edu/ karypis/metis
7. Миряха В.А., Санников А.В., Петров И.Б. Численное моделирование динамических процессов в твердых деформируемых телах разрывным методом Галеркина // Математическое моделирование. 2015. V. 27, № 3. C. 96-108.
8. Petrov I.B., Favorskaya A.V., Khokhlov N.I., Miryakha V.A., Sannikov A.V., Golubev V.I. Monitoring the state of the moving train by use of high performance systems and modern computation methods // Mathematical Models and Computer Simulations. 2015. V. 7, N 1. P. 51-61.
9. Kaser M., Dumbser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - I. The two-dimensional isotropic case with external source terms // Geophysical Journal International. 2006. V. 166, N 2. P. 855-877.
10. Hesthaven J.S., Warburton T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer. 2008.
11. Ладонкина М.Е., Неклюдова О.А., Тишкин В.Ф. Использование разрывного метода Галеркина при решении задач газовой динамики // Математическое моделирование. 2014. T. 26, № 1. C. 17-32.
12. LeVeque, R.J. Finite Volume Methods for Hyperbolic Problems // Cambridge University Press, 2002.
13. Wilcox L.C., Stadler G., Burstedde C., Ghattas O. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media // Journal of Computational Physics. 2010. V. 229, N 24. P. 9373-9396.
14. Saenger E.H., Shapiro S. Effective velocities in fractured media: A numerical study using the rotated staggered finite-difference grid // Geophysical Prospecting. 2002. V. 50, N 2. P. 183-194.
15. Saenger E.H., Kriiger O.S., Shapiro S. Effective elastic properties of randomly fractured soils: 3D numerical experiments // Geophysical Prospecting. 2004. V. 52, N 3. P. 183-195.
References
1. Kvasov, I.E. and Petrov, I.B. Numerical study of the anisotropy of wave responses from a fractured reservoir using the grid-characteristic method. Mathematical Models and Computer Simulations. 2012. V. 4, N 3. P. 336-343.
2. Leviant V.B., Petrov I.B. and Muratov M.V. Numerical simulation of wave responses from subvertical macrofractures system. Seismic Technologies. 2012. N 1. P. 5-21. (in Russian).
3. Muratov M.V., Petrov I.B. Simulation of wave responses from subvertical macrofracture systems using grid-characteristic method. Matem. Mod. 2013. V. 25, N 3. P. 89-104. (in Russian).
4. Leviant V.B., Petrov I.B., Golubev V.I. and Muratov M.V. 3D modelling of seismic responses from large vertical fractures. Seismic Technologies. 2012. N 2. P. 5-21. (in Russian).
5. Shewchuk J.R. A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator. https://www.cs.cmu.edu/ quake/triangle.html
6. George Karypis METIS, Serial graph partitioninge. 1998. http://www.cs.umn.edu/ karypis/metis
7. Miryaha V.A., Sannikov A.V., Petrov I.B. Discontinuous Galerkin method for numerical simulation of dynamic processes in solids. Matem. Mod. 2015. V. 27, N 3. P. 96-108. (in Russian).
8. Petrov I.B., Favorskaya A.V., Khokhlov N.I., Miryakha V.A., Sannikov A.V., Golubev V.I. Monitoring the state of the moving train by use of high performance systems and modern computation methods.Mathematical Models and Computer Simulations. 2015. V. 7, N 1. P. 51-61.
9. Kaser M., Dumbser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - I. The two-dimensional isotropic case with external source terms. Geophysical Journal International. 2006. V. 166, N 2. P. 855-877.
10. Hesthaven J.S., Warburton T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer. 2008.
11. Ladonkina M.E., Neklyudova O.A., Tishkin V.F. Application of the RKDG method for gas dynamics problems. Matem. Mod. 2014. V. 26, N 1. P. 17-32. (in Russian).
12. LeVeque R.J. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. 2002. (in Russian).
13. Wilcox L.C., Stadler G., Burstedde C., Ghattas O. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. Journal of Computational Physics. 2010. V. 229, N 24. P. 9373-9396.
14. Saenger E.H., Shapiro S. Effective velocities in fractured media: A numerical study using the rotated staggered finite-difference grid. Geophysical Prospecting. 2002. V. 50, N 2. P. 183-194.
15. Saenger E.H., Kriiger O.S., Shapiro S. Effective elastic properties of randomly fractured soils: 3D numerical experiments. Geophysical Prospecting. 2004. V. 52, N 3. P. 183-195.
Поступила в редакцию 29.10.2015.
