CC BY
32
Вычислительные технологии
Том 16, № 3, 2011
О численном исследовании автоколебательного обтекания затупленных конусов неоднородным потоком
В. И. Пинчуков
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: pinchvi@ict.nsc.ru
Исследуются неоднородные автоколебательные течения сжимаемого газа возле затупленных конусов. Двумерные уравнения Рейнольдса решаются неявной схемой Рунге — Кутты третьего порядка. Приводятся результаты численных исследований.
Ключевые слова: автоколебательные течения, схема Рунге — Кутты, методы высокого разрешения, уравнения Рейнольдса.
Введение
В настоящее время в исследованиях автоколебательных течений доминируют экспериментальные подходы. Математическое моделирование имеет вторичный характер и основывается на экспериментальных данных. Однако развитие численных методов и рост мощности ЭВМ позволяют повысить роль численных исследований в изучении и поиске новых нестационарных течений. Этот поиск осуществляется посредством выдвижения и проверки гипотез о механизмах автоколебаний. В данной работе на основе гипотезы, предложенной в [1, 2], проводится математическое моделирование обтекания затупленного конуса неоднородным потоком.
Для исследования автоколебательных течений используются как уравнения Эйлера, так и уравнения вязкого газа с дифференциальными и алгебраическими моделями турбулентности (см. обзор [3]). В [3] используется алгебраическая RNG-модель подсеточной вязкости, определяющая турбулентную вязкость по свертке тензора скоростей деформации с вычислением масштаба длины на основе размеров локальной сеточной ячейки. Таким образом, в турбулентной вязкости представлены лишь мелкомасштабные пульсации. Следует отметить, что в данном подходе при проверке сходимости решений по числу узлов сетки возникают проблемы, так как при изменении этого числа меняется и турбулентная вязкость.
В настоящей работе применяется родственный подход с вычислением алгебраической вязкости по формулам типа Себечи — Смита, в которых масштаб длины определяется по толщине локального сдвигового слоя, причем, поскольку алгебраические формулы дают приемлемые результаты лишь для пульсаций локального характера, этот масштаб ограничивается некоторой малой длиной. Естественно, данная модель, как и другие алгебраические модели турбулентности, позволяет лишь приближенно описывать течения (хотя и дифференциальные модели могут давать довольно большую ошибку [3]). Тем не менее она позволяет обнаружить автоколебательные решения
и проверить их сходимость по числу узлов сетки. Далее найденные решения могут быть исследованы на основе более совершенных моделей.
На основе численных исследований [1, 2, 4, 5] в [1, 2] сформулировано предположение о том, что генерация автоколебаний может происходить при взаимодействии возмущений, генерируемых контактным разрывом в результате неустойчивости Кельвина — Гельмгольца, с течением в окрестности точек пересечения двух ударных волн или ударной волны и контактного разрыва. Особая чуетвительноеть течений возле таких точек констатирована в [6]. Таким образом, контактный разрыв играет роль генератора возмущений, точка пересечения разрывов — роль усилителя, и если складываются условия для обратного влияния усиленных возмущений на течение в окрестности контактного разрыва, то могут возникнуть автоколебания потока.
Вероятность реализации этого механизма стала стимулом исследования натекания неоднородного потока на коническое тело. Оказалось: несмотря на то, что течение носит более сложный, чем предполагалось, характер, обнаружены автоколебательные режимы взаимодействия сверхзвукового потока, содержащего дозвуковую струю, с затупленным конусом,
1. Математическая модель
Следует отметить, что автоколебательные течения имеют определенную специфику, для учета которой при моделировании, кроме достаточно больших ресурсов ЭВМ, требуются надежные и эффективные численные методы. Так, большинство автоколебательных течений содержат контактные разрывы, неустойчивость которых существенно влияет на реализацию нестационарных режимов. Поэтому вычислительный алгоритм должен давать хорошее разрешение контактных разрывов, обеспечивающее моделирование их неустойчивости, для чего целесообразно применение методов высоких порядков, Кроме того, для надежного различения стационарных и автоколебательных режимов требуются хорошие свойства устойчивости. Наконец, опыт исследований показал, что при использовании уравнений Эйлера для моделирования автоколебаний возникают проблемы, так как с измельчением сетки возникают дополнительные мелкомасштабные возмущения, усиливающиеся в результате неустойчивости контактного разрыва и влияющие на поток, что приводит к отсутствию сходимости по числу узлов сетки,
В данной работе для численного интегрирования по времени определяющих уравнений применяется неявная консервативная схема Рунге — Кутты [7, 8] третьего порядка по времени и четвертого по пространственным переменным (вязкие слагаемые аппроксимируются со вторым порядком), В [8] показана ее абсолютная устойчивость для многомерного уравнения переноса, В схеме используются искусственные диффузионные члены с шаблоном из семи узлов по каждой переменной, В этих диффузионных слагаемых, так же как и в слагаемых, связанных с давлением и конвективным переносом, применена процедура адаптации. Формулы адаптации предназначены, во-первых, для демпфирования осцилляции возле разрывов, во-вторых, для повышения нелинейной устойчивости схемы, под которой подразумевается следующее. Формулы адаптации в случае модельного уравнения переноса в пределе асимптотически малого временного шага обеспечивают нелинейную энергетическую оценку при произвольно меняющихся коэффициентах адаптации (как известно, быстро меняющиеся коэффициенты разностных уравнений являются источником неустойчивости численных решений), В итоге
метод обладает хорошей сходимостью по времени, что в совокупности с высоким порядком позволяет исследовать течения, которые могут иметь как стационарные, так и автоколебательные режимы. Метод программно реализован для произвольных криволинейных координатных преобразований х = х(С, п), У = У (С, п) имеющих необходимую гладкость и отображающих единичный квадрат на плоскости переменных С, П в криво-
х, У
Для расчета вязкости используется алгебраическая турбулентная модель типа Себе-чи — Смита, Вычисление вязкости начинается с определения слоев смешения, что производится на основе расчета завихренности. Текущая точка принадлежит слою смешения, если тБ = Е(мгДгг) > еЕг | щ || Дг где т — завихренность, Б — площадь четырех элементарных ячеек разностной сетки, прилегающих к данному узлу сетки, Е — знак суммирования по сторонам ячеек, щ — вектор скорости в узлах сетки, Дг — вектор между узлами сетки, е — малая константа, для которой принято значение 0,01, Сумма слева приближенно представляет циркуляцию скорости по контуру четырехугольника. Для вычисления турбулентной вязкости внутри слоя смешения для вычисления турбулентной вязкости используется формула Прандтля ^ =| т | рг2, где р — плотность, г — масштаб длины, который рассчитывается по формуле Кармана
г = 0.4 Ь,
здесь Ь — расстояние до границы слоя смешения. Применяется также модифицированная формула
г = 0.4шт(Ь, д),
где д — параметр-ограничитель, позволяющий представить в турбулентной вязкости лишь коротковолновые пульсации, в то время как остальные описываются сеточным решением. При расчетах эта вязкость появляется в сдвиговых слоях и предназначена для подавления наиболее коротковолновых компонент решения, порождаемых в результате неустойчивости контактного разрыва,
В качестве граничных используются следующие условия:
а) на поверхности конуса — условие непротекания и экетраполяционные соотношения для давления, плотности и касательной компоненты скорости;
б) на границах, соответствующих оси симметрии, — равенство нулю радиальной компоненты скорости и экетраполяционные соотношения для других параметров;
в) на подветренной части внешней границы — экетраполяционные соотношения;
г) наконец, на наветренной части внешней границы задаются все параметры,
2. О возможной реализации автоколебаний вблизи конических тел
По-видимому, известные автоколебательные течения могут быть разделены на два семейства, Первое соответствует взаимодействию однородного потока е телом [1, 2, 4, 9, 10], второе — натеканию недораеширенной струи на тело [1, 3, 5, 11, 12], В обоих случаях экспериментальная реализация автоколебаний связана е созданием емкости с повышенным давлением, вытекание газа из которой при взаимодействии его с телом приводит к автоколебательному течению.
Рис. 1. Варианты относительного движения самолетов: a — без автоколебаний, б— автоколебания возможны
В работе рассматривается ситуация, по двум причинам не вписывающаяся в указанные случаи. Во-первых, рассматривается струя с давлением, не превышающим давления в окружающем потоке, а равным ему, во-вторых, натекающий поток имеет две ненулевые скорости, а не одну, как в этих семействах. Экспериментальное моделирование такой ситуации стандартными методами сопряжено с определенными проблемами.
Однако в принципе возможна натурная реализация данного случая, а именно, при попадании самолета в струю от реактивного двигателя другого самолета, но лишь при движении их в разных направлениях (рис. 1). Учитывая, что струя от двигателя достаточно быстро диссипируется и что самолету для попадания в струю на умеренном расстоянии необходимо проведение резкого маневра, реализация автоколебательного режима при обратном движении самолетов, по-видимому, не наблюдалась. Естественно, для установления вероятности нестационарных режимов при попадании самолета в струю двигателя другого, движущегося в противоположном направлении, самолета требуются дополнительные исследования и здесь констатируется лишь принципиальная возможность такого рода эффектов.
Определенные трудности, связанные с экспериментальными исследованиями и реализацией автоколебаний возле конусов, стимулируют применение для данных целей численных методов. Эти методы позволяют получать объемную, полную и относительно дешевую информацию и без дополнительных трудностей могут быть использованы для исследования нестационарного обтекания затупленных конусов.
3. Численные примеры автоколебательного обтекания конусов
Ввиду большого количества имеющихся в исследуемой задаче входных параметров ограничимся исследованием течений, которые могут реализоваться при попадании самолета в струю от реактивного двигателя другого самолета (при их противоположном движении). Иными словами, давления в струе и в основном потоке принимаем одинаковыми, температуру струи задаем выше температуры потока, скорость струи — меньше скорости потока (поскольку движение самолетов разнонаправленно) и меньше скорости звука (так как при сверхзвуковых скоростях струи и потока получаются лишь стационарные режимы, которые в данном случае не рассматривается).
Оказалось, что увеличение угла полураствора асоп затупленного конуса способствует переходу от стационарных режимов к автоколебательным. Рассмотрим течения возле конуса с углом полураствора асоп = 30°, для которого при указанных выше условиях все варианты оказались автоколебательными. На рис. 2, a приведено распределение
Рис. 2. Обтекание конуса: а — распределение плотности, б— векторы скорости
плотности для чисел Маха основного потока М81г = 2 и струи М^ = 0.9. Соотношение радиусов струи и затупления конуса принято равным ^е1/гсоп = 5/6, соотношение плотности струи и потока Р^/Рбъ = 0.8. Используется сетка 490 х 544, параметр ограничивающий размер представленных в турбулентной вязкости пульсаций, равен Гсоп/60, угол полураствора конуса асоп = 30о. На рис. 2, б для этого же варианта приведены векторы скорости в подобласти, примыкающей к затуплению. Они изображены в каждом десятом (по преобразованным переменным £ и п) узле сетки, т. е. в общей сложности в каждом сотом узле.
По векторам скорости видно, что вблизи затупления имеет место вихревая зона. В целом течение похоже на автоколебательное осесимметричное, реализующееся вблизи цилиндра с передним выступом (см., например, [1, 4, 9]). Априори предполагалось, что усилителем возмущений, генерируемых контактным разрывом, является точка пересечения головного скачка и границы между струей и основным потоком. Однако оказалось, что головная ударная волна в зоне, примыкающей к струе, трансформируется в веер волн сжатия. По-видимому, усилителем возмущений является точка пересечения контактного разрыва и второго скачка, располагающегося ниже циркуляционной зоны. Таким образом, механизм автоколебаний сходен с таковым в случае обтекания цилиндра с передним выступом.
Для исследования сходимости по количеству узлов вариант, отличающийся от рассмотренного выше лишь радиусом струи (он принят равным 0.666 от радиуса затупления, т^/гсоп = 2/3), рассчитан на сетках 327 х 363, 490 х 544, 653 х 725. Сравнение позволяет констатировать приемлемое согласие решений. Вместе с тем интенсивность и периодичность автоколебаний оказались более чувствительными к числу узлов. На рис. 3 для этого же варианта приведена динамика плотности в точке торможения на затуплении конуса для сеток 490 х 544 (ромбики) и 653 х 725 (квадратики). В норме С за счет локальных выбросов имеется различие, но в целом согласие удовлетворительное. Однако на сетке 327 х 363 результаты оказались неприемлемы. Поэтому далее обсуждаются расчетные результаты, полученные на сетке 490 х 544. Сетка 327 х 363 используется для получения начальных данных и оценки необходимой дистанции между наветренной границей расчетной области и поверхностью конуса, которая различается для разных вариантов более чем в два раза, сетка 653 х 725 применяется эпизодически для контроля сходимости.
Все рассмотренные варианты при дозвуковой скорости струи оказались "слабонестационарными", а именно, при расширении диапазона масштабов пульсаций, включенных в турбулентную вязкость (при увеличении параметра d в шесть раз до значения d = гсоп/10), течения становятся стационарными. Амплитуды автоколебаний могут существенно различаться. Например, при увеличении радиуса струи до г^ = 2гсоп колебания плотности в точке торможения значительно уменьшаются по амплитуде. При этом оказалось, что эпицентр автоколебаний смещается на боковую поверхность затупленного конуса. На рис. 4 показана динамика плотности в точке торможения при 0 < £ < 10 и 10 < £ < 20. На рис. 5 приведено распределение плотности в расчетной области. Входные параметры, кроме радиуса струи, соответствуют варианту на рис. 2 и указаны выше. Стрелкой обозначена точка, для которой изображена динамика плотности на рис. 4. Из рис. 5 видно, что циркуляционная зона сместилась вниз по потоку, веер волн сжатия, в который трансформируется головная ударная волна вблизи струи, значительно больше, чем на рис. 2, и на контактном разрыве ниже циркуляционной зоны имеют место достаточно выраженные колебания.
Расчеты проведены в диапазоне входных параметров 1/3 < г^/гсоп < 2, 0.5 < Р^/^т < 1.3 0.5 < М^ < 0.9. Некоторые параметры были фиксированны: = 2, ^соп = 30°. Во всех случаях структура течения имеет характерные черты: циркуляционная зона возле затупления или несколько ниже по потоку и веер волн сжатия, в который
р
1.25 1.00 0.75 0.50 0.25
2 4 6 8 1 0 12 14 16 18 1
Рис. 3. Динамика плотности для сеток Рис. 4. Динамика плотности на затуплении (при 490 х 544 и 653 х 725 0 < г < 10) и на боковой поверхности (при 10 <
г < 20) конуса
Рис. 5. Распределение плотности в расчетной области
трансформируется головная ударная волна вблизи струи. Хвостовой скачок может как присутствовать, так и исчезать, В последнем случае пропадают и автоколебания в потоке.
Таким образом, в работе обнаружены автоколебательные режимы обтекания затупленных конусов сверхзвуковыми потоками, содержащими дозвуковую изобарическую струю. Предложенный механизм автоколебаний объясняет случаи наличия или отсутствия нестационарных режимов.
Список литературы
[1] Пинчуков В.И. Численное моделирование нестационарных течений с переходными режимами // Журн. вычисл. математики и матем. физики. 2009. Т. 49, № 10. С. 1865-1873.
[2] Пинчуков В.И. Численное моделирование автоколебаний и поиск новых автоколебательных решений // Материалы VIII Междунар конф. NPNJ-2010. М.: МАИ-ПРИНТ, 2010. С. 226-229.
[3] Волков К.Н. Взаимодействие круглой турбулентной струи с плоской преградой // ПМТФ. 2007. Т. 48, № 1. С. 55-67.
[4] Пинчуков В.И. Численное исследование автоколебательных течений неявной схемой четвертого порядка // Вычисл. технологии. 2005. Т. 10, № 2. С. 114-126.
[5] Пинчуков В.И. Математическое моделирование колебаний, возникающих при втекании струи в полость // Там же. 2007. Т. 12, № 6. С. 73-80.
[6] Алексеев А.К. К оптимальному управлению переходом между регулярным и маховским взаимодействием ударных волн // Материалы VIII Междунар. конф. NPNJ-2010. М.: МАИ-ПРИНТ, 2010. С. 212-224.
[7] Пинчуков В.И. О численном моделировании нестационарных течений на больших интервалах по времени с использованием неявных схем высоких порядков // Матем. моделирование. 2004. Т. 16, № 8. С. 59-69.
[8] Пинчуков В.И. О численном решении уравнений вязкого газа неявной схемой Рунге — Кутты третьего порядка // Журн. вычисл. математики и матем. физики. 2002. Т. 42, № 6. С. 896-904.
[9] Антонов А.Н., Елизарова Т.Г., Павлов А.Н., Четверушкин А.Н. Математическое моделирование колебательных режимов при обтекании тела с иглой // Матем. моделирование. 1989. Т. 1, № 1. С. 13-23.
[10] Базыма Л.А. О влиянии вдува на стабилизацию течения при обтекании осесимметрич-ного тела с полостью, обращенной навстречу сверхзвуковому набегающему потоку // ПМТФ. 1998. Т. 39, № 4. С. 84-90.
[11] Адрианов А.Л., Безруков A.A., Гапоненко Ю.А. Численное исследование взаимодействия сверхзвуковой струи газа с плоской преградой // Там же. 2000. Т. 41, № 4. С. 106-111.
[12] Бочарова О.В., Лебедев М.Г. Моделирование нестационарного взаимодействия звуковой струи с преградой // Матем. моделирование. 2007. Т. 19, № 8. С. 31-36.
Поступила в редакцию 26 июля 2010 г., с доработки — 17 января 2010 г.
