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

УЧЕНЫЕ ЗАПИСКИ Ц АГ И Т о м XV 198 4

№ 2

УДК 533.6.011

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

С. М. Босняков, В. В. Коваленко, А. Н. Минайлос

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

1. Результаты исследования плоских сверхзвуковых воздухозаборников различной конфигурации широко используются при проектировании перспективных летательных аппаратов. Эти результаты в настоящее время получаются как экспериментальным путем, так и численно [1—6]. Расчетные исследования развивались путем последовательного усложнения формы конфигураций, такой путь обеспечивал преемственность и физическую наглядность результатов. Большое значение для подтверждения достоверности расчетов имеют серии специальных экспериментов [2 — 6], проводимых на моделях, полностью соответствующих телам, исследуемым расчетным путем. Получено хорошее соответствие между данными расчетов и экспериментов с учетом в ряде случаев влияния толщины вытеснения пограничного слоя.

Метод сквозного счета использует декартову прямоугольную систему координат, связанную с клином (рис. 1,а) так, что на нижней, наветренной стороне клина расчетная сетка имеет дробные ячейки. Расчет в этих ячейках проводится нестандартным способом. Газодинамические разрывы в численном решении не выделены, а получаются в виде областей больших градиентов параметров потока. При этом за разрывами параметры определяются с точностью первого порядка аппроксимации. При таком подходе расчетная сетка не зависит от решения, и при малых значениях х отсутствует его аппроксимация (между размытой ударной волной и поверхностью тела оказывается недостаточное количество счет-

Рис. 1

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

Более точные результаты в области малых значений д: можно получить с помощью метода установления [7], однако они ограничены только областью коничности течения от угловых точек клина и неприемлемы для длинных клиньев, представляющих практический интерес. Получение более точных численных результатов требует перехода к разностным схемам высокого порядка аппроксимации и представления ударных волн в поле течения в виде разрывов.

При переходе к схемам второго порядка аппроксимации по сравнению со схемами первого порядка, использующими „распады разрывов" на боковых гранях счетных ячеек (схемы [8] и [9]), возможна, кроме повышения точности, также значительная (в 6—8 раз) экономия машинного времени в счетном узле. Это объясняется тем, что итерации при расчете распадов разрывов требуют больших временных затрат. Так как для оптимизации форм аэродинамических объектов необходимо наличие большого числа результатов численных расчетов, вопрос об экономии машинного времени чрезвычайно важен и актуален в настоящее время.

2. Ниже используется разработанная авторами программа расчета обтеканйя сложных тел и компоновок сверхзвуковым потоком невязкого газа. Программа использует маршевую схему второго порядка аппроксимации типа [10 — 11] с выделением внешней ударной волны.

Выбор рабочей системы координат (^, Ь, 5) определяется особенностями геометрической формы поперечных сечений исследу-

емого тела. Поперечные сечения аппарата (сечения х — const связанной с телом декартовой системы координаты х, у, z) преобразуются в плоскости x — q = const в сечении, по форме близкие к кругу, с помощью ряда последовательных преобразований Кармана— Треффтца [13]:

?/+i

* ( hjj

i+i \b + hn

j—U 2, . ... J, (1)

здесь \j—-комплексная переменная в плоскости поперечного сечения (при у = 1 tl = z-\-iy, = f(x) = zp -г iyp — комплексная коор-

дината угловой точки поперечного сечения), h* — точка, сопряженная точке h, = <fj(x) — внешний угол угловой точки, J—число угловых точек поперечного сечения и число преобразований.

Соотношение (1) преобразует плоскость ^ в E/+i так, что образ точки hjj переходит в точку S/+I = l, и угловая точка исчезает. Последовательное преобразование контура поперечного сечения от ^ до 5/+1 позволяет получить гладкий контур, близкий к окружности. Все особые точки преобразования (1) размещаются внутри контура тела и якобиан преобразования g = d(x, у, z)/d(q, t, s) не имеет особенностей в расчетной области (вне контура). Однозначность преобразования обеспечивается проведением разрезов, соединяющих пары сопряженных особых точек h и h* таким образом, чтобы разрезы проходили внутри контура тела.

Затем в плоскости S/+i осуществляется переход к полярным координатам г, ф (с центром системы в точке %j+1 = 0 и осью ф = 0, направленной горизонтально) и, наконец, к переменным

т rs(q, ф)-/•„(?, Ф) v ’

где гв и rs — радиусы сечения тела и головной волны в плоскости

£/+1-

Проведя замену q = х, получим систему координат q, t, s. При этом равномерное распределение счетных узлов вдоль поверхности тела в сечении q = const соответствует в плоскости х — const сгущению узлов в окрестности острых выпуклых углов сечения тела. Преобразование (х, у, z) -*(q, t, s) является неособым в исследуемой области и переводит эту область в параллелепипед.

В конкретном рассматриваемом примере поперечное сечение тела в физической плоскости представляет собой прямоугольник, у которого горизонтальные стороны постоянны, а вертикальные увеличиваются пропорционально значению координаты л:. Для трансформации счетного поля используется два преобразования Кармана — Треффтца. Особые точки преобразований А и В (рис. 1, б) располагаются вблизи угловых точек, и преобразования выполняются в таком порядке, чтобы в последней плоскости %2 после преобразований ось ф = 0 проходила через точку В. Это позволяет сохранить на нижней поверхности, представляющей наибольший интерес, постоянное число (а именно половину) счетных узлов, расположенных на поверхности тела. Остальные узлы на теле располагаются на сторонах СА и АВ, при этом с ростом значения ко-

ординаты х сторона АВ увеличивается, и расчетные узлы сдвигаются в сторону точки В и могут переходить со стороны СА на сторону АВ. В плоскости после обезразмеривания расстояния между узлами сетки одинаковы, и расчет проводится на равномерной сетке.

Уравнения Эйлера в системе (<7, Ь, 5) записываются в дивергентном виде:

Здесь

Е =

+ Р /+^5 = 0;

р ■■= Р (1 — И2 — 'V* — ни2).

(3)

рия р и1 р и?

рии* + крцх р — рии* -\-kptj. 0 = щи? -)- крьх

рг»и» + крцу > ' — уии* + кр1у 1 V р юи? + к,р8у

ртаи4 + крдг ртюи* р пои? + крз2

2%

х — показатель адиабаты,

и,4 = Щх-\-ЮЦу + Ч0Ц„

и'■ = иЬх 4- vty + ио1г, и5 = тх + ту + те>зг,

(4)

Р£,

Р = Рё, £

д(х, у, г) д(д, Ь, 5)

Входящие в (4) производные определяются дифференцированием системы (2).

Граничными условиями являются:

1) на поверхности тела (* = 0):

I *=0 = ■

0

(5)

— условие непротекания;

2) условие симметрии течения (вектор скорости набегающего потока лежит в плоскости симметрии):

а> = и2 = *>г=Л = р, = 0

при = 1)

(6)

3) на головной ударной волне (2 =

— условия Рэнкина—Гюгонио.

4) на острой боковой кромке тела при перетекании потока через кромку

— условие Жуковского—Кутта образования тангенциального разрыва на острой кромке.

Задача ставится для ^-гиперболической системы уравнений и применяется маршевый подход, поэтому как и в методе сквозного счета основным требованием применимости методики является условие соблюдения во всем поле течения сверхзвуковой величины компонента скорости и4. В качестве начальных данных в некоторой

плоскости g — q0 задается известное поле течения. Обычно это или рассчитанное ранее поле, или поле, заполненное параметрами невозмущенного набегающего потока.

В конкретном случае исследуемого тела на некотором начальном отрезке оси х имеются области конического и плоского течения, поэтому возможно использование процесса „установления" конического течения по координате х от условий равномерного набегающего потока без предварительного расчета области в окрестности передней кромки тела. Выбор положения начального сечения q0 — x0 определяется путем численных экспериментов так, чтобы коническое поле устанавливалось на некотором малом начальном участке оси х и изменение величины х0 в два раза не влияло на результаты расчетов.

Опишем кратко применяемый численный метод.

Во внутренних точках расчетной области используются схемы, учитывающие направление вектора скорости и ориентацию конуса влияния возмущений в исследуемой точке. Обозначим точку, в которой ищется решение, индексами я+1, J (индексы соответствуют координатам q, t, s). Очевидно, что в зависимости от ориентации конуса влияния в точке л + 1, /, j решение в этой точке зависит от параметров в различных точках слоя п. Проиллюстрируем выбор шаблона для координаты t при s = const. Рассмотрим четыре типа явных двухшаговых схем второго порядка точности:

а)

Е, = Е?

£?+1 =

El + 4г (2 - tf-i + ^-2)

м

(7)

б)

Е?+1 =

в)

г)

%= ^--w(Fl+i-F?)-,

-i- [El+Et - (2 (Fi+l - ?,) - FUi + E1+ 2)

Et = E?----^ (F"+i — Fl)\

ег\= 4- [e7+% - ;

Et — El----^J~(F"—Ei—i)',

Г7Я+1 ____

ti = -7Г

Эти схемы применяются в трех различных вариантах программы. В первом из них во всем счетном поле используется одна однородная схема „в“ или одна схема „г“. Это метод Мак-Кормака [10]. Во втором варианте по каждому направлению ^ и « используются обе схемы „в“ и „г“ выбор между ними определяется наклоном ком-

понентов скорости и* и и? по отношению к линиям t = const и ^ = const соответственно. При 0 используется схема „в“, при и*^>0 — схема „г“ (аналогично для направления s). Эта схема „флюгер" на основе схемы Мак-Кормака.

Наконец, схемы „а“ и „6“ совместно со схемами „в“ и „г“ применяются в третьем варианте алгоритма. Схемы „а“ и „б“ используются в случаях, когда конус влияния проходит ниже (в случае „а"}, или выше (в случае ,,б“) линии t = const. Для выбора расчетной схемы в узле рассчитываются наклоны характеристического конуса в точке в плоскостях t = const и s = const и по их значениям заполняются две матрицы условных чисел, определяющих выбор схемы. Прием оказался возможным благодаря факторизации схемы счета по координатным направлениям f и s. Односторонние схемы „а“ и „6“ ■близки к схеме, предложенной в работе [11].

В рассмотренных ниже примерах применяется третий вариант алгоритма.

Полученный вектор Еп+1 обрабатывается с помощью моното-низатора, близкого к предложенному в работе [12].

Обычно на поверхности тела (£ = 0) используется для счета схема „б“, на головной ударной волне (£=1)—схема *а“, в плоскости симметрии течения —схемы „в“ и „г“. При этом величины вектора G в точках, расположенных за плоскостью симметрии течения, вычисляются с учетом условия (6).

Для корректировки значений параметров на поверхностях тела и ударной волны используются процедуры, разработанные в [14, 15], алгоритм образования на кромке тангенциального разрыва создан в соответствии с рекомендациями работы [1].

3. В качестве исследуемого тела взят клин с углом раствора 14°, установленный в сверхзвуковом потоке так, что угол атаки нижней его поверхности а = 15°. За характерный размер, равный единице, взята половина ширины клина. Расчеты проведены при числе Моо = 2,5 До значений л; = 6, в методе второго порядка точности получены результаты как без применения, так и с применением алгоритма принудительного образования тангенциального разрыва [1] на острой кромке поперечного сечения тела (вблизи особой точки В).

На рис. 2 для двух расчетов, проведенных на разных сетках с выделением головной волны, сопоставляются линии постоянных

Рис. 2

значений скосов arctg (vju) и arctg (wju). Штриховые линии соответствуют сетке с числом узлов 22X16, а сплошные линии—сетке с числом узлов 42 X 32, при этом в первом случае на нижней поверхности располагалось И узлов, а во втором—21 узел. Изолинии проводились с шагом 2°. Поля изолиний параметров очень чувствительны к изменению размеров счетных узлов, и полученные отличия положения изолиний для двух вариантов в целом, следует считать малыми, а сходимость по сетке в поле течения хорошей. Исключение составляет окрестность угловой точки тела,

Рис. 3

особенно у боковой поверхности клина. Эти результаты получены без применения алгоритма образования отрыва на кромке и свидетельствуют о том, что число узлов у боковой поверхности клина,, несмотря на их сгущение в окрестности угловой точки, недостаточно для детального описания решения (известно, что за угловой точкой течение имеет сложную структуру с областью перерасши-рения и „висячей" ударной волной).

На рис. 3 сравнивается распределение давления на нижней поверхности тела при значении х = 0,5, полученное тремя различными численными методами. При х = 0,5 еще не нарушена конич-ность течения, и, следовательно, применима программа, использующая методику [7]. Данные на рисунке соответствуют области влияния угловой точки. Координата г определяется выражением:

* = (** — *)/ (* 1 ). где гк—координата кромки (при принятом обезразмеривании гк—\). Сплошной линией изображены результаты, полученные методом сквозного счета, штриховой—методом второго порядка и штрих-пунктирной—методом первого порядка с выделением головной волны [7]. В окрестности кромки (при г <0,5) все три метода дают близкие результаты; в области г»0,7 результаты, полученные с помощью методов первого порядка, отличаются друг от друга на 5%. Причины такого отличия авторам не ясны.

Рис. 4 иллюстрирует поведение давления в возмущенном поле течения на вертикальных линиях 2 = 0,95 и 1,05, проходящих в окрестности боковой кромки. Для сравнения с программой, использующей методику [7], результаты представлены при малом значении х — 1,6, в котором течение еще остается коническим. Вертикальная координата у на графике соответствует принятой в методе [2]. При этом верхняя поверхность клина располагается в горизонтальной плоскости у = сопб1 = 5; нижняя поверхность пересекает плоскость я = 1,6 при 3; ~ 4,61.

Результаты расчетов, полученные с помощью метода [2] на грубой и мелкой сетке, показаны на графике штриховыми и сплошными линиями соответственно. Головная волна в этом методе не выделена, она размыта и представлена в виде области больших градиентов в нижней части графика. Горизонтальные линии, проведенные посредине области больших градиентов, указывают приближенно положение ударной волны после перехода к представлению ее в виде разрыва. Определенное таким образом положение волны не совпадает для двух сеток, а на густой сетке—совпадает

х=1,В

»У

2 =0,95

777777777777

ТТ

450-

4,25-

4,0 -

О

О

о

3,15

с результатами расчета методом второго порядка на сетках 22X16 и 42X32 узла.

Результаты, полученные методом второго порядка на сетке 22X16 узлов показаны квадратами черными — без учета отрыва потока в окрестности точки В и белыми — с учетом отрыва. Черные квадраты со штрихом —результаты без учета отрыва на сетке 42X32 узла. При 2 = 0,95 показаны только черные квадраты, так как все три результата ложатся практически на одну кривую, это справедливо и при г =1,05 в области значений .у <4,61. У боковой поверхности клина, как уже отмечалось при анализе полей течения (см. рис. 2), сходимость результатов при уменьшении размеров ячеек на применяемых сетках отсутствует. Это объясняется сильным разрежением за угловой точкой и образованием затем висячей ударной волны. Эта ударная волна в методе второго порядка не выделялась и представлена в виде области резкого увеличения давления в диапазоне 4,65 <.у<4,85. Мельчение узлов сетки приводит к росту разрежения в окрестности угловой точки и к увеличению интенсивности ударной волны. В случае введения области

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

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

Белыми кружками изображены результаты расчетов по методу [7]. При приближении снизу к угловой точке они расходятся

с данными, полученными двумя другими методами; особенно велики расхождения при 2=1,05. В области за угловой точкой методы [2] и [7] дают примерно одинаковые результаты.

На рис. 5 показано изменение величины коэффициента прост-1

ранственности ЛГпр = J(ср!ср пл) йг на нижней поверхности клина

о

в зависимости от значений х. Здесь ср — коэффициент давления, српл—коэффициент давления в плоском течении на клине с углом в 15°. Белые кружки—результаты обработки экспериментальных данных [2] по распределению давления. Расчетные результаты работы [2] для сеток, содержащих 10 и 20 узлов на полуразмахе, изображены штриховой и сплошной линиями, данные, полученные по методике [7] — пунктирной линией, результаты метода второго порядка—черными и белыми квадратами и штрихпунктирной линией для сеток с числом узлов 10X8, 22 X 16 и 42X32 соответственно. На двух последних сетках результаты на графике практически совпадают. Они совпадают также с экспериментальными данными в диапазоне х>1. Близки к экспериментальным данным также результаты, полученные при применении метода [7]. Результаты метода второго порядка на сетке 10x8 узлов, когда на нижней поверхности клина расположено только 5 узлов, имеют низкую точность. Метод сквозного счета во всем диапазоне изменения значений х дает завышение значения Кпр по сравнению с результатами эксперимента и метода второго порядка. Сгущение узлов

дает сдвиг в нужную сторону, а применение экстраполяции Ричардсона на нулевой размер шага существенно приближает результаты к экспериментальным. Для схем первого порядка экстраполяция линейна и при отношении шагов hjh2 = 0,5 она изменяет результаты расчетов на густой сетке на значения разности между грубым и более точным результатом в сторону, противоположную результатам, полученным на грубой сетке.

После применения экстраполяции отличия от экспериментальных значений при х<3 не превосходят 2,6%.

Авторы благодарят Н. Н. Захарова, предоставившего им материалы расчета, полученные по программе работы [7].

ЛИТЕРАТУРА

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

МЖГ, 1978, № 1.

2. Босняков С. М., Минайлос А. Н., Р е м е е в Н. X. Обтекание клина конечной ширины сверхзвуковым потоком газа.—Ученые записки ЦАГИ, 1977, т. VIII, № 6.

3. Босняков С. М., Минайлос А. Н., Ремеев Н. X. Исследование пространственного обтекания двухступенчатых клиньев конечной ширины сверхзвуковым потоком газа.— Ученые записки ЦАГИ*, 1980, т. XI, № 1.

4. Босняков С. М., Ремеев Н. X. Исследование пространственного обтекания клина конечной ширины сверхзвуковым потоком газа при наличии углов атаки и скольжения.—Ученые записки ЦАГИ,

1981, т. XII, № 6.

5. Босняков С. М., Ремеев Н. X. Исследование пространственного обтекания плоского воздухозаборника с боковыми щеками сверхзвуковым потоком газа.—Ученые записки ЦАГИ, 1980, т. XI, № 5.

6. Босняков С. М., Ремеев Н. X. Исследование пространственного обтекания плоского воздухозаборника в компоновке с треугольным крылом сверхзвуковым потоком газа.—Ученые записки ЦАГИ, 1983, т. XIV, № 4.

7. Дуганов В. В., Иванов М. Я. Сверхзвуковое обтекание боковой кромки половины клина.—Ученые записки ЦАГИ, 1977, т. VIII, № 6.

8. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики.—Матем. сб., 1959, т. 3,

№ 47.

9. Колган В. П. Применение принципа минимальных значений производной к построению конечно-разностных решений газовой динамики,—Ученые записки ЦАГИ, 1972, т. III, № 6.

10. М с С о г ш a k R. W. The effect of viscosity in hypervelocity impact cratering.—AIAA Paper 1969, N 64—354.

11. Warming R. F.,'Beam R. M. Upwing second-order difference schemes and applications In aerodynamic flows.—AIAA J., 1976, vol. 14,

N 9.

12. Ж м а к и н А. И., Фурсенко А. А. Об одном классе монотонных разностных схем сквозного счета.—Препринт ФТИ АН СССР им. А. Ф. Иоффе, 1979, № 623.

13. Moretti G. Conformal mapping for computations of steady, threedimensional, supersonic flows.—Numerical Laboratory Computer Met- ■ hods in Fluid Mechanics. ACME, 1976.

14. Warming R. F., Kutler P., L o m a x H. Second and third-order noncentered difference schemes for nonlinear hyperbolic equations.—

AIAA J„ 1973, vol. 11, N 2.

15. T h о m a s P. D.. Vinokur М.. В a s t i a n о n R., С о n t i R. J. Numerical solution for the three-dimensional inviscid supersonic flow of a blunt body.—AIAA J., 1974, vol. 10, № 7.

Рукопись поступила 26jVII 1982 г.