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

Том XЬїї

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

№ 4

УДК 532.516.5: 533.596.8

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ПРОСТРАНСТВЕННЫХ ТЕЧЕНИЙ

ВЯЗКОГО ГАЗА В СОПЛАХ

А. П. МАЗУРОВ

Разработана методика численного моделирования турбулентных течений газа в трехмерных соплах произвольной формы на основе осредненных по Рейнольдсу уравнений Навье —

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

Для вычисления турбулентной вязкости применяется двухпараметрическая к-є модель турбулентности, учитывающая влияние стенки. Численное интегрирование уравнений движения производится с помощью неявной разностной схемы третьего порядка аппроксимации конвективных членов. Приведены примеры расчетов течений в трехмерных соплах с разной формой поперечных сечений (круглой, прямоугольной и треугольной). Результаты численных расчетов сопоставляются с экспериментальными данными.

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

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

Ранее численные исследования течений в трехмерных соплах проводились в основном на базе уравнений Эйлера. Значительный объем информации был получен о пространственных течениях в сверхзвуковых эллиптических соплах [4 — 6] и в соплах с несимметричной дозвуковой частью [7]. С появлением быстродействующих компьютеров и эффективных разностных схем стало возможным создание вычислительных алгоритмов для расчета сложных пространственных течений газа на основе полных уравнений Навье — Стокса или их приближений. Численное исследование турбулентного течения в сверхзвуковом сопле прямоугольного сечения, проведенное в работе [14], показало, что результаты, полученные с помощью таких алго-

МАЗУРОВ Анатолий Павлович

кандидат физикоматематических наук, ведущий научный сотрудник ЦАГИ

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

В настоящей работе представлен численный метод расчета пространственных турбулентных течений вязкого газа в сверхзвуковых соплах с произвольной формой поперечных сечений. В качестве математической модели течений используются осредненные по Рейнольдсу нестационарные уравнения Навье — Стокса, записанные в приближении тонкого слоя относительно криволинейных координат. Для вычисления турбулентной вязкости применяется двухпараметрическая к-8 модель турбулентности [8]. Численное интегрирование уравнений движения выполняется по неявной разностной схеме, записанной в виде законов сохранения массы, импульса и энергии для ячейки конечного объема. Конвективные члены в уравнениях движения аппроксимируются с помощью разностей против потока третьего порядка точности по пространственным переменным, члены с вязкостью аппроксимируются центральными разностями второго порядка точности. Решение системы трехмерных разностных уравнений производится методом релаксаций Гаусса — Зейделя по поперечным плоскостям, при котором за одну итерацию выполняются прямой и обратный обходы расчетной области вдоль продольной координаты. Итерации выполняются до тех пор, пока численное решение не станет установившимся. При проведении расчетов контроль сходимости решения осуществляется посредством мониторинга интегральных характеристик сопла (коэффициентов расхода и тяги).

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

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

1. ОСНОВНЫЕ УРАВНЕНИЯ ДВИЖЕНИЯ ГАЗА

Рассматривается турбулентное течение вязкого совершенного газа в пространственном сверхзвуковом сопле произвольной формы, имеющем одну плоскость симметрии. Общий вид сопла и оси декартовой системы координат х, у, г показаны на рис. 1, а. Плоскость хг принимается в качестве продольной плоскости симметрии сопла, и течение рассматривается при у > 0 в области, ограниченной этой плоскостью и поверхностью сопла. В расчетной области физического пространства вводится система криволинейных координат £,, п, ? при помощи преобразования общего вида

% = £(х, у, г), 11 = Г|(х,у, г), С = £(х, у, г) (1)

д(Е л С)

в предположении, что его якобиан J =- —-— нигде не обращается в нуль. Конкретный вид

д(х,у,г)

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

Рис. 1. Схема пространственного сопла (а) и система криволинейных координат (б, в)

В данной работе применяется преобразование пространственных переменных, при котором поверхность сопла является координатной поверхностью г) = const, а координаты с и С направлены в продольном и окружном направлениях соответственно (рис. 1, б, в). Для построения сетки криволинейных координат г) = const, С = const в поперечных плоскостях расчетной области применяется алгоритм, основанный на численном решении уравнения Лапласа [9].

Для описания течения используются осредненные по Рейнольдсу нестационарные уравнения Навье — Стокса сжимаемого газа, записанные в приближении тонкого слоя относительно криволинейных координат г), в следующем виде:

, дО dF dG дН J — + — + — + -

і да.

V

dt д£, дц дС, Re дц

(2)

где

Q =

р и pv

р W

Е

F = F£), G = Р(ц), Я = F(Q, F(P) = J~l

( Р^р > puU^+p$x pvC/p +P$y pwt/p +pfiz {Е + рЩ

Gv =J

0

\x{nuл +тцх) \x{nv^+mx\y) \х.(гтц + mr\z)

(2 2 2 \

U +V +w

+m(T\xu + r\vv + T\zw)

■ Kne_

= PxM + P.yv + Pzw> Р = ^Л,С, Ц =

K = 7

V-L

Vpt

Hr

Prr

Л

п = ЧІ+Ч2у+ЧІ, т = ^(цхиц +ЦуУц +r\zwn).

Здесь и, V, м? — составляющие скорости в декартовых координатах х, у иг; р — плотность; р — давление; е — внутренняя энергия единицы массы; Е — полная энергия единицы объема; у — отношение удельных теплоемкостей; \1Ь и — коэффициенты ламинарной и турбулентной вязкости; рх, Ру, [Зг — метрические производные, где Р представляет с,, 1] и с;: / — время; Рг/ = 0.72 и PгJ^ = 0.9 — ламинарное и турбулентное числа Прандтля соответственно.

Величины, входящие в уравнения (2), представлены в безразмерном виде. Составляющие скорости и, V и м> отнесены к критической скорости в сопле а*, плотность р — к критической

плотности р*, давление р и полная энергия Е — к р;;а?, внутренняя энергия — к а* . В качестве характерной температуры Тя, по которой определяется характерная вязкость |ад , принята вели-

чина

Тк = а}/с,- , где Су — удельная теплоемкость при постоянном объеме. Время I отнесено

к Х/а*, где /, — характерный линейный размер; Яс = ^*а*^ — число Рейнольдса. Давление определяется из уравнения состояния для совершенного газа

р = (у-1)ре .

п

2

Зависимость безразмерного коэффициента ламинарной вязкости от безразмерной температуры принимается в соответствии с формулой Сатерленда:

гр 3/2 1 + С 110К

Рь = т тг—, с = -

Т + с

Тт,

2. МОДЕЛЬ ТУРБУЛЕНТНОСТИ

Величина турбулентной вязкости \1Т определяется с использованием двухпараметрической дифференциальной к-в модели турбулентности [8]. Модель позволяет проводить расчеты в непосредственной близости от поверхности тела и поэтому широко используется при численном решении задач аэродинамики. Уравнения этой модели в приближении тонкого слоя относительно криволинейных координат могут быть записаны в виде

Vе/

до а дР дО дН 1 ,дСт

J -?—!- + — + — +------------= — (——+ £г),

5? 52, дц дС, Яе дц

'Рир1^

^ = ^), а = ад, я = я о, ЯР) = J

(

\

Ст =ги/

I I

У г~г г) )

=Г

.Р-Кере-^^к

у1

2

У

СХ-Р-С2 р —Яе-^е 2

к к у

У П

(3)

Мт Мт

Н*=Ц£ +—, — >

о*

се

Щ =Рхм + Р^ + Рг^5 Р = ^л,С-

Здесь к — кинетическая турбулентная энергия, отнесенная ка*; 8 — скорость диссипации турбулентной энергии, отнесенная к а1 /ь ; Р — скорость генерации турбулентной энергии;

у = — безразмерная координата «закона стенки»; уп — расстояние от стенки;

IIх - ^т1Г/р1Г — динамическая скорость (ти. ри.. |1И. — напряжение трения, плотность и вязкость на стенке).

Эмпирические константы и функции, используемые в данной модели, имеют вид

а£=1.3, су* =1.0,

С1 = 1.35,

с0 =1.8

1-0.22ехр<^ -

Яс,

,РК

где Яс7 = Яе------турбулентное число Рейнольдса.

№

Безразмерная турбулентная вязкость определяется по формуле:

Мт =Кес / р

к

-1

2

6

2

где с(| — коэффициент, равный 0.09; функция /ц =1-ехр(-0.015у+) учитывает демпфирующее влияние стенки (из-за молекулярной вязкости). На большом расстоянии от стенки (в переменных стенки) функция / равна 1.

Скорость генерации турбулентной энергии в приближении тонкого слоя может быть записана в следующем виде:

Р — \хт

9 О 2 1 2

и(ыл +Уп+м>п) + - (Г| хиц + Г1 ,,УЛ + Т1 _М/Л)

ХУ п

\2 ггГ1>

--Яерк(гігил +л>,^1 +лги’л).

3. РАЗНОСТНАЯ СХЕМА

Численное интегрирование уравнений движения (2) осуществляется с помощью неявной разностной схемы, записанной в интегральной форме для ячейки конечного объема (рис. 2). Целочисленные индексы /,7, к связываются с направлениями с. ц. С соответственно и приписываются к центру ячейки, а полуцелые индексы / ± 1/2 , ] ± 1/2, к ± 1/2 относятся к граням ячейки. Уравнения (2) интегрируются по конечному объему ячейки на временном слое Тп + 1 = {п + 1) \/ с

„ „ дО

заменой производной по времени — разностью

дґ

Оп+1-Оп

Дґ

— (0'1+1 -Оп) + Ьг$п+1 + 5 (§'1+1 + 5С#7+1 = — 8 ,

(4)

\п+1

АГ~ ~ ь 11 ь Яе 11

гле нпичинм 8 Э1+1-Э,+1 -Э,+1 8^,г+1-^,г+1 -&п+1 8^,г+1-^,г+1 -&п+1

где величины 0*1* --П+1/2 “п1* _1У/+1/2 °7-1/2> _0К/+1/2 °Г../-1/2>

Ъ-&"+] = &1-+\ 2 ~^к-\12 представляют собой приращения в ячейке конечного объема разностных потоков вдоль направлений г| и V— объем ячейки в физическом пространстве; А/ — временной шаг.

При таком подходе можно отказаться от вычисления метрических производных р.т, Рг. рг (Р = г), С), а использовать векторы площадей граней ячейки 8. Я и Т соответственно в направ-

лениях г) и С, (рис. 2). Тогда конвективные потоки через эти грани могут быть записаны в виде:

(

#Ь) =

Л

рА

р и#ь + рЬх +рЬу + рЬ-(Е + р)ёь у

где $1 = и1,х + уЬу + жЬ: = (и • Ь) — скалярное произведение вектора скорости и и вектора площади Ь = 8, Ы, Т.

Вектор &у вязких потоков через площадь грани Яу + 1/2 в направлении г| может быть записан в виде:

V

./+1/2

0

\х,(Ш^и + МКХ) )

+ МЯ-)

(и и)

Л

2

+ /в(К-и)] + т$Жс/Г|е

1

/€=(!*• II), /ё = —(К-^и),

Рис. 2. Расчетная ячейка конечного объема

1

где V,! | 2 =^(^у+1 + ^"/) — среднее значение объемов соседних ячеек, примыкающих к рассматриваемой грани; с/ц11 — | ~иг с/ц= \,/ | -\,/. и' = | -—разности составляющих ско-

рости; с/ и = (б/п», б/п б/п VI') — разность вектора скорости на грани / + 1/2 ячейки.

Уравнения (4) представляют собой систему нелинейных уравнений относительно искомых составляющих вектора <2Г/£ в ячейке с номерами /, / и к на (п + 1)-м временном слое. Общим

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

уравнений в окрестности решения. Пусть 0* — известное приближение к О"11, тогда для нахождения следующего приближения О' 1 =(2* + AQS используется система линейных уравнений:

(5)

—/ +— 5е#+§ $+5 $--------------§ (§

АГ д(){ 5 11 С Яе

где Я1=-

—03я -0:п) + \р +ьс&

&

1 АІ

—Ь^Оґу представляет собой невязку уравнений (4)

Яе

на 5-ой итерации; I — единичная матрица.

Разностные потоки через грани ячеек при вычислении Щ определяются по методу Ройя [10].

В частности, для потока вдоль направления £, на грани / + 1/2 ячейки с номером / используется выражение

1

т^,2)+т^п)-

лііоє

(бг+1/2 бг+1/г)

/у Н.

где бг+1/2 и (А+і 2 — значения вектора зависимых переменных слева и справа от грани і + 1/2, ЛКое — матрица, элементы которой вычислены по специальной процедуре осреднения параметров потока на грани ячейки [10].

I; К

Векторы Єг+1/2 и Єг+1/2 вычисляются с использованием вектора примитивных переменных

Т

q = (р. и. V. и1, р) , величины которого слева и справа от грани ячейки находятся по формулам [11]:

їч-і/2 = <7, +1[(! - ф)^г + (1 + ф)Аг ], д*+1,2 = Я1+1 - 7 [(! + ф)^г+1 + (1 - ф) Аг+1 ],

где Уг =тіпто<і (Vq,bAq)i

тіп

Аг = тіп тос1(Д<7.ЪЩ\, Ас], = qi+l Vql = ql - ql , ,

3 -ф

іп тос1(х, у ) - 8І§п(х) тах[0, тіп {|х|, у 8І§п(х)} ], Ь =---, -1 < ф < 1.

1 1 1-ф

Ф

Параметр ф определяет порядок точности разностной схемы, например, при ф = -1 получаем второй порядок точности. В настоящей работе бралось значение ф = 1/3, при котором конвективные члены в уравнениях движения аппроксимировались с третьим порядком точности на гладких решениях. Для аппроксимации вязких членов в уравнениях (5) применялись центральные разности второго порядка точности.

Система линейных разностных уравнений (5) относительно приращений А£?/ ■ к решается

методом релаксаций Гаусса — Зейделя по плоскостям, соответствующим индексу /. Этот метод позволяет перейти от решения трехмерной системы разностных уравнений к решению двумерной системы [12]. В каждой поперечной плоскости / решение системы двумерных разностных урав-

2

нений производится методом приближенной факторизации с усиленным диагональным преобладанием:

(А + )А(А + )А£?* — Щ1,

_ V ~ V

А — —I + А = —I + + £"2^

А? ^ Л? 5 л ^

Ц= =0.5(1А11+112 + 1А11_112), Ц =0.5(/В/М12 + /В/^12), Ц. =0.5(/С7*+1/2 + /С7*_1/2),

3^6 <3(§- 0/$-

где А = В = ^ = — матрицы Якоби невязких потоков; — соответствующие

разностные операторы вдоль направлений г| и С,-

Аналогичный подход использовался при численном решении уравнений модели турбулентности (3). Только в этом случае для аппроксимации конвективных потоков кинетической энергии турбулентности и диссипации через грани расчетных ячеек применялись односторонние разности против потока первого порядка точности. Кроме того, при решении разностных уравнений методом релаксаций Гаусса — Зейделя выполняется только прямой ход в направлении Решение разностных уравнений модели турбулентности было сдвинуто на один шаг по времени вперед относительно решения разностных уравнений движения.

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

Начальное распределение параметров течения в сопле задавалось по одномерной теории в невязком ядре и по закону степени 1/7 в пограничном слое на стенках сопла.

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

4. АНАЛИЗ РЕЗУЛЬТАТОВ РАСЧЕТОВ

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

На рис. 3 дано сравнение результатов настоящего расчета с результатами экспериментального исследования [13] поля течения в трансзвуковой области осесимметричного сопла Лаваля с углами наклона образующих в дозвуковой и сверхзвуковой частях, равными соответственно 45° и 15°, и относительным радиусом кривизны контура в критическом сечении Н/г,, = 0.625 (г* — радиус критического сечения). Хотя рассматриваемое сопло является осесимметричным, расчеты были проведены на трехмерной сетке с целью проверки численного алгоритма на предмет сохранения окружной симметрии потока. Сетка содержала I х J х К = 75 х 45 х 36 ячеек (I, J, К — количество ячеек соответственно в осевом, радиальном и окружном направлениях). Анализ показал, что численное решение очень хорошо сохраняет симметрию потока. Например, статическое давление на поверхности сопла в критическом сечении остается постоянным в азимутальном направлении с точностью до четырех знаков, за исключением одной точки при к = 1.

Р%с

б)

0.8

0.6

0.4

0.2

стоящий расчет

эксперимент [13]: (2) ось симметрии —

V А ст енк а

А \

\ 4 \

ч к

кА к 1

0

-0.5

0.5

1

х-х,

1,5

Рис. 3. Сравнение расчетных и измеренных распределений числа Маха (а) и относительного давления (б) в осесимметричном сопле

Рассчитанные (сплошные линии) и измеренные (значки) распределения уровней числа Маха в меридиональной плоскости сопла сравниваются на рис. 3, а. В окрестности минимального сечения сопла наблюдается хорошее согласование расчета и эксперимента. Некоторое расхождение результатов, заметное в сверхзвуковой части сопла, связано, по-видимому, с разным нарастанием пограничного слоя в расчете и эксперименте. На рис. 3, б сравниваются рассчитанные и измеренные распределения относительного статического давления р/рос (р0с — полное давление в сопле) на стенке и на оси симметрии сопла. Видно, что результаты численного расчета хорошо согласуются с экспериментальными данными.

Для проверки влияния размеров ячеек сетки на численное решение были выполнены расчеты на трех сетках с разным количеством ячеек: I х 3 х К = 52 х 30 х 24, 75 х 45 х 36 и 102 х 60 х 48. Они показали, что зависимости распределения числа Маха, а также распределения статического давления на стенке и оси сопла, полученные на разных сетках, практически сливаются на графиках. Рассчитанные на этих сетках значения коэффициента расхода сопла составили ц = 0.9807,

0.9803 и 0.9799 соответственно, что хорошо согласуется с измеренной величиной ц = 0.983 ±0.006.

На рис. 4, 5 приведено сравнение результатов расчета настоящим методом с расчетными и экспериментальными данными работы [14] для сверхзвукового сопла прямоугольного сечения. Рассмотренное сопло имело параллельные боковые стенки, поэтому ширина проточной части сопла была постоянной (рис. 4, а). Верхняя и нижняя образующие располагались симметрично относительно горизонтальной плоскости. Контур сверхзвуковой части сопла представлял собой прямолинейный отрезок, расположенный под углом 2.56° к оси сопла и плавно сопряженный с дугой окружности относительным радиусом Я/А* = 0.24 (А* — высота критического сечения, равная 2.08 см; все размеры на рис. 4, а даны в сантиметрах).

Расчет параметров потока в сопле проводился в области, ограниченной вертикальной плоскостью симметрии и стенками сопла, на разностной сетке с числом ячеек I х 3 х К = 105 х 52 х 48. На рис. 4, б сравниваются рассчитанные в настоящей работе и в работе [14] распределения числа

б)

Рис. 4. Уровни числа Маха в вертикальной плоскости симметрии трансзвуковой области прямоугольного сопла

Маха в вертикальной плоскости симметрии трансзвуковой области сопла (ввиду симметрии течения относительно горизонтальной плоскости результаты нанесены для верхней половины сопла). Видно, что оба расчета хорошо согласуются друг с другом. На рис. 4, б можно заметить, что вниз по течению от критического сечения на стенке сопла формируется волна сжатия, которая затем отражается от горизонтальной плоскости симметрии и образует (вместе с симметричной волной на нижней стенке) волновую структуру течения в сверхзвуковой части сопла. На существование такой структуры указывает, в частности, волнообразная форма приведенных на рис. 5 распределений статического давления на поверхности и на оси симметрии сопла. Вследствие неоднородности сверхзвукового потока в окрестности критического сечения сопла, когда скорости вблизи стенки выше, чем на оси, возникает волна сжатия, которая затем пересекается на горизонтальной плоскости симметрии с такой же волной сжатия, возникшей на противоположной стенке. В результате последовательного отражения этих волн от стенок сопла образуется периодическая ударно-волновая структура. На рис. 5 видно, что максимумы давления на верхней стенке в плоскости симметрии (ряд 1) и около боковой стенки (ряд 2) располагаются примерно в сечениях, где образуются минимумы давления на оси симметрии (ряд 0) и на боковой стенке в плоскости симметрии (ряд 3). При приближении к срезу сопла интенсивность волн сжатия понижается, и волновая структура практически исчезает.

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

Ниже представлены результаты численного расчета пространственного течения газа в сверхзвуковых соплах с круглыми входным и критическим сечениями и с треугольным выходным сечением. Расчеты выполнены для нескольких вариантов сопел, имеющих одинаковую геометрию дозвуковой части и разные длины сверхзвуковой части при одном и том же значении относительной площади среза =7^//^ = 2.05 (Т7* — площадь критического сечения). В качестве

Рис.

5. Распределение давления по поверхности прямоугольного сопла

примера на рис. 6 показана геометрия сопла с относительной длиной сверхзвуковой части

1.с = /.с= 2.44. Все размеры на этом рисунке отнесены к радиусу критического сечения сопла. Отношение площади входного сечения сопла к площади критического сечения составляло /'|)х= 4.3, угол наклона образующей дозвуковой части сопла равен 20.31°. Выходное сечение имело форму равностороннего треугольника и располагалось таким образом, чтобы ось симметрии дозвуковой части сопла проходила через центр треугольника (точку пересечения медиан). Отношение длины стороны треугольника к радиусу критического сечения сопла равно 3.856. Основные геометрические параметры сопла, изображенного на рис. 6, взяты из работы [15], в которой аналогичное сопло (только со скругленными углами выходного сечения) исследовалось экспериментально. Кроме того, были проведены расчеты течения в сверхзвуковом сопле треугольного сечения с таким же распределением площади по длине, что и у сопла с круглым критическим и треугольным выходным сечениями.

Расчеты течения в сопле проводились при перепаде давлений в сопле пс = рос/р® = 11.06 (рда — давление в окружающем пространстве), температуре торможения Тос = 300 К и показателе адиабаты у = 1.4. Число Рейнольдса, вычисленное по параметрам торможения и радиусу критического сечения сопла, бралось равным 106. На рис. 7 показана разностная сетка в плоскости симметрии (а) и в поперечном сечении (б), которая использовалась при численном решении. Сетка содержала 96 ячеек в продольном направлении, из которых 40 приходилось на сверхзвуковую часть сопла, 30 — в радиальном и 45 — в окружном направлениях.

Д.-4 хс~6А4

Рис. 6. Г еометрия сопла с круглым критическим и треугольным выходным сечениями

Рис. 8. Уровни числа Маха в плоскости симметрии сопла с треугольным выходным сечением (а) и в эквивалентном осесимметричном сопле (б)

Расчеты показали, что поле течения в сверхзвуковой части сопла с треугольным выходным и круглым критическим сечениями имеет сложную структуру, обусловленную продольной и поперечной неоднородностью газодинамических параметров. На рис. 8, а представлено распределение уровней числа М в вертикальной плоскости симметрии сопла с треугольным выходом. Темные области на этом рисунке относятся к пограничному слою на поверхности сопла. Из-за сильного разгона потока в окрестности угловой точки на верхнем контуре сопла формируется локальная пространственная волна сжатия. Следует отметить, что такая же волна сжатия формируется около другой образующей, которая проходит через нижнюю вершину треугольника. В дозвуковой части сопла, вплоть до критического сечения, течение является осесимметричным, и распределение параметров в нем имеет такой же характер, как и в большинстве осесимметричных сопел с конической дозвуковой частью при умеренных углах сужения.

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

Рис. 9. Уровни числа Маха в поперечных сечениях сопла с круглым критическим и треугольным выходными сечениями

сопел одинаковая, то и поля течения в них практически одинаковые. Распределения уровней числа М в сверхзвуковых частях существенно различаются.

Изменение формы поперечных сечений по длине сверхзвуковой части сопла приводит к тому, что в каждом сечении формируется различная структура неоднородности газодинамических параметров. Характерное распределение уровней числа М в поперечных сечениях х = 0, 0.5, 1 и 2.4 показано на рис. 9. В сечении, расположенном на расстоянии примерно х = 0.5 от критического сечения, около вершин треугольника зарождаются конусоподобные волны сжатия. На поверхности сопла эти волны формируют весьма неоднородное распределение статического давления (рис. 10, а). В случае треугольного сопла, когда все сечения имеют треугольную форму, качественная картина распределения статического давления на поверхности примерно такая же, но без четко выраженных волн сжатия (рис. 10, б). Следует отметить, что приведенные на рис. 10 распределения давления существенно отличаются от случая осесимметричного сопла, в котором распределение давления на поверхности представляет собой семейство окружностей.

На рис. 11 приведены распределения относительного статического давления р/р0с вдоль образую-

настоящий расчет:

- круглое критическое сечение _ осесимметричное

' сопло (эквивалентное)

- треугольное критическое сечение

эксперимент [15): круглое критическое сечение осесимметричное сопло (эквивалентное)

0.4

рр*

0.3

0.1

/ ч N 'Л

]° .а/ /

\у’ 1 у 1 А Д 1 \ А \ , — ■

N В

0.5

1.5

2 2.5

(х-х,) /;

Рис. 10. Изобары на поверхности сопла с круглым (а) и треугольным (б) критическими сечениями

Рис. 11. Распределение давления в соплах с треугольным выходным сечением и в эквивалентном осесимметричном сопле

щих сверхзвуковой части сопла, проходящих посередине стороны (линия А) и через вершину (линия В) треугольника в выходном сечении. Сплошные кривые относятся к соплу с круглым критическим и треугольным выходным сечениями, штриховые кривые — к соплу, у которого критическое и выходное сечения имеют форму треугольника. В обоих случаях после понижения давления за угловой точкой на стенке сопла происходит его повышение до некоторого максимума, после чего давление снова понижается. При этом на образующей А максимум давления в сопле с круглым критическим сечением составляет р/рос = 0.34 и расположен примерно на расстоянии х = 0.4 от критического сечения, в то время как в сопле с треугольным критическим сечением максимум давления равен р/рос = 0.4 и расположен на том же расстоянии х = 0.4. На образующей В максимум давления в сопле с круглым критическим сечением достигает величины р/р0с = 0.195 в точке х = 0.9, а в сопле с треугольным критическим сечением р/р0с = 0.15 в точке х = 1.8. По мере приближения к срезу сопла давление на образующей А становится ниже, чем на образующей В, однако разница между ними, после некоторого возрастания на небольшом участке, становится незначительной на срезе (< 0.02 р0с).

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

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

Для практики представляет интерес сопоставление расходных и тяговых характеристик пространственных сопел и эквивалентных осесимметричных сопел с таким же распределением площади поперечных сечений. В результате численных расчетов были получены значения коэффициента расхода ц = ('/С/('7ИЛ. потерь импульса А/с =1-/с//р и потерь тяги АРС = 1 - Рс /Рс ид.

где Ос — действительный расход через сопло; Gид — расход, вычисленный по одномерной теории; 1с — импульс сопла; 1р — расчетный импульс, вычисленный по одномерной теории; Рс — тяга сопла; Рс ид — идеальная тяга при полном изоэнтропическом расширении струи.

В таблице приведены значения ц. А/с и АРС для сопел с треугольным выходным сечением и эквивалентного осесимметричного сопла с относительной площадью среза Рс = 2.05 и относительной длиной сверхзвуковой части Ьс = 2.44. Из расчетов следует, что коэффициент расхода принимает одинаковое значение, равное 0.967, как для сопел с треугольным выходом, имеющих круглое или треугольное критическое сечение, так и для эквивалентного сопла. Расчеты показали, что, несмотря на существенное различие в распределениях давления по поверхности сверхзвуковой части пространственных сопел с круглым и треугольным критическим сечениями

Расчет Эксперимент [15]

Характе- ристики сопла Тип сопла

Круглое критическое сечение треугольный выход Треугольное критическое сечение треугольный выход Круглое критическое сечение — круглый выход Круглое критическое сечение треугольный выход Круглое критическое сечение — круглый выход

ц 0.967 0.967 0.967 0.968 0.968

0.0213 0.0196 0.0114 0.0147 0.0091

ЛРс 0.0247 0.0227 0.0133 0.0191 0.0111

_ 0.06

Л/с.

0.05

0.04

0.03

0.02

0.01 0

1,5 2 2.5 3 3.5 _ 4

Рис. 12. Зависимости потерь импульса и тяги от длины сопла:

1 — круглое критическое и треугольное выходное сечения; 2 — эквивалентное осесимметричное сопло

(см. рис. 10), тяговые характеристики этих сопел отличаются незначительно. Например, по потерям тяги это различие составляет 0.002. Следует отметить, что тяга рассмотренных сопел направлена практически вдоль центральной оси. Угол отклонения вектора тяги не превышает 0.14° у сопла с круглым сечением и 0.02° у сопла с треугольным сечением. Потери тяги сопел с треугольным выходным сечением с РС = 2.05 и ЬС = 2.44 выше, чем у эквивалентного осесимметричного сопла, примерно на 1% идеальной тяги.

Для сравнения в таблице представлены результаты измерений [15]. Видно, что результаты настоящих расчетов хорошо согласуются с экспериментальными данными. Расхождение значений потерь импульса относительно расчетного составляет примерно 0.2% для круглого сопла и 0.7% для сопла с треугольным выходом. Более высокая разница в потерях импульса трехмерного сопла обусловлена, по-видимому, некоторым отличием геометрии сопел в расчетах и в эксперименте (при расчетах ребра сопла были острые, а в эксперименте скругленные). Коэффициенты расхода равны 0.968 в эксперименте и 0.967 в расчетах.

С помощью разработанного метода была проведена серия расчетов пространственных сопел с треугольным выходным и круглым критическим сечениями при одинаковой относительной площади среза /^ = 2.05, но с различными длинами сверхзвуковой части в диапазоне /,с = 1.6 -ь 3.8. Для каждого варианта пространственного сопла был выполнен расчет соответствующего осесимметричного сопла с таким же распределением площади поперечных сечений. Результаты расчетов представлены на рис. 12, где показаны зависимости потерь тяги АРс и потерь импульса А/с от относительной длины сопла /с. При увеличении /,с потери тяги и импульса сопла уменьшаются, причем наиболее заметное влияние длины сопла наблюдается при изменении ЬС от 1.6 до 2.5. Например, на этом отрезке значений /,с величина Д/с трехмерного сопла изменяется примерно на 2.7% от расчетного импульса (с 0.048 до 0.021), в то время как при увеличении

Гс _

от 2.5 до 3.8 величина Д/с изменяется всего на 0.4% (с 0.021 до 0.017). Примерно такой же характер наблюдается при изменении потери тяги АРс.

На рис. 12 показано также сравнение тяговых характеристик трехмерных сопел с треугольным выходом и эквивалентных круглых сопел. Наибольшая разница в их характеристиках наблюдается при меньших значениях ЬС, и по мере увеличения длины сопла характеристики трехмерных сопел приближаются к характеристикам круглых сопел. В частности, при ЬС = 1.6 потери импульса сопла с треугольной формой выходного сечения выше, чем у эквивалентного осесим-

метричного сопла, примерно на 2.3% от расчетного импульса. При Zc> 2.5 величина Д/с сопел с треугольным выходом примерно на 0.5 — 1% выше, чем у осесимметричных сопел. Приведенные на рис. 12 данные показывают, что влияние пространственности течения на потери импульса и тяги сопла наиболее заметно проявляется при небольших длинах (Lc < 2) сверхзвуковой части и значительно в меньшей степени при больших длинах (Lc > 3).

ЗАКЛЮЧЕНИЕ

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

Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект 10-01 -00208-а).

ЛИТЕРАТУРА

1. Пиру мов У. Г., Росляков Г. С. Газовая динамика сопел. — М.: Наука, 1990,

3б8 с.

2. Л а в р у х и н Г. Н. Газодинамика реактивных сопел. Т. 1. Внутренние характеристики сопел. — М.: Физматлит, 2002, 376 с.

3. Mazurov A. P. Computation of viscous gas flow field in a lobed nozzle // XlV 1SABE. —

Florence, ltaly, 1999.

4. Иванов М. Я., Крайко А. Н. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. 11 // ЖВМ и МФ. 1972. № 3, с. 805 — S13.

5. Иванов М. Я., Рылько О. А. К анализу трансзвукового течения в эллиптических соплах // Изв. АН СССР. МЖГ. 1972. № 3, с. 161 — 1б3.

6. Пирумов У. Г. Пространственные до- и сверхзвуковые течения в соплах и каналах переменного сечения // ПММ. МЖГ. 1972. Т. 36, № 2, с. 239 — 247.

7. Д в о р е ц к и й В. М., И в а н о в М. Я. К расчету смешанных течений в соплах с несимметричной дозвуковой частью // Ученые записки ЦАГИ. 1974. Т. 5, № 5, с. 39 — 45.

S. Chien K.-Y. Prediction of channel and boundary-layer flows with low-Reynolds-number turbulence model // AlAA J. 19S2. V. 20, N 1, p. 33 — 3S.

9. Годунов С. К., Забродин А. В., Иванов М. Я., Крайко А. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики / Под ред. С. К. Годунова. — М.: Наука, 197б, 400 с.

10. R o e P. Approximate Riemann solvers, parameter vectors and difference schemes //

J. Comput. Phys. 19S1. V. 43, p. 357 — 372.

11. Obayashi S. Numerical simulation of underexpanded plumes using upwind algorithms // A1AA Paper SS-4360 CP. 19SS.

12. Чакравати С. Р., Жем К.-Й. Расчет трехмерных сверхзвуковых течений с дозвуковыми зонами на основе уравнений Эйлера // Аэрокосмическая техника. 1987. № 11, с. 22 — 35.

13. Cuffel R. F., Back L. H., Massier P. F. Transonic flowfield in a supersonic nozzle with small throat radius of curvature // A1AA J. 19б9. V.7, N 7, p. 13б4 — 13бб.

14. Compton W. B., 111. Comparison of turbulence models for шжк-айегЬ^у flows with propulsive jets // NASA TP-3592. 199б.

15. Курочкин Н. А., Лаврухин Г. Н. Характеристики сверхзвуковых сопл с различными выходными сечениями // Труды ЦАГИ. 1974, вып. 7074.

Рукопись поступила 21/IX 2010 г.