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

________УЧЕНЫЕ ЗАПИСКИ и А Г И

Т о м XI 19 80

М 5

УДК 533.66J.013

РАСЧЕТ СВЕРХЗВУКОВОГО ОБТЕКАНИЯ ИНТЕРФЕРИРУЮЩИХ КРЫЛЬЕВ

Д. Н. Минайлос

Метод сквозного счета [1] использован для расчета обтекания системы крыльев сверхзвуковым потоком. Решены две задачи. В первой из них два тонких треугольных крыла расположены друг за другом по схеме .утка*; представлено сравнение с результатами расчетов, использующих линеаризованные уравнения газовой динамики. Во второй задаче под крылом носителя симметрично расположены два тонких треугольных крыла меньшего удлинения.

1. При использовании сверхзвуковой интерференции в связи с большой стоимостью экспериментальных исследований и недостаточной точностью приближенных теоретических подходов все большее значение приобретают численные методы решения полных уравнений Эйлера, требующие применения мощных ЭВМ. Такими методами к настоящему времени решено несколько задач обтекания сложных тел с интерференцией. Исследованы, например, течение в двугранном угле (наиболее полное исследование выполнялось в [2]), взаимодействие течения около кругового конуса с ударной волной [3], течение около самолетных компоновок различных типов. Последней задаче посвящены несколько работ, но все они не содержат анализа интерференции частей самолета, а ограничиваются получением ноля течения. Отметим первую в этой группе работу [4] и сборник [5], содержащий несколько других работ.

Численный метод [1] обладает рядом достоинств, позволяющих использовать его при решении задач интерференции: простотой задания и изменения форм тел, широким диапазоном возможного задания входных данных задачи, относительно небольшим потребным временем счета и возможностями довольно просто изменить размеры областей счета и счетных ячеек. Расчет проводится в декартовой системе координат х, у, г (и, V, а1)- Порядок аппроксимации используемой схемы заключен между первым и вторым и на сетке в поперечной плоскости у, г, содержащей — 3200 узлов.

дает ошибки в значениях локальных функции около 2—3%. Точность определения интегральных характеристик должна быть значительно выше; однако в окрестности вершины крыла, когда на его размахе расположено только несколько узлов сетки, аппроксимация течения отсутствует, что приводит к ошибкам в определении локальных аэродинамических характеристик, достигающим 10—20%. Локальными аэродинамическими характеристиками здесь названы характеристики части крыла, расположенной до некоторого сечения л: = const, отнесенные к площади этой части крыла.

лг*

Рис. 1

В диапазоне изменения углов атаки, углов стреловидности и чисел Мос, рассматриваемых в настоящей работе, область с низкой точностью локальных характеристик Ьх имеет размеры — 10А (Л = ку = Иг — размер ячеек счетной сетки в направлении осей у и г, рис. 1). При счете в направлении оси х с увеличением числа счетных узлов на размахе крыла точность возрастает и решение устанавливается, суммарная относительная погрешность уменьшается до 1—2%. Но в окрестности носовой части крыла всегда есть некоторая область Дх, где точность аппроксимации низка. Это особенно важно в тех случаях, когда одно из интерферирующих крыльев мало относительно основного масштаба течения и область Дл занимает значительную часть малого крыла. В этом случае для получения высокой точности необходимо после расчета всего поля провести дополнительные расчеты части поля в окрестности маломасштабного крыла на мелкой сетке в неравномерном поле потока, параметры которого получены интерполяцией из расчета всего течения.

Как только возмущения течения достигают некоторой части границы счетной области, на ней включается алгоритм „срезки" возмущенного поля, основанный на экстраполяции.

Для повышения точности на начальном участке расчета при решении второй рассмотренной задачи используется алгоритм „сжатия” поля; начальное поле, занимающее небольшую физическую область, рассчитывается на сетке, размеры ячеек которой вдвое меньше основных. Это позволяет вдвое уменьшить зону плохой аппроксимации поля Дл. В качестве размеров счетного поля второй

задачи на начальном участке взяты величины у1( г,, показанные на рис. 1, б. Как только область возмущений численного решения достигает границы счетного ноля (координата х*, рис. 1, б), происходит автоматическое укрупнение счетных ячеек или в размерах счетного поля — уменьшение вдвое линейных размеров крыла и возмущенной области в плоскости у, г — „сжатие1* поля. Для такого пересчета интерполяция не используется, просто из поля выбрасываются все узлы, имеющие по направлениям у или г нечетную индексацию:

/(1В№+ 1 - /, ^=/(ГВАг +1-2/, 2Д

///</</ЯЛ /(/ВА2 + /, У) = /ивлг-т 21, 2J),

1ВА/.+ 1</</К,

/= \Р, Р. «. V, ■гг»);

здесь р и — давление и плотность, IBAZ, 1Н и /А — индексы слоя ячеек, расположенных непосредственно под крылом, нижнего и верхнего слоев возмущенной области, / и J— индексы слоев ячеек в направлении осей у и г.

Внешние узлы поля, не входящие в новом масштабе в область возмущения, заполняются параметрами набегающего потока. Так как в методе центры ячеек расположены на расстояниях Ау(/ —0,5) и Аг(У— 0,5) от оси х, то при пересчете новое поле в плоскости у, г образуется со сдвигом в сторону начала координат на размер А/4.

В результате работы такого алгоритма полосы разрывов в иоле решения, размытые при счете на 2—4 интервала счетной сетки, сжимаются по ширине вдвое относительно размеров сетки (до 1 — 2 интервалов). Такое сужение зон разрыва не соответствует порядку аппроксимации схемы, но ввиду ее монотонности дальнейший счет проводится без каких-либо затруднений, и уже через несколько шагов разрывы размываются до обычных размеров на заданной сетке. Численные эксперименты показывают, что метод может быть использован для расчета разрывов большой интенсивности, описанных с любой степенью размытости, допускаемой дискретным заданием начального поля.

Кроме алгоритма „сжатия* поля, применяется счет на более мелкой сетке, позволяющий рассчитать только начальную часть поля интерференции. Во второй задаче такие расчеты проводятся с применением аппарата „сжатия“ ноля.

2. В первой задаче два тонких плоских треугольных крыла расположены друг за другом но схеме „утка“ и имеют общую плоскость симметрии. Для сравнения с результатами, основанными на линейной теории, геометрические характеристики крыльев взяты из (6]. Оба крыла имеют одинаковый угол стреловидности у=653. Как и в [6], варьировались относительные размеры крыльев и высота Н расположения второго крыла относительно первого. Длина центральной хорды первого крыла с1 принята равной 1 (рис. 1, я), число Мо,. = 2, угол атаки первого крыла а, = 4°, угол атаки второго крыла относительно набегающего потока а, = 0;2; 4; 8°. Шаг счетной сетки А = 0,02. В расчетной плоскости у, г располагается 72Х X 42 счетных узла. Задняя кромка первого крыла и вершина второго лежат в одной плоскости х=1. После того как проведен расчет обтекания первого крыла и получено поле течения в сечении 1, производится поворот системы координат х, у, г на

угол Да = а,— а,, чтобы при расчете обтекания второго крыла оно располагалось в плоскости у = const. При этом значения компонентов скорости и и v пересчитываются, а сдвигом точек поля плоскости х = 1 в направлении оси х пренебрегается. Оценки показывают, что погрешности в локальных параметрах поля при таком пересчете не превосходят 0,6%. Первое крыло располагается в плоскости У=Ух =0,6 (или yt =0,8).

Обтекание одного первого крыла исследовано и сопоставлено с результатами [6J. Течение коническое, крыло обтекается в режиме Z? — 1 — I [7]. Ударная волна отсоединена от кромок (в линейной теории кромки дозвуковые), линия растекания потока расположена на кромке крыла.

Сопоставим численное решение с точным решением линейной теории для конического течения. Максимальное отличие в локальных значениях компонентов возмущенной скорости в декартовой системе координат на нижней поверхности крыла составляет 60%

для "~м" и 50% для - , а на верхней поверхности — соответст-

° со а со

венно 160% и 265%. При этом значения wB имеют в некоторой области разные знаки. Сравнение аэродинамических характеристик можно провести, рассматривая их в поле параметров а|3 и £>. Г? = (М^ — 1, >. = 4ctgz — удлинение крыла]; определяющим является параметр jft. Рассматриваемый вариант задачи соответствует области параметров, где отличие результатов, полученных разными методами (численным и по линейной теории), должно быть небольшим. Действительно, оно составляет ~ 1,0%. Этот результат позволяет надеяться на хорошее соответствие и в случае оценок интерференции.

За первым крылом расположены веер течения разрежения, контактный разрыв и хвостовой скачок уплотнения, которые в окрестности плоскости симметрии близки к двумерным, и течение в этой области может моделироваться взаимодействием двух плоских потоков [8]. Течение около второго крыла, расположенного в возмущенном потоке, уже не является коническим. Если второе крыло расположено выше первого но оси у(Н<^0), оно интерферирует с хвостовым скачком первого крыла. Если крылья лежат в одной плоскости (// = 0), то носовая часть второго крыла лежит в плоскости луча веера течения разрежения, и параметры потока, набегающего на это крыло, меняются слабо. В случае же //>• 0 второе крыло с некоторого значения х=^х0 пересекает веер разрежения, а хвостовой скачок не взаимодействует со вторым крылом.

Это отличие в обтекании приводит к качественному отличию в распределении давления по второму крылу при различных значениях Н (рис. 2). В случае /У>0 давление над крылом при *>-*<> резко падает. При Я = 0 изменения давления при росте координаты х невелики. При /У<0 хвостовой скачок падает на нижнюю поверхность второго крыла в окрестности д: = дс0> резко повышая давление в диапазоне 1,3<jc<1,6. Эти особенности сказываются на интегральных характеристиках крыла (рис. 3 и 4). На рисунках представлены локальные значения коэффициентов подъемной силы Cyi и момента тг, второго крыла в зависимости от координаты х. Индекс i означает, что коэффициенты взяты с учетом интерференции. В окрестности х= \ расположена область Ьх, поте-аппроксимации. В качестве характерной площади взята площадь

Ю

второго крыла до данного сечения х, в качестве характерной длины — расстояние от носка этого крыла до данного сечения -с. Черными треугольниками показаны результаты расчетов, использующих густую сетку с шагом А, = А/1,6.

Штриховыми линиями показаны характеристики при значении Н = 0,22, сплошными — при /7 = 0 и штрихпунктирными — при //=—0,22. Светлыми треугольниками отмечены характеристики первого крыла, кружками — характеристики второго крыла в невозмущенном набегающем потоке. К этим последним значениям должны стремиться величины Су/ и тг1 при неограниченном росте ко-

распространяется далеко вниз по потоку. С ростом угла х, отличие результатов от предельных возрастает. Вдоль оси х на рис. 3 и 4 можно выделить три участка; на каждом из них поведение зависимостей С,, и тг1 (при Нф 0) значительно отличается от поведения на соседних участках.

В области / (1,1 <х < 1,3) характеристики имеют ярко выраженное „плато**. Затем в области 11 влияние веера разрежения или ударной волны (при //= + 0,22 и —0,22 соответственно) приводит к резкому увеличению модулей производных функций. Наконец, в области 111 при д:>1,7 наклоны кривых значительно уменьшаются, и результаты начинают асимптотически стремиться к предельным значениям. Характерно, что зависимости, соответствующие значениям Н = — 0,22 и 0,22, близки по своему поведению и численным значениям друг к другу до величин а„~8°, хотя нелинейности в этих случаях вызваны разными причинами. При //<0 поведение кривых определяется падением скачка уплотнения от задней кромки первого крыла на нижнюю поверхность второго крыла, а при //>0 — взаимодействием веера разрежения от задней кромки первого крыла с верхней поверхностью второго. Близость этих характеристик свидетельствует о том, что в поле за первым крылом интеграл изменения давления соответствует соотношениям линейной теории (число М.х невелико). Прирост значений Су , в области нелинейностн в обоих случаях составляет около 0,065—0,07 при всех значениях угла а;.

В случае # = 0 поведение коэффициентов Су1 и тг1 иное: они слабо зависят от координаты х. Это объясняется тем, что в этом случае крыло не пересекает резких неоднородностей в иоле течения.

Полученные результаты позволяют построить зависимости Су1 (а„) и тг1(а2) при заданных значениях х = const. Сопоставим эти характеристики с характеристиками первого крыла. Функции Су1(а,) и /кг,(а,) в диапазоне изменения угла атаки а, <8° практически линейны. Характеристики второго крыла в этом диапазоне углов атаки имеют некоторую нелинейность. При этом знак второй производной зависимости Су1(а,) совпадает со знаком величины Н С уменьшением значения х нелинейность характеристик возрастает.

Рассмотренные зависимости Су1(я2) позволяют также оценить для второго крыла заданной длины осредненный по всей поверхности крыла скос потока. Его величина определяется значением угла а0 = а, (Су| = 0). Естественно, что значения a0<ai. они различны для разных величин Н и убывают с ростом координаты х. В пределе при х -* оо значения а0 должны стремиться к нулю.

Следуя [6], введем для оценок линейной теории коэффициенты интерференции. Пусть и К, —значения подъемной силы неинтерферирующих крыльев; )\2 — часть подъемной силы второго крыла, определяемая интерференцией; 5, и S,- площади крыльев; q — скоростной напор; К,= К, + К,,. Тогда производные коэффициента подъемной силы по углам атаки определяются формулами:

к, >з „ к12

а‘— <7$,*, * а2— аП — qsi а, •

Локальные аэродинамические коэффициенты системы, отнесенные к площади 5,, можно приближенно записать так:

Значения х с индексами — координаты точек приложения соответствующих подъемных сил. Момент рассчитан относительно точки ;с = 0, величиной Н по сравнению со значениями х пренебрегаем.

В рассматриваемом случае а1 = а2 = а = 43; Су , = Су2 = Су=0,143; .<1 = 2/3, х.., = 1 + 2 3 (х — 1); тг , = ^^ тг, = — 0,0955.

В численных расчетах определены локальные значения Су1 и тг1 с учетом интерференции (см. рис. 3 и 4). Они позволяют найти коэффициенты интерференции:

йц = — [Су1 — Су1 ;

л*,, = 1 + | —|тг1 + Су (х2 — 1)||.

Коэффициенты отнесены к площади 5-.., момент *тг1 рассчитан относительно вершины второго крыла х=\. Сравнение этих величин с полученными в [6) (рис. 5) показывает, что линейная теория качественно, а при //> 0 и количественно, правильно отражает поведение значе-ний а,2. Отличия от результатов численных расчетов не превосходят 10% (для вари- V антов со значениями // = 0,22 и 0,23; при совпадении значений Н отличие было бы еще меньше) и 23% (Н — 0). В окрестности точки 52/5, ■= 1 результаты линейной теории имеют в значениях а,,(// = 0) значительно более резкий экстремум. Перегибы кривой .г,, (при Н — 0,22) в численном решении значительно более резкие. Отличия в значениях л12 достигают 40% (/7 = 0) и 71%

(//=*0,22 и 0,23). Таким образом, момент-ные характеристики определяются на основе линейной теории с погрешностями примерно в четыре раза бйльшими, чем в зна- 1,5

чениях коэффициента подъемной силы.

Пунктирной кривой показано положение центра давления второго крыла хг в невозмущенном потоке.

3. Во второй задаче под тонким плос- и

ким треугольным крылом с углом стреловидности Х:= 70° симметрично расположены

два треугольных крыла с углами у2 = 75°. Схемы и расположени крыльев даны на рис. 1,6. Число" Мао —2,5, углы атаки обои крыльев а = 5° или 10°.

В расчетной области по у, г расположено 54X57 ячеек. Значение Ха — 0,6; Хв= 1: у а = 0,1; гд = 0,25; _у1 = 0,27, у* = 0,54. На начальном участке счет проводится на мелкой сетке с шагом Л = 0,005; вершина крыла-носителя (крыла /) расположена в точке (0; 0,16; 0). „Сжатие” поля осуществляется в области, определяемой значением х* = 0,224 (в варианте с углом а = 5°) и соответственно х* = 0,193 (при г = 10°). На сетке с шагом Л = 0,01 вершина крыла / лежит в точке (0; 0,32; 0), крыла подвесного груза (крыла 2) — в точке (0,6; 0,22; 0,25). В сечении х= 1, где лежит задняя кромка крыла 2, на его размахе размещается 22 счетных ячейки, а на полуразмахе крыла /—37 ячеек.

В условиях равномерного набегающего потока режимы обтекания крыльев определяются с помощью диаграммы 6 — а (7). Для крыла 2 это режим £> — 1, для крыла / — £> — 1 или С—1. Линии растекания располагаются на кромках, вихревые жгуты над кромками отсутствуют. Над верхней поверхностью крыла / осуществляется течение без отрыва (а = 5°) или с отрывом (а = 10°). Во втором случае (режим И) отрыв образуется под висячей ударной волной [7]. Над крылом 2 течение без отрывов (режим I).

Вершина крыла 2 расположена в невозмущенном набегающем потоке, ниже крыло пересекает основная волна, и оно оказывается в возмущенном от крыла / поле потока. Вдоль размаха крыла 2 локальные характеристики этого поля при а =10° изменяются при увеличении координаты х от 0,6 до 1 в следующих пределах: давление 1,34 <р!(рх < 1,48; плотность 1,33<р/рсо< 1,44; число М 2,3<М<2,18; угол скоса потока в вертикальной плоскости (местный угол атаки крыла 2) 6,4° > агс^ ~ > 2,0°; угол скоса потока в горизонтальной плоскости (местный угол скольжения крыла 2) 4,1° < аг^ ■— < 6,Г. Таким образом, число М около крыла 2 значительно меньше, чем в набегающем потоке; но это не приводит к смене режима обтекания кромок крыла с первого режима на второй (т. е. к образованию на кромках отрывного течения), так как существенно уменьшаются локальные углы атаки.

Характерным для обтекания крыла 2 является значительное скольжение. Конические течения со скольжением исследованы и классифицированы в [7]. Эту классификацию можно использовать в анализе обтекания крыла 2, учитывая, что оно не является строго коническим. Анализ поля проводится в конических переменных, имеющих центр в вершине крыла 2. Режим течения Б— 1—1 — — 0—1 —I. Скольжение приводит к асимметричной форме головной волны и возникновению момента М,.

На рис. 6 показано поле изобар р рх в сечении л: = 1,0. Цифрами в кружочках обозначены: / — основная волна и 2 — висячий скачок у крыла /, 3—основная волна и 4, 5—области разрежения у крыла 2.

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

Рис. 6

крыльев. Эти сгущения и определяют форму и положение сечений крыльев. Над крылом / изобары не выводились на графопостроитель с целью экономии времени построения чертежа. Положение висячего скачка 2 определено из других чертежей. Изобары построены с постоянным шагом ^p|pv = 0,05. Ломаные линии снизу и справа ограничивают область возмущенного течения в сеточном иоле. Как показывают расчеты, волны 1 и 3 в результате взаимодействия сливаются в одну волну. По крайней мере до сечения а'=1,4 разрешающая способность сетки не позволила ' выявить разделение волн после пересечения.

Наличие крыла 2 влияет на распределение давления по нижней поверхности крыла /: образуется расширяющаяся область разрежения, давление в которой постепенно выравнивается. Распределение давления по верхней поверхности крыла 1 слабо зависит от наличия второго крыла. За крылом 2 возникает течение, включающее в себя веер разрежения, контактный разрыв и хвостовой скачок. Это течение первоначально близко к плоскому, но хвостовой скачок имеет по размаху разную интенсивность: в окрестности внутренней кромки крыла 2, т. е. ближайшей к плоскости симметрии течения, перепад давления в хвостовом скачке 1р!р:с = 0,15, а в окрестности внешней кромки Ьр\р?э = 0,03 (в сечении х=1,05). Это вызвано тем, что над внешней кромкой крыла хвостовой скачок проходит область больших значений давлений

Рис. 7

и плотностей, чем над внутренней. Область разрежения, образованная над внутренней кромкой крыла 2, быстро „всплывает* при х> 1 к нижней поверхности крыла 7, образуя на ней область пониженного давления. Эта область затем смещается к плоскости симметрии и размывается. Область разрежения над внешней кромкой не „всплывает”, так как над ней расположена зона повышенного давления в окрестности линии растекания на крыле 1. Картина течения в следе за крылом 2 является асимметричной относительно плоскости г = 0,25.

Локальные аэродинамические коэффициенты (рис. 7) Су1, тг, для крыла 1 и Су 2, т,2 для крыла 2 отнесены к площадям соответствующих крыльев. В качестве характерной длины взято локальное значение х. Значения тг) рассчитаны относительно оси х = 0, у — 0,32; тг, — относительно оси х = 0,6, у = 0,22; значения тхг — относительно оси .у = 0,22, г = 0,25. Кружками показаны значения для крыльев, обтекаемых равномерным набегающим потоком. Область потери аппроксимации накладывается на зону входа крыла 2 в возмущенную область, что снижает точность определения характеристик в этой зоне.

1. М и н а й л о с А. Н. Расчет сверхзвукового обтекания крыльев с учетом схолящих с кромок тангенциальных разрывов в рамках модели, использующей систему уравнений Эйлера. „Изв. АН СССР, МЖГ*. 1978, № 1.

2. Anderson D. A., Nang! a R. К. Comparison of numerical and experimental .conical* flow fields In supersonic corners wilh compression and/or expansion. Aeron. Quart., 1977, november.

3. К u 11 e r P., S a k e 11 L. Three-dimensional, shoclk-on shock interaction problem. A1AA. Paper N 75—49, 1975.

4. D'Attorre L., Bilyk M. A., Serglant R. J. Three-dimensional supersonic flow field analysis of the B-l airplane by a finite difference technique and comparision with experimental data. AIAA Paper N 74— 189, 1974.

5. Aerodynamic analysis requiring advanced computers. NASA SP-347, Langly Research Center, march 4 -6, 1975.

6. Purshouse М., Nangia R. K. Application of linearised supersonic wing theory to the calculation of some aircraft interference flows. In .Computational methods and problems in aeronautical fluid dynamics*, 1976.

7. M и ii a ii л о с A. H. Расчет обтекания крыльев сверхзвуковым потоком газа. Труды 6-й международной конференции по численным методам в гидродинамике (Тбилиси, 1978), изд-во АН СССР, 1979.

8. Минай л ос А. Н. Расчет поля сверхзвукового течения с вихрями за тонким прямоугольным крылом. .Ученые записки ЦАГИ*, т. 9, Д6 5. 1978.

Рукопись поступила 23j IV 1979

«Ученые записки» .V* 5