Документация Engee
Notebook

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

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

Рассмотрим трубопровод длиной 10 метров с круглым поперечным сечением и диаметром 0.1 метра. Массовый расход в трубе постоянен и составляет 0.4 кг/с. Вода поступает из резервуара при температуре 293.15 К и при атмосферном давлении. Изначально эффекты динамической сжимаемости жидкости не учтены.

image.png

Зададим массовый расход, длину и диаметр трубы.

In [ ]: 
mdot = 0.4;       # [кг/с] Массовый расход   
pipe_length = 10; # [м] Длина трубы 
D = 0.1;          # [м] Диаметр трубы

Падение давления из-за перепада высот

Как будет вести себя модель этой трубы если перепад высоты между ее входом и выходом будет равен 3 метра? Пренебрежем влиянием трения и рассчитаем падение давления аналитически, а потом используем модель чтобы проверить результат.

image.png

Перепад давления в трубе — сумма потерь давления из-за трения о стенки трубы и из-за перепада высот.

$$p_A - p_B = \Delta p_f + \Delta p_h $$

где $p_h$ - падение давления из-за изменения высоты по трубе, $p_f$ — падение давления из-за трения о стенки трубы. Поскольку мы пренебрегаем трением, можно считать $\Delta p_f$ равным нулю и использовать определить $\Delta p_h$ через уравнение Бернулли.

Уравнение Бергулли связывает скорость, давление и высоту жидкости в разных точках гидродинамической системы, и мы используем его чтобы рассчитать свойства жидкости в точках A и B модели трубы:

$$p_A + \frac{1}{2} \rho \nu_A^2 + \rho g h_A = p_B + \frac{1}{2} \rho \nu_B^2 + \rho g h_B $$

где $p$ — давление в точке A или B, $\rho$ — плотность жидкости, $\nu$ — скорость потока в точке A или B, $g$ — ускорение свободного падения, $h$ — высота жидкости в точке A или B.

При моделировании несжимаемой жидкости масса сохраняется и $\dot m_a + \dot m_b = 0$, где $\dot m_a$ - расход в точке A, $\dot m_b$ — расход в точке B.

Расход — произведение плотности, скорости и площади поперечного сечения трубы. Если предположить, что плотность жидкости и площадь поперечного сечения трубы постоянны по всей ее длине, то в силу сохранения расхода $\dot m_a^2 = \dot m_b^2$ и сохранения импульса можно использовать равенство $\nu_a^2 = \nu_b^2$, благодаря чему термины с этой переменной по обеим сторонам уравнения Бернулли можно исключить и выражение можно сократить до:

$$p_A - p_B = \Delta p_h = \rho g \left ( h_B - h_A \right ) = \rho g \Delta z $$

где $\Delta z$ — перепад высоты между точками A и B.

Выполним аналитический расчет ожидаемой потери давления из-за перепада высоты. Поменять перепад высот в модели можно через изменение переменной delta_z:

In [ ]: 
delta_z = 3;     # [м] Перепад высоты
g = 9.81;        # [м/с^2] Ускорение свободного падения
rho = 998;       # [кг/м^3] Плотность воды при 293.15 К и 101325 Па

deltaP_h = rho*g*delta_z   # [Па] Аналитический асчет падения давления из-за перепада высот
Out[0]:
29371.140000000003

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

image.png

In [ ]: 
# engee.close_all()
model = engee.open("$(@__DIR__)/PressureLossAndFlowRateInTLPipe_Insulated.engee")
delta_z = 3;
engee.run(model);

Построим график и изучим результаты расчета:

In [ ]: 
t = pressure_drop.time;
delta_P = pressure_drop.value;
plot( t, delta_P, ylimits = (2.93e4,2.94e4),
      xlabel="Время, с", ylabel="Перепад давления, 10⁴ Па",
      yformatter = x->"$(x/10000)", legend=false, lw=2 )
Out[0]:

Полностью отклюить трение в блоке "Труба" невозможно, поэтому результаты включают потери на трение, но они существенно меньше чем потеря из-за перепада высоты. Результаты показывают постоянную потерю давления, составляющую около $2,938 \cdot 10^4$ Па, что согласуется с аналитическими результатами (ошибка в три сотых от процента).

Потери давления из-за трения

Потери на трение, как правило, гораздо меньше, чем потери из-за перепада высоты. Однако, если перепад высоты $\Delta z$ равен нулю, потери на трение станут основной причиной потери давления в трубе.

Рассмотрим латунную трубу,у которой шероховатость поверхности внутренней стенки равна $1,52 \cdot 10^{-6}$ м. Проведем аналитический расчет потерь давления из-за трения в трубе, а затем проверим модель Engee на согласованность с результатами.

image.png

Потери давления из-за трения прежде всего зависят от числа Рейнольдса жидкости. Число Рейнольдса — безразмерная величина, используемая в гидродинамике для прогнозирования характера течения. Оно представляет собой отношение инерционности жидкости к ее вязкости. При низких числах Рейнольдса жидкость, как правило, имеет скорее ламинарный поток (инерционность намного превосходит вязкость), с увеличением числа Рейнольдса поток становится более турбулентным.

Число Рейнольдса обычно выражается следующим образом:

$$Re = \frac{\rho \nu D}{\mu} $$

где:

  • $\nu$ — скорость потока,

  • $\mu$ — динамическая вязкость жидкости,

  • $D$ — диаметр трубы.

Используя выражение для расхода $\dot m = \rho \nu A$ можно переписать выражение:

$$Re = \frac{\dot m D}{\mu A} $$

где $A$ — площадь поперечного сечения трубы.

Рассчитаем число Рейнольдса для воды в трубе:

In [ ]: 
S = pi/4 * D^2;      # [м^2] Площадь поперечного сечения трубы 
mu = 0.00102;        # [Па*с] Динамическая вязкость воды при 293 К и 0,101325 МПа

Re = mdot*D/(S*mu)   # число Рейнольдса
Out[0]:
4993.096253863384

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

$$\Delta p_f = \frac{f L}{2 \rho D A^2} \dot m | \dot m | $$

где $L$ — длина трубы, $f$ — коэффициент трения Дарси.

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

$$f = \left [ -1.8 log_{10} \left ( \frac{6.9}{Re} + \left ( \frac{\epsilon}{3.7 D} \right )^{1.11} \right ) \right ]^{-2} $$

где $\epsilon$ — мера абсолютной шероховатости внутренней поверхности.

In [ ]: 
epsilon = 1.52e-6;    # [м] Абсолютная шероховатость внутренней поверхности для коммерческих латунных, свинцовых, медных или пластиковых труб

f = 1 / (-1.8 * log10(6.9/Re + (epsilon/3.7)^1.11))^2;          # Коэффициент трения Дарси 
deltaP_f = f * pipe_length * mdot * abs(mdot) / (2*rho*D*S^2)   # [Па] Потеря давления из-за трения
Out[0]:
4.905190947467246

Запустим модель PressureLossAndFlowRateInTLPipe_Insulatedмодель и проверим величину потери давления из-за трения.

In [ ]: 
model = engee.open("$(@__DIR__)/PressureLossAndFlowRateInTLPipe_Insulated.engee")
delta_z = 0;       # [м] Высота трубы
engee.run(model);

Построим график результатов и изучим их.

In [ ]: 
t = pressure_drop.time;
delta_P = pressure_drop.value;
plot( t, delta_P, ylimits = (4.85,4.95),
      xlabel="Время, с", ylabel="Перепад давления, 10⁴ Па",
      legend=false, lw=2 )
Out[0]:

Модель показывает постоянное падение давления 4.878 Па, что соотносится с аналитическим результатом 4.905 Па (ошибка около половины процента).

Мы можем оценить влияние материала стенок трубы с разными коэффициентами потерь на трение на потерю давления:

In [ ]: 
Material  = ["Латунь/Свинец/медь/Пластик", "Сталь или железо (кованое)", "Чугун", "Чугун и ржавчина"];
Roughness = [1.52e-6, 4.6e-5, 2.59e-4, 1.5e-3]; # [м] Абсолютная шероховатость внутренней поверхности

f         = 1 ./ (-1.8 .* log10.(6.9/Re .+ (Roughness./3.7).^1.11)).^2;    # Коэффициент трения Дарси
Pressure_Loss  = f .* pipe_length .* mdot .* abs(mdot) ./ (2*rho*D*S^2);   # [Па] Перепад давления из-за трения
using DataFrames
T = DataFrame("Материал" => Material, "Шероховатость" => Roughness, "Перепад давления" => Pressure_Loss)
Out[0]:
4×3 DataFrame
RowМатериалШероховатостьПерепад давления
StringFloat64Float64
1Латунь/Свинец/медь/Пластик1.52e-64.90519
2Сталь или железо (кованое)4.6e-54.90897
3Чугун0.0002594.93133
4Чугун и ржавчина0.00155.08436

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

Сохранение массы и сжимаемость жидкости

Теперь исследуем, как на расход в трубе скажутся учёт эффекта динамической сжимаемости и добавление источника тепла в модель.

Если отключить параметр "Динамическая сжимаемость" (Fluid dynamic compressibility), то модель трубопровода не будет учитывать динамические эффекты в массе жидкости внутри трубы и расход на выходе всегда будет равен расходу на входе:

$$\dot m_a + \dot m_b = 0 $$

где $\dot m_a$ - расход через сечение A, $\dot m_b$ - расход через сечение B.

Поскольку рассматриваемая до сих пор модель устроена именно так, расход на обоих концах трубы будет одинаковым.

In [ ]: 
model = engee.open("$(@__DIR__)/PressureLossAndFlowRateInTLPipe_Insulated.engee")
delta_z = 0;       # [м] Высота трубы
engee.run(model);

Вот результаты расчета для модели без учета сжимаемости.

In [ ]: 
t = pressure_drop.time;
mdot_a_incompressible = mdot_a.value;
mdot_b_incompressible = mdot_b.value;
plot( t, [mdot_a_incompressible -mdot_b_incompressible],
      xlabel="Время, с", ylabel="Массовый расход, (кг/с)", lw=2,
      ls=[:solid :dash], label=["mdot_a" "- mdot_b"], ylimits=(-0.5,0.5) )
Out[0]:

Массовый расход в портах А и В равен по величине и противоположен по знаку.

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

$$\dot m_a + \dot m_b = V_p \left ( \frac{1}{\beta} \frac{dp}{dt} - \alpha \frac{dT}{dt} \right ) $$

где:

  • $V$ — объем жидкости в трубе,

  • $\rho$ — плотность,

  • $p$ — давление,

  • $T$ — температура,

  • $\beta$ — изотермический модуль объемной упругости,

  • $\alpha$ — изобарический коэффициент теплового расширения.

Поскольку труба идеально изолирована, а потеря давления из-за трения незначительна, $\frac{dp}{dt}$ и $\frac{\partial T}{\partial m}$. Оба значения можно пренебречь. В этом случае, поскольку правая часть уравнения баланса массы по-прежнему приблизительно равна нулю, труба ведёт себя так же, как и в случае несжимаемой среды, и массовый расход между портом A и портом B не меняется.

Включим динамическую сжимаемость жидкости в блоке "Трубопровод":

In [ ]: 
engee.set_param!( "PressureLossAndFlowRateInTLPipe_Insulated/Pipe (Advanced) (TL)",
   "dynamic_compressibility" => true )

Запустим модель и построим графики:

In [ ]: 
engee.run(model);

t = pressure_drop.time;
mdot_a_compressible = mdot_a.value;
mdot_b_compressible = mdot_b.value;
plot( t, [mdot_a_compressible -mdot_b_compressible],
      xlabel="Время, с", ylabel="Массовый расход, (кг/с)", lw=2,
      ls=[:solid :dash], label=["mdot_a" "- mdot_b"], ylimits=(-0.5,0.5) )
Out[0]:

Массовые расходы в портах А и В по-прежнему равны по величине и противоположны по знаку.

Чтобы пронаблюдать существенное различие в расходе, добавим учет теплоотдачи, подключив к модели источник тепла. Тогда компонент уравнения, содержащий термин $\frac{\partial T}{\partial m}$ достигнет заметного размера. Рассмотрим модель, в которой температура источника тепла, подключённого к трубе, изменяется линейно во времени.

image.png

image.png

Построим график температуры при наличии источника тепла в модели.

In [ ]: 
t = 0:0.1:120
temperature = 200 .+ 1.5t;
plot( t, temperature, xlabel="Время (с)", ylabel="Температура (К)")
Out[0]:

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

In [ ]: 
model = engee.open("$(@__DIR__)/PressureLossAndFlowRateInTLPipe_HeatSource.engee");
engee.set_param!( "PressureLossAndFlowRateInTLPipe_HeatSource/Pipe (Advanced) (TL)",
   "dynamic_compressibility" => true )
engee.run(model);

И построим график с результатами

In [ ]: 
t = pressure_drop.time;
mdot_a_compressible = mdot_a.value;
mdot_b_compressible = mdot_b.value;
plot( t, [mdot_a_compressible -mdot_b_compressible],
      xlabel="Время, с", ylabel="Массовый расход, (кг/с)", lw=2,
      ls=[:solid :dash], label=["mdot_a" "- mdot_b"], ylimits=(0.39,0.41) )
Out[0]:

Заключение

В начале моделирования температура источника равнялась 200 К, что меньше чем температура воды в трубе, которая составляет 293.15 К. После истечения 55 секунд, температура источника тепла превысила температуру воды, втекающей в трубу. В первой половине эксперимента скорость изменения температуры в трубе была отрицательной, то есть выполнялось $\dot m_a < \dot m_b$. После того, как температура источника тепла превысила температуру воды, поступающей в трубу, тепло начинает передаваться в трубу. В этом случае скорость изменения температуры стала положительной, что приводит к $\dot m_a > \dot m_b$.

Блоки, использованные в примере