Моделирование газовых потоков около поверхности гиперзвуковых летательных аппаратов методом начального аналитического приближения

Размещено на http://www.allbest.ru/

На правах рукописи

МОДЕЛИРОВАНИЕ ГАЗОВЫХ ПОТОКОВ ОКОЛО ПОВЕРХНОСТИ ГИПЕРЗВУКОВЫХ ЛЕТАТЕЛЬНЫХ АППАРАТОВ МЕТОДОМ НАЧАЛЬНОГО АНАЛИТИЧЕСКОГО ПРИБЛИЖЕНИЯ

05.13.18 - Математическое моделирование, численные методы и комплексы программ

05.07.01 - Аэродинамика и процессы теплообмена летательных аппаратов

АВТОРЕФЕРАТ

диссертации на соискание ученой степени доктора технических наук

Владимир Пантелеевич Котенев

Москва - 2008

Работа выполнена на кафедре вычислительной математики и математической физики МГТУ им. Н.Э. Баумана.

Официальные оппоненты: доктор технических наук, профессор Феоктистов Владимир Васильевич

доктор технических наук

Финченко Валерий Семенович

доктор физико-математических наук, профессор

Формалев Владимир Федорович

Ведущая организация: ФГУП ЦАГИ имени профессора Н.Е. Жуковского

Защита состоится «___» _____________ 2008 г. в __ ч. на заседании диссертационного совета Д 212.141.15 в Московском государственном техническом университете имени Н.Э. Баумана по адресу: 105005, г. Москва, 2-ая Бауманская ул., д.5.

С диссертацией можно ознакомиться в библиотеке МГТУ им. Н.Э. Баумана по адресу: 105005, г. Москва, 2-ая Бауманская ул., д.5.

Автореферат разослан «___» ____________ 2008 г.

Ученый секретарь диссертационного совета

кандидат технических наук, доцент А.В. Аттетков

Общая характеристика работы

Актуальность проблемы. При проектировании летательных аппаратов (ЛА) необходимы предварительные сведения о полях распределения газодинамических величин около их поверхности. В случае гиперзвуковых скоростей полета в ударном слое происходит изменение химического состава газа и характеристики течения должны уточняться путем учета влияния физико-химических процессов. Экспериментальные исследования дороги, а во многих случаях в наземных условиях моделирование физико-химических процессов, сопровождающих полеты ЛА в атмосфере Земли, принципиально невозможно. Поэтому получение необходимой информации может быть осуществлено посредством решения на ЭВМ соответствующим образом поставленных математических задач.

Нестационарные течения описываются системой уравнений, принадлежащей к гиперболическому типу при любых числах Маха, что позволяет применять метод установления по времени для получения искомого стационарного решения. Его появление связано с работой, основанной на конечно-разностной схеме, предложенной С.К. Годуновым в 1959 году. Различные модификации этого метода успешно применяются при расчете смешанных течений, возникающих при исследовании обтекания тел, в том числе для расчета химически неравновесного ударного слоя. С целью сокращения временных затрат в этом случае используется специальная процедура раздельного решения нестационарных газодинамических и стационарных релаксационных уравнений [Савинов К.Г., Шкадова В.П.,1975]. В данной работе для решения задач сверхзвукового обтекания тел с учетом физико-химических превращений в потоке используется метод установления на основе конечно-разностной схемы Мак-Кормака. Предложено использовать дивергентную форму записи уравнений газовой динамики, что позволяет ввести единственный коэффициент , связывающий эти две системы уравнений [10,11]. Эффективный показатель адиабаты , плавно меняясь в поле течения многокомпонентной смеси, является удобным параметром для организации итерационного процесса при раздельном интегрировании газодинамической и релаксационной систем уравнений. В целом, численные методы характеризуются простотой применения, слабыми упрощающими допущениями, возможностью включения в процедуры оптимизации, полным описанием течения, отсутствием ограничений по числам Маха и Рейнольдса, умеренными затратами средств. При этом в качестве недостатков можно отметить не всегда хорошую точность разностных методов, недостаточное быстродействие и объем памяти ЭВМ.

Приближенные методы позволяют получать решение в замкнутой форме при минимальных временных затратах. Например, по формуле Ньютона давление на элемент поверхности тела определяется в основном углом встречи элемента с направлением набегающего потока. Однако приближенные теории не дают возможности найти давление на участки поверхности тела с малыми углами встречи с направлением набегающего потока (<30O), так как давление на них по этим теориям близко к нулю или отрицательно, что не всегда соответствует действительности. При применении метода установления по времени успех решения задачи обтекания тел сверхзвуковым потоком газа во многом зависит от начальных данных, задаваемых с помощью приближенных решений. Поэтому по-прежнему актуально создание такого метода для определения давления и других параметров на поверхности тел, который сочетал бы простоту приближенных формул, точность строгих методов и был пригоден для сферы, эллипсоида, параболоида, гиперболоида и других тел этого класса, характерного для носовой части ЛА.

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

Цель диссертационной работы. Целью работы является создание новых методов моделирования обтекания ЛА как в двумерной, так и в пространственной постановке потоком невязкого или вязкого газа. Основу диссертации составляют:

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

2. Создание экономичных алгоритмов определения контуров ударной волны для использования их в качестве начальных данных в методе установления при решении задач сверхзвуковой аэродинамики.

3. Анализ влияния различных термохимических моделей на параметры течения газа, выбор экономичной модели с точки зрения учета основных факторов.

4. Разработка программных комплексов для решения задач обтекания тел сверхзвуковым потоком химически реагирующего газа.

5. Определение с помощью разработанных методик аэродинамических характеристик и параметров многокомпонентной смеси газов около поверхности ЛА.

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

Научная новизна. В данной работе разработан новый метод начального аналитического приближения (НАП) для расчета параметров установившихся газовых потоков на поверхности обтекаемых тел, сочетающий простоту приближенных подходов и точность строгих численных и аналитических решений. Применение метода НАП для моделирования обтекания гиперзвуковых ЛА в рамках уравнений Навье-Стокса с неравновесными химическими реакциями показало сокращение в несколько раз времени расчета установившихся режимов по сравнению с методами, использующими традиционный подход к заданию исходных данных.

Разработан программный комплекс, реализующий метод НАП для расчета обтекания осесимметричных затупленных тел и для специальных типов пространственных компоновок ЛА.

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

Предложены новые формы записи уравнений гидродинамики в строго консервативной форме.

Практическая ценность. Рассматриваемый в диссертации подход к решению задачи пространственного сверхзвукового обтекания, разработанные аналитические методы, алгоритмы и метод НАП численного решения представляют собой математический аппарат для моделирования течения газа около поверхности современных и перспективных ЛА. Созданный метод НАП позволяет проводить параметрические расчёты для различных режимов пространственного течения с учётом реальных свойств газа. При этом предложенный в работе метод может быть использован для решения как полных уравнений Навье-Стокса, так и упрощённых уравнений Навье-Стокса в приближении тонкого слоя. С помощью разработанной программы численного решения уравнений химически неравновесного вязкого ударного слоя можно проводить систематическое моделирование неравновесных течений около неразрушающихся поверхностей с различными каталитическими свойствами. Полученные в работе результаты используются в практике инженерных расчетов аэродинамических, тепловых нагрузок ЛА и влияния плазменных образований на работу различных приборов, а также в учебном процессе в МГТУ им Н.Э.Баумана.

На защиту выносятся следующие положения:

1. Для генерации начальных данных при применении метода установления по времени построена математическая модель двумерных течений газа в переменных давление - функция тока.

2. Предложены новые аналитические зависимости для определения давления и других параметров течения на поверхности обтекаемых тел, а также на произвольной линии тока с известной формой.

3. Установлена возможность перехода к трехмерным течениям с помощью осесимметричной аналогии.

4. Получены новые формы записи уравнений гидродинамики в строго консервативной форме для применения эффективных и экономичных конечно-разностных схем сквозного счета.

5. Разработан метод начального аналитического приближения для моделирования пространственного обтекания гиперзвуковых ЛА в рамках уравнений Эйлера и Навье-Стокса с учетом физико-химических превращений в потоке.

Апробация работы. По теме диссертации опубликовано 25 работ, из них 7 - в рецензируемых научных журналах перечня ВАК.

Основные результаты диссертации докладывались: на научных конференциях «Ломоносовские чтения» МГУ им. М.В. Ломоносова (1980,1981); на Всесоюзных школах-семинарах по механике реагирующих сред(Междуреченск, 1980;1981); на II Всесоюзной школе-семинаре «Современные проблемы аэрогидродинамики»(Махачкала,1982); на научной конференции «Задачи космического образования в XXI веке» МАИ (2001); на Всероссийской научно-технической конференции «Аэрокосмические технологии» МГТУ им. Н.Э. Баумана (2002); на VI Международной научно-практической конференции «Человек и космос»( Днепропетровск, 2004); на I Международной научно-технической конференции «Аэрокосмические технологии», посвящённой 90-летию со дня рождения академика В.Н.Челомея (Москва-Реутов, 2004); на XXIX академических чтениях по космонавтике(Москва, 2005), а также на ряде семинаров в МГУ, МГТУ и отраслевых институтах.

Структура и объем работы. Диссертация состоит из введения, трех глав, заключения, списка литературы и приложения. Диссертация изложена на 222 страницах машинописного текста, содержит 96 иллюстраций. Библиография включает 142 наименования.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

летательный аппарат аэродинамика газ

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

Первая глава посвящена разработке метода начального аналитического приближения (НАП) для моделирования идеальных газовых потоков в окрестности ЛА. Особенности моделирования гиперзвуковых течений с учетом неравновесных физико-химических превращений в потоке описаны в 1.1. В 1.2 дается общая концепция и основные этапы реализации метода НАП. Стационарное решение задачи получается как предельное решение нестационарной системы уравнений газовой динамики, когда время , при этом эволюция во времени промежуточных решений не важна.

Важным аспектом с точки зрения получения стационарного решения является удачное задание начальных данных. Классический подход состоит в использовании простых зависимостей типа формулы Ньютона для вычисления газодинамических параметров на поверхности тела. Для тел сложной формы отдельным элементам конструкции ставятся в соответствие простые тела с тем же, что и на них распределением параметров газовой динамики. Вид образующей головной ударной волны, как правило, принимается параболическим. Между поверхностями ударной волны и тела осуществляется линейная интерполяция для определения всех искомых функций в узлах разностной сетки. Такой приближенный подход сопряжен с большими временными затратами получения стационарного решения, а иногда приводит к срыву процедуры интегрирования нестационарных уравнений газовой динамики. Поэтому задача определения начальных данных с хорошей точностью является актуальной, а ее решение является основой метода начального аналитического приближения (НАП). В 1.3 рассматривается преимущественно аналитический подход к решению этой задачи. Необходимая программная реализация не требует сколько-нибудь значительных временных затрат. Исходя из записи стационарной системы газовой динамики для плоских и осесимметричных течений газа в динамических переменных, в 1.3.1 получено уравнение относительно искомых функций. Это уравнение эквивалентно условию равенства нулю полной кривизны плоскости, в которой происходит течение газа в двумерном случае. Его решение в 1.3.2 используется для определения давления на тех участках, где угол встречи набегающего потока с поверхностью тела не превышает . При больших углах удовлетворительные результаты получаются по формуле Ньютона. Плотность и скорость находятся из условия постоянства энтропии на теле и интеграла Бернулли. Распределение давления на поверхности тела используется для построения контура ударной волны, а значит, и нахождения параметров на ней (1.3.3). Для определения параметров на известной линии тока, не обязательно совпадающей с поверхностью тела, предлагается также использовать специальный ряд скалярных произведений (1.3.4), применение которого рассматривается в 1.3.6. Для расчета начальных данных на линии тока в окрестности точки торможения проводится аналогия систем уравнений газовой динамики и несжимаемой жидкости (1.3.5). Построение начального приближения в трехмерном случае осуществляется с использованием осесимметричной аналогии (1.3.7). В 1.4 описывается программная реализация принципа установления нестационарной системы уравнений газовой динамики в области эллиптичности стационарной системы уравнений, а в 1.5 - численное интегрирование стационарных уравнений в сверхзвуковой части потока.

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

Рис.1

Здесь - компоненты вектора скорости, - единичные базисные векторы. Осесимметричные стационарные течения газа описываются следующей системой уравнений:

(1)

где - плотность, - энтропия газа.

Введем функцию тока , т.е. такую функцию , что . В переменных уравнение неразрывности удовлетворяется тождественно, условие сохранения энтропии вдоль линии тока: , авторое и третье уравнения системы (1):

(2)

где - контравариантная компонента вектора скорости.

Как показывает анализ численных решений, линии тока и изобары пересекаются вблизи поверхности тела в одной точке, поэтому введение координат корректно в этой области (рис.2).

Рис. 2. ударная волна, - - - изобары, линии тока

Ковариантные и контравариантные базисные векторы, якобиан преобразования и коэффициенты метрической матрицы имеют вид:



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

Отметим, что контравариантная составляющая скорости , т.е. на поверхности тела (линии тока) выполнено условие непротекания. Поэтому

. (3)

Здесь и - ковариантные составляющие скорости, причем . Система уравнений (2) может быть записана в виде проекций векторного уравнения на контравариантные базисные вектора: Индексу соответствует , индексу - , - ковариантная производная по аргументу ковариантных компонент вектора скорости. Умножая эти уравнения на , а затем, складывая их (по повторяющимся немым индексам производится суммирование), получим . Откуда, с учетом , Если - расстояния вдоль линии тока и по нормали к ней, - радиус кривизны линии тока, то:

. Следовательно, и откуда, с учетом выражения (3), получаем: . (4)

Т.к. набегающий поток однородный, то .

В координатах : . (5)

Течение газа происходит в плоскости, поэтому в переменных :

Подставляя сюда метрические коэффициенты и используя (4),(5):

.

Вместо аргументов можно использовать координаты , тогда для функции имеем: . (6)

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

Если , то для оценки давления на телах, обтекаемых сверхзвуковым потоком газа, применяется формула Ньютона:

. (7)

Индексами обозначены значения параметров перед скачком, в точке торможения потока (критической точке) и на поверхности тела соответственно. При этой формулой пользоваться нельзя, т.к. здесь распределение давления существенно зависит от числа Маха. Иногда используется приближенный подход, при котором распределение давления по местному углу наклона поверхности тупого тела берется таким же, как и для сферы. При этом относительная погрешность может достигать 60% и более, т.к. разность между величинами часто имеет такой же порядок, как и сами значения давления. Формула (7) дает еще более высокую погрешность.

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

. (8)

По формуле (8) определяем на заданном теле . Используя таблицы газодинамических функций, находим . Для проверки (8) исследованы тела (контуры представлены на рис. 3) и режимы их обтекания, рассмотренные в таблицах Любимова А.Н., Русанова В.В.,1970:

1.Эллипс при =2,4,6,20;

2.Эллипс при =4,6,20;

3.Парабола (y2=2px, p=1) при =4,6,20;

4.Парабола (p=0.25) при =2;

5.Парабола (p=0.5) при =2,4,6,20. Рис.3

Во всех вариантах относительная погрешность, как правило, не превышает 5%, иногда локально достигая 10%. Распределение давления по поверхности тела в сочетании с законами сохранения используется для «быстрого» приближенного определения формы ударной волны в 1.3.3. Некоторые результаты сравнения даны на рис.4-6 (крестики - табличные данные).



Рис.4 Рис.5



Рис.6. Отход ударной волны

a - сфера; b - эллипсоид b/a=1/2;

с - эллипсоид b/a=3/2; d - параболоид p=1;

e - параболоид p=1/4; f - параболоид p=1/2.

В 1.3.4 для уравнений газовой динамики, записанных в переменных «давление - функция тока», получено разложение решения в ряд по приращению давления и коэффициенты ряда. Этот ряд использован для построения эффективного алгоритма интегрирования решения вдоль линии тока:

(9)

Используя формулу , из уравнения последовательно получаем:

, (10)

, (11)

(12)

т.е.. (13)

В общем случае, Поэтому,

.

Если, например, ,то есть выражения в квадратных скобках из (12). Откуда для :

(14)

Для определения как функции выразим правые части (11),(13),(14) через параметры газовой динамики. Выше получено

. (15)

Сделаем подстановки в первое уравнение движения из (15), тогда, после преобразований, получаем:

. (16)

Поэтому, . (17)

Из уравнений движения Эйлера (15) следует:

. (18)

Используя (13),(16),(18) получаем

. (19)

Пусть давление отнесено к давлению торможения, скорость - к максимально возможной скорости, тогда

. (20)

Аналогично, из (14),(15), (16),(18) следует:

(21)

Полагая в формулах (10),(17),(20),(21) , автоматически получаем коэффициенты искомого разложения в ряд для модели несжимаемой жидкости. В 1.3.5 показано, что разложение (9) в переменных будет таким, как для несжимаемой жидкости в переменных

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

Применяя разложение (9) и учитывая (10), (17), (19), получаем:

(22)

(23)

Для конечно-разностной аппроксимации производных , применяются центральные разности. Формулы (22), (23) эффективны для пошагового с фиксированным продвижения решения вдоль линии тока (поверхности тела) ненулевой кривизны. Для примера они были использованы при определении давления на поверхности сферы, цилиндра, эллипсоида вращения с полуосями a=2b=2c. В качестве начальных данных задавались , , . Здесь давление соответствует сферическому (полярному) углу =200 . При этом проведена параболическая аппроксимация таблиц Любимова А.Н., Русанова В.В. в окрестности критической точки. Для дозвуковых режимов обтекания использована аппроксимация Покровского А.Н., Фролова Л.Г.,1985. В общем случае для сверхзвуковых режимов обтекания начальные данные в окрестности точки торможения потока достаточно точно задаются с помощью формул Ньютона и Буземана или их совместного применения. Для дозвуковых течений можно воспользоваться тем, что уравнения газовой динамики в переменных по форме совпадают с «несжимаемой» системой уравнений в переменных , т.е. использовать в окрестности критической точки известные «несжимаемые» распределения скоростей бесциркуляционного обтекания сферы или цилиндра.

Значения получены с помощью формул (22), (23). Реализован простой алгоритм. При местных числах применялась формула (23) . Вместо , использовались , , т.е. ошибка на одном шаге возрастала до . При местных числах на шаге предиктор применялась формула (22), при этом решалось квадратное уравнение для неизвестной , а на шаге корректор - формула (23).

Рис.7

На рис. 7 приведены данные Любимова А.Н., Русанова В.В., причем сплошная линия соответствует сфере, пунктирная - цилиндру, штрихпунктирная - эллипсоиду. Кружочками обозначены значения давления, полученные по (22),(23). Потенциальность течения, когда , заранее не предполагалась. Сравнение свидетельствует о хорошей точности представленных соотношений.

Из построения решения по (22),(23) следует, что функциядолжна быть одинаковой для семейства линий тока, имеющих общую начальную часть и сопряжение достаточной степени гладкости. Действительно, на рис. 8 изображен набор осесимметричных тел, полученных путем сопряжения без разрыва кривизны сферы с параболоидом вращения. Здесь даны условные знаки для обозначения давления на графиках. Зависимости заметно отличаются друг от друга (-сферический угол). Зависимость же хорошо описывается одной кривой (рис. 9).

Размещено на http://www.allbest.ru/

Рис.8 Рис.9



Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Рис.10 Рис.11

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

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

Для определения начального распределения давления на поверхности ЛА, обтекаемого трехмерным газовым потоком, в 1.3.7 предлагается использовать метод осесимметричной аналогии. Согласно данной аналогии расстояние вдоль линии тока рассматривается как расстояние вдоль оси эквивалентного осесимметричного тела, а коэффициент Ламе, перпендикулярного к линии тока координатного направления как цилиндрический радиус. Каждая невязкая линия тока на поверхности ЛА соответствует телу вращения, обтекаемого под нулевым углом атаки.

На примере эллиптического параболоида, обтекаемого сверхзвуковым потоком газа под ненулевым углом атаки показано, что формула Ньютона дает удовлетворительные результаты в сравнении с эталонными данными только на наветренной стороне. Построено эквивалентное тело вращения, соответствующее подветренной стороне рассматриваемого параболоида в плоскости его симметрии. Установлено, что распределение давления на теле вдоль линии тока, обтекаемого трехмерным газовым потоком, совпадает с распределением давления на поверхности осесимметричного аналога, обтекаемого под нулевым углом атаки, с приемлемой точностью.

В п. 1.4 описан разностный метод решения исходной нестационарной трехмерной задачи с учетом неравновесных физико-химических превращений в потоке около лобовой поверхности ЛА. Начальное приближение задается с использованием подходов, изложенных выше. Уравнения газовой динамики включают в себя уравнения движения Эйлера, неразрывности и энергии. Для их решения применяется явная конечно-разностная схема Мак-Кормака, основанная на дивергентной форме записи системы уравнений записанная в нормированной сферической системе координат. Реализация граничных условий на ударной волне и поверхности тела для совершенного газа аналогична той, что использована в работе [11].

Учет неравновесных физико-химических превращений в потоке сводится к следующей процедуре. Решение газодинамической части системы уравнений при заданной в любой точке расчетной сетки функции ищем на основе метода установления, как и для совершенного газа. Начальное распределение эффективного показателя адиабаты в ударном слое берем постоянным, обычно =1.3333. Используя найденные из решения уравнений газовой динамики давление и скорость, интегрируем вдоль линий тока стационарные релаксационные уравнения. Для этого из расчетных точек ближайшей к ударной волне координатной поверхности (-- координата, нормализующая расстояние между телом и волной) выпускаются линии тока до пересечения с головным скачком. Определив координаты точки пересечения, находим в ней значения газодинамических функций , и Т и концентраций ( для атомов на скачке) при помощи интерполяции. Затем по неявной схеме рассчитываем концентрации компонентов вдоль выстроенной линии тока (плотность, входящая в выражение для , исключается с помощью интеграла энергии, а концентрации -- из условия сохранения элементарного состава многокомпонентной смеси и ее квазинейтральности):

.

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

В п.1.5 рассматриваются вопросы построения решения в сверхзвуковой части потока. В сверхзвуковой расчетной области система внутренних ударных волн, волн разряжения, контактных поверхностей, образующихся на выступах и на изломах обтекаемой поверхности, может быть достаточно сложна. Здесь применяются методы сквозного счета, автоматически удовлетворяющие условиям на разрывах и позволяющие переходить через скачок, размазывая его в пределах нескольких ячеек расчетной области [Годунов С.К., 1957; Kutler P., Warming R., Lomax H., 1973; Рицци В., Клавинс Э., Мак-Кормак У., 1977]. Система уравнений газовой динамики при этом записывается в дивергентной форме для того, чтобы соответствующая ей система разностных аналогов удовлетворяла законам сохранения с точностью до машинного округления. При использовании криволинейных систем координат в уравнениях движения появляются, вообще говоря, слагаемые, не стоящие под знаком производной. Наличие этих членов аналогично действию фиктивных массовых сил, возникающих от использования неинерциальной системы координат. При численном решении задачи присутствие "источников" препятствует достижению полного сохранения импульса, делает неопределенным вид разностной аппроксимации уравнений.

Пусть координаты относительно произвольной системы, а - относительно декартовой системы. Базисные ковариантные и контравариантные векторы:

, , ().

Компонентам векторов, тензоров, скалярам и всем параметрам, характеризующим течение газа, будем приписывать индекс "0" в произвольной фиксированной точке евклидового пространства . Для любой точки , принадлежащей окрестности , введем два набора двухточечных скалярных произведений базисных векторов по следующему закону:

, (24)

. (25)

Величины , при фиксированной точке являются функциями координат . Из определений (24), (25) вытекают следующие свойства:

, , (26)

, (27)

, (28)

, . (29)

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

, ,

где ; ; - контравариантные компоненты скорости; - ковариантные компоненты скорости ( и - составляющие вектора скорости в базисах и ). Используя свойства (26)-(29) и формулу получаем:

, (30)

. (31)

Итак, рассматривая в произвольной точке проекции векторного уравнения импульсов на оси, определяемые векторами или векторами, получаем строго консервативные формы (30-31).

Пусть - декартов базис. Проектируя векторное уравнение импульсов на направления, определяемые векторами , вместо (30) и (31) получаем

(32)

( - компоненты в декартовом базисе, т.е. , а ).

Приведенные рассуждения носят общий характер и в равной степени относятся также к уравнениям Навье-Стокса для многокомпонентного газа, поскольку последние можно записать в инвариантной форме. Строго консервативные формы записи уравнений гидродинамики (30),(31) являются оригинальными и впервые получены автором в [3,4]. Отметим, что существуют и другие способы записи уравнений гидродинамики в строго консервативной форме [Lax P.D., Wendroff, 1960; Vinokur M., 1973; Hicks, Wooten, 1980]. Выбор в пользу той или иной формы записи уравнений зависит от способа введения криволинейной системы координат и конкретной физической задачи.

Удобным с точки зрения уменьшения объема вычислений в эллиптико-гиперболической зоне является такое задание поверхности Коши, когда нормальная составляющая скорости к ней сверхзвуковая, в то время как проекции скорости на направление оси тела может быть еще дозвуковой. Используемая в работе искомая поверхность в процессе счета может деформироваться, и при необходимости трансформируется в плоскость для продолжения счета в цилиндрических координатах. Наличие сложной формы рассчитываемых тел и больших углов атаки делает разумным требование близости линий тока к одному из координатных направлений. Предельными случаями построенной системы являются сферическая и цилиндрическая системы, широко применяемые для расчета соответствующих конфигураций.

Дифференциальные уравнения в частных производных, представляющие собой уравнения неразрывности и импульсов, записываются в строго консервативной форме. При расчете течений совершенного газа они замыкаются интегралом Бернулли: . При рассмотрении течений с учетом неравновесных физико-химических превращений в потоке к этим уравнениям добавляются релаксационные уравнения, записанные вдоль линий тока, уравнения состояния и зависимости констант скоростей химических реакций от плотности, температуры и концентраций компонентов. Подробные сведения об уравнениях химической кинетики и входящих в них переменных приведены в главах 1-3 и таблицах 2,3 приложения. Интегрирование газодинамической части системы осуществляется при помощи явной двухшаговой конечно-разностной схемы второго порядка точности, предложенной Мак-Кормаком. Для интегрирования обыкновенных дифференциальных уравнений химической кинетики, относящихся к классу "жестких", применяется неявная конечно-разностная схема, описанная в 1.4.

На первом этапе решение продвигается вдоль линий тока, оканчивающихся в расчетных узлах нового координатного слоя, где уже известны "предикторные" значения «консервативных» газодинамических комплексов. Координаты точек пересечения линий тока со старым слоем и значения комплексов в них находятся при помощи интерполяции.

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

.

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

Процедура удовлетворения граничным условиям и расчет контура ударной волны осуществляется при помощи соотношений Рэнкина-Гюгонио. Граничные условия на теле: . На поверхности тела используется схема Аббетта.

Процессы конвективного переноса могут проходить медленнее, чем химические процессы, особенно при приближении к равновесию, и для правильного учета последних осуществляется несколько дробных шагов на каждом крупном шаге, определяемом из условия Куранта-Фридрихса-Леви.

Во второй главе приводятся результаты численного моделирования газовых потоков методом НАП. Модель совершенного газа с постоянным показателем адиабаты была использована для проведения расчетов по определению внешнего давления (коэффициент ) по длине и в поперечных сечениях затупленного тела с изломами образующей. Расчет физико-химических реакций в равновесной постановке требует вдвое и более раз большего времени и памяти ЭВМ по сравнению со случаем расчета в рамках совершенного газа. Эта задача значительно упрощается, если определенным образом аппроксимировать таблицу термодинамических функций воздуха с достаточной для практики точностью [Синченко С.Г., 1968]. В расчетах с учетом неравновесных химических процессов временные затраты ЭВМ численного решения задачи возрастают на порядок. При больших скоростях полета из-за высоких температур в ударном слое происходит распад молекул и ионизация. Это приводит к снижению температуры в потоке, увеличению плотности, оказывает влияние, правда в меньшей степени, и на давление газа. Применение различных моделей возбуждения колебательных степеней свободы (Лайтхилла и равновесной) не сказывается на распределении давления в ударном слое, однако влияет незначительно на температуру и концентрации компонентов диссоциированного воздуха. При использовании равновесной модели колебаний более низкие значения температуры и степени диссоциации молекул и объясняются большими энергетическими затратами потока на возбуждение колебаний в этом случае. На рис. 12-13 приведены результаты расчетов в неравновесной постановке распределений поперек ударного слоя температуры и концентраций компоненты в плоскости симметрии течения при обтекании эллипсоида. Режим обтекания: угол атаки , , км. Кривые 1,3 соответствуют наветренной, а кривые 2,4-подветренной частям тела. Квадратики и треугольники - результаты [Иванова В.Н., Радвогин Ю.Б., 1981].

Рис.12 Рис.13

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

Исследование влияния учета различных моделей физико-химических процессов на обтекание классических форм важно для прогнозирования аэродинамики более сложных аппаратов в реальных условиях обтекания. Рассмотрено обтекание конуса с углом полураствора (радиус сферического затупления Rn =5см): M? = 20, H = 50 км, . Показано, что учет неравновесных химических процессов уменьшает на 2-4% значения осевой силы по сравнению с равновесным режимом протекания химических реакций. Значения для совершенного газа также незначительно ниже, чем для равновесного воздуха. В случае равновесия химических реакций коэффициент нормальной силы для малых длин конуса меньше, чем для совершенного газа. При x>5и x>4значения этого коэффициента становятся большими (до 15% при x=10-15) для равновесно реагирующего воздуха. Это связано с разным характером в поведении давления на наветренной стороне конуса. Для равновесного газа минимум давления на поверхности тела в «ложке» вследствие меньших размеров ударного слоя достигается ближе к критической точке, а нарастание давления в дальнейшем происходит более интенсивно и до больших значений, чем для совершенного газа и неравновесного воздуха. Значения в последнем случае лежат выше, чем для равновесного воздуха и совершенного газа с г=1.4. Это объясняется более сильным расширением потока с подветренной стороны конуса. Отличие в может достигать 15% для химически равновесного и 20% для неравновесного режимов при всех углах атаки. Аналогичное поведение наблюдается в распределении . Для определенных длин конусов учет равновесных химических реакций приводит к увеличению до 15% значений Сd по сравнению с совершенным газом. Учет кинетики химических процессов увеличивает значения Сd в зоне минимума (x=4) до 3-4% по сравнению с равновесным режимом. При это различие сохраняется для более длинных конусов, а значения Сd лежат между соответствующими значениями для совершенного и равновесного воздуха.

Численное моделирование методом НАП течений газа около поверхности ЛА в рамках уравнения Навье-Стокса рассматривается в третьей главе. Классический подход к решению задач обтекания состоит в разделении течения на внешнюю невязкую область и вязкий пограничный слой с присущей ему проблемой, связанной с нахождением условий на внешней границе пограничного слоя. B рамках уравнений Навье-Стокса расчёт проводится для всей области течения сразу, и указанные трудности не возникают. Хотя для расчета обтекания затупленных тел вязким газом разработано множество программ, часто они не могут быть использованы в качестве быстрого и надежного инструмента определения характеристик обтекания хотя бы наиболее часто встречающихся форм. В частности, в универсальных методах сквозного счёта используются достаточно трудоёмкие процедуры построения качественных сеток в области интегрирования для каждого расчётного случая, а в некоторых случаях и для режима обтекания. Это приводит порой к очень большим затратам получения решения. Расчеты при помощи соответствующих алгоритмов и программ наряду с экспериментальными данными являются эталонными для оценки достоверности результатов, получаемых другими методами. С другой стороны, упрощенные, в той или иной степени, уравнения Навье-Стокса имеют существенные ограничения на их применение по характеристикам режимов обтекания. Поэтому одной из задач данной работы было создание эффективного численного метода решения полных уравнений Навье-Стокса, позволяющего моделировать пространственные гиперзвуковые течения с учётом физико-химических превращений в газовой фазе за приемлемое для практических расчётов время.

Для получения стационарного решения при обтекании тела равномерным сверхзвуковым потоком в данной работе используется метод установления нестационарной системы, состоящей из уравнений неразрывности, Навье-Стокса, и диффузии продуктов химических реакций. Набегающий сверхзвуковой поток рассматривается как смесь 23.3% кислорода O2 и 76.7% азота N2. Обтекающий тело газ состоит из семи компонентов (O, N, NO, O2, N2, NO+, e-), а модель химических реакций, протекающих в ударном слое при высоких температурах, является общепринятой [Rakich J.V., Bailey H.E, Park C., 1983]. Коэффициент вязкости неравновесной смеси газов рассчитывается по формуле Буденберга-Уилки. Тепловой поток к поверхности тела обусловлен теплопроводностью и переносом энергии диффундирующими компонентами. В работе [Blottner F.G., 1964] показано, что концентрации компонентов слабо зависят от степени упрощения диффузионной модели для рассматриваемых режимов обтекания, поэтому используется модель бинарной диффузии. Распределение энергий по поступательным и вращательным степеням свободы принимается равномерным. Колебания молекул считаются возбужденными равновесно, или же, по модели Лайтхилла, равны половине максимально возможного равновесного значения. Замыкают систему уравнений сведения о константах равновесия и скоростях химических реакций и уравнения состояния.

Естественные граничные условия для системы заданы на бесконечном удалении от тела (в невозмущённом потоке) и на самом теле. На поверхности обтекаемого тела для компонентов скорости задавались условия прилипания . Режим теплообмена - охлаждение стенки или изменение температуры по закону, близкому к линейному. Для концентраций компонент атомов O и N, молекул NO граничные условия имеют вид , как для некаталитической поверхности стенки. В отношении рекомбинации заряженных частиц используются условия идеально каталитической поверхности . При вычислении концентраций используются условия сохранения атомарного состава и квазинейтральности смеси. Для обеспечения устойчивости расчёта накладывалось ещё одно условие на стенке:.

Переход через ударную волну происходит «замороженным» образом, т.е. на ней концентрации совпадают с соответствующими значениями набегающего потока. Для задания «внешней» границы области интегрирования необходимы априорные данные о структуре поля течения. Для больших чисел Рейнольдса (Re?>103) можно пренебречь влиянием структуры тонкой головной ударной волны на течение вниз по потоку и принять, что ударная волна является поверхностью разрыва газодинамических параметров, на которой выполняются нестационарные условия Рэнкина-Гюгонио.

Основой предлагаемого метода расчета является расщепление задачи обтекания по физическим процессам.

На первом этапе с помощью процедуры, описанной в 1.4, интегрируется «невязкая» часть системы. Она включает в себя уравнения движения Эйлера, неразрывности и энергии в невязкой постановке.

На втором этапе решение уточняется с использованием «вязкой» части уравнений Навье-Стокса и энергии. «Вязкая» подсистема интегрируется при помощи метода прогонки.

Наличие узких областей с большими градиентами параметров потока, таких как пограничный слой вблизи тела, зона размазанной ударной волны при низких числах Рейнольдса, не позволяет получать надёжные количественные результаты, используя равномерные сетки. Для адекватного учёта поведения параметров в этих областях больших градиентов вводится сгущение в радиальном направлении.

Для расчёта «вязкой» подсистемы уравнений используется неявная разностная схема. Для её описания представим «вязкую» часть уравнений количества движения и энергии в виде



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

, ,

где - любая из неизвестных функций; и заданы или определяются из граничных условий.

Таким образом, в «вязкой» части последовательно решаются разностные уравнения движения, записанные неявно относительно искомых параметров. Затем, с использованием полученного поля скоростей решаются уравнения энергии и диффузии. Для интегрирования уравнений диффузии метод прогонки применяется в сочетании со специальной неявной разностной схемой, используемой для аппроксимации источниковых членов и описанной в 1.4. При этом сначала решается система разностных уравнений с источниковыми членами, но без учета диффузии. Затем учитывается диффундирование компонентов смеси путем последовательного решения уравнений этой подсистемы методом прогонки.

Для исследования течений около удлиненных тел используется разделение всей области интегрирования на ряд взаимно перекрывающихся подобластей и последовательный расчёт в каждой из них. Такое разделение области возможно в силу слабой передачи возмущений вверх по потоку при обтекании тел сверхзвуковым набегающим потоком вязкого газа [Кокошинская Н.С., Павлов Б.М., Пасконов В.М., 1980], и позволяет проводить расчёт длинных затупленных тел до 100 калибров и более.

Из-за большой сложности уравнений Навье-Стокса невозможно получить аналитическое выражение устойчивости для описанной схемы расчета. Однако практическое применение алгоритмов показало возможность использования, без потери устойчивости счёта, эмпирической формулы , где - коэффициент запаса (=0.8); - определяется по критерию Куранта-Фридрихса-Леви для невязких линейных гиперболических уравнений в частных производных.

Предусматривается возможность применения упрощенных уравнений Навье-Стокса для некоторых режимов обтекания или для генерации начальных данных. Расщепление задачи по физическим процессам, использование смешанных явно-неявных численных алгоритмов позволяет значительно сократить трудоемкость получения практических результатов. Затраты машинного времени при расчете пространственного течения около затупленных конусов примерно в 10 раз меньше, чем при использовании полностью неявных схем во всей рассматриваемой области без выделения головной ударной волны. Программа расчета пространственного неравновесного гиперзвукового течения вязкого газа позволяет исследовать все режимы протекания химических реакций: “замороженный”, неравновесный и равновесный. Приводятся результаты расчётов сверхзвукового обтекания затупленных тел в рамках полных уравнений Навье-Стокса, а также сравнения с экспериментальными и расчётными данными других авторов по давлению, температуре, концентрациям химических компонент и другим параметрам течения.

Сегментально-коническая форма тел состоит из лобовой части в виде шарового сегмента с сужающейся конической кормовой части и используется при проектировании спускаемых космических аппаратов. Характер обтекания таких тел является очень сложным, так как связан не только с наличием дозвуковых и сверхзвуковых зон течения, но и с областями отрыва потока. При осесимметричном сверхзвуковом обтекании (угол атаки ) расчет распределения давления по лобовой поверхности не вызывает затруднений. Вдоль боковой поверхности действует практически постоянное давление, свидетельствующее об открытом, разомкнутом отрыве с угловой кромки в месте соединения лобовой и боковой поверхностей. Величина давления может быть на боковой поверхности меньше, чем в набегающем сверхзвуковом потоке, а ее определение является не простой задачей. Представлено полученное по предложенному в диссертационной работе методу распределение давления, и характер поведения линий тока, начинающихся в зоне затупления аппарата. Угол конической поверхности . Режим обтекания: , число вычислено по диаметру миделя. Полученные в расчете коэффициенты давления на боковой конической поверхности практически совпадают с данными эксперимента [К.П. Петров, 1998].

Рассмотрим некоторые выявленные особенности неравновесных пространственных течений на примере расчетов при помощи предложенного в данной работе численного метода обтекания вязким газом затупленных по сфере конусов. Характерной особенностью неравновесных режимов течения вязкого газа является явно выраженная «ложка», которая отсутствует, если физико-химические процессы в потоке не рассматриваются [17]. Различия тепловых потоков здесь достигают 20% в сторону уменьшения для модели, в которой учитываются химические реакции. Отмечено, что при обтекании тел под большими углами атаки ионизация и диссоциация газа происходит преимущественно с подветренной стороны тела. Для подтверждения этого факта рассмотрено обтекание затупленного конуса и приведено сравнение результатов с данными [Bhutta B.A., Lewis C.H., 1989]. В этой работе принималось, что концентрации компонент на стенке равны концентрациям в невозмущённом потоке, а температура тела . Угол полураствора конуса 7°, радиус затупления равен 5.08 см, общая длина тела - 50.8 см, высота полета H=53.3 км, число , угол атаки . Максимальная концентрация с подветренной стороны тела примерно в два раза превосходит её наибольшее значение с наветренной стороны (рис. 15). Линии равных концентраций на рис.14 соответствуют данным работы [Bhutta B.A., Lewis C.H., 1989] и подтверждают этот вывод. Для атомарного азота максимальная концентрация с подветренной стороны превышает её максимальное значение с наветренной стороны тела примерно в 60 раз. Распределение концентраций и также подтверждает интенсивную диссоциацию молекулярного кислорода и молекулярного азота с подветренной стороны тела. И хотя максимальная концентрация в три раза меньше, чем с наветренной стороны, ионизованный слой с подветренной стороны имеет большую толщину примерно в 13 раз.

...