Акустические волны в пористых флюидонасыщенных средах: компьютерное моделирование на мезоуровне Текст научной статьи по специальности «Физика»

Акустические волны в пористых флюидонасыщенных средах: компьютерное моделирование на мезоуровне

В.В. Крючкова

Томский филиал Института геологии нефти и газа ОИГГМ СО РАН, Томск, 634055, Россия

В статье приведен метод моделирования на мезоуровне среды с регулируемой пористостью. Модель создается посредством генерирования равномерно распределенных случайных чисел. Проводятся оценки теоретических и расчетных скоростей распространения волн в пористых средах. Теоретические значения были получены по теории Гассмана, а расчетные — из модельного эксперимента. Приводятся волновые поля и сейсмограммы.

1. Введение

Исследованию различных свойств пористых сред, насыщенных газом или жидкостью, в последнее время уделяется большое внимание. В том числе интерес вызывает акустика таких сред, то есть особенности распространения в них упругих волн. Это, с одной стороны, вызвано необходимостью отличать нефтегазосодержащие породы от пород, насыщенных другими флюидами. С другой стороны, для таких явлений, как оползни, актуальной является проблема прогноза, не в последнюю очередь опирающегося на сейсмические методы [1]. При анализе (интерпретации) волновых полей в средах, имеющих различную пористость и насыщенных различным флюидом, исследователи опираются на существующие упрощенные теории флюидонасыщенных пород.

Одной из самых распространенных является теория Гассмана [2]. Эта теория, как, впрочем, и другие, опирается на лабораторные измерения свойств сухой породы (сухого скелета). Зная упругие свойства сухой породы, на основе теории Гассмана можно предсказать значение модуля плоского деформирования при насыщении данной породы любыми флюидами. (Значение этого модуля определяет скорость распространения продольной волны в среде.) Нами предлагается определять свойства сухого скелета из численного эксперимента. Кроме того, необходимость численных расчетов возникает и в связи с изучением деформирования пористых сред в случае сложной геометрии и неоднородности среды.

© Крючкова В.В., 2000

В достаточно широком диапазоне скважинных и лабораторных акустических исследований используется частота от 104 до 106 Гц. Длина волны при этом будет для большинства горных пород заключена в пределах 0.5 мм - 50 см. Для длины волны порядка 50 см и более можно пренебречь размером отдельной поры и выполнить осреднение упругих свойств. Далее пористую среду считают однородной, и задачу решают на макроуровне. Но в сантиметровом диапазоне мы неизбежно оказываемся на мезоуровне и поры становятся отдельными структурными элементами среды. С точки зрения вычислительной для сантиметрового диапазона мы сталкиваемся с необходимостью задавать поры как элементы расчетной области. Это делает актуальной задачу создания (синтезирования) исходной модели пористой среды с регулируемой пористостью.

Расчеты нами будут проводиться в сравнении с теорией Гассмана. Изложим поэтому вкратце основные моменты этой теории.

2. Теория Гассмана

В теории Гассмана [2] основным предположением является то, что относительное движение жидкости и скелета настолько мало, что не оказывает влияния на распространение сейсмических волн во флюидонасыщенных средах.

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

хого скелета, константы флюида и константы флюидонасыщенной породы. Однородный изотропный упругий материал, из которого состоит скелет, имеет плотность р0, модуль всестороннего сжатия k0 и модуль сдвига ц0. Сухой скелет (т. е. тот же упругий материал, но пористый и без флюида) имеет пористость Ф = Уп/У (где Уп — объем порового пространства; V — весь объем образца), среднюю плотность pdry = (1 -Ф)р0, модуль всестороннего сжатия kdry, модуль сдвига ^dry и модуль плоского деформирования M ^ = kdry + (4/3)^^. Флюид (газ или жидкость), насыщающий поровое пространство, имеет плотность рf и модуль всестороннего сжатия kf; при этом цf = 0. Средние свойства флюидонасыщенной породы имеют плотность psat; модуль всестороннего сжатия ksat; модуль сдвига ^sat и модуль плоского деформирования Msat. Необходимо свойства флюидонасыщенной породы выразить через свойства флюида и скелета.

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

Psat =^f + (1 -Ф)Р0. (1)

Следующее предположение, что модуль сдвига в флюидонасыщенной среде не отличается от модуля сдвига в сухом скелете (но не в материале скелета!), вследствие малого воздействия флюида на твердую фазу:

M"sat = цdry • (2)

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

ksat k dry +

(І - k^/ко)2

Ф/кf + (І - Ф Vк0 - kdry /k0

(3)

Учитывая, что в общем случае М = к + (4/3)ц, можем записать аналогичную (3) формулу для модуля плоского деформирования:

M sat = M dry +

(І - k^/ к о)2

Ф/к f + (І -Ф)/ко - kd^/ ко2'

(4)

При низких частотах флюидонасыщенный материал ведет себя как однородное тело, поэтому плоские продольные и поперечные волны будут распространяться со скоростями:

(5)

Vp =4 M sat IP,

Vs = лКй/P . (6)

В формулах (3) и (4) есть неизвестное отношение kdry /k0 (собственно говоря, неизвестны величины kdry

или Магу). Не останавливаясь пока подробно на смысле этого отношения, заметим лишь, что в [2] предлагается для определения этого отношения использовать эмпирическую связь, полученную по измерениям в сухих песчаниках:

к^/ ко = (1 + 50Ф)-1. (7)

Расчеты по формулам Гассмана проводят следующим образом. Задаются некоторой пористостью; опираясь на какие-либо соображения (или на формулу (7)), находят модуль всестороннего сжатия сухого скелета кагу. Затем подставляют найденное значение в формулу (3). Предполагая, что модуль сдвига при насыщении флюидом не меняется (формула (2)), находят модуль плоского деформирования М5Л по формуле (4).

3. Некоторые особенности расчетов по теории Гассмана

Итак, для выполнения расчетов необходимо знать значение модуля всестороннего сжатия в сухом скелете кагу, а при переходе от формулы (3) к формуле (4) обязательно нужно знать значение модуля сдвига в сухом скелете Цйгу. Это составляет основную проблему в использовании теории Гассмана.

3.1. Для оценки величины кагу в работе [3] используется понятие “эффективной сжимаемости порового пространства сухого скелета” и приводится следующая зависимость:

І

kdry

=—+—

ко кФ

Здесь величина (к^)-1 имеет смысл сжимаемости сухого скелета, а (кф )-1 — эффективная сжимаемость пустого порового пространства сухого скелета. После преобразования получим выражение:

k dry

ко

І+^ Ф

кФ

(8)

Сравнение (7) и (8) показывает, что формула (7) из [2] есть лишь частный случай более общего соотношения (8), которое справедливо при любой пористости и для любой геометрии пор. Формулу (7) можно использовать только для песчаников и нельзя обобщать на остальные исследования с другими материалами.

3.2. Что касается определения величины цагу, то можно оценивать ее из лабораторных опытов по возбуждению в сухой породе поперечных волн. Авторами [4] предлагается альтернативный способ — оценивать цагу, не измеряя скорость поперечных волн. Ими предполагается, что в широком диапазоне значений пористости и для большого класса пород имеет место соотношение

M

dry

M

№

к й

#: А

■ ■■■■ ШШШ

й? £

д: £

Я й В В

Рис. 1.

что позволяет после определения (или измерения) модуля плоского деформирования М легко получать значение цЗатем по формуле kdry = М- (4/3)ц подсчитывается модуль всестороннего сжатия сухого скелета, что и дает возможность использовать формулы Гассмана (3) и (4).

4. Описание численного эксперимента

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

Расчет выполняется следующим образом: сначала строится модель мезоструктуры пористой среды (шаг 1), затем осуществляется собственно расчет волнового поля при возбуждении плоской волны в такой среде (шаг 2), после чего результаты расчета могут быть представлены графически в формате, удобном для последующей обработки (шаг 3).

4.1. Создание моделей пористых сред с заданной пористостью

Простейший способ задания пористой среды, использованный нами ранее, — располагать поры в шахматном порядке (рис. 1).

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

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

В качестве достаточно удобной модели мезострукту-ры пористой среды нами предложена модель среды со случайно распределенными порами. Алгоритм в общих чертах таков. Когда генератор случайных чисел доставляет очередное значение в интервале от 0 до 1, оно, будучи масштабированным соответствующим образом, определяет абсциссу “поры”, следующее значение определяет ее ординату. Заполнение матрицы заканчивается тогда, когда общая пористость достигает заданного значения. Под общей пористостью в численном методе понимается отношение количества расчетных ячеек, соответствующих порам, к общему числу расчетных ячеек всего материала. На рис. 2 приведена модель мезоструктуры среды с пористостью 63 %. Черным цветом обозначено поровое пространство, белым — скелет. В верхней части модели находится сплошная среда, такой однородный участок необходим для формирования импульса. Частоту и форму импульса можно тогда регулировать, об этом будет сказано ниже.

Растровый образ модели мезоструктуры среды используется далее для преобразования в формат входных данных для расчетной программы (см., например, [6]).

4.2. Задание исходного акустического сигнала

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

Рис. 2. Модель мезоструктуры среды с пористостью 63 %. Поры распределены случайным образом

сигналы приведены на рис. 3. Приведем аналитические выражения для этих сигналов [7, 8]: импульс Рикера:

у = 2к/4~е(-(/ - Т)) е-2(п/('-Т))2,

импульс с колоколообразной огибающей:

у = е-Р 2('-Т )2

sm(2я/(t - /о)),

импульс Берлаге:

у = /пе/ вт(2/), затухающая синусоида: у = е/ вт(2/).

Выше говорилось, что модель мезоструктуры среды — это входные данные для расчетной программы. Одним из основных параметров при этом является скорость распространения продольных колебаний в материале. Задавая значение этой скорости Кр и частоту сигнала /, мы тем самым определяем длину волны X = = ¥р! /. Например, при частоте 400 кГц и скорости звука в среде 6 000 м/с длина волны получится равной 1.5 см. Если шаг по пространству равен 0.01 см, то на одну длину волны придется 150 расчетных точек. Это обусловливает высокую точность расчетов, так как согласно работе [9] для модели упругой дискретной среды такой “протяженный” по количеству расчетных ячеек сигнал будет длительно распространяться, сохраняя свою форму и амплитуду.

4.3. Обсуждение результатов

Для расчетов был выбран кварц с константами: k0 = 37.9 ГПа, ц 0 = 44.3 ГПа, М0 = k0 + (4/3)ц 0 = = 97.0 ГПа, р0 = 2 650 кг/м3. Было рассмотрено два значения пористости 11 и 20 %.

Волновое поле рассчитывалось для сигнала с частотой 400 кГц. Этот сигнал в виде плоской продольной волны скорости смещения частиц среды возбуждался на верхней границе модели ЙВ, в однородной среде

Рис. 4. Геометрическая постановка задачи для модели мезоструктуры среды при пористости 11 %. Цифрами показаны: 1 - середина модели; 2 — датчики приема сигнала; 3 — плоская продольная волна

EDCF. Распространяясь вниз со скоростью Кр

Рис. 3. Виды исходных сигналов. Цифрами обозначены: 1 — импульс Рикера; 2 — импульс с колоколообразной огибающей; 3 — импульс Берлаге; 4 — синусоидальный импульс с затуханием

/М0/р0, плоская волна достигает границы раздела EF “однородная среда - пористая среда” в некоторый момент времени / . Отраженная волна распространяется вверх также со скоростью Кр0. Волна, идущая вниз, испытывает множественные отражения-преломления на встречающихся порах и является преломленной волной, распространяющейся с некоторой средней скоростью. На рис. 5 показаны волновые поля (так называемые численные снимки) в моменты времени / * (слева) и 2/* (справа). На левом фрагменте — неискаженная форма исходного сигнала, на правом фрагменте хорошо заметно изменение волны в результате дифракции на порах. В целом форма импульса сохраняется.

Начиная с момента времени / и до момента, когда волна достигнет нижней границы среды АВ, в программе определялась средняя скорость распространения продольной волны в пористой среде. Этот расчет проводился следующим образом. Вдоль оси 1 в каждой расчетной точке расставлялись (мысленно) датчики 2 (рис. 4). На уровне линий у1 и у2 фиксировалось время прихода экстремумов продольной волны. Затем по формуле V = (где £ — расстояние от точки у1 до у2,

£ — время прихода волны) находили скорость распространения продольной волны. Эта процедура производилась как для средней линии 1, так и для нескольких линий справа и слева от нее. Скорость распространения вычисляли как среднее арифметическое найденных значений. Такой способ определения скорости применялся и для сухого скелета, и для флюидонасыщенной породы. В случае сухого скелета считалось, что поры заполнены воздухом kшг = 0.00014 ГПа, ршг = 1.2 кг/м3. Для насыщенной среды в качестве флюида бралась вода kf = = 2.25 ГПа, рг = 1 000 кг/м3.

t = t* t = 2t*

Рис. 5. Фрагменты волновых полей для разных моментов времени. Пунктиром показано местоположение границы раздела

На рис. 6 приведены сейсмограммы вертикальной компоненты скорости смещения. На сейсмограммах по горизонтали отложено время. Каждая отдельная кривая (трасса) строится по результатам приема сигнала в каждом датчике 2 (рис. 4). Здесь хотелось бы обратить внимание на отличие волновых полей от сейсмограмм. Волновые поля (рис. 5) строятся в фиксированный момент времени во всех точках пространства. Сейсмограммы же состоят из отдельных трасс, каждая из которых строится во все моменты времени в одной точке пространства, и отражают эволюцию волновой картины. Способ построения сейсмограмм соответствует наблюдениям in situ с помощью приемников (датчиков). Способ построения волновых полей соответствует “фотографированию” поля.

Найденное значение скорости позволяло подсчитать значения M и k для сухой и насыщенной пород и провести затем сравнение с теорией Гассмана. Приведем подробнее результаты расчетов для пористости 11 и 20 %. Значения скоростей и модулей, оцененные в программе, будем называть расчетными, а значения, рассчитанные по формулам Гассмана, — теоретическими.

4.3.1. Пористость 11 %

Сначала выполнялся расчет волнового поля для продольной волны в сухом скелете. По сейсмограмме вертикальной компоненты скорости смещения была подсчитана скорость распространения продольной волны (см. рис. 6, а). Полученное значение скорости Fpdiy = = 5 080 м/с позволило оценить модуль плоского деформирования в сухом скелете для данной пористости M dry = 60.86 ГПа (здесь использовано значение плотности сухого скелета pdry = 2 358 кг/м3 , вычисленное по формуле (1)).

Имея значение модуля M dry, получаем по формуле (9) модуль сдвига ц dry = 27.79 ГПа. При этом значении цdry скорость поперечной волны будет: Vsdiy = 3 430 м/с. Это теоретическая оценка значения скорости поперечной волны.

С другой стороны, значение скорости поперечной волны в сухом скелете и, стало быть, значение модуля сдвига в нем ц ^ можно также получить из расчетов. Из сейсмограммы для поперечной волны получаем удгу _ з 640 м/с, а затем по формуле ц^ = (V^ )2 р^ находим ц_ 31.24 ГПа.

Отношение расчетного значения поперечной скорости V*/ = 3 430 м/с к теоретическому значению Vsdry _ 3 640 м/с составляет 5.8 %. Таким образом, для этой пористости предпочтительнее, по-видимому, пользоваться расчетными значениями ц ^.

Далее по формуле (8) получаем отношение: кдгу /к0 _ 0.507. А по формуле (7) это выражение получается кдгу /к0 _ 0.154. Такое несовпадение результатов показывает неприменимость формулы (7) к расчетам с кварцем. Формула (8) поэтому является необходимой для расчета скоростей в пористых средах по формулам теории Гассмана.

4.3.2. Пористость 20 %

По сейсмограмме вертикальной компоненты скорости смещения рассчитали скорость распространения

Р-волна S-волна

Рис. 6. Сейсмограммы для продольной (а) и поперечной (б) волн. Линиями помечены: сплошной — падающая волна; штриховой — преломленная волна. Наклоны линий соответствуют скорости распространения волны

продольной волны У^ = 3 740 м/с. Плотность сухого скелета для такого материала рйгу = 2 120 кг/м3. По формуле Мйгу = (Кргу)2рйгу получили модуль плоского деформирования в сухом скелете для данной пористости М йгу = 29.65 ГПа. Зная эту величину, находим модуль сдвига ц йгу = 13.54 ГПа. Можно получить теоретическую скорость поперечной волны по формуле (9): У^ = = 2 527 м/с.

Теперь для нахождения той же величины цйгу проводим компьютерный расчет. При этом скорость поперечной волны в сухом скелете получается из сейсмограммы такая: у^ = 2 508 м/с. Модуль сдвига равен ц*у = 13.33 ГПа.

Для пористости 20 % по расчетам получили отношение: к^у/к0 = 0.306. По формуле (7) это выражение получается кдг/1 к0 = 0.091.

Отношение расчетного значения поперечной скорости = 2 508 м/с к теоретическому значению = 2 527 м/с составляет 0.75 %. Анализируя результат, можно сделать вывод о том, что для значений пористости порядка 20 % лучше пользоваться теоретическим значением цагу, полученным по формуле (9).

5. Заключение

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

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

Литература

1. Голъдин C.B., Колесников Ю.И., Полозов C.B. Раснространение акустических волн в грунтах в условиях изменяющегося сдвигового нанряжения (внлоть до разрушения образцов) // Физ. мезо-мех. - 1999. - Т. 2. - № б. - С. 105-113.

2. Уайт Дж. Э. Возбуждение и раснространение сейсмических волн. - М.: Недра, 198б. - 2б1 с.

3. Mavko G., Mukerji T. Seismic pore space compressibility and Gass-mann‘s relation // Geophysics. - 1995. - V. 6G. - No. б. - P. 17431749.

4. Mavko G., Chan C., Mukerji T. Fluid substitution: estimation changes in Vp without knowing Vs // Geophysics. - 1995. - V. 6G. - No. б. -P. 175G-1755.

5. Крючкова B.B. Комньютерное моделирование сейсмических волновых нолей и возможности структурных ностроений // Структурный анализ в геологических исследованиях. Материалы Междунар. науч. семинара. - Томск: ЦНТИ, 1999. - С. 242-244.

6. Крючкова B.B., Немирович-Данченко M.M. Численное моделирова-

ние раснространения акустических волн в анизотронных средах // Физ. мезомех. - 1999. - Т. 2. - № 1-2. - С. 43-4S.

7. Nielsen P. Numerical modelling of seismic waves: on the elimination of grid artifacts // Norwegian Research Council, Department NTNF, Reference No. 1GG463/41G.

S. Cnравочник геофизика. - М.: Недра, 19бб. - Т. 4. - 750 с.

9. Cепян Л.И. Нестационарные унругие волны. - Л.: Судостроение, 1972. - 374 с.