CC BY
6
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: [email protected]
Цыдыпов Севан Гуро-Цыренович, преподаватель кафедры информационных технологий Бурятского государственного университета, е -mail: [email protected]
Bubenchikov Aleksey Mikhalovich, doctor of physical and mathematical sciences, professor, head of the department of theoretical mechanics, MMF Tomsk State University, е -mail: [email protected]
Tsyrenova Valentina Babasanovna, doctor of pedagogical sciences, head of the department of geometry, Buryat State University, е-mail: [email protected]
Tsydypov Sevan Guro-Tsyrenovich, teacher, department of information technologies, Buryat State University, е -mail: [email protected]
© Ж.Г. Дамбаев
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВЗАИМОДЕЙСТВИЯ ВОЛН НАПРЯЖЕНИЙ ПРИ ВЗРЫВЕ СИСТЕМЫ ЗАРЯДОВ, РАСПОЛОЖЕННЫХ В РЯД
Представлены результаты численного решения задачи взрыва системы зарядов, расположенных в ряд. Демонстрируется технология щадящего взрывания модельных и натурных экспериментов.
Ключевые слова: математическое моделирование, взрыв системы зарядов.
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.
Постановка задачи
Для решения задачи взрыва системы зарядов, расположенных по линии, применяется численный метод решения уравнений динамической теории упругости. Рассмотрим взаимодействие волн напряжений между
