Динамика газожидкостной среды в угольном пласте Текст научной статьи по специальности «Математика»

5. Математическое моделирование и обработка данных

УДК 532.546

© А.М. Бубенчиков, В.Б. Цыренова, С.Г. Цыдыпов ДИНАМИКА ГАЗОЖИДКОСТНОЙ СРЕДЫ В УГОЛЬНОМ ПЛАСТЕ

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

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

A.M. Bubenchikov, V.B. Tsyrenova, S.G. Tsydypov DYNAMICS OF GAS AND LIQUID MEDIA IN A COAL SEAM

The problem in this work represented by the equations (1)-(2) is complex, including for numerical methods, as a type of governing equations changes in subareas of physical area which are surrounding areas S=0 and S=1. So, the most universal iterative calculating technology has been chosen, it has allowed to obtain the concrete result without appeal to the type of initial equations.

Keywords: porous volume, anisotropic media, two-phase filtration, moisture saturation, capillary pressure, iterative differential technology.

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

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

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

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

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

др0т (1 - 5 )

Ы

-

( 1Кк0

V

р0

(-Р°Г)] = ; (1)

др^т , (1Кк^ р

а

(( -р^г)] = др ; (2)

дг

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

V / —

I с\ V

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

5 - эффективная влагонасыщенность (0-1);

1К - тензор проницаемости, в общем случае полнозаполненный;

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

о

ц - вязкость газа;

д0 - источниковый член в балансе газовой фазы;

зы;

др - источниковый член в балансе жидкости;

р0 рр

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

кРг =45 (1 - 41 - 52) - относительная проницаемость для жидкой фа-ц^ - вязкость флюида;

52

рр = р0 - рс (5) ; рс (5) = ра--капиллярное давление;

р

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

а

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

Одномерное приближение

Предельно упростим ситуацию. Примем, что источники отсутствуют, то есть ц_=^=0. Будем рассматривать одномерное движение, в этом случае У = /—, = — . Кроме того, примем, что движение происходит в

дх дх

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

др_т (1 - л) =д_( К_ _ _ р

dt

dpFms д dt дх

дх

fdpG W

Ц

дх

f KkL

ц F

V г- V /у

f F W

др дх

(3)

(4)

/у

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

G ,0 F ,0

G р G F р F

Р =—гр , р =~гр •

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

дрG aG -д- V дрСл Л

дt дх дх у

F F др a д f eF V дрF Л

дt дх дх у

(5)

(6) (7)

-G,° kG

G ,0 F ,0

aG = m(1 - s, aF = ms, pG = K

V V р° р0

рG, вF = K Р

р0 Ц

F рР ,

(8)

р - атмосферное давление, р_0 - плотность газа при атмосферном давлении, рр0 - плотность жидкости при атмосферном давлении.

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

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

„_„_ я ( Л

(9)

дt

д_

дх

PG др

дх

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

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

аА -аРР° + +ь = о, (10)

2Р° 1 2Р°1

1-2 1+2

где ае =т-V?-^ а* = 7-V?-V (11)

(х1+1- х1-1) (х1- х1 -1) (х1+1- х1-1) (х1+1- х1)

(11)

а° (< )п-1 (р; Г

а = а + а +-1, Ь = —\-JJ- . (12)

Р е ™ АГ АГ У '

Коэффициенты в этом разностном уравнении определены на новом слое по времени, но с использованием значений на предыдущей итерации.

Выражая из (10) р . , найдем

а рй, + а р°+, + Ь ро = ^+1-, ( = 2,..., N -1). (13)

«р

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

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

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

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

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

4) делаем порядка N глобальных итераций. После этого считается законченным расчет на очередном слое по времени;

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

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

О = аеР°-1Л + а*.Р°+1Л + а*Р^-1 + апРС1к+1 + Ь (14)

^ ар '

Ь {аР,к Г (Р^ Г (15)

где а = ае + + ах + ап + 1, Ь = ' ; / 1 ' , (15)

р АГ АГ

а также на пространственный случай

с _ ^Р-Щ + ^Л+Щ + 1,1 + Щщ + аьРщ-х + а>Рщ+1 + Ь (16)

Рм ар '

Здесь

К к1 Г (*%, )п 1 (р-и V

ар _ ае + ак + аз + ап + аь + аг + у '' , Ь = А . (17)

р & АГ

При этом ае и будут определяться формулами (11), а для ай ап, аь, а(

будем иметь

2рл 1 2рл 1 к — к+— - _ 2 ап _--^--, (18)

(ук+1- Ук-1) (ук- Ук-1У n {Ук- Ук-1) (Ук- Ук У

2вг; 1 2pG i l— l+— 2 at =---2--. (19)

(+1- -1)( - -1)' ' (+1--1 +1-^)

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

Рк-1 _ 2 .

2 ^

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

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

Тестирование алгоритма при постоянной насыщенности

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

Ь = 6 м - длина трубы;

р^ = 1000 кг/м3; рЛ0 = 1,27 кг/м3;

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

/ = 8,9 ■ 10-4 Пас; = 1,78 ■ 10-5 Пас.

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

5 = 0,8;

Рв°х = 1,2 • 105 Па;

5 = 0,3; Рв0ых = 105 Па.

Рис. 1. Одномерная область двухфазной фильтрации

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

ф° = КкО (5) д

(

дГ ц°т (1 - 5 ) дх

дР

дх

(20)

На рис. 2 представлены результаты решения этого уравнения рассматриваемым численным методом. В случае установившегося течения, когда

дРо

дг

■ = 0, решение легко находится аналитически. На графике аналитиче-

ское решение практически совпадает с численным, что говорит о корректности применения рассматриваемого численного метода для решения данных уравнений.

рв, кПа

аналитическое решение численное решение

4.5 5 5.5 X, М

Рис. 2. Распределение давления газа по длине трубы в разные моменты времени

В рассматриваемом примере влага поступает с левого конца трубы за счет скачкообразного увеличения насыщенности в этой части канала от значения 5=0,3 до значения 5=0,8. Поэтому профиль давления газа, изна-

чально имеющий вид «левого нижнего уголка», постепенно деформируется в почти линейное распределение.

Численный эксперимент

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

Рис. 3. Распространение волны насыщенности по объему пористого пространства

Как видно из рис. 3, наблюдается волновой механизм распространения насыщенности по пористому пространству трубы. Причем зона изменения насыщенности с течением времени увеличивается. Характерная выпукло-вогнутая форма ее распределения со временем становится более простой и к моменту, определяющему установившееся состояние течения, кривая s = s(x) имеет вполне определенный знак кривизны.

Литература

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 / M. Th. van Genuchten // Soil Sci. Soc. Am. J. - 1980. - Vol. 44. - P. 892-898.

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

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

Бубенчиков Алексей Михайлович, доктор физико-математических наук, профессор, заведующий кафедрой теоретической механики ММФ Томского государственного университета, е-mail: Aleksy121 @mail.ru

Цыренова Валентина Бабасановна, доцент, доктор педагогических наук, заведующий кафедрой геометрии Бурятского государственного университета, е-mail: v.ts@mail.ru

Цыдыпов Севан Гуро-Цыренович, преподаватель кафедры информационных технологий Бурятского государственного университета, е -mail: sivan77@mail.ru

Bubenchikov Aleksey Mikhalovich, doctor of physical and mathematical sciences, professor, head of the department of theoretical mechanics, MMF Tomsk State University, е -mail: Aleksy121@mail.ru

Tsyrenova Valentina Babasanovna, doctor of pedagogical sciences, head of the department of geometry, Buryat State University, е-mail: v.ts@mail.ru

Tsydypov Sevan Guro-Tsyrenovich, teacher, department of information technologies, Buryat State University, е -mail: sivan77@mail.ru

© Ж.Г. Дамбаев

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ ВОЛН НАПРЯЖЕНИЙ ПРИ ВЗРЫВЕ СИСТЕМЫ ЗАРЯДОВ, РАСПОЛОЖЕННЫХ В РЯД

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

Ключевые слова: математическое моделирование, взрыв системы зарядов.

Zh.G. Dambaev

MATHEMATICAL MODELING OF VOLTAGES WAVES INTERACTION AT EXPLOSION OF CHARGES SYSTEM ARRANGED IN A ROW

In this work the results of numerical solution of the problem of explosion of charges system arranged in a row are presented. A technology of sparing explosion of model and field experiments has been shown.

Keywords: mathematical modeling, explosion of charges system.

Постановка задачи

Для решения задачи взрыва системы зарядов, расположенных по линии, применяется численный метод решения уравнений динамической теории упругости. Рассмотрим взаимодействие волн напряжений между