...

Motivation

...

Definition

...

...

...

...

...

...

...

...

и с учетом:

...

One of the key problems in designing the pipelines and controlling the pipeline fluid transport is to predict the temperature and pressure losses during the stationary fluid transport.

Pipeline flow simulator is addressing this problem. It should account for the varying pipeline trajectory, gravity effects, fluid friction with pipeline walls and varying heat exchange with surroundings.

Definition

...

Given 

• space coordinates are 

  LaTeX Math Inline
  body \{ x, \, y, \, z \}

   with 

  LaTeX Math Inline
  body z

  -ccordinate facing down to the Earth Centre
  
  
• inflow pipeline coordinates 

  LaTeX Math Inline
  body \{ x_s = 0, \, y_s = 0, \, z_s = 0 \}

  
  
  
• pipeline trajectory 

  LaTeX Math Inline
  body \{ x_w(l), \, y_w(l), \, z_w(l) \}

  , where 

  LaTeX Math Inline
  body l = \int_0^l \sqrt{dx^2 + dy^2 + dz^2} = \int_0^l \sqrt{\dot x^2 + \dot y^2 + \dot z^2} dl

  ,  is pipeline length from inflow point 

  LaTeX Math Inline
  body \{ x_s = 0, \, y_s = 0, \, z_s = 0 \}

  
  
  
• pipeline cross-section area 

  LaTeX Math Inline
  body A(l)

  
  
  
• earth gravity vector 

  LaTeX Math Inline
  body {\bf g} = (0, \, 0, \, g)

   where 

  LaTeX Math Inline
  body g = 9.81 \ \rm m/s^2

  
  
  
• inflow temperature 

  LaTeX Math Inline
  body T_s

  , inflow pressure 

  LaTeX Math Inline
  body p_s

  , inflow rate 

  LaTeX Math Inline
  body q_s

  
  
  
• PVT properties of water 

  LaTeX Math Inline
  body \rho(T, p)

  , 

  LaTeX Math Inline
  body \mu(T, p)

  
  
  
• surroundings initial temperature  

  LaTeX Math Inline
  body T_g(l)

  , thermal diffusivity 

  LaTeX Math Inline
  body a_e(l)

  , thermal conductivity 

  LaTeX Math Inline
  body \lambda_e(l)

   of surrounding media
  
  
• heat exchange coefficient 

  LaTeX Math Inline
  body U(l)

   based on pipeline schematics
  
  

Simulate

• along-pipe distribution of stabilized pressure 

  LaTeX Math Inline
  body p(l)

  , flow rate

  LaTeX Math Inline
  body q(l)

   and average flow velocity 

  LaTeX Math Inline
  body u(l)

   
  
  
• along-pipe distribution of fluid flow temperature 

  LaTeX Math Inline
  body T(t, l)

   after a flow period of time 

  LaTeX Math Inline
  body t

   when the flow stabilization achieved
  
  

 A(l) \, \rho(l) \, v(l) = \rm const

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

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

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

...

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

Уравнение

 A(l) \, \rho(l) \, v(l) = \rm const

\frac{dp}{dl} = \rho \, g \, \sin \theta - \rho \, v \, \frac{dv}{dl} - \frac{ f \, \rho \, v^2 \, }{2 d}

\frac{dp}{dl} =  \bigg(rho \frac{dp}{dl} \bigg)_g + \bigg( \frac{dp}{dl} \bigg)_v + \bigg(, g \, \sin \theta - \rho \, v \, \frac{dpdv}{dl} \bigg)_f

где

- \frac{

 f \

, \rho \, 

v^2 \, 

}{2 d}

где

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

Уравнение

\bigg( \frac{dp}{dl} \bigg)_v = - \rho \, v \, \frac{dv}{dl},

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

\frac{dp}{dl} =  \bigg( \frac{dp}{dl} \bigg)_

g 

+ 

\bigg( \frac{

dp}{dl} \bigg)_v + \

bigg( \

frac{dp}{

фрикционная компонента вариации давления, формируемая трением флюида со стенкой скважины

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

dl} \bigg)_f

где

\bigg( \

Для несжимаемой жидкости

frac{dp}{dl} \bigg)_g = \rho \, g \, \sin \theta

,

и может быть явно проинтегрировано:

\bigg( \frac{dp}{dl} \bigg)_v = - \rho \, 

v \,

 \

frac{

dv}{

dl}

,

и называется уравнением Бернулли.

\bigg( \frac{

dp}{

dl} 

\bigg)_f = - \frac{

 f \, \rho v^2}{2 d},

Для несжимаемой жидкости

 \rho \, v \

откуда и следует формула

\frac{dp}{dl} = \rho \, g \, \sin \theta -

Для вывода уравнения

Рассмотрим стационарное (то есть установившееся во времени) течение потока по трубе.

Движение поперек трубы отсуствует и, следовательно, сумма проекций всех сил на трансверсальное направление

dF_p \bigg |_{l_{\perp}} + dF_g \bigg |_{l_{\perp}} + dF_f \bigg |_{l_{\perp}}+ dF_N \bigg |_{l_{\perp}} =0

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

Уравнение движения флюида вдоль оси трубы

 dF_p \bigg |_l + dF_g \bigg |_l + dF_f \bigg |_l+ dF_N \bigg |_l = \frac{d I}{dt}\bigg |_l

где

Изменение импульса c учетом стационарности скорости потока

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

Проекция гравитационной силы

Сила трения со стенками трубы дается феноменологическим уравнением Дарси-Вейсбаха:

Аксиальная компонента реакции опоры труб по определению отсутствует

Подставляя вышеприведенные выражения в уравнение

- A dp + \rho \, A\, \delta l \, g \, \sin \theta - \frac{f \, \rho \, v^2}{2 d} \, A \, \delta l = \rho \, A \, v \, dv

Разделив уравнение на бесконечно малый объем элемента

Если дебит скважины на устье составляет

A \, \rho \, v = \rho_s \, q_s

откуда можно выразить явно профиль скорости потока по стволу:

v(l) = \frac{\rho_s \, q_s}{\rho(p) \, A(l)}

Подставляя 

, \frac{dv}{dl}

и может быть явно проинтегрировано:

p(l) - \rho \, g \, l \, \sin \theta +  \frac{1}{2} \rho \, v^2 = \rm const

и называется уравнением Бернулли.

\frac{\partial (\rho A)}{\partial t} + \frac{\partial}{\partial l} \big( A \, \rho \, v \big) = 0

\frac{\partial (\rho A)}{\partial t} = - \frac{\partial}{\partial l} \big( A \, \rho \, v \big) = 0

dF_p \bigg |_{l_{\perp}} + dF_g \bigg |_{l_{\perp}} + dF_f \bigg |_{l_{\perp}}+ dF_N \bigg |_{l_{\perp}} =0

 dF_p \bigg |_l + dF_g \bigg |_l + dF_f \bigg |_l+ dF_N \bigg |_l = \frac{d I}{dt}\bigg |_l

- A dp + \rho \, A\, \delta l

 \, g \, \sin \theta - \frac{f \, \rho

 \, 

v^2}{

2 d} \

, A \, \

delta l = \

rho

 \, 

A \, v \, dv

Если дебит скважины на устье составляет \theta может быть выражен через абсолютные отметки глубин  

\sin \theta = \frac{dz}{dl}

\frac{dp}{dl} =A \, \rho \, g \, \frac{dz}{dl} -  \frac{v = \rho_s^2s \, q_s^2}{A}  \frac{d}{dl} \bigg( \frac{1}{A \, \rho} \bigg) -s

откуда можно выразить явно профиль скорости потока по стволу:

v(l) = \frac{\rho_s^2s \, q_s^2 s}{2 A^2 d} \frac{f}{\rho}

\rho(p) \, A(l)}

Подставляя 

\frac{dp}{dl} = \rho \, g \, \sin \frac{dz}{dl}theta - \frac{\rho_s^2 \, q_s^2}{A^2A}  \frac{d}{dl} \bigg( \frac{1}{A \, \rho} \bigg) - \frac{\rho_s^2 \, q_s^2 }{2 A^2 d} \frac{f}{\rho}

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

\frac{d}{dl}sin \bigg(theta = \frac{1dz}{\rho} dl}

и уравнение для давление примет вид:

\bigg) = -\frac{1dp}{\rho^2dl} \frac{d= \rho}{ dl} 
= -\, g \, \frac{1dz}{\rho^2}dl} -  \frac{d\rho_s^2 \rho, q_s^2}{dpA}  \frac{dpd}{ dl}
=- \bigg( \frac{c1}{A \, \rho} \bigg) - \frac{dp}{ dl}

где

\rho_s^2 \, q_s^2 }{2 A^2 d} \frac{f}{\rho}

Диаметр труб, вдоль которых идет движение воды, остается постоянным на долгом протяжении и меняется редко (например, километр НКТ и потов выход потока в колонну), и это позволяет решать задачу нахождения профиля давления на кусках постоянного диаметра 

\

\bigg( 1 -  \frac{c(p) \, \rho_s^2 \, q_s^2}{A^2}   \bigg )  \frac{dp}{dl} = \rho(p) \, g \, \frac{dz}{dl}  - \frac{\rho_s^2 \, q_s^2 }{2A^2} A^2 \frac{d}{dl} \bigg( \frac{f(p)1}{\rho(p)}} \bigg) - \frac{\rho_s^2 \, q_s^2 }{2 A^2 d} \frac{f}{\rho}

Процесс движения воды вдоль трубы происходит в состоянии термодинамического равновесия и плотность воды является функцией только давления Функция

Как будет показано ниже коэффициент трения

Если предположить постоянство коэффициента трения

, следовательно:

\frac{d}{dl} \bigg( \frac{1}{\rho} \bigg) = -\frac{1}{\rho^2} \frac{d \rho}{ dl} 
= - \frac{1}{\rho^2}\frac{d \rho}{dp} \frac{dp}{ dl}
=- \frac{c}{\rho} \frac{dp}{ dl}

где

\bigg( 1 -  \frac{c(p) \, \rho_s^2

p(l) = p_s + \rho \, g \, z(l) - \frac{\rho_s \, q_s^2 }{2A^2} A^2 d} \, f_s \, l = p_s + \rhobigg )  \frac{dp}{dl} = \rho(p) \, g \, z(l) - \frac{dpdz}{dhdl}  - \frac{\Bigg|_{loss}rho_s^2 \, l

Первые два слагаемых описывают гидростатический столб неподвижного флюида, а последнее слагаемое выражает потери на трение при движении флюида:

\frac{dp}{dl} \Bigg|_{loss} =  \frac{\rho_s \, q_s^2 }{2 A^2 d} \, f_s

В калькуляторе Well Flow Performance Calculator можно оценить величину потерь на трения для различных сценариев диаметров труб и дебитов скважин.

Коэффициент трения

Коэффициент трения Дарси

Для гладкой трубы

f = 64 \, \rm Re^{-1}

переходной режим течение

f = 0.32 \, \rm Re^{-0.25}

турбулентный режим течения

f = 0.184 \, \rm Re^{-0.2}

сильно турбулентный поток режим течения

 q_s^2 }{2 A^2 d} \frac{f(p)}{\rho(p)}

Функция

Как будет показано ниже коэффициент трения

Если предположить постоянство коэффициента трения

p(l) = p_s + \rho \, g \, z(l) - \frac{\rho_s \, q_s^2 }{2 A^2 d} \, f_s \, l

Pressure gradient will be:

\frac{dp}{dl} = \rho \, g \, \cos \theta(l) - \frac{\rho_s \, q_s^2 }{2 A^2 d} \, f_s 

where

The first term defines the hydrostatic column of static fluid while the last term defines the friction losses under fluid movement:

\frac{dp}{dl} \Bigg|_{loss} =  \frac{\rho_s \, q_s^2 }{2 A^2 d} \, f_s 

В калькуляторе Well Flow Performance Calculator можно оценить величину потерь на трения для различных сценариев диаметров труб и дебитов скважин.

Коэффициент трения

Коэффициент трения Дарси

Для гладкой

где

число Рейнольдса

Для переходных и турбулентных режимов течения коэффициент трения удовлетворяет эмпирической модели Колбрука-Уайта (Colebrook–White), которая учитывает шероховатость внутренней поверхности трубы

f

 = 

64 \, \

rm Re^{-1}

Существует множество явных аппроксимаций решения уравнения

f = 0.25 \, \bigg[ \log \bigg( \frac{\epsilon / d}{3.7065} - \frac{5.0272}{\rm Re} \log \Lambda \bigg)   \bigg]^{-2}

где

\Lambda = \frac{(\epsilon/d)}{3.827} - \frac{4.657}{\rm Re} \log \Bigg[  \bigg( \frac{\epsilon/d}{7.7918} \bigg)^{0.9924} + \bigg( \frac{5.3326}{208.815+Re} \bigg)^{0.9345} \Bigg]

Однако, в пределах измерительной погрешности (< 2 %) можно пользоваться универсальной корреляцией (Churchil) для всех режимов течения, от ламинарного до сильно турбулентного:

f = \frac{64}{\rm Re} \, \Bigg [ 1+ \frac{\big(\rm Re / 8 \big)^{12} }{ \big( \Theta_1 + \Theta_2 \big)^{1.5} }  \Bigg]^{1/12}

где

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

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

Тем не менее, зависимость от дебита является слабой. Из формулы

Зависимость коэффициента трения от давления формируется только через число Рейнольдса:

{\rm Re} = \frac{ d \, \rho_s \, q_s}{A \, \mu(p)}

отсюда следует, что зависимость коэффициента трения от давления формируется вязкостью

δμ/μ = 25 % при вариации μ = 2.4·10^-5 Па · с для p = 1 атм до μ = 3.0·10^-5 Па · с для 300 атм (cм. Свойства воды).

Для оценки числа Рейнольдса для нагнетаемой по 2.5 " НКТ воды можно пользоваться формулой

Отсюда видно, что при дебитах более 18 м³/сут число Рейнольдса становится больше 4,000 и режим течения является турбулентным и коэффициент трения можно считать практически постоянным вдоль ствола нагнетательной скважины.

А учитывая, что рост давления с глубиной сопровождается увеличением температуры, что компенсирует рост вязкости воды, то для большинства практических реализаций ППД можно полагать, что вариация коэффициента трения вдоль ствола не превышает 2-3 % и в оценках потери напора на трение принимать коэффициент трения постоянным

Профиль температуры 

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

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

В этом случае говорят о квазистационарном распределении температурного поля.

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

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

Температурный профиль

\rho \, c \, \frac{\partial T}{\partial t} = \frac{d}{dl} \, \bigg( \lambda \, \frac{dT}{dl} \bigg)  - \rho \, c \, v \, \frac{dT}{dl}

с начальным условием:

T(t=0, l) = T_g(l)

и граничным условием на поверхности:

T(t, l=0) = T_s(t)

Распределение температуры в массиве горных пород

\rho_e \, c_e \, \frac{\partial T_e}{\partial t} = \nabla ( \lambda_e \nabla T_e)

с начальным условием:

T_e(t=0, l, r) = T_g(l)

и граничным условием на бесконечном удалении от скважины:

T_e(t, l, r \rightarrow \infty) = T_g(l)

Геотермическое распределение температуры (также называемое геотермой) вдоль ствола скважины

T_g(l) = T_{0e} + \int_{z_0}^{z(l)} G_T(z) dz = T_{0e} + \int_{l_0}^l G_T(z(l)) \sin \theta dl 

геотермический градиент задается отношением регионального теплового потока из недр Земли

G_T(z(l)) = \frac{j_e}{\lambda_e(l)}

где

В регионах, где геотермический градиент остается постоянным

T_g(l) = T_{0e} + G_T \, z 

Однако в большом количестве практических случаев это не так и применение среднего по всему разрезу значения геотермического градиента для оценки геотермического распределения температур по формуле

f = 0.32 \, \rm Re^{-0.25}

f = 0.184 \, \rm Re^{-0.2}

где

Для переходных и турбулентных режимов течения коэффициент трения удовлетворяет эмпирической модели Колбрука-Уайта (Colebrook–White), которая учитывает шероховатость внутренней поверхности трубы

\frac{1}{\sqrt{f}} = -2 \, \log \Bigg( \frac{\epsilon}{3.7 \, d}  + \frac{2.51}{{\rm Re} \sqrt{f}} \Bigg)

Типичное значение шероховатости труб

Существует множество явных аппроксимаций решения уравнения

f = 0.25 \, \bigg[ \log \bigg( \frac{\epsilon / d}{3.7065} - \frac{5.0272}{\rm Re} \log \Lambda \bigg)   \bigg]^{-2}

где

\Lambda = \frac{(\epsilon/d)}{3.827} - \frac{4.657}{\rm Re} \log \Bigg[  \bigg( \frac{\epsilon/d}{7.7918} \bigg)^{0.9924} + \bigg( \frac{5.3326}{208.815+Re} \bigg)^{0.9345} \Bigg]

Однако, в пределах измерительной погрешности (< 2 %) можно пользоваться универсальной корреляцией (Churchil) для всех режимов течения, от ламинарного до сильно турбулентного:

f = \frac{64}{\rm Re} \, \Bigg [ 1+ \frac{\big(\rm Re / 8 \big)^{12} }{ \big( \Theta_1 + \Theta_2 \big)^{1.5} }  \Bigg]^{1/12}

где

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

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

Тем не менее, зависимость от дебита является слабой. Из формулы

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

Зависимость коэффициента трения от давления формируется только через число Рейнольдса:

При этом число Рейнольдса

{\rm Re} = \frac{ d \, \rho_s \, q_s}{A \, \mu(p)}

отсюда следует, что зависимость коэффициента трения от давления формируется вязкостью

δμ/μ = 25 % при вариации μ = 2.4·10^-5 Па · с для p = 1 атм до μ = 3.0·10^-5 Па · с для 300 атм (cм. Свойства воды).

Это приводит к 25 % вариации коэффициента трения для ламинарного потока (в котором сила трения минимальна) и порядка 4.5 % для турбулентного потока (и максимальным вкладом трения).

Для оценки числа Рейнольдса для нагнетаемой по 2.5 " НКТ воды можно пользоваться формулой

Отсюда видно, что при дебитах более 18 м³/сут число Рейнольдса становится больше 4,000 и режим течения является турбулентным и коэффициент трения можно считать практически постоянным вдоль ствола нагнетательной скважины.

А учитывая, что рост давления с глубиной сопровождается увеличением температуры, что компенсирует рост вязкости воды, то для большинства практических реализаций ППД можно полагать, что вариация коэффициента трения вдоль ствола не превышает 2-3 % и в оценках потери напора на трение принимать коэффициент трения постоянным

Профиль температуры 

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

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

В этом случае говорят о квазистационарном распределении температурного поля.

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

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