Application laplace transformation to flow problems in reservoirs (1949)

В 1949 году выходит статья «О применении преобразования Лапласа для решения задач течения флюида в пластах» за авторством A. F. Van Everdingen (Shell) и W. Hurst, которая демонстрирует технику применения операционного исчисления для решения уравнения диффузии.

Статья ссылается на фундаментальную книгу «Conduction of Heat in Solids» H.S.Carslaw и J.C.Jaeger (1947), посвященную решению задач теплопроводности твердых тел. Как упоминается в книге, дифференциальное уравнение течения сжимаемой жидкости через пористую среду полностью соответствует уравнению теплопроводности, поэтому неудивительно, что методы решения были перенесены и на задачи, возникающие при разработке месторождений.

Раздел посвященный переходу к безразмерным координатам выписан из приложения книги «Pressure Transient Testing», J.Lee, J.B. Rollins и J.P.Spivey, также для операции с функциями Бесселя, как и у авторов, используется книга «Mathematical Methods in Engineering», T.Karman и M.A.Biot, а также «Treatise on the theory of bessel functions», G.N.Watson (1945).

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

***

Родственные теплопроводности дифференциальные уравнения обобщено называют уравнением диффузии, которые в цилиндрических (радиальных) координатах записываются в виде,

$\mathtt{ \frac{\partial^2 p}{\partial r^2} + \frac{1}{r} \frac{\partial p}{\partial r} = \frac{1}{\eta} \cdot \frac{\partial p}{\partial t} }$

где коэффициент диффузии,

$\mathtt{ {\eta} = \frac{k}{\phi \cdot \mu \cdot c_t} }$

Уравнение диффузии традиционно решается при помощи рядов Фурье-Бесселя, однако используя преобразование Лапласа можно получить более элегантное решение. Умножив обе части уравнения на $\mathtt{e^{-st}}$ и проинтегрировав по $\mathtt{t}$ , получим изображение уравнения диффузии,

$\mathtt{ \int_0^\infty e^{-st} \, \frac{\partial^2 p}{\partial r^2} \, dt + \frac{1}{r} \, \int_0^\infty e^{-st} \, \frac{\partial p}{\partial r} \, dt = \frac{1}{\eta} \, \int_0^\infty e^{-st} \, \frac{\partial p}{\partial t} \, dt }$

Правая часть уравнения представляет собой преобразование производной функции,

$\mathtt{ \int_0^\infty e^{-st} \, \frac{\partial p}{\partial t} \, dt = s\cdot P(s) - p(0) }$

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

Рассмотрим перепад давления и его частные производные,

$\mathtt{ \Delta p(r, t) = p_i - p(r, t) } }$

$\mathtt{ \frac{\partial \Delta p(r, t)}{\partial t} = -\frac{\partial p(r, t)}{\partial t} }$

$\mathtt{ \frac{\partial \Delta p(r, t)}{\partial r} = -\frac{\partial p(r, t)}{\partial r} }$

Заменяя в уравнении диффузии давление на величину перепада давления, получим такое же по форме уравнение диффузии,

$\mathtt{ \frac{\partial^2 \Delta p}{\partial r^2} + \frac{1}{r} \frac{\partial \Delta p}{\partial r} = \frac{1}{\eta} \cdot \frac{\partial \Delta p}{\partial t} }$

Применяя подобные преобразования стараются перейти к удобным для вычислений значениям в новых переменных. Например, заменив давление на перепад давления и учитывая, что давление в нулевой момент времени равняется начальному пластовому давлению $\mathtt{p(r,0) = p_i}$ , получим что,

$\mathtt{ \int_0^\infty e^{-st} \, \frac{\partial \Delta p}{\partial t} \, dt = s\cdot P(s) - \Delta p(0) = s\cdot P(s) }$

Далее знак дельты будет опускаться и под давлением всегда понимается разница текущего давления от начального пластового.

Также нам понадобится масштабный множитель для радиуса, который выбирается так, чтобы при $\mathtt{r = r_w}$ безразмерное расстояние равнялось единице,

$\mathtt{ r_D = \frac{r}{r_w} }$

$\mathtt{ \frac{1}{r_w^2} \cdot \frac{\partial^2 p}{\partial r_D^2} + \frac{1}{r_w^2 \cdot r_D} \cdot \frac{\partial p}{\partial r_D} = \frac{1}{\eta} \cdot \frac{\partial p}{\partial t} }$

$\mathtt{ \frac{\partial^2 p}{\partial r_D^2} + \frac{1}{r_D} \cdot \frac{\partial p}{\partial r_D} = \frac{r_w^2}{\eta} \cdot \frac{\partial p}{\partial t} }$

Множитель справа, перед производной по времени, будет мешать в дальнейшем получить уравнение Бесселя, поэтому спрячем его под производную давления,

$\mathtt{ \partial \left (\frac{t \cdot \eta}{r_w^2} \right ) = \partial t_D }$

Делая таким образом переход к безразмерному времени,

$\mathtt{t_D = t \cdot \frac{\eta}{r_w^2} }$

Уравнение диффузии с двумя новыми безразмерными переменными примет следующий вид,

$\mathtt{ \frac{\partial^2 p}{\partial r_D^2} + \frac{1}{r_D} \cdot \frac{\partial p}{\partial r_D} = \frac{\partial p}{\partial t_D} }$

Для левой части изображения уравнения диффузии вспомним формулу Лейбница,

$\mathtt{ u(r) = \int_0^\infty f(t, r) \, dt }$

$\mathtt{ \frac{d u(r)}{dr} = \frac{d}{dr} \left [ \int_0^\infty f(t, r) \right ] dt = \int_0^\infty \frac{\partial f(t, r)}{\partial r} dt }$

По определению, изображение функции давления,

$\mathtt{ P(s) = \int_0^\infty e^{-st} \cdot p(t, r) \, dt }$

Первая производная запишется как,

$\mathtt{ \frac{d P(s)}{dr} = \int_0^\infty e^{-st} \cdot \frac{\partial p(t, r)}{\partial r} dt }$

Вторая производная,

$\mathtt{\frac{d^2 P(s)}{dr^2} = \int_0^\infty e^{-st} \cdot \frac{\partial^2 p(t, r)}{\partial r^2} dt }$

Преобразование Лапласа свело дифференциальное уравнение диффузии с частными производными к обыкновенному дифференциальному уравнению, которое называется вспомогательным уравнением,

$\mathtt{ \frac{d^2 P(s)}{dr_D^2} + \frac{1}{r_D} \, \frac{d P(s)}{dr_D} = s\cdot P(s) }$

Совершая очередную замену,

$\mathtt{x = \sqrt{s}\cdot r_D }$

$\mathtt{\frac{d^2 P(s)}{dr_D^2} = s \cdot \frac{d^2 P(s)}{dx^2} }$

$\mathtt{\frac{d P(s)}{dr_D} = \sqrt{s} \cdot \frac{d P(s)}{dx} }$

Получим модифицированное уравнение Бесселя нулевого порядка,

$\mathtt{P'' + \frac{1}{x} \cdot P' - P = 0}$

Решение которого записывается в виде

$\mathtt{P(r_D, s) = A \cdot I_0 \left (\sqrt{s}\cdot r_D ) + B \cdot K_0 (\sqrt{s}\cdot r_D )}$

Для определения коэффициентов $\mathtt{A}$ и $\mathtt{B}$ , потребуется два граничных условия, заданных на забое скважины или на контуре питания.

Решение для бесконечного пласта

Известно, что для модифицированных функций Бесселя при стремлении аргумента $\mathtt{\sqrt{s}\cdot r_D \to \infty}$ , значение $\mathtt{I_0}$ становится очень большим, а значение $\mathtt{K_0}$ стремится к нулю. В начальный момент времени во всех точках пласта перепад давления равен нулю и изображение перепада давления также равно нулю, поэтому при любых $\mathtt{r_D}$ должно выполнятся,

$\mathtt{0 = A \cdot I_0 \left (\sqrt{s}\cdot r_D ) + B \cdot K_0 (\sqrt{s}\cdot r_D )}$

Равенство нулю возможно только при $\mathtt{A = 0}$ , поэтому в общем виде решение для бесконечного пласта следующее,

$\mathtt{P(r_D, s) = B \cdot K_0 (\sqrt{s}\cdot r_D )}$

Конкретный вид будет зависеть уже от условия на забое скважины, при $\mathtt{r_D = 1}$ .

Первым рассмотрим режим поддержания постоянного дебита скважины.

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

$\mathtt{q \cdot B = \upsilon \cdot A = \frac{k}{\mu} \frac{\partial p}{\partial r} \cdot 2 \pi r h }$

$\mathtt{q = 2 \pi \cdot \frac{kh}{B \mu} \cdot \left (r \frac{\partial p}{\partial r} \right )}$

Вспоминая, что в уравнении диффузии мы перешли от давления к перепаду давлений,

$\mathtt{ \frac{\partial \Delta p(r, t)}{\partial r} = -\frac{\partial p(r, t)}{\partial r} }$

Запишем уравнение Дюпюи иначе,

$\mathtt{q = 2 \pi \cdot \frac{kh}{B \mu} \cdot \left (-r \frac{\partial \Delta p}{\partial r} \right )}$

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

$\mathtt{ \left (-r \frac{\partial \Delta p}{\partial r} \right )_{r = r_w} = q \cdot \frac{B \mu}{ 2 \pi kh } }$

$\mathtt{ \left (-r_D \frac{\partial \Delta p}{\partial r_D} \right )_{r_D = 1} = q \cdot \frac{B \mu}{ 2 \pi kh } }$

Здесь безразмерный радиус помогает избавится от радиуса скважины,

$\mathtt{ \left ( \frac{\partial \Delta p}{\partial r_D} \right )_{r_D = 1} = -q \cdot \frac{B \mu}{ 2 \pi kh } }$

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

$\mathtt{ \frac{d P(1, s)}{dr_D} = -q \cdot \frac{B \mu}{ 2 \pi kh} \cdot \frac{1}{s} }$

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

$\mathtt{ K_0'(z) = -K_1(z) }$

$\mathtt{ \frac{\partial P(r_D, s)}{\partial r_D} = \frac{\partial}{\partial r_D} \left ( B \cdot K_0 (\sqrt{s}\cdot r_D ) \right ) = -B \cdot \sqrt{s} \cdot K_1(\sqrt{s} \cdot r_D ) }$

Зная, чему равно значение производной из заданного дебита скважины при $\mathtt{r_D = 1}$ ,

$\mathtt{ q \cdot \frac{B \mu}{ 2 \pi kh} \cdot \frac{1}{s} = B \cdot \sqrt{s} \cdot K_1(\sqrt{s} ) }$

$\mathtt{ B = q \cdot \frac{B \mu}{ 2 \pi kh} \cdot \frac{1}{s \cdot \sqrt{s} \cdot K_1(\sqrt{s} ) } }$

Подставляя найденный коэффициент в общее решение, получим изображение изменения перепада давления при постоянном дебите для бесконечного пласта,

$\mathtt{P(r_D, s) = q \cdot \frac{B \mu}{ 2 \pi kh} \cdot \frac{K_0 (\sqrt{s}\cdot r_D) }{s^{3/2} \cdot K_1(\sqrt{s})} }$

Для изображения перепада давления на забое скважины,

$\mathtt{P(1, s) = q \cdot \frac{B \mu}{ 2 \pi kh} \cdot \frac{K_0 (\sqrt{s}) }{s^{3/2} \cdot K_1(\sqrt{s})} }$

можно заметить, что форма решения не зависит от свойств пласта и с каким дебитом работает скважина. Поэтому разумно получить единственное общее решение в виде,

$\mathtt{ \frac{P(1, s)}{q \cdot \frac{B \mu}{ 2 \pi kh}} = \frac{K_0 (\sqrt{s}) }{s^{3/2} \cdot K_1(\sqrt{s})} }$

Ещё раз вспоминая определение преобразования Лапласа,

$\mathtt{ P(s) = \int_0^\infty e^{-st} \cdot p \, dt }$

Учтем множитель, сформировав новую безразмерную переменную,

$\mathtt{ \frac{P(r_D, s)}{q \cdot \frac{B \mu}{ 2 \pi kh}} = \int_0^\infty e^{-st} \cdot \frac{\Delta p}{q \cdot \frac{B \mu}{ 2 \pi kh}} \, dt }$

$\mathtt{ p_D(r_D,t_D) = 2 \pi \cdot \frac{kh}{q B \mu } \cdot (p_i - p(r_D, t_D) )}$

Уравнение диффузии теперь с учетом трех безразмерных переменных принимает свою самую известную форму,

$\mathtt{ \frac{\partial^2 p_D}{\partial r_D^2} + \frac{1}{r_D} \cdot \frac{\partial p_D}{\partial r_D} = \frac{\partial p_D}{\partial t_D} }$

Внутреннее граничное условие упрощается,

$\mathtt{ \left ( \frac{\partial p_D}{\partial r_D} \right )_{r_D = 1} = -1 }$

Решение для бесконечного пласта при режиме поддержания постоянного дебита скважины,

$\mathtt{P(1, s) = \frac{K_0 (\sqrt{s}) }{s^{3/2} \cdot K_1(\sqrt{s})} }$

В этой заметке для получения обратного преобразования Лапласа мы будем использовать несложный в реализации численный алгоритм Stehfest.

Полученное решение представляет скважину как полый цилиндр, на границе которого создается забойное давление. Известно также решение уравнении диффузии методом точечного источника, при котором скважина представлена бесконечно тонкой линией нулевого радиуса (line-source well) $\mathtt{r_w \to 0}$ .

$\mathtt{p_D(r_D, t_D) = -\frac{1}{2} \cdot Ei \cdot \left[ -\frac{r_D^2}{4t_D} \right ] }$

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

Для безразмерного времени $\mathtt{t_D > 10}$ , решения дают одинаковые результаты.

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

Для режима поддержания постоянного забойного давления, перепад давления на забое скважины равен $\mathtt{\Delta p = p_i-p_w = const}$ . Найдем преобразование Лапласа,

$\mathtt{L[(p_i-p_w)] = (p_i-p_w) \cdot \frac{1}{s} }$

Теперь можно выразить коэффициент $\mathtt{B}$ ,

$\mathtt{P(1,r_D) = (p_i-p_w) \cdot \frac{1}{s} = B \cdot K_0 (\sqrt{s})}$

и получить преобразование Лапласа для распределения давления по всему пласту при поддержании $\mathtt{p_w = const}$ ,

$\mathtt{P(r_D,t_D) = (p_i-p_w) \cdot \frac{K_0(\sqrt{s}\cdot r_D )}{s \cdot K_0(\sqrt{s})} }$

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

$\mathtt{ \left ( \frac{\partial P(r_D,t_D)}{\partial r_D} \right)_{r_D = 1} = \frac{(p_i-p_w)}{s \cdot K_0(\sqrt{s})} \cdot \frac{\partial }{\partial r_D} \left [K_0(\sqrt{s}\cdot r_D) \right ] = - (p_i-p_w) \cdot \frac{K_1(\sqrt{s})}{\sqrt{s} \cdot K_0(\sqrt{s})} }$

Повторяя размышления о безразмерном давлении, постоянная разница давлений вносится под интеграл преобразования Лапласа. В итоге безразмерное давление записывается как,

$\mathtt{p_D = \frac{p_i - p}{p_i-p_w} }$

Уравнение диффузии с тремя безразмерными параметрами записывается также как и раньше,

$\mathtt{ \frac{\partial^2 p_D}{\partial r_D^2} + \frac{1}{r_D} \cdot \frac{\partial p_D}{\partial r_D} = \frac{\partial p_D}{\partial t_D} }$

Внутреннее граничное условие упрощается,

$\mathtt{ \left ( \frac{\partial p_D}{\partial r_D} \right )_{r_D = 1} = 1 }$

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

$\mathtt{ \left ( \frac{\partial p}{\partial r_D} \right )_{r_D = 1} = q \cdot \frac{B \mu}{ 2 \pi kh } }$

Найдем производную безразмерного давления,

$\mathtt{ \frac{\partial p_D}{\partial r_D} = -\frac{1}{p_i-p_w} \frac{\partial p}{\partial r_D}}$

Следовательно,

$\mathtt{ \left ( (p_i-p_w) \cdot \frac{\partial p_D}{\partial r_D} \right )_{r_D = 1} = -q \cdot \frac{B \mu}{ 2 \pi kh } }$

Сформируем комплекс безразмерного дебита скважины,

$\mathtt{ q_D = q \cdot \frac{B \mu}{ 2 \pi kh } \cdot \frac{1}{p_i - p_w}}$

Найдем изображение безразмерного дебита,

$\mathtt{ L[q_D](s) = - \left (\frac{\partial P}{\partial r_D} \right )_{r_D = 1} = \frac{K_1(\sqrt{s})}{\sqrt{s} \cdot K_0(\sqrt{s})}}$

Численное решение обратного преобразования Лапласа позволяет найти изменение безразмерного дебита скважины,

Таким образом при сохранении постоянного забойного давления, темпы падения дебита в безразмерных координатах имеют единый вид для всех вертикальных скважин .

Выразим накопленную добычу,

$\mathtt{Q = \int_0^t 2 \pi \cdot \frac{kh}{B \mu} \cdot \left ( r \frac{\partial p}{\partial r} \right )_{r=r_w} \, dt} =$

$\mathtt{ = 2 \pi \cdot \frac{kh}{B \mu} \cdot \frac{\phi \mu c_t \cdot r_w^2 }{k} \cdot \int_0^{t_D} \left ( -(p_i-p_w)\cdot\frac{\partial p_D}{\partial r_D} \right )_{r_D=1} \, dt_D}$

Сформируем комплекс безразмерной накопленной добычи,

$\mathtt{ Q_D = Q \cdot 2 \pi \cdot \frac{kh}{B \mu} \cdot \frac{\phi \mu c_t \cdot r_w^2 }{k} \cdot (p_i-p_w) }$

$\mathtt{ Q_D = -\int_0^{t_D} \left (\frac{\partial p_D}{\partial r_D} \right )_{r_D=1} \, dt_D}$

Найдем преобразование Лапласа,

$\mathtt{ L[Q_D] = -\int_0^\infty e^{-st} \left [ \int_0^{t_D} \left (\frac{\partial p_D}{\partial r_D} \right )_{r_D=1} \, dt_D} \right ] dt =$

$\mathtt{ = -\int_0^\infty \left [ \int_0^{t_D} \left (\frac{\partial p_D}{\partial r_D} \right )_{r_D=1} \, dt_D} \right ] d \left(-\frac{1}{s} \cdot e^{-st} \right ) =$

применим метод интегрирования по частям,

$\mathtt{ \int u \, dv = uv - \int v \, du }$

$\mathtt{ = \left( \frac{1}{s}e^{-st} \int_0^{t_D} \left (\frac{\partial p_D}{\partial r_D} \right )_{r_D=1} \, dt_D} \right ) \Bigg|_0^\infty - \frac{1}{s} \int_0^\infty e^{-st} d \left( \int_0^{t_D} \left (\frac{\partial p_D}{\partial r_D} \right )_{r_D=1} \, dt_D} \right ) = }$

левый интеграл при границе $\mathtt{t\to\infty}$ равен нулю, а при $\mathtt{t\to 0}$ накопленная добыча в нулевой момент равен нулю, поэтому

$\mathtt{ L[Q_D] = -\frac{1}{s} \cdot \frac{\partial P}{\partial r_D }}$

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

$\mathtt{ L[Q_D]= \frac{K_1(\sqrt{s})}{s^{3/2} \cdot K_0(\sqrt{s})} }$

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

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

Решение для замкнутого пласта

Замкнутый пласт представлен условием на внешней границе как отсутствие притока,

$\mathtt{\left(\frac{\partial p_D}{\partial r_D} \right )_{r_D = R} = 0 }$

Из общего решения,

$\mathtt{P(r_D, s) = A \cdot I_0 \left (\sqrt{s}\cdot r_D ) + B \cdot K_0 (\sqrt{s}\cdot r_D )}$

используя свойство $\mathtt{K_0'(z)=-K_1(z)}$ и $\mathtt{I_0'(z)= I_1(z)}$ найдем производную по расстоянию,

$\mathtt{ \frac{\partial P}{\partial r_D} = A \sqrt{s} \cdot I_1 \left (\sqrt{s}\cdot r_D) - B \sqrt{s} \cdot K_1 (\sqrt{s}\cdot r_D )}$

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

$\mathtt{L\left[ \left(\frac{\partial p_D}{\partial r_D} \right )_{r_D = 1} = - 1\right] = -\frac{1}{s} }$

$\mathtt{L\left[ \left(\frac{\partial p_D}{\partial r_D} \right )_{r_D = R} = 0\right] = 0 }$

которые позволяют выразить коэффициенты $\mathtt{A}$ и $\mathtt{B}$ и записать общий вид уравнения,

$\mathtt{ P(r_D,s) = \frac{1}{s^{3/2}} \cdot \frac{K_1(\sqrt{s}\cdot R_D) \cdot I_0(\sqrt{s} \cdot r_D) + I_1 (\sqrt{s} \cdot R_D) \cdot K_0(\sqrt{s}\cdot r_D)}{I_1(\sqrt{s} \cdot R_D) \cdot K_1 (\sqrt{s}) - I_1(\sqrt{s})\cdot K_1(\sqrt{s} \cdot R_D) } }$

Получим серию решений для различных значений безразмерного радиуса ограниченной залежи $\mathtt{R_D = R_e / r_w}$ ,

Для режима поддержания постоянного забойного давления, безразмерный постоянный перепад давления на скважине равен единице,

$\mathtt{L[p_D = 1] = \frac{1}{s} }$

На внешней границе,

$\mathtt{L\left[ \left(\frac{\partial p_D}{\partial r_D} \right )_{r_D = R} = 0\right] = 0 }$

Из общего решения,

$\mathtt{P(r, s) = A \cdot I_0 \left (\sqrt{s}\cdot r ) + B \cdot K_0 (\sqrt{s}\cdot r )}$

можно записать систему из двух уравнений,

$\mathtt{ \frac{1}{s} = A \cdot I_0 \left (\sqrt{s} ) + B \cdot K_0 (\sqrt{s})}$

используя свойство $\mathtt{K_0'(z)=-K_1(z)}$ и $\mathtt{I_0'(z)= I_1(z)}$

$\mathtt{ \frac{\partial P}{\partial r_D} = A \sqrt{s} \cdot I_1 \left (\sqrt{s}\cdot r_D) - B \sqrt{s} \cdot K_1 (\sqrt{s}\cdot r_D )}$

$\mathtt{ 0 = A \cdot I_1 \left (\sqrt{s}\cdot R_D) - B \cdot K_1 (\sqrt{s}\cdot R_D)}$

Получим следующее общее решение,

$\mathtt{P(r_D, s) = \frac{1}{s} \cdot \frac{K_1(\sqrt{s} \cdot R_D) \cdot I_0(\sqrt{s} \cdot r_D) + I_1 (\sqrt{s}\cdot R_D) \cdot K_0(\sqrt{s}\cdot r_D)}{I_0(\sqrt{s}) \cdot K_1 (\sqrt{s}\cdot R_D ) + I_1(\sqrt{s} \cdot R_D)\cdot K_0(\sqrt{s}) } }$

Чтобы найти безразмерный дебит и накопленную добычу, найдем производную давления,

$\mathtt{-\frac{\partial P}{\partial r_D} = \frac{1}{s^{1/2}} \cdot \frac{K_1(\sqrt{s}) \cdot I_1(\sqrt{s} \cdot R_D) - I_1 (\sqrt{s}) \cdot K_1(\sqrt{s}\cdot R_D)}{I_0(\sqrt{s}) \cdot K_1 (\sqrt{s}\cdot R_D ) + I_1(\sqrt{s} \cdot R_D)\cdot K_0(\sqrt{s}) } }$

Ранее мы использовали переход от производной давления,

$\mathtt{ L[q_D](s) = - \left (\frac{\partial P}{\partial r_D} \right )_{r_D = 1}}$

$\mathtt{ L[Q_D] = -\frac{1}{s} \cdot \left( \frac{\partial P}{\partial r_D } \right) _{r_D=1} }$

Получим серию решений для различных значений безразмерного радиуса ограниченной залежи $\mathtt{R_D = R_e / r_w}$ ,

$\mathtt{ L[q_D] = \frac{1}{s^{1/2}} \cdot \frac{K_1(\sqrt{s}) \cdot I_1(\sqrt{s} \cdot R_D) - I_1 (\sqrt{s}) \cdot K_1(\sqrt{s}\cdot R_D)}{I_0(\sqrt{s}) \cdot K_1 (\sqrt{s}\cdot R_D ) + I_1(\sqrt{s} \cdot R_D)\cdot K_0(\sqrt{s}) } }$

$\mathtt{ L[Q_D] = \frac{1}{s^{3/2}} \cdot \frac{K_1(\sqrt{s}) \cdot I_1(\sqrt{s} \cdot R_D) - I_1 (\sqrt{s}) \cdot K_1(\sqrt{s}\cdot R_D)}{I_0(\sqrt{s}) \cdot K_1 (\sqrt{s}\cdot R_D ) + I_1(\sqrt{s} \cdot R_D)\cdot K_0(\sqrt{s}) } }$

Решение для поддержания начального пластового давления на границе контура

Поддержание постоянного пластового давления на границе представляется как отсутствие перепада давления на внешней границе,

$\mathtt{ (\Delta p)_{r_D = R} = 0 }$

Из общего решения,

$\mathtt{P(r_D, s) = A \cdot I_0 \left (\sqrt{s}\cdot r_D ) + B \cdot K_0 (\sqrt{s}\cdot r_D )}$

получим,

$\mathtt{0 = A \cdot I_0 \left (\sqrt{s}\cdot R_D ) + B \cdot K_0 (\sqrt{s}\cdot R_D )}$

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

$\mathtt{P(r_D, s) = \frac{1}{s^{3/2}} \cdot \frac{K_0(\sqrt{s} \cdot r_D) \cdot I_0(\sqrt{s} \cdot R_D) - I_0 (\sqrt{s}\cdot r_D) \cdot K_0(\sqrt{s}\cdot R_D)}{I_1(\sqrt{s} \cdot r_D) \cdot K_0 (\sqrt{s}\cdot R_D ) + I_0(\sqrt{s} \cdot R_D)\cdot K_1(\sqrt{s} \cdot r_D) } }$

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

$\mathtt{p_D = ln(R_D)}$

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

$\mathtt{\frac{1}{s} = A \cdot I_0 \left (\sqrt{s}\cdot r_D ) + B \cdot K_0 (\sqrt{s}\cdot r_D )}$

$\mathtt{0 = A \cdot I_0 \left (\sqrt{s}\cdot R_D ) + B \cdot K_0 (\sqrt{s}\cdot R_D )}$

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

$\mathtt{ L[q_D] = \frac{1}{s^{1/2}} \cdot \frac{K_0(\sqrt{s} \cdot R_D) \cdot I_1(\sqrt{s}) + I_0 (\sqrt{s}\cdot R_D) \cdot K_1(\sqrt{s})}{I_0(\sqrt{s} \cdot R_D) \cdot K_0 (\sqrt{s}) - I_0(\sqrt{s})\cdot K_0(\sqrt{s} \cdot R_D) } }$

Накопленная добыча нефти,

$\mathtt{ L[Q_D] = \frac{1}{s^{3/2}} \cdot \frac{K_0(\sqrt{s} \cdot R_D) \cdot I_1(\sqrt{s}) + I_0 (\sqrt{s}\cdot R_D) \cdot K_1(\sqrt{s})}{I_0(\sqrt{s} \cdot R_D) \cdot K_0 (\sqrt{s}) - I_0(\sqrt{s})\cdot K_0(\sqrt{s} \cdot R_D) } }$

Исходные данные для графиков приведены в файле VanEverdingen

Алгоритм Stehfest

В 1970 году в ежемесячном журнале ассоциации вычислительной техники (ACM) под номером 368 опубликован алгоритм, позволяющий найти обратное преобразование Лапласа любой функции. Реализация на C++(Qt) получилась такой,

double factorial(int N) // factorial
{
double x = 1.0;
for (int i = 2; i <= N; i++) {
x = x * i;
}
return x;
}

QList StehfestCoefficients(int N)
{
int N2 = static_cast(N / 2);
int NV = 2 * N2;

QList H(NV);

int sign = 1;

if ((N2 % 2) != 0) {
sign = -1;
}

for (int i = 0; i < NV; ++i)
{
int kmin = (int) ((i + 2) / 2);
int kmax = i + 1;
if (kmax > N2) {
kmax = N2;
}

H[i] = 0;
sign = -sign;

for (int k = kmin; k <= kmax; ++k) {
H[i] += (pow(k, N2) / factorial(k)) * (factorial(2 * k) / factorial(2 * k - i - 1))
/ factorial(N2 - k) / factorial(k - 1) / factorial(i + 1 - k);
}

H[i] *= sign;
}

return H;
}

double StehfestAlgorithm(double t, double (*func)(double))
{
const double ln2 = 0.69314718056;
const int stehfest_steps = 18;
double ln2t = ln2 / t;

QList stehfest_coeff = StehfestCoefficients(stehfest_steps);
double s = 0.0;
double y = 0.0;

for (double c : stehfest_coeff) {
s += ln2t;
y += c * func(s);
}

return ln2t * y;
}

Для примера решим задачу определения безразмерного давления для ограниченного пласта при $\mathtt{q = const}$ . Сначала определим функцию изображения,


double Fs5(double s) // Ограниченный пласт, постоянный дебит скважины
{
    double R = 30;
    double sqs = sqrt(s);

    return (1 / pow(s, 3.0 / 2.0))
           * (std::cyl_bessel_k(1, sqs * R) * std::cyl_bessel_i(0, sqs)
              + std::cyl_bessel_i(1, sqs * R) * std::cyl_bessel_k(0, sqs))
           / (std::cyl_bessel_i(1, sqs * R) * std::cyl_bessel_k(1, sqs)
              - std::cyl_bessel_i(1, sqs) * std::cyl_bessel_k(1, sqs * R));
}

Затем для заданных $\mathtt{t_D}$ найдем обратные преобразования,


int main(int argc, char *argv[])
{
    double scale = 0.01;
    for (int M = 1; M < 7; M++) {
        for (int T = 1; T <= 9; T++) {
            qInfo() << "T = " << T * scale << " F(t)=" << StehfestAlgorithm(T * scale, Fs5);
        }
        scale = scale * 10;
    }
    return 0;
}

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

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

Комментарии

Архивы