Известия РАН. Механика жидкости и газа, 2021, № 6, стр. 125-135

ПОСТАНОВКА ЗАДАЧИ И ЧИСЛЕННАЯ МОДЕЛЬ

Рассматривается взаимодействие равномерного сверхзвукового потока вязкого совершенного газа с бесконечным проницаемым экраном в виде плоской периодической решетки параллельных цилиндрических стержней. Плоскость решетки перпендикулярна вектору скорости V₁ набегающего сверхзвукового потока. Показатель адиабаты газа γ = 1.4. Число Рейнольдса, вычисленное по параметрам невозмущенного потока и диаметру D цилиндрического элемента решетки Re = 10⁵. Число Маха M₁ набегающего потока варьируется в диапазоне от 2.4 до 4.2. Период решетки h = 5D обеспечивает достаточно большую степень проницаемости экрана σ = 1 – D/h = = 0.8, при которой скачки уплотнения перед элементами решетки остаются локализованными и не образуют глобальной ударной волны перед проницаемым экраном [12]. При этом позади экрана формируется вторичный сверхзвуковой поток [2, 10] как результат слияния дискретных струек, образовавшихся между элементами решетки.

Течение газа описывается системой двумерных уравнений Навье–Стокса с граничными условиями прилипания на изотермической поверхности элементов решетки. Решение строится численно в полосе шириной h между двумя продольными плоскостями, на которых ставится условие периодичности решения по поперечной координате Y. Внутри полосы находится один элемент решетки (рис. 1a).

Результаты сравнительных расчетов [13] обтекания системы элементов конечной решетки и обтекания одного элемента с условиями периодичности для бесконечной редкой решетки подтверждают отсутствие качественных отличий, в частности, сохраняется разнообразие схем обтекания бесконечной решетки. При этом граничное условие периодичности позволяет существенно улучшить разрешение расчетной сетки.

Входная граница расчетной области располагается в невозмущенном потоке с числом Маха M₁ ≥ 2.4. Выходная граница располагается далеко вниз по потоку в области сверхзвукового следа и на ней ставятся граничные условия свободного вытекания. В дальнейшем все линейные размеры в области течения нормируются на радиус элемента решетки, в частности, D = 2.

При γ = 1.4 и σ = 0.8 необходимое условие существования стационарного сверхзвукового потока перед бесконечной решеткой [12]

сводится к условию M₁ > 2.15. По расчетам [11] переход от режима с ветвящейся системой локальных ударных волн перед элементами решетки к режиму с гладкой отошедшей ударной волной происходит на интервале 2.35 < M₁ < 2.4. Перечисленные факты были учтены при выборе нижней границы M₁ = 2.4 интервала варьирования параметра M₁ в настоящем исследовании. Верхняя граница M₁ = 4.2 выбрана из условия отсутствия альтернативы схеме регулярного пересечения головных скачков уплотнения от соседних цилиндрических элементов решетки [11]. При M₁ > 4.2 угол 2α регулярного пересечения этих скачков меньше, чем теоретический предел 2 α_N существования тройной маховской конфигурации в этой точке. Здесь α_N = f(M₁, γ) – угол Неймана [15, 16].

Решение поставленной краевой задачи строится на основе многоблочной вычислительной технологии [7, 12, 17] с использованием локальных адаптированных к поверхности тел криволинейных сеток “B”, имеющих конечные области перекрытия с глобальной прямоугольной сеткой “A” для всей расчетной области (рис. 1б,в). Вязкие пограничные слои разрешаются с использованием уравнений Навье–Стокса на локальных сгущающихся к поверхности тел сетках, а эффекты интерференции ударно-волновых структур в сверхзвуковом следе описываются в рамках уравнений Эйлера. В областях перекрытия сеток применяется интерполяция функций до границ перехода от одной сетки к другой. Математическая модель и метод интегрирования детально описаны в [17]. В основной серии настоящих расчетов шаг глобальной сетки “A” (рис. 1) составляет 0.0125 (в долях радиуса цилиндрического элемента решетки); минимальный шаг локальной сетки “B” вблизи границ элементов решетки – 0.0005; общее число расчетных узлов – 3 270 000. Отладочные и контрольные расчеты выполнялась на сетках вдвое меньшей и вдвое большей густоты.

Численные решения строятся методом установления по времени. Для диагностики возможных гистерезисных режимов применяется тактика [11, 12] последовательного изменения числа Маха M₁ на входной границе расчетной области, при которой в качестве начального газодинамического поля берется решение с предыдущего малого шага по M₁. Для визуализации результирующих полей течения около решетки рассчитываются изолинии числа Маха, давления, модуля градиента плотности, а также характерные линии тока. Коэффициенты С_x, С_y аэродинамической силы, действующей на элемент решетки, определяются путем интегрирования давления по поверхности этого элемента. Отличия С_x от значений C_x0 коэффициента сопротивления одиночного цилиндра в неограниченном сверхзвуковом потоке характеризуют влияние аэродинамической интерференции между элементами решетки. Полученные расчетные величины C_x0 на интервале 2.4 ≤ M₁ ≤ 4.2 достаточно хорошо описываются аппроксимационной зависимостью C_x0 = = 1.39–0.097ln(M₁ – 1).

РЕЗУЛЬТАТЫ ВЫЧИСЛЕНИЙ

На рис. 2а,б представлены картины течения для значений M₁ на концах интервала варьирования (густота оттенков серого цвета пропорциональна модулю градиента плотности газа). При M₁ = 4.2 возмущения от соседних элементов решетки взаимодействуют между собой в сверхзвуковой части аэродинамического следа, не оказывая влияния ни на ударный слой перед цилиндром, ни на ближний след (рис. 2a). Головная ударная волна перед цилиндром и основная часть ближнего следа точно такие же, как при обтекании цилиндра в неограниченном сверхзвуковом потоке. Совпадают коэффициенты аэродинамического сопротивления элементов решетки и одиночного цилиндра: C_x = C_x0 = 1.277. В области сверхзвукового течения позади решетки формируется сложная ударно-волновая структура. Головные ударные волны “1” от соседних элементов сначала пересекаются между собой и с хвостовыми скачками уплотнения “2” по схеме регулярного взаимодействия, а затем образуют тройную маховскую конфигурацию в области сверхзвукового следа за промежуточным элементом решетки (рис. 2в). При этом за маховской ножкой “3” в окрестности оси следа образуется локализованная вторичная дозвуковая область, а уходящая ветвь тройной маховской конфигурации “4” догоняет преломленный фронт “5” хвостового скачка уплотнения и сливается с ним, образуя сверхзвуковую λ-конфигурацию скачков уплотнения со слабым контактным разрывом [18]. Следует отметить, что в расчетах [11] внимание уделялось только зонам первичного взаимодействия головных ударных волн, идущих от смежных элементов решетки, а все дальнейшее многообразие вторичных ударно-волновых взаимодействий позади решетки не воспроизводилось и не изучалось. В частности, отрывное донное течение и хвостовые скачки уплотнения в ближнем следе за элементами решетки не могли быть получены без учета вязкости газа. Соответственно, не воспроизводилась ячеистая конфигурация ударных волн в сверхзвуковом следе (рис. 2a).

Существенно иная картина наблюдается при M₁ = 2.4 (рис. 2б). Головные скачки уплотнения от соседних элементов решетки пересекаются по схеме маховского взаимодействия с образованием явно выраженной маховской ножки. При этом отраженные ударные волны непосредственно воздействуют на дозвуковую область ближнего следа, вызывая почти трехкратное увеличение ее размеров и значительное повышение донного давления на элементах решетки. Соответствующее значение коэффициента аэродинамического сопротивления элемента решетки C_x = 1.02 существенно меньше, чем C_x0 = 1.36 для одиночного цилиндра при M = 2.4.

В обоих случаях для значений M₁ на концах интервала варьирования течение взаимодействия сверхзвукового потока с решеткой (рис. 2) однозначно определяется величиной числа Маха M₁. Однако внутри интервала 2.4 < M₁ < 4.2 существуют области неоднозначной зависимости течений от этого параметра.

На рис. 3 представлены два различных режима обтекания решетки при одном и том же М₁ = 3.2, но при различной предыстории установления этого значения. Течение на рис. 3a получено по сценарию с постепенным уменьшением М₁, начиная от М₁ = 4.2. В главном оно аналогично стационарному течению, показанному на рис. 2a. Возмущения от соседних элементов решетки не влияют на ближнее поле течения. Поэтому все основные характеристики обтекания элементов решетки и одиночного цилиндра в неограниченном потоке совпадают. В частности, совпадают такие признаки, как конфигурация дозвуковой области ближнего следа и значение коэффициента аэродинамического сопротивления (в данном случае C_x = C_x0 = 1.31). В противоположность этому течение на рис. 3б,получено путем постепенного увеличения М₁, начиная от М₁ = 2.4. На этом пути происходит переход от стационарного к нестационарному режиму течения в аэродинамическом следе решетки. Вторичные скачки уплотнения от соседних элементов решетки взаимодействуют с неоднородной трансзвуковой областью ближнего следа, в результате чего возникают висячие отрывные структуры, которые, в свою очередь, взаимодействуют с дозвуковой частью ближнего следа. Это сложное взаимодействие имеет ярко выраженный нестационарный характер. Аэродинамическое сопротивление и подъемная сила элемента решетки совершают квазипериодические колебания. Среднее значение коэффициента сопротивления C_x = 1.25 значительно меньше величины C_x0 = 1.31 одиночного цилиндра.

Рисунок 4 иллюстрирует двойственность режимов обтекания решетки при М₁ = 2.6. Течение на рис. 4a получено путем постепенного увеличения числа Маха, начиная от М₁ = 2.4. Оно аналогично стационарному течению, показанному на рис. 2б. Отраженные скачки уплотнения в составе тройной маховской конфигурации при нерегулярном взаимодействии головных ударных волн от соседних элементов решетки непосредственно воздействуют на дозвуковую область ближнего следа, вызывая увеличение ее размеров и повышение донного давления на элементах решетки. При этом результирующее течение остается стационарным, а коэффициент аэродинамического сопротивления элемента решетки C_x = 1.03 намного меньше, чем значение C_x0 = 1.34 для одиночного цилиндра при M₁ = 2.6. В противоположность этому течение на рис. 4б получено путем постепенного уменьшения числа Маха. Оно имеет ярко выраженный нестационарный характер, при котором аэродинамическое сопротивление C_x и подъемная сила C_y элемента решетки совершают квазипериодические колебания. По оси абсцисс на графике колебаний C_y (рис. 4в) отложено i = 0.01t/Δt – количество сотен временных шагов Δt в ходе численного моделирования. При этом среднее значение C_x = 1.25 остается меньше коэффициента сопротивления C_x0 = 1.34 одиночного цилиндра.

Отношение C_x/C_x0 характеризует вклад эффектов аэродинамической интерференции в аэродинамическое сопротивление решетки. На рис. 5a построены расчетные данные, полученные по сценариям непрерывного увеличения и уменьшения числа Маха M₁ в диапазоне от 2.4 до 4.2. Зависимость C_x/C_x0 = f(M₁) распадается на три изолированных ветви. Нижняя ветвь I отвечает стационарным режимам течения как на рис. 2б и рис. 4a. На этом нижнем режиме аэродинамическое сопротивление элемента в составе решетки на 20–27% меньше, чем сопротивление одиночного элемента в неограниченном потоке. Верхняя ветвь III соответствует стационарным режимам течения как на рис. 2а и рис. 3a. На верхнем режиме аэродинамическое сопротивление элемента решетки и сопротивление такого же одиночного элемента в неограниченном потоке совпадают. Промежуточная ветвь II соответствует нестационарным режимам течения по типу как на рис. 3б и рис. 4б. Переходы с одной ветви на другую происходят пороговым образом так, как показано стрелками на рис. 5a. Очевидно, имеется две петли гистерезиса и, соответственно, два интервала значений параметра M₁, для которых решение задачи обтекания решетки неединственно. На интервалах неединственности отбор одного из возможных режимов течения определяется не только конкретным значением параметра M₁, но еще и предысторией установления этого значения. Это свойство двойственности решения является причиной гистерезисного поведения зависимости C_x/C_x0 от M₁.

Если интерпретировать решетку вместе с прилегающим слоем неоднородного диссипативного течения как поверхность сильного газодинамического разрыва, то коэффициент сопротивления единицы площади решетки C_x1 можно вычислить по теоретической формуле [9]

где M₁, M₂ – значения числа Маха до и после поверхности разрыва.

Расчетные газодинамические поля использовались для определения осредненных по ширине полосы периодичности (рис. 1) значений газодинамических параметров течения. Отношение коэффициента сопротивления проницаемой поверхности (1) к фактическому коэффициенту сопротивления решетки (1 – σ)C_x

характеризует точность модели проницаемой поверхности [9] применительно к решетке цилиндрических стержней. На рис. 5б построена расчетная зависимость K = f(M₁) для стационарных режимов течения I и III, показанных на рис. 5a. При вычислении C_x1 по формуле (1) использовалось осредненное расчетное значение числа Маха M₂ в поперечном сечении X = 20 в следе решетки, а при вычислении коэффициента сопротивления C_x элемента решетки использовалось фактическое расчетное давление на поверхности этого элемента. Незначительное отклонение значений параметра (2) от 1 характеризует высокую точность формулы (1) при стационарных режимах обтекания решетки. В случае нестационарных режимов II (рис. 5a) вопрос о моделировании решетки поверхностью газодинамического разрыва остается открытым.

На рис. 6 показаны поля числа Маха M при последовательном уменьшении значений M₁ с переходом от стационарного состояния течения на верхней ветви III кривой C_x/C_x0 = f(M₁) к нестационарному состоянию на промежуточной ветви II (рис. 5a)¹¹. Маховское взаимодействие вторичных скачков уплотнения в сверхзвуковой части следа при M₁ = 4.0 порождает локализованную дозвуковую область (рис. 6a). С уменьшением M₁ эта дозвуковая область увеличивается (рис. 6б). При этом наличие локального минимума числа Маха в ядре неоднородного сверхзвукового следа за элементом решетки порождает предвестник в виде вытянутого навстречу потоку локального выступа на фронте ударной волны, пересекающей ось этого следа (рис. 6б). По мере уменьшения M₁ предвестник быстро разрастается и проникает в область ближнего следа за элементом решетки (рис. 6в). В результате образуется длинный дозвуковой “язык”, взаимодействующий с отраженными и преломленными ударными волнами в примыкающих сверхзвуковых областях течения. Неустойчивость и колебания этой сложной ударно-волновой и отрывной структуры (рис. 6г) являются причинами формирования нестационарной системы вихревых сгустков, дрейфующих вниз по сверхзвуковому следу позади решетки (рис. 3б). Детали соответствующего автоколебательного процесса иллюстрирует рис. 7, где приблизительно на полупериоде автоколебаний визуализированы нестационарные ударно-волновые конфигурации и дозвуковые зоны (M < 1)²² для последовательных моментов времени с интервалом 500Δt (Δt – безразмерный шаг по времени в расчете).

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Три главных механизма определяют сложность взаимодействия сверхзвукового потока с редкой решеткой, не вызывающей запирания и возникновения отошедшей ударной волны.

Первый механизм – это воздействие локальных скачков уплотнения от соседних элементов решетки непосредственно на донную отрывную область ближнего следа за промежуточными элементами (по типу, как показано на рис. 2б), в результате чего происходит значительное увеличение размеров донной области и увеличение донного давления на элементах решетки. Связанная с этим механизмом двойственность структуры течения и соответствующий параметрический гистерезис возникают в тех случаях, когда указанные локальные скачки смещены вниз по потоку настолько, что проходят мимо короткой невозмущенной донной области, но способны попадать в зону предварительно возмущенной увеличенной донной области. В настоящей работе такой гистерезис идентифицирован по параметру M₁. В работе [12] аналогичный гистерезис был выявлен при варьировании расстояния между элементами решетки при фиксированном M₁ = 6. Очевидно, такой же механизм может приводить к гистерезису и по другим параметрам задачи, включая показатель адиабаты и число Рейнольдса. Этот же механизм может быть причиной параметрического гистерезиса характеристик сверхзвукового обтекания не только решеток, но и иных систем тел.

Второй механизм связан с существованием диапазона неоднозначности схем пересечения и отражения ударных волн от соседних элементов решетки. Изменение числа Маха M₁ приводит к изменению угла наклона и интенсивности скачков уплотнения на оси симметрии между элементами решетки. Существует известный диапазон этих параметров, при котором отражение ударных волн может происходить по схеме как регулярного, так и маховского взаимодействия [15, 16]. Определяемый этим механизмом гистерезис ударно-волновых конфигураций между элементами решетки описан в [11].

Третий механизм – образование предвестника в сверхзвуковом следе за элементом решетки при взаимодействии отраженных ударных волн с ядром неоднородности следа (подобно тому, как показано на рис. 6б,в). Наличие в сверхзвуковом потоке ядра поперечной неоднородности типа локального минимума числа Маха является известной причиной возникновения предвестника в виде выступающей вперед отрывной структуры в области пересечения ударной волны с ядром неоднородности [19]. Соответствующим нестационарным аналогом является возникновение предвестника перед фронтом ударной волны, движущейся вдоль локализованного канала пониженной плотности в покоящемся газе [20]. Проникание предвестника в область ближнего следа за элементом решетки приводит к образованию протяженного дозвукового канала (рис. 6в), через который происходит взаимодействие отрывной донной области с ударно-волновыми структурами в сверхзвуковой части следа. Неустойчивость такой конфигурации порождает автоколебания (рис. 6г, рис. 7), сопровождающиеся образованием дискретной дорожки вихревых сгустков, дрейфующих в сверхзвуковом следе (рис. 3б).

Одновременное действие второго и третьего механизмов может приводить к вторичным гистерезисным явлениям. Это связано с тем, что переход от регулярного к маховскому отражению происходит с внезапным образованием маховской ножки конечной высоты, что приводит к скачкообразному смещению (в сторону решетки) отраженных ударных волн, взаимодействующих с аэродинамическим следом. В частности, для решетки с проницаемостью 80% такой переход неизбежен в диапазоне уменьшения M₁ от 3.2 до 3.1 из-за достижения на этом интервале теоретического предела по углу поворота потока по схеме регулярного отражения ударных волн [16]. Представленные на рис. 6c регулярные отражения ударных волн при M₁ = 3.15 близки к указанному предельному состоянию. Дальнейшее уменьшение M₁ приводит к внезапному появлению тройных конфигураций ударных волн с маховской ножкой конечной высоты (рис. 6г). Однако возвратное движение по параметру M₁ в сторону увеличения не приводит к восстановлению схемы регулярного отражения при M₁ = 3.15. По мере роста M₁ происходит плавное уменьшение размеров маховской ножки вплоть до нуля. Детально этот процесс был исследован в [11] с помощью специальных мер по измельчению расчетных сеток в окрестности точки отражения.

Следует также отметить, что взаимодействие отраженных ударных волн с хвостовыми скачками и между собой может также происходить по схемам маховского или регулярного отражения, что, в свою очередь, может влиять на условия их взаимодействия с ядром неоднородного течения в следе. В конечном счете это может приводить к множественному параметрическому гистерезису при обтекании решеток и систем тел. Полученный в настоящих расчетах некоторый разброс точек на средней ветви II режимной кривой (рис. 5a) можно объяснить проявлением вторичных гистерезисов.

ЗАКЛЮЧЕНИЕ

На примере сверхзвукового обтекания бесконечной решетки цилиндрических стержней, имеющей проницаемость 80%, дана классификация режимов течения при непрерывном увеличении и уменьшении числа Маха набегающего потока воздуха в диапазоне от 2.4 до 4.2. Воздействие ударных волн от соседних элементов решетки на область ближнего аэродинамического следа за промежуточными элементами порождает неоднозначные структуры течений, а также нестационарные явления в следе. Идентифицированы три основных режима течений (нижний I, промежуточный II, верхний III, рис. 5a) и объяснены механизмы, определяющие многообразие соответствующих течений в следе, включая проявления параметрического гистерезиса при переходах от одного режима к другому. Интерпретация решетки как поверхности газодинамического разрыва дает функциональную зависимость между коэффициентом сопротивления и числами Маха осредненных сверхзвуковых потоков по обе стороны от проницаемого экрана. Результаты настоящих расчетов подтверждают высокую точность такого моделирования в случае стационарных режимов I и III (рис. 5б).

Работа выполнена при частичной поддержке РФФИ (проект 19-01-00242). Расчеты проводились на МВС-100К МСЦ РАН.