ПРИМЕНЕНИЕ 3D КОДА РАСЧЕТА ТУРБУЛЕНТНОГО ТЕЧЕНИЯ К АНАЛИЗУ "ГИСТЕРЕЗИСА" ДАВЛЕНИЯ В ВАКУУМНОЙ КАМЕРЕ ГАЗОВОГО ЭЖЕКТОРА Текст научной статьи по специальности «Физика»

Труды МАИ. Выпуск № 103

http://trudymai.ru/

УДК 621.431.75

Применение 3D кода расчета турбулентного течения к анализу «гистерезиса» давления в вакуумной камере газового эжектора

Ларина Е.В.1*, Ципенко А.В.1**, Афанасьев А.А.1***, Крюков И.А.2****

1 Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993, Россия Институт проблем механики им. А.Ю. Ишлинского Российской академии наук, проспект Вернадского, 101, корп.1, Москва, 119526, Россия *e-mail: larinaelenav@gmail.com

**e-mail: tsipenko_av@mail.ru ***e-mail: tempero@tempero.com ****e-mail: kryukov@ipmnet.ru

Аннотация

В статье представлены результаты анализа экспериментальных данных о

работе газового эжектора как вакуумного насоса и результаты численных

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

трехмерного течения по неявной схеме по времени с третьим порядком точности по

пространству. Код тестирован на течениях с однородным сдвигом, турбулентным

пограничным слоем вдоль плоской пластине и автомодельной струе. Выполнено

сравнение с экспериментом и результатами расчетов сверхзвуковой

перерасширенной струи по стандартной k-s модели, k-s модели с поправкой Chen,

k-s модели RNG. Получено качественное соответствие эксперименту по ударно-

волновой картине течения и приемлемое количественное согласие по давлению

Пито. Численным экспериментом показано, что петля «гистерезиса» в зависимости

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

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

Введение

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

а) б)

в) г)

Рис. 1. Зависимость давления в вакуумной камере от давления на входе в сопло для

различных опытов

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

другая модель, учитывающая временную задержку в изменении напряжений Рейнольдса при внезапном изменении скоростей деформации [6].

Анализ экспериментальных данных

Рассматривается одноступенчатый осесимметричный эжектор. Рабочий газ -азот. Параметры канала: выходное отверстие осесимметричного сопла рабочего газа соответствует числу Маха 3.31; около выходного сечения из сопла среднее число Маха 4.035, а среднее давление 0,063 бар. Подробные данные и эскиз эжектора приведены в работах [2, 4]. Особенностью переходного режима потока являются низкочастотные пульсации давления, связанные с периодическим смещением точки отрыва псевдо-скачка. Анализ экспериментальных кривых (рис. 2) показывает, что расход потока изменяется в переходном режиме. Колебания давления в вакуумной камере и на выходе сопла аналогичны колебаниям осциллятора переменной (убывающей) массы с затуханием (для сравнения рис.2 и рис.3).

Рис.2. Вертикальные шкалы в psi (1 psi = 0.06804596390998 физическая атмосфера). Графику давления на входе в сопло соответствует правая шкала, остальным - левая.

Рис. 3. Амплитуда колебаний осциллятора переменной (убывающей) массы с

затуханием как функция времени (шкалы условные, время - по горизонтальной оси).

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

сохраняться при меньшем давлении в предкамере сопла. Это явление представляет собой «гистерезис» давления эжектора (рис. 1-б, в, г).

Это явление даже на одном и том же устройстве наблюдается не всегда (сравни рис. 1-а и 1-г). Явный «гистерезис» наблюдался при малом времени выхода эжектора на минимальное давление. Результаты, показанные на рис. 1-а,б,в,г, сведены в Таблицу. Видно, что большая петля «гистерезиса» (рис.1-б,в,г) связана с коротким временем запуска, то есть при увеличении времени запуска и плавном повышении давлении петля должна уменьшиться, что подтверждает опыт (см. рис. 1-а). Аналогичный «гистерезис» отмечен в работе [1], но без низкочастотных пульсаций, возможно, это связано с осреднением экспериментальных данных для однотипных запусков.

Таблица

Связь «гистерезиса» со временем запуска и остановки эжектора

рис. давление на входе в сопло растет давление на входе в сопло падает примечания

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

1-б 12 4 5 3 нет «гистерезис»

1-г 3 - нет 60 18 «истинный гистерезис»

1-в 8 3 нет

20 8 3 30 10 нет «гистерезиса»

Небольшое отличие в давлениях в вакуумной камере на рис. 4-г при запуске и остановке эжектора в области стационарного течения при большом давлении эжектирующего газа иллюстрирует «истинный гистерезис» (на рис. 1г - hysteresis), но это единичный опыт, который требует повторения и более детального изучения.

Тем не менее, на основе имеющихся данных можно сделать вывод, что «гистерезис» при запуске-остановке эжектора связан с нестационарным изменением параметров в процессе запуска-остановки эжектора. То есть при стационарных граничных условиях параметры в канале эжектора и вакуумной камере будут определяться однозначно. Поэтому простые одномерные модели (например, [7,8,9]) будут давать, при соответствующей настройке, вполне приемлемые результаты. Также 2D и 3D численное моделирование можно проводить методом установления с использованием локального шага по времени. Однако при эксплуатации реальных установок время выхода на оптимальный режим при плавном повышении давления заметно больше, чем при сначала резком выходе на завышенные параметры и затем небольшом плавном дросселировании давления рабочего газа. Именно такой процесс отражает петля «гистерезиса».

Этот вывод проверялся в численном эксперименте, о котором рассказано

ниже.

Математическая модель и численный метод

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

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

дt(pk) + V ■ (puk) = V ■ (ß + ß /ак)Vk)+?k -ps(1 + M2t) 51(ps)+V ■ (piie) = V ■ (ß + ßJaE)Vs)+CeiPk elk - CE2ps2/k, дtр)+ V\puVt) = CrPs(ytE -ytVk>

Здесь и - скорость газа, p - плотность, - порождение кинетической энергии

турбулентности, заданное формулой

Г Г 2 ^ 2 ^

Pk = \д и+ди- - stj д тит у- PkSy Jd и,

турбулентное число Маха Mt =422к/а , где a - скорость звука, равновесная турбулентная вязкость vtE = C pkVs и константы: Csi = 1.44,Cs2 = 1.92 ,ok = 1 ,*e = 13,Q = 0.75 .

Для численного решения используется метод Годунова и неявная схема Гаусса-Зейделя (LU-SGS scheme) [10] по времени. Для нахождения потоков на грани ячейки сетки используется точный решатель задачи о распаде разрыва.

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

f (x) = /о + Z a (Х - x)+1Z bj (Х - xXxj - x°),

i=1 2 i,j=1

5f ( ü\ -7 52f

где a = Iх ) - первые производные в центре ячеики, b = —

5х J йхх.

(х0) - вторые

производные в центре ячейки, x0 = (x( ).=1 3 - геометрический центр ячейки, x = (xt ).=1 3 - произвольная координата внутри ячейки, функция f (x) обозначает одну из неизвестных величин (плотность, компонента скорости, давление, кинетическая энергия турбулентности, скорость диссипации, турбулентная вязкость).

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

~ пЪ ( пЪ (Л/ ( пЪ (Л

~ = (x, - x Imax (x - x ),

где xnb = (xnb

(xf) - геометрический центр какоИ-либо соседней ячеИки. Соседними

считаются все ячейки, имеющие с данной общую грань или общую вершину. Рассматриваем среднее квадратичное отклонение Ф

3 13 \ С 3 13

ф=Е f - fnb+£ a, ~nb+- х ь, = x f - ft+£ a a"b a+1£ ь„ a

nb

2 Л

nb

2 tit-

nba nb

x

U i J

nb,i

где ~ = а/Д - новая неизвестная, Д = шах*/16 - максимальный модуль разности

координат.

Ищем критическую точку функционала Ф ^ Arg min . Тогда имеем систему

линейных алгебраических уравнений

1 5Ф

2

С 3 1 3

= Е anb f - fn+У a anbA+-£ b

/ 1 i J ü J ü / 1 mm ~ / 1 i

nb у m=1 2 m,p=1

3

1 _5Ф

2 Ж"

anbanb

mp m p

nb

m, p=1 3

С 3 3

= £ anbanb f - fnb+y a anbA+1 £ b

/ 1 i J Jü Jü / i mm ~ / 1 I

\

1 5Ф 1

ü ü ' / ^"m"m у m=1

3

anbanb

mp m p 2 m,p=1

=1 y(xnbj2 f - fnb + Ya a "bA +1 £b

О Я/i О ' ' i ' ' mm ~ / . I

2 5bii 2 nb у m=1 2

anbanb

mp m p

' m,p=1

i = 1,...,3

i * J

i = J

2

i=1

i=1

:

Полученная система решается методом Гаусса с выбором ведущего элемента.

Течение с однородным сдвигом.

Задача о течении с однородным сдвигом рассматривается в упрощенной одномерной постановке, и для тестирования алгоритма рассчитывается численно. Параметры моделирования соответствуют эксперименту [11]. Начальные условия для упрощенной системы и входные граничные условия соответствуют экспериментальным значениям в точке x/h=7.5. Среднее течение задается скоростью и = {и(у) 0 о), где скорость линейно растет вдоль координаты у с наклоном

dU ( y) dy

= 46.8s 1. Давление в области задавалось постоянным и соответствующим

нормальным условиям. Градиент статической температуры dT (y) = 9.5 Ks 1 .

dy

9 9 ^ 9 ^

Начальные значения турбулентных величин k=0.268m /s , s=1.94-10 m /s , ^t/^L=107.

В случае рассматриваемого течения взяты средние параметры течения вдоль

центральной линии. В этом случае параметры турбулентных величин можно найти

из упрощенной системы уравнений. В случае k-s-^t модели система следующая

Dk/Dt = vtS2 -s,

Ds/Dt = CslvtS2 s/k - Ce2 s2/k,

DvjDt = Qs (С, k V s-v t )/ k,

для lag модели [12] имеем:

Dk/Dt = vtS2 - Ceka,

Da IDt = Cv 2alk - C«2a: DvjDt = CTtt)a(CM k/ a-Vt)

C = 0.09, С = 1.44, С. = 1.92, С = 0.35

s ? al ? a 2 ? та

<

В случае k-ю и lag моделей значение ю на входе задается по формуле с = Cste.

Полученные результаты одномерного приближения и численного моделирования представлены на рис. 4-6. Видно, что трехмерный код (3D case) дает те же результаты, что получаются по упрощенной модели. Аналогичная проверка в случае других моделей также дает идеальное соответствие результатов (на рисунках не показано). Наличие граничного условия для турбулентной вязкости позволяет использовать все исходные средние параметры турбулентности, заданные в эксперименте. В случае стандартной k-s и k-ю моделей без изменения используемых констант можно задать только два параметра турбулентности на левой границе области. Градиенты турбулентных величин в k-s-^t модели лучше других моделей согласуется с экспериментальными значениями. Lag-модель дает заметно завышенные значения турбулентных величин. Градиент турбулентной вязкости, описываемый моделью vt-92 [13] , очень точно согласуется с экспериментом. Чуть менее точный градиент вязкости получен и для k-ю модели, хотя для нее имеются заметные расхождения с экспериментом по остальным турбулентным величинам.

Рис. 4: Профили кинетической энергии турбулентности в течении с однородным

сдвигом

Рис. 5: Скорость диссипации кинетической энергии турбулентности в течении с

однородным сдвигом

3

А

107р

2,5 2

1,5 1

~~89 X/h 10 11

Рис. 6: Профили турбулентной вязкости в течении с однородным сдвигом.

Штриховая линия соответствует стандартной k-s модели, сплошная линия: k-s-штрих-пунктир: lag модели [12], пунктир: k-ю модели Wilcox[14]. Символы: круги -эксперимент [11], квадраты -расчет с k-s- моделью, треугольники - расчет с

моделью vt-92 [13]. X0 / h = 7.5, h = 0.31 м.

Турбулентный пограничный слой на плоской пластине

Проведено тестирование программного кода в случае течения в пограничном слое на плоской пластине. Параметры течения выбраны следующие: M=0.5, длина пластины L=10м, T=300K, стенка изотермическая с температурным фактором 1, ламинарная вязкость определена по Сазерленду. На рис. 7 изображены коэффициент трения (рис. 7A) и скорость (рис. 7Б) в логарифмических координатах в сравнении с известными аналитическими решениями для данного случая. При моделировании сетка выбиралась таким образом, чтобы первая ячейка попадала в логарифмический подслой для высокорейнольдсовой модели. Из рисунков видно, что соответствие

расчета с известными решениями хорошее. Расчет проводился на последовательности из трех сеток 50x70x10, 100x140x20, 200x280x40, различие результатов по модулю скорости и коэффициенту трения на стенке для этих сеток составил менее 1%.

Рис. 7: Продольная скорость в безразмерных координатах (А) и коэффициент трения (Б) в течении на плоской пластине. Символы - расчет, кривые - аналитические

кривые.

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

Рассматривалось автомодельное течение в плоской струе. Жидкость несжимаемая. Вводим функцию тока р ( и = дур , V = -дхр ), и автомодельные

координаты £ = х ,] = у/£(х), где 5{х) - толщина пограничного слоя. Предполагаем, что все функции в системе в безразмерном виде зависят только от поперечной координаты г. Тогда

р = *(Йи(Й®(?) ,г12 = т0(£)г]г), к = к №(?) ,£ = ^Мг) = ^Мг),

где и (£) - значение скорости на оси струи, размерные функции

г0 £) = и2 {£), к, (£) = и2, ^ (£) = и3 , ^ £) = и(^).

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

Д(ф ' )2 -ДФФ ' + Т' = 0,

Т + ДФТ/Ы + Д(Ф ')2 = 0,

2Д Ф Т - ДФР ' - (р N + Р N + Е + ТФ ' = 0,

(4Д - Д )Ф' Е - ДФЕ ' - (ЕN ' + ЕN)/ае + Се1 ТЕФ '/Р - Се2 Е2/Т = 0,

Д (#Ф ' - ') - С Е/Е (см Р2/Е - N)= 0.

где параметры Д = би/и , Д=б{би/5 + и)/и , являются постоянными. Учитывая сохранения импульса по сечению струи, получаем Д = -Д.

Для полученной системы ОДУ построена конечно-разностная схема. Решение ищется итерационно и на каждом шаге корректируется исходное значение параметров. В результате (рис.8) имеем хорошее согласие с экспериментальными данными Фертмана [15].

0,4 0,8 У/У с 1,2 1,6

Рис. 8: Продольная скорость в безразмерных координатах. Ш - скорость на оси струи, Уе - полутолщина струи. Сравнение моделирования (кривые) с

экспериментом Фертмана [15].

<

Сверхзвуковая перерасширенная струя

Проведено моделирование перерасширенной струи воздуха, истекающей в затопленное пространство с расчетным числом Маха 3.005. Данное течение соответствует экспериментальной работе [16]. Параметры течения следующие: температура торможения Т0=287К в струе, статическое давление и температура окружающего воздуха РИ=1.008 бар, ТИ=294К, отношение давления торможения в струе к давлению в окружающем пространстве Po/Ph=21.8. Использовалась последовательность неструктурированных сеток (3 млн., 7 млн., 10 млн. ячеек), содержавших шестигранные и призматические ячейки. Расчетная сетка имела плавные сгущения к стенкам сопла.

Рис. 9: Поле числа Маха перерасширенной осесимметричной сверхзвуковой струи

Mach Number

Яз.; з.:

3.80 3.23 2.74 2.33 1.98 1.68 1.42 1.21 1.03 0.87 0.74 0.63 0.53 0.45 0.38 0.33 0.28 0.24 0.20

Z

О 20 Х/Х* 40

Рис. 10: Давление Пито на оси перерасширенной осесимметричной сверхзвуковой

струи

Поле числа Маха, полученное при моделировании, представлено на рис. 9. Давление Пито на оси струи, полученное в расчетах с использованием четырех моделей турбулентности, сопоставлено на рис. 10 с экспериментальными данными. Падение давления Пито, начиная от среза сопла до отражения падающего скачка уплотнения от оси, все модели предсказывают одинаково, так как влияние вязкости в данной области пренебрежимо мало. В расчете по «стандартной» к-е модели не образовался диск Маха, поэтому течение осталось сверхзвуковым и давление снизилось незначительно. Давление Пито за отражением косой волны возрастает резко потому, что при больших числах Маха небольшое изменение числа Маха заметно отражается на изменении давления Пито. В расчете по к-е-^ модели давление Пито после взаимодействия хорошо соотвествует эксперименту. Снижение давления Пито связано с влиянием волны разрежения, а последующее увеличение с влиянием волн

сжатия. На рис. 10 представлены результаты, полученные по моделям Chen и RNG. Здесь под моделью Chen понимается модель, содержащая поправку на неравновесность турбулентности из работы [17] с моделью сжимаемой диссипации [18]. под RNG моделью понимается поправка на неравновесность [19] с той же моделью сжимаемой диссипации [18]. В k-s моделях Chen и RNG величина диска Маха немного завышена и течение за областью отражения дольше остается дозвуковым, поэтому давление Пито занижено по сравнению с экспериментальными данными. Течение внутри сопла моделируется с использованием пристеночных функций с учетом градиента давления, одинаковых для всех четырех моделей турбулентности. В дальнейших бочкообразных структурах k-B-^t модель начинает завышать давление Пито, хотя соответствие эксперименту остается приемлемым. «Стандартная» модель завышает давление Пито, а затем внутренняя граница внешнего слоя смешения попадает на ось раньше, чем в эксперименте. В результатах по модели Chen смыкание внешнего слоя смешения происходит дальше по потоку по сравнению с экспериментом, а по модели RNG немного ближе по потоку.

Анализ результатов численного моделирования течения в эжекторе

Поток в эжекторе от старта до установления стационарного течения

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

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

основании приведенных выше результатов применялась k-B-^t модель

турбулентности. Численное моделирование показало, что как при повышении, так и

при понижении давления рабочего газа в эжекторе реализуются одинаковые стационарные или пульсационные режимы течения (рис. 11).

Рис. 11. Результаты численного эксперимента (красные ромбы) в сравнении с экспериментом (рис. 1-а). Для режимов с пульсирующей точкой отрыва показан

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

Выводы

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

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

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

Вычислительный эксперимент с использованием к-в-^ модели турбулентности показал, что «гистерезис» при запуске-остановке эжектора связан с нестационарным изменением параметров в процессе запуска-остановки эжектора. То есть при стационарных граничных условиях параметры в канале эжектора и вакуумной камере будут определяться однозначно. Однако при эксплуатации реальных установок время выхода на оптимальный режим при плавном повышении давления заметно больше, чем при сначала резком выходе на завышенные параметры и затем небольшом плавном дросселировании давления рабочего газа. Именно такой процесс отражает петля «гистерезиса».

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

Работа выполнена при финансовой поддержке гранта РФФИ № 16-38-60185.

Библиографический список

1. Соболев А.В., Запрягаев В.И., Мальков В.М. Одноступенчатый эжектор большой степени сжатия // Теплофизика и аэромеханика. 2005. Т. 12. № 1. С. 149 - 158.

2. Ларина Е.В., Ципенко А.В. Экспериментальные данные о потоке в газовом эжекторе для верификации моделей турбулентности // Труды МАИ. 2017. № 97. URL: http: //trudymai. ru/published. php?ID=87135

3. Verma S.B., Stark R., Genin C., Haidn O. Cold gas dual-bell tests in high-altitude simulation chamber // Shock Waves, 2011, vol. 21, issue 2, pp. 131 - 140, doi 10.1007/s00193-011 -0302-6

4. Kartovitskiy L., Lee J.H., Tsipenko A. Numerical and Experimental Investigation of Non-Stationary Processes in Supersonic Gas Ejector // 29th Congress of the International Council of the Aeronautical Sciences. ICAS 2014, ISBN CD: 3-932182-80-4.

5. Исаев А.И., Скоробогатов С.В. Гидродинамическая верификация и валидация численных методов расчета течения в камере сгорания газотурбинного двигателя // Труды МАИ. 2017. № 97. URL: http://trudymai.ru/published.php?ID=87336

6. Ivanov I.E., Kryukov I.A., Larina E.V. Effect of the turbulnet viscosity relaxation time on the modelling of nozzle and jet flows // Fluid Dynamics, 2014, vol. 49, issue 5, pp. 694 - 702, https://doi.org/10.1134/S00154628

7. Письменный В.Л. Математическая модель звукового газового эжектора с

цилиндрической камерой смешения в системе турбоэжекторного двигателя // Труды

МАИ. 2003. № 12. URL: http://trudymai.ru/published.php?ID=34456

8. Письменный В.Л. Турбоэжекторный двигатель // Труды МАИ. 2003. № 11. URL: http://trudymai.ru/published.php?ID=34476

9. Голубев В.А., Монахова В.П. Методы исследования эжекторных усилителей тяги (ЭУТ) // Труды МАИ. 2006. № 22. URL: http://trudymai.ru/published.php?ID=34102

10. Borisov V.E., Davydov A.A., Kudryashov I.Yu., Lutsky A.E., Men'shov I.S. Parallel implementation of an implicit scheme based on the LU-SGS method for 3D turbulent flows // Mathematical Models and Computer Simulations, 2015, vol. 7, issue 3, pp. 222 - 232, https://doi.org/10.1134/S2070048215030035

11. Tavoularis S., Corrsin S. Experiments in a nearly homogeneous shear flow with a uniform mean temperature gradient, Part 1 // Journal of Fluid Mechanics, 1981, no. 104, pp. 311 - 347.

12. Olsen M.E., Coakley T.J. The Lag Model, a Turbulence Model for Non Equilibrium Flows // 15th AIAA Computational Fluid Dynamics, 2001, pp. 2001 - 2564.

13. Gulyaev A.N., Kozlov V.E., Secundov A.N. A Universal One-Equation Model for Turbulent Viscosity // Fluid Dynamics, 1993, vol. 28, no. 4, pp. 485 - 494.

14. Wilcox D.C. Turbulence Modeling for CFD, DCW Industries, Inc., Griffin Printing, Glendale, California, 1994, XIX, 460 p.

15. Abramovich G.N. The theory of turbulent jets, The M.I.T. Press, Cambridge, Massachusetts, 1963, 671 p.

16. Zapryagaev V.I., Kudryavtsev A.N., Lokotko A.V. An experimental and numerical

study of a supersonic jet shockwave structure // West East High Speed Flow Fields -

CIMNE, Barcelona, Spain, 2002, pp. 346 - 351.

17. Chen Y.S., Applications of a new wall function to turbulent flow computations // AIAA Pap. 86-0438, 1986, 11 p.

18. Sarkar S., Erlebacher G., Hussaini M.Y., Kreiss H.O. The analysis and modeling of dilatational terms in compressible turbulence // Journal of Fluid Mechanics, 1991, vol. 227, pp. 473 - 493.

19. Yakhot V., Orszag S.A., Thangam S., Gatski T.B., Speziale C.G. Development of turbulence models for shear flows by a double expansion technique // Physics of Fluids A, 1992, no. (4)7, pp. 1510 - 1520.