CC BY
47
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2013
Математика и механика
№ 6(26)
УДК 629.78.001
Д.О. Диль, А.М. Бубенчиков
ДВУХФАЗНАЯ ФИЛЬТРАЦИЯ В АНИЗОТРОПНОМ ПРОСТРАНСТВЕ
Представлено решение задачи двухфазной фильтрации в анизотропной по проницаемости среде, полученное методом контрольного объёма. Для аппроксимации потоков на гранях контрольного объёма использовалась нелинейная двухточечная схема. Решение дискретных уравнений осуществлялось итерационным методом Ньютона.
Ключевые слова: двухфазное течение, метод контрольного объёма, нелинейная двухточечная аппроксимация, метод Ньютона.
Решение задач двухфазной фильтрации имеет большое значение для повышения эффективности добычи полезных ископаемых [1-3]. Для математического моделирования процесса двухфазной фильтрации воспользуемся моделью, в основу которой положен закон Дарси. Два основных уравнения для независимых давления газа р° и насыщенности флюида 5 представляют собой законы сохранения массы газа (1) и флюида (2) и имеют следующий вид:
тензор проницаемости К - зависящим только от пространственных координат, а
мое для того, чтобы несмачивающая жидкая фаза начала вытеснять смачивающую фазу. В данной работе будем считать его равным 10 кПа.
В качестве расчётной области рассмотрим квадрат со стороной Ь с непроницаемыми стенками, входом на верхней стенке и выходом на правой стенке (см рис. на с. 76, Ь = 12 м). Обозначим множество всех граничных точек через Г,
(1)
(2)
Пористость т, вязкость и плотность р^ флюида, а также вязкость газа будем считать неизменными, плотность газа р° - линейной функцией давления,
относительные проницаемости для каждой из фаз (5) и к.^ (5), а также капил-
лярное давление рс (5) - зависящими от насыщенности по формулам Ван-Генухтена - Муалема [4]:
РС (5) = ^ , РР = Р° -РС (5)• (4)
5
(4)
Здесь ра - пороговое давление - минимальное капиллярное давление, необходи-
множество точек входа - через Геп = {(х, у), Ь < х < Ь2, у = Ь}; выхода - через
Гех = {(х,у), х = Ь, Ь3 < у < Ь4} (Ь = Ь3 = 2 м, Ь2 = Ь4 = 8 м). В качестве гранич-
ных условий на входе и выходе будем использовать условия Дирихле:
р° = pen, 5 = ^ при (х у) 6 Геп; (5)
р° = рех.5 = ^ при (x, у) 6 Гех.
На остальной границе поток через стенку будем считать нулевым, что математически может быть выражено как равенство нулю производной по нормали к стенке:
д5
= 0,~^ = а при (x, у) 6 Г\(Геп и Гех). (6)
дп дп
При численном расчёте методом контрольного объёма такие граничные условия реализуются простым приравниванием к нулю потока через грань контрольного объёма, лежащую на границе. В начальный момент времени давление газа и влагонасыщенность внутри расчётной области примем постоянными и равными значениям на выходе.
Исходные уравнения (1), (2) как система, имеющая замыкания (3), (4), глобально имеет параболический тип. При установившихся режимах фильтрации задача может быть решена без эволюционных членов, то есть как глобально эллиптическая. Если в какой-либо точке рассматриваемой области 5 ^ 1, то порядок стремления к нулю эволюционного члена будет (1 - 5), а диффузионного (1 - 5)32. Это означает, что при таком стремлении относительный вклад диффузионной части будет исчезающее малым. В результате асимптотически вместо (1) получим
дР°т = 0 д/
С математической точки зрения это проявление гиперболических свойств системы уравнений в рассматриваемой точке. С физических же позиций последнее означает, что плотность газа, а в рассматриваемой постановке и его давление, перестают меняться в данной точке пространства и времени.
То же самое можно было бы сказать и в отношении уравнения (2), исследовав его состояние при 5 ^ 0 , если бы представленные замыкания (3), (4) не требовали бы существования предельной влагонасыщенности 5пр > 0 , определенной величиной порогового давления ра.
Таким образом, для определенных краевых условий в расширенной области определения задачи (физическая область, дополненная координатой времени) могут существовать точки, линии, поверхности и т.д., в которых система исходных уравнений меняет тип. Если такие подобласти существуют, то предсказать их расположение до решения задачи мы не можем.
Для многих специализированных методов смена типа интегрируемой системы при продвижении от точки к точке расширенной области представляет собой неразрешимую проблему. Используемый в настоящей работе метод контрольного объема не является критичным к указанным свойствам системы.
Аппроксимация потоков на гранях контрольного объёма
Для построения дискретного аналога определяющих уравнений использовался метод контрольного объёма с нелинейной двухточечной аппроксимацией потока на многогранных ячейках, предложенной в [3]. Рассмотрим эту методику.
При интегрировании основных уравнений по контрольному объёму объёмный интеграл от дивергенции переходит в поверхностный, который имеет физический смысл суммы потоков флюида или газа через грани контрольного объёма. Для каждой грани с номером/такой поток будет иметь следующий вид (а = О, Е):
Здесь через п/ обозначена внешняя нормаль к грани/, а через Б/ - её площадь. Далее рассматривается так называемый вектор конормали, равный произведению тензора проницаемости на вектор нормали к текущей грани КП/ (обозначим его
через С/ ). Скалярное произведение конормали и градиента давления можно выразить через производную по её направлению:
Чтобы осуществить дискретизацию - выразить потоки через значения давлений в центрах контрольных объёмов - для каждой конормали вводится базис, состоящий из векторов, направленных из центра текущего объёма к центрам соседних. При этом выбираются те векторы, для которых конормаль будет лежать внутри трёхгранного угла, образованного ими (рис. 1).
у
(7)
(8)
-Ч
Рис. 1. Выбор базиса для конормали
є/ , і
Используя разложение конормали в выбранном базисе (обозначим его Є/ і , і = 1,3), перепишем выражение для производной по направлению конормали:
дт / = &,.,£*. (9)
дс/ 1=1 деГ,1
где а у і - коэффициенты разложения конормали по нормированному базису
/И,1=1,3.
Заменяя частные производные разностными, будем иметь следующее выражение для скалярного произведения конормали и градиента давления:
— - 3 ра - ра
С/^р°=£ а/,1-^. (10)
і=1 Iе/, і
Обозначим полученное выражение через Ну+ и, используя в качестве основного соседний контрольный объём, также содержащий грань / , повторим ту же процедуру. В итоге будем иметь две различных аппроксимации одного и того же
выражения. Избавимся от слагаемых, в которые входят значения давления в цен-
трах контрольных объёмов, не содержащих грань / Это необходимо, чтобы итоговая аппроксимация была двухточечной. Для этого вводятся веса (п+, П- ), сумма которых равняется единице:
с/-Ура Н/ +п-Ну- = ОуТр- -Оу+ р+, (11)
и выражение для потока газа или флюида через грань контрольного объёма / приобретает следующий вид:
К (5 г)
Є“ =J-/ Ра (О/-ра - Оу+ ра)Бу, (12)
где 5/ - значение влагонасыщенности на поверхности / взятое против потока,
ра - значение давления с внешней стороны грани / а р+" - с внутренней. Их
среднее арифметическое используется для вычисления плотности газа, входящей в формулу для потока.
Метод Ньютона
Интегрируя определяющие уравнения по контрольному объёму и заменяя частные производные по времени разностными, будем иметь 2Ы дискретных уравнений (Ы - число контрольных объёмов):
^т!Р^;211-51])П+^:(£!(£О)Ы): = £е°, (13)
/=1
и+1 : И к __
Угтрр г г =ХЯ, і = 1,N. (14)
А' /=1
Для решения полученной системы дискретных уравнений воспользуемся методом Ньютона решения систем нелинейных уравнений. Выпишем нелинейные невязки уравнений для 1-го приближения к величине, изменяемой на временном шаге п+1 в г-м контрольном объеме:
Дифференцируя невязки по влагонасыщенности и давлению газа, запишем систему линейных алгебраических уравнений с разреженной матрицей Якоби. Решая данную систему, будем определять искомые приращения на каждой итерации, пока процесс не сойдётся с заданной точностью.
Вычислительный алгоритм тестировался для одномерных и двумерных задач при постоянной влагонасыщенности. Полученные результаты обнаружили совпадение распределений давления в одномерном случае с простейшим аналитическим решением, а в двумерном - с численным решением уравнения однофазной фильтрации газа, полученным с помощью метода Гаусса - Зейделя. На рис. 2 приведены распределения давления газа для двухфазного течения при постоянной влагонасыщенности и однофазного течения газа через пористое пространство с тем же значением проницаемости через 10 ч. Распределения практически идентичны.
0123456789 10 11 12
Рис. 2. Распределения давления (в Па) для однофазного течения (-)
и двухфазного течения с постоянной влагонасыщенностью (- - -)
в области 12x12 м. Расчеты отвечают установившемуся режиму фильтрации
(15)
Г=1
к
(16)
г=1
■Р
11-
10-
2-
/
I
\
/
1-
Ч
о-
Для анализа аппроксимационной сходимости были проведены расчёты на равномерной сетке с различными шагами по пространству. В качестве параметра, с помощью которого проводился анализ, была использована средняя скорость фильтрации как сумма фильтрационных потоков через каждую из граней всех контрольных объёмов, делённая на количество граней. Каждый из таких потоков является объёмным расходом фильтрующейся среды и является ключевым показателем динамики фильтрационного течения. На рис. 3 приведены зависимости средней скорости фильтрации флюида и газа от времени для сеток с шагом 1, 0,5 и 0,25 м.
Ш, мкм/с
Ш, мкм/с
Рис. 3. Зависимости средней скорости для флюида (вверху) и газа (внизу) от времени при расчётах на сетках с шагом 1 (—), 0,5 (- - -) и 0,25 м (сплошная линия)
Численный эксперимент
В каменноугольных породах проницаемость вдоль трещин может в 3-10 раз превышать проницаемость в других направлениях. В данной работе проводился расчёт течения смеси газа и флюида через область, содержащую пласт с анизотропной проницаемостью. На рис. 4 схематично изображён пласт с поперечными трещинами (перпендикулярными границам пласта), проницаемость вдоль которых в 10 раз превышала проницаемость в продольном направлении и вне пласта; вход (вверху) и выход (справа), на которых задавались граничные условия Дирихле. На левой границе ордината точек, лежащих на пласте, удовлетворяет условию 4,5 <у < 8,5, на правой - 7,5 <у < 11,5. Угол наклона пласта по отношению к оси х: атС^(1/3), при 0 < х < 8, и ат^(1/5), при 8 < х < 12.
Рис. 4. Расположение анизотропной части среды
Численный эксперимент проводился при следующих значениях определяющих параметров:
I = 12 м; т = 0,05; рр = 1000 кг/м3; р0,0 = 1,07 кг/м3;
К = 9,869233 • 10-13 м2 = 1 Дарси - абсолютная проницаемость изотропных частей среды; проницаемость вдоль трещин бралась 10 Дарси;
= 8,9 • 10-4 Пас; = 1,78 • 10-5 Пас - вязкости сред. Сетка бралась равно-
мерной с шагом 0,5 м. Шаг по времени - 10 мин.
Давление газа и влагонасыщенность на входе задавались постоянными и равными 120 кПа и 0,9 соответственно; на выходе, а также в начальный момент времени во всей расчётной области давление газа принималось равным 10 кПа, вла-гонасыщенность - остаточной насыщенности 0,2.
На рис. 5 представлена динамика развития течения в разнородной среде, а именно распределение влагонасыщенности в разные моменты времени.
Рис. 5. Распределения влагонасыщенности в моменты времени 150, 250, 450, 525, 700 и 1000 часов
Полученные результаты позволяют анализировать течение и характер заполнения пористого пространства, содержащего трещины. Так вдоль трещин скорость течения намного превосходит скорость течения в изотропном пространстве, что обуславливает скопление жидкости на границе между породами. Кроме того, метод контрольного объёма позволяет легко оценить количество жидкости или газа, проходящих через различные сечения, поскольку потоки, вычисляемые при решении уравнений имеют физический смысл массовых потоков через грани выделенного объёма.
ЛИТЕРАТУРА
1. Lama R.D. A model for determination of methane flow parameters in coal from desorption tests / R. D. Lama, V. U. Nguyen // Materials of Int. Symposium on the Applications of Computers and Mathematics in the Mineral Industries, Johannesburg, SAIMM. 1987. V. 1. P. 275-282.
2. 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.
3. Никитин К.Д. Метод конечных объемов для задачи конвекции-диффузии и моделей двухфазных течений: дис. ... канд. физ.-мат. наук. М., 2010. 105 с.
4. Van GenuchtenM.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated Soils // Soil Sci. Soc. Am. J. 1980. V. 44. P. 892-898.
Статья поступила 19.09.2013 г.
Dill D.O., Bubenchikov A.M. TWO-PHASE SEEPAGE IN AN ANISOTROPIC MEDIUM. A solution of the two-phase seepage problem by the FVM is presented. To approximate flows through finite volume faces, the nonlinear two-point scheme was used. The discrete equations were solved by the Newton method.
Keywords: two-phase flow, finite volume method, nonlinear two-point approximation, Newton's method.
Dill Denis Olegovich (Tomsk State University)
E-mail: gradpower@list.ru
Bubenchikov Aleksey Mikhailovich (Tomsk State University)
E-mail: alexy121@mail.ru
