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

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

Том XIII

198 2

№ 2

УДК 533.6.011.55.011.6

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

А. Г. Зарубин

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

1. Метод С. К. Годунова для решения систем квазилинейных дифференциальных уравнений гиперболического типа [2] основан на использовании точных решений уравнений с кусочно-ностоянными начальными условиями. В случае стационарных сверхзвуковых течений для этого используется решение задачи о взаимодействии равномерных сверхзвуковых нотоков (ЗВРСП).

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

В работе [3] для реального газа предложен метод решения ЗВРСП, в котором для ускорения вычислений используется система таблиц и аппроксимаций термодинамических функций [4]. Однако время расчетов в этом случае заметно увеличивается по сравнению с расчетами течений совершенного газа с теми же граничными условиями по методу [1]

Дальнейшая оптимизация алгоритма вычислений может быть достигнута, если в итерационных нроцедурах численного метода использовать такие термодинамические функции газа, которые слабо меняются в условиях термодинамического равновесия. Так, вместо энтальпии к (р, р) более удобна функция 7 (р, р), определяемая из соотношения

и/ Ч 7 (Рг р) Р ...

7 (Р> Р) — 1 ? ’

где р — давление газа, р — плотность.

При этом удается построить алгоритм решения ЗВРСП, совпадающий с [1] в случае течений совершенного газа и требующий нри совнадающих граничных

условиях в среднем на 15—25% времени больше, чем в случае расчета течений

реального газа.

Для описания термодинамических свойств газа будем, использовать еще одну функцию р), входящую в уравнение для скорости звука

(2)

В термодинамических уравнениях, описывающих реальный газ, функции 1 (Р> Р) и Х(Р>Р) выполняют роль, аналогичную показателю адиабаты совершенного газа. Для вычисления %(р, р) и Л(р,р) используется система таблиц и аппроксимаций [4].

С учетом соотношения (1) адиабату Гюгонио можно представить в виде

7 (Ра. Рг) + I _ , Л Р2 = Рі 7 (Ръ> Рг) ~ 1 Рг 7(Л>Рі)+ 1 .

ї(л.р1>-іЛ+Л

(3)

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

2, При рассмотрении плоской ЗВРСП как в совершенном, так и в реальном газе будем предполагать, что решение относится к одному из четырех основных типов, схематически показанных на рис. I. Здесь двойные линии обозна-

а)

<Г)

чают ударные волны, штриховые — тангенциальные разрывы, пунктирные — границы веера характеристик.

Указанные йинии разбивают плоскость течения на ряд областей, пронумерованных на рис. I цифрами от 1 до 6. В областях 1—4 параметры течения постоянны. Область 1 занимает нижний невозмущенный поток, область 2—верхний. В областях 3, 4 давление и угод наклона скорости газа однн и те же. Область 5 занимает веер характеристик, если происходит разрежение нижнего потока. Область 6 имеет аналогичное назначение для верхнего потока. Параметры течения в областях будем обозначать величинами с индексом, равным номеру области.

Итерационный метод решения ЗВРСП несколько видоизменяется в зависимости от того, какая из показанных на рис. 1 возможностей реализуется, одиако

во всех случаях начальные значения давления р и тангенса угла наклона скорости Cг=tgв в областях 3 и 4 вычисляются, как и в [I], в рамках первого приближения по формулам

Р — iPi + Ps)№ + (£1 — Сг)/2/га, |

С = С! - т (р - Pi) » С2 + тп (р — р2), /

где

—4-^ • (5>

2 1 Pi 9? 1 Рз?2 2 i

pi — давление, рi — плотность, qt—модуль скорости» С* — тангенс угла наклона скорости, М* — число М для нижнего (/=1) и верхнего (/=2) невозмущенных потоков.

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

Ь+2 =(р» р»+2> = х (А> Pi) = */ = const> (6)

Тi+2 = 1 (р> Р»Ч2) = 7 (яь Pi) = = const. (7)

Способ вычисления плотности газа в областях 3 и 4 зависит от соотношения между давлением газа, пайденным по (4), и давлением газа в невозмущенном потоке. Если (случай волны разрежения), плотность газа в г-{~2-й

области определяется из (2) с учетом (6):

1

Pi+2 = Pi(PlPi) *

При p^>pi (случай ударной волиы) для вычисления плотности газа нужно использовать выражение (3), которое при условиях (7) принимает вид

_ (Tt + 1 )р + (7g — 0 Pi р,‘+2 р/ (7 г + 1 )pi + (и —1)р

Модуль q и компопенты вектора скорости и, v, а также число М в областях 3 и 4 вычисляются по формулам

Ч,+2 ~ У 2 («Г- r-^ї ~) ■ М/+2=ї<+2 V■ й*+2 ~ ” ^*+ 2’ г" ~ 1 ’2’ у 1 + са

(8)

где #0 — энтальпия торможения, а при вычислении параметров нулевого приближения используются (б), (7).

Если относительные отклонения р и С от аналогичных величин с индексом 1 и 2 великн, итерационный процесс следует продолжить. В этом случае уточненные значения р и С можно представить в виде

Щ С} — т\ Р\ Ч~ Рг-

р = т{ + т2 ’

С = «1 — т1 (Р “ Р\) = "9 ^2 + т1 (р ~ />»)•

Коэффициенты л* в (9) вычисляются следующим образом: в случае волн разрежения

. и2

(9)

I Pi Ь Рі+2 Чі+2

щ — 1, і = I, 2;

в случае ударных волн

щ ~ п" Лі = 1 ± Сі tg of > * = l’2>

*

где верхний знак соответствует / = I, нижний —г = 2, 5г —угол отклонения /-го потока при прохождении через ударную волну, gi определяется из соотношения

Из законов сохранения массы и импульса на скачке можно получить следующую форму для gi:

где отношение р//р/_)_2 определяется нз адиабаты (3) следующим образом:

Для совершенного газа формула (10) переходит в выражение, полученное в [1].

Для вычисления плотности газа в областях 3 и 4 в случае волн разрежения проинтегрируем приближенно уравнение (2), заменив %(р, р) постоянной величиной х/:

В случае ударных волн плотность газа вычисляется кз (11).

Другие параметры газа в областях 3, 4 определяются по формулам (8), в которых тенерь */+2 — * (р, Р/+2)» 7*4* = т(Р> р|+2)«

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

Соответствующие изменения в формулах численного метода вполне аналогичны описанным в [I].

3. Изложенный численный метод проиллюстрируем на примере расчета конического течения воздуха около нижней поверхности треугольного крыла эллиптического поперечного сечения (отношение большой и малой полуосей ajb = 10) с углом стреловидности у =50°. 5

Параметры невозмущенного ' потока: число Мто — 20, угол атаки а = 30°, давление и плотность соответствуют высоте Н = 50 км, показатель адиабаты воздуха % = 1,4-

Обтекание рассматривается в декартовой системе координат (х, у, г), центр которой совмещен с вершиной крыла. Ось крыла совпадает с осью х, а большая ось эллипса поперечного сечения лежит в плоскости у = 0. Используется разностная схема [5].

Пространственная расчетная область Q совнадает с одной из двух симметричных относительно плоскости половин возмущенного течения. В качестве

поверхностей начальных данных используется семейство плоскостей х = const, поэтому х — компонент скорости газа — должен быть сверхзвуковым всюду в области Q.

tg Si - g (р, Pu ?и Qi) (Р — Pi)>

(Ю)

Pi _______________________________

(И)

f"« Щ*±±Р+Р1

Т/+2—1

'Ч=~2 + х

где рг+2 берется из предыдущего приближения.

Это дает следующее выражение для плотности:

1

х-1

где

1

* = Т (Л’ pl) + * р2>} •

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

В начальном сечении х = const вблизи вершины крыла расчетная область задается в виде узкой полосы вокруг крыла, газодинамические параметры совпадают с параметрами невозмущенного потока. Расчет ведется методом установления по переменной х.

Поверхность крыла задается в параметрическом виде:

у = bx cos г — ах sin ЧГ, |

0^ЧГ<те, tf = ctgx, b = 0,1 a, j где *£'—угол относительно вертикальной оси в плоскости х*= const.

Рис. 3 Рис. 4

Пространственная разностная сетка строится путем разбиения поперечного сечения Q на элементарные четырехугольники (см., например, [2]). Одним семейством сеточных линий служат лучи, исходящие из точек на поверхности крыла, задаваемых формулами (12), в которых угол l* принимает значения

W=k(nlK)> k — 0, 1,... , К.

Нанравления лучей совпадают с нормалью к эллипсу поперечного сечения.

Одна из сеточных линий второго семейства совмещена с поверхностью крыла, другая—с головным скачком. Остальные сеточные линии второго семейства равномерно разбивают отрезок каждого луча между внутренней и внешней границами на заданное число N интервалов. Точки разбиения являются узлами сетки. В приводимых ниже расчетах числа А/ и К равны соответственно 8 и 32.

На рис. 2 и 4 сравниваются расчетные характеристики при обтекании крыла воздухом с учетом реальных свойств (сплошные линии) и при обтекании совершенным газом при %=1,4 (штриховые линии).

На рис. 2 показано положение головного скачка на наветренной стороне крыла в сечении х = const (масштаб по оси у вдвое больше, чем по оси г, величина z отнесена к местному значению полуразмаха крыла). В совершенном газе угол наклона головного скачка больше, чем в реальном. Причиной является уменьшение эффективного показателя адиабаты в реальном газе, как это видно на рис. 3, где представлено распределение функций х и 7 по нижней поверхности крыла в плоскости х — const.

Как следует из рис. 4, почти на всей нижней поверхности крыла коэффициент давления в реальном газе приблизительно на 3% ниже, чем в совершенном. Исключением является область пнка давлении в окрестности кромки, где в случае реального газа давление нарастает более интенсивно и достигает большей величины, чем при обтекании совершенным газом.

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

Аналогичный результат при несколько отличающихся параметрах невозмущенного потока и для эллиптических конусов с другим отношением полуосей получен в работе [6].

ЛИТЕРАТУРА

1. Иванов М. Я., КрайкоА. Н., Михайлов Н. В. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. „Ж. вычисл. матем. и матем. физ.*, 12, № 2, 1972-

2. Годунов С. К., Забродин А. В., Иванов М. Я., Крайко А. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики. М., „Наука', 1976.

3. ГолубинскнйА. А., К о с ы х А. П. К расчету стационарных сверхзвуковых течений с учетом реальных свойств газа. „Ученые записки ЦАГИ“, т. IX, № 4, 1978.

4. Дьяконов Ю. Н. Пространственное обтекание затупленных тел с учетом равновесных физико-химических реакций. ДАН СССР, т. 157, № 4( 1964.

5. 3 а р у б и н А. Г. О точности метода Годунова в различных системах координат. „Ученые записки ЦАГИ", т. XIII, № 4, 1977.

6. Б а з ж и н А. П,} Трусова О. Н., Челышева И. Ф. Влияние реальных свойств воздуха на параметры течения около эллиптического конуса. Аэродинамические характеристики эллиптических конусов при больших углах атаки. „Ученые записки ЦАГИ*, т. I, № 2, 1970.

Рукопись поступила Щ1Х 1980 г.