УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XVIII 19 87
№ 1
УДК 629.735.33.015.4 :533.6.013.42
РАСЧЕТ ОКОЛОЗВУКОВОГО ОБТЕКАНИЯ КОМБИНАЦИИ КРЫЛО—ФЮЗЕЛЯЖ С УЧЕТОМ влияния УПРУГОСТИ КОНСТРУКЦИИ КРЫЛА
В. И. Савицкий
На основе метода коэффициентов влияния [1] в рамках околозвуковой теории малых возмущений с применением метода последовательной релаксации разработан итерационный метод учета влияния упругости конструкции крыла на околозвуковое обтекание комбинации крыло — фюзеляж. Проведены исследования сходимости итерационного процесса решения задачи, найдены его оптимальные параметры. Приведены результаты расчетов характеристик статической аэроупругости крыла по разработанной методике.
В настоящее время для определения влияния упругости конструкции на суммарные аэродинамические характеристики самолета в околозвуковом диапазоне чисел М<*, в практических целях, как правило, проводятся эксперименты на упругоподобных моделях. В очень редких случаях удается провести испытания дренированных упругоподобных моделей. Также используется полуэмпирический метод расчета: распределение коэффициентов давления на жесткой модели определяется экспериментально, а поправки на влияние упругости находятся из расчета с помощью линеаризованных методов решения задач статической аэроупругости [1]. Необходимость в получении экспериментальных данных для выполнения такого расчета существенно ограничивает возможности его использования.
Успешное развитие численных методов анализа пространственных околозвуковых течений в последние годы явилось основой разработки методов расчета характеристик статической аэроупругасти элементов самолета на околозвуковых режимах полета. Первые результаты в этом направлении были получены за рубежом. В работе [2] предложен метод расчета околозвукового обтекания изолированного крыла с учетом влияния упругости конструкции, основанный на итерационном соединении -двух этапов вычислений: расчета околозвукового обтекания изолированного крыла в рамках нелинейной околозвуковой теории малых возмущений и расчета упругих деформаций крыла по линейному методу. В работе [3] для расчета околозвукового обтекания изолированного упругого крыла подобная итерационная процедура составляется
с использованием метода расчета околозвукового обтекания крыла на основе численного решения полного уравнения для потенциала по программе РЬО-22 [4], имеющей в СССР аналог-[5], и с использованием методики расчета упругих деформаций крыла на основе линейного метода коэффициентов влияния.
В настоящей статье впервые предлагается метод расчета околозвукового обтекания комбинации упругое крыло — жесткий фюзеляж. Так же, как и в двух работах, упомянутых выше, здесь используется метод итераций. В качестве аэродинамической теории предлагается использование нелинейного метода расчета околозвукового обтекания комбинации крыло—фюзеляж в рамках околозвуковой теории малых возмущений [6], а для расчета упругих деформаций крыла применяется линейный метод коэффициентов влияния. Использование околозвуковой теории малых возмущений, как показывают результаты работы [6], обеспечивает необходимую для решения задачи статической аэроупругости точность расчета околозвукового течения, значительно упрощает вычисления и сокращает время счета на ЭВМ по сравнению с решением полного уравнения для потенциала.
1. Основные предположения. Рассматривается установившееся околозвуковое обтекание комбинации упругое крыло — жесткий фюзеляж. Упругое крыло предполагается жестко заделанным по борту фюзеляжа.
При описании упругих свойств конструкции крыла принимается существование базовой плоскости, которая выбирается близкой к срединной поверхности крыла. Предполагается, что упругие смещения точек конструкции направлены по нормали к базовой плоскости, деформации малы, а матрицы коэффициентов упругого влияния можно рассматривать постоянными при изменении формы крыла. При выполнении расчета упругих деформаций срединная поверхность крыла разбивается на панели трапециевидной формы в плане. Считается, что суммарные аэродинамические силы на каждой панели направлены по нормали к базовой плоскости и приложены в геометрических центрах тяжести площадей панелей. Инерционная разгрузка крыла для расчетной перегрузки учитывается при задании исходной геометрии крыла.
При выполнении расчетов околозвукового обтекания комбинации крыло—фюзеляж предполагается, что скачки уплотнения имеют малую интенсивность и изменением энтропии в частицах, пересекающих их фронты, можно пренебречь. Это предположение позволяет описывать течение газа уравнением для потенциала скорости. В рамках околозвуковой теории малых возмущений условие непротекания на поверхности крыла выполняется в плоскости у = О (см. рис. 1 , а). Эта плос-
кость совпадает с базовой плоскостью, вводимой для расчета упругих деформаций крыла. При задании граничного условия на поверхности вихревой пелены, сходящей с задней кромки крыла, сворачивание вихревой пелены не учитывается, считается, что вихревая пелена лежит в плоскости у=0.
2. Описание метода расчета. Основное уравнение околозвуковой теории малых возмущений имеет вид
А 9 XX + В (?уу 4* ¥ гг) — 2М^о <РХ <?хг = 0 , (1)
где
А =
1—М*о —(т+ 1)mL*,-!±^M2o<p2*| ; £ = [ 1 — (т — 1) ML cpj ;
Moo — число М набегающего невозмущенного потока; у — показатель адиабаты, а ф = ср(л:, у, г) — потенциал возмущений скорости. Это уравнение интегрируется численно при наложении соответствующих граничных условий (см. [6]) с помощью метода последовательной верхней линейной релаксации с использованием неконсервативной конечно-разностной схемы [7]. В процессе выполнения расчета одного режима околозвукового обтекания (жесткой) комбинации крыло—фюзеляж проводится 150 итераций метода последовательной релаксации на грубой расчетной сетке (31X21X15) соответственно направлениям Ox, Оу, Oz (см. рис. 1, а) и 100 итераций — на удвоенной расчетной сетке 62Х X42X30). Учет влияния упругости конструкции на околозвуковое обтекание комбинации крыло — фюзеляж в настоящей статье проводится одновременно с выполнением, итерационного процесса работы [6] на этапе осуществления расчетов на мелкой (удвоенной) расчетной сетке посредством итераций, которые будем называть глобальными. Каждая глобальная итерация включает несколько итераций метода последовательной верхней линейной релаксации для уточнения аэродинамического решения, описывающего потенциальное течение около комбинации крыло — фюзеляж и задающего распределение аэродинамических нагрузок на поверхности крыла, и расчет упругих деформаций крыла под воздействием этих аэродинамических нагрузок.
Для учета влияния упругости в граничное условие непротекания на поверхности крыла при расчете его обтекания совместно с фюзеляжем включены члены, описывающие изменения наклонов участков поверхности крыла вследствие его упругих деформаций. Эти члены определяются по результатам расчета упругих деформаций крыла на предыдущей глобальной итерации. Глобальные итерации повторяются до достижения сходимости решения задачи об обтекании комбинации упругое крыло — жесткий фюзеляж с одновременным окончательным определением геометрии этой комбинации..
Сходимость итерационного процесса оценивалась величиной Rescp=<p<n+1)—ф(п), где ф<п) — значение потенциала возмущений на п-й итерации, а черта сверху обозначает усреднение по расчетным точкам на поверхности крыла. Итерационный процесс считался сошедшимся к решению задачи при выполнении условия Res(p<10~5.
Граничное условие непротекания, записанное в приближении око-
дозвуковой теории малых возмущений с учетом влияния упругих деформаций крыла, имеет вид
где г/=Уп(х,г)—уравнение исходной поверхности жесткого крыла; а — угол атаки крыла, задаваемый углом атаки его бортового сечения;. А|а=Да(д:, г) —изменение угла атаки средней линии крыла вследствие его упругих деформаций; верхние индексы « + » и «—» относятся к верхней и нижней поверхностям крыла соответственно; верхний индекс К указывает номер глобальной итерации.
В предположении, что упругие деформации крыла малы, линейный метод коэффициентов влияния позволяет выписать матричные соотношения
{Д«} = [01{<?}; (3)
{^“[СЧКП, (4)
которые предполагают разбиение расчетной схемы крыла на п элементарных участков. Суть обозначений в формулах (3) и (4) будет ясна из дальнейшего изложения.
Рассмотрим в качестве расчетной схемы проекцию срединной поверхности крыла на плоскость хОг, разделенную на п трапециевидных элементов, основания которых параллельны оси симметрии комбинации крыло—фюзеляж (см. рис. 1, б). Аэродинамическая нагрузка (¿т , действующая на т-й элемент расчетной схемы крыла при проведении К-й глобальной итерации учета влияния упругости на околозвуковое обтекание комбинации крыло — фюзеляж имеет вид:
с& = я $ ьр* аэ, (5>
здесь интегрирование ведется по т-му элементарному участку расчетной схемы крыла; q — значение скоростного напора; Арк=Арк(х, г) — значение перепада давления в точке с координатами (х, г) на т-и расчетном элементе перед проведением расчета упругих деформаций крыла на /С-й глобальной итерации. Как следствие матричных формул (3) и (4), получим соотношения
Д 4, 2 ст, I дГ; (б>
/-1
= £ с1,1 <2?, (7)
г=1
где Аато — приращение местного угла атаки т-го элемента, №т — прогиб т-го элемента расчетной схемы крыла, компонент матрицы упругого влияния Ст, I — приращение местного угла атаки т-го элемента,
обусловленное единичной поперечной силой, приложенной в контроль-
ной точке /-го элемента, а компонент матрицы упругого влияния Ст, /— прогиб т-го элемента расчетной схемы крыла, обусловленный единичной поперечной силой, приложенной в контрольной точке /-го элемента. Использование формул (6), (7) позволяет завершить выполнение К-й глобальной итерации и с видоизмененными граничными условиями (2) перейти к осуществлению (/С+1)-й глобальной итерации. Для обеспечения надежной сходимости глобальных итераций предусмотрена возможность использования нижней релаксации. В формуле (2) вместо изменения угла атаки (Да)к, рассчитанного на предыдущей глобальной итерации, используется линейная комбинация вида
(1 - “о) (Да)*-1 + ша (Да)АГ , (8)
где сое —0,8—0,9.
При численной реализации изложенного метода расчета соотношение (2) использовалось в конечно-разностном виде
1дл\* 1 К '(дЛ)
С II II >> /\ к \ дх /<х=о_ /, к
а— (Да),-, к , (9)
где индексы & — номера узлов расчетной сетки метода расчета [6] по направлению осей координат Ох и Ог соответственно.
По узлам расчетной сетки (/,£), описывающим плоскость крыла, проводится численное интегрирование соотношения (5). Так как расчетная схема крыла линейного метода аэродинамических коэффициентов влияния отличается от сеточного представления поверхности крыла в методе расчета [6], — в первом из них крыло разбивается на п=100 участков, а во втором — поверхности крыла соответствует г= (г'0Х&о) = = (34X22) =748 узлов расчетной сетки, — проводилась сплайновая интерполяция результатов расчетов по формулам (6) и (7). Найденные в результате интерполяции значения (Д«)£ к в расчетном узле (г, £) затем использовались в соотношении (9).
Основной целью первых расчетов по предлагаемой методике и разработанной программе, результаты которых представлены далее, является выбор определяющих параметров итерационной процедуры для обеспечения наиболее быстрой сходимости и экономии времени счета. Эти параметры следующие: — коэффициент релаксации для глобальных итераций в формуле (8); — число итераций релаксацион-
ного метода [6], выполняемых на мелкой расчетной сетке на одной глобальной итерации; Ыв — число глобальных итераций. Из опыта использования разработанной программы расчета [6] известно, что для получения сошедшегося решения уравнения трансзвуковой теории малых возмущений при расчете околозвукового обтекания жесткой комбинации крыло — фюзеляж требуется выполнить 100 итераций на мелкой сетке. Решено было и при расчете с учетом влияния упругости конструкции сохранить общее число итераций релаксационного метода М0, связанное с параметрами и Ыв соотношением
N0 = NF■No, (Ю)
равным 100. Таким образом, осталось два независимых определяющих параметра м<з и оптимальный выбор которых можно сделать на основе численных экспериментов. В качестве критерия сходимости при этом использовалось как требование уменьшения поправки потенциала
возмущений Res ф до заданного значения, так и требование стабилизации текущего значения подъемной силы крыла.
Во всех приведенных ниже примерах результатов расчета по разработанной методике также выполнялась нулевая глобальная итерация — первый расчет упругих деформаций крыла проводился после выполнения первой итерации релаксационного метода [6] на мелкой расчетной сетке. Остальные глобальные итерации выполнялись в соответствии с соотношением (10).
Матрицы коэффициентов упругого влияния, входящие в формулы (3), (4) и (6), (7), предполагались заданными и в процессе расчета считывались с магнитного носителя (диска). Эти матрицы были предоставлены автору Д. Д. Евсеевым.
Разработанная программа реализована на языке FORTRAN-IV на ЭВМ типа ЕС 10-60 и требует затрат машинного времени на 10—15% больше, чем исходная программа расчета околозвукового обтекания жесткой комбинации крыло — фюзеляж.
3. Результаты расчетов с использованием итерационного метода. Далее приводятся результаты расчетов по новой программе для комбинации упругое крыло — жесткий фюзеляж, близкой к компоновке современного пассажирского самолета. Будем рассматривать две зависимости крутки крыла ф(г) по его размаху: крейсерскую крутку крыла, задаваемую при получении чисто аэродинамических характеристик жесткого крыла, и стапельную крутку или крутку реального крыла самолета, полученную из крейсерской с некоторым учетом влияния упругости крыла. В дальнейшем изложении крыло с первым законом крутки будем называть крылом «К», а крыло со стапельной круткой будем обозначать крылом «С».
Динамика сходимости итерационного процесса разработанного метода расчета показана на рис. 2 и 3. На рис. 2 показана зависимость
Res?
0
25
50
15
Nít
— жесткое крыло „С" :Nr =0 ;cv=B531
— - •• „К" \ 0; 0.W
— упругое
Рис. 2
■ тестнве крылв„К”\Ы£=0]Ыр=Н10; сух=0,441
------упругое
5\ 20 \
20
5
Gyynp 0,349, цЧВОООН/м'-
Рис. 3
поправки потенциала возмущений Res <р от числа выполненных итераций метода последовательной верхней линейной релаксации [6] для различных параметров глобальных итераций учета влияния упругости конструкции на околозвуковое обтекание комбинации крыло — фюзеляж. Из рисунка видно, что учет влияния упругости крыла улучшает сходимость итерационного процесса при расчете обтекания комбинации крыло — фюзеляж. Скорость сходимости итерационного процесса практически не зависит от параметров глобальных итераций NG, NF. По-видимому, этот вывод будет справедлив для любой комбинации упругое крыло — жесткий фюзеляж с крылом прямой стреловидности.
На рис. 3, а показана зависимость текущего значения коэффициента подъемной силы крыла от числа итераций метода последовательных релаксаций [6] для различных параметров глобальных итераций. Следует отметить, что увеличение Ng соответствует увеличению числа обращений к программе расчета упругих деформаций крыла, т. е. к возрастанию времени счета. Приведенные сравнения указывают на целесообразность выбора Ng = 5, Nf = 20, так как при этом расчет выполняется с наименьшими затратами машинного времени без каких-либо потерь в точности получаемых результатов. Во всех дальнейших приведенных расчетах число глобальных итераций учета влияния упругости конструкции крыла на его обтекание равняется пяти. Как показывают результаты расчета, приведенные на рис. 3, б, скорость сходимости итерационной процедуры предлагаемого метода расчета для комбинации крыло — фюзеляж рассматриваемого типа слабо зависит от коэффициента релаксации ©о- Поэтому в дальнейших приведенных расчетах cog = 0,8.
На рис. 4, а приведено сравнение распределений коэффициента подъемной силы сечений крыла по его размаху суСеч(г) для упругого крыла «С» при значениях скоростного напора q=> 10 000 Н/м2 и 15000 Н/м2 и для жесткого крыла «К» при М = 0,8, с,уа = 0,5=const. Значение скоростного напора q= 10 000 Н/м2 соответствует высоте крейсерского режима полета //=11 000 м. Следует отметить существенное перераспределение циркуляции по размаху при переходе от жесткого крыла к упругому.
-----местное крыло,,К";д=01!/мг ; а = 5,5°;сХ1 = 0,0112
----- ” " „С”
-----ипригое » •• ,д=10000н/мг',<х=5,ЗЯ°\сх;=0,0115
----- ~ •• ■■ ; 15000Н/м2] 6,39°] 0,0117
Рис. 4
М„~0,б
О 0,25 0,50 су„
Рис. 5
Как показывают результаты расчетов, деформация средней линии профилей сечения крыла с учетом влияния упругости составляет не более 0,01% от местной хорды и, таким образом, влияние упругости конструкции крыла на околозвуковое обтекание рассматриваемой комбинации крыло — фюзеляж, в основном, сводится к изменению закона крутки крыла б (рис. 4,6). Следует отмети^ несоответствие закона крутки упругого крыла крейсерской крутке жесткого крыла, что указывает на необходимость проведения точного учета влияния упругости конструкции крыла при разработке и проектировании его аэродинамической компоновки.
Расчетные зависимости суммарных аэродинамических характеристик обтекания упругого крыла тг(сУа) и сУа (а) для разных значений скоростного напора ц при Моо = 0,8 приведены на рис. 5. В частности,
из рис. 5, б видно, что значения коэффициентов подъемной силы крыла «С» при <7 = 0 Н/м2 (жесткое крыло) и q= 15 ООО Н/м2 при а = 5° = const различаются примерно на 30%.
ЛИТЕРАТУРА
1. Евсеев Д. Д. Расчет некоторых аэродинамических характеристик упругого самолета методом коэффициентов влияния. — Ученые записки ЦАГИ, 1978, т. 9, № 6.
2. Chipman R., Waters С., Mackenzie D. Numerical computation of aeroelastically corrected transonic loads. — AIAA Paper, 1979,
N 79—0766.
3. Whitlow W. Jr., Bennet R. M. Application of a transonic potential flow code to the static aeroelastic analysis of three — dimensional wings.—AIAA Paper, 1982, N 82—0689.
4. Jameson A., С a u g h e у D. A. Numerical calculation of the transonic flow past a swept wing. — N. — J. University ERDA Report COO 3077—140, 1977.
5. Владимирова H. А. Исследование обтекания прямых и стреловидных крыльев большого удлинения при околозвуковых скоростях. — Ученые записки ЦАГИ, 1983, т. 14, № 4.
6. Савицкий В. И. Расчет невязкого околозвукового обтекания крыла с фюзеляжем. — Ученые записки ЦАГИ, 1983, т. 14, № 2.
7. Mu г man Ё. М., Cole J. D. Calculation of plane steady transonic flow. — AIAA J., 1971, vol. 9, N 1.
Рукопись поступила З/Х 1985 г.