г. Москва и Московская область, Россия
Статья посвящена математическому моделированию распределения потоков в гидравлических сетях. Расчеты гидравлических сетей проводятся на стадии их проектирования и эксплуатации. Результаты численного моделирования используются при управлении работой гидравлической сети в режиме реального времени. Математическая модель распределения потоков в гидравлической сети представляет собой систему нелинейных уравнений. Метод узловых давлений, используемый для численного решения системы уравнений, представляет собой n-мерный метод Ньютона. Для обеспечения стабильной и быстрой сходимости итерационного процесса предлагается использовать начальное приближение, учитывающее топологию сети и параметры ее объектов, нижний релаксационный множитель, оптимизировать структуру матрицы Максвелла. Представленные в статье алгоритмы позволяют существен-но уменьшить размерность решаемой системы нелинейных уравнений.
гидравлическая сеть, метод узловых давлений, ленточная матрица, эквивалентирование
Введение
Основные положения теории гидравлических сетей были заложены В.Я. Хасилевым и А.П. Меренковым в шестидесятые годы прошлого века. Их идеи развиваются в работах их соратников и учеников. Исторически для решения задачи распределения потоков в гидравлических сетях предпочтение отдавалось методу контурных расходов (МКР). Даже для полностью разомкнутых сетей контуры создавались искусственно. Бесспорным преимуществом МКР перед методом узловых давлений (МУД) является значительно меньшая размерность системы уравнений. В МКР размерность равна количеству главных контуров, в МУД – количеству узлов с неизвестными давлениями. Преимуществом МУД перед МКР является то, что метод позволяет моделировать работу гидравлической сети любой конфигурации (с контурами, без контуров, с любым количеством узлов с заданным давлением) без дополнительного искусственного зацикливания сети.
Переход на «умные» технологии требует моделирования работы гидравлических сетей в реальном времени. Это означает, что пользователь должен иметь возможность оперативно изменять режим эксплуатации сетью (открывать/закрывать запорные устройства, переключать насосы, подключать/отключать дополнительные источники и т.д.) и быстро получать результаты моделирования. Эти задачи равноценны по значимости. Первая задача – разработка интерфейса. Пользователь готов затратить много времени и сил на сбор и первичный ввод информации, но при этом корректировать данные оперативно. Вторая задача – собственно моделирование работы гидравлической сети. Основное требование - быстродействие для моделирования сетей реальных размеров. В статье рассматривается только вторая задача.
Methods
При известной топологии и параметрах гидравлической сети ее математическая модель строится на основании законов Кирхгофа:
где A – матрица инцидентности; В – матрица главных контуров; S – диагональная матрица сопротивлений связей; |X| - диагональная матрица расходов на участках; 
- вектор расходов в узлах; 
- вектор действующих напоров на участках; 
- вектор расходов на участках; 
- вектор разностей давлений на участках; 
- вектор давлений в узлах.
После упрощения системы уравнений (1-4) узловая система уравнений приобретает вид:
Итерационный процесс строится с использованием формулы:
матрица Максвелла, вычисленная при текущем векторе расходов 
- невязка в узле.
Метод узловых давлений – это метод Ньютона, который, как известно, имеет высокую скорость сходимости только в окрестности точки решения. Выбор начального приближения для итерационного процесса очень важен. Рекомендации по решению этой проблемы не существуют.
Для поиска начального приближения к решению системы уравнений (5-6) предлагается заменить зависимость (6) на линейную зависимость:
Тогда начальное приближение 
может быть найдено в результате решения системы линейных уравнений:
где 
матрица Максвелла, вычисленная при 
Таким образом, при вычислении начального приближения используется вся имеющаяся информация о топологии сети и текущих параметрах объектов. Не смотря на это, в некоторых случаях скорость сходимости не достаточно высока. Для обеспечения стабильной монотонной сходимости в формулу (7) вводят релаксационный множитель 
, значение которого вычисляют на каждом шаге итерации. Такой подход существенно увеличивает время, затрачиваемое на одну итерацию. В теоретически обосновывается, что для положительно определенного Якобиана, для обеспечения монотонной сходимости достаточно выбирать <1. Практический опыт автора позволил установить следующий ряд значений =0,5 для k < 6, = 0,6 для k < 8, = 0,7 для k < 10, =1 для 
. Данный подход обеспечивает сходимость процесса за максимум 15 итераций.
На каждом шаге итерации необходимо решать систему линейных уравнений или обращать матрицу (8). По своей структуре матрица Максвелла (8) совпадает с матрицей смежности, которая для графов, моделирующих гидравлические сети, является разреженной. Для уменьшения времени, затрачиваемого на одну итерацию, узлы графа переномеровываются в соответствии с алгоритмом, что позволяет привести матрицу к ленточной форме с минимальной шириной ленты. Следующий шаг в оптимизации времени расчета – использование метода Холецкого для решения систем линейных уравнений с положительно определенными ленточными матрицами.
Размерность решаемых задач стала такой большой, что предложенных способов оптимизации времени расчета уже не достаточно. В работах предлагаются разные способы уменьшения размерности решаемой задачи, одни из них не вносят дополнительной погрешности в результаты, другие способы вносят минимальную погрешность. В работах излагается идея упрощения задачи, но алгоритмы ее решения не приводятся.
При разработке алгоритмов важна организация данных. В описанных ниже алгоритмах информация о топологии сети хранится в двумерном массиве, каждая строка которого соответствует участку сети, а столбцы содержат номера узлов, инцидентных участку.
В процессе эксплуатации гидравлической сети происходит перекрытие задвижек. В результате чего, некоторые части сети могут оказаться отключенными от источника. Выявление подобных ситуаций необходимо делать автоматически.
Высокоуровневое описание алгоритма выявления множества узлов, связанных с источником, представлено ниже.
Connectivity_checkg – выделение узлов сети, связанных с источником (выявление связных компонент графа сети).
- Маркировать узел «Источник».
- Кол-во_маркированных_узлов = 1.
- Повторять пока Кол-во_маркированных_ узлов >0:
- Кол-во_маркированных_ узлов = 0
- Для каждого участка:
- Найти участок, у которого маркирован один из инцидентных ему узлов;
- Маркировать смежный узел;
- Кол-во_маркированных_ узлов увеличить на 1.
В дальнейшем в расчетах рассматриваются только узлы, связанные с источником.
На втором этапе происходит поиск и удаление из рассмотрения ненагруженных деревьев. Под ненагруженным деревом понимается граф, не имеющий циклов, для узлов которого вектор G=0. Высокоуровневое описание алгоритма представлено ниже.
Unloaded_trees - поиск ненагруженных деревьев.
search for unloaded trees
- Вычислить валентность всех узлов графа.
- Кол-во_маркированных_связей = 1.
- Повторять пока Кол-во_маркированных_связей >0:
- Кол-во_маркированных_ связей = 0
- Для каждого немаркированного участка:
- Если валентность одного из инцидентных связи узлов равна 1, он не имеет нагрузки и не является источником, то:
- валентность узлов, инцидентных связи, уменьшить на 1;
- связь маркировать;
- Кол-во_маркированных_связей увеличить 1.
Дальнейшее упрощение расчетной схемы происходит за счет вычленения и удаления нагруженных деревьев из графа сети. Процесс начинается с терминальных узлов. В результате определяются расходы на всех участках и нагрузка корневого узла, которая увеличивается на суммарную нагрузку всех узлов, вошедших в дерево (рисунок 1). Индексы всех узлов и связей, вошедших в нагруженные деревья, равны 1. После того, как найдено давление в корневом узле, находятся давления во всех узлах, принадлежащих дереву (модуль Unfurl). Высокоуровневое описание алгоритмов представлено ниже.
Рисунок 1 Последовательность эквивалентирования нагруженных деревьев гидравлической сети a) – исходное состояние сети; b) перенос нагрузки в узлы 3 и 4; отсекаемые узлы - 1 и 2; c) перенос нагрузки в узел 4; отсекаемый узел 3; d) перенос нагрузки в узел S; отсекаемый узел 4.
Furl – скручивание нагруженных деревьев.
- Вычислить валентность всех узлов графа.
- Кол-во_маркированных_связей = 1.
- Повторять пока Кол-во_маркированных_связей >0:
- Кол-во_маркированных_связей = 0
- Для каждой немаркированной связи j (i;k):
- Если валентность одного из инцидентных связи узлов равна 1, и он не является источником, то:
- Валентность обоих узлов уменьшить на 1;
- Связь маркировать;
- Присвоить индексу узла I значение 1;
- Кол-во_маркированных_связей увеличить на 1;
- Если отсекаемый узел – конец связи (k), то:
- присвоить значению расхода значение нагрузки отсекаемого узла
- увеличить значение нагрузки смежного узла на значение нагрузки отсекаемого узла
- Если отсекаемый узел – начало связи (i), то:
- присвоить значению расхода инверсное значение нагрузки отсекаемого узла
- уменьшить значение нагрузки смежного узла на значение нагрузки отсекаемого узла
Unfurl – раскручивание свернутых нагруженных деревьев.
Кол-во_маркированных_связей = 1.
- Повторять пока Кол-во_маркированных_связей >0:
- Кол-во_маркированных_связей = 0.
- Для каждой немаркированной связи j (i;k) , принадлежащей дереву:
- Если индекс одного из инцидентных связи узлов равен 1 и он не маркирован, то:
- Вычислить давление в смежном узле по формуле:
- Маркировать смежный узел;
- Присвоить значение 1 индексу связи;
- Связь маркировать;
- Кол-во_маркированных_связей увеличить на 1.
Эффективность применения модуля Furl зависит от типа сети и решаемой задачи. Например, наладочный расчет систем теплоснабжения производится при заданной тепловой нагрузке абонентов, которая пересчитывается в расчетные расходы теплоносителя. В этом случае большую часть сети, а в некоторых случаях и всю сеть, удается «свернуть».
Графы, моделирующие расчетные схемы гидравлических сетей, могут содержать большое количество последовательных и параллельных связей. Эквивалентирование таких связей позволяет существенно уменьшить размерность решаемой системы уравнений.
Последовательные ребра – ребра, инцидентные одному и тому же узлу, который не имеет нагрузки, степень которого равна 2. При замене двух последовательных связей одной (рисунок 2 a-b) на единицу уменьшается количество узлов и добавляется новая связь. Гидравлическое сопротивление новой связи равно сумме сопротивлений замененных связей.
Параллельными связями называются связи, у которых совпадают пары инцидентных им узлов. Порядок перечисления узлов в паре не имеет значения. При замене двух параллельных связей одной не меняется количество узлов и добавляется новая связь, гидравлическая проводимость которой равно сумме проводимостей замененных связей, соединяющая те же концевые узлы (рисунок 2 c-d).
Процесс упрощения схемы начинается с поиска и эквивалентирования последовательных связей, затем происходит поиск и эквивалентирование параллельных связей, затем снова поиск последовательных связей, затем параллельных и так до тех пор, пока не останется последовательных и параллельных связей. Связи, инцидентные узлам, моделирующим насосы или регуляторы не эквивалентируются.
Рисунок 2 Эквивалентирование связей сети. a-b) последовательное соединение; c-d) параллельное соединение; a, c ) - исходное состояние; b-d) результат эквивалентирования.
Обратный ход состоит в присвоении значений расходов дочерних связей значениям расходов родительских связей (последовательные связи) и в вычислении значений расходов родительских связей (параллельные связи). Высокоуровневое описание алгоритмов представлено ниже.
Equivalenting_forward
- Вычислить валентность всех узлов графа.
- Определить количество узлов без нагрузки, валентность которых равна 2 (number_serial).
- Определить количество параллельных связей (number_parallel).
- Количество связей new=nl
- Повторять пока есть последовательные и параллельные связи (number_serial+ number_parallel>0):
- Если есть последовательные связи (number_serial>0), то
- Для каждой немаркированного узла i, валентность которого равна 2, и у которого отсутствует нагрузка:
- Найти 2 немаркированные связи j1 (i; k1) и j2 (i; k2), инцидентные узлу i;
- Маркировать найденные связи j1 и j2;
- Маркировать узел i;
- Создать новую связь new=new+1 (k1;k2);
- Вычислить гидравлическое сопротивление связи new по формуле Snew =Sj1+Sj2;
- Вычислить гидравлическую проводимость связи new Rnew=Snew-2;
- Сохранить номер связи new как дочерний subsidiary affiliate child
для связей j1 и j2.
- Определить количество параллельных связей (number_parallel).
- Если есть параллельные связи (number_parallel>0), то
- Для каждой немаркированной связи j (i; k):
- Найти параллельную ей немаркированную связь j1 (i; k);
- Маркировать найденные связи j и j1;
- Создать новый участок new=new+1 (i; k);
- Вычислить гидравлическую проводимость связи new по формуле Rnew =Rj+Rj1;
- Вычислить гидравлическое сопротивление связи new по формуле Snew = Rnew-0.5;
- Сохранить отрицательный номер связи new как дочерний для связей j и j1.
- Определить количество узлов без нагрузки, валентность которых равна 2 (number_serial).
Equivalenting_back
- Для каждого маркированной связи j, включая новые, начиная с большего номера с шагом -1:
- Найти узлы (i;k), инцидентные связи j;
- Вычислить
; - Для каждой немаркированной связи j1, начиная с 1 до j:
- Найти узлы (i1;k1), инцидентные связи j1;
- Найти родительскую связь (по совпадению номера дочернего участка j1 c номером j);
- Если номер дочерней связи >0 (т.е. последовательное соединение), то вычислить расходы на родительских связях:
, и
.
- Если номер дочерней связи <0 (т.е. параллельное соединение), то
, и
.
- Вычислить
- Маркировать связь.
Выводы
- Анализ состояния гидравлических сетей требует проведения расчетов в реальном времени. Этот факт повышает требования к надежности и быстродействию вычислительных процессов.
- Эквивалентные преобразования сети позволяют существенно уменьшить размерность систем нелинейных уравнений, а в отдельных случаях ликвидировать итерационный процесс, заменив его на последовательные вычисления.
- Использование в методе узловых давлений нижнего релаксационного множителя обеспечивает быструю и монотонную сходимость итерационного процесса.
Приведение матрицы Максвелла к ленточной форме и оптимизация ширины позволяют сократить время проведения одной итерации.
1. Меренков А.П., Хасилев В.Я. Теория гидравлических цепей М.: Наука, 1985. 279 с.
2. Новицкий Н.Н., Шалагинова З.И., Михайловский Е.А. Объектно-ориентированные модели тепловых пунктов теплоснабжающих систем Вестник ИрГТУ 21(9) 2017 с.157-172
3. Шалагинова З.И., Новицкий Н.Н., Токарев В.В., Гребнева О.А. Многоуровневое моделирование теплогидравлических режимов больших систем теплоснабжения Энергетика России в XXI веке. Инновационное развитие и управление 2015 Иркутск с. 389-398
4. Китайцева Е.Х. Обобщенные методы расчета воздушного режима зданий и факторов, влияющих на качество внутреннего воздуха. Автореферат кандидата технических наук МГСУ Москва 1995 18 с.
5. Дж. Ортега, В. Рейнболдт Итерационные методы решения нелинейных систем уравнений со многими неизвестными М.: Мир, 1975. 558 с.
6. Review of optimization models for the design of polygeneration systems in district heating and cooling networks Ortiga J, Bruno J C, Coronas A and Grossman I E 17th European symposium on Computer Aided Process Engineering. ESCAPE 17 Elsevier 2007
7. М. Свами, К. Тхуласираман Графы, сети и алгоритмы М.: Мир, 1984. 455 с.
8. Тьюарсон Р. Разреженные матрицы М.: Мир, 1977. 189 с.
9. Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке Алгол. Линейная алгебра М.: Машиностроение, 1976. 389 с.
10. Логинов К.В. Эквивалентирование гидравлических схем при моделировании крупных районов теплосетей Математические структуры и моделирование 2004 13 с.62-71
11. С.В. Якшин Применение метода расщепления графа при оптимизации параметров тепловой сети Вестник ИрГТУ 2018 22(10) с. 129-140
12. Щербаков В.И., Нгуен Х. К. Принцип энергетического эквивалентирования для расчета сетей водоснабжения с множеством участков. 2016 Вестник РУДН 2016 4 с. 27-34
13. С.В. Якшин Метод расщепления графа и принцип аддитивности тепловой сети Вестник ИрГТУ 017 21(4) с. 127-138
14. Сидорова В.Т., Карчин В.В. Перераспределение потока мощностей в сложно-замкнутых воздушных сетях 10 кВ для уменьшения потерь и улучшения качества электроэнергии 2016 Проблемы энергетики 2016 11 с. 51-54
15. Todini E, PIlati S A gradient algorithm for the analysis of pipe network Computer Applications in Water Supply London John Wiley & Sons 1 с. 1-20
16. Wang Jinda, ZhouZhigang and Zhau Jianing A method for the steady-state thermal simulation of district heating systems and model parameters calibration Energy Conversion and management 2016 120 с. 294-305
17. Static study of traditional and ring networks and use of mass flow control in district heating applications Mauna Kuosa, Kaisa Kontu, Tapio Makila, Markku Lampinen and Risto Lahdelma Applied Thermal Engineering 2013 54(2) с. 450-459
18. Lipovka A Y, Lipovka Y L Application of “Gradient” algorithm to modeling Thermal pipeline networks with pumping stations J Siberian Federal University/ Engineering and technologies 2013 6 с. 28-35
