CC BY
90
УЧЕНЫЕ ЗАПИСКИ Ц АГИ
Т о м XX 198 9 М3
УДК 533.6.011.32:629.7.025.1
РАСЧЕТ ИНДУКТИВНОГО СОПРОТИВЛЕНИЯ КРЫЛА С ПРОИЗВОЛЬНОЙ ДЕФОРМАЦИЕЙ
Ф. И. Ганиев, Нгуен Дык Кыонг
Предложена методика расчета индуктивного сопротивления с подсасывающей силой и без нее крыла с гладкой деформацией, движущегося в среде идеальной несжимаемой жидкости. Расчет основан на методе дискретных вихрей в нелинейной постановке с выполнением условия Чаплыгина— Жуковского на задней и боковых кромках. Вихревая пелена и циркуляция вихревых отрезков определяются последовательными приближениями. Подсасывающая сила рассчитывается по первым вихревым отрезкам. Результаты расчета сравниваются с точными решениями в предельных случаях и с результатами расчета крыльев, представляющих аналитически заданные криволинейные поверхности, а также с экспериментальными данными для моделей крыльев малого удлинения. Показана возможность получения аэродинамического качества выше его максимального значения, определяемого по линейной теории.
Расчет индуктивного сопротивления несущей поверхности с помощью метода дискретных вихрей рассмотрен в линейной постановке в работах [1, 2], в нелинейной постановке — в [3]. В работе [2] показано, что расчет индуктивного сопротивления искривленных поверхностей в линейной постановке с непосредственным интегрированием проекций нормальных нагрузок (т. е. по «ближнему полю» скоростей у поверхности крыла), а не с помощью теоремы количества движения, (т. е. по «дальнему полю» скоростей в плоскости Треффтца), сильно зависит от количества вихрей по хорде N. При этом не было ясно, достаточно ли довести N до 32, как было сделано в этой работе. В зарубежных работах также отмечается медленная сходимость расчета сопротивления по «ближнему полю» по сравнению с расчетом по «дальнему полю» [4]. Это объясняется тем, что индуктивное сопротивление представляет собой, как правило, небольшую разность близких величин. Если даже они получены с приемлемой точностью, то их разность может иметь значительные относительные погрешности. Результаты расчета по «дальнему полю» зависят только от распределения суммарной циркуляции по размаху, поэтому имеет место сходимость при увеличении числа вихрей.
При решении задачи в нелинейной постановке расчет индуктивного сопротивления по «дальнему полю» существенно усложняется, так как при этом необходимо учитывать реальную форму вихревой пелены. По
этой причине в настоящее время расчет этих характеристик проводится только по «ближнему полю». При решении нелинейных задач количество дискретных вихрей, моделирующих обтекание несущей поверхности, ограничивается возможностями ЭВМ. Поэтому получаемые значения индуктивного сопротивления в нелинейных задачах часто имеют значительные погрешности.
В настоящей статье предлагается методика расчета индуктивного сопротивления с подсасывающей силой и без нее в нелинейной постановке по «ближнему полю», которая позволяет получить удовлетворительную точность при существующих ограничениях на возможности ЭВМ. Задача рассматривается в рамках нелинейной теории стационарного безотрывного обтекания тонкой несущей поверхности несжимаемым потоком идеальной жидкости при помощи метода дискретных вихрей [3].
1. Методика расчета. Задано тонкое крыло с криволинейной срединной поверхностью а в связанной системе координат Охуг (рис. 1). Поверхность а может иметь произвольную кривизну и крутку, но должна быть гладкой. Крыло движется равномерно в безграничной жидкости с поступательной скоростью ио под некоторым постоянным углом атаки а без скольжения. Предполагается, что обтекание плавное, на задней и боковых кромках имеет место тангенциальный сход потока (выполняется условие Чаплыгина — Жуковского). Требуется найти коэффициенты сопротивления среды с подсасывающей силой Сх+ и без нее
С лг_, а также коэффициент подъемной силы Су.
Задача решается методом дискретных вихрей. Из системы линейных алгебраических уравнений определяются т неизвестных циркуляций вихревых отрезков на одной половине поверхности а от плоскости хОу. Положение свободных вихревых отрезков, моделирующих пелену, определяется совместно с циркуляциями путем последовательного уточнения [3]. По найденным значениям циркуляций определяются безразмерные нагрузки Др [3].
Из рис. 1 видно, что коэффициент нормальной силы, направленной по оси Оу, рассчитывается без подсасывающей силы по интегралу:
уУ
х
Рис. 1
(1)
(а)
где 5 — характерная площадь крыла.
Коэффициент продольной силы без подсасывающей силы определяется по формуле
СТ=ЛЯ Л,-
cos (ti, x) do . ^2)
(з>
Переходя в скоростную систему координат, получим коэффициенты Су_ и Сх_ без подсасывающей силы:
Сх- = Сп sin а + Ст cos а , | ^
Су_ = Сп cos а — Ст sin а . J
Подсасывающая сила Q рассчитывается по циркуляциям первых (ближайших к передней кромке) вихревых отрезков [1]. Хотя в [1] методика разрабатывалась только для линейных задач, но основное соотношение, связывающее элемент подсасывающей силы dQ и характер изменения интенсивности вихревого слоя на данном участке dL передней кромки L, получено без линеаризации. Поэтому воспользуемся им также и для решения задач в нелинейной постановке. Коэффициенты соответствующих составляющих подсасывающей силы в скоростной системе координат будут:
(4)
1 г / \ Г / \
Сп =~7 [cosa cos IdQ, л;) dQ-f sin a J cos (dQ, y) aQl,
* щ <£>
CQy = —[— sin a J* cos (dQ, x) dQ + cos a J cos (dQ, y) dQ], I
(L) (I.) )
где q — скоростной напор, а вектор dQ направлен перпендикулярно к элементу передней кромки dL и к нормали поверхности о в данной точке.
Коэффициенты подъемной силы и сопротивления с подсасывающей силой будут
Су,=Су -j- Cq ,
Сх+ = СХ_ + Сдх (
Как правило, | С<эу [ <С! | Су+ | , поэтому в дальнейшем величины Cv- и Су+ не будем различать и обозначим через Су.
Как было выше указано, расчет по соотношениям (1) — (5) может привести к значительным погрешностям. Для повышения точности расчета в предлагаемой методике приняты следующие меры.
Во-первых, численное решение соответствующего сингулярного интегрального уравнения при помощи схемы расположения вихрей и контрольных точек по хорде по закону косинуса «COS» равномерно сходится к точному решению [5], а при расположении по схеме «1/4» численное решение тоже сходится, но неравномерно (в окрестности передней кромки). Поэтому была принята схема «COS».
Во-вторых, вместо интегрирования по формуле прямоугольников при расчете по соотношениям (1) — (2) принято интегрирование по хорде с помощью параболического сплайна с аналитическим вычислением несобственного интеграла у передней кромки, а по размаху шаг инте-
грирования уменьшен на порядок по сравнению с расстоянием между соседними продольными вихрями.
В-третьих, более точное определение направляющих косинусов нормалей в контрольных точках и в точках расчета нагрузок (в центрах отрезков) также существенным образом уточняет численный расчет коэффициентов Сх+ и Сх... При аналитическом задании формы поверхности а получение точных значений направляющих косинусов не вызывает затруднений. Однако в этом случае возникают трудности в подготовке исходных данных для решения большинства практических задач. Поэтому целесообразно задание криволинейной поверхности 0 в виде совокупности числовых координат (опорных точек). Характерные точки вихревой решетки, в общем случае, не совпадают с опорными. Необходимо также вычислить направляющие косинусы нормалей. Возникает потребность в аппроксимации дискретно заданной поверхности.
Принята следующая аппроксимация криволинейной поверхности. В каждом сечении кривая дужка аппроксимируется одномерным кубическим сплайном. Далее, через каждые шесть соседних точек вихревой решетки строится одна «порция» криволинейной поверхности второго порядка:
у = Ах2 + Bz2 + Cxz 4- Dx -f Ez +F ,
где неизвестные коэффициенты находятся по условиям прохождения через эти шесть точек. Вычисляются направляющие косинусы нормали:
/\ . —ду/дх /\ 1 / ч —dyjdz
cos (в, я) =-------- , cos (п, у) = — , cos (п, z) =-------- ,
II ЯII II л II я
м= у 1+ (-£)+(-£)’• -£=2Лх+с*+°- -£-™г+сх + Е.
Здесь значения А, В, С, Д Е соответствуют той «порции» поверхности, которой принадлежит точка искомой нормали.
2. Исследование методики расчета. Время расчета варианта крыла для одного значения угла атаки с числом моделирующих вихрей т — 90 на ЭВМ ЕС-1060 составляет около 10-5-15 мин.
Сопоставление результатов численного расчета с известным точным решением для дужки окружности с кривизной /=15% приведено в табл. 1. Буквами «а» и «б» обозначено соответственно интегрирование по прямоугольникам и уточненное интегрирование.
Таблица 1
Характеристика с,+ су
Точное решение [6] 0 — 1,882
Численный расчет „1/4“ а —0,0114 -0,0062 1,839
б —0,0284 -0,0232 1,867
У"*»ч II о о о ,COS“ а 0,0012 0,0012 1,753
б 0,0025 0,0025 1,766
Здесь прямоугольное крыло с удлинением Х= 1000 моделировалось вихревой решеткой (5 вихрей по полуразмаху и 9 вихрей по хорде) по схемам «1/4» и «COS». По значению Су схема «COS» несколько хуже, чем «1/4», но для величины Сх+ Принятая схема дает повышение точности на порядок.
У крыльев большого удлинения при плавном и стационарном обтекании решение по нелинейной теории с уменьшением деформации крыла и угла атаки должно приближаться к результатам линейной теории. Расчеты показывают практическое совпадение распределенных и суммарных характеристик, полученных по линейной и нелинейной теориям, у крыльев с удлинением Х>5, со стреловидностью от нуля до 45° и Су<с0,4. Следовательно, в указанных пределах точные решения линейной теории могут быть использованы для проверки численных расчетов по нелинейной теории.
Таблица 2
Профиль Плоская пластинка, а = 2,5° Параболическая дужка а = 0
Характеристика Сх х+ с,_ су Сх х+ С,
Линейная теория [5] 0 0,0119 0,2733 0 0 . 0,392
5X9 вихрей X = 1000 „1/4“ а б 0,0002 0,0008 0,0119 0,0128 0,2712 0,2939 —0,0012 —0,0012 -0,0009 —0,0009 0,398 0,406
• COS* а б 0,0003 0,0003 0,0117 0,0119 0,2678 0,2743 0,00005 0,0001 0,00005 0,0001 0,379 0,382
В табл. 2 сравниваются точные решения линейной теории и результаты численного расчета тонких профилей в случае безударного входа потока (симметричная параболическая дужка с небольшой кривизной / = 3,1%) и при отсутствии безударного входа потока (плоская пластинка под небольшим углом атаки). Видно, что схема «1/4» дает даже отрицательный коэффициент сопротивления значительной величины у параболической дужки, в то время как схема «COS» при принятой методике интегрирования дает результаты весьма близкие к точным значениям. Что касается значения Су, то наиболее точным оказывается расчет по схеме «1/4» при интегрировании по формуле прямоугольника. Известно [1], что при этом получается точное значение Су для пластины при любом числе вихрей. В работе [7] доказано, что при схемах «1/4» и «COS», если интерполировать полиномом определенной степени с особенностью, можно получить точную распределенную нагрузку для пластины в однородном и неоднородном потоке. В настоящей работе определение распределенной нагрузки в нелинейной постановке выполнено по методике [3] в точках, где расположены дискретные вихри. С целью обеспечения гибкости программы расчета интерполяция нагрузок и их интегрирование производится с помощью параболического сплайна. Поэтому наблюдаются отличия полученных значений Су от его точных значений по линейной теории, особенно по схеме «1/4». В табл. 3 по-
казана сходимость численного расчета фактора индуктивного сопротивления В+ = іс X Сх+/С2у для плоского эллиптического крыла в зависимости от числа вихрей.
Таблица 3
Линейная теория В+ = 1,0
Численный текст N — & т 40 48 56 64
В+ 1,06 1,03 1,02 1,01
Отношение поперечной оси у эллипса к его продольной оси равно 4. Угол атаки а = 2°. В численном расчете эллипс был подрезан на концах на величину 1% поперечной оси. Число вихрей по хорде N = 8 и фиксировалось. Видно, что для плоского крыла практическая сходимость достигается уже при т = 50-^60.
Были проведены также расчеты слабоизогнутых поверхностей различных форм в плане с эллиптическим законом распределения циркуляции по размаху. При этом исследовались случаи с безударным входом потока у передней кромки и без него. Во всех расчетах наблюдалась удовлетворительная сходимость к точному решению линейной теории при т = 90-г-110. В качестве примера на рис. 2 приведены зависимости численных значений фактора индуктивного сопротивления без подсасывающей силы В_=лАС*-/ Су у слабоизогнутых поверхностей (Су-^0) с эллиптическим законом распределения циркуляции по размаху и с безударным входом потока у передней кромки. В данном случае по линейной теории В+=В-— 1,0. Видно, что схема «COS» вместе с уточненным интегрированием дает удовлетворительные результаты уже при т = = 80-^100. Уточненное интегрирование, вскрывая погрешность расчета распределенных нагрузок по схеме «1/4», еще больше отдаляет численные результаты от точного решения.
Наряду со сравнением с линейной теорией для слабоизогнутых поверхностей необходимо исследовать практическую сходимость численных расчетов, проводимых по аппроксимированным и точным геометрическим данным у поверхностей с достаточно большой кривизной и крут-
095
0,90
0,8$
50
—
f
_ 1
J ~
t к
/ Л
'l 1 ' 1 і I. —
100
О
50
100
Рис. 2. Зависимость величины В_ от принятой схемы расчета: /—схема «COS» и уточненное интегрирование; 2—схема «COS» и интегрирование по прямоугольникам; 3—схема «1/4» и интегрирование по прямоугольникам; 4— схема «1/4» и уточненное интегрирование
кой.В качестве примера на рис. 3 приведены результаты такого исследования для параболоидальной поверхности.
(У -Н Уо)2 г* _х+х0 № ‘ с2 х, ’
у2
где Ь= 1, с = 2, *1 = 3, г/о=0,2, Хо=*1 — . Поверхность ограничена плос-
ь2
костями 2— ±0,85 с, х = х1 и у = у0. Она представляет собой тонкое крыло со значительной кривизной и круткой. Удлинение Я =1,43, сужение г] =
/и0
Рис. 3. Результаты расчета:
1—по точным геометрическим даннцм; 2—по аппроксимированным исходным данным
= 3,42, угол между осью параболоида и вектором невозмущенной скорости составляет 24,4°. Из рисунка видна удовлетворительная сходимость результатов при т» 100.
3. Сопоставление с физическим экспериментом. Для подтверждения достоверности методики необходимо сопоставление с результатами продувок моделей крыла с заведомо значительным влиянием боковых вихревых образований. С этой целью в аэродинамической трубе был проведен эксперимент с 4 моделями прямоугольного крыла с удлинением А, = 0,5. Размеры моделей 300X600 мм. Скорость потока равна 34 м/с, что соответствует числу Рейнольдса Ке=1,4-106. Начальная турбулентность е = 0,9%. Сечение рабочей части представляет собой восьмигранник с диаметром вписанной окружности 2,25 м. Из четырех моделей две плоские и две изогнутые. Плоская модель выполнена с острыми передней и задней кромками и относительной толщиной с=1,5%. После продувок на ее переднюю часть был установлен обтекатель с заполнителем, что придавало носку профилированную форму с с = 3%. Изогнутая модель
выполнена в виде «корыта» с толщиной с = 0,4% и острыми передней и задней кромками. После продувок на ее переднюю часть также был установлен аналогичный обтекатель для получения профилированной формы носка, при этом с = 3%. Во всех случаях боковые торцы крыльев были плоскими (незаостренными).
В программу эксперимента включались весовые испытания и исследование физической картины обтекания.
На рис. 4, а показаны экспериментальные данные для плоских моделей с острым и профилированным носком. Здесь приведены также расчетные кривые с подсасывающей силой и без нее. Видно, что экспериментальные данные у модели с острым носком весьма близки к кривой без подсасывающей силы вплоть до су~0,6 (а*»20°). На больших углах атаки наблюдался сильный отрыв с передней кромки, что не моделировалось в расчетах. Поэтому при Су>0,6 имеет место расхождение
Рис. 4. Сопоставление экспериментальных и расчетных данных:
1—с учетом подсасывающей силы; 2—без подсасывающей силы
эксперимента и расчета. У плоской модели с профилированным носком на небольших углах атаки (а<10° или Су<0,3) экспериментальные точки находятся близко к кривой' с подсасывающей силой, далее они приближаются к кривой без подсасывающей силы. Из рис. 4, а также видно, что экспериментальные точки лежат выше, чем кривая с подсасы-
На рис. 4, б приводятся результаты продувок изогнутых моделей с профилированным и острым носком. За индуктивное сопротивление при обработке опытных данных принималась разность С, —{Сх0)у_0, где (^1)7= о было взято из эксперимента у плоских моделей, которые вмес-
те с сопротивлением узлов крепления модели составило 0,011. Срединная поверхность у этих моделей была одной и той же и имела аэродинамическую крутку: кривизна в центральных сечениях / = 5,8%, а на концах уменьшается до нуля. Поверхность была специально подобрана таким образом, чтобы на углах а» 10° обеспечивался безударный вход потока. Исследование поля скоростей вблизи передней кромки модели с острым носком с помощью щупа с волосинками подтвердило, что на а=Ю—11° наблюдалось примерно одинаковое разветвление струек по обеим сторонам острого носка по всей передней кромке. На рис. 4, б также приведены расчетные зависимости Сх^—/{су) полученные по геометрическим данным реальных моделей (со всеми погрешностями изготовления). Видно хорошее согласование результатов расчета и эксперимента, особенно вблизи углов атаки а~10° (Су^0,4). На больших углах атаки а>18° (Су>0,7) у модели с профилированным носком теряется подсасывающая сила и ее экспериментальные точки приближаются к кривой без подсасывающей силы. В отличие от плоских моделей, у изогнутых моделей как с профилированным носком, так и с острым носком экспериментальные данные мало отличаются от расчетных вплоть до Су*&\,0 (а~25°). Это объясняется тем, что полный отрыв потока с передней кромки у изогнутых моделей наблюдался в эксперименте лишь при а>30° и обтекание этих моделей хорошо описывается безотрывной схемой в широком диапазоне углов атаки. На малых углах атаки (Су<0,4) у модели с профилированным носком подсасывающая сила сохранялась вплоть до а—7° (С^О). В этом диапазоне углов атаки у модели с острым носком наблюдается некоторое расхождение теории с экспериментом, особенно при а = 0-^4° (С„ = 0,1—0,2). Это, по-видимому, объясняется наличием местного отрыва на нижней стороне острого носка и уменьшением пика разрежения.
На рис. 4,6 показана оптимальная (огибающая) поляра, получен-
ная по линейной теории |СЛ;+= . В экспериментах удалось у моде-
лей тонкого крыла при Су = 0,4-ь0,6 реализовать величину индуктивного сопротивления на 30н-40%! меньшую, чем минимально возможная по линейной теории.
Таким образом, разработана эффективная методика расчета индуктивного сопротивления в нелинейной постановке, позволяющая выявить дополнительный резерв улучшения поляры крыла на дозвуковых скоростях полета.
Авторы выражают благодарность Караску А. А. за методическое руководство экспериментом.
вающей силой, полученная по линейной теории
1. Белоцерковский С. М., Скрипач Б. К. Аэродинамические производные летательного аппарата и крыла при дозвуковых скоростях.— М.: Наука, 1975.
1. Г'л ушков Н. Н. О точности расчета аэродинамических характеристик тонких крыльев и профилей методом дискретных вихрей. — Ученые записки ЦАГИ, 1982, т. 13, № 3.
3. Белоцерковский С. М., Ништ М. И. Отрывное и безотрывное обтекание тонких крыльев идеальной жидкостью. — М.: Наука, 1978.
4. Miranda L. R. Application of computational aerodynamics to airplane design. — J. of Aircraft, 1984, vol. 21, N6.
5. Белоцерковский С. М., Л и ф а н о в И. К. Численные методы в сингулярных интегральных уравнениях и их применение в аэродинамике, теории упругости, электродинамике. — М.: Наука, 1985.
6. Бетчелор Д. Введение в динамику жидкости. — М.: Мир, 1973.
7. Баскин В. Э., Каплан В. С. Дискретные вихревые схемы плоской пластины в неоднородном потоке. — Ученые записки ЦАГИ, 1986, т. 17, № 1.
Рукопись поступила 20/VII 1987
