Математическое моделирование колебаний, возникающих при втекании струи в полость Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 12, № 6, 2007

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ КОЛЕБАНИЙ, ВОЗНИКАЮЩИХ ПРИ ВТЕКАНИИ СТРУИ В ПОЛОСТЬ

В. И. Пинчуков

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: pinchvi@ict.nsc.ru

Unsteady gas flows in Hartman device are investigated. Two-dimensional Euler or Reynolds equations are solved using an implicit third order Runge—Kutta scheme. Numerical results reveal that flows in Hartman device are unsteady both in turbulent or laminar approach under the assumption that the flow leaks through the bottom of a cavity.

Введение

С целью исследования механизмов автоколебаний рассматривается влияние эффектов турбулентности и проницаемости стенок на существование квазипериодических режимов в устройстве Гартмана. В отличие от “слабо нестационарных” автоколебательных течений, имеющих место при взаимодействии солнечного ветра с межзвездным потоком и при натекании сверхзвукового потока на цилиндр с выдвинутым вперед фрагментом, в которых учет турбулентной диссипации позволяет получить стационарный режим, в “свистке Гартмана” имеет место “сильно нестационарное” течение, которое не удается стабилизировать ни учетом турбулентной диссипации, ни включением проницаемости дна полости устройства.

В исследованиях автоколебательных течений наиболее популярен экспериментальный подход [1-4]. Однако если, например, методы измерения давления в динамике на поверхности аэродинамических конструкций хорошо развиты, то визуализация нестационарных течений и получение о них детальной информации представляет собой трудную проблему. Преимущества математического подхода — в полноте создаваемой на его основе картины течения, а также в эффективности параметрических исследований близких по геометрии вариантов или исследований с использованием модификаций исходной модели.

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

коротковолновых компонент сеточных функций. По-видимому, именно наличие коротковолновых компонент численного решения обусловило использование в расчетах автоколебательных течений методов невысоких порядков [5-9], в основном первого порядка.

Как известно, схемы первого порядка дают хорошие результаты при наличии ударных волн, однако чрезмерно (значительно больше, чем ударные волны) размазывают протяженные контактные разрывы, для которых отсутствует эффект “градиентной катастрофы”, заключающийся в тенденции “укручения” волны сжатия. Поэтому, в частности, они плохо описывают течения с циркуляционными зонами.

В данной работе определяющие уравнения интегрируются на основе неявной схемы третьего порядка [10, 11]. Для более корректного описания течений с зонами коротковолновых колебаний применяется учет их в рамках алгебраической модели турбулентной вязкости. Численно моделируется осесимметричное течение, возникающее при втекании сверхзвуковой недорасширенной струи в цилиндрическую полость умеренной глубины. Изучается влияние эффектов турбулентности и проницаемости дна полости на интенсивность автоколебаний.

1. Численная модель

Решаются двумерные осесимметричные уравнения совершенного вязкого газа, записанные для случая произвольных криволинейных неортогональных координат. Вязкость вычисляется с использованием алгебраической турбулентной модели типа Себечи— Смита. Расчет вязкости начинается с определения слоев смешения, что производится на основе расчета завихренности. Текущая точка принадлежит слою смешения, если тБ = ^(иДг) > £^ |и||Дг|, где т — завихренность; Б — площадь четырех элементарных ячеек разностной сетки, прилегающих к данному узлу сетки; ^ — знак суммирования по сторонам ячеек; и — вектор скорости в узлах сетки; Дг — векторы, соединяющие узлы сетки; £ — малая константа, для которой принято значение 0.01. Сумма слева приближенно представляет циркуляцию скорости по контуру фигуры, состоящей из четырех ячеек. Внутри слоя смешения используется формула Прандтля ^ = р|т|г2, где р — плотность; г — масштаб длины, который рассчитывается по формуле Кармана

Ь — расстояние до границы слоя смешения. Используется также модифицированная формула:

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

Для численного интегрирования по времени определяющих уравнений используется неявная консервативная схема Рунге—Кутты [10] третьего порядка по времени и четвертого — по пространственным переменным. В [11] показана абсолютная устойчивость ее для многомерного уравнения переноса. В схеме используются искусственные диффузионные члены с шаблоном из семи узлов по каждой переменной. В этих диффузионных слагаемых, как и в слагаемых, связанных с давлением и конвекционным переносом, используется процедура адаптации. Формулы адаптации построены с

г = 0.4 Ь,

(1)

г = 0.4шіп(Ь, і),

(2)

учетом требования, чтобы для модельного уравнения переноса в пределе асимптотически малого временного шага они обеспечивали нелинейную энергетическую оценку при произвольно меняющихся коэффициентах адаптации (как известно, быстро меняющиеся коэффициенты разностных уравнений являются источником неустойчивости численного метода). Это свойство адаптационной процедуры (в совокупности с линейной устойчивостью при фиксированном шаблоне) обеспечивает хорошую сходимость по времени. Схема имеет приемлемое разрешение разрывов, хотя и несколько хуже, чем у ЕМО — схем высоких порядков (см. обзор [12]). Реализация схемы связана с обращением пятиточечных операторов с матричными коэффициентами посредством скалярных прогонок и локальных матричных преобразований. Метод программно реализован для произвольных криволинейных координатных преобразований х = х(С,п), У = У(С,п), имеющих необходимую гладкость и отображающих единичный квадрат на плоскости переменных £, п в криволинейный четырехугольник на плоскости физических переменных х,у. В рамках такого подхода получить хорошие сетки для сложных физических областей бывает затруднительно. Поэтому создана также версия основной программы для случая, когда функции х = х(С,п), У = у(С,п) отображают квадрат с вырезами (С < £0,п < По}, (С ^ СъП — П1} на криволинейный четырехугольник с двумя криволинейными же вырезами. Именно эта версия используется в расчетах, причем делается лишь один вырез (С < £0,П < П0}, соответствующий боковой стенке полости.

На рис. 1 схематически изображена расчетная область в плоскостях физических х, у и преобразованных С,П переменных. Точки, отмеченные одинаковыми буквами в физической и преобразованной плоскостях, отображаются друг в друга. Жирной линией отмечены границы, совпадающие с твердыми поверхностями. Криволинейная сетка в плоскости физических переменных получается отображением прямоугольной равномерной сетки на плоскости преобразованных переменных.

В качестве граничных условий используются:

1) на границах АВ, ЕЕ и частично СБ, соответствующих твердым поверхностям, — условие непротекания и экстраполяционные соотношения;

2) на дне полости может применяться соотношение

ир = (р — р0 )а, (3)

где р0 — давление за проницаемой поверхностью дна; а — коэффициент, моделирующий ее “проницаемость”;

Рис. 1. Расчетная область в физических и преобразованных координатах

3) на свободной границе ВС задается давление атмосферы при нормальных условиях и, если поток втекает внутрь области через эту границу, энтропия и касательная компонента скорости (она равна нулю). Другие параметры определяются из экстраполяционных соотношений;

4) на границе БЕ, соответствующей оси симметрии, используются экстраполяционные соотношения и равенство нулю радиальной компоненты скорости;

5) на части границы СБ, соответствующей втекающей струе, задаются коэффициент нерасчетности п, нормальная и касательная к границе компоненты скорости (касательная принимает нулевое значение) и температура, принятая равной значению температуры атмосферы при нормальных условиях.

2. Результаты расчетов

В качестве иллюстрации “слабо нестационарного” автоколебательного течения рассмотрим кратко один из вариантов решения задачи о взаимодействии солнечного ветра, т. е. сферически симметричного гиперзвукового течения от точечного источника, с межзвездным сверхзвуковым потоком, который полагается однородным и постоянным.

Используется модель идеального одноатомного газа с показателем адиабаты к = 5/3. Расчетная область и типичное распределение плотности в потоке изображены на рис. 2. В качестве единицы длины выбран радиус орбиты Земли. Внутренней границей является полуокружность радиуса 240, т. е. она расположена в 240 раз дальше Земли от Солнца. Поскольку параметры солнечного ветра задаются на расстоянии орбиты Земли, для их определения на внутренней границе производится решение алгебраических уравнений одномерной газовой динамики в предположении стационарности и сферической симметричности течения в окрестности Солнца. Более детально постановочная задача и метод решения изложены в [13].

В ламинарном приближении решение этой задачи получается нестационарным для различных параметров солнечного ветра и межзвездного потока, опубликованных в научной литературе, турбулентное приближение позволяет получить стационарное решение. Рассмотрим вариант с параметрами межзвездного потока Ц» = 2 км/с, п» = 0.1 в 1 см3, Т» = 100 К, параметрами солнечного ветра Це = 360 км/с, пе = 3 в 1 см3, Те = 140 000 К на расстоянии орбиты Земли. В рамках турбулентной модели получено стационарное решение. На рис. 2 изображено распределение плотности в логарифми-

Рис. 2. Стационарное распределение плотности. Турбулентный режим

Те = 140 000 К на расстоянии орбиты Земли. В рамках турбулентной модели получено стационарное решение. На рис. 2 изображено распределение плотности в логарифмическом масштабе (т. е. рисуется функция log1oр) для сетки 301 х 304. Имеются слабая внешняя ударная волна (слева у границы области), интенсивная внутренняя ударная волна и контактный разрыв между ними.

Далее полученное решение используется в качестве начальных данных при моделировании ламинарного течения. То есть используется формула (2), в которой положено d =10 (отметим, что расстояние от Солнца до правой границы расчетной области составляет 7000). Постепенно в решении вблизи контактного разрыва появляются возмущения, они возрастают, распространяются на более протяженные области. Со временем

д

Рис. 3. Распределения плотности в моменты времени с интервалом в 1740 лет. Ламинарный режим: а — £ = £0, б — £ = £0 + Д£, в — £ = £0 + 2Д£, г — £ = £0 + 3Д£, д — £ = £0 + 4Д£

устанавливается квазипериодический процесс, стадии которого изображены на рис. 3, содержащем распределения плотности в логарифмическом масштабе с интервалом по времени Д£ = 1740 лет.

Таким образом, причина автоколебательного режима, вероятно, — в неустойчивости контактного разрыва, так как стабилизация его турбулентной вязкостью (которая отлична от нуля лишь в сдвиговых зонах) приводит к получению стационарного режима. Отметим, аналогичное положение дел наблюдалось в наших расчетах натекания сверхзвукового потока на цилиндр с выдвинутым вперед фрагментом.

Ситуация с режимами течения в “свистке Гартмана” принципиально иная. Как включение турбулентной диссипации, так и задание граничных условий на дне полости, соответствующих проницаемости дна, приводит лишь к уменьшению амплитуды автоколебаний.

Рассмотрим результаты расчета втекания в полость струи с числом Маха Мс = 2, нерасчетностью п = 1.5, температурой, равной температуре атмосферы при нормальных условиях. Сетка содержит 391 х 211 узлов.

Радиус струи составляет 70 % от радиуса г цилиндрической полости, длина полости равна ее радиусу г. Срез струи расположен на расстоянии Ь от полости, которое вычисляется по формуле

Ь = г 1.72 Мс /к п, (4)

где к = 1.4 — показатель адиабаты. Эта эмпирическая формула дает длину первой бочки струи в отсутствие препятствий (см. [1]). То есть здесь расстояние от среза струи до полости принято равным длине первой бочки. Такой выбор обусловлен тем, что он обеспечивает существование автоколебаний для разных условий эксперимента [1]. В процессе расчета выводилось давление на дне полости в центре. Это давление в развитом течении носит колебательный характер. Рассмотрим структуру развитого ламинарного течения для случая дна без протекания, т. е. для а = 0. На рис. 4 и 5 приведены распределения плотности в моменты времени, соответствующие максимуму (рис. 4) и минимуму (рис. 5) давления на дне в центре его. На этих рисунках виден характер процесса в “свистке Гартмана”. Он заключается в периодическом значительном сжатии газа в полости и последующем отражении сжатой массы с образованием разреженной области вблизи дна.

1

Рис. 4. Втекание струи в полость. Распределе- Рис. 5. Втекание струи в полость. Распреде-

ние плотности в момент максимального дав- ление плотности в момент минимального давления на дне полости ления на дне полости

Рис. 6. Давление в точке Е Рис. 7. Втекание струи в полость с проницаемым

дном. Распределение плотности

Включение турбулентной диссипации приводит лишь к уменьшению интенсивности колебаний давления. Для иллюстрации степени влияния этого эффекта был проведен расчет турбулентного течения, далее решение использовалось в качестве начальных данных в расчете ламинарного течения. График давления на рис. 6 соответствует турбулентному (до момента времени £ = 24) и далее ламинарному течению в устройстве Гартмана.

Рассмотрим также результаты расчета втекания в полость струи с числом Маха Мс = 1.04, нерасчетностью п = 2.1. Расстояние от среза струи до полости снова вычисляется по формуле (4), остальные параметры не изменились. Сетка содержит 341 х 211 узлов. Рассматривается турбулентный режим с протеканием на дне полости, коэффициент а в формуле (3) принят равным 0.9, давление за дном полости р0 выбрано равным среднему по периоду давлению в точке Е, которое определяется в предварительном расчете. На рис. 7 показано распределение плотности в развитом течении. Следует отметить, что сжатая масса газа постоянно находится в полости, тем не менее колебания давления имеют место, однако значительно меньшей интенсивности, нежели в предыдущем варианте (см. рис. 6).

Таким образом, расчеты не выявляют единого механизма автоколебательных течений. Так, для взаимодействия солнечного ветра с межзвездным потоком и для натекания сверхзвукового потока на цилиндр с выдвинутым вперед фрагментом основной причиной нестационарности является неустойчивость контактного разрыва, поскольку неустойчивость исчезает при его стабилизации турбулентной вязкостью (напомним, она отлична от нуля лишь в сдвиговых зонах). В то же время при втекании струи в полость невозможность подавления нестационарности посредством даже двух факторов, указанных в работе, заставляет искать другие причины автоколебаний.

Список литературы

[1] ГлАЗНЕВ В.Н., КОРОБЕЙНИКОВ Ю.Г. Эффект Гартмана. Область существования и частоты колебаний // ПМТФ. 2001. Т. 42, № 4. С. 62-67.

[2] ГОРШКОВ Г.Ф., Усков В.Н. Особенности автоколебаний, возникающих при обтекании ограниченной преграды сверхзвуковой недорасширенной струей // ПМТФ. 1999. Т. 40, № 4. С. 143-149.

[3] Набережнова Г.В., Нестеров Ю.Н. Неустойчивое взаимодействие расширяющейся сверхзвуковой струи с преградой // Труды ЦАГИ. 1976. Вып. 1765.

[4] Антонов А.Н., Шалаев С.П. Экспериментальное исследование нестационарного течения в полостях, обтекаемых сверхзвуковым потоком // Изв. АН СССР. Механика жидкости и газа. 1979. № 5. С. 183-188.

[5] Базыма Л.А. Аэродинамика тел вращения с передней срывной зоной. Харьков, 1988. Деп. в ВИНИТИ. 3.10.1988, № 7262-В88.

[6] Базыма Л.А. О влиянии вдува на стабилизацию течения при обтекании осесимметричного тела с полостью, обращенной навстречу сверхзвуковому набегающему потоку // ПМТФ. 1998. Т. 39, № 4. С. 84-90.

[7] Антонов А.Н., Елизарова Т.Г., Павлов А.Н., Четверушкин А.Н. Математическое моделирование колебательных режимов при обтекании тела с иглой // Математическое моделирование. 1989. Т. 1, № 1. С. 13-23.

[8] Гапоненко Ю.А., Адрианов А.Л., Безруков А.А., Фаворский В.С. Расчетноэкспериментальное исследование взаимодействия сверхзвуковой струи газа с различными преградами // Мат. моделирование в механике: Сб. науч. тр. Красноярск: ВЦ СО РАН. 1997. С. 108-124. Деп. в ВИНИТИ 7.11.97, № 3357-В97.

[9] Адрианов А.Л., Безруков А.А., Гапоненко Ю.А. Численное исследование взаимодействия сверхзвуковой струи газа с плоской преградой // ПМТФ. 2000. Т. 41, № 4. С. 106-111.

[10] Пинчуков В.И. О численном моделировании нестационарных течений на больших интервалах по времени с использованием неявных схем высоких порядков // Матем. моделирование. 2004. Т. 16, № 8. С. 59-69.

[11] Пинчуков В.И. О численном решении уравнений вязкого газа неявной схемой Рунге— Кутты третьего порядка // ЖВМ и МФ. 2002, Т. 42, № 6. С. 896-904.

[12] Пинчуков В.И., Шу Ч.-В. Численные методы высоких порядков для задач аэрогидродинамики. Новосибирск: Изд-во СО РАН, 2000. 232 с.

[13] Пинчуков В.И. Численное исследование автоколебательных течений неявной схемой четвертого порядка // Вычисл. технологии. 2005. Т. 10, № 2. С. 114-126.

Поступила в редакцию 23 июля 2007 г., в переработанном виде —10 октября 2007 г.