CC BY
36
УДК 622.276
Исследование влияния трансформации двухфазной фильтрации на формирование зон невыработанных запасов нефти
С.И.ГРАЧЕВ, В.А.КОРОТЕНКО, Н.П.КУШАКОВАИ
Тюменский индустриальный университет, Тюмень, Россия
В статье с целью исследования процесса фильтрации флюидов при заводнении нефтяного месторождения используется модель Рапопорта - Лиса непоршневого вытеснения нефти водой. При плоскорадиальной фильтрации в однородном пласте определены радиусы зон возмущений с учетом и без учета концевого эффекта. Выявлено влияние изменения величины градиента капиллярного давления на распределение коэффициента водонасыщенности в зоне непоршневого вытеснения для высоко- и низкопроницаемых коллекторов. Применение модели элемента пятиточечной системы размещения нагнетательной и добывающих скважин показало, что при традиционной технологии заводнения разработки плоскорадиальная фильтрация жидкости трансформируется в прямолинейно-параллельную. При решении уравнения водонасыщенности использован метод интегральных соотношений Баренблатта, позволяющий определить время трансформации. Решением уравнения насыщенности для прямолинейно-параллельной фильтрации определено изменение величины коэффициента водонасыщенности на забое добывающей скважины для неограниченной и замкнутой залежи. Показано, что увеличение коэффициента обводненности добывающей скважины возможно только для замкнутого пласта. Для определения коэффициента водонасыщенности в замкнутой залежи получено дифференциальное уравнение с переменными коэффициентами, предложен итерационный метод решения. В элементе пятиточечной системы выявлены нефтенасыщенные зоны, не охваченные разработкой. Для каналов низкого фильтрационного сопротивления установлены условия их размещения в горизонтальной и вертикальной плоскостях. Показано, что при поддержании пластового давления в пласте существует линия изобар, соответствующая начальному пластовому давлению, расположение которой определяет направление скоростей перетоков флюидов. Интенсивность перетоков влияет на эффективность применения гидродинамических, физико-химических, тепловых и других методов увеличения нефтеотдачи.
Ключевые слова: модель Рапопорта - Лиса; метод Баренблатта; коэффициент водонасыщенности; трансформация фильтрационного процесса; остаточные подвижные запасы нефти
Как цитировать эту статью: Грачев С.И. Исследование влияния трансформации двухфазной фильтрации на формирование зон невыработанных запасов нефти / С.И.Грачев, В.А.Коротенко, Н.П.Кушакова // Записки Горного института. 2020. Т. 241. С. 68-82. DOI: 10.31897/РМ1.2020.1.68
Введение. Модель вытеснения нефти водой хорошо изучена в классических работах по подземной гидродинамике, тем не менее исследование процесса извлечения нефти при заводнении является одной из сложнейших задач проектирования и контроля разработки нефтяных месторождений. Это обусловлено тем, что при закачке воды в призабойной зоне пласта (ПЗП) нагнетательной скважины имеет место концевой эффект (КЭ). Вода накапливается в области радиусом гк, в которой коэффициент водонасыщенности возрастает до определенного значения sk. При превышении этой величины вода поступает в пласт, начинается процесс двухфазной фильтрации и непоршневого вытеснения нефти водой. На внутренней границе области КЭ вода поступает, на внешней - вода не вытекает, а вытекает нефть. Исследование фильтрационного процесса осложняется нелинейным граничным условием на внешней области КЭ [14]. В связи с этим концевым эффектом пренебрегают, коэффициент водонасыщенности на забое нагнетательной скважины принимает предельное (максимальное) значение s . Фильтрационно-емкостные свойства (ФЕС) пласта, зависящие от коэффициентов насыщенности, геологическое строение залежи, размещение скважин влияют на изменение характера течения флюидов в пласте и, следовательно, на выработку подвижных запасов нефти.
Постановка проблемы. Проблема формирования невыработанных остаточных подвижных запасов (ОПЗ) нефти обусловлена физическими свойствами флюидов, реализованной системой разработки, макронеоднородностью фильтрационных параметров по латерали и толщине нефтяного пласта [2, 8, 9, 12, 15, 16, 18]. Поверхностные явления, возникающие на границе раздела фаз, градиент капиллярного давления, гидрофильность или гидрофобность коллектора (микронеоднородность), существенно влияют на отклонение фактических показателей разработки от проектных технологических [7, 11]. Для устранения несоответствия между показателями используются
технологии (методы) повышения нефтеотдачи (МПН). При выборе и оценке эффективности соответствующего метода необходимо учитывать изменение (трансформацию) фильтрационных процессов и влияние степени насыщенности флюидами емкостного пространства коллектора.
Методология. Физические модели многофазной фильтрации непоршневого вытеснения нефти рассмотрены в работах [4, 6, 10, 13]. Из уравнений неразрывности, движения и реологических уравнений получена система дифференциальных уравнений в частных производных, решение которой сводится к определению распределения коэффициентов насыщенности и давления фаз в области возмущения - двухфазной фильтрации. Для связи коэффициента водонасыщенно-сти с давлением вводится дополнительное замыкающее уравнение. Все рассматриваемые модели можно разделить на две группы: учитывающие и неучитывающие градиент капиллярного давления. Решение систем дифференциальных уравнений, как правило, осуществляется приближенными численными методами. В предлагаемой работе для решения уравнений Рапопорта - Лиса используется метод интегральных соотношений Баренблатта [1]. Рассмотрена трансформация зон фильтрации в процессе разработки эксплуатационного объекта (ЭО). Полученные приближенные аналитические решения значений коэффициента водонасыщенности позволяют выделить зоны существования двухфазной фильтрации и остаточной нефти.
Обсуждение. Результаты исследования. Для исследования трансформации фильтрационных процессов рассмотрим модель заводнения однородного пласта. С вводом в эксплуатацию нагнетательной скважины происходит плоскорадиальное непоршневое вытеснение нефти водой. Давление на забое нагнетательной скважины больше начального пластового и забойного давлений добывающей скважины. С течением времени плоскорадиальная фильтрация трансформируется в прямолинейно-параллельную фильтрацию между нагнетательной и добывающей скважинами [17]. Для определения значений коэффициента водонасыщенности в модели, предложенной Рапопортом - Лисом, частным случаем которой является модель Бакли - Леверетта, пласт и флюиды считаются несжимаемыми. Задача сводится к решению дифференциального уравнения относительно коэффициента насыщенности с учетом капиллярного давления на границе раздела фаз и соответствует жесткому упруговодонапорному режиму:
т0 * + Ж (г )дЙ + ^ div
дг дх ( 2
♦ I др д5
к* (я)/ (01 —-— + Дpg вт а + Дgo
дя дх
= 0, (1)
где т0 - коэффициент открытой пористости; я - коэффициент водонасыщенности в зоне двухфазной фильтрации; Ж(г) - суммарная скорость фильтрации фаз; х - координата, для плоскорадиальной фильтрации х заменяется на г; /(я) - функция Бакли - Леверетта; к0 - коэффициент абсолютной проницаемости; к2 (я) - относительная фазовая проницаемость (ОФП) коллектора по нефти; рк - капиллярное давление; Ар = р2 - рь р2 и р1 - плотности нефти и воды; g - ускорение свободного падения; Аg0 = g2 - gl, g2 и gl - начальные градиенты давления при вытеснении нефти водой в зоне двухфазной фильтрации.
Функция Бакли - Леверетта определяется экспериментально:
/(*)= к.( к' (Л}к*( ,, (0 (2) М5)+(0к2 ) (2
где (1, (2 - коэффициенты динамической вязкости воды и нефти; к*(5), к2(5) - ОФП коллектора по воде и нефти.
Уравнение (1) содержит производную функции 5(х,г) по времени, и оператор div соответствует уравнению параболического типа. Следовательно, для решения (1) можно применить метод интегральных соотношений Баренблатта [1].
Начало координат совместим с забоем нагнетательной скважины. Радиус зоны возмущения двухфазной фильтрации обозначим Я(г, г). Поскольку изменение коэффициента водонасыщенности соответствует изменению давления нагнетаемой в пласт воды, то решение находим в виде [6]
з(г,г)= ^0 + л(г; 5 = 5* -50, (3)
/1(rс, Щ))
где s0 - коэффициент остаточной водонасыщенности (связанной воды); s - предельное значение коэффициента водонасыщенности; гс - радиус скважины. Единичная функция Хэвисайда ^(0 равна нулю при t = 0 и единице при t > 1:
fi (r, R(t )) = -n ln
(
R(t I
+1; fi(rc, R(t))=-nln
R(t) -
'c v
R (t),
+1.
(4)
Отметим, что функция/\(г,Я({)) < 0 для гс < г < R(t).
Если концевым эффектом пренебречь, то функция s(r,t) удовлетворяет следующим краевым условиям:
ds
s(r,0) = i0, s(r,, t) = s(R(t) t) = s0 —
dr
= 0.
(5)
r=R(t)
Давление на границе зоны двухфазной фильтрации R(t) равно начальному давлению в залежи р0.
При плоскорадиальном вытеснении нефти водой уравнение (1) примет вид
m,
= -qJ_ df - k1 d_ dt r dr д 2 r dr
rk2 (s)f (s) [^ Is + Apg sin a + Ago
I ds dr
(6)
Q(t)
где =-, Q(t) - расход воды на нагнетательной скважине; h - толщина пласта.
2пИ
Для определения радиуса зоны двухфазной фильтрации R(t) проинтегрируем (6) по г в интервале от гс до R(t). Получим
Asn2 d ( R2 Л
4(n + 2) dt { fi(rc, R)
4 ()f (s)
Ид
kr.
m0^2
rk* (s) [%Is + Apg sina + Ag0
I ds dr
(7)
Из граничных условий (5) следует, что при г = Я({) коэффициент водонасыщенности s = s0. ОФП коллектора по нефти принимает максимальное значение к * (50) = к 2тах , а функция Бакли -
Леверетта /(s0) = 0. При г = гс, s = 5*, к*( 5 *) = 0 функция Бакли - Леверетта /(5*) = 1. В результате имеем обыкновенное дифференциальное уравнение:
d ( Я2 4(п + 2)
dt
fi (rc, R(t))
n
n 2m0As
4 s(t).
(8)
Если qs = const после разделения переменных и интегрирования, получим трансцендентное уравнение для определения R(t):
R(t) =
r -
4(n + 2)
m0 Asn2
fi k, R{t)hJ
1
r +-
4(n + 2) (]n Rt) -1
Л
m0 Asn ^ rc n j
4st.
(9)
Следовательно, при заданных граничных условиях (5) капиллярное давление и ОФП на величину радиуса зоны двухфазной фильтрации в явном виде не влияет.
Однако, если не пренебрегать концевым эффектом радиуса гк, то второе граничное условие (5) на забое нагнетательной скважины примет вид
^к, 0=*к <гк > гс.
При г = гк и 5 = 5к ОФП коллектора по нефти к* () Ф 0, а функция Бакли - Леверетта
/Ы < 1.
После интегрирования в интервале [гк, Я(0] правая часть уравнения (7) примет вид
n
Г „ V
r
R
R
r.
r
()+МЫ г/ (^)
т т0р 2
gAp sin а + Ag +
дРк
дs
д1 Л
1 N к гк ) т0
А _ к о к 2( я к )
Ак - '
2
gAp sin а + Ag +
дРк
дя
д1 Л
1 N к гк)
При концевом эффекте s(r,t) и ее производная равны
')=10+л(')А*'ЖЖ;
к ~ 50;
дя
пА
( гп Л 1 - к
дг г/1 (ге,Я(г)){ Я" (г)
и
Для а = 0 и Дg = 0 получим
А =
к0к2 (^к) пА1к дРк
№ /1 (Гк,Я) д1
1к
(10)
После подстановки (10) в правую часть уравнения (7), получим дифференциальное уравнение, отличающееся от (8) правой частью:
йг
(
Я2 (г)
л
/(гс, Я(г))
- 4(п + 2)
п2 т0 Ая
/ 1 к (г )(1 + А)к.
(11)
Из формулы (10) видно, что параметр Ак зависит от технологического параметра ц фазовой проницаемости нефти, производной капиллярного давления по 1 в точке % радиуса зоны двухфазной фильтрации Я. По приведенным в работе [3] данным, производная капиллярного давления по 1 равна тангенсу угла наклона к оси 1 ^а) в интервале 10 < 1 < 1 . Наибольшие значения производной капиллярного давления наблюдаются в окрестностях значений коэффициента водо-насыщенности 10 соответствующему фронту вытеснения нефти водой, удаленного от забоя нагнетательной скважины (рис.1). Значения производной изменяются в диапазоне -10 < tgа < -8. На забое нагнетательной скважины изменение капиллярного давления для высокопроницаемого коллектора (кривая 1) практически не влияет на вытеснение нефти: tga « -0,1. Для низкопроницаемого коллектора (кривая 2) tga « -1.
Действительно, рассмотрим интервал изменения радиуса зоны возмущения 4 < Я < 200 м. Для высокопроницаемого коллектора при следующих значениях: Q = 300 м3/сут; к0 = 200 мД; к * = 0,6; р2 = 2 мПа-с; tga = -0,1 МПа; Д1к = 0,72 значения параметра Ак изменяются в интервале 0,01< Ак < 0,06. Для низкопроницаемого коллектора при заданных значениях: Q = 8 м3/сут; к0 = 4 мД; к2* = 0,2; р2 = 2 мПа-с; tga = -1 МПа; Д1к = 0,27 значения параметра Ак изменяются в интервале 0,03 < Ак < 0,16. С увеличением радиуса Я(г) параметр Ак убывает. Поскольку Ак < 1, то для определения радиуса зоны возмущения можно использовать формулу (9).
Пример 1. Дано: гс = 0,1 м; h = 10 м; 10 = 0,2;
= 0,8; т0 = 0,2; Q = 300 м/сут.
рк, МПа -
0,2
0,15
0,1
0,05
0,2
0,4
0,6
0,8
Рис. 1. Зависимость капиллярного давления от коэффициента водонасыщенности
1 - для высокопроницаемого коллектора; 2 - для низкопроницаемого коллектора
1
к
к
Определение значений радиуса зоны возмущения по формуле (9) для разных значений п
Таблица1 Результаты расчетов приведены в табл.1.
Рассмотрим элемент эксплуатационного объекта, состоящего из нагнетательной и добывающих скважин (рис.2). Пока выполняется неравенство Я(0 + Я2(0 < L, где L - расстояние между скважинами; Я2(0 - радиус зоны возмущения 2 добывающей скважины, в которой происходит только фильтрация нефти, имеет место плоскорадиальная фильтрация. Поскольку давление нагнетания больше начального пластового давления и давления на забое добывающей скважины, основная часть закачиваемой воды будет перемещаться к зоне перфорации добывающей скважины. Фильтрация становится прямолинейно-параллельной. Давление на границе раздела зон 1 и 2 (точки С,) равно начальному пластовому давлению р0. Другая часть воды будет перетекать со скоростью У в зону 1. В области 2 возмущения добывающей скважины в линейный канал фильтрации будет поступать нефть со скоростью У2. Причем скорости перетоков с течением времени будут убывать. Отметим, что при заданных граничных условиях точки С, неподвижны.
Время трансформации t определится из условия [6, 7]
t, сут R(t), м
n = 1 n = 2 n = 3 n = да
0,01 3 3 3 2,2
0,1 14 12 11 8
1 50 42 39 30
2 73 61 56 44
10 176 147 135 105
100 607 506 464 361
365 1210 1006 923 718
R(t*)+ R2 (t*) = L.
(12)
По известным формулам
R2 (t) ;
X 2 =
k 0k 2 Д 2Р2
где С = const.
Подставляя R2 и выражение (9) в (12), получим время трансформации из трансцендентного уравнения:
t * =
L
CX 2
1 + . 1 + -
4(n + 2)qs [ R(t*) -1 CX 2m0nAs l гс n
(13)
B1
R
n
D
D4
B4
B2
R2
D2
R2
B3
D3
2
б
а
Рис.2. Схема трансформации фильтрационного процесса для элемента обратной пятиточечной системы заводнения: а: Я - радиус зоны 1 двухфазной фильтрации нагнетательной скважины А; Я2, - радиусы зон 2
возмущения добывающих скважин В; 3 - зона (область), не охваченная заводнением; б: У\, У2 - скорости перетоков флюидов в зонах 1 и 2; 1 - область влияния нагнетательной скважины А!
Пример 2. Дано: гс = 0,1 м; h = 10 м; 50 = 0,2; = 0,8; то = 0,2; Q = 300 м/сут; L = 400 м; *0 = 100 мД; k2* = 0,6; ц2 = 2 мПа-с; р2 = 10-3 1/МПа;
Таблица 2
Определение времени начала трансформации по формуле (13) для разных значений п
* сут Я, м
п = 1 п = 2 п = 3
3 3.3 3.4 92 80 75
0 „2 L = 400 м; С = 12.
Результаты расчетов t по формуле (13) приведены в табл.2.
Для п = 2 расстояние от нагнетательной скважины до границы трансформации R = 80 м, давление на границе равно р0. В зоне 1 имеет место двухфазная
фильтрация - вытеснение нефти водой по модели Рапопорта - Лиса. Радиус зоны возмущения добывающей скважины, в которой осуществляется однофазная фильтрация нефти, равен 320 м. С момента начала эксплуатации нагнетательной скважины время изменения характера фильтрационного процесса для данных ФЕС и заданных технологических параметров не превышает 4 сут. Поскольку давление на забоях добывающих скважин меньше давления в других точках залежи, то плоскорадиальная фильтрация трансформируется в прямолинейно-параллельную. Обводнение продукции добывающих скважин начнется по каналам А\Ви
Для описания линейного вытеснения нефти водой воспользуемся уравнением (1). Начало координат совпадает с забоем нагнетательной скважины. Краевые условия в этом случае изменятся, так как на отрезке [хс, R] в пласте уже произошло изменение коэффициента водонасыщен-ности, а на забое нагнетательной скважины для х = хс = гс коэффициент водонасыщенности 5 = 5 . С течением времени на границе х = R коэффициент водонасыщенности будет возрастать: 5 > 50. Поэтому рассмотрим 5(х, О на двух интервалах: 51(х,0 в интервале хс < х < 1(0 < R и 52(х,0 в интервале R < х < 1(0.
Для интервала хс < х < 1(0 < R решение находим в виде
(х, t )= эЯ (х, Я ) + <
х - х<=. + Ь (х - хс)2 .
/ ^) /2 (;) '
,(х, Я )= 50 + А
А (х, Я) . А (хс, Я);
А (х, Я )= -п 1п Я-( Я ]П +1; / (хс, Я )= -п 1п Я-И +1 .
(14)
Для Я < х < l(t)
( , х - Я (х - Я )п
52 (х, t ) = а0 + а1 ч + а„ _ , ч ; п > 2.
/ (t) п 1п (t)
(15)
Здесь Я - константа, соответствующая п (табл.2); с, Ь, а0, а^ ап - неизвестные параметры, определяемые из граничных условий и условий сшивки функций 51 и 52.
Граничные условия для функции 52 зададим на границе фронта вытеснения:
,(/(t), t) = ^0; ~~~
дх
= 0.
(16)
х=/^ )
Условия сшивки функций 51 и 52 для х = Я примут вид
Ч \ д5, д52 д2 5, д2
(Я,0 = 52(Я,t); -1 = 2
дх дх ' дх2 дх2 Из совместного решения системы (16), (17) для п = 2 параметры а, Ь, с равны
(, Я Л2 „ (, Я Л 2А5
а0 = 50 + а21 1--I ; а1 = -2а21 1--I; а2 =---г-;
0 0 2 1 / У I / / 2 А (хс, Я)
(17)
2А5
/
2
Ь = а + / (х., я к я )•
— I ; с = -2а -
4А5 /
А (хс, Я) Я'
(18)
5
После подстановки (18) в (15) получим
(I (О - х)2
52 0 = 50 + а2 /2 (^ •
(19)
Для практической апробации определения коэффициентов водонасыщенности и обводнения на забое добывающей скважины рассмотрим случай Я < х < I(Х), 5 = 52.
После подстановки 52 в уравнение (1) и последующего интегрирования по х в интервале от Я до /(Х) получим дифференциальное уравнение для определения линейного размера зоны двухфазной фильтрации /(Х):
(/ - я)3
/2
= ^ / (5)
1 + А
(/ - я)
к/ /2
А = 4Алк 0к 2*(О дрк
к/ 1^(0/ (хс, Я)
(20)
Здесь 5 значение коэффициента водонасыщенности, соответствующее х = Я. Заметим, что 5 > 50 величина переменная, которая с течением времени будет возрастать. Функция Бакли - Леве-ретта изменяется в интервале 0 < / (5) < 1 и также возрастает.
Проверочный расчет параметра Ак/ показывает, что для различных значений входящих в формулу величин, его порядок не превышает 10-4, и из (20) с ростом /(Х) второе слагаемое в квадратной скобке стремится к нулю. Поэтому уравнение для определения /(Х) имеет вид
Л
(/(О - Я)3
/ (0 2
-
а2 т0
(/(х) - Я)3 /(х)2
/ (Х)
Ж(Х) /(5— Ф(Х).
(21)
Я а 2т0
Откуда получим кубическое уравнение
/3 -(3Я + Ф)/2 + 3Я2/ - Я3 = 0.
Применяя метод Кардано, найдем /(Х):
/ (х) = Я + 3 ф(х)
1 + 3
1+
Я
.1,5
+ 3
ф(х) j 1
1 -
Я
Ф0,
1,5
(22)
(23)
Значение 5 определяется из формулы (19), которая для х = Я примет вид
(/(х) - Я)2
5 — 50 I ^а 2
/(х)2 '
/ > Я.
(24)
Для получения приближенного аналитического решения оценим функцию Ф(Х), входящую в правую часть (21). Пусть суммарная скорость фильтрации фаз Ж(Х) постоянна. Тогда, применяя теорему о среднем, можно записать
ф(х) — ¡/ № (5)х; — *
п т ~ ' —
15 < 5 < 5 < 5 .
а2т0 0
а 2 т0
(25)
Из (24) следует, что значение 5 зависит от времени, следовательно, 5 также возрастает с течением Х. Предельное значение
* -. (26)
5 = 50 + а2 = 50 + Я
1п--0,5
х„
а2 т0
Для 50 = 0,2; Д5 = 0,6 и Я = 80 м полу- Таблица 3
чим 5 = °,24-85. Таким образом, функция Зависимость изменения размера зоны линейной фильтрации 1(() Бакли - Леверетта /(5) в подынтегральном от времени в бесконечном пласте
выражении изменяется незначительно (слабоменяющаяся функция) и ее вынос за знак интеграла оправдан. В дальнейшем
примем 5 — 5.
Отметим, что при Х ^ ж формула (23) примет простой вид:
/ (х) — Я + Ф(х) — Я + -3^/(5)Х. (27) а2 т0
Из (23) и (27) следует, что после трансформации размер зоны двухфазной фильтрации /(Х) линейно зависит от времени. Формула (23) справедлива для малых времен, а для Я с Ф(Х) можно использовать (27).
Значения функции Бакли - Леверетта /(5) определяются, если известны формулы зависимости ОФП от коэффициента водонасыщенности. Согласно исследованиям Чень Чжун-сяна [16]
Х, сут 1(0, м (23) /(0, м (27) х > Я 52(х,0 (19)
0,1 80,1 80,03 81 0,2924
1 80,9 80,3 100 0,2913
2 82 81 300 0,2803
10 89 83 361 0,2771
100 169 117 500 0,2701
365 405 394 800 0,2560
3650 3333 3332 3333 0,2000
к* (5) —
'т- 0,2 ^
V 0,8 у
; 0,2 < 5 < 1; к2* (5) —
'0,85 - 5 ^2,8
0,85
(1 + 2,45); 0 < 5 < 0,85.
(28)
В этом случае для 5 = 0,2485 функция Бакли - Леверетта равна /(5) —1,81-10 4.
В табл.3 приведены результаты расчета /(Х) для разных значений времени.
Во втором и третьем столбцах приведены расчетные значения /(Х) в интервале изменения времени от 0,1 сут до 10 лет, рассчитанные по формулам (23) и (27). В четвертом столбце приведены значения координаты х от 81 до 3333 м, в пятом столбце - расчетные значения коэффициента водонасыщенности на этом интервале. Отметим, что если х = L - расстоянию между скважинами, то из формулы (19) следует, что при /(Х) ^ ж, = 50 + а2, х - постоянная величина. Это противоречит практике, поскольку обводненность увеличивается, следовательно, на забое добывающей скважины доля воды в потоке жидкости должна возрастать. Соответственно, размер зоны фильтрации должен быть ограничен значением /к = /(Хк), которое определяется из геометрического строения залежи.
Распределение давления в ограниченных залежах рассматривалось в работах [1, 16]. По аналогии изменение коэффициента водонасыщенности в замкнутой ограниченной залежи находим в виде
,(х, Х) — D(t >2 (х, ¡к ) — D(t)
50 1 а2
/к - х
V /к у
; Я < х < ¡к; D(0)— 1.
(29)
Неизвестную функцию D(t) определим после подстановки (29) в уравнение (1) и после интегрирования в интервале [Я,/к] получим дифференциальное уравнение:
^ — С - BD, ЛХ
(30)
где
С — У 5)-/ 5*; Н — (/к - Я )Т50 + 3 а,
т0 Н I 3
1 2а2 к0к2 (5Я )/ (5Я ) (/к - Я) дР
Ч - ^ 2
В —
т0 Н
Д 2
/2
д5
V /к у ; В < 0.
(31)
Здесь 5Я соответствует значению коэффициента водонасыщенности при х = Я, 5к соответствует х = /к, 5 > 5Я > 5к > 50. Для функций Бакли - Леверетта выполняется соотношение
/5) </5) ^ 1.
5
Я
Из (29) следует, что функция D(t) ограничена сверху, ее максимальное значение Дтах=а*/а0. Для постоянной суммарной скорости фильтрации W(t) = W и ограниченных слабозависящих от времени функций G, В решение (30) примет вид
D(t ) — Б [G-(G - Б) e Bt ].
Раскладывая экспоненту в ряд, получим приближенное решение:
D(t ) = 1 + (G - B)t.
Для расчета предлагается следующая итерационная схема:
1. По известным величинам lk и R из первой формулы (31) находим параметр H.
2. Зададим значения времени t и параметр D(t) = 1,01.
3. Из (29) определим
(32)
(33)
«з, (L, t ) — D(t >2 (L, lk ) — D(t )
S0 ^ a2
Su — s
— s„ (R, t ) — DS2 (R, h )— D(t )
' i-L,' '
V lk y " f
; R < L < lk ;
S0 ^ a2
lk - R
V lk
Sk — S31 (lk , t) — DS2 (lk , lk ) — D(t)S0 ,
где L - расстояние между скважинами.
4. Рассчитаем значения коэффициентов ОФП и функции Бакли - Леверетта по формулам (28) и (2) для и
дРк
5. Из рис.1 или функциональной зависимости находим
ds
, по формулам (31) определим
параметры В, G.
6. Подсчитаем новое значение коэффициента а32. Если [а31 - а32 ] > 10 3, то повторим процесс
с пункта 3) итерационной схемы.
Задавая новое значение времени, определяем возрастание коэффициента водонасыщенности
на забое добывающей скважины. Для определения коэффициента обводненности а по каналу А1В1 воспользуемся формулой
V иг
а = —; W = V + V,,
где V - скорость фильтрации воды; v2 - скорость фильтрации нефти; W- суммарная скорость потока жидкости на забое добывающей скважины. Если W постоянна, то изменяются слагаемые, скорость фильтрации нефти убывает, а скорость фильтрации воды возрастает.
Скорость фильтрации воды выражается через суммарную скорость и зависит от градиента капиллярного давления [15]:
— f (s )
ТГГ 1 k 2* (s )dpk ds W + k0— p 2 ds dx
Коэффициент обводненности
a — f (s )(1 + A); A — k(
k 2*(s ) dPk Wp 2 ds
SL
ds dx
— - D(t )k 0
k2 (s) 2a2 (lk- - L) dpk
WP2 lk ds
sL
. dp. ds
< 0.
(34)
2
s
R
v
L
При малых значениях 5, близких к 50, ОФП к2*(5) принимает значения, близкие к максимальным, а функция Бакли - Леверетта мала. При 5 ^ 5* функция /(5) ^ 1, а к*(5) ^ 0. Следовательно, при больших значениях коэффициента водонасыщенности на забое добывающей скважины функция Бакли - Леверетта соответствует коэффициенту обводненности продукции. Таким образом, экспериментальные работы по определению значений ОФП играют важную роль при моделировании технологических показателей и выборе режима эксплуатации скважин.
Время частичного обводнения продукции ХЬ скважины В1 по каналу А1В1 соответствует значению коэффициента водонасыщенности, определенного по формуле (29), при Х = ХЬ и х = Ь. После подстановки в (33) получим
Ь — ^; D{tL)—-2. (35)
С '/к -Ьл
50 1 а 2
/к
Параметр С определяется из формулы (31) при/(5Я) = 1, к*(5*) = 0.
Пример 3. Дано: Q = 200 м3/сут; h = 10 м; Я = 80 м; Ь = 400 м; /к = 1000 м; для значения функции /(5к) = 0,5 время обводнения продукции, подсчитанное по формуле (35), равно ХЬ = 6,3 сут. Если значение функции Бакли - Леверетта принять равным /(5к) = 0,9, то ХЬ = 31,7 сут. Если размер зоны фильтрации /к = 3000 м, то для значений функции Бакли - Леверетта /(5к) = 0,5 и /(5к) = 0,9 значения времен обводнения ХЬ = 16,1 сут и ХЬ = 90 сут соответственно.
Размер зоны фильтрации /к зависит от геологического строения залежи и существенно влияет на время обводнения. В модели Рапопорта - Лиса коэффициент пористости постоянен, флюиды и коллектор несжимаемы, поэтому линейными каналами фильтрации можно пренебречь. Следовательно, в областях 3 (рис.2, а) останутся не вовлеченные в разработку запасы нефти.
Для определения не охваченных разработкой запасов нефти в зонах 3 определим Х22 - время соприкосновения зон возмущения добывающих скважин ВВ2 в точке D1.
Расстояние между добывающими скважинами В1В2 = Ь2 =Ьл/2. Радиус зоны фильтрации нефти Я22 для однородного коллектора
Я22 — ^ — 1Т Ь. (36)
Обозначим через Х22 время, соответствующее достижению точки D1, - границе раздела зон фильтрации добывающих скважин:
Я22 — >/С%2Х22 , (37)
где Х22 - коэффициент пьезопроводности нефти; С - числовой коэффициент. Радиус Я2 и время начала трансформации Х удовлетворяют соотношению
Я2 — Ь - Я — дС/. (38)
Из формул (36) и (37) получим
Г>2 Т2
* _/ 22 *
Х22 — Х — Х--Г^ (39)
22 Я22 2(Ь - Я)2
Воспользуемся данными примера 2 и результатами расчета табл.2. Для Ь = 400 м, Я = 80 м, Х = 3,3 сут и Я2 = 320 м из (38) и (39) получим Х22 = 2,6 сут, Я22 = 283 м. Следовательно, зоны возмущения двух добывающих скважин соприкасаются в точке D1 раньше, чем соединяются в точке С1 зоны возмущения нагнетательной и добывающей скважин.
Сделаем существенное замечание. Пластовое давление в точке D1 при Х22 = 2,6 сут равно р0, которое в дальнейшем будет убывать. Следовательно, в интервале времени Х22 < Х < Х в зоне 3 начнется фильтрация нефти в направлении точки D1. Если из-за «близости» величин временем Х22 и Х пренебречь, то можно определить запасы нефти в зонах, не охваченных заводнением.
B
п
D
B2
Рис.3. Фрагмент элемента эксплуатационного объекта для определения площади зоны 3. Произвольная точка Da е кривой DA. Угол 0 < а < 45°.
Итак, на момент времени £ в элементе ЭО пятиточечной системы образуются три зоны с разным характером насыщения (рис.2).
1. В окрестности нагнетательной скважины 1-я зона имеет форму круга радиуса Я. Характер насыщения - вода и нефть.
2. Четыре зоны 2, каждая из которых имеет форму замкнутой области, ограниченной кривыми линиями С^О. Характер насыщения - нефть.
3. Четыре зоны 3, ограниченные кривыми СС-Оъ не охваченные разработкой, насыщенные нефтью, давление в которых равно начальному давлению р0.
Обозначим площадь зоны 3 через Я3. Рассмотрим равнобедренный прямоугольный треугольник В1В2А1, катеты которого равны L. Из рис.2 определяем площадь зоны 3:
1 %R 2
S3 = SA S A1C1C2 -2Sa = 2 L 4--2Sa
(40)
где Яа - площадь треугольника А1В1В2; ¿А1с1с2 - площадь сектора А1С1С2; Яа - площадь фигуры В1О1С1 с центром в точке В1, ограниченной дугой О1С1 с центральным углом а0 = 45°. Радиус кривизны О1С1 изменяется от Я22 до Я2. Площадь зоны Я3 зависит от точности определения (расчета) площади Яа.
Формула для определения Яа для переменного радиуса кривизны имеет вид
1 "-0
Sa= ^ J Г2da .
(41)
Рассмотрим несколько вариантов расчета площади фигуры В1О1С1 (рис.3). Как видно из рисунка, дуга О1С1 выпукла по отношению к центру В1. При изменении угла в интервале 0 < а < л/4 радиус Яа возрастает от Я22 = 283 м до Я2 = 320 м. Такому условию удовлетворяет функция
BDa= Ra= L-
sin I 45 -a
+ (L - R )tga.
(42)
cos a
Для интервала 0 < a < л/4 расстояние от нагнетательной скважины А\ до точки Da определяется соотношением
sin
ADa= L-
л
-a
cos a
+ Rtg a.
(43)
Из формул (42) и (43) следует, что радиус кривизны Яа и расстояние А1Ба до нагнетательной скважины А1 зависят от угла а. Результаты расчета приведены в табл.4.
Пример 4. Дано: Я = 80 м; L = 400 м; Я2 = 320 м; Я22 = 283 м; 0 < а < л/4.
Таблица 4
Расчетные значения радиусов кривизны фигуры BjDjCj и расстояния ^4iDa, м
Параметры Номер расчета
1 2 3 4 5 6 7
a, рад 0 0,087 0,175 0,349 0,524 0,698 0,785
Ra, (42) 283 286 289 296 304 314 320
ADa, (43) 283 265 247 209 166 113 80
4
После подстановки (42) в (41) площадь B1D1C1 равна
5 »= 2
4 2
(
К -
л/2
Л
+1
ь
^1п2 + 2
(
К -
л/2
+1
ь
1
= 28682 м2. (44)
Если для определения площади 5а взять средний радиус гср, равный полусумме радиусов
К2 = 320 м и К22 = 283 м, то
5 <4--'
пг
ср
8
= 35661м2.
Таким образом, аппроксимация участка BlDlCl на сектор постоянного среднего радиуса дает завышенное значение площади, и, как следует из формулы (40), - меньшее значение площади зоны 3, что приведет к неверной оценке запасов.
В элементе эксплуатационного объекта (ЭО) В1В2В3В4 (см. рис.2) содержатся четыре зоны 3, каждая площадью 53. Объем нефти (геологические запасы)
V = 45 3^ (1 - ).
(45)
Для определения объема извлекаемых запасов нефти в зоне 3 примем коэффициент вытеснения нефти водой
Л: =
"н.н "ос.н _ " "0
1 - "п
где "н.н = 1 - "0 - начальный коэффициент нефтенасыщенности; "ос.н = 1 - - остаточный коэффициент нефтенасыщенности; "0 - коэффициент остаточной водонасыщенности; " - предельный коэффициент водонасыщенности, при котором фильтрация нефти прекращается. Объем извлекаемых запасов нефти
Уъь = Л1Л2^ = 4Л253кт0 * - "0 )
(46)
где л2 - коэффициент охвата пласта заводнением.
Пример 5. Пусть т0 = 0,2; к = 10 м; "0 = 0,2; " = 0,8; л2 = 1. Остальные данные возьмем из примера 3.
Из формул (40), (45), (46) получим 53 = 17613 м2, ¥3 = 12721 м3, = 84540 м3. Если входной дебит скважины, вскрывшей зоны 3, равен 50 м3/сут, то время, потраченное на извлечение подвижных запасов ¥3г2, равно 4 годам, при условии, что границы зоны 3 не изменятся. Объем геологических запасов нефти в элементе ЭО равен 512000 м3. Доля неизвлеченных запасов составит 0,22.
Однако с течением времени произойдет вторжение воды в зону 3, в которой возникнет двухфазная фильтрация, дебит нефти будет падать.
Рассмотрим, как изменятся фильтрационные процессы с течением времени при одновременной эксплуатации скважин. Целесообразно следующее замечание. В силу заданных граничных условий в элементе ЭО выполняется неравенство:
Рв1 < РD1 < PDa < Рс1 = Ра, = Р0 < Рлх
(47)
где произвольная точка Оа е зоне 3.
Условие (47) выполняется с момента времени начала формирования зон. Следовательно, заводнение зон 2 и 3 будет происходить неравномерно. С ростом времени t > t > t22 в зоне 3 начнется двухфазная фильтрация, граница С1Оа которой ограничена изобарой с начальным давлением р0. Зона однофазной фильтрации 2 деформируется, через границу по каналу C1Da начнет поступать закачиваемая в пласт вода.
Таким образом, в однородном по проницаемости в плоскости Оху коллекторе возникнет линейная фильтрация в зонах 3 и 2. Величины коэффициентов насыщенности и давлений в каналах фильтрации будут отличаться от первоначальных значений в начале эксплуатации.
2
2
2
2
*
н.н
Представим процесс фильтрации в элементе ЭО в виде системы каналов, соединяющих нагнетательную Л\ и добывающую В1 скважину.
В статье [5] показано, что обводнение продукции скважин осуществляется при наличии в пластах каналов низкого фильтрационного сопротивления (НФС). Рассмотрены случаи расположения каналов НФС по толщине пласта и в плоскостях Оху. Для определения длины каналов введены параметры ш7, зависящие от концентрации пиков индикатора и времени фиксации, 7 -порядковый номер фиксации индикатора. Для рассмотренного выше элемента пятиточечной системы разработки имеем
Lj = ш L = AxDa В, = — cos a
Откуда
2sinI —-a 1 + sin a
14 J
(48)
2sin I —-a 1 + sin a
ш =-^-J-. (49)
j cos a
Угол a = 45° соответствует прямой A1B1 = L - кратчайшему расстоянию между нагнетательной и добывающей скважинами, параметр Шj = 1. Угол a = 0° соответствует наибольшему расстоянию A1D1B1 =lV2 , параметр шj■ = 42 . Следовательно, каналы НФС, удовлетворяющие условиям
L < L} < la/2; 1 < ш ■ <42, (50)
лежат в горизонтальных плоскостях Оху. Если условия (50) не выполняются, то каналы НФС расположены в вертикальной плоскости.
В табл.3-5 статьи [5] рассмотрены случаи горизонтального расположения каналов НФС Lj для трех заданных параметров ш,.
В табл.5 приведены значения шj■ и Lj из работы [5, см. табл.3] и зависимости шj■ и A1DaB1 от угла a.
Таблица 5
Сравнение параметров Шу, Lj для L = 505 м [5] с расчетными значениями (48) (49)
Пики j fflj [5] Lj = rajL, м [5] a, град. fflj (49) A1D0B1 (48)
1 1 505 45 1 505
2 1,44 707 30 1,17 593
3 2,6 1313 20 1,26 638
4 3,8 1919 10 1,34 677
5 5,4 2828 0 1,41 714
Очевидно, что при заданном параметре ш 7 = (7 / ^ в горизонтальной плоскости находятся
два канала Lf. 7 = 1 и 7 = 2. Остальные каналы НФС располагаются в вертикальной плоскости. Заметим, что методика расчета объемов поровых каналов, необходимая для применения выравнивания профиля приемистости, приведена в статье [5].
Каждая модель имеет свои ограничения, свою область применения, в том числе и модель Бакли - Леверетта. Подтверждением достоверности является соответствие расчетных и фактических промысловых показателей [7]. Известно, что обводнение добывающей скважины происходит не «скачком», а постепенно увеличивается с возрастанием коэффициента водонасыщенности.
Если пласт слоисто-неоднородный, состоит из прослоев (пропластков) разной проницаемости, то для каждого /-го пропластка образуются свои зоны 1, 2, 3. Время трансформации фильтрации для каждого прослоя будет свое, определяемое по формулам (9), (12), (13). При отсутствии перетоков между пропластками (в рамках принятой модели Рапопорта - Лиса) возможны следующие граничные условия: расход закачиваемой в пласт воды следует распределить пропорционально коэффициентам проводимости пропластков [8, 15] или принять, что забойное давление постоянно во всех пропластках [11]. Отметим, что каждый случай требует отдельного исследова-
ния. Каналы НФС соответствуют высокопроницаемым пропласткам (ВП). Для стационарной фильтрации установлено наличие «языков» воды в наиболее проницаемых участках пласта, которое осуществляется по ВП. Изменение коэффициента водонасыщенности определяется по формулам (29), (31), (33).
Оригинальный метод исследования распределения остаточных подвижных запасов нефти в неоднородном по латерали пласте предложен М.И.Максимовым в монографии [8]. Показано, что в зависимости расположения зон проницаемостью 40, 200, 1000 мД в неоднородном по площади пласте имеет место «саморегулирование разработки». На наш взгляд, это следствие перетоков флюидов (см. рис.2, б), индивидуальное для каждой зоны.
В 50-х годах ХХ в. В.Н.Щелкачевым [18] сформулированы условия (критерии) применения форсированного отбора жидкости (ФОЖ), рекомендованные для эксплуатационных объектов (ЭО), находящихся на поздней стадии разработки: обводнение продукции не менее 90 %, высокие коэффициенты продуктивности, устойчивость коллектора и т.д. Считаем, что большое значение коэффициента обводненности не является решающим критерием определения поздней четвертой стадии ЭО.
Как известно, технология ФОЖ заключается в поэтапном изменении режимов работы добывающих и нагнетательных скважин, в создании высоких градиентов давлений в пласте, позволяющих вовлечь в разработку низкопроницаемые прослои, тупиковые, застойные зоны. Это приводит к изменению направлений фильтрационных потоков, а наличие гидродинамической связи между разнопроницаемыми пропластками осуществляет вовлечение в разработку низкопроницаемых разностей, т.е. циклическое заводнение. Следовательно, имеет место комбинированное нестационарное заводнение. Существование областей 3, сформированных в результате трансформации фильтрационных потоков и содержащих подвижные неизвлеченные запасы нефти, свидетельствует о необходимости применения методов увеличения нефтеотдачи (МУН), ранее не задействованных на данном ЭО. Существуют методики, определяющие плотность запасов нефти в областях, не охваченных разработкой, в частности способ Вороного. В многочисленных научных работах последних лет [15, 16] рекомендуется ввод в эксплуатацию вторых стволов с пологим и горизонтальным окончанием и многостадийный гидроразрыв пласта. Длина ствола в продуктивном пласте, как следует из приведенных выше расчетов, не должна превышать радиуса зоны дренирования К2 добывающей скважины (рис.2, б). Методы расчета горизонтальной дрены в пласте достаточно полно рассмотрены в литературе [15, 16].
Таким образом, изменение направления фильтрационных потоков ведет к формированию линейных каналов фильтрации между эксплуатационными скважинами и образованию зон, не охваченных разработкой (рис.3). Очевидно, что это необходимо учитывать при проектировании и расположении вторых стволов для эффективной выработки запасов нефти.
Поскольку забойное давление нагнетательных скважин превышает начальное пластовое р0, а забойное давление добывающих скважин, наоборот, меньше р0, то в пласте будет существовать линия изобар, соответствующая ему - 00. Причем существование линии 00 не зависит от выбранной физической модели. Расположение этой линии в пространстве зависит от граничных условий, влияет на формирование потоков фильтрации. С другой стороны, изобара С0 ограничивает области эффективного воздействия всех существующих МУН: гидродинамических, физико-химических, газовых, тепловых, ГРП, зон эффективной эксплуатации горизонтальных окончаний скважин и т.д. Наличие перетоков флюидов в горизонтальной плоскости (см. рис.2, б), а также перетоки в гидродинамически связанных пропластках в слоисто-неоднородных пластах влияют на эффективность применяемых технологий.
Выводы
1. Исследованием традиционной технологии заводнения с применением модели Рапопорта -Лиса непоршневого вытеснения нефти водой для плоскорадиальной фильтрации с учетом и без учета концевого эффекта установлено, что изменение градиента капиллярного давления практически не влияет на распределение коэффициентов насыщенности в зоне двухфазной фильтрации.
2. Моделированием одновременной работы нагнетательной и добывающей скважин выявлено, что плоскорадиальная фильтрация трансформируется в прямолинейно-параллельный поток. Определено время этого процесса и распределение коэффициентов водонасыщенности при прямолинейно-параллельной фильтрации в бесконечном и конечном (замкнутом) пласте.
3. На примере рассмотренной модели установлены зоны пласта, содержащие 22 % неизвле-ченных запасов нефти. В реальных условиях разработки эксплуатационного объекта могут быть получены более негативные результаты.
4. Для слоисто-неоднородного пласта без учета перетоков между пропластками предложена схема определения коэффициентов водонасыщенности в каждом прослое.
5. Установлены критерии возможного размещения каналов низкого фильтрационного сопротивления по латерали и толщине пласта.
ЛИТЕРАТУРА
1. Баренблатт Г.И. Движение жидкостей и газов в природных пластах / Г.И.Баренблатт, В.М.Ентов, В.М.Рыжик. М.: Недра, 1984. 211 с.
2. Желтое Ю.П. Разработка нефтяных месторождений. М.: Недра, 1998. 365 с.
3. Зозуля Г.П. Физика нефтяного и газового пласта / Г.П.Зозуля, Н.П.Кузнецов, А.К.Ягафаров; ТюмГНГУ. Тюмень, 2006. 244 с.
4. Каневская Р.Д. Математическое моделирование гидродинамических процессов разработки месторождений углеводородов. Москва-Ижевск: Институт компьютерных исследований, 2002. 140 с.
5. Коротенко В А. Интерпретация результатов трассерных исследований с учетом конвективного массопереноса / В.А.Коротенко, С.И.Грачев, А.Б.Кряквин // Записки Горного института. 2019. Т. 236. С. 185-193. DOI: 10.31897/PMI.2019.2.185
6. Коротенко В.А. Особенности фильтрации и вытеснения нефти из аномальных коллекторов / В.А.Коротенко, Н.П.Кушакова; Тюменский индустриальный университет. Тюмень, 2018. 150 с.
7. Крэйг Ф.Ф. Разработка нефтяных месторождений при заводнении: Пер. с англ. / Под ред. В.Л.Данилова. М.: Недра, 1974. 192 с.
8. МаксимовМ.И. Геологические основы разработки нефтяных месторождений. М.: Недра, 1975. 534 с.
9. Маскет М. Физические основы технологии добычи нефти. Москва-Ижевск: Институт компьютерных исследований, 2004. 606 с.
10. Механика насыщенных поровых сред / В.Н.Николаевский, К.С.Басниев, А.Т.Горбунов, Г.А.Зотов. М.: Недра, 1970. 336 с.
11. Мирзаджанзаде АХ. Физика нефтяного и газового пласта / А.Х.Мирзаджанзаде, И.М.Аметов, А.Г.Ковалев. Москва-Ижевск: Институт компьютерных исследований, 2005. 280 с.
12. Михайлов Н.Н. Остаточное нефтенасыщение разрабатываемых пластов. М.: Недра, 1992. 270 с.
13. Николаевский В.Н. Механика пористых и трещиноватых сред. М.: Недра, 1984. 232 с.
14. Подземная гидравлика / К.С.Басниев, А.М.Власов, И.Н.Кочина, В.М.Максимов. М.: Недра, 1986. 303 с.
15. ТелковА.П. Гидромеханика пласта применительно к нефтегазопромысловым задачам разработки месторождений наклонно-направленными и горизонтальными стволами / А.П.Телков, С.И.Грачев. СПб: Наука, 2012. 160 с.
16. Телков А. П. Гидромеханика пласта применительно к прикладным задачам разработки нефтяных и газовых месторождений: В 2 ч. / А.П.Телков, С.И.Грачев; ТюмГНГУ. Тюмень, 2009. Ч. 1. 240 с.
17. Трансформация процессов фильтрации при разработке залежей углеводородов / В.А.Коротенко, С.И.Грачев, Н.П.Кушакова, Р.Р.Сабитов // Успехи современного естествознания. 2017. № 2. С. 86-93.
18. Щелкачев В.Н. Разработка нефтеносных пластов при упругом режиме. М.: Гостоптехиздат, 1959. 467 с.
Авторы: С.И.Грачев, д-р техн. наук, профессор, grachevsi@tyuiu.ru (Тюменский индустриальный университет, Тюмень, Россия), В.А.Коротенко, канд. техн. наук, доцент, korotenkova@tyuiu.ru (Тюменский индустриальный университет, Тюмень, Россия), Н.П.Кушакова, канд. техн. наук, доцент, kushakovanp@tyuiu.ru (Тюменский индустриальный университет, Тюмень, Россия).
Статья поступила в редакцию 30.05.2019.
Статья принята к публикации 12.09.2019.
