Метод характеристик для волнового уравнения на

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

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

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

Введение

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

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

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

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

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

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

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

где — плотность, — скорость, — давление, — потери на трение, — перепад давлений, вызванный гравитационной силой.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Метод характеристик при решение задачи коши для уравнений гиперболического типа

Стерлитамакский филиал Башкирский государственный университет

NovaInfo58, с. 11-15
Опубликовано 25 января 2017
Раздел: Физико-математические науки
Просмотров за месяц: 77
CC BY-NC

Аннотация

В статье рассматривается решение задачи Коши для уравнения гиперболического типа. Продемонстрировано решение данного уравнения методом характеристик.

Ключевые слова

Текст научной работы

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

где u=u(x,y,z,t), a — скорость распространения волны в данной среде. В одномерном случае это уравнение примет вид

которое является уравнением вынужденных колебаний однородной струны [1, 12].

В одномерном случае рассмотрим уравнение струны [2, 26]:

Задача Коши: Найти решение u(x,y) данного уравнения, удовлетворяющее начальным условиям:

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

Приведем уравнение (1) к каноническому виду. Для этого составим уравнение характеристик

где A=0, 2B=e y , C=-1. Вычислим D=B^2-AC=\frac><4>>0

. Следовательно, уравнение (1) является уравнением гиперболического типа.

Подставляя в уравнение характеристик наши значения, получим:

Волновое уравнение

Одним из наиболее распространенных в инженерной практике уравнений с частными производными второго порядка является волновое уравнение, описывающее различные виды колебаний. Поскольку колебания — процесс нестационарный, то одной из независимых переменных является время t. Кроме того, независимыми переменными в уравнении являются также пространственные координаты х, у, z. В зависимости от их количества различают одномерное, двумерное и трехмерное волновые уравнения.

Одномерное волновое уравнение – уравнение, описывающее продольные колебания стержня, сечения которого совершают плоскопараллельные колебательные движения, а также поперечные колебания тонкого стержня (струны) и другие задачи. Двумерное волновое уравнение используют для исследования колебаний тонкой пластины (мембраны). Трехмерное волновое уравнениеописывает распространение волн в пространстве (например, звуковых волн в жидкости, упругих волн в сплошной среде и т.п.).

Рассмотрим одномерное волновое уравнение, которое можно записать в виде

(2.63)

Для поперечных колебаний струны искомая функция U(x,t) описывает положение струны в момент t. В этом случае а2 = Т/ρ, где Т — натяжение струны, ρ — ее линейная (погонная) плотность. Колебания предполагаются малыми, т.е. амплитуда мала по сравнению с длиной струны. Кроме того, уравнение (2.63) записано для случая свободных колебаний. В случае вынужденных колебаний в правой части уравнения добавляют некоторую функцию f(x,t), характеризующую внешние воздействия, при этом сопротивление среды колебательному процессу не учитывается.

Простейшей задачей для уравнения (2.63) является задача Коши: в начальный момент времени задаются два условия (количество условий равно порядку входящей в уравнение производной по t):

(2.64)

Эти условия описывают начальную форму струны и скорость ее точек .

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

(2.65)

Рассмотрим некоторые разностные схемы для решения задачи (2.63)-(2.65). Простейшей является явная трехслойная схема типа крест (шаблон показан на рис. 2.21). Заменим в уравнении (2.63) вторые производные искомой функции Uпо tи х их конечно-разностными соотношениями с помощью значений сеточной функции в узлах сетки :

Рис. 2.21. Шаблон явной схемы

Отсюда можно найти явное выражение для значения сеточной функции на ( j + 1)-ом слое:

(2.66)

Здесь, как обычно в трехслойных схемах, для определения неизвестных значений на (j + 1)-ом слое нужно знать решения на j-ом и (j — 1)-ом слоях. Поэтому начать счет по формулам (2.66) можно лишь для второго слоя, а решения на нулевом и первом слоях должны быть известны. Их находят с помощью начальных условий (2.64). На нулевом слое имеем

(2.67)

Для получения решения на первом слое воспользуемся вторым начальным условием (2.64). Производную заменим конечно-разностной аппроксимацией. В простейшем случае полагают

(2.68)

Из этого соотношения можно найти значения сеточной функции на первом временном слое:

(2.69)

Отметим, что аппроксимация начального условия в виде (2.68) ухудшает аппроксимацию исходной дифференциальной задачи: погрешность аппроксимации становится порядка , т.е. первого порядка по τ, хотя сама схема (2.66) имеет второй порядок аппроксимации по hи τ. Положение можно исправить, если вместо (2.69) взять более точное представление:

(2.70)

Вместо нужно взять . А выражение для второй производной можно найти с использованием исходного уравнения (2.63) и первого начального условия (2.64). Получим

Тогда (2.70) примет вид:

(2.71)

Разностная схема (2.66) с учетом (2.71) обладает погрешностью аппроксимации порядка

При решении смешанной задачи с граничными условиями вида (2.65), т.е. когда на концах рассматриваемого отрезка заданы значения самой функции, второй порядок аппроксимации сохраняется. В этом случае для удобства крайние узлы сетки располагают в граничных точках (х0=0, xI = l). Однако граничные условия могут задаваться и для производной.

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

(2.72)

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

Рассмотренная разностная схема (2.66) решения задачи (2.63) — (2.65) условно устойчива. Необходимое и достаточное условие устойчивости:

(2.73)

Следовательно, при выполнении этого условия и с учетом аппроксимации схема (2.66) сходится к исходной задаче со скоростью O(h2+τ2). Данная схема часто используется в практи-ческих расчетах. Она обеспечивает приемлемую точность получения решения U(x,t), которое имеет непрерывные производные четвертого порядка.

Рис. 2.22. Алгоритм решения волнового уравнения

Алгоритм решения задачи (2.63)-(2.65) с помощью данной явной разностной схемы приведен на рис. 2.22. Здесь представлен простейший вариант, когда все значения сеточной функции, образующие двумерный массив, по мере вычисления хранятся в памяти компьютера, а после решения задачи выводятся результаты. Можно было бы предусмотреть хранение решения лишь на трех слоях, что сэкономило бы память. Результаты в таком случае можно выводить в процессе счета (см. рис. 2.13).

Существуют и другие разностные схемы решения волнового уравнения. В частности, иногда удобнее использовать неявные схемы, чтобы избавиться от ограничений на величину шага, налагаемых условием (2.73). Эти схемы обычно абсолютно устойчивы, однако алгоритм решения задачи и программа для компьютера усложняются.

Построим простейшую неявную схему. Вторую производную по tв уравнении (2.63) аппроксимируем, как и ранее, по трехточечному шаблону с помощью значений сеточной функции на слоях j— 1, j, j + 1. Производную до х заменяем полусуммой ее аппроксимации на (j + 1)-ом и (j— 1)-ом слоях (рис. 2.23):

Рис. 2.23. Шаблон неявной схемы

Из этого соотношения можно получить систему уравнений относительно неизвестных значений сеточной функции на (j+ 1)-ом слое:

(2.74)

Полученная неявная схема устойчива и сходится со скоростью . Систему линейных алгебраических уравнений (2.74) можно, в частности, решать методом прогонки. К этой системе следует добавить разностные начальные и граничные условия. Так, выражения (2.67), (2.69) или (2.71) могут быть использованы для вычисления значений сеточной функции на нулевом и первом слоях по времени.

При двух или трех независимых пространственных переменных волновые уравнения принимают вид

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