CC BY
35
Vol. 19, № 04, 2016
Civil Avition High TECHNOLOGIES
УДК 533.6.011.72
ВЕРИФИКАЦИЯ ГИБРИДНОЙ ЧИСЛЕННОЙ СХЕМЫ ДЛЯ ЗАДАЧИ НАТЕКАНИЯ СЖИМАЕМОЙ СТРУИ НА ТВЕРДУЮ ПРЕГРАДУ
М.В. КРАПОШИН, C.B. СТРИЖАК
В статье рассматриваются вопросы математического моделирования истечения свободной сжимаемой турбулентной струи из модельного сопла и натекания струи на твердую преграду при различных степенях нерасчетно-сти. Для решения задачи используется разработанный авторами статьи решатель pisoCentralFoam на базе гибридной численной схемы Kurganov-Tadmor, алгоритма PISO и метода контрольного объема. Для расчета сжимаемой струи используется модель на базе нестационарных уравнений Рейнольдса и k-omega SST модели турбулентности с пристеночными функциями. Приведена постановка задачи для расчета натекания струи на преграду. Расчетная область представляла собой прямоугольник. Для упрощения постановки задачи рассматривалось только половина сопла. Для случая свободной струи на задания давления на выходе расчетной области использовалось смешанное граничное условие 3-го рода. На входе расчетной области задавалось специальное табличное условие для давления, которое позволяло постепенно поднимать абсолютное значений для давления. Значение степени нерасчетности струи выбиралось равным n = 2,5 и n = 5,0. Проведен анализ сеточной сходимости на сетках от 100 тысяч до 500 тысяч ячеек. Среднее значение величины y+ составило 270. Расчеты проводились для конечного времени Tend = 0,01 секунды. Получены результаты распределения поля модуля скорости, давления на оси симметрии. Проведено сравнение результатов расчета распределения давления по длине сопла на разных расчетных сетках с результатами эксперимента. Получено совпадение с результатами эксперимента с точностью в 5 %.
Ключевые слова: Сжимаемая струя, твердая преграда, бочка Маха, устойчивость, сетка, численная схема, нестационарный расчет.
ВВЕДЕНИЕ
Изучение струйных течений имеет большое практическое значение для решения задач в ракетно-космической технике, газокислородной резке металла и в вакуумной технике.
Задача расчета газодинамических параметров и акустического поля для случая натекания сжимаемой турбулентной струи на твердую преграду, расположенную перпендикулярно потоку, имеет большое фундаментальное и прикладное значение. К такой задаче можно отнести задачу об изучении газоструйной колебательной системы с соплом и резонатором Гельмгольца, задачу об изучении колебательных свойств струи при лазерной резке металла, задачу об изучении автоколебаний струи при натекании на преграду при различной степени нерасчетности струи. Известно, что перерасширенные струи характерны для случая работы первой ступени ракеты-носителя [1-3] и для случая лазерной резки металла [4]. Недорасширенные струи, в том числе при большой степени нерасчетности, характерны для случая разделения ступеней ракеты- носителя [5, 6]. Подобные течения характеризуются числом Маха на срезе сопла и степенью нерасчетности n, где n есть отношение давлений на срезе сопла к давлению в окружающей среде. Данные газодинамические параметры определяют структуру ударных волн и волн разряжения, наличие бочек с дисками Маха, неравномерное распределение давления, температуры в различных сечениях струи. Решение подобных задач позволяет оценивать силовые, тепловые и акустические нагрузки на элементы конструкции ракеты-носителя. В данной работе рассматривается задача численного моделирования взаимодействия недорасширенной струи, истекающей из модельного сопла, с твердой преградой при различных степенях нерасчетности с целью верификации выбранной численной схемы.
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Для решения подобных задач можно воспользоваться математическим моделированием и данными эксперимента с целью верификации математической модели. Известно, что для случая больших градиентов давления необходимо использовать специальные численные схемы на базе
Civil Avition High TECHNOLOGIES
Vol. 19, № 04, 2016
решения задачи о распаде произвольного разрыва. К таким схемам можно отнести схемы Годунова, Roe, HLLE, HLLE+, HLLE++ [7]. Данные схемы сложны для реализации в программном коде и не всегда дают устойчивый счет. Авторами данной статьи был предложен новый подход на базе гибридной численной схемы для неструктурированных сеток. В основе разработанного решателя используется гибридная схема, объединяющая возможности численной схемы Kurganov-Tadmor и численной схемы алгоритма PISO. На базе открытого пакета OpenFOAM, в основе которого лежит метод контрольного объема, была разработана новая программа на языке С++ для расчета течений при изменении числа Маха в широком диапазоне и новое граничное условия для учета особенностей задания величины давления для струйных течений на выходе расчетной области [8]. Исходная система уравнений включает в себя уравнения неразрывности, количества движения и энергии:
IP + V- (pü) = 0; (1)
Р+V ■ (pÜÜ) = V • П+i^; (2)
Р+ V^ (pije) = V^ (П-Ü) + Fb ■ Ü-V^ q, (3)
где p - плотность, Ü - скорость, Fb - главный вектор массовых сил, П - тензор нормальных и вязких напряжений, e - полная энергия потока (сумма кинетической и внутренней энергий), q - турбулентный поток тепла. В приближении ньютоновской жидкости тензор напряжения
вычисляется по формуле: П= - p + 2 V ■Üj/ + ¡VÜ +(VÜ) J, где p - давление, I - единичный тензор, 1 - динамическая вязкость. Для вычисления теплового потока используется закон Фурье: q = -XVVT, где T - температура среды, Л - коэффициент теплопроводности.
^ ~ R
Предполагается, что среда является совершенным газом, т. е. p = pRT, R =--индивидуаль-
M
ная газовая постоянная, R - универсальная газовая постоянная, M - молярная масса (вес) рассматриваемого газа. Идея гибридной численной схемы состоит в введении переключателя, который при достижении дозвуковой скорости переключит схему для определения потоков к стандтартному варианту PISO алгоритма:
\ V\Ü¥)V = \ dS ■ (Üy) = ^Wf (ü Л ) = ^VfVf, (4)
«-» direction «+» direction
P N
о — О
X,- X \J + 1 / 2 X./+1
X
где у/ - некоторая произвольная переносимая величина, / - номер грани на поверхности ячейки, - произведение вектора нормали к рассматриваемой грани на ее площадь. Величины с индексом / интерполируются на поверхность грани. Для сопоставления КТ/ККР (Киг§аиоу-Таёшог / Киг§аиоу-Кое11е-Ре1гоуа) схемы, предложенной в работе [9] со стандартными обозначениями, используемой в методе контрольного объема, обратимся к рис. 1. Точкой Р обозначен центр произвольного контрольного объема, для которого сформулированы балансовые уравнения для массы, количества движения и энергии. В одномерном случае этот объем имеет координату Х}.. Выражения
для потоков записываются на правой границе этого объёма, причём нормаль направлена вдоль
Рис. 1. Конечно-объемная вычислительная молекула
Vol. 19, № 04, 2016
Civil Avition High TECHNOLOGIES
оси X от точкиР к точке N, положение центра грани - X]+1/2, положение центра соседнего объема N - X]+1. Одним из важных свойств этой схемы является независимость аппроксимирующего выражения для потоков физических величин от характеристик исследуемой системы уравнений. В соответствии с КТ/ККР схемой поток величины поля у через грань может быть представлен как весовая сумма потоков в положительном и отрицательном направлениях. При этом в потоки входит не только макроскопическая скорость движения среды, но также акустические (звуковые) скорости распространения возмущений:
у((+1/2)((+1,2) = 1Д-'(У+1/2"1/21 (У1/2), (5)
a j+1/2 a j+1/2
aj+1/2 aj+1/2 aj+1/2 — aj+1/2
(j+1/2 -¥j+1/2 ).
Для перехода к системе обозначений, используемой в методе контрольного объем, введем следующие обозначения:
max a mm
a f м и 1
^f aj+1/2' af aj+1/2' af max min ' ^f „max , „min
ap + af J af + af
VPf =Wj+xn, VNf = w)+\n, Pf = f [wj+x/2), Pf = f [wj+x/2I ( = afa
Тогда, выражение для вычисления конвективного потока через грань принимает вид:
yrfpf = af pfyf + aNf p - (Qf (f -yf ), (7)
или после перестановки слагаемых:
WfPf =¥f \afPf +afaf ) + yf \af pf -afaf ), (8)
Здесь верхний индекс указывает ячейку, значение которой используется для интерполяции на грань f, если f то используется ячейка, по отношению к которой нормаль внешняя, если N то используется ячейка, по отношению к которой нормаль внутренняя.
Для вычисления коэффициентов af, а1^ используются объемные потоки af3^ и af1", относительные к локальной скорости звука в среде:
a^ах = max
№1 )+P +pP >.
a;in = min (-cff |Sf I) + pP , -of |Sf I + pNf ),
of
=Yf, ■ of=YRf, (9)
с
где у = —-— показатель адиабаты газа. Для расчета значений величин ур и у^ на грани ис-су
пользуются интерполяционные схемы второго порядка точности с ограничителями. Обычно в
Civil Avition High TECHNOLOGIES
Vol. 19, № 04, 2016
качестве ограничителя используются симметричные TVD схемы - vanLeer, Minmod, vanAlbad и другие. Данный подход реализован в решателе pisoCentralFoam для случая турбулентного нестационарного течения.
ПОСТАНОВКА ЗАДАЧИ И РЕЗУЛЬТАТЫ РАСЧЕТА
В качестве тестового примера выбиралась задача с результатами экспериментов, доступных в открытой печати: эксперимент NASA с плоским соплом и истечением сжимаемой турбулентной струи [10]. Диаметр выходного сечения модельного сопла da = 0,0493 м, угол наклона внутренней поверхности сопла равен 11 градусов, геометрическая степень расширения сопла равна 1,797. В эксперименте варьировался перепад давления n в сопле, т. е. отношение давления газа на входе в сопло к давлению окружающего пространства. Расчетная область представляла собой прямоугольник с размерами 0,85 x 0,65 x 0,01 м. Для упрощения постановки задачи рассматривалась только половина сопла.
В данной работе в качестве математической модели выбирался подход на базе URANS и модели турбулентности k-omega S ST. На входе расчетной области, совпадающей с входом в камеру сгорания сопла, задавалось специальное табличное условие для давления, которое позволяло постепенно поднимать абсолютное значений для давления. Во внешней области величина давления варьировалась. На выходе расчетной области задавалось смешанное граничное условие для давления для случая свободной струи, на стенках - условие непротекания, для величин k и omega на стенках выбирались условия для пристеночных функций. На оси симметрии расчетной области задавалось граничное условие симметрии. В ходе расчета варьировалось значение перепада давления n. Анализировались размеры бочек, длина струи, распределение давления, температуры в различных сечениях. Проведен анализ сеточной сходимости на сетках от 100 тысяч до 500 тысяч ячеек, сравнение с экспериментальными данными. Блочно-структурированная сетка и неструктурированная сетка были построены в пакете Salome.
Ранее проводились расчеты для случая натекания струи на твердую преграду с использованием различных моделей [11, 12]. Постановка задачи с использованием метода крупных вихрей (LES) накладывает серьезные требования на качество сетки, что не всегда бывает оправдано при проведении серии расчетов. Так, в работе [11] получен достоверный результат для случая натекания струи на преграду на сетке 22 млн ячеек. В другой работе, с использованием LES, расчеты проведены на сетке в 200 млн ячеек [12]. Для инженерной оценки иногда бывает достаточно выбрать подход на базе URANS с выбранной моделью турбулетности. Расчет истечения свободной турбулентной струи проведен для случая n = 2,5 и n = 5. Количество ячеек для многоблочной сетки составило 225900. Среднее значение величины y+ составило 270. На рис. 2 представлены результаты расчета поля модуля скорости в моменты времени от t = 0,001 до Tend = 0,01 секунды для случая свободной струи при n = 5.
Рис. 2. Значения поля модуля скорости для свободной струи
Vol. 19, № 04, 2016
Civil Avition High TECHNOLOGIES
Рис. 2 (продолжение). Значения поля модуля скорости для свободной струи
Рис. 3. Распределение давления и скорости на оси симметрии
p/piii ехр p/pin mesh No. 1 p/pín mesh No, 2
. 0.1 0,2 0.3 0.4 0.5 0.6 0.7 0.8 03 1 1.1 1,2 1,3 1.4 1.5 1.6 1.7 1.8 1.9 _ Non dimensional distance from inlet
Рис. 4. Сравнение величины безразмерного давления по длине сопла с результатами эксперимента [10]. Расчет проведен для двух вариантов сеток: 1) неструктурированная на базе тетраэдров; 2) блочно-структурированная на базе гексаэдров
На рис. 5 представлены результаты расчета иатекаиия струи на твердую преграду при n = 2,5. Представлено поле модуля скорости в моменты времени от t = 0,001 до Tend = 0,01 секунды.
Рис. 5. Результаты расчета натекания струи на твердую преграду
Civil Avition High TECHNOLOGIES
Vol. 19, № 04, 2016
Рис. 5 (продолжение). Результаты расчета натекания струи на твердую преграду
Рис. 6. Распределение давления и скорости вдоль оси симметрии
На рис. 7 представлены результаты расчета натекання струи на твердую преграду при n = 5. На рис. 8 представлены результаты расчета давления и скорости в момент времени Tend = 0,01 секунды для оси симметрии.
Рис. 7. Значения полей модуля скорости
Vol. 19, № 04, 2016
Civil Avition High TECHNOLOGIES
На рис. 3, 6, 8 представлено распределение давления и скорости на оси симметрии. Можно отметить, что параметры струи (количество бочек Маха, абсолютная величина давления) значительно меняются при изменении степени нерасчетности струи.
ЗАКЛЮЧЕНИЕ
Проведен расчет и визуализация значений полей скорости для случая натекания на твердую преграду. Расчеты проводились с использованием ресурсов web-лаборатории UniHUB.
Работы выполнялись в Федеральном государственном бюджетном учреждении науки «Институт системного программирования Российской академии наук» при финансовой поддержке Министерства образования и науки Российской Федерации (соглашение № 14.607.21.0090 о предоставлении субсидии от 24 ноября 2014, уникальный идентификатор соглашения RFMEFI60714X0090).
СПИСОК ЛИТЕРАТУРЫ
1. Антонов А.Н., Купцов В.М., Комаров В.В. Пульсации давления при струйных и отрывных течениях. М.: Машиностроение, 1990. 272 с.
2. Горшков Г.Ф., Усков В.Н. Автоколебания в сверхзвуковых перерасширенных импакт-ных струях // ПМТФ. 2002. Т. 43. N 5. С. 49-54.
3. Бирюков Т.П. и др. Газодинамика стартовых комплексов. М.: РЕСТАРТ, 2010. 364 с.
4. Ковалев О.Б., Фомин В.М. Физические основы лазерной резки толстых листовых материалов. М.: Физматлит, 2013. 256 с.
5. Горшков Г.Ф., Усков В.Н., Фаворский B.C. Особенности нестационарного обтекания безграничной преграды недорасширенной струей // ПМТФ. 1993. Т. 43. N 4. С. 58-65.
6. Varnier J., Requenet W. Experimental Characterization of the Sound Power Radiated by Impinging Supersonic Jet // AIAA Journal. 2002. Vol. 40. N 5. P. 825-831.
7. Куликовский А.Г., Погорелов H.B., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. 2-е изд. М.: Физматлит, 2012. 656 с.
8. Kraposhin M., Bovtrikova A., Strijhak S. Adaptation of Kurganov-Tadmor Numerical Scheme For Applying in Combination With the PISO Method in Numerical Simulation of Flows in a Wide Range of Mach Numbers // Procedia Computer Science. Vol. 66. 2015. P. 43-52.
9. Greenshields C.J. et al. Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for high-speed viscous flows // IJNMF. 63:1-21. 2010.
10. Asbury S.C., Hunter C.A. Static Performance of a Fixed-Geometry Exhaust Nozzle Incorporating Porous Cavities for Shock-Boundary Layer Interaction Control. NASA/TM-1999-209513. P. 135.
11. Dauptain A. et al. Large Eddy Simulation of Stable Supersonic Jet Impinging on Flat Plate // AIAA Journal. Vol. 48. No. 10. 2010.
12. Uzun A. et al. Simulation of Tonal Noise Generation for Supersonic Impinging Jets // AIAA Journal. Vol. 51. No. 7. 2013.
VERIFICATION OF HYBRID NUMERICAL SCHEME FOR THE CASE OF COMPRESSIBLE JET IMPINGIMENT ON FLAT PLATE
Kraposhin M.V., Strijhak S.V.
The article deals with the questions of mathematical modeling of compressible jet outflow from model nozzle and jet impingiment on flat plate at various values of n. pisoCentralFoam solver which is based on the Kurganov-Tadmor hybrid numerical scheme, PISO algorithm and finite volume method, is used for the solution of this problem. The model, based on unsteady Reynolds equation and K-omega SST turbulence model with boundary functions is used for compressible jet calculation. The problem definition for calculation of jet impingiment on flat plate is given. The simulation domain
Civil Avition High TECHNOLOGIES
Vol. 19, № 04, 2016
was selected as a rectangle. Only a half of the nozzle was considered for simplification. The mixed boundary condition for pressure setting in case of free jet was used on the outlet of simulation domain. The special condition for the pressure with table data, allowed to increase the value of pressure gradually, was used on the inlet of simulation domain. The value of the jet pressure degree was selected as n = 2.5 and n = 5.0. The results of distribution of the velocity magnitude, field pressure, upon symmetry axes were received. The simulations were done with grids 100 000-500 000 cells. The average value of y+ was equal to 270. The calculations were done for the end time Tend = 0.01 s. Comparison of the results of pressure distribution calculation based on nozzle length on different grids with the results of the experiment is carried out. The coincidence to engineering accuracy of 5 % is received.
Key words: Compressible Jet, flat plate, Mach barrel, stability, grid, numerical scheme, unsteady simulation.
REFERENCES
1. Antonov A.N. et al. Pulsatsii davleniya pri stryinykh i otryvnykh techeniyach. Moscow, Mash-inostroenie, 1990. 272 p. [In Russian].
2. Gorshkov G.F., Uskov V.N., Favorskii V.S. Avtokolebaniya v sverzvykhovykh perera-shirenykh impaktnykh stryakh. PMTF, 1993. T. 43, N 4. Pp. 58-65. [In Russian].
3. Byuryukov G.P. et al. Gazodinamika startovykh kompleksov. Moscow, RESTART, 2010. 364 p. [In Russian].
4. Kovalev O.B., Fomin V.M. Physizheskie osnovy lazernoy tolstykh listovych metallov. Moscow, Fizmatlit, 2013. 256 p. [In Russian].
5. Gorshkov G.F., Uskov V.N. Osobennosti nestasionarnogo obtekaniya bezgranichoy pregrady nedorashirenoy stryei. PMTF, 2002. T. 43, N 5. Pp. 49-54. [In Russian].
6. Varnier J., Requenet W. Experimental Characterization of the Sound Power Radiated by Impinging Supersonic Jet. AIAA Journal, 2002. Vol. 40, N 5. Pp. 825-831.
7. Kulikovsky A.G., Pogorelov N.V., Semenov A.Y. Matematicheskie voprosy chislennogo resheniya giperbolicheskih sistem uravnenii, 2-ed. Moscow, Fizmatlit, 2012. 656 p. [In Russian].
8. Kraposhin M., Bovtrikova A., Strijhak S. Adaptation of Kurganov-Tadmor Numerical Scheme For Applying in Combination With the PISO Method in Numerical Simulation of Flows in a Wide Range of Mach Numbers. Procedia Computer Science. Vol. 66, 2015. Pp. 43-52.
9. Greenshields A.J. et al. Implementation of semi-discrete, non-staggered central schemes in a collocated, polyhedral, finite volume framework, for high-speed viscous flows. IJNMF, 2010; 63: 1-21.
10. Asbury S.C., Hunter C.A. Static Performance of a Fixed-Geometry Exhaust Nozzle Incorporating Porous Cavities for Shock-Boundary Layer Interaction Control. NASA/TM-1999-209513, p. 135.
11. Dauptain A. et al. Large Eddy Simulation of Stable Supersonic Jet Impinging on Flat Plate. AIAA Journal, Vol. 48. No. 10, 2010.
12. Uzun A. et al. Simulation of Tonal Noise Generation for Supersonic Impinging Jets. AIAA Journal. Vol. 51. No. 7, July 2013.
СВЕДЕНИЯ ОБ АВТОРАХ
Крапошин Матвей Викторович, старший научный сотрудник ИСП РАН.
Стрижак Сергей Владимирович, инженер ИСП РАН, к.т.н.
