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

Профилирование сопла с центральным телом

Ким Ч.В fcwkim@mail.ru )

Московский физико-технический институт

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

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

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

Когда энтропия S = const, уравнения состояния задаются в виде :

T

E (T) = J Cv (T )dT (1)

о

T

p(T) = p(T )*(T), K(T) = J R(T)dT (2)

о

Несовершенный газ с уравнениями состояния (1, 2) определен как «квази-совершенный». Зависимости Cp(T), Cv(T) заданы формулами : Теплоёмкость смеси : Cp, Cv, к

Cp 1 = A1 + B1 • T, Cp2 = A2 + B2 • T , Cp 3 = A3 + B3 • T, Cp4 = A4 + B4 • T • 0.37 Cp = У1 • Cp 1 + У2 • Cp 2 + У3 • Cp 3 + У4 • Cp 4

С1 = ср 1 - 1, С 2 = СР 2 - 1, С 3 - Ср 3 - 1, Су 4 - Ср 4 -1

С = У1 ■ С 1 + У2 ■ С 2 + У3 ■ С 3 + У4 ■С 4

к = Ср / Су

А1 = 5.30886 В1 - 0.29682 У1 - 0.1838

А2 = 3.60820 В2 - 0.35166 У2 - 0.02

Аз = 3.78380 В3 - 0.11131 У3 - 0.0102

А4 = 3.35202 В4 - 0.14020 У4 - 0.786

Т = Т/273 °К Ср = Ср / 7 Ъо С = Су / 7 Ъо К =8.314 1/шоЬК

Уравнения осесимметричного безвихревого течения идеального газа в ортогональной системе координат (у,ф), связанной с линиями тока, представлены в [1]. Обозначив 1п Q(Л) = | (М 2 - 1)(1п Л)ф dф , получим окончательное представление уравнений движения

^АГ^ , М<АТ)^ (3)

Л ^ А,То) ■ Ъъ * у^ А,То) ■ Ъъ

где Л = |К|, в- угол наклона вектора скорости по отношению к оси О, М = Л / а,

Н3 = 1, у для плоского и осесимметричного случаев соответственно.

Итак, уравнения стационарного изэнтропического движения несовершенного газа инвариантны относительно уравнений состояния (1, 2).

В уравнениях (3) можно совершить преобразование Чаплыгина [2] в плоскость годографа скорости точно так же, как и в уравнениях для совершенного газа [3, 4], путём умножения якобианов (д( ф, у) /д(Л, в)).

Перекрёстно дифференцируя выражения фЛ и фв получим окончательное уравнение для функции тока ¥(Л, в) в несовершенном газе.

д у (1 + М2) ду (1 -М2) д> дЛ + Л ~дЛ+ Л

n<

y Ш2 + (1 - M2 ) y W2 _ СОЗв + 2 Sin в

У у* Л ^ л2 Уу в

yjrt_ ssne D

в

(4)

yA * p yA y Л ^ p уЛ

где D =(д( ф, у) /д(Л, в)) и n = 0, 1 для плоского и осесимметричного случаев соответственно.

В плоском течении при подстановке Л^Х, Q(A,T0) ^ q(X) система (3) переходит в систему уравнений для совершенного газа. (Х - коэффициент скорости, газодинамическая функция q(X) выражает удельный расход в совершенном газе). В осесимметричном течении при соотношениях т =

\V\2/\V\2 = ^/2CpTa, т* = \V* f/\V\2 и Q^o^ = Q^t^^ □ P*/ Po уравнение (4)

/ max J f max 1 1

переходит в обобщенное уравнение Чаплыгина для совершенного газа.

ДОЗВУКОВАЯ ЧАСТЬ СОПЛА

Корректная постановка задачи профилирования сопла Лаваля методом годографа [4] не претерпевает никаких изменений в случае несовершенного газа.

$

в ¥=1 С

А D

Ч>=0

Рис. 1 Рис. 2

Прямоугольнику АБСБ (Рис. 1) в плоскости годографа соответствует в физической области контур сопла abcd (Рис. 2). Задача Дирихле в прямоугольнике (Рис. 1) для уравнения (4) решена численным методом. После решения для получения координат сопла проводилось интегрирование вдоль

линий ВС, СВ и ВА в плоскости годографа. Результаты расчётов представлены на рис. 3.

На рис. 3. Контуры сопел для разной ширины щели в критическом сечении при фиксированной наружной обечайке.

На рисунке сплошной линией показан контур, полученный при То -1000°К, пунктиром - при То - 3000°К, а треугольником - при То -5000 °К. Все

контуры почти совпадают.

т о т* скорость (ш^) соотношение удельного расхода на входе и выходе соотношение высоты на входе и выходе (^/Ь*)

на входе на выходе

1000 859.97 36.025 97.449 0.3462 1.8

2000 1745.14 50.631 137.066 0.3485 1.8

3000 2649.65 61.672 167.091 0.3499 1.8

4000 3569.30 70.893 192.177 0.3519 1.8

5000 4500.90 78.932 214.091 0.3530 1.8

Таб. 1

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

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

Рис. 4. Контур дозвуковой части осесимметричного сопла при То = 3000°К. Разные кривые соответствуют различным отрезкам звуковой линии.

СВЕРХЗВУКОВАЯ ЧАСТЬ СОПЛА

Схема вычисления течения для несовершенного газа не отличается от схемы для совершенного газа [3], однако, алгоритмы на каждом этапе разные. Сверхзвуковая часть контура в несовершенном газе вычисляется в четыре этапа, как и в совершенном газе.

А' А''' Б

\\\\\\\\\\\\\ч\\\

Е

Рис. 5 Область определения течения в физической плоскости х, у.

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

На втором этапе возникают некоторые изменения от случая совершенного газа в алгоритме вычисления характеристики второго семейства А "Б в течении Прандтля-Майера вблизи угловой точки.

Используя условия совместности вдоль характеристики второго семейства можно найти число Маха М, в точке с известным числом Маха М^ предыдущей точки на характеристике. Число Маха М, должно удовлетворять следующему уравнению.

I

V-1

1

V2 ^V

1--атсзт

Г 1 Л

Ср(Т)п„у V

Я(Т)

V ,

+ агсзт

' л.л

М-1У

в =0 (5)

С„(Т)

где в = (а + в ¿.¡) - (а , + в ¡), а - угол Маха. Уравнение (5) решается методом проб и ошибок. При этом возникает трудности с выбором шага интегрирования. По мере уменьшения шага интегрирования увеличивается точность, с другой стороны объём вычислений значительно возрастает. Такая же трудность возникнет при вычислении методом характеристик на каждом этапе. При проведении расчёта в области А "'А "БЭ (между двумя характеристиками, полученными на первых двух этапах) сначала вычисляется в характеристическом четырехугольнике А "'А "БС, а затем в треугольнике А'" СЭ, ограниченном цилиндрической стенкой.

На третьем (вычисление решения в области А "'А "БЭ) и четвёртом этапах (вычисление контура центрального тела - участка АЕ) расчёты проводились методом характеристик с использованием следующих формул :

=

У-1- У-и + & -1- в- -1 X -1 +{а1-1/ + Д.-,,, ъ (а.-1- в-1)+ # к-и+Д-и)

У/,.. = Ум,; + & (а/-1,. + Д-1,/ Xх/,;- - х-и )

Вычисление числа Маха М'. и угла наклона скорости в'. в каждой точке можно осуществлять двумя подходами.

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

в = 2

V-

А-,, +Р,;-1 + |

V

-1

1

V2

С(Т) С(Т) ' 7

, (IV у 1-+

V

I

Ч;-1

1

V2

С(Т) СУ(Т) ' 7

сIV ~У~

в;-1 -в-1. = I

-1

1

V2

Ср(Т)

, С^ 'ч 1 — +

ЩТ)

V

I

1

V2

Ср(Т)

1

С(Т)

ЩТ)

су

~У~

(7)

С(Т)

При решении двух уравнений основная погрешность в расчёте связана с величинами шага V';, Т,- и шагов сеток.

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

Во втором подходе используется соотношение между касательной к характеристике первого семейства, проходящей в плоскости (х, у) через одну точку М и касательной к характеристике второго семейства, проходящей через соответствующую точку М' в плоскости (Уx, Уy), (ось параллельна оси Уx).

При этом получается Уij, решая уравнения Уx Уy у :

1

Ухг, =

УУг-1, / - УУг , /-1 + КА-1,/ - К2У\/-1 + (х! , / - /-1) (х! ,, / - Х,--1, / )

2 \ лг ,/

Уу,/ = УУ-1,/ - К2(Ухи - Ух<-1,/) К =

' Л-1,] 2' г,] г-1,/7

1 Ух2 - а2

К1 - К2 (хг,/ - хг-1,/)

1 Ус1 УхУу + а^Ух2 + Уу2 -1 Ух2 - а2

2 2 2 а

К2 =

01 =

02 =

а2УУ1УС /

У1(УУ12 - а2 ) а2Уу2 Ус^

(8)

УхУу - а^Ух2 + Уу2 - а2 2 У2 (Уу1 - а2 )

Ус2

где £ = 0, 1 относятся к плоскому и осесимметричному случаям, а индексы С1, С2 - к характеристикам первого семейства и второго семейства

соответственно.

Координаты контура (хг/, у/ центрального тела могут быть получены вычислением пересечения линии характеристики первого семейства от точки ху,-.] с линией (от точки хг-],уг.]/) с углом наклона в.] / вектора скорости относительно оси х.

|\ х,ч,, у,^ /

\vall |\ /

1 СЦ-Ц | / ^у' У у

Рис. 6

Число Маха и угол наклона вектора скорости найдётся интерполяцией со значениями М1/.], и М, в в точке х , у на линии воображённой характеристики второго семейства от хг-]уг.]/ (Рис. 6).

Теперь рассмотрим влияние температуры торможения на контур сопла. На рис. 7 показаны результаты расчёта при разных температурах торможения. В сверхзвуковой части сопла увеличение температуры торможения приведёт к росту скорости на выходе сопла, а к уменьшению числа Маха.

Линией показан контур сопла при температуре торможения Т0 = 1000°К, пунктиром - при Т0 =3000 °К, треугольником - при Т0 =5000 °К. В этих случаях

скорости потока, скорости звука и число Маха на выходе сопла написаны в таблице 2.

Рис. 7. Контуры сверхзвуковой части осесимметричных сопел с центральным телом.

Т 1 о скорость потока Уo скорость звука ао число Маха Мо

1000 171.82 78.86 2.179

3000 299.82 140.41 2.135

5000 389.12 184.76 2.106

Таб. 2

Как показаны расчёты на рис. 7 и таб. 2, течение сверхзвуковой части зависит от температуры торможения Т0. Поэтому сопло, спроектированное под какую-то температуру торможения, при её изменении окажется в нерасчётном режиме. Иначе говоря, если сверхзвуковой часть оказывается короче, чем нужно, то по теореме А.А.Никольского возникнет скачок уплотнения, что приводит к дополнительным потерям [5].

Дальше проводилось сравнение контура для несовершенного газа с контуром для совершенного газа.

На рис. 8. Линией показан контур сопла при температуре торможения Т0 = 3000°К для несовершенного газа, пунктиром - для совершенного газа с показателем адиабаты к = 1.3. В последнем случае число Маха на выходе сопла состоится 2.062.

СПИСОК ЛИТЕРАТУРЫ

[1] Шифрин Э. Г. Потенциальные и вихревые трансзвуковые течения идеального газа. М: Физматлит, 2001.

[2] Кочин Н. Е., Кибель И. А., Розе Н. В. Теоретическая гидромеханика, часть вторая. Государственное издательство физико-математической литературы, Москва 1963.

[3] Шифрин Э. Г, Ким Ч. В. Профилирование сопла с центральным телом методом Чаплыгина. Доклады Академии Наук, Т. 401 № 1, 2005.

[4] Подсыпанина Н. А. Использование плоскости годографа при профилировании численным методом осесимметричного сопла Лаваля. Изв. АН СССР, МЖГ, № 1, 1977.

[5] Никольский А.А. О течениях газа вблизи остроконечных задних кромок тел вращения. Сборник теоретических работ по аэродинамике, М., оборонгиз. 1957.