CC BY
150
Том XXXV
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 200 4
№ 1 — 2
УДК 533.6.011.3/55 629.782.015.3
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННОГО ОБТЕКАНИЯ СВЕРХЗВУКОВЫХ ЛЕТАТЕЛЬНЫХ АППАРАТОВ И ИХ ЭЛЕМЕНТОВ НА ОСНОВЕ МНОГОЗОННОЙ ТЕХНОЛОГИИ
А. П. КОСЫХ, Г. Г. НЕРСЕСОВ, И. Ф. ЧЕЛЫШЕВА, В. Л. ЮМАШЕВ
Создана и апробирована новая программная система АРГОЛА-2 расчета трехмерных стационарных невязких течений около летательных аппаратов (ЛА) сложной формы на основе многозонной технологии с применением метода Годунова — Колгана и принципа установления по времени. В случае необходимости АРГОЛА-2 позволяет также проводить моделирование нестационарных газодинамических процессов. По сравнению с ранее разработанным программным комплексом АРГОЛА система АРГОЛА-2 существенно расширяет класс решаемых задач как по степени сложности компоновки летательного аппарата и топологии течения, так и по диапазону скоростей и углов атаки. При этом для гиперзвуковых скоростей учитываются не моделируемые в АДТ факторы — влияние теплофизических свойств воздуха.
Программная реализация созданной системы содержит ряд специальных технологических решений. К их числу относятся: использование динамических структур данных, управляемых с помощью дескрипторов; процедура автоматического сопряжения расчета в смежных областях; средства подачи заданий на расчет и др. Все это облегчает использование системы в исследовательском коллективе при проведении массовых расчетов.
Следует отметить, что в случае необходимости многозонная технология расчета естественным образом вписывается в многопроцессорную архитектуру ЭВМ.
В данной работе осуществлена верификация системы АРГОЛА-2 как путем сопоставления с экспериментальными результатами, так и с расчетными данными, полученными другими методами. По новой технологии проведено численное моделирование различных режимов обтекания элементов ЛА.
При решении разнообразных задач аэрогазодинамики на ЭВМ широкое распространение получили комплексы и пакеты прикладных программ (111111). рассматриваемые как средства исследования и автоматизации вычислительного эксперимента. В ЦАГИ в 80-е годы был создан ППП АРГОЛА для проведения расчетов сверх- и гиперзвукового обтекания летательных аппаратов (ЛА) [1]. [2]. За прошедшие годы эксплуатации этой программной системы получен ряд фундаментальных результатов в области аэрогазодинамики ЛА различной формы. Так. например. были проведены уникальные численные исследования обтекания и аэродинамических характеристик ЛА типа «Буран». «Бор». «Гермес». идеализированных компоновок гиперзвуковых летательных аппаратов (ГЛА) с воздушно-реактивным двигателем (ВРД) и их элементов.
Как показал опыт численного моделирования с помощью ППП АРГОЛА. проведение расчетов является достаточно трудоемким делом и требует большой подготовительной работы. В частности. необходимы: задание достаточно гладкой поверхности вычислительной модели ЛА; выделение головной ударной волны и разбиение расчетной области на несколько подобластей в соответствии с типом течения и используемыми численными методами (методы Бабенко — Русанова. Маккормака. Годунова — Колгана и др.); интерфейс между подобластями; графическая обработка полей течений и т. д. [3] — [6]. ППП АРГОЛА был ориентирован на
существовавший в то время парк ЭВМ сравнительно невысокой производительности. Вследствие этого в пакете программ для получения приемлемых численных решений на достаточно «грубых» сетках применялись методы повышенного порядка точности с выделением головной ударной волны как поверхности разрыва газодинамических функций и выделением подобластей расчета, где образуются местные дозвуковые зоны, внутренние скачки уплотнения и другие особенности в течении. При использовании такого рода технологии требовалось прерывать расчеты на ЭВМ и организовывать интерфейс на границах между смежными подобластями (сращивание решений). Следует отметить, что даже на «грубых» сетках расчет одного режима обтекания ЛА, например, орбитальных самолетов «Буран» или «Гермес», требовал значительных временных затрат (десятки часов процессорного времени).
За прошедшие после создания ППП АРГОЛА годы мощность ЭВМ значительно выросла, объем экспериментальных исследований из-за их дороговизны резко уменьшился, и увеличилась потребность в проведении расчетных исследований ЛА самой сложной формы в широком диапазоне изменения скоростей полета от дозвуковых до гиперзвуковых. В связи с изложенным выше возникла необходимость в разработке новой программной системы АРГОЛА-2, в которой реализована более универсальная технология математического моделирования обтекания ЛА с использованием нестационарных уравнений газодинамики. Этот программный комплекс основан на методе Годунова сквозного счета с кусочно-линейной аппроксимацией в элементарных ячейках, многозонной технологии разбиения расчетной области на ряд подобластей и принципе установления по времени для всего поля течения. Преимуществом многозонной технологии является гибкость, надежность, относительная простота управления процессом численного моделирования и возможность расчета топологически сложных пространственных течений. Расчет в каждой из подобластей (зон) проводится по единому вычислительному алгоритму, при этом для получения решения повышенной точности в некоторых из них необходимо использовать адаптивные сетки.
1. Многозонное представление расчетной области. Сложная структура течения около ЛА делает затруднительным, а чаще невозможным, построение единой регулярной сетки во всей расчетной области. Поэтому приходится либо вводить нерегулярную сетку, либо разбивать область течения на ряд подобластей с регулярной сеткой. Каждый из этих подходов имеет свои преимущества и недостатки. В настоящей работе выбрано разбиение области течения на подобласти и использование регулярной сетки. Регулярная сетка позволяет лучше контролировать качество аппроксимации, кроме того, такой подход соответствует опыту и традициям авторов.
Каждая отдельная область топологически эквивалентна кубу, но в рамках этого может быть как угодно деформированной и повернутой относительно базовой системы координат.
Опишем суть методологии создания многозонной сетки, покрывающей область возмущенного течения. По исходным данным, имеющим вид аналитических зависимостей или таблиц,
а также конструкторским чертежам строится математическая модель ЛА и его элементов. Для проведения расчетов обтекания в интерактивном режиме задается внешняя граница (оболочка) области возмущенного потока. Далее проводится разбиение области течения на непересекающиеся подобласти (зоны), топологически эквивалентные кубам, полностью заполняющим расчетную область. Затем каждая подобласть разрезается на элементарные ячейки, т. е. строится пространственная расчетная сетка.
42
г
х
Построение ведется в декартовой системе координат х, у, г . Три измерения сетки пронумерованы в произвольно выбранном порядке и обозначены на рисунке 41, 42, 4з . По направлениям 42, 4з сетка содержит N1, N2, N3 интервалов соответственно. Значения обобщенных криволинейных координат 4!, 42, 4з определяются просто как номера узлов вдоль соответствующих измерений и пробегают значения от 0 до N, N2, N3.
В реализованной версии системы АРГОЛА-2 дополнительно делается предположение о регулярном сопряжении областей друг с другом. В понятие регулярного сопряжения входят следующие условия:
всякая область с каждой из 6 сторон может граничить только с одной какой-либо областью. Другими словами, не допускаются никакие перехлесты и лоскуты, когда боковая граница одной области делится между несколькими другими областями;
узлы сеток двух областей на общей границе обязаны строго совпадать. Это предполагает, в частности, одинаковое количество интервалов узлов сетки в обеих областях по соответствующим измерениям.
Сформулированные условия облегчают практическую реализацию вычислительных алгоритмов, но несколько сужают свободу деления всей области течения на подобласти. Со временем предполагается освободиться от этих ограничений.
2. Модульная структура программы и динамические массивы данных. Программа расчета течений состоит из головного модуля и набора стандартных программ, выполняющих отдельные части расчета. Головной модуль обеспечивает повторение основного расчетного цикла по времени и обращается к другим модулям по мере необходимости в зависимости от числа областей и их взаимного расположения. Началу расчета предшествует подготовительная стадия. Программа получает формализованное задание на расчет с помощью специального языка, удобного и понятного для пользователя. Задание содержит все параметры расчета, а также ссылки на дополнительные источники информации: файлы с данными о расчетной сетке, начальном поле течения и т. д. На основе задания формируются все массивы и структуры данных, необходимые для расчета.
Для дискретного представления поля течения в программе должен быть предусмотрен обширный набор многомерных числовых массивов. Эти массивы содержат геометрические характеристики сетки — декартовы координаты узлов, объемы ячеек, площади и координаты векторов нормалей к граням ячеек; газодинамические параметры в ячейках сетки — декартовы компоненты вектора скорости, давление, плотность и другие характеристики, представляющие интерес. Потребные размеры числовых массивов могут изменяться в широких пределах в зависимости от конфигурации расчетной сетки в конкретных вариантах расчета. Поэтому крайне неудобно использовать для этих целей статические массивы, предоставляемые языками верхнего уровня и присущие традиционному стилю программирования. Ситуация еще больше осложняется при многозонном представлении поля течения, поскольку в этом случае требуется свой набор числовых массивов для каждой области.
Перечисленных трудностей удается избежать при использовании методики динамического распределения памяти. Согласно этой методике в головной программе определяется одномерный статический массив достаточной длины-рабочая память, в которой в дальнейшем размещаются все динамические массивы и структуры данных. Всякий динамический массив (или структура данных) характеризуется базой — индексом своего начала в рабочей памяти. Доступ к внутренним элементам динамического массива производится путем вычисления смещения — приращения индекса относительно базы.
Каждая подобласть поля течения характеризуется своими индивидуальными размерами расчетной сетки и набором динамических числовых массивов соответствующей длины. Указатели размеров сетки и базы всех массивов размещаются в особой числовой таблице, которая называется дескриптором области. Подобласть имеет свой дескриптор, полностью определяющий ее массивы данных, необходимые для проведения расчета.
Количество подобластей в поле течения будет разным в различных задачах. Очевидно, количество дескрипторов тоже будет разным. Естественно размещать дескрипторы в динамической памяти, несмотря на то, что все они имеют один и тот же фиксированный размер. В головной программе каждый дескриптор представлен своей базой в динамической памяти, которая является простой скалярной переменной и служит для ссылки на дескриптор и далее на всю совокупность данных соответствующей области течения. Тем самым структура головной программы предельно упрощается.
3. Физическая и численная модель течения. Рассматривается течение невязкого сжимаемого газа, численное моделирование которого проводится с использованием нестационарного метода сквозного счета Годунова [3] и его уточняющих модификаций [7], [8]. В общем случае термодинамическая модель газа не конкретизируется. Это может быть совершенный газ с заданным значением показателя адиабаты либо воздух в состоянии термодинамического равновесия. Наконец, может рассматриваться течение с неравновесными физико-химическими процессами. Стационарное поле течения рассчитывается установлением по времени путем задания некоторого начального распределения параметров потока газа и расчета его эволюции к конечному стационарному состоянию. Принцип установления по времени позволяет единым образом рассчитывать течения в широком диапазоне изменения скоростей потока от дозвуковых до гиперзвуковых.
При классическом подходе для невязкого течения обычно записывается система дифференциальных уравнений Эйлера, формулируются начальные и граничные условия и далее поставленная математическая задача решается тем или иным методом. Однако в случае метода Годунова данный формальный этап является излишним, поскольку уравнения Эйлера в методе Годунова нигде в явном виде не фигурируют. Вместо этого проводится прямое дискретное моделирование газодинамического процесса обтекания на основе интегральных законов сохранения.
Метод Годунова можно рассматривать как частный случай метода конечных объемов. Вся область течения разбивается на большое число элементарных ячеек, для которых рассматриваются законы сохранения массы, импульса и энергии. Решаются две отдельные задачи:
для каждой пары соседних ячеек определяются потоки массы, импульса и энергии из одной ячейки в другую в зависимости от текущих значений параметров газа в ячейках с использованием модели распада разрыва;
для каждый ячейки определяется изменение содержащихся в ней массы, импульса и энергии за малые промежутки времени At в зависимости от потоков из соседних ячеек.
Согласно сказанному выше течение рассматривается на регулярной сетке в подобластях, топологически эквивалентных кубам. Очевидно, каждая ячейка сетки также оказывается топологически эквивалентной кубу и имеет шесть боковых граней, по которым она граничит с соседними ячейками. Для определенности присвоим граням ячейки номера от 1 до 6 в порядке возрастания индекса по обобщенным сеточным координатам д2, д3. При этом номера 1 и 2 присваиваются граням, поперечным измерению дь номера 3 и 4 — измерению д2, номера 5 и 6 — поперечным д3, как показано на рисунке:
Пусть V — объем ячейки, £ — площадь боковой грани, р — давление, р — плотность, Ух, Уу, у2 — декартовы компоненты вектора скорости, е — внутренняя энергия, И — энтальпия, к — кинетическая энергия газа. Обозначим через В — вектор массы, импульса и энергии газа в ячейке, ¥ — вектор потоков массы, импульса и энергии через боковую грань:
В = V
р
РУх
РУу
Рvz р(е + к)
¥ = £
О
Оух + РПх
вУу + рПу
Оу2 + рп2 О(к + к)
Здесь О — поток массы через боковую грань; пх, Пу, п2 — компоненты вектора нормали к боковой грани. Условимся считать положительным направлением потоков и нормалей направление возрастания индекса по обобщенной сеточной координате. Тогда законы сохранения массы, импульса и энергии газа в ячейке принимают вид:
ДВ=Д?(¥1 - ¥2 + ¥3 - ¥4 + ¥5 - ¥6).
Здесь ДВ — изменение параметров в ячейке за промежуток времени Д?, индексы потоков соответствуют номерам граней ячейки. Это основная расчетная формула метода Годунова на регулярной сетке.
Аппроксимационные свойства исходной схемы Годунова со ступенчатым распределением параметров газа в ячейке хорошо изучены. Она относится к классу явных схем, является условно устойчивой в соответствии с критерием Куранта — Фридрихса — Леви и имеет первый порядок точности. Кусочно-постоянная аппроксимация обеспечивает безотказность метода Годунова, однако вносит в численное решение заметную схемную вязкость. Поэтому появились модификации схемы Годунова, повышающие точность аппроксимации приблизительно до второго порядка. Это схема Колгана [7] для одномерного случая и ее обобщение на многомерный случай [8]. Отличительная особенность модификации метода Годунова заключается в способе вычисления потоков ¥ через грани ячеек. В этих модификациях также используется модель распада разрыва, однако по обе стороны разрыва берутся не кусочно-постоянные в ячейке распределения параметров газа, а кусочно-линейные с учетом местных градиентов. Существуют также зарубежные аналоги этих методов, известные как схемы TVD.
В общем случае между ячейками имеет место разрыв — скачкообразное изменение газодинамических параметров. Предполагается, что этот разрыв распадается с течением времени в соответствии с законами газовой динамики. Хорошо известная задача о распаде произвольного разрыва является автомодельной и имеет точное решение, из которого находятся местные значения параметров газа на границе между ячейками. Именно эти значения, не зависящие от времени в силу автомодельности, используются для вычисления потоков массы, импульса и энергии по приведенным выше формулам. Наряду с точным решением задачи о распаде произвольного разрыва нередко используются различные упрощенные варианты решения: использование линейного приближения, методика [9] и т. д. Применение приближенных методик дает существенный выигрыш в вычислительном быстродействии и обычно слабо влияет на качество решения.
Задача о распаде произвольного разрыва обычно формулируется и решается в одномерной постановке. В случае пространственного течения между соседними ячейками имеет место разрыв векторных величин и необходим переход к локальной одномерной задаче о распаде разрыва. Рассмотрим схематично две смежные ячейки 1 и 2:
Ч\
Здесь ^1 и ^2 — векторы скорости газа в ячейках, п — единичный вектор нормали к общей грани между ячейками. Как было принято ранее, этот вектор направлен в сторону возрастания сеточного индекса и определяет локальное положительное направление отсчета. Вектор скорости с каждой стороны от границы может быть разложен на нормальную и тангенциальную составляющие
Vп =(v, п )• п , ^= V - ^п .
Величина (V ,П) представляет собой проекцию вектора скорости на нормальное направление и участвует в качестве скорости газа в одномерной задаче о распаде произвольного разрыва. Из решения этой задачи находится результирующая скорость vR по нормали к границе раздела. Поскольку тангенциальная составляющая вектора скорости не изменяется в процессе распада разрыва, результирующий вектор скорости на границе будет равен
^ • п + ух .
Разумеется, тангенциальную составляющую надо брать либо из первой, либо из второй ячейки в соответствии с направлением движения контактной границы при распаде разрыва. Аналогичная операция проводится с плотностью, температурой и другими параметрами физической природы.
4. Методика учета реальных свойств газа. Из анализа публикаций по вычислительной аэродинамике нетрудно убедиться, что большинство аэродинамических расчетов в мире выполняются для модели совершенного газа. В этом случае действуют простые соотношения между термодинамическими параметрами среды, которые изначально закладываются как в исходные уравнения, так и в последующую программную реализацию. В результате становится проблематичным обобщение вычислительной технологии на случай реальных термодинамических свойств среды. Типичным примером такой ситуации могут служить популярные приближенные соотношения для задачи о распаде произвольного разрыва [9].
В противоположность этому при разработке нашей многозонной вычислительной системы не делалось никаких предположений о конкретном виде термодинамических соотношений. Принято, что термодинамические свойства среды определяются функциями общего вида e = e(p, р), С = ^, р).
Здесь е — внутренняя энергия, с — скорость звука. В используемой технологии аэродинамического расчета переменные p и р считаются первичными переменными поля
течения. Исходя из этого, термодинамические функции e(p, р), c(p, р) позволяют находить недостающие параметры и вычислять потоки массы, импульса и энергии между ячейками расчетной сетки. Однако далее на каждом шаге интегрирования законов сохранения в ячейках получаются новые значения внутренней энергии, по которым надо восстанавливать соответствующие значения давления. Для этого требуется обратная термодинамическая функция Р = Р(Р, e).
В системе АРГОЛА-2 обе «прямые» термодинамические функции реализованы в виде одного программного модуля, где давление р и плотность р служат входными параметрами, а
внутренняя энергия e = e(p, р) и скорость звука c = ^, р) — выходными. Все параметры определяются в системе СИ: давление — в Па (Н/м ), плотность — кг/м , внутренняя энергия — дж/кг, скорость звука — м/с. Внутреннее содержание модуля термодинамических функций может варьироваться в зависимости от рассматриваемой среды. Выбор термодинамической модели среды для конкретного варианта расчета производится в формализованном задании на расчет наравне с прочими параметрами.
В настоящее время реализованы три варианта модели среды:
совершенный газ с заданными значениями показателя адиабаты и молекулярного веса;
равновесно-диссоциирующий воздух, состояние которого описывается при помощи числовых таблиц на основе аппроксимационных зависимостей [10];
расчет течений с учетом неравновесных физико—химических превращений в газе [11].
Технология расчета неравновесных процессов позволяет варьировать состав среды и модель химических реакций в широких пределах, причем, это делается без внесения изменений в программу.
Заметим, что общность описания термодинамических зависимостей достигается ценой увеличения вычислительных затрат. Эти затраты связаны с передачей формальных параметров, дополнительными уровнями вложенности модулей, итерационным процессом (который автоматически осуществляется даже в случае совершенного газа). Однако здесь авторы сознательно пошли на снижение чисто вычислительной эффективности программы, поскольку, с другой стороны, упрощается эксплуатация и модификация программы, а в итоге достигается выигрыш с точки зрения эргономики всего цикла расчетных исследований. В принципе существует возможность сократить вычислительные затраты, связанные с термодинамикой, если применить более сложную, но эффективную и универсальную методику представления термодинамических моделей [12].
5. Результаты численного моделирования обтекания элементов летательных аппаратов. Вычислительная система АРГОЛА-2 создавалась, прежде всего, для математического моделирования стационарных газодинамических процессов в широком диапазоне изменения определяющих параметров — числа Мш, угла атаки а и угла скольжения р. В частности, число Мш может варьироваться от дозвуковых до сверхзвуковых значений. При гиперзвуковых скоростях полета учитывается влияние на аэродинамику ЛА и его элементов реальных термодинамических свойств газа. Новая вычислительная система всесторонне верифицируется и находит применение в фундаментальных и прикладных исследованиях. С ее помощью проведено численное моделирование обтекания разнообразных аэродинамических форм: носовые затупления, удлиненные тела вращения, эллиптические конусы,
схематизированные аэрокосмические ЛА и т. д.
Для получения решения на основе принципа установления по времени в качестве начального поля во всех ячейках расчетной области задавались параметры невозмущенного набегающего потока (px, рот vxx, vyx, vzx ). На теле выполнялось условие непротекания. На
внешней границе всей расчетной области в зависимости от режима обтекания выполнялись асимптотические условия — задавались параметры невозмущенного потока либо рассчитанные на предыдущем шаге по времени параметры потока экстраполировались на границу области.
Возможности созданной системы АРГОЛА-2 продемонстрированы при решении ряда задач. Были проведены расчеты обтекания совершенным газом (у = 1,4) модели «затупленный конус — цилиндр» под нулевым углом атаки и числах Мш = 0,9 и 1,5. Данная форма обтекаемого тела выбрана в соответствии с экспериментом, проведенным в аэродинамической трубе Т-108 ЦАГИ [13]. Полуугол раствора конуса равен 18°, радиус затупления ~ 20 мм, радиус цилиндрической части ~ 60 мм, длина конической части ~ 140 мм, длина модели ~ 380 мм. Расчетные сетки для до- и сверхзвукового обтекания показаны на рис. 1, где для наглядности изображения число ячеек уменьшено в несколько раз.
Количество подобластей в расчетах варьировалось. Так, например, для случая Мш = 0,9 расчетная область состояла из 7 подобластей с общим числом ячеек сетки ~2-104. При необходимости ячейки сетки сгущались вблизи тела и в области больших градиентов газодинамических параметров. Расчетная область для Мш = 1,5 включала в себя 5 подобластей. На рис. 2, 3 представлены некоторые результаты проведенного численного моделирования. Поля изобар p = const [Па] показаны на рис. 2, где видна сложная топология течений, соответствующих различным режимам обтекания. Вблизи излома образующей тела формируется волна разрежения, при этом в случае Мш = 0,9 в сверхзвуковой зоне формируется замыкающий скачок уплотнения.
У
а)
Рис. 2. Поля изобар течения для модели «затупленный конус-цилиндр»:
а) М„ = 0,9, а = 0; б) М„ = 1,5, а = 0
На рис. 3 приведены распределения коэффициента давления ср вдоль тела (здесь и далее буквой Ь обозначается длина модели). Сплошной линией изображены результаты моделирования по модифицированному методу Годунова сквозного счета (АРГОЛА-2), а пунктирной линией показаны расчетные данные, полученные с помощью маршевого метода Бабенко — Русанова (АРГОЛА) с выделением головной ударной волны [14]. Здесь же для сравнения приведены результаты эксперимента [13]. Видно согласование результатов настоящего расчета при Мш = 1,5 с экспериментом и данными, полученными по маршевому методу 2-го порядка точности. Отметим, что отличия в расчетных значениях коэффициента сопротивления сх, полученных по программам АРГОЛА-2 и АРГОЛА, невелики и не превосходили ~3%. В работе [13] аэродинамические характеристики рассмотренной модели не исследовались.
Была проведена серия расчетов обтекания модели затупленного по сфере тонкого кругового конуса с полууглом при вершине 5° (радиус затупления — 2,5 мм, длина модели — 545 мм) при различных числах Мш и углах атаки. Расчетная область содержала 7 подобластей и около 105 ячеек. Результаты сравнивались с данными расчетов по методу Бабенко — Русанова (АРГОЛА) и с экспериментом [15]. На рис. 4 показаны распределения давления в продольном сечении (2 = 0) на наветренной и подветренной поверхностях конуса при Мш = 4 и а = 5°. Отличие результатов расчета по программе АРГОЛА-2 от «эталонного» решения по методу 2-го порядка точности не превосходит 3 — 4%. Проведенное сравнение аэродинамических характеристик свидетельствует
а)
1,5
1
0.5
О
-0,5
-1
*1,5
!
1
А •
1) 0 1 0 2 0 3 0 4 0 5 • 0 6 0 7 0 8 0
Уш
1
£
б)
, 9 •- —
0 0 1 0 2 0 3 ^ ■ ■у *4? 8 0 1- •9
х
Ь
0,6
0,8
X 11
Рис. 3. Распределение коэффициента давления су вдоль модели «затупленный конус — цилиндр»:
а) М„ = 0,9, а = 0; б) М„ = 1,5, а = 0;----АРГОЛА-2, метод Годунова;
......— АРГОЛА, метод Бабенко — Русанова; • • • — эксперимент, Т-108
ЦАГИ
об удовлетворительном согласовании результатов расчетов коэффициента подьемной силы су, полученных по различным вычислительным технологиям и в эксперименте. Однако отличие между расчетным и экспериментальным значениями коэффициента сопротивления сх достигало 10%, что объясняется, по-видимому, неучетом в расчетном исследовании эффектов вязкости для рассматриваемого удлиненного тела.
piПа] 250000--
200000.
0 ' fte ' гЬо ' ifeo ' Йо 1 5&0 ^*
Рис. 4. Распределение давления по поверхности затупленного конуса в плоскости симметрии
течения, 0к = 5°, М = 4, а = 5°
1 — АРГОЛА; 2 — АРГОЛА-2; 3 — наветренная сторона; 4 — подветренная сторона
На рис. 5 с помощью линий M = const изображена картина более сложного течения около затупленного эллиптического конуса с отношением полуосей 1:2 в поперечном сечении (y, z). Носовая часть представляла собой трехосный эллипсоид с отношением полуосей 0,25:1:2 в системе координат x, y, z. Вся расчетная область разбита на 15 подобластей с общим числом ячеек ~ 4-104.
Наряду с простыми аналитически заданными телами были рассмотрены также тела, близкие по форме к реальным объектам. Было проведено исследование аэрогазодинамики обтекания передней части схематизированного ЛА. Моделирование осуществлено по единому численному алгоритму (АРГОЛА-2) в широком диапазоне изменения входных параметров, в частности, число Мш варьировалось от дозвуковых до гиперзвуковых значений. Рассматривались режимы обтекания, включая переход через скорость звука, в том числе околозвуковые и закритические течения.
Для проведения численного моделирования конструировалась пространственная сетка, охватывающая возмущенное поле течения и содержащая до 1 млн. ячеек. Топология расчетной сетки для сверхзвуковых режимов обтекания (Мш > 1) показана на рис. 6. При этом для наглядности графического изображения перед построением сетка была существенно укрупнена по всем трем направлениям. Для решаемой задачи расчетная область разрезалась на 10 подобластей, форма и размеры которых выбирались исходя из геометрических особенностей модели ЛА и реализуемых режимов трехмерных течений. В случае Мш < 1 расчетная область содержала дополнительно 4 подобласти за кормовым срезом.
Структура течений детально изучалась путем построения изолиний газодинамических функций в различных сечениях возмущенного потока, а также построением распределений этих функций по поверхности аппарата. В расчетах и на рисунках линейные размеры взяты в миллиметрах, а газодинамические параметры — давление р, плотность р, компоненты скорости vx, vy, vz отнесены соответственно к рж рот у]р^/р^, . Возможности численного моделирования
разработанной программной системы при до- и сверхзвуковых скоростях полета иллюстрируют рис. 7, 8. Используемая единая универсальная технология расчета позволяет корректно
моделировать до- и трансзвуковые течения при обязательном выполнении основного условия, чтобы искомое решение задачи обтекания не зависело от формы и положения внешней границы.
Рис. 5. Поле течения около эллиптического конуса, М = 0,8, а = 16°, 5э = 0,5:
а) линии М = const в кормовом сечении; б) линии M = const в продольном сечении
На рис. 7 показана структура течения в плоскости симметрии для трансзвукового режима обтекания (Мш = 0,8) с местной сверхзвуковой зоной и замыкающим скачком уплотнения у подветренной поверхности тела. Рис. 8 соответствует другому режиму обтекания (Мш = 2), когда образуется головная ударная волна с дозвуковым течением около носового затупления и сверхзвуковым течением далее вниз по потоку. В окрестности изломов образующих тела формируются волны разрежения с последующим торможением потока без образования скачков уплотнения. Как показали расчеты, при сверхзвуковых числах Маха у поверхности рассматриваемого тела всегда формируется энтропийный слой газа, который может утоньшаться на наветренной поверхности (а > 0).
Рис. 6. Пространственная расчетная сетка для сверхзвукового обтекания передней части схематизированного ЛА
В последние годы пристальное внимание многих исследователей за рубежом сосредоточено на проблеме разработки гиперзвукового летательного аппарата интегральной компоновки с ВРД. Создание такого рода новой техники связано с решением ряда фундаментальных задач аэротермодинамики. При этом используется комплексный подход к решению поставленной задачи, в котором одновременно находят применение теоретические, экспериментальные и численные методы.
С помощью программной системы АРГОЛА-2 проведено расчетное исследование схематизированной модели ГЛА, в том числе с целью оптимизации характеристик потока на входе в воздухозаборник ВРД. На рис. 9 — 14 представлены результаты численного
моделирования пространственного течения около передней части ГЛА — клиновидных тормозящих поверхностей и течения внутри профилированного канала. Расчеты были выполнены для набегающего потока с числом Мш = 6 и углов атаки а = 0
и 10°.
Общий вид течения для одного из вариантов расчета показан на рис. 9, 10 в виде изобар (p/pa = const) в плоскости симметрии и векторов скорости на поверхности модели. По поведению изолиний можно наблюдать сложную картину интерференции скачков уплотнения и волн разрежения. На рисунках также изображены поля направлений линий тока внутри профилированного канала. Особенности стока с боковых кромок обусловлены прямоугольной формой поперечного сечения модели, при этом боковая поверхность модели оказывается в области разрежения.
Рис. 7. Изобары и распределение давления в плоскости симметрии, М = 0,8, а = 10°:
------на нижней поверхности;.......— на верхней поверхности
У
-°іооо ' о1 ' 1 Ьоо ' гЬоо ьЬоо ' ^Ьоо ' эооо
Рис. 8. Изобары и распределение давления в плоскости симметрии,
М = 2, а = 10°:
---------на нижней поверхности;.....— на верхней поверхности
Как показал анализ результатов расчетов, течение имеет существенно пространственный характер, что выражается в поперечной неоднородности газодинамических параметров. Для случая а = 10° по сравнению с а = 0 на участке 0 < х < 12 наблюдается простое количественное усиление неоднородности потока. Однако в канале (12 < х < 20) резко изменяется качественная картина течения. Следы падения скачков на верхнюю поверхность канала приобретают сложную форму, а векторы скорости показывают, что происходит поворот потока от боковой стенки к центру канала и возникает область возвратного течения (рис. 11, 12). Структура полей течения сви-
Рис. 11 Изолинии давления (а) и числа М (б) в продольных сечениях канала
Рис. 12 Изолинии давления (а), числа М (б), продольной скорости (в) и векторы скорости (г) на верхней стенке канала
Р РОС
во 60
6)
40
20 О
0.60 0.70 0.80 0,90 1.00
Рис. 13. Сравнение с экспериментом по структуре течения и распределению
давления при а = 0:
а) распределение давления вдоль поверхности А; б) распределение давления вдоль поверхности В; ------------- — расчет; • • • • — эксперимент
детельствует о наличии здесь обширной дозвуковой зоны. На рис. 13, 14 дано сравнение расчетных и экспериментальных распределений давления по поверхности ГЛА вне и внутри канала ВРД (Мш = 6, а = 0 и 10°, Яеш ~ 3-106). Экспериментальное исследование проведено в аэродинамической трубе Т-121 ЦАГИ коллективом авторов под руководством В. Н. Гусева.
Как следует из проведенного сравнения, до входа в воздухозаборник расчетные и экспериментальные распределения давления удовлетворительно согласуются между собой. Однако внутри профилированного канала топологии невязкого и вязкого течения в деталях заметно различаются. Рассмотрим структуру течения газа в канале у поверхности А при а = 0. При невязком обтекании на входе в канал (х > 0,6) образуется центрированная волна разрежения, тогда как в эксперименте реализуется иная картина течения. По-видимому, скачок от обечайки достигает поверхности А и происходит отрыв пограничного слоя от этой поверхности (0,6 <х < 0,7) с образованием отрывной зоны с повышенным давлением.
Далее вниз по течению (0,7 <х < 0,8) в расчете и в эксперименте наблюдается течение расширения с последующим торможением потока и качественное согласование распределения давления. На участке 0,8 <х < 1 видно не только количественное, но и качественное расхождение в поведении давления. При численном моделировании формируется скачок уплотнения от излома поверхности, а в эксперименте из-за очередного отрыва потока образуется изобарическая область небольших размеров, и давление в вязком потоке на рассматриваемом участке изменения х повышается непрерывным образом. Существенное различие в области выходного сечения гондолы можно объяснить неполным соответствием условий на свободной границе в задаче невязкого обтекания условиям трубного эксперимента.
Р Р* а)
•
•
. 1—•
г • —1 и
0,00 0.20 0,40 0,60 0.80 1,00
Р Рх
150
100
б)
50 •
1 •
х/Ъ
0,60 0,70 0,80 0,90 1,00
Рис. 14. Сравнение с экспериментом по структуре течения и распределению давления при а = 10°: а) распределение давления вдоль поверхности А; б) распределение давления вдоль поверхности В; ----------- — расчет; • • • • — эксперимент
На внутренней поверхности В обечайки на участке 0,6 < х < 0,75 получено качественное соответствие в распределениях давления. Ниже по течению (0,75 <х < 1) заметное расхождение в значениях давления объясняется тем, что в расчете отраженный от поверхности А скачок из-за его ослабления в волне разрежения не достигает обечайки, тогда как в эксперименте отраженный интенсивный скачок приходит на обечайку и соответственно повышает уровень давления в области (0,75 <х < 1).
На основании накопленного опыта в заключение отметим, что в широких диапазонах изменения входных параметров численное моделирование в рамках уравнений Эйлера позволяет с приемлемой для практики точностью определять распределенные и интегральные силовые нагрузки на ЛА (элементы ЛА), а также характеристики возмущенного потока, если влияние эффектов вязкости невелико. В противном случае получаемые численные решения можно рассматривать как некоторый газодинамический фон, на котором разыгрываются сложные реальные процессы обтекания.
Работа выполнена при финансовой поддержке РФФИ, грант № 01-01-00633.
ЛИТЕРАТУРА
1. Михайлов Ю. Я., Юмашев В. Л. Комплекс АРГОЛА: Автоматизированный расчет гиперзвукового обтекания летательного аппарата/Материалы VII Всесоюзного семинара по комплексам программ математической физики. — Новосибирск. — 1982.
2. Михайлов Ю. Я., Савин И. В., Челышева И. Ф., Юмашев В. Л. Комплекс APГОЛA: Aвтоматизированный расчет гиперзвукового обтекания летательного аппарата // Труды ЦДГИ. — 1993. Вып. 2478.
3. Годунов С. К., Забродин A. В., Иванов М. Я., Крайко A. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики. — М.: Наука. — 197б.
4. Базжин A. П., Пирогова С. В. Aлгоритм расчета трехмерных смешанных течений газа//Труды ЦДГИ.— 1974. Вып. 1604.
5. Базжин A. П., Михайлов Ю. Я., Нерсесов Г. Г. Специализированный алгоритм расчета сверхзвуковых трехмерных течений газа//Труды ^AFH.— 1984. Вып. 2248.
6. Голубинский A. A., Косых A. П., Савин И. В., Челышева И. Ф.
Численное моделирование сверхзвукового пространственного обтекания идеализированных компоновок ВКС совершенным газом и равновесно-диссоциирующим воздухом/Сб. докладов научной Школы-семинара «Механика жидкости и газа», ч. II.— 1992.
7. Колган В. П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики//Ученые записки ^AFH.— 1972. Т. III, № 6.
8. Тилляева Н. И. Обобщение модифицированной схемы С. К. Годунова на произвольные нерегулярные сетки//Ученые записки ^AFH.— 1986. Т. XVII, № 2.
9. Roe P. L. Aproximate Riemann solvers, parameter vectors and difference schemes//J. of Comp. Phys.— 1981. Vol. 43.
10. Дьяконов Ю. Н. Пространственное обтекание затупленных тел с учетом равновесных физико-химических реакций//ДДН.— 1964. Т. 157, № 4.
11. Kireev A. Yu., Yumashev V. L. Numerical simulation of reacting and radiating flow in HSCT engines / Lecture Notes in Physics, 16th Int. Conference on Numerical Methods in Fluid Dynamics.— 1998. Vol. 515.
12. Ворожцов И. И., Юмашев В. Л. Система стандартного представления теплофизических свойств моделей реального газа в численных задачах газовой динамики//Труды ^Ara.— 1993. Вып. 2490.
13. Петров К. П. Aэродинамика ракет.— М.: Машиностроение.— 1977.
14. Бабенко К. И., Воскресенский Г. П., Любимов A. Н., Русанов В. В. Пространственное обтекание гладких тел идеальным газом.— М.: Наука.— 1964.
15. Aртонкин В. Г., Леутин П. Г., Петров К. П., Столяров Е. П. Aэродинамические характеристики острых и затупленных конусов при дозвуковых и сверхзвуковых скоростях//Труды ^AFH.— 1972. Вып. 1413.
Рукопись поступила 26/XII2002 г.
