Главная » Хабрахабр » [Из песочницы] Расчёт волновых процессов в гидравлической линии методом характеристик

[Из песочницы] Расчёт волновых процессов в гидравлической линии методом характеристик

В этой статье я расскажу про создание математической модели длинного трубопровода для CAE-программы SimulationX на языке Modelica. Привет, Хабр! Несмотря на то, что этот метод довольно старый, в рунете довольно мало информации о его применении для решения прикладных задач. Речь пойдёт о расчёте волновых процессов (пульсации давления, гидроудар и т.п.) в гидравлической линии методом характеристик.

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

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

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

С тех пор, формулу для расчёта скачка давления при быстром закрытии задвижки, во всём мире называют формулой Жуковского:
Гидроудар с необходимой для решения практических задач точностью описал в конце XIX века Николай Жуковский, решив таким образом проблему аварий на Московском водопроводе.

$\Delta p=\rho \cdot c\cdot \Delta v$

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

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

их расчёт требовал углубления в матан и численные методы, к которым на протяжении всей учёбы я относился со снисходительным пренебрежением. В своей практике я долго пытался отмахиваться от волновых процессов, т.к. Тем более к моему стыду, за меня уже начал писать модель трубопровода на С++ мой научный руководитель. Но когда однажды я своими глазами увидел, что стандартные советы (поставить везде РВД, гидроаккумулятор, организовать подпор на всасе насоса) не помогают ни избавиться от пульсаций на стенде, ни, тем более, приближают к пониманию происходящих процессов, пришлось-таки углубляться в матан.

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

$\frac{\partial t}+\frac{\partial(\rho v)}{\partial x}=0,$

$\frac{\partial(\rho v)}{\partial t}+\frac{\partial(\rho v^2+p)}{\partial x}=\rho\cdot(G+F),$

где $\rho$ — плотность, $v$ — скорость, $p$ — давление, $F$ — потери на трение, $G$ — перепад давлений, вызванный гравитационной силой.

интегрировать теперь нужно не только по времени $t$, но и по пространственной координате $x$.
В случае с жидкостями, можно ещё немного упростить себе жизнь, если переписать уравнения из консервативных переменных в примитивные переменные (скорость и давление):
Т.е.

$\frac{\partial p}{\partial t}+v\frac{\displaystyle\partial p}{\displaystyle\partial x}+\rho\;c^2\frac{\partial v}{\partial x}=0$

$\frac{\partial v}{\partial t}+\frac1\rho\frac{\displaystyle\partial p}{\displaystyle\partial x}+v\frac{\displaystyle\partial v}{\displaystyle\partial x}=F+G$

где $c$ — скорость звука.

Теперь, если принять, что скорость звука существенно больше скорости движения жидкости $v\ll c$ (что справедливо при отсутствии кавитации), то уравнения станут ещё немного проще:

$\frac{\partial p}{\partial t}+\rho\;c^2\frac{\partial v}{\partial x}=0$

$\frac{\partial v}{\partial t}+\frac1\rho\frac{\displaystyle\partial p}{\displaystyle\partial x}=F+G$

Чтобы решить эти уравнения, нужно тем или иным способом избавиться от дифференцирования по пространственной координате $x$. В лоб это можно сделать, если заменить пространственный дифференциал конечно-разностной схемой, а в случае с временем тогда просто перейти к полному дифференциалу, сказав, что в рамках одной ячейки, параметры состояния не зависят от координаты:

$\frac{dp}{dt}=-\rho\;c^2\frac{\triangle v}{\triangle x}$

$\frac{dv}{dt}=F+G-\frac1\rho\frac{\displaystyle\triangle p}{\displaystyle\triangle x}$

Теперь эти уравнения можно решить как обыкновенные дифференциальные уравнения, разбив длину трубы на множество конечных объёмов. Так это и делается, например, в пакете Simscape, в MATLAB Simulink, так проблема решалась до недавнего времени в SimulationX.
Что-то таким образом, конечно, можно посчитать, но сильно мешают возникающие при этом численные колебания:

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

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

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

Википедия по запросу “Метод характеристик” рекомендует:

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

Это как философский камень, только вместо превращения металлов в золото мы превращаем уравнения в частных производных в обыкновенные, и наоборот. Возникает вопрос: “как же применить это на практике?”, причём желательно более эффективно, чем это делали средневековые алхимики…

В нашем распоряжении в начальный момент времени имеется какое-то распределение давлений и скоростей по длине трубы. Для начала, разберёмся с постановкой задачи. Первым делом мы разобьём трубу на конечное число элементов и для каждой грани присвоим своё значение давления $p_i$ и скорости $v_i$.

Перенесёмся в пространство-время и расположим состояние трубы в будущем выше начального состояния: Интересует нас то, как изменятся значения в этих точках через момент времени $\triangle t$.

Рабоче-крестьянское объяснение заключается в том, что все изменения в трубе происходят со скоростью звука. Вот здесь-то нам и пригодятся “магические” характеристики! Иллюстрируется это следующим образом: Давление и скорость в текущий момент времени в точке $i$ зависят от давления и скорости в тех точках трубы, где звуковая волна была (бы) $\triangle t$ секунд назад.

Это и есть те самые характеристики, вдоль которых уравнения в частных производных превращаются в обыкновенные дифференциальные уравнения. Из какой-либо точки проводятся две симметричные линии, угол наклона которых определяется скоростью звука. Если мы назовём точки, в которых характеристики пересекаются с состоянием трубы в прошлом как $c$ и $f$, уравнения запишутся следующим образом:

$dv+\frac1{\rho c}dp-(F+G)dt=0\;при\;\frac{dx}{dt}=+c$

$dv-\frac1{\rho c}dp-(F+G)dt=0\;при\;\frac{dx}{dt}=-c$

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

Для этого шаг времени должен удовлетворять критерию Куранта — Фридрихса — Леви (CFL):
Важно учитывать, что эти точки всегда должны находиться в пределах соседних ячеек!

$\triangle t<\frac{\triangle x}c$

Теперь к этим уравнениям можно применить хоть самую простую разностную схему:

$(v_{i,\;j+1}-v_c)+\frac1{\rho c}(p_{i,\;j+1}-p_c)-(F+G)\triangle t=0$

$(v_{i,\;j+1}-v_f)-\frac1{\rho c}(p_{i,\;j+1}-p_f)-(F+G)\triangle t=0$

В получившейся системе из двух уравнений две неизвестных: давление $p_{i,\;j+1}$ и скорость $v_{i,\;j+1}$. Можно решать её численно, но нет особых проблем получить и аналитическое решение. Тогда, если принять постоянство скорости звука, получится полностью явная разностная схема.

Для закрепления приведу анимацию работы метода характеристик:

На самом деле...

В этом случае характеристики, строго говоря уже не будут прямыми линиями, а для того, чтобы найти давление, нужно уже будет знать скорость звука, которая от этого давления зависит. … скорость звука зависит от давления жидкости. схема будет уже неявной.
При создании модели, я принял допущение, что скорость звука слабо меняется от шага к шагу. Т.е. Чтобы быть уверенным в результате, модель лучше использовать при давлениях от 10 бар. Для жидкости это справедливо в случае низкого газосодержания и при отсутствии кавитации.

Мне представилась возможность окончательно довести до ума модель уже когда я начал работать в фирме ESI ITI GmbH в Дрездене. Как-то раз, я получил тикет в Helpdesk, где инженеры фирмы URACA жаловались, что не могут добиться сходимости с экспериментом с нашей “старой” трубой.
Они изготавливают водяные плунжерные насосы высокого давления, эдакие огромные “Керхеры” и хотели бы уметь предсказывать возможные резонансные эффекты, обусловленные в т.ч. волновыми процессами в трубопроводе. Проблема заключается в том, что у таких насосов, как правило, довольно мало плунжеров и работают они на низких оборотах (250-500 об/мин):

Из-за этого, а также из-за влияния сжимаемости жидкости, на выходе получается очень неравномерная подача:

Кроме того, у них уже были модели в SimulationX, где они учитывали динамику механической части насоса, упругости рамы, характеристику электродвигателя, поэтому интересно было бы посмотреть как на это повлияет трубопровод. Разрывы и нелинейности делают затруднительным процесс линеаризации и анализ модели в частотной области, а расчёты CFD для такой задачи — стрельба из пушки по воробьям.

Схема испытательного стенда довольно простая:

В начале трубопровода установлен датчик давления pd1, на расстоянии 22 метра от него — датчик давления pd2. Имеется простой трубопровод общей длиной примерно 30 метров. В конце трубопровода клапан, которым настраивается давление в системе.

Я предложил протестировать бета-версию своей модели, в результате в SimulationX была собрана такая модель:

Результаты даже меня приятно удивили:


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

Там пришлось использовать метод на основе метода Годунова, который в свою очередь базируется на решении задачи Римана о распаде произвольного разрыва, ну да об этом уже как-нибудь в другой раз… Этот опыт позволил оперативно запустить в релиз SimulationX новую модель гидравлической линии, а я погрузился в эту тему и не заметил как вместе со студентом-практикантом запилил и модель пневматической линии, где всё было гораздо интереснее.

  1. В отечественной литературе метод характеристик для инженерного применения лучше всего описан в книге «Гидромеханика», Д. Н. Попов, С. С. Панаиотти, М.В. Рябинин.
  2. В своей публикации

    Pipeline simulation by the method of characteristics for calculating the pressure pulsation of a high-pressure water plunger pump

    Uwe Grätz and Dipl.-Ing. «Dr.-Ing.(Rus) Maxim Andreev, Dipl.-Ing. IFK, March 19-21, 2018, Aachen, Germany, за текстом обращайтесь в личку (FH) Achim Lamparter», The 11th International Fluid Power Conference, 11.

    я более подробно рассмотрел проблемы стыковки метода характеристик и решателя ОДУ.

  3. У кого есть доступ к немецким библиотекам, лучший обзор методов решения гиперболических уравнений в применении к гидравлическим линиям, который я встречал, содержится в следующей диссертации: Beck, M., Modellierung und Simulation der Wellenbewegung in kavitierenden Hydraulikleitungen, Univ. Stuttgart, Germany, 2003.
  4. Классика жанра по гиперболическим уравнениям в целом: Randall J. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, Cambridge, United Kingdom, 2002.


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

Ваш email нигде не будет показан
Обязательные для заполнения поля помечены *

*

x

Ещё Hi-Tech Интересное!

[Из песочницы] LVEE — самая неформальная и душевная ИТ-тусовка

О событии Под Минском завершилась ежегодная конференция LVEE 2018 (Linux Vacation/Eastern Europe), организованная белорусской группой пользователей Linux еще в 2005 году. В мероприятии принимали участие более сотни людей из 7 стран, а программу составили 23 доклада, лайтнинги, воркшоп, круглый стол, ...

Губительная ошибка новичков в геймдеве

Перед началом любого дела необходимо составить план, сделать «пробы пера», одним словом — черновик. Именно это помогает определить стартовую точку и понять направление движения. Не хотите тратить тонны усилий впустую? Хотите делать быстрее и качественней остальных? 90% начинающих разработчиков этого ...