Моделирование характеристики гидравлической сети

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

$$\Delta p = RQ^2,$$

где $R$ - коэффициент, задающий крутизну параболы. Чем больше значение $R$, тем круче парабола. Назовём $R$ гидравлическим сопротивлением (по аналогии с электрическим сопротивлением). Для определения значения $R$ необходимо найти расчётным или опытным путём потери давления в трубопроводе при каком либо расходе.

Под потерями давления подразумевается снижение значения полного давления, вызванное потерями энергии при движении среды в связи с наличием трения и местными потерями. Нивелирная составляющая (давление, вызванное разницей в высотных положениях входного и выходного сечения трубопровода) не должна входить в $\Delta p$. Размерности величин расхода и давлений не важны, главное - не менять выбранные размерности во время дальнейшего расчёта. В нашем случае потери давления будут выражаться в метрах столба жидкости $\Delta H \equiv \Delta p$, а расход в $м^3/ч$.

Один раз определив для трубопровода значение $R$, мы сможем находить значение потерь давления для любого расхода, а так же по потерям давления определять расход. Следует заметить, что значение коэффициента трения при при малых скоростях движения среды зависит от скорости движения среды. В связи с этим, при низких скоростях движения среды, значение $R$ при разных расходах будет отличаться.

In [1]:
from scipy.interpolate import interp1d  #класс для интерполяции для моделирования 
                                        #напорной характеристики насоса
import numpy as np #работа с массивами
from math import sqrt
In [2]:
class Pipeline(object):
    """Объект трубопровод"""
    def __init__(self, Q, dH):
        """
        Q, dH: значение расхода и потерь давления при заданном расходе (единицы измерения не важны)
        """
        self.Q = Q
        self.dH = dH
        self.R = dH / Q / Q # гидравлическое сопротивление

    def dH_Q(self, Q, H_start = 0):
        """
        Q: расход через трубопровод
        H_start: - начальная ордината, откуда выходит характеристика сети
        return: значение потерь давления при заданном расходе + H_start
        """
        return self.R * Q * Q + H_start

    def Q_dH(self, dH, H_start = 0 ):
        """
        dH: потери давления
        return: расход, при котором потери давления составляют dH - H_start
        """
        return sqrt((dH - H_start) / self.R)

    def get_curve (self, Qmin, Qmax, H_start=0., n=51):
        """
        Возвращает характеристику трубопровода
        Qmin: начальная абсцисса
        Qmax: конечная абсцисса
        H_start: значение H при Q = 0
        n: число точек в графике
        return: кортеж (Q, H), где Q - список значений расхода,
                                   H - список значений потерь давления для заданных Q
        """
        H = []
        Q = []
        dq = (Qmax - Qmin) / (n-1)
        q = Qmin
        for _ in range(n):
            Q.append(q)
            H.append(self.dH_Q(q) + H_start)
            q += dq
        return Q, H

В приведённом выше коде класса Pipeline требует пояснения переменная H_start - начальная ордината характеристики трубопровода. Данная переменная нужна для приведения системы координат характеристики трубопровода к системе коодинат напорной характеристики насоса.

Предположим, что насос качает жидкость из бака, уровень жидкости в котором относительно оси всаса насоса составляет 4 м. Уровень жидкости в баке поддерживается постоянным вне зависимости от величины расхода из него. Высотное положение конца трубопровода такое же как и у его начала (всаса насоса), жидкость из трубопровода выливается в приёмный бак, уровень жидкости в котором расположен ниже оси трубопровода.

В данном случае часть работы по преодолению сопротивления трубопровода движению жидкости совержается за счёт давления столба жидкости высотой 4 метра. Данный столб жидкости можно рассматривать как насос с напором 4 м, установленный на всасе исследуемого насоса. Расход через трубопровод будет больше, чем в случае отсутствия "помощника". Чтобы отразить в системе координат напорной характеристики исследуемого насоса наличие "помощника" и увеличения в связи с этим расхода, начало характеристики трубопровода необходимо сместить на 4 метра вниз, т.е. значение нивелирной составляющей потери напора составляетH_start = -4.

Теперь рассмотрим случай, когда уровень жидкости в баке-источнике ниже уровня жидкости в баке-приёмнике на те же 4 м. В этом случае исследуемый насос, кроме сопротивления трубопровода, должен совершить работу по подъёму жидкости на высоту 4 метра. Расход жидкости в трубопроводе будет меньше, чем при отсутствии необходимости подъёма жидкости. Чтобы отразить в системе координат напорной характеристики исследуемого насоса уменьшение расхода в связи с необходимостью совершения насосом дополнительной работы по подъёму жидкости, точка начала характеристики трубопровода сдвигается на 4 метра вверх, т.е. значение нивелирной составляющей потери напора составляет H_start = 4.

Пусть имеется трубопровод потери давления в котором при расходе $Q = 3\ м^3/ч$ составляют $\Delta H = 14\ м$

In [3]:
Q_nom = 3; dH_nom = 14
pl = Pipeline(Q_nom, dH_nom) #создаём объект Трубопровод

Так же нам нужен насос.

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

In [4]:
#Возьмём несколько точек с напорной характеристики (включая две крайних)
Qs = np.array([1.19, 2, 3, 4, 4.43])
Hs = np.array([18.01, 16.55, 13.64, 8.91, 6.14])
#Напорная характеристика насоса
f_pump = interp1d(Qs, Hs, kind = 'cubic')
Qs_int = np.linspace(Qs[0], Qs[-1], 51) #значения подач для построения напорной характеристики насоса
In [5]:
Qs_pl, Hs_pl = pl.get_curve(0, 4) #создаём набор точек для построения характеристики трубопровода
Hs_pl = np.array(Hs_pl)
height = 4 #нивелирная составляющая, м
In [6]:
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq
from scipy.optimize import brentq #решение нелинейных уравнений
In [7]:
#Функция для поиска решения. Решение - значение Q при котором функция возвращает ноль.
#Данная функция передаётся в качестве первого параметра в brentq,
#который ищет Q при котором функция возвращает ноль,
#второй и третий параметр brentq - границы диапазона, внутри которого ищется решение. 
def f(Q):
    return pl.dH_Q(Q) + H_start - f_pump(Q)
In [8]:
#Варианты значений нивелирных составляющих
Hs_start = [-height, 0, height]
In [9]:
#С помощью brentq найдём значения Q точек пересечения характеристик трубопроводов с напорной
#характеристикой насоса (поиск корня нелинейного уравнения pl.dH_Q(Q) + H_start - f_pump(Q) = 0)
Qs_pl_pump = []
for h in Hs_start:
    H_start = h 
    Qs_pl_pump.append(brentq(f,Qs[0], Qs[-1]))
In [10]:
#Найдём значения H точек пересечения характеристик трубопроводов с напорной характеристикой насоса
Hs_pl_pump = []
for i in range(len(Hs_start)):
    Hs_pl_pump.append(pl.dH_Q(Qs_pl_pump[i], Hs_start[i]))
In [11]:
import matplotlib.pyplot as plt #библотека для построения графиков
plt.rcParams['figure.figsize'] = [14, 7]
plt.rcParams.update({'font.size': 14})
plt.plot(Qs_int, interp1d(Qs, Hs, kind='cubic')(Qs_int), label = "Напорная характеристика насоса")
plt.plot(Qs_pl, Hs_pl, label = "Нивелирная составляющая отсутствует")
plt.plot(Qs_pl, Hs_pl + height, label = "Подъём жидкости на 4 м")
plt.plot(Qs_pl, Hs_pl - height, label = "Подпор 4 м столба жидкости на всасе ")
plt.plot(Qs_pl_pump, Hs_pl_pump, 'o', label = "Найденные с помощью brentq решения")
plt.xlabel('$Q, м^3/ч$'); plt.ylabel('$H, м$')
plt.xlim(Qs[0],Qs[-1]); plt.ylim(Hs[-1],Hs[0])
plt.xlim(0, 4.6); plt.ylim(-4, 20)
plt.xticks(np.arange(0, 4.61, 0.2)); plt.yticks(np.arange(-4, 21, 2))
plt.legend(); plt.grid()

Инженерные расчёты на Python, С.В. Медведев, 2020
Использование Python и Jupyter Notebook для инженерных расчётов, С.В. Медведев, 2020