CC BY
97
УДК 550.834, 539.3
Численное моделирование распространения сейсмических волн в среде, содержащей одиночную пору или пористый слой
М.М. Немирович-Данченко, А.А. Шатская
Национальный исследовательский Томский политехнический университет, Томск, 634050, Россия
В работе рассмотрено влияние пористости сред на упругие волновые поля. На основе численного моделирования построена дифракционная картина при прохождении плоской волны через одиночную пору в карбонатных породах. Свойства матрицы (кальцит и доломит) и флюида (вода) моделируются по изображению тонкого среза керна. Проводится качественное сравнение с имеющимися расчетными данными. Изучено влияние пористой среды на сейсмическое поле для случая ансамбля пор. Рассмотрены значения пористости от 2 до 20 %. Для сопоставления с экспериментальными данными рассмотрена модель пористого спеченного алюминия Al-6061. Обработка результатов численного моделирования позволила оценить средние скорости в модели пористого алюминия и сопоставить их с данными физического моделирования. Показано, что расчеты в целом дают хорошее качественное (для одиночной поры) и количественное (для ансамбля пор) соответствие результатов.
Ключевые слова: пористые среды, сейсмические волны, численное моделирование
Numerical modeling of seismic wave propagation in a medium with
a pore or porous layer
M.M. Nemirovich-Danchenko and A.A. Shatskaya
National Research Tomsk Polytechnic University, Tomsk, 634050, Russia
The influence of the porosity of media on elastic wave fields was studied. Numerical modeling was performed to construct a diffraction pattern during passage of a plane wave through a single pore in carbonate rock. The properties of the matrix (calcite and dolomite) and fluid (water) were modeled using the image of a thin drill core section. A qualitative comparison to available computational data was carried out. The influence of a porous medium on the seismic field was studied for an ensemble of pores. The considered porosity range varied from 2 to 20%. A model of porous sintered aluminum Al-6061 was used for comparison to experimental data. When processing the numerical modeling results, we estimated average velocities in the porous aluminum model and compared them to physical modeling data. The calculations provide a good qualitative (for a single pore) and quantitative (for an ensemble of pores) correspondence of results.
Keywords: porous media, seismic waves, numerical modeling
1. Введение
Интерес к пористым средам и их влиянию на поле упругих волн не ослабевает уже несколько десятилетий. Это связано и с задачами материаловедения, и с разведкой недр при поисках углеводородов. Две эти различные области знаний становятся родственными при физическом и математическом моделировании. Для лабораторного моделирования сейсмических волн в пористых средах в качестве скелета используются различные материалы: алюминий [1, 2], натриево-известковое стекло [3], гипс [4] и другие материалы [5]. В лабораторных экспериментах оцениваются как статические, так и динамические характеристики, основными из которых яв-
© Немирович-Данченко М.М., Шатская A.A., 2015
ляются скорости продольной и поперечных волн. Создание одной лабораторной модели пористой среды позволяет решить целый ряд задач [5], при этом удается создавать модели со сравнительно несложной заменой одного флюида другим [4]. В то же время присущие методам лабораторного моделирования ограничения, прежде всего по выполнению требований подобия, определяют актуальность математического моделирования [6].
Существует широкий спектр методов математического моделирования — от полуаналитических (лучевой метод и методы, основанные на лучевом) до конечно-разностных (подробный обзор методов, их достоинств и недостатков приведен в [7]).
Нами используется конечно-разностный метод, позволяющий вычислять полное волновое поле с использованием полной системы уравнений механики сплошных сред [8]. Применяется явная условно-устойчивая схема второго порядка, позволяющая при использовании ряда усовершенствований [9] проводить расчеты для большого числа длин волн (>1001).
Численное моделирование требует соотнесения результатов расчетов с имеющимися экспериментальными или теоретическими результатами. Анализируя численные результаты при расчете волновых полей в пористых средах, необходимо оценивать изменения скоростей распространения волн (количественный анализ). Для целей качественного анализа подходят картины дифракции, сравнение фронтов волн, поляризации колебаний и т.п.
2. О постановке задач
Моделирование полного волнового поля выполняется для двух задач, двух разных моделей сред и в двух постановках.
В первой постановке рассмотрена модель одной крупной поры, заполненной жидкостью, и распространение упругих волн в такой модели. Для этой постановки целью расчета является сравнение общего волнового поля с имеющимися расчетными результатами (т.е. сравнение качественных картин).
Во второй постановке рассматривается падение плоской волны на пористый слой со случайным распределением пор. В этой модели в ходе расчета оценивается, как за счет пористости изменится скорость распространения продольной волны. Расчет проводится для нескольких значений пористости, для двух из которых имеются данные лабораторного моделирования.
Для одиночной поры расчет сравнивается с результатами конечно-разностного моделирования [10]. Для расчета волнового поля в случае среды с ансамблем пор была сгенерирована модель, эквивалентная пористому алюминию (спеченный А1-6061), для которого имеются количественные данные физического моделирования [2].
Обе задачи решаются линеаризованным конечно-разностным методом Уилкинса [7-9].
2.1. Задача об одиночной поре
Было рассмотрено влияние одиночной поры на проходящую сейсмическую волну. Одна из недавних работ [10] посвящена расчету дифракции сейсмического поля на одиночной поре для различных частот. В ней приведен срез изображения образца карбонатной породы и упрощенная матрица упругих свойств для проведения расчетов.
Нами решена задача о распространении плоской упругой волны в модели образца из работы [10]. На рис. 1 дана геометрия первой задачи, показаны модель, содержащая одну крупную пору, и вид падающего сиг-
нала. Сигнал (упругая волна) взят в виде затухающей синусоиды с длиной волны, примерно в 10 раз меньшей характерной толщины поры L (показана на рисунке). В расчете размер поры L составил 200 расчетных ячеек, длина волны 1 = 20 расчетных ячеек.
В такой постановке основной особенностью начала дифракционного процесса является обтекание фронта падающей волны вокруг поры. Материал поры (флюид) обладает существенно более низким импедансом. Совокупность преломлений-отражений на отдельной поре приводит к хорошо различаемому эффекту локального замедления фронта волны.
На более поздней стадии (рис. 2, а), когда волна прошла через пору, отразилась от границы модели и вновь прошла через пору, волновая картина в целом становится почти хаотичной. Остаются небольшие участки с цугами продольных волн (показаны на рис. 2 литерами Р). Для сравнения на рис. 2, б приводится рисунок из работы [10].
Видно, что пора в обоих случаях насыщена коротко-периодными колебаниями (чередование черных и белых тонов на рисунках — это экстремумы волн), берега поры отчетливо прорисованы сейсмическим полем. На рис. 2, а момент времени характерен для начала формирования того псевдохаотичного поля, которое наблюдается в работе [10]. В целом решение первой задачи показало, что отдельная пора отчетливо проявляется в волновом поле вследствие контрастности упругих свойств, при этом волновому полю присущи явно выраженные картины дифракции. Тем самым было выполнено качественное тестирование расчетной программы, которая далее применена для решения следующей задачи.
2.2. Задача о пористом слое
Рассмотрим задачу о падении плоской упругой волны на среду, содержащую слой, в котором случайно заданы поры. Рассмотрен пористый алюминий А1-6061, т.к. именно для этого материала в литературе имелись данные лабораторного моделирования [2]. Геометрия задачи показана на рис. 3. На границу АВ падает плоская
Рис. 1. Матрица скоростей продольных волн в образце [6], стрелкой показано направление падения плоской волны в виде затухающей синусоиды, L — характерный размер поры
р
Рис. 2. Волновая картина: результат расчета (а), расчет [10, рис. волн в обоих расчетах
. Буквой Р обозначены явно выраженные цуги продольных
продольная волна скорости смещения. Среда до линии Е¥ однородная, с параметрами алюминия А1-6061: Vp = = 6260 м/с, V = 3080 м/с, р = 2700 кг/м3, после линии Е¥ задается случайная пористость, свойства флюида (воздух при комнатной температуре): Ур = 331 м/с, V = = 0 м/с, р = 1.225 кг/м3.
В отличие от предыдущей задачи, необходимо оценить перемещение экстремума волны вдоль модели, поэтому в качестве исходного сигнала взят импульс Рикера
и&у I^ = РЦ) = -( - ^ )в~2п/_
где f — частота; to = 1/ / — период импульса; и у — ^-компонента скорости смещения.
Поскольку цель данного расчета — количественное сравнение с результатами физического моделирования [2], то это определило выбор размеров пор и характерных длин волн. Размер поры порядка 20 • 10-6 м, частота акустического сигнала 5 МГц, длина волны —0.1 см. Таким образом, размер поры всего в 50 раз меньше длины волны, такие поры можно описывать несколькими расчетными ячейками. На рис. 4 приведено распределение пор по размерам, минимальный размер — 1 расчетная ячейка, максимальный для пористости 17 % — 10 ячеек.
При распространении плоской волны в такой среде поведение отдельной поры будет подобно поведению поры из предыдущей задачи при подобных частотах.
Групповой эффект от ансамбля пор очевиднее всего сказывается на кинематическом параметре — скорости распространения упругих волн. Для спеченного алюминия опубликованы данные физического моделирования [2] и рассмотрены два случая пористости — 9 и 17%. Для этих значений получены следующие лабораторные оценки скорости продольной волны: Ур (9 %) = = 5125 м/с, Ур (17 %) = 4120 м/с. Для сравнения скорость продольной волны для материала А1-6061 ^ (0 %) = 6260 м/с.
В расчетах скорость распространения волны определялась по углу наклона годографа преломленной волны в плоскости (I, у) (I — время, у — вертикальная координата). Необходимо отметить, что верхняя часть модели (рис. 3) и скелет нижней части заданы одним материалом — А1-6061. Нижняя часть отличается только наличием пор и флюидом в этих порах. Очевидно, что изменение скорости распространения волн, определяемое по углу наклона годографа, вызвано только значением пористости и свойствами флюида.
Была проведена серия расчетов для одного и того же флюида (воздух при комнатной температуре, константы приведены выше) и нескольких значений пористости (от 2 до 20 %). Результаты расчетов приведены на рис. 5. На рисунке квадратами выделены значения из экспериментальной работы [2], ромбами — результаты данной работы.
Рис. 3. Геометрия задачи о волновом поле в пористой среде
Рис. 4. Распределение пор по размерам
Рис. 5. Зависимость скорости продольной волны от пористости
Отметим, что в расчетах были получены следующие значения: при пористости 9 % значение скорости продольной волны Vp = 5123 м/с, а при пористости 17 % значение скорости продольной волны Vp = 4074 м/с. Таким образом, при численном моделировании значения скоростей в пористом алюминии отличались от скоростей из эксперимента при пористости 9 % всего на 0.04 %, а при пористости 17 % отличие составило 1.1 %.
3. Обсуждение результатов
Одним из технологичных методов геофизики является сейсморазведка. Она проводится в диапазоне частот, редко превышающим 100 Гц (длина волны составляет первые десятки метров). При обработке данных сейсморазведки используются, кроме всего прочего, результаты исследований в скважинах — данные акустического каротажа. Эти данные получаются на частотах звукового (0.5-15 кГц) и ультразвукового (2050 кГц, 0.3-2.0 МГц) диапазонов (длина волны в последнем случае — несколько миллиметров). По результатам измерений строят геоакустические модели разрезов скважин для интерпретации данных сейсморазведки, проводят оценку пористости продуктивных_пластов, определяют упругие модули горных пород и т.д. (см., например, [11]).
Таким образом, частоты и результаты при физическом моделировании сопоставимы с частотами и данными акустического каротажа, а последний служит для анализа результатов сейсморазведочных работ. С другой стороны, упругие свойства керна и породы в скважине определяются почти исключительно по скоростям распространения продольной и поперечной волн. Численное моделирование позволяет, меняя свойства скелета и замещая один флюид другим, строить модели, наиболее близко отвечающие как данным каротажа, так и сейсморазведки.
В работе решены две задачи о сейсмическом поле в пористых средах: задача о дифракции волнового поля
на одиночной поре и задача о прохождении упругой волны через пористый слой со случайно распределенными порами. Дифракционная картина в случае одиночной поры качественно хорошо соответствует расчетным результатам, опубликованным другими авторами. В случае пористого слоя по результатам численного моделирования были оценены средние скорости продольной волны в слое, значения скоростей сопоставлены с данными физического моделирования.
Представленный нами подход может использоваться на этапе изучения геологического разреза в тех случаях, когда для оценки петрофизических параметров требуется знать геометрию пористой среды и величину пористости. Обычно эти значения определяются по данным керна и на основании исследований в скважине с учетом атрибутов сейсморазведки. Численный расчет позволяет воссоздать поле упругих колебаний по заданным параметрам пористой модели и сравнить с реальным сейсмическим разрезом, проверяя таким образом корректность определения породы-коллектора на качественном и количественном уровне.
Литература
1. Pyrak-Nolte L.J. Fracture anisotropy: the role of fracture-stiffness gra-
dients // The Leading Edge. - 2007. - V. 26. - No. 9. - P. 1124-1127.
2. Bonnan S., Hereil P.-L., Collombet F. Experimental characterization of quasi static and shock wave behavior of porous aluminum // J. Appl. Phys. - 1998. - V. 83. - No. 11. - P. 5741-5749.
3. Li X., Zhong L., Pyrak-Nolte L.J. Physics of partially saturated porous media: Residual saturation and seismic-wave propagation // Annu. Rev. Earth Planet. Sci. - 2001. - V. 29. - P. 419-460.
4. Ошкин A.H. Изучение распространения упругих волн в средах с цилиндрической симметрией методами лабораторного моделирования / Дис. ... канд. физ.-мат. наук. - М.: МГУ, 2009. - 118 с.
5. Ивакин Б.Н. Методы моделирования сейсмических волновых явле-
ний. - М.: Наука, 1969. - 287 с.
6. Жукова Т.М. Методика моделирования волновых полей в неоднородных средах на основе численного решения прямой динамической задачи сейсморазведки / Дис. ... канд. техн. наук. - М.: УлГТУ, 1984. - 193 с.
7. Немирович-Данченко М.М. Модель гипоупругой хрупкой среды и ее применение в сейсмике / Дис. ... докт. физ.-мат. наук. - Новосибирск: ИГ СО РАН, 2004. - 217 c.
8. Немирович-Данченко М.М., СтефановЮ.П. Применение конечно-
разностного метода в переменных Лагранжа для численного расчета волновых полей в сложнопостроенных средах // Геология и геофизика. - 1995. - № 11. - С. 96-105.
9. Немирович-Данченко М.М. Численное моделирование трехмерных
динамических задач сейсмологии // Физ. мезомех. - 2002. - Т. 5. -№ 5. - С. 99-106.
10. Nurhandoko B.E.B., Wardaya P.D., Adler J., Siahaan K.P. Seismic wave propagation modeling in porous media for various frequencies: A case study in carbonate rock // AIP Conf. Proc. - 2012. - V. 1454. -P. 109-112. - doi 10.1063/1.4730699.
11. Горная энциклопедия / Под ред. Е.А. Козловского. - М.: Сов. энцикл., 1984. - Т. 1. - 560 с.
Поступила в редакцию 23.03.2015 г., после переработки 21.04.2015 г.
Сведения об авторах
Немирович-Данченко Михаил Михайлович, д.ф.-м.н., проф. ТПУ, [email protected] Шатская Александра Андреевна, ассист. ТПУ, [email protected]
