Фильтрация газо-жидкостной среды в плоской горизонтальной области

2014

Математика и механика

№ 6

МЕХАНИКА

УДК 532.546

А.М. Бубенчиков, В.Б. Цыренова, С.Г. Цыдыпов

ФИЛЬТРАЦИЯ ГАЗО-ЖИДКОСТНОЙ СРЕДЫ

В ПЛОСКОЙ ГОРИЗОНТАЛЬНОЙ ОБЛАСТИ

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

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

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

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

Эти результаты позволяют распространить теорию двухфазной фильтрации на течения в анизотропных средах и указать метод проведения и интерпретации экспериментов для определения коэффициентов в тензорных связях.

В настоящей работе используется математическая модель двухфазной фильтрации в пористых пластах, сформулированная на основе анализа работ [1-4]:

Зр°т

д/

дрртя д/

--

-

Г К?

Ъ РО (О-р°й )) = дО; Р? ( -Ррй)) = д?;

О О ( О \ ъ 0 Р 0

Р = Р (Р ) = Р ' —г , Р - атмосферное давление; Р

V

т = V - пористость ( Vп - объем пор, V- общий объем);

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

к<Ъ = л/ 1 - 5 (1 - 52) - относительная проницаемость для газовой фазы;

дО - вязкость газа; дО - источниковый член в балансе газовой фазы; Ц -источниковый член в балансе жидкости;

Р =Р

?,0 Р

р , - плотность флюида при атмосферном давлении;

к? =\[5 (1 - V1 - 52) - относительная проницаемость для жидкой фазы;

д - вязкость флюида;

? О с / \ с / \ ^

р = р - р (5); р (5) = р

л/Т-

5

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

, Р я

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

а

Плоская задача изотропной фильтрации

Примем, что источники отсутствуют, то есть дО = д? = 0. Будем решать пло-

- - д д д д скую задачу, в этом случае V = /--+ у— , шу =--1--. Кроме того, примем,

дх ду дх ду

что движение происходит в горизонтальном направлении: (я)х = 0, (я) = 0. В этом случае уравнения (1), (2) можно переписать следующим образом:

дРОт (1 - 5) = д

д/ дх

дР?т5 = д

Кк!

г РО

д/

дх

V Ц

Кк?

V дх у/

д 'ду

Кк

г РО

( дрО ^

\

Ц

д/ V дх У/

ду

V Ц

Кк?

ду

V ^ у/

Л Л

Ц

др? ду

Здесь К - коэффициент проницаемости, в настоящих расчетах скалярная характеристика.

Ъ

0

В уравнениях (3), (4) перейдем к давлениям рр. Для этого заменим рО и рЕ,

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

О,0 Е,0

„О Р „О Р р Е Р = ~Р , Р = ~Р .

Р Р

В результате получим два уравнения следующего вида:

дРО аО д/ Ч £ = ( РО V дРО Л дх у д + дУ ( вО V дРО \ дУ У

ъ Е Е дР а д/ д дх ( РЕ V дРЕ 1 дх у д + дУ ( РЕ V дРЕ ^ дУ У

0,0 Е,0 0,0 / О

„О____\р „,Е_„,„Р дО _ ТУ- р кг „О иЕ

= т(1 -5, аЕ = тя^——, р° = К-

Р0 ЦО

= К-

Р0

Р° - атмосферное давление; рО,° - плотность газа при атмосферном давлении; рЕ 0 - плотность жидкости при атмосферном давлении.

Итерационно-разностный алгоритм

Вернемся теперь к уравнениям (6), (7). Решение этой системы разберем на примере уравнения (6):

д/

_д_

дх

дР

о\

дх

д

"дУ

Рс

дР

о\

ду

Из зависимости р° от ро (см. (8)) уравнение (9) является нелинейным, поэтому решение его будет строиться с использованием итерационного подхода. Предварительно строим полностью неявную разностную схему для уравнения (9), которая, по опыту решения задач теплопроводности и диффузии, должна быть абсолютно устойчивой.

Аппроксимируя пространственные производные простейшими разностями на неравномерной сетке, а производную по времени разностями назад, найдем

аеР°1-\,к + ам,Р%\,к - аРР<О,к + а*Р1к-1 + апР°к+\ + Ь = 0 >

где

2Р°

1 - ?к

а„, = -

2Р°

1 + ~2'к

2рОк ! м - 2

ап = -

2вОк 1

1 + 2

а

а = а + а + а + а +

Р "е ^ ^ "я ^ п ^ д/

1,к .

ь =

(О,))1 1

д/

п—1

ае =

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

во = в % + Р у,к-1

в 1 " о ' М- 2

Коэффициенты в этом разностном уравнении определены на новом слое по времени, но с использованием значений на предыдущей итерации. Выражая из (10) рО, найдем

1 а„

О "еР°-1, к + а^Р } +1, к + а5Р у , к-1 + апР }, к+1 + Ь

Рук =---------.

а

Р

Заметим, что при вычислении источника Ь необходимо помнить значения величин на предыдущем слое по времени (в формулах (13) эти величины обозначены индексом п - 1 вверху).

Расчет с использованием итерационно-разностной технология (ИРТ) выполняется следующим образом:

1. Задаем начальные распределения давлений (газовой и жидкой фаз) по всей области.

2. С помощью граничных условий вычисляем значения плотностей на концах трубы на новом слое по времени.

3. Используя (14), последовательно перевычисляем давления во внутренних узлах расчетной области. Когда перебор точек по индексам 1 и к закончен, считаем одну глобальную итерацию завершенной.

4. Делаем порядка М = Ы-К глобальных итераций. Здесь Ы,К - количества шагов по отдельным координатным направлениям. После этого считается законченным расчет на очередном слое по времени.

5. Далее осуществляется переход к новому слою по времени, то есть возвращаемся к пункту 2.

С помощью (11) - (14) можно решать все плоские задачи двухфазной фильтрации в однородных анизотропных пластах.

Тестирование алгоритма

С использованием описанной математической модели были решены различные тестовые задачи, в частности одномерная задача.

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

р? = 1000 кг/м3; р°,° = 1,27 кг/м3;

К = 9,869233 • 10-13 м2;

= 8,9 • 10-4 Пас; дО = 1,78 • 10-5 Пас.

Было рассмотрено движение газо-жидкостной среды в трубе, заполненной пористым материалом, в частности песком (т = 0,4).

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

Результаты вычислений

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

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

Рх = 105 Па , рх = 6,9 -104 Па,

Рвь1х = 9-104 Па, = 5-104 Па .

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

Задача решалась как эволюционная с итерациями на каждом шаге по времени. Как видно из рис. 1, к моменту времени t = 106 с распределение давления газа становится симметричным, что при симметричном расположении окон проницаемости говорит о правильности работы вычислительной программы.

0 1,2 2,4 3,6 4,8 х

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

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

Рис. 2. Распространение давления флюида в плоской области с двумя окнами проницаемости

0 1,2 2,4 3,6 4,8 х

Рис. 3. Распределение влагонасыщенности в плоском горизонтальном пласте

ЛИТЕРАТУРА

1. Shikuo C. Displacement mechanism of the two-phase flow model for water and gas based on adsorption and desorption in coal seams / Chen Shikuo, Yang Tianhong, Wei Chenhui // 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 GenuchtenM.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils / M. Th. van Genuchten // Soil Sci. Soc. Am. J. 1980. Vol. 44. P. 892-898.

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

4. Никитин К.Д. Метод конечных объемов для задачи конвекции-диффузии и моделей двухфазных течений: дис. ... канд. физ.-мат. наук. М., 2010. 105 с.

Статья поступила 18.07.2014 г.

BubenchikovA.M., Tsyrenova V.B., TsydypovS.G. FILTRATION OF A GAS-LIQUID MEDIUM IN A PLANE HORIZONTAL REGION

In describing two-phase filtration flows, along with problems of constructing a model of a physical process, there are difficulties of the computational character as well, which are associated with the manifestation of different-scale filtration processes of phases. The paper presents a universal computing technology that implements the van Genuchten mathematical model. An example of two-phase filtration in a plane region is presented.

Keywords: porous volume; anisotropic medium; two-phase filtration; moisture saturation; capillary pressure; iterative-difference technique.

BUBENCHIKOV AlexeyMikhailovich (Doctor of Physics and Mathematics, Prof., Tomsk State University, Tomsk, Russian Federation) E-mail: alexy121@mail.ru

TSYRENOVA Valentina Babasanovna (Candidate of Physics and Mathematics, Prof., Buryat State university, Ulan-Ude, Russian Federation) E-mail: v.ts@mail.ru

TSYDYPOV Sevan Guro-Tsyrenovich (Buryat State university, Ulan-Ude, Russian Federation) E-mail: sivan77@mail.ru

REFERENCES

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, 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, vol. 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. Vadose Zone J., 2006, vol. 5, pp. 27-34.

4. Nikitin K.D. Metod konechnykh ob"emov dlya zadachi konvektsii-diffuzii i modeley dvukhfaznykh techeniy. Dis. kand. fiz.-mat. nauk. Moskow, 2010. 105 p.