ОФП ещё раз про любовь — часть четвертая
Казалось бы, знание теории BL уже позволяет нам строить фронт вытеснения, однако для расчета показателей работы скважины возникает необходимость численного вычисления интегральных функций.
В книге «Waterflooding» (1986), в приложении А приведен листинг программы на FORTRAN66, однако мне показалось, что можно переписать проще. В книге также приведен пример расчета, который использовался для контроля правильности написания.
Переписанный код на C# githib
Относительные фазовые проницаемости, описаны функцией Кори
k_{ro}=\alpha_1(1-S_{wd})^m k_{rw}=\alpha_2 (S_{wd})^n
public double GetKro(double Swd)
{
return a1 * Math.Pow((1 - Swd), m);
}
public double GetKrw(double Swd)
{
return a2 * Math.Pow(Swd, n);
}
где нормирование насыщенности проводится в интервале от начальной до остаточной нефтенасыщенности,
S_{wd}=\frac{S_w-S_{iw}}{1-S_{or}-S_{iw}}
public double GetSwd(double Sw)
{
return (Sw - Swi) / (1 - Sor - Swi);
}
Доля воды в потоке (обводненность),
f_w=\frac{(S_{wd})^n}{(S_{wd})^n+ A \cdot (1-(S_{wd})^m}где,
A = \frac{\alpha_1 \cdot \mu_w}{\alpha_2 \cdot \mu_o}
public double Getfw(double Swd) // Функция fw в точке Swd
{
double A = (a1 * viscw) / (a2 * visco);
return Math.Pow(Swd, n) / (Math.Pow(Swd, n) + A * Math.Pow((1 - Swd), m));
}
Значение производной от обводненности
\frac{\partial f_w}{\partial S_w}=\frac{A}{1-S_{or}-S_{wi}} \cdot \frac{n (S_{wd})^{n-1}(1-S_{wd})^m + m (S_{wd})^n (1-S_{wd})^{m-1}}{((S_{wd})^n + A (1-(S_{wd})^m)^2}
public double GetdfwdSw(double Swd) // Производная fw` в точке Swd
{
double A = (a1 * viscw) / (a2 * visco);
return (1 / (1 - Sor - Swi)) * A * (n * Math.Pow(Swd, (n - 1)) * Math.Pow((1 - Swd), m) + m * Math.Pow(Swd, n) * Math.Pow((1 - Swd), (m - 1))) / (Math.Pow((Math.Pow(Swd, n) + A * Math.Pow((1 - Swd), m)), 2));
}
Определение значения насыщенности на фронте вытеснения происходит простым итеративным методом половинного деления.
Зная, что решение находится в интервале нормированной насыщенности от нуля до единицы, задается начальное приближение Swd = 0.5 и строится уравнение касательной к графику функции fw в этой точке. Если насыщенность подобрана верно, то касательная упадет в точку Swd = 0. В зависимости от того, перелетела ли или не долетела касательная, делается вывод об изменении насыщенности на следующей итерации.
public double GetSwf(double Sw)
{
double L = 0;
double R = 1;
double C = 0;
double Swd; // Нормированная насыщенность
double dfwdSw;
double Y = 1;
double eps = 1e-5;
while (Math.Abs(Y) > eps)
{
C = (L + R) * 0.5;
Swd = GetSwd(C);
dfwdSw = GetdfwdSw(Swd);
// Проводим линию касательной к функции fw, если в точке Swi она не пересекает ноль, значит решение пока не найдено
Y = Getfw(Swd) + dfwdSw * (Sw - C);
if (Y < -eps)
{
L = C;
}
if (Y > +eps)
{
R = C;
}
}
return C;
} // Получить насыщенность на фронте вытеснения
Это всё что касается вспомогательных подпрограмм.
Рассмотрим технологический режим разработки поддержания постоянного перепада давления между нагнетательной и добывающей скважиной. Дебит жидкости в элементарном сечении пласта можно записать как,
q_t=\left (\frac{k_{ro}}{\mu_o} + \frac{k_{rw}}{\mu_w} \right ) \cdot k \cdot A \cdot \frac{dp}{dx}Приведем к привычному уравнению для дебита,
q_t= \frac{k}{\mu^*} \cdot A \cdot \frac{dp}{dx}Запишем в целом для пласта,
q_t \int_0^L \mu^* dx = k \cdot A \int_{p_i}^{p_p} dpИзменение дебита связано с изменением вязкости жидкости,
q_t= \frac{k}{\mu_t} \cdot A \cdot \frac{\Delta p}{L}Значение вязкости жидкости в интервале насыщенности от заданной Sw до остаточной нефтенасыщенности, находится интегрированием методом трапеции.
public double GetLaAverage(double Sw) // Определение средней вязкости в интревале насыщенности от Sw до (1 - Sor)
{
// Решение интеграла методом трапеции
int N = 49;
double dSw = (1 - Sor - Sw) / N;
double Swd;
double Kro, Krw;
double La1, La2;
double dfwdSw1, dfwdSw2;
double dfwdSwi;
double LaSum;
double h;
Swd = GetSwd(Sw);
Kro = GetKro(Swd);
Krw = GetKrw(Swd);
La1 = 1 / (Kro / visco + Krw / viscw);
if (dSw == 0) return La1;
dfwdSw1 = GetdfwdSw(Swd);
dfwdSwi = dfwdSw1;
LaSum = 0;
for (int i = 0; i < (N); ++i)
{
Sw = Sw + dSw;
if (Sw > (1 - Sor)) Sw = 1 - Sor;
Swd = GetSwd(Sw);
Kro = GetKro(Swd);
Krw = GetKrw(Swd);
La2 = 1 / (Kro / visco + Krw / viscw);
dfwdSw2 = GetdfwdSw(Swd);
h = dfwdSw1 - dfwdSw2;
LaSum = LaSum + 0.5 * (La1 + La2) * h;
La1 = La2;
dfwdSw1 = dfwdSw2;
}
return LaSum / dfwdSwi;
}
}
На первом этапе расчета, фронт вытеснения нагнетаемой воды не достиг добывающей скважины. Наблюдается период безводной добычи нефти.
Дебит жидкости изменяется обратно пропорционально вязкости системы, по мере прохождения фронта вытеснения,
q_t= \frac{k}{\mu_t \cdot x/L + \mu_o \cdot (1 - x/L)} \cdot A \cdot \frac{\Delta p}{L}Вязкость жидкости в заводненной части остается постоянной и определяется нахождением вязкости от насыщенности на фронте вытеснения Swf до значения остаточной нефтенасыщенности (1-Sor).
В примере используются значения,
\mu_o=2,\mu_w=1Для заданных параметров Кори, насыщенность на фронте вытеснения,
S_{wf}=0.6639И вязкость жидкости,
\mu_t=3.2642Сравнивая полученную вязкость с вязкостью воды, можно отметить, кажущейся вязкость нагнетаемой воды в более чем 3.2 раза выше вязкости воды.
В советской литературе (Борисов Ю.П), использовалось понятие коэффициент увеличения фильтрационных сопротивлений, который показывал, во сколько раз возрастает сопротивление в зоне замещения нефти водой.
\alpha=\frac{\mu_w}{\mu_o} \cdot (1.7+8 \cdot S^*+25 \cdot (S^*)^2) \frac{\mu_o}{\mu_w} \leqslant 10где
S^*=1-S_{wf}-S_{or}Определив который, можно рассчитать и среднюю вязкость,
\mu_t=\alpha \cdot \mu_oПодстановка значений в формулу дает такое же значение вязкости 3.2.
Коэффициент безводной нефтеотдачи или же доля отбора нефти из пласта до начала обводнения от порового объема,
Q_{ibt}=\frac{1}{\displaystyle \left (\frac{\partial f_w}{\partial S_{wf}} \right )}Вязкость системы будет линейно изменятся от вязкости нефти в начальный момент времени, до вязкости жидкости в заводненной части в момент прорыва воды,
\mu=\mu_o + (\mu_t-\mu_o)\cdot \frac{Q_i}{Q_{ibt}}В тексте программы безводный период делится на 10 равных частей, для каждой расчетной точки определяется дебит жидкости, который равен дебиту нефти, накопленная добыча нефти в процентах от порового объема и накопленная добыча нефти.
Последняя неопределенная величина, это время.
В каждой расчетной точке нам известен мгновенный дебит жидкости и накопленные добычи жидкости. Сначала определим средний дебит жидкости за период времени,
\tilde q = \frac{q^{n+1} + q^{n}}{2}Затем добычу жидкости за рассмотренный период времени,
Q = (Q_i^{n+1} - Q_i^{n}) \cdot V_pгде поровый объем обозначен как,
V_p=A \phi Lнакопленная закачка воды, выраженная в поровых объемах записана через отношение,
Q_i=\frac{W_i}{V_p}Время которое проходит при переходе от одной расчетной точки к другой,
\Delta t = \frac{Q}{\tilde q}Далее просто накапливаем время,
t^{n+1}=t^n+\Delta tСобственно, всё готово для расчета показателей разработки в безводный период.
double Qibt = 1 / BL.GetdfwdSw(BL.GetSwd(Swf));
double dQibt = Qibt / 10;
double Labt = BL.GetLaAverage(Swf);
double time = 0;
double Qi_t_prev = 0;
double qt_t_prev = 0;
double qw = 0;
double qo = 0;
double WCT = 0;
for (int J = 0; J <= 10; ++J) // Десять расчетных шагов до начала обводнения
{
double Qi_t = dQibt * J;
double La_t = BL.visco + (Labt - BL.visco) * J / 10;
double qt_t = BL.GetLiquidRate(La_t);
if (J > 0)
{
time = time + 2 * (Qi_t - Qi_t_prev) * 160285 / (qt_t + qt_t_prev); // PV = 160285
}
qo = qt_t;
if (J == 10) // Прорыв в добывающую скважину
{
var fw = BL.Getfw(BL.GetSwd(Swf));
qw = qt_t * fw;
qo = qt_t - qw;
WCT = fw;
}
text.WriteLine($"{time:F1}\t{qo:N2}\t{qw:N2}\t{qt_t:N2}\t{WCT:N3}\t{Qi_t:N3}\t{Qi_t:N3}\t{Qi_t*160.285:N3}\t{La_t:N4}");
Qi_t_prev = Qi_t;
qt_t_prev = qt_t;
}
Прорыв фронта вытеснения происходит при отборе 30.2% от запасов нефти.
Второй этап расчета, это период обводнения, когда в добываемой продукции обводненность постоянно возрастает по мере отбора запасов нефти.
Обводненность продукции зависит от значения насыщенности в прискважинной зоне. В момент прорыва фронта вытеснения в скважину, обводненность определяется насыщенностью на фронте вытеснения. Для условий нашего примера, это значение
f_w = 0.896Таким образом, обводненность изменяется мгновенно от 0% до 89.6%. В дальнейшем водонасыщенность продолжит расти, до достижения значения остаточной нефтенасыщенности. Для расчета интервал изменения нефтенасыщенности делится на 50 равных частей.
Накопленная добыча нефти, в поровых объемах определяется для каждой точки насыщенности на фронте вытеснения,
Q_i=\frac{1}{\displaystyle \left (\frac{\partial f_w}{\partial S_{w}} \right )}И вязкость жидкости, определяется в интервале от текущей насыщенности до остаточной нефтенасыщенности.
Это позволяет построить зависимость вязкости системы от накопленной закачки воды, что позволяет определить дебит жидкости. Зная обводненность нетрудно определить и дебит нефти.
До прорыва воды в добывающую скважину, вязкость системы увеличивалась. Однако после прорыва воды, вязкость начинает снижается и после двух поровых объемов снижается ниже вязкости нефти 2 cP. Таким образом в нашем примере, появление воды в продукции скважины увеличивает продуктивность скважины, что сказывается на росте дебита жидкости.
N = 51;
dSw = (1 - BL.Sor - Swf) / (N - 1);
Sw = Swf;
for (int i = 0; i < N - 2; ++i)
{
Qi = 1 / BL.GetdfwdSw(BL.GetSwd(Sw));
var La = BL.GetLaAverage(Sw);
var qt_t = BL.GetLiquidRate(La);
var fw = BL.Getfw(BL.GetSwd(Sw));
var Swav = Sw + (1 - fw) / BL.GetdfwdSw(BL.GetSwd(Sw));
qw = qt_t * fw;
qo = qt_t - qw;
WCT = fw;
if (WCT > 0.999) break;
if (i > 0)
{
time = time + 2 * (Qi - Qi_t_prev) * 160285 / (qt_t + qt_t_prev);
text.WriteLine($"{time:F1}\t{qo:N2}\t{qw:N2}\t{qt_t:N2}\t{WCT:N3}\t{Qi:N3}\t{Qi:N3}\t{160.285 * (Swav - BL.Swi):N3}\t{La:N4}");
}
Qi_t_prev = Qi;
qt_t_prev = qt_t;
Sw = Sw + dSw;
}
Объединяя показатели для безводного периода и периода после прорыва воды, получим показатели работы скважины.
Характеристика вытеснения нефти водой,
Расчет накопленной добычи нефти производится через изменение средней водонасыщенности в пласте,
\frac{1}{Q_i} = \frac{1-f_w}{\bar{S}_w-S_{wi}} \bar{S}_w = S_{wi} + Q_i \cdot (1-f_w)