CC BY
42
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Лобачевского, 2010, № 2 (1), с. 140-150
УДК 532.54+531.011+519.853
ПРИМЕНЕНИЕ ПРИКЛАДНОЙ АНАЛИТИЧЕСКОЙ ГИДРОМЕХАНИКИ И МЕТОДОВ ПРИНЯТИЯ ОПТИМАЛЬНЫХ РЕШЕНИЙ В ЗАДАЧЕ НАХОЖДЕНИЯ ПОТОКОРАСПРЕДЕЛЕНИЯ В ГИДРОСИСТЕМАХ
© 2010 г. Л.В. Смирнов 1, В.А. Гришагин 1, Д.Н. Добряев 1, Н.В. Данилова 1,2
1 Нижегородский госуниверситет им. Н.И. Лобачевского 2 Нижегородский филиал Института машиноведения им. А.А. Благонравова РАН
Поступила в редакцию 15.07.2009
Показано, что задачу нахождения стационарного потокораспределения в произвольной гидросистеме можно свести к конечномерной многоэкстремальной задаче безусловной оптимизации. На примере четырехпетлевой системы циркуляции теплоносителя ядерной энергетической установки показана работоспособность предложенного алгоритма.
Ключевые слова: гидросистема, стационарное потокораспределение, аналитическая гидромеханика.
Задача нахождения потокораспределения в произвольной гидросистеме
Системы различных аппаратов и устройств, соединенные между собой трубопроводами, называются гидросистемами или гидравлическими сетями (ГС). Это системы тепло- и водоснабжения, системы перекачки рабочей среды в различных производствах, системы циркуляции теплоносителя ядерных реакторов. При решении технических и экономических проблем, возникающих при проектировании и эксплуатации ГС, встают два вида задач: задача нахождения равновесных режимов (задача стационарного потокораспределения), состоящая в определении стационарных значений расходов жидкости в различных участках ГС; и более общая задача исследования динамических свойств ГС, включая устойчивость состояний равновесия, структуру фазового пространства и их зависимости от параметров. Эту задачу можно назвать просто задачей потокораспределения. Геометрическая структура ГС представляется в виде ориентированного графа, ребрами которого являются участки ГС, а вершинами (или узлами) -элементы соединения и разделения потоков различных участков. Узлами могут быть тройники, коллекторы и смесительные камеры, в которых происходит соединение или разделение потоков.
В основе традиционной методики математического моделирования потокораспределения в ГС лежит используемый в трубопроводной гидравлике подход, состоящий в следующем [1]. Течение жидкости на участках предполагается
одномерным, осредненным по поперечному сечению потока и турбулентным пульсациям, а потери на трение определяются с помощью выражений, основанных на обработке и обобщении большого количества экспериментальных данных и эмпирических упрощенных соображений (формулы Дарси - Вейсбаха). При этом не рассматриваются быстрые волновые процессы, связанные со сжимаемостью жидкости, и течение со свободной поверхностью. Основным уравнением нестационарного движения вязкой жидкости при этом является уравнение Бернулли, а в качестве гидродинамических характеристик узлов используются давления и уравнения неразрывности в узлах.
При таком подходе математическая модель ГС из N участков, N из которых содержат насосы с угловыми скоростями вращения рабочих колес ъъ1, имеет вид:
= р1 - р2 + Р&(^1 - Ъ2) -АР, (бг,ъ1 X
і = 1, N
Хійі = Рі1 - Рі2 + р§(zi1 - zi2) -^Рі (0-і)
і = N +1, N;
(1)
мк
£а }б} = 0, к = 1, Ь,
У=1
Jlъг = м;к,ег) -мш,шг), / = 1,N1,
где Ql - расход жидкости на 1-м участке; х',2 СХ
т1 =р I -------; х1 Ха, х12 - координаты попе-
Sl (х)
Х1 14 '
речного сечения потока на участке, на его входе
и выходе соответственно; р - плотность жидкости; Б1(х1) - площадь сечения прохода; g - ускорение силы тяжести; Р11, г11 - соответственно давление и высота центра сечения прохода на входе (I = 1) и выходе (I = 2) участка;
АР, (Ql, ъ) = AP/(Ql) + АР'Щ., ъ) - РQl х
; i = 1, N; ДPi' (Qi) - по-
Sf(xa) Sf( Xi 2)
уравнения вращения твердого тела относительно неподвижной оси.
Для гидродинамических процессов воспользуемся уравнениями (1) при постоянных или меняющихся квазистатически медленных механических переменных ъ,:
TiQi = (P-1 -P2) + pg(zi1 - zi2) -ДР (Qi,®i )
i = 1, N;
(2)
Mk
тери на трение; APf”(Q,, ю,) - характеристика
насоса на участке, если он есть; M'(&i, 0,) -движущий момент, в общем случае зависящий от внешних воздействий 0, на привод насоса; M "(Qi, ю,) - момент сопротивления, определяющийся трением и расходом перекачиваемой среды; J - момент инерции рабочего колеса насоса и вращающихся элементов, связанных с ним механически. Последнее слагаемое выражения для AP, (Q,, ю,) представляет собой разность удельных кинетических энергий потока на входе и выходе участка.
Если имеются внешние воздействия на привод, то необходимо задать закон изменения переменных 0, или соответствующие уравнения (например, это могут быть уравнения электрических процессов в сетях питания электродвигателей насосов). Параметр а, - постоянный множитель, равный +1, если жидкость поступает в узел, и -1, если жидкость выходит из узла. Знаки этих множителей можно определить, зная геометрическую структуру и назначение ГС. При определении коэффициентов т, интегрирование выполняется по направлению течения от входа до выхода участка. Точкой здесь и ниже обозначается дифференцирование по времени.
Первые N уравнений системы (1) - это уравнения Бернулли [1] для участков, содержащих и не содержащих насосы. Следующие L уравнений - это уравнения неразрывности для L узлов, в каждом из которых соединяются Mk участков. Последние N уравнений - это уравнения угловых скоростей вращения рабочих колес насосов. Если на некоторых участках ГС имеются запорно-регулирующие устройства, то их гидравлические характеристики могут быть учтены, так же как и для насосов.
Система (1) содержит уравнения двух видов взаимодействующих друг с другом процессов -гидродинамических и механических. Первые N + L уравнений описывают гидродинамическую подсистему, а последние N уравнений представляют собой известные из механики
^а jQj = 0, k = 1, L.
J=1
При решении задачи стационарного потоко-распределения предположение о постоянстве переменных ъ, является естественным. Что же касается задач динамики, то это либо соответствует предположению об идеальности регулятора оборотов и типично для традиционных подходов к решению задачи, либо основывается на разделении движений, когда процессы изменения оборотов являются медленными. Последний случай типичен для систем циркуляции ядерных энергетических установок (СЦ ЯЭУ). Более подробно эти вопросы обсуждаются в литературе (см. например [2], где имеются указания на известные публикации).
Функции АР,, входящие в (1) и (2), суть суммарные гидравлические характеристики соответствующих участков системы. Аналитические выражения этих функций для участков с насосами и участков без насоса, согласно данным прикладной гидромеханики и теории гидравлических машин [1, 3], имеют вид:
АТЗ(п ) {аи01 -Ьг®гй -Сг®2, ^ - 0;
АР (й, Ъ Н 2 2 (3)
I- а2& - Ь1 ®гй - С1 ®г > й < 0;
iQi2, Qt > 0; i2
(4)
ЛР(&)=' 1г /
і- ^2ійі , йі < 0; где а1і, а2і, Ьі, сі, й\і, ^2і - положительные коэффициенты, вид которых зависит от гидравлических сопротивлений участков и гидравлических характеристик насосов. Качественный вид графиков характеристик участков без насоса и с насосом представлен соответственно на рис. 1 (а и б).
ГС может быть замкнутой или разомкнутой, как, например, система городского водоснабжения. В последнем случае на границах вводятся внешние узлы, в которых задаются предполагающиеся постоянными расходы или давления. Этот случай является частным.
Исследование системы нелинейных уравнений (2), которая, как правило, имеет высокий порядок, обычно проводится численно. В боль-
1
1
X
шинстве публикаций, посвященных исследованию ГС, рассматривается только задача стационарного потокораспределения и решаются численно уравнении статики, следующие из (2) при
(0, = 0, I = 1, N. Решение такой задачи порождает проблему выбора метода расчета, начального приближения и сходимости итерационного процесса, а проблемы неоднозначности состояний равновесия, устойчивости системы и существования других режимов работы системы [4, 5] даже не обсуждаются [6].
Предлагаемый в данной работе подход к исследованию системы уравнений (2) (в том числе и нахождению стационарных режимов) опирается на методы аналитической механики [7], теории нелинейных колебаний [8, 9] и теории принятия оптимальных решений [10, 11]. Система (2) может быть представлена в форме системы частного вида уравнений Лагранжа с интегрируемыми связями обобщенных скоростей, которыми в данном случае являются расходы
Яг, , = :
С дТ дR —
-------=-------, , = 1, N;
С дЯ, дЯ,
мк ___ (5)
Д. = 0, к = 1, Ь;
J=1
N Т О2
Т = 1^;
,=1 2
N
К = X [р2 - р1 + рg(2 - zil)]Qi +
,=1
N 0
+ Х|АР- (^ + С1,
,=1 0
Т - кинетическая энергия системы, учитывающая движение жидкости во всех N участках системы. Функция К, определенная с точностью до постоянного слагаемого Сь - образующая
функция, дифференцирование которой дает правые части дифференциальных уравнений системы (2).
Из этого представления следует, что обобщенные координаты являются скрытыми и в явном виде в уравнения не входят. Число обобщенных скоростей, входящих в уравнения (5), избыточно и может быть сокращено с помощью уравнений связи. Такими уравнениями являются алгебраические уравнения системы, среди которых могут быть уравнения-следствия.
После исключения избыточных переменных система (5) принимает вид [2]:
С дТ дК —
------=-------, , = 1, п;
Ж д0, О ’ ’
N т О2 N 0 (6)
Т = Е-0-; К = 1{АР(^+С1.
,=1 2 ,=1 0
Отличие системы (6) от системы (5) состоит в уменьшении числа дифференциальных уравнений и обращении в нуль первой суммы в выражении для функции К. При дифференцировании Т и К для подстановки в уравнения Лагранжа уже для независимых переменных учитываются следующие из уравнений связи выражения избыточных переменных через независимые. Размерность фазового пространства системы (1) равна п = N - L , где Ь - число независимых
уравнений связи (Ь < Ь).
Исследование динамических свойств ГС проведено с помощью прямого метода Ляпунова с использованием К(01,...,0п) в качестве функции Ляпунова [8, 9]. На основании определяющих вид зависимостей АР, (0,),, = 1, N , выражений (3), (4) можно утверждать, что эта функция ограничена снизу и неограниченно растет при удалении от начала координат фазового пространства 01,.,0п . Элементарные преобразования показывают, что
f= _^ 2’ dt i=1
(7)
то есть производная dR|dt отрицательна и равна нулю в состояниях равновесия.
Из выражений для К и dt на основании теорем прямого метода Ляпунова следует дис-сипативность системы и другие важные выводы, дающие полную информацию о структуре фазового пространства. В частности, диссипа-тивность системы означает, что все движения независимо от начальных условий заканчиваются в одном из устойчивых состояний равновесия, в которых К имеет минимум. Таким образом, все фазовые траектории идут из бесконечности в некоторую ограниченную область фазового пространства, в которой находится либо одно устойчивое состояние равновесия, и тогда К имеет единственный минимум, либо несколько состояний равновесия (устойчивые и неустойчивые), в этом случае К имеет несколько минимумов. Это позволяет обосновать принципиально новый метод нахождения состояний равновесия ГС, описываемой системой уравнений (2). Метод состоит в поиске минимумов
функции К(01,...,0п). В окрестности каждого из состояний равновесия 0°,...,0° = 00 имеем:
R = R
Q0,-, Q0)+ Ci
n n
(Qi - Q0) + LX
+
d 2 R
,=i j=i
dQdQj
x (8)
x(Q, -Q0)(Q, -Q0).
Первые два слагаемых в (8) выбором C1 можно обратить в нуль. Коэффициенты первой суммы равны нулю в каждом состоянии равновесия системы. Это условие соответствует обращению в нуль правых частей системы (2) с учетом содержащихся в этой системе уравнений неразрывности. Оно обычно и используется для поиска состояний равновесия без оценки устойчивости. Таким образом, функция R(QjQn)
в малой окрестности состояния равновесия аппроксимируется второй суммой ряда (8). Коэффициенты двойной суммы положительны в устойчивых состояниях равновесия, где R имеет минимум.
Таким образом, задачу нахождения стационарного потокораспределения можно свести к конечномерной задаче безусловной оптимизации [10, 11]:
R(Q) ^ min, Q е D,
(9)
где величины Ai, Bi, i = 1, n, - константы, задающие границы изменения координат вектора
Q = (Q1,k,Qn)T .
Как правило, данная задача многомерна и многоэкстремальна, поэтому целесообразно применять хорошо зарекомендовавшие себя схемы редукции размерности наряду с характеристическими методами поиска глобального экстремума.
Многошаговая схема редукции размерности
Рассмотрим многоэкстремальные задачи оптимизации, т.е. задачи, в которых целевая функция в области D может иметь несколько локальных экстремумов. На сложность решения данных задач существенное влияние оказывает размерность. К примеру, для класса многоэкстремальных функций, удовлетворяющих условию Липшица, имеет место экспоненциальный рост вычислительных затрат при увеличении размерности: если в одномерной задаче для достижения точности в требуется p вычислений функции,
то в задаче с размерностью N для решения с той же точностью необходимо осуществить количе-
N
ство испытании порядка p .
Многие подходы к конструированию численных методов анализа многомерных многоэкстремальных задач основаны на идее редукции сложности: решение исходной задачи заменяется решением системы более простых задач.
Одним из наиболее эффективных является подход, связанный с редукцией размерности, когда решение многомерной задачи сводится к решению одной или нескольких одномерных подзадач. В качестве примеров можно привести следующие схемы такого рода: многошаговая схема (схема вложенной оптимизации) [10] и редукция размерности на основе кривых, заполняющих пространство (кривых Пеано) [11].
Рассмотрим в качестве исходной задачу (9). В силу отсутствия функциональных ограничений исходное соотношение можно представить в виде [10]:
min R(Q) = min min ...
QeD Q1e[ A1,B1] Q2e[ A2,B2]
••• minn nR(QbQ2’ -’Qn). (10)
Qn e[ An A ]
Введем семейство функций:
R1 (Q1,...,Qi) = iinRi+1(Q1,...,Qi, Qi+1): Qi+1 е [A+1, Bf+1]}
1 < i < n -1.
= mini
0
0
Q
Q
Для решения исходной задачи достаточно решить одномерную задачу:
r1(Q1) ^ min Q1 е [ Ab B1] •
При этом каждое вычисление функции R1 (Q1) в некоторой фиксированной точке Q1 е [A1, B1 ] представляет собой решение одномерной задачи:
R 2(Q15 Q2) ^ min> Q2 е [A2’ B2] .
Эта задача является одномерной задачей минимизации по Q2 , так как Q1 фиксировано. Каждое вычисление значения функции R (Q1, Q2) при фиксированных Q1Q2 требует решения одномерной задачи, и т.д. вплоть до решения задачи
Rn (Qb-,Qn) -
= R(Q1,Q2, .,Qn) ^ min, Qn е [An,Bn] при фиксированных Q1,..., Qn _1.
То есть решение задачи (10) сводится к решению семейства вложенных одномерных подзадач вида:
Rl (Q1v, Qi ) ^ min, Qi е [Ai, Bi ], (11)
где Qlv.., Qi _1 - фиксированный вектор.
Решение исходной задачи (9) через решение системы одномерных рекурсивно связанных подзадач называется многошаговой схемой редукции размерности.
Характеристические алгоритмы глобального поиска
Для решения подзадач (11), полученных в результате использования многошаговой схемы редукции размерности, можно использовать одномерные алгоритмы глобального поиска. Рассмотрим класс численных методов поиска глобального экстремума, называемый классом характеристических алгоритмов [10, 12].
Рассмотрим следующую одномерную задачу: f (х) ^ min, х е [A, B] . (12)
Алгоритм решения задачи (12) называется характеристическим, если, начиная с некоторо-
/X -\Jk+1
го шага поиска, выбор координаты х очередного испытания заключается в выполнении следующих действий.
1) Задать набор
Лk = (х0,..., хт } конечного числа т + 1 = т(£) + 1 точек области G = [A, B], полагая, что A еЛ£, B еЛ£ , все
координаты предшествующих испытаний расположены по возрастанию:
А = х0 < х1 <... < хт-1 < хт = В .
2) Каждому интервалу (х,-1, х,), 1 << т, поставить в соответствие число К(,), называемое характеристикой интервала.
3) Определить интервал (х,-ьх,), которому соответствует максимальная характеристика:
К^) = тах{К(,') :1 < < т}.
4) Провести очередное испытание в точке
хк+1 = ) е (х -1, х).
Для проведения нового испытания отрезок [А, В] разбивается на т интервалов. Далее выбирается интервал, у которого характеристика наилучшая. Точка очередного испытания размещается внутри этого интервала в соответствии с правилом ).
В вычислительную схему характеристического алгоритма можно ввести условие оста-
новки вида
xt - xt-1 < S ,
где в > 0 - заданная точность поиска (по координате), то есть следует прекращать вычисления, когда длина интервала с максимальной характеристикой станет меньше заданной точности в. Тогда процесс поиска будет остановлен через конечное число шагов.
В качестве примера характеристического метода можно привести информационностатистический алгоритм глобального поиска (АГП) [10, 11]:
2
H(i) = m(xt - xt-1) +
(zi - zi-1)
m(x, - x,-1)
- 2(zi + г,-1)
„k+1
= 0.5(xt + xt-1) - (zt - zt-1)/(2m),
где z, - результаты испытании в точках x,. Величина m > 0 вычисляется:
[rM, M > 0
m =
1, M = 0
I z- - z- | где M = max —l--------------------l—^ , а r > 1
1<<T xi - xi-1
параметр
метода.
Индексная схема учета функциональных ограничений
В силу того, что задача (9), как правило, многоэкстремальна, а в результате поиска с использованием указанного выше подхода мы получаем только глобальный минимум, то возникает необходимость поиска остальных минимумов, чтобы найти все устойчивые состояния
равновесия. Одним из способов для нахождения всех минимумов является исключение уже найденных экстремумов путем ввода функциональных ограничений вида
(Q1 _Q*} + (Q2 _Q2} + ••• + Q _Q*} ^5,
/ * * * \ где Q1,Q2,. .,Qn) - точка найденного минимума, а 5 - радиус сферы, исключающей найденный минимум из дальнейшего рассмотрения.
После ввода вышеуказанных неравенств получаем многомерную задачу оптимизации с функциональными и координатными ограничениями, которая может быть представлена в следующем виде:
R(Q) ^ min, Q е G с Rn, (13)
G = {Q е D : g; (Q) < 0,1 < j < m},
D = {Q е Rn : Qi е [A;,Bt ], 1 < i < n}.
Таким образом, возникает необходимость учёта функциональных ограничений в процессе поиска. Многошаговую схему (10) можно обобщить на класс задач (13) [10] таким образом, что одномерные подзадачи (11) будут иметь следующий общий вид:
f (х) ^ min, х е [A, B ], gj(х) < 0,1 < j < m. (14)
Для решения этой задачи был предложен индексный алгоритм [11, 12], который может быть реализован в рамках характеристической вычислительной схемы и многошаговой схемы (в случае нескольких переменных). Отметим, что каждая итерация предлагаемого алгоритма включает определение индекса 1 < v(х1) < m + 1, равного номеру первого нарушенного ограничения. Если v(хi) < m +1, то точке хi соответствует значение zi = gv (х{). Если v(х{) = m + 1 (то есть все ограничения выполнены), то z{ = f (х{). Далее, согласно правилам R(i) и D(t), осуществляется процесс поиска в рамках характеристической схемы.
Индексный метод имеет ряд преимуществ перед классическим методом штрафных функций, так как не требует подбора каких-либо параметров (константы штрафа, нормирующих коэффициентов при ограничениях).
Система глобальной оптимизации PSoDI
В рамках подхода с использованием многошаговой схемы редукции размерности наряду с характеристическими методами и индексной схемой учёта ограничений была реализована
программная система PSoDI [13], позволяющая осуществить следующую последовательность действий, необходимую для решения многомерных задач многоэкстремальной оптимизации.
1) Постановка задачи и выбор методов поиска. Целевую функцию, координатные и функциональные ограничения можно задавать аналитически. Минимизируемая функция может быть либо выбрана из предложенного набора, либо введена пользователем. Также предусмотрена возможность задания оптимизируемых функций в виде подключаемых DLL-библиотек. Далее по каждой переменной необходимо сконструировать характеристический метод. Алгоритм поиска может быть как выбран из списка, предлагаемого программой, так и определен самостоятельно путем задания правил вычисления характеристики и точки следующего испытания аналитически. Способность конструирования алгоритмов дает широкие возможности как для эффективного решения задачи, так и для исследований.
2) Решение поставленной задачи. Процесс поиска оптимального значения осуществляется на основе многошаговой схемы редукции размерности с использованием сконструированных алгоритмов для решения одномерных подзадач. Для учета ограничений используется индексный метод.
3) Визуализация исходной задачи и процесса поиска. Средства программной среды позволяют строить целевую функцию (сечения функции, если количество переменных больше двух), линии уровня и отображать точки испытаний. Визуализация дает возможность наглядно представить минимизируемую функцию и процесс поиска, а также качественно оценить правильность работы вычислительной схемы.
Исследование стационарных режимов работы гидросистемы на примере системы циркуляции теплоносителя ядерной энергетической установки
Применение описанной выше методики продемонстрировано на примере решения задачи нахождения равновесных режимов системы циркуляции теплоносителя ЯЭУ, являющейся частным видом ГС. Нарушение нормальной, предусмотренной проектом, работы системы циркуляции теплоносителя недопустимо с точки зрения безопасности реактора. Результаты ее решения также имеют самостоятельное значение.
Рис. 2
На рис. 2 представлена расчетная модель типовой системы циркуляции. Теплоноситель, нагревающийся в активной зоне реактора 1, поступает по трубопроводам в параллельно работающие теплообменные петли 2 и после охлаждения в теплообменниках 3 снова подается в активную зону. Течение теплоносителя обеспечивается работой циркуляционных насосов 4, снабженных электродвигателями 5. Стрелками показано направление течения теплоносителя в нормальном эксплуатационном режиме. Эта схема при числе петель, равном четырем, соответствует установке типа ВВЭР-1000. Изменение плотности теплоносителя при нагревании и охлаждении невелико, и оценки, имеющиеся в литературе, показывают, что влиянием этих изменений в рассматриваемых ниже гидромеханических процессах можно пренебречь [5]. Это дает основание изучать гидромеханические процессы независимо от теплофизических.
Полная математическая модель динамики гидромеханических процессов, соответствующая представленной расчетной схеме для четырехпетлевой системы циркуляции, имеет вид:
Тгйг = (Р2 - Р1) + Pg (^2 - Zl) -АР, , а, ),
I = 1,4;
15^5 = (Р - Р2) + Р£(^ - Z2) - АР^); (15)
Jlа, = м;(а,, 0,) - МЩ-,а,), , = 1,4;
4
I в] =
]=1
где обозначения соответствуют обозначениям для системы (1). Переменные Ql, , = 1, 5 , и ю,, , = 1, 4 , имеют смысл обобщенных скоростей, а
обобщенные координаты, получающиеся интегрированием по времени, являются скрытыми и в уравнения движения не входят.
Математическая модель (15) с учетом зависимости плотности от температуры является частью более сложной математической модели, используемой для изучения не только гидромеханических, но и теплофизических процессов. При пренебрежении изменением плотности теплоносителя (о чем было упомянуто выше) эта система уравнений является автономной и пригодна для изучения гидромеханических процессов и общих динамических свойств системы циркуляции. Данная модель отражает динамику двух взаимодействующих процессов разной природы. Это, соответственно, процессы течения теплоносителя и вращения рабочих колес циркуляционных насосов и связанных с ними масс (редуктор, ротор электродвигателя и маховик). В ряде работ [5, 14] аналитически и численным счетом для параметров типового реактора ВВЭР-1000 показано, что гидродинамические процессы являются быстрыми и могут быть изучены в предположении, что механические переменные ю,, = 1, 4 , постоянны или
изменяются квазистатически по отношению к быстрым переменным Ql, , = 1, 5. Малым параметром, позволяющим разделить гидродинамические и механические процессы, является отношение кинетических энергий гидродинамических и механических процессов в стационарном режиме [2, 5]. Это отношение в соответствии с приведенными оценками является малым:
5
тт0
1Г
^ = ТГГ
ТМ
2
Ji (а
а I
<< 1.
2
Таким образом, оказывается возможным раздельно изучать быстрые и медленные процессы, то есть понизить порядок исследуемой системы уравнений (15) вдвое.
Для быстрых процессов после исключения Q5 с помощью последнего из уравнений системы (15) имеем систему:
4 / \ Г 4 ^
т5] + тiQi = -Ар , а° )-АР5 IQ
]=1
(16)
где значения медленных переменных ю , предполагаются постоянными либо медленно меняющимися в соответствии с решением системы уравнений медленных процессов, следующей из (15), при т, = 0, , = 1, 5.
т
Нелинейные функции АР, ^, ю ,), , = 1, 4, и АР5®5) определяются гидравлическими характеристиками участков (3), (4).
Функция Я приобретает вид:
5 Q¡
Я = ЦЩ (^ + С1, (17)
г=1 0
и ее аналитическое выражение получается интегрированием соответствующих зависимостей
АР, , , = 175, (3), (4).
Построение и анализ функции Я дают наглядное представление о структуре фазового пространства и об изменениях этой структуры в результате бифуркаций.
Для нахождения состояний равновесия системы (16) при традиционном подходе имеем:
Г 4 ^ _
, = 1,4. (18)
и=! )
Решение этой системы уравнений соответствует нахождению стационарного потокораспре-деления, то есть нахождению стационарных
значений расходов Qi0,, = 1, 4 .
Задача (18), в отличие от общего случая, имеет частный вид и решается с помощью несложного графического построения или соответствующего ему разработанного численного алгоритма [5].
Схематический способ решения системы (18) для общего случая, когда число петель равно не четырем, а п, состоит в следующем. Каждое уравнение системы (18) разрешим относительно Ql:
Г п У
— ЛР,(о,,ш»)=ЛР5 Юу
ЛР.
I в,
V ,=1
і = 1, п.
Сложим все п уравнений (19) и обозначим выражение в правой части через функцию Ф:
їв =1^
ЛР
( п \
[і О
V ,=1
= ф
ЛРп+11ІО,
\
Далее, найдем обратную для ф(і,=1в,)
функцию и вместо уравнений (19) получим од-
ІП=1 в,:
но уравнение относительно
(19)
ф-1 [і в, і=др„11:;*; в, і. (2°)
После решения уравнения (20) и подстановки результата в (18) получаем систему уже несвязанных уравнений для нахождения стационарных значений переменных в0, і = 1,п .
На рис. 3 представлены результаты такого построения для случая двух петель (п = 2) и
одинаковых зависимостей — ЛР, (в,,ю°), і = 1, 2 (кривая 4). Решение задачи стационарного по-токораспределения соответствует точкам пересечения характеристики общего участка АР3 О + 02) (кривые 1, 2, 3) и получающейся графически или численно суммарной характеристики системы петель Ф-1(01 + в2) (кривая 5). Последняя имеет двойную петлю, отмеченную жирной линией, и каждая точка на ней соответствует двум симметричным режимам, так как индексы , = 1, 2 входят симметрично.
При изменении нелинейных характеристик АР, (в,, ю0,) и характеристики общего участка, обусловленном изменением параметров, например гидравлического сопротивления реак-0
тора или величин ю ,, стационарные значения расходов в, изменяются, что соответствует изменению координат особых точек в пространстве состояний системы уравнений (18). Могут
г=1
,=1
происходить бифуркации в виде попарного рождения и исчезновения некоторых из этих точек с нулевым суммарным индексом Пуанкаре. Например, двухпетлевая система может иметь одно, пять или три состояния равновесия, что определяется возможным числом точек пересечения кривых. На рис. 3 это продемонстрировано для случая увеличения гидравлического сопротивления реактора АP3(Ql+Q2) (соответственно кривые 1, 2, 3).
Важным качественным результатом, представляющим интерес с точки зрения общих динамических свойств системы циркуляции как ГС частного вида, является возможность существования нескольких состояний равновесия. Проектному расчету отвечает состояние равновесия, соответствующее рабочей точке на участке ММ характеристики каждой петли на рис. 3. Работа системы циркуляции в стационарных режимах, соответствующих другим состояниям равновесия, с технической точки зрения недопустима, так как при этом в одной из петель расход теплоносителя либо мал, либо отрицателен.
Алгоритм графического решения системы (18) был реализован числено для двух-, трех- и четырехпетлевой системы циркуляции. При этом функция Ф-1, так же как и гидравлические характеристики участков АР,(в,, ю,), , = 1, п , и общего участка, были построены для задаваемых с терминала коэффициентов, определяющих вид гидравлических характеристик.
Расчеты при различных значениях сопротивлений реактора, так же как и при различных значениях ю,, позволяют получить наглядное представление об изменениях структуры фазового пространства при бифуркациях в виде рождения и исчезновения состояний равновесия. На рис. 4 слева направо представлены: результаты решения задачи потокораспределения, структура фазовой плоскости и вид функции R для случая пяти состояний равновесия (число петель системы циркуляции п = 2).
На рис. 5 представлены результаты решения задачи стационарного потокораспределения для случая четырех петель и одинаковых зависимостей — Лр (в,, ш0),, = 1, 4 .
Рис. 4
Расчеты в случае четырех петель проводились при следующих заданных параметрах характеристик ЛР,, , = 1, 5, системы (16) (см. за-
0000 висимости (3) и (4)): ш 1 = ш 2 = ш 3 = ш 4 = 1;
аъ = 2; а2, = - 4; Ь, = 3; с, = 6; ^ ^ = 0.4. Гид-
равлические характеристики для всех петель полагались одинаковыми, поэтому коэффициенты а1,, а2,, Ь,, с, зависимостей ЛР,, , = 1, 4, совпадают.
Далее, задача нахождения стационарного потокораспределения четырехпетелевой системы циркуляции была решена путем нахождения экстремумов 4-мерной функции Ляпунова R с помощью описанной выше многошаговой схемы редукции размерности. Были найдены 5 состояний равновесия:
в1 = -0.72, 02 = 1.52, 03 = 1.52, в4 = 1.52;
RmIn= -19.14; в1 = 1.52, в2 = -0.72, в3 = 1.52, вл = 1.52;
Rmm= -19.14; в1 = 1.52, в2 = 1.52, 03 = -0.72, 64 = 1.52;
Rmm= -19.14; в1 = 1.52, 02 = 1.52, в3 = 1.52, 04 = -0.72;
RmIn= -19.14; в1 = 1.04, 02 = 1.04, в3 = 1.04, 04 = 1.04.
Rmm = -18.85.
Результаты численного поиска экстремумов функции R (в1, в2, в3, в4) для указанного набора параметров системы имеют полное соответствие с результатами, полученными путем графического решения системы (18) и представленными на рис. 5. Функция R (в1, в2, в3, в4) для численных расчетов выбиралась в виде (17), где в5 = в1 + в2 + 63 + в4 , а нелинейные
функции АР, (в ю,), , = 1,4 , и АР5(в5) определялись соответственно выражениям (3) и (4) при указанных выше значениях параметров а1,; а2,, Ь,, с,, , = 1, 5, d1 и d2.
Заключение
Таким образом, на примере системы циркуляции теплоносителя ЯЭУ показана работоспособность разработанного алгоритма поиска минимумов функции R для решения задачи нахождения стационарного потокораспределения. Следует отметить, что в общем случае ГС использование указанного выше алгоритма исключено.
Приведенные выше теория и пример расчета стационарных режимов ГС относятся к методике нахождения расходов среды на участках. Эта методика в общем случае сводится к поиску
минимумов функции Ляпунова, представляющей собой сумму интегралов нелинейных гидравлических характеристик участков. Необходимое для этого исключение избыточных переменных приводит к отсутствию зависимости этой функции от давлений в узлах ГС. Однако стационарные значения избыточных переменных и давлений в узлах при необходимости могут быть определены с помощью соответствующих алгебраических уравнений.
В современных системах теплоснабжения используются так называемые замкнутые схемы. Рассмотренный выше пример системы циркуляции теплоносителя ядерного реактора можно считать также и примером системы теплоснабжения, в которой источником тепла является не ядерный реактор, а теплообменник, в котором тепло от теплоносителя центральной сети передается теплоносителю локальной замкнутой сети. Петель с насосами в этой локальной сети может быть несколько. В рассмотренном примере их четыре.
Интересующее практику количество передаваемого на каждом участке тепла, кроме расхода теплоносителя, определяется задаваемыми или рассчитываемыми поверхностями теплообмена, коэффициентами теплопередачи и так называемым температурным напором, то есть разностью температур теплоносителя и внешней среды. Таким образом, для решения проблемы теплоснабжения результат рассмотренного выше решения задачи стационарного по-токораспределения является главным [15, 16].
Для проектирования и эксплуатации систем теплоснабжения представляет также интерес оценка напряженно-деформированного состояния трубопроводов, теплообменного и другого оборудования. Для решения связанных с проектированием и эксплуатацией проблем обеспечения прочности и отсутствия кипения теплоносителя необходимо знать не только расходы на участках, но и давление в узлах. Для нахождения этих давлений в узлах используются найденные на предыдущем этапе значения расходов на участках 0°, , = 1, N , и уравнения статики, получающиеся из системы (1) при (2, = 0, , = .
Работа выполнена в рамках ведомственной целевой программы «Развитие научного потенциала высшей школы» (проект 2.1.2/3863) и федеральной целевой программы «Научные и научнопедагогические кадры инновационной России» (гос-контракт № 02.740.11.5018).
Список литературы
1. Чугаев Р.Р. Гидравлика: Учебник для вузов. 4-е изд., доп. и перераб. Л.: Энергоиздат, 1982. 672 с.
2. Смирнов Л.В., Данилова Н.В. Основы прикладной аналитической гидромеханики напорного течения несжимаемой жидкости: Учеб.-метод. пособие. Н. Новгород: Нижегор. гос. ун-т, 2009. 65 с.
3. Степанов А.И. Центробежные и осевые насосы. М.: Физматгиз, 1960. 375 с.
4. Смирнов Л.В. Динамика гидромеханических процессов в гидросистемах. Основы прикладной аналитической гидромеханики // 5 Межд. конф. «Проблемы колебаний» (ICOVP-2001), Москва, ИМАШ, 8-10 октября 2001 г. Сборник докл. М.: ИМАШ, 2002. С. 416-420.
5. Смирнов Л.В. Математические модели динамики и устойчивость систем принудительной циркуляции теплоносителя. М.: Энергоатомиздат, 1992. 128 с. (Физика и техника ядерных реакторов. Вып. 44).
6. Меренков А.П., Хасилев В.Я. Теория гидравлических цепей. М.: Наука, 1985. 230 с.
7. Лурье А.И. Аналитическая механика. М.: Гос. изд. физ.-мат. лит., 1961. 825 с.
8. Андронов А.А., Витт А.А., Хайкин С.Э. Теория колебаний. М.: Наука, 1959. 917 с.
9. Горяченко В.Д. Элементы теории колебаний. Учеб. пособие для вузов. 2-е изд., перераб. и доп. М.: Высш. шк., 2001. 395 с.
10. Стронгин Р.Г., Гергель В.П., Городецкий С.Ю. и др. Современные методы принятия оптимальных решений: Учебник. Н. Новгород: Изд-во Нижег. ун-та, 2002. 189с.
11. Strongin R.G., Sergeyev Ya.D. Global Optimization with Non-Convex Constraints Sequential and Parallel Algorithms. Dordrecht. Kluwer Academic Publishers, the Netherlands, 2000. 728 p.
12. Городецкий С.Ю., Гришагин В.А. Нелинейное программирование и многоэкстремальная оптимизация: Учеб. пособие. Н. Новгород: Изд-во ННГУ, 2007. 489 с.
13. Добряев Д.Н., Гришагин В.А. Программная система конструирования и исследования характеристических методов многоэкстремальной оптимизации // Матер. конф. «Технологии Microsoft в теории и практике программирования», Н. Новгород, 3-4 апреля 2007 г. С. 68-71.
14. Смирнов Л.В., Яскеляин А.В., Овчинников В.Ф. и др. Динамические свойства системы циркуляции теплоносителя первого контура ЯЭУ // Вопросы атомной науки и техники. Сер. Физика ядерных реакторов. 1991. Вып.3. С. 25-31.
15. Соколов Е.Я. Теплофикация и тепловые сети: Учебник для вузов. 5-е изд., перераб. М.: Энергоиз-дат,1982. 360 с.
16. Ионин А.А., Хлыбов Б.М., Братенков В.Н. Теплоснабжение. М.: Стройиздат, 1982. 336 с.
APPLICATION OF ANALYTICAL HYDROMECHANICS AND OPTIMAL DECISION-MAKING METHODS IN THE PROBLEM OF FLOW DISTRIBUTION IN HYDRAULIC SYSTEMS
L. V. Smirnov, V.A. Grishagin, D.N. Dobryaev, N. V. Danilova
The problem of a stationary flow distribution in an arbitrary hydraulic system has been shown to be reduced to a finite-dimensional multiextremal problem of the unconstrained optimization. The proposed algorithm’s efficiency has been shown by the example of the four-looped circulation system of a nuclear power plant coolant system.
Keywords: hydraulic system, stationary flow distribution, analytical hydromechanics.
