CC BY
41
2013
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Математика и механика
№ 5(25)
МЕХАНИКА
УДК 532.546
Д.О. Диль, А.М. Бубенчиков
ДВУХФАЗНАЯ ФИЛЬТРАЦИЯ В ТРУБЕ, ЗАПОЛНЕННОЙ ПОРИСТЫМ МАТЕРИАЛОМ
На основе анализа современного состояния работ по фильтрации составлена математическая модель двухфазной фильтрации в анизотропной пористой среде. Решение задачи строится с использованием метода контрольного объема и итерационной технологии Ньютона. В качестве примера рассмотрено одномерное развивающееся течение газожидкостной среды в однородном пористом материале.
Ключевые слова: пористый объем, двухфазная среда, влагонасыщенность, капиллярное давление, контрольный объем, нелинейные невязки уравнений, метод Ньютона.
Математическая модель двухфазной фильтрации в пористых пластах сформулирована на основе анализа работ [1-4]:
(1)
(2)
/ о\ V
т = т (р ) - пористость (, ¥п - объем пор);
5 - эффективная влагонасыщенность (в принципе меняется от 0 до 1); ІК - тензор проницаемости, в общем случае полно заполненный;
ко = 1 - 52) - относительная проницаемость для газовой фазы;
д° - вязкость газа; ц0 - источниковый член в балансе газовой фазы; Ц -источниковый член в балансе жидкости;
рр
рр = рр 0 —— , р^0 - плотность флюида при атмосферном давлении;
т = т
р
кГ = л/5 (1 -V1 - 52) - относительная проницаемость для жидкой фазы;
- вязкость флюида (цр » );
pF = pG - pc (s) ; pc (s) = p‘
s
sd
- капиллярное давление; p - пороговое
F
d P g
давление, причем p =--------, где а - параметр, характеризующий пористую сре-
а
ду. Таким образом независимыми переменными будут давление газа и влагона-сыщенность. В качестве граничных условий для них будем использовать условия Дирихле. В начальный момент времени давление газа и влагонасыщенность будем считать постоянными и равными значениям на выходе.
Одномерное приближение
Предельно упростим ситуацию. Примем, что источники отсутствуют, то есть
G F ~ ■* д
q = q = 0. Будем рассматривать одномерное движение, в этом случае V = i —,
dx
div = —. Кроме того, примем, что движение происходит в горизонтальном на-
dx
правлении: (g)x = 0 . В этом случае уравнения (І), (2) можно переписать следующим образом:
dPGm (І - s) д
dt dx
dPF ms = d
KkG
f dpG 11
dt
dx
V Ц
KkF
dx
V ил //
F
V
dp
dx
(3)
(4)
Здесь К - скалярный коэффициент проницаемости.
Плотность флюида будем считать неизменной, а плотность газа - линейной функцией давления:
G,0
PG =P-0“ PG.
В качестве независимых переменных будем рассматривать давление газа и влагонасыщенность. Поэтому давление флюида будем определять по формуле
PF = PG - pc (s ),
где pc (s) = p‘
л/і - s2
. Эта зависимость представлена на рис. І.
В результате получим два уравнения следующего вида:
dpGm (І - s) = d f Kk<G (s) G dp
dt
dx
dms df KkF (s) d(P° -Pc
dt dx
dx
(s ))
dx
(5)
s
Рис. 1. Характер зависимости капиллярного давления от насыщенности
Метод контрольного объема
Для построения дискретных аналогов исходных дифференциальных уравнений разобьем область решения на контрольные объемы вдоль оси х и проинтегрируем эти уравнения, заменяя объемные интегралы поверхностными. Поскольку стенки трубы непроницаемы, при интегрировании по поверхности рассматриваются только две ее составляющие - два круга, через которые осуществляется течение. Рассматривая значения давления газа и влагонасыщенности в центрах контрольных объемов, а также учитывая все предполагаемые выше упрощения исходных дифференциальных уравнений для случая одномерного приближения, будем иметь разностные уравнения следующего вида:
>>((|-5)> а-((1 -5>> = Й{Г'• еог-■ ей;,*1. (7)
И+1 - п
= еш*1 - 0^+1. ()
0,п+1 _0,п+1
V С ьО {7,п+1 \ ^0,п+1 1
К ■ Ьс ■ кг (5] ) Р/ - Р-,
где е°п1 =-с '0^ -, ( = 1ей, г1мИ,);
ц0 А/
,П+1
к■ Ьс ■ к,р (•?;+>) р0-п+1 -Рс ()-Р0,п+1 + Р (¿Т1)
А/ ’
- объем /-го контрольного объема, Ьс - площадь сечений трубы, через которые осуществляется течение, ро’п+1 = (р0 п1 + р0,п+1)/2 - осредненное значение дав-
ления газа на/-й поверхности контрольного объема, pi ,п+ - давление газа в центре /-го контрольного объема, §П+1 - значение насыщенности на /-й поверхности,
взятое против потока.
Для решения полученной системы 2Ы дискретных уравнений (Ы - число контрольных объемов) воспользуемся технологией Ньютона, используемой для решения систем нелинейных уравнений. Выпишем нелинейные невязки уравнений для 1-го приближения к величине, изменяемой на временном шаге п+1 в /-м контрольном объеме:
right
4, = v,m[pG(1 -s\)-pG,n(1 -s,n)]-Ai X pG’lQTl;
j=lefi
(9)
right
<, = Vm[sl-s,n]-A X ef' .
j=lefi
(10)
Определяя производные от невязок по значениям насыщенности и давления газа, запишем для каждого контрольного объема систему линейных алгебраических уравнений из двух уравнений для двух неизвестных приращений рассматриваемых переменных:
,,
dpl dR
Л
F ,i
dR
F ,i
ds,
^p\ ^ v8si /
f-R0?
- R
(11)
Благодаря тому, что для каждой ячейки записывается независимая система уравнений, легко решаемая с помощью правила Крамера, вычислительный алгоритм становится достаточно простым. Недостатком данного метода являются некоторые ограничения для шага по времени.
Вычисления проводились при следующих значениях определяющих параметров:
Ь = 6 м - длина трубы;
рр = 1000 кг/м3; р0,0 = 1,27 кг/м3;
К = 9,869233 • 10-13 м2;
^ = 8,9 • 10-4 Пас; = 1,78 • 10-5 Пас.
Далее будем рассматривать движение газо-жидкостной среды в трубе, заполненной пористым материалом, например песком (т = 0,4).
s = 0,8 pBGx = 1,2 -105 Па
s = 0,3
105 Па
Рис. 2. Одномерная область двухфазной фильтрации
Тестирование алгоритма при постоянной насыщенности
В случае постоянной насыщенности уравнение для давления газа принимает вид квазилинейного уравнения теплопроводности:
др?
дґ
Кк? (5) д
(
(1 - 5) дх
др
о\
дх
(12)
На рис. 3 представлены результаты решения этого уравнения рассматривае-
мым численным методом. В случае установившегося течения, когда
др
дґ
= о, ре-
шение легко находится аналитически. На графике аналитическое решение практически совпадает с численным, что говорит о корректности применения рассматриваемого численного метода для решения данных уравнений.
Рис. 3. Распределения давления газа по длине трубы в моменты времени 36, 108, 180, 360 и 36 000 с
В рассматриваемом примере профиль давления газа, изначально имеющий вид «левого нижнего уголка», постепенно деформируется в почти линейное распределение, что характерно для решения уравнения теплопроводности.
Численный эксперимент
Для анализа аппроксимационной сходимости были проведены вычисления с различными параметрами сетки, а именно шага вдоль оси х. На рис. 4 можно наблюдать характер сходимости решения разностной задачи. Заметной особенностью более точного решения является большая кривизна. В точках перегиба решения совпадают.
Рис. 4. Распределение насыщенности в моменты времени 10, 100, 360 часов для различных шагов вдоль оси х (1 - 0,5 м, 2 - 0,2 м, 3 - 0,1 м)
При указанных выше параметрах проводился численный эксперимент, отслеживающий динамику развития течения в трубе, заполненной пористым материалом. Шаг сетки вдоль оси х равнялся 0.1 м. В этом случае рассматриваемая область решения была покрыта 60-ю контрольными объёмами. Были получены распределения давления флюида и газа, а также распределение насыщенности флюидом пористого пространства с течением времени, представленное на рис. 5. Полученные результаты позволяют анализировать течение и характер заполнения пористого пространства.
Рис. 5. Распределение насыщенности по объему пористого пространства в моменты времени 2, 10, 30, 60, 100, 150, 220, 360 ч
Как видно из рис. 5, наблюдается волновой механизм распространения насыщенности по пористому пространству трубы. Причем зона изменения насыщенности с течением времени увеличивается. Характерная выпукло-вогнутая форма
ее распределения со временем становится более простой и к моменту, определяющему установившееся состояние течения, кривая 5 = 5(х) имеет вполне определенный знак кривизны.
В таблице приводятся данные итерационной сходимости при различных шагах по времени. Точность расчётов - 10-2. Для давления газа столь высокая точность является излишней и приводит к значительному увеличению итераций. Однако эта точность необходима для значений влагонасыщенности, которые колеблятся в пределах от 0 до 1.
Данные итерационной сходимости при расчёте течения с общим временем 15 дней
Шаг по времени, мин (число шагов) Общее число итераций Среднее число итераций на шаг Время выполнения программы, с
6 (3600) 590 534 768 315
60 (360) 276 477 164 142
ЛИТЕРАТУРА
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 Genuchten M.Th. A Closed-form equation for predicting the hydraulic conductivity of unsaturated soils // Soil Sci. Soc. Am. J. l980. V. 44. P. 892-898.
3. Schaap M.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. V. 5. P. 27-34.
4. Никитин К.Д. Метод конечных объемов для задачи конвекции-диффузии и моделей двухфазных течений: дис. ... канд. физ.-мат. наук. М., 2010. 105 с.
Статья поступила 30.03.2013 г.
Dil’ D.O., Bubenchikov A.M. TWO-PHASE FILTRATION IN A PIPE FILLED WITH A POROUS MATERIAL. Based on the analysis of the present-day state of works on the seepage theory, a mathematical model of two-phase flow in anisotropic porous medium is established. A numerical solution using the control volume method and Newton's iteration algorithm is proposed. As an example, one-dimensional developing motion of a gas-fluid medium in a homogeneous porous material is considered.
Keywords: porous volume, two-phase medium, water saturation, capillary pressure, control volume, nonlinear residuals, Newton's method.
DIL’Denis Olegovich (Tomsk State University)
E-mail: [email protected]
BUBENCHIKOV AlekseyМikhailovich (Tomsk State University)
E-mail: [email protected]
