Применение стандартного многосеточного метода к задачам с условием Дирихле на внутренних границах. Приложение для моделирования нефтяных месторождений Текст научной статьи по специальности «Математика»

УДК 519.633+004.021 +532.546

И. Т. Усманов, Д. В. Шевченко

ПРИМЕНЕНИЕ СТАНДАРТНОГО МНОГОСЕТОЧНОГО МЕТОДА К ЗАДАЧАМ С УСЛОВИЕМ ДИРИХЛЕ НА ВНУТРЕННИХ ГРАНИЦАХ. ПРИЛОЖЕНИЕ ДЛЯ МОДЕЛИРОВАНИЯ НЕФТЯНЫХ МЕСТОРОЖДЕНИЙ

Ключевые слова: многосеточный метод, скорость сходимости, моделирование нефтяных месторождений.

Предложеныпреобразования разностных уравнений для задания условия Дирихле на внутренних границах при решении задачи геометрическим многосеточным методом. В приложении к моделированию нефтяных месторождений преобразования позволяют задавать значения забойных давлений скважин во внутренних узлах ко-нечноразностной сетки и распределения давления во внутренних границах.

Keywords: multigrid method, convergence rate, modeling of oilfields.

Proposed conversion of differential equations to specify Dirichlet conditions at internal borders in solving the problem of the standard multigrid method. In application to the modeling of oil fields conversions allow you to specify the values of the bottom-hole pressure of the wells in the inner nodes of the finite-difference grid and pressure distribution in the inner borders.

Введение

Спецификой нефтяного месторождения являются его большие геометрические размеры и значительное время эксплуатации, и, вместе с тем, наличие скважин, размер которых несоизмеримо мал по сравнению с размерами пласта. Поэтому гидродинамический расчет для него требует достаточной детализации и большого объема вычислений. Одним из лучших вариантов организации численного расчета в этом случае является применение многосеточной технологии [1]. Использование нескольких уровней детализации позволяет избежать конфликтов между компонентами решения, имеющими разный масштаб [2].

Наличие источников или стоков (скважин) приводит к особенностям в решении уравнения пье-зопроводности. В работе [3] описан способ подправки коэффициентов задачи, при использовании которого разностная схема строится без аппроксимации контуров скважин. При этом мощность источника (стока) хорошо согласуется со значением функции давления в узле, которое интерпретируется как «забойное давление». В результатезабойные давления записываются в узлы сетки, коэффициен-тызадачи подправляются. Дополнительные особенности скважин как правило моделируются с использованием понятия скин-фактора или эффективного радиуса скважины [4,5].

Задачи с внутренними и криволинейными границами часто приводят к решению задачи с условием Дирихле. Видразностного уравнения в узле конечноразностной сетки с заданным значением искомой функции отличается от вида уравнений в других узлах. Решитьсистему разностных уравнений напрямую с использованием стандартного геометрического многосеточного метода в данном случае нельзя. Связано это с проблемой переноса сеточного оператора на грубые сетки, а получившаяся матрица коэффициентов системы разностных уравнений на мелкой сеткеявляется несимметричной.

В данной работе предлагается способ преобразования разностных соотношений для задания

известных значений искомой функции. Приведено соответствующее преобразование матрицы системы. Эффективность предложенного подхода продемонстрирована на примере задачи определения давления в нефтяном пласте.

Постановка задачи

Рассмотрим задачу изотермической фильтрации двухфазной жидкости в природном пласте, вскрытом системой добывающих и нагнетательных скважин [4]. Фильтрация происходит по обобщенному линейному закону Дарси. В двумерной постановке (при осреднении по толщине пласта [6]) решение задачи сводится к определению приведенного пластового давления Р (х, у, /), подчиняющегося

следующему уравнению в частных производных

дР

(а вгаё Р )"3— = 0 (1)

где а = а( х, у, /) - гидропроводность пласта, 3 = х, у, /) - упругоемкость пласта и жидкостей.

Уравнение справедливо в многосвязной области, ограниченной внешним контуром питания Г и внутренними границами: скважинами ук,

(к = 1..Ы), а также внутренними контурами питания Гг, (I = 1,...,М ). Последние образуются, например, за счет слияниянефтяного пласта с мощным водоносным пластом.

На внешнем и внутреннем контурах питания известно распределение давления

Р = р( х, у), (х, у )еГиГг, (( = 1,..., М). (2) Для скважин задано забойное давление либо

дебит:

Р = Рк (г), (х,у) еГк, (к = 1,...,Мр),

= Чк ^), (к = Мр +1,...,N). (3)

Гк

На начальный момент времени известно распределение давления Р(х, у,0).

Применение многосеточного метода в случае, когда во всех скважинах задан дебит, рассмотрено в [1]. Здесь рассмотрим только случаи заданных забойных давлений (условие (3)) и заданного давления во внутренних подобластях (2). Объединение всех случаев не представляет сложности.

Построение разностной схемы

Введем сетку

аь ={(*,,У.): X = хо + гН, У = Уо + А г = О,-.,Кх> } = 0,..., Ку} с равномерным шагом Н , где (Ых +1), (Му +1) - число узлов по направлениям х и у .

Будем предполагать, что траектории скважин и границы контуров питания проходят через узлы сетки. Для точек с координатами х, ± Н/2 будем использовать обозначение г ±12, аналогично ] ± 1/2 для направления У .

Как было отмечено выше, подправка коэффициента гидропроводности [3] отождествляет давление в узле, через которое проходит скважина, с забойным давлением. Коэффициент гидропроводно-сти в полузлах сетки умножается на поправочный коэффициент. Значение поправочного коэффициента отлично от нуля лишь в полуузлах на границах ячеек скважин.

В итоге для уравнения пьезопроводности (1) консервативная разностная схема [7] будет иметь вид

1 (( " Ри ) " * 1 . (] " Рг-1,] ) +

+с

=Н1

Т

+1 (Рг, ]+1 - Ри ] У*,. -1 ((] - Ри ]-1 ) =

(4)

где Рг. - давление с прошлого временного шага.

Разностное уравнение справедливо во всех точках сетки (оН кроме граничных уН и узлов скважин с заданным забойным давлением. В этих точках давление известно:

р, , ] = ъ, ], С, ]) еп; Рг,, Л = Рк, к = l,..., ыр.

В итоге приходим к системе, в которой наряду со «стандартными» уравнениями типа (4) имеются «особенные» - в узлах скважин и граничных узлах. При стандартном применении операторов сужения [2], данная особенность не переносится на грубую сетку. Кроме того, матрица системы линейных уравнений становится несимметричной. В итоге сходимость многосеточного алгоритма существенно портится. Далее предлагается способ преобразования разностных коэффициентов задачи, приводящий к однотипности условий во всех узлах сетки и симметризации матрицы, что позволяет использовать все преимущества многосеточной технологии.

Модификация коэффициентов задачи для расчета в режиме заданных забойных давлений

Способ модификации коэффициентов, для простоты изложения, покажем на одномерном слу-

чае. Тогда разностные соотношения в узлах г и ( г +1) примут вид

*+12 (Р,+1 - Р, ) - *-12 (Р, - Р,-1 ) =

-НН-Р( - РО)

ст,-

(Рг+ 2 - Р,+1 ) - *+12 (Р,+1 - Рг ) =

32

Н2 „ ( о (+1 - Рг+1

-Р+1 (

(5)

(6)

Осуществим преобразование уравнений для задания известных значений давления. Пусть в узле г известно давление Р,. В данном узле, во-первых, должно выполняться разностное уравнение (5), во-вторых,

Рг = Рг .

(7)

Для выполнения обоих условий введем новые значения коэффициентов гидропроводно-сти,упругоемкости и значения функции на предыдущем временном слое *г+1/2, Р, Р°+1 и т.д. Здесь и

далее новые значения сеточных величин обозначены чертой сверху. Приняв в узле с известным значением

*+1/2 = 0 , *-1/2 = 0 , (8) изуравнения (5) получим

Д( Р, - Р, ) = 0. (9)

Следующим шагом является преобразование уравнений в соседних с г узлах. Сами разностные соотношения относительно неизвестных значений давления при этом должны сохранить первоначальный вид, меняются только известные коэффициенты. Тогда уравнение (6) с учетом (8) примет вид:

*,+3/2 ( Р,+ 2 - Р,+1 ) - 0 ' (Р,+1 - Р, ) =

Н^ / Р0

Рг+1 - Р0+1

Р+1 (

(10)

Соотношение (10) будет в точности совпадать с (6) при

Р+1 =Р+1 +

*,+1/2Т Н2 '

Р0+1 _ Р Р0+1 + , 2 1 Рг. Р+1 Н Рг+1

(11)

В узле (г -1) разностное уравнение имеет

аналогичный вид.

Новое значение коэффициента упругоемко-сти Р является любым отличным от нуля числом, но от величины Р зависит точность выполнения условия (7): с учетом (9) для невязки е имеем оценку! Р, -Р,\ <еР .

Для двумерного случая формулы получаются аналогичным образом.

Рассмотрим физический смысл преобразо-ваний.Наиболее наглядная интерпретация получается для двумерного случая (рис. 1).

Рис. 1 - Физический аналог модификации коэффициентов

Обнуление гидропроводностей равнозначно вводу непроницаемых стенок на границах ячейки с заданным значением давления. Потерянные потоки компенсируются в соседних ячейках изменением объема жидкости.

Для рассмотрения алгебраического смысла преобразования необходимо выписать систему разностных уравнений для одномерной задачи. Для известного значения функции в точке х. система будет иметь вид:

(,

0.

0

Л

а1, 1 - а1, 1 а1, 1+1

ап -1, п-2 ап-1, п -1

{ Р 1 {1 ^

Л

V рп-1 у V ^п-\

где

(12)

а1.1 ±1 =-а1 ±1/2 , (( * /), а>, 1 = \, а1,1 ±1 = 0, 1 = Р, ( = /).

Для межсеточных переходов необходимо, чтобы разностным соотношениям во всех точках сетки соответствовали строки с коэффициентами вида (12), условие(5) при этом также должно выполняться. Обоим требованиям удовлетворяем после умножения слева левой и правой части матричного уравнения на матрицу вида

{10.................................01

0.0 1 0........................0

0......0 1 - а1-и 0 0......0

0.........0 3 0.........0

0......0 0 - а,+и 10......0

0........................0 10 ... 0

0.

.0 1

V ................................. /

Выбор значения ¡3 был описан выше, зависит от требования к точности выполнения (5).

Преобразование разностных уравнений позволяет обрабатывать точку с заданным значением

при операциях многосеточного метода также, как и остальные. Фактически меняются коэффициенты задачи и правые части на самой мелкой сетке.

В результате, для решения системы разностных уравнений применяется стандартная геометрическая многосеточная технология.

Распределение давления на внутренних контурах питания задаются аналогично. Подобласти внутренних контуров питания при операциях сглаживания на мелкой сетке могут быть пропущены.На грубых сетках при операциях сглаживания могут быть пропущены узлы, полученные сужением узлов только данной подобласти (Рис. 2).

Рис. 2 - Работа с точками в областях выклинивания (■ -граничные точки, О - точки, не участвующие в операциях многосеточного метода)

Тестирование метода

В качестве тестов для проделанной модификации коэффициентов были проведены следующие численные эксперименты. Во-первых, проведено сравнение решений задачи при заданных дебитах и при заданных забойных давлениях. Для этого сделано следующее. Задача решена при заданных деби-тах, получено распределение давления.Далее полученное для режима заданных дебитов значение в узле скважины задавалось в качестве забойного давления, и задача решалась с данным условием. Полученные решения отличались не более чем на величину допустимой погрешности, с которой проводился расчет.

Во втором тесте задача вначале решалась при заданных дебитах. Далее полученные значения давления в произвольных узлах за исключением граничных задавались в качестве известных и задача решалась повторно. Распределения давления также отличаются не более чем на заданную погрешность расчета. Данный тест демонстрирует, что предложенное преобразование не вносит дополнительных возмущений в решение.

Анализ скорости сходимости полученного многосеточного метода показал, что после преобразования разностных соотношений коэффициент понижения невязки за итерацию не ухудшился.

Было проведено сравнение скорости сходимости многосеточного метода после преобразования разностных соотношений при разных полях коэффициента гидропроводности. Для гладкого поля со значением единица и поля с резкими градиентами гидропроводности, где значения на отдельных линиях вырождались с единицы до 1е-8, получены следующие коэффициенты понижения невязки (таблица 1).

Таблица 1 - Коэффициент понижения невязки при различных полях гидроповодности и размерностях сетки

Размер Коэффициент понижения невязки

сетки для гладкого для поля

поля с резкими

градиентами

200х200 0,1 0,3

1000х1000 0,1 0,35

Наконец, было проведено сравнение численного и точного решений. Для задачи определения давления в круговом пласте с единичной скважиной, на которой задан дебит, известно [4] аналитическое решение

Р (х, у)= ¿1п [(( - Х0 ) + (У - У0 )) . (13)

Особенностью сравнения с данным точным решением является то, что численное решение-должно отличаться от аналитического на некоторую величину, являющейся неустранимой погрешностью. Величина данной погрешности равнаразности между фундаментальным решением оператора Лапласа и главным членом асимптотического разложения данного решения [8].

В тесте в узле скважины и в узлах в окрестности круга радиуса К = 0.4 задавались значения, вычисленные по формуле (13). По результатам теста численное решение отличается от точного на величину неустранимой погрешности.

Проведенные тесты показали, чтопреобра-зования разностных уравнений для задания известных значений функции позволяют сохранить все

преимущества многосеточного алгоритма. Предложенный алгоритм применим для задания значений на внутренних и криволинейных границах при решении задач стандартным геометрическим многосеточным методом.

Литература

1. Д.В. Шевченко Применение многосеточных методов для расчета давления в нефтяном пласте. // Математическое моделирование. 2002. Т.14, №8. - С. 61-66.

2. K. Stuben, U. TrottenbergMulti-grid methods: fundamental algorithms, model problem analysis and applications // Multi-grid methods, Springer lecture notes in mathematica. New York, Springer Verlag. № 960. 1982. - P. 1-176.

3. А.Н. Чекалин Численные решения задач фильтрации в водонефтяных пластах. Казань, Изд-воКГУ, 1982. -208 с.

4. Г.И. Баренблатт, В.М. Ентов, В.М. Рыжик Движение жидкостей и газов в природных пластах. М., Недра, 1984. - 211 с.

5. Д.В. Шевченко, О.В. Панченко Учет в численных алгоритмах последствий применения методов значительного улучшения свойств призабойной зоны скважин. // Вестник Казанского технологического университета. 2014. Т.17, №6. - С. 267-269.

6. В.Я. Булыгин Гидромеханика нефтяного пласта. М.: Недра, 1974. - 232 с.

7. А.А. Самарский, Е.С. Николаев Методы решения сеточных уравнений. М.: Наука, 1978. - 592 с.

8. В.Б. Андреев, С.А. Кряквина О функции источника сеточного оператора Лапласа. // ЖВМ и МФ. 1972г., Т.12, №2. -С.364-373.

© И. Т. Усманов - канд. физ.-мат. наук, доцент, Елабужский институт (филиал) К(П)ФУ, [email protected]; Д. В. Шевченко - канд. физ.-мат. наук, доц. каф. интеллектуальных систем и управления информационными ресурсами КНИТУ, [email protected].

© I. T. Usmanov - candidate of physical-mat. sciences, Elabuga Institute, Kazan Federal University, [email protected]; D. V. Shevchenko - candidate of physical-mat. sciences, Associate Professor of the Department of Intelligent Systems & Information Systems Control, KNRTU, [email protected].