CC BY
4«Труды МАИ». Выпуск № 82
www.mai.ru/science/trudy/
УДК 629.735.45
Методические особенности численного моделирования в рамках сеточных методов поля течения около несущего винта на режиме висения с учетом вихревой структуры
Вершков В.А.*, Воронич И.В.**, Вышинский В.В.***
Московский физико-технический институт (государственный университет), Институтский пер., 9, Долгопрудный, 141700, Россия *e-mail: vershkov. va@gmail. com **e-mail: voronich. iv@mipt. ru ***e-mail: vyshinsky@falt.ru
Аннотация
В работе рассмотрена методика численного расчета обтекания несущего винта на режиме висения в рамках комплекса Numeca FINE/Turbo. Для ее валидации использованы данные тестовых экспериментов Карадонны-Танга для жесткого двухлопастного винта. На основе обзора работ сделана попытка систематизации факторов, определяющих точность и возможности аналогичных методик. В расчетах использовались генерируемые в полуавтоматическом режиме структурированные многоблочные расчетные сетки, адаптированные в пограничном слое и в области вихревой системы. Для оценки индустриальной пригодности методик предложен критерий, основанный на погрешности поверхностных распределений коэффициента давления. Полученные численные результаты удовлетворительно соответствуют экспериментальным данным.
Ключевые слова: несущий винт, аэродинамические характеристики, численное моделирование, Numeca FINE/Turbo.
Список обозначений b - хорда лопасти в сечении на радиусе r [м] С - число Куранта Cp - коэффициент давления
\ÔCp\=\Cpрасчет - Срэксперимент\ - погрешность коэффициента давления
Cy - коэффициент нормальной силы
n - частота вращения винта [оборот/мин]
p, p0 - статическое, полное давление [Па]
R - радиус несущего винта [м]
r - радиус сечения лопасти несущего винта [м]
ReK - число Рейнольдса в концевом сечении
T, T0 - статическая, полная температура [К]
VK - скорость конца лопасти [м/с]
ф0 - угол установки лопасти (угол "общего шага") [град.] Y - угловая длина концевого вихря [град.]
Введение
Задача расчета аэродинамических характеристик несущего винта является основной в вертолетостроении. Построению вихревых теорий винта и методов расчета, основанных на дискретных вихревых моделях, посвящен ряд основополагающих трудов [1-4] и работ, дополняющих и развивающих базовые модели с использованием эмпирических данных [5-7]. До настоящего времени методы этого класса и по-
2
строенные на их основе программные комплексы [7, 8] заслуженно являются ведущими для практических приложений. Ограничениями применимости дискретных вихревых моделей являются предположения о существенно дозвуковом и безотрывном режиме обтекания и гипотеза плоских сечений при использовании экспериментальных аэродинамических характеристик профилей.
На первом этапе применения методы расчета обтекания несущего винта в рамках подхода вычислительной аэродинамики сталкивались с необходимостью предварительного задания формы вихревой пелены и влекли существенные вычислительные затраты [9]. Применение численных моделей на основе уравнений Навье-Стокса и Рейнольдса на рубеже 1980-1990 гг. по мере роста производительности вычислительной техники дало возможность проводить расчет поля течения и вихревой системы как его части «из первых принципов». До 2000 г. такой подход к решению задач обтекания несущего винта не был вполне успешным из-за недостаточной точности моделирования вихревой пелены и концевых вихрей [10]. Современный уровень инструментов вычислительной аэродинамики позволяет рассчитывать на прогресс в решении этой проблемы [10-16].
В настоящей работе рассматривается методика численного расчета обтекания несущего винта на режиме висения в рамках Numeca FINE/Turbo. Необходимыми являются верификация, проверка сходимости и корректности получаемых решений, и валидация, проверка в задачах, имеющих достоверные и достаточно подробные решения. К таким задачам относятся тестовые эксперименты Карадонны-Танга [17] для жесткого двухлопастного винта радиусом R=1.143 м с профилем NACA0012.
Наряду с интегральными характеристиками несущего винта детальное представление о поле течения дают локальные аэродинамические характеристики, прежде всего, распределения коэффициента давления Ср на поверхности лопасти. Для получения корректного решения необходимо разрешение глобальной вихревой системы в рамках самосогласованного поля течения. Представляется необходимым формирование простого критерия, позволяющего судить об индустриальной пригодности методики. В качестве такого критерия в настоящей работе предложен уровень погрешности распределений коэффициента давления на лопасти: |^Ср|<0.05.
Анализ работ показал, что для применения в данной области расчетная методика в рамках полных уравнений Навье-Стокса и Рейнольдса должна включать модели вращающихся и неподвижных блоков с их корректным сопряжением, набор гибких граничных условий с фиксацией полных или статических параметров и экстраполяцией поля скорости [10-14]. Эффективность методики определяется возможностями автоматизированного построения адаптированных к вихревой системе расчетных сеток и ускорения расчета.
Работа имеет целью систематизацию факторов, определяющих точность численного решения задачи обтекания несущего винта, и тестирование расчетной методики, имеющей потенциал к достижению индустриальной пригодности.
1. Краткий обзор работ по сеточным методам В работах [10-11] представлен обзор методов, рассматриваемых на момент 2000 г. как перспективные для расчета несущих винтов. Существенно, что наиболее важным фактором признана адаптация расчетных сеток. Повышение порядка точности разностных схем также рассматривается как необходимое направление разви-
тия. Эти выводы актуальны по настоящее время. Работа [12] подчеркивает роль технологии построения и начальной адаптации расчетных сеток, в том числе с учетом последующего взаимного движения блоков в процессе расчета. В этом отношении одни авторы считают технологию перекрывающихся расчетных сеток (overset grids) наиболее универсальной для адаптации к вихревой системе и учета движения лопасти [10,13], другие считают технологию деформируемых расчетных сеток (deformed grids) предпочтительной [12]. Общим является использование адаптированных к геометрии (body fitted) пристенных блоков, при этом для способа их стыковки с внешней сеткой имеются упомянутые расхождения в подходах. Успешные примеры применения решателей, использующих декартовы сетки, также показывают необходимость построения адаптированных к поверхности тела сеточных слоев.
Типичной для текущего состояния общеупотребительных методов является методика, изложенная в работе [14] и посвященная валидации кода SPARC. Алгоритмы этого кода основаны на методе конечного объема в рамках блочно-структурированных расчетных сеток, численный метод - на схеме центральных разностей 2-го порядка точности с искусственной вязкостью по методу SLIP и явном многошаговом методе Рунге-Кутты интегрирования по времени. Программа решает полные осредненные по Рейнольдсу уравнения Навье-Стокса с замыканием по различным моделям турбулентности. В расчетах использовалась низкорейнольдсовая модель турбулентности Спаларта-Аллмараса, которая доказала свою эффективность для подобных задач [10, 15]. Для ускорения сходимости применялись технологии локальных шагов по времени и неявного сглаживания невязки, которые включены в явную схему. В коде реализован многосеточный подход с V-циклами, который дает
дополнительное ускорение сходимости. Расчетная область была построена на основе С-топологии вокруг хорды лопасти и Н-топологии вдоль ее размаха. В отличие от О-И топологии, С-Н топология успешнее разрешает вихревую систему винта. Течение на режиме висения является стационарным во вращающейся системе координат и пространственно периодическим, поэтому расчетная область строится для одной лопасти. Адаптированная структурированная сетка была создана в приложении IGG (Кишееа) в полуавтоматическом режиме с использованием языка сценариев. Полученная сетка включала 5.7 миллионов ячеек на периодическую область с разрешением пограничного слоя на уровне у =1 на всем размахе лопасти и сгущениями в областях пограничного слоя, следа, концевых вихрей.
Граничные условия в [14] были следующие: на поверхности лопасти использовалось условие прилипания с нулевым тепловым потоком, на поверхности вала -условие непротекания (Эйлеровская стенка), на внешних границах - комбинированное условие входа/выхода с фиксацией полного давления для втекающего потока и статического давления для вытекающего, на поверхностях периодичности - условие периодичности. Помимо этого на границах втекания фиксировалось направление скорости по нормали к границе и отношение турбулентной к ламинарной вязкости, равное 1. При моделировании процесса использовался шаг по времени, соответствующий 1/62 оборота (5.8о за итерацию).
Полученные в [14] распределения коэффициента давления Ср показывают разумное соответствие с экспериментом. Погрешность в расчете коэффициента тяги винта составила ~1%. Проблемными для разрешения областями являются области
разгона потока на верхней и нижней поверхностях лопасти и область замыкающего сверхзвуковую зону скачка. В этих зонах критерий точности |^Ср|<0.05 не выполняется. Коэффициент нормальной силы Cy в сечениях лопасти согласуется с данными эксперимента при использования одинакового метода интегрирования распределений давления. Радиальное и вертикальное положение концевого вихря хорошо согласуется с данными эксперимента вплоть до углау=450 Это подтверждает возможность получения корректных результатов в рамках описанной расчетной методики при использовании адаптированной расчетной сетки.
Основной вывод по рассмотренным работам заключается в приоритете адаптации расчетной сетки в области вихрей. Разрешение концевых вихрей на протяженности одного оборота оказывается достаточным для получения достоверных значений распределений коэффициента давления и интегральных характеристик на режиме висения. По данным работы [10] очевидно, что начальная адаптация расчетной сетки не позволяет гарантировать разрешение особенностей для режимов с боковым обдувом. Решением этой проблемы является динамическая адаптация расчетной сетки к получаемому решению [10], обеспечивающая не менее 10 ячеек на диаметр вихря. Такая технология в настоящее время не вполне устоялась и имеет свои варианты для перекрывающихся и деформируемых сеток. Не вызывает сомнений вывод о том, что численное решение менее зависимо от подробности сетки при использовании методов 3-го и выше порядка пространственной аппроксимации [10,11]. Потенциал методов высокого порядка может быть использован в полной мере только на адаптированных и достаточно равномерных сетках, так как точность
решения падает на грубых и сильно неравномерных сетках. Отмеченное условие часто не выполняется в индустриальных расчетах, которые при этом требуют устойчивости и предсказуемости алгоритмов моделирования. Методы конечного элемента теоретически слабо зависят от качества сетки и позволяют длительно сохранять без искажения завихренность как основную переменную [11]. Тем не менее, они также подчинены указанному ограничению разрешения вихря и имеют проблемы в реализации равномерного порядка аппроксимации для конвективных, диффузионных и источниковых членов, особенно при нестационарном моделировании.
2. Постановка задачи В настоящей работе рассматривается один режим тестового эксперимента [17], служащего источником данных для тестирования и валидации программ, применяемых для расчета несущих винтов.
В экспериментах исследовались аэродинамические характеристики жесткого винта с лопастями прямоугольной формы в плане на основе профиля КЛСЛ0012, без геометрической крутки. Радиус винта ^=1.143 м, хорда лопасти ¿=0.1905 м. Угол установки лопасти (общий шаг винта) ф0=8°. В настоящей работе рассматривается дозвуковой режим вращения лопастей с частотой вращения «=1750 об/мин с соответствующими скоростью в концевом сечении лопасти Кк=209 м/с и числом Рейнольдса Reк=2 880 000. Схема установки приведена на рисунке 1, взятом без изменений из работы [17].
Рисунок 1 Схема экспериментальной установки [17].
Установка работала в помещении, где обеспечивалось отсутствие паразитных потоков за счет отсоса воздуха. Это позволяет рассматривать вычислительную задачу без учета влияния внешних границ - винт в бесконечном пространстве. Стоит обратить внимание на наличие шайбы в корневом сечении лопасти для предотвращения образования комлевого вихря.
В эксперименте измерялись распределения коэффициента давления ^ на поверхности лопастей в сечениях r/R=0.5^; 0.68; 0.8; 0.89; 0.96. Помимо этого проводились измерения поля скорости в области под винтом для идентификации положения
концевых вихрей.
3. Расчетная модель
В качестве модели, описывающей движение газа, принята система уравнений Навье-Стокса, осредненная по Рейнольдсу, с замыканием по моделям турбулентно-
сти Спаларта-Алмараса и SST. Вся расчетная область, содержащая лопасть винта,
рассматривается во вращающейся системе отсчета с учетом членов, отвечающих за действие центробежной и кориолисовой сил. Рабочей средой является воздух, рассматриваемый как вязкий сжимаемый идеальный газ, коэффициенты динамической вязкости и теплопроводности которого полагаются зависящими от температуры по закону Сазерленда с опорными значениями, взятыми при нормальных условиях.
В работе использовались генерируемые в полуавтоматическом режиме структурированные расчетные сетки Н-типа с О-слоем вокруг лопастей, рисунок 2. О-слой был локализован так, чтобы он вносил минимальные искажения в решение при сходе пограничного слоя с задней кромки. Особенностью получившейся сетки является достаточно малая толщина О-слоя, вследствие чего на расстоянии от лопасти порядка 10-4 м ячейки уже ориентированы по потоку, рисунок 3.
Сгущение узлов сетки осуществлялось в областях пограничного слоя, следа, концевых вихрей. Количество элементов в них значительно варьировалось, для расчета была использована расчетная сетка с 10 млн. узлов на периодическую область. Разрешение сеткой пограничного слоя обеспечивалось на уровне у+<1 на всем размахе лопасти (в среднем у+~0.1).
Рисунок 2 Структура сетки в сечении
лопасти.
Рисунок З Структура сетки около задней кромки лопасти.
Для моделирования с помощью FINE/Turbo использовались следующие граничные условия. На входе в расчетную область (рисунок 4) втекает втягиваемый винтом поток. Фиксировались значения полного давления p0 = 103 027 Па, полной температуры T0 = 2S9.75 K и параметров турбулентности - отношения турбулентной и ламинарной вязкости на уровне 1. На выходе из расчетной области фиксировалось ос-редненное по поверхности статическое давление, равное невозмущенному значению. Площадь этой границы должна быть достаточно большой, чтобы граничное условие не влияло на струю. Внешняя граница расчетной области может быть определена как Эйлеровская стенка с условием непротекания, если она расположена достаточно далеко. Аналогичное условие ставится на поверхности вала. На поверхностях лопастей использовалось условие прилипания с нулевым тепловым потоком.
Ввиду характера течения используются периодические граничные условия для связи сторон расчетной области, которые совмещаются при вращении вокруг оси на угол 180°.
Рисунок 4 Расчетная сетка.
При построении расчетной сетки учитывается, чтобы узлы расчетных сеток на различных сторонах периодической границы совпадали при вращении на угол периодичности. Поле течения для режима висения является нестационарным в неподвижной системе координат, тогда как во вращающейся системе его можно рассматривать как стационарное при отсутствии неустойчивостей, связанных с взаимодействием вихрей. В FINE/Turbo такой подход позволяет получить корректное решение с использованием методов ускорения сходимости.
В решателе реализована многосеточная технология ускорения расчета за счет использования вложенных сеток различной подробности с переходом между уровнями. В качестве основного численного метода для стационарных задач используется схема центральных разностей с искусственной вязкостью с применением локального шага по времени. Для гарантии устойчивости решения используется ограничение на число Куранта С<1. Для проверки полученных результатов проводились расчеты в нестационарной постановке задачи в рамках метода дуального шага по времени и с использованием противопотоковых разностных схем.
Начальные условия выбирались в соответствии с параметрами невозмущенного потока при наличии малой осевой скорости для корректной работы условий входа
и выхода.
В стационарной постановке задачи расчет велся на протяжении 10 000 итераций. Критериями сходимости численного решения являются значения невязок (скорости изменения параметров поля течения) и интегральные характеристики. В расчетах использованы критерии сходимости для тяги винта |£У|~2 Н и крутящего момента |£МК|~0.5 Н*м. После этого расчет проводился в нестационарной постановке на протяжении 5-10 периодов вращения с шагом по времени, соответствующим по-
вороту лопасти на угол 2°.
4. Анализ результатов
Сопоставление полученных численных решений с экспериментальными данными показало их удовлетворительное соответствие. На рисунке 5 представлена общая картина линий тока в плоскости, перпендикулярной плоскости винта. Их структура соответствует теоретическим представлениям.
На рисунке 6 показаны распределения коэффициента давления Ср в различных сечениях лопасти в сравнении с экспериментальными значениями. Критерий точности |^Ср|<0.05 выполняется неравномерно: в сечении г/Я=0.8 на верхней поверхности лопасти наблюдается расхождение. На рисунке 7 представлено распределение коэффициента нормальной силы Су вдоль размаха лопасти.
Рисунок 5 Структура течения Расхождения наблюдаются в концевых сечениях около лопасти. г/К=°.8 и °.89.
Ср -1.2
1
Г К,
>
>
1/ —
/ > } N
1 1
1
!
Распределений Ср в сечен им 0.5
0 ф ф Эксперимент -|---1---1- Расчет
1
а)
Ср -1-е
В)
Ср -1.5
х/Ь
л.
к.
к
л-
—/ / / /
}
Распределение Ср о сечении 0.89 Ф ф ф Эксперимент -{----{- Расчет
1
х/Ь
Ср И.5
А.
и К
1 I
/ г*- -4—
(
1
! * |
Распределение Ср в сече мии 0.8
0 9 9 Эксперимент + + + Расчет
0 г с 4 б) 06 0 8 X/
1 *
N ч
* 1 >1
} — I
I у-*-1......
/ N
/
1
Ряситредепений Ср н оечпнии ОЛМЗ Ф 0 ф Эксперимент ---Расчет
1 1 '
ог о.а
х/Ь
Г)
Рисунок 6 Распределения Ср в сечениях.
Зависимости вертикальной и радиальной координат концевого вихря от угла поворота у представлены на рисунке 8 в сравнении с экспериментальными координатами [17] и моделью [18]. Ошибки координат центра вихря, полученного в расчете, определялись следующим образом. Для каждого угла у было построено уникальное изображение, на котором изолировался нужный фрагмент вихря. Зона с наи-
большей интенсивностью имеет конечный размер, который и принимался за диапазон ошибки при определении положения центра концевого вихря. В полученных результатах спуск и поджатие вихря соответствуют модели [18] до угла у=360°.
Рисунок 7 Распределение Су.
Рисунок 8 Координаты вихря.
Типичная форма поверхности равной завихренности изображена на рисунке 9. Концевые линии тока во вращающейся системе представлены на рисунке 10. Точность моделирования концевых вихрей сказывается в большей степени на концевых сегментах лопасти [14], рисунок 7.
Рисунок 9 Изоповерхность завихренности.
15
Рисунок 10 Концевые линии тока.
Использование модели турбулентности SST не дало принципиальных отличий в численных решениях [10, 15]. При достаточном разрешении поля течения улучшение результата может достигаться за счет моделей турбулентности с поправками на кривизну линий тока.
Следует отметить, что в экспериментальных исследованиях встречаются режимы, при которых число Рейнольдса в сечениях лопасти не превышает 1 млн. Это делает необходимым использование расчетных моделей, основанных на подходе прямого моделирования в рамках уравнений Навье-Стокса или на подходах моделирования ламинарно-турбулентного перехода [19]. Оба варианта являются ресурсоемкими, так как требуют нестационарной постановки задачи и подробной сетки [20].
Полученные результаты согласуются с выводами о роли адаптации расчетной сетки и возможности получения корректных результатов с помощью методов 2-го порядка аппроксимации.
Выводы
Для оценки индустриальной пригодности методик численного расчета обтекания несущего винта предложен критерий, основанный на погрешности поверхностных распределений коэффициента давления: (¿^<0.05. Данный критерий соответствует типовой погрешности экспериментальных данных и обеспечивает разумную точность интегральных характеристик. Полученные численные результаты в основном соответствуют этому критерию. Превышение заданного уровня погрешности ^ в концевых сечениях обуславливает отличие значений коэффициента нормальной силы ^ Анализ вертикальной и радиальной координат концевого вихря показал,
что спуск и поджатие вихря соответствуют модели [18] до угловой протяженности у=360°.
Основным фактором, определяющим точность, является адаптация расчетной сетки в области вихрей, предпочтительно динамическая. Для снижения сеточной зависимости численных решений необходимо повышение точности разностной схемы.
Работа выполнена при финансовой поддержке и в интересах проекта "Разработка программно-аппаратного комплекса реалистичного восприятия летчиком сложных режимов полета и оценки его психофизиологического состояния" (Договор № 02.С25.31.0017 между ОАО «РСК «МиГ» и Министерством образования и науки РФ об условиях предоставления и использования субсидии на реализацию комплексного проекта по созданию высокотехнологичного производства, выполняемого с участием МФТИ) с целью валидации методик моделирования вихревых структур, используемых в проекте.
Библиографический список
1. Теория несущего винта / Под ред. А.К. Мартынова. - М.: Машиностроение, 1973. -364 с.
2. Джонсон У. Теория вертолета. Т.1, 2. - М.: Мир, 1983.- 503 с., - 529 с.
3. Белоцерковский С.М., Локтев Б.Е., Ништ М.И. Исследование на ЭВМ аэродинамических и аэроупругих характеристик винтов вертолетов. - М.: Машиностроение, 1992. - 218с.
4. Leishman J.G. Principles of Helicopter Aerodynamics 2nd Edition. New York: Cambridge University Press, 2008.
5. Игнаткин Ю.М., Макеев П.В., Гревцов Б.С., Шомов А.И. Нелинейная лопастная вихревая теория винта и ее приложения для расчета аэродинамических характеристик несущих и рулевых винтов вертолета // Вестник Московского авиационного института. 2009. Т. 16. № 5. С. 24-31.
6. Головкин В.А., Миргазов Р.М. Метод расчета аэродинамических характеристик крыла и несущего винта на основе обратной процедуры использования "гипотезы плоских сечений" // Научный вестник МГТУ ГА. 2010. № 154. С. 34-41.
7. Игнаткин Ю.М., Макеев П.В., Шомов А.И. Программный комплекс для расчета аэродинамических характеристик несущих и рулевых винтов вертолетов на базе нелинейной лопастной вихревой теории // Электронный журнал «Труды МАИ», 2010, выпуск №38: http : //www. mai. ru/science/trudy/published. php?ID= 14148 (дата публикации 23.06.2010).
8. Головкин М.А., Кочиш С.И., Крицкий Б.С. Методика расчета аэродинамических характеристик комбинированной несущей системы летательного аппарата // Электронный журнал «Труды МАИ», 2012, выпуск №55: http://www.mai.ru/science/trudy/published.php?ID=30023 (дата публикации 16.05.2012).
9. Ramachandran K., Tung C., Caradonna F.X. Rotor Hover Performance Prediction Using A Free Wake Computational Fluid Dynamics Method // Journal of Aircraft. 1989. Vol. 26. No. 12. P. 1105-1110.
10. Hariharan N., Sankar L. A Review of Computational Techniques for Rotor Wake Modeling // AIAA Paper 2000-0114.
11. Boelens O.J., van der Ven H., Oskam B. and Hassan A.A. Accurate and Efficient Vortex-Capturing for a Helicopter Rotor in Hover // NLR-TP-2000-420.
12. Steijl R., Barakos G., Badcock K. A Framework for CFD Analysis of Helicopter Rotors in Hover and Forward Flight // Int. Journal for Numerical Methods in Fluids. 2006. V. 51. P. 819-847.
13. Dong-Kyun Im, Seong-Yong Wie, Kim Eugene et al. Aerodynamic Analysis of Rotor Blades using Overset Grid with Parallel Computation // Lecture Notes in Computational Science and Engineering. V. 74. Parallel CFD 2008. P. 101-110.
14. Doerffer P., Szulc O. Numerical Simulation of Model Helicopter Rotor in Hover // TASK Quarterly. 2008. V. 12. N 3. P. 227-236.
15. Игнаткин Ю.М., Константинов С.Г. Исследование аэродинамических характеристик профиля и законцовок лопасти несущего винта вертолета методами CFD // Электронный журнал «Труды МАИ», 2012, выпуск 57: http://www.mai.ru/science/trudy/published.php?ID=30874 (дата публикации 30.06.2012).
16. Игнаткин Ю.М., Константинов С.Г. Исследование аэродинамических характеристик несущего винта вертолета методом CFD // Электронный журнал «Труды МАИ», 2012, выпуск №57:
19
http://www.mai.ru/science/trudy/published.php?ID=30875 (дата публикации 30.06.2012).
17. Caradonna F.X., Tung C. Experimental and Analytical Studies of a Model Helicopter Rotor in Hover // NASA TM 81232. 1981. 60 p.
18. Kocurek J.D., Tangler J.L. A Prescribed Wake Lifting Surface Hover Performance Analysis // 32nd Annual National Forum of the American Helicopter Society. 1976. Preprint 1001.
19. Garcia A.J., Barakos G.N. Hover Predictions on the S-76 Rotor using HMB2 // 53rd AIAA Aerospace Sciences Meeting. 2015.
20. Hirsch Ch. Numerical Computation of Internal and External Flows. Vol. 1, 2. John Wiley & Sons, 1994.
