Метод псевдохарактеристик для расчета сверхзвуковых пространственных течений газа Нурбаев Улукбек Джаныбекович

Содержание к диссертации

Введение

Глава 1. Метод псевдохарактеристик 17

Глава 2. Задача сверхзвукового обтекания пространственных тел 26

Глава 3. Метод псевдохарактеристик для задачи обтекания пространственных тел 39

3.1 Расчет искомых функций во внутренней точке области течения .40

3.2 Расчет искомых функций в точке на обтекаемом теле 48

3.3 Расчет искомых функций в точке на ударной волне 50

Глава 4. Устойчивость метода псевдохарактеристик 55

Глава 5. Примеры расчетов 70

5.1 Осесимметричное обтекание 70

5.2 Пространственное обтекание 77

Заключение 85

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

• Задача сверхзвукового обтекания пространственных тел
• Расчет искомых функций в точке на обтекаемом теле
• Расчет искомых функций в точке на ударной волне
• Пространственное обтекание

Введение к работе

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

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

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

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

Неявные сеточные конечно-разностные методы (см., например: Бабенко К.И., Воскресенский Г.П., Любимов А.Н., Русанов В.В., 1964 г.) имеют простую реализацию, но требуют для получения нужной точности много машинного времени. В методах такого типа при решении задачи о сверхзвуковом трехмерном течении газа сложно удовлетворять краевым условиям на поверхности тела и на ударной волне. В 1960-е годы Годунов С.К. предложил явный сеточный метод решения нестационарных задач газовой динамики, записывая уравнения движения в интегральной форме и используя автомодельные задачи о распаде разрыва. В 1972 г. Иванов М.Я. и Крайко А.Н. распространили метод на решение пространственных стационарных задач газовой динамики. (Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П., 1976 г.).

В обратных характеристических методах (Кацкова О.Н., Чушкин П.И., 1968 г.) необходима двумерная интерполяция искомых величин на предыдущем расчетном слое. Прямые характеристические методы (Борисов В.М., Михайлов

И.Е., 1978 г.) обладают высокой точностью и быстродействием, но имеют плавающую сетку, что требует дополнительных вычислений для определения газодинамических функций в фиксированных точках.

Сеточно-характеристические методы (Магомедов К.М., Холодов А.С., 1988 г.) опираются на достоинства сеточных и характеристических методов, однако для его использования необходимы одномерные интерполяции газодинамических величин на предыдущем расчетном слое.

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

Цели и задачи работы. Основная цель диссертационной работы состоит в разработке численного метода для решения задач обтекания пространственных тел сверхзвуковым потоком газа. Для достижения поставленной цели были сформулированы следующие задачи:

1. Разработка численного метода, обладающего фиксированной сеткой, быстродействием и точностью

2. Исследование устойчивости разработанного метода

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

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

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

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

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

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

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

4. Научно-практический семинар "Теория, численные методы и математический эксперимент в газовой динамике", ЦИАМ, Москва, 2009 г.

5. 53 научная конференция МФТИ, Москва-Долгопрудный, 2010 г.

6. Седьмой международный аэрокосмический конгресс, Москва, 2012 г.

7. Научно-практический семинар в Институте автоматизации проектирования, Москва, 2013 г.

В диссертационной работе использованы результаты, полученные автором в ходе исследований, проводимых в рамках научно-исследовательской работы по проекту РФФИ № 11-01-00835-а.

Публикации и личный вклад. По теме диссертации опубликованы 5 работ [1-5], в том числе 2 работы в издании, входящем в перечень ведущих журналов и изданий, рекомендованных ВАК для публикации основных результатов диссертаций на соискание ученой степени кандидата наук. Получено свидетельство о государственной регистрации программы для ЭВМ [6]. Личный вклад автора заключается в разработке и реализации численного метода, разработке комплекса программ, адаптированного для расчета течения газа при пространственном обтекании тел, анализе результатов расчетов, построении математической модели для исследования устойчивости метода, а также подготовке публикаций.

Структура и объем диссертации. Диссертационная работа состоит из введения, обзора литературы, пяти глав, заключения и списка использованной литературы. Диссертация изложена на 93 страницах, включает 10 таблиц, 37 рисунков и 77 наименований использованной литературы.

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

Итак, с помощью метода характеристик можно определить искомую функцию в точке 4 области течения, которая не совпадает с фиксированным узлом сетки 3. Метод характеристик имеет высокую точность и быстродействие, поскольку основывается на внутренних свойствах уравнений. Однако при использовании метода характеристик расчетная сетка получается плавающей, что не всегда является удобным и практичным для решения конкретных прикладных задач.

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

Сеточно-характеристические методы имеют фиксированную сетку. Поэтому, шаблон для расчета выглядит следующим образом (см. Рис. 1.3): г

Расчет искомых функций сеточно-характеристическим методом Из точки 3, совпадающей с узлом фиксированной сетки, выпускается характеристика под углом наклона — = а и пересекает ось г в некоторой точке /. Как правило, точка / не совпадает с фиксированными узлами сетки 1 и 2. Поэтому с помощью интерполяции определяют значение искомой функции в точке / по известным значениям в узлах сетки 1 и 2. Далее, используя условие совместности, расситывают значение искомой функции в точке 3. Разностные уравнения для определения искомой функции в фиксированной точке 3 расчетной сетки выглядят следующим образом (при линейной интерполяции):

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

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

Рассмотрим идею метода псевдохарактеристик. Как и прежде, будем решать уравнение переноса

Уравнения (1.2) и (1.3) представляют собой альтернативную запись уравнения переноса и поэтому позволяют определить искомые функции в области течения газа. Далее введем сетку в рассматриваемой области течения газа. Шаблон для расчета искомых функций состоит из трех точек 1, 2 и 3, совпадающих с узлами сетки (см. Рис. 1.4).

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

Отметим, что левая часть уравнений не содержит выводящих производных с линий 1-3 и 2-3, а правая часть - содержит (в данном случае иг). Поэтому введем дополнительную искомую переменную R = ur. Тогда получим следующую систему, где искомыми функциями являются вектор (u{x,r),R)\

Идея рассматривать выводящие производные в качестве новых дополнительных искомых функций для случая метода характеристик была предложена Борисовым В.М. [62]. Идея использовать производную искомой функции в качестве новой искомой функции для случая конечно-разностных схем рассматривалась Грудницким В.Г. и Прохорчуком Ю.А. в 1977 г. [63]. В данной работе авторы предложили способ построения явных разностных схем на трехточечном (минимальном) шаблоне для решения гиперболических уравнений дивергентного вида. Кроме того, авторы выписывали дополнительное дифференциальное уравнение для новой искомой функции. В работе [64] приведен расчет двумерного нестационарного сверхзвукового обтекания источника энергии. Система уравнений, полученная авторами с помощью дифференцирования исходных уравнений, сохраняет характеристические свойства гиперболических уравнений и в литературе называется продолженной или расширенной системой уравнений [65].

Расчет искомых функций в точке на обтекаемом теле

Соединив между собой точки 1-3, получим псевдохарактеристику 1-3 (первого семейства). Данную линию можно записать в параметрической форме х = х, г = г(х), Р = РІ, взяв в качестве параметра независимую переменную х. Тогда на этой линии любая произвольная функция g(x,r,(pt) = g{x,r{x),(pi) является функцией одной переменной х. В таком случае, производная по х от этой функции вдоль псевдохарактеристики 1-3 примет вид где gx и gr - частные производные по соответствующим переменным, а индекс 1 указывает, что равенство записано на псевдохарактеристике 1 -3. Рассмотрим теперь некоторый венок, лежащий в плоскости xs. Уравнение венка имеет вид г = r(xs, ф), таким образом, его можно записать в параметрической форме х = xs, г = г(ф), ср = ср. Тогда на венке любая произвольная функция g(xs,r,y) = g(xs,r((p),(p) является функцией одной переменной (р. Полная производная по р вдоль венка примет вид

Введем новую искомую вектор-функцию R = Ur. В таком случае мы получим в точке 3 следующую систему уравнений для функций (U,R) \

Покажем, что функция U, являющаяся решением системы (3.6) в точке 3, также является решением исходной системы (2.8). Вычтем из первого уравнения второе уравнение системы (3.6), получим:

Таким образом, из (3.6) следует, что Ur = R на псевдохарактеристиках 1-3 и 2-3 в некоторой меридиональной плоскости. Теперь сложим два уравнения системы (3.6), получим:

В результате, получим исходную систему (2.8). Поэтому, функция U, являющаяся решением системы (3.6) в точке 3 расчетного шаблона, также является решением системы (2.8). Системы (3.6) и (2.8) являются эквивалентными. Систему (3.6) будем аппроксимировать на двух псевдохарактеристических поверхностях 1-3 и 2-3. Псевдохарактеристические поверхности не являются, вообще говоря, характеристическими, однако, как и условия совместности на характеристических поверхностях, система (3.6) не содержат производных в направлениях, пересекающих эти поверхности. Поэтому аппроксимация уравнений на псевдохарактеристических поверхностях осуществлялась аналогично аппроксимации канонической характеристической системы в [44]. А именно, производные по х будем аппроксимировать центральными разностями, а производные по р находить дифференцированием кубического сплайна, построенного по всем точкам соответствующего венка.

Выбранный способ обеспечивает второй порядок аппроксимации по х и ср. Действительно, запишем значения функции U в точках шаблона, используя разложение функции в ряд Тейлора в окрестности средней точки, получим из=Щхср+Ьх/2) = Щхср) + dU dx Так, производные по переменной х имеют второй порядок аппроксимации на средней точке отрезка [х17х3].

Производные по р находятся с помощью дифференцирования кубического сплайна. Тогда для первой производной некоторой достаточно гладкой функции fix) на отрезке [хк,хш] справедливы следующие оценки погрешности [66]

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

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

N. Искомыми величинами являются вектор-функция (U,R) в точке 3 шаблона. Поэтому мы получаем 8N уравнений с 8N неизвестными и таким образом определяем все искомые величины сразу по всему венку.

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

Расчет искомых функций в точке на ударной волне

Систему (3.6) будем аппроксимировать на двух псевдохарактеристических поверхностях 1-3 и 2-3. Псевдохарактеристические поверхности не являются, вообще говоря, характеристическими, однако, как и условия совместности на характеристических поверхностях, система (3.6) не содержат производных в направлениях, пересекающих эти поверхности. Поэтому аппроксимация уравнений на псевдохарактеристических поверхностях осуществлялась аналогично аппроксимации канонической характеристической системы в [44]. А именно, производные по х будем аппроксимировать центральными разностями, а производные по р находить дифференцированием кубического сплайна, построенного по всем точкам соответствующего венка.

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

Так, производные по переменной х имеют второй порядок аппроксимации на средней точке отрезка [х17х3]. Производные по р находятся с помощью дифференцирования кубического сплайна. Тогда для первой производной некоторой достаточно гладкой функции fix) на отрезке [хк,хш] справедливы следующие оценки погрешности [66]

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

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

N. Искомыми величинами являются вектор-функция (U,R) в точке 3 шаблона. Поэтому мы получаем 8N уравнений с 8N неизвестными и таким образом определяем все искомые величины сразу по всему венку.

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

Расчет искомых функций в точке на обтекаемом теле На линии пересечения плоскостей xs и cpi введем точки 1 и 2, а на линии пересечения xs+l и cpt -точку 3 (см. Рис. 3.3).

Шаблон для расчета искомых функций на обтекаемом теле в некоторой меридиональной плоскости Координаты точек 1, 2 и 3 задаются по формулам (3.1) при условии, что jx = О, j2 = 1 и Уз = 0. Таким образом,

Шаблон для расчета искомых функций на обтекаемом теле состоит из пересечения псевдохарактеристических поверхностей первого и второго семейств и меридиональной плоскости р = %. Поэтому, соединив точки 1 -3 и 2-3, получим минимальный шаблон из трех точек (см Рис. 3.3). По известным газодинамическим величинам в точках 1 и 2 необходимо определить искомые газодинамические величины в фиксированной точке 3 шаблона. Система дифференциальных уравнений, описывающие течение газа в точке на поверхности тела, совпадают с системой (3.6). где пт - вектор единичной нормали к поверхности тела.

Таким образом, мы имеем систему из 9 уравнений для определения 8 искомых величин (U, R). Разностная система уравнений для нахождения искомых функций на обтекаемом теле совпадает с (3.7). Записывая вместо одного из уравнений системы (3.7) условие непротекания (2.10), получим модифицированную систему из 8N уравнений с 8N неизвестными, которую также решаем с помощью итераций. Опыт расчетов показал, что наилучшая точность достигается при замене первого уравнения вдоль псевдохарактеристики 1-3.

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

Сложность нахождения искомых величин в точке на ударной волне заключается в том, что угловой коэффициент наклона — и положение ударной волны rj{xs+l,(pt) в искомой точке заранее неизвестны и также подлежат определению. Поэтому, расчет углового коэффициента наклона и искомых газодинамических функций в точке на ударной волне проводился не с помощью метода псевдохарактеристик, а с помощью итераций как в методе характеристик. Подробный вывод формул для определения искомых функций в точке на ударной волне можно найти в литературе [см., например 37 ]. Приведем здесь лишь основные результаты и логику вычисления искомых функций в точке на ударной волне.

Итак, в качестве начального приближения для угла наклона ударной волны в плоскости xs+1 принимается угол наклона в плоскости с известными данными xs. По известному углу наклона ударной волны с помощью соотношений Ренкина-Гюгонию (2.13) и (2.16) определяются искомые величины в точке 3 на ударной волне. Ниже на Рис. 3.4 изображен шаблон для расчета искомых величин в точке 3 на ударной волне в некоторой меридиональной плоскости

Пространственное обтекание

Поскольку коэффициенты при старшем и младшем членах уравнения (4.10) равны единице, то ХЪХ6? \ = 1. Отметим, что аналогичное условие имеет место при расчете методом псевдохарактеристик стационарных и нестационарных сверхзвуковых течений в декартовых координатах [69, 70]. По-видимому, оно отражает общее свойство метода псевдохарактеристик [70].

Найдем собственные числа уравнения (4.10). Решая данное уравнение с помощью замены у = Л + —, получим Таким образом, мы получили значения для оставшихся четырех собственных чисел. Определим соотношения между параметрами течения, при которых выполняется необходимое условие устойчивости (4.8).

Для дальнейшего удобства введем следующие коэффициенты 0 и к2 О, связывающие шаги по пространственным переменным между собой:

Рассмотрим течение газа в сверхзвуковой области, то есть такой, что в любой точке области выполняется неравенство —г -\=М2лок 1, где Млок а а0 локальное число Маха. Значения собственных чисел, являющихся решениями уравнения (4.10) при различных значениях входных параметров, вычислялись на компьютере. Анализ полученных результатов позволил сделать следующие выводы относительно устойчивости метода псевдхарактеристик:

1. Метод псевдохарактеристик удовлетворяет условию устойчивости при выполнении определенных ограничений на шаги hx и hr. При выполнении этих ограничений, дополнительные ограничения на шаги по h9 не требуются.

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

Графически это означает, что область устойчивости метода псевдохарактеристик находится под графиком k2msK(Мжж). Ниже на Рис. 4.2 показан график к2тах(Млок) ПРИ 0(hx) = 0,05. «Квадратами» указаны значения k2max 5 вычисленные автором с помощью компьютера.

Исследование устойчивости метода псевдохарактеристик было проведено с помощью аппарата спектральной устойчивости. Спектральный признак оценивает устойчивость только по начальным данным. Хотя данный аппарат не дает точного ответа на вопрос об устойчивости, тем не менее схемы, признанные на его основе устойчивыми, как правило, на самом деле являются таковыми, и в более или менее общем случае из такой устойчивости следует устойчивость по правой части [71]. Анализ устойчивости алгоритма по краевым условиям выходит за рамки данной диссертации.

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

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

Газ будем считать идеальным, то есть k = 1,4 . Угол атаки набегающего потока равен нулю, число Маха набегающего потока Мю равно 2, Мю = 2. Начальные данные брались из таблиц [4] на линии х0 = 1,5. Число венков J = 9, шаг по маршевой переменной hx=0,\. Ниже на рисунках и таблицах представлены расчетные и табличные значения положения ударной волны, тангенса угла наклона ударной волны к оси х, а также значения искомых газодинамических функций в точке на ударной волне, во внутренней точке области течения и в точке на обтекаемом теле. Результаты осесимметричного обтекания были представлены автором на различных конференциях и научных семинарах [73-75].

Газ будем считать идеальным, то есть к = 1,4 . Число Маха набегающего потока Мю равно 4, Мю = 4. Начальные данные брались из таблиц [4] на плоскости х0 = 2. Число венков J = 9, число меридиональных плоскостей N = 10, шаг по маршевой переменной hx=0,\. Расчеты проводились при двух различных углах атаки набегающего потока: 0 и 5. При расчете внешнего обтекания эллиптического параболоида сверхзвуковым потоком стационарного газа методом псевдохарактеристик были получены газодинамические функции, которые представлены в нижеследующих таблицах и рисунках. Результаты пространственного обтекания были опубликованы автором в работе [76].