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

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

Введение

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

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

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

$\$\$display\$\$Delta\; p=rho\; cdot\; ccdot\; Delta\; v\$\$display\$\$$

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

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

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

1. Одномерная модель гидравлической линии в распределённых параметрах

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

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

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

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

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

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

$\$\$display\$\$frac\{partial\; v\}\{partial\; t\}+frac1rhofrac\{displaystylepartial\; p\}\{displaystylepartial\; x\}+vfrac\{displaystylepartial\; v\}\{displaystylepartial\; x\}=F+G\$\$display\$\$$

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

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

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

$\$\$display\$\$frac\{partial\; v\}\{partial\; t\}+frac1rhofrac\{displaystylepartial\; p\}\{displaystylepartial\; x\}=F+G\$\$display\$\$$

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

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

$\$\$display\$\$frac\{dv\}\{dt\}=F+G-frac1rhofrac\{displaystyletriangle\; p\}\{displaystyletriangle\; x\}\$\$display\$\$$

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

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

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

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

2. Метод характеристик

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

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

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

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

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

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

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

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

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

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

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

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

$\$\$display\$\$(v\_\{i,;j+1\}-v\_c)+frac1\{rho\; c\}(p\_\{i,;j+1\}-p\_c)-(F+G)triangle\; t=0\$\$display\$\$$

$\$\$display\$\$(v\_\{i,;j+1\}-v\_f)-frac1\{rho\; c\}(p\_\{i,;j+1\}-p\_f)-(F+G)triangle\; t=0\$\$display\$\$$

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

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

3. Эксперимент

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

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

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

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

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

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

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

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

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

Литература

1. В отечественной литературе метод характеристик для инженерного применения лучше всего описан в книге «Гидромеханика», Д. Н. Попов, С. С. Панаиотти, М.В. Рябинин.
2. В своей публикации
   Pipeline simulation by the method of characteristics for calculating the pressure pulsation of a high-pressure water plunger pump

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

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

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.