Молекулярно-динамическое моделирование фильтрации Текст научной статьи по специальности «Физика»

М.Н. Овчинников

Казанский государственный университет Marat.Ovchinnikov@ksu.ru

МОЛЕКУЛЯРНО-ДИНАМИЧЕСКОЕ МОДЕЛИРОВАНИЕ

ФИЛЬТРАЦИИ

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

Введение

Целью данного исследования является определение соотношения между скоростью фильтрации жидкости в пористой среде и градиентом давления. Для стационарных фильтрационных потоков оно известно как закон Дарси:

цг =--Ур

И '

(1)

где № - скорость фильтрации, р - давление, к - проницаемость среды, и - вязкость жидкости. А как будет выглядеть ситуация для нестационарных потоков? Ответ на этот вопрос актуален как для интерпретации результатов натурных исследований нестационарными гидродинамическими методами насыщенных жидкостями пористых пластов, так и при анализе фильтрациижидкостей в мезоскопических системах.

Ряд исследователей (Христианович, 2000; Молокович и др., 1980) для описания нестационарных потоков используют феноменологические модели фильтрации с уравнением эволюции в виде:

(2)

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

Модель

Исследуемая нами система представляла собой прямоугольный «ящик», состоящий из трех отсеков. На границах этого ящика располагались жестко зафиксированные частицы. Отсек II (центральный) представлял собой, собственно, модель двумерной пористой среды, в которой размещались ячейки пористого скелета и поровые ячейки с частицами жидкости. Пример нескольких таких ячеек, располагающихся на границе между первым и вторым отсеками, показан на рис.1. Черным цветом снизу и сверху показаны непроницаемые границы. Точки в пределах голубых ячеек (поровых каналов) показывают местоположение частиц жидкости, темные прямоугольники обозначают элементы пористого скелета. Первый и третий отсеки заполнены только жидкостью. Отсек III соприкасается, с одной стороны, с пористой средой (отсеком II), а с противоположной от него стороны находилась жесткая стенка. Первый же отсек, с одной стороны, также соприкасается с порис-

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

Взаимодействие частиц жидкости между собой осуществлялось посредством потенциала вида

и= а (г - р)2

(3)

где параметр р принимал обычно значение 1, а параметр а изменялся в пределах от 0 до100 (обычно а =72). Частицы пористого скелета взаимодействовали с частицами жидкости посредством потенциала Леннарда-Джонса

17 = 4е

Г, N12

а

(4)

где г - расстояние между частицами, а параметры потенциала (4) принимались равными е =1, а =2"1/6. Массы всех частиц принимались равными 1.

В пределах пористого отсека II мы выделяли прямоугольные подсистемы (пример которой показан на рис. 1 красным цветом) с площадью Ах^-к, где Ах* - ширина к-го прямоугольника в направлении х, которое совпадает с направлением движения поршня (показано стрелками) и, в соответствии с геометрией задачи, направлением фильтрационного потока, к - высота прямоугольной подсистемы. Обычно Ах* включало в себя от 1 до 8 фильтрационных ячеек. Сами ячейки представляли собой квадратные прямоугольники с размерами от 2х2 до 10х10.

При молекулярно-динамическом моделировании методом частиц (Хокни и др., 1987; Хеерман, 1990), в результате решения уравнений Гамильтона, мы вычисляем в любой момент времени координаты каждой /-ой частицы и компоненты ее скорости (ух,уу). Временной шаг интегрирования был равен 2-10 3. Далее мы рассчитываем

N /

V* = Ы* - среднюю скорость частиц жидкости в к-

1=1 / к ой подсистеме Ах*-к, где И, - полное число частиц в выбранном интервале Ах*-к в момент времени t (обычно Ю2в начальный момент времени t=0). V* принимается далее за скорость фильтрационного потока в применяемом нами мезоскопическом подходе. Отметим, что всюду в данном изложении используется истинная скорость фильтрационного потока V , а не применяемая, зачастую, для удобства фиктивная скорость фильтрации Ж=V т, где т - пористость блока II.

научно-технический журнал

1 (15) 2004 I еоресурсы

Аналогом давления выступала концентрация частиц рР (площадная плотность частиц жидкости) в каждом к-ом элементе Ахк-к и выиислялась как рк = ¡N1, где Ы* - текущее число частиц в элементе Ах^-й, а Ыд - начальное невозмущенное (нормировочное) число частиц. Градиент концентрации частиц (градиент "давления") выгаислялся как Vр к=

.-ТТЛ (Рш-Ры +Г хк-1*- Началь-

I ные значения координат час-

Рис. 1. Пример модельной тиц жидкости и их скоростей пористои среды (импульсов) выбирались та-

ким образом, чтобы суммарные значения компонент скоростей V*, У* в пределах каждой ячейки были равны нулю,

ЛГ» _ __'

также, как и моменты импульсов X * у<, а суммарные зна-

г

чения кинетической и потенциальной энергии были бы одинаковы для каждой ячейки, что исключало возникновение «нефильтрационных» внутренних перетоков частиц между элементами системы в начальный момент времени. В процессе эволюции системы, по мере продвижения поршня, У* - компонента приобретает ненулевое значение, компоненты же V* колеблются вблизи нулевых значений.

Отметим, что в исследуемой системе движение поршня с постоянной скоростью Уп не приводит к постоянному расходу жидкости на границе I и II отсеков: число частиц жидкости, пересекающих в единицу времени границу между отсеками, изменяется. Некоторое время такой расход равен нулю, затем возмущение концентрации, создаваемое в I блоке, достигает границы I и II блоков, далее расход непрерывно нарастает.

Результаты

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

0 2 4 6 8 10 12 14 16 1

Рис. 3. Изменение во времени скорости V* для прямоугольныгхпод-систем Ах^-И, где ширина прямоугольника Ах охватымает ячейки с координатами 0 < х <4(а), 4 < х<8(б),8 < х <12 (в), 12 < х <16(г).

0 2 4 6 8 10 12 14 16 I

Рис. 2. Изменение во времени концентрации р* для прямоугольные подсистем Ахк-И, где ширина прямоугольника Ах охва-ты1вает ячейки с координатами 0 < х <4 (кривая а), 4 < х < 8

(б), 8 < х< 12 (в), 12 < х < 16 (г).

тамиО < х <4, 4 <х <8, 8 <х < 12,12 < х < 16.Мывидим,что рост концентрации частиц наблюдается, начиная с некоторого порогового момента Х ~ 4, когда возмущение концентрации частиц в отсеке I достигает границы с отсеком II.

Аналогично выплядит (рис. 3)и картина для х-компонент средних скоростей У* этих участков. Поначалу распространяется возмущение типа ударной волны с уменьшающейся амплитудой по мере распространения вглубь отсека II. Скорость распространенияэтойволны У^ ~10, чтосуществен-но больше скорости фильтрационного потока У*.

Со временем устанавливается относительно стабильный фильтрационный поток, возникает распределение концентрации частиц, убывающее по мере удаления от границы раздела !иП блоков. Отметим, что по мере продвижения вглубь пористого образца ударные волны рассеиваются, и их амплитуда достаточно быстро спадает на расстояниях Ь ~20. В целом, наблюдаемая картина соответствует эффектам рассеяния упругих колебаний на неоднородностях среды.

Таким образом, в модельной пористой системе возникает чередование изменений концентраций и скоростей потоков, связанный друг с другом, при этом между изменениями этих величин в пределах выделенной ячейки пористого пространства возникает временная задержка между V* и V/?*. Это видно из рис. 4, несмотря на значительные флуктуации у* и V/?*, что связано с относительной малостью числа частиц и малостью размеров исследуемой системы. При этом можно выщелить мелкомасштабные флук-

1,2"

1,0

х 0,8

•з Й

0,6

> 0,4

0,2-

0,0

7 8 9 10 11 12 13 14 1

Рис. 4. Временная зависимость скорости фильтрации Ук и градиента концентрации Vдля прямоугольной подсистемы,I при 0< х<4.

научно-технический журнал

Георесурсы 1 (15) 2004

0,0-

■ЯА.

1,6

Ж-1'4 р 1,2

2, 1,о *

£ °.8

0,6 ** 0,4

0,2 0,0-|"1

0,0

V—Vp (кривая а) и К + т" т*=0.3 (б).

0,2 0,4 0,6 0,8 1,0 1,2 {гга<1р

Рис. 5. Зависимость скорости фильтрации от градиента концентрации для прямоугольной подсистемы Ахк-И при 4 < х < 8.

туации рассматриваемых величин с резкими изменениями их значений в диапазоне времен Д^ < 10"2 и крупномасштабные флуктуации, охватывающие значительные изменения рассматриваемых величин на временах Д^ ~ 0.1 — 0.5.

На рисунке 5 показаны корреляции значений расчетных скоростей Ух в определенные моменты времени и соответствующих этим моментам времени значений градиентов концентрации Vрх- Мы видим, что, хотя и с учетом довольно больших разбросов значений вследствие флуктуаций, бульшим значениям градиентов соответствуют большие значения скоростей потоков, и можно предположить, что после статистического усреднения по значительному количеству вычислительных экспериментов будет получена зависимость между Ух и Vркх, близкая к линейной, соответствующей закону Дарси.

Рассмотрим подробнее временные корреляции между значениями скоростей Ух и градиентами концентрации Уркх для прямоугольной подсистемы АXе-к при 4<х<8. Для этого рассчитаем корреляционные функции следующего вида:

о,о од 0,4 0,6 0,8 l.o gradp Рис. 7. Линии трендов зависимостей .8V_ dt

Рассмотрение корреляционных функций, аналогичных (5), но для компонент скоростей фильтрационного потока, измеряемых в соседних прямоугольных подсистемах, показывает, что такие корреляции существуют на расстояниях Ь <10 и исчезают на больших расстояниях.

Теперь вернемся к рис. 5 и построим новый график не только в переменных

дУ

<->Ур при У&Ур, но и в переменных^ + г*— <-» Ур

1

N;

(5)

N -

где Ур\ = Ур* (,) - Урх У* ({) = У* (/) - У* нормировочный множитель, угловые скобки означают усреднение по времени t, а за средние значения величин

градиентов плотностей Урх (?) и скоростей Ух (?) в момент времени t принимаются значения линейных трендов этих величин, построенных на протяжении временных интервалов АtT~1. Мы видим, что максимальная корреляция наблюдается на временах т~0.2 - 0.4.

Аналогичная картина наблюдается и для уменьшенного участка Ах1, охватывающего прямоугольную подсистему с координатами 4< х <6. Максимум корреляций наблюдается при т~0.3. Данное значение интервала времени можно принять за характерное время релаксационного процесса т*.

Заметим, что в макроскопических системах могут существовать более острые пики корреляционных функций

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

8V(x,t)

8t

dt

легко вычисляется по из-

при т*= 0.3. Производная вестным значениям У(х^). Чтобы полученная картина не выглядела как хаотический набор точек, подобно рис. 5, изобразим на рисунке только линии трендов для зависимо' дУ

стей У^Ур иУ + т* — <-> Ур (рис. 7). Ширина линий трен-

8t

да:

показанных на рисунке, пропорциональна величине дУ

флуктуаций V и V + т*—для фиксированных значений У р. ' ' dt

Мы видим, что в новых координатах корреляционная дУ

зависимость Ур <-> V + для случая т*=0.3 значительно более близка к линейной, нежели для случая т*=0, особенно при больших значениях Ур. Данное обстоятельство, на наш взгляд, может служить дополнительным подтверждением необходимости учета релаксационных процессов в условиях нестационарной фильтрации, по крайней мере, на малых временах наблюдения в мезоскопических системах.

Литература

Христианович С.А. Избранные работы. М.: изд-во МФТИ, 2000. Молокович Ю.М., Непримеров H.H., Пикуза В.И., Штанин А.В. Релаксационная фильтрация. Казань: изд-во Казанского ун-та, 1980.

Соболев С.Л. Локально-неравновесные процессы модели процессов переноса. Успехи физических наук. 167. 1997. №10. 1095-1106.

Хокни Р., Иствуд Дж. 'Численное моделирование методом частиц. М.: Мир, 1987.

Хеерман Д.В. Методы компьютерного эксперимента в теоретической физике. М.: Наука, 1990.

М.Н. Овчинников

ДИНАМИКА ЖИДКОСТЕИ1 И КОНТРОЛЬ РЕСУРСОВ ПОДЗЕМНОЙ ГИДРОСФЕРЫ

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

специалистов в области подземной гидродинамики.

ISBN 5-98180-1 10-7 Казань: Казанский государственный университет, 2004. - 140 с

■— научно-технический журнал

i (15) 2004 I еоресурсы