ОФП ещё раз про любовь — часть четвертая — Практическое моделирование

Казалось бы, знание теории 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;
                }

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