CC BY
78
__________УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XXIX 1998 ~
№3—4
УДК 629.735.33.015.3.025.33
МЕТОД РАСЧЕТА НЕСТАЦИОНАРНОГО ОБТЕКАНИЯ ПРОФИЛЯ С ИНТЕРЦЕПТОРОМ И ЕГО АЭРОДИНАМИЧЕСКИЕ ХАРАКТЕРИСТИКИ
В. В. Богатырев
Предложен метод расчета нестационарной картины течения и аэродинамических характеристик для профиля с интерцептором. Решается нестационарная задача в рамках идеальной среды с поправкой на толщину вытеснения пограничного слоя. Для моделирования вихревого следа используется метод дискретных вихрей.
Сравнение с экспериментальными данными показало возможность оценки распределенных и суммарных аэродинамических характеристик, в том числе шарнирных моментов интерцептора, а также характерных частот вихреобразования.
Системы поперечного управления современных самолетов отличаются как механической, так и аэродинамической сложностью. Интерцепторы в настоящее время обеспечивают высокую поперечную управляемость, выполняя, кроме того, и функции уменьшения подъемной силы и торможения. Многофункциональность интерцепторов увеличивает сложность системы управления. Поэтому при проектировании интерцепторов необходимо правильно понимать их аэродинамику, учитывать производимый ими вихревой след и его влияние на горизонтальное оперение.
Аэродинамические характеристики крыльев с интерцепторами исследовались экспериментально ранее во многих работах [1]—[6]. Был разработан ряд стационарных методов расчета для профиля с интерцептором [7]—[И]. Однако нестационарная картина течения и вихревые структуры, образующиеся за крылом с отклоненным интерцептором, исследованы меньше как экспериментально, так и теоретически [2]—[6].
Обтекание профилей с отклоненными на большие углы интерцепторами характеризуется существованием обширных отрывных зон и образованием развитого вихревого следа, т. е. такое обтекание по своей природе является нестационарным. Обычно вихревой след является периодическим и поэтому важно знать характерные частоты вихреобра-
зования, так как вихревой след позади интерцептора может вызывать бафтинг горизонтального оперения.
Нестационарное обтекание профиля с интерцептором в настоящей работе рассчитывается с помощью панельного метода распределенных особенностей с учетом толщины вытеснения пограничного слоя. Привлекается ряд эмпирических закономерностей. Вихревой след при этом моделируется дискретными вихрями. Метод дискретных вихрей нашел широкое применение при расчете отрывных течений в методиках школы С. М. Белоцерковского [12].
Описанный в данной работе метод и созданная на его основе программа расчета позволяют оценивать распределение давления по профилю с интерцептором. По распределению давления рассчитываются суммарные аэродинамические характеристики, в том числе шарнирные моменты интерцептора. Кроме того, расчет картины обтекания позволяет оценивать характерные частоты вихреобразования. Результаты расчета сравнивались с экспериментальными данными, содержащимися в [1]— [3], [Ю].
1. Метод расчета. Предлагаемый метод расчета отрывного обтекания профиля с интерцептором основан на решении нестационарной задачи в рамках идеальной среды с поправкой контура профиля на величину толщины вытеснения пограничного слоя. Вихревой след моделируется дискретными вихрями.
Пограничный слой при проведении расчетов считался квазиста-ционарным, т. е. в каждый момент времени рассматривался как стационарный, но в качестве определяющих параметров брались их мгновенные значения. Для расчета пограничного слоя (ламинарного и турбулентного) использовался метод описания пограничного слоя А. А. Дородницына в сочетании с методом конечных элементов [13].
Для определения турбулентной вязкости использовались теория пути смешения Прандтля с поправкой на демпфирующий фактор Ван-Дриста во внутренней, пристеночной подобласти турбулентного пограничного слоя и гипотеза Клаузера с модификацией Ю. В. Лапина и М. X. Стрельца во внешней подобласти, кроме того, учитывался коэффициент перемежаемости Клебанова [14], [15]. .
Точки перехода ламинарного пограничного слоя в турбулентный в приводимых далее расчетах фиксировались в соответствии с положением турбулизаторов в эксперименте. В качестве критерия отрыва турбулентного пограничного слоя используется условие превышения поперечной составляющей скорости в пограничном слое 10% от величины продольной составляющей скорости в пограничном слое на расстоянии от стенки, равном толщине вытеснения пограничного слоя.
При обтекании профиля с отклоненным на достаточно большой угол интерцептором перед осью вращения интерцептора обычно происходит отрыв пограничного слоя вследствие большого положительного градиента давления. Если угол атаки профиля не слишком велик, то оторвавшийся пограничный слой может присоединиться к поверхности интерцептора. В результате в области оси вращения интерцептора образуется зона возвратного течения. Приближенно это явление моделируется с помощью участка квадратичной параболы (рис. I). Парабола стро-
Рис. 1. Модификация контура профиля с учетом толщины вытеснения пограничного слоя и моделирования отрывного пузыря перед интерцептором
ится в системе координат, начало которой расположено в месте отрыва пограничного слоя перед интерцептором, а оси направлены соответственно по касательной и перпендикулярно к поверхности профиля в точке отрыва пограничного слоя (точка 5). Точка присоединения пограничного слоя (точка К) определяется параметром /д, характеризующим расстояние от точки присоединения до кончика интерцептора. В конкретных расчетах величина 1ц задавалась равной 2% хорды профиля. Ранее моделирование зоны возвратного течения перед интерцептором с помощью участка квадратичной параболы использовалось Зумвольтом [10].
Обтекание профиля с интерцептором рассчитывалось с помощью панельного метода с кусочно-линейным распределением вихрей у и кусочно-постоянным распределением источников д (рис. 2). Вихри распределялись по контуру, получаемому с учетом добавления толщины вытеснения пограничного слоя на большей части поверхности профиля и моделирования отрывного пузыря перед интерцептором. На оставшейся части верхней поверхности профиля позади оси вращения интерцептора распределялись источники. В результате при представлении рассматриваемых контуров с помощью панелей образуются две ломаные — первая соответствует большей части поверхности профиля и верхней поверхности интерцептора, вторая — части верхней поверхности профиля, находящейся в отрывной зоне позади интерцептора.
Рис. 2. Схема моделирования обтекания профиля с интерцептором
К крайним панелям первой ломаной добавлялись по одной добавочной панели (на рис. 2 на добавочных панелях нет крестика посередине), которые моделировали самый начальный участок вихревых пелен, сходящих с кончика интерцептора и нижней поверхности задней кромки профиля. Плотность интенсивности вихря вдоль добавочных панелей была постоянной. Вихревой след моделируется сходящими на каждом шаге расчета по времени с кончика интерцептора и нижней поверхности задней кромки профиля дискретными вихрями. Каждая из добавочных панелей используется для определения циркуляции первого сходящего из соответствующей точки дискретного вихря, который в
дальнейшем уносится потоком, не меняя своей циркуляции. В результа-
1 л
те образуются две системы Г(‘ и Г. вихрей, сходящих соответственно с
нижней поверхности задней кромки профиля и с кончика интерцептора (рис. 2). Каждая из добавочных панелей направлена так же, как и соответствующая примыкающая спереди панель. Длина добавочной панели определяется как произведение скорости схода вихрей на шаг расчета по времени. В предположении малости скорости на внутренней стороне панелей за скорость схода вихрей принимается половина плотности интенсивности вихря в середине примыкающей к добавочной спереди панели. Эта же скорость применяется для определения положения первого сходящего дискретного вихря, который сходит в поток из середины добавочной панели вдоль ее направления.
Таким образом, на первой ломаной с добавочными панелями непрерывно распределены вихри, плотность интенсивности которых у линейно меняется вдоль каждой из панелей и не терпит разрывов при переходе к соседней панели. При этом вдоль добавочных панелей плотность интенсивности вихря постоянна. На второй ломаной распределены источники, плотность интенсивности которых 9 постоянна вдоль каждой из панелей, но меняется при переходе к соседней панели.
В контрольных точках (серединах) каждой из панелей с распределенными вихрями, кроме добавочных, ставится условие непротекания. Суммарная интенсивность вихрей определяется условием выполнения на каждом шаге по времени теоремы Томсона о постоянстве циркуляции по замкнутому контуру, охватывающему тело и его вихревой след.
Для панелей с распределенными источниками налагается дополнительное условие о равенстве нулю суммарной обильности источников. Это условие требует введения дополнительной ‘неизвестной, необходимой для замыкания линейной системы уравнений. Эта неизвестная появляется после замены для панелей с распределенными источниками условия непротекания на условие одинаковости нормальной составляющей скорости в контрольных точках панелей. Величина этой нормальной составляющей скорости и служит дополнительной неизвестной. Величина дополнительной неизвестной получается малой (сотые-десятые доли процента от скорости набегающего потока) вследствие замкнутости контура профиля. Отсутствие суммарного источника исключает появление заметного «вытекания» изнутри или «втекания» внутрь контура профиля в областях в начале и в конце ломаной, составленной из панелей с распределенными источниками.
В результате нормальная м>„ и тангенциальная м>(. компоненты
скорости и потенциал <р/ в /-й контрольной точке будут определяться следующим образом:
к2 *3-1 2 Мт
*4 = ЛАууУ У + * ЕХ^г” + Ащ
./'=*1 +1 /'=*2+1 «=1у=1
*2 *з-1 2 N т
ч, = 2Х-ъ- + *12 + Вщ
У=*1+1 У=*2+1 т=\]=\
к2 Аг3-1 2 Мт
*< = + 1ч^ + 1Ес?#г"+с»(-
у'=*]+1 у=Лг2 +-1 /и= 1У=1
где Ау.., ВУу, Су — соответствующие коэффициенты влияния от распределенных вихрей (уу — плотность интенсивности вихря в у'-й узловой точке, у = £] +1, £2 (см. рис. 2)); Ац.., Вд.., Сд.. — коэффициенты влияния от распределенных источников Ц- — плотность интенсивности источника на у'-й панели, у=&2+1, &з -1 (см. рис. 2)); АТ.., Вг , Сг.. — коэффициенты влияния от дискретных вихрей с циркуляциями
Гу”, т = \, 2, у = 1, Ыт (N1 и N2 — числа дискретных вихрей, сошедших с задней кромки профиля и с кончика интерцептора соответственно); Аи., Вщ , Сщ — коэффициенты влияния от набегающего потока.
Линейная система уравнений для определения неизвестных плотностей интенсивности вихрей в узловых точках и плотностей интенсивности источников на панелях на заданном шаге по времени будет содержать следующие уравнения:
Здесь первая группа уравнений задает равномерность распределения плотности интенсивности вихря вдоль добавочных панелей; вторая группа уравнений определяет условия непротекания для панелей с распределенными вихрями; третья — условие одинаковости нормальной составляющей скорости для панелей с распределенными источниками — величина этой одинаковой нормальной составляющей скорости); четвертое уравнение определяет условие выполнения теоремы Томсона (/,• — длина у'-й панели, на концах которой (узлах) плотности интенсивности вихрей равны уу и уу-ц)> величина константы Го определяется начальными условиями; последнее уравнение определяет равенство нулю суммарной обильности источников (/у — длина у-й панели с равномерно распределенным с плотностью интенсивности #у источником).
У к\ +1 — У*] + 2> У к2 ~Ук2-Ь
м>„. =0, / = к{ + 2, к2 - 2;
*з-1
./ = *2 + 1
При определении абсолютной величины скорости, генерируемой дискретным вихрем с циркуляцией Г на расстоянии г от него, использовалось выражение
где гв — эффективный радиус вихря (вихревого жгута). Таким образом обычно исключается появление слишком больших значений скорости вблизи вихря. В конкретных расчетах использовался гв =1,5% хорды профиля.
Для движущихся цепочек дискретных вихрей характерна неустойчивость к возмущениям малого порядка. Для повышения устойчивости к возмущениям с длиной порядка расстояния между соседними вихрями в цепочках с примерно равномерным расположением дискретных вихрей применялась следующая процедура. Для цепочки вихрей последовательно, начиная от точки схода вихрей, рассчитывались расстояния между соседними вихрями на заданном шаге по времени. При нахождении положения дискретного вихря на следующем шаге по времени скорость этого вихря, определяемая суммарным влиянием набегающего потока и распределенных и дискретных особенностей, получала дополнительное приращение. Вектор дополнительной скорости н>доп направлен от заданного вихря к более удаленному из двух вихрей, между которыми находится заданный. Величина waon определялась выражением н'дол = к ■ M/At, где AL — разница в расстояниях между заданным вихрем и каждым из двух соседствующих с разных сторон вихрей; At — шаг расчета по времени; к=0,3 — эмпирический коэффициент (& = 0 соответствует отсутствию дополнительной скорости; в случае если бы сумма расстояний между заданным вихрем и двумя соседствующими не изменилась на следующем шаге по времени, то к=0,5 соответствовал бы полному выравниванию расстояний между заданным дискретным вихрем и соседствующими).
Рассмотренная процедура применялась до начала разрушения цепочки дискретных вихрей, когда расстояние между соседними вихрями превышало заданное (в конкретных расчетах 1,5% хорды профиля).
Влияние процедуры повышения устойчивости цепочек дискретных вихрей проиллюстрировано на рис. 3, на котором представлена картина течения для профиля GA(W)-2 при угле атаки а =12° . Расчет с применением процедуры повышения устойчивости проводился до значе-
Рис. з. Влияние процедуры повы- ния параметра х = til /с = 4,1 (безразмер-.
шения УСТОЙЧИВОСТИ цепочек ДИС- ч г -
кретных вихрей. Профиль ное время), после чего расчет был про-
GA(W)-2, а=12° должен без применения рассматриваемой
процедуры. На рис 3 можно видеть, что для момента времени т = 4,5 процесс разрушения цепочек вихрей происходит значительно быстрее.
Для уменьшения времени расчета применялись две процедуры сокращения количества вихрей. А именно, допускалось одно слияние двух дискретных вихрей, сходящих с разных точек отрыва, при достаточном сближении вихрей (в конкретных расчетах при расстоянии между вихрями меньше 1% хорды профиля). При этом новый вихрь располагался посередине между заменяемыми вихрями, циркуляции которых суммировались. Вторая процедура заключалась в слиянии пар вихрей, сошедших с одной точки отрыва, на достаточном удалении вниз по потоку (в конкретных расчетах на расстоянии от задней кромки профиля, большем одной хорды профиля). Вихрь, полученный при слиянии пары вихрей, располагался в центре масс (если под массой понимать циркуляцию). Циркуляции вихрей суммировались. Данная процедура могла применяться неоднократно при достаточном сближении вихрей (при расстоянии между вихрями, меньшем 40% расстояния поперек набегающего потока между точками схода вихрей).
Распределение давления по профилю находилось с помощью интеграла Лагранжа—Коши. Коэффициент давления в /-й контрольной точке
, 2 ~ Зф,-
с„. = 1 - и- 2—^.
Р,
Для упрощения расчета в отрывной зоне позади интерцептора, где скорость потока близка к нулю, производная потенциала по времени
считалась постоянной по зоне. Величина определялась с под/ <Эг
мощью методики, описанной в [12], исходя из величины циркуляции Г1
первого сходящего с кончика интерцептора дискретного вихря:
дФт Г,
~ д1 А/ ’
Зфт „
где — производная потенциала по времени на верхней поверхно-
сти интерцептора вблизи его кончика, А? — шаг расчета по времени.
Суммарные аэродинамические характеристики находились интегрированием полученных эпюр распределения давления.
2. Примеры расчета. Начальные данные при расчете соответствовали безотрывному обтеканию профиля без интерцептора. На первом шаге по времени на профиле появлялся отклоненный интерцептор. Таким образом, рассматривается «мгновенное» отклонение интерцептора.
Отсчет времени проводился в единицах безразмерного времени х = Ш/с, где и — скорость набегающего потока, с — хорда профиля, г — размерное время.
Шаг по времени между сходом очередных дискретных вихрей составлял Ат = 0,01.
На рис. 4 представлены примеры расчета картины течения и эпюры распределения давления для профиля ОА(\У)-2 относительной тол-
х
—— расчет о эксперимент о Со
100 х,°/с
п
.чА|
‘•ЧУ •
Рис. 4. Картина течения и эпюра распределения давления. Профиль
СА(\У)-2, сх=8°, 8„
=60°
щины 13% с интерцептором с относительной хордой 10%, расположенным на 77,5% хорды профиля от его передней кромки. Угол отклонения интерцептора 8ИНТ = 60°. Угол атаки а = 8°. Число Рейнольдса Ле = 2,2 ■ 106. Эпюра распределения давления сравнивается с эпюрой, полученной в стационарном эксперименте [1]. Видно, что позади профиля образуется развитая вихревая дорожка. В распределении давления наблюдается характерная «полка» на верхней поверхности профиля позади интерцептора. : , '
На основе рассчитанных эпюр распределения давления были получены картины установления (при переходе от начального приближения) коэффициентов подъемной силы су, сопротивления сх и момента тангажа ст относительно передней кромки профиля при трех углах отклонения интерцептора 5ИНТ= 10, 40, 60° (рис. 5). Для сравнения на рис. 5 отмечены значения су, сх и ст, полученные в стационарном эксперименте [1]. Соответствие рассчитанных и экспериментальных значений су, сх и ст получается удовлетворительным.
С помощью интегрирования рассчитанного распределения давления по интерцептору был вычислен шарнирный момент интерцептора в зависимости от утла отклонения интерцептора. Рассмотренная модель отрывного обтекания не предполагает возможности исследования случаев малых углов отклонения интерцептора, когда возможно присоединение оторвавшегося с кончика интерцептора потока к верхней поверхности профиля за интерцептором. Однако сравнение с экспериментом [1],
12
-0,5 J
Рис. 5. Установление по времени коэффициентов подъемной силы Су, сопротивления сх и момента тангажа с,„. Профиль ОА(\У)-2, а=8°:
расчет:----— 8ИНТ=]0°:
!
эксперимент: О -•
іт=40°; ■
т~60°;
^инг” * ■
■ ^инх- {*0°
■ 8ИНТ 40°
показало возможность оценки шарнирного момента начиная примерно с 8ИНТ = 10°
(рис. 6).
Для оценки возможности исследования особенностей вихревого следа позади профиля с интерцептором были использованы имеющиеся экспериментальные данные по частотным спектрам мощности турбулентных пульсаций скорости в следе позади профиля Воеіп§ относительной толщины 11,3% с интерцептором с относительной хордой 15,7% с осью вращения, расположенной на 73,3% хорды профиля от его передней кромки, для двух чисел Рейнольдса — И.е = 2,2 • Ю6 [2] и Яе = 0,29 • 106 [3] при бинт = 60°. Пики частот, связанные с прохождением вихревых сгустков, полученные при этих числах Рейнольдса, отличаются мало. Соответствующие числа Струхаля— $1і= 1,16 и 1,13.
Для определения характерных частот по результатам расчетных исследований применялась следующая процедура. На расстоянии от 0,5 до 0,6 хорд: .г профиля от его задней кромки выделялось «окно» (рис. 7), в котором велся подсчет числа Мт) наблюдаемых дискретных вихрей, сходящих либо с кончика интерцептора, либо с задней кромки профиля, в зависимости от времени. После чего находилась величина
эксперимент
расчет
Рис. 6. Коэффициент шарнирного момента интерцептора. Профиль СА(\\0-2, <х=8°
Рис. 7. Определение частоты вихреобразования. Профиль Воете, а=0 , 6ИНТ=60°:
------сход вихрей с задней кромки профиля;-------сход вихрей с кончика
интерцептора; Д — БИэксп =1,16
параметра /?(Sh), характеризующего амплитудно-частотную характеристику (АЧХ) сигнала N(t):
р( Sh) =
где Sh = fc/U — число Струхаля по хорде профиля с; -с, = Дт • / — безразмерное время на i-м шаге по времени; i0 и iy ограничивают диапазон времени, для которого определяется АЧХ.
На рис. 7 представлены зависимости N(x) и p(Sh) для вихрей, сходящих с задней кромки профиля и с кончика интерцептора при 8ИНТ = 60°. Величина /?(Sh) определялась на интервале безразмерного времени от 3 до 10. На АЧХ />(Sh) отмечено число Sh= 1,16, соответствующее частоте вихреобразования, определенной в эксперименте [2].
Таким образом, разработанный метод позволяет оценивать распределение давления при отрывном обтекании профиля с интерцептором при умеренных углах атаки (когда интерцептор не погружается полностью в отрывную зону на верхней поверхности профиля), а также суммарные аэродинамические характеристики и шарнирные моменты интерцептора. Кроме того, с помощью моделирования вихревого следа позади профиля с интерцептором можно оценивать характерные частоты вихреобразования при больших углах отклонения интерцептора.
ЛИТЕРАТУРА
1. Wentz W. Н., Jr. Wind-tunnel tests of the GA(W)-2 airfoil with 20% aileron, 25% slotted flap and 10% slot lip spoiler // Wichita State University,
NASA CR-145139. - 1977.
2. W e n t z W. H., Jr., Ostowari C., Seetharam H.C. Effects of design variable on spoiler control effectiveness, hinge moments and wake turbulence//AIAA-81-0072. — 1981.
3. Bodapati S., Karamcheti K., Mack M. D. Basic stadies of the flow fields of airfoil-flap-spoiler system//AIAA-82-0173. — 1982.
4. Consigny H., Graveile A., Molinaro R. Aerodynamic characteristics of a two-dimensional moving spoiler in subsonic and transonic flow //
J. Aircraft.— 1984. V. 21, N 9.
5. Jaeggy В. C., Czichowsky W., Meyer P., Koerber G., Kauffmann R. Contribution a l'etude experimentale et theorique de 1 ecoulement autour d'un profil d'aile muni d'un spoiler//L'Aeronautique et l'Astronautique. — 1985, N 112.
6. Costes М., Graveile A., Philippe J. J. Steady and unsteady investigation of spoiler and flap aerodynamics in two-dimensional subsonic flows // ICAS-86-2.9.2.— 1986.
7. Wo о d s L. C. Theory of airfoil spoilers // British R&M 2969.—1956.
8.Barnes C. S. A developed theory of spoilers on airfoils// British ARC CP887. - 1965.
9. Parkinson G. V., Tam Do о P., Bernier R. A prediction method for spoiler performance//SAE 770459. — 1977.
10. P f e i f f e r M. J., Z u m w a 11 G. W. A computational model for low speed flows past airfoils with spoilers// AIAA Paper N 81-0253. — 1981.
11. Той H. B., Hancock G. J. Inviscid theory of two-dimensional airfoil/spoiler configurations at low speed// Aeronautical Journal.— October, 1987.
jV(t, )cOs(2tI Sh-C;)
t—10
Ar(t,)sin(2n Sh-t/)
l=10
12. Белоцерковский С. М. Математическое моделирование плоскопараллельного отрывного обтекания тел.— М.: Наука. — 1988.
13. Флетчер К. Вычислительные методы в динамике жидкости.— М.: Мир. - 1991.
14. Л о й ц я н с к и й Л. Г. Механика жидкости и газа.— М.: Наука. —
1987.
15. Бондарев Е. Н., Дубасов В. Г., Рыжов Ю. А., Свир-щевский С. Б., Семенчиков Н. В. Аэрогидродинамика,— М.: Машиностроение. — 1993.
Рукопись, поступила 3/11996 г.
