CC BY
116
Вычислительные технологии Том 14, № 2, 2009
Многорежимная оптимизация формы рабочего колеса гидротурбины*
Д.В. Банников Новосибирский государственный университет, Россия e-mail: denis. bannikov@gmail. com
С. Г. Черный, Д. В. Чирков Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: [email protected]
В. А. Скороспелов, П. А. Турук Институт математики им. С.Л. Соболева СО РАН, Новосибирск, Россия
Представлено обобщение системы автоматического проектирования формы рабочих колес гидротурбин [1] на случай многорежимной оптимизации. В ходе проектирования одновременно улучшаются характеристики 3Б обтекания рабочих колес на нескольких режимах работы, рассмотрены их новые целевые функционалы и параметризация. Представлены результаты двухрежимной оптимизации формы рабочих колес Сангтудинской ГЭС.
Ключевые слова: оптимизация формы, несжимаемая жидкость, гидравлические турбины, многоцелевой генетический алгоритм.
Введение
В [1-7] предложен и развивается подход к оптимизационному проектированию проточных частей гидротурбин, основанный на последовательности расчетов течений для многих тысяч вариаций их геометрий и выборе такой формы, которая обеспечивает минимум одного или нескольких заданных целевых функционалов. Отличительная особенность предложенного метода автоматического проектирования от существующих — применение в нем для расчета трехмерных течений воды в проточных частях гидротурбин уравнений Эйлера несжимаемой жидкости. Классическая модель идеальной жидкости не может быть использована для расчета потерь энергии, а следовательно, и для определения коэффициента полезного действия гидротурбины. Диссипация энергии отсутствует в описываемых течениях жидкости. Процесс сводится к комбинации поступательного движения и вращения жидкости как твердого тела. Однако диссипативные механизмы, отсутствующие в уравнениях Эйлера, неизбежно появляются в их разностных
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 08-01-00364).
© ИВТ СО РАН, 2009.
аналогах (численных моделях уравнений Эйлера), используемых для приближенного решения исходных дифференциальных уравнений. Указанные механизмы играют важную роль при сходимости итераций предельных решений от начальных данных. Причем для отбора требуемого решения важен не конкретный вид “диссипативных” членов, появляющихся в разностных уравнениях и зависящих от выбранной разностной схемы, а само наличие “механизма диссипации”. Таким образом, численная модель уравнений Эйлера, являясь чрезвычайно экономичной, позволяет без использования условий Жуковского—Чаплыгина выделять близкое к действительности решение задачи обтекания решеток профилей, характерной для определения течения в гидротурбине. Однако корректное вычисление КПД с ее помощью невозможно, хотя качественное поведение этой важной характеристики передается данной моделью хорошо. Поэтому при ее использовании в методе автоматического проектирования проточных частей особо актуальной становится задача выбора подходящих целевых функционалов, которые за счет улучшения кинематических свойств потока косвенно уменьшали бы потери энергии в гидротурбине.
Отметим ряд публикаций [8-12], идейно близких к использованному авторами подходу автоматизации проектирования формы рабочего колеса (РК) — основного элемента проточного тракта гидротурбины. Во всех работах расчетная область ограничена частью проточного тракта между соседними лопастями РК. Моделирование стационарного течения выполняется в предположении, что в остальных межлопастных каналах течение циклически повторяется. Для моделирования потока обычно используются уравнения Навье—Стокса совместно с (к—е)-моделью турбулентности.
Авторам известна только одна работа [8], в которой при оптимизации формы РК также используются уравнения Эйлера несжимаемой жидкости, а выбор оптимальной геометрии проводится по косвенным критериям, которые учитываются с разными весами в одном целевом функционале. Недостатком подхода одноцелевой оптимизации является то, что значения весовых коэффициентов заранее неизвестны, и приходится многократно решать оптимизационную задачу с различными комбинациями весов.
В указанных выше собственных работах авторов для поиска геометрий, улучшенных одновременно по нескольким критериям качества, формулируется задача многоцелевой оптимизации, решение которой ищется с помощью генетического оптимизационного алгоритма (ГА). Предложены новые критерии качества работы гидротурбины. С их помощью удается обеспечить отсутствие противотоков на входной кромке лопасти, минимизировать область возможной кавитации, создать требуемое распределение параметров потока. Вариация формы лопасти осуществляется 16 угловыми координатами ее срединной поверхности.
В [9, 10] при оптимизации формы РК гидротурбины предлагается не только рассматривать гидравлические характеристики проточного тракта, но и проводить прочностной анализ каждой геометрии. Движение жидкости моделируется в рамках уравнений Навье—Стокса, а для анализа напряженного состояния лопастей используется упрощенная модель упругого равновесия. Предположения о непрерывной зависимости рассматриваемых критериев от варьируемых параметров и близости оптимальной геометрии к исходному прототипу, созданному опытными инженерами, позволяют использовать градиентный метод для быстрого нахождения экстремума. В работах решается задача оптимизации прочностных характеристик лопасти при сохранении гидравлических характеристик как у прототипа. Для вариации формы используется не более 10 свободных геометрических параметров. Единственным критерием оптимизации служит взвешен-
ная сумма прочностных и гидравлических функций качества. К недостаткам данных работ можно отнести использование одноцелевой оптимизации, малое число варьируемых параметров и локальный поиск оптимальной геометрии.
В [11] рассматривается оптимизация формы РК поворотно-лопастной гидротурбины. При этом в отличие от предыдущих работ с помощью многоцелевого генетического алгоритма отыскивается геометрия, одновременно улучшенная по двум критериям качества на одном из режимов работы гидротурбины. Критерии качества — КПД рабочего колеса и коэффициент восстановления давления. Для моделирования движения жидкости используются уравнения Навье—Стокса с (к—е)-моделью турбулентности. Число варьируемых параметров равно 10. Отмечено, что точность определения КПД на сетке с 32 000 узлов около 3 % по сравнению с экспериментом. Общее число анализируемых геометрий в оптимизационном расчете достигает 120. После нахождения оптимальной геометрии проводится ее анализ на других режимах работы гидротурбины. В [11], как ив [9, 10], при поиске оптимальной геометрии исследуется небольшое число вариаций исходной геометрии. Применяемая грубая сетка не позволяет с достаточной точностью оценивать КПД рабочего колеса.
В [12] в отличие от указанных выше работ при поиске оптимальной формы РК рассматриваются три режима работы гидротурбины: оптимального КПД, максимальной и неполной загрузок. На каждом режиме минимизируется один функционал, равный сумме потерь энергии в межлопастном канале и кинетической энергии вращающегося потока на выходе из РК. Расчетная сетка содержит 112 000 ячеек. При варьировании 30 геометрических параметров проанализировано 2000 различных модификаций при помощи ГА. В результате получено множество решений оптимизационной задачи, составляющих фронт Парето. В статье при оптимизационном проектировании не учитываются кавитационные характеристики РК, существенно воздействующие на работу всей турбины в режиме максимальной загрузки. Также на режиме неполной загрузки за РК возникает прецессирующий вихревой жгут, который не может быть описан в используемой в работе постановке. В то же время его влияние на улучшаемые характеристики значительно.
В настоящей работе предлагается дальнейшее развитие предложенного ранее авторами подхода. В частности, расширена вариация формы РК за счет параметризации поверхностей обода и ступицы, входной и выходной кромок лопасти. Представлены формулировки новых целевых функционалов и ограничений. В предыдущих работах авторов поиск оптимальных геометрий проводился для одного предварительно заданного режима работы гидротурбины. Однако одним из требований при проектировании проточной части гидротурбины является получение приемлемых характеристик в широком диапазоне режимов ее работы, каждый из которых обладает своими особенностями. Так, режим максимальной мощности характеризуется сильной кавитацией на тыльной стороне лопасти, размер которой необходимо минимизировать. На режиме максимального КПД следует обеспечить отсутствие сильных поперечных течений на лопасти и равномерность поля скоростей на входе в отсасывающую трубу. В данной работе рассмотрена постановка задачи многорежимной оптимизации, позволяющая находить форму РК, удовлетворяющую нужным критериям качества одновременно на каждом из заданных режимов работы гидротурбины. С помощью представленного алгоритма проведено оптимизационное проектирование формы радиально-осевого РК на основе решения двух- и трехцелевой однорежимной, а также трехцелевой двухрежимной оптимизационных задач.
1. Обобщенная постановка задачи автоматизированного проектирования и особенности ее решения
1.1. Общая схема решения задачи оптимизационного проектирования с учетом особенностей нескольких режимов
Здесь формулируется постановка задачи многорежимной оптимизации (МРО), обобщающая предложенные ранее авторами постановки многоцелевой оптимизации (МО). Описываются детали алгоритма, присущие задаче МРО.
Пусть имеется геометрическая форма проточной части некоторого прототипа гидротурбины, которая подлежит оптимизационному проектированию. Под этим в статье подразумевается такая модификация формы прототипа, после которой работа гидротурбины на каждом из заданных Ь режимов будет максимальным образом удовлетворять наборам критериев качеств, сформулированных для каждого из режимов. Варьирование формы проточного тракта при поиске ее оптимальной модификации осуществляется путем ее параметризации N параметрами
X =(Х1,...,ХМ ) (1)
и их изменения. Набор параметров (1) однозначно определяет геометрию проточного тракта гидротурбины и будет отождествляться в статье с самой формой.
Процесс МРО состоит из следующих шагов.
1. Случайным образом формируется N форм проточного тракта гидротурбины начального поколения в = 0:
Хт , Хк, ..., X N в
N в ^ор
Форма прототипа проточного тракта включена в это множество. Будем также называть форму индивидуумом.
2. Для каждой формы хк автоматически строится сетка, проводится Ь гидродинамических расчетов течений на ней, каждое из которых соответствует одному из Ь заданных режимов работы гидротурбины. По полям течения в форме хк 1-го режима (I = 1, ...,Ь) вычисляются М1 функционалов:
рм,...., рм, (2)
являющихся критериями качества работы гидротурбины на этом режиме. Общее количество функционалов М для одного индивидуума равно
ь
М = ^ М1. (3)
1=1
Множество функционалов
РМ ,...,РМ (4)
характеризует геометрию хк на Ь режимах ее работы.
3. Оптимизационный алгоритм по значениям N • М функционалов всего текущего
поколения в определяет лучшие формы этого поколения и, если не выполнено условие
остановки процесса, создает новое (в + 1)-е поколение
Х+1 Х+1 Х+1 (5)
Х1 ,..., хк ,..., хма+1. (5)
^фор
Число индивидуумов в в-м поколении N при переходе к (в + 1)-му поколению ^Ор1, вообще говоря, не сохраняется. Оно может как увеличиваться, так и уменьшаться. Далее осуществляется переход на шаг 2 с переприсвоением в = в + 1 до тех пор, пока не будет рассчитано требуемое число поколений.
После завершения оптимизационного процесса получается множество оптимальных решений, каждое из которых не хуже, чем прототип.
Для вычисления М функционалов, характеризующих одну геометрию, необходимо провести Ь гидродинамических расчетов. Время оптимизации при этом увеличивается в Ь раз по сравнению с однорежимным проектированием, но позволяет более полно учесть требования к проектируемой гидротурбине.
Ниже сформулирована задача многорежимной оптимизации и приведены отличительные особенности ее решения.
1.2. Математическая постановка задачи многорежимной оптимизации
Математически формулировка задачи МРО совпадает с постановкой задачи МО и состоит в нахождении вектора х, обеспечивающего минимальные значения М функционалов
тт Р1(х),..., тт РМ (х) (6)
при фазовых:
Хь,г < Хг < Хд,г, г =1,...,^ (7)
геометрических:
"фг(х) < 0, г = 1,...,/, (8)
гидродинамических:
(х) < 0, э = 1,...,^, I = 1,.., Ь, (9)
ограничениях.
Фазовые ограничения определяют пределы варьирования каждого из параметров. Геометрические ограничения позволяют исключить из рассмотрения формы, не приемлемые с инженерной точки зрения. Гидродинамические ограничения определяют допустимые требования к рассчитанному потоку жидкости, при их нарушении геометрия проточного тракта признается недопустимой. Решением задачи (6)-(9) является семейство точек х, называемое оптимальным фронтом Парето.
1.3. Отличительные особенности генетического алгоритма поиска решения задачи многорежимной оптимизации
Операторы ГА кодирования, рекомбинации, мутации и клонирования для решения задачи МРО совпадают с соответствующими операторами для МО, детально описанными в [1]. Далее рассматриваются особенности оператора селекции, характерные для задачи МРО. Предлагается новый способ выполнения ограничений оптимизационной задачи.
1.3.1. Оператор селекции
Оператор селекции выбирает из текущего в-го поколения определенную долю Тг • N индивидуумов, имеющих наименьшее значение критерия качества. Критерий качества — это функция
С(Рь...,Рм): ЯМ ^ Я, (10)
которая характеризует степень пригодности проверяемой геометрии. Формы, лучшим образом удовлетворяющие заданным критериям, должны иметь наименьшие значения функции (10) и, следовательно, участвовать в построении нового поколения.
1.3.2. Критерии качества в многорежимной оптимизации
Построение функции (10) может быть основано на различных принципах. В данной работе для решения задачи МРО предлагаются следующие два.
Критерий качества, основанный на ранге
В [1] на основе значений рассчитанных целевых функционалов каждому индивидууму Xfc сопоставляется скалярная величина, называемая рангом:
L Ji
rank(xk) = 1 + ak + NopE ЕПЦ(xk))> (11)
i=i j=1
где ak — число индивидуумов текущей популяции лучших, чем Xfc, одновременно по всем функционалам Ff,..., F_M. Последнее слагаемое в правой части (11) позволяет учитывать число невыполненных ограничений. Функция П определяется как
= 1'^> 0’ (12)
= 0, р < 0. (12)
При таком определении П за каждое невыполненное ограничение к рангу добавляется штраф, равный размеру поколения Np0p. Из двух индивидуумов всегда хуже тот,
у которого не выполнено большее число ограничений. Однако функция штрафа, учитывающая, насколько сильно не выполнены ограничения, будет лучше, чем функция, зависящая только от числа нарушенных ограничений.
Поэтому в настоящей работе критерий качества, основанный на ранге, модифицирован. Сначала в значениях вычисленных функционалов (2) с помощью метода штрафных функций учитываются невыполненные ограничения
Fk = Fk
1 + 13 Cj (^j (xk) + (xk)|)
j=1
i =1,..,Ml, (13)
где С^ — коэффициенты штрафа. Здесь, в отличие от (11), на критерии качества /-го режима влияют только ограничения этого же режима. Затем на основе всех функционалов р0,..., РМ вычисляются значения рангов по формуле
гапк(хм) = 1 + ам. (14)
Очевидно, что индивидуумы, принадлежащие фронту Парето, имеют ранг 1. Чем больше ранг, тем больше расстояние до фронта Парето. Окончательно критерий качества, основанный на ранге, формулируется как
гапк(хк) —1
Сг (хм) = 1 + п(г), (15)
г=1
где п(г) — число особей текущей популяции, имеющих ранг г.
Дальнейший отбор особей можно проводить по величине 0Г (хм), однако в этом случае может наблюдаться преждевременная сходимость алгоритма к небольшому подмножеству оптимального фронта Парето. Это происходит вследствие того, что какая-нибудь наиболее приспособленная особь и ее потомки вытесняют менее приспособленных на данный момент, но, возможно, более перспективных, индивидуумов и популяция локализуется вблизи нескольких лучших особей. Такая ситуация также может приводить к сходимости алгоритма к локальному минимуму.
Критерий качества, основанный на методе разделения
Для того чтобы избежать локального скопления популяции, а также получить равномерное распределение точек вдоль фронта Парето, используется метод разделения. В основе метода лежит следующая идея: на индивидуумы, расстояние между которыми в пространстве значений функционалов меньше, чем некоторое а8ьаге, налагаются “штрафы”. Для этого критерий качества 0Г, определенный в (15), умножается на нишевое число тм, вычисляемое для каждого индивидуума. Нишевое число тм является оценкой количества особей популяции, располагающихся в окрестности индивидуума хм. Это число рассчитывается по всей популяции, включая к-й индивидуум:
N 5
^рор
тм = 53 5Ь(р(М)). (16)
г=1
Здесь р(к, г) — расстояние между индивидуумами хм и хг, определяемое как
р(к,г) = тах Р,м — Р ; (17)
4 ' 1<1<М 1 1 47
£^(р) — функция разделения, задается как
б'ВД = (1 — р/азЬагеГ, £й(0) = 1, £Л,(р > ОЛаге) = 0, (18)
где а^ге — параметр разделения, подбираемый эмпирически. При таком определении
функции следует, что тм > 1. Индивидуумы, находящиеся на расстоянии меньшем
а^^ге, ухудшают (увеличивают) функции качества друг друга. Исследование влияния вида функции БН на работу алгоритма показало, что линейная зависимость (а = 1) позволяет равномерно распределять формы вдоль фронта Парето. В итоге отбор индивидуумов из рассчитанного поколения происходит на основе функции качества:
(гапк(х^ )-1 \
1+ 53 П(г)1 тм. (19)
1.3.3. Клонирование в многорежимной оптимизации
Для сохранения лучших форм проточной части 5-го поколения в следующем (в + 1)-м поколении в него копируются без изменения все индивидуумы ранга 1 — обозначим их р|. Таким образом, популяцию (в + 1)-го поколения
х«+1 х«+1 х«+1 х«+1 (20)
х1 >•••> хр , хр+1,--ч хр+р| (20)
составляют р рекомбинированных и мутированных индивидуумов, а также р^ особей, взятых с фронта Парето. Общее количество индивидуумов (в + 1)-го поколения есть
ЛОр1 = р + р|. (21)
1.3.4. Задание параметров генетического алгоритма МРО
Операторы ГА МРО определяются следующими параметрами [1]: селекции — р, Тг, рекомбинации — мутации — ц и 8. Выбор этих параметров существенно влияет на эффективность его работы. Размер рекомбинированных индивидуумов р влияет на скорость и точность поиска фронта Парето. Если значение р слишком мало, то пространство решений не может быть исследовано должным образом и возникают сложности с определением положения оптимального фронта Парето. Если р слишком большое, то ГА будет расходовать вычислительные ресурсы на излишних особей, что приведет к неоправданно долгому времени оптимизации. Параметр Тг характеризует относительную часть популяции, информация о которой будет использоваться при построении нового поколения. Задание слишком маленького значения Тг (Тг < 0.2) приводит к тому, что все потомки получаются похожими на родителей и алгоритм очень быстро сходится к локальному минимуму. Слишком большое значение Тг (Тг > 0.8) приводит к тому, что выбор родительских индивидуумов будет происходить случайным образом из всей популяции. В этом случае алгоритм очень долго сходится, так как нет явного направления поиска и развитие популяции похоже на случайное блуждание. Параметр d определяет максимальный размер области, в которой может лежать потомок. При выборе значений d < 0 все получаемые значения параметров лежат “между” родительскими, и тем самым алгоритм имеет тенденцию очень быстро стягиваться в точку. Значения параметра d > 1 приводят к слишком сильному “разбрасыванию” потомков и замедляют сходимость алгоритма к оптимуму. Параметры мутации ц, 8 показывают, как далеко и с какой вероятностью могут измениться параметры потомка, полученного в результате рекомбинации двух родительских особей. Величину мутации ц следует задавать достаточно большой (±10%), а вероятность 8 — относительно небольшой (менее 0.1). Таким образом, случайные большие изменения параметров дают возможность получать принципиально новые решения. Параметр а8ьаге характеризует размер шара в пространстве функционалов, в котором не желательно нахождение других индивидуумов. Существуют попытки теоретически оценить величину а^ге в зависимости от размерности пространства решений и числа найденных решений [13]. Параметр а8ьаге может задаваться до начала оптимизации либо корректироваться в процессе счета на каждом поколении. Обычно выбирается первый вариант. Поскольку целевые функционалы могут принимать значения на отрезках, имеющих разные масштабы, то для каждого функционала необходимо задать свой а8ьаге. Чтобы этого избежать, в данной работе все функционалы предварительно нормируются на отрезок [0, 1].
Генетический оптимизационный алгоритм не гарантирует сходимости к глобальному оптимальному решению, однако правильное задание параметров алгоритма увеличивает шансы на успех. Влияние параметров р, Тг, d, ц, 8 на работу построенного алгоритма в случае многоцелевой оптимизации исследуется в [1].
Приведем значения параметров, используемых в данной работе при оптимизации формы РК: N = 23; Л£ор = 200; р = 70 ... 100; Тг = 0.3; d = 0.7; ц = 0.1; авЬаге = 0.01.
2. Конкретизация предложенной постановки
к оптимизационному проектированию формы рабочего колеса
2.1. Параметризация рабочего колеса
Во всех предыдущих работах авторов варьирование геометрии проточного тракта РК осуществлялось только посредством вариации угловой координаты срединной поверхности лопасти. В настоящей работе к этой параметризации добавляется параметризация меридиональной проекции РК.
2.1.1. Параметризация меридиональной проекции рабочего колеса
Параметризация меридиональной проекции РК (рис. 1) позволяет изменять форму осесимметричных поверхностей обода и ступицы, а также кривизну и положение входной и выходной кромок лопасти в меридиональной проекции.
Входная кромка лопасти варьируется двумя параметрами — р1 и р4, задающими касательные к ней в точках С1 и С4 векторы Н1 и Н4 соответственно (рис. 1). В предельном случае она может совпадать с отрезком С1С4 или должна быть выгнута в направлении потока с постоянным знаком кривизны.
Выходная кромка лопасти также варьируется двумя параметрами — р2 и р3, задающими касательные к ней в точках С2 и С3 векторы Н2 и Н3 соответственно.
При описании формы поверхности обода фиксируются положения точек А и С2, а также значения диаметров ^1 и ^2. В точках А и С2 задаются и не меняются направления касательных к ободу векторов. В частности, касательный вектор в точке С2 определяет угол конуса отсасывающей трубы а. Варьируется только один параметр — р5, задающий касательное направления Н5 к ободу в точке С1.
Форма поверхности ступицы зависит от одного параметра — р6, задающего касательный к ней в точке В вектор Н6.
Диаметры ^3 и ^4 определяются параметрами р7 и р8 соответственно. Они характеризуют положение точек С3 и С4 на линии ступицы.
Для примера рассмотрим, как происходит варьирование формы входной кромки РК. Вводится локальная прямоугольная система координат Хм Ум, ось Хм которой
Рис. 1. Меридиональная проекция рабочего колеса
совпадает с прямой С1С4. Уравнение модифицированной линии входной кромки С1С4 в системе координат Хм Ум выглядит так:
где Хм, Ум — координаты исходной линии ступицы; х*м, у*м — координаты модифицированной линии в системе координат Хм Ум; 0(и) — линейная функция, удовлетворяющая условиям: "0(0) = рі, 0(1) = р4, где — варьируемые параметры. Таким образом,
параметр р1 определяет направление касательной к ступице в точке С1 (направление вектора Н1), параметр р4 — направление касательной в точке С4. Очевидно, значение параметров р1 = 1, р4 = 1 (0(и) = 1) соответствует исходной форме входной кромки. На параметры р1 и р4 накладывается ограничение вида
гарантирующее выпуклость модифицированной линии входной кромки. Параметры а и Ь зависят от геометрии прототипа. Например, а = 0.6 и Ь = 2.
2.2. Гидродинамические ограничения
Гидродинамические ограничения должны обеспечивать выполнение режимных и кавитационных требований к потоку жидкости в РК.
2.2.1. Режимные ограничения
Режим работы РК задается тремя параметрами: расходом ^, числом оборотов в минуту п и напором Нс в рабочем колесе.
Расход выполняется за счет задаваемого во входном сечении поля скоростей:
где V — вектор абсолютной скорости; Ж — вектор нормали к элементарной площадке входного сечения, величина которого совпадает с ее площадью.
Частота вращения п учитывается в источниковом члене уравнения количества движения в виде кориолисова 2ш х V и центростремительного ш х ш х г ускорений, где ш = (0, 0,ш) — вектор угловой скорости, ш = 2п/п, г — радиус-вектор точки от оси вращения колеса.
Напор определяется как разность полных удельных энергий на входе Е1п и выходе
где р — давление, г — положение элементарной площадки по оси Z, 7 = рд, д — ускорение свободного падения. Вообще говоря, после расчета поля течения напор может не соответствовать режимному значению. Оценка потерь в РК дает
(22)
ар1 < р4 < Ьр1,
(23)
1П
(24)
Еоиі из РК:
Нс Е1п Еоиі,
(25)
где Н — полный напор гидротурбины. Для сохранения заданного напора станции ставится обязательное ограничение:
|Нс - 0.95Н| < е, (26)
где е — малый параметр.
2.2.2. Кавитационное ограничение
Кавитационные явления, наступающие при понижении давления на поверхности лопасти ниже давления насыщенного пара рсау:
Р < Рсау, (27)
существенно влияют на срок эксплуатации РК. Если область кавитации занимает более 15 % площади всей тыльной стороны лопасти, наблюдается резкое падение КПД турбины. Поэтому необходимо минимизировать область кавитации либо добиваться ее смещения как можно ближе к выходной кромке лопасти. Первое обеспечивается ограничением
5^ < 0.15, (28)
^БИС
где 5сау — площадь зоны выполнения (26), 58ис — площадь тыльной стороны лопасти. Второе условие достигается ограничением
5*ау = I ^са^5 < 0, (29)
-Фэис
где $сау = 1 при выполнении условия (26) и $сау = 0 в противном случае. Весовая
функция ,ш = и>(7), зависящая от расстояния до входной кромки, управляет положением
зоны кавитации.
2.3. Целевые функционалы
Основные целевые функционалы, рассматриваемые в настоящей работе, предложены и обоснованы ранее в [1, 2, 5]. Здесь предлагается и исследуется дополнительно к ним новый функционал, минимизирующий силовое воздействие потока на входную кромку лопасти. Поскольку при описании результатов оптимизационного проектирования будут упоминаться как ранее введенные, так и новый функционалы, то для удобства изложения приведем их краткие описания.
2.3.1. Кинетическая энергия в выходном сечении рабочего колеса
Кинетическая энергия потока в выходном сечении 5ои4 расчетной области вычисляется по следующей формуле:
^ = 2дЬ / м'2 (у ■ж)' (30)
где V — вектор абсолютной скорости; Ж — вектор нормали к элементарной площадке ^5 выходного сечения, равный по модулю площади этой площадки и направленный вовне
расчетной области; Q — расход. Минимизация функционала эквивалентна выравниванию профиля скорости в сечении 5ои4 и направлена, таким образом, на повышение эффективности работы отсасывающей трубы. Уменьшение закрутки потока приводит к уменьшению циркуляционных потерь.
2.3.2. Отклонение линий тока от осредненного по окружному направлению потока
Для уменьшения “профильных” потерь энергии в лопастной системе РК предлагается приближать предельные линии тока на поверхности лопасти к линиям тока осесимметричного потока. Для этого необходимо потребовать, чтобы |в| ^ 0, где в — угол
между предельной линией тока на поверхности лопасти и линией тока осесимметрично-
го потока. Условие отсутствия возвратных течений в терминах угла в можно выразить неравенством |в| < п/2.
Указанные требования оцениваются с помощью следующего функционала:
F2 = 1 / (1 - ст(в)в)dS, (31)
S
где S — площадь поверхности лопасти; весовая функция а имеет вид
а(в)= {0;в > п/2: (32)
Параметр а0 = const > 1 характеризует степень выполнения требования по отсутствию возвратных течений на лопасти.
2.3.3. Силовое воздействие потока на входную кромку лопасти
Для уменьшения “ударных” потерь на входе в лопастную систему предлагается добиваться безударного угла натекания потока на лопастную систему. Для этого необходимо минимизировать силовое воздействие потока на входную кромку лопасти в направлении, перпендикулярном входному элементу срединной поверхности. Сформулированные условия учитываются в функционале ^3:
*3 = У |Др„| (33)
Звх
где Дрп — проекция перепада давления на нормаль к входному элементу срединной поверхности лопасти, £вх — поверхность входного элемента (рис. 2). Длина входного элемента срединной поверхности лопасти определяется как /0 = $/, где / — длина срединной поверхности лопасти.
2.3.4. Расположение области кавитации
Если постановкой кавитационных ограничений (28), (29) не удается получить лопасть, имеющую область кавитации желаемого размера, можно добиваться минимально возможного размера этой области либо пытаться смещать ее как можно ближе к выходной
Рис. 2. К определению целевого функционала
кромке. В этом случае кавитационные ограничения (28), (29) формулируется в виде целевых функционалов:
ОС
О.
О:
О,
*
еау
О.,
(34)
(35)
где £сау, 5^ определены в п. 2.2.2. Функционал ^4 позволяет минимизировать область кавитации. Функционал ^5 с использованием весовой функции позволяет смещать область кавитации к выходной кромке лопасти.
3. Результаты оптимизационного проектирования формы рабочего колеса
Логика расчетов, результаты которых приводятся ниже, следующая. Проводились многоцелевые, но однорежимные оптимизационные расчеты. В отличие от ранее сделанных расчетов [1], использовалась расширенная параметризация РК. Анализировались полученные результаты и отмечались их недостатки, связанные с проектированием при учете одного режима. После этого для их устранения осуществлялась двухрежимная оптимизация геометрии.
Моделирование потока жидкости проводилось для гидротурбины Сангтудинской ГЭС с диаметром РК О = 1 ми напором Н = 1 м. Рассмотрены два наиболее важных режима работы гидротурбины:
I) максимального КПД (ф = 0.97, п = 77);
II) номинальной мощности (ф = 1.18, п = 78.8).
Конечно-разностная сетка, выстраиваемая в одном межлопастном канале РК, имела размеры 50 х 20 х 20 в продольном, от обода к ступице и окружном направлениях соответственно. Влияние способа постановки граничных условий, размеры расчетной области и сгущение сетки детально исследованы в [5]. Следуя этой работе, поле скорости
на входе в расчетную область бралось из предварительного расчета течения в направляющем аппарате совместно с РК. Во время всего оптимизационного процесса эти параметры полагались неизменными. В выходном сечении задавалось условие радиального равновесия для давления. Варьировалось 23 геометрических параметра: 8 параметров, изменяющих меридиональную проекцию РК, и 15 угловых параметров, характеризующих кривизну срединной поверхности лопасти. Время решения задачи двухрежимной оптимизации на кластере ИВТ СО РАН (4х2 Intel Xeon 3.0 GHz) с использованием восьми процессоров составило около 60 ч.
3.1. Двухцелевая оптимизация на режиме I
Моделирование потока в проточном тракте прототипа на режиме I показало наличие области возвратного течения на входной кромке лопасти и сильных поперечных течений по рабочей поверхности лопасти (рис. 3, а). В то же время наблюдается достаточно равномерное поле скорости на выходе из РК и отсутствие кавитации на поверхности лопасти.
В качестве целевых функционалов задачи (6)-(9) рассматривались F (30) и F2 (31). Дополнительно было установлено ограничение на отсутствие кавитации Scav = 0. Данная постановка задачи рассматривалась как попытка устранить возвратные течения на входной кромке лопасти, приблизить поток к осесимметричному и, если возможно, улучшить профиль скорости на выходе из РК, а также сохранить кавитационные качества, как у прототипа.
Рис. 3. Линии тока на рабочей стороне лопасти (вверху) и меридиональные сечения РК с профилем скорости в выходном сечении (внизу): а — прототип; б — ОрЬ-1; в — Орі-2; г — Орі-3; 1 — прототип; 2 — после оптимизации
14
Орі-1
Орі-2
б
11.8
12
12.2
12.4
1
Рис. 4. Лучшие индивидуумы для различных поколений
На рис. 4 в плоскости минимизируемых функционалов показаны текущие фронты Парето для 0-, 5-, 10-, 20-, 30-, 60-го поколений. Визуальная сходимость фронта была получена на 40-м поколении. В дальнейшем фронт Парето перестал двигаться, и алгоритм более тщательно исследовал найденный фронт. Так, после 40-го поколения фронт Парето содержал 25 решений, а после 60-го — 75. Отметим равномерность расположения найденных решений вдоль оптимального фронта Парето. Это позволяет выбирать решения, в различной мере удовлетворяющие заданным требованиям. Выбор компромиссной геометрии из всего найденного множества решений остается за проектировщиком.
Для анализа полученного множества решений рассматривались РК, соответствующие крайним точкам фронта Парето. Геометрия ОрЬ-1 соответствовала наименьшему значению ^1, ОрЬ-2 — ^2. На рис. 3 приведены предельные линии тока по поверхности лопасти и меридиональные проекции РК. Уменьшение функционала ^2 привело к приближению линий тока к осесимметричному потоку и совместило линию растекания с входной кромкой (рис. 3, в), в отличие от лопасти, оптимальной по (рис. 3, б). Геометрия, оптимальная по функционалу , имела наиболее равномерный профиль скорости на выходе из РК.
3.2. Трехцелевая оптимизация на режиме I
При проведении трехцелевой оптимизации на режиме I одновременно минимизировались функционалы (30), ^2 (31) и ^3 (33). Как и при двухцелевой оптимизации, ставилось ограничение на величину напора и отсутствие кавитации. Число рекомбинируемых индивидуумов р в ГА было увеличено до 100. После расчета 60 поколений получена визуальная сходимость фронта Парето, содержащего 223 решения.
На рис. 5 представлен полученный фронт Парето в плоскости функционалов и ^2. На нем функционал ^3 изображен оттенками серого. Отметим, что все три функционала являлись конфликтующими, т. е. уменьшение одного из них приводило к росту хотя бы одного из других функционалов. Там же представлено положение начальной геометрии и для сравнения кругами приведен фронт Парето, полученный при двухцелевой оптимизации (п. 3.1). Видно, что найденные фронты достаточно близки и имели
Рис. 5. Сравнение фронтов Парето, полученных при двухцелевой (кружки) и трехцелевой (закрашенная область) оптимизациях
практически одинаковый охват максимальных и минимальных значений функционалов и
В качестве оптимальной была выбрана лопасть 0рЬ-3, которая существенно улучшена по второму и третьему критериям, а величина кинетической энергии (значение ) такая же, как у прототипа (рис. 5). На рис. 3, г представлены результаты моделирования потока в РК 0рЬ-3.
3.3. Оптимизация на режимах I и II
На рис. 6, а представлены результаты расчета течения в прототипе РК на режиме
II, который показал, что в этом случае на тыльной стороне лопасти возникает зона кавитации площадью 19%. Геометрии, полученные при однорежимной оптимиза-
Рис. 6. Обтекания рабочей (слева) и тыльной (справа) стороны лопастей прототипа (а) и оптимального РК ОрЬ~4 (б): область кавитации обозначена черным цветом
Рис. 7. Сравнение фронтов Парето, полученных при двухрежимной оптимизации (точки) и двухцелевой однорежимной (кружки) оптимизации
ции, также имеют области кавитации на режиме II, например, геометрия 0рЬ-1 — 8%, 0рЬ-2 — 24, 0рЬ-3 — 12 %.
Для того чтобы получить РК без областей кавитации и возвратных течений на обоих режимах, была проведена двухрежимная оптимизация. На режиме I минимизации подвергались функционалы (30) и Г2 (31) при ограничении на напор и размер области кавитации ScalV < 0. На режиме II минимизировался функционал ^5 (35). Было рассчитано 60 поколений с числом рекомбинируемых индивидуумов р = 100 в каждом. На рис. 7 представлен фронт Парето в плоскости функционалов — ^2. Отметим, что зависимость от функционала ^5 слабая, т. е., немного ухудшив геометрию относительно функционалов — Г2, можно добиться отсутствия кавитации на тыльной стороне
лопасти.
Для сравнения на рис. 7 приведено множество решений, полученное при однорежимной двухцелевой оптимизации. И здесь мы наблюдаем достаточно близкое расположение фронтов Парето. Строго говоря, фронт Парето, полученный при двухцелевой оптимизации, должен быть частью фронта Парето трехцелевой оптимизации. Несовпадение может быть связано с недостаточно большим числом рекомбинируемых индивидуумов р. С увеличением числа критериев и ограничений оптимизационной задачи возрастает сложность определения оптимального фронта Парето.
В качестве оптимальной была выбрана геометрия 0рЬ-4, у которой отсутствовала кавитация на обоих режимах; геометрия 0рЬ-4 существенно улучшена по критерию Г2 и не хуже прототипа по критерию (см. рис. 6, б).
Заключение
В работе рассмотрен один из возможных подходов к формализации процесса проектирования элементов гидротурбин путем постановки и решения соответствующей многорежимной оптимизационной задачи. Для решения этой задачи в системе автоматизированного проектирования применена обобщенная версия ГА многоцелевой оптимизации. ГА обычно оказываются существенно эффективнее, чем традиционные градиентные методы, при решении многокритериальных практических задач. Больший объем вычис-
лений, требуемый для получения оптимального решения в случае однокритериальной оптимизации, компенсируется легким и естественным обобщением ГА на многоцелевой случай и возможностью применять данный алгоритм на многопроцессорных системах. Проведено двухрежимное оптимизационное проектирование формы РК относительно заданных критериев.
В настоящий момент выполнено распараллеливание ГА, что позволило существенно сократить время оптимизационного расчета. В будущем планируется рассмотреть возможность оптимизации элементов гидротурбины в рамках одновременных расчетов течения в направляющем аппарате, РК и отсасывающей трубе. А также осуществлять расчет течения в рамках турбулентных моделей [14] и выбирать оптимальную форму непосредственно по виду зависимости КПД гидротурбины от режима ее работы.
Список литературы
[1] Лобарева И.Ф., Черный С.Г., Чирков Д.В. и др. Многоцелевая оптимизация формы лопасти гидротурбины // Вычисл. технологии. 2006. Т. 11, № 5. С. 63-76.
[2] Skorospelov V.A., Turuk P.A., Aulchenko S.M. et AL. Solution of the 3D optimization problem of the aerohydrodynamic shape of turbine components // Proc. of Xll-th Intern. Conf. on the Methods of Aerophysical Research, Novosibirsk, 28 June—3 July, 2004. Vol. 2. P. 172-177.
[3] Лобарева И.Ф., Скороспелов В.А., Турук П.А. и др. Об одном подходе к оптимизации формы лопасти гидротурбины // Вычисл. технологии. 2005. Т. 10, № 6. C. 52-73.
[4] Cherny S., Chirkov D., Lapin V. et al. 3D Euler flow simulation in hydro turbines: unsteady and automatic design // Springer Series: Notes on Numerical Fluid Mechanics and Multidisciplinary Design. 2006. Vol. 93. P. 33-51.
[5] Черный С.Г., Чирков Д.В., Лапин В.Н. и др. Численное моделирование течений в турбомашинах. Новосибирск: Наука, 2006. 202 с.
[6] Bannikov D.V., Lobareva I.F., Cherny S.G. et al. Multiobjective optimization in problems of hydrodynamics of turbomachines // Proc. of XIII-th Intern. Conf. on the Methods of Aerophysical Research, Novosibirsk, 5-10 February, 2007. Part 1. P. 22-27.
[7] Пылев И.М., Малышев А.К., Черный С.Г., Скороспелов В.А. Оптимизационное проектирование проточных частей гидротурбин // Тяжелое маш-ние. 2007. № 4. С. 10-13.
[8] Sallaberger I., Fisler M., Michaud M. et al. The design of Francis turbine runners by 3D Euler simulations coupled to a breeder genetic algorithm // Proc. of 20th IAHR Symp. on Hydraulic Machinery and Systems. 2000.
[9] Mazzouji F., Couston M., Ferrando L. et al. Multicriteria optimization: viscous fluid analysis — mechanical analysis // Proc. of 22nd IAHR Symp. on Hydraulic Machinery and Systems. 2004.
[10] Mazzouji F., Francois M., Tomas L. et al. Refinements in Francis turbine design // Hydropower & Dams, Issue One. 2004. P. 53-58.
[11] Lipey A., Poloni C. Design of Kaplan runner multiobjective genetic algoritm optimization // J. of Hydraulic Research. 2000. Vol. 38, N 1. P. 73-79.
[12] Enomoto Y., Kurosawa S., Suzuki T. Design optimization of Francis turbine runner using multi-objective genetic algorithm // Proc. of 22nd IAHR Symp. on Hydraulic Machinery and Systems. 2000.
[13] Fonseca C.M., Fleming P.J. Genetic algorithm for multiobjective optimization: formulation, discussion and generalization // Proc. 5th Intern. Conf. on Genetic Algorithms. 1993. P. 416-423.
[14] Cherny S.G., Chirkov D.V., Lapin V.N. et al. Numerical simulation of a turbulent flow in Francis hydroturbine // Russ. J. Numer. Anal. Math. Modeling. 2006. Vol. 21, N 5. P. 425-446.
Поступила в редакцию 3 сентября 2008 г., в переработанном виде —10 ноября 2008 г.
