ОФП ещё раз про любовь — часть четвертая
Казалось бы, знание теории BL уже позволяет нам строить фронт вытеснения, однако для расчета показателей работы скважины возникает необходимость численного вычисления интегральных функций.
В книге «Waterflooding» (1986), в приложении А приведен листинг программы на FORTRAN66, однако мне показалось, что можно переписать проще. В книге также приведен пример расчета, который использовался для контроля правильности написания.
Переписанный код на C# githib
Относительные фазовые проницаемости, описаны функцией Кори
k_{ro}=\alpha_1(1-S_{wd})^m k_{rw}=\alpha_2 (S_{wd})^npublic 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)