Двухфазная фильтрация в трещиновато-пористой среде Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

УДК 519.6

DOI: 10.18101/2304-5728-2019-2-104-115

ДВУХФАЗНАЯ ФИЛЬТРАЦИЯ В ТРЕЩИНОВАТО-ПОРИСТОЙ СРЕДЕ

© Цыдыпов Севан Гуро-Цыренович

старший преподаватель,

Бурятский государственный университет имени Доржи Банзарова Россия, 670000, г. Улан-Удэ, ул. Смолина, 24а E-mail: sivan77@mail.ru

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

Ключевые слова: пористый объем; анизотропная среда; двухфазная фильтрация; деформация; теория упругости; итерационно-разностная технология; математическая модель.

Для цитирования:

Цыдыпов С. Г.-Ц. Двухфазная фильтрация в трещиновато-пористой среде // Вестник Бурятского государственного университета. Математика, информатика. 2019. № 2. С. 104-115.

Введение

Деформация каркаса сильно влияет на течение жидкости и газа в трещиновато-пористой среде. Для построения вычислительной технологии мы опирались на математическую модель фильтрации газовой и жидкой фаз в анизотропных средах, построенную на основе трудов Ван Генухтена (Van Genuchten):

др^ш (1 - 5 )

дл

я " др ш5

дл

-

(^ г1 г г г рг (ург-ргЯ

т

г

- div

г

рр (?рр -р "я

т" 1

г

(1)

(2)

г г г\ г,0 р 0 ,

Р =Р Р = Р —тт, Р — атмосферное давление;

ш = ш

( г\ V

(р ) — пористость (-^ , У„ — объем пор);

5 — эффективная влагонасыщенность;

1К — тензор проницаемости, в общем случае полно заполненный; к<г =л! 1 - 5 (1 - 52) — относительная проницаемость для газовой фазы;

цг — вязкость газа; цг — источниковый член в балансе газовой фазы;

ц" — источниковый член в балансе жидкости;

рР

р" = р" 0 —от , р" 0 — плотность флюида при атмосферном давлении; Р

к" = л/5 (1 - V1 - 52) — относительная проницаемость для жидкой фазы; ц" — вязкость флюида;

р" = рг - р (5); рс (5) = р'

капиллярное давление; Р

Р я

пороговое давление, причем рй =-, где а — параметр, характеризую-

а

щий пористую среду.

Для построения численной модели будем использовать итерационно-разностный алгоритм и формулы расчета давлений газа, которые подробно расписаны в работе [5]:

рг,»= ■

а , + а , + а рг\,, + а рг,,, + а^р\,, + а р\,, + Ь

е± у-1,к ,1 ] +1,к ,1 ] ,к-1,1 п± ] ,к+1,1 Ь± ] ,к ,1 -1 ] ,к ,1 +1

(3)

5

а

Здесь

а = а + а + а + а + а + а +

(амX Ь = (аX (рУ

АЛ

АЛ

(4)

При этом будем иметь

2р°1 2р°1

^ (*+1 - *-1 )(*" ) ' ^ (*+1 - ^-1 )(*+1 - ' (5)

2р° 1 2р°+1

а (л+1-у,-1 )(У, -у«)' а (у,+1-У,-1 )(У,+1-У»)' ()

2РГ1 2р° 1

' (3+1 -3-1 )(3-3-1) ' ' (3+1 -3-1 )(з+1 -з)

Для упрощения записи у коэффициента в опускаем целые индексы. Значения в в дробной точке определяем через значения в целых точках:

РГ : ^

к- 2 2 ^

Деформация пористого пласта

Отметим, что предложенная разностная схема обладает большой универсальностью, поскольку позволяет решать задачи в предельно широком диапазоне изменения влагонасыщенности. Опытные зависимости для коэффициентов относительной проницаемости для жидкой фазы, так и для газовой среды имеют нулевые значения при = 0 и при = 1. Это приводит к многообразию ситуаций по влиянию окружения на разностное решение в точке. Из-за этого специализированные методы, ориентированные на конкретный тип уравнения, не справляются с решением задачи двухфазной фильтрации. Рассмотрим возникающие ситуации, характерные для плоской задачи фильтрации. В то же время пласт будем считать линейно деформируемой средой, следующей закону Гука. В случае отсут -ствия влияния массовых сил перемещения точек пористой среды будут зависеть только от коэффициента Пуассона. Этот коэффициент, в свою очередь, будет зависеть от влагонасыщенности пласта. Однако из-за отсутствия данных о влиянии влаги, находящейся в порах, на характер деформации пористого каркаса будем исходить из постоянной величины коэффициента Пуассона. При этом предлагаемая вычислительная технология без труда обобщается на случай переменной величины этого коэффициента. Таким образом, мы включаем в расчет напряжения, которые возникают в пористой среде из-за воздействия слоев, окружающей пласт среды.

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

В декартовых координатах уравнения равновесия в перемещениях выглядят так [8]:

V2и 50/ = 0, (8)

1 - 2у сх О х

VV + -±-50 +-В. / = 0, (9)

1 - 2п5у

V2у +—-+-Р / = 0, (10)

в = ди+ду + д^, (П)

Сх Су д3

где и, V, у — проекции смещения точки; х, у, 3 — декартовы координаты; р — плотность; /х, /, / — проекции массовой силы на оси выбранной системы координат; 0 — объемная деформация; О = Е—- — модуль уп-

2 (1 + п)

ругости при сдвиге; V — коэффициент Пуассона; Е — модуль Юнга;

52 52 52 сх2 ду2 д32

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

/ = / = / = 0.

•Ух ¿у 3

Подставляя (11) в (8)-(9) и раскрывая оператор Лапласа, получим следующую систему уравнений в перемещениях:

52и 52и 52и , ч 52V , ч 52у

Ь — + —-+ — + (Ь -1)-+ (Ь -1)-= 0, (12)

СХ2 ду2 д32 у ' дхду у ' дхдз

— + Ь-С^ + + (Ь -1) д2и +(Ь -1) д2у = 0 (13)

СХ2 + Су2 +С32 +( )дХду+( )Суд3" '

д2у д2у , д2у п л д2и п 1Ч д2v

—г + ^ + Ь-- + (Ь -1)-+ (Ь -1)-= 0, (14)

дх2 ду2 дз2 у дхд3 ' дуд3

, 2 - 2У где Ь =

V2 = +--2 +--2 — оператор Лапласа.

1 - 2у

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

Разностная схема

Умножая каждое из определяющих уравнений на г2 и аппроксимируя производные в уравнениях (12)—(14) симметричными разностями на неравномерной сетке, получим следующие соотношения:

а (и,]+1,к- и,],к)- ап (и,],к- и,]-1,к)+Ьау, (и+1,],к- и,],к)-

~Ьае (и,],к - и-1,],к) + аь (и,],к+1 - и,],к)- (15)

-ал (",■,],к - "/,],к-1 ) - ,],к + 5р1 = 0

ьа (V, ]+1,к- Ч] ,к)- Ьап (V ] ,к- V, ]-1,к)+ау (V+1,],к- V, ],к)-

(\],к- V-]) +а (V ,],к+1- \],к)- а (\],к- )- (16)

-Ч, ],к + 5р 2 =0,

а (у, ]+1,к- у, ] ,к)- ап (у ], к- у, ] -1,к)+ау (у,+1, ] ,к- у, ] ,к)-

-ае (У,] ,к - У-1,],к ) + Ьаь (уг,],к +1 - У ] ,к ) - (17)

-ьа (у ] ,к- у, ], к-1)+5р з = 0

где

22

а =7-V-V, ау = 7-V-\, (18)

(х - хм) (х+1 - хм) (х+1 - х) (х+1 - X-1)

2 2

п (У]- У]-1)(У]+1 - У]-1Г 5 (У]+1 - У])(У]+1 - У]-1)'

22

( 71 к 7к-1 )( 7к+1 2к-1 ) ( 7к+1 7к )(7к+1 7к-1 )

(20)

^ 1 =

р1

(Ь - 1)(V

■ + 1, ] + 1,к Vi+1, ]-1,к Vi-1, ] + 1,к

+ Vi-1, ]-1,к )

(Ь - 1)(у

(х+1 - х-1)(У]+1 - У]-1)

■ +1, ], к + 1 - у ■ + 1, ], к-1 - у-1, ], к + 1 + у ■ -1, ], к-1)

( +1 -1)( 7 к +1 7 к -1)

(21)

2 =

(Ь - 1)(и

■ +1, ] + 1,к 1, ]-1,к и -1, ]' + 1,к

+ -1,]-1,к )

(х+1 - х-1)(У]+1 - У]-1) (Ь -1)(,]+1,к+1 - ,]+1,к-1 - ,]-1,к+1 + wi,]-1,к-1) (У]+1 - У]-1)(7к+1 - 7к-1)

(22)

(Ь - 1)(и , + - и , , - и , . , + и , , ,)

^ _ V /V 1 +1, у, к + 1 1 +1, у, к-1 1 -1, у, к +1 1 -1, у, к-1 /

( Х1 +1 - Х1 -1)( 3к +1 - 3к-1) (23)

+ (Ь - 1)(V,у +1,к + 1 - V,у +1,к-1 - V,у-1,к +1 + V,у-1,к-1)

( у у+1 - у, -1)( 3к+1 - 3к-1)

5

Шаблон разностной сетки (стандартный для регулярной сетки)

Стрелками показаны направления возрастания индексов. Из разностных уравнений (15) - (17) получаем следующие рекуррентные формулы для определения значений перемещений во внутренних точках:

и

а и1 , ,, + ауи1 ,,, , + Ьапи1 , , к + Ьаи,, ,, + аи ,, , + и,. ,, ,, + я ,

е 1, у-1,к у 1, у+1, к п 1 -1, у ,к $ 1+1, у ,к ? 1, у ,к-1 Ь 1, у ,к +1 р1

1,у,к ар + ап (Ь-1) + а$ (Ь-1)

(24)

V

IV ,, + Ьа V + а V , ., + а V. + ., + ау. ., , + а V ., + + $ 0

е 1, у-1,к у 1, у +1,к п 1 -1, у ,к $ 1 +1, у, к ? 1 ,у, к-1 о 1, у ,к+1 р 2

1, У ,к

, у-1,к 1 , у +1,к 1 апу1 -1, у ,к 1 +1, у ,к 1 ЬаГу1, у ,к-1 1 ЬаЬу1, у ,к+1 1 3

ар + ае (Ь - 1)+ ^ (Ь - 1)

а у ,. ,, + а,.,у + а^,. , ,■,, + ау,, ,■,, + Ьаж ,, . + Ьау ,,,, + $„

у ,, _ -

1 ,У,к ар + а( (Ь -1) + аЬ (Ь -1)

, (25)

, (26)

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

Здесь

ар _ ае + ау + ап + + а + аЬ .

Формулы (24)-(26) реализуют решения пространственных уравнений Ляме. Используя их вместе с разностными граничными условиями можно решить все задачи линейной теории упругости в декартовых координатах на неравномерных сетках. Таким образом, с помощью итерационного метода мы решаем задачу упругого деформирования пласта.

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

1. Если все ае, ап, ам> 0, задача является эллиптической в точке р.

О

р

о

п

-о И'

2. Если один из коэффициентов равен нулю, например ап=0, задача параболическая в точкер, причем ось параболичности вертикальная.

О

о

-а

-О и'

Ж

3. Ш-эллиптичность. По одному из направлений отсутствует связанность задачи, в то время как по другому направлению присутствуют двусторонние связи.

4. Гиперболичность с характеристиками, совпадающими с направлением координатных осей

5. Ультрагиперболичность. Характеристики совпадают между собой.

s

о

а—Ж11

6. Полное отсутствие влияния окружения. Стабилизация состояния в точкер при продвижении процесса по итерационной оси.

S

о

Ж

Ж

Я-(V—еЖ-9р-Ж1'

Ж

ж

п

e

Результаты расчетов

Описанная ранее вычислительная технология была применена к решению задачи об упругом вдавливании базальтовой плиты в прямоугольную область нефтеносного пласта. Предполагается, что пористая структура пласта, заполненная нефтью, при небольших деформациях ведет себя как упругая среда, подчиняющаяся закону Гука. Квадратный пласт 100 х 100 м2 и толщиной 4 м в центральной части подвергается действию плиты, которая вдавливается на глубину 0,1 м.

На рис. 1, 2 представлены распределения нормальной и касательной компоненты (одной из двух касательных, которые из-за глобальной симметрии условий деформирования ведут себя схожим образом). Как видно из рис. 1, продавливание породы в центральной части приводит к глобальному выталкиванию всего остального фона и локальному несимметричному выходу материала с левой стороны рисунка. Причем эта асимметрия при использованном направлении перебора внутренних узлов всегда появлялась в одном месте.

0.07 -р"

О.Об 4

О.Ов-О.ОУ 0.05-0.06 I 0.04-0.05 I о.оз-о.оа I 0.02-0.03 0,01-0.02 I 0-0.01 I -О-ОТ-О

Рис. 1. Нормальная компонента перемещений точек пласта в сечении 3 = 2 м

Область по вертикали(в направлении 3 координаты) составляет 4 м, поэтому 3 = 2 м — это среднее сечение пласта.

I 0.02-0.03

I 0.01-0.02 10-0.01 1 -0.01-0

I -0.02-0.01 I -0.03-0.02

Рис. 2. Касательная компонента перемещений точек пласта в сечении 3 = 2 м

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

2300 2200 2100 2000 1900 1800

■ 2200-2300

■ 2100-2200 ■ 2000-2100

■ 1900-2000

■ 1800-1900

Рис. 3. Изменение плотности при вдавливании квадратной базальтовой

области в пласт

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

Рис. 4. Изолинии давления газа при вдавливании квадратной базальтовой области по центру.

На западной и восточной стенках находятся окна проницаемости

Из-за разности проницаемостей индивидуальных компонентов в задачах просачивания флюида и выдавливания метана наблюдается относительно быстрое прохождение газовой фазы. Это отражается в том, что даже в стационарных задачах быстрее происходит установление параметров газовой фазы. На рис. 4, 5 показаны распределения давления в среднем сечении пласта в двух формах: в виде изолиний (рис. 4) и поверхности (рис. 5). Почти симметричная форма графиков определяется симметрией в постановке задачи и свидетельствует о завершении процесса установления.

Рис. 5. Поверхность, характеризующая распределение давления газа

Заключение

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

Литература

1. Shikuo C., Tianhong Y., Chenhui W. Displacement Mechanism of the Two-Phase Flow Model for Water and Gas Based on Adsorption and Desorption in Coal Seams // Materials of Int. Symposium on Multi-field Coupling Theory of Rock and Soil Media and Its Applications. Chengdu City, China. 2010. P. 597-603.

2. Van Genuchten M. Th. A closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils // Soil Sci. Soc. Am. J. 1980. Vol. 44. P. 892-898.

3. Schaap M. G., van Genuchten M. Th. A modified Mualem-van Genuchten formulation for improved description of the hydraulic conductivity near saturation // Va-dose Zone J. 2006. Vol. 5. P. 27-34.

4. Бубенчиков А. М., Цыдыпов С. Г. Построение математической модели двухфазной фильтрации в анизотропных средах // Геометрия многообразий и ее приложения. Улан-Удэ: Изд-во Бурят. гос. ун-та, 2014. С. 75-81.

5. Хан Х. Теория упругости. Основы линейной теории и ее применения. М.: Мир, 1988. 344 с.

TWO-PHASE FILTERING IN A FRACTURED POROUS MEDIUM

Sevan G.-Ts. Tsydypov Senior Lecturer,

Dorzhi Banzarov Buryat State University 24a Smolina St., Ulan-Ude 670000, Russia E-mail: sivan77@mail.ru

There are two main problems in describing two-phase filtration flows: the construction of physical processes model and the complexity of accounting for the processes of gas- and liquid-phase filtration. The article presents a universal computational technology that implements an arithmetical model of two-phase filtering of the processes of gas and fluid filtration in a fractured porous medium with consideration of the horizontal seam deformation. The Van Genuchten model for fluid and methane movement in coal seams is used as a basic model. A model of linear elastic deformation is constructed for accounting the influence of deformations. The universality of the obtained computing technology allows us to generalize the arithmetical model further. This technology takes into account the effect of water saturation on the elastic parameters, as well as the appearance of new fractures in the porous framework. The used differential scheme allows solving problems in an extremely wide range of changes in water saturation.

Keywords: porous volume; anisotropic medium; two-phase filtering; deformation; elasticity theory; iterative differential technology; arithmetic model.

References

1. Chen Shikuo, Yang Tianhong, Wei Chenhui. Displacement Mechanism of the Two-Phase Flow Model for Water and Gas Based on Adsorption and Desorption in Coal Seams. Materials of Int. Symposium on Multi-field Coupling Theory of Rock and Soil Media and Its Applications. Chengdu City, 2010. Pp. 597-603.

2. Van Genuchten M. Th. A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J. 1980. V. 44. Pp. 892-898.

3. Schaap M. G., van Genuchten M. Th. A Modified Mualem-Van Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation. Va-dose Zone J. 2006. V. 5. Pp. 27-34.

4. Bubenchikov A. M., Tsydypov S. G. Postroenie matematicheskoi modeli dvukhfaznoi filtratsii v anizotropnykh sredakh [Construction of an Arithmetic Model of Two-Phase Filtering in Anisotropic Media]. Geometriya mnogoobrazii i ee priloz-heniya. Ulan-Ude: Buryat State University Publ., 2014. Pp. 75-81.

5. Khan Kh. Teoriya uprugosti. Osnovy lineinoi teorii i ee primeneniya [The Theory of Elasticity. Fundamentals of Linear Theory and Its Applications]. Moscow: Mir Publ., 1988. 344 p.