Теплофизика высоких температур, 2021, T. 59, № 5, стр. 730-736

ВВЕДЕНИЕ

Исследование течений газа, реализующихся в результате взаимодействия ударной волны (УВ) с препятствиями различной природы и формы, представляет значительный интерес и рассматривается в различных постановках в большом количестве работ [1–9]. В настоящей работе проводится моделирование течения в замкнутой полости, имеющей форму полусферы, после вхождения в нее ударной волны. Данная задача содержит особенность, связанную с неограниченным ростом давления и температуры за ударной волной в момент ее отражения от плоского донышка. При моделировании течения газа, возникающего при фокусировке ударной волны, как правило, используются уравнения Эйлера и Навье‒Стокса, при этом параметры газа определяются в середине расчетной ячейки, примыкающей к оси симметрии. Такой способ делает решение сеточно-зависимым. В ряде работ считается, что в реальности в каналах всегда имеется пусть небольшое, но плоское донышко, на котором ставится условие непротекания [3, 10, 11]. Целью настоящей работы является анализ влияния координаты (радиуса) плоского донышка и термодинамическое описание воздействия свойств высокотемпературного воздуха на реализуемые в канале параметры течения.

ФИЗИКО-МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Исследуется течение, возникающее в канале с полусферическим торцом (рис. 1), после входа в него плоской УВ. Канал представляет собой трубу радиуса R, заканчивающуюся полусферическим торцом с плоским донышком. Радиус полусферического торца меняется от R (при х = R) до r – плоское донышко при х = X₀ (рис. 1). Исследуется влияние координаты плоского донышка на параметры течения за падающей и отраженной ударными волнами. В начальный момент времени в канале находится покоящийся воздух при нормальных условиях (P₀= 101325 Па, Т₀= 298.15 К). Параметры воздуха за ударной волной рассчитываются из соотношений Ренкина–Гюгонио при заданном числе Маха М ударной волны:

Здесь ρ, u, p, h – плотность, скорость, давление, удельная энтальпия газа соответственно, D – скорость ударной волны. Индексом “0” отмечено начальное состояние.

Рассматриваются три модели протекания физико-химических процессов.

Первая модель – “идеальная” – предполагает, что теплоемкость воздуха С_p и молекулярный вес M постоянны:

Вторая – “замороженная”, где считается, что воздух представляет собой смесь кислорода (O₂), азота (N₂) и аргона (Ar), для анализа термодинамических свойств которых используется модель термодинамики, описываемая удельным термодинамическим потенциалом Гиббса [12] для смеси совершенных газов:

где R – универсальная газовая постоянная, Т – температура, P₀ = 101 325 Па – стандартное давление, $G_{i}^{0}(T)$ – температурная часть стандартных молярных потенциалов Гиббса отдельных компонент, аппроксимационные функции для которых приведены в [12]. Соответствующие (3) термическое и калорическое уравнения состояния имеют вид

В третьей – “равновесной” [13] для описания воздуха используется смесь совершенных газов (3), (4), включающая 13 компонентов (O₂, N₂, NO, ${\text{O}}_{2}^{ + }$, ${\text{N}}_{2}^{ + }$, NO⁺, O⁺, N⁺, Ar⁺, O, N, Ar, e – электронный газ) и удовлетворяющая условиям термодинамического равновесия

Здесь N_e – число элементов в системе (в рассматриваемом случае N_e = 4: O, N, Ar, e), $A_{K}^{i}$ – матрица состава, $\gamma _{K}^{0}$ – заданные значения мольно-массовых концентраций элементов ($\gamma _{{{{N}_{e}}}}^{0} = 0$ – условие электронейтральности), ${{\mu }_{i}}$ – химический потенциал i-го компонента, z_k – неизвестные параметры, число которых равно числу элементов.

Необходимо отметить, что для описания термодинамических свойств компонентов при использовании второй и третьей моделей используются полиномы, приведенные в [12] для диапазона температур от 100 до 20000 К. В случае выхода за данные границы свойства веществ (3), (4) рассчитываются с фиксированными теплоемкостями, значения которых равны их значениям на границе данного диапазона. Соответственно, полученные результаты расчетов при температурах выше 20 000 К носят качественный характер.

Для предварительного анализа характеристик течения, реализуемых за ударными волнами, может быть рассчитана ударная адиабата (УА) воздуха (рис. 2) и прямая Михельсона, являющиеся следствием (1) и дополненные моделями термодинамики (2)–(5). Зависимости параметров течения вдоль адиабаты как функции скорости ударной волны представлены на рис. 3:

На рис. 2, 3 показано множество значений параметров, которые могут быть достигнуты за ударной волной, распространяющейся по покоящемуся воздуху при стандартных условиях. Ударные адиабаты, рассчитанные по всем трем моделям, совпадают на фазовой плоскости p–V (рис. 2) справа от точки v ~ 0.2 м³/кг, p ~ 1.35 × 10⁶ Па. Данной точке соответствуют скорость УВ D ~ 1180 м/с и температура T ~ 940 К. Далее идеальная УА располагается правее замороженной и равновесной. Равновесная и замороженная УА расходятся при v ~ 0.156 м³/кг, p ~ 3 × 10⁶ Па, D ~ 1740 м/с и T ~ ~ 1650 К. Затем равновесная УА проходит левее замороженной до точки их пересечения на плоскости p–v при v ~ 0.074 м³/кг, p ~ 1.36 × 10⁸ Па, D ~ ~ 11 200 м/с. При этом на равновесной УА T ~ ~ 17 300 К, а на замороженной – 34 800 К (рис. 3а). Взаимное расположение адиабат и различие значений параметров течения вдоль них определяются существенным отличием реальных термодинамических свойств веществ от идеальных, а также протеканием процессов диссоциации и ионизации при использовании “равновесной” модели (рис. 3б). Так, процессы диссоциации начинают оказывать заметное влияние при температуре, превышающей 3500 К (D ~ 3900 м/c), ионизации – при T > 11 000 К (D > 8000 м/c). Необходимо отметить, что давление и температура являются монотонно возрастающими функциями скорости ударной волны, а плотность монотонно возрастает на идеальной и замороженной адиабатах и имеет максимум на равновесной при D = 9200 м/c. Различие температур на идеальной, замороженной и равновесной адиабатах в случае высоких скоростей ударной волны может достигать нескольких десятков процентов, при этом давления различаются существенно меньше (рис. 3а).

Газодинамическое течение в канале считается квазиодномерным нестационарным, вязкость, теплопроводность и диффузия не учитываются. Для описания течения в областях непрерывности используются уравнения физической газовой динамики в дифференциальной форме [11, 14]

где e – удельная внутренняя энергия газа; F – площадь канала

На ударной волне выполняются соотношении Ренкина–Гюгонио (1).

Для замыкания системы (6) применялись три вышеописанные модели термодинамики: “идеальная”, “замороженная” и “равновесная”.

ОСОБЕННОСТИ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

Для моделирования рассматриваемого течения использовался оригинальный сеточно-характеристический метод, позволяющий рассчитывать квазиодномерные нестационарные течения многокомпонентного реагирующего газа с явным выделением произвольного числа взаимодействующих разрывов (ударных волн, контактных разрывов, характеристик семейств ${{C}^{ \pm }}$) [10, 11, 13–17]. При численном моделировании падающая и отраженные ударные волны явно выделялись, их координаты, скорости и параметры слева и справа рассчитывались с высокой точностью. В момент отражения УВ от стенки решалась задача об отраженной ударной волне с учетом используемой модели течения. Например, при моделировании течения воздуха с равновесным протеканием химических превращений [13] численно решалась система нелинейных алгебраических уравнений, включающая соотношения Ренкина–Гюгонио (1), условия термодинамического равновесия (5), уравнения состояния (4) и равенство нулю скорости потока за отраженной УВ. Особенностью моделирования течения с выделением падающей ударной волны является то, что область течения между плоским донышком и УВ может не рассчитываться вплоть до прихода УВ на донышко.

Расчетная сетка включала фиксированные узлы, которые на отрезке [X₀, X₁] сгущались к центру по геометрической прогрессии. На отрезке [X₁, X_k] использовалась равномерная сетка

Для расчета q и N₂ использовался следующий алгоритм. Вначале задавались X₁, h₁ и желаемый шаг вблизи плоского донышка $h_{2}^{0}$ ≤ X₀. Вычислялись $q = 1 + {{\left( {{{h}_{1}} - h_{2}^{0}} \right)} \mathord{\left/ {\vphantom {{\left( {{{h}_{1}} - h_{2}^{0}} \right)} {\left( {{{X}_{1}} - {{X}_{0}}} \right)}}} \right. \kern-0em} {\left( {{{X}_{1}} - {{X}_{0}}} \right)}}$, ${{{\text{N}}}_{2}} = \operatorname{int} \,\, \times $ $ \times \left( {\ln {{\left( {{{{{h}_{1}}} \mathord{\left/ {\vphantom {{{{h}_{1}}} {h_{2}^{0}}}} \right. \kern-0em} {h_{2}^{0}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{{{h}_{1}}} \mathord{\left/ {\vphantom {{{{h}_{1}}} {h_{2}^{0}}}} \right. \kern-0em} {h_{2}^{0}}}} \right)} {\ln \left( q \right)}}} \right. \kern-0em} {\ln \left( q \right)}}} \right) + 1$ и уточнялся шаг ${{h}_{2}} = $ $ = \frac{{\left( {{{X}_{1}} - {{X}_{0}}} \right)\left( {q - 1} \right)}}{{{{q}^{{{{{\text{N}}}_{2}}}}} - 1}}$.

РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

Для тестирования используемой методики моделирования была рассмотрена автомодельная задача об отражении ударной волны от центра симметрии, аналитическое решение которой было найдено в [6, 7]. В табл. 1 приведены результаты сравнения при Х₀= 10^–6, Х₁= 0.1, h₁= 0.005, Х_k = 1.1, начальное число Маха УВ – M_УВ= 4. Используются обозначения из работы [7], $n = {{\tau {{D}_{{{\text{УВ}}}}}} \mathord{\left/ {\vphantom {{\tau {{D}_{{{\text{УВ}}}}}} {{{X}_{K}}}}} \right. \kern-0em} {{{X}_{K}}}}$, D_УВ – начальная скорость УВ при x = X_k, τ – время, за которое УВ достигает левой границы.

При численном моделировании не использовалось предположение о том, что падающая УВ является сильной. Отличие расчетных данных от аналитического решения (табл. 1) составило не более 5%.

Исследовалось влияние координаты плоского донышка и числа Маха падающей УВ на параметры течения за падающей и отраженной УВ. Для определенности радиус канала и полусферического торца равнялся 1 м, координата плоского донышка X₀ варьировалась от 10^–6 до 10^–2 м, что соответствует радиусу плоского донышка от 0.0014 до 0.14 м, X_k = 1.1 м. Начальное число Маха падающей УВ варьировалось от двух до пяти. Рассматривались три вышеописанные модели воздуха: “идеальная”, “замороженная” и “равновесная”. Результаты численного моделирования приводятся на рис. 4, 5, а также в табл. 1 и 2. Из анализа рис. 4а видно, что координата левой стенки в диапазоне от 10^–2 до 10^–6 м незначительно сказывается на временно́й развертке процесса – траектории ударной волны. В приведенном масштабе они практически неразличимы. При этом используемая модель течения оказывает заметное влияние. В увеличенном масштабе (рис. 4б) видно, что процесс отражения УВ от стенки происходит в отличающиеся моменты времени и с существенно различными скоростями отраженных ударных волн.

В табл. 1, 2 даны значения параметров течения в характерных точках канала в момент прохождения через них ударной волны, рассчитанные по “идеальной” (и), “замороженной” (з) и “равновесной” (р) моделям, при координате плоского донышка X₀ = 10^–3–10^–6, числе Маха падающей УВ, равном трем (табл. 2) и пяти (табл. 3).

Цифрами 1 и 2 отмечены параметры в момент прихода падающей ударной волны на стенку (1 – параметры за падающей УВ, 2 – за отраженной). Цифрой 3 отмечены параметры воздуха за отраженной ударной волной в момент времени, когда она находится в сечении X = 0.05 м, цифрой 4 – Х = 1 м. Из табл. 1, 2 видно, что “идеальная” модель дает существенно завышенные значения давления и температуры за падающей ударной волной в момент ее прихода на стенку. За отраженной УВ отличие давления уменьшается, а “равновесная” температура отличается от “идеальной” от двух до шести раз при числе Маха падающей УВ, равном трем, и от 3.5 до пяти раз при числе Маха, равном пяти. В этом случае “равновесная” температура ниже “замороженной” почти в два раза. Пиковые значения давления и температуры достигаются, как и следовало ожидать, при X₀ = 10^–6 м. Необходимо отметить, что параметры газа, достигаемые в момент отражения, слабо влияют на значения параметров течения на некотором удалении от плоского донышка. Данный факт оправдывает постановку граничных условий на небольшом расстоянии от точки фокусировки. Как следует из проведенного исследования, вполне достаточно размещать границу расчетной области в точке с координатой, равной 10^–4R.

На рис. 5 приводятся мольно-массовые концентрации атомарных кислорода, азота и электронного газа при продольной координате, равной X = 10^–6, 0.05, 0.3 и 1 м, полученные при расчете по “равновесной” модели течения при начальном числе Маха, равном четырем. Концентрации в рассматриваемых сечениях начали заметно меняться после прохождения через них падающей УВ. Так, например, график концентрации атомарного кислорода (кривая 10), содержит два характерных скачка при t ~ 0.0001 c и t ~ 0.0022 c, связанных с прохождением падающей и отраженной ударной волны. Аналогичный характер имеют и остальные графики.

В окрестности плоского донышка наблюдаются высокие концентрации атомарных азота, кислорода и электронов, причем концентрация азота превышает концентрацию кислорода. Это связанно с тем, что рассчитанное значение температуры за отраженной ударной волной составляет 41556 К, соответственно, молекулярный азот практически полностью диссоциирует и при этом молекулярный вес смеси составляет ~11 г/моль.

ЗАКЛЮЧЕНИЕ

В работе рассмотрена задача о фокусировке и отражении от стенки ударной волны, распространяющейся в канале с полусферическим торцом радиуса R с плоским донышком. Варьировались число Маха падающей ударной волны (M = 2–5) и координата плоского донышка (X₀ = 10^–6–10^–2). Использовались три модели течения: “идеальная” – теплоемкости и молекулярный вес воздуха постоянны; “замороженная” – теплоемкости зависят от температуры, молекулярный вес – постоянный; “равновесная” – воздух моделируется находящейся в состоянии термодинамического равновесия смесью совершенных газов, в состав которой входят 13 компонентов. Показано, что для каждой модели процесса координата плоского донышка слабо влияет на параметры течения за отраженной ударной волной уже на расстоянии 0.05R от плоского донышка, что может быть использовано при постановке граничных условий. При X = R параметры течения, реализуемые за ОУВ, рассчитанные по “замороженной” и “равновесной” моделям, близки (отличие менее 2%). Использование “идеальной” модели в данном случае приводит при M = 3 к ошибкам по температуре и плотности 6%, а при M =  5–20%. Вблизи плоского донышка, как и следовало ожидать, наблюдается максимальное отличие параметров течения, рассчитанных по различным моделям. Так, температура за падающей ударной волной в момент отражения, рассчитанная по “равновесной” модели при M = 3 и X₀ = 10^–6_, ниже “идеальной” в пять раз, “замороженной” – в два раза, а при M = 5 и X₀ = 10^–6 ниже “идеальной” в восемь раз и “замороженной” в 2.5 раза. Отличия по давлению при этом существенно ниже и составляют соответственно 1.7, 1.1 и 1.9, 1.04. Результаты расчетов по “идеальной” модели существенно отличаются от расчетов, выполненных по “замороженной” и “равновесной” моделям. Параметры течения, полученные по “замороженной” модели, близки к “равновесным” на удалении от торца, но дают существенно завышенную температуру вблизи плоского донышка.

Работа выполнена по государственному заданию № FSFF-2020-0013.