CC BY
36
Вестник ПНИПУ. Аэрокосмическая техника. 2019. № 57
Б01: 10.15593/2224-9982/2019.57.06 УДК 62-135.1
В.И. Милешин1, ИХ. Орехов1, C.К. Щипин2
1 Центральный институт авиационного моторостроения им. П.И. Баранова, Москва, Россия 2 АО «Российская самолетостроительная корпорация "МиГ"», Москва, Россия
ПРОЕКТИРОВАНИЕ ТРАНСЗВУКОВОГО РОТОРА КОМПРЕССОРА НИЗКОГО ДАВЛЕНИЯ НА ОСНОВЕ РЕШЕНИЯ 3D-ОБРАТНОЙ ЗАДАЧИ С ЦЕЛЬЮ ОБЕСПЕЧЕНИЯ ЕГО РАБОТОСПОСОБНОСТИ В ШИРОКОМ ДИАПАЗОНЕ ЗНАЧЕНИЙ ЧАСТОТЫ ВРАЩЕНИЯ
Рассматриваются особенности проектирования трансзвукового ротора 1-й ступени компрессора низкого давления без входного направляющего аппарата на основе решения ЗЭ-обратной задачи с помощью программного комплекса ЗО-ИЧУЕВвЕ.ЕХВ!. с целью обеспечения его работоспособности в широком диапазоне значений частоты вращения.
Обратная задача строится на задании предпочтительного распределения статического давления на поверхности лопатки, заданной толщине лопатки и разности значений статического давления (называемой нагруженностью) в соответствующих точках спинки и корыта лопатки. Входные и выходные газодинамические параметры берутся из решения прямой задачи о течении в многоступенчатом компрессоре и остаются неизменными в процессе решения обратной задачи.
Решение обратной задачи строится на использовании подвижной расчетной сетки (изменение потоков расхода, импульса и энергии, обусловленное движением граней ячейки, учитывается). Ввиду этого каждый шаг по времени начинается с определения изменения поверхностей лопатки с последующим построением новой вычислительной сетки.
В данной работе обратная задача обеспечивает эффективность и работоспособность ротора 1-й ступени КНД в широком диапазоне значений частоты вращения (ппр = 70...100 %) для случая отсутствия входного направляющего аппарата. Для этого перепроектирование лопатки ротора разбивается на две части. Сначала производится модификация лопатки в окрестности передней кромки для ппр = 70 % с целью минимизировать или совсем убрать срыв потока из-за больших углов атаки. Затем производится модификация только 2-й половины лопатки, включая окрестность задней кромки, чтобы обеспечить параметры ротора на проектных оборотах ппр = 100 %.
Ключевые слова: трансзвуковой ротор, проектирование, компрессор низкого давления, ЗЭ-обратная задача, широкий диапазон значений частоты вращения.
V.I. Mileshin1, I.K. Orekhov1, S.K. Shchipin2
1 Central Institute of Aviation Motors named after P.I. Baranov (CIAM), Moscow, Russian Federation 2 Russian Aircraft Corporation "MiG", Moscow, Russian Federation
3D-INVERSE DESIGN OF TRANSONIC FAN ROTORS EFFICIENT FOR A WIDE RANGE OF RPM
The paper is dealing with the peculiarities in design of first rotor of multistage low pressure compressor (LPC) in case of absence of intel guide vane (IGV) on base of in-house software package "3D-INVERSE.EXBL" user to solve 3D invers problem. Presented in the paper inverse solution provides effectiveness and operability of LPS first rotor for a wide range of RPM.
Inverse problem is based on desired static pressure distribution on suction side of blade, given blade thickness and pressure difference (named loading) in corresponding points of suction and pressure sides of blade. Inlet and outlet gas-dynamic parameters (pressure, density and flow velocity vector) are taken from direct solution of flow within multistage compressor and remains unchanged during inverse problem solution.
Solution of inverse problem is determined using moving grid. Normal speed of face of grid cell adjacent to blade surface is determined using given static pressure (inverse mode) with the aid of relationships, which are the elements of Godunov scheme applied for integration of flow equations.
In the paper inverse solution provides effectiveness and operability of first rotor of multistage low-pressure compressor (LPC) for a wide range of rpm (70+100 %) in case of absence of inlet guide vane (IGV).
Keywords: transonic fan rotor, design, low-pressure compressor, 3D-Inverse problem, wide range speed.
Введение
ЗБ-обратная задача исследована достаточно глубоко в работах [1-14]. В этих работах ЗБ-обратная задача применена к проектированию роторов в точке одной выбранной характеристики при выбранной частоте вращения ротора. В работе [3] ротор, улучшенный на 100 % оборотов, дал улучшение также и на 90 % оборотов.
В разработанном методе решения ЗБ-обратной задачи [5] поверхность лопатки является непроницаемой и на ней выставляются граничные условия прилипания или другие условия, совместимые с уравнениями Навье - Стокса. Решение обратной задачи строится на использовании подвижной расчетной сетки, которая перестраивается на каждом шаге итерационной процедуры установления по времени (изменение потоков расхода, импульса и энергии, обусловленное движением граней ячейки, учитывается). Скорость перемещения границы ячейки, лежащей на поверхности лопатки по нормали к ней, находится из условия заданного распределения статического давления по соотношениям, которые являются элементом схемы Годунова [15]. Граница ячейки, лежащая на поверхности лопатки, в процессе перемещения остается непроницаемой. Координаты узлов поверхности лопатки на новом временном слое находятся по специальному алгоритму, который обеспечивает самостабилизацию новой поверхности профиля (отсутствие изломов, пилы и т.д.) без использования процедуры сглаживания [5, 16].
В случае вязкого потока расчетная сетка должна быть адаптивной и густой в областях больших градиентов параметров (пограничные слои, ударные волны, передние и задние кромки профилей лопаток и т.д.). Алгоритм построения поверхностей спинки и корыта лопатки должен надежно работать в случае неравномерной конечно-разностной сетки при решении ЗБ-обрат-ной задачи.
Хорошо известно, что ЗБ-обратная задача при проектных оборотах ппр = 100 % используется для снижения интенсивности ударной волны путем перепрофилирования поверхности лопатки с целью уменьшения числа Маха перед ударной волной либо конструирования двухскач-ковой системы (косой скачок уплотнения плюс прямой) путем расщепления одной ударной волны высокой интенсивности.
Для современных вентиляторов необходимо обеспечить эффективность и работоспособность ротора 1-й ступени компрессора низкого давления (КНД) в широком диапазоне оборотов вала (70-100 %), особенно для случая отсутствия входного направляющего аппарата (ВНА).
На 70 % оборотов отсутствие ВНА приводит к большим углам атаки на передней кромке ротора 1-й ступени КНД [17]. Ввиду этого периферийные сечения лопатки ротора работают в срыве. Как правило, в реальных конструкциях эта проблема решается применением надро-торного устройства. Хотя надроторное устройство (НРУ) и увеличивает запас устойчивой работы, но может вызывать ухудшение КПД на проектных оборотах. Исходя из этого важно улучшить течение на 70 % оборотов (при ппр = 70 %) только за счет изменения формы лопатки, без использования НРУ.
3Б-обратная задача для дозвукового предсрывного потока сама по себе непроста. Стараясь уменьшить угол атаки, приходится увеличивать кривизну лопатки на передней кромке. Однако в этом случае максимальный массовый расход через венец на 100 % оборотов имеет тенденцию к уменьшению.
Как известно, для увеличения запаса устойчивой работы на 100 % оборотов межлопаточный скачок должен быть сдвинут к задней кромке лопатки при уменьшении числа М перед скачком для увеличения КПД. Для достижения этого лопатка ротора должна быть перепрофилирована в окрестности задней кромки.
Та же самая задача на 70 % оборотов достигается путем подстройки передней кромки лопатки к течению с большими углами атаки. Это требует перепроектирования лопатки в окрестности передней кромки с использованием уменьшенных нагрузок.
Эта работа посвящена развитию и применению метода решения 3Б-обратной задачи для проектирования венцов компрессора. Подправленное распределение нагрузки на 70 % оборотов используется для усовершенствования геометрии лопатки ротора, чтобы повысить КПД и запасы устойчивой работы современного трехступенчатого КНД [17].
Затем решается прямая задача для всего компрессора на 70 и 100 % оборотов. Трехмерные расчеты всего трехступенчатого КНД в целом показали, что КПД ротора 1-й ступени увеличился на 4,85 %, степень повышения давления возросла на 4 на 70 % оборотов. Чтобы КПД КНД на 100 % оборотов не уменьшился, интенсивность ударной волны должна быть снижена путем решения ЗБ-обратной задачи. Обычно требуется 2 или 3 коррекции по обратной задаче, чтобы найти распределение нагрузки, обеспечивающее требуемые параметры КНД в широком диапазоне оборотов (70-100 %).
Система уравнений для ЗО-обратной задачи в трехмерной постановке
и элементы численной схемы
Нестационарные уравнения Навье - Стокса могут быть записаны в безразмерной форме в декартовой системе координат, вращающейся с угловой скоростью ю относительно оси х, в относительном движении следующим образом:
где е - удельная внутренняя энергия; к - коэффициент теплопроводности; тхх, тху, тх;, ...т;; -компоненты тензора вязких напряжений; ю - частота вращения ротора.
Декартовы составляющие скорости (и, V, ж) отнесены к критической скорости звука, посчитанной по параметрам на входе, плотность р - к критической плотности на входе, линейный размер - к величине периферийного радиуса передней кромки лопатки, вязкость ц - к ламинарной вязкости воздуха при критической температуре набегающего потока, давление р - к величине скоростного напора, подсчитанного по значениям величин критической плотности и скорости звука на входе. Величина скорости звука а определяется из уравнения
ЕV = ( |(х + иХ ) + ^ Ж |(иу + Vx ) + Жх ) иТхх + ^ху + WT XZ + fx ^ ^^
К = ( |(у + ^ )|(у + Vy ) + ^ Ж + Жу ) итух + ^уу + Жу; + ^ ) Ке-1;
Gv =(0,+ ж),+ Wy),+ Ж; ) + ХШуж,ит^ + vтZy + жтгг + fz) Ие-1;
fx = к (к-1)-1 Рг-1 (а2 )х; fy = к (к-1)-1 Рг-1 (а2 )у; fz = к (к-1)-1 Рг-1 (а2);
д = (р, ри, рv, рж, е)'; 5 = (0,0,р(ю2у + 2юж),р(ю2; - 2юv),рю2^у + ж;)) ;
е = р(р /(к-1)р + (и2 + V2 + ж2) /2); Х = -2/3ц,
Параметр к - отношение значений удельной теплоемкости для газа. При проведении расчетов принималось, что ламинарное число Прандтля Р^ = 0,72, а турбулентное Ргг = 0,9.
Для описания турбулентных течений использована алгебраическая двухслойная модель вихревой вязкости, предложенная Болдуином и Ломаксом [18].
Для интегрирования уравнений используется модифицированная схема С.К. Годунова [15, 19, 20]. В соответствии с модифицированным методом Годунова, нестационарные уравнения Навье - Стокса записываются в форме интегральных законов сохранения для элемента объема ёт с учетом скорости О движения границ этого элемента:
Здесь Л(г) - произвольно изменяющийся по времени г объем; Дг) - его граница; ео8(я,х), ео8(я,у), ео8(и,г) - косинусы углов внешней нормали к поверхности границы Ь({) с соответствующими осями х, у, z; Ох, Оу и Ог - компоненты вектора скорости О границы Ь(г), движущейся по нормали к себе.
Стационарные распределения функций течения вырабатываются в процессе установления по времени при достаточно произвольных начальных и граничных условиях. Установление по времени проводится с использованием подвижной конечно-разностной сетки. Деформация расчетной сетки определяется движением соответствующей поверхности лопатки, строящейся на каждом временном слое. Интегрирование системы уравнений (1) проводится по явной схеме Годунова повышенного порядка точности [15, 19, 20] в обратной постановке и по неявной схеме в прямой постановке.
На границах рассчитываемой области функции течения должны удовлетворять заданным граничным условиям.
На входной границе задается распределение по радиусу полного давления, полной температуры в абсолютном движении, угол между окружной компонентой скорости в абсолютном движении и меридиональной плоскостью и угол между радиальной компонентой скорости и осью вращения компрессора. Из расчетной области на входную границу вдоль характеристики С- сносится левый инвариант Римана. На выходной границе задается величина статического давления на втулке, а распределение статического давления вдоль радиуса определяется из уравнения радиального равновесия. На периодических границах задаются условия периодичности. На твердых стенках - условие прилипания (или проскальзывания в случае использования пристеночных функций) в случае решения прямой задачи или распределение статического давления по спинке и корыту лопатки в случае решения обратной задачи. Корректное задание статического давления по спинке и корыту лопатки, удовлетворяющее условию отсутствия самопересечения поверхности спинки и корыта, заранее трудно получить. Исходя из этого удобно использовать способ, позволяющий избежать проблемы самопересечения профилей. Суть этого способа заключается в том, что статическое давление задается либо по спинке, либо по корыту лопатки, а другая сторона поверхности лопатки находится из условия заданного распределения толщины профиля лопатки по высоте.
В системе многоступенчатого компрессора граничные условия на входе и выходе из расчетной области берутся в сечениях mixing plane после решения прямой задачи. На рис. 1 показаны первые две поверхности МР1 и МР2, в которых были использованы граничные условия в приближении mixing plane. 3Б-обратная задача решалась для отдельного лопаточного венца
(E - Ev - qDx)cos(n,x) +
+ (G - Gv - qDz )cos(n, z)
(1)
Граничные условия
КНД, в котором нужно было улучшить течение и, соответственно, его интегральные параметры. Для обеспечения согласования венцов в системе КНД граничные условия для решения ЗБ-обратной задачи брались из решения прямой задачи сквозного расчета ЗБ-вязкого течения потока в КНД.
MPI
а б
Рис. 1. Сетка типа Н: а - меридиональное сечение; б - сечение в окружном направлении
от лопатки к лопатке
Например, для решения 3Б-обратной задачи с целью оптимизации РК1 левая граница расчетной области совпадала с поверхностью mixing plane МР1, а правая с поверхностью МР2 (см. рис. 1).
При решении обратной задачи граничные условия для потока на поверхностях МР1 и МР2 совпадают с соответствующими граничными условиями в приближении mixing plane при решении прямой задачи.
Конечно-разностная сетка
Для простоты реализации алгоритма решения обратной задачи используется сетка типа Н со сгущением узлов возле поверхностей лопатки в окружном направлении и возле передней и задней кромок в меридиональном направлении. В окружном направлении к поверхности лопатки сетка сгущается по экспоненциальному закону. Сетка строится по аналитическим формулам, поэтому при ее пересчете на каждом новом временном слое тратится минимум вычислительного времени (см. рис. 1).
Алгоритм построения профиля
Как известно, в схеме Годунова по параметрам в соседних ячейках сначала решаются вспомогательные задачи о распаде произвольного разрыва. В частности, по параметрам в ячейках, примыкающих к поверхности лопатки, и по заданной величине статического давления на поверхности лопатки, которая должна соответствовать этому давлению, находится нормальная
к поверхности лопатки составляющая скорости D. Для этого при решении задачи о распаде произвольного разрыва на грани ячейки разностной сетки, прилегающей к поверхности лопатки, положим, что самой поверхности лопатки будет соответствовать контактный разрыв. Таким образом, построение формы лопатки при решении обратной задачи сводится к выделению поверхности контактного разрыва при рассмотрении стандартной задачи для схемы Годунова -задачи о распаде произвольного разрыва на поверхности лопатки. Решение данной задачи используется для построения поверхности лопатки как для уравнений Эйлера, так и для уравнений Навье - Стокса. Только для уравнений Навье - Стокса используется модифицированная
схема Годунова. Также подчеркнем, что в соответствии с определением контактного разрыва через него не происходит протекание потока, т.е. на каждом временном слое г + ёг на движущемся элементе поверхности лопатки выполняется условие твердой стенки - условие непротекания. Скорость контактного разрыва будет определять скорость движения элемента поверхности лопатки. Например, для решетки лопаток, приведенной на рис. 1, скорость перемещения спинки лопатки Б.. по нормали к себе определяется по формуле
= (РР.. - р.)
О.. — - ,
Л0..
а скорость перемещения корыта лопатки Ор. - по формуле
( - Рр:: )
Бр. — --
Л0 р8
здесь РР:: и РРр: - заданные распределения статического давления по спинке и по корыту лопатки соответственно; р:: и рр: - величины статического давления в центрах граней ячеек, прилегающих к лопаткам; А0 - массовая скорость.
Величина А0 определяется по формуле, которая для спинки лопатки может быть представлена в следующем виде: Если РР
Если РР** < р* то
Л0 —. 1^ РР„ + 1 р.
1 - РР
Л0 = 5 1 Р а _р
Л0 = 0 Р
25
1-
(РР ^
5-1 ' 25
р*
Аналогичная формула используется для расчета скорости контактного разрыва со стороны поверхности корыта лопатки.
После того как найдены все нормальные к элементам поверхности лопатки, скорости их перемещения Б.. или Ор., можно определить смещение поверхности лопатки по нормали к себе за временной шаг ёг. Поскольку конфигурация сетки в момент г + ёг определяется положением узлов на направляющих линиях сетки, то ее построение состоит в определении смещений узлов, а не смещений центров граней ячеек, составляющих поверхность лопатки. На рис. 1 дана схема сетки и расчетной области. Сеточные линии, выходящие из одной поверхности лопатки к другой («направляющие рельсы»), фиксируются, а сеточные линии, пересекающие их, являются подвижными.
Определение скорости перемещения элемента поверхности лопатки - важный этап алгоритма построения поверхности лопатки, но недостаточный для его успешного построения. Это обусловлено тем, что начальное распределение функций течения в ячейках расчетной области может быть достаточно произвольным: на поверхность лопатки могут приходить сильные ударные волны или влиять другие особенности течения, приводящие к большим градиентам параметров. Большие неоднородности функций течения вдоль поверхности лопатки могут вызвать появление пилообразной структуры поверхности лопатки, что может привести в конечном счете к неустойчивости решения обратной задачи. Для того чтобы такой неустойчивости не возникало, алгоритм построения поверхности лопатки должен обладать свойством самостабилизации. Это свойство должно обеспечивать ограниченность скорости перемещения элемен-
та лопатки вдоль направляющих рельсов, даже в случае образования пилообразной структуры на поверхности лопатки. Для лопатки с изломами поверхности «геометрическая скорость» движения узлов граней ячейки по направляющей может существенно превышать скорость движения грани по нормали к себе. Для ограничения этой скорости воспользуемся аналогией с известным принципом Гюйгенса [5, 17]. Напомним, как по указанному принципу строится фронт акустической волны. Каждая точка фронта в момент Х рассматривается как источник возмущений, распространяющихся относительно газа со скоростью звука а во всех направлениях и к тому же сносящихся со скоростью потока W. Фронт волны получается как граница возмущенной области.
В схеме Годунова, в которой по параметрам на слое Х из решения задачи о распаде произвольного разрыва находятся нормальные компоненты скорости центров всех элементов контактного разрыва, построение контура поверхности контактного разрыва можно вести аналогичным образом. При этом роль скорости звука играет своя для каждого элемента нормальная скорость О движения этого элемента.
Таким образом, построение элемента поверхности контактного разрыва (элемент поверхности лопатки) на временном слое Х + йХ начинается путем смещения узлов этого элемента Ь1(Х), Ых), Ь3(Х), Ь4(Х) на временном слое Х по нормали к себе на величину А = ОйХ. Но в отличие от слоя Х на временном слое Х + йХ к границам этого элемента будут прилегать четыре сферы с радиусами А = ОйХ и с центрами в узлах этого элемента Ь1(Х), Ь^Х), Ь3(Х) и Ь4(Х). Кроме того, к границам элемента Ь1(Х + йХ) Ь2(Х + йХ) Ь3(Х + йХ) Ь4(Х + йХ) на слое Х + йХ будут также прилегать четыре цилиндра с радиусами А = ОйХ и с осями вдоль соответствующего ребра элемента контактного разрыва на временном слое Х. Например, на рис. 2 показано, что к узлам Ь1(Х + йХ) и Ьг(Х + йХ) прилегают участки сферических поверхностей М1(Х + йХ) Ь1(Х + йХ) М8(Х + йХ) и М3(Х + йХ) Ь2(Х + йХ) М2(Х + йХ) с центрами сфер в узлах Ь1(Х) и ^(Х) соответственно и с радиусами А = ОйХ. Между этими двумя участками сферических поверхностей расположен участок цилиндрической поверхности М1(Х + йХ) М2(Х + йХ) Ь2(Х + йХ) Ь1(Х + йХ) с осью цилиндра ^(Х)^(Х) и с радиусами А = ОйХ. Аналогично на рис. 2 показаны участки сферических и цилиндрических поверхностей для остальных точек и ребер элемента поверхности контактного разрыва Ь1(Х + йХ) ^(Х + йХ) Ь3(Х + йХ) Ь4(Х + йХ) на слое Х + йХ.
Рис. 2. Схема построения элемента поверхности лопатки на временном слое Х + йХ с учетом участков сферических и цилиндрических поверхностей, примыкающих к ребрам грани конечно-разностной ячейки
Л<0
М-,
м.
I
2 t-
м,
\
\
л/й
В соответствии с вышесказанным на рис. 3 показана проекция границы «возмущенной» области, образованной от перемещения элемента контактного разрыва Ь1(г)Ь2(г)Ь3,(г)Ь4(г) с временного слоя г на временной слой г + ёг. Важно подчеркнуть, что максимальное удаление точек границы возмущенной области от границ четырехугольника Ь1(г)Ь2(г)Ь3(г)Ь4(г) не превышает величину А = Бёг.
На следующем этапе полученная кон- м,
фигурация элемента контактного разрыва ь1(г + ёг) Д(г + ёг) ь3(г + ёг) ь4(г + ёг) на слое
/ 1 ¿л» ' ,
> > /I \ М-
I + ёг сносится на вектор И7,., ск. Здесь
есть вектор скорости потока в центрах ячеек, прилегающих к поверхности спинки лопатки. В общем случае величина вектора Wss отлична от нуля и может достигать значительных значений, когда используется модель турбулентности с пристеночными функциями либо когда обратная задача решается в невязкой постановке для уравнений Эйлера.
На рис. 4 показана схема выбора узлов конечно-разностной сетки ТС1ТС2Тс3Та4 по направляющей рельсе Р3 для построения поверхности лопатки на временном слое г + ёг.
Для простоты анализа на этой фигуре показана только направляющая рельса, выходящая из узла Ь3(г). Направляющие рельсы являются сеточными линиями 3-го семейства, которые выходят из узлов, расположенных на поверхности одной лопатки, и входят в узлы, расположенные на поверхности соседней лопатки. Подчеркнем, что направление рельс не изменяется в процессе построения поверхности лопатки.
/ 1 Оск \ \
/ \ / \ __'______« --- "Г/
Ш I ' 1
_и4(0 !
Ш ¡\ ш? 1
\ \ ч ч К М7\ V
Рис. 3. Проекция границы «возмущенной» области, образованной от перемещения элемента контактного разрыва Ь1(г)Ь2(г)Ь3(г)Ь4(г) с временного слоя г на временной слой г + ёг
Рис. 4. Схема выбора узлов конечно-разностной сетки ТС1 Та2 Та3 Та4 по направляющей рельсе Р3 для построения поверхности лопатки на временном слое г + ёг
В соответствии с вышеуказанным принципом Гюйгенса построения «возмущенной» области на временном слое t + dt построим образы четырех элементов поверхности контактного разрыва G1, G2, G3 и G4 на слое t окружающих рельс P3. Более отдаленные грани ячеек от рельса P3 в анализе участвовать не будут в силу ограничения шага по времени dt - поверхность лопатки выделяется по явной схеме. На рис. 4 индекс n соответствует продольной сеточной линии, индекс m соответствует радиальной сеточной линии и индекс o - направляющим рельсам.
Первоначально каждая поверхность элементов G1, G2, G3 и G4 переносится на соответствующий вектор (Dssdt) с образованием необходимых участков сфер и цилиндров. После этого полученные элементы поверхностей на временном слое t + dt переносятся на соответствующий каждому элементу вектор (Wssdt). Далее находится пересечение смещенных элементов поверхностей контактного разрыва G1, G2, G3 и G4 (t + dt) с рельсой P3. Анализ этой фигуры показывает, что максимальное удаление от узла L3(t) имеет точка TG2,
полученная пересечением элемента G2 (t + dt) после его смещения на вектор (Wssdt) .
\ ' G2
Именно эта точка TG2 будет новым узлом поверхности лопатки на временном слое t + dt на рельсе P3.
Оптимизация рабочего колеса осевой ступени на оборотах n = 70 %
Анализ исходного варианта рабочего колеса R1 в точке рабочей линии на промежуточных оборотах n = 70 % показывает, что при высоте лопатки 70 % < h < 100 % образуется срыв потока с передней кромки лопатки, что хорошо видно на рис. 5. Указанный срыв потока существенно снижает запасы газодинамической устойчивости и может явиться источником возбуждения срывного флаттера. На основе трехмерной обратной задачи рабочее колесо R1 было перепроектировано с целью устранения срыва потока с передней кромки лопатки R1. С этой целью были модифицированы распределения статического давления на стороне разрежения лопатки во всех сечениях лопатки R1 от втулки до периферии. На рис. 6 показан пример такой модификации для трех сечений лопатки R1 - 10, 60 и 90 % высоты лопатки h . На этом рисунке видно, что для модифицированного распределения статического давления обрезан пик давления в окрестности передней кромки. Далее линия модифицированного распределения статического давления проведена подобно линии для исходного давления на спинке лопатки, но несколько ниже исходной линии. Это необходимо для сохранения суммарной нагрузки в каждом сечении лопатки. Здесь под нагрузкой понимается площадь между двумя кривыми, соответствующими статическому давлению на корыте и спинке.
Далее для каждой ячейки конечно-разностной сетки при x = const вычисляется разность значений статического давления на корыте и спинке лопатки, которая затем держится постоянной во время решения 3Б-обратной задачи и используется в этой задаче в качестве граничного условия на поверхности лопатки. 3Б-обратная задача решается на основе программного комплекса 3D-INVERSE.EXBL, разработанного в ЦИАМ. Некоторые результаты решения этой задачи приведены на рис. 7, на котором показаны профили исходной и модифицированной лопаток для трех сечений по высоте: h = 10, 60 и 90 %. Анализ этого рисунка показывает, что у модифицированных профилей существенно увеличился изгиб профиля в окрестности передней кромки с целью устранения срыва потока.
Рис. 5. Исходный вариант Ш, п = 70 %. Распределение чисел Маха в трех сечениях по высоте лопатки: к = 90 %, периферия (а); к = 60 %, среднее сечение (б) и к = 10 %, втулка (в) на поверхности вращения и в меридиональном сечении вдоль стороны разрежения (г) и вдоль стороны давления (д)
а
г
б
д
Р* ' у .....//'' .......1.............:
...................... ............\ Лт
■ 1+ /г4-р; : / Модифицировании? ........./--/-.............;...............1.............1} —Г Р ...... ......
0,0130162976 0,111648626 0,2164624
0,583532 0,16910145
а)периферия
б) среднее сечение
в) втулка
Рис. 6. Пример задания модифицированной нагрузки для решения 3Б-обратной задачи для трех сечений по высоте лопатки при п = 70 %: а - к = 90 %, б -Н = 60 %, в - к = 10 %; Ря - распределение статического давления по спинке лопатки, Рр5 - распределение статического давления по корыту лопатки
Рис. 7. Сравнение профилей исходной и модифицированных лопаток Ш после решения 3Б-обратной задачи для трех сечений по высоте лопатки при п = 70 %: а - к = 90 %, периферия; б - Н = 60 %, среднее сечение; в - Н = 10 %, втулка
а
б
в
То, что поставленная цель достигнута, будет продемонстрировано при анализе течения при п = 70 % для окончательного варианта лопатки Ш после решения ЗБ-обратной задачи при п = 100 %. Анализ рис. 7 показывает, что модифицированная лопатка будет иметь некоторые проблемы при п = 100 %. Прежде всего из-за уменьшения горла проходных сечений лопатки в средних и периферийных сечениях уменьшится на несколько процентов величина расхода воздуха, что неприемлемо для реальных условий проектирования. Кроме того, возможен некоторый переразгон потока в окрестности задней кромки, что вызовет увеличение интенсивности замыкающего скачка уплотнения и, соответственно, снижение КПД колеса. Для устранения двух последних замечаний была решена ЗБ-обратная задача при п = 100 % в рабочие точки характеристики.
Оптимизация рабочего колеса осевой ступени на оборотах п = 100 %
Сначала решается прямая задача для сквозного расчета ЗБ-вязкого течения во всем компрессоре с модифицированным Ш на основе программного комплекса ЦИАМ ЗБ-1МР-МиЬТ1. Результаты этого расчета показаны на рис. 8.
Во всех сечениях можно отменить хорошее течение, определяющее высокий КПД этого колеса (на 1 % выше, чем у исходного варианта п = 100 %). Однако, как уже отмечалось, у этого варианта колеса имеется заметный недобор расхода, что приводит к снижению числа Маха для потока в Я1. Для устранения этого недостатка решается ЗБ-обратная задача для модификации Ш в рабочей точке при п = 100 %. Для задания модифицированной нагрузки теперь распределение статического давления в окрестности передней кромки профиля лопатки не изменяется, а изменяется нагрузка в сечениях, больших, чем 25 % длины хорды лопатки от передней кромки. В основном эти изменения давления касаются окрестности вблизи задней кромки лопатки, как это показано на рис. 9. На этом рисунке приведены не только распределения давления по спинке и корыту лопатки в сечении Ь по высоте лопатки, но также и в соседних сечениях с номерами Ь - 1 и Ь + 1.
На рис. 10 показан пример окончательной модификации Ш при п = 100 % для трех сечений лопатки - 10, 60 и 90 % высоты лопатки. В качестве положительного результата, полученного после решения ЗБ-обратной задачи, можно отметить увеличение горла проходных сечений в периферийной части лопатки Ш, что хорошо видно на рис. 10.
После ЗБ-обратной задачи полученное окончательное решение для Ш при п = 100 % было вставлено в трехступенчатый компрессор низкого давления. На основе программного комплекса ЗБ-1МР-МиЬТ1 был сделан сквозной расчет ЗБ-вязкого течения при п = 100 %. На рис. 11 показаны результаты расчета по прямой задаче в модифицированном варианте Ш при п = 100 %. Анализ полученных результатов показывает, что примерно на 0,1 увеличилось число Маха в средних сечениях по сравнению с вариантом на рис. 8. Это вызвано увеличением расхода воздуха через Ш на 0,З % по сравнению с исходным вариантом. Однако из-за увеличения числа Маха в средних сечениях Ш уровень КПД снизился на 1 % и вернулся к уровню исходного варианта.
После этого был проведен сквозной расчет ЗБ-вязкого течения в трехступенчатом КПД при п = 70 %. Результаты этого расчета приведены на рис. 12. Сразу можно отметить главный результат решения ЗБ-обратной задачи - отсутствие срыва потока с передних кромок Я1. Такое течение в Ш позволило увеличить запасы ГДУ на п = 70 % почти в 4 раза и увеличить КПД на 4,85 %.
Отсутствие срывных режимов обтекания передних кромок Ш в диапазоне п = 70.. .100 % дает возможность использовать ротор Ш, полученный на основе решения ЗБ-обратной задачи, без применения надроторного устройства.
Рис. 8. Промежуточный ротор (после решения ЗБ-обратной задачи, п = 70 %). Распределение чисел Маха в трех сечениях по высоте лопатки: к = 90 %, периферия (а); к = 60 %, среднее сечение (б) и к = 10 %, втулка (в) на поверхности вращения и в меридиональном сечении вдоль стороны разрежения (г) и вдоль стороны давления (д)
б
д
в
-0,0050271763 0,0897276255 0,188950525 0,2640397
0,0385025108 0.139333875 0,2283228
-0,0224752603 0,06538801 0,14948785 0,22526705
0,0177023508 0,111805875 0,18860355
Рис. 9. Пример задания модифицированной нагрузки Рис. 10. Сравнение профилей исходной и мо-в окрестности задней кромки для решения ЗБ-обратной дифицированной лопатки Ш после решения задачи при п = 100 % для трех сечений по высоте ЗБ-обратной задачи при п = 100 % для трех лопатки: а - к = 90 %, втулка; б - к = 60 %, среднее сечений по высоте лопатки: а - к = 90 %, втул-значение; в - к = 10 %, периферия ка; б - к = 60 %, среднее значение; в - к = 10 %,
периферия
б
б
в
в
б
а
в
д
Рис. 11. Окончательное решение ЗБ-обратной задачи при п = 100 %. Распределение чисел Маха в трех сечениях по высоте лопатки: к = 90 %, периферия (а); к = 60 %, среднее сечение (б) и к = 10 %, втулка (в) на поверхности вращения и в меридиональном сечении вдоль стороны разрежения (г),
и вдоль стороны давления (д)
г
г д
Рис. 12. Окончательное решение ЗБ-обратной задачи при п = 70 %, демонстрирующее отсутствие срыв-ного обтекания передних кромок лопатки Ш. Распределение чисел Маха в трех сечениях по высоте лопатки: к = 90 %, периферия (а); к = 60 %, среднее сечение (б) и к = 10 %, втулка (в) на поверхности вращения и в меридиональном сечении вдоль стороны разрежения (г) и вдоль стороны давления (д)
Выводы
ЗБ-обратная задача применена к проектированию ротора 1-й ступени КНД. Применение перепроектирования на 70 % оборотов позволило улучшить картины течения и интегральные газодинамические параметры.
Трехмерные расчеты всего трехступенчатого компрессора в целом показали, что КПД ротора 1-й ступени увеличился на 4,85 %, степень повышения давления возросла на 4 % на 70 % оборотов. При этом КПД на 100 % оборотов не уменьшился.
Улучшенное течение в роторе 1-й ступени дало четырехкратное увеличение запаса устойчивой работы на 70 % оборотов. КПД всего трехступенчатого компрессора увеличился на 5 %.
Отсутствие срыва на передней кромке ротора 1-й ступени в широком диапазоне оборотов от 70 до 100 % позволяет использовать этот ротор, спрофилированный по ЗБ-обратной задаче, без надроторного устройства.
Библиографический список
1. Meauze G. An inverse time marching method for the definition of cascade geometry // J. of Turbo-machinery. - 1982. - Vol. 11. - P. 650-656.
2. Dang Т., Damie S., Qiu X. Euler-based inverse method for turbomachine blades. P. 2. Three -demensional flows // AIAA J. - November 2000. - Vol. 38, no. 11. - P. 2007-2013.
3. Medd A.J., Dang T.Q., Larosiliere L.M. 3D inverse design loading strategy for transonic axial compressor blading, Larosiliere // ASME TURBO EXPO 2003, Atlanta, Georgia, USA. - 2003. - Paper no. GT2003-38501. - P. 489-495. - 7 p. DOI: 10.1115/GT2003-38501
4. Watanabe H., Zangeneh M. Design of the blade geometry of swept transonic fans by 3D inverse design // ASME TURBO EXPO 2003, Atlanta, Georgia, USA. - 2003. - Paper no. GT2003-38770. - P. 603-612. -10 p. DOI: 10.1115/GT2003-38770
5. New 3D inverse Navier - Stokes based method used to design turbomachinery blade rows / V.I. Mileshin, I.K. Orekhov, S.K. Shchipin, A.N. Startsev // ASME Heat Transfer/Fluids Engineering, Charlotte, USA. - 2004. - Paper no. HF-FED04-56436. - 9 p.
6. 3D viscous invers design of turbomachinery using one-equation turbulence model / Jinguang Yang, Xiaofang Wang, Yan Liu, Hu Wu // ASME Turbo Expo 2016: Turbomachinery Technical Conference and Exposition, Seoul, South Korea. - 2016. - Paper No. GT2016-56332. - P. V02CT39A014. - 9 p. DOI: 10.1115/GT2016-56332
7. 3D inverse design of transonic fan rotors efficient for a wide range of RPM / V.I. Mileshin, I.K. Orekhov, S.K. Shchipin, A.N. Startsev // ASME TURBO EXPO 2007 Power for Land, Sea and Air. Montreal, Canada. - 2007. - Paper No. GT2007-27817. - P. 341-352. - 12 p. DOI: 10.1115/GT2007-27817
8. Effect of tip clearance on flow structure and integral performances of six-stage HPC / V.I. Mileshin, I.K. Orekhov, V.A. Fateyev, S.K. Shchipin // Proc. of ISABE Conf., ISABE-2007-1179-Paper, Beijing, China. -2007. - 18 p.
9. Hield P. Semi-inverse design applied to an eight stage axial flow compressor // ASME TURBO EXPO 2008. - Paper No. GT2008-50430. - 2008. - P. 293-303. - 11 p. DOI: 10.1115/GT2008-50430
10. Qiu. X., Ji. M., Dang. T. Three-dimensional viscous inverse method for axial blade design // J. of Inverse Problems in Sci. and Eng. - 2009. - No. 17(8). - P. 1019-1036.
11. Experimental and numerical study of two first highly-loaded stages of compressors as a part of HPC and separate test unit / V.I. Mileshin, P.G. Kozhemyako, I.K. Orekhov, VA. Fateev // Proc. of 28th Inter. Cong. of Aeron. Sci. - 2012. - Paper no. ICAS2012-474. - P. 2675.
12. Arbabi A., Ghaly W., Medd Aerod A. Inverse blade design of axial compressors in three-dimensional flow using a commercial CFD program // Proc. of ASME Turbo Expo 2017, Charlotte, NC, USA. - 2017. -Paper No. GT2017-65194. - P. V02BT41A055. - 12 p. DOI: 10.1115/GT2017-65194
13. Arbabi A., Ghaly W. Inverse design of turbine and compressor stages using a commercial cfd program // Proc. of ASME Turbo Expo 2013. - 2013. - Paper No. GT2013-96017. - P. V06BT37A045. -14 p. DOI: 10.1115/GT2013-96017
14. Development of direct-driven and geared fan stages with reduced tip speeds / S.V. Pankov, V.I. Mileshin, I.K. Orekhov, V.A. Fateev // Proc. of ASME Turbo Expo Conference, Charlotte, USA. - 2017. -Paper No. GT2017-64585. - P. V02AT39A036. - 11 p. DOI: 10.1115/GT2017-64585
15. Численное решение многомерных задач газовой динамики / С.К. Годунов, А.В. Забродин, М.Я. Иванов, А.Н. Крайко, Г.П. Прокопов. - М.: Наука, 1976. - 400 c.
16. Милешин В.И. Расчет сверхзвукового пространственного обтекания воздухозаборников на режимах с выбитой ударной волной // Журнал вычислительной математики и математической физики. -1986. - Т. 26, № 11. - C. 1704-1718.
17. Stall margin improvement in three-stage low pressure compressor by use of slot type casing treatments / F.Sh. Gelmedov, V.I. Mileshin, P.G. Kozhemyako, I.K. Orekhov // Proc. of ASME Turbo 2014. - 2014. -Paper No. GT2014-26298. - P. V02AT37A036. - 11 p. DOI: 10.1115/GT2014-26298
18. Baldwin B.W., Lomax H. Thin layer approximation and algebraic model for separated turbulent flows // AIAA. - 1978. - Paper No. 78-257. - 9 р.
19. Gouskov O.V., Kopchenov V.I., Nikiforov D.A. Flow numerical simulation in the propulsion elements of aviation space system within full Navier - Stokes equations // Int. Conf. on the Methods of Aeroph. Research. Proc. - Novosibirsk, 1994. - Part 1. - P. 104-109.
20. Иванов М.Я., Крупа В.Г., Нигматуллин Р.З. Неявная схема С.К. Годунова повышенной точности для интегрирования уравнений Навье - Стокса // Журнал вычислительной математики и математической физики. - 1989. - Т. 29. - C. 888-901.
Reference
1. Meauze, G. An Inverse Time Marching Method for the Definition of Cascade Geometry // Journal of Turbomachinery, 1982, Vol. 11, pp. 650-656.
2. Т. Dang, S. Damie, and X. Qiu. Euler-Based Inverse Method for Turbomachine Blades. Part 2: Three-Demensional Flows // AIAA Journal, Vol. 38, No. 11, November 2000, pp. 2007-2013.
3. Medd, A.J., Dang, T.Q. and Larosiliere, L. M. 3D inverse design loading strategy for transonic axial compressor blading, Larosiliere // ASME TURBO EXPO 2003, Atlanta, Georgia, USA, 2003, Paper no. GT2003-38501, pp. 489-495, 7 p. DOI: 10.1115/GT2003-38501
4. Watanabe, H. and Zangeneh, M. Design of the blade geometry of swept transonic fans by 3D inverse design // ASME TURBO EXPO 2003, Atlanta, Georgia, USA, 2003, Paper no. GT2003-38770, pp. 603-612, 10 p. DOI: 10.1115/GT2003-38770
5. Mileshin, V.I., Orekhov, I.K., Shchipin, S.K., and Startsev, A.N. New 3D Inverse Navier-Stokes Based Method Used to Design Turbomachinery Blade Rows // ASME Heat Transfer/Fluids Engineering, Charlotte, USA, 2004, Paper no. HF-FED04-56436, 9 p.
6. Jinguang Yang, Xiaofang Wang, Yan Liu, Hu Wu. 3D Viscous Invers Design of Turbomachinery Using One-Equation Turbulence Model // ASME Turbo Expo 2016: Turbomachinery Technical Conference and Exposition, Seoul, South Korea, 2016 Paper No. GT2016-56332, P. V02CT39A014, 9 p. DOI: 10.1115/GT2016-56332
7. Mileshin, V.I., Orekhov, I.K., Shchipin, S.K., and Startsev, A.N. 3D inverse design of transonic fan rotors efficient for a wide range of RPM // ASME TURBO EXPO 2007 Power for Land, Sea and Air. Montreal, Canada, 2007, Paper No. GT2007-27817, P. 341-352, 12 p. DOI: 10.1115/GT2007-27817
8. Mileshin V.I., I.K. Orekhov, V.A. Fateyev, S.K. Shchipin. Effect of tip clearance on flow structure and integral performances of six-stage HPC // Proceedings of ISABE Conference, ISABE-2007-1179-Paper, Beijing, China, 2007, 18 p.
9. Hield, P., 2008. Semi-inverse design applied to an eight stage axial flow compressor // ASME TURBO EXPO 2008, Paper No. GT2008-50430, pp. 293-303, 11 p. DOI: 10.1115/GT2008-50430
10. Qiu. X., Ji. M., and Dang. T., 2009. Three-dimensional viscous inverse method for axial blade design // Journal of Inverse Problems in Science and Engineering, 2009, no. 17(8), pp. 1019-1036.
11. V.I. Mileshin, P.G. Kozhemyako, I.K. Orekhov, VA. Fateev. Experimental and numerical study of two first highly-loaded stages of compressors as a part of HPC and separate test unit // Proceedings of 28th International Congress of Aeronautical Sciences, 2012, paper no. ICAS2012-474, pp. 2675.
12. A. Arbabi, W. Ghaly, A. Medd Aerod. Inverse Blade Design of axial compressors in three-dimensional Flow using a commercial CFD Program // Proceedings of ASME Turbo Expo 2017, Charlotte, NC, USA, Paper No. GT2017-65194, P. V02BT41A055, 12 p. DOI: 10.1115/GT2017-65194
13. A. Arbabi, and W. Ghaly, 2013. Inverse design of turbine and compressor stages using a commercial cfd program // Proceedings of ASME Turbo Expo 2013, 2013, Paper No. GT2013-96017, P. V06BT37A045, 14 p. DOI: 10.1115/GT2013-96017
14. S.V. Pankov, V.I. Mileshin, I.K. Orekhov, V.A. Fateev. Development of direct-driven and geared fan stages with reduced tip speeds // Proceedings of ASME Turbo Expo Conference, Charlotte, USA, 2017, Paper No. GT2017-64585, P. V02AT39A036, 11 p. DOI: 10.1115/GT2017-64585
15. S.K. Godunov, A.V. Zabrodin, M.Ya. Ivanov, A.N. Kraiko, G.P. Prokopov. Chislennoye resheniye mnogomernykh zadach gazovoy dinamiki [Numerical Solution of Multi-Dimensional Problems Gas Dynamics]. Moscow: Nauka, 1976, 400 p.
16. V.I. Mileshin. Raschet sverkhzvukovogo prostranstvennogo obtekaniya vozdukhozabornikov na rezhimakh s vybitoy udarnoy volnoy [Computation of 3D Supersonic Flow Engine Intel at the Detached Shock Flow Regime]. Journal of Computational Mathematics and Mathematical Physics, 1986, vol. 26, no. 11, pp. 1704-1718.
17. F.Sh. Gelmedov, V.I. Mileshin, P.G. Kozhemyako, I.K. Orekhov. Stall margin improvement in three-stage low pressure compressor by use of slot type casing treatments. Proceedings of ASME Turbo 2014, 2014, Paper No. GT2014-26298, P. V02AT37A036, 11 p. DOI: 10.1115/GT2014-26298
18. Baldwin B.W. and Lomax H. Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows // AIAA Paper 78-257, 1978, 9 р.
19. O.V. Gouskov, V.I. Kopchenov, D.A. Nikiforov. Flow numerical simulation in the propulsion elements of aviation space system within full Navier-Stokes equations. International Conference on the Methods of Aerophysical Research. Proceedings. Part 1. Novosibirsk, 1994, pp. 104-109.
20. Ivanov M.Ya., Krupa V.G., Nigmatullin R.Z.. Neyavnaya skhema S.K. Godunova povyshennoy tochnosti dlya integrirovaniya uravneniy Navye-Stoksa [Implicit S.K. Godunov's scheme of Raised Accuracy for Integration of Navier-Stokes Equations]. Journal of Computational Mathematics and Mathematical Physics, 1989, vol. 29, pp. 888-901.
Об авторах
Милешин Виктор Иванович (Москва, Россия) - кандидат физико-математических наук, начальник отделения ФГУП ЦИАМ им. П.И. Баранова (111116, г. Москва, ул. Авиамоторная, д. 2; е-mail: mileshin@ciam.ru).
Орехов Игорь Константинович (Москва, Россия) - старший научный сотрудник ФГУП ЦИАМ им. П.И. Баранова (111116, г. Москва, ул. Авиамоторная, д. 2; е-mail: orekhov@ciam.ru).
Щипин Сергей Константинович (Москва, Россия) - кандидат технических наук, заместитель главного конструктора АО «Российская самолетостроительная корпорпция "МиГ"» (125284, г. Москва, 1-й Боткинский проезд, д. 7; е-mail: Sershchipin@mail.ru).
About the authors
Viktor I. Mileshin (Moscow, Russian Federation) - CSc in Physical and Mathematical Sciences, Head of Division, Central Institute of Aviation Motors named after P.I. Baranov (CIAM) (2, Aviamotornaya st., Moscow, 111116, Russian Federation; e-mail: mileshin@ciam.ru).
Igor K. Orekhov (Moscow, Russian Federation) - Senior Researcher, Central Institute of Aviation Motors named after P.I. Baranov (CIAM) (2, Aviamotornaya st., Moscow, 111116, Russian Federation; e-mail: orekhov@ciam.ru).
Sergey K. Shchipin (Moscow, Russian Federation) - CSc in Technical Sciences, Deputy Chief Designer, Russian Aircraft Corporation "MiG" (7, 1-st Botkinsky drive, Moscow, 125284, Russian Federation; e-mail: Sershchipin@mail.ru).
Получено 23.04.2019
