CC BY
18
Вычислительные технологии
Том 24, № 5, 2019
Математическое моделирование газодинамических процессов в импульсной аэродинамической установке и расчет некоторых параметров потока в рабочей части*
Н. П. СкивинА
Томский государственный университет, Россия Контактный e-mail: uss.skibina@gmail.com
Описывается математическое моделирование газодинамических процессов, протекающих в основных функциональных частях импульсной аэродинамической установки. Осуществляется расчет локальных и интегральных параметров обтекания модели конуса, расположенной в рабочей части установки. Приводится качественное и количественное сравнение результатов математического и физического моделирования.
Ключевые слова: вычислительная гидродинамика, математическое моделирование, аэродинамическая установка, коэффициент силы лобового сопротивления.
Библиографическая ссылка: Скибина Н.П. Математическое моделирование газодинамических процессов в импульсной аэродинамической установке и расчет некоторых параметров потока в рабочей части // Вычислительные технологии. 2019. Т. 24, № 5. С. 38-48. DOI: 10.25743/ICT.2019.24.6.004.
Введение
Разработка сверх- и гиперзвуковых летательных аппаратов как военного, так и гражданского назначения приводит к необходимости изучения аэродинамики, тепломассообмена и физической механики полета на высоких скоростях. Актуальность исследований в данном направлении подтверждается национальными и международными программами, реализуемыми на территории Европы, США, Индии, Японии и Российской Федерации. Результаты научно-исследовательских и опытно-конструкторских работ все чаще встречаются в публикациях зарубежных и отечественных авторов, в новостных заголовках и научно-популярных статьях.
Одной из задач наземных испытаний летательных аппаратов и их частей является воспроизведение натурных и близких к ним условий эксперимента при исследовании аэродинамических характеристик модельных объектов. Соблюдение соответствующих критериев подобия, например чисел Маха и Рейнольдса, позволяет моделировать
* Title translation and abstract in English can be found on page 48.
© ИВТ СО РАН, 2019.
различные эффекты. Экспериментальные исследования с использованием аэродинамических установок — один из наиболее распространенных способов определения локальных и интегральных характеристик моделей и отдельных частей летательного аппарата. Классические эксперименты направлены:
— на определение сил и моментов, действующих на тело (весовые испытания);
— определение давлений на поверхности тела (дренажные испытания);
— получение картины обтекания (визуализация).
Импульсная аэродинамическая установка Томского государственного университета используется для проведения аэродинамических и аэрофизических исследований в сверх- и гиперзвуковых потоках с числом Маха М = 2... 7. Основными достоинствами установки являются широкий диапазон режимов работы, простая конструкция, низкие расходы сжатого газа и электроэнергии.
Осуществление измерений в высокоскоростных потоках — довольно трудоемкий процесс, так как его организация и проведение требуют учета различных аспектов. Например, необходимо учитывать длительность процесса обтекания и подбирать соответствующие регистрирующие приборы. Затем выбранные устройства следует разместить так, чтобы их влияние на течение было минимальным — не возникали нежелательные искажения и срывы потока. Полученные данные подвергаются обязательной статистической обработке, что приводит к необходимости проведения серий экспериментов и влечет за собой значительные временные, энергетические и финансовые затраты. Менее ресур-созатратным способом исследований является численный расчет. При этом критерием истинности остается эксперимент: сравнительный анализ результатов расчета и результатов, полученных опытным путем, позволяет судить о достоверности данных математического моделирования и применимости подхода в целом.
Цели работы — математическое моделирование процессов, протекающих в основных частях импульсной аэродинамической установки [1], и изучение обтекания тел, закрепленных в ее рабочей части.
1. Физическая постановка задачи
В процессе подготовки к численным расчетам проведен обзор технической документации, в том числе инструкции по эксплуатации импульсной аэродинамической установки [2], выполнены обработка данных и перенос в цифровой формат геометрических размеров форкамеры, набора осесимметричных сопел, рабочей части и диффузора. Цикл работы аэродинамической установки представляет собой движение воздуха внутри перечисленных элементов. В форкамере происходит выравнивание потока, который затем поступает в профилированное сопло, где разгоняется до необходимой скорости и уже в рабочей части встречается с закрепленным телом, инициируя процесс сверх- или гиперзвукового обтекания. Для изучения обтекания тел в рабочей части использованы модели в виде клина и конуса [3]. Выбор тел данной геометрии обусловлен тем, что для них ранее были проведены серии экспериментов с использованием импульсной аэродинамической установки. Модели имели одинаковое сечение с углом полураствора 15 град., но структура течения вблизи тел различна.
Получены результаты дренажных испытаний [3], в ходе которых осуществлялось измерение статического давления в сходных точках на поверхности клина и конуса. Для конуса проводились весовые испытания, направленные на измерение аэродинамического коэффициента силы лобового сопротивления. Результаты весовых испытаний
Сх 0.50
2 3 4 5
Рис. 1. Коэффициент сопротивления модели в виде конуса: I — результаты экспериментов, • — данные [4], ▲ — данные [5]
Рис. 2. Фрагменты визуализации процесса сверхзвукового обтекания конуса, закрепленного в рабочей части (M = 3): а — процесс установления, б — стационарный режим обтекания
в сравнении с данными других авторов приведены на рис. 1. При проведении всех экспериментов осуществлялась съемка процесса обтекания высокоскоростной видеокамерой, два фрагмента приведены на рис. 2.
На рис. 2, а представлено начало процесса обтекания конуса потоком воздуха со скоростью, соответствующей M = 3. Можно видеть равномерный поток заданной скорости, выходящий из сопла, и формирование ударной волны в головной части тела, при этом процессы за ударной волной еще не установившиеся. Кадр на рис. 2, б относится к уже полностью установившемуся процессу: отчетливо видна зона стационарного течения, называемая ромбом Маха, где все параметры потока имеют постоянные значения.
2. Математическая постановка задачи
Для решения задачи использовался пакет вычислительных программ Ansys Fluent, широко применяемый для гидро- и газодинамических расчетов. Приняты следующие допу-
щения: силы инерции значительно преобладают над силами вязкости, характер течения турбулентный, отсутствует теплообмен между потоком рабочего газа и моделью. В рамках поставленной задачи осуществлялось решение системы нестационарных осреднен-ных уравнений Навье — Стокса [6], включающей в себя уравнение неразрывности
ж + v с и) = 0-
уравнение сохранения импульса (количества движения)
9 (рИ)
т
уравнение сохранения энергии
+ V (pU 0 U) = —Vp + Vт,
^^¡Г -1 + V (рUhtot) = V (XeffVT) + V (U r)
^VU + (VU)T - V U j
J 2
где т = peff ( VU + (VU) — -SV U J — тензор вязких напряжений; peff = p + pt;
Срр^ к 1 2
Ле// = Л + ——, ^ = р—; = к + -Ц , к = СрТ; р — давление, кш — полная гГ( ш 2
энтальпия, к — энтальпия, р — плотность газа, и — вектор скорости, Т — температура. Уравнение состояния имеет вид
Р
R ' T
м„
где Mw — молекулярная масса; R — универсальная газовая постоянная.
Для описания турбулентного характера течения применялся RANS-подход, в качестве модели турбулентности, замыкающей исходную систему уравнений, выбрана SST-мо-дель [7] (Shear Stress Transport model). Особенностью данной модели является эффективное сочетание точной формулировки к — ш-модели турбулентности в пристеночной области с независимым описанием поля удаленного течения с помощью к — е-модели. Стандартная к—ш-модель и преобразованная к—е-модель одновременно умножаются на некоторую функцию смешивания и складываются друг с другом. Функция смешивания сконструирована таким образом, что в пристеночной области принимает значение 1 и активирует стандартную к — ш-модель, а при значении 0, вдали от поверхности, активирует преобразованную к—е-модель. Приведенные выше особенности делают SST-модель более точной и надежной для широкого класса потоков [8], например для потоков с неблагоприятными градиентами давления, аэродинамических поверхностей, трансзвуковых течений и ударных волн.
Система уравнений для SST-модели имеет следующий вид:
д д д ( дк\ - (Рк) + — (Р кчг) = — (Tfc —J —Yk,
д д д ( дк\ - (Рш) + — (РшЩ) = ^ (Г* д^) —Y + D,
где Г k, Г* и Yk, Y* — эффективная диффузия и диссипация величин в турбулентном потоке для к и ш соответственно.
Турбулентная вязкость в ББТ-модели рассчитывается как
рк 1
№
ш
тах
1 ЗД
а* а\ш
где $ — модуль тензора деформаций [7];
1
Ок
1
^ 1 - ^
о.
К 1- П
+
+
Функции смешивания или весовые функции имеют следующую формулировку:
= 1апЬ (Ф1) ,
/к 500^'
Ф1 = тт
д+
тах
тах
(
4рк
0.09шу ру2шу оШ)2Щу2 1 1 дк дш
2р
аш,2 ш дху дх3
10
-10
= 1апЬ (Ф2)
Ф2 = тах 2
(2
/к 500^
0.09шу ру2ш
)
где у — расстояние до ближайшей поверхности и И+ — положительная часть кросс-диффузионного члена [7]. Функция Р1 отвечает за переключение между моделями турбулентности — вблизи поверхностей к — ^-модель, вдали от поверхностей к — е-модель. Функция Г2 активирует ограничитель в формуле для определения турбулентной вязкости. Слагаемое представляет собой диссипацию кинетической энергии турбулентности и определяется аналогично стандартной к — ^-модели:
П = Р/Зш2.
ББТ-модель основана на стандартной к — ^-модели и к — е-модели, чтобы объединить их вместе, к — е-модель была преобразована в уравнения, основанные на к и е. Это привело к возникновению термина кросс-диффузии:
2(1 — р
1 дк дш шаш>2 дх^ дх^
Значения модельных констант: ак,1 = 1.176, аш,1 = 2.0, ак,2 = 1.0, аш,2 = 1.168, а1 = 0.31, Рг,1 = 0.075, = 0.0828. Все дополнительные модельные константы а^, а^, а0, Д*,, Яр, Як, Яш, С * и М^ имеют те же значения, что и в стандартной к — ^-модели [7, 8].
В качестве граничных условий на входе в расчетную область задаются значения полного давления и полной температуры.
Полное давление
Ро = ps( 1 + ^M2) 7-1,
где ps — статическое давление на входе, 7 — показатель адиабаты, M — число Маха.
Полная температура
То = Ts (1 + M^ .
Здесь Ts — статическая температура.
Условие прилипания задается на всех твердых стенках и поверхности тела:
их = 0, иу = 0,
где Ux и Uy — компоненты вектора скорости в декартовой системе координат (х, у). На выходе из рабочей области задаются мягкие граничные условия и соответствующие граничные условия для оси симметрии [7, 8].
3. Численный расчет
Решение системы уравнений производилось методом конечных элементов. Для этого расчетная область, полностью дублирующая внутреннюю геометрию импульсной аэродинамической установки, разбивалась на конечное число расчетных элементов, после чего путем дискретизации система дифференциальных уравнений сводилась к системе алгебраических уравнений, решаемой итерационно. Для создания сеточных моделей использовался программный пакет Gambit 2.4.6. Каждая сеточная модель в среднем содержит от 150 000 до 500 000 расчетных элементов. Такое число элементов выбрано в процессе проверки сеточной сходимости [9] — дальнейшее уменьшение размеров расчетных элементов и увеличение их числа не оказывают влияния на точность получаемого решения. Расчетные элементы имеют форму четырехугольников, соотношение длин сторон которых находится в пределах, соответствующих хорошему и отличному качеству сеточной модели. Сама сетка является структурированной и сгущается вблизи стенок аэродинамической установки и исследуемого тела.
На первом этапе работы рассмотрено движение воздуха внутри форкамеры, сопла (M = 2), рабочей части и диффузора. В результате анализа полученных картин обтекания и распределения газодинамических параметров выявлено, что описанная выше математическая модель работает правильно: движение потока воздуха направлено от форкамеры к рабочей части, при прохождении через критическое сечение сопла скорость возрастает до величины, соответствующей M = 2.08, что с погрешностью 0.5% отвечает скорости, получаемой в экспериментах.
Затем исследовалось влияние расстояния между срезом сопла и входным сечением диффузора, которое может изменяться в диапазоне от 100 до 250 мм в зависимости от габаритов модели, закрепленной в рабочей части. В экспериментах с клином и конусом это расстояние составляет 145 мм — такое относительное расположение сопла и диффузора приводит к тому, что вихри, образующиеся на выходе из сопла, не контактируют со стенками конструкции и выходят в диффузор, не оказывая существенного влияния на характер потока.
4. Обзор результатов
Дальнейшая работа над топологией расчетной области и генерацией сеток была направлена на численный расчет аэродинамического коэффициента силы лобового сопротивления для тела, закрепленного в рабочей части, на определение величины статического давления в фиксированной точке на теле и получение визуализации обтекания.
Численный расчет производился для импульсной аэродинамической установки с соплами 2... 5 (номер сопла приблизительно соответствует числу Маха на входе в рабочую часть). На рис. 3 приведены фрагменты визуализации, полученные в ходе расчета для варианта с соплом 3. Фрагмент на рис. 3, а соответствует началу процесса обтекания — представлены встреча потока заданной скорости с телом и формирование ударной волны. Фрагмент на рис. 3, б соответствует полностью установившемуся процессу обтекания тела [10] в рабочей части — отчетливо видна зона стационарного течения, в которой газодинамические параметры имеют постоянные значения. Относительное рассогласование по величине скорости потока между расчетом и экспериментом составило 1.72%. Следует отметить, что фрагменты, представленные на рис. 2, а и 3, а, описывают одинаковые процессы, как и на рис. 2, б и 3, б. Таким образом, можно говорить, что результаты численного расчета не только количественно, но и качественно согласуются с данными, полученными в ходе эксперимента. Аналогичные результаты получены для сопел 2, 4, 5. Контроль сходимости численного расчета осуществлялся путем наблюдения за невязками. Кроме того, при достижении коэффициентом силы лобового сопротивления установившегося значения с погрешностью 0.5 % решатель автоматически оповещал о достижении сходимости. На рис. 4 представлены результаты численного расчета для сопел 4 и 5.
На рис. 5 приведено сравнение экспериментальных [2] и расчетных данных по величине статического давления. Видно, что характер его изменения в зависимости от скорости потока сохраняется, но результаты математического моделирования обтекания имеют погрешность порядка 10 % при сопоставлении с экспериментальными данными для клина и конуса. В таблице приведены значения чисел Маха, полученные
Рис. 3. Фрагменты визуализации процесса обтекания, полученные при численном расчете (Ы = 3): а — процесс установления; б — стационарный режим обтекания
а
Contours OFMach Number (Tlms=1.7530fr+{I1) Sep 26.2010
AN5YS Fluem 15.0 farf. dp, dt>ns окр, Irans-sst, Lranslenli
Рис. 4. Поле распределения скоростей потока в числах Маха по аэродинамической установке с соплами 4 (а) и 5 (б)
Рис. 5. Сравнение результатов численного расчета с экспериментальными значениями статического давления для моделей в виде клина и конуса
Относительное рассогласование результатов физического и численного экспериментов
№ Вычисленное Экспериментальное Ш,%
сопла значение M значение M
2 2.15 2.09 2.87
3 2.96 2.91 1.72
4 3.91 4.03 3.19
5 4.94 5.02 1.59
в результате численных расчетов и серий экспериментов, также их относительное рассогласование, составившее 2-3 %.
Выдвинуто предположение, что полученное рассогласование по величине статического давления на поверхности модели обусловлено спецификой реализации решения двумерной задачи. Для учета структуры течения и распределения давлений при обтекании моделей плоской (клина) и осесимметричной (конуса) формы в дальнейшем будет осуществляться решение задачи в трехмерной постановке, что позволит получить значения, приближенные к экспериментальным.
Таким образом, проведено математическое моделирование газодинамических процессов в импульсной аэродинамической установке. Отлажена методика, позволяющая производить численный расчет параметров обтекания тела в рабочей части при заданной скорости набегающего потока. Сформированная база знаний открывает возможности размещения в рабочей части установки тел более сложной формы с целью изучения их обтекания в сверх- и гиперзвуковых потоках. Грамотная состыковка эксперимента и математического моделирования позволяет получать достоверные данные для широкого диапазона аэро- и термогазодинамических параметров и дает информативную картину их изменения в исследуемом процессе.
Благодарности. Работа выполнена при финансовой поддержке Программы повышения международной конкурентоспособности Томского государственного университета на 2013-2020 гг.
Список литературы / References
[1] Звегинцев В.И. Газодинамические установки кратковременного действия. Ч. I. Установки для научных исследований. Новосибирск: Параллель, 2014. 551 с. Zvegintsev, V.I. Gas dynamic devices of short duration. Pt I. Devices for scientific research. Novosibirsk: Parallel, 2014. 551 p. (In Russ.)
[2] Модельная аэродинамическая установка. Техническое описание. Новосибирск: ИТПМ им. С.А. Христиановича СО РАН, 2009. 48 c.
Model aerodynamic devices. Technical description. Novosibirsk: ITPM im. S.A. Khristianovich SO RAN, 2009. 48 p. (In Russ.)
[3] Маслов Е.А., Жарова И.К., Фарапонов В.В., Мацкевич В.В., Чижов С.Ю.
Физическое моделирование обтекания тел плоской и осесимметричной формы сверхзвуковым потоком // V Междунар. молодеж. науч. конф. "Актуальные проблемы современной механики сплошных сред и небесной механики", 25-27 ноября 2015 г. Томск: Изд-во Том. ун-та, 2016. С. 84-88. DOI: 10.17223/9785751124199/17.
Maslov, E.A., Zharova, I.K., Faraponov, V.V., Matskevich, V.V., Chizhov, S.Yu.
Physical modeling of supersonic flows around plane or axisymmetric bodies // V Mezhdunar. Molodezh. Nauch. Konf. "Aktual'nye problemy sovremennoy mekhaniki sploshnykh sred i nebesnoy mekhaniki", 25-27 noyabrya 2015 g. Tomsk: Izd-vo Tom. Un-ta, 2016. P. 84-88. DOI: 10.17223/9785751124199/17. (In Russ.)
[4] Артонкин В.Г., Леутин П.Г., Петров К.П. Аэродинамические характеристики острых и притупленных конусов при дозвуковых и сверхзвуковых скоростях // Труды ЦАГИ. 1972. Вып. 1413. 92 с.
Artonkin, V.G., Leutin, P.G., Petrov, K.P. Aerodynamic characteristics of sharp and dulled cones at subsonic and supersonic speeds // Trudy TsAGI. 1972. Iss. 1413. 92 p. (In Russ.)
[5] Ш^елудько Ю.В. Измерение донного давления осесимметричных тел малого удлинения // Физико-газодинамические баллистические исследования / Под ред. Г.И. Мишина. Л.: Наука, 1980. С. 68-77.
Sheludko, Yu.V. Measurement of the bottom pressure of slightly elongated axisymmetric bodies // Physical-gas-dynamic ballistic studies / Ed. G.I. Mishina. Leningrad: Nauka, 1980. P. 68-77. (In Russ.)
[6] Скибина Н.П., Савкина Н.В., Фарапонов В.В. Численное определение аэродинамических характеристик цилиндрических тел с разным видом геометрии при сверхзвуковых скоростях // Матер. XVII Всерос. конф. молодых ученых по матем. моделированию, 30 октября — 3 ноября 2016 г. Новосибирск, Россия. Новосибирск: ИВТ СО РАН, 2016. С. 67-68.
Skibina, N.P., Savkina, N.V., Faraponov, V.V. Numerical determination of the aerodynamic characteristics of cylindrical bodies with different types of geometry at supersound speeds // Proc. of the XVII All-Russ. Conf. of Young Scientists on Math. Modeling. October 30 — November 3, 2016, Novosibirsk, Russia. Novosibirsk: ICT SB RAS, 2016. P. 67-68. (In Russ.)
[7] ANSYS Fluent User's Guide. ANSYS, Inc., 2013. 2692 p. Available at: http://www.pmt.usp. br/academic/martoran/notasmodelosgrad/ANSYS7o20Fluent%20Users%20Guide.pdf
[8] Снегирев А.Ю. Высокопроизводительные вычисления в технической физике. Численное моделирование турбулентных течений: Учеб. пособ. СПб.: Изд-во Политехн. ун-та, 2009. 143 с.
Snegirev, A.Yu. High performance computing in engineering physics. Numerical modeling of turbulent flows: Ucheb. posob. SPb.: Izd-vo Polytechn. Univ., 2009. 143 p. (In Russ.)
[9] Скибина Н.П., Савкина Н.В., Фарапонов В.В. Численное моделирование обтекания цилиндрического тела сверхзвуковым потоком // Фундамент. и прикл. пробл. совр. механики (ФППСМ-2016): Сб. тр. IX Всерос. науч. конф., 21-25 сентября 2016 г., Томск. Том. гос. ун-т, 2016. С. 151-153.
Skibina, N.P., Savkina, N.V., Faraponov, V.V. Numerical simulation of supersonic flow around a cylindrical body // Fundamental and Applied Problems of Modern Mechanics (FPSPM-2016): Proc. of the IX All-Russ. Sci. Conf., Sept. 21-25, 2016, Tomsk. Tom. Gos. Univ., 2016. P. 151-153. (In Russ.)
[10] Альбом сверхзвуковых течений / Под ред. П.И. Ковалева, Н.П. Менде. СПб.: Изд-во Политехн. ун-та, 2011. 251 с.
Album of supersonic flows / Eds P.I. Kovalev, N.P. Mende. SPb.: Izd-vo Polytechn. Univ., 2011. 251 p. (In Russ.)
Поступила в 'редакцию 26 сентября 2019 г.
48
H. n. CKHÔHHa
Mathematical modeling of gas-dynamic processes in a pulsed aerodynamic facility and the numerical calculation of flow parameters in the test chamber
Skibina, Nadezhda P.
Tomsk State University, Tomsk, 634050, Russia
Corresponding author: Skibina, Nadezhda P., e-mail: uss.skibina@gmail.com
Purpose. The aim of this paper is the mathematical modelling of gas-dynamic processes and the incoming airflow around the model fixed in the test chamber of the main functional elements of a pulsed aerodynamic facility.
Methodology. Physical and mathematical modelling methods are jointly used for comprehensive description of the processes under study. The experimentally obtained results of drainage and balance measurements combined with visualization of the flow around the model using the shadow method are used to validate the numerical solution of the problem. System of unsteady Reynolds averaged Navier — Stokes (RANS) equations and SST turbulence model are used for mathematical description of the movement of working gas flow from pre-chamber through the nozzle to the test chamber where the body was fixed. Finite elements method implemented in Ansys Fluent computational complex is used for the numerical solution.
Findings. Mathematical modelling of the flow around the cone model by the incoming flow with Mach numbers M = 2 ... 5 was carried out. The distribution fields of gas-dynamic parameters (p, P, T, M) over the entire internal volume of the pulsed aerodynamic facility (pre-chamber, axisymmetric nozzles set, test chamber) are obtained. The values of static pressure at the model surface and Mach number at nozzle exit are calculated for each velocity. Relative error between numerical calculation and experimental data is obtained to be of the order of 10 and 3 %, respectively. For flow structure near the model we obtain a qualitative agreement that the density distribution field is identical to the shadowgraph flow visualization pattern.
Value. A numerical method for calculation of local and integral parameters of supersonic and hypersonic flow around models in test chamber of aerodynamic facility was implemented. Information about the gas-dynamic flow parameters at the nozzle exit allowing the study of the flow around the model for the cases with more complex geometry using mathematical modelling based on the results of experiments.
Keywords : computational fluid dynamics, mathematical modeling, aerodynamic facility, aerodynamic drag coefficient.
Cite : Skibina, N.P. Mathematical modeling of gas-dynamic processes in a pulsed aerodynamic facility and the numerical calculation of flow parameters in the test chamber // Computational Technologies. 2019. Vol. 24, No. 5. P. 38-48. (In Russ.) DOI: 10.25743/ICT.2019.24.5.004.
Acknowledgements. This work has been done in the framework of TSU Competitiveness Enhancement Program 2013-2020.
Received September 26, 2019
© ICT SB RAS, 2019
