CC BY
39
Известия Института математики и информатики УдГУ
2015. Вып. 2 (46)
УДК 519.6, 533
© И.М. Кузьмин, Л.Е. Тонкое, A.A. Чернова
МОДЕЛИРОВАНИЕ ТУРБУЛЕНТНОСТИ В СОПЛАХ С ВЫСОКОЙ СТЕПЕНЬЮ РАСШИРЕНИЯ И ОЦЕНКА ПАРАМЕТРОВ БОКОВОЙ СОСТАВЛЯЮЩЕЙ ТЯГИ1
Рассматривается численное решение задачи о нарушении осевой симметрии течения и возникновения боковой составляющей тяги на начальном этапе работы сверхзвукового сопла с высокой степенью геометрического расширения. Приводится сравнение результатов, полученных с использованием k—ш SST модели турбулентности для осредненных по Рейнольдсу уравнений Навье^Стокса, и результатов моделирования методом отсоединенных вихрей SST DDES.
Ключевые слова: метод конечных объемов, метод отсоединенных вихрей, сверхзвуковое сопло, боковая составляющая тяги.
Введение
Проблема нарушения осесимметричности течения и возникновения боковых нагрузок на начальном этапе работы присуща в той или иной степени большинству конструкций ракетных двигателей первой ступени, что обусловлено прежде всего высокой геометрической степенью расширения сопла. С этими явлениями в разное время сталкивались и американские исследователи (ЖРД J-2, J-2S, J-2X), и европейские (Vulcain), и японские (LE-7A) [1].
При рассмотрении процессов, протекающих в начальный момент работы сверхзвукового сопла ракетного двигателя, можно выделить несколько характерных особенностей, присущих достаточно широкому классу конструкций. Прежде всего, это формирование системы скачков уплотнения в сверхзвуковой части сопла, обусловленное перерасширенным режимом течения. При значительной величине градиента давления вблизи стенки сопла возникает характерная А-конфигурация скачков уплотнения и отрыв потока. Как экспериментально, так и теоретически было установлено [2], что существует два хорошо различимых типа отрывного течения: неограниченный отрыв, при котором поток отделен от стенки развитой областью возвратного течения, и ограниченный отрыв, когда после отрыва пограничного слоя происходит его присоединение на некотором расстоянии ниже по потоку, и затем возможны новые отрывы и присоединения.
Неизбежно существующая в любой реальной газодинамической системе неоднородность параметров потока приводит к тому, что на некоторых режимах реализуется особый тип отрыва, имеющий существенно трехмерный и нестационарный характер, что и приводит к возникновению динамических высокочастотных боковых нагрузок при возрастании давления в камере сгорания на этапе запуска двигателя [3]. Данные боковые нагрузки, в свою очередь, приводят не только к колебаниям конструкции, но и к изменению формы сопла, нарушению его осесимметричности, что дополнительно сказывается на характере истечения газа.
Применяемые при численном исследовании начального этапа работы сверхзвукового сопла математические модели и численные схемы, очевидно, должны адекватно учитывать каждый из перечисленных механизмов, участвующих в формировании асимметричного течения и боковой составляющей тяги. Весьма важным здесь является исследование влияния вида граничных условий во входном сечении и способа их численной реализации, так как в силу специфики моделируемого явления именно распространение возмущений вниз по потоку от входной границы области определяет пространственно-временные характеристики процессов формирования зон отрывного течения и, как следствие, потери симметрии течения.
1 Работа выполнена при частичной поддержке РФФИ, проект № 14-08-00064-а.
№ б)
Рис. 1. Расчетная сетка: а) образующая плоская сетка с обозначением вида границ; б) общий вид расчетной сетки
§ 1. Постановка задачи и метод численного решения
Рассмотрим модель нестационарного внутреннего течения совершенного газа по соплу в виде следующих уравнений сохранения:
dp/dt + V • pU = 0,
dpU/dt + V • pUU = -Vp + Va, (1.1)
dpE/dt + V • pUE = -p(V • U) - V • q - V • (a • VU),
где a — тензор вязких напряжений, E = ^U2 + cvT — полная энергия, q — тепловой поток. Остальные обозначения общепринятые.
Расчетная область включает внутреннее пространство сопла и область свободного течения. Профиль стенки сверхзвуковой части сопла соответствует приведенному в [4] (ЖРД J-2S). Построение вычислительной сетки осуществлялось вращением относительно оси симметрии сопла плоской сетки, содержащей около 28 ООО четырехугольных ячеек (рис. 1, а). Выбранный способ построения позволяет при проведении вычислений полностью исключить влияние сетки как механизма потери решением осевой симметрии. Все дальнейшие результаты получены на сетке, состоящей в общей сложности из 4 242 ООО элементов, 46 500 из которых — треугольные призмы (рис. 1, б"). Аналогично [1] на входной границе 1 (рис. 1, а) задается зависимость полной энергии от времени, на стенках сопла 2 — условия прилипания и отсутствия тепловых потоков. На внешней границе области 3 — условия дозвукового истечения. Все используемые в расчетах параметры среды совпадают с принятыми в [4]. В начальный момент времени область заполнена неподвижным газом с параметрами окружающей среды. Более подробно особенности постановки начальных и граничных условий приведены в [6].
Система (1.1) аппроксимировалась методом конечных объемов, что позволило естественным образом использовать понятие обобщенного решения для течений, содержащих разрывы в виде ударных волн. Численное интегрирование уравнений сохранения выполнялось при помощи неявного трехшагового алгоритма с коррекцией давления PISO, который сохраняет свои аппроксимационные свойства и устойчивость при расчете до- и сверхзвуковых течений в широком диапазоне чисел Маха [5]. Для обеспечения повышенного порядка аппроксимации по пространству применялась кусочно-полиномиальная реконструкция функций внутри дискретной ячейки, что, в свою очередь, потребовало введения функции-ограничителя, снижающей порядок аппроксимации до первого вблизи поверхностей разрывов. Во всех расчетах использовался ограничитель minmod Ван-Лира.
§ 2. Моделирование турбулентности
Очевидно, что решение поставленной задачи интегрированием только системы (1.1) на рассматриваемой сетке без привлечения дополнительных моделей турбулентного переноса не поз-
воляет надеяться на получение физически содержательного решения, и прежде всего это относится к областям отрывного течения, определяющим характер потери течением симметрии. Однако, как будет показано далее, решение в приближении ламинарного течения может быть полезным с точки зрения оценки влияния численной диссипации на динамику боковой составляющей силы тяги.
Одним из перспективных методов моделирования турбулентных течений является метод отсоединенных вихрей (DES), основная идея которого состоит в том, что «точно» разрешаются только крупные вихревые структуры, а присоединенные пограничные слои описываются полуэмпирическими моделями турбулентности (RANS). Для рассматриваемой задачи, где структура течения в равной степени определяется как крупными вихрями в области за прямым скачком уплотнения, так и малыми в пристенных областях отрывного течения, применением моделирования отсоединенных вихрей возможно достичь существенного сокращения вычислительных затрат по сравнению с методом крупных вихрей (LES) и, тем более, прямым численным моделированием турбулентности (DNS).
В работе применялся метод отсоединенных вихрей в модификации SST DDES [7] с фильтром Гаусса шириной А, равной кубическому корню из объема ячейки дискретизации. С формальной точки зрения система (1.1) дополняется двумя уравнениями переноса — кинетической энергии турбулентности и удельной скорости диссипации:
dpk/dt + v ■ pUk = v ■ ((u + Ckßsos)vk) + Pk - pk3/2/e, (2.1)
дрш/öt + V • pUw = V • ((/i + CufisGsWu) + a—-—Pk - Âpw2 + (1 - F^Dk», (2.2)
ßSGS
при этом в системе (1.1) тензор напряжений примет вид
<T = -(fi + ßsGs)(V и + (VU)T) + + ßsGs)(VV)I,
Usgs ~ pk/ш — «подсеточная» (вихревая) вязкость. Линейный масштаб турбулентности I = min (k1/2 / s А), Cdes = const разделяет области LES и RANS моделирования
турбулентности; значения констант можно найти в [7].
Таким образом, имея программную реализацию численной процедуры совместного решения (1.1), (2.1), (2.2), можно практически без изменения кода получить решение исходной зада-
k w SST модели
Ментера для осредненных по Рейнольдсу уравнений Навье-Стокса и турбулентного по модели SST DDES.
§ 3. Результаты численного моделирования
Основной целью серии проведенных вычислительных экспериментов являлось исследование симметрии получаемых решений. Для каждого из трех рассматриваемых приближений было выполнено несколько повторений расчета параметров течения для одних и тех же начальных и граничных условий, так как в силу природы моделируемого процесса осреднение по испытаниям не может быть заменено осреднением по времени. Все расчеты выполнялись на вычислительной системе АПК' 11)11 Межведомственного суперкомпьютерного центра РАН.
Прежде всего следует отметить, что потеря течением симметрии наблюдалась в каждом из расчетов, в том числе и в приближении ламинарного течения при стационарных граничных условиях. На рисунке 2 показана типичная мгновенная рассчитанная картина течения, соответствующая отрезку времени, когда диск Маха колеблется вблизи среза сопла. Несимметричное относительно оси сопла расположение отрывных зон, взаимодействующих с системой скачков, приводит к ярко выраженному нарушению симметрии течения и, как следствие, возникновению боковой составляющей тяги.
На рисунке 3, б показана одна из рассчитанных в ламинарном приближении реализаций проекции годографа вектора тяги на плоскость, нормальную к оси симметрии в сравнении с численными результатами [4] (рис. 3, а). Максимум амплитуды боковой составляющей хорошо
Рис. 2. Изолинии числа Маха (t = 0.15 с)
-4е+04 -2е+04
2е+04 4е+04
-4е+04 -2е+04
2е+04 4е+04
-4е+04 -2е+04 0 2е+04 4е+04
а.)
-4е+04 -2е+04 0е+00 2е+04 4е+04
-4Р+04 -2Р+04 Ор+ОО 2р+04 4е + 04
4е+04 Зе+04 2е+04 1е+04 О
-1е+04 -2е+04 -Зе+04 -4е+04
б)
-4е+04 -Зе+04 -2е+04 -1е+04 О
-1е+04 -2е+04 -Зе+04 -4е+04
-4е+04 -2е+04 0 2е+04 4е+04
-4е+04 -2е+04 0е+00 2е+04 4е+04
—2е+04 —2е+04-
-4Р+04 -2Р+04 Ор+ПП 2Р+04 4Р+04
г)
Рис. 3. Годограф вектора боковой составляющей силы тяги Fxy (Н): а) результаты численного моделирования [4]; б) расчет в ламинарном приближении; в) расчет по модели k—w SST; г) расчет по модели SST DDES
согласуется с данными [4], но во всех выполненных расчетах годограф имеет характерный вид (рис. 3, б) с точностью до поворота относительно начала координат.
На рисунке 3, в показана проекция годографа, полученная решением осредненной по Рей-нольдсу системы уравнений Навье-Стокса (1.1), дополненной k w SST моделью турбулентности. Несмотря на ярко выраженные колебания боковой составляющей силы тяги, их амплитуда весьма мала, кроме того, проведенные вычисления показали, что изменения даже на несколько порядков величин k и w на входной границе не позволяют получить картину, близкую к показанной на рисунке 3, а. По-видимому, в данном случае необходимо задание нестационарных граничных условий для U и p, что требует отдельного дополнительного исследования.
Весьма обнадеживающие результаты (рис. 3, г) позволяет получить модель отсоединенных вихрей SST DDES. Полученный годограф боковой составляющей тяги качественно значительно ближе к результатам [4] (рис. 3, а), но частота колебаний очевидно выше и в процессе доминируют отклонения средней амплитуды. Здесь необходимо отметить, что применение в расчетах турбулентных течений моделей отсоединенных вихрей и их модификаций (DES, DDES, IDDES) уже само по себе требует весьма значительного объема тестовых вычислений с целью подбора параметров модели, оценки влияния расчетной сетки и способа аппроксимации исходных уравнений [8]. Но, несмотря на значительные вычислительные затраты, есть все основания считать, что моделирование отсоединенных вихрей при решении рассматриваемой задачи наиболее перспективно.
Заключение
Единственным первичным механизмом нарушения осевой симметрии численного решения в каждой из рассмотренных здесь моделей является погрешность округления при выполнении арифметических операций с плавающей точкой. Следовательно, можно считать, что изначально осесимметричное течение возмущается некоррелированным по времени случайным процессом с нулевым средним значением, то есть белым шумом, амплитуда которого имеет порядок машинного нуля. Это дополнительно подтверждает, что наблюдаемые в численном эксперименте явления имеют физическую природу.
Количество выполненных испытаний мало для получения статистически устойчивых оценок, но предварительно можно заключить, что основные параметры несимметричного течения (максимум амплитуды, момент времени, когда он достигается) воспроизводятся от расчета к расчету даже в такой простой постановке без учета деформаций стенок сопла.
Предварительные результаты, полученные применением метода отсоединенных вихрей в модификации SST DDES [7], продемонстрировали целесообразность его использования как с точки зрения эффективности вычислений, так и с позиций наиболее полного соответствия получаемых результатов экспериментальным данным среди рассмотренных ранее подходов.
Список литературы
1. Tomita Т., Sakamoto Н., Onodera Т., Sasaki М., Takahashi М., Tamura Н. Experimental evaluation of side load characteristics on TP, CTP, and TO Nozzles // AIAA Paper. 2004. № 04. 3678. 9 p.
2. Zmijanovic V., Rasuo В., Chpoun A. Flow separation modes and side phenomena in an overexpanded nozzle // FME Transactions. 2012. Vol. 40. №3. P. 111-118.
3. Wang T.-S. Transient three-dimensional startup side load analysis of a regeneratively cooled nozzle // Shock Waves. 2009. Vol. 19. №3. P. 251-264. DOI: 10.1007/s00193-009-0201-2
4. Zhao X., Bayyuk S., Zhang S. Aeroelastic response of rocket nozzles to asymmetric thrust loading // Computers and Fluids. 2013. Vol. 76. P. 128-148. DOI: 10.1016/j.compfluid.2013.01.022
5. Demirdzic I., Lilek Z., Peric M. A collocated finite volume method for predicting flows at all speeds // International Journal for Numerical Methods in Fluids. 1993. Vol. 16. № 12. P. 1029-1050.
6. Копысов С.П., Тонкое Л.Е., Чернова А.А. Постановка граничных и начальных условий при моделировании процесса запуска сопла // Химическая физика и мезоскопия. 2013. Т. 16. № 2. С. 216-222.
7. Paik J., Sotiropoulos F., Porte-Agel F. Detached eddy simulation of flow around two wall-mounted cubes in tandem // Int. J. Heat and Fluid Flow. 2009. Vol. 30. №2. P. 286-305.
DOI: 10.1016/j.ijheatfluidflow.2009.01.006
8. Gritskevich M.S., Garbaruk A.V., Menter F.R. Fine-tuning of DDES and IDDES formulations to the k w shear stress transport model // Progress in Flight Physics. 2013. Vol. 5. P. 23-42.
DOI: 10.1051/eucass/201305023
Поступила в редакцию 01.10.2015
Кузьмин Игорь Михайлович, младший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: imkuzmin@gmail.com
Тонков Леонид Евгеньевич, к. ф.-м. н., старший научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: tnk@udman.ru
Чернова Алена Алексеевна, к.т.н., научный сотрудник, Институт механики УрО РАН, 426067, Россия, г. Ижевск, ул. Т. Барамзиной, 34. E-mail: alicaaa@gmail.com
I. M. Kuz'min, L. E. Tonkov, A. A. Chernova
Turbulence modeling approaches to prediction of the overexpanded supersonic nozzles side load
Keywords: finite volume method, detached-eddy simulation, the side load, the supersonic nozzle. MSC: 76H05, 76F65, 65M08
A simulation of gas flow in a overexpanded axisymmetric nozzle at a start-up stage is performed. The calculations exposed in this paper point out a side load phenomena. Computations are carried out with the use of the Reynolds-averaged Navier-Stokes к-ш SST model and Delayed Detached-Eddy Simulations (SST DDES). The results are discussed with respect to the turbulence modelling technique.
REFERENCES
1. Tomita T., Sakamoto H., Onodera T., Sasaki M., Takahashi M., Tamura H. Experimental evaluation of side load characteristics on TP, CTP, and TO Nozzles, AIAA Paper, 2004, no. 04, 3678 9 p.
2. Zmijanovic V., Rasuo B., Chpoun A. Flow separation modes and side phenomena in an overexpanded nozzle, FME Transactions, 2012, vol. 40, no. 3, pp. 111-118.
3. Wang T.S. Transient three-dimensional startup side load analysis of a regeneratively cooled nozzle, Shock Waves, 2009, vol. 19, no. 3, pp. 251-264. DOI: 10.1007/s00193-009-0201-2
4. Zhao X., Bayyuk S., Zhang S. Aeroelastic response of rocket nozzles to asymmetric thrust loading, Computers and Fluids, 2013, vol. 76, pp. 128-148. DOI: 10.1016/j.compfluid.2013.01.022
5. Demirdzic I., Lilek Z., Peric M. A collocated finite volume method for predicting flows at all speeds, Internationa! Journal for Numerical Methods in Fluids, 1993, vol. 16, no. 12, pp. 1029-1050.
6. Kopysov S.P., Tonkov L.E., Chernova A.A. Positing of boundary and entry conditions for modeling progress of start up of a nozzle, Khimicheskaya fizika i mezoskopiya, 2013, vol. 16, no 2. pp. 216-222 (in Russian).
7. Paik J., Sotiropoulos F., Porte-Agel F. Detached eddy simulation of flow around two wall-mounted cubes in tandem, Int. J. Heat and Fluid Flow, 2009, vol. 30, no. 2, pp. 286-305.
DOI: 10.1016/j.ij heatfluidflow .2009.01.006
8. Gritskevich M.S., Garbaruk A.V., Menter F.R. Fine-tuning of DDES and IDDES formulations to the k-ш shear stress transport model, Progress in Flight Physics, 2013, vol. 5, pp. 23-42.
DOI: 10.1051/eucass/201305023
Received 01.10.2015
Kuz'min Igor' Mikhailovich, Junior Researcher, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: imkuzmin@gmail.com
Tonkov Leonid Evgen'evich, Candidate of Physics and Mathematics, Senior Researcher, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: tnk@udman.ru
Chernova Alena Alekseevna, Candidate of Engineering, Researcher, Institute of Mechanics, Ural Branch of the Russian Academy of Sciences, ul. T. Baramzinoi, 34, Izhevsk, 426067, Russia. E-mail: alicaaa@gmail.com
