Практическое моделирование

и другие вопросы разработки нефтяных месторождений
frac

Расчет скважины с трещиной бесконечной проводимости

10.09.2026

Продолжим исследовать приток флюида к плоскости трещины.

Функция мгновенного точечного источника, расположенного в точке (x_w,y_w)

    \begin{equation*} S(x,t) = \frac{1}{2 \sqrt{\eta}} \cdot \frac{1}{\sqrt{\pi t }} \cdot \exp \left [ -\frac{(x-x_w)^2}{4 \eta t} \right ] \end{equation*}

    \begin{equation*} S(y,t) = \frac{1}{2 \sqrt{\eta}} \cdot \frac{1}{\sqrt{\pi t }} \cdot \exp \left [ -\frac{(y-y_w)^2}{4 \eta t} \right ] \end{equation*}

Представим трещину как объединение точечных источников вдоль оси x, где под x_f будем понимать теперь полудлину трещины,

    \[ \Delta p(x,y,t) = \frac{1}{\phi c_t} \, \int\limits_0^t q_f(\tau) \int\limits_{-x_f}^{+x_f} S(x,t-\tau) \, dx_w \cdot S(y,t-\tau) d\tau \]

Получим,

    \[ \Delta p(x,y,t) = \frac{1}{4 \eta \pi \phi c_t } \, \int\limits_0^t q_f(\tau) \left ( \int\limits_{-x_f}^{+x_f} \exp \left [ -\frac{(x-x_w)^2 + (y-y_w)^2}{4 \eta (t - \tau)} \right ] dx_w \right ) \frac{d\tau}{t - \tau} \]

Трещину разместим в начале координат (y_w=0) и давление будет определяться только вдоль оси x, там где y=0, поэтому

    \[ \Delta p(x,t) = \frac{1}{4 \eta \pi \phi c_t } \, \int\limits_0^t q_f(\tau) \left( \int\limits_{-x_f}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] dx_w \right ) \frac{d\tau}{t - \tau} \]

Разделим время на N произвольных интервалов, в течение которых дебит равен некоторому среднему дебиту. Тогда к моменту t, изменение давления согласно принципу суперпозиции, будет складываться из суммы,

    \[ \Delta p(x,t) = \frac{1}{4 \eta \pi \phi c_t } \, \sum_{l=1}^{N} q_{f}^{l} \, \int\limits_{t_{l-1}}^{t_l} \left ( \int\limits_{-x_f}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] dx_w \right ) \frac{d\tau}{t - \tau} \]

Здесь предполагается, что приток жидкости равномерен по всей плоскости (uniform flux), поэтому дебит скважины q_f выносится за знак интеграла,

Внутренний интеграл можно разложить на положительное и отрицательное направление трещины,

    \[ \int\limits_{-x_f}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] dx_w = \int\limits_{0}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] dx_w + \int\limits_{-x_f}^{0} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] dx_w = \]

    \[=\int\limits_{0}^{x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t - \tau)} \right ] dx_w \]

Распишем сумму к моменту времени t_N,

    \[ \Delta p(x,t_N) = \frac{1}{4 \eta \pi \phi c_t } \Bigg [q_{f}^{1} \int\limits_{0}^{t_1} \bigg( \int\limits_{0}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t_N - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t_N - \tau)} \right ] dx_w \bigg) \frac{d\tau}{t_N - \tau} + \]

    \[+ q_{f}^{2} \int\limits_{t_1}^{t_2} \bigg ( \int\limits_{0}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t_N - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t_N - \tau)} \right ] dx_w \bigg ) \frac{d\tau}{t_N - \tau} + ... \]

    \[+ q_{f}^{N} \int\limits_{t_{N-1}}^{t_N} \bigg ( \int\limits_{0}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t_N - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t_N - \tau)} \right ]dx_w \bigg ) \frac{d\tau}{t_N - \tau} \Bigg ] \]

Обозначим коэффициентом G множитель после дебита,

    \[ \Delta p(x,t_N) = q_{f}^{1} G^{N,1} + q_{f}^{2} G^{N,2} + \cdots + q_{f}^{N} G^{N,N} \]

    \[G^{N,l} = \frac{1}{4 \eta \pi \phi c_t } \int\limits_{t_{l-1}}^{t_l} \bigg( \int\limits_{0}^{+x_f} \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t_N - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t_N - \tau)} \right ]dx_w \bigg) \frac{d\tau}{t_N - \tau} \]

Для всех временных шагов можно составить матрицу,

    \begin{equation*} \begin{bmatrix} {\Delta p}^1 \\ {\Delta p}^2 \\ {\Delta p}^3 \\ \vdots \\ {\Delta p}^N \end{bmatrix} = \begin{bmatrix} {G}^{1,1} & 0 & 0 & \dots & 0 \\ {G}^{2,1} & {G}^{2,2} & 0 & \dots & 0 \\ {G}^{3,1} & {G}^{3,2} & {G}^{3,3} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {G}^{N,1} & {G}^{N,2} & {G}^{N,3} & \dots & {G}^{N,N} \end{bmatrix} \begin{bmatrix} {q_f}^1 \\ {q_f}^2 \\ {q_f}^3 \\ \vdots \\ {q_f}^N \end{bmatrix} \end{equation*}

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

    \[ \mathbf{\Delta p}^k = \sum_{l=1}^{k-1} \mathbf{G}^{k,l} \cdot \mathbf{q_f}^l + \mathbf{G}^{k,k} \cdot \mathbf{q_f}^k \]

И так шаг за шагом получить решение системы.

***

Поделим одно крыло трещины на M сегментов.

Длина одного сегмента составит L=x_f/M. Каждый сегмент работает со своим равномерно распределенным дебитом q_{f,m} где (m=1,M). Первый сегмент располагается от 0 до L, второй сегмент от L до 2L, последний сегмент от (M-1)L до x_f.

    \[ \Delta p(x,t) = \frac{1}{4 \eta \pi \phi c_t } \, \sum_{l=1}^{N} \sum_{m=1}^{M} (q_f)_m^l \int\limits_{t_{l-1}}^{t_l} \int\limits_{(m-1)L}^{mL} \left ( \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t - \tau)} \right ] dx_w \right ) \frac{d\tau}{t - \tau} \]

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

    \[ (q_f)_m^l = (q_w)_m^l \cdot \frac{M}{x_f h} \]

    \[ q_w^l = 2 \sum_{m=1}^{M} (q_w)_m^l \]

    \[ q_w^l = \sum_{m=1}^{M} (q_f)_m^l \frac{2x_f h}{M} \]

Из последнего выражения записывается уравнение связывающее удельные дебиты сегментов с дебитом скважины,

    \[ M = 2 x_f h \sum_{m=1}^{M} \frac{(q_f)_m^l}{q_w^l} \]

Дебит скважины связан с удельным дебитом,

    \[ q_f^l = q_w^l \cdot \frac{1}{2 x_f h} \]

Удобней перейти к отношению дебита сегмента к дебиту трещины,

    \[ \Delta p(x,t) = \frac{q_f^l}{4 \eta \pi \phi c_t } \, \sum_{l=1}^{N} \sum_{m=1}^{M} \frac{(q_f)_m^l}{q_f^l} \int\limits_{t_{l-1}}^{t_l} \int\limits_{(m-1)L}^{mL} \left ( \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t - \tau)} \right ] dx_w \right ) \frac{d\tau}{t - \tau} \]

Введем следующие безразмерные параметры,

    \[ (q_D)_m^l = \frac{(q_f)_m^l}{q_f^l} \]

    \[ p_D = 2 \pi \cdot \frac{kh}{\mu} \cdot \frac{\Delta p }{q_w} \]

Что позволит записать,

    \[ (p_D)^l = \frac{1}{4 x_f} \, \sum_{l=1}^{N} \sum_{m=1}^{M} (q_D)_m^l \int\limits_{t_{l-1}}^{t_l} \int\limits_{(m-1)L}^{mL} \left ( \exp \left [ -\frac{(x-x_w)^2}{4 \eta (t - \tau)} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta (t - \tau)} \right ] dx_w \right ) \frac{d\tau}{t - \tau} \]

Напоследок введем ещё одно обозначение: \Delta t = t - \tau.

***

Пришло время разобрать внутренний интеграл по расстоянию,

    \[ \int\limits_{(m-1)L}^{mL} \bigg ( \exp \left [ -\frac{(x-x_w)^2}{4 \eta \Delta t} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4 \eta \Delta t} \right ] \bigg) dx_w \]

Здесь нам поможет следующее решение,

    \[ \int \exp \bigg( -\frac{\xi^2}{a^2} \bigg ) d\xi = \frac{a \sqrt{\pi} }{2} \operatorname{erf} \frac{\xi}{a}\]

    \[\int\limits_{(m-1)L}^{mL} \bigg ( \exp \left [ -\frac{(x-x_w)^2}{4\eta\Delta t} \right ] + \exp \left [ -\frac{(x+x_w)^2}{4\eta \Delta t} \right ] \bigg) dx_w = \]

    \[\sqrt{\pi\eta\Delta t} \left( \operatorname{erf} \frac{m L - x}{2\sqrt{\eta\Delta t}} - \operatorname{erf} \frac{(m-1)L - x}{2\sqrt{\eta\Delta t}} + \operatorname{erf} \frac{m L + x}{2\sqrt{\eta\Delta t}} - \operatorname{erf} \frac{(m-1)L + x}{2\sqrt{\eta\Delta t}} \right) \]

Для произвольного сегмента j, центр сегмента имеет координату x = jL-L/2,

    \[=\sqrt{\pi\eta\Delta t} \left( \operatorname{erf} \frac{x_f (m-j+0.5)}{2M\sqrt{\eta\Delta t}} - \operatorname{erf} \frac{x_f(m-j-0.5)}{2M\sqrt{\eta\Delta t}} + \operatorname{erf} \frac{x_f(m+j-0.5)}{2M\sqrt{\eta\Delta t}} - \operatorname{erf} \frac{x_f(m+j-1.5)}{2M\sqrt{\eta\Delta t}} \right) \]

Можно немного внести косметических изменений и вынести j перед m,

    \[=\sqrt{\pi\eta\Delta t} \left( -\operatorname{erf} \frac{x_f (j-m-0.5)}{2M\sqrt{\eta\Delta t}} + \operatorname{erf} \frac{x_f(j-m+0.5)}{2M\sqrt{\eta\Delta t}} + \operatorname{erf} \frac{x_f(j+m-0.5)}{2M\sqrt{\eta\Delta t}} - \operatorname{erf} \frac{x_f(j+m-1.5)}{2M\sqrt{\eta\Delta t}} \right) \]

Сделаем следующие замены,

    \[ \alpha = \frac{j-m+1/2}{2M} \]

    \[ \beta = \frac{j-m-1/2}{2M} \]

    \[ \gamma = \frac{j+m-1/2}{2M} \]

    \[ \delta = \frac{j+m-3/2}{2M} \]

    \[=\sqrt{\pi\eta\Delta t} \left( \operatorname{erf} \frac{x_f\alpha}{\sqrt{\eta\Delta t}} -\operatorname{erf} \frac{x_f \beta}{\sqrt{\eta\Delta t}} + \operatorname{erf} \frac{x_f\gamma}{\sqrt{\eta\Delta t}} - \operatorname{erf} \frac{x_f\delta}{\sqrt{\eta\Delta t}} \right) \]

Здесь напрашивается безразмерное давление,

    \[ \Delta t_D = \Delta t \cdot \frac {\eta}{x_f^2} \]

    \[=x_f \sqrt{\pi \Delta t_D} \left( \operatorname{erf} \frac{\alpha}{\sqrt{\Delta t_D}} -\operatorname{erf} \frac{\beta}{\sqrt{\Delta t_D}} + \operatorname{erf} \frac{\gamma}{\sqrt{\Delta t_D}} - \operatorname{erf} \frac{\delta}{\sqrt{\Delta t_D}} \right) \]

Получим промежуточный результат изменения давления в середине произвольного сегмента,

    \[ (p_D)_j^l = \frac{\sqrt{\pi}}{4} \, \sum_{l=1}^{N} \sum_{m=1}^{M} (q_D)_m^l \int\limits_{t_{l-1}}^{t_l} \frac{1}{\sqrt{\Delta t_D}}\left( \operatorname{erf} \frac{\alpha}{\sqrt{\Delta t_D}} -\operatorname{erf} \frac{\beta}{\sqrt{\Delta t_D}} + \operatorname{erf} \frac{\gamma}{\sqrt{\Delta t_D}} - \operatorname{erf} \frac{\delta}{\sqrt{\Delta t_D}} \right) d\tau_D \]

***

Разберем внутренний интеграл по времени, при помощи такого вот решения,

    \[ \int \frac{1}{\sqrt{t-u}} \operatorname{erf} \frac{a}{\sqrt{t-u}} du = -2\sqrt{t-u}\operatorname{erf}\bigg(\frac{a}{\sqrt{t-u}} \bigg) + \frac{2a}{\sqrt{\pi}} \operatorname{Ei}\bigg(-\frac{a^2}{t-u}\bigg) \]

Получим,

    \[ (p_D)_j^l = \frac{\sqrt{\pi}}{4} \, \sum_{l=1}^{N} \sum_{m=1}^{M} (q_D)_m^l \Bigg ( \]

    \[ \bigg(2\sqrt{\Delta t}\operatorname{erf}\bigg(\frac{\alpha}{\sqrt{\Delta t}} \bigg) - 2\sqrt{\Delta t}\operatorname{erf}\bigg(\frac{\beta}{\sqrt{\Delta t}} \bigg) + 2\sqrt{\Delta t}\operatorname{erf}\bigg(\frac{\gamma}{\sqrt{\Delta t}} \bigg) - 2\sqrt{\Delta t}\operatorname{erf}\bigg(\frac{\delta}{\sqrt{\Delta t}} \bigg) \Bigg |_{t_{l}}^{t_{l-1}} + \]

    \[ + \bigg(- \frac{2\alpha}{\sqrt{\pi}} \operatorname{Ei}\bigg(-\frac{\alpha^2}{\Delta t}\bigg) + \frac{2\beta}{\sqrt{\pi}} \operatorname{Ei}\bigg(-\frac{\beta^2}{\Delta t}\bigg) - \frac{2\gamma}{\sqrt{\pi}} \operatorname{Ei}\bigg(-\frac{\gamma^2}{\Delta t}\bigg) + \frac{2\delta}{\sqrt{\pi}} \operatorname{Ei}\bigg(-\frac{\delta^2}{\Delta t}\bigg) \bigg)\Bigg |_{t_{l}}^{t_{l-1}} \Bigg ) \]

Введем два крупных комплекса величин,

    \[ X_{j,m}^{N,l} = 2\sqrt{\Delta t_{N,l}} \left [ \operatorname{erf}\bigg(\frac{\alpha_{j,m}}{\sqrt{\Delta t_{N,l}}} \bigg) - \operatorname{erf}\bigg(\frac{\beta_{j,m}}{\sqrt{\Delta t_{N,l}}} \bigg) + \operatorname{erf}\bigg(\frac{\gamma_{j,m}}{\sqrt{\Delta t_{N,l}}} \bigg) - \operatorname{erf}\bigg(\frac{\delta_{j,m}}{\sqrt{\Delta t_{N,l}}} \bigg) \right ] \]

    \[ Y_{j,m}^{N,l} = - \frac{2}{\sqrt{\pi}} \bigg( \alpha_{j,m} \operatorname{Ei}\bigg(-\frac{\alpha_{j,m}^2}{\Delta t_{N,l}}\bigg) - \beta_{j,m} \operatorname{Ei}\bigg(-\frac{\beta_{j,m}^2}{\Delta t_{N,l}}\bigg) + \gamma_{j,m} \operatorname{Ei}\bigg(-\frac{\gamma_{j,m}^2}{\Delta t_{N,l}}\bigg) - \delta_{j,m} \operatorname{Ei}\bigg(-\frac{\delta_{j,m}^2}{\Delta t_{N,l}}\bigg) \bigg) \]

Итоговая формула определения падения давления в центре j сегмента трещины,

    \[ (p_D)_j^l = \frac{\sqrt{\pi}}{4} \, \sum_{l=1}^{N} \sum_{m=1}^{M} (q_D)_m^l\cdot (X_{j,m}^{N,l-1} - X_{j,m}^{N,l} + Y_{j,m}^{N,l-1} - Y_{j,m}^{N,l} ) \]

Разложение трещины на сегменты утяжеляет итоговую расчетную матрицу,

    \begin{equation*} \begin{bmatrix} (p_D)_1^k \\ (p_D)_2^k \\ \vdots \\ (p_D)_M^k \end{bmatrix} = \sum_{l=1}^k \begin{bmatrix} A_{1,1}^{k,l} & A_{1,2}^{k,l} & \dots & A_{1,M}^{k,l} \\ A_{2,1}^{k,l} & A_{2,2}^{k,l} & \dots & A_{2,M}^{k,l} \\ \vdots & \vdots & \ddots & \vdots \\ A_{M,1}^{k,l} & A_{M,2}^{k,l} & \dots & A_{M,M}^{k,l} \end{bmatrix} \begin{bmatrix} (q_D)_1^l \\ (q_D)_2^l \\ \vdots \\ (q_D)_M^l \end{bmatrix} \end{equation*}

Для текущего временного шага,

    \[ \begin{bmatrix} (p_D)_1^k \\ (p_D)_2^k \\ \vdots \\ (p_D)_M^k \end{bmatrix} = \sum_{l=1}^{k-1} \begin{bmatrix} A_{1,1}^{k,l} & A_{1,2}^{k,l} & \dots & A_{1,M}^{k,l} \\ A_{2,1}^{k,l} & A_{2,2}^{k,l} & \dots & A_{2,M}^{k,l} \\ \vdots & \vdots & \ddots & \vdots \\ A_{M,1}^{k,l} & A_{M,2}^{k,l} & \dots & A_{M,N}^{k,l} \end{bmatrix} \cdot \begin{bmatrix} (q_D)_1^l \\ (q_D)_2^l \\ \vdots \\ (q_D)_M^l \end{bmatrix} + \begin{bmatrix} A_{1,1}^{k,k} & A_{1,2}^{k,k} & \dots & A_{1,M}^{k,k} \\ A_{2,1}^{k,k} & A_{2,2}^{k,k} & \dots & A_{2,M}^{k,k} \\ \vdots & \vdots & \ddots & \vdots \\ A_{N,1}^{k,k} & A_{N,2}^{k,k} & \dots & A_{M,M}^{k,k} \end{bmatrix} \cdot \begin{bmatrix} (q_D)_1^k \\ (q_D)_2^k \\ \vdots \\ (q_D)_M^k \end{bmatrix} \]

Сумму предыдущих известных шагов обозначим матрицей B,

    \[ \begin{bmatrix} (p_D)_1^k \\ (p_D)_2^k \\ \vdots \\ (p_D)_M^k \end{bmatrix} = \begin{bmatrix} B_1^l \\ B_2^l \\ \vdots\\ B_M^l \end{bmatrix} + \begin{bmatrix} A_{1,1}^{k,k} & A_{1,2}^{k,k} & \dots & A_{1,M}^{k,k} \\ A_{2,1}^{k,k} & A_{2,2}^{k,k} & \dots & A_{2,M}^{k,k} \\ \vdots & \vdots & \ddots & \vdots \\ A_{M,1}^{k,k} & A_{M,2}^{k,k} & \dots & A_{M,M}^{k,k} \end{bmatrix} \cdot \begin{bmatrix} (q_D)_1^k \\ (q_D)_2^k \\ \vdots \\ (q_D)_M^k \end{bmatrix} \]

Решение задачи строится на добавлении замыкающего уравнения,

    \[\sum_{m=1}^{M} (q_D)_m = M \]

Отсутствие потерь давления при течении жидкости внутри трещины, при переходе от элемента к элементу, позволяет заменить собственные падения давления на одно, общее падение,

    \[(p_D)_1^k = (p_D)_2^k = \cdots = (p_D)_w^k \]

Теперь можно расширить матрицу до (M+1)(M+1),

    \[ \begin{bmatrix} A_{1,1}^{k,k} & A_{1,2}^{k,k} & \dots & A_{1,M}^{k,k} & -1\\ A_{2,1}^{k,k} & A_{2,2}^{k,k} & \dots & A_{2,M}^{k,k} & -1\\ \vdots & \vdots & \ddots & \vdots & -1 \\ A_{M,1}^{k,k} & A_{M,2}^{k,k} & \dots & A_{M,M}^{k,k} & - 1\\ 1 & 1 & \dots & 1 & 0 \end{bmatrix} \begin{bmatrix} (q_D)_1^k \\ (q_D)_2^k \\ \vdots \\ (q_D)_M^k \\ (p_D)_w^k \end{bmatrix} = \begin{bmatrix} -B_1^l \\ -B_2^l \\ \vdots\\ -B_M^l \\ M \end{bmatrix} \]

Сравнение для случая M=100 с данными опубликованными в 1974 показывает хорошее совпадение,

Распределение дебита вдоль полуоси трещины при t_D=1000,

Сравнение притока к трещине бесконечной проводимости и равномерного притока,

И реализация в коде,

     public void Calculate()
     {
         const int M = 100;
         double[,] qD = new double[tD.Length, M];
         double[] pD = new double[tD.Length];

         int tindex = 1;

         double[,] A = new double[M + 1, M + 1];
         double[] B = new double[M + 1];

         while (tindex < tD.Length)
         {
             // Заполнение матрицы A текущего шага

             for (int j = 1; j <= M; j++)
                 for (int m = 1; m <= M; m++)
                 {
                     double Xj = X(tD[tindex], tD[tindex - 1], j, m, M);
                     double Yj = Y(tD[tindex], tD[tindex - 1], j, m, M);
                     A[j - 1, m - 1] = Math.Sqrt(Math.PI) / 4.0 * (Xj + Yj);
                 }

             // Заполнение исторической матрицы B

             for (int j = 1; j <= M; j++)
             {
                 double Bi = 0;

                 for (int it = 1; it < tindex; ++it)
                 {
                     for (int m = 1; m <= M; m++)
                     {
                         double Xj = X(tD[tindex], tD[it - 1], j, m, M) - X(tD[tindex], tD[it], j, m, M);
                         double Yj = Y(tD[tindex], tD[it - 1], j, m, M) - Y(tD[tindex], tD[it], j, m, M);

                         Bi = Bi + qD[it, m - 1] * (Xj + Yj);
                     }
                 }
                 B[j - 1] = -Bi * Math.Sqrt(Math.PI) / 4.0;
             }


             // Заполнение матрицы для расчета изменения дебита

             for (int i = 0; i < M; i++)
             {
                 A[M, i] = +1; A[i, M] = -1;
             }

             A[M, M] = 0; B[M] = M;

             var Am = MathNet.Numerics.LinearAlgebra.Matrix.Build.DenseOfArray(A);
             var Bm = MathNet.Numerics.LinearAlgebra.Vector.Build.Dense(B);
             var Cm = Am.Solve(Bm);

             for (int m = 1; m <= M; m++)
             {
                 qD[tindex, m - 1] = Cm[m - 1];
             }
             pD[tindex] = Cm[M];
             tindex++;
         }
     }
 }

***

Восстановление порядка вывода итоговой формулы стало возможным благодаря изданным позже следующим статьям,

SPE-4051 «Unsteady-State pressure distribution created by a well with a single infinty-conductivity vertical fracture» (Aug 1974).
SPE-8281-PA «Non-Darcy Flow in Wells With Finite-Conductivity Vertical Fractures» 1982
SPE-9344-PA «Effect of Non-Darcy Flow on the Constant-Pressure Production of Fractured Wells» 1981
SPE-104004-MS «Effect of Pressure in a Well With a Vertical Fracture With Variable Conductivity and Skin Fracture» 2006

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Архивы

Яндекс.Метрика