Analytical solution of the boundary value problem of mathematical modeling for non-stationary oil flows through the trunk pipeline in the presence of internal pressure sources

Full Text

Введение

Магистральные трубопроводы (МТП), предназначенные для транспортировки нефти и нефтепродуктов, имеют большую протяженность, ввиду чего с позиции задачи управления их следует рассматривать как объекты с распределенными параметрами (ОРП) [1–4]. Для решения широкого круга прикладных задач оптимального управления системами с распределенными параметрами, таких как поиск оптимальных программ управления и синтез систем регулирования с обратной связью, востребованными остаются аналитические модели управляемых процессов [3–7]. Такие проблемно-ориентированные модели могут быть получены путем решения соответствующих краевых задач математической физики [8–11]. Проблемная ориентированность моделей для целей синтеза систем оптимального управления подразумевает применение специальных методов решения краевой задачи, в частности метода функций Грина [2, 9–11].

Движение условно несжимаемой нефти по трубопроводу постоянного диаметра описывается уравнениями математической физики гиперболического типа [8, 9]. В частности, в работах [1, 5] на основе системы одномерных уравнений нестационарного движения жидкости И.А. Чарного [12], широко используемой для описания гидродинамики магистральных нефтепроводов [13, 14], предложена краевая задача математического моделирования процесса трубопроводного транспорта нефти при наличии внутренних сосредоточенных источников давления и расхода. Численное решение задачи получено в работе [1], в работах [4, 5] с использованием метода функций Грина получено решение линеаризованной задачи для пространственно-временного распределения скорости потока при наличии в трубопроводе внутренних сосредоточенных источников давления. В то же время более востребованными для решения задач оптимального управления режимами работы магистральных нефтепроводов являются модели, описывающие пространственно-временное распределение внутренних давлений в трубопроводе [5, 12–14]. В статье будет представлено аналитическое решение описанной выше краевой задачи математического моделирования процесса трубопроводного транспорта нефти для пространственно-временного распределения давлений в МТП при наличии внутренних источников давления, сосредоточенных в некоторых точках расположения нефтеперекачивающих станций (НПС).

 

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

Взаимосвязь внутреннего избыточного давления P и средней по сечению скорости ω потока нефти плотностью ρ, движущейся по трубопроводу постоянного диаметра D и длиной L, в любой точке x, $x\in \left[\mathrm{0,}L\right]$по направлению движения потока и в любой момент времени t, $t\ge 0$ описывается согласно [1, 10–14] системой двух одномерных дифференциальных уравнений в частных производных:

$\frac{\partial P\left(x,t\right)}{\partial x}=-\rho \cdot \left(\frac{\partial \omega \left(x,t\right)}{\partial t}+2\cdot \overline{\upsilon }\cdot \omega \left(x,t\right)+{\overline{\upsilon }}_0+g\cdot \sin \alpha \left(x\right)\right)+u\left(x,t\right)$,     (1)

$\frac{\partial P\left(x,t\right)}{\partial t}=-c^2\cdot \rho \cdot \frac{\partial \omega \left(x,t\right)}{\partial x}$,      (2)

где $\alpha \left(x\right)$ – угол наклона оси трубопровода к произвольной горизонтальной поверхности; $g$ – ускорение свободного падения; c – скорость распространения волн в жидкости, текущей в стальной трубе с толщиной стенки $d$, определяется по формуле Жуковского [1, 12, 14].

Система уравнений (1), (2) дополняется начальными условиями

$\omega \left(x\mathrm{,0}\right)={\omega }_0$   (3)

или соответствующими условиям (3) начальными условиями

$P\left(x\mathrm{,0}\right)=P_0-\rho \cdot \left({\overline{\upsilon }}_0\cdot x+2\cdot \overline{\upsilon }\cdot {\omega }_0\cdot x+g\cdot z\left(x\right)\right)$,     (4)

описывающими исходное стационарное состояние в трубопроводе, которое сохраняется до момента ${\tau }_{s_k}\ge 0$ появления в некоторой внутренней точке с координатой $x_k$ источника давления, величина которого во времени меняется согласно зависимости $u_{s_k}\left(t\right)$. В (1) и (4) ${\overline{\upsilon }}_0$ и $\overline{\upsilon }$ – коэффициенты линеаризации. Тогда в уравнении (1) функция распределения $u\left(x,t\right)$ по длине трубопровода внутренних источников давления, приложенных в точках $x_k\in \left[\mathrm{0,}L\right]$ расположения НПС, имеет вид

$u\left(x,t\right)=\sum _{k=1}^K\sum _{s_k=1}^{S_k}{u}_{s_k}\left(t\right)\cdot \delta \left(x-x_k\right)$,     (5)

где $\delta \left(x-x_k\right)$ – функции Дирака; $k=\overline{\mathrm{1,}K}$ – индекс работающих НПС; $s_k=\overline{\mathrm{1,}{S}_k}$ – индекс работающих на каждой k-ой НПС насосных агрегатов.

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

$P\left(\mathrm{0,}t\right)=P_0$, $P\left(L,t\right)=P_L$ (6)

или соответствующие условиям (6) граничные условия

$\frac{\partial \omega \left(\mathrm{0,}t\right)}{\partial x}=0$, $\frac{\partial \omega \left(L,t\right)}{\partial x}=0$.    (7)

Совместно уравнения (1), (2) с начальными условиями (3), (4) и граничными условиями (6), (7) составляют краевую задачу математического моделирования нестационарного течения нефти по магистральному трубопроводу при наличии внутренних источников давления с пространственно-временной функцией распределения в форме (5).

Схемы линеаризации уравнения (1) описаны в [10, 12] и основаны на замене произведения нелинейной функции гидравлических потерь $\lambda \left(\omega \right)$ и квадрата средней скорости потока $\omega$ в трубопроводе линейной зависимостью

$\frac{\lambda \left(\omega \left(x,t\right)\right)\cdot {\omega }^2\left(x,t\right)}{2\cdot D}=2\cdot \overline{\upsilon }\cdot \omega \left(x,t\right)+{\overline{\upsilon }}_0$,   (8)

как показано на рис. 1.

 

Рис. 1. Схема выбора коэффициентов линеаризации уравнения (1)

 

Схема выбора коэффициентов линеаризации (см. рис. 1) содержит две точки пересечения прямой $2\cdot \overline{\upsilon }\cdot \omega +{\overline{\upsilon }}_0$ и исходной нелинейной зависимости $\frac{\lambda \left(\omega \right)\cdot {\omega }^2}{2\cdot D}$, которая может быть построена по формулам гидравлических сопротивлений для зоны гидравлически гладких труб и зоны смешанного трения турбулентного режима течения [15], характерных для большинства режимов транспортировки нефти и нефтепродуктов.

Система уравнений (1), (2) приводится в работах [10, 11] к каноническому виду в форме 

$\frac{{\partial }^2\omega \left(x,t\right)}{\partial t^2}-{с}^2\cdot \frac{{\partial }^2\omega \left(x,t\right)}{\partial x^2}+2\cdot \overline{\upsilon }\cdot \frac{\partial \omega \left(x,t\right)}{\partial t}=\frac{U\left(x,t\right)}{\rho }$,   (9)

где $U\left(x,t\right)$ определяется как

$U\left(x,t\right)=\sum _{k=1}^K\sum _{s_k=1}^{S_k}\frac{d{u}_{s_k}\left(t\right)}{dt}\cdot \delta \left(x-x_k\right)=\sum _{k=1}^K\sum _{s_k=1}^{S_k}{u^{\text{'}}}_{s_k}\left(t\right)\cdot \delta \left(x-x_k\right)$.   (10)

В качестве типовой программы пуска (останова) насосного агрегата может рассматриваться программа с постоянной скоростью роста ($+u_{s_k}^{\max }$) или снижения ($-u_{s_k}^{\max }$) перепада давления на насосе, приведенная в [5], функциональная зависимость от времени которой имеет вид

$u_{s_k}^{}\left(t\right)=\left\{\begin{array}{l}0,t<t_{s_k}^{},\\ \pm u_{s_k}^{\max }\cdot t,t_{s_k}^{}\le t\le \left(t_{s_k}^{}+\Delta \right),\\ \pm u_{s_k}^{\max }\cdot \Delta ,t>\left(t_{s_k}^{}+\Delta \right),\end{array}\right.$ (11)

${u^{\text{'}}}_{s_k}^{}\left(t\right)=\left\{\begin{array}{l}0,t<t_{s_k}^{},\\ \pm u_{s_k}^{\max },t_{s_k}^{}\le t\le \left(t_{s_k}^{}+\Delta \right),\\ 0,t>\left(t_{s_k}^{}+\Delta \right)\end{array}\right.$(12)

и приведена на рис. 2.

 

Рис. 2. Программа пуска и останова s_k-го насосного агрегата k-й НПС

 

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

Аналитическое решение задачи относительно скорости потока жидкости в трубопроводе получено в [11] с помощью метода функций Грина. Представим краевую задачу (9) относительно пространственно-временной функции распределения скоростей потока $\omega (x,t)$ с начальными условиями (3), граничными условиями (7) и функцией распределения внутренних источников давления (10) в стандартной форме. Под стандартной формой будем понимать эквивалентную задаче (9), (3), (7), (10) краевую задачу с нулевыми начальными и однородными граничными условиями [2–4, 9]. Характеристикой краевой задачи в стандартной форме или ее импульсной переходной функцией будет являться функция Грина

$G_{\omega }\left(x,\xi ,t-\tau \right)$,   (13)

такая, что позволяет получить аналитическую зависимость для $\omega (x,t)$ в виде интеграла

$\omega \left(x,t\right)=\underset{0}{\overset{t}{\int }}\underset{0}{\overset{L}{\int }}{G}_{\omega }\left(x,\xi ,t-\tau \right)\cdot w_{\omega }\left(\xi ,\tau \right)d\xi d\tau$,   (14)

где $w_{\omega }(\xi ,\tau )$ – стандартизирующая функция, позволяющая представить исходную краевую задачу (9), (3), (7), (10) с ненулевыми начальными и, в общем случае, неоднородными граничными условиями в виде эквивалентной краевой задачи с нулевыми начальными и однородными граничными условиями.

С учетом аналитического решения рассматриваемой краевой задачи относительно скорости потока жидкости в трубопроводе $\omega (x,t)$, полученной в [11], и правил, установленных в [2–4], получим следующие выражения для $G_{\omega }\left(x,\xi ,t-\tau \right)$:

$\begin{array}{c}{G}_{\omega }\left(x,\xi ,t-\tau \right)=\frac{1}{L}\cdot \left(\frac{e^{2\alpha \cdot \left(t-\tau \right)}-1}{2\cdot \alpha }+\right.\\ +2\cdot \sum _{n=1}^{N^{\ast }}\cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \cos \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot \left(\frac{e^{\alpha \cdot \left(t-\tau \right)}\cdot sh\left({\beta }_n^{}\cdot \left(t-\tau \right)\right)}{{\beta }_n^{}}\right)+\\ \left.+2\cdot \sum _{n=N^{\ast }+1}^{\infty }\cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \cos \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot \left(\frac{e^{\alpha \cdot \left(t-\tau \right)}\cdot \sin \left({\beta }_n^{\ast }\cdot \left(t-\tau \right)\right)}{{\beta }_n^{\ast }}\right)\right),\end{array}$(15)

$\alpha =-\overline{\upsilon },{\beta }_n=\sqrt{{\overline{\upsilon }}^2-{\left(\frac{\pi \cdot n\cdot c}{L}\right)}^2},{\beta }_n^{\ast }=\sqrt{{\left(\frac{\pi \cdot n\cdot c}{L}\right)}^2-{\overline{\upsilon }}^2},N^{\ast }=⌊\frac{L\cdot \overline{\upsilon }}{\pi \cdot c}⌋$ (16)

где $N^{*}$– целое число, полученное округлением в меньшую сторону отношения $\frac{L·\overline{\nu }}{\pi ·c}$. Выражение для стандартизирующей функции $w_{\omega }(\xi ,\tau )$ принимает вид

$w_{\omega }\left(\xi ,\tau \right)=\left(2\cdot \overline{\upsilon }\cdot \delta \left(\tau \right)+{\delta }^{\text{'}}\left(\tau \right)\right)\cdot {\omega }_0+\frac{U\left(\xi ,\tau \right)}{\rho }$.   (17)

Проинтегрировав в (14) по пространственной координате ξ, с учетом свойств δ-функции и линейности интеграла свертки получим выражение для $\omega (x,t)$ в более простом виде:

$\omega \left(x,t\right)={\omega }_0+\frac{1}{\rho }\cdot \sum _{k=1}^K\sum _{s_k=1}^{S_k}\underset{0}{\overset{t}{\int }}{G}_{\omega }\left(x,x_k,t-\tau \right)\cdot {u^{\text{'}}}_{s_k}\left(\tau \right)d\tau$.   (18)

Аналитическое решение исходной краевой задачи относительно пространственно-временного распределения давлений $P\left(x,t\right)$ можно получить в виде функции Грина $G_p\left(x,\xi ,t-\tau \right)$ путем подстановки (15) в уравнение (2) и последующего интегрирования результата подстановки

$-\frac{1}{c^2\cdot \rho }\cdot \frac{\partial G_p\left(x,\xi ,t-\tau \right)}{\partial t}=\frac{\partial G_{\omega }\left(x,\xi ,t-\tau \right)}{\partial x}=$

$=-\frac{2}{L}\cdot \sum _{n=1}^{N^{\ast }}\frac{\pi \cdot n}{L}\cdot \cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \sin \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot \underset{f_{1n}\left(t\right)}{\underbrace{\left(\frac{e^{\alpha \cdot \left(t-\tau \right)}\cdot sh\left({\beta }_n^{}\cdot \left(t-\tau \right)\right)}{{\beta }_n^{}}\right)}}-$ (19)

$-\frac{2}{L}\cdot \sum _{n=N^{\ast }+1}^{\infty }\frac{\pi \cdot n}{L}\cdot \cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \sin \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot \underset{f_{2n}\left(t\right)}{\underbrace{\left(\frac{e^{\alpha \cdot \left(t-\tau \right)}\cdot \sin \left({\beta }_n^{\ast }\cdot \left(t-\tau \right)\right)}{{\beta }_n^{\ast }}\right)}}$

по схеме, приведенной далее, и составления новой стандартизирующей функции ${\omega }_p(\xi ,\tau )$, удовлетворяющей начальным условиям (4) и граничным условиям (6).

Интегрирование по t правой части (19) сводится к получению двух интегралов:

$I_{1n}\left(t\right)=\int f_{1n}\left(t\right)dt$,   (20)

$I_{2n}\left(t\right)=\int f_{2n}\left(t\right)dt$.   (21)

Введем новые переменные интегрирования:

${\phi }_1={\beta }_n^{}\cdot \left(t-\tau \right)$,   (22)

${\phi }_2={\beta }_n^{*}\cdot \left(t-\tau \right)$.   (23)

Тогда (20), (21) примут вид

$I_{1n}\left({\phi }_1\right)=\frac{1}{{\beta }_n}\int f_{1n}\left({\phi }_1\right)d{\phi }_1$,   (24)

$I_{2n}\left({\phi }_2\right)=\frac{1}{{\beta }_n^{*}}\int f_{2n}\left({\phi }_2\right)d{\phi }_2$,   (25)

Получение интеграла (24) не представляет особой сложности после раскрытия выражения для $\text{sh}\left({\phi }_1\right)$. Результат интегрирования (24) имеет вид

$I_{1n}\left({\phi }_1\right)=\frac{e^{k_1\cdot {\phi }_1}}{2\cdot {\beta }_n^2}\cdot \left(\frac{e^{{\phi }_1}}{k_1+1}-\frac{e^{-{\phi }_1}}{k_1-1}\right)+C_{1n}$,   (26)

$k_1=\frac{\alpha }{{\beta }_n}$.   (27)

Для интегрирования (25) выделим интеграл:

$J_{1n}\left({\phi }_2\right)=\int e^{k_2\cdot {\phi }_2}\cdot \sin \left({\phi }_2\right)d{\phi }_2$,   (28)

$k_2=\frac{\alpha }{{\beta }_n^{*}}$.   (29)

Для получения решений в (28) выполним двойное интегрирование по частям и сведем интеграл к «самому себе». Для интегрирования по частям в (28) положим:

$\vartheta =e^{k_2\cdot {\phi }_2}\Rightarrow d\vartheta =k_2\cdot e^{k_2\cdot {\phi }_2}d{\phi }_2$,   (30)

$d\sigma =\sin \left({\phi }_2\right)d{\phi }_2\Rightarrow \sigma =-\cos \left({\phi }_2\right)$,  (31)

тогда после первого этапа интегрирования по частям получим:

$\begin{array}{c}{J}_{1n}\left({\phi }_2\right)=\int \vartheta d\sigma =\vartheta \cdot \sigma -\int \sigma d\vartheta =\\ =-e^{k_2\cdot {\phi }_2}\cdot \cos \left({\phi }_2\right)+k_2\cdot \int e^{k_2\cdot {\phi }_2}\cdot \cos \left({\phi }_2\right)d{\phi }_2.\end{array}$ (32)

Для повторного интегрирования по частям (32) положим вместо (31)

$d\sigma =\cos \left({\phi }_2\right)d{\phi }_2\Rightarrow \sigma =\sin \left({\phi }_2\right)$,    (33)

тогда после второго этапа интегрирования по частям (32) примет вид

$\begin{array}{c}{J}_{1n}\left({\phi }_2\right)=\int \vartheta d\sigma =\\ =-e^{k_2\cdot {\phi }_2}\cdot \cos \left({\phi }_2\right)+k_2\cdot e^{k_2\cdot {\phi }_2}\cdot \sin \left({\phi }_2\right)+k_2^2\cdot \underset{J_{1n}\left({\phi }_2\right)}{\underbrace{\int e^{k_2\cdot {\phi }_2}\cdot \sin \left({\phi }_2\right)d{\phi }_2}}.\end{array}$ (34)

Как видно из выражения (34), интеграл $J_{1n}\left({\phi }_2\right)$ свелся к «самому себе», и после незначительных преобразований получим

$J_{1n}\left({\phi }_2\right)=\frac{e^{k_2\cdot {\phi }_2}}{1+k_2^2}\cdot \left(k_2\cdot \sin \left({\phi }_2\right)-\cos \left({\phi }_2\right)\right)+{С}_{21n}$.   (35)

Выражение для (21) примет вид

$I_{2n}\left({\phi }_2\right)=\frac{1}{{{\beta }_n^{*}}^2}\cdot \frac{e^{k_2\cdot {\phi }_2}}{1+k_2^2}\cdot \left(k_2\cdot \sin \left({\phi }_2\right)-\cos \left({\phi }_2\right)\right)+C_{2n}$.   (36)

После подстановки (22) и (27) в (26) получим:

$I_{1n}\left(t\right)=\frac{e^{\alpha \cdot \left(t-\tau \right)}}{2\cdot {\beta }_n}\cdot \left(\frac{e^{{\beta }_n^{}\cdot \left(t-\tau \right)}}{\alpha +{\beta }_n}-\frac{e^{-{\beta }_n^{}\cdot \left(t-\tau \right)}}{\alpha -{\beta }_n}\right)+C_{1n}$.   (37)

После подстановки (23) и (28) в (36) получим:

$I_{2n}\left(t\right)=\frac{e^{\alpha \left(t-\tau \right)}}{{\alpha }^2+{{\beta }_n^{*}}^2}\cdot \left(\frac{\alpha }{{\beta }_n^{*}}\cdot \sin \left({\beta }_n^{*}\cdot \left(t-\tau \right)\right)-\cos \left({\beta }_n^{*}\cdot \left(t-\tau \right)\right)\right)+C_{2n}$.   (38)

Тогда результат интегрирования (19) окончательно примет вид:

$\begin{array}{c}{G}_p\left(x,\xi ,t-\tau \right)=\frac{2\cdot c^2\cdot \rho }{L^2}\cdot \left(\sum _{n=1}^{N^{\ast }}\pi \cdot n\cdot \cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \sin \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot I_{1n}\left(t\right)+\right.\\ +\left.\sum _{n=N^{\ast }+1}^{\infty }\pi \cdot n\cdot \cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \sin \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot I_{2n}\left(t\right)\right).\end{array}$ (39)

Функция Грина (39) является решением задачи (1), (2) при нулевых начальных условиях, обращается в ноль при t=0, начальные условия являются статичными и не зависят от τ (следовательно, можно принять τ=0), тогда постоянные интегрирования $C_{1n}$и $C_{2n}$ должны удовлетворять тождеству

$\begin{array}{c}{G}_p\left(x,\xi \mathrm{,0}\right)=\\ =\frac{2\cdot c^2\cdot \rho }{L^2}\cdot \left(\sum _{n=1}^{N^{\ast }}\pi \cdot n\cdot \cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \sin \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot \left(-\frac{1}{{\alpha }^2-{\beta }_n^2}+C_{1n}\right)+\right.\\ +\left.\sum _{n=N^{\ast }+1}^{\infty }\pi \cdot n\cdot \cos \left(\frac{\pi \cdot n\cdot \xi }{L}\right)\cdot \sin \left(\frac{\pi \cdot n\cdot x}{L}\right)\cdot \left(-\frac{1}{{\alpha }^2+{{\beta }_n^{*}}^2}+C_{2n}\right)\right)=0.\end{array}$ (40)

Отсюда $C_{1n}$и $C_{2n}$ определяются как

$C_{1n}^{}=\frac{1}{{\alpha }^2-{\beta }_n^2}$, $C_{2n}^{}=\frac{1}{{\alpha }^2+{{\beta }_n^{\ast }}^2}$.   (41)

Соответствие граничных условий (7) условиям (6) и начальных условий (3) условиям (4) позволяет провести замену ${\omega }_0$ в (17) следующим выражением:

${\omega }_0=\frac{P_0-P_L}{2\cdot \rho \cdot \overline{\upsilon }\cdot L}+g\cdot \frac{z\left(0\right)-z\left(L\right)}{2\cdot \overline{\upsilon }\cdot L}-\frac{{\overline{\upsilon }}_0}{2\cdot \overline{\upsilon }}$,   (42)

после чего вместо стандартизирующей функции $w_{\omega }\left(\xi ,\tau \right)$ для граничных условий (7) и начальных условий (3) можно использовать функцию ${\omega }_p(\xi ,\tau )$ для граничных условий (6) и начальных условий (4) в форме

$w_p\left(\xi ,\tau \right)=\left(2\cdot \overline{\upsilon }\cdot \delta \left(\tau \right)+{\delta }^{\text{'}}\left(\tau \right)\right)\cdot \left(\frac{P_0-P_L}{2\cdot \rho \cdot \overline{\upsilon }\cdot L}+g\cdot \frac{z\left(0\right)-z\left(L\right)}{2\cdot \overline{\upsilon }\cdot L}-\frac{{\overline{\upsilon }}_0}{2\cdot \overline{\upsilon }}\right)+\frac{U\left(\xi ,\tau \right)}{\rho }$,   (43)

а аналитическое решение краевой задачи (1), (2) относительно $P\left(x,t\right)$ находится в виде интеграла:

$P\left(x,t\right)=\underset{0}{\overset{t}{\int }}\underset{0}{\overset{L}{\int }}{G}_p\left(x,\xi ,t-\tau \right)\cdot w_p\left(\xi ,\tau \right)d\xi d\tau$.   (44)

Проинтегрировав в (44) по пространственной координате ξ аналогично (18), получим решение для $P\left(x,t\right)$ в виде

$P\left(x,t\right)=p_0\left(x\right)+\sum _{k=1}^K\sum _{s_k=1}^{S_k}\underset{0}{\overset{t}{\int }}{G}_p\left(x,x_k,t-\tau \right)\cdot {u^{\text{'}}}_{s_k}\left(\tau \right)d\tau$,   (45)

где $p_0\left(x\right)=P_0\cdot \left(1-\frac{x}{L}\right)+P_L\cdot \frac{x}{L}+\rho \cdot g\cdot \left(z\left(0\right)\cdot \left(1-\frac{x}{L}\right)-z\left(x\right)+z\left(L\right)\cdot \frac{x}{L}\right)$.   (46)

Сравнительный анализ результатов аналитического и численного решения краевой задачи

Полученные интегральные формы (18) и (45) позволяют получить в аналитическом виде в любой момент времени функции распределения давлений и средних скоростей потока нефти по длине магистрального трубопровода при наличии внутренних источников давления, которые в достаточно полной мере описывают гидродинамику изотермических нефтепроводов. Методика численного решения в пакете MathCAD нелинейной краевой задачи математического моделирования гидродинамики магистрального нефтепровода представлена в работе [1].

В данном разделе приводится сравнение результатов аналитического решения в форме (18), (45) и численного решения по методике [1] линеаризованной краевой задачи (1) – (7) на примере технологического участка магистрального трубопровода, имеющего протяженность L = 200 км, внутренний диаметр D = 1200 мм. Транспортируемый продукт – нефть плотностью $\rho =850кг/{м}^3$, вязкостью $\nu =25сСт$, для условий течения которой определены коэффициенты линеаризации ${\overline{\upsilon }}_0=-7.14\times {10}^{-3}м/{с}^2$, $\overline{\upsilon }=0.016{с}^{-1}$ согласно схеме (см. рис. 1).

В начальный момент времени режим течения нефти в трубопроводе стационарный со средней по сечению скоростью потока ${\omega }_0=0.94м/с$, давления в начале и в конце трубопровода составляют соответственно $P_0=4\times {10}^6Па$, $P_L=1\times {10}^5Па$, профиль трассы ровный. В момент времени ${\tau }_{pk}^{}=0$ в точке расположения НПС $x_{pk}^{}=100$ км происходит пуск насосного агрегата по программе (11), (12) с параметрами $u_{pk}^{\max }=2\times {10}^5Па/с$ и $\Delta =10с$. Для расчетов учитывались первые 400 членов ряда в (15) и (40), при N^*=1 численный расчет по методике [1] реализован на сетке размером 400×400.

Результаты расчетов распределения по длине трубопровода средних по сечению скоростей потока представлены на рис. 3, распределения давлений – на рис. 4 для четырех моментов времени (10 с, 35 с, 70 с, 200 с).

 

Рис. 3. Распределение скорости потока по длине трубопровода в различные моменты времени: а – численное решение; б – аналитическое решение

 

Рис. 4. Распределение давления по длине трубопровода в различные моменты времени: а – численное решение; б – аналитическое решение

 

Сравнительный анализ результатов численного и аналитического решения линеаризованной краевой задачи (1) – (7) показывает:

– высокую степень совпадения результатов расчетов во временной области: характер распределения функций $P\left(x,t\right)$ и $\omega \left(x,t\right)$ в частности положение фронта волны давления (скорости потока), полученные как по численной, так и по аналитической модели, совпадают в моменты времени, отличные друг от друга менее чем на 1 %;

– высокую степень совпадения результатов расчетов по пространственной координате: абсолютные значения функций $P\left(x,t\right)$ и $\omega \left(x,t\right)$, полученные как по численной, так и по аналитической модели и совпадающие по характеру в некоторые моменты времени, отличаются друг от друга менее чем на 1 %;

– зависимость точности совпадения результатов расчетов от количества узлов сетки численной модели (шага расчета по пространственной и временной координате) и количества членов ряда в (15) и (40). Указанная выше точность совпадения результатов расчетов достигается для следующих параметров численной и аналитической моделей: шаг сетки по пространственной координате численной модели $\Delta x\le 0.5км$, шаг сетки по временной координате численной модели $\Delta t\le 0.5c$, количество членов ряда в (15) и (40) $n\ge 400$.

Заключение

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

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

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

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

About the authors

Alexander A. Afinogentov

Samara State Technical University

Email: Pondex@yandex.ru
ORCID iD: 0000-0003-3498-6587
SPIN-code: 9243-2840
Scopus Author ID: 56695164500
ResearcherId: N-4453-2014

Ph.D. (Techn.), Associate Professor

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

Yulia A. Tychinina

Samara State Technical University

Author for correspondence.
Email: ytychinina@list.ru
SPIN-code: 5315-1230

Ph.D. (Techn.), Associate Professor

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

References

Supplementary files