УДК 519.6:621
Метод контрольного объема для расчета гидравлических сетей
© О.В. Белова1, В.Ю. Волков2, А.П. Скибин3
1 МГТУ им. Н.Э. Баумана, Москва, 105005, Россия 2 ОАО ВНИИАЭС, Москва, 141980, Россия 3 ОАО ОКБ «ГИДРОПРЕСС», Московская область, Подольск, 142103, Россия
Решение задачи потокораспределения является неотъемлемой и крайне важной проблемой при моделировании работы нефте- и газопроводов, а также водопроводных и тепловых сетей (далее — гидравлические сети). Обычно методики расчета подобных систем базируются на первом и втором законах Кирхгофа, законе сохранения массового баланса и законе сохранения энергии. На этих же законах основаны и численные методы решения гидравлических сетей вплоть до современных программных продуктов. Авторами был совершен переход от методов, основанных на совместном решении законов Кирхгофа, к численному решению распределения потоков в гидравлической сети с помощью дискретизации уравнения неразрывности, которое приводится к разностному аналогу дифференциального уравнения второго порядка относительно давления. Данный метод успешно апробирован на примере решения нескольких тестовых задач.
Ключевые слова: гидравлическая сеть, метод контрольного объема, поле давления, расход, уравнение неразрывности.
Введение. Принято считать, что первые водопроводы появились приблизительно в 700-х гг. до н. э. Вода с горных вершин под действием силы тяжести по открытым каналам поступала в города, расположенные неподалеку. Одни из самых знаменитых водоводов древности — древнеримские водопроводы, или акведуки. Протяженность самого большого из них составляла более 90 км. Система трубопроводов строилась из камня, щебня, кирпича, грубого бетона, а трубы, как правило, были выдолблены из деревянных бревен или камня, однако встречались и образцы, изготовленные из глины и свинца. С ростом числа и сложности трубопроводов увеличивалась и сложность гидравлических систем, а также росло количество различной используемой арматуры, что привело к проблемам оценки распределения давления и определения скоростей в различных частях системы. Поиск метода анализа всей сети в целом привел к появлению отдельного раздела гидравлики, называемого «теория гидравлических сетей», или «теория гидравлических цепей» [1-3].
В настоящее время накопленные данные, рост компьютерных мощностей и применение современных измерительных приборов, имеющих возможность получать и обрабатывать данные в реальном времени, расширили границы использования теории гидравлических сетей. Так, помимо различных переходных режимов на станциях, гидравлические расчеты позволяют определить нелегальное использование воды, уменьшить затраты электроэнергии для насосных станций. Также расчеты и анализ переходных режимов позволяют уменьшить возможные разрывы трубопроводов и загрязнение воды. Одними из наиболее актуальных проблем для моделирования процессов потокораспределения являются пожары, при которых требуется определить изменение работы гидравлической сети при заборе жидкости из водяных сетей города. Многие современные программные продукты позволяют выполнять калибровку систем, чтобы помочь пользователю определить «тупиковые» ветки гидравлической сети и провести оптимизацию гидравлической сети для устранения ее недостатков.
Основы гидравлических расчетов были заложены в ХУШ-Х1Х вв. в работах Бернулли, Пуазейля, Рейнольдса. Однако и сейчас огромное количество работ посвящено решению вопросов задач потоко-распределения в гидравлических цепях, так как уравнения, описывающие даже изотермическое течение однофазной жидкости, нелинейны и задача их решения является нетривиальной и требует индивидуального подхода к каждой проблеме.
В общем виде для решения задачи потокораспределения в произвольной гидравлической сети должны быть записаны уравнения сохранения импульса и неразрывности для каждого элемента сети (трубы, насосы, клапаны и т. д.). Число уравнений, которые должны быть решены даже для малой гидравлической сети, является достаточно большим, что делает данный процесс решения крайне трудоемким.
Наиболее эффективными методами расчета потокораспределения в трубопроводных сетях стали итерационные, или так называемые «увязочные» методы, впервые предложенные в 30-х гг. XX в. Анд-рияшевым и Кроссом [1, 4, 5]. Данные методы являются покоординатными релаксационными методами последовательных приближений. «Увязка» сети с заданными сопротивлениями и напорами имеет своей целью найти такие значения расходов на всех участках и давлений в узлах, которые с наперед заданной точностью удовлетворяли бы законам Кирхгофа [6]. Основными недостатками «увязочных» ме-
тодов являются необходимость хорошего начального приближения и низкая скорость сходимости.
Этих недостатков лишен градиентный метод (Global Gradient Algorithm, GGA), предложенный Эзио Тодини в 1987 г. [2]. Данный алгоритм сочетает в себе быструю сходимость и высокую точность решения. Наиболее полно описание градиентного метода изложено в руководстве пользователя программного комплекса EPANET, написанного Льюисом Россманом [7].
Однако для гидравлических сетей со сложной топологией и сопротивлениями участков контура, отличающимися более чем на 10 порядков, используемые методы не обеспечивают сходимости итерационного процесса, особенно для сетей, на которых установлены регуляторы расхода или давления [3, 8]. Поэтому авторами предложен новый метод численного моделирования гидравлических сетей со сложной топологией, основанный на дискретизации уравнения неразрывности. Данный метод аналогичен методу расчета гидравлических сетей при помощи метода конечных элементов и описанного в работах [9, 10].
Постановка задачи. Рассмотрим постановку задачи потокорас-пределения для гидравлической сети, показанной на рис. 1. Гидравлическая сеть состоит из гидравлических связей, которые представляют собой каналы (трубы) заданной длины с постоянной площадью сечения, и узлов, в которых происходит объединение гидравлических связей.
1 2 3 4 5 6 7
♦ * Л --1 ■ 1 | < ■ < ■ • ь
8 9 10 11 12
Рис. 1. Схема гидравлической сети
На рис. 1 узлы отмечены цифрами и показаны в виде ромбов, а гидравлические связи представляют собой отрезки между узлами. При этом значения давления и источники массы задаются только в узлах, а коэффициенты гидравлического сопротивления и гидравлического трения — только для гидравлических связей.
Математическая модель. Движение жидкости в гидравлической сети считается стационарным и изотермическим. Тогда в общем случае описание изотермического движения жидкости происходит при
помощи системы уравнений Навье—Стокса для ламинарного движения и уравнений Рейнольдса для турбулентного.
--дP
ё1у(рмм - ц§гаё(м)) =--, (1)
дx
дР
div(pvu - ц grad(v)) = - —, (2)
dy dP
div(pwU - p,grad( w)) =--, (3)
dz
div(pu) = q. (4)
В данном случае в уравнениях (1-3) вместо давления p применяется полное давление P, определяемое выражением
P = Р - pg(eg, г), (5)
где р — плотность движущейся среды, кг/м3; r — радиус-вектор текущей точки пространства, м; (e , r) — скалярное произведение
единичного вектора силы тяжести на радиус-вектор, м; u — вектор скорости, м/c; u, v, w — компоненты скорости по оси х, y или z соответственно, м/c; x, у, z — координаты, м; д — динамическая вязкость среды, Па • с; g — ускорение свободного падения, м/с ; q — массовый источник, кг/м3 с.
Граничные условия на входе в гидравлическую сеть:
P = Pinlet или U = Uinlet > (6)
Граничные условия на выходе из гидравлической сети:
P = Poutlet или U = Uoutlet> (7)
где inlet и outlet — обозначения сечений на входе и выходе из гидравлической системы
Метод расчета. В применяемом методе решения системы нелинейных дифференциальных уравнений, составляющих математическую модель процессов в гидравлической сети (1-7), используется конечно-разностный вариант метода контрольного объема (МКО) на совмещенной сетке [11-13]. Для описания задачи стационарного движения в фрагментах гидравлической сети можно сделать следующее упрощение: если направить координатную ось OX по направлению оси гидравлической связи, то уравнения сохранения количества движения по осям OY и OZ рассматривать нецелесообразно. Поэтому из рассматриваемой системы уравнений Навье—Стокса можно
исключить уравнения (3) и (4). В данной постановке уравнения движения для компоненты вектора скорости, совпадающей с направлением движения среды в канале (1), рассматриваются как уравнения для определения поля скорости, а уравнение неразрывности (4) — как уравнение для определения поля давления.
Дискретные аналоги уравнений (1) и (4) получаются их интегрированием по каждому контрольному объему, на который разбита расчетная область. Особенностью разрабатываемого метода является непосредственная подстановка дискретного аналога для уравнения движения в уравнение неразрывности для получения дискретного аналога уравнения для определения поля давления.
Дз
■ф-
и
-О
" л1
рр
ех/2
65 /2
0Х
■ у;
Ре
2 -3
♦-1-*
е Е
№ V/ Р
5 1 Ре
б
Рис. 2. Контрольные объемы для уравнения движения (а) и уравнения неразрывности (б)
а
Рассмотрим интегральную формулировку уравнения движения (1), для чего проинтегрируем уравнение (3) по типичному контрольному объему (рис. 2, а) и применим теорему Остроградского— Гаусса:
| ё1у(рмм - [ц] егаё(м))dV = | (- — ) dV, (8)
V V { дх)
|((рии -[ц]¡^(и)),п^ = \{-^~]dV, (9)
5 V { дХ)
где п — единичный вектор внешней нормали к поверхности КО; 5 — площадь боковой поверхности КО.
После вычисления интегралов (9) дискретный аналог уравнения количества движения (8) принимает следующий вид:
{к + ^] иР = К№и№ - \JPdV, (10)
Р - Р
= (11)
и =-
V
25х дР
= -
V о.
ОХ
Рм, К1 £ дх
Р№
где К№ и КР — соответственно расходы через грани «V» и «Р» контрольного объема
К =| Ри^ = рЛ^ = рР =| Ри<*5 = РРиР5Р, (12)
Sw 5Р
где X — распределенный коэффициент гидравлического трения; Б — гидравлический диаметр, м; д — коэффициент гидравлического сопротивления; 5х — длина гидравлической связи, м.
Для получения дискретного аналога для определения давления проинтегрируем уравнение неразрывности (4) по контрольному объему, изображенному на рис. 2, б. Согласно формуле Остроградского—Гаусса, для левой части уравнения (4) имеем:
| ё1у(ри) dv = | ((ри), п =
V 8
= | ((ри), п + | ((ри), п + | ((ри), п)dS.
Следовательно, дискретный аналог уравнения неразрывности (4) имеет вид:
к - К+= е, (13)
где е — массовый источник в контрольном объеме, ассоциированном с узлом Р, кг/с.
Определим конвективные потоки ¥е и К№ на гранях ^ и е аналогично (11):
Ъ — (ри =-Р^А
дР Рв^в
дх РЕ 8хв
дР МА
дх РЯ 5х*
■(РЕ -РР), (14) (Р3 -Рр ). (15)
Введя обозначения
п = РАЯ± п = п (16)
пв с- , с- , с- , V1
5хе
получим дискретный аналог уравнения неразрывности в следующем виде:
аррр = ажр + аЕРЕ + + ь (17)
где
а8 = Пя ,
а — п ,
аЕ — пЕ,
ь —б,
ар — аш аЕ ^ а^.
Для решения системы уравнений (1-4) предлагается следующая итерационная процедура:
1. Вводятся предполагаемые поля скорости и давления.
2. Рассчитываются значения с1 для всех гидравлических связей по формуле (11).
3. Рассчитываются коэффициенты дискретного аналога поля давления (16) и определяются поля давления и градиентов давления. Определяются массовые потоки через грани КО.
4. Рассчитываются значения скорости для каждой гидравлической связи, используя коэффициенты дискретного аналога уравнения движения и полученный градиент поля давления.
5. Возврат к п. 2 до тех пор, пока не будет достигнута сходимость.
Тестирование метода расчета. В качестве одной из тестовых задач рассмотрена задача о моделировании потокораспределения участка тепловой сети, представленная в [14]. В данной работе задача решалась методом контурных расходов Андрияшева—Лобачева—Кросса [1, 3]. Расчетная и принципиальная схемы рассматриваемого участка приведены на рис. 3. Расчетный участок состоит из десяти связей, в одной из
которых задан постоянный напор, и восьми узлов, в двух из которых заданы расход и приток массы. В ходе гидравлического расчета требуется определить объемные расходы для каждого элемента теплосети. В данной схеме коэффициенты трения в гидравлических связях являлись величинами постоянными и не зависящими от числа Рейнольдса.
а
Рис. 3. Гидравлическая схема первой тестовой задачи: а — принципиальная; б — расчетная
Задача, расчетная схема которой приведена на рис. 1, была решена при помощи представленного выше метода контрольного объема. Исходные данные для тестового гидравлического расчета представлены в табл. 1.
Сравнение полученных объемных расходов для каждого элемента представлено в табл. 2.
Сравнение результатов, полученных при помощи метода контрольного объема, и результатов работы [14] определялось по следующей зависимости:
5 = 100 (Им-И МКО )
Няв/
где НКв^ — значения напоров в узлах гидравлической сети из [14]; НМКО — значения напоров в узлах, полученные с помощью МКО.
Таблица 1
Исходные данные для гидравлического расчета
№ связи Гидравлическое сопротивление Действующий напор, м
0 1е-12 1,5684е-8 115
1 0е+0 0,13190284 0
2 8,41е-6 0,13190284 0
3 8,41е-6 0,13190284 0
4 8,41е-6 0,13190284 0
5 8,41е-6 0,13190284 0
6 8,41е-6 0,13190284 0
7 8,41е-6 0,13190284 0
8 5,50е-4 8,62622615 0
9 5,50е-4 8,62622615 0
10 5,50е-4 8,62622615 0
№ узла Расход / Приток, мм3/ч
0 0
1 0
2 -100
3 100
4 0
5 0
6 0
7 0
8 0
Таблица 2
Объемный расход для каждого элемента теплосети
№ участка м3/ч м3/ч 5, %
1 100 100 0,00000
2 1264,21 1264,15 0,00438
3 759,31 759,27 0,00529
4 376,8 376,78 0,00637
5 -1164,21 -1164,15 0,00475
6 -759,31 -759,27 0,00529
7 -376,8 -376,78 0,00637
8 404,9 404,88 0,00376
9 382,51 382,49 0,00422
10 376,8 376,78 0,00637
Авторами проведено исследование скорости сходимости задачи от числа итераций. Изменение невязки по итерациям представлено на рис. 4.
1.Е+00
1,Е-01
1,ЕпШ
я 1.Е-03 й
№
1.Е-05 1.Е-06 1.Е-07
—•—МКО
6 г ю № итерации
Рис. 4. Скорость сходимости: МКО — метод контрольного объема; Loginov — данные статьи [14]
Для определения невязки Ф использовалось следующее выражение:
Е \аг1 - аг\
ф , =
т+1
I=1, N
Е
I=1, N
где От — объемный расход на 1-м элементе расчетной схемы на т-й итерации, м3/с.
Как видно из рисунка, скорость сходимости предлагаемого метода контрольного объема выше по сравнению с вариантом метода контурных расходов Андрияшева—Лобачева—Кросса, применямого в [14]. При этом можно отметить, что при численном решении с помощью МКО невязка уменьшается монотонно, тогда как для «увязочных» методов сходимость задачи является немонотонной.
В качестве другой тестовой задачи для проверки работоспособности предложенного метода контрольного объема была рассмотрена тестовая задача из статьи [15], описывающая потокораспределение в гидравлической сети, с коэффициентами трения, являющимися функцией числа Рейнольдса. Топология расчетной области представлена на рис. 5. Гидравлическая сеть состоит из 18 связей, номера которых показаны на рис. 5 в середине каждой связи, и 13 узлов, номера которых также приводятся на рис. 5. В узле 1 задан постоянный напор, схематически изображенный на рис. 5. В результате гидравлического расчета требуется определить напор в узлах сети. Полная постановка данной задачи приводится в [15], где также приводится об-
3
суждение решения задачи, полученное при помощи двух вариантов градиентного метода, и распределение напоров по узлам гидравлической сети.
Авторами использовались те же зависимости для определения коэффициента гидравлического сопротивления трения, что и в работе [15]. Данные зависимости [16] аналогичны формулам Блазиуса для ламинарного и Альтшуля для турбулентного режима течения, используемых в отечественной литературе [17].
Сравнение работы подпрограммы, вычисляющей коэффициент гидравлического сопротивления трения, с используемыми в работе [15] и данными, полученными авторами, представлено на рис. 6. Как видно из анализа рисунка, гидравлическое сопротивление трения рассчитывается корректно во всем диапазоне рассматриваемых чисел Рейнольдса. Сравнение результатов, полученных при помощи МКО, и результатов статьи [15] сведены в табл. 3.
Анализ сравнения результатов работы [15] и результатов разработанного метода (МКО) показывает, что максимальная относительная погрешность не превышает 0,03%. Однако следует отметить, что скорость сходимости МКО незначительно уступает скорости сходимости градиентного алгоритма Тодини, а общее число итераций МКО превышает количество итераций стандартного программного обеспечения ЕРАКЕТ [7] всего на 4 итерации.
£ 0,2
К
га 0,05
.и о
Si mpson
N jmerical
20000 30000
Число Рейнолъдса.
Рис. 6. Гидравлические сопротивления: Simpson — данные статьи [15]; Numerical — данные, используемые авторами
Сравнение результатов
Таблица 3
HFdM, м HRefi м 5, %
240,00000 240,000 0,00000
203,38794 203,403 0,00740
200,49390 200,509 0,00753
223,58027 223,588 0,00346
202,34512 202,360 0,00735
200,48419 200,499 0,00739
197,03134 197,048 0,00845
191,80189 191,818 0,00840
191,04145 191,058 0,00866
118,35004 118,384 0,02869
139,38592 139,416 0,02158
188,43973 188,456 0,00863
Заключение. Авторами был совершен переход от стандартной формы построения численного метода для решения задачи потоко-распределения, основанного на уравнениях Кирхгофа, к численному решению задачи потокораспределения для гидравлической сети с помощью построения дискретного аналога для уравнения неразрывности, которое приводится к дифференциальному уравнению второго порядка относительно давления.
Данный метод успешно протестирован на примере решения задач участка тепловой сети и сети водоснабжения как с постоянными коэффициентами сопротивлениями, так и с коэффициентами сопротивления, являющимися функцией от числа Рейнольдса. При этом скорость сходимости предлагаемой методики выше, чем у «увязочных» методов и близка к скорости сходимости градиентного алгоритма Тодини, уступая всего на несколько итераций.
ЛИТЕРАТУРА
[1] Cross H. Analysis of flow in networks of conduits or conductors. Engineering Experiment Station, University of Illinois, 1936, 38 p.
[2] Todini E., Pilati S. A gradient method for the analysis of pipe networks. International Conference on Computer Applications for Water Supply and Distribution. Leicester Polytechnic, UK, September 8—10, 1987.
[3] Меренков А.П., Хасилев М.Ю. Теория гидравлических цепей. Москва, 1985, 279 с.
[4] Андрияшев М.М. Техника расчета водопроводной сети. Москва, ОГИЗ, 1932, 62 с.
[5] Лобачев В.Г. Новый метод увязки колец при расчете водопроводных сетей. Санитарная техника, 1934, № 2, с. 8-12.
[6] Kirchhoff G. Ueber die Auflo- sung der Gleichungen, auf welche man bei Untersuchung der linearen Verthei- lung, galvanische Strome gefuhrt wird. Leipzig.
Annalen der Physik und Chemie (Poggendorf), 1847, Bd. 72, N 12, ss. 497-508.
[7] Lewis A. Rossman, EPANET 2, USERS MANUAL. Water Supply and Water Resources Division National Risk Management Research Laboratory Cincinnati, OH 45268, September, 2000, 200 p.
[8] Файзуллин Р.Т. Расчет и оптимизация больших гидравлических сетей. Труды Международной конференции RDAMM-2001, т. 6, ч. 2, Спец. выпуск. Новосибирск, с. 638-641.
[9] Норри Д., де Фриз Ж. Введение в метод конечных элементов. Москва, Мир, 1981, 299 с.
[10] Расчет параметров микроклимата с учетом конденсации влаги в рудничной вентиляционной сети. Аэрология, метан, безопасность: Сб. статей. Отдельный выпуск Горного информационно-аналитического бюллетеня. Москва, Изд-во «Горная книга», № ОВ7, 2011, с. 331-344.
[11] Применение совмещенной сетки для численного решения трехмерных задач гидродинамики и теплообмена методом контрольного объема. Инженерно-физический журнал, 1998, вып. 71, № 4, с. 740-745.
[12] Patankar S.V. Numerical heat transfer and fluid flow. Taylor and Francis, 1981, 196 p.
[13] Скибин А.П. Вариант конечно-элементного метода контрольного объема для решения задач тепломассообмена. Дис. ... канд. техн. наук. Москва, 1993, 222 с.
[14] Логинов К.В. Модели и алгоритмы расчетов режимов работы сложных гидравлических сетей. Дис. ... канд. техн. наук: Омск, 2004, 137 с.
[15] Jacobian Matrix for Solving Water Distribution System Equations with the Darcy-Weisbach Head-Loss Model. J. Of Hydraulic Engineering, 2011, vol. 137 (6), pp. 696-700.
[16] Swamee P., Jain A. Explicit equations for pipe flow problems. J. Hydraul. Div., 1976, vol. 102 (5), pp. 657-664.
[17] Идельчик И.Е. Справочник по гидравлическим сопротивлениям. Москва, Машиностроение, 1975, 559 с.
Статья поступила в редакцию 31.05.2013
Ссылку на эту статью просим оформлять следующим образом: Белова О.В., Волков В.Ю., Скибин А.П. Метод контрольного объема для расчета гидравлических сетей. Инженерный журнал: наука и инновации, 2013, вып. 5. URL: http://engjournal.ru/catalog/machin/vacuum/764.html
Белова Ольга Владимировна родилась в 1971 г., окончила МГТУ им. Н.Э. Баумана в 1995 г. Канд. техн. наук, доцент кафедры «Вакуумная и компрессорная техника» МГТУ им. Н.Э. Баумана. Автор более 20 научных работ в области компьютерного моделирования инженерных систем. e-mail: ovbelova@yandex.ru
Волков Василий Юрьевич родился в 1989 г., окончил МГТУ им. Н.Э. Баумана в 2012 г. Инженер ОАО ВНИИАЭС. Специализируется в области вычислительной газодинамики и тепломассообмена.
Скибин Александр Петрович родился в 1963 г., окончил МВТУ им. Н.Э. Баумана в 1986 г. и МГУ им. М.В. Ломоносова в 1988 г. Канд. техн. наук, начальник бюро ОАО ОКБ «ГИДРОПРЕСС». Автор более 80 научных работ в области вычислительной гидрогазодинамики и тепломассообмена. e-mail: askibin@yandex.ru