CC BY
89
Вычислительные технологии
Том 15, № 4, 2010
Оптимизационное проектирование формы проточной части гидротурбины и анализ течения в ней*
Д. В. Банников Новосибирский государственный университет, Россия e-mail: denis .baimikov@gmail. com
С. Г. Черный, Д. В. Чирков e-mail: [email protected]
Институт вычислительных технологий СО РАН, Новосибирск, Россия
В. А. Скороспелое, П. А. Турук Институт математики СО РАН, Новосибирск, Россия
Представлен обзор современных методов проектирования лопастных систем гидротурбин. Проводится обобщение предложенной ранее системы оптимизационного проектирования включением в нее методики параметризации распределения толщины лопасти. Рассматривается вопрос получения оптимальных характеристик на предварительно заданном режиме работы станции. Решена задача выбора формы лопасти рабочего колеса, дающей высокие энергетические характеристики на оптимальном режиме работы и малые пульсации давления, вызванные вихревым жгутом, на режиме неполной загрузки. Выполнены расчеты прецессии вихревого жгута в исходной и оптимальной геометриях.
Ключевые слова: оптимизационное проектирование, рабочее колесо гидротурбины, параметризация поверхности, невязкая несжимаемая жидкость, прецессия вихревого жгута.
Введение
В настоящей работе получено дальнейшее развитие методов оптимизационного проектирования проточных частей гидротурбин (ГТ) [1, 2]. Расширен вектор параметров, определяющих форму проточной части ГТ. Впервые варьируются максимальная толщина лопасти рабочего колеса (РК) и линия ее расположения. Модифицированный комплекс автоматизированного проектирования применяется для совершенствования проточной части РК строящейся Гоцатлинекой ГЭС на реке Аварское Койсу в Дагестане. С целью увеличения сработки напора на РК и выравнивания профиля скорости на входе в отсасывающую трубу (ОТ) проводится автоматизированный выбор формы проточного тракта РК, обеспечивающей минимумы интегральных кинетических энергий в сечениях, расположенных сразу за РК и во входном сечении в ОТ. При этом с фронта Парето выбирается геометрия, гарантирующая некоторое увеличение расходной составляющей скорости в центре ОТ, что при уменьшении расхода ниже оптимального препятствует развитию прецессии приосевого вихревого жгута, сходящего с втулки РК.
*Работа выполнена при финансовой поддержке РФФИ (грант № 08-01-00364).
Тем самым было снижено динамическое воздействие прецессирующего вихревого жгута на режиме неполной загрузки турбины.
Выбранные с фронта Парето формы анализируются на возможность обеспечения устойчивой работы ГТ в режиме неполной загрузки. Для этого проводится моделирование нестационарного потока в ГТ и оцениваются динамические характеристики прецессирующего вихревого жгута.
1. Прямые задачи гидродинамики турбин
Целью прямых задач расчета гидродинамики турбин является определение пространственного течения в заданной проточной части ГТ и гидродинамических характеристик потока. Для расчета течений в неподвижных элементах проточного тракта используется абсолютная декартова система координат (СК), В рабочем колесе, вращающемся с постоянной угловой скоростью ш, рассматривается относительное течение в системе координат, вращающейся вместе с РК,
1.1. Уравнения Эйлера движения несжимаемой жидкости
Для моделирования пространственного течения в абсолютной СК (х1; х2, х3) применяются уравнения Эйлера несжимаемой жидкости
дху
др
дЬ дху дхг
0, (1)
= 9г,г =1, 2, 3, (2)
где ,ш1, и>2, и>3 — декартовы компоненты вектор а скорости, р — давление, деленное на плотность жидкости, 91, 92, 93 — компоненты вектора силы тяжести.
Во вращающейся с угловой скоростью ш СК связь между вектором абсолютной скорости ^ ^ ^^^^^^^^^ ^^^^еттельной скорости и = (п1,п2,п3) задается как
w = и + ш х г, (3)
где г — радиус-вектор точки.
Уравнения (1), (2) переписываются в относительной СК в виде
дху
(4)
дпг дщпз др д^ дху дхг
¡г, г =1, 2, 3, (5)
где ¡1, ¡2, ¡3 = — ш х ш х г — 2ш х и, Здесь ш х ш х г центростремительное ускорение, 2ш х и кориолисово ускорение.
1.2. Полная нестационарная постановка
Наиболее общей является постановка прямой задачи гидродинамики турбин, в которой моделирование нестационарного течения проводится во всем проточном тракте ГТ. Нестационарная постановка позволяет моделировать весь диапазон режимов работы ГТ, в том числе и режимы неполной загрузки, учитывать взаимодействие ротора и статора, описывать пульсации сил и моментов на лопатках, моделировать прецессирующий вихревой жгут за РК и т, д.
Для расчета потока строится многосвязная блочно-етруктурированная сетка, покрывающая весь проточный тракт ГТ, состоящий из спиральной камеры, статора, направляющего аппарата (НА), РК и отсасывающей трубы. На каждой итерации расчет проводится во всех блоках с последующим обменом данных между соседними блоками. Итерации повторяются до тех пор, пока не будет найдено решение на текущем временном слое, затем происходит переход на следующий шаг по времени. Для реализации полной нестационарной постановки необходим значительный объем вычислительных ресурсов. Так, расчет периодически нестационарного течения с прецессирующим вихревым жгутом на одном периоде (около трех оборотов РК) требует нескольких дней работы процессора Соге2Био 2,6 ГГц,
1.3. Полная стационарная постановка (приближение замороженного колеса)
В данной постановке фиксируется положение РК относительно лопаток направляющего аппарата и отыскивается стационарное решение во всем проточном тракте. При этом вращение РК учитывается в уравнениях количества движения подстановкой угловой скорости ш в (5) в расчете течения в межлопастных каналах РК,
Выполнив вычисления для нескольких углов поворота РК относительно статора, можно получить оценку высокочастотных динамических нагрузок на элементы ГТ, Расчеты квазинестационарного течения в приближении замороженного колеса достаточно экономичны. Например, моделирование на процессоре Соге2Био 2,6 ГГц течения во всем проточном тракте на сетке с 0,8 млн узлов требует 2,5 ч,
1.4. Циклическая постановка
Расчеты потока в ГТ на режимах, близких к оптимальным, показывают, что поток на выходе из статора имеет слабую окружную неравномерность, не оказывающую заметного влияния вниз по потоку, а течение за РК близко к осесимметричному. При таких условиях может быть рассмотрена циклическая постановка задачи, в которой предполагается, что течения во всех межлопаточных каналах направляющего аппарата и рабочего колеса имеют осевую цикличность, т, е, повторяются от канала к каналу, В этом случае расчет проводится только в одном межлопаточном и одном межлопастном каналах.
Расчет течений в НА и в РК в циклической постановке на сетке с 60 000 узлов занимает около 5 мин процессорного времени, что позволяет использовать данную постановку для построения системы автоматизированного оптимизационного проектирования форм направляющего аппарата и РК,
1.5. Краевые и начальные условия
Рассматриваемая в настоящей работе расчетная область ГТ состоит из направляющего аппарата, РК и конуса отсасывающей трубы (диффузора). На входе в направляющий аппарат задаются угол входа $сп равномерного потока и расход жидкости Q. В выходном сечении диффузора задается давление, определяемое из условия радиального равновесия
= (в)
дг г
где г — расстояние от оси РК до точки сечения, Си — окружная компонента абсолютной скорости в этой точке. Условие радиального равновесия (6) является следствием радиальной проекции уравнения количества движения при предположении, что течение в диффузоре близко к осесимметричпому и радиальная компонента скорости Сг мала. Распределение давления в выходном сечении получается путем интегрирования уравнения (6) вдоль радиуса от стенки диффузора (где задается постоянное давление рои^ к его оси. На твердых стенках ставится условие непротекания. Если в расчете есть границы цикличности, то на них ставится условие периодичности течения. Угловая скорость ш
Решения, получаемые в стационарных постановках, не зависят от выбранных начальных данных. Рассматриваемая в работе динамика вихревого жгута имеет периодическую нестационарность, которая также не зависит от начальных данных. Однако процесс достижения периодического решения для грубого начального приближения может быть очень долгим, поэтому начальное поле для нестационарного расчета было получено при расчете потока в приближении замороженного колеса.
2. Метод решения прямой задачи 2.1. Концепция искусственной сжимаемости
Численный алгоритм решения уравнений (4), (5) основан на методе искусственной сжимаемости, который состоит во введении в уравнение неразрывности производной по псевдовремени от давления, а в уравнения количества движения — производных по псевдовремени от соответствующих компонент скорости. Модифицированная система принимает вид
др дпз
Вт-^г (7)
дпг дпг дпгщ др
дт дЬ дхз дхг
¡г, г = 1, 2, 3, (8)
где в ~ коэффициент искусственной сжимаемости. Система уравнений (7), (8) решается численно. При этом на каждом шаге по физическому времени Ь проводится уста-
т
решения методом установления как стационарных, так и нестационарных задач. Если задача нестационарна, то маршевый ход по физическому времени проводится до тех пор, пока не будет получена устойчивая периодичность, для стационарной задачи — пока не занулитея производная по физическому времени.
2.2. Неявный метод конечных объемов
Для построения консервативной разностной схемы применяется неявный метод конечных объемов. При этом уравнения (7), (8) записываются в виде интегральных законов сохранения
(тс-^ + К4^) I + у К • ¿Б = I тау, (9)
^ ' V дУ V
где V — объемы элементарных ячеек, на которые разбивается расчетная область,
К* = ^(0,1,1,1), Кт = ^(1,1,1,1),
д
(V \
Щ
щ \ из )
К
( вщ ви2 и2 + р и\и2 и1и2 и1из
виз
и1из и'2 + V и2из
Е
0 /1 /2 V /з /
и2из и3 + р )
Дискретизация уравнения (9) для каждой ячейки Vijk дает выражение
^(<3"+')'^-(СГ+'Г + + „ = (кнз„+1),+1, (10)
где Ат, А£ — шаги по псевдовремени и физическому времени соответственно, 5 — номер итерации по псевдовремени, п — номер слоя по времени. Дискретизация (10) имеет
т
по физическому времени ¿. Правая часть имеет структуру
ЯШ = - (К ■ ^+1/2 - (К ■ S)i—1/2 + (К ■ 8),+1/2 - (К ■ 8),-—
(К ■ 8)к+1/2 - (К ■ 8)к—1/Н + FV,
где (К ■ 8)т+1/2 — разностный поток через общую г рань ячеек т и т+1, вектор 8 имеет величину площади грани т + 1/2 и направлен то нормали к ней, наружу ячейки т. Вычисление потоков проводится по формуле
(К • 8)т+1/2 = - (К(Цт) + 1ОДт+1))8т+1/2 -
А
т+1/2
Ат+1/2д
- ТО".
т+1/2,
(Н)
где Ат+1/2д = дт+1 - (Зт, W — добавочный член, предложенный в работе [3], повы-
А
= А = А+ + А", |А| = А+-А". (12)
А+ А—
ченпя соответственно.
2.3. Приближенная Ш-факторизация
Для решения системы нелинейных уравнений (10) последние линеаризуются с использованием метода Ньютона
\Дт
V - (^гЫК 2Д£ ) 1<9д ,
(д^1)^1 — (д^1)
га+1\ я
-К'
з(дга+1)5 - 4сг + дга_1
2Аг
V + (КШ^1)5.
(13)
При построении неявного оператора в левой части (13) предполагается, что на всех гранях поток аппроксимируется с первым порядком. Применяя формулу (11) при вы-д
чпсленпн ——Н.11Н с \¥ = 0, получим следующую систему линейных уравнений:
дд
(Д7КТ + 2Д*К') У + К-1/2^-1/2 + А+_ 1/2А,-1/2+
+А-+1/2Д3+1/2 + А+_1/2Д3-1/2 + +А-+1/2Дй+1/2 + А+_1/2 Дк_1^ Д*+1д
(14)
где Д^д = (д^1)^1 — (д^1)5. Перепишем (14) в виде
В + Аг+1/2Тг+ А+-1/2Тг + А_+1/2Тз+ А+_1/2Т3 +
+Ак+1/2Тк А+_1/2Тк
Д5+1д = + V + (кн8га+1г,
(15)
где Т± — оператор сдвига па один узел вперед (+) или назад (—) по индексу т,
В = (Д^^ + V + ~ Аг+1/2 + А^-1/2 - ~А/+1/2 + Ак-1/2 ~ Ай+1/2"
Теперь оператор в левой части (15) может быть приближенно факторизован:
В ^-1/2^ А+ —1/2Tг А+"_1/2ТА:
+Ak_+l/2Tfc+] Д*+1д = — к
и последовательно обращен:
В_1
я Л глп I г\п_1
В + А +
3(дп+1)* — 4дп + дп
Д*+1/2дгз к = В _1
К
2ДЬ
з(дга+1)5 - 4дга + д™-1 2Аг
-V + (КШ^1)5
(16)
V + (КНБ^1)*
+А+1/2Д*+1/2дг_ 13 к + А+_1/2Д*+1/2дгЗ _1 к + А+_1/2 Д*+1/2дгз к _1
Д*+1дгзк = Д*+1/2дгзк — в _1 [А"+1/2Д*+1дг+1зк+
(17)
(18)
+А_+1/2Д + дг3+1 к + А _+1/2Д + дг3 к+1
На первом шаге по формулам бегущего счета (17) совершается разовый обход области в направлении возрастании всех индексов и определяется вспомогательная величина Д5+1/2дгзк; на втором шаге (18) в направлении убывания индексов определяется величина Д*+1дгз к, по которой находится вектор неизвестных па (в + 1)-итерации:
(дз 1)*+1 = (дз 1)5+Д*+1дгз к.
Итерации по в повторяются до сходимости (дп+1)5 ^ дга+1.
3. Обратные задачи гидродинамики турбин — задачи проектирования
Обратные задачи гидродинамики турбин заключаются в нахождении формы проточной части Г Г. удовлетворяющей заданным критериям на заданных режимах работы станции. Критерии включают в себя требования к энергетическим показателям, режимным параметрам, прочностным, технологическим, кавитационным и другим характеристикам. Обратные задачи гидродинамики турбин называют в литературе задачами проектирования. Методы решения задач проектирования можно разделить на два класса: оптимизационного и прямого проектирования.
3.1. Методы оптимизационного проектирования
Обзор работ, в которых применяются методы оптимизационного проектирования проточных частей ГТ, приведен в [1, 4]. В них варьируется форма проточной части. Для каждой вариации решается прямая задача гидродинамики и выбирается форма-вариация, максимально удовлетворяющая заданным критериям.
Критерии формулируются в виде целевых функционалов, которые могут характеризовать локальные свойства потока, такие как распределение скоростей или давления в заданном сечении, либо интегральные характеристики — КПД, мощность, стоимость установки, запас кавитации и прочности. Если целевой функционал один, то решается задача одноцелевой оптимизации, и экстремальное значение функционала соответствует геометрии, неулучшаемой по данному критерию. Обычно требуется улучшить сразу несколько критериев по разным качественным характеристикам. В этом случае задача сводится к поиску экстремума одного взвешенного функционала либо решается задача многоцелевой оптимизации и строится оптимальный фронт Парето.
Вариация формы проточной части осуществляется посредством изменения параметров, определяющих ее геометрию.
В зависимости от числа варьируемых параметров и целевых функционалов выбирается стратегия поиска оптимальных решений.
Классические градиентные методы спуска и анализа чувствительности [5, 6] подходят для решения задач с непрерывным функционалом и малым числом варьируемых параметров. Если информации о гладкости решения нет, то применяются детерминированные прямые методы Хука—Джнвса или Нелдера—Мида [4, 5, 7]. Методы, использующие функции отклика [4, 8], позволяют аппроксимировать исходную модель более простой полиномиальной, для которой затем решается задача оптимизации. Такие модели обычно строятся путем планирования численного эксперимента либо обучением нейронных сетей. Удачно построенная аппрокеимационная модель позволяет значительно сократить количество анализируемых геометрий для нахождения оптимального решения. Существенным ограничением данного метода является требование гладкости у исходных функциональных зависимостей.
В настоящее время большую популярность имеют стохастические оптимизационные алгоритмы, такие как метод моделируемого отжига [7] и эволюционные алгоритмы [1, 5, 7], не требующие гладкости от целевых функций, обладающие способностью не останавливаться в локальных экстремумах и применяемые для решения многокритериальных задач. Необходимость большого числа вычисления целевых функций в стохастических методах компенсируется возможностью их параллельного расчета.
Метод моделируемого отжига определяет экстремум функционала путем случайных вариаций текущего вектора геометрических параметров. При этом худшая вариация может приниматься как основная с некоторой вероятностью, уменьшающейся в процессе итераций алгоритма. Чем меньше скорость убывания этой вероятности, тем выше возможность найти глобальный оптимум, но и больше время работы алгоритма. Эволюционные алгоритмы осуществляют поиск оптимальных решений, одновременно анализируя множество текущих геометрий, называемых популяциями, которые эволюционируют на протяжении многих поколений в соответствии с предписанными правилами, Наиболее известные эволюционные алгоритмы — генетический оптимизационный алгоритм и метод колонии частиц.
Формализация схемы оптимизационного проектирования в виде математической задачи поиска экстремумов функционалов по параметрам задачи позволяет автоматизировать ее решение на ЭВМ, При этом создаются комплексы программ, сочетающие в себе блоки оптимизации, алгоритмы восстановления параметризованной поверхности проточного тракта и построения расчетных сеток, гидродинамических солверов, пре- и постпроцессоров для подготовки и обработки расчетных данных.
Чтобы оптимальная геометрия имела рассчитанные характеристики, необходимо адекватное моделирование процессов, происходящих в ГТ, Непосредственное решение трехмерных уравнений Навье—Стокса во всей проточной части требует значительных затрат времени, особенно при решении нелинейных многокритериальных задач с большим числом варьируемых параметров. Поэтому при построении автоматизированной системы выбирается компромисс между количеством исследуемых форм в поиске наилучшей геометрии, подробностью сеток и размером расчетной области, сложностью используемых моделей движения жидкости, точностью учета и выявления принципиально важных локальных особенностей течения и соответствующим временем, затрачиваемым на такой поиск,
3.2. Методы прямого проектирования
В методах прямого проектирования устанавливается непосредственная связь между формой проточной части ГТ и характеристиками течения в ней. Эти методы отличаются используемыми моделями течения жидкости, методиками представления геометрии, уравнениями связи форм и параметров потока, а также итерационными алгоритмами, которые находят форму, соответствующую заданным для потока условиям,
В настоящее время методы прямого проектирования ГТ [9-21] созданы только для лопастных систем (ЛС), Во всех этих работах фиксируются формы обода и ступицы и варьируется угловая координата лопасти, В [11, 12, 14] дополнительно подбирается распределение толщин профиля,
В работах [9-12] задача проектирования профиля ЛС решается с применением двумерных уравнений потенциального потока, в [13-17] решаются квазитрехмерные уравнения идеальной жидкости, в [18, 19] — трехмерные уравнения Эйлера, Для выделения действительного течения в [9-19] на выходной кромке ставится условие Кутты-Жуковского. В [21] рассмотрена двумерная вязкая модель течения жидкости, а работа [20] является одной из первых, где предложен метод прямого проектирования формы Л С с использованием трехмерных уравнений движения вязкой жидкости.
Во всех указанных работах задается один из двух критериев, которому должна удовлетворять проектируемая геометрия: первый — распределение закрутки гСи в меж-
лопастном канале, второй — распределение по поверхности профиля давления p либо нагрузки Ар, равной разности давлений между рабочей и тыльной сторонами профиля. Для построения уравнения связи между заданными критериями и формой профиля в работах [9, 10] используется метод конформных отображений, в [11, 12] — метод годографа скорости, В [13-15] для плоских решеток в слое переменной толщины на основе уравнений движения жидкости выводится соотношение, связывающее угловую координату профиля А< и закрутку потока rCu. Интегрируя полученное соотношение вдоль хорды профиля, находится его угловая координата <р. В [16] уравнение связи между вектором нормали к поверхности бесконечно тонкой лопасти и вектором локальной скорости записывается в виде дифференциального уравнения первого порядка, В работе [17] связь моделируется добавлением сил в источниковый член уравнений количества движения, В [18] впервые предложена концепция проницаемых стенок, успешно примененная в [19, 20],
Существует несколько итерационных алгоритмов нахождения формы лопасти, определяющей требуемые характеристики потока, В [13, 15, 16, 18] на первом шаге решается прямая задача течения жидкости для фиксированной геометрии ЛС, на втором шаге из условия связи модифицируется форма профиля. Шаги повторяются до тех пор, пока не перестанет изменяться геометрия, В разных работах для нахождения решения требуется от трех до нескольких десятков итераций алгоритма, В [17, 19-21] стационарное решение уравнений движения жидкости отыскивается методом установления, при этом форма лопасти корректируется после каждой итерации по псевдовремени. Такой подход позволяет получать форму профиля одновременно с решением уравнений движения жидкости, В работе [19] указывается, что время решения обратной задачи на 20% больше времени решения прямой задачи тем же методом, но при фиксированной форме лопасти. Увеличение времени связано с необходимостью перестраивать расчетную сетку после каждой модификации профиля и с возрастанием числа итераций, требуемых для сходимости, В [21] время решения обратной задачи методом установления на двумерной сетке с 4000 узлами на процессоре Pentium 2 GHz заняло 3 ч,
В работе [21] эмпирически исследуются существование и единственность решения обратной задачи для сжимаемой жидкости в двумерном случае. Получено, что используемый здесь метод установления не всегда сходится к стационарному решению, особенно при нате капни потока на профиль под большим углом атаки, В [9] были выведены три необходимых условия существования решения для двумерного несжимаемого потенциального потока, однако для трехмерных уравнений сжимаемой жидкости такие условия неизвестны [20], В работах [11, 12] указывается, что по заданным распределениям скоростей на профиле в результате решения обратной задачи может получаться геометрия лопасти с неприемлемым характером распределения толщин, В [20] также отмечается, что постановка условий в виде нагрузки на профиль может давать решения с нереалистичной геометрией лопасти,
К недостатку методов прямого проектирования относится следующее. Основные характеристики ЛС, такие как КПД, коэффициенты потерь энергии и давления, напрямую не связаны с распределениями скорости или давления, задаваемыми в качестве целевых критериев. Поэтому нерешенными остаются вопросы о том, насколько заданный профиль закрутки в межлопастном канале или нагрузки на лопасть соответствуют оптимальным энергетическим характеристикам всей установки и в какой мере возможны улучшение ее КПД или снижение потерь в проточной части. Часто требуемые
распределения скоростей или давлений априори неизвестны. Кроме того, в данном случае нет гарантии, что для заданных распределений существуют решения.
Несомненным достоинством указанного подхода, позволяющим применять его в ежедневной инженерной практике, является малое время решения задачи проектирования — существенно меньшее затрачиваемого при оптимизационном проектировании.
4. Метод оптимизационного проектирования проточного тракта рабочего колеса радиально-осевой гидротурбины
Для решения задачи оптимизационного проектирования РК авторами разработана параметризация его формы, содержащая 28 параметров. Варьируются формы обода и ступицы, входной и выходной кромок лопасти, распределение толщин лопасти и угловая координата ее срединной поверхности. Сформулированы целевые критерии, характеризующие энергетические и кавитационные качества лопастной системы, определяемые на основе решений прямой задачи в приближении пространственных уравнений Эйлера. Приведена математическая постановка задачи многоцелевой оптимизации с фазовыми, геометрическими и гидродинамическими ограничениями, для решения которой используется многокритериальный генетический алгоритм (МКГА). Представлена общая схема решения задачи оптимизационного проектирования.
4.1. Параметризация формы рабочего колеса
В работах [1, 2] была предложена методика построения параметрического семейства поверхностей, задающих геометрию РК гидротурбины для решения задачи выбора оптимальной формы рабочего колеса. В этих публикациях поверхность лопасти (рис. 1) представляется в цилиндрической системе координат в виде
где г(и, у) = (Я(и,у(и, у), Ф(и, у)) — срединная поверхность, п — орт нормали к срединной поверхности, ^(и, у) — функция распределения толщины.
И(и, у) = г(и, у) + й(и, у) • п, и, у Е [0,1]
(19)
Ступиц;
Выходная кромка "
а
Линия максимальной толщины
Вариация формы РК проводится за счет изменения угловой координаты срединной поверхности лопасти Ф(п, у) и меридиональной проекции (К(п,у), Z(и, у)) линий обода, ступицы, входной и выходной кромок. Для параметризации срединной поверхности лопасти г (и, у) используется подход с относительным представлением ее угловой координаты Ф(и, у) в виде бикубического полинома с 16 свободными параметрами (<р1,..., ^16), Вариация меридиональной проекции РК осуществляется восемью свободными параметрами, изменяющими форму обода (р1,р2), ступицы (р3,р4), входной (р5,р6) и выходной (р7,р8) кромок лопасти,
В настоящей работе проводится развитие предложенной ранее методики [1, 2] добавлением параметризации функции распределения толщины лопасти ¿(и, у). Параметры ^ и ¿2 характеризуют положение линии максимальной толщины на ободе и ступице соответственно, Параметры ¿3 и ¿4 показывают, во сколько раз изменяется максимальная толщина лопасти на ободе и ступице по сравнению с прототипом. Примеры вариации толщины профиля лопасти приведены на рис, 2,
Вектор варьируемых параметров, задающих форму РК, имеет вид
4.2. Формулировка целевых функционалов и ограничений задачи оптимизационного проектирования
В [13, 15] сформулированы основные требования к рабочему колесу. При проектировании формы РК необходимо обеспечить следующие характеристики:
— сохранение заданных параметров оптимального режима, т, е, величин расхода, угловой скорости вращения РК и напора;
— максимально возможный КПД на оптимальном режиме;
— хорошие кавитационные показатели;
— высокую надежность и дол го временность работы гидромашины.
Для выполнения последнего требования необходимо, чтобы статические и динамические нагрузки не превышали заданных предельных значений.
Для нахождения геометрии РК, удовлетворяющей этим требованиям, ниже формулируются целевые функционалы и ограничения оптимизационной задачи,
4-2.1. Сохранение заданных режимных параметров
Режим работы ГТ характеризуется расходом воды Q, угловой скоростью вращения РК разностью энергий потоков Н между входом в спиральную камеру и выходом из отсасывающей трубы, называемой напором станции. Расход Q связан с распределением скоростей во входном сечении расчетной области 50 соотношением
и сохраняется в процессе решения задачи в силу задания этих скоростей в 50 в качестве краевого условия и выполнения уравнения неразрывности (4),
Угловая скорость вращения РК входит в выражения для центробежных и корио-лиеовых сил в источниковых членах уравнений (5) и поэтому в силу явного задания также сохраняется в процессе решения задачи оптимизационного проектирования.
X = (Х1,...,ХМ) = (^1,...,^16,Р1,...,Р8,^1,...,4).
(20)
(21)
Напор станции Н имеет структуру
Н = Нд + Н, (22)
где Нд — полезная энергия станции, вырабатываемая в виде электроэнергии, Н — гидродинамические потери энергии, обусловленные вязкостью воды и представляемые как сумма потерь в каждом элементе ГТ:
Н = ДНСП + ДНСт + ДНна + ДНрк + ДНот. (23)
Н
обще говоря, может отклоняться от заданного режимного значения. Поэтому его постоянство должно формулироваться в виде ограничения оптимизационной задачи, В процессе оптимизационного проектирования рассчитывается невязкое течение в одном РК, По насчитанному полю Нд определяется как разность полных энергий
(24)
Н
Н
Н
ваетея ограничение на напор Нд в виде неравенства
|Нд - 0.95Н|< е, (25)
где е — малый параметр, принимаемый обычно как е = 0.005,
4-2.2. Энергетические характеристики потока в рабочем колесе
КПД станции п определяется в виде
П = (26)
Н
геометрию, дающую максимально возможный КПД, можно найти, приняв в качестве
Н
задачу ее минимизации. Однако в настоящей работе для гидродинамических расчетов применяется модель невязкой жидкости, не позволяющая прямо определять потери энергии. Поэтому вслед за работами [1, 2, 4, 9-21] здесь используются критерии, формулируемые на основе кинематики невязкого потока, косвенно позволяющие уменьшать потери энергии, К таким критериям относятся отсутствие резких градиентов давлений и скоростей, монотонность профиля распределения давления на поверхности лопасти, равномерность поля скорости и малая закрутка потока на выходе из ЛС, В частности, минимизация кинетической энергии потока в сечении $1, расположенном сразу за ЛС рабочего колеса (рис, 3),
1 Г М2
^ = (27)
0 0.2 0.4 0.6 Я
Рис. 3. Меридиональное сечение расчетной области
при фиксированном значении энергии перед ЛС приводит к выравниванию профиля меридиональной скорости и, как следствие, к уменьшению потерь энергии в РК, То же самое относится к сечению S2 на входе в отсасывающую трубу. Минимизация кинетической энергии здесь приводит к выравниванию профиля скорости в52ик уменьшению гидравлических потерь энергии в ОТ,
Минимальное значение кинетической энергии при фиксированном расходе Q имеет равномерный поток, скорость которого V = Q/S^.
1 г V2 V2 Q2
= = (28)
В силу сказанного в качестве целевых функционалов рассматриваются выражения
Я = (29)
Ш1П Ькг
где индекс г относит их к сечению Si,
4-2.3. Кавитационные характеристики потока
Явление кавитации возникает в том месте, где локальное давление становится меньше давления насыщенного пара
Р < Рса^ (30)
Дня обеспечения требуемых кавитациоппых качеств РК ставится ограничение па относительную площадь кавитации в виде
Д:ау = М [ 5слЧ8<5- 100%, (31)
»-'вис / •эвис
где — площадь тыльной стороны лопасти, 6са¥ = 1 при выполнении условия (30) и 6са¥ = 0 — иначе, 6 — допустимый уровень кавитации, обычно составляющий 15% всей площади тыльной стороны .попасти.
4.3. Математическая постановка задачи оптимизационного проектирования
Задача оптимизационного проектирования состоит в нахождении вектора х в (20), обеспечивающего минимальные значения M функционалов
min Fi(x),..., min FM (x) (32)
при фазовых
XL,i < Xi < xRyi, i = 1,..., N, (33)
геометрических
^i(x) < 0, i = 1,...,/, (34)
и гидродинамических
(x) < 0, j = 1,...,J, (35)
ограничениях. Решением задачи (32)-(35) является семейство точек х, называемое оптимальным фронтом Парето,
4.4. Многокритериальный генетический алгоритм решения задачи оптимизационного проектирования
Для поиска геометрий, определяющих минимальные значения функционалов, используется многокритериальный генетический алгоритм (МКГА), Подробно данный алгоритм для случая многоцелевой оптимизации представлен в монографии [1], где также проведено исследование влияния параметров МКГА на эффективность работы алгоритма.
Решение задачи оптимизационного проектирования состоит из следующих шагов,
1, На основе общих статистических данных и имеющихся результатов экспериментального исследования ГТ подбирается геометрия проточной части, наиболее удовлетворяющая требованиям к проектируемой. Устанавливаются значения ^0пт
И ^ОПТ 5 ^
пределах которых должны быть получены заданные характеристики. Выбранный прототип соответствует начальному набору геометрических параметров хо.
2, Формируется начальная, соответствующая в = 0 и р = р0, популяция
хъ xp, (36)
р
метров (20), определяющий геометрию РК, Начальная популяция создается случайным образом так, чтобы каждый индивидуум удовлетворял ограничениям (33), (34), Прототип х0, построенный на первом шаге решения задачи, включен в популяцию (36),
3, Для каждой формы проточной части, задаваемой вектором х*, строится расчетная сетка, на ней проводится гидродинамический расчет течения для заданного режима, по которому вычисляются значения функционалов
Р* = (Я(х*),..,*М (х*)). (37)
Течения в различных формах рассчитываются одновременно на многопроцессорных системах, при этом получаемое ускорение счета близко к линейному. По значениям функционалов и ограничений вычисляется критерий качества индивидуума х*
С(^1(х*),...,^м(х*)) : Дм ^ Д. (38)
4, Производится отбор Тг •р индивидуумов, имеющих наименьшие значения функционала (38), по которым с помощью операторов ГА строится новая популяция х!+1,..., хр+1
р
5, Производится переход на третий шаг с переприсвоением в = в + 1 до тех пор, пока не будет рассчитано поколений, обеспечивающих нахождение оптимального фронта Парето,
5. Результаты решения задачи оптимизационного проектирования формы рабочего колеса и выбор оптимальной геометрии с фронта Парето
Изложенная выше методика применена дня оптимизационного проектирования РК турбины Гоцатлипской ГЭС, Сначала рассчитывается течение в прототипе РК, выявляются его недостатки. Затем ставится оптимизационная задача, направленная па их устранение. Анализируется полученный фронт Парето и выбирается геометрия РК, наиболее полно соответствующая поставленным критериям.
5.1. Геометрия, расчетные сетки и режимные параметры
Проточный тракт гидротурбины включает в себя 20 лопаток направляющего аппарата и 13 лопастей рабочего колоса. Сетка покрывает один межлонаточный канал НА, один межлонаетной канал РК и часть конуса ОТ, Общее число ячеек в НА 36 000, в РК и конусе ОТ — 49 000, На рис, 4 изображены изометрии лопаток направляющего аппарата и рабочего колеса и сеток в каналах дня расчета в циклической постановке.
Расчеты течений выполнялись для приведенной ГТ с диаметром РК D = 1 ми напором H = 1 м. Режимные параметры были следующими: расход Q = Qon-r; угловая скорость вращения ш = шопт, открытие НА составляет 75 % от а0,тах. Параметры численного алгоритма: в = 86.9, Дт = 0.1, Для получения установившегося решения в циклической постановке потребовалось выполнить 4000 итераций но нсевдовремепи, что заняло около 5 мин работы CPU Intel 2,6 ГГц,
Рис. 4. Изометрии лопаток направляющих) аппарата и рабочих) колеса и сеток в каналах для расчета в циклической постановке; вид со стороны спиральной камеры: а сверху, б снизу
5.2. Результаты оптимизационного проектирования
Выполнен поиск геометрий, определяющих минимумы значений функционалов и (29), Форма входной кромки была зафиксирована, как у прототипа, а варьировались распределение толщин лопасти и геометрия выходной кромки лопасти, так как именно эта геометрия в значительной мере определяет профиль скорости за РК,
Для получения геометрии с кавитационными характеристиками не хуже, чем у прототипа, задавалось ограничение на размер области кавитации (31) 6 = 0.15, Для сохранения режимных параметров ставилось ограничение на напор (25),
На рис, 5 в плоскости значений функционалов — ^2 представлен фронт Парето, состоящий из 31 геометрии, взятый после 25 поколения 60 геометрий-индивидуумов. Видно, что функционалы ^ и ^2 являются конфликтующими, так как уменьшение одного из них приводит к росту другого.
Для анализа были выбраны три геометрии, соответствующие средней (Ор^16) и крайним (Ор^1 и Ор^31) точкам фронта Парето, В табл. 1 приведены интегральные характеристики выбранных геометрий и прототипа, В целях удобства сравнения там же указаны наименьшие возможные значения функционалов , ^2 для гипотетической геометрии х*. Дополнительно были проведены расчеты в приближении уравнений Рей-нольдса [1] в выбранных геометриях, которые позволили оценить вязкие потери энергии Д^рк- Из таблицы видно, что минимизация функционала приводит к уменьшению
ДЛРК ■
1 ...........................
1.6 1.7 1.8 1.9 2 г
1 1
Рис. 5. Фронт Парето, полученный при двухцелевой оптимизации: 1 — прототип, 2 — Ор1>1, 3 - ОрМб, 4 - Ор^31
Таблица 1. Значения интегральных величин в различных геометриях РК
Геометрия А-сау> % Д^рк, %
x* 1.00 1.00 _ _
ОрМ 1.66 1.47 0 2.2
ОрМб 1.72 1.14 9 2.3
Ор^31 1.92 1.08 12 2.5
Прототип 1.98 1.08 17 2.6
С-, м/с
1.5
0.5
г
0.4 0.2
-0.2
0.4 1-, ■ ■ , ■
0.1 0.2 0.3 0.4
0.1 0.2 0.3 0.4 Л9м
Рис. 6. Распределение осевой С (а) и окру ж ной Си (б) компонент скорости в сечении Б2 для различных геометрий с фронта Парето: 1 — прототип, 2 — Ор^!, 3 — Ор^16, 4 — Ор^З!
Рис. 7. Профили лопасти на ступице (а), средней поверхности (б) и ободе (в) у прототипа (сплошная линия) и геометрии Ор^16 (штрих)
На рис. 6 представлены профили скоростей Сх и Си в сечении Б2 для выбранных геометрий. Геометрия Ор!-1 имеет наименьшее значение но существенно неравномерный профиль скорости на входе в отсасывающую трубу, что приведет к увеличению потерь в ней. Геометрия Ор1-31 наиболее близка к начальному РК по распределениям скоростей и в пространстве целевых функционалов (см. рис. 5). Геометрия Ор1-16
Си
оптимизированной по функционалам и ¥2. На рис. 7 приведено сравнение профилей лопасти РК прототипа и Ор1-16 на трех поверхностях: ободе, ступице и средней поверхности между ними.
В результате решения оптимизационной задачи не удалось найти геометрию с еще меньшим значением функционала чем у исходной геометрии. Это объясняется тем, что прототип имел практически однородное распределение скорости по компоненте Сх (см. рис. 6). У всех геометрий, лежащих на фронте Парето, ограничение по кавитации выполнено, причем если двигаться вдоль фронта Парето от точки Ор1-1 к Ор1-31, то область кавитации увеличивается: у геометрии Ор!-1 6 = 0, Ор!-16 6 = 0.09, Ор^З! 6 = 0.12.
6. Анализ выбранной оптимальной формы
Большое внимание при оптимизационном проектировании ГТ на ее оптимальном режиме работы уделяется тому, как спроектированная ГТ поведет себя на других режимах работы. В частности, на режимах частичной загрузки за РК возникает нестационарный прецессирующий вихревой жгут (рис. 8), обусловливающий пульсации давления. Если размах колебаний больше радиуса диффузора, то возникают сильные биения по стенкам диффузора, которые приводят к снижению надежности и долговечности конструкции. Наличие вихревого жгута также уменьшает площадь поперечного сечения на входе в ОТ, что приводит к возрастанию расходной компоненты скорости и, как следствие, к увеличению вязких гидравлических потерь энергии. Поэтому на стадии проектирования необходимо оценивать характеристики вихря, добиваясь по возможности меньших его интенсивности и размаха колебаний.
Возрастание осевой компоненты скорости к оси РК у выбранной геометрии Ор^16 может способствовать улучшению стабильности работы ГТ на режимах частичной загрузки. Для проверки этой гипотезы были проведены нестационарные расчеты в полной постановке для начальной и оптимальной геометрий на режиме неполной загрузки. Найдены интенсивности, радиус и размах колебаний осевого вихря, проведен амплитудно-частотный анализ.
6.1. Геометрия, расчетные сетки и режимные параметры
Расчетная область включала все межлопаточные каналы НА и РК и диффузор отсасывающей трубы. Общее число ячеек в НА 674 000, в РК 672 ООО и в конусе ОТ 323 000. По сравнению с расчетной областью (см. раздел 5.1) конус диффузора был удлинен на три диаметра РК с^ = 0.8 до 4.0.
Режимные параметры были следующими: расход Q = 0.75^опт, угловая скорость ш = шопт, открытие НА составляло 60% от а0,тах. Шаг по псевдовремени Ат = 0.1, параметр псевдосжимаемости в = 15.9.
6.2. Расчет течения в постановке замороженного колеса
Расчет течения в приближении замороженного колеса не позволяет выйти на стационарное решение. За РК формируется вихревой жгут, положение которого в процессе итераций по псевдовремени изменяется.
Рис. 8. Вихревой жгут за рабочим колесом
Си,м/с 2г}
-0.4 -0.2
х,М
-0.4 -0.2 0 0.2 0.4 х,м
Рис. 9. Картина течения в сечении Б2 в прототипе: распределение давления и линии тока (а) и профиль окружной скорости Си при у = 0 (б); постановка замороженного колеса
Таблица 2. Амплитуда, радиус и интенсивность вихревого жгута в диффузоре, полученные в расчете в постановке замороженного колеса
Сечение Геометрия Rb n г
Z = 0.7 Прототип 0.118 0.106 0.065
Opt-16 0.094 0.078 0.042
Z = 1.0 Прототип 0.210 0.103 0.065
Opt-16 0.152 0.078 0.048
Z = 1.2 Прототип 0.253 0.116 0.065
Opt-16 0.210 0.082 0.049
На рис. 9 приведена картина течения в сечении S2\ а — распределение давления и линии тока на текущей итерации, б — распределение окружной компоненты скорости Cu вдоль линии y = 0, проходящей через центры диффузора и вихря. Хорошо виден участок графика, где Cu линейно зависит от x. На этом участке жидкость в цилиндре радиуса rb вращается как твердое тело. Диаметр вихря равен максимальному расстоянию, на котором сохраняется линейная зависимость Cu от x. Амплитуда вихря Rb определяется расстоянием от центра диффузора до центра вихря.
В табл. 2 приведены параметры вихревого жгута, найденные в разных сечениях Z — const. Интенсивность вихревого жгута определяется по формуле
Г = rb Cu. (39)
Из таблицы видно, что rb и Г слабо зависят от Z, а радиус Rb растет при удалении от РК.
6.3. Расчет течения в нестационарной постановке
Чтобы оценить динамику вихревого жгута, было проведено моделирование потока в нестационарной постановке. Использовалась та же сетка, что и в постановке замороженного колеса. Так как течение в направляющем аппарате не изменяется со временем, то для сокращения времени расчета он был исключен из расчетной области. На входе
Таблица 3. Амплитуда, радиус и интенсивность вихревого жгута в диффузоре, полученные в расчете в нестационарной постановке
Сечение Геометрия Rb п Г
Z = 0.7 Прототип 0.159 0.141 0.114
Opt-16 0.124 0.108 0.082
Z = 1.0 Прототип 0.221 0.140 0.114
Opt-16 0.182 0.108 0.088
Z = 1.2 Прототип 0.257 0.182 0.151
Opt-16 0.220 0.112 0.089
в РК задавалось поле скорости, полученное при расчете в постановке замороженного колеса. За один шаг по времени РК поворачивается на 3.08 град, а полный оборот делает за 117 шагов. Подытерации выполняются до достижения невязки Err = 10-3. Время расчета одного оборота РК на шести процессорах Xeon 2.66 ГГц занимало 12 ч. Было рассчитано 30 оборотов РК. Полученная картина течения близка к течению, найденному в постановке замороженного колеса.
В табл. 3 представлены параметры вихревого жгута, рассчитанные в нестационарной постановке. Видно, что rb и Г не зависят от Z, а радиус как и в постановке замороженного колеса, растет при удалении от РК. В нестационарном расчете увеличились параметры вихря, т. е. вихрь в этом случае был более ярко выражен.
30 t, с
А/H, % 21.5-1 1
0.5 -I
0-1 ю]
л/
10"
Ю1 ///;.
Рис. 10. Пульсации давления в сечении конуса Z = 1.2 (слева) и их спектр (справа) для геометрии РК прототипа (а) и Ор^16 (б)
На рис, 10 изображены пульсации давления на стенке конуса ОТ в сечении Z = 1.2, где наблюдается их максимальная амплитуда, и соответствующее преобразование Фурье. Период прецессии вихревого жгута у прототипа tb = 2.95 с, основная частота 0.28 fr (fr — частота вращения РК), у оптимизирован ной геометрии tb = 2.78 с, основная частота пульсаций 0,3fr, Амплитуда пульсаций давления у исходной геометрии составляла 2.5%, у оптимизированной 1.8%. Таким образом, на режиме неполной нагрузки пульсации давления уменьшились. Второй пик на обоих графиках спектра пульсаций
fr
Таким образом, в работе проведено оптимизационное проектирование рабочего колеса гидротурбины. Задача решалась с двумя функционалами при ограничениях на площадь кавитации и напор станции. В результате получен фронт Парето, содержащий множество неулучшаемых геометрий, из которого выбирается та, что в наибольшей степени удовлетворяет исходным требованиям. Расчеты показали, что спроектированная геометрия обладает улучшенными энергетическими и кавитационными характеристиками. Проведены расчеты течения в полной стационарной и нестационарной постановках на режиме неполной загрузки турбины, и сопоставлены максимальные амплитуды и частоты пульсаций давления в конусе отсасывающей трубы.
Список литературы
[1] Черный С.Г., Чирков Д.В., Лапин В.Н. и др. Численное моделирование течений в турбомашинах. Новосибирск: Наука, 2006. 202 с.
[2] Банников Д.В., Черный С.Г., Чирков Д.В. и др. Многорежимная оптимизация формы рабочего колеса гидротурбины // Вычисл. технологии. 2009. Т. 14, № 2. С. 32-50.
[31 Chakravarthy S.R., Osher S. A New Class of High Resolution TVD Schemes for Hyperbolic Conservation Laws. AIAA Paper 85-0363. 1985.
[4] Аэродинамический расчет и оптимальное проектирование проточной части турбома-шин / А.В. Бойко, Ю.Н. Говорущенко, С.В. Ершов, А.В. Русанов, С.Д. Северин. Харьков: НТУ "XIII Г. 2002. 356 с.
[5] Eisinger R., Ruprecht A. Automatic shape optimization of hydro turbine components based on CFD // TASK Quarterly. 2001. Vol. 6. P. 101-111.
[61 Mazzouji F., Couston XL. Ferrando L. et al. Multicriteria optimization: Viscous fluid analysis — mechanical analysis // Proc. of 22nd IAHR Svmp. "On Hydraulic Machinery and Systems". Sweden, Stockholm, 2004.
[71 Rogalsky Т., Derksen R.W., Kocabiyk S. Differential Evolution in Aerodynamic Optimization // Proc. of 46th Annual Conf. of the Canadian Aeronautics and Space Institute. Canada, Montreal, 2007. P. 29-36.
[81 Marjavaara B.D., Lundstrom T.S. Response surface-based shape optimization of a Francis draft tube // Intern. J. Numerical Methods Heat k, Fluid Flow. 2007. Vol. 17, iss. 1. P. 34-45.
[91 Lighthill J. A New Method of Two Aerodynamics Design / Aeronautical Res. Council's Rep. and Memoranda. No. 2112. London, R&M, 1945. 45 p.
[101 Лесохин А.Ф. Расчет лопастей рабочих колес осевых турбин (решетка профилей конечной толщины) // Тр. .'IIIII. 1953. № 5. С 49-55.
[11] Степанов Г.Ю. Гидродинамика решеток турбомашин. М.: Физматгиз, 1962. 512 с.
[12] Garabedian P., Korn D. A systematic methods for computer design of supercritical airfoil in cascade // Comm. Pure Appl. Math. 1976. Vol. 29. P. 369-382.
[13] Этинберг И.Э., Раухман B.C. Гидродинамика гидравлических турбин. Л.: Машиностроение (Ленингр. отд-ние), 1978. 280 с.
[14] Викторов Г.В. Гидродинамическая теория решеток: Уч. пособие. М.: Высшая шк., 1969. 368 с.
[15] Топаж Г.И. Расчет интегральных гидравлических показателей гидромашин. Л.: Изд-во Ленингр. гос. ун-та, 1989. 208 с.
[16] borges J.E. A Three-dimensional inverse method for turbomachinerv. Part I // J. Turbo-machinery. 1990. Vol. 112. P. 346-354.
[17] Dang Т., isgro V. Euler-based inverse method for turbomachine blades. Part 1: Two dimensional cascades // AIAA J. 1995. Vol. 33, No. 12. P. 2309-2315.
[18] Demeulenaere A., Van den Braembussche R.A. Three-dimensional inverse method for turbomachinerv blading system //J. Turbomachinerv. 1998. Vol. 120. P. 247-255.
[19] Dang Т., Damle S., QlU X. Euler-based inverse method for turbomachine blades. Part 2: Three-dimensional flows // AIAA J. 2000. Vol. 38, No. 11. P. 2007-2013.
[20] QlU X., Dang T. Three-dimensional viscous inverse method for axial blade design // Inverse Problems Sci. and Eng. 2009. Vol. 17, No. 8. P. 1019-1036.
[21] pascoa J.C., Mendez A.C., Gato L.M.C. A fast iterative method for turbomachinerv blade design // Mech. Res. Commun. 2009. Vol. 36. P. 630-637.
Поступила в редакцию 7 мая 2010 г.
