CC BY
91
Том XLV
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
2014
№ 2
УДК 532.517:4
МОДЕЛИРОВАНИЕ ПЕРИОДИЧЕСКИХ ВИХРЕВЫХ СТРУКТУР
В СЛЕДЕ ЗА ПРОФИЛЕМ
П. А. БАРАНОВ, С. В. ГУВЕРНЮК, С. А. ИСАЕВ, А. Г. СУДАКОВ, А. Е. УСАЧОВ
Параметрические расчеты нестационарного обтекания профиля NACA 0012 при фиксированном числе Рейнольдса 40 000 выполнены на нескольких многоблочных пересекающихся сетках разного масштаба и густоты, охватывающих совместно ближний и практически весь дальний след вплоть до начального разгонного вихря. Решения, полученные с помощью полуэмпирических моделей, сравниваются с имеющимися оценочными и экспериментальными данными. Верификация двумерной модели выполнена при сравнении численных прогнозов поперечного обтекания толстой пластинки с экспериментальными данными Игараши.
Ключевые слова: численное моделирование, турбулентность, аэродинамический профиль, процедура коррекции давления, многоблочные вычислительные технологии, пакет VP2/3.
ВВЕДЕНИЕ
Нестационарные вихревые поля являются источником акустических возмущений. В данной работе численно анализируются периодические гидродинамические процессы генерирования и распространения крупномасштабных вихревых структур в ближнем и дальнем турбулентном следе за профилем NACA 0012 при варьировании углом атаки от 0 до 180° (круговая обдувка).
Проблематика, связанная с обтеканием профилей, по-прежнему привлекает внимание аэродинамиков. В последние годы интерес к ней проявлялся при численном моделировании трансзвукового нестацио-
ГУВЕРНЮК Сергей Владимирович
кандидат физико-математических наук, заместитель директора ИМех МГУ
ИСАЕВ Сергей Александрович
доктор физико-математических наук, профессор СПбГУГА
СУДАКОВ Александр Григорьевич
доктор технических наук, главный конструктор ОАО АК «Ригель»
БАРАНОВ Павел Андреевич
научный сотрудник ОАО АК «Ригель»
УСАЧОВ Александр Евгеньевич
кандидат физико-математических наук, ведущий научный сотрудник ЦАГИ
нарного режима течения [1, 2], а также при постановке экспериментов для анализа возможностей управления обтеканием с помощью выдува струи [3].
Известно существенное количественное отличие численных решений двумерных нестационарных осредненных по Рейнольдсу уравнений Навье — Стокса (URANS), замкнутых с помощью полуэмпирических моделей турбулентности, от теоретико-экспериментальных оценок осредненных интегральных характеристик профиля NACA 0012, приведенных в [4, 5]. Оно послужило основанием для вывода в [6] о неприемлемости подхода Рейнольдса к моделированию нестационарных отрывных течений и преимуществе альтернативного подхода на базе модели отсоединенных вихрей. Однако двумерный подход URANS оказался вполне приемлемым для моделирования нестационарного обтекания кругового цилиндра в до- и закритических режимах [7—10], а также в отдельных расчетах профиля NACA 0015 [11]. Чтобы прояснить ситуацию, в данной работе проводятся детальные расчеты на различных сетках с использованием многоблочных вычислительных технологий [12, 13] нестационарных отрывных течений около профиля NACA 0012 при фиксированном числе Рейнольдса 40 000 и сравниваются результаты, полученные с использованием нескольких полуэмпирических моделей дифференциального типа, с вышеуказанными экспериментальными оценками. Важное значение для оценки применимости двумерного подхода к описанию турбулентных отрывных течений имеет сопоставление полученных численных прогнозов с данными физического эксперимента Игараши [14] для турбулентного обтекания поперечной толстой пластины.
ОБЪЕКТ ИССЛЕДОВАНИЯ
Профиль NACA 0012 является классическим симметричным профилем с аналитическим заданием координат точек верхней и нижней дужек контура [15]. Длина хорды L выбирается в качестве линейного масштаба обезразмеривания. Решение задачи строится в поточной системе координат (x, y), имеющей начало в середине хорды, а ось x ориентирована по направлению невозмущенного потока. Наклон хорды профиля по отношению к вектору скорости невозмущенного потока U задается углом атаки а. Задача обтекания несжимаемой вязкой жидкостью профиля NACA 0012 решается в прямоугольной расчетной области, границы которой располагаются на достаточном расстоянии от профиля, причем входная граница перпендикулярна скорости U. В начальный момент времени считается, что профиль, движущийся со скоростью U, внезапно тормозится, и в дальнейшем происходит процесс постадийного формирования его вихревого обтекания в зависимости от числа Рейнольдса, определенного по U, L и кинематической вязкости среды (см., например, [12]).
В отличие от использования для решения задачи моноблочных или неструктурированных сеток, как, например, в [6], в данной работе развивается альтернативный подход на основе многоблочных сеток, объединяющий вводимые априори разнообразные по типу и масштабу структурированные сетки простой топологии (близкие к ортогональным), накладываемые друг на друга с пересечением. Каждая из сеток нацелена на отображение соответствующего элемента гидродинамической структуры обтекания профиля. Так, на рис. 1 внешняя декартовая сетка А согласовывается с границами расчетной области, расположенными достаточно далеко (на несколько сотен хорд) от профиля (его даже не видно на рис. 1,а). Нестационарный дальний след за профилем детально воспроизводится на внутренней декартовой сетке В (рис. 1, б). Весьма подробная цилиндрическая сетка, согласованная с омываемой поверхностью профиля, нацелена на отображение сравнительно тонких пограничных и сдвиговых слоев. Наложение этой криволинейной сетки С на декартовые сетки А и В разной густоты показано на рис. 1, в.
Важным методическим этапом численного исследования является сравнительный анализ результатов, полученных с использованием дифференциальных полуэмпирических моделей турбулентности различного типа (одно- и двухпараметрических). Наряду со стандартной версией модели переноса сдвиговых напряжений MSST Ментера в варианте 2003 г. [16] применяется ее модификация MSST SM, учитывающая влияние кривизны линий тока на характеристики турбулентности в рамках подхода Смирнова — Ментера [17]. Для сравнения также используется одно-параметрическая модель переноса вихревой вязкости Спаларта — Аллмараса SA [18].
Помимо представленных сеточных методов численного моделирования, отдельные расчеты обтекания профиля NACA 0012 были выполнены с помощью бессеточного метода решения урав-
Рис. 1. Фрагменты многоблочной сетки: внешняя декартовая А (а), внутренняя декартовая В для отображения дальнего следа (б) и наложение декартовых сеток А, В и цилиндрической сетки С в ареале профиля под углом атаки 45°
нений Навье — Стокса — метода вязких вихревых доменов (ВВД) [19], в котором вообще отсутствуют ограничения на размеры расчетной области, а степень разрешения вихревых структур определяется лишь принятым уровнем начальной дискретизации поверхности профиля. Применение такого комплексного подхода к моделированию нестационарного обтекания профиля диктуется нацеленностью на отображение дальнего вихревого следа за ним.
Большое методологическое значение имеет сохранение в расчетной области разгонного вихря, особенно интенсивного при больших углах атаки. Именно это обстоятельство предопределило выбор значительных размеров расчетной области. Таким образом, время решения задачи подбирается так, чтобы практически весь вихревой след вплоть до разгонного вихря располагался внутри сетки А.
Фиксирование числа Рейнольдса в расчетах не связано с принципиальными методическими ограничениями, а также со стремлением экономить вычислительные ресурсы. Главным образом оно обусловлено необходимостью тщательной обработки огромного объема информации по
нестационарным процессам обтекания профиля. Хотя нельзя не отметить, что для числа Яе = 4 • 104 существенно снижаются требования к выбору величины пристеночного шага.
Как уже указывалось, обработка результатов расчета нестационарных характеристик течения представляет особую проблему, причем она связана не только со значительными массивами данных. В настоящей работе продолжается развитие подхода, ранее использованного в [8 —11], заключающегося в достижении квазипериодического режима изменения характеристик обтекания профиля (ранее это был цилиндр). Обычно расчеты нестационарных течений приводят к хаотическим колебательным процессам, обработка которых проводится статическими методами. Здесь наблюдается изменение характеристик на периоде колебаний интегрального параметра, например, подъемной силы. Осреднение процесса происходит по этому периоду. Изучение поведения решения на этом временном отрезке позволяет детально проанализировать нестационарный режим течения с акцентом на динамику вихревых структур.
ПОСТАНОВКА ЗАДАЧИ И ТЕХНОЛОГИЯ РЕШЕНИЯ
Задача турбулентного нестационарного обтекания несжимаемой вязкой жидкостью профиля под углом атаки решается в расчетной прямоугольной области с фиксированными входными условиями.
В качестве масштабов обезразмеривания выбираются скорость потока на входе и и длина хорды Ь. Продольная компонента скорости и =1, а поперечная — V = 0, степень турбулентности 1% (как в рабочей части аэродинамической трубы), масштаб турбулентности 1. На выходных проточных границах задаются условия продолжения решения (мягкие граничные условия). На омываемой поверхности профиля ставится условие прилипания. Расчет характеристик пристеночной турбулентности проводится аналогично методологии [13].
Решение нестационарных осредненных по Рейнольдсу уравнений Навье — Стокса (URANS) проводится в рамках метода коррекции давления, основанного на концепции расщепления по физическим процессам и связанного с заменой уравнения неразрывности на уравнение для поправки давления (SIMPLEC-метод).
Нестационарный подход Рейнольдса к описанию турбулентности позволяет, как и в стационарном случае, перейти к осредненному течению и описывать корреляции составляющих пульсаций скорости полями вводимых турбулентных характеристик. При этом осредненное течение может изменяться во времени.
Как уже отмечалось, ищется периодическое решение задачи с циклическим изменением интегральных характеристик, в частности, действующей на профиль подъемной силы. Подобно [12, 13], обобщенное уравнение переноса записывается в приращениях зависимых переменных, причем явная часть дискретизируется по схемам высокого (не ниже второго) порядка аппроксимации (противопоточная схема с квадратичной интерполяцией Леонарда [20] для конвективных членов уравнений количества движения, схема TVD [21] для конвективных членов уравнений характеристик турбулентности и центрально-разностная схема для диффузионных членов), а неявная часть представляется с помощью противопоточной схемы с односторонними разностями. Решение на каждом временном шаге находится глобальными итерациями [13]. Использование центрированных сеток с размещением зависимых переменных в центрах расчетных ячеек приводит к необходимости монотонизации поля давления за счет введения коррекции Рхи — Чоу со специально подобранным из тестов [12, 13] сомножителем, равным 0.1. Расчеты проводятся на многоблочных сетках (см., например, рис. 1). Внешняя сетка построена в области размером 300 х 200 с центром координатной системы, удаленным от входной границы на 100 хорд профиля. Ячейки сгущаются по мере приближения в зоне расположения профиля. Минимальные шаги внешней сетки равны 0.1. Количество узлов сетки 672 х 201. Внутренняя сетка близка к равномерной с шагом 0.03 по обеим координатам, содержит 633 х 164 ячейки и охватывает пространство протяженностью порядка 60 хорд профиля. Цилиндрическая сетка согласовывается с профилем. Шаг сетки в окрестности носка профиля — 0.001, шаг у задней кромки — 0.001 (радиус скругления задней кромки — 0.002). Шаг около срединной части профиля — 0.01. Пристеночный
шаг — 8 х 10-5. Сетка содержит 329 х 62 ячейки. Временной шаг выбирается равным 0.01.
МЕТОДИЧЕСКИЕ ИССЛЕДОВАНИЯ. ВЕРИФИКАЦИЯ.
АНАЛИЗ ИНТЕГРАЛЬНЫХ СИЛОВЫХ ХАРАКТЕРИСТИК
Выполненная серия расчетов обтекания профиля при варьировании углом атаки от 15 до 150° характеризуется выходом на автоколебательный режим с периодическим поведением погрешности, интегральных продольных и поперечных нагрузок, как демонстрируется на рис. 2 для а = 60 и 90°. Аналогично турбулентному обтеканию кругового цилиндра [12, 13], начальная фаза взаимодействия потока с профилем протяженностью 10—20 безразмерных временных единиц характеризуется повышенными пульсирующими нагрузками, постепенно снижающимися по мере развития следа.
Анализ поведения поперечных нагрузок на автоколебательном режиме обтекания профиля позволил выделить участки их периодического изменения (рис. 3, 4). В отдельных случаях (при углах атаки 135 и 150°) приходится рассматривать сдвоенные участки. Для сопоставления зависимостей интегральных нагрузок X (t), Y (t) на участках их периодического изменения вводится приведенное время t*, начало отсчета которого берется в момент достижения минимума Y.
Сравнение зависимостей X (t*) и Y (t*) на рис. 3 носит также методический характер, так
как результаты получены в расчетных областях с существенно различающимися размерами 40 х 20 и 300 х 200. Как уже указывалось, для больших размеров расчетных областей практически весь вихревой след вплоть до разгонного вихря находится внутри области, а для области сокращенных размеров разгонный вихрь покидает область. Отмечается, что периодический характер процесса вихреобразования затруднительно выявить в областях сокращенных размеров, что повлияло на проведение расчетов обтекания профиля при больших углах атаки (свыше 60°). В то же время, в пределах углов атаки до 60° (см. рис. 3) согласие результатов расчетов нестационарных интегральных характеристик профиля NACA 0012 вполне удовлетворительное (см. кривые 7,11 и 8,12).
12 г
Y .
Рис. 2. Поведение во времени лобовойX (a, в) и подъемной Y (б, г) силы, действующей на профиль NACA 0012
под углами атаки:
a, б — 60°; в, г — 90°
Рис. 3. Сравнение зависимостей от приведенного времени Г продольной (а) и поперечной (б) нагрузки на профиль при углах атаки от 15 до 60° для расчетных областей с размерами 50 x20 (кривые 1 — 8, полученные с использованием модификации МББТ Смирнова — Ментера) и 300 х 200 (кривые 9—12):
1 — а =15°; 2 — 20°; 5 — 25°; 4 — 30°; 5 — 35°; 6 — 40°; 7 — 45°; 8 — 55°; 9 — 15°; 10 — 30°; 11 — 45°; 12 — 60°
/
на профиль при углах атаки от 75 до 157.5°:
1 — а=75°; 2 — 85°; 5 — 90°; 4 — 105°; 5 — 120°; 6 — 135°; 7 — 150°; 8 — 157.5°
Следует отметить, что по мере увеличения угла атаки периодические распределения интегральных нагрузок во времени становятся все более отличными от синусоидальных, характерных для малых углов (до 25° и свыше 155°).
В [6] на основе рассогласования результатов численных прогнозов, полученных с использованием моделей Ментера MSST (вариант 1993 г.) [22] и Спаларта — Аллмараса SA (Re = 105), а
также данных Хернера [4], сделан вывод о неприемлемости использования двумерного подхода URANS для расчета обтекания профиля под большими углами атаки. К указанным выше сведениям можно добавить оценки [5] (рис. 5). Однако полученные при Re = 4 -104 в настоящей работе результаты хорошо согласуются с данными [6]. Более того, они также коррелируют с результатами расчетов обтекания профиля бессеточным методом вязких вихревых доменов [19]. Полученные результаты указывают на ограниченность вывода [6].
Рассмотрим вариант расположения профиля под углом атаки 90°, этот случай мало отличается от поперечно обтекаемой пластинки. Экспериментальные данные [14] показывают, что
cx = 2.8 для пластинки с относительной толщиной 0.2 при Re = 4 -104. По данным [20] для тонкой пластинки cx ~ 2.5. Для сравнения с экспериментальными данными была решена специальная тестовая задача обтекания поперечно ориентированной толстой пластины при равном экспе-
Рис.5. Зависимость от угла атаки осредненных по периоду колебаний коэффициентов лобового сопротивления (а) и подъемной силы (б):
1 — расчет иКА^-М88Т [3]; 2 — расчет ИЯАШ^А [3]; 5 — расчет иКА^-М88Т (данная работа); 4, 5 — данные [10], [11] соответственно; 6 — метод ВВД;-----[4]
¿jjj^PTl^_I_._._L_i_,_L
uULbS-J 20° ¿r Í:
Рис. 6. Сравнение зависимостей от умеренных (а, б) и больших (в, г) углов атаки осредненных по периоду колебаний cx (а) (а, в) и cy (а) (б, г), полученных с помощью URANS для расчетных областей с размерами 300 х 200 (1) и 50 х 20 (2), а также в рамках метода ВВД (3), с оценками (4) [11]
риментальному числе Re. Расчетный прогноз составил cx = 2.76. Таким образом, рассогласование данных по сопротивлению профиля при а = 90° в [6] нельзя отнести к неадекватности двумерного подхода URANS. Скорее всего, причина отмеченных выше случаев рассогласования заключается в различии двумерного и пространственного характера обтекания профиля, наблюдаемого и в физических экспериментах [4, 5, 14].
Более детальное сравнение расчетных и оценочных интегральных силовых нагрузок на профиль в диапазонах умеренных углов атаки представлено на рис. 6. При малых а (до 7—8°) численные прогнозы по сеточным и бессеточным методам практически совпали с оценочными данными [4]. Этот диапазон углов атаки соответствует безотрывному обтеканию профиля NACA 0012.
При определении характеристик отрывных режимов обтекания профиля наблюдаются некоторые рассогласования численных прогнозов и их прогрессирующее отклонение от оценочных данных. Хотя зависимости для cx (а) хорошо сочетаются для сеточных и бессеточных методов (рис. 6, а), причем примерно до 40° они подчиняются близкому к квадратичному закону, результаты расчетов Су (а) при умеренных а (порядка 15—25°) демонстрируют различие в отображении «ложки» локального падения Су при формировании отрывного характера течения. При этом
прогнозы по URANS для различных расчетных областей довольно близки друг к другу.
На больших углах атаки, близких к 180°, численные прогнозы cx практически не отличаются от оценочных данных [4], а прогнозы cy, хотя и несколько отличаются от данных [4], но воспроизводят мелкую «ложку» cy. Расслоение результатов при а< 157.5°, скорее всего, свиде-
* Расчет бессеточным методом ВВД [19] выполнил Я. А. Дынников.
Рис. 7. Зависимости от угла атаки пульсационных продольных (а) и поперечных (б) нагрузок на профиль
Рис. 8. Зависимости от угла атаки числа Струхаля (а) и приведенного числа Струхаля (б):
1 — 300 X 200; 2 — 50 х 20
тельствует о пределах применимости двумерной модели для описания существенно пространственного обтекания профиля.
Уровень пульсирующих нагрузок на профиль довольно велик (рис. 7). Если пульсации сх достигают 1 при максимуме сх порядка 3 (т. е. составляют примерно треть), то пульсации су доходят до 1.8 при наибольшей по модулю величине су, равной 2 (т. е. оказываются одного порядка). При этом максимальные пульсации су соответствуют развороту профиля острым концом к набегающему потоку. Снижение пульсаций су при положении профиля, близком к нормальному, связано с сокращением проекции хорды профиля на направление набегающего потока.
Обработка данных по периодам колебаний на рис. 3, 4 позволяет построить зависимости определенного по хорде профиля числа Струхаля (8Ь) от угла атаки (рис. 8, а), причем прогнозы для расчетных областей различного размера удовлетворительно согласуются между собой. Значительный разброс от 0.7 до 0.1, по-видимому, является нефизическим и связан с тем, что проекция профиля на нормаль к потоку сокращается при приближении углов атаки к нулю и к 180° соответственно. Определение приведенного по переменному размеру — проекции хорды профиля на направление, перпендикулярное набегающему потоку, существенно сокращает диапазон варьирования числа — от 0.11 до 0.19 (рис. 8, б). Таким образом, угол атаки существенно влияет на характер периодических колебаний в следе, так как число при а = 90° оказывается порядка 0.1, что более чем вдвое ниже, чем для кругового цилиндра. Во многом это определяется изменением топологии вихревых структур в дальнем следе, которые анализируются ниже.
АНАЛИЗ ЭВОЛЮЦИИ ВИХРЕВЫХ СТРУКТУР И ЛОКАЛЬНЫХ ХАРАКТЕРИСТИК
В БЛИЖНЕМ И ДАЛЬНЕМ СЛЕДЕ
На рис. 9 сравниваются картины вихревых структур в следе за профилем в различные моменты времени 40(А), 60(5) и 100(С) при варьировании в широких пределах углом атаки от 15 до 165°. Визуализация структур осуществляется изолиниями вихревой вязкости, причем для каждого из рассчитанных вариантов число изолиний фиксировано и равно 21, а максимум вихревой вязкости различный.
Проводится качественное сравнение формирующихся вихревых дорожек по форме сгенерированных вихрей, их концентрации и протяженности. По мере выхода интегральных характеристик профиля на автоколебательный режим сокращается число вариантов анализируемых дорожек, т. е. сокращается диапазон варьирования углом атаки. На малых и очень больших, близких к 180°, углах атаки выход на автоколебательный режим происходит на временах порядка 40 (рис. 9, А), причем протяженность следа, в целом, превосходит протяженность вихревых дорожек при углах атаки в диапазоне 60—120°.
Вихревая дорожка за профилем при малых углах атаки (порядка 15°) имеет характер сцепленных попарно вихрей разной ориентации, образующихся при периодическом срыве потока поочередно со срединной части контура профиля и с его острой кромки [11]. Как видно из рис. 9, А, в момент г = 40 при а = 15 и 165° структура следа получается упорядоченной и однородной, причем волнистая в ближнем следе граница следа постепенно выглаживается и в дальнем следе образуется вихревая полоса постоянной ширины, на конце которой сохраняется разгонный вихрь, а точнее то, что от него осталось к этому моменту времени. Выход на автоколебательный режим обтекания профиля происходит гораздо раньше отмеченного момента времени. Следует подчеркнуть, что эволюция вихревого следа, в целом, не является периодическим процессом.
С увеличением угла атаки до 30° наблюдается более интенсивное смешение в следе, чем при малых углах атаки (см. рис. 9). Более рельефно происходит спаривание вихрей разной ориентации с формированием полосчатой дорожки постоянной ширины, которая становится со време-
Рис. 9. Картины вихревой вязкости для а от 15 до 150° с шагом 15° (а — /), 157.5° (к) и 165°(/)
при г = 40(А), 60(5), 100(С)
нем (рис. 9, В, С) все более протяженной. Обращает на себя внимание, что конфигурация разгонного вихря меняется незначительно за 20 единиц безразмерного времени от ^ = 40 до 60.
Расположение профиля под углом атаки 157.5° сопровождается разворотом задней острой кромки навстречу набегающему потоку. Как результат, сходящие вихри в начальные моменты времени (примерно до 30 единиц) собираются в отдельные вихревые сгустки (см. рис. 9). В дальнейшем (см. рис. 9, С) наблюдается более упорядоченное поведение крупномасштабных вихрей с формированием дорожки, аналогичной случаю угла атаки 30°.
Следующая в систематизации группа вихревых дорожек соответствует углам атаки 45, 135 и 150°. Во всем временном диапазоне от 40 до 100 ее отличает формирование изолированных вихрей при сходе их с поверхности профиля. Правда, позже эти вихри все же объединяются в дорожку, которая имеет полосчатый «протекторный» характер.
Последний тип вихревых дорожек соответствует углам атаки от 60 до 120°. Для них характерна довольно резкая генерация вихрей в ближнем следе, сопровождающаяся расширением вихревой дорожки. Топология дорожек однородна и имеет указанный выше полосчатый «протекторный» вид.
Осредненные по периоду колебаний поперечной силы (см. рис. 3, 4) картины линий тока около профиля под различными углами атаки в двух масштабах показаны на рис. 10 и 11. В определенной степени они подтверждают выводы, сделанные при анализе нестационарных вихревых структур. Действительно, судя по рис. 10, отрывные течения около профиля можно разбить на группы: первая — 30, 45, 135 и 150°, вторая — от 60 до 120°. Для первой группы характерна небольшая протяженность циркуляционной зоны в следе, причем образовавшиеся крупномасштабные вихри примыкают к профилю. Во второй группе за профилем формируются множественные протяженные зоны отрыва, причем они простираются довольно далеко вниз по потоку.
Более детально характер присоединенных к профилю крупномасштабных вихрей представлен на рис. 11. По мере увеличения угла атаки вихрь с точкой отрыва в окрестности передней сглаженной кромки, первоначально занимающий всю верхнюю дужку профиля, начинает сокращаться, т. е. точка присоединения перемещается в сторону от задней острой кромки. При этом продольный (вдоль профиля) размер вихря оказывается не меньше половины хорды, за исключением угла атаки 90°, при котором фактически этого вихря нет, т. е. он является частью общей вихревой системы. В то
Рис. 10. Осредненные картины обтекания профиля от 30 до 150° через 15° (а — г)
Рис. 11. Осредненные по периоду колебаний су картины линий тока около профиля для углов атаки от 15 до 165°
с шагом 15° (а—к)
же время поперечный размер вихря остается примерно одинаковым, несмотря на то, что вихрь при больших углах атаки оказывается погруженным в вихревой след.
При а = 90° структура вихревого следа в окрестности профиля симметрична относительно центральной линии, проходящей через начало координат в середине профиля. При развороте профиля острой кромкой навстречу набегающему потоку присоединенный вихрь с точкой отрыва на кромке сначала несколько увеличивается в размерах (продольный размер оказывается длиннее полухорды профиля, а поперечный превосходит хорду), затем, по мере увеличения угла атаки, уменьшается. При а = 165° структура присоединенного вихря примерно такая же, как и при а = 15°. Наблюдается образование вторичных малогабаритных вихрей в окрестности кромок при больших углах атаки, а также отрывных пузырей при а = 135 и 150° (рис. 11, г, у).
Следующая группа расчетных материалов носит, в основном, методический характер и связана с анализом поведения локальных характеристик обтекания профиля при угле атаки 90°, причем в расчетах используются различные модели турбулентности: стандартная МББТ и с поправкой Смирнова — Ментера (МББТ БМ), а также модель Спаларта — Аллмараса БЛ.
Донное давление за профилем при угле атаки 90° на периоде колебаний поперечной силы изменяется в достаточно большом диапазоне от -3 до нуля (в обезразмеренном виде эта величина имеет смысл половины коэффициента донного давления). При осреднении процесса получается, что осредненное давление довольно сильно уменьшается в центральной части профиля, достигая уровня -1.5, в то время как вблизи кромок давление увеличивается до -1.1 в окрестности перед-
Рис. 12. Сравнение вихревых структур, рассчитанных МББТ (а), МББТ БМ (б) и БЛ (в), а также распределений продольных (г) и поперечных (д) нагрузок на профиль при
а =90°:
1 — МББТ; 2 — МББТ БМ; 3 — БЛ (показаны два периода колебаний, Г = 100)
ней кромки и до -0.9 вблизи задней острой кромки. Пульсации донного давления высоки, составляя минимум 0.9 в центральной части и имея два локальных максимума, равных 1.5 и 1.4 в верхней и нижней половине верхней дужки профиля. На нижней дужке, обращенной навстречу потоку, наибольшие пульсации давления реализуются в окрестности кромок и достигают 1 —1.2.
Судя по рис. 12, вихревой след распадается на характерные группы. Разгонный вихрь вместе с системой беспорядочных мелких комкообразных вихрей соединяется перемычкой со следующей группой, напоминающей по неупорядоченной структуре первую. Затем через переходный слой следует упорядоченная группа вихрей протекторной формы, которые, по-видимому, и определяют автоколебательный режим интегральных характеристик профиля.
На рис. 12, a — в, помимо прогнозов по модели Ментера, представляются картины вихревой вязкости, отображающие вихревые структуры в следе за профилем под углом атаки 90°, полученные с помощью модели Спаларта — Аллмараса SA. Обращает на себя внимание гораздо меньший (в 4 раза) уровень максимальной вязкости, генерируемой моделью SA. Однако, как показано на рис. 12, г, д, распределения интегральных силовых нагрузок на профиль отличаются незначительно. Упорядоченный характер вихревых структур на рис. 12, в сохраняется.
ВЫВОДЫ
1. Численно обоснован автоколебательный режим изменения во времени интегральных нагрузок на профиль в широком диапазоне углов атаки (от 15 до 165°) при моделировании в рамках решения осредненных по Рейнольдсу уравнений Навье — Стокса с помощью многоблочных вычислительных технологий и полуэмпирических моделей турбулентности на протяженных сетках с разрешением в расчетной области разгонного вихря. Определены периоды автоколебаний по протяженности периодического изменения поперечной силы, действующей на профиль.
2. Обоснована приемлемость двумерного подхода к моделированию вихревого следа за профилем на примере обтекания поперечной толстой (0.2 хорды) пластины. Рассогласование данных по сопротивлению профиля при а = 90° в [6] не является свидетельством неадекватности двумерного подхода URANS. Причина, скорее всего, заключается в различии двумерного и пространственного характера обтекания профиля, наблюдаемого, в том числе, в физических экспериментах [4, 5, 14].
3. В ходе численного исследования показано, что прогнозы осредненных интегральных и локальных нагрузок на профиль, полученные с использованием в подходе URANS двухпарамет-рической стандартной модели переноса сдвиговых напряжений Ментера (2003), модифицированной версии этой модели в рамках подхода Смирнова — Ментера, однопараметрической модели вихревой вязкости Спаларта — Аллмараса, а также при решении уравнений Навье — Стокса бессеточным методом вихревых доменов, хорошо согласуются между собой. При этом уровень максимальной вихревой вязкости в следе для модели Спаларта — Аллмараса в четыре раза ниже, чем для модели Ментера, в том числе в модифицированном варианте. Следует подчеркнуть, что причины хорошего согласия прогнозов cx для профиля при решении уравнений URANS сеточным методом и неосредненных уравнений Навье — Стокса бессеточным методом требуют дополнительного анализа.
4. Анализ эволюции вихревых структур в следе за профилем в моменты времени 40, 60 и 100, а также осредненных по периоду колебаний поперечной силы картин отрывного обтекания профиля позволил выделить две структурные группы с образованием компактных (на углах атаки 30, 45, 135, 150°) и протяженных (на углах атаки от 60 до 120°) циркуляционных зон. Отмеченные группы разделяются также по поведению турбулентных характеристик: для первой характерно смещение максимумов вихревой вязкости на значительное расстояние от профиля в область дальнего следа, а для второй — локализация максимума энергии турбулентности и вихревой вязкости в ближнем следе. Также следует отметить, что второй группе соответствует ярко выраженный участок изменения пульсаций cx .
5. Как и ожидалось, при безотрывном обтекании профиля при малых углах атаки (до 7°) прогнозы cy хорошо согласуются с оценками [4, 5]. На отрывных режимах обтекания с образованием отрывных пузырей и низкоинтенсивных циркуляционных зон при углах атаки до 20° и со 157.5 до 180° получено вполне удовлетворительное согласование осредненных по периоду
автоколебаний расчетных результатов и оценочных данных [1] по cx и несколько худшее по cy ,
однако при этом качественно и количественно воспроизводится «ложка» на кривой cy. Не вполне
понятны причины рассогласования данных, связанные, возможно, с неясностью оценивания экспериментальных нестационарных нагрузок.
6. При анализе осредненных нагрузок на профиль показано, что в диапазоне углов атаки от 40 до 135o пульсации продольной силы составляют примерно 30% от максимальной нагрузки на профиль, а максимальные пульсации поперечной силы (при 135o) одного порядка с наибольшей поперечной силой, действующей на профиль.
7. На основе анализа автоколебаний поперечной силы, действующей на профиль, показано, что определенное по проекции хорды на нормаль к вектору скорости набегающего потока приведенное число Струхаля в диапазоне углов атаки от 40 до 150o меняется от 0.11 до 0.14.
Работы выполнены при финансовой поддержке РФФИ (проекты M 11-01-00039, 12-0S-90001 и 12-01-009S5) и ФЦП «Научные и педагогические кадры инновационной России» M 14.740.11.0SS0.
ЛИТЕРАТУРА
1. Головина Н. В. Сравнение численных расчетов методом, основанным па разпо-стной схеме Годунова — Колгана — Родионова, с экспериментальными данными для случая трансзвукового обтекания профиля RAE 2S22 // Ученые записки ЦАГИ. 2009. Т. XL, M 5, с. 41—52.
2. Пилипепко А. А., Полевой О. Б., Приходько А. А. Численное моделирование влияния числа Маха и угла атаки на режимы трансзвукового турбулентного обтекания аэродинамических профилей // Ученые записки ЦАГИ. 2012. Т. XLIII, M 1, с. 3—31.
3. Богатырев В. В., Петров А. В. Численное моделирование воздействия поперечной струи па обтекание крылового профиля // Ученые записки ЦАГИ. 2012. Т. XLIII, M 2, с. 3 —9.
4. H o e r n e r S. F. Fluid-dynamic lift. Hoerner Fluid Dynamics. — Brick Town. NJ, 1975, p. 4—19.
5. Sheldahl R. E., Klimas P. C. Aerodynamic characteristics of seven symmetrical airfoil sections through 1S0-degree angle of vertical axis wind turbines // Sandia National Lab. Energy report. SANDS0-2114. Ш1, 120 p.
6. Гарбарук А. В., Стрелец М. Х., Шур М. Л. Моделирование турбулентности в расчетах сложных течений. Уч. пособие. — СПб: Изд. Политехн. ун-та, 2012, SS с.
7. Isaev S. A., Zhdanov V.L.Niemann H.-J. Numerical study of the bleeding effect on the aerodynamic characteristics of a circular cylinder // J. Wind Eng. Ind. Aerodyn. 2002. V. 90, Issue 11, p. 1217—1226.
S. Исаев С. А., Баранов П. А., Кудрявцев Н. А., Жукова Ю. В. Численное моделирование нестационарного теплообмена при турбулентном обтекании кругового цилиндра. Ч. 1. Методическое исследование // Теплофизика и аэромеханика. 2005. Т. 12, M 1, с. 27—39.
9. Исаев С. А., Баранов П. А., Кудрявцев Н. А., Жукова Ю. В. Численное моделирование нестационарного теплообмена при турбулентном обтекании кругового цилиндра. Ч. 2. Анализ автоколебательного режима // Теплофизика и аэромеханика. 2005. Т. 12, M 2, с. 271—2S3.
10. Исаев С. А., Баранов П. А., Кудрявцев Н. А., Лысенко Д. А., Уса-чов А. Е. Сравнительный анализ пакетов VP2/3 и FLUENT при расчете нестационарного обтекания кругового цилиндра с использованием моделей турбулентности Спаларта — Аллмараса и Мептера // ИФЖ. 2005. Т. 7S, M 6, с. 14S—162.
11. Исаев С. А., Судаков А. Г., Баранов П. А., Усачов А. Е., Стри-жак С. В., Лохапский Я. К., Гуверпюк С. В. Разработка, верификация и применение основанного на многоблочных вычислительных технологиях распараллеленного пакета открытого типа VP2/3 для решения фундаментальных, прикладных и эксплуатационных задач аэромеханики и теплофизики // Вестник ЮУрГУ. Серия «Математическое моделирование и программирование». 2009. M 17 (150), вып. 3, с. 59—72.
12. Управление обтеканием тел с вихревыми ячейками в приложении к летательным аппаратам интегральной компоновки (численное и физическое моделирование) / Под ред. А. В. Ермишина и С. А. Исаева. — М.: МГУ, 2003, 360 с.
13. Быстров Ю. А., Исаев С. А., Кудрявцев Н. А., Леонтьев А. И. Численное моделирование вихревой интенсификации теплообмена в пакетах труб. — СПб: Судостроение, 2005, 39S c.
14. Igarashi T., Terachi N. Drag reduction of flat plate normal to airstream by flow control using a rod // Expanded abstracts of 4th Int. Colloquium on Bluff Body Aerodynamics & Application. — Bochum: Ruhr-University Bochum, 2000, p. 635—638.
15. Флетчер К. Вычислительные методы в динамике жидкостей. Т. 2. — М.: Мир, 1991, 552 с.
16. Menter F., Kuntz M., Langtry R. Ten years of industrial experience with the SST turbulence model // Turbulence, Heat and Mass Transfer 4. Ed. K. Hajalic, Y. Nogano, M. Tummers. Begell House, Inc, 2003, 8 p.
17. Smirnov P., Menter F. Sensitization of the SST turbulence model to rotation and curvature by applying the Spalart — Shur correction term // Proc. ASME Turbo Expo Conf., 2008. N GT2008-50480, 10 p.
18. Spalart P. R., A 11mare s S. R. A one-equation turbulence model for aerodynamic flows // AIAA Paper 92-0439, 1992, 22 p.
19. Андронов П. Р., Гувернюк С. В., Дынникова Г. Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. — Изд. МГУ, 2006, 184 с.
20. L e o n a r d B. P. A stable and accurate convective modeling procedure based on quadratic upstream interpolation // Comp. Meth. Appl. Mech. Eng. 1979. V. 19, N 1, p. 59—98.
21. L i e n F. S., L e s c h z i n e r M. A. Approximation of turbulence convection in complex flows with a TVD-MUSCL scheme // Proc. 5th Int. Symp. Refined flow modelling and turbulence measurements. — Paris: 1993, p. 183 —190.
22. Menter F. R. Zonal two equation k-m turbulence models for aerodynamic flows // AIAA Paper 93-2906, 1993, 21 p.
23. Белоцерковский С. М., Ништ М. И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. — М.: Наука, 1978, 352 с.
Рукопись поступила 15/II2013 г.
