Об одном подходе к оптимизации формы лопасти гидротурбины Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 10, № 6, 2005

ОБ ОДНОМ ПОДХОДЕ К ОПТИМИЗАЦИИ ФОРМЫ ЛОПАСТИ ГИДРОТУРБИНЫ*

И.Ф. ЛОБАРЕВА Новосибирский государственный университет, Россия

В.А. Скороспелов, П.А. Турук Институт математики СО РАН, Новосибирск, Россия

С.Г. Черный, Д.В. Чирков Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]

An approach to automatic blade shape optimization of hydroturbine runner is presented. It is based on parameterization of blade surface, 3D inviscid flow simulation in the runner, and genetic algorithm (Breeder GA) for finding a global minimum of the prescribed objective function. Three different objective functions are considered. Capabilities of the proposed method are demonstrated in application to runner blade shape optimization for Bratskaya hydropower station refurbishment project.

Введение

В настоящее время повышение качества гидротурбин стало весьма актуальной проблемой, необходимость решения которой обусловлена острой конкуренцией на мировом рынке среди их производителей, а также тенденцией роста требований к техническим характеристикам. Так, КПД, требуемый заказчиками, превышает 93.5 % для радиально-осевых турбин и 92 % для поворотно-лопастных. Значительно возросли требования к кавитационным характеристикам, допустимым напряжениям, вибрации, шуму и др.

Проектирование форм компонентов турбомашин, при которых выполняются предъявляемые к турбомашинам требования, в настоящее время, как правило, проводится опытными проектировщиками вручную методом возмущения формы некоторого известного прототипа с использованием инженерных формул, а также численного и экспериментального исследования пространственного потока в проточном тракте. Такой подход к совершенствованию проточной части турбомашин чрезвычайно трудоемок, так как требует перебора большого количества комбинаций геометрических параметров и анализа соответствующих гидродинамических полей.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-000246) и Интеграционным проектом СО РАН (№ 27).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

В последнее время в мире появляется все больше и больше работ, посвященных проблеме автоматизации процесса проектирования и оптимизации форм элементов гидротурбин. Это связано со значительным прогрессом в развитии методов вычислительной гидродинамики, методов решения оптимизационных задач и повышением производительности вычислительной техники. Общее в подходе в этих работах состоит в предварительной параметризации рассматриваемой поверхности и последующем автоматическом переборе различных комбинаций геометрических параметров, численном моделировании течения в полученных проточных трактах и отыскании (с помощью некоторой стратегии) такой комбинации параметров, которая будет давать минимум заданной целевой функции.

Большинство публикаций посвящены автоматической оптимизации формы основного элемента — рабочего колеса радиально-осевой [2-4, 6-8] и поворотно-лопастной [1, 5] гидротурбины. В [6] рассмотрена оптимизация формы лопаток направляющего аппарата. Работы [9-11] посвящены оптимизации формы отсасывающей трубы.

Для параметризации поверхностей активно используется аппарат кривых Безье и сплайн-функций, при этом, как правило, варьирование формы осуществляется путем вариации нескольких десятков параметров. Так, в [2] форма лопасти, образующих обода и ступицы рабочего колеса определяется 27 параметрами. В [3] оптимизируется только форма срединной поверхности лопасти, которая задается 35 параметрами. В [5] для варьирования поверхности лопасти поворотно-лопастной турбины используются девять параметров, задающих углы входной и выходной кромки, и семь параметров, задающих распределение толщины. В работе [10] форма симметричной отсасывающей трубы определяется 10 свободными параметрами.

При оптимизации рабочего колеса во всех указанных выше работах расчетной областью является проточная часть между двумя соседними лопастями рабочего колеса (межлопастной канал), захватывающая некоторое пространство перед рабочим колесом и часть диффузора за ним. Расчет стационарного пространственного течения в межлопастном канале проводится в предположении, что течения во всех остальных межлопастных каналах циклически повторяются. Такой подход позволяет существенно сократить объем вычислений и для большинства режимов работы гидротурбины вполне приемлем.

Режим работы гидротурбины определяется заданием расхода Q и напора Н или заданием расхода Q и частоты вращения п. Перед началом процесса оптимизации фиксируются рабочий режим турбины и, как правило, некоторый исходный вариант рабочего колеса. В работах [1-3, 5, 6, 9, 10] оптимизация проводится для теоретически предсказанного режима наибольшего КПД (оптимального режима). В [7] кроме оптимального учитываются также режимы максимальной и неполной загрузки. После задания рабочего режима турбины проводится расчет течения в рабочем колесе с учетом направляющего аппарата и отсасывающей трубы, необходимый для определения граничных условий на входе в рабочее колесо. В дальнейшем, в процессе модификации формы рабочего колеса, эти условия остаются неизменными.

Гидродинамические расчеты перебираемых конфигураций в работах [1, 3-7, 9-11] проводятся путем решения уравнений Навье — Стокса несжимаемой жидкости. При этом для моделирования турбулентности, как правило, используется к — е-модель [1, 7, 9, 10, 11]. В работе [2] расчет течения в рабочем колесе проводится в рамках уравнения Эйлера.

Вообще говоря, основными показателями качества гидротурбины являются КПД, с которым гидравлическая энергия преобразуется в механическую, а также ее кавитацион-ная устойчивость. Однако для корректного расчета КПД необходимо проводить расчет во всей гидротурбине с учетом направляющего аппарата и особенно отсасывающей тру-

бы с использованием достаточно подробных сеток и адекватных моделей турбулентности. Подобные расчеты требуют больших вычислительных затрат и пока не могут быть использованы для нужд оптимизации. Поэтому в большинстве работ для построения целевого функционала используются косвенные данные, такие как, например, потери полного давления в рабочем колесе (фактически КПД одного рабочего колеса), качество поля скорости на выходе из рабочего колеса, т. е. на входе в диффузор отсасывающей трубы. Кроме того, в качестве функционалов, характеризующих качество гидротурбины, рассматриваются вращающий момент, связанные с кавитационной устойчивостью размер и положение зоны низкого давления на тыльной стороне лопасти и т. д.

Среди работ по оптимизации рабочего колеса стоит выделить [2], поскольку здесь гидродинамические расчеты проводятся в рамках уравнений Эйлера идеальной несжимаемой жидкости. При использовании такой модели расчет КПД невозможен. Поэтому здесь для оценки качества лопасти используются такие косвенные параметры, как характер распределения давления от ступицы к ободу на рабочей стороне лопасти, монотонность понижения давления от входной кромки лопасти к выходной, характер распределения меридиональной и окружной компонент скорости на входе в отсасывающую трубу. Запас устойчивости к кавитации характеризуется в [2] превышением давления потока на тыльной стороне лопасти над давлением насыщенного пара.

В работах [9-11], посвященных оптимизации отсасывающей трубы, в качестве целевого функционала, подлежащего максимизации, используется основной показатель качества диффузора — коэффициент восстановления давления.

Множество требований, предъявляемых к рабочему колесу (как правило, противоречивых) , диктует необходимость использования специальных методов при решении задачи оптимизации. Здесь можно выделить несколько подходов. Первый состоит в рассмотрении в качестве целевого функционала, подлежащего минимизации, взвешенной суммы всех интересующих функций качества. Такой подход использован, в частности, в работе [2]. Недостаток его очевиден: для адекватного выбора весовых коэффициентов для каждой новой турбины необходимы накопление статистики оптимизационных расчетов с различными комбинациями весов и тщательное сопоставление результатов. Второй подход предполагает при минимизации одного функционала накладывать дополнительные ограничения на величину всех остальных интересующих функционалов. Однако наиболее перспективным, с точки зрения авторов, является подход многоцелевой оптимизации, в котором улучшение формы лопасти проводится сразу по нескольким функционалам [1, 6, 7, 12]. В этом случае результатом решения задачи оптимизации является не одна форма, а множество форм, образующих фронт Парето.

Среди различных методов поиска глобального минимума целевого функционала применительно к задаче оптимизации форм элементов турбомашин наибольшее распространение получили различные модификации генетического алгоритма (Genetic Algorithm, GA). Он использован в работах [1-3, 5, 7, 8]. Из достоинств этого алгоритма следует отметить основное — он допускает параллельный расчет нескольких десятков конфигураций, образующих поколение, что позволяет в десятки раз сократить общее время проведения оптимизации. В [4, 6, 9] используются как GA, так и метод градиентного спуска и проводится их сравнение. В работе [11], посвященной оптимизации отсасывающей трубы, поиск минимума осуществляется методом градиентного спуска, а в [10] используется так называемый "метод поверхности отклика" (Response Surface Methodology).

В работах [1, 7, 12], где проводится многоцелевая оптимизация, используется Multi-Objective Genetic Algorithm (MOGA), являющийся модификацией GA. При осуществлении

многоцелевой оптимизации оптимум ищется у нескольких функционалов Fi, ... , Fm одновременно. Отметим, что если функционалы не являются независимыми, то решение такой задачи не единственно, а образует множество, внутри которого невозможно улучшить значение любого из функционалов без ухудшения значений каких-либо остальных. Эта группа решений называется множеством Парето. При использовании классических оптимизационных алгоритмов, таких как градиентные методы, задача нахождения множества Парето сводится к проведению серии из k оптимизационных расчетов, в каждом из которых минимизируется сумма с весами Wj интересующих проектировщика функционалов

M

Fi,... , FM, называемая функцией качества: Fk = WjFj. Структура GA позволяет

j=l,...,M

проводить многоцелевую оптимизацию без рассмотрения последовательности функций качества и получать множество Парето после проведения одного оптимизационного расчета [12].

В настоящей работе описывается подход к оптимизации формы лопасти радиально-осевого рабочего колеса, идейно близкий работе [2], где оптимизация проводится в рамках уравнений Эйлера. Основным достоинством модели идеальной жидкости является быстрота расчета: с использованием разработанного авторами алгоритма расчет трехмерного поля течения для одной конфигурации требует нескольких минут. Это обстоятельство, а также тот факт, что влияние вязких эффектов в рабочем колесе гидротурбины мало, подтолкнуло авторов настоящей работы к проведению оптимизации лопасти рабочего колеса в приближении уравнений Эйлера, направленной на улучшение кинематических свойств потока. При этом, как отмечалось выше, особо актуальной становится задача выбора подходящего целевого функционала. Рассмотрены следующие варианты целевых функционалов, подлежащих минимизации: кинетическая энергия потока, покидающего рабочее колесо, степень отклонения линий тока на поверхности лопасти от "осесимметричных" линий тока, а также размер области низкого давления на тыльной стороне лопасти. Необходимо отметить, что на данный момент оптимизация проводится только относительно одного из них. Важное место в работе занимают вопросы выбора метода поиска глобального минимума целевого функционала. На тестовых задачах проведено сравнение двух оптимизационных алгоритмов: алгоритма ПОИСК [18], который стоит ближе к детерминированным алгоритмам, и генетического алгоритма Breeder Genetic Algorithm [2]. В качестве объекта оптимизации выбран один из вариантов нового рабочего колеса для Братской ГЭС.

1. Общее описание алгоритма оптимизации

Постановка задачи. Пусть задана некоторая начальная лопасть, подлежащая модификации (оптимизации). Зафиксирован режим работы гидротурбины посредством задания: расхода Q, напора Н и частоты вращения п. Требуется модифицировать форму лопасти так, чтобы она была лучше начальной с точки зрения некоторого критерия качества (целевого функционала Г) и удовлетворяла заданным ограничениям (например, отсутствие кавитации).

Схема решения. Вариация формы лопасти в данной работе осуществляется путем варьирования значений нескольких геометрических параметров х\,...,х^, определяющих кривизну срединной поверхности лопасти.

Схематически процесс оптимизации состоит из следующих шагов (рис. 1):

1. Берутся начальные значения параметров XI,...,х^, соответствующие начальной

Рис. 1. Схема процесса оптимизации.

форме лопасти рабочего колеса.

2. По значениям Х\,... , Xn восстанавливается срединная поверхность лопасти, к которой добавляются фиксированные распределения ее толщины d(u,v). После этого в межлопастном канале рабочего колеса строится регулярная конечно-разностная сетка.

3. Полученная сетка подается на вход модуля расчета невязкого 30-течения в рабочем колесе.

4. По рассчитанному полю скорости и давления вычисляется значение целевого функционала F.

5. Оптимизационный алгоритм получает значение функционала F и вырабатывает новый набор параметров Х\,... , Xn. Далее повторяется шаг 2.

Указанный процесс продолжается до тех пор, пока оптимизационный алгоритм не найдет минимум функционала F в допустимой области изменения параметров. Схема работы Breeder Genetic Algorithm несколько отличается от указанной выше, о чем будет подробно сказано в разд. 6.3.

2. Параметризация лопасти

Варьирование формы лопасти осуществляется путем варьирования числовых значений конечного числа параметров Х\,..., хN, задающих ее форму. Параметризация формы лопасти должна быть, с одной стороны, достаточно гибкой, с другой — содержать по возможности меньшее число параметров. Ниже описана использованная в данной работе "относительная" параметризация срединной поверхности лопасти бикубической функцией с 16 свободными параметрами.

С рабочим колесом радиально-осевой гидротурбины связана декартова система координат OXYZ, ось OZ которой совпадает с осью вращения колеса и направлена по потоку

Рис. 2. RZ-проекция лопасти.

в колесе. Поверхность лопасти представлялась в виде ее срединной поверхности

ф,^) = {X(и,^),У(и^)}, и,^ е [0,1], (1)

и функции распределения толщин й = ^(и,^). Здесь (и,^) — так называемая "естественная" параметризация поверхности лопасти (рис. 2). Уравнение срединной поверхности (1) может быть записано в цилиндрической системе координат:

Я = Я(и,и) = у/Х (и,^)2 + У (и,^)2, £ = £ (и^),

N У (и,^) Ф = Ф(и, V) = аг^^-^.

X (и, V)

Характер формы срединной поверхности в направлении от входной кромки к выходной во многом определяется зависимостью Ф(и, V), поэтому в данной работе вариация срединной поверхности осуществлялась путем вариации только функции Ф(и,^) при сохранении Я£-проекции и заданного изначально распределения толщин й(и, V). Пусть Ф0(и, V) — функция угловой координаты исходной лопасти, а Ф(и, V) — функция угловой координаты возмущенной лопасти. В данной работе для варьирования формы лопасти использовался "относительный" подход, т.е. параметризации подлежала не сама зависимость Ф(и,^), а отклонение Ф'(и, V) = Ф(и, V) — Ф0(и, V). Вариация Ф'(и, V) осуществлялась путем вариации 16 коэффициентов бикубического полинома, представленного в форме Бернштейна:

3 3

Ф' М ) = ^ В<(и)В (V), (2)

¿=0 у=0

3! _

где Б&(эд) = ку(з—к)!^(1 — . Достоинством рассмотренной "относительной" параметризации является тот факт, что исходная лопасть, которая, вообще говоря, имеет сложную форму, всегда лежит в множестве допустимых параметризацией форм (при = 0). Для удобства геометрические параметры в дальнейшем обозначаются . Для реальных лопастей функция угловой координаты Ф(и,^) обычно является монотонно возрастающей по переменной V: (и, V) < 0 (небольшое нарушение условия монотонности может иметь место лишь вблизи входной кромки). Для выполнения этого условия на геометрические параметры накладываются соответствующие ограничения.

Кроме того, для каждого параметра х^- задается диапазон его изменения (так называемые фазовые ограничения):

3. Расчет 3В-поля течения в рабочем колесе

Для каждого набора геометрических параметров расчеты течения в полученном рабочем колесе проводились в невязком приближении в рамках уравнений Эйлера несжимаемой жидкости. Расчеты выполнялись в одном межлопастном канале рабочего колеса при предположении, что течения в остальных межлопастных каналах циклически повторяются (рис. 3).

3.1. Основные уравнения

Во вращающемся с постоянной угловой скоростью ш рабочем колесе рассматривается относительное течение в системе координат у1,у2,у3, вращающейся вместе с колесом. Связь между векторами скорости V в абсолютном движении и w = (/ш\,/ш2,/ш33) — в относительном задается выражением

где ш — вектор угловой скорости, ш = (0,0,ш); г — радиус-вектор точки. Уравнения Эйлера для относительного движения записываются в виде

хь,з < Х3 < хк,з

V = w + ш X г,

(3)

w = 0;

(4)

(5)

Здесь ¥ = (/ь/2,/3) = -ш хш хг - 2ш xw + g.

Рис. 3. Расчетная область и сетка в ней.

Рис. 4. Меридиональная проекция расчетной области.

3.2. Краевые и начальные условия

В качестве входных данных для расчета течения помимо сетки должны быть заданы частота вращения рабочего колеса n, распределение вектора скорости во входном сечении AA', а также распределение давления или условие радиального равновесия — в выходном сечении BB' (рис. 4). Это условие получается из радиальной составляющей уравнения количества движения

с dcr + с« dcr + c dCr - Cu = - dP (6)

r dr r z dz r dr

при предположении о малости радиальной составляющей скорости cr ~ 0 и осесиммет-ричности течения. В этом случае уравнение (6) при условии z = const, соответствующем выходному сечению расчетной области, превращается в обыкновенное дифференциальное уравнение

dp = ^u (7)

dr r

На боковых границах выше и ниже лопасти ставятся условия периодичности. На твердых стенках задается условие непротекания w • n = 0, где n — единичная нормаль к твердой стенке.

Во время всего оптимизационного процесса частота вращения рабочего колеса, распределение вектора скорости во входном сечении и закон распределения давления в выходном сечении расчетной области предполагались неизменными и соответствующими заданному (Q, H, n) режиму работы гидротурбины с исходным рабочим колесом. В качестве начального приближения для 2-го и всех последующих расчетов бралось установившееся поле скорости и давления, полученное при расчете предыдущей конфигурации. Такой подход позволил сэкономить около 50 % времени счета.

3.3. Численный метод

3.3.1. Метод искусственной сжимаемости и конечных объемов решения уравнений Эйлера несжимаемой жидкости

Одна из главных проблем при разработке алгоритмов численного решения уравнений Эйлера несжимаемой жидкости с использованием переменных скорость — давление связана с выбором способа расчета давления, обеспечивающего соленоидальность поля скорости. Для ее решения в настоящей работе используется метод искусственной сжимаемости, заключающийся в добавлении в уравнение неразрывности (4) производной по времени

dP + в div w = 0, (8)

где в есть коэффициент искусственной сжимаемости. Для построения консервативного численного алгоритма модифицированные уравнения (8), (5) записываются в виде интегральных законов сохранения для произвольного фиксированного объема V

д J QdV + у H • dS = J FdV, (9)

V dV V

где dV — замкнутая поверхность объема V; dS = n • dS — элемент поверхности, умноженный на единичную внешнюю нормаль к ней; Q = (p, wi,w2,w3)T; H = (вw,wiw + pei,w2w + pe2,w3w + pes)T; F = (0, / ,/2,/з)Т.

Рис. 5. Ячейка расчетной сетки.

Обозначим средние значения д и Е в ячейке расчетной области с индексами (г,], к) и объемом Уф* через

Щ* = I Р" йУ, Пк = ± ( ЕйУ

фк

V]

Уцк

1]к «

Уг^к

и отнесем их к центру ячейки (рис. 5). Верхний индекс п означает номер слоя по времени. Дискретизируем уравнение (9)

д^1 - д"

А

У = БН8

"+1

(10)

где А* — шаг по времени

ИН8 = — (Нг+1/2 — Нг—1/2 + Нф+1/2 — Нф-1/2 + Н*+1/2 — Н*-1/^ + ЕУ

Разностные потоки через грани ячейки Н^±1/2 вычисляются с использованием ТУБ-схемы

Чакраварти — Ошера [14, 15]:

Н

1

т+1/2

(Нт+1 + Нт) • 8т+1/2 — |А|т+1/2 Ат+1/2д

- W

т+1/2 5

(11)

А =

^ = А+ + А-, |А|= А+ - А-.

Матрица Якоби А представляется в виде А = ИОЬ, где И и Ь — матрицы правых и левых собственных векторов соответственно, ИЬ = ЬИ = I; О = diag (Л1, Л2, Л3, Л4), Л1)2 = и,

Л3,4 = и ± д/ и2 + в 8 • 8, и = w • 8, А± = - И (О ± |О|) Ь; Wm±1/2 содержит дополнитель-

2

ные члены, повышающие порядок аппроксимации до второго или третьего [14, 15].

3.3.2. Решение нелинейных уравнений

Для линеаризации системы уравнений (10) применяется метод Ньютона

"У ( д I^ - ^ИН8 а* V д д I

(д"+1 - д") = ИН8

(12)

2

Используя при вычислении дИНЯ/дQ в разностных потоках только направленные разности первого порядка, получим из (12) систему линейных уравнений

V - +

1 А + А-+1/2 Лг+1/2 + А+-1/2А^-1/2 +

+ Aj+1/2AJ+1/2 + A+-1/2Aj-i/2 +

+ A-+1 /2 Ak+1/2 + A+-1/2 Ak-1/2

An+1Q = RHSn, (13)

где An+1Q = Qn+1 - Qn.

Система (13) приближенно LU-факторизуется и разрешается за два дробных шага

An+1/2Qjfc = B-1 ( RHSn+A+-1/2An+1/2Qi-1jfc+A+_1/2An+1/2Qij_1fc+A^ An+1/2Qijfc_1 An+1Qijfc = A"+1/2 Qijfc — B 1 ^Ai+1/2Ara+1Qj+1jfc + A-+1/2An+1Qjj+1fc + Afc+1/2An+1Qjjfc+1 где

V _ + - + _ +

B = 1 At — Ai+1/2 + Ai_1/2 — Aj+1 /2 + Aj_1/2 — Afc+1 /2 + Afc_1/2"

Уравнения каждого из дробных шагов решаются по формулам бегущего счета в направлении возрастания и убывания индексов. Описанный выше метод аналогичен использованному ранее авторами в работах [15-17].

3.3.3. Выбор коэффициента искусственной сжимаемости в

Для обеспечения высокой скорости сходимости итерационного процесса и точности пространственной аппроксимации коэффициент искусственной сжимаемости в нужно выбрать так, чтобы собственные значения матрицы A в (11) были одного порядка. Оценивая величины собственных значений матрицы A

Л1,2 = U, Лз,4 = U U2 + вS ■ S, U = w ■ S,

получаем, что коэффициент искусственной сжимаемости должен быть порядка |w|2.

Аналогичный результат получается из условия хорошей обусловленности матрицы R, т. е. из условия, что число ее обусловленности близко к единице:

K = ||R|| ■ ||L|| « 1.

Для выполнения этого условия в работе [13] установлено, что коэффициент искусственной сжимаемости должен вычисляться по формуле

в = max (0.3, r|w|2) , 1 < r < 5.

В настоящей работе коэффициент искусственной сжимаемости в полагался

в = 8Wx2

хар'

где 8 = 5 ... 10; Wxap — величина характерной скорости задачи.

4. Ограничения

Помимо линейных ограничений на геометрические параметры ставились также ограничение на напор в рабочем колесе и кавитационное ограничение.

Ограничение на напор. Чтобы полученное в результате оптимизации колесо соответствовало заданному изначально режиму работы, необходимо, чтобы разность полных энергий потока на входе в спиральную камеру гидротурбины и на выходе из отсасывающей трубы совпадала с напором Н, обеспечиваемым станцией для данного режима работы. Однако в процессе оптимизации рассчитывается течение только в рабочем колесе. Поэтому для сохранения заданного напора станции ограничение накладывалось на так называемый напор в рабочем колесе Нс, который исходя из оценки потерь в элементах реальной гидротурбины составляет около 95 % от полного напора. Таким образом, ограничение на напор имеет вид

\НС - 0.95Н\< е, (14)

где е — малый параметр.

Кавитационное ограничение. Важнейшими в работе гидротурбины являются ее кавитационные характеристики. Явление кавитации возникает в случае, если величина статического давления в точке потока становится меньше, чем величина давления насыщенного пара при данной температуре:

Р < Pv. (15)

При расчете течения в рабочем колесе по уравнениям Эйлера поле давления в потоке определяется с точностью до постоянного давления р0, заданного в выходном сечении расчетной области. Для получения абсолютного давления в потоке, используемого в условии (15), необходимо оценить абсолютную величину давления р0. Это осуществлялось с помощью уравнения Бернулли, записанного для линии тока, идущей от точки "0" до точки "*" в отсасывающей трубе (рис. 6).

Известно, что наиболее подвержены кавитации тыльные поверхности лопастей. В случае, когда область кавитации занимает более 15 % площади всей тыльной поверхности лопасти, наблюдается резкое падение мощности и КПД турбины, что крайне нежелательно. Предварительные оптимизационные расчеты рабочего колеса Братской ГЭС показали,

Рис. 6. К оценке давления ро-

что путем варьирования только угловой координаты Ф(и, у) в рамках 16-параметрического ее задания полностью избавиться от кавитации на заданном режиме невозможно, однако удается снизить площадь области кавитации до значений, меньших допустимого 0.15Б8ИС, где Б8ИС — площадь тыльной поверхности лопасти. Поэтому ограничение по кавитации ставилось в следующем виде:

< 0.15, (16)

где Бсау — площадь зоны на тыльной поверхности лопасти, где выполняется (15).

5. Целевые функционалы

Как уже говорилось, в идеале оптимизация геометрии турбомашины должна сводиться к максимизации КПД на заданном режиме работы или, что то же, к минимизации потерь энергии в проточном тракте. В случае расчета одного рабочего колеса в рамках модели невязкой жидкости явный расчет потерь невозможен, поэтому необходимо сформулировать целевой функционал, косвенно учитывающий потери энергии в проточном тракте. При этом неоценимую помощь оказывает опыт специалистов-проектировщиков. Ниже представлены формулировки некоторых целевых функционалов.

1. Ех — кинетическая энергия в выходном сечении рабочего колеса. Предположим, что вся кинетическая энергия потока на выходе из рабочего колеса теряется. Тогда наибольшая эффективность работы гидротурбины достигается в том случае, когда поток покидает рабочее колесо с минимальной кинетической энергией. В соответствии с этим в качестве минимизируемого функционала логично взять кинетическую энергию потока в выходном сечении 50 расчетной области:

Р = 1/ М2(у■ ¿8). (17)

Здесь V — вектор абсолютной скорости; ¿8 — вектор нормали к элементарной площадке ¿Б выходного сечения, равный по модулю площади этого сечения и направленный вовне расчетной области. При фиксированном расходе минимум функционала (17) соответствует равномерному распределению осевой компоненты сх и нулевой окружной компоненты си абсолютной скорости в выходном сечении 50.

2. Е2 — относительный размер области кавитации. В качестве минимизируемого функционала берется относительная площадь области кавитации на тыльной стороне лопасти:

р=. (18)

3. Е3 — отклонение линий тока от "осесимметричных". Интуиция опытных проектировщиков подсказывает, что более качественной является та лопасть, на поверхности которой предельные линии тока ближе к линиям тока "осесимметричного" потока. Соответствующий целевой функционал можно сформулировать следующим образом:

Р = 1^(1 - а (в) сое в) ¿Б, (19)

5

Рис. 7. Предельные линии тока на рабочей стороне лопасти.

где Б — площадь поверхности лопасти; в — угол между предельной линией тока и линией тока осесимметричного потока (рис. 7), а весовая функция а имеет вид

а(в) Г 1, в < п/2,

а(в} П ао, в > п/2.

За счет выбора подходящего весового коэффициента а0 минимизация этого целевого функционала позволяет также выполнить одно из основных требований к лопастной системе — отсутствие угла атаки при обтекании лопасти (совпадение линии растекания жидкости с линией входной кромки лопасти).

6. Алгоритмы поиска минимума функционала 6.1. Постановка оптимизационной задачи

Математически задача оптимизации формы лопасти формулируется следующим образом: найти

min F(x), x = (xi,... , xN) £ X

при ограничениях

X = {x : XL,i < Xi < XR,i} , (20)

<Pj(x) < 0, j = l,...,m, (21)

^fc(x) < 0, k = m + 1,...,Q, (22)

где x — вектор параметров, определяющих геометрию лопасти; F — целевой функционал; (20) — фазовые, (21) — гидродинамические, (22) — геометрические ограничения.

В данной работе исследовались два оптимизационных алгоритма: ПОИСК, разработанный в ИТПМ СО РАН, и генетический алгоритм Breeder Genetic Algorithm (BGA). Ограничения при решении оптимизационной задачи в обоих случаях учитывались методом штрафных функций путем введения составного функционала Fc.

6.2. ПОИСК

Данный алгоритм является симбиозом детерминированных и случайных методов поиска с элементами самообучения. Он реализован в виде специализированного комплекса программ ПОИСК, синтезированного на основе модифицированных методов вращающихся координат, направляющего конуса и матричного спуска [18]. Составной функционал имеет следующий вид:

Я

Рс = Р 1 + Е К

3 = 1

j +

\ 2"|

j

) _

(23)

Здесь £ — задаваемая точность выполнения ограничений; К3- — коэффициенты штрафа, начальные значения которых задаются и уточняются в процессе работы программы. Для настройки программы ПОИСК пользователю предоставляется возможность задать точность выполнения ограничений, точность определения минимума по параметрам, точность по функционалу. Так как и целевой функционал, и ограничения нелинейным образом зависят от определяющих параметров задачи (х\,... , х^), решение, получаемое алгоритмом ПОИСК, вообще говоря, зависит от выбранной начальной точки в пространстве варьируемых параметров, требуемой точности выполнения ограничений, критериев прекращения работы программы и т. д.

6.3. Breeder Genetic Algorithm

В основу BGA положена теория естественного отбора Дарвина, согласно которой популяция индивидуумов меняется в течение нескольких поколений путем рекомбинации (скрещивания) и мутации отдельных индивидуумов, подчиняясь при этом законам естественного или искусственного отбора относительно предписанного критерия. Применительно к проектированию лопасти рабочего колеса гидротурбины индивидуумом является рабочее колесо с лопастью какой-то конкретной формы. Также индивидуумом мы будем называть набор параметров x = (xi,... , Xn), определяющих эту лопасть. Вероятность перехода данного индивидуума в следующее поколение зависит от того, насколько хорошо он удовлетворяет требованиям проектировщика. В данной работе реализован вариант BGA, сходный с описанным в работе [2]. Составной функционал имеет следующий вид:

Fc = F

1 + 5] C (j + I jI)

(24)

где С3 — весовые коэффициенты, которые задаются при запуске программы и в процессе оптимизации не изменяются.

Процесс оптимизации состоит из следующих шагов (рис. 8):

1. Формируется начальная популяция, состоящая из р0 индивидуумов. Каждый индивидуум — это набор параметров х = (х1,... , х^), определяющих объект, подлежащий оптимизации, при этом хг £ [х^.г, хд,г]. Начальная популяция создается случайным образом применительно ко всем параметрам так, чтобы каждый индивидуум удовлетворял геометрическим ограничениям. При этом в начальную популяцию включается индивидуум, соответствующий начальной геометрии лопасти.

2. Для каждого индивидуума в популяции вычисляется соответствующее значение составного функционала Рс.

m

Рис. 8. Схема алгоритма BGA.

3. Отбор (selection). На этом этапе выживает определенная доля Tr (0 < Tr < 1) наилучших индивидуумов. В данном случае это будут индивидуумы, которым соответствуют наименьшие значения составного функционала. Остальные индивидуумы отбрасываются.

4. Построение новой популяции. На данном этапе строится новая популяция, состоящая из индивидуумов, согласно следующей схеме.

A. Рекомбинация (recombination). Из выживших Trp индивидуумов случайным образом выбираются два индивидуума-родителя — x = (xi,...,XN ) и y = ), которые производят новый индивидуум w = (w1, :,Wn). Это реализуется следующим образом: wi = xi + аДу — Xi), где ai — случайное число между — d и 1 + d. Параметр d называется параметром рекомбинации. Необходимо отметить, что ai может быть больше единицы и меньше нуля, так что рекомбинация включает в себя элемент экстраполяции. На этапе рекомбинации число индивидуумов восстанавливается до p.

Б. Мутация (mutation). На этапе мутации все индивидуумы слегка меняются следующим образом: wi = wi ± ^(wL;i — wR;i)5, где ^ G [0,1] — мутационный параметр, 5 = 2-16y, а y — случайное число между 0 и 1.

B. Клонирование (cloning). Для того чтобы наилучший индивидуум в каждом поколении не потерялся в результате мутации и рекомбинации, он копируется неизменным в следующее поколение.

5. Переход на шаг 2, пока не будет рассчитано Noen поколений.

В задаче оптимизации формы лопасти для нахождения глобального минимума функционала обычно требуется рассчитать число поколений Noen ~ 30.

Выполнение геометрических ограничений легко проверить, так как они могут быть

сформулированы явно в терминах x = (xi, :,Xn)• Поэтому новая популяция строится так, чтобы эти ограничения заведомо выполнялись. Ниже описывается один из способов реализации этого подхода. Пусть D G X — область параметров, в которой выполняются геометрические ограничения. При невыполнении хотя бы одного из геометрических ограничений для индивидуума x из вновь построенной i-й популяции (т. е. при x G X\D) он заменяется индивидуумом w, лежащим на границе области D, который строится следующим образом: w = x + a(xm — x) G ÖD, a G [0,1]. В качестве xm используется индивидуум из (i—1)-го поколения, соответствующий минимальному значению составного функционала Fc и лежащий в допустимой геометрическими ограничениями области D.

6.4. Методическое исследование оптимизационных алгоритмов

Главным критерием качества оптимизационного алгоритма является гарантия поиска глобального, а не локального минимума функционала. Другим ценным качеством является скорость поиска, т. е. число комбинаций параметров, которые требуется перебрать, чтобы найти минимум целевого функционала.

В данном разделе приводятся результаты сравнения алгоритмов ПОИСК и BGA на трех тестовых функционалах. Одной из целей исследования было отыскание такого набора внутренних параметров алгоритма BGA (d, р, р0, Tr, ß), который позволяет осуществлять поиск глобального минимума функционала с максимально возможной скоростью.

Предварительные исследования показали, что на работу алгоритма BGA оказывают существенное влияние два параметра: параметр рекомбинации d и число индивидуумов в популяции р. Поэтому в данном исследовании значение d варьировалось в диапазоне [0, 2], параметр р принимал значения 40, 50, 60, 70, а остальные параметры в силу их слабого влияния фиксировались и принимали значения р0 = 60, ß = 0.1, Tr = 0.3. Необходимо однако отметить, что в случае сложного вида целевого функционала с множеством локальных минимумов количество индивидуумов в начальной популяции р0 также будет играть существенную роль.

Поскольку алгоритм BGA является стохастическим, для исследования его работоспособности проводился статистический анализ серии запусков. Так, на каждом тестовом функционале для каждого набора параметров (d, р, р0, Tr, ß) проводилось 70 запусков, после чего определялась вероятность нахождения алгоритмом глобального минимума.

Результат работы программы ПОИСК, хоть и имеющей элемент случайности, при фиксированных внутренних параметрах алгоритма, как правило, не меняется от запуска к запуску. В этом смысле ПОИСК близок к детерминированным алгоритмам. Предварительные расчеты показали также, что при наличии в функционале локальных минимумов результат алгоритма ПОИСК зависит от выбора начальной точки. Поэтому при исследовании программой ПОИСК варьировалась также начальная точка.

Тестовые функционалы зависят от 20 независимых переменных, принимающих значения в области X = {x : —10 < Xj < 10}:

20

(25)

i=1

20

i=1

20

FT3 (x) = 1+ g(xi)+ g(x2 ) + жД (27)

i=3

где функция $(ж) имеет следующий вид (рис. 9):

{(ж + 3)2, x G [-10, 0],

- (ж - 3)2 + 18, ж G [0, 5], (ж - 7)2 + 10, ж G [5,10].

Функционалы (25) и (26) имеют единственный глобальный минимум. Функционал (27) имеет глобальный минимум в точке (-3, - 3,0,..., 0) и три локальных минимума в точках (-3, 7, 0,..., 0), (7, -3, 0,..., 0), (7, 7, 0,..., 0) (см. рис. 9).

В табл. 1-3 приведены результаты поиска минимумов тестовых функционалов по алгоритмам ПОИСК и BGA при различных значениях d и p для BGA и различных стартовых точках для алгоритма ПОИСК. В случае исследования алгоритма ПОИСК на функционалах Ft1 , Ft2 использовалась одна стартовая точка, поскольку в данном случае этот алгоритм находит минимум независимо от стартовой точки, что ясно из описания метода и вида функционалов.

Как видно из табл. 1, 2, в случае наличия у функционала единственного ярко выраженного минимума ПОИСК уверенно находит этот минимум менее чем за 1500 итераций, при этом (в таблицах не указано) определяет его с высокой точностью. Таблица 3 иллюстрирует склонность алгоритма ПОИСК к нахождению локальных минимумов при неудачном выборе стартовой точки. Генетический алгоритм BGA с большой вероятностью сходится к глобальному минимуму при d = 0.7 и при этом требует существенно большего количества вычислений целевого функционала, нежели ПОИСК. Для функционала Ft1 было рассчитано 40 поколений, для Fy2 и Fy3 — 80 поколений. Необходимо однако отметить, что при использовании алгоритма BGA зачастую требуется рассчитать лишь 20-30 поколений, с дальнейшим увеличением числа поколений результат практически не улучшается.

В работе [2] рекомендован следующий набор параметров BGA: d = 0.25, p = 40, p0 = 60, ^ = 0.1, Tr = 0.3, при этом функционал зависит от 27 переменных. Однако, как видно из табл. 1-3, ни для одного из рассмотренных тестовых функционалов такой набор не позволяет найти глобальный минимум. Отсутствие сходимости к глобальному минимуму

Рис. 9. Функция g(x) и изолинии функции 1+g(xi) + g(x2).

Таблица1

Сравнение алгоритмов ПОИСК и BGA при минимизации функционала Ft1

Алгоритм ПОИСК BGA

Внутренние параметры алгоритма Ро = 70, d = 0.25 Ро = 70, d = 0.7 Ро = 40, d = 1.5

Стартовая точка (10,..., 10) — —

Результат Найден минимум Процесс не сходится к минимуму Минимум найден в 97% случаев Минимум найден в 97 % случаев

Число 1364 — 2860 1260

итераций

Сравнение алгоритмов ПОИСК и Таблица2 BGA при минимизации функционала Fy2

Алгоритм ПОИСК BGA

Внутренние параметры алгоритма Ро = 40, 50, 60, 70; d = 0.25 Ро = 70, d = 0.7 Ро = 40, 50, 60, 70; d = 1.5

Стартовая точка (10,..., 10) — — —

Результат Найден минимум Процесс не сходится к минимуму Минимум найден в 97 % случаев Процесс не сходится к минимуму

Число 1404 — 5660 —

итераций

ТаблицаЗ

Сравнение алгоритмов ПОИСК и BGA при минимизации функционала Fy3

Алгоритм ПОИСК BGA

Внутренние — Ро = 40, 50, Ро = 70, Ро = 40, 50,

параметры 60, 70; d = 0.7 60, 70;

алгоритма d = 0.25 d = 1.5

Стартовая (8, 8, (-8, -5, (4, 2, — — —

точка 10,...,10) 10,...,10) 10,..., 10)

Результат Найден Найден Найден Процесс не Глобальный Процесс не

локальный глобальный локальный сходится ни минимум сходится ни

минимум минимум минимум к локальному, ни к глобальному минимуму найден в 97% случаев к локальному, ни к глобальному минимуму

Число 1520 1325 1640 — 5660 —

итераций

при малых d можно объяснить следующим образом. Параметр d определяет величину "отскока" индивидуума, переходящего в следующее поколение, от родителей-индивидуумов из текущего поколения. Таким образом, при малых значениях параметра рекомбинации, например при d = 0.25, поиск осуществляется в окрестности оставшихся после отбора индивидуумов начальной популяции. Как показали дополнительные расчеты, с увеличением начальной популяции до p0 = 200, 500,5000 при d = 0.25 достичь сходимости к глобальному минимуму не удалось. Согласно результатам таблицы, большие значения параметра рекомбинации (d = 1.5) также не обеспечивают сходимость эволюционного процесса к глобальному экстремуму.

На основании проведенных тестовых расчетов можно заключить, что хоть число оптимизационных шагов, необходимых для нахождения минимума функционала у алгоритма ПОИСК в большинстве случаев меньше, чем у BGA, генетический алгоритм более применим к задаче оптимизации лопасти, так как при сложном виде целевого функционала он с большей вероятностью осуществляет поиск глобального минимума.

7. Результаты оптимизационных расчетов

В качестве объекта оптимизации было выбрано радиально-осевое рабочее колесо гидротурбины Братской ГЭС (см. рис. 3). Сетка в межлопастном канале имела размеры 40 х 20 х 20. Для получения стационарного решения с точностью div(w) < 10-3 на этой сетке требуется от 1000 до 4000 шагов по времени, что составляет соответственно от 1 до 4 мин времени работы процессора Athlon 2500+. Поиск минимума целевой функции осуществлялся с помощью алгоритма BGA с параметрами d = 0.7, p = 70, p0 = 200, ß = 0.1, Tr = 0.3. Нахождение оптимального решения требует построения Аоеп ~ 30 поколений, что соответствует перебору и расчету 2-3 тысяч конфигураций. Таким образом, в целом процесс поиска оптимальной формы занимает около трех суток. Ниже представлены результаты оптимизационных расчетов, полученные при минимизации каждого из рассмотренных в разд. 5 функционалов.

На рис. 10 представлены начальная лопасть и лопасть, полученная путем минимизации функционала (17) при ограничении на напор (14) и кавитационном ограничении (16). Здесь же показаны соответствующие распределения меридиональной проекции вектора скорости на выходе из расчетной области. Как и следовало ожидать, для оптимальной лопасти

V

Рис. 10. Минимизация кинетической энергии на выходе из расчетной области: вид лопасти со стороны входной кромки и распределение скоростей на выходе; слева — начальная лопасть, справа — оптимальная.

поле скорости существенно более равномерное, при этом кинетическая энергия потока в выходном сечении была уменьшена в полтора раза.

На рис. 11 приведены распределения давления на тыльной поверхности для начальной лопасти и для лопасти, полученной в результате минимизации функционала (18), отвечающего за кавитацию, при единственном ограничении на напор (14). Заштрихованная область — зона кавитации, где выполняется условие (15). В результате оптимизации площадь зоны кавитации уменьшилась с 5 % до примерно 1 % от площади всей тыльной поверхности лопасти.

Результаты минимизации отклонения линий тока от "осесимметричных" показаны на рис. 12. При решении задачи оптимизации ставились ограничения на напор (14) и размер области кавитации (16). Проведено два расчета на минимизацию функционала в первом а0 = 1, во втором а0 = 10. Видно, что при а0 = 1 линии тока на поверхности лопасти стали существенно ближе к осесимметричным линиям тока, однако при этом линия растекания удалилась от пересечения срединной поверхности с входной кромкой лопасти. При а0 = 10 линии тока практически не изменились по сравнению с исходной лопастью, но линия

Рис. 11. Минимизация кавитации на тыльной стороне лопасти: слева — начальная лопасть, справа — оптимальная.

Рис. 12. Линии тока на рабочей стороне лопасти: а — начальная, б — оптимальная при ао в — оптимальная при ао = 10.

а

в

а б в

Рис. 13. Линии тока в окрестности входной кромки: а — начальная, б — оптимальная при сто = 1, в — оптимальная при сто = 10; 1 — линия пересечения срединной поверхности и входной кромки лопасти, 2 — линия растекания.

растекания полностью совпала с линией пересечения срединной поверхности с входной кромкой лопасти. На рис. 13, где в другом ракурсе показана окрестность входной кромки, видно, что при ст0 = 10 удалось полностью устранить перетекание жидкости с рабочей стороны на тыльную в окрестности входной кромки и тем самым обеспечить нулевой угол атаки и снизить ударные потери в рабочем колесе. Таким образом, выбором константы ст0 в функционале можно влиять на результат оптимизации: при малых значениях ст0 будут выравниваться линии тока на всей поверхности лопасти, но при этом положение линии растекания будет учитываться слабо, и наоборот.

Итак, по результатам настоящего раздела можно судить, что созданный инструмент позволяет достаточно эффективно улучшать качество лопасти с точки зрения какого-то одного фиксированного критерия. Однако улучшение какой-то одной характеристики лопастной системы, как правило, влечет ухудшение других характеристик. Поэтому для практического применения системы автоматического проектирования необходимо научиться оптимизировать лопасть сразу по нескольким критериям, что составляет направление дальнейших исследований.

Заключение

Рассмотрен подход к формализации процесса проектирования формы лопасти гидротурбины путем постановки и решения соответствующей оптимизационной задачи. На основе разработанного подхода создана система автоматической оптимизации, которая позволяет достаточно эффективно улучшать качество лопасти с точки зрения какого-то одного фиксированного критерия.

Остается открытым вопрос выбора целевых функционалов для решения задачи оптимизации в рамках уравнений Эйлера, наиболее отвечающих требованиям проектировщиков.

Показано, что выбранная параметризация с 16 свободными параметрами позволяет достаточно эффективно варьировать форму лопасти. Однако для оптимизации всего рабочего колеса этого недостаточно. В дальнейшем планируется увеличение числа геометрических параметров с целью обеспечения возможности варьировать распределение толщины на лопасти, форму обода и ступицы, а также положение входной и выходной кромок в

RZ-проекции рабочего колеса. Среди направлений дальнейшей работы можно также отметить распараллеливание BGA, переход на модель турбулентного течения и исследование

методов многоцелевой оптимизации, таких как MOGA.

Список литературы

[1] Lipej A., Poloni C. Design of Kaplan runner using multi-objective genetic algorithm optimization // J. of Hydraulic Research. 2000. Vol. 38. P. 73-77.

[2] Sallaberger I., Fisler M., Michaud M. et al. The design of Francis turbine runners by 3D Euler simulations coupled to a breedergenetic algorithm // Proc. of 20th IAHR Symp. on Hydraulic Machinery and Systems. 2000.

[3] Tomas L., Pedretti C., Chiappa T. et al. Automated design of a Francis turbine runner using global optimization algorithms // Proc. of 21st IAHR Symp. on Hydraulic Machinery and Systems. 2002.

[4] Mazzouji F., Francois M., Tomas L. et al. Refinements in Francis turbine design // Hydropower & Dams, Issue one. 2004. P. 53-58.

[5] Kueny J.-L., Lestriez R., Helali A. et al. Optimal design of a small hydraulic turbine // Proc. of 22nd IAHR Symp. on Hydraulic Machinery and Systems. 2004.

[6] Mazzouji F., Couston M., Ferrando L. et al. Multicriteria optimization: viscous analysis-mechanical analysis // Proc. of 22nd IAHR Symp. on Hydraulic Machinery and Systems. 2004.

[7] 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. 2004.

[8] Ferrando L., Kueny J.-L., Avellan F. et al. Surface parametrization of a Francis turbine for optimum design // Proc. of 22nd IAHR Symp. on Hydraulic Machinery and Systems. 2004.

[9] Eisinger R., Ruprecht A. Automatic shape optimization of hydro turbine components based on CFD // Seminar "CFD for Turbomachinery Applications", Gdansk, Sept. 2001.

[10] Puente L.R., Reggio M., Guibault F. Automatic shape optimization of a hydraulic turbine draft tube // Proc. of Intern. Conf. CFD2003, Vancouver, May 28 — June 1, 2003.

[11] Marjavaara B.D., Lundstrom T.S. Shape optimization of a hydropower draft tube // Proc. of 22nd IAHR Symp. on Hydraulic Machinery and Systems. 2004.

[12] Favaretto C.F.F., Funazaki K.-I., Tanuma T. The development of a genetic algorithm code for secondary flow injection optimization in axial turbines // Proc. of the Intern. Gas Turbine Congress, Tokyo, Nov. 2-7, 2003.

[13] Rizzi A.W., Ericsson L.-E. Computation of inviscid incompressible flow with rotation //J. Fluid. Mech. 1985. Vol. 153. P. 275-312.

[14] Chakravarthy S.R., Osher S. A New Class of High Resolution TVD Schemes for Hyperbolic Conservation Laws. AIAA Paper 85-0363, 1985.

[15] Грязин Ю.А., Черный С.Г., Шаров С.В.Численное моделирование течений несжимаемой жидкости на основе метода искусственной сжимаемости // Вычисл. технологии: Сб. науч. тр. / РАН. Сиб. отд-ние. ИВТ. 1995, Т. 4, № 13. C. 180-203.

[16] Грязин Ю.А., Черный С.Г., Шаров С.В., Шашкин П.А. Об одном методе численного расчета решения трехмерных задач динамики несжимаемой жидкости // Докл. РАН. 1997. Т. 353, № 4. С. 478-483.

[17] Cherny S.G., Sharoy S.V., Skorospeloy V.A., Turuk P.A. Methods for three-dimensional flows computation in hydraulic turbines // Russ. J. Numer. Anal. Math. Modelling. 2003. Vol. 18, N 18. P. 87-104.

[18] Латыпов А.Ф., Никуличев Ю.В. Специализированный комплекс программ оптимизации. Препр., 1985 (АН СССР. Сиб. отд-ние. ИТПМ; №15-85.)

Поступила в редакцию 22 июля 2005 г.