УЧЕНЫЕ ЗАПИСКИ Ц АГ И Том VIII 1977
№ 4
УДК 533.6.011.5
ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ОБ ОБТЕКАНИИ ОСЕСИММЕТРИЧНОГО ВОЗДУХОЗАБОРНИКА СВЕРХЗВУКОВЫМ ПОТОКОМ ИДЕАЛЬНОГО ГАЗА НА РАЗЛИЧНЫХ УГЛАХ АТАКИ
В. Е. Макаров
Приведены результаты расчета обтекания равномерным сверхзвуковым потоком идеального (невязкого, нетеплопроводного) газа изолированного осесимметричного воздухозаборника смешанного сжатия при нулевом и ненулевом углах атаки. Стационарная система уравнений газовой динамики, записанная в виде интегральных законов сохранения, решалась с использованием конечноразностного метода сквозного счета первого порядка точности. В рассмотренных примерах получены значения коэффициентов расхода и дополнительного сопротивления воздухозаборника на нулевом угле атаки и углах атаки, отличных от нулевого.
Несмотря на прогресс в развитии численных методов расчета течений идеального и вязкого газа, решение большинства газодинамических проблем, возникающих при создании высокоэффективных входных устройств (ВРД), осуществляется, в основном, экспериментальным путем, что связано с большими затратами средств и времени. В силу исключительной сложности течения в воздухозаборнике используемые расчетные методы носят чаще всего, эмпирический и полуэмпирический характер. Торможение потока в сверхзвуковых диффузорах смешанного сжатия осуществляется в системе косых скачков уплотнения, создаваемых центральным телом, и в замыкающем прямом скачке, который расположен в канале воздухозаборника вблизи сечения минимальной площади. На расчетном режиме работы все косые скачки сходятся на передней кромке обечайки. Используя известные методы расчета сверхзвуковых течений идеального газа, можно найти параметры во всей возмущенной области вплоть до замыкающего скачка уплотнения, определить расход и дополнительное сопротивление воздухозаборника, возникающее из-за отклонения режима его работы от расчетного. В тех случаях, когда замыкающий скачок находится вне канала вблизи сечения входа и перепуск воздуха
через кромку обечайки невелик, результаты, полученные на основе расчета сверхзвукового течения,также могут быть использованы для оценки характеристик входного устройства.
Для исследования осесимметричных многоскачковых диффузоров целесообразно применять конечноразностные методы расчета. В этом случае удается проанализировать пространственное течение, возникающее при обтекании воздухозаборника сверхзвуковым потоком на углах атаки, отличных от нулевого. Несмотря на то, что этот режим работы входного устройства является типичным, в настоящее время нет методики, которая позволяла бы достаточно быстро и точно определять расходные характеристики и дополнительное сопротивление осесимметричного воздухозаборника в диапазоне углов атаки, встречающихся на практике. Более того, поскольку в эксперименте трудно определить величину дополнительного сопротивления, возникающего из-за нерасчетности течения во входном устройстве, в литературе отсутствуют конкретные данные о величине и поведении этой важной характеристики воздухозаборника на различных углах атаки.
В данной работе исследовалось обтекание однородным сверхзвуковым потоком идеального газа осесимметричного диффузора со смешанным сжатием. Система уравнений Эйлера интегрировалась методом сквозного счета сверхзвуковых течений [1, 2], обладающего, применительно к рассматриваемой задаче, рядом преимуществ. Получившие большое распространение различные варианты метода характеристик и неявного конечноразностного метода [3] весьма эффективны при расчете сверхзвукового обтекания гладких тел, однако их трудно использовать при изучении течений, которые содержат особенности (ударные волны, тангенциальные разрывы, изломы стенок и т. п.). Выделение этих особенностей ведет к значительному усложнению вычислительного алгоритма и часто требует предварительной информации об их положении. Методы сквозного счета, основанные на дивергентной форме записи уравнений газовой динамики, позволяющие не заботиться о разрывах, широко используется для решения сложных задач внутренней и внешней аэродинамики. Так например, в работе [4] метод, предложенный в работах [1, 2], был успешно применен для анализа сверхзвуковых конических течений. Возможности этого метода при расчете сложных пространственных течений были показаны в работе [5] на примере исследования обтекания конфигурации типа летательного аппарата. В работе [6] метод с успехом использовался при изучении течения в плоских и кольцевых каналах при наличии в потоке косых скачков уплотнения и центрированных волн разрежения.
1. Расчет стационарного течения идеального газа с постоянными теплоемкостями проводился в цилиндрической системе координат х, г, ср, ось Ох которой совпадала с осью симметрии исследуемой конфигурации. Предполагалось, что всюду в потоке выполнено условие и^>а, где и — проекция вектора скорости в рассматриваемой точке на ось Ох, а — местная скорость звука. Исходная система уравнений, система конечноразностных соотношений и метод ее решения, который является стационарным аналогом известного метода С. К. Годунова [7], подробно описаны в работах [1,2]. В процессе расчета по известному распределению газодинамических параметров в некоторой заданной области плоскости х =х0 при определенных граничных условиях находятся границы возмущенной области и распределение параметров в сечении х =х0 1г, где
к — шаг в направлении Ох, выбираемый из условия устойчивости системы конечноразностных уравнений. Границами рассчитываемой области по координате г могут быть поверхности обтекаемых тел, поверхность скачка уплотнения и свободная поверхность (граница струи). Достоинством применяемого метода является то, что он позволяет без усложнения вычислительной процедуры рассматривать все типы граничных условий.
Программа, с помощью которой проводились расчеты в данной работе, была составлена на основе алгоритма, использованного в работе[4] для анализа конических течений. Основная и вспомогательные программы написаны на языке АЛГОЛ-бО применительно к ЭВЦМ БЭСМ-6.
2 Центральное тело исследуемого воздухозаборника, схематически изображенного на фиг. 1, представляет собой тело вращения, образующая которого состоит из нескольких прямолинейных
участков. Соседние участки с разными углами наклона к оси симметрии либо пересекаются с образованием угловой точки (точки В, С и В', С' на фиг. 1), либо плавно сопрягаются посредством дуги окружности заданного радиуса (участки О, Е и Ь', Е' на фиг. 1). Начальный участок центрального тела, носок которого помещен в начало координат, представляет собой круговой конус ■с полууглом раствора бк=15°. Участки ВС (В'С') и СО(С'Ь') ■образуют с осью х углы 20 и 25° соответственно. Точки В, С и В’, С’ служат местами возникновения косых скачков уплотнения. Образующие внутренней и внешней поверхностей обечайки имеют на ее передней кромке углы с осью симметрии в 14 и 20° и состоят из ряда плавно сопрягающихся прямолинейных участков. За характерный размер задачи принимался диаметр ИР' входа в канал воздухозаборника. Все линейные размеры отнесены к этой величине. Расстояние д:вх от носка центрального тела до сечения входа в канал менялось от 1,09 до 1,4 с целью изменения режима работы входного устройства при фиксированном числе М«, в набегающем потоке и заданном угле атаки а. Здесь и в дальнейшем индекс оо указывает на то, что соответствующие параметры берутся в невозмущенном потоке. Все полученные ниже результаты соответствуют значению показателя адиабаты х=1,4.
И'
Фиг. 1
3. Течение на начальном участке центрального тела является коническим и находилось в соответствии с работой [4] путем установления решения по координате л;. В качестве исходных параметров брались параметры в невозмущенном потоке. Расчетная область в сечении х = const была заключена по координате г между телом и ударной волной и по координате ср — между наветренной (ср = 0) и подветренной (<р = 180°) полуплоскостями симметрии. Степень установления решения определялась по максимальной разности между газодинамическими параметрами на поверхности тела при одном и том же значении ср в соседних по л: расчетных сечениях. Результаты расчетов сравнивались с данными работы [3] при одних и тех же значениях Мао, 6К и а, а также между собой при различном числе разбиений расчетной области. При угле атаки не равном нулю и сетке, содержащей 20 слов по координате г и 32 сектора по координате <р, максимальное отличие параметров на поверхности тела от взятых из работы [3] не превышало 2%, причем ошибка по давлению и продольному компоненту скорости была меньше 1%. Время установления решения в этом случае составляло 25—30 мин. При нулевом угле атаки при том же числе расчетных слоев между телом и ударной волной в направлении ср бралась всего одна ячейка. Указанная выше точность определения параметров на поверхности конуса в этом случае достигалась установлением решения в течение 1,5 мин. По известному в каждой меридиональной плоскости f= const значению угла головного скачка с осью д: возмущенная область из конечного сечения установления пересчитывалась в сечение х= 1,09, которое являлось исходным для дальнейшего расчета обтекания центрального тела. Распределение газодинамических параметров при этом, в силу коничности течения, остается неизменным. В процессе счета, который велся до значения х = хвх, определялась внешняя граница возмущенной области, являвшаяся продолжением головного скачка уплотнения. Внутренние скачки, исходящие от мест излома образующей центрального тела, не выделялись. Время, потребное для расчета этой области, не превосходит 5 мин при сетке, содержащей 20X32 = 640 точек. Для продолжения расчета в канале или на внешней поверхности обечайки воздухозаборника параметры в сечении входа посредством линейной интерполяции пересчитывались на сетку, связанную с границами соответствующей области. В первом случае исходной областью является кольцевой зазор между передней кромкой обечайки и центральным телом. Во втором случае кромка обечайки определяет внутреннюю границу исходной расчетной области. Внешней границей служила линия пересечения головного скачка уплотнения с плоскостью х = хвх. При таком способе задания начальной области для расчета обтекания мотогондолы необходимо, чтобы скачок целиком охватывал обечайку и проходил на достаточно большом расстоянии от ее передней кромки. Если граница возмущенной области попадает в канал или находится вблизи кромки, в качестве внешней границы исходной области используется окружность, охватывающая область возмущенного течения в плоскости х = хвх. Расчет на этом этапе для такой же разностной сетки, как и на предыдущем этапе, от х = хвх до значения х = \,65 занимал не больше 3 мин. Таким образом, расчет пространственного обтекания центрального тела и внешней поверхности обечайки, а также сверхзвукового течения в канале воздухозаборника занимает примерно 40 мин на ЭВЦМ БЭСМ-6.
Фиг. 2
Соответствующее время для двумерного случая (а = 0) равно 2 мин.
В качестве примера расчета на фиг. 2 показано распределение статического давления, отнесенного к статическому давлению в невозмущенном потоке, вдоль образующих обтекаемых поверхностей. Число М в набегающем потоке было равно 2,35, сечение входа в канал располагалось на расстоянии хвх=1,34 от начала координат. Кривые 1, 2 и 3 относятся соответственно к центральному телу, внутренней и внешней поверхностям обечайки. Сплошными линиями дано распределение давления при нулевом угле атаки, штриховыми — в наветренной (|? = 0), штрихпунктирными — в подветренной (<р=180°) полуплоскостях симметрии при угле атаки, равном 10°. Характер изменения давления свидетельствует о наличии в рассчитываемой области многочисленных внутренних скачков уплотнения и зон ускорения потока. Это подтверждается представленной на фиг. 3 картиной течения в меридиональной плоскости, совпадающей с плоскостью симметрии задачи. Сплошными линиями указаны линии постоянства числа М. Их сгущение определяет области резкого изменения параметров в потоке (скачки уплотнения, волны разрежения).
4. Рассмотрим начинающуюся на бесконечности поверхность тока, которая приходит на переднюю кромку обечайки. При установившемся течении и фиксированном режиме работы входного устройства в канал попадает только тот газ, который находится в объеме, ограниченном этой поверхностью тока. Выделим плоскостью л; = 0 часть объема полученной таким образом трубки тока. Внешняя 2, и внутренняя 2а поверхности, ограничивающие этот объем, являются поверхностями тока, причем, совпадает с поверхностью центрального тела при 0 х-< .хвх. Поверхности £оо и Евх, имеющие площади Fco и являются сечениями рассматриваемой трубки тока плоскостями х = 0 и х — хвх. На фиг. 1 площадь, ограниченная контурами НРйА и Н'Е'СА представляет собой сечение выделенной части трубки тока меридиональной плоскостью,
а линии HF и Н' F', AG и AG', НН\ FF' есть следы поверхностей Е,, Е2, Seo и Евх в этой плоскости.
Обозначим через и, v, w проекции вектора скорости на оси декартовой системы координат, связанной с исходной цилиндрической таким образом, что ось Ох и начало у них являются общими, а координатная плоскость ху лежит в плоскости симметрии задачи.
Пусть V = (и2 -4- V2 + т2)2 —модуль вектора скорости, р и р— давление и плотность газа. Введем вектор Р с компонентами X, У, Z следующим равенством:
где п■—вектор внешней нормали к поверхности Заметим, что из-за наличия в задаче плоскости симметрии Z = 0. Проекции X и
У вектора Р на оси Ох и Оу можно определить, применяя закон сохранения импульса к выделенному ранее объему трубки тока, заключенному между плоскостями х — 0 и х = л:вх:
X = § р(1у [ {р -г ри2) йу йг — § (р + ри2) йу йг + / р йу йг,
Интегралы по Евх вычисляются непосредственно по известным параметрам в сечении л: = хвх. Последние слагаемые в выражениях
Y — J р dx dz = I puv dy dz — J puv dy dz — f p dx dz.
S1 SBX Sco S2
для X и У представляют собой волновое сопротивление Л'ц. т и подъемную силу Уц.т части центрального тела прил;-<л:вх, которые определяются в процессе расчета его обтекания. Интегралы по Ноо могут быть вычислены по параметрам в невозмущенном потоке. С учетом сказанного выражения для X и У можно переписать в другом виде:
X — | (р + Р«2) йу йг — (роэ +- роэ и*) /г00 + Хт,
2ВЛ
У = J рUV dy dz — peo Moo ®oo /Гсо + *Vr,
где
Ноо= Voo cos a, Z>co = l^ooSina.
Для определения Foa запишем условие постоянства расхода в трубке тока:
j рtidy dz = \ ри dy dz.
SBX *00
Значение интеграла в левой части и, тем самым,величина расхода G могут быть определены непосредственно. Второй интеграл определяется через параметры невозмущенного потока:
J Щ dy dz = poo l^oo cos a^oo.
^co '
Обозначая через Fo6 площадь сечения, определяемую контуром передней кромки обечайки, и вводя коэффициент расхода воздухозаборника / по формуле
/=_______-_____
Рос ^0= Роб ’
непосредственно служащей для его определения, можно выразить отношение Fсо IFоб через / и а:
Fco = / .
Fob cos а
Коэффициенты дополнительного сопротивления и дополнительной подъемной силы (или коэффициенты сопротивления и подъемной силы поверхности тока 2^ определяются следующим образом:
- О _____У____
' Poo ^ Foo ’ У fo3Vl Fo6 ■
При угле атаки, отличном от нуля, дополнительные сопротивление и подъемная сила в системе координат x'y'z', ось Ох' которой направлена по вектору скорости невозмущенного потока,а ось Oz’ совпадает с осью Oz, определяются коэффициентами с'х и с , которые равны
с'х = сх cos a + су sin я, c'y = — cx sin a + cy cos a.
На фиг. 4 показано изменение величин сх и / для исследуемого воздухозаборника при нулевом угле атаки, различных числах в невозмущенном потоке и различных положениях сечения входа относительно носка центрального тела. Особенность рассматриваемой конфигурации состоит в том, что ни на одном из рассчитанных режимов головной скачок не приходит на переднюю кромку обечайки и, поэтому, при всех значениях М«, и л:вх величина сх
больше нуля. Таким образом, несмотря на то, что предлагаемая методика предназначена прежде всего для анализа пространственных течений, простота и возможность расчета входных устройств сложной формы на нерасчетных режимах работы делает ее привлекательной и для осесимметричных течений.
Более интересны в практическом отношении кривые, иллюстрирующие поведение относительного коэффициента расхода /=///0
и относительного коэффициента сопротивления Сх=='с'х/сх,, где / и с* — коэффициенты при угле атаки а, /о и сх при а = 0, при раз-
0,8
' X \ \ \ \ \
\ \ \ N ^-1,15- ~-1,19’ Л28-Л’Зй' ^1,32, — 134 — N43
г
X і,гг«Ґ
9 ; №
Ч С, с ■СГ /,» * — 116 1,09
2,0
2.5
ом
'1,30 ,т ¿¡г
,1,26
\26
\ 1,28 _ 1/
4 1,30
т
ос
Фиг. 4
Фиг. 5
личных углах атаки и различных положениях сечения входа. Для М = 2,35 такие кривые, полученные при помощи развитого в работе подхода, приведены на фиг. 5.
В заключение автор благодарит А. Н. Крайко, М. Я- Иванова и Т. В. Никитину за полезные советы и внимание к работе.
ЛИТЕРАТУРА
1. Иванов М. Я , Крайко А. Н., Михайлов Н. В. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. „Журн. вычислит, матем. и матем. физ.“. 1972, т. 12, № 2.
2. И в а н о в М. Я., К р а й к о А. Н., Михайлов Н. В. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. „Журн. вычислит, матем. и матем. физ“. 1972, т. 12, № 3.
3. Бабенко К.И., Воскресенский Г. П., Любимов А. Н., Русанов В. В. Пространственное обтекание гладких тел идеальным газом. М., „Наука“, 1964.
4. И в а н о в М. Я., Крайко А. Н. К расчету сверхзвукового обтекания конических тел. „Журн. вычислит, матем. и матем. физ.“ 1973, т. 13, № 6.
5. Иванов М. Я., Никитина Т. В. К расчету пространственного обтекания сверхзвуковым потоком тел сложной формы. „Ученые записки ЦАГИ“, т. 4, № 4, 1973.
6. Виноградов В. А., Захаров Н. Н., Иванов М. Я. Расчет торможения двумерного сверхзвукового потока невязкого газа в канале и возможность реализации такого течения. „Ученые записки ЦАГИ“, т. 6, № 2, 1975.
7. Г о д у н о в С. К. Разностный метод численного расчета разрывных решений уравнений газовой динамики. „Матем. сб.“, т. 47 (89), вып. 3, 1959.
Рукопись поступила 9ІІХ ¡976 гг..