Метод расчета обтекания летательного аппарата дозвуковым потоком идеального газа Текст научной статьи по специальности «Механика и машиностроение»

_______УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Том XVI 1985

№ 4

УДК 629.7.015.3

МЕТОД РАСЧЕТА ОБТЕКАНИЯ ЛЕТАТЕЛЬНОГО АППАРАТА ДОЗВУКОВЫМ ПОТОКОМ ИДЕАЛЬНОГО ГАЗА

А,, Н. Любимов, Ю. С. Сорокин

Излагается метод расчета обтекания установившимся потенциальным дозвуковым потоком идеального газа летательного аппарата, состоящего из фюзеляжа и несущих элементов. Метод включает методику представления формы ЛА в ЭВМ и методику определения параметров потока газа на поверхности ЛА. Первая методика основана на замене элементов ЛА совокупностью отсеков, поверхность которых определяется аналитической зависимостью. Вторая методика основана на представлении потенциала возмущенных скоростей в виде суммы потенциалов от источников, распределенных на поверхности ЛА, и диполей, распределенных по срединным поверхностям несущих элементов и пелене. Интенсивность источников и диполей определяется после решения системы линейных алгебраических уравнений, к которой сводится интегральное уравнение для производной по нормали к поверхности ЛА от потенциала скоростей. Приводятся примеры расчетов и сравнение с экспериментальными данными.

1. Линеаризация нелинейного уравнения для потенциала скоростей основана на предположении о малости возмущений, возникающих в пространстве около тела при обтекании его дозвуковым потоком газа.

Известное линеаризованное уравнение для потенциала скоростей <р(х, у, г) сводится к уравнению Лапласа с помощью преобразования Прандтля — Глауэрта — Х — х/У\ — ; У = у; 2 = г,

где Моо — число Маха невозмущенного потока.

Отметим, что некоторые аналитические решения задач о дозвуковых потенциальных течениях идеального газа даны в [1], а некоторые вопросы конечно-разностных методов решения уравнения Лапласа изложены в [2]. Однако решение задач обтекания дозвуковым потоком трехмерных тел типа ЛА конечно-разностными методами связано с большими трудностями. Это объясняется сложностью построения корректных конечно-разностных алгоритмов для уравнений эллиптического типа и прежде всего трудностями, связанными с заданием условий на свободных границах и обеспечением построения устойчивых разностных схем. Кроме того, затраты машинного времени при численном решении трехмерных задач итерационными методами чрезмерно велики даже при использовании мощных современных ЭВМ.

Вследствие отмеченных выше причин в дозвуковой аэродинамике получили распространение методы, в которых потенциал скорости <р записывается в виде суммы потенциалов от особенностей (источников, вихрей, диполей), распределенных на некоторых плоских базовых поверхностях, схематически представляющих форму ЛА [3]. При этом граничное условие непротекания выполняется не на самой поверхности ЛА, а сносится на плоские базовые поверхности. Такой метод соответствует принципу линеаризации исходного нелинейного уравнения для потенциала скоростей.

В последние годы разрабатываются методы, в которых особенности распределяются непосредственно на поверхности ЛА или на его срединной поверхности, а граничное условие непротекания выполняется на самой поверхности Л А [4, 5]. Такой подход не согласуется с принципом линеаризации, если при его реализации рассматриваются ЛА с большой относительной толщиной элементов компоновки. Эти методы позволяют определить не только разность коэффициентов давления на элементах ЛА, но и значения скоростей и коэффициентов давления р в точках на поверхности ЛА.

Характерной особенностью методов этого направления является необходимость достаточно точной аппроксимации поверхности ЛА, имеющей в общем случае сложную форму. Для решения этой задачи наибольшее распространение получили так называемые панельные методы, в которых поверхность ЛА заменяется совокупностью плоских или неплоских поверхностей — панелей. Например, в [4—6] поверхность ЛА, заданную совокупностью точек (Хг, у и 2г), аппроксимируют четырехугольными плоскими панелями. Другой подход представления трехмерной поверхности ЛА предложен в [7].

Для решения задачи обтекания ЛА источники обычно распределяются на панелях, аппроксимирующих форму ЛА, а диполи или вихревой слой распределяются на срединных поверхностях и пелене, сходящей с задних кромок несущих элементов ЛА. При этом интенсивность особенностей обычно принимается постоянной на всей панели [4, 6]. В работе [8] показано, что результаты расчетов [9] обтекания газом тел сложной формы с постоянной и переменной интенсивностью особенностей в пределах панели практически совпадают.

Для расчета обтекания ЛА потоком газа, как правило, интегральное уравнение, определяющее потенциал скоростей или его производную по нормали к поверхности, сводится к системе линейных алгебраических уравнений, которая решается численно. В [10, 11] численно решается непосредственно интегральное уравнение.

Разработанные методы аппроксимации поверхности ЛА панелями и методы аэродинамического расчета требуют больших затрат машинного времени и занимают значительные объемы оперативной памяти ЭВМ. Поэтому их реализация в системах автоматизированного проектирования ЛА затруднительна. Для использования в таких системах необходима разработка специальных методик аппроксимации поверхности ЛА, разбиения ее на панели и методик аэродинамического расчета, позволяющих рассчитывать аэродинамические характеристики всего ЛА сложной формы с относительно необльшими затратами машинного времени. В настоящей работе излагаются такие методики, специально разработанные для включения их в систему автоматизированного проектирования ЛА.

2. Рассмотрим в линейной постановке задачу об обтекании ЛА установившимся дозвуковым потенциальным потоком идеального сжимаемого газа. Летательный аппарат состоит из фюзеляжа и нескольких

несущих элементов, расположенных по отношению к фюзеляжу произвольно, но имеющих небольшие углы установки и небольшую относительную толщину. Угол атаки ЛА будем считать малым, а угол скольжения ¡3 = 0.

Потенциал скоростей У*, 2*) в произвольной точке

пространства Прандтля—Глауэрта РУ*, 2Г.*) представим в виде суммы интеграла, взятого от источников, распределенных по поверхности ЛА, и интеграла, взятого от диполей, распределенных по срединным поверхностям несущих элементов ЛА и пелене, и запишем для точки Р* производную по нормали к поверхности от потенциала скоростей, устремив точку Р* к поверхности ЛА

где Г, — поверхность фюзеляжа; М — количество несущих элементов; ¿2 — поверхность т-то несущего элемента; Е3 — поверхность, состоящая из срединной поверхности яг-го несущего элемента и поверхности пелены Е5, сходящей с задней кромки т-го несущего элемента; — внешняя нормаль к поверхности Е3; /? =

= V (X — Х%)2 + (X — ^*)2 + (2 — где X, У, Е — координаты по-

верхности ЛА, являющиеся переменными интегрирования; 5 — интенсивность источников; О — интенсивность диполей.

Граничное условие непротекания поверхности Л А потоком газа запишем в виде

где иоо — скорость набегающего потока, Ых — косинус угла между нормалью к поверхности Л А и осью X.

Таким образом, задача об обтекании ЛА дозвуковым потоком идеального сжимаемого газа сведена к решению интегрального уравнения (1) с граничным условием (2) и условием Жуковского — Чаплыгина, заключающимся в равенстве нулю перепада давления на задних кромках и пелене.

3. Для представления формы ЛА при проведении аэродинамического расчета методом, изложенным в п. 2, использована методика, которая является развитием работы [12].

Летательный аппарат считается состоящим из фюзеляжа и нескольких несущих элементов (рис. 1). Фюзеляж и несущие элементы мо-

(1)

(2)

о

У

10

гут иметь один или несколько отсеков, следовательно, ЛА представляется совокупностью отсеков. С каждым отсеком свяжем систему координат ОХ]У]2з, где / — номер отсека, и введем систему координат Охуг, связанную с ЛА.

Задав поверхность /-го отсека в системе координат, связанной с отсеком, с помощью преобразования

где А, и/?з — соответственно матрица и вектор преобразования, можно определить поверхность отсека в системе координат ЛА. Задав поверхности всех отсеков, определим поверхность всего ЛА.

Поверхность отсека представим линейчатой поверхностью, уравнение которой в параметрическом виде запишем следующим образом:

где Ь и но — параметры, а кривые г1 = г1(Ь) и = г21 (¿) являются

направляющими линейчатой поверхности. Параметр * изменяется от 0 до 2тг, а но изменяется от 0 до 1.

Уравнение линейчатой поверхности, представляющей поверхность /-го отсека фюзеляжа (рис. 2) в декартовой системе координат Ох ¡у ¡г], запишем в виде

где ¿/ф — длина отсека фюзеляжа, ус1 и ус2 — расстояния от начала систем координат Оург® и Оу£ г§ (рис. 2) до оси Охі.

Функции у*(£) и £*(£) определяют форму сечений отсека плоскостью Xj = 0 (>і=1) и плоскостью хі; = ЬФ (п — 2) и имеют вид

где 5П — ширина п-го сечения отсека; Нп — толщина сечения; Рп — кривизна срединной линии сечения; йщ■—коэффициенты, которые позволяют изменять форму верхней и нижней дужки сечения отсека.

(3)

Гу (і, т) = г1{і)(\—тю)+ г21 (*) но,

(4)

х; = и, 0 < * < 2*, 0 ^u = w■Lf^CLf,

У І = ІУ? (О + Ус 1І (1 - + [У2ф (*) + у, 2І(5)

г, = г? (О (1 - «/¿?) + (*) «/1?,

(6)

л=г у/

Рис. 2

Таким образом, поверхность /-го отсека фюзеляжа будем считать заданной уравнениями (5—6), если заданы параметры отсека: Lj, Sn, Нп> Fп у din•

Уравнение линейчатой поверхности, представляющей поверхность ¿-го отсека несущего (управляющего) элемента (рис. 3), в декартовой системе координат OxííjíZí запишем в виде

х[ = х1 (0(1 — u/L4)-j-[x2(t)cos'zl-{-y2(t)smxí] и,Щ + uígxai, У]=У\ (0 0 — u/L*) + [_y2(¿)cosx¿ — х2 (t) sin хг] tifUi + и tg Y0 z¡ = и, 0< г<2тс, 0-<u = w-LK¡ LKi,

(7)

где ¿Г — длина отсека несущего элемента; Хо I — угол стреловидности передней кромки отсека; у01— угол поперечного V отсека; ^ — угол поворота сечения отсека плоскостью = относительно сечения отсека плоскостью г( = 0 (рис. 3).

Рис. 3

Функции xn(t) и yn(t) определяют форму сечений отсека несущего элемента плоскостью z¡ — 0 (я == 1) и плоскостью z¡ = Lf {fi = 2) и имеют вид

х„ (t) = Ьп/2 + bn¡2 {cos (t -f aln sin t) +

+ a2n[aSn + (\ — лвя) sin í] sin (2£ + a4nsin2¿)}, (8)

Уп (t) = bn/2 (cn sin t,+ 2fn sin21),

где bn — хорда профиля в ti-м сечении; сп — максимальная относительная толщина профиля; ]п — максимальная относительная кривизна срединной линии профиля в п-м сечении; ain — коэффициенты, которые позволяют изменять кривизну верхней и нижней дужки профиля в п-м сечении.

Поверхность i-го отсека несущего элемента J1A будем считать заданной уравнениями_ (7)— (8), если заданы параметры отсека:

Li ) Zo b То i' Х(> Ьт Сп> fn> в/я-

Представленная методика, основанная на аналитическом описании поверхности отсеков, позволяет достаточно просто задавать в ЭВМ JIA •сложной конфигурации и разбивать его поверхность на панели. Однако главное достоинство методики состоит в том, что она удобна при проектировании и проведении аэродинамических расчетов, позволяя легко и целенаправленно изменять форму JIA путем изменения традиционно применяемых в аэродинамике геометрических параметров, входящих в уравнения поверхности отсеков (5) — (8).

4. Кратко отметим основные особенности численного решения сформулированной в п. 2 задачи. При построении алгоритма предполагается, что пелена от каждого т-то несущего элемента ЛА образована прямыми вихревыми линиями, сходящими с задней кромки и параллельными оси х. Поверхность фюзеляжа разбивается на А^ф малых четырехугольных панелей. Поверхность каждого т-го несущего элемента разбивается на А^т панелей, а его срединная поверхность — на Ыт/2 панелей. Поверхность пелены разбивается на полосы, количество которых соответствует количеству панелей вдоль размаха т-то несущего элемента.

Считается, что значения £>;, 5 и (ду/д1\) постоянны в пределах каждой панели и равны значениям ^ и {д^!дЩ1 в центре панелей.

Тогда на основе интегрального уравнения (1) для всех центров Р* (Хк, Ук, Z¿) панелей поверхности ЛА, учитывая граничное условие (2) и условие Жуковского — Чаплыгина и заменяя операцию интегрирования суммированием, можно записать систему линейных алгебраических уравнений для определения неизвестных 5, и 0$

и- (Ли =4- + £(Л^ VI,) 5г

м

тп~ 1 (.1 = 1

/=1 лт,2

где

/=1

4-1,2,... (лГф+^ЛГ..),

^‘“-тИЫтг)^,

{Ук1)т = (- 11 V* (-£-) Й22/)т .

(9)

(У*)т

/ т

(10)

(П)

(12)

(У*/)/п = 0 для всех панелей срединной поверхности т-го несущего элемента, кроме панелей, имеющих общую сторону с задней кромкой;

(13)

для панелей срединнои поверхности т-го несущего элемента, имеющего общую сторону с задней кромкой;

25з- — слой (полоса) пелены, ограниченный двумя линиями тока, начинающимися от панелей срединной поверхности т-го несущего элемента, имеющих общую сторону с задней кромкой.

Для удовлетворения условия Жуковского — Чаплыгина в уравнения (9) значения интенсивности диполей на задней кромке т-то несущего элемента £>з. к. ш приняты равными значениям т в центре панелей срединной поверхности, имеющих общую сторону с задней кромкой.

м

Таким образом, имеется система из ( Д^ф 4- Е-'Чп) уравнений

Используем условие симметрии источников [13] на поверхности т-го несущего элемента, т. е. ¿>гт = 5(г+Л^т/2)т, для г от 1 до Л^/2. Тогда систему (9) можно записать в виде

Отметим, что слой диполей не заканчивается на поверхности фюзеляжа, а продолжается внутрь фюзеляжа, причем интенсивность диполей в подфюзеляжной части постоянна вдоль размаха и равна значению интенсивности диполей в бортовом сечении несущего элемента.

Решая систему линейных алгебраических уравнений (14), найдем значения интенсивности источников 5 и диполей В в центрах панелей поверхности ЛА, далее определим значения скоростей и коэффициентов давления р в этих точках. Интегрируя, легко определить суммарные аэродинамические характеристики ЛА от сил давления.

Зная параметры потока на поверхности ЛА — скорости и давления, можно определить плотность, температуру, вязкость и местные значения числа Рейнольдса, рассчитанные по характерному размеру, равному расстоянию от рассматриваемой точки поверхности ЛА до носика фюзеляжа, если точка принадлежит поверхности фюзеляжа, или до носика профиля сечений т-го несущего элемента, если точка лежит на поверхности несущего элемента. Используем принцип локального подобия и будем считать, что параметры потока в точке на поверхности такие же, как если бы эта точка принадлежала пластине, с длиной, равной характерному размеру, и обтекаемой со скоростью, равной местной скорости в рассматриваемой точке поверхности ЛА. Тогда по известным соотношениям для плоских пластин можно определить местные значения коэффициентов трения, а далее и суммарные аэродинамические характеристики от сил трения, в том числе и коэффициент лобового сопротивления при нулевой подъемной силе — сх .

Следует заметить, что интегральное уравнение (1) записано в общем виде для ЛА сложной конфигурации. Однако его можно применять и для решения задач обтекания потоком газа отдельных элементов ЛА. В этом случае в уравнении (1) остаются только интегралы, взятые по соответствующим поверхностям того или иного элемента ЛА, характеристики которого рассчитываются.

(

"Ф

- £/оо (ЛОО* = 4“ + Е (Л^ VI) 5г + 1=1

Была исследована сходимость результатов расчетов по описанной выше методике по количеству панелей на отдельных элементах ЛА — фюзеляже и крыле. Результаты некоторых из этих исследований приведены в [14].

5. Изложенный метод применен для определения обтекания несжимаемым и сжимаемым газом ЛА различной формы с несколькими несущими элементами. Ниже приводятся некоторые результаты расчетов двух компоновок фюзеляж — крыло, для которых известны экспериментальные данные, позволившие провести сравнение расчетных и экспериментальных результатов.

На рис. 4 и 5 представлены коэффициенты давления для компоновки фюзеляж — крыло при Моо —0,113 и а = 3,77°. Фюзеляж является эллипсоидом вращения с большой осью Ь и малой осью 0,2 Ь. Крыло имеет размах Ь и прямоугольную форму в плане, удлинение крыла 1 = 5, хорда Ь = 0,2 Ь, максимальная относительная толщина профиля

с = 0,10, профиль NACA 65010. Крыло установлено на фюзеляже по схеме среднеплан. Передняя кромка крыла находится на расстоянии 0,3 L от носка фюзеляжа. Фюзеляж задается в ЭВМ десятью отсеками, а крыло — одним отсеком. Общее количество панелей 208, из них 88 панелей представляет фюзеляж и 120 панелей — поверхность крыла. Время расчета одного варианта на ЭВМ ЕС1033 составило около 40 мин.

На рис. 4 дано распределение коэффициента давления р вдоль фюзеляжа без крыла (пунктирная линия) и с крылом (сплошная линия) в меридиональном сечении фюзеляжа плоскостью 0 = 45°, которая показана на схеме. На рис. 5 приведено распределение коэффициента давления на верхней и нижней поверхности крыла на расстоянии от фюзеляжа, равном 0,1 L. Сплошной линией показаны значения р компоновки фюзеляж — крыло, а пунктирной линией — изолированного крыла. Крестиками на рис. 4 и 5 нанесены экспериментальные данные работы [15] для компоновки фюзеляж — крыло.

На рис. 6 представлены результаты расчетов при Моо = 0,6 коэффициента подъемной силы су (а), коэффициента продольного момента тг (су) и коэффициента лобового сопротивления сх (Моо) при нулевой подъемной силе компоновки фюзеляж'—крыло [16]. Фюзеляж имеет носовую часть оживальной формы удлинения Хи.ч=2, переходящую в цилиндр диаметром d= 1 и удлинения Л.ц= 10,3. Крыло имеет удлинение Я = 6,54 и сужение г| = 1,63, угол стреловидности передней кромки %о=44° и установлено на расстоянии 5,4 диаметра от вершины фюзеляжа. Профиль крыла чечевицеобразной формы с максимальной относительной

Рис. 4

Рис. 5

толщиной профиля с = 0,10. Общее количество панелей—180, из них 80 представляют поверхность фюзеляжа и 100 — поверхность крыла. Фюзеляж для расчетов задан пятью отсеками, а крыло — одним. Время расчета одного варианта (одно число Маха) составляет около 35 мин.

На рис. 6 сплошной линией показаны расчеты, а крестиками— экспериментальные данные [16].

Приведенные примеры расчетов демонстрируют удовлетворительное согласование экспериментальных и расчетных данных.

ЛИТЕРАТУРА

1. Тихонов А. Н., Самарский А. А. Уравнения математической физики.— М.: Наука, 1977.

2. Роуч П. Вычислительная гидродинамика. — М.: Мир, 1980.

3. Белоцерковский С. М., Скрипач Б. К. Аэродинамические производные летательного аппарата и крыла при дозвуковых скоростях. —

М.: Наука, 1975.

4. Н е s s J. L. Review of integral-equation techniques for solving potential flows problems with emphasis on the surface source method. — Computer Methods in Applied Mechanics and Engineering. Ill, 1975, vol. 5, N 2.

5. Woodward F. A. An improved method for the aerodynamic analysis of wing-body-tail configurations in subsonic and supersonic flow.

P. 1. Theory and Applications, NASA—CR—2228, 1973.

6. M о r i n о L., Chen L. T. and S u с i u E. O. Steady and oscillatory subsonic and supersonic aerodynamics around complex configurations. —

AIAA Journal, 1975, vol. 13, N 3.

7. Давыдов Ю. В. Уравнение одной поверхности зависимых сечений. —В кн.: Кибернетика графики и прикладная геометрия поверхностей. МАИ, 1976, вып. 349.

8. М a s k е w В. Prediction of subsonic aerodynamic characteristics — a case for low-order panel methods. — AIAA Paper, N 81—0252, 1981.

9. J о h n s о n F. Т., R u b b e r t P. E. Advanced panel-type influence methods applied to subsonic flows. — AIAA Paper, N 75—50, 1975.

10. Вернигора В. H., Ир а кл ионов В. С., Павлове ц Г. А. Расчет потенциальных течений около крыльев и несущих конфигураций крыло — фюзеляж. — Труды ЦАГИ, 1976, вып. 1803.

11. Маслов JI. А., Тимербулатов А. М. Расчет давлений на поверхности произвольной комбинации фюзеляжа с несущим крылом при малых скоростях. — Труды ЦАГИ, 1979, вып. 2005.

12. Лебедь В. Г., Любимов А. Н., Русанов В. В. Метод представления и визуализации формы в диалоговых системах машинного проектирования. — Автометрия, 1982, № 4.

13. Р е t г i е J. А. Н. A surface source and vorticity panel method.— Aeronaut. Quart., vol. 39, N 4, 1978.

14. Любимов A. H., С орокин Ю. С. Методика расчета обтекания нетонких крыльев дозвуковым потоком сжимаемого идеального газа.— В сб.: Экспериментальное и теоретическое исследование аэродинамических характеристик ЛА и его частей.— МАИ, 1983.

15. Shin ju Suzuki, Kyiuchiro Washizu. Calculation of wing-body pressures in incompressible flow using Green’s function method. —

J. Aircraft, 1980, vol. 17, N 5.

16. Потапова Л. А., Штейнберг P. И. Волновое сопротивление крыльев с прямой и обратной стреловидностью при околозвуковых скоростях.— Ученые записки ЦАГИ, 1980, т. XI, № 3.

М~=0,6

о расчет

х экспериментов]

0,15 0,5 су 0,6 0,1 0,8 0,9 1,0

\Ч,М;т\ =Ш/

Рис. 6

Рукопись поступила 17/VI1 1984 г.