Для исследования характера многофазных фильтрационных течений в неоднородных пластах широко применяются методы математического моделирования. Решение практических задач многофазной фильтрации возможно только с применением численных методов. Проблемы построения эффективных вычислительных алгоритмов для данного достаточно сложного класса нелинейных задач хорошо известны и далеко не преодолены [2], [3], [6].
Рассмотрим совместное течение двух несжимаемых жидкостей в среде с постоянной пористостью и проницаемостью. Полагая в качестве неизвестных s –насыщенность порового пространства первой фазой и P – давление во второй фазе, система двух уравнений относительно определяемых полей насыщенности и давления в безразмерном виде имеет следующий вид
Система (1), (2) дополняется стандартными граничными условиями для давления и граничными и начальным условиями для насыщенности, определяющимися конкретной постановкой задачи.
Подходы, лежащие в основе предлагаемого численного алгоритма, рассмотрим на примере решения уравнения вида
, (3)
содержащего все дифференциальные операторы, входящие в систему (1), (2).
Пусть уравнение (3) решается в области 
в некоторой системе координат (ξ, η). Допускается, что границы областей ξ = z0(η,τ) и ξ = zN(η,τ) могут перемещаться с течение времени.
Для дискретизации расчетной области используется полукриволинейная подвижная сетка 
, составленная совокупностью узлов, соответствующих точкам пересечения семейства фиксированных координатных линий η = ηk и семейства подвижных линий вида ξ = zn(η,τ) . Для последних должно выполняться требование 
для любого n и для каждой пары значений (η,τ). Фрагмент такой расчетной сетки представлен на рис. 1.
На каждый момент времени с любым узлом (n, k) связана тройка величин (ξnk(τ), ηk, σnk(τ)), где первая из координат ξnk(τ) = zn(ηk,τ) и значение решения в этом узле σnk(τ) = σ(ξn(τ), ηk) могут изменяться во времени. Положение узла на координатной линии h = hk может определяться заданием одной из величин: координаты ξnk(τ) или значения узловой величины σnk(τ). Таким образом, расчетная сетка может включать узлы двух типов: перемещающиеся вдоль данной координатной линии по заданному закону и прослеживающие перемещение точки с заданным изменением рассчитываемой узловой характеристики. В частном случае величины ξnk и σnk могут быть фиксированными.
Рисунок 1. Фрагмент расчетной сетки
В рамках сделанных допущений можно определить расчетные сетки самой разнообразной структуры, сочетающей узлы двух описанных типов.
Построение системы уравнений, представляющих дискретный аналог исходного дифференциального уравнения (3), основывается на интегральных тождествах, справедливых для окрестности каждого узла расчетной сетки.
Пусть (n, k) – произвольный внутренний узел. Его сеточная окрестность D(n, k) представляется объединением элементарных подобластей
Интегральное тождество для узла получается интегрированием уравнения (3) по каждой из подобластей, составляющих его окрестность, со своими весовыми функциями и последующим суммированием результатов интегрирования.
Веса интегрирования определяются с помощью набора функций следующего вида:
где нижний индекс совпадает с индексом узла, а верхний определяет ориентацию подобласти относительно данного узла. Для подобласти Dnk узла (n, k) весовая функция есть произведение 
, для подобласти Dn-1,k этого же узла – произведение 
и т.д.
Далее будем считать, что координатная система (ξ, η) является ортогональной, так что развернутое представление уравнения (3) имеет вид
Результат интегрирования уравнения (4) по подобласти Dnk узла (n, k) приводится к равенству
Для остальных подобластей окрестности узла получаются равенства, аналогичные (5). Заметим, что правые части равенств (5) представлены потоковыми членами по общим границам подобластей и при суммировании взаимно уничтожаются. В итоге интегральное тождество по окрестности произвольного внутреннего узла будет иметь вид
. (6)
Для узлов, лежащих на границе расчетной области, тождества интегрального баланса получается аналогично за исключением того, что окрестность граничного узла состоит из одной или двух подобластей, а члены, выражающие потоки на внешних границах области заменяются заданными функциями из граничных условий.
Система численных уравнений относительно неизвестных узловых величин выводится на основе тождеств (6) последовательным раскрытием интегралов по обеим координатам. Для неизвестного распределения σ(ξ, η) в пределах каждой подобласти используется линейная аппроксимация по координате x вдоль каждого фиксированного значения координаты η:
При получении аппроксимирующих уравнений входящие в (4) коэффициенты в пределах каждой элементарной подобласти Dnk заменяются их осредненными по подобласти значениями.
Производные от узловых величин по времени аппроксимируются неявными разностными отношениями, гарантирующими безусловную устойчивость вычислительной схемы от величины временного шага, что является принципиально важным при использовании расчетных сеток с большой пространственной неравномерностью. В конечном итоге расчет одного временного шага приводит к девятиточечной системе алгебраических уравнений относительно неизвестных величин в точках расчетной сетки, входящих в единую узловую окрестность. В случае явного задания расположения узлов xnk получаемая система уравнений линейна. В общем случае, когда сетка включает узлы с заданными значениями определяемого решения snk, система является нелинейной относительно координат xnk. Для ее решения используется линеаризация по методу Ньютона, приводящая на каждой итерации к линейным девятиточечным уравнениям с блочно-трехдиагональной матрицей. К последним применяется матричная прогонка [5]. Отметим, что численная реализация алгоритма не зависит от конкретной структуры расчетной сетки. При этом могут использоваться различные координатные системы, а также разбиение расчетной области на конечное число подобластей, геометрия которых для подходящей системы координат не противоречит оговоренным соглашениям.
В качестве примера рассматривается процесс вытеснения нефти водой в однородном пласте на элементе пятиточечной системы скважин. С учетом симметрии задачи исследуемая область представляет единичный квадрат с круговыми вырезами радиуса r0 с центрами в двух противоположных углах, соответствующих контурам нагнетающей и добывающей скважин:
.
Постановка задачи включает уравнения двухфазной фильтраци с граничными условиями:
Для численного решения использовалась координатная система ξ = 1nρ, η — φ где ρ и φ – радиус и угол полярных координат. Расчет поля давлений и поля насыщенности производился на отдельных сетках.
Для учета больших градиентов давления в прискважинных зонах и получения пространственной дискретизации, близкой к конфигурации изобар и линий тока удобнее разбить исследуемую область на две подобласти D1 и D2 диагональю проходящей через вершины, в которых нет скважин, и перейти для каждой подобласти к своим полярным координатам. В каждой из двух подобластей полярный полюс располагаем в той вершине, где находится ось скважины. Таким образом, 
, где 
. В каждой области задавалась равномерная по координатам ξ и η расчетная сетка, дающая автоматическое сгущение узлов в узких прискважинных зонах, где поле давлений имеет большие градиенты
Поле насыщенности рассчитывалось в системе координат (ξ, η) с полюсом в центре нагнетающей скважины на подвижной сетке, соответствующей совокупности линий ξ = zn(η,τ) с одинаковыми значениями насыщенности sn в интервале snÎ[0; 1]. Таким образом, расчет распределения σ(ξ, η, τ) сводился к прослеживанию перемещения вдоль прямых h = hk точек линий изонасыщенности с заданными значениями 
. Для интервала изменения величин σ и η задавалось равномерное разбиение. Положение изолинии ξ = z0(η,τ) с насыщенностью σ0=1 фиксировано и совпадает с контуром скважины. Изолиния ξ = zN(η,τ) с насыщенностью σN=0 представляет подвижную границу области двухфазного течения с соответствующим граничным условием
.
Связанная система уравнений (5), (6) является внешне нелинейной относительно определяемых полей давления и насыщенности. Такой же нелинейностью характеризуется и система численных уравнений при использовании полностью неявной аппроксимации по времени. Поэтому расчет каждого шага по времени требует внешних итераций. Последовательность действий для перехода на новый временной слой включала следующие этапы:
- Для заданного распределения насыщенности из дискретного аналога уравнения (2) находится поле давлений. Уравнение (2) является линейным относительно давлений и решается на фиксированной сетке, поэтому система численных уравнений также линейна.
- Исходя из найденного поля давлений из дискретного аналога уравнения (1) определяется поле насыщенности. Так как для нелинейного уравнения (1) используется подвижная сетка с фиксированными значениями насыщенности, то система численных уравнений нелинейна относительно координат узлов и решается итерациями по методу Ньютона [5].
Оба этапа итерируются до достижения сходимости внешних итераций по насыщенности. Так как каждое из определяемых полей рассчитывается на разных пространственных сетках, то перед каждым этапом производится интерполяция давления и насыщенности с одной сетки на другую.
Предлагаемый численный алгоритм был протестирован на рассмотренной задаче с использованием последовательности сгущающихся сеток для различных значений безразмерных величин 
. Тестирование показало, что хорошее пространственное разрешение особенностей рассчитываемых полей достигается на сетках с небольшим числом узлов. В качестве примера на рис. 2 представлены поля давления и насыщенности, рассчитанные с помощью изложенного алгоритма на сетках 
.
Рисунок 2. a) Распределение давления b) Распределение насыщенности
Численное моделирование задачи двухфазной фильтрации позволяет использовать созданный в ней математический, алгоритмический и программный продукт для решения задач прогнозирования многофазных фильтрационных течений, в частности в пластах сложной структуры. Результаты численного исследования дают ряд новых выводов об особенностях вытеснения нефти из слоисто-неоднородной среды, которые могут найти практическое применение при совершенствовании технологий разработки нефтяных пластов с послойной неоднородностью.
Список литературы:
- Влияние межпластовых перетоков и капиллярных сил на процесс вытеснения нефти в слоисто-неоднородном пласте [Текст] / Ю. А. Сигунов, Г. Р. Усманова // Известия РАН. Механика жидкости и газа. — 2007. — N 6. — С. 85-92 : 4 рис. — Библиогр.: с. 92 (7 назв. )
- Математическое моделирование пластовых систем : пер. с англ. / Х. Азиз, Э. Сеттари . – 2-е изд., стереотип . – М. : Ин-т компьют. исслед., 2004 . – 416 с. – (Современные нефтегазовые технологии) . — ISBN 5-939723-55-1
- Коновалов А.Н. Задачи фильтрации многофазной несжимаемой жидкости. – Новосибирск: Наука, 1988. — 166 с. — ISBN 5-02-028569 -2.
- Сигунов Ю.А. Методы решения классической задачи Стефана – Сургут: РИО Сургутского государственного педагогического университета, 2009. – 140с.).
- Самарский А.А., Гулин А.В. Численные методы: Учебное пособие для вузов. – М.: Наука. Гл. ред. физ.-мат. лит., 1989. – 432 с.
- Флетчер К. Вычислительные методы в динамике жидкостей. В 2–х томах. – М.: Мир. 1991. – 504 c. – ISBN: 5-03-001880-8, 5-03-001881-6 (т.1), 5-03-001881-4 (т.2)[schema type=»book» name=»МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧ ДВУХФАЗНОЙ ФИЛЬТРАЦИИ НА СВЯЗАННЫХ ПОДВИЖНЫХ СЕТКАХ» description=»В данной статье рассматривается решение двумерных задач двухфазной фильтрации на подвижных сетках, пространственная структура которых определяется самими значениями рассчитываемого решения. Для задач нелинейной теплопроводности применение таких сеток, представляющих собой совокупность подвижных изотерм, дает эффективный способ выделения и прослеживания таких особенностей, как движущиеся границы фазовых переходов и бегущие температурные волны. В статье дан алгоритм численного метода к решению двумерных задач фильтрации на подвижных сетках, привязанных к значениям рассчитываемого поля температур, с обобщением на класс задач многофазной фильтрации.» author=»Прозорова Гульшат Ринатовна» publisher=»БАСАРАНОВИЧ ЕКАТЕРИНА» pubdate=»2016-12-26″ edition=»euroasia-science.ru_26-27.02.2016_2(23)» ebook=»yes» ]