CC BY
14
УДК 519.6 Б01: 10.15350/17270529.2020.2.17
МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ВЯЗКОГО ГАЗА СКВОЗЬ СУЖАЮЩИЕСЯ КОНИЧЕСКИЕ СОПЛА
1ЛИПАНОВ А. М., 2КАРСКАНОВ С. А.
1 Институт прикладной математики им. М.В. Келдыша РАН, 125047, Москва, Миусская пл., д. 4
Удмуртский федеральный исследовательский центр Уральского отделения РАН, 426067, г. Ижевск, ул. Т. Барамзиной, 34
АННОТАЦИЯ. Методом прямого численного моделирование рассчитывается нестационарное течение вязкого сжимаемого газа сквозь сужающиеся сопла. Расчеты выполнены в осесимметричной постановке. Угол сужения сопла является варьируемым параметром. Показано, что с ростом угла сужения расход в сопле падает. Определено, что интенсивные вихревые возмущения в сужающемся сопле затухают; причем, чем сильнее сопло сужается, тем быстрее гасятся возмущения.
КЛЮЧЕВЫЕ СЛОВА: прямое численное моделирование, сужающееся коническое сопло, вязкий поток, высокий порядок аппроксимации.
ВВЕДЕНИЕ
Возрастание требований к характеристикам ракетных и авиационных двигателей требует разработки сопловых блоков с улучшенными свойствами. Важным этапом при создании и тестировании современных двигательных установок является математическое моделирования физико-химических процессов, протекающих при работе устройства в сопле.
Численному исследованию процессов течения газа в соплах посвящено достаточное количество работ. Например, в работах [1, 2] численно исследуется возникновение отрыва в плоских и круглых соплах Лаваля. Моделирование проводится на основе решения двумерных нестационарных уравнений Навье-Стокса, описывающих динамику сжимаемого вязкого газа совместно с к-8 моделью турбулентности, применяется метод Годунова повышенного порядка точности. Рассматриваются течения в плоских, конических и круглых профилированных соплах, проводится детальное изучение структуры потоков и отдельных его элементов.
В [3] представлены результаты моделирования аэрогазодинамических характеристик эллиптического сопла смешанного расширения. Дано сравнение характеристик эллиптического и плоского сопел смешанного расширения с одинаковой относительной площадью миделя и входа в сопло. Приведены результаты по влиянию перепада давления в струе и формы поперечного сечения на тяговые характеристики сопла. Расчеты основаны на усредненных уравнениях Навье-Стокса в приближении тонкого слоя также с привлечением указанной выше модели турбулентности.
В работе [4] проводится математическое моделирование пространственного течения в осесимметричном сопле с кососрезанным выходным сечением при несимметричном перекрытии критического сечения с использованием модели Навье-Стокса вязкого теплопроводного газа. Используется консервативный численный метод конечных объемов, реализованный на вычислительных системах параллельной архитектуры.
Стоит отметить, что среди множества численных исследований большинство основывается на использовании тех или иных моделей турбулентности. Исследования, основанные на прямом численном моделировании, встречаются редко, так как сопряжены с вычислительными сложностями.
Настоящая работа посвящена моделированию методом DNS процессов течения газа сквозь сужающиеся сопла.
Целью работы является определение влияния угла сужения сопла на характер течения и распределение параметров потока.
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Решались уравнения Навье-Стокса в осесимметричной постановке, записанные для вязкого сжимаемого газа:
Ф + Фи*+ Фur , Pur = q Q)
dt dx dr r
dpUx , diPuX + P ~X xx) | d(puxur xr) | pUxUr -X xrr = Q (2)
dt dx dr r
Фи + d(puxllr -Xrx ) + d(puX + P -Xrr) + puj -Xrr + Xee = 0 (3)
dt dx dr r
dpE d (ux (pE + P) + qx - uxXxx - urXrx ) d (ur (pE + P) + 4r - uxXxr - urXrr )
dt dx dr
, ur (pE + P) + qr - uxXxr - urX
r
= 0. (4)
„ dux X (du„ duv
Xxx = -~PdlvU , Trx =Xxr = p|
f""", "rx -xr r^l ^ ^ '
dx 3 I dx dr j
xrr = Xp-^^-—p div u, Xee = Xp — - — p div u,
dr 3 r 3
xdL, dx qr = - -xdL, dr
_dux + dur + ^.
dx dr r
x
Для энергии можно записать
zr ul + u2r 1 P
E = —-r- +--. (5)
X k -1 p
Здесь х, г, I - цилиндрические координаты и время. Коэффициенты молекулярной вязкости - ц и теплопроводности - X полагаются постоянными.
Давление определяется из уравнения (5), а температура в уравнениях исключается с помощью уравнения Клапейрона-Менделеева:
Р/р = ЯТ.
При выполнении расчетов использовались безразмерные гидромеханические параметры (ГМП). В качестве масштабов брались: для линейных величин - величина радиуса Я сопла на входе, для скоростей - скорость звука с, для давления и плотности -величины давления РоШ и плотности роШ окружающей среды. В качестве окружающей среды выступал воздух, давление в окружающей среде равнялось одной атмосфере.
Граничные условия формулируются в соответствии с расположением вычислительной области, показанной на рис. 1. На входе задавались постоянными давление (р = 1.05 • Рш )
и температура (Тпп = 293 К). На выходе задавалось давление окружающей среды, остальные
параметры определялись с помощью условия «сноса» (равенство нулю первой производной на границе). На твердых обтекаемых поверхностях задавалось условие прилипания, кроме того, поверхности считались адиабатическими. На оси симметрии, совпадающей с осью г, задавались условия симметрии.
Счет велся на равномерной прямоугольной сетке, вблизи границ сопла сетка считалась неравномерной. Радиус цилиндрической части сопла на входе (1 единицы длины) составлял 1.0, затем сопло резко расширялось до радиуса 1.2, далее постепенно сужаясь на протяжении девятнадцати единиц (до л = 20). Угол сужения обозначен как а. Число типа Рейнольдса посчитанное по скорости звука (масштабу для скорости), равнялось Яеж = 2.06 • 10б.
Характерное число Рейнольдса на входе в канал было Яе = 0.5 • 106. 'к
0 1 20
Рис. 1. Область интегрирования
Мощность сетки составляла 11200x960 = 10 752 000 узлов. Вычислительный процесс распараллеливался с помощью технологии MPI. Расчеты выполнены на суперкомпьютере «Уран» (ИММ УрО РАН, г. Екатеринбург).
Для аккуратного интегрирования уравнений (1 - 4) с адекватным воспроизводством физических эффектов присущих течениям данного класса была применена следующая методика. Производные по пространству аппроксимировались с высоким порядком с помощью WENO-схем; по времени уравнения интегрировались со вторым порядком на основе TVD схемы Рунге-Кутты.
Итак, для расчета конвективных слагаемых (первых производных по пространству) использовалась конечно-разностная WENO-схема пятого порядка [5]. Поток в точке (i +1 / 2) записывался в виде суперпозиции потоков
3
f _ q(v) f(v )
Ji+1/2 Ji+1/2,
v=1
где Q(1) = 1/10, Q(2) = 6/10, Q(3) = 3/10 являются весовыми коэффициентами, а J - обобщенная переменная, соответствующая уравнениям (1 - 4), компоненты которой вычисляются через значения в соседних узлах:
/•(V) =
Л +1/2 =
Д1) = И / _ 7 / + 2 /•
&г+1/2 = 6 6 ■/г_1 + 6 2
+21)/2 = 2 /'+1 + 5 _ 1 _ 1
/1+1/ 2 = _ 1 2 + 5 ¿г+1 + 2 6 6 6
Схема WENO устанавливает веса линейной комбинации малыми для областей, содержащих разрыв и близкими к оптимальным весам для гладкого решения
3
¿■+1/2 = 2ю V ¿г+1/2,
V=1
(V) ^
Ю а(1) +а(2) +а(3), а (V) =. ^
(е +18 (V)) Р
Здесь 18- индикаторы гладкости [6], р = 2, а е - малая величина, предотвращающая
деление на 0, принималась равной 10_6.
По времени уравнения (1 - 4) интегрировались со вторым порядком точности. TVD-схема Рунге-Кутта второго порядка, используемая в работе, выглядит следующим образом [7]
W(1) = W"+Лt Ь(^) ,
W(и+1) =1 Wи +1 W(1) +1М L(W(1)), 2 2 2
где L - конечно-разностная аппроксимация. Для явной схемы значение шага для интегрирования по времени определяется с помощью предварительных тестовых вычислений, обеспечивающих стабильность счета. Требование к шагу по времени оказывается более строгим, чем у условия Куранта-Фридрихса-Леви, которое является необходимым, но не достаточным условием.
РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ
Было проведено прямое численное моделирование дозвукового протекания вязкого газа сквозь сужающиеся сопла. Выполнено четыре серии расчетов. Угол сужения изменялся от 0 до 1.5 градусов с шагом 0.5 градуса.
Поля распределения мгновенной и осредненной продольной компоненты вектора скорости потока в соплах с различными углами сужения показаны на рис. 2 и рис. 3 соответственно. На рис. 2 четко визуализируется вихревая дорожка за уступом в виде положительных и отрицательных экстремумов. Видно, что в сопле с углом сужения 1.5 градуса в окрестности выхода нет зон с отрицательными скоростями, то есть вихри, зародившиеся на входе, успели полностью диссипировать.
В цилиндрическом канале и сопле с малым углом сужения (0.5 град.) «медленные» приграничные зоны визуализируются почти по всей длине канала. В сопле с углом сужения 1.5 градуса поток заметно разгоняется, при этом приграничные области занимают значительно меньшее пространство, вытесняясь ядром потока.
х их/с
—1-0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 -0.05 -0.10 -0.15
Рис. 2. Распределение безразмерной мгновенной продольной компоненты вектора скорости в соплах с различными углами сужения: а) - 0°, б) - 0.5°, в) - 1.0°, г) - 1.5°
Рис. 3. Распределение безразмерной осредненной продольной компоненты вектора скорости в соплах с различными углами сужения: а) - 0°, б) - 0.5°, в) - 1.0°, г) - 1.5°
Данный тезис еще больше подтверждает распределение кинетической энергии турбулентности кг, показанное на рис. 4 для каналов с различными углами сужения. В цилиндрическом сопле с нулевым углом сужения ненулевая кинетическая энергия распределена по всей длине сопла. Видно, что уже при угле сужения 0.5 градуса кинетическая энергия турбулентности значительно падает. Это свидтельствует о снижении уровня пульсаций скорости во всей рассматриваемой области. В сопле с углом сужения в 1.5 градуса кт имеет лишь редкие всплески, оставаясь близкой к нулю почти во всей области.
Таким образом, вихревые возмущения, генерируемые на входе в сопло, при нулевом угле сужения сопла распространяются по всему потоку. С возрастанием угла конусности, то есть в более узких соплах, вихревые возмущения интенсивно гасятся, и уже при а = 1.0° практически не проникают вглубь сопла.
—г 0 01
0.006 X Ы- 0 005 0.004
Рис. 4. Распределение кинетической энергии турбулентности в соплах с различными углами сужения: а) - 0°, б) - 0.5°, в) - 1.0°, г) - 1.5°
Далее, рассмотрим распределение генерации энергии турбулентности
"к Н | ■
что описывает потери энергии для осредненного потока на подпитку пульсаций. В осессиметричной постановке формула для Рк приобретает вид:
"к=-< -< ^+% ]-<
дт дх
дт
Поля распределения генерации энергии турбулентности в различных соплах показаны на рис. 5. Видим, что в сопле с малым углом сужения и в цилиндрическом канале генерация происходит по всей длине сопла. При увеличении угла сужения сопла Р имеет значения отличные по модулю от нуля лишь в окрестности входа в сопло, где происходит активное вихреобразование за счет «разрывной» геометрии. Вглубь сопла вихревые возмущения не проникают, происходит их активное затухание.
а)
б)
в)
г)
Рис. 5. Распределение генерации энергии турбулентности в соплах с различными углами сужения: а) - 0°, б) - 0.5°, в) - 1.0°, г) - 1.5°
Численные расчеты показывают, что расход в сопле с ростом угла сужения падает. Данные таблицы наглядно свидетельствуют об этом. Площадь сечения (5) и расход (0) в цилиндрическом канале взяты за эталон. Падение расхода связано с уменьшением площади критического сечения сопла. Однако прямой зависимости между падением расхода и уменьшением площади критического сечения не наблюдается. Так, в сопле с углом сужения 1.5 градуса площадь критического сечения уменьшилась почти в три раза до 0.343*5, при этом расход упал только в два раза до 0.497*0.
Таблица
Относительные площадь критического сечения и расход
Угол сужения Площадь критического сечения Расход
0.0 1.000*5 1.000*0
0.5 0.743*5 0.859*0
1.0 0.524*5 0.658*0
1.5 0.343*5 0.497*0
ЗАКЛЮЧЕНИЕ
Проведено исследование влияние геометрических параметров сопла (угла сужения) на характер распределения гидромеханических параметров. Показано, что угол сужения сильно влияет на распространение вихревых возмущений, которые перестают проникать вглубь сопла.
Имеет место реламинаризация высокоскоростного завихренного течения за счет «разрывной» геометрии потока. С ростом угла сужения процессы реламинаризации только интенсифицируются.
Следует отметить, что расход в соплах с меньшим критическим сечением оказывается также меньше, то есть, чем больше угол сужения, тем меньше расход. Однако прямой зависимости между падением площади критического сечения и падением расхода не выявлено.
СПИСОК ЛИТЕРАТУРЫ
1. Глушко Г. С., Иванов И. Э., Крюков И. А. Численное моделирование отрывных течений в соплах // Физико-химическая кинетика в газовой динамике. 2010. Т. 9. С. 1-8.
2. Глушко Г. С., Иванов И. Э., Крюков И. А. Метод расчета турбулентных сверхзвуковых течений // Математическое моделирование. 2009. Т. 21, № 12. С. 103-121.
3 Мазуров А. П. Численное исследование характеристик эллиптического сопла смешанного расширения при сверхзвуковом обтекании // Ученые записки ЦАГИ. 2019. Т. 50, № 1. С. 57-72.
4. Бабаков А. В., Белошицкий А. В., Гайдаенко В. И., Дядькин А. А. Моделирование пространственного течения в осесимметричном сопле с асимметричным критическим и выходным сечениями // Математическое моделирование. 2018. Т. 30, № 11. C. 13-26.
5. Jiang G. S., Shu C. W. Efficient Implementation of Weighted ENO Schemes // Journal of Computational Physics, 1996, vol. 126, no. 1, pp. 202-228.
6. Shu C. W. High Order ENO and WENO Schemes for Computational Fluid Dynamics // High-Order Methods for Computational Physics. Lecture Notes in Computational Science and Engineering, 1999, vol. 9, pp. 439-582.
7. Gottlieb S., Shu C. W. Total Variation Diminishing Runge-Kutta Schemes // Mathematics of Computation, 1998, vol. 67, pp. 73-85.
Simulation of viscous gas flow through tapered conic nozzles
:Lipanov A. M.,,2Karskanov S. A.
1 Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences, Moscow, Russia
2 Udmurt Federal Research Center, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. The unsteady viscous gas compressible flow through the tapered nozzle is calculated. The calculations are performed in an axisymmetric formulation. The computational process was parallelized using MPI technology. Data was received using «Uran» supercomputer of IMM UB RAS. For accurate integration of Navier-Stokes equations with adequate reproduction of physical effects high-order approximation WENO scheme was used for calculation spatial derivative and second-order TVD Runge-Kutta scheme was used for time integration. The nozzle taper angle was a variable parameter. Fields of instantaneous and average velocity were analyzed. It is shown, that the taper angle increasing cause the flow rate decreasing. In a nozzle with a large tapering angle, the flow is noticeably accelerated, while the border regions occupy a much smaller space being displaced by the flow core. In a cylindrical nozzle with a zero taper angle nonzero kinetic energy is distributed over the entire length of the nozzle. In a nozzle with a taper angle of 1.5 degrees the kinetic energy of turbulence has only rare bursts, remaining close to zero in almost the entire region. The nozzle flow rate decreases when taper angle increasing. There is no direct relationship between the drop in flow and the decrease in the critical section area.
KEYWORDS: direct numerical simulations, tapered conic nozzle, viscous flow, high-order approximation.
REFERENCES
1. Glushko G. S., Ivanov I. E., Kryukov I. A. Numerical simulation of separated flow in nozzles. Physical-Chemical Kinetics in Gas Dynamics, 2010, vol. 9, pp. 1-8. http://chemphys.edu.ru/issues/2010-9/articles/143/
2. Glushko G. S., Ivanov I. E., Kryukov I. A. Computational method for turbulent supersonic flows. Mathematical Models and Computer Simulations, 2010, vol. 2, no. 4, pp. 407-422. https://doi.org/10.1134/S2070048210040010
3. Mazurov A. P. Numerical Investigation of Elliptical Mixed Expansion Nozzle Performance in Supersonic Flow. TsAGI Science Journal, 2019, vol. 50, iss. 1, pp. 71-89. https://doi.org/10.1615/TsAGISciJ.2019030084
4. Babakov A. V., Beloshitskiy A. V., Gaydaenko V. I., Dyadkin A. А. Simulation of a Spatial Flow in an Axisymmetric Nozzle with Nonsymmetrical Critical and Exit Sections. Mathematical Models and Computer Simulations, 2019, vol. 11, no. 3, pp. 457-468. https://doi.org/10.1134/S2070048219030050
5. Jiang G. S., Shu C. W. Efficient Implementation of Weighted ENO Schemes. Journal of Computational Physics, 1996, vol. 126, no. 1, pp. 202-228. https://doi.org/10.1006/jcph.1996.0130
6. Shu C. W. High Order ENO and WENO Schemes for Computational Fluid Dynamics. High-Order Methods for Computational Physics. Lecture Notes in Computational Science and Engineering, 1999, vol. 9, pp. 439-582. https://doi.org/10.1007/978-3-662-03882-6 5
7. Gottlieb S., Shu C. W. Total Variation Diminishing Runge-Kutta Schemes. Mathematics of Computation, 1998, vol. 67, pp. 73-85. https://doi.org/10.1090/S0025-5718-98-00913-2
Липанов Алексей Матвеевич, академик РАН, главный научный сотрудник Института прикладной математики им. М.В. Келдыша РАН, e-mail: aml35@yandex. ru
Карсканов Сергей Андреевич, кандидат физико-математических наук, старший научный сотрудник УдмФИЦ УрОРАН, е-mail: [email protected]
