Расчет скважины с МГРП

Собственно, ничего не мешает рассмотреть случай $n$ трещин, расположенных параллельно друг другу, вскрывающих пласт на полную высоту и обладающие бесконечной проводимостью.

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

$\Delta p(x,y,t) = \frac{q_f}{\phi c_t} \, \int\limits_0^t \frac{1}{4 \sqrt{\eta \pi \tau}} \left [ \mathtt{erf} \, \frac{y - (y_w - y_f/2)}{2 \sqrt{\eta \tau}} - \mathtt{erf} \, \frac{y - (y_w + y_f/2)}{2 \sqrt{\eta \tau}} \right ] \cdot \exp \left [ -\frac{(x-x_w)^2}{4 \eta \tau} \right ] d\tau$

Удельный дебит на площадь трещины, выражается через дебит трещины,

$q_f = \frac{q \cdot B}{y_f \cdot h}$

Суперпозиция позволяет получить падение давления на забое первой трещины, как сумму падений давления трех трещин. Чтобы не загромождать и так громоздкие формулы, номер трещины вынес перед переменной,

$\Delta \textbf{1}p(x,y,t) =$

$\frac{\textbf{1}q_{f}}{\phi c_t} \, \int\limits_0^t \frac{1}{4 \sqrt{\eta \pi \tau}} \left [ \mathtt{erf} \, \frac{\textbf{1}y_w - (\textbf{1}y_w - \textbf{1}y_f/2)}{2 \sqrt{\eta \tau}} - \mathtt{erf} \, \frac{\textbf{1}y_w - (\textbf{1}y_w + \textbf{1}y_f/2)}{2 \sqrt{\eta \tau}} \right ] \cdot \exp \left [ -\frac{(\textbf{1}x_w-\textbf{1}x_w)^2}{4 \eta \tau} \right ] d\tau +$

$\frac{\textbf{2}q_{f}}{\phi c_t} \, \int\limits_0^t \frac{1}{4 \sqrt{\eta \pi \tau}} \left [ \mathtt{erf} \, \frac{\textbf{1}y_w - (\textbf{2}y_w - \textbf{2}y_f/2)}{2 \sqrt{\eta \tau}} - \mathtt{erf} \, \frac{\textbf{1}y_w - (\textbf{2}y_w + \textbf{2}y_f/2)}{2 \sqrt{\eta \tau}} \right ] \cdot \exp \left [ -\frac{(\textbf{1}x_w-\textbf{2}x_w)^2}{4 \eta \tau} \right ] d\tau +$

$\frac{\textbf{3}q_{f}}{\phi c_t} \, \int\limits_0^t \frac{1}{4 \sqrt{\eta \pi \tau}} \left [ \mathtt{erf} \, \frac{\textbf{1}y_w - (\textbf{3}y_w - \textbf{3}y_f/2)}{2 \sqrt{\eta \tau}} - \mathtt{erf} \, \frac{\textbf{1}y_w - (\textbf{3}y_w + \textbf{3}y_f/2)}{2 \sqrt{\eta \tau}} \right ] \cdot \exp \left [ -\frac{(\textbf{1}x_w-\textbf{3}x_w)^2}{4 \eta \tau} \right ] d\tau$

или заменяя переменными,

$\Delta \textbf{1}p(x,y,t) = \textbf{1}q \, a_{11} + \textbf{2}q \, a_{12} + \textbf{3}q \, a_{13}$

Коэффициенты $a$ зависят от взаимного положения трещин и времени и представляют собой величину обратную продуктивности. Индексация коэффициентов, например «12», отвечает за влияние на работу первой трещины второй трещины и так далее.

Записывая те же выражения для второй и третьей трещины,

$\Delta \textbf{2}p(x,y,t) = \textbf{1}q \, a_{21} + \textbf{2}q \, a_{22} + \textbf{3}q \, a_{23}$

$\Delta \textbf{3}p(x,y,t) = \textbf{1}q \, a_{31} + \textbf{2}q \, a_{32} + \textbf{3}q \, a_{33}$

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

$p_w = \textbf{1}p_w + \textbf{2}p_w + \textbf{3}p_w$

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

Дебит скважины представляет собой сумму дебитов трещин.

$q = \textbf{1}q + \textbf{2}q + \textbf{3}q$

что позволяет составить систему уравнений,

$\begin{equation*} \begin{cases} q_1 \cdot a_{11} + q_2 \cdot a_{12} + q_3 \cdot a_{13} = \Delta p_1 \\ q_1 \cdot a_{21} + q_2 \cdot a_{22} + q_3 \cdot a_{23} = \Delta p_2 \\ q_1 \cdot a_{31} + q_2 \cdot a_{32} + q_3 \cdot a_{33} = \Delta p_3 \\ q = q_1 + q_2 + q_3 \end{cases} \end{equation*}$

которую подготовим к матричному представлению,

$\begin{equation*} \begin{cases} q_1 \cdot a_{11} + q_2 \cdot a_{12} + q_3 \cdot a_{13} - \Delta p = 0 \\ q_1 \cdot a_{21} + q_2 \cdot a_{22} + q_3 \cdot a_{23} - \Delta p = 0 \\ q_1 \cdot a_{31} + q_2 \cdot a_{32} + q_3 \cdot a_{33} - \Delta p = 0\\ q_1 \cdot 1 + q_2 \cdot 1 + q_3\cdot 1 - 0 = q \end{cases} \end{equation*}$

$\begin{equation*} \left[ \begin{array}{cccc} a_{11} & a_{12} & a_{13} &-1 \\ a_{21} & a_{22} & a_{23} &-1\\ a_{31} & a_{32} & a_{33} &-1 \\ 1 & 1 & 1 & 0 \\ \end{array} \right] \begin{bmatrix} q_1 \\ q_2 \\ q_3 \\ \Delta p \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ q \end{bmatrix} \] \end{equation*}$

Полученная система уравнений решается методом Гаусса (Gauss–Jordan elimination).

Пробный расчет для пласта низкой проницаемости,

m = 0.15;
ct = 1e-09;
mu = 0.328 * 0.001;
k = 1 * 0.986923 * 1e-15;
h = 8.5;
B = 1.019;

Задал три стадии ГРП и дебит скважины 100 m3/d.

xw = 100, yw = 100, xf = 70 m
xw = 100, yw = 150, xf = 50 m
xw = 100, yw = 300, xf = 20 m

Решение системы дает следующее,

q1 = 49.9 m3/d
q2 = 35.7 m3/d
q3 = 14.4 m3/d

Коэффициенты продуктивности трещин изменяются во времени неравномерно, что приводит к появлению различия забойных давлений между трещинами.

Для выравнивания забойных давлений нужно изменить дебит трещины. Согласно принципу суперпозиции, изменение дебита каждой трещины, эквивалентно появлению источника нагнетания или отбора. В итоге на поздних временах, забойное давление складывается из всех прошлых источников.

Следующая система уравнений позволяет найти изменение дебита каждой трещины и выравненое забойное давление скважины.

$\begin{equation*} \begin{cases} \Delta q_1 \cdot a_{11} + \Delta q_2 \cdot a_{12} + \Delta q_3 \cdot a_{13} = p_i - \sum \Delta p - p_w \\ \Delta q_1 \cdot a_{21} + \Delta q_2 \cdot a_{22} + \Delta q_3 \cdot a_{23} = p_i - \sum \Delta p - p_w \\ \Delta q_1 \cdot a_{31} + \Delta q_2 \cdot a_{32} + \Delta q_3 \cdot a_{33} = p_i - \sum \Delta p - p_w \\ 0 = \Delta q_1 + \Delta q_2 + \Delta q_3 \end{cases} \end{equation*}$

Сравнение результатов полуаналитического расчета проводится с моделью созданной в tNavigator.

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

TNAVCTRL
FRACTURE_BUILD_LOGIC UNSTRUCTURED_GRID /
UNSTRUCTURED_FRACTURE_GRID_LEVELS 9 /
/

Сравнение изменения дебита трещин во времени с расчетом симулятора.

В начале расчета симулятор дает численный дребезг — черные точки отлетают от аналитического расчёта. После отметки в 100 часов распределение дебитов стабилизируется.

Начиная с некоторой большой величины вложенности, расчет начинается рассыпаться, забойное давление может неожиданно вырасти в несколько раз, что можно заметить оценивая зависимость динамики забойного давления от размера ячейки. В итоге достигается эмпирический предел возможности симулятора провести расчет, однако являются ли полученные результаты валидными могут ответить полуаналитические модели. А для отладки алгоритма расчета полуаналитической модели, хорошо бы иметь примеры расчета симулятора.

Вот так рука об руку происходило развитие аналитического фундамента расчета МГРП.

Сравнение динамики забойного давления.

Результаты определенно демонстрирует схожую динамику. Разница между расчетами вызвана тем, что в аналитической модели предполагается равномерный приток к плоскости трещины. В статье 1974 года Alan Gringarten предлагает добавить псевдо-скин или другими словами дополнительный перепад давления, мы же сделаем переход к модели трещины конечной проводимости, где этот недостаток будет устранен.

Есть и второй вопрос. Численное интегрирование интеграла при определении значения коэффициента $a$ , очень не стабильно при малых временах $\Delta t \to 0$ . Рекомендации о решении интеграла приводятся в нескольких работах, мы будем придерживаться работ Erdan Ozkan, который с помощью преобразования Лапласа расширил возможность решения задачи притока к трещинам ГРП на все известные случаи граничных условий на скважине и внешней границы.

Его диссертация есть логическое продолжение работы 1949 года Application laplace transformation…, изучив которую я попробую улучшить качество расчета модели.

Исходные тексты расчетов выложены в открытом доступе на https://github.com/romanio/wRate.

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

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

Комментарии

Архивы