Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 3 (2), с. 168-172
УДК 681.324
ЧИСЛЕННЫЕ МЕТОДЫ ТЕПЛОГИДРАВЛИЧЕСКОГО РАСЧЕТА СЕТЕЙ ИНЖЕНЕРНЫХ КОММУНИКАЦИЙ
© 2011 г. Р.А. Штыков
Муромский институт Владимирского государственного университета
Поступила в редакцию 20.01.2011
Рассматривается новая методика проведения гидравлических расчетов в автоматизированных системах управления инженерных коммуникаций большой размерности.
Ключевые слова: инженерные сети, теплоснабжение, матрица, математическая модель, алгоритм расчета.
Г идравлические расчеты лежат в основе анализа режимов тепловых, газовых, водопроводных и напорных канализационных сетей. Любые информационные системы по инженерным сетям, не предусматривающие проведения гидравлических расчетов, имеют крайне ограниченные возможности применения и потому вряд ли могут рассматриваться всерьез.
Методов решения задач гидравлического расчета вполне счетное количество, и они также хорошо известны; но проблема состоит в его более или менее приличном изготовлении. Поэтому на первый план выступает качество и алгоритмов, и программной реализации гидравлического расчета.
Рассмотрим схему с установившимся движением жидкости, состоящую из и участков, 7 узлов и к линейно независимых контуров, к = и - г +1.
Разработку линейно независимых циклов для выполнения расчета потокораспределения в гидравлических сетях можно осуществить при помощи основного дерева графа гидравлической системы для определения матриц соединений и матриц контуров. Построение основного дерева выполняется с помощью алгоритма Краскала. На схеме выделяется некоторое дерево, связывающее все её 7 узлы. В результате все участки разобьются на (г -1) участок дерева и к участков, не вошедших в это дерево, которые называются хордами. Каждая хорда замыкает определенную последовательность участков дерева и однозначно определяет контур, который фиксируется соответствующей строкой матрицы контуров.
Определим эти две основополагающие матрицы.
Введем полную матрицу А соединений 7 узлов и и ветвей сети, однозначно описывающую ее конфигурацию, безотносительно к конкретным длинам ветвей и фактическому месторасположению узлов. В этой матрице на пересечении строки у, соответствующей узлу у и столбца
I, соответствующего ветви I, помещается элемент а-. Этот элемент принимает следующие
значения: а- = 0 - если ветвь I не соединена с
—
узлом у, а- = -1 - если ветвь I исходит из узла у, а- = 1 - если узел у является для ветви I конечным. Таким образом, в каждом столбце матрицы А только два не нулевых элемента, это 1 и -1, поэтому сумма всех ее строк дает нулевую строку, что означает их линейную зависимость. Поэтому ранг матрицы, то есть максимальное число ее линейно независимых строк или столбцов составит 7 - 1. Для расчетов будем использовать усеченную матрицу соединений с размерами (г -1) • и , имеющую только линейно
независимые строки, ее обозначим А . Она получается из матрицы А вычеркиванием любой из строк, как правило этой строкой является строка 7.
Введем полную матрицу В контуров. В матрице контуров на пересечении строки г, соответствующей контуру к и столбца I, стоит элемент Ъ .. Этот элемент принимает следующие
значения: Ь = 0, если ветвь I не принадлежит
контуру г, Ъ„. = 1, если ветвь I входит в контур к
и ее ориентация совпадает с направлением его обхода, Ь . = -1, если ветвь I входит в контур г,
но ее ориентация противоположна направлению его обхода.
Для гидравлических сетей имеет место квадратичный закон гидравлического сопротивления:
^ + Ні = Г х2 >
(1)
где х. - неизвестный расход, гг - сопротивление участка, - напор равный разности пьезометрических отметок на концах участка (линейное падение напора), Н. - активный (действующий)
напор.
Зная длину трубопровода и его внутренний диаметр, данное сопротивление найдем по формуле:
А Оча) + 4)
гі (к) - ■
і(к) + 4, (к) , (к)
(2)
где /э - эквивалентная длина местных участков
(сопротивлений), I - длина участка сети, 1 -внутренний диаметр трубопровода, g - ускорение свободного падения, к - коэффициент, характеризующий хорду, Л8 - постоянный коэффициент, зависящий от фэ,
(3)
где ф - абсолютная эквивалентная шероховатость,
к =-
к
(4)
сать в виде
С использованием матрицы соединений А, балансы могут быть записаны сразу для всей схемы: Ах = V.
Слева в формуле стоит алгебраическая сумма расходов по всем участкам, имеющим общий узел у, а справа расход в узле: нагрузка у_ > 0,
если в узле у находится потребитель, приток V- < 0, если в узле у находится источник и
V = 0, если узел является точкой разветвления
потоков по схеме.
Значение V. < 0, должно быть задано таким
образом, чтобы имел место их общий нулевой баланс по всем 7 узлам схемы:
2 г-1
X V = 0, Кг =-£г7. . (8)
7=1 М
Второй закон Кирхгофа требует суммарного нулевого изменения напоров й,- для любого
контура схемы - для этого достаточно, чтобы выполнялось равенство
X *, = 0. (9)
с
С использованием матрицы контуров В, напоры могут быть записаны сразу для всей схемы: Бк = 0.
С учетом (6) и (9) уравнения, соответствующие второму закону Кирхгофа, имеют вид
где I * - сумма коэффициентов местных сопротивлений, dм - внутренний диаметр местных участков, X - коэффициент гидравлического трения.
В случае, когда нет данных о местных сопротивлениях, эквивалентная длина принимается приближенно в долях от линейной длины:
1э = а I. (5)
Поскольку по знаку величины гг необходимо судить о направлении потока на участке и соответственно о знаке й., то (1) можно запи-
X я,- Х| Х| = нк,
(10)
где Я., Хі - гидравлическое сопротивление и расходы на всех участках і контура к, Ик - алгебраическая сумма действующих напоров на всех участках, входящих в контур к .
Таким образом, математическая модель по-токораспределения для гидравлических сетей сводится к системе уравнений, состоящей 7 - 1 линейных уравнений вида (7) и к нелинейных уравнений вида (11):
Xх = V> ■/ = 1>->т-1>
X Ъ х,\х\ = нк.
(11)
(6)
Для любого потокораспределения должны выполнятся два сетевых закона Кирхгофа. Во-первых, в каждом узле выполняется материальный баланс:
(7)
Решение задачи потокораспределения сводится к решению системы нелинейных алгебраических уравнений вида (11). Как показала практика расчетов больших гидросистем (крупных городов или мегаполюсов), для систем уравнений с матрицами, размерность которых превосходит 500-700, начальное приближение
к
к
должно лежать настолько близко к решению, что практическая ценность вышеуказанных алгоритмов теряется - технические требования к точности решения на порядки слабее, чем требования увязочных методов к начальному приближению.
В качестве расчетного метода мы используем предложенный в работе [4] метод последовательных приближений, который обеспечивает сходимость на системах с числом неизвестных до 10000 при небольшом числе итераций.
Рассмотрим метод последовательных приближений применительно к системе (11).
Значение V- должно быть задано таким образом, чтобы имел место их общий нулевой баланс по всем 7 узлам схемы.
Будем решать полученную систему уравнений методом последовательных приближений. На каждой итерации решается система линейных уравнений вида:
Xі = sign(X)Я^\X,.|.
а)
в*
\п \і
да
У./1
' Ц " і .и :1 ,1.1111 ;у1' 11 .и 11' її ні
б)
і а ним №№0
(12)
|С
1
0Л
0Д1
(■.«і
о
V * .
*
1 і"),"1 і'і'Цііч. мц
*5
Нпри
где а; - коэффициент линейной части системы (а = 1 если энергоресурс в трубопроводе входит в узел, д. = -1 если энергоресурс трубопроводе выходит из узла и д. = 0 если участок I инцидентен узлу); й. - коэффициенты соответствующие нелинейной части системы, получаемые следующим образом:
В работах Файзулина, Денисова и Меренко-ва (см., например [4]) показано, что в качестве начального приближения целесообразно принять X0 = 1, (г = 1,2,..., п), а коэффициенты
а = 0.3, в = 0.7 , итерации продолжаются до тех пор, пока величина невязки между левой и правой частью уравнений нелинейной части не будет меньше заданной величины.
С целью обеспечения более быстрого затухания колебаний невязки и уменьшения числа итераций была проверена и предложена модификация, в которой после первой итерации применяется поправочный коэффициент ю = 0.48
Рис. 1. Расчет последовательных приближений:
а) стандартным методом, б) модифицированным методом
Докажем, что матрица соединений А и выбранное на схеме цепи дерево однозначно определяет матрицу контуров В, что может быть использовано для ее автоматического построения при расчетах на компьютере.
Перенумеруем переменные х; таким образом, чтобы расходы оказались в конце матрицы, то есть хг = (Хв, хк), где Х0 = (х1..х2_1) - расходы на участках дерева, хк = (..хи) - расходы на хордах.
Матрицу В можно записать в следующем виде:
В =
Ь11 К-1
(14)
(13)
То есть В = (Вб, Ек), А = ( Аб , А к), где Е - единичная матрица.
Таким образом,
Это позволило приблизить порядок получения вектора X к решению, сократив число итераций, колебания невязки быстро затухают. В качестве примера рассмотрим гидравлическую сеть, состоящую из 1000 ребер графа и 960 узлов (см. рис. 1).
Аэ хэ + Ак хк = V, следовательно,
х г = А г) (V — А к хк).
(15)
(16)
Матрица А в , обратная к , здесь всегда существует, поскольку в случае графа-дерева
к
1
1
его матрица соединений (без последней строки) будет иметь определитель, отличный от нуля, матрица не особенная.
Выявим связь между матрицами соединений А и контуров В. Рассмотрим однородную систему уравнений первого закона Кирхгофа Ах = 0, и покажем, что любая строка матрицы В является ее решением. Действительно, пусть а] = (а]Л а]п) - у'-я строка А(у = 1г-1),
= (Ьп агп) - г-я строка В(г = 1к), где
элементы а^ и Ьгг. (/ = 1,..., и) принимают значение 0, 1 или -1.
Когда контур г не проходит через узел у ненулевые элементы Оу. и Ъп- имеют обязательно различные номера I и, поэтому, скалярное произведение аТЬг = 0. Если же, простой контур г проходит через узел у, то ему могут принадлежать лишь две ветви 1Х, 12, инцидентные данному узлу, и только для них одновременно
не равны нулю соответствующие элементы а ■■
№
и Ън , так что скалярное произведение фактически будет сводится к сумме двух слагаемых:
ЛБ1 = 0.
(19)
а ■ Ь = а Ь + а Ь .
] Г ]1\ Г1\ ]12 Г1г .
Таким образом, если рассматривать матрицы блочной структуры то:
АВ1 =
= \Ав А
:]в Ек I
[а в А к ][в В к I
= А ВВВ + А кЕТк = 0, (20)
то есть
Ав Бп + А к = 0,
откуда следует, что
Бв = — А в А к.
(21)
(22)
(17)
Как видно из рис. 2, в любом случае прохождения контура через узел у, произведения в правой части равенства (17) обязательно имеют разные знаки, так что их сумма равна 0 и поэтому всегда аТ Ьг = 0 :
Рис. 2. Варианты прохождения контура г через узел j
Так как это справедливо для любых у и г, то:
ЛЬГ = 0, (г = 1,..., к), (18)
то есть любая строка матрицы В (берется здесь как вектор-столбец) удовлетворяет системе (14). Это может быть записано сразу для всех строк, что приведет к искомому произведению матриц:
Это выражение связывает матрицы Вв, А в, А к и отражает тот факт, что матрица соединений А и выбранное на схеме цепи дерево однозначно определяет матрицу контуров В.
Указанные действия позволили определить следующие данные:
1) методику гидравлического и теплового расчета инженерных коммуникаций, позволяющую провести расчет дополнительных параметров (температура, скорость) в рамках гидравлического расчета, что не требует проведения теплового и конструкторско-гидравлического расчета и экономит время и ресурсы;
2) при расчете стационарного потокорас-пределения в многокольцевых гидравлических сетях для обеспечения более быстрого затухания колебаний невязки и уменьшения числа итераций учесть поправочный коэффициент ю = 0.48;
т ——1 —
3) выражение Во = -Ао Ак, связывающее
матрицы БП, А б , А к и отражающее тот факт,
что матрица соединений А и выбранное на схеме цепи дерево однозначно определяет матрицу контуров В.
Список литературы
1. Андрияшев М.М. Техника расчета водопроводных сетей. М.: Сов. законодательство, 1932. 64 с.
2. Андрианов Д.Е., Макаров К.В., Штыков Р.А. Системы оперативного управления пространственно распределенными объектами. М.: Радио и связь, 2005. 211 с.
3. Бакластов А.М., Бродянский В.М., Голубев Б.Н. Промышленная теплоэнергетика и теплотехника: Справочник. М.: Энергоатомиздат, 1983. 551 с.
4. Денисов Е.Е. Математическое моделирование динамики жидкости с использованием теории графов // Математическое моделирование. 1996. № 2. С. 91-105.
NUMERICAL METHODS OF THERMOHYDRAULIC CALCULATION OF UTILITY NETWORKS
R.A. Shtykov
A new method of hydraulic calculations for automated control systems of large-scale utility networks is considered.
Keywords: utility networks, heat supply, matrix, mathematical model, calculation algorithm.