Численный Расчет течения вязкого газа в плоском сопле ГПВРД Текст научной статьи по специальности «Физика»

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

Том XXIV 1993

№ 3

УДК 629.7.015.3.036:533.697.4 629.7.036.22.016.55

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

А. П. Мазуров

На основе неявной разностной схемы ТУО разработан численный метод расчета параметров турбулентного течения в плоском несимметричном сопле ГПВРД при больших степенях нерасчетносги истечения реактивной струи. Для описания движения газа используются осредненные по числу Рейнольдса уравнения Навье — Стокса, записанные в приближении тонкого слоя. Изменением термодинамических свойств и химического состава продуктов сгорания вдоль тракта сопла пренебрегается. Приведены примеры расчетов течения в плоском сопле ГПВРД с профилированным клином при Ми »6 и 10. Проанализированы влияние располагаемой степени понижения давления на картину течения в сопле и его тяговые характеристики.

Для гиперзвуковых летательных аппаратов (ГЛА) в качестве основной двигательной установки рассматривается гиперзвуковой прямоточный воздушно-реактивный двигатель (ГПВРД) в компоновке с плоским несимметричным соплом (см., например, [1]). Большие сложности моделирования условий гиперзвукового полета в аэродинамических трубах существенно ограничивают возможности экспериментальных исследований на моделях. В связи с этим значительно возрастает роль численных методов при определении аэродинамических характеристик ГЛА и его двигательной установки.

С ростом скоростей полета происходит заметное увеличение вклада тяговых характеристик реактивного сопла в суммарную эффективность двигательной установки с ВРД. Поэтому одним из важных требований, предъявляемых к численным методам расчета течения в соплах ГПВРД, является высокая точность определения тяговых характеристик сопла.

Расчеты и анализ характеристик плоских несимметричных сопл ВРД в приближении невязкого газа выполнены в работах [2 — 4]. В работе [5] получено решение вариационной задачи о построении контура оптимального плоского сопла с учетом аэродинамических характеристик летательного аппарата.

В последнее время для численного моделирования течения в плоских соплах ГПВРД стали применяться разностные методы решения уравнений Навье — Стокса (см., например, [6, 7]). Результаты, полученные с помощью таких методов, показывают, что использование уравнений движения вязкого газа позволяет рассчитывать очень сложную структуру потока в сопле ГПВРД на различных режимах истечения, включая отрыв пограничного слоя.

В настоящей статье для численного интегрирования уравнений движения газа применяется неявная разностная схема [8]. Турбулентная вязкость в пристеночных пограничных слоях и в слое смешения рассчитывается по формулам вихревой модели турбулентности [9].

1. Течение в сопле рассматривается в предположении, что термодинамические свойства и химический состав продуктов сгорания ГПВРД не изменяются вдоль тракта сопла, т. е. рассматривается течение совершенного газа с постоянным показателем адиабаты *. При этом осред-ненные по числу Рейнольдса уравнения Навье — Стокса, записанные в приближении тонкого слоя в преобразованных пространственных переменных £ = ¿¡(х), п=‘п(х,у), имеют вид

ÔU âF дС ât д$ drj

1 âW Re ât]

(1)

где

U = J-'Q, Q =

/ p / \ РиУц Г pK 1

pu , F = (pu2 + p)yr) » G = puV -yçp

po 7 pouyn puV + х^р

\E J (E + p)u y , (E + p)V y

W =

0 ^

" a2un a3o4-a2un

{a4un-a50n + a6en)

р = (х-1)(Е-ри2+?и2), J-I=xçy4,

V —Xç — uy

Здесь и, о — продольная и поперечная составляющие скорости в декартовых координатах х, у; р — плотность; р — давление; е — внутренняя энергия единицы массы; Е—полная энергия единицы объема; Re — число Рейнольдса; xç,y^yn — метрические производные; t—время.

Коэффициенты a; (i = 1 - 6) , входящие в выражение вектора вязких членов W, имеют вид

Oj = J; а2 - -^nJxçyÿ а3 = + у\ j:

a4 =0]U-a2y; a5=a3v-a2U; ae=kJ(x^+yp,

где ц и к — коэффициенты вязкости и теплопроводности. Для турбулентного течения коэффициенты вязкости и теплопроводности берутся

В виде сумм Ц = Цт+Цт И к = + где Цт —коэффициент

Рг1 Ргт

ламинарной вязкости; цт—коэффициент турбулентной вязкости; Рг£ =0,72 — число Прандтля для ламинарного течения; Рг,. =0,90 — для турбулентного. Коэффициент турбулентной вязкости цт вычисляется по формулам алгебраической модели турбулентности [9]. Все величины, входящие в уравнение (1), записаны в безразмерном виде и отнесены к параметрам невозмущенного потока и в качестве харак-

терного линейного размера взята высота миделя сопла Ам.

Граничные условия на твердых стенках при решении уравнений (1) ставятся с учетом условия прилипания при заданной температуре стенок. Стационарное решение находится с помощью численного метода в процессе установления при стремлении времени / к бесконечности. Метод разработан на основе разностной схемы ТУБ. Схемы данного класса правильно отслеживают распространение возмущений в потоке и являются монотонными без добавления искусственной диссипации. Это их свойство особенно важно при расчете течений в соплах ГПВРД, в которых истечение реактивных струй происходит, как правило, на нерасчетных режимах, в частности при очень больших перепадах давлений.

Численное решение уравнений (1) выполняется с помощью приближенно факторизованной разностной схемы

где — разностные операторы нестационарной части уравнений (1)

соответственно вдоль направлений £ и 77; АС} = (}п+] - (}п —приращение вектора искомых переменных за временной шаг Аґ; Ь,^п—разностный аналог стационарной части уравнения (1) на п-м временном слое:

(2)

/

^ (

\

КеАт/

V

Здесь целые индексы /, у относятся к центру ячеек разностной сетки; полуцелые индексы / ±^ и ] ±~ относятся к граням ячеек; А£, Ат) —раз-

65

меры ячеек в расчетной области соответственно вдоль направлений £ и т/;

Л Л Л _

Рп, Сп и |уп обозначают разностные аналоги векторов Рп, С” и 'ЧРп.

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

ячеики I + - используется выражение

. 1

1+~

2

г+- г+-

. 2 2

\ . (<?+, -<?-,)

1+- ¿+- й—

2 2 2 .

Здесь

А . »+2

матрица, которая с помощью преобразования по-

добия может быть приведена к диагональной матрице с элементами, равными абсолютным значениям собственных чисел якобиевой матрицы А = дР / д(?; (?+ , <3~ — значения вектора (} соответственно спра-

¡+- »+2 2

ва и слева от 1рани * + ~, вычисленные с помощью соотношений [8].

Стационарное решение уравнения (2) находится в процессе установления по времени с использованием векторных прогонок при обращении нестационарных операторов и

2. Разработанный метод был применен для расчета течения в плоском сопле ГПВРД с учетом взаимодействия реактивной струи с внешним потоком. Схема течения и геометрия сопла в декартовых координатах х, у показаны на рис. 1, а. Расчеты выполнены при следующих значениях параметров: длина сопла хс =2,50; абсцисса точки излома нижнего контура х0 =0,30; длина обечайки хоб = 0,60; ординаты образующих нижней и верхней стенок сопла /^ = 0,90 и Ас = 0,98; высота входа в сопло квх = 0,08 (размеры даны в долях высоты миделя сопла /гм); контур клина сопла образован полиномом второй степени, причем

угол наклона образующей клина в точке излома равен 30°. Отношение

высоты среза сопла к высоте входа составляло Лс = Ас / йвх =12,25.

При расчетах принимался одинаковый в струе и во внешнем потоке показатель адиабаты х= 1,2 (поскольку методика применима только для однокомпонентной среды); температура невозмущенного потока составляла Тм = 218К; число Рейнольдса, определенное по параметрам

невозмущенного потока и высоте миделя сопла, Ие = 1,5 • 106. При числе

о

6) 14

0,5

1.0

1,5

2,0

1П

1Е

Рис. 1. Геометрия плоского сопла ГПВРД: а — схема течения; б — разностная сетка

Мю = 6 температура торможения в струе равнялась Т0с = 2800К, число Маха на входе в сопло Мвх = 1,50; при М«, = 10 соответственно 7’0е= 3800К, Мвх=1,92. Относительные температуры наружной (Тк =Т„ / Тж) и внутренней (ТКс = Тш / Тт) поверхностей сопла соответствовали 4,4; 5,4 и 6,4. Располагаемая степень понижения давления я'ср = р0с / рх (р0с —полное давление в реактивной струе, рж —давление невозмущенного внешнего потока) изменялась от 100 до 600 при Мю = 6 и от 100 до 1200 при Мте = 10. Толщины пограничного слоя во входном сечении как внутри сопла, так и снаружи составляли ¿=<?с=0,02 (см. рис. 1, а).

Расчеты выполнялись на неравномерной разностной сетке, связанной с контуром сопла (рис. 1, б). Общее число узлов сетки равнялось

Шх. Ш = 73 х 59, при этом в поперечном сечении проточной части сопла содержался ж = 31 узел, в продольном направлении до кромки обечайки Ш = 32 узла. Во внешнем потоке распределение узлов в поперечном направлении задавалось с помощью соотношения

л] = л(6) + & = ¿¿¡У-1), V] = А чи ~ 1).

где >>„(#,•) — ордината образующей наружной поверхности обечайки в точке £,•; Н—вертикальный размер расчетной области; функция <р(т)) имеет вид [11]

_ /?+1-(/?-!)ехр[с(1- 7)]

1 + ехр[с(1 - 7)]

Здесь р — коэффициент сгущения узлов, коэффициент с находится из условия <р(т]л) = 0 и равен

-- 1 -ш-^1

1-*д Р~1

где т)л - Лт^ЗЬ-1)—координата 77, соответствующая наружной образующей обечайки и номерам узлов _/ = л,.

Распределение узлов сетки поперек проточной части сопла задавалось следующим образом:

Ух} = УЛ&) + 0,5[ун(й) - Ус(6)]Р1(*7/)> 0 ^ Ч) ^ *7°,

Уу = Ун(6) + 0,5[ун(й) - >'С(^)]<Р2(^), тр гик,

гае ус(&), Ун(&) — ординаты образующих клина и внутренней поверхности обечайки в точке £,•; 77О — значение поперечной координаты 77, соответствующее линии стыковки узлов; функции ^(77) И <р2(Т]) имеют вид

, V Рс + 1 - (^С - 1)ехр[с!(7° - 7)1 1 р + 1

91(4)=-------------Г----п --1---- ’ С1 = “(Г ---->

1 + ехр[с! (Тр - 7)] 7° Рс-1

■ Рс + 1-(/?с -1)ехр[с2(7- 7°)1 1 я+1

<р2{ Т}) =---------=----!=--=----с2=--------------1п-2----.

1 + ехр|с2(7- 7°^ Чж - V Рс-'

В расчетах принималось тР =0,Ъг]Ж и соответственно

С1 = с2 - —— 1п — +1. Распределение узлов в продольном направлении

Члс Рс -1

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

3. Рассмотрим влияние на картину течения и тяговые характеристики сопла располагаемой степени понижения давления лср. Для сопла

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

На рис. 2 для примера показана картина течения в сопле при

М«, = 6 и лср = 300. Сплошными кривыми здесь нанесены уровни чисел М, штриховой кривой обозначена условная граница реактивной струи, проходящая примерно посередине слоя смешения. При возрастании

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

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

; %ср-т , Тиг'Тгс-*,*

поэтому в рассмотренном течении нет четкого разделения на области невязкого и вязкого потоков.

Дополнительную информацию о структуре течения в плоском сопле с косым срезом дает рис. 3, на котором показано поле векторов скорости при Мда =6 и д-ср = 300. Поскольку давление в реактивной струе выше, чем в набегающем потоке, то за кромкой обечайки линии тока отклоняются в сторону внешнего потока. Возрастание ?гср приводит к увеличению углов наклона векторов скорости в области взаимодействия и оставляет неизменным поле векторов скорости в области, прилегающей к поверхности клина. При Мда =10 качественная картина течения имеет такой же характер. Отличие состоит лишь в том, что из-за разных значений числа М на входе в сопло при одних и тех же значениях 7гср возмущенная область в спутном потоке при Мю = 10 имеет меньшие поперечные размеры, чем при Мм =6.

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

на рис. 4. Здесь для Мж = 6 и Т„ = Т„с =5,4 дано сопоставление профилей продольной скорости и / м«, и температуры Т / Тж в сечениях х = 0,69; 1,25 и 1,85 при яср = 100, 300 и 600. Основное отличие профилей наблюдается в слое смешения и в ударном слое над ним. Вблизи поверхности клина по всей длине сопла профили и / и„ и Т / Тх при разных значениях тгср практически сливаются. Отметим также, что при увеличении лср происходит снижение температуры в слое смешения. Выше слоя смешения температура в ударном слое превышает температуру невозмущенного потока вследствие сжатия в ударной волне.

Рис. 4. Распределения профилей продольной скорости и

температуры

На рис. 5 приведены распределения относительного статического давления руу / р0с вдоль образующих клина и внтуренней стенки обечайки, полученные в расчетах при Мда =6, лср = 300 (сплошные кривые) и Мм = 10, • 2гср = 900 (штриховые кривые). Отметим, что в диапазоне значений лср = 100 -*■ 1200 статическое давление на стенках сопла, отнесенное к полному давлению в реактивной струе, практически одно и то же. За точкой излома контура клина наблюдается резкое падение давления, после чего существует участок плавного снижения давления (при Мм = 10 и Мвх = 1,92 за точкой излома наблюдается плато давления); примерно с половины длины клина давление на нем изменяется довольно слабо. Различие давлений на стенках сопла при Мте = 6 и 10 связано с неодинаковыми значениями числа Маха на входе. В случае Мм = 6 статическое давление вблизи среза клина равно давле-

ГІГ'.РжЬн 0,1*

1,11 в,т до до до до

О V 1,9 1,5 Іо *

Рис. 5. Распределения давления вдоль образующих клина и обечайки сопла ГПВРД

нию невозмущенного потока примерно при 7гср = 110, а для Мте =10 при яср =! 170, Однако для сопла с косым срезом, в отличие от обычных симметричных сопл, указанные режимы не являются расчетными, поскольку давление на обечайке с внутренней стороны заметно отличается от давления невозмущенного потока. Так, относительное давление Ру, / р0с на срезе обечайки равно 0,020 при М«, = 6 и 0,014 при М,,, = 10. Отсюда следует, что в первом случае давление на срезе обечайки при тгср =110 превышает давление невозмущеннош потока в 2,20 раза, а во втором при = 170 в 2,38 раза. Таким образом, для несимметричного сопла ГПВРД довольно проблематичным становится выбор режима истечения, при котором давления на срезе клина и на срезе обечайки были бы одновременно близки к давлению невозмущенного потока.

Особенностью расчета течения вязкого газа в сопле ГПВРД является повышение статического давления на начальном участке проточной части сопла (см. рис. 5). Например, для исходных данных, соответствующих числу М«, =6, это повышение составляет примерно 10% от давления на входе в сопло. Поскольку проходные сечения на этом отрезке сопла не изменяются, то полученный результат связан непосредственно с влиянием трения и теплопередачи на стенках сопла.

Контрольные расчеты, выполненные на трех сетках с разным количеством ячеек поперек проточной части сопла (ЇК= 21, 31 и 41) при одинаковом суммарном числе ячеек в расчетной области (/£ = 59), показали хорошую сходимость распределений давления вдоль образующих сопла, а также профилей скорости в поперечных сечениях, за исключением узкой области, прилегающей к стенкам.

Рассмотрим поведение тяговых характеристик сопла в зависимости от располагаемой степени понижения давления в реактивной струе. На рис. 6 приведены результаты расчета коэффициента продольной тяга

Рис. 6. Зависимости угла отклонения и коэффициента продольной составляющей вектора тяги сопла ГПВРД от располагаемой степени понижения давления

сопла РХс = РХс / Рсил (Рхс — продольная тяга сопла, Рс ид — идеальная

Рус

тяга) и угла отклонения вектора тяги у/с = агс^-р- (Рус — поперечная

Р% с

тяга сопла) при числах М*, = 6 и 10 для температуры стенок сопла, равной 5,4 Т^. Видно, что при повышении яср значение Рхс изменяется немонотонно: сначала возрастает, затем, после достижения максимума, снижается. При Мш = 6 максимальное значение Рхс составляет примерно 0,960 при /гср = 200, при М«, = 10 этот максимум равен 0,95 при *гср = 300.

Потери продольной тяги сопла АРХС =1- Рхс в точках максимума

кривых Рх с = /(^ср) обусловлены главным образом потерями из-за рассеяния потока на срезе сопла (с учетом взаимодействия струи с внешним потоком). Потери тяги левее максимума вызваны перерасши-рением сопла, правее — недорасширением.

В отличие от коэффициента продольной тяги угол отклонения вектора тяги не имеет экстремума в рассмотренном диапазоне значений лср.

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

случая Мж = 6 при тгср = 100 угол отклонения вектора тяги равен 5,0°,

а при 7гср = 600 равен 3,8°. Это указывает на то, что поперечная тяга

сопла ГПВРД в рассмотренном диапазоне яср может быть как положительной, так и отрицательной. (Отметим, что результаты приведены для случая нижнего расположения клина, для верхнего расположения клина знак угла следует поменять на обратный.)

Расчеты тяговых характеристик при различных значениях температуры стенок сопла показали, что охлаждение стенок •приводит к некоторому снижению коэффициента продольной тяги. В частности, при Мм = 6 и яср = 400 значение Рхс соответствует 0,952 при относительной

температуре стенки Т = 6,4 и 0,947 при Тш = 4,4. Причем изменение

Рхс в рассмотренном диапазоне 4,4 ^ Тш й 6,4 близко к линейному.

Наконец, рассмотрим одно следствие, вытекающее из того обстоятельства, что на режимах недорасширения безразмерное давление Ру» / Ро с на поверхности клина практически не зависит от располагаемой степени понижения давления в сопле. Оно состоит в том, что,

зная коэффициент продольной тяги при одном значении я£\ можно рассчитать его значение при другом я^К Формула пересчета коэффициентов продольной тяги Р^1 и Р^ при значениях я^ и я™ получается из условия независимости от яср интеграла сил безразмерного давления Рк / Рос, взятого по поверхности клина, и имеет следующий ввд:

І” Ж(2)

р(2) _ р(1) ИД | кс Ф _ |

Х° ХС Л(2) Р(2) я-(,)

(3)

ид с.ид \ ср у

где р(2) =х

2

— идеальная тяга при

с.ИД

Х+1

\

Р ср > ИД , ж_ і

1-

п,

1 *

идеальная приведенная скорость;

д (Авх) — функция приведенного расхода на входе в сопло.

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

Пример пересчета коэффициентов продольной тяги по формуле (3) показан на рис. 6 пприхпунктириыми линиями. При Мда = 6 пересчет

производился при значении я^ = 400, при котором = 0,949, а при

Мда = 10 я^ = 600 и = 0,944. Видно, что на режимах недорасшире-

ния результаты пересчета хорошо согласуются с результатами численных расчетов. Таким образом, ігри оценке потерь продольной составляющей тяги сопла ГПВРД при больших располагаемых степенях понижения давления достаточно знать потери при одном-двух значениях я^.

ЛИТЕРАТУРА

1. К у р з и н е р Р. И. Реактивные двигатели для больших сверхзвуковых скоростей полета,—М.: Машиностроение. — 1989.

2. Бер ля нд А. Т., Я гуд ин С. В. Расчет истечения невязкого газа из плоского сопла с учетом влияния внешнего потока//Труды объединенных научных чтений по космонавтике/Двигатели летательных аппаратов. — М.: Изд. ИИЕТ АН СССР. — 1980.

3.Тагиров' Р. К., Левин М. П. Расчет характеристик плоских несимметричных сопл//Ученые записки ЦАГИ. — 1981. Т. 12, Mi 6.

4.Ягудин С. В. Влияние насимметрии плоского сверхзвукового течения невязкого газа на характеристики плоского сопла в статических условиях// Ученые записки ЦАГИ. — 1983. Т. 14, № 5.

5. Рыло в А. Н. О построении компактных несимметричных сопл максимальной тяпс при заданной подъемной силе//Изв. АН СССР, МЖГ. — 1981, №6.

6.Harloff G. J., Lai Н. Т., Nelson Е. S. Two-dimensional viscous flow computations of hypeisonic sc ramjet nozzle flowfields at design and ofF-design conditions//AIAA Paper, 88 — 3280.

7.Ruffin S. М., Venkatapathy E., Keener E. R., Nagaraj N. Computational design aspects of a NASP nozzle/afteibody experiment//AIAA Paper, 89 — 0446.

8. Woan C. J., Chakravarthy S. R Transonic Euler calculations of a wing-body configuration using a high-accuracy TVD scheme//AIAA Paper, 88 — 2547 CP.

9.Baldwin B. S., Lomax H. Thin-layer approximation and algebraic model for separated turbulent flow//AIAA Paper, 78 — 257.

10.Roe P. L. Approximate Riemann solvere, parameter vectors and difference schemes//Journal of Computational Physics. — 1981. V. 43.

11.Holst T. L. Numerical solution of axisymmetric boattail fields with plume simulators//AIAA Paper, 77 — 224.

Рукопись поступила 3/IX 1991г.