CC BY
29
2014
Математика и механика
№ 6(32)
МЕХАНИКА
УДК 532.546
А.М. Бубенчиков, В.Б. Цыренова, С.Г. Цыдыпов
ФИЛЬТРАЦИЯ ГАЗО-ЖИДКОСТНОЙ СРЕДЫ
В ПЛОСКОЙ ГОРИЗОНТАЛЬНОЙ ОБЛАСТИ
При описании двухфазных фильтрационных течений, наряду с проблемами построения модели физического процесса, существуют сложности вычислительного характера, связанные с проявлением разномасштабности процессов фильтрации фаз. В работе представлена универсальная вычислительная технология, реализующая математическую модель Van Genuchten'a. Дан пример двухфазной фильтрации в плоской области.
Ключевые слова: пористый объем, анизотропная среда, двухфазная фильтрация, влагонасыщенность, капиллярное давление, итерационно-разностная технология.
Проблемы охраны окружающей среды выдвигают в качестве важнейшей задачу исследования миграции загрязнённой солесодержащей воды (флюида) в угольных пластах. Так как угольные пласты являются, как правило, неоднородными и имеют трещинноватую структуру, то требуются новые численные модели, позволяющие учитывать неоднородность и анизотропию пористого скелета, а также пространственный характер движения двухфазной среды в пласте.
Теория двухфазной фильтрации базируется на обобщенном законе Дарси, справедливом для медленной стационарной фильтрации несмешивающихся сред. Согласно обобщенному закону Дарси, неподвижная пористая среда и одна из подвижных фаз рассматриваются как некая фиктивная пористая среда, в которой происходит фильтрация другой фазы. Такая схематизация предполагает, что при медленном стационарном течении формируется равновесное распределение фаз, которое в процессе двухфазной фильтрации сохраняется статистически постоянным. В процессе двухфазной фильтрации формируется фиктивная пористость, состоящая из активных, соединяющихся между собой пор, образующих каналы, по которым происходит движение фаз из застойных зон, где фазы неподвижны, или находятся в состоянии медленного циркуляционного движения. Смачивающая и несмачивающая фазы движутся каждая по своей системе каналов. Движение каждой из фаз происходит под действием своего фазового давления, а проницаемость фиктивной пористой среды определяется своей фазовой проницаемостью.
Эти результаты позволяют распространить теорию двухфазной фильтрации на течения в анизотропных средах и указать метод проведения и интерпретации экспериментов для определения коэффициентов в тензорных связях.
В настоящей работе используется математическая модель двухфазной фильтрации в пористых пластах, сформулированная на основе анализа работ [1-4]:
Зр°т (1 - 5)
д/
дрртя д/
--
-
Г К?
Ъ РО (О-р°й )) = дО; Р? ( -Ррй)) = д?;
(1) (2)
О О ( О \ ъ 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 ^ у/
Л Л
Ц
др? ду
(3)
(4)
Здесь К - коэффициент проницаемости, в настоящих расчетах скалярная характеристика.
Ъ
0
В уравнениях (3), (4) перейдем к давлениям рр. Для этого заменим рО и рЕ,
входящие в левые части этих уравнений, по следующим формулам:
О,0 Е,0
„О Р „О Р р Е Р = ~Р , Р = ~Р .
Р Р
В результате получим два уравнения следующего вида:
дРО аО д/ Ч £ = ( РО V дРО Л дх у д + дУ ( вО V дРО \ дУ У
ъ Е Е дР а д/ д дх ( РЕ V дРЕ 1 дх у д + дУ ( РЕ V дРЕ ^ дУ У
0,0 Е,0 0,0 / О
„О____\р „,Е_„,„Р дО _ ТУ- р кг „О иЕ
= т(1 -5, аЕ = тя^——, р° = К-
Р0 ЦО
= К-
Р0
(5)
(6)
(7)
(8)
Р° - атмосферное давление; рО,° - плотность газа при атмосферном давлении; рЕ 0 - плотность жидкости при атмосферном давлении.
Итерационно-разностный алгоритм
Вернемся теперь к уравнениям (6), (7). Решение этой системы разберем на примере уравнения (6):
д/
_д_
дх
дР
о\
дх
д
"дУ
Рс
дР
о\
ду
(9)
Из зависимости р° от ро (см. (8)) уравнение (9) является нелинейным, поэтому решение его будет строиться с использованием итерационного подхода. Предварительно строим полностью неявную разностную схему для уравнения (9), которая, по опыту решения задач теплопроводности и диффузии, должна быть абсолютно устойчивой.
Аппроксимируя пространственные производные простейшими разностями на неравномерной сетке, а производную по времени разностями назад, найдем
(10)
аеР°1-\,к + ам,Р%\,к - аРР<О,к + а*Р1к-1 + апР°к+\ + Ь = 0 >
где
2Р°
1 - ?к
а„, = -
2Р°
1 + ~2'к
((1 х1-1 )(х1 Х]-1) ((1 х1-1 )(х1+1 х1 )
2рОк ! м - 2
(Ук+1- Ук-1 )(Ук- Ук-1)
ап = -
2вОк 1
1 + 2
(Ук+1 - Ук-1 )(Ук+1 - Ук)
а
а = а + а + а + а +
Р "е ^ ^ "я ^ п ^ д/
1,к .
ь =
(О,))1 1
д/
п—1
(11)
(12)
(13)
ае =
Значения диффузионного коэффициента в в дробной точке определяются как полусумма значений в целых точках, например
во = в % + Р у,к-1
в 1 " о ' М- 2
Коэффициенты в этом разностном уравнении определены на новом слое по времени, но с использованием значений на предыдущей итерации. Выражая из (10) рО, найдем
1 а„
О "еР°-1, к + а^Р } +1, к + а5Р у , к-1 + апР }, к+1 + Ь
Рук =---------. (14)
а
Р
Заметим, что при вычислении источника Ь необходимо помнить значения величин на предыдущем слое по времени (в формулах (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: [email protected]
TSYRENOVA Valentina Babasanovna (Candidate of Physics and Mathematics, Prof., Buryat State university, Ulan-Ude, Russian Federation) E-mail: [email protected]
TSYDYPOV Sevan Guro-Tsyrenovich (Buryat State university, Ulan-Ude, Russian Federation) E-mail: [email protected]
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. (in Russian)
