МОДЕЛИРОВАНИЕ РАСПРЕДЕЛЕНИЯ ДАВЛЕНИЯ ПРИ ГИДРАВЛИЧЕСКОМ УДАРЕ В ТРУБОПРОВОДЕ С ИСПОЛЬЗОВАНИЕМ COMSOL MULTIPHYSICS

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

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

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

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

представлен очень детальный обзор теории и практики гидравлического удара на 2005 год. В последующие годы в литературе имеются по крайней мере два разных подхода к моделированию трения на внутренней поверхности трубы во время моделирования гидравлического удара. Первая группа состоит из моделей, основанных на мгновенных изменениях в локальной и конвективной производных скорости, а вторая группа — модели, основанные на интеграле свертки и полной истории течения. Более популярные модели — это модели из первой группы, их использование требует эмпирических коэффициентов. Вторая группа все еще недооценена, даже если основана на хороших теоретических законах и не требует любых эмпирических коэффициентов. Это несомненно связано со сложностью вычислений интеграла свертки. А в последние несколько лет для примера представлены в хронологическом порядке следующие работы. В статье проведено численное исследование кавитационного течения при гидравлическом ударе. Анализ нестационарных моделей трения, используемых в техническом (инженерном) программном обеспечении для анализа гидравлического удара, выполнен в программе WANDA . Точная и эффективная схема, включающая нестационарное трение для переходного течения (потока) в трубе, предложена в . Эффекты времен и законов закрытия в трубопроводе шарового клапана при гидравлическом ударе рассмотрены в работе . Моделирование гидравлического удара, используя упрощенную основанную на свертке нестационарную модель трения, предлагается в . В этой работе предложено новое улучшенное эффективное решение интеграла свертки (применяя интегральное преобразование Лапласа), которое характеризуется использованием упрощенной весовой функции, состоящей только из двух экспоненциальных членов. Такой подход значительно содействует численным расчетам базовых параметров течения (давления и скорости). Далее представлен энергетический анализ квазидвумерной модели трения для моделирования переходных течений (потоков) в вязкоупругих трубах. Гидродинамика сглаженных (гладких) частиц с нестационарной моделью трения для течения в трубах при гидравлическом ударе рассмотрена в . Источники затухания при гидравлическом ударе и сравнение различных методов для моделирования закрытия клапана при гидравлическом ударе с программой CFD представлены в

[10]

, соответственно. А новые успехи в проблемах гидравлического удара (обзор) подробно изложены в статье .

Физическое моделирование гидравлического удара в реальных условиях очень трудное. Из-за размеров систем трубопроводов проведение исследования в реальных масштабах невозможно. Однако благодаря развитию численных методов изучение гидравлического удара и его эффектов может быть представлено с использованием программы математического моделирования, например программного комплекса Comsol Multiphysics.

Поэтому рассмотрим краевую задачу о распространении скачка давления в трубопроводе с начальной скоростью жидкости u0 в исходном состоянии. Допустим, что в начальный момент времени t = 0 в сечении x = 0 и в трубопроводе начальное давление p0, а в сечении x = l произошло скачкообразное изменение давления ΔPJ (например, при мгновенном закрытии клапана), которое действует в течение всего времени процесса вплоть до установления стационарного состояния. Требуется найти распределение давления при гидравлическом ударе по длине трубопровода во времени

, решая численно систему двух нелинейных уравнений гидродинамики.

Исходя из изложенного, цель работы:

1) использовать Comsol Multiphysics

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

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

Полная математическая постановка задачи по определению распределения давления в трубопроводе при гидравлическом ударе в виде системы двух нелинейных уравнений гидродинамики (уравнения непрерывности и уравнения импульсов), используя COMSOL MULTIPHYSICS, будет

Здесь p(x, t) — давление; x — продольная координата; t — время;

— вектор тангенциальной скорости жидкости ; — дифференциально-векторный оператор набла; = const — плотность жидкости; A — площадь поперечного сечения трубы; fD — коэффициент трения Дарси (или коэффициент гидравлического сопротивления); d — внутренний диаметр трубы; — вектор ускорения свободного падения; p0 – начальное давление в трубе и в точке x = 0; u0 – начальная скорость жидкости в трубе; — увеличение (скачок) давления, приложенного в точке x = l и действующего в течение всего времени процесса вплоть до установления стационарного состояния (); l — длина трубопровода.

Скорость звука в капельной упругой жидкости, текущей в трубе с упругими стенками, т.е. скорость распространения ударной волны — волны гидравлического возмущения c дается выражением

, :

где k — модуль упругости жидкости;

— толщина стенки трубы; E — модуль упругости материала стенки трубы.

Представленная формула для определения скорости звука в капельной упругой жидкости (c) справедлива при u/c < 1 и

, где u — скорость жидкости.

Дополнительно (только при определенных условиях — модель коэффициента трения Стокса

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

где

— фундаментальная формула Н.Е. Жуковского

которая определяет увеличение (скачок) амплитуды давления гидравлического удара, например, при мгновенном закрытии клапана в конце трубопровода x = l в начальный момент времени t = 0, а u0 может быть и средней скоростью жидкости до закрытия клапана. В целом уравнение (6) математически описывает затухающие гармонические колебания

.

Величина

для ламинарного режима, учитывая формулу Пуазейля для коэффициента гидравлического сопротивления, приводится к виду ,

где

— кинематическая вязкость жидкости.

По физическому смыслу коэффициент

(8) определяет коэффициент затухания колебаний (6) .

Для примера в

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

Круговая или циклическая частота затухающих колебаний

, принимая во внимание (6), находится по формуле

где безразмерный параметр For = const, характеризующий гидравлическое сопротивление жидкости, учитывающий ее вязкость, скорость звука в ней, коэффициент гидравлического сопротивления, а также длину трубопровода, определяется так

Тогда согласно формуле

, период затухающих колебаний равен

При необходимости можно определить логарифмический декремент затухания

и добротность колебательной системы .

В соответствии с видом функции (6) затухающие колебания можно рассматривать как гармонические колебания частоты

с амплитудой, изменяющейся по закону

Таким образом, для описания зависимости амплитуды давления от времени в определенном сечении трубопровода предложена более общая формула (12), которая включает фундаментальную формулу Н.Е. Жуковского

(7), описывающую увеличение (скачок) амплитуды давления гидравлического удара, например, при мгновенном закрытии клапана в конце трубопровода x = l в начальный момент времени t = 0.

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

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

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

(8), зависимости амплитуды давления от времени в определенном сечении трубопровода (12), периода колебаний ударной волны давления T (11).

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

: u0 = 14,8573 м/с — начальная скорость жидкости в трубе; p0 = 1,0133ּ⋅105 Па — начальное давление в трубе и в точке x = 0; = 1100 кг/м3 — плотность жидкости; k = 1100 МПа — модуль упругости жидкости; = 0,006 м — толщина стенки трубы; d = 0,207 м — внутренний диаметр трубы; E = 2⋅105 МПа — модуль упругости материала стенки трубы; = 5,75⋅10-6 м2/с — кинематическая вязкость жидкости; l = 200 м — длина трубопровода (или расстояние, где к трубопроводу присоединен, например, клапан, регулятор расхода жидкости или регулятор давления); e = 0,15 мм — шероховатость внутренней поверхности трубы; коэффициент трения Дарси fD (или коэффициент гидравлического сопротивления) для примера — модель коэффициента трения Черчиля (англ. Churchill) ,

где

В уравнениях (13)–(15) Re — это безразмерное число Рейнольдса, которое определяется так

Для ньютоновских однофазных жидкостей модель Черчиля для коэффициента трения Дарси fD может быть использована для полной области значений Re (ламинарной, переходной и турбулентной) и полной области e/d.

По представленным данным скорость звука в капельной упругой жидкости, текущей в трубопроводе с упругими стенками будет равна c = 916,7948 м/с (задача (1)–(5)), величина

согласно формуле (8) для ламинарного режима будет равняться a = 0,0021 1/с, а безразмерный параметр For, используя формулу (10), будет равен For = 1,1395⋅106.

Далее для представленной вязкой жидкости в трубопроводе проведено численное решение системы двух нелинейных уравнений гидродинамики (уравнения непрерывности (1) и уравнения импульсов (2)) с начальными (3) и граничными (4) условиями методом конечных элементов, используя модель коэффициента трения Черчиля (13)–(16). Поэтому на рисунке 1 представлено распределение давления p-p0 от времени t (уравнения (1)–(4)) при For = 1,1395⋅106 для x = l (или по-другому в конце трубопровода) в виде шаговой, ступенчатой функции. Кроме этого, на рисунке 1 в виде горизонтальной линии представлена величина скачка давления по формуле Н.Е. Жуковского (7) при t = 0 и здесь также ясно виден процесс затухания (уменьшения) амплитуды давления со временем.

Поэтому расстояние c⋅dt должно быть меньше, чем типичный размер сетки dx. При численном решении в Comsol Multiphysics было использовано условие Куранта-Фридрикса-Леви CFL = c⋅dt/dx = 0,1. Это числовое условие требует, чтобы размер временного шага был меньше времени прохождения волной пространственного шага сетки dx и оно необходимо для согласования шагов по времени и по пространственной переменной под технические возможности компьютера. Таким образом, выполняется необходимое условие устойчивости разностных схем. Осцилляции Гиббса могут быть ограничены дополнением некоторого количества высокочастотного демпфирования в

– обобщенный решатель, зависимый от времени, что и было сделано в модели и описано в инструкциях Comsol Multiphysics .

Кроме модели Черчиля COMSOL MULTIPHYSICS имеет еще несколько встроенных моделей коэффициентов трения для ньютоновских и неньютоновских жидкостей

. Например, для ньютоновской жидкости в ламинарном режиме (Re<2000) имеется коэффициент трения Дарси fD независимый от шероховатости внутренней поверхности трубы, который дается моделью (формулой) Стокса (англ. Stokes)

Решение нелинейных уравнений (1), (2) возможно лишь путем численного интегрирования. Для упрощения нелинейного уравнения (2) И.А. Чарный

для ламинарного режима предложил способ линеаризации второго слагаемого в правой его части, приняв множитель fD u/2d постоянным и равным его среднему значению по длине трубы и времени, используя модель Стокса (17), (16)

Таким образом, из уравнения (18) после соответствующих сокращений получаем формулу (8) для определения величины

. Тогда используя (18), уравнение (2) становится линейным. В этом случае с целью упрощения получения аналитического решения уравнения (1) и (2) сводятся к одному гиперболическому уравнению относительно, например давления

[15]

, что позволяет предложить упрощенный вариант для первого приближенного аналитического описания в трубопроводе основных физических характеристик гидравлического удара в виде уравнения (6), (7) и соответствующих формул (8)–(12) .

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

), используя модель коэффициента трения Стокса (17), (16) в COMSOL MULTIPHYSICS . Поэтому на рисунке 4 представлено соответствующее распределение давления p-p0 от времени t (уравнения (1)–(4)) при For = 1,1395⋅106 для x = l в виде шаговой, ступенчатой функции.

По числовым данным для задачи ((1)–(4), модель коэффициента трения Стокса) (рис. 4) приближенно определен период колебаний T = 0,88 с, а для уравнения (6) по формуле (11) аналитически (просто) при For = 1,1395⋅106 период затухающих колебаний T = 0,8726 с. В результате, при сравнении на рисунке 4 численного решения уравнений (1)–(4) (ступенчатая синяя линия) с первым приближенным аналитическим решением более простого (упрощенного, «инженерного») варианта (затухающая косинусоида — красная линия) для периода колебаний получена относительная погрешность ε = 0,85%, а также практически получено визуальное совпадение с предсказанной (уравнение (12)) амплитудой давления гидравлического удара (p-p0)1max от времени t (верхняя экспоненциальная зеленая линия).

Таким образом, подтверждается справедливость более простого (упрощенного, «инженерного») варианта для первого приближенного аналитического описания в трубопроводе (уравнения (6), (12)) основных физических характеристик гидравлического удара для ламинарного режима, если возможно использовать линеаризованную модель в виде одного дифференциального уравнения в частных производных второго порядка гиперболического типа.

Для математического описания гидравлического удара в трубопроводе использована модель в виде системы двух нелинейных уравнений гидродинамики (уравнения непрерывности (1) и уравнения импульсов (2)) с начальными (3) и граничными (4) условиями.

Если возможно для ламинарного режима использовать линеаризованную модель для описания гидравлического удара в виде одного уравнения гиперболического типа, предлагается более простой (упрощенный, «инженерный») вариант для первого приближенного аналитического описания (6) в трубопроводе основных физических характеристик гидравлического удара (формулы): например, коэффициента затухания a (8), зависимости амплитуды давления от времени в определенном сечении трубопровода (p-p0)1max (12), периода колебаний ударной волны давления T (11).

Для конкретной вязкой жидкости в трубопроводе проведено математическое (численное) моделирование распределения давления при гидравлическом ударе в трубопроводе методом конечных элементов, используя COMSOL MULTIPHYSICS и две модели коэффициента трения: Черчиля и Стокса.

Проведено сравнение (рис. 4) численного решения системы двух уравнений гидродинамики (1)–(4), используя модель коэффициента трения Стокса, методом конечных элементов с первым приближенным аналитическим решением (6) для более простого (упрощенного, «инженерного») варианта.

Для описания зависимости амплитуды давления от времени в определенном сечении трубопровода предложена более общая формула (12), которая включает фундаментальную формулу Н.Е. Жуковского

(7).

Получено совпадение максимальной величины давления полного решения p-p0 от времени t (уравнения (1)–(4) для модели коэффициента трения Стокса) с предсказанной амплитудой давления (p-p0)1max (12) от времени.