Метод псевдохарактеристик для расчета математической модели сверхзвукового пространственного течения газа тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Нурбаев, Улукбек Джаныбекович

• Нурбаев, Улукбек Джаныбекович
• кандидат наук
• 2017

Оглавление диссертации кандидат наук Нурбаев, Улукбек Джаныбекович

Содержание

Введение

Глава 1. Моделирование обтекания пространственных тел сверхзвуковым газом

1.1 Математическая модель движения газа, описываемая уравнениями Эйлера

1.1.1 Цилиндрическая система координат

1.1.2 Декартовая система координат

1.2 Начальные данные и граничные условия

1.3 Выводы из главы 1

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

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

2.2 Введение в метод псевдохарактеристик

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

2.3.1 Расчетная сетка

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

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

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

2.3.5 Алгоритм расчета методом псевдохарактеристик

2.4 Выводы из главы 2

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

3.1 Линейное приближение модели обтекания пространственных тел сверхзвуковым газом

3.2 Анализ устойчивости метода псевдохарактеристик

3.3 Выводы из главы 3

Глава 4. Примеры расчетов типовых тел

2

4.1 Конус

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

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

4.2 Круговой параболоид

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

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

4.3 Эллиптический параболоид (а1=3, а2=5)

4.4 Эллиптический параболоид (а1=5, а2=3)

4.5 Выводы из главы 4

Заключение

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

Приложение А. Акт внедрения (копия)

Приложение Б. Свидетельство о государственной регистрации программы для

ЭВМ

Приложение В. Структура и код программы метода псевдохарактеристик

Введение диссертации (часть автореферата) на тему «Метод псевдохарактеристик для расчета математической модели сверхзвукового пространственного течения газа»

Введение

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

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

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

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

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

В соответствии с поставленной целью были сформулированы следующие задачи:

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

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

3. Тестирование предложенного метода на расчетах сверхзвуковых течений газа около различных пространственных тел;

4. Разработка комплекса программ на основе предложенного метода для расчета внешнего сверхзвукового обтекания пространственных тел. Научная новизна результатов. В диссертационной работе предлагается

новый маршевый численный метод для расчета нелинейных уравнений в частных

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

Предложенный метод имеет регулярную расчетную сетку по маршевой

переменной, второй порядок аппроксимации и не требует промежуточных

интерполяций. Кроме того, в отличие от характеристических методов он не

использует характеристическую поверхность тока и условий совместности на ней.

На модели линеаризованного течения сверхзвукового газа в цилиндрических

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

по пространственным переменным, обеспечивающие устойчивость метода.

Дополнительно показано, что метод сохраняет устойчивость на оси симметрии.

Алгоритм метода реализован в виде комплекса программ, написанного на языке

программирования Ма^аЬ. Комплекс программ адаптирован к расчету обтекания

5

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

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

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

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

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

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

2. Ограничения на шаги по пространственным переменным, обеспечивающие устойчивость метода. Сохранение устойчивости метода на оси симметрии.

3. Комплекс программ, адаптированный к расчету обтекания стационарным сверхзвуковым потоком газа пространственных тел под различными углами атаки с помощью метода псевдохарактеристик. Комплекс программ написан с использованием пакета прикладных программ Ма^аЬ.

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

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

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

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

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

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

Соответствие диссертации паспорту специальности. Диссертационная работа соответствует научной специальности 05.13.18 «Математическое моделирование, численные методы и комплексы программ» в части п.п. 2, 3 и 4 области исследования паспорта специальности.

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

Структура и объем диссертации. Диссертационная работа состоит из введения, четырех глав, заключения, списка литературы и трех приложений. Работа изложена на 136 страницах, включает 28 таблиц, 63 рисунка и 99 наименований использованной литературы.

Глава 1. Моделирование обтекания пространственных тел сверхзвуковым

газом

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

Эйлера

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

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

Рисунок 1.1 - Внешнее обтекание пространственных тел сверхзвуковым газом

г

В дозвуковой области течение газа описывается уравнениями в частных производных эллиптического типа, а в сверхзвуковой - гиперболического типа [2].

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

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

-5

неприемлем не только из-за очень большого их числа (в 1см воздуха при нормальных условиях содержится 2,7-1019 молекул), но также ввиду невозможности указать точные начальные данные. Поэтому в математической модели движения газа используется осредненное описание движения и взаимодействия. Тем не менее, такая математическая модель всесторонне апробирована практикой и в большинстве случаев с достаточной для практики точностью описывает обтекание тел и параметры газа. Причиной совпадения математической модели и опыта служит дополнительное допущение, связанное с непрерывным распределеннием средних величин по занимаемому объему (принцип непрерывности) [1].

Принцип непрерывности приходится нарушать в особых случаях, например

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

или, другими словами, отсутствие протекания газа сквозь поверхность тела.

Однако на практике газ не скользит вдоль тела, скорость частиц, граничащих с

телом, равна нулю и газ прилипает к телу. Тем не менее, эта скорость резко

возрастает при удалении от поверхности тела и на границе довольно тонкого

пограничного слоя достигает значений, соответствующих модели непротекания

10

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

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

1.1.1 Цилиндрическая система координат

Будем рассматривать в качестве обтекаемых тел - тела плавной вытянутой формы с затупленной носовой частью и/или конусообразной формы (параболоиды вращения, конусы, фюзеляж самолета и др.). Для таких тел «естественной» является цилиндрическая система координат.

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

— — —

векторов (е*, er, eq>), причем под маршевой переменной будем принимать переменную х. В системе координат, движущейся вместе с телом, последнее неподвижно. Газ будем считать идеальным, то есть уравнение состояние газа имеет вид [1]:

Р

-V = const рк

где p - давление газа, р - плотность газа, к - показатель адиабаты.

Исходная система уравнений математической модели, описывающей течение идеального газа, представляет собой законы сохранения массы, импульса и энергии (записанные в дифференциальной форме), и выглядит следующим образом [1]:

pt+ div(pV) = 0

V+(V-V)v+— = о p

st + (V -V)s + — div(V) = 0 P

где s - удельная внутренняя энергия газа, V = (u, v, w)' - скорость газа (верхний индекс штрих - знак транспонирования), V - оператор Гамильтона.

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

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

™ = о ф= о, ф = о, ^ = о dt dt dt dt

Следствием уравнений сохранения импульсов и энергии при стационарном течении газа является интеграл Бернулли [1], который имеет вид:

V2 k p

--1---= const

2 k-1 p

В задачах обтекания константа, стоящая в правой части интеграла Бернулли, независит от траектории движения частиц газа и поэтому одинакова во всей области течения. Такие течения называются изоэнергетическими [3].

Систему уравнений Эйлера, описывающую движение газа, можно записать в следующем виде:

сИу(р\) = 0

(V-У)У + ^ = 0 (1.1)

р

2 к-1 р

V2 к р

- н---= еож1

Искомыми функциями являются пять газодинамических функций: скорость V, давление р и плотность газа р. При этом система (1.1) содержит пять уравнений и является замкнутой системой уравнений.

Раскроем скобки и приведем подобные в системе (1.1). Используя известные соотношения в цилиндрической системе координат:

div(рV) = V(Vр) + р ■ (Иу(У)

VР = РХ +РГ +1 Рр г

V 1

^(У) = их + V +- + - * г г

и учитывая, что

дех _ дх _дех дг = дех = 0 дф

дег _ к м де = 0, г = дф

дх дг

деФ _ дх II де = 0, ф = дф

—е

г

запишем первое и второе уравнения системы (1.1) в проекциях на оси (х,г,ф):

1 Р Р

иРх + VРr+-*Р<р+Рих +РЛ?г + — * ф =--V

г г г

* Рх <л

иих + vur +—иф+— = 0

г ф р

* Р, 2 (!.2)

^ + vv н--V н—- = *

Х г ф

г р

* 1 рф V ■ * и* + V* н— * н---— =--

Х г ф

г г р г

ф

Дифференцируя интеграл Бернулли из системы (1.1), получим выражение для полного дифференциала плотности газа в зависимости от его скорости и давления

VdV + -L Рdp-Р ' dP = 0 k -1 р2

dp= (k - 1)Р VdV + Аф (1.3)

а а

2 1 Р

где а = - местная скорость звука. Р

В проекциях на оси (х,г,ф) выражение (1.3) принимает вид

к -1 .к

Рх =~Г Р(UUx + Wx + ) + — Рх а а

к -1 к

рг Р(ииг + Wr + wwr ) + ~ Рг (14)

п п V • /

а а

к -1 *

Р?> = —г Р(UU^ + vvV + wwV) + ~ Р^

а а

Подставив (1.4) в первое уравнение системы (1.2) и используя остальные выражения в системе (1.2), получим окончательную запись уравнения неразрывности:

и2. UV uw v2. vu vw

(1--r)Ux--г vx--Г wx + (1--Г )vr--Г Ur--Г wr +

а а а а а а (15)

,, w2. 1 WU 1 wv 1 v

+(1--- ---Г ~UV--Г -v^= —

а гаг а r r

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

„ Р Р v V

рб = —, Рб = —, V6 = / .

Р» Р» ^Р» / Рос

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

V2

кРб

1

2 (к — 1)ре рх / рх

V кр 2 (к —1)р

р„ / рх

крх

У2

2 н (к—1)рК

У2

2(рт /рт) (к — 1) 2

к 2 к ■ = — М2 н -

к—1

где Мш - число Маха набегающего потока.

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

аб = к

1 н—М2

к — 1

(1.6)

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

Таким образом, уравнение (1.5) и второе, третье и четвертое уравнения системы (1.2) формируют систему из четырех уравнений для нахождения четырех газодинамических величин: скорости V = (и,V,*)' и давления р. Данная система дифференциальных уравнений в частных производных выглядит следующим образом:

и2. Ш и* /1 V2 . vu V*

(1--2)их--2 V--2 *х н (1--б)vг--2 иг--2 *г н

а а а а а а

2 '"ф 2 иф 2 'ф а г а г а г г

.. 1 *и 1 wv 1 V

н(1--7 )*ф---5Т ~иф--2" " ^ = — "

* рх п

иих н vuг н — ифн — = 0 г р

* рг 2 ш н vv н--V н—- = *

х г ф

г р

* 1 рф V ■ *

и*х нvwг н — * = ■

г г р г

(1.7)

Систему (1.7) можно представить в матричном виде:

ли + ви + си = f

(1.8)

где А, В, С - квадратные матрицы размера 4х4, и - вектор искомых переменных, и, и и и - частные производные по соответствующей координате.

л =

( 2 - ^ а иу uw Л 0 г иу 1 - vw Л 0

- а2 - а2 а2 а2

кр 2 и 0 0 1 кр 2 У 0 0 0

а2 0 кр ~и а 0 0 , в = а2 0 кр а 0 1

0 0 кр ^и 0 0 0 кр 2 V 0

V а ) V а )

с = -

и4 vw . - 4 а \ 0

- а - а

кр — w а 0 0 0

0 кр — w а 0 0

0 0 кр — 4 а 1 )

и=

Г V

[ и Л Г 0

V , f = Р42

V р У Г

Ур4

v г у

Плотность газа р определяется через остальные газодинамические величины

2 , Р

через выражения для местной скорости звука а = к—, полученной для интеграла

Р

Бернулли (1.6).

1.1.2 Декартовая система координат

На практике также часто встречаются обтекаемые тела, для которых «естественной» является декартовая система координат (например, крыло

самолета). Рассмотрим декартовую систему координат (х,у,г). Учитывая, что в декартовой системе координат

VР = РХ + РУ +Рг ¿НЮ = их + Уу + м

используя следующие выражения

де де„

де„

дх ду дг

де де де

у у у

дх ду дг

дег _ _ дег _ _ дег

дх ду дг

0 = 0 0

и выполняя те выкладки, которые были указаны для цилиндрической системы координат, систему (1.1) можно представить следующей системой в матричном виде

Аи +Ви +Си =0

где А, В, С - квадратные матрицы размера 4х4, и - вектор искомых переменных, и, И и И - частные производные по соответствующей координате. Вид матриц А, В и С представлен ниже:

А =

1 м2 IV им

- а2 - а

кр 2 11 0 0

а2

0 кр —и а2 0

0 0 кр —и а2

0

В

IV 1 - а2 УМ

а2 а2

кр 2 V 0 0

а2

0 кр —у а 0

0 0 кр — У а

0

V 0

C

Г 2 л

uw vw „ w

-— - — 1 - — O a a a

kP n —^r w O

a

kp

OO

O -Z-w O O a

O O kpw 1

a

U=

u

v w

v P J

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

a2 = kp, полученной для интеграла Бернулли (1.6). Р

1.2 Начальные данные и граничные условия

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

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

одной точке. Угол полураствора конуса Маха определяется уравнением sin а = a,

где а - угол полураствора, a - скорость звука и V - скорость газа. Отметим, что если конус Маха пересекает поверхность в нескольких точках, то такая поверхность называется времениподобной. Если задавать начальные данные на таких поверхностях, то задача является некорректной [98].

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

(V, п„) = 0

(1.10)

и на ударной волне условия Ренкина-Гюгонио (в размерных величинах):

^ = К, РзУзп = РхКп

Р3 + Рз^п =

к р 1 тг2 к п 1 2 -— + — У3 =-^ + — V

к -1 Рз 2 п к -1 рх 2 оои

3п

(111)

где V3 - скорость в точке непосредственно за ударной волной, V ^ - скорость в точке непосредственно перед ударной волной, V«, , V, и - соответственно проекция векторов на нормаль п к ударной волне и касательный вектор t к ударной волне, пи - вектор единичной нормали к поверхности тела.

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

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

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

Vзt=^¡kMaa

РзУЪп =4кМх пх

Рз + РУ1 = 1 + кМТО Пх2 (112)

к Рз ,1 т^ _ к , 1 7-»Г2„2

+ - У 2 =-+ - кМ2 п2

к -1 Рз 2 3п к -12 " х

где пх - проекция вектора нормали п на единичный орт ех оси х.

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

Рз 4кМя пх

У3п

рз = 1 + км> Пх2 -4кмвпхузп

Подставив выражения для давления и плотности в четвертое уравнение системы (1.12) и приводя подобные члены, получим следующее квадратное уравнение для определения нормальной составляющей вектора скорости У3п:

12 к 1+^/кМл)\/ + _п

— У 1„---;=-+--_ 0

2 зп к +1 4кМ п зп к +1

* то х

Уравнение имеет два корня. Первый корень уравнения уя =4кМагпх не удовлетворяет условию переходу газа через ударную волну, так как он предполагает сохранение нормальной составляющей газа, а не ее изменение [1]. Второй корень уравнения имеет вид:

ТЛ 2к 1 к - 1 г-Л/Г

Узп -+ Т~Т ^кМ » Пх

к +1 V кМгхпх к +1

Подставляя второй корень уравнения в систему (1.12) и приводя подобные слагаемые, получим следующие выражения для газодинамических величин после ударной волны:

V = ———+^4км п 3п к+14кмп к+\

(1.13)

1 - к 2

(#МЛ )2

Р3 =-+-

3 к +1 к +1

Р3 =

Из системы (1.13) видно, что газодинамические величины за ударной волной зависят не только от соответствующих значений набегающего потока, но и от положения самой ударной волны в пространстве. Система (1.13) позволяет определить искомые газодинамические величины в точке за поверхностью ударной волны: нормальную составляющую вектора скорости ^, касательную составляющую вектора скорости , давление р3 и плотность р3.

Выведем уравнения для определения значений проекций скорости V - и, V, w - в точке за поверхностью ударной волны на оси координат. Рассмотрим точку 3 на поверхности ударной волны (см. Рисунок 1.2). Из этой точки выпустим три вектора - единичную внутреннюю нормаль п к поверхности ударной волны УВ и два перпендикулярных ей вектора 1 и т, лежащие на поверхности ударной волны.

Рисунок 1.2 - Поверхность ударной волны

Положение векторов задается с помощью следующих формул

/

1' =

Ыт

\

1; л;0,

у ЫЛ J

г дх _ дт л дф дф

П ' = 8[1X т]

(1.14)

где 8 - нормировочный коэффициент и такой, что вектор нормали п является единичным.

Введем вектор 1 = а\ + [5т. Тогда любой вектор в точке 3 на поверхности ударной волны можно представить как сумму векторов п и 1. В частности, для вектора скорости V имеем следующее соотношение:

V = Vn П + V1 = (V, п)П + (V, 1)1

(1.15)

откуда можно выразить вектор 1

^ V - П(У, П) IV»- п(Ую , п)|

где V,' = (иш, о, 0) = (у[кмх, о, 0) при нулевом угле атаки. Из соотношения (1.15) также следует, что

и3ех + ^т + ^ = У3 = + V3f 1

Приравнивая коэффициенты при ортах е^е^е^,, получим выражения для нахождения проекций скорости V в точке 3 на поверхности ударной волны.

и3 = + УъК

V3 = ^ + У3К

= У3п% + ^

(1.16)

Таким образом, дифференциальная система в матричном виде (1.8), граничные условия на теле (1.10), ударной волне (1.13) и (1.16) и начальные данные на пространственноподобной поверхности П (1.9) позволяют определить искомые величины во всей рассматриваемой сверхзвуковой области течения газа.

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

1.3 Выводы из главы 1

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

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

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

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

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

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

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

Список литературы диссертационного исследования кандидат наук Нурбаев, Улукбек Джаныбекович, 2017 год

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

1. Овсянников, Л.В. Лекции по основам газовой динамики / Л.В. Овсянников: Москва-Ижевск. Институт компьютерных исследований. 2003. 368 с.

2. Курант, Р. Уравнения с частными производными / Р. Курант: М. Мир. 1964. 823 с.

3. Лойцянский, Л.Г. Механика жидкости и газа / Л.Г. Лойцянский: Москва -Ленинград. Государственное издательство технико-теоретической литературы. 1950. 678 с.

4. Бабенко, К.И. Численный метод расчета пространственного обтекания тел сверхзвуковым потоком газа / К.И. Бабенко, Г.П. Воскресенский // Ж. вычисл. матем. и матем. физ.. 1961. Т.1. № 6. С.1051-1060.

5. Бабенко, К.И. Пространственное обтекание гладких тел идеальным газом / К.И. Бабенко, Г.П. Воскресенский, А.Н. Любимов, В.В. Русанов М.: Наука. 1964 г. 505 с.

6. Любимов, А.Н. Течение газа около тупых тел /А.Н. Любимов, В.В. Русанов. М.: Наука. 1970. Т.1. Метод расчета и анализ течений. 287с.

7. Русанов, В.В. Пространственное обтекание затупленного тела сверхзвуковым потоком газа / В.В. Русанов // Ж. вычисл. матем. и матем. физ.. 1968. Т.8. № 3. С.616-633.

8. Любимов, А.Н. Течение газа около тупых тел /А.Н. Любимов, В.В. Русанов. М.: Наука. 1970. Т.2. Таблицы газодинамических функций. 379 с.

9. Дьяконов, Ю.Н. Сверхзвуковое обтекание затупленных конусов и тел сложной формы / Ю.Н. Дьяконов, Н.С. Кокошинская, В.М. Лукашин и др.: Некоторые применения метода сеток в газ.динамике. Вып. V. Обтекание затупленных тел газом. М. Моск. ун-т. 1974. 4. 139.

10.Ткаченко, А.С. Численное исследование тяговых характеристик и структура пространственных течений в соплах / А.С. Ткаченко // Изв. АН СССР. МЖГ. 1981. № 5. С. 168-172.

11.Иванова, В.П. Трехмерное сверхзвуковое обтекание гладких тел неравновесно реагирующим воздухом / В.П. Иванова, Ю.Б. Радвогин: Препринт № 92. М.: Ин-т прикл. матем. им. М. В. Келдыша. АН СССР. 1981.

12. Бородин, А. И. Численное исследование сверхзвукового обтекания затупленных тел сложной формы под углами атаки и скольжения / А. И. Бородин, С. В. Пейгин // ТВТ. 38:3. 2000. С.468-476

13. Леонтьева, Н. В. Численное моделирование сверхзвукового асимметричного обтекания острого кругового конуса / Н. В. Леонтьева, Ю. П. Головачев // Ж. вычисл. матем. и матем. физ. Т.40. №4. 2000. С.638-646

14. Шевелев, Ю. Д. О влиянии многокомпонентной диффузии на сверхзвуковое обтекание тел потоком углекислого газа / Ю. Д. Шевелев, Н. Г. Сызранова // ТВТ. 47:6. 2009. С.863-869

15.Годунов, С.К.Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной / С.К. Годунов, А.В. Забродин, Г.П. Прокопов // Ж. вычисл. матем. и матем. физ.. 1961. Т.1. № 6. С.1020-1050.

16.Годунов, С.К. Численное решение многомерных задач газовой динамики / С.К. Годунов, А.В. Забродин, М.Я. Иванов, А.Н. Крайко, Г.П. Прокопов. М.: Наука. 1976. 400 с.

17.Иванов, М.Я. Метод сквозного счета двумерных и пространственных течений / М.Я. Иванов, А.Н. Крайко, Н.В. Михайлов // Ж. вычисл. матем. и матем. физ.. 1972. Т.12. № 2. С.441-463.

18.Иванов, М.Я. Метод сквозного счета двумерных и пространственных течений / М.Я. Иванов, А.Н. Крайко // Ж. вычисл. матем. и матем. физ.. 1972. Т.12. № 3. С.805-813.

19. Родионов, А.В. Численный метод решения уравнения Эйлера с сохранением аппроксимации на деформируемой сетке / А.В. Родионов // Ж. вычисл. матем. и матем. физ.. 1996. Т.36. № 3. С.117-129.

20. Шевелев, Ю.Д. Сверхзвуковое обтекание тел вязким теплопроводным газом при разделении ступеней ракет / Ю.Д. Шевелев, Н.А. Михалин, Н.Г. Сызранова // Матем. моделирование. 2002. Т.14. № 7. С.3-14.

21.Кузьмина, А. В. Расчет полей течения на входе в воздухозаборник сверхзвукового самолета (СПС) маршевой схемой второго порядка аппроксимации / А. В. Кузьмина // Матем. моделирование. 2000. Т.12. № 6. С.79-84

22.Бабарыкин, К. В. Автоколебательные процессы при обтекании сверхзвуковым потоком тел с образованием передней срывной зоны / К. В. Бабарыкин, В. Е. Кузьмина, Е. А. Угрюмов, А. И. Цветков // Матем. моделирование. 2003. Т.15. №6. С. 27-32

23.Луцкий, А. Е. Влияние вязкости на параметры сверхзвукового обтекания профиля крыла / А. Е. Луцкий, А. С. Черногузов: Препринты ИПМ им. М. В. Келдыша. 2004. 026. 27 с.

24.Васильев, Е.И. Трехмерное обобщение для W-модификации метода Годунова / Е.И. Васильев, А.С. Демин // Ж. вычисл. матем. и матем. физ.. 2008. Т.48. № 9. С.1659-1672.

25.Глушко, Г. С. Метод расчета турбулентных сверхзвуковых течений / Г. С. Глушко, И. Э. Иванов, И. А. Крюков // Матем. моделирование. 2009. Т.21. №12. С. 103-121

26.Пахомов, Ф. М. Теоретико-экспериментальное исследование сверхзвукового обтекания тел выпукло-вогнутой конфигурации / Ф. М. Пахомов, В. А. Антонов, Г. Ф. Костин, Н. В. Чурилов // Вестн. Томск.гос. ун-та. Матем. и мех. 2009. № 4(8). С.93-97

27. Голованов, А. Н. Сверхзвуковое обтекание затупленного тела при наличии вдува газа-охладителя и продольных потоку колебаний стенки / А. Н. Голованов, В. Д. Гольдин // Вестн. Томск.гос. ун-та. Матем. и мех. 2010. № 1(9). С. 87-95

28.Lax, P.D. Difference schemes for hyperbolic equations with high order of accuracy / P.D. Lax, B.Wendroff // Communs Pure and Appl. Math. 1964. 17. № 3. P. 381-398.

29.MacCormack, R.W. The effect of viscosity in hypervelocity impact cratering / R.W. MacCormack // AIAAPaper. 1969.№ 364.

30.Погорелов, Н.В. Маршевый явно-неявный метод расчета сверхзвукового обтекания тел / Н.В. Погорелов, Ю.Д. Шевелев // Ж. вычисл. матем. и матем. физ.. 1985. Т.25. № 9. С.1391-1400.

31.Глушко, А. И. Об одном численном методе решения смешанных задач для гиперболических систем уравнений на нерегулярной сетке / А. И. Глушко, И. И. Нещеретов // Ж. вычисл. матем. и матем. физ. 2001. Т. 41. №1. С. 114125

32. Куликовский, А.Г. Математические вопросы численного решения гиперболических систем уравнений / А.Г. Куликовский, Н.В. Погорелов, А.Ю. Семенов. М.: Физматлит. 2001. 608 с.

33.Георгиевский, П. Ю. Режимы сверхзвукового обтекания тел в условиях энерговклада различной геометрической конфигурации в набегающий поток / П. Ю. Георгиевский, В. А. Левин // ТВТ. 48. дополнительный выпуск. 2010. С. 77-84

34. Максимов, Ф. А. Сверхзвуковое течение в осесимметричном канале / Ф. А. Максимов // Мат. моделир. и числ. Методы. 2015. № 5. С. 109-120

35.Яненко, Н.Н. Численное решение задач механики жидкости / Н.Н. Яненко // Тр. III Всес. семинара по моделям мех. сплош. среды. Новосибирск. 1976. С. 177-199.

36.Яненко, Н.Н. Метод дробных шагов решения многомерных задач математической физики /Н.Н. Яненко: Наука. Новосибирск. 1967. 195с.

37.Evans, M.W. Calculation of unsteady supersonic flow past a circular cylinder / M.W. Evans, F. H. Harlow //ARS Journal. 1959. 29. № 1. P.46-48

38.Белоцерковский, О.М. Численное исследование современных задач газовой динамики / О.М. Белоцерковский, Ю.П. Головачев, В.Г. Грудницкий, Ю.М. Давыдов и др. М.: Наука, 1974. 398 с.

39.Белоцерковский, О.М. Метод крупных частиц в газовой динамики: Вычислительный эксперимент / О.М. Белоцерковский, Ю.М. Давыдов М.: Наука. 1982. 392 с.

40.Ковеня, В.М. Методы расщепления в задачах газовой динамики / В.М. Ковеня, Н.Н. Яненко: Наука. Новосибирск. 1981. 304 с.

41.Карпенко, А. П. Численное моделирование нестационарного трехмерного течения газа с ударными волнами и отрывом потока от поверхности / А. П. Карпенко, В. Н. Ляхов, И. Н. Протасов, В. Е. Фортов // Матем. Моделирование. 7:8. 1995. С. 36-59

42.Ковеня, В. М., Алгоритмы оптимального расщепления в задачах аэро и гидродинамики / В.М. Ковеня // Матем. моделирование. 16:6. 2004. С. 23-27

43.Ковеня, В. М. Метод предиктор-корректор для численного решения уравнений Эйлера и Навье-Стокса / В. М. Ковеня, А. А. Еремин // Вестн. НГУ. Сер. матем., мех., информ. 15:2. 2015. С.22-37

44.Подладчиков, Ю. Н. Метод характеристик для расчета пространственных сверхзвуковых течений газа / Ю. Н. Подладчиков // Изв. АН СССР. Механика. 1964. №4. С. 3-12.

45.Борисов, В.М. Численный метод характеристик для трехмерных стационарных течений газа / В.М. Борисов, И.Е. Михайлов // Сб. теорет. работ.погидромех. М.: ВЦ АН СССР. 1970. С.6-29.

46.Борисов, В.М. Расчет установившихся сверхзвуковых пространственных течений с головной ударной волной методом характеристик / В.М. Борисов, Ю.В. Куриленко // Ж. вычисл. матем. и матем. физ.. 1998. Т.38. № 4. С. 659668.

47.Ахметов, Г.Б. Об одной схеме численного метода пространственных

характеристик для расчета безвихревых течений газа / Г.Б. Ахметов, И.Е.

Михайлов // Изв. АН КазСССР. Сер. физ.-матем.1975. №3. С.29-35.

112

48.Кацкова, О.Н. Об одной схеме численного метода характеристик / О.Н. Кацкова, П.И. Чушкин //Докл. АН СССР. 1964. 154. №1. С. 26-29.

49.Магомедов, К.М. Метод характеристик для численного расчета пространственных течений газа / К.М. Магомедов // Ж. вычисл. матем. и матем. физ. 1966. Т.6. № 2. С.313-325.

50.Чушкин, П.И. Метод характеристик для пространственных сверхзвуковых течений / Чушкин П.И.: Тр. ВЦ АН СССР. 1968. 121 с.

51.Кацкова, О.Н. Трехмерное сверхзвуковое равновесное течение газа около тел под углом атаки / О.Н. Кацкова, П.И. Чушкин // Ж. вычисл. матем. и матем. физ. 1965. Т. 5. №3. С. 503-518.

52.Григорьев, Е.И. Об одной прямой схеме метода характеристик для расчета пространственных течений газа / Е.И. Григорьев, К.М. Магомедов // Ж. вычисл. матем. и матем. физ. 1969. Т.9. № 6. С. 1413—1419.

53.Борисов, В.М. Метод характеристик для расчета пространственных сверхзвуковых безвихревых течений / В.М. Борисов, И.Е. Михайлов // Ж. вычисл. матем. и матем. физ. 1978. Т.18. № 5. С.1243-1255.

54.Крайко, А. Н. К построению головной ударной волны при обратном расчете сверхзвукового течения методом характеристик / А. Н. Крайко, В.Е. Макаров, Д.Е. Пудовиков // Ж. вычисл. матем. и матем. физ. 1999. Т.39. №11. С.1889-1894

55.Крайко, А. Н. Метод характеристик и полухарактеристические переменные в задачах профилирования сверхзвуковых частей осесимметричных и плоских сопел / А. Н. Крайко, Н. И. Тилляева // Ж. вычисл. матем. и матем. физ. 1996. Т.36. №9. С.159-160

56.Магомедов, К.М О построении разностных схем для уравнений гиперболического типа на основе характеристических соотношений / К.М. Магомедов, А.С. Холодов // Ж. вычисл. матем. и матем. физ. 1969. Т.9. №2. С.373-386.

57.Магомедов, К.М. Сеточно-характеристический метод для решения задач

газовой динамики / К.М. Магомедов // Тр. II Междунар. коллоквиума по

113

газовой динамике взрыва и реагирующих систем. Новосибирск. 1969. 19-23 августа. М.: ВЦ АН СССР. 1971

58.Магомедов, К.М. Сеточно-характеристические численные методы / К.М. Магомедов, А.С. Холодов: М.: Наука, 1988. 290 с.

59.Дородницын, А.А. Об одном методе численного решения некоторых нелинейных задач аэрогидродинамики / А.А. Дородницын // Тр. III Всес. матем. съезда. 1956 г. Т.Ш.М. Изд-во АН СССР. 1958. С.447-453.

60.Белоцерковский, О.М. Численный метод интегральных соотношений / О.М. Белоцерковский, П.И. Чушкин // Ж. вычисл. матем. и матем. физ. 1962. Т.2. №5. С.731-759.

61. Чушкин, П. И. Метод интегральных соотношений для сверхзвуковых пространственных течений / П. И. Чушкин // Ж. вычисл. матем. и матем. физ. 1968. Т.8. №4. С.853-864

62.Белоцерковский, О. М. Исследование нестационарных течений газа со сложной внутренней структурой методами интегральных соотношений / О. М. Белоцерковский, В. Г. Грудницкий // Ж. вычисл. матем. и матем. физ. 1980. Т.20. №6. С. 1400-1415

63.Белоцерковский, О.М. Исследование влияния сильного вдува газа с поверхности на гиперзвуковое обтекание тел / О.М. Белоцерковский, М.М. Голомазов, А.В. Шабалин // Ж. вычисл. матем. и матем. физ.. 1981. Т.21. №4. С.1010-1030.

64. Голомазов, М.М. Исследование сверхзвукового обтекания сегментальных и конических тел с учетом химических реакций / М.М. Голомазов // Ж. вычисл. матем. и матем. физ. 2005. Т.45. №11. С.2070-2076.

65.Белоцерковский, О.М. Метод интегральных соотношений для расчета обтекания затупленных тел с отошедшей ударной волной / О.М. Белоцерковский, Г.А. Цветков // Ж. вычисл. матем. и матем. физ. 2011. Т.51. №1. С.120-130.

66.Pulliam, T. Implicit Finite-Difference Simulations of Three-Dimensional Compressible Flow / Т. Pulliam, J. Streger // AIAA Journal. 1980. v18. №2. P.159-167.

67.Shu, C.-W. High order ENO and WENO schemes for computational fluid dynamics in High-Order Methods for Computational Physics / Ed. By T.J. Barth and H. Deconinck. Lecture Notes in Computational Science and Engineering. v.9. Springer. 1999. P.439-582

68. Рогов, Б.В. Маршевый расчет ударной волны при невязком сверхзвуковом обтекании затупленных тел / Б. В. Рогов, И. А. Соколова // Матем. моделирование. 13:5. 2001. С. 110-118

69.Chou, C.-S. High order residual distribution conservative finite difference WENO schemes for steady state problems on non-smooth meshes / C.-S. Chou, C.-W. Shu // Journal of Computational Physics. 241. 2006. P.698-724

70. Кудрявцев, А. Н. Применение схем высокого порядка точности при моделировании нестационарных сверхзвуковых течений / А.Н. Кудрявцев, Т.В. Поплавская, Д.В. Хотяновский // Матем. моделирование. 19:7. 2007. С.39-55

71. Журавлева, Г. С. Численное моделирование турбулентных течений при сверхзвуковом обтекании тел / Г. С. Журавлева // Изв. Иркутского гос. унта. Сер. Математика. 3:3. 2010. С.67-79

72. Давыдов, А.А. Численное моделирование сверхзвукового течения в следе за крылом и его взаимодействия с пересекающимися ударными волнами / А.А. Давыдов, А.Е. Луцкий: Препринты ИПМ им. М. В. Келдыша. 2014. 098. 16 с.

73.Королев, С.А. Математическое моделирование обтекания тела вращения сверхзвуковым потоком газа / С. А. Королев, С. А. Карсканов // Вестн. Удмуртск. ун-та. Матем. Мех. Компьют. Науки. 2014. № 3. С. 123-133

74.Knox, E. A comparison of experimental and theoretically predicted pressure distributions and force and stability coefficients for a spherically blunted cons at

Мш □ 18 and angles of attack / E. Knox, C. Lewis // AEDC. 1966. 11 N TR-65-234.

75.Гридин, А.Ю. Расчетно-экспериментальное исследование сверхзвукового обтекания затупленного тела с иглой при наличии электрического разряда в его головной части / А.Ю. Гридин, Г.Б. Ефимов, А.В. Забродин: Препр. Инт проблем механики. 1995. № 19.

76.Третьяков, П.К. Управление сверхзвуковым обтеканием тел с использованием мощного оптического пульсирующего разряда / П.К. Третьяков, А.Ф. Гаранин, Г.Н. Грачев и др. // Докл. РАН. 1996. Т.31. №3. С. 339-340.

77.Пилюгин, Н.Н. Аэротермодинамические характеристики сопутствующего тела при сверхзвуковом обтекании / Н. Н. Пилюгин, В. С. Хлебников // ТВТ. 39:4. 2001. С.620-628

78.Пахомов, Ф.М. Теоретико-экспериментальное исследование сверхзвукового обтекания тел выпукло-вогнутой конфигурации / Ф.М. Пахомов, В.А. Антонов, Г.Ф. Костин, Н.В. Чурилов // Вестн. Томск.гос. ун-та. Матем. и мех. 2009. № 4(8). С.93-97

79.Аксенов, В. С. Сверхзвуковое обтекание воздухом профиля крыла при инициировании скользящего разряда на его поверхности / В.С. Аксенов, В.В. Голуб, С.А. Губин, А.С. Савельев, В.А. Сеченов, Э.Е. Сон // ТВТ. 48. дополнительный выпуск. 2010. С.93-101

80.Рябенький, В.С. Введение в вычислительную математику: Учеб. Пособие. 2-е изд., исправл. / В.С. Рябенький: М. ФИЗМАТЛИТ. 2000. 296 с.

81.Холодов, А.С. Численные методы решения уравнений и систем гиперболического типа / А.С. Холодов: М. 2007. 55 с.

82.Борисов, В.М. О методе пространственных характеристик в решении задач вычислительной аэродинамики / В.М. Борисов// Методы аэрофиз. иссл. Новосибирск: ИТПМ СО АН СССР. 1990. С. 16-23.

83.Грудницкий, В.Г. Один прием построения разностных схем с произвольным

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

116

производных / В.Г. Грудницкий, Ю.А. Прохорчук // Докл. АН СССР. 1977. Т. 234. № 6. С. 1249-1252.

84.Белоцерковский, О.М. Разностная схема второго порядка точности на минимальном шаблоне для гиперболических уравнений / О.М. Белоцерковский, В.Г. Грудницкий, Ю.А. Прохорчук // Ж. вычисл. матем. и матем. физ. 1983. Т.23. №1. С.119-126.

85.Азарова, О.А. Численные эксперименты по моделированию стационарных структур в задачах сверхзвукового обтекания с несимметричным подводом энергии / О. А. Азарова // Ж. вычисл. матем. и матем. физ. 2010. 50:10. С. 1840-1853

86.Азарова, О.А. Стационарное обтекание тел сверхзвуковым потоком газа, содержащим бесконечный тонкий разреженный канал / О. А. Азарова, В.Г. Грудницкий, Ю.Ф. Колесниченко // Матем. моделирование. 18:1. 2006. С.79-87

87.Азарова, О.А. Воздействие тонкого разреженного канала на сверхзвуковое обтекание цилиндрического тела с полостью / О. А. Азарова, Ю. Ф. Колесниченко // Матем. моделирование. 20:4. 2008. С. 27-39

88. Рождественский, Б.Л. Системы квазилинейных уравнений и их приложения к газовой динамике / Б.Л. Рождественский, Н.Н. Яненко: М.: Наука, 1979. 546 с.

89.Кацкова, О.Н. Опыт расчета плоских и осесимметричных сверхзвуковых течений газа методом характеристик / О.Н. Кацкова, Н.И. Наумова, Ю.Д. Шмыглевский, Н.П. Шулишнина: М.: ВЦ АН СССР. 1961. 59 с.

90.Годунов, С.К. Разностные схемы / С.К. Годунов, В.С. Рябенький: М.: Наука. 1977. 440 с.

91.Федоренко, Р.П. Введение в вычислительную физику / Р.П. Федоренко Учеб.пособие. М.: Изд-во Моск. физ.-техн. ин-та. 1994. 528 с.

92.Нурбаев, У.Д. Об устойчивости метода псевдохарактеристик / У.Д. Нурбаев

// Вестник Тверского Государственного Университета. Серия: Прикладная

математика. Тверь. 2014. №2. С. 21-28.

117

93.Михайлов И.Е. Метод псевдохарактеристик для расчета стационарных пространственных сверхзвуковых течений газа / И.Е. Михайлов, У.Д. Нурбаев // Теория, численные методы и математический эксперимент в газовой динамике. Материалы Научно-практического семинара. Москва. 2526 августа 2009 г. М.: ЦИАМ, 2009. С. 91 - 92.

94. Нурбаев, У.Д. Об одном методе расчета стационарных пространственных сверхзвуковых течений газа / У.Д. Нурбаев // Современные проблемы фундаментальных и прикладных наук. Часть VII. Управление и прикладная математика. Труды LIII Научной конференции МФТИ. М. Долгопрудный, 2010. С. 104 - 106

95.Mikhailov, I. On method of pseudocharacteristics for solving 3D stationary supersonic flow of ideal gas / I. Mikhailov, U. Nurbaev // VII International Aerospace Congress IAC'12. August 25-31. 2012. Moscow. Russia. Proceedings. P. 224-229. Электронный вид. Зарегистрировано в ВГУП НГЦ в ИНФОРМ-РЕГИСТР. Гос. рег. № 0321303652. 2013.

96.Михайлов, И.Е. Метод псевдохарактеристик для расчета стационарных пространственных сверхзвуковых течений газа / И.Е. Михайлов, У.Д. Нурбаев // Вестник Тверского Государственного Университета. Серия: Прикладная математика. Тверь, 2013. № 15. С.19 - 27.

97.Нурбаев, У.Д. Метод псевдохарактеристик для расчета пространственных сверхзвуковых течений газа / У.Д. Нурбаев // Свидетельство о государственной регистрации программы для ЭВМ № 2015610150 от 12 января 2015 г.

98.Михайлов, И.Е. О численном решении пространственного аналога задачи Гурса / И.Е. Михайлов // Ж. вычисл. матем. и матем. физ. 1983. 23:5. С. 1238-1240

99.Пирумов У.Г. Газовая динамика сопел / Пирумов У.Г., Г.С. Росляков: М.: Наука. Гл. ред. физ.-мат. лит, 1990. 368с

Приложение А. Акт внедрения (копия)

Приложение Б. Свидетельство о государственной регистрации программы

для ЭВМ

Приложение В. Структура и код программы метода псевдохарактеристик

Структура программы:

1. Загрузка начальных данных

2. Расчет вспомогательных > параметров /

3. Метод псевдохарактеристик

>

УВ

^ ИТ ^ ОТ ^

4. Преобразование полученных расчетных данных

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

5. Построение графиков, таблиц

УВ - ударная волна ПТ - поле течения ОТ - обтекаемое тело

Код программы метода псевдохарактеристик (в части расчета искомых функций на ударной волне, в поле течения и на обтекаемом теле):

%расчет ударной волны for it=1:Qu for fi=1:N if y==1; r0(y,it,fi)=Fx0(y,fi); end; if y>=2; ro(y,it,fi)=Fx l(y- 1,fi); end;

if ((it==1) && (y==1)); T(y,it,fi)=T0(y,1,fi); end; if ((it==1) && (y>=2)); T(y,it,fi)=T(y-1,Qu,fi); end; if it==2; T(y,it,fi)=T(y,it-1,fi)-0.001; end;

if it>=3; T(y,it,fi)=T(y,it- 1,fi)-Et(y,it- 1,fi)*(T(y,it- 1,fi)-T(y,it-2,fi))/(Et(y,it- 1,fi)-Et(y,it-2,fi));

end;

r3u(y,it,fi)=r0(y,it,fi)+T(y,it,fi)*hx; end;

%сплайн для r в 3 точке for i=2:(N-1);

if i==2; Alpha_ch_ud(i)=0; Beta_ch_ud(i)=0; end; Alpha_ch_ud(i+1)=-1/(Alpha_ch_ud(i)+4);

Beta_ch_ud(i+1)=Alpha_ch_ud(i+1)*(Beta_ch_ud(i)-3*(r3u(y,it,i+1)-r3u(y,it,i-1))/hf); end;

Mr3(N)=0; for i=N:-1:2;

Mr3(i-1)=Alpha_ch_ud(i)*Mr3(i)+Beta_ch_ud(i); % производные для r3 end;

%конец сплайна

for fi=1:N %l3

lx3=1; lr3=T(y,it,fi); lf3=0;

l3=[lx3;lr3;lf3];

%tau3

taux3=0; taur3=Mr3(fi); tauf3=r3u(y,it,fi);

tau3=[taux3;taur3;tauf3];

%n3

cn3=1/((T(y,it,fi)*r3u(y,it,fi))A2+r3u(y,it,fi)A2+(taur3-T(y,it,fi)*taux3)A2)A(1/2); nx3=cn3*T(y,it,fi)*r3u(y,it,fi); nr3=-cn3*r3u(y,it,fi); nf3=cn3*(taur3-T(y,it,fi)*taux3); n3=[nx3;nr3;nf3]; % n3_v=cross(l3,tau3);

% nx3=n3_v(1)*(n3_v(1)A2+n3_v(2)A2+n3_v(3)A2)A(1/2); % nr3=n3_v(2)*(n3_v(1)A2+n3_v(2)A2+n3_v(3)A2)A(1/2); % n^=n3_v(3)*(n3_v(1)A2+n3_v(2)A2+n3_v(3)A2)A(1/2); % n3=[nx3;nr3;nf3]; %Vbesk

Vbesk_vector = Vbesk*[cos(ugol); -sin(ugol)*cos((fi-1)*hf); sin(ugol)*sin((fi-1)*hf)]; %t3

txx3=Vbesk_vector(1)-nx3*dot(Vbesk_vector,n3);

trr3=Vbesk_vector(2)-nr3*dot(Vbesk_vector,n3);

tff3=Vbesk_vector(3)-nf3*dot(Vbesk_vector,n3);

ct=1/(txx3A2+trr3A2+tff3A2)A(1/2);

tx3=ct*txx3; tr3=ct*trr3; tf3=ct*tff3;

t3=[tx3;tr3;tf3];

h_ud(y,it,fi) = dot(Vbesk_vector,n3);

V3n(y,it,fi) = (2*kappa/h_ud(y,it,fi) + (kappa-1)*h_ud(y,it,fi))/(kappa+1); V3t(y,it,fi) = dot(Vbesk_vector,t3); V3(y,it,fi)=(V3n(y,it,fi)A2+V3t(y,it,fi)A2)A(1/2); p3(y,it,fi) = (1-kappa+2*h_ud(y,it,fi)A2)/(1+kappa); po3(y,it,fi)=h_ud(y,it,fi)/V3n(y,it,fi);

a3(y,it,fi) = kappa*(1+(kappa-1)*(MahA2)/2)-(kappa-1)*(V3(y,it,fi)A2)/2;

u3(y,it,fi)=V3n(y,it,fi)*nx3+V3t(y,it,fi)*tx3; v3(y,it,fi)=V3n(y,it,fi)*nr3+V3t(y,it,fi)*tr3; w3(y,it,fi)=V3n(y,it,fi)*nf3+V3t(y,it,fi)*tf3; V3vector=[u3(y,it,fi);v3(y,it,fi);w3(y,it,fi)];

%наклон характеристики %nh3

V3_x_tau3=cross(V3vector/(a3(y,it,fi)A(1/2)),tau3);

nh3 = (cross(tau3,V3_x_tau3) + V3_x_tau3*(dot(V3_x_tau3,V3_x_tau3)-dot(tau3,tau3))A( 1/2))/dot(V3_x_tau3,V3_x_tau3); rx3u(y,it,fi)=-nh3(1)/nh3(2); %lh3

lh3=[1,rx3u(y,it,fi),0];

if it==1; r1u(y,it,fi)=r3u(y,it,fi)-rx3u(y,it,fi)*hx; end;

if it>=2; r1u(y,it,fi)=r3u(y,it,fi)-(rx3u(y,it,fi)+rx1u(y,it-1,fi))*hx/2; end;

num(y,it,fi) = round((r1u(y,it,fi)-Gx0(y,fi))/hy(y,fi));

if num(y,it,fi)<=8; num_1(y,it,fi)=10-num(y,it,fi); end; if num(y,it,fi)>8; num_1(y,it,fi)=0; end; % if num(y,it,fi)==9; num_1(y,it,fi)=0; end;

% if num(y,it,fi)==10; num_1(y,it,fi)=0; end;

% if num(y,it,fi)==11; num_1(y,lt,fl)=0; end;

% end;

u_1(fi) = Input(fi,num_1(y,lt,fi)*4+1); v_1(fi) = Input(fi,num_1(y,lt,fi)*4+2); w_1(fi)=Input(fi,num_1(y,lt,fi)*4+3); p_1(fi) = Input(fi,num_1(y,lt,fi)*4+4); r_1(fi) = num(y,lt,fi)*hy(y,fi)+Gx0(y,fi);

u_2(fi) = Input(fi,num_1(y,lt,fi)*4+5); v_2(fi) = Input(fi,num_1(y,lt,fi)*4+6); w_2(fi)=Input(fi,num_1(y,lt,fi)*4+7); p_2(fi) = Input(fi,num_1(y,lt,fi)*4+8); r_2(fi) = (num(y,lt,fi)-

1)*hy(y,fi)+Gx0(y,fl);

u_3(fi) = Input(fi,num_1(y,lt,fi)*4+9); v_3(fi) = Input(fi,num_1(y,lt,fi)*4+10); w_3(fi)=Input(fi,num_1(y,lt,fi)*4+11); p_3(fl) = Input(fi,num_1(y,lt,fi)*4+12); r_3(fi) = (num(y,lt,fi)-

2)*hy(y,fi)+Gx0(y,fi);

u1_u(fi) = u_3(fí)+(u_2(fi)-u_3(fi))*(r1u(y,lt,fi)-r_3(fi))/hy(y,fí)+(u_3(fí)-2*u_2(fi)+u_1(fi))*(r1u(y,lt,fi)-r_3(fi))*(r1u(y,lt,fi)-r_2(fi))/(2*hy(y,fi)*hy(y,fi));

v1_u(fi) = v_3(fí)+(v_2(fi)-v_3(fi))*(r1u(y,lt,fi)-r_3(fi))/hy(y,fí)+(v_3(fí)-2*v_2(fí)+v_1(fí))*(r1u(y,lt,fí)-r_3(fi))*(r1u(y,lt,fi)-r_2(fí))/(2*hy(y,fi)*hy(y,fí));

w1_u(fi) = w_3(fl)+(w_2(fi)-w_3(fi))*(r1u(y,lt,fl)-r_3(fi))/hy(y,fi)+(w_3(fi)-2*w_2(fi)+w_1(fi))*(r1u(y,lt,fi)-r_3(fi))*(r1u(y,lt,fi)-r_2(fi))/(2*hy(y,fi)*hy(y,fi));

p1_u(fi) = p_3(fí)+(p_2(fi)-p_3(fi))*(r1u(y,lt,fi)-r_3(fi))/hy(y,fí)+(p_3(fí)-2*p_2(fí)+p_1(fí))*(r1u(y,lt,fí)-r_3(fi))*(r1u(y,lt,fi)-r_2(fí))/(2*hy(y,fi)*hy(y,fí));

V1(fi)=((u1_u(fi)л2+v1_u(fi)л2+w1_u(fi)л2))л(1/2);

a1(fi) = kappa*(1+(kappa-1)*(Mahл2)/2)-(kappa-1)*(V1(fl)л2)/2; po1(fi) = kappa*p1_u(fi)/a1(fi); V1vector=[u1_u(fi);v1_u(fi);w1_u(fi)] ;

%tau1

if num(y,lt,fi)>=11; num(y,lt,fi)=11;

df_taur = (-3*Mr(fi,num(y,lt,fi)-2)+4*Mr(fi,num(y,lt,fi)-1)-Mr(fi,num(y,lt,fi)))/(-2*hf); end;

if num(y,lt,fi)<11;

df_taur = (Mr(fi,num(y,lt,fi)+1)+Mr(fi,num(y,lt,fi)))/(2*hf); end;

taux1=0; taur1=df_taur; tauf1=r1u(y,it,fi); %taur1=Mr(fi,11);

tau1=[taux 1;taur1;tauf1] ;

%nh1

V1_x_tau1=cross(V1vector/(a1(fl)^(1/2)),tau1);

nh1 = (cross(tau1,V1_x_tau1) + V1_x_tau1*(dot(V1_x_tau1,V1_x_tau1)-dot(tau1,tau1))A( 1/2))/dot(V 1_x_tau1,V1_x_tau1); rx 1u(y,lt,fi)=-nh1 (1)/nh1(2); %lh1

lh1=[1;rx1u(y,it,fi);0]; %условие совместности

tau_x_n3=cross(tau3,nh3); n_x_l3=cross(nh3,lh3); tau_x_n1=cross(tau 1,nh1); n_x_l1=cross(nh1,lh1);

q13 = cross(tau3,V3vector)/(a3(y,lt,fi)^(1/2)); q23 = dot(V3vector,tau_x_n3)/a3(y,lt,fi); q33 = dot(V3vector,n_x_l3)/a3(y,lt,fi);

q43 = cross(V3vector,lh3)/(a3(y,it,fi)A(1/2));

q11 = cross(tau1,V1vector)/(a1(fi)A(1/2)); q21 = dot(V1vector,tau_x_n1)/a1(fi); q31 = dot(V1vector,n_x_l1)/a1(fi); q41 = cross(V1vector,lh1)/(a 1(fi)A(1/2));

E4=(q41(1)+q43(1))*(u3(y,it,fi)-u1_u(fi))/(2*hf)+(q41 (2)+q43 (2))*(v3(y,it,fi)-v1_u(fi))/(2*hf)+(q41 (3)+q43(3))* (w3(y,it,fi)-w1 _u(fi))/(2*hf)-(w1_u(fi)*q41(2)+w3(y,it,fi)*q43(2))/2+(v1_u(fi)*q41(3)+v3(y,it,fi)*q43(3))/2;

E1=(q11(1)+q13(1))*(u3(y,it,fi)-u1_u(fi))/(2*hx)+(q11(2)+q13(2))*(v3(y,it,fi)-v1_u(fi))/(2*hx)+(q11(3)+q13(3))*(w3(y,it,fi)-w1_u(fi))/(2*hx); E2=(q23/po3(y,it,fi)+q21/po1(fi))*(p3(y,it,fi)-p1_u(fi))/(2*hx); E3=(q33/po3(y,it,fi)+q31/po1(fi))*(p3(y,it,fi)-p1_u(fi))/(2*hf); if (fi>=6)&&(fi<10) E4=-(q41(1)+q43(1))*(u3(y,it,fi)-u1_u(fi))/(2*hf)-(q41(2)+q43(2))*(v3(y,it,fi)-v1_u(fi))/(2*hf)-(q41 (3)+q43(3 ))*(w3(y,it,fi)-w1 _u(fi))/(2*hf)-(w1_u(fi)*q41(2)+w3(y,it,fi)*q43(2))/2+(v1_u(fi)*q41(3)+v3(y,it,fi)*q43(3))/2; E3=-(q33/po3(y,it,fi)+q31/po1(fi))*(p3(y,it,fi)-p1_u(fi))/(2*hf); end;

Et(y,it,fi)=E1+E2+E3+E4; end; %fi end; %it

%fid2 = fopen(ss1, 'a'); for fi=1:N;

fprintf(fid2,'%5.4f ', u3(y,it,fi), v3(y,it,fi), w3(y,it,fi), p3(y,it,fi)); end;

fclose(fid2);

% конец ударная волна for fi=1:N

Fx 1(y,fi)=Fx0(y,fi)+hx*(T(y,Qu,fi)+T(y,1,fi))/2;

Gx1(y,fi)=(2*(x_start+(kt+1)*hx)/(5*cos((fi-1)*hf)A2+3*sin((fi-1)*hf)A2))A(1/2); %Gx1(y,fi)=(2*(x_start+(kt+1)*hx)/(3*cos((fi-1)*hf)A2+5*sin((fi-1)*hf)A2))A(1/2); hy 1(y,fi)=(Fx 1(y,fi)-Gx 1(y,fi))/10; end;

for fi=1:N

r1(y,fi)=Fx0(y,fi)-2*hy(y,fi); r2(y,fi)=Fx0(y,fi); r3(y,fi)=Fx1(y,fi)-hy1(y,fi); end;

%сплайны для r при x=x+hx. Мг3_ро1е-матрицы производных for fi=1:N %углов

for j=1:11 %точек на каждом угле

rr3_pole(fi^j)=Gx1(y,fi)+(j-1)*hy1(y,fi); end; end;

for j=1:11

for i=2:(N-1)

if i==2 Alpha_ch(i,j)=0; Beta_ch(i,j)=0; end; Alpha_ch(i+1,j)=-1/(Alpha_ch(i,j)+4);

Beta_ch(i+1,j)=Alpha_ch(i+1,j)*(Beta_ch(i,j)-3*(rr3_pole(i+1,j)-rr3_pole(i-1,j))/h^; end; end;

for j=1:11;

Mr3_pole(N,j)=0; for i=N:-1:2;

Mr3_pole(i-1,j)=Alpha_ch(i,j)*Mr3_pole(i,j)+Beta_ch(i,j); % производные для r end; end;

%конец сплайны для r при х=х+hx

%расчет точек в поле течения for j=9:-1:0;

if j==9; t=0; end; if j==8; t=0; end; if j<=1; t=24; end; if j==0; st=1; end;

j_2=j+2; j_1=j;

j_3=j+1;

if j==0; j_1=1; end; fid2 = fopen(ss1, 'a'); for l=1:3; fid4=fopen('splain.txt', 'w+'); for fi=1:N;

u01=Input(fi,t+1); v02=Input(fi,t+2); w03=Input(fi,t+3); p04=Input(fi,t+4);

u11=Input(fi,t+5); v12=Input(fi,t+6); w13=Input(fi,t+7); p14=Input(fi,t+8); %1

u31=Input(fi,t+9); v32=Input(fi,t+10); w33=Input(fi,t+11); p34=Input(fi,t+12);

u21=Input(fi,t+13); v22=Input(fi,t+14); w23=Input(fi,t+15); p24=Input(fi,t+16); %2

u41=Input(fi,t+17); v42=Input(fi,t+18); w43=Input(fi,t+19); p44=Input(fi,t+20);

u1=u21; v1=v22; w1=w23; p1=p24; u2=u11; v2=v12; w2=w13; p2=p14; if j==9; u1=u31;v1=v32;w1=w33;p1=p34; u2=u01;v2=v02;w2=w03;p2=p04; end; if j==1; u1=u41;v1=v42;w1=w43;p1=p44; u2=u31;v2=v32;w2=w33;p2=p34; end;

if st==1;

u1=u41; v1=v42; w1=w43; p1=p44; u2=u21; v2=v22; w2=w23; p2=p24; r1(y,fi)=Gx0(y,fi); r2(y,fi)=Gx0(y,fi)+hy(y,fi); r3(y,fi)=Gx 1(y,fi); end;

rx1(y,fi)=(r3(y,fi)-r1(y,fi))/hx; rx2(y,fi)=(r3(y,fi)-r2(y,fi))/hx;

if l==1;

u0(fi)=(u1 +u2)/2; v0(fi)=(v1+v2)/2; w0(fi)=(w1+w2)/2; p0(fi)=(p1+p2)/2; end;

a 1(y,fi)=kappa*(1+(kappa- 1)*(Mahл2)/2)-(kappa- 1)*(v1*v1+u1*u1+w1*w1)/2; a2(y,fi)=kappa*(l+(kappa- l)*(мahл2)/2)-(kappa- l)*(v2*v2+u2*u2+w2*w2)/2;

R1 = [(u31-u41)/(2*hy(y,fi)); (v32-v42)/(2*hy(y,fi)); (w33-w43)/(2*hy(y,fi)); (p34-p44)/(2*hy(y,fi))];

if j==1; R1=[(-3*u41+4*u21-u31)/(2*hy(y,fi));(-3*v42+4*v22-v32)/(2*hy(y,fi));(-

3*w43+4*w23-w33)/(2*hy(y,fi)); (-3*p44+4*p24-p34)/(2*hy(y,fi))]; end;

if j==9; R1=[(u11-u21)/(2*hy(y,fi));(v12-v22)/(2*hy(y,fi));(w13-w23)/(2*hy(y,fi));(p14-p24)/(2*hy(y,fi))]; end;

if st==1; R1=[(-3*u41+4*u21-u31)/(2*hy(y,fi));(-3*v42+4*v22-v32)/(2*hy(y,fi));(-3*w43+4*w23-w33)/(2*hy(y,fi)); (-3*p44+4*p24-p34)/(2*hy(y,fi))]; end;

R2 = [(u01-u31)/(2*hy(y,fi)); (v02-v32)/(2*hy(y,fi)); (w03-w33)/(2*hy(y,fi)); (p04-p34)/(2*hy(y,fi))];

if j==9; R2=[(-3*u01+4*u11-u31)/(-2*hy(y,fi));(-3*v02+4*v12-v32)/(-2*hy(y,fi)); (-3*w03+4*w13-w33)/(-2*hy(y,fi)); (-3*p04+4*p14-p34)/(-2*hy(y,fi))]; end;

if j==1; R2=[(u11-u21)/(2*hy(y,fi)); (v12-v22)/(2*hy(y,fi)); (w13-w23)/(2*hy(y,fi)); (p14-p24)/(2*hy(y,fi))]; end;

if st==1; R2=[(u31-u41)/(2*hy(y,fi)); (v32-v42)/(2*hy(y,fi)); (w33-w43)/(2*hy(y,fi)); (p34-p44)/(2*hy(y,fi))];end;

a0(l,fi)=kappa*(1+(kappa-1)*(Mahл2)/2)-(kappa-1)*(v0(fi)л2+u0(fi)л2+w0(fi)л2)/2;

%C

c11_1=-u1*w1/(a1(y,fi)*r1(y,fi));

c12_1=-v1*w1/(al(y,fi)*rl(y,fi));

c13_1=(1-w1Л2/a1(y,fi))/rl(y,fi);

c21_1=kappa*p1*w1/(a1(y,fi)*r1(y,fi));

c32_1=kappa*p1*w1/(al(y,fi)*rl(y,fi));

c43_1=kappa*p1*w1/(a1(y,fi)*r1(y,fi));

c44_1=1/r1(y,fi);

c11_2=-u2*w2/(a2(y,fi)*r2(y,fi));

c12_2=-v2*w2/(a2(y,fi)*r2(y,fi));

c 13_2=(1-w2Л2/a2(y,fi))/r2(y,fi);

c21_2=kappa*p2*w2/(a2(y,fi)*r2(y,fi));

c32_2=kappa*p2*w2/(a2(y,fi)*r2(y,fi));

c43_2=kappa*p2*w2/(a2(y,fi)*r2(y,fi));

c44_2=1/r2(y,fi);

c11_3=-u0(fi)*w0(fi)/(a0(l,fi)*r3(y,fi));

c12_3=-v0(fi)*w0(fi)/(a0(l,fi)*r3(y,fi));

c13_3=(1-w0(fi^2/a0(l,fi))/r3(y,fi);

c21_3=kappa*p0(fi)*w0(fi)/(a0(l,fi)*r3(y,fi));

c32_3=kappa*po(fi)*wo(fi)/(ao(l,fi)*r3(y,fi));

c43_3=kappa*po(fi)*wo(fi)/(ao(l,fi)*r3(y,fi));

c44_3=1/r3(y,fi);

%1

d11=2*(1-u^2/a 1(y,fi)+1 -u0(fi)л2/a0(l,fi)); d12=2*(-u1*v1/a1(y,fi)-u0(fi)*v0(fi)/ a0(l,fi)); d13=2*(-u1*w1/a1(y,fi)-u0(fi)*w0(fi)/a0(l,fi)); d14=0;

d15=hx*(-u1*v1/a1(y,fi)-rx1(y,fi)*(1-u1Л2/a1(y,fi))-u0(fi)*v0(fi)/a0(l,fi)-rx1(y,fi)*(1-u0(fi)л2/a0(l,fi))-c11_1*Mr(fi,j_1)-c11_3*Mr3_pole(fi,j_3));

d16=hx*(1-v1Л2/a1(y,fi)+rx1(y,fi)*u1*v1/a1(y,fi)+ 1-

v0(fi)л2/a0(l,fi)+rx1(y,fi)*u0(fi)*v0(fi)/a0(l,fi)-c12_1*Mr(fi,j_1)-c12_3*Mr3_pole(fi,j_3));

d17=hx*(-v1*w1/a1(y,fi)+rx1(y,fi)*u1*w1/a1(y,fi)-v0(fi)*w0(fi)/a0(l,fi)+rx1(y,fi)*u0(fi)*w0(fi)/a0(l,fi)-c13_1*Mr(fi,j_1)-c13_3*Mr3_pole(fi,j_3)); d18=0; %2

d21=2*kappa*(u1*p1/a1(y,fi)+u0(fi)*p0(fi)/a0(l,fi));

d22=0;

d23=0;

d24=2*(1+1);

d25=hx*kappa*((v1-rx1(y,fi)*u1)*p1/a1(y,fi)+(v0(fi)-rx1(y,fi)*u0(fi))*p0(fi)/a0(l,fi)-c21_1*Mr(fi,j_1)/kappa-c21_3*Mr3_pole(fi,j_3)/kappa); d26=0; d27=0;

d28=hx*(-2*rx1(y,fi)); %3 d31=0;

d32=2*kappa*(u1*p1/a1(y,fi)+u0(fi)*p0(fi)/a0(l,fi));

d33=0;

d34=0;

d35=0;

d36=hx*kappa*((v1-rx1(y,fi)*u1)*p1/a1(y,fi)+(v0(fi)-rx1(y,fi)*u0(fi))*p0(fi)/a0(l,fi)-c32_1*Mr(fi,j_1)/kappa-c32_3*Mr3_pole(fi,j_3)/kappa); d37=0;

d38=hx*(1+1); %4 d41=0; d42=0;

d43=2*kappa*(u1*p1/a1(y,fi)+u0(fi)*p0(fi)/a0(l,fi));

d44=0;

d45=0;

d46=0;

d47=hx*kappa*((v1-rx1(y,fi)*u1)*p1/a1(y,fi)+(v0(fi)-rx1(y,fi)*u0(fi))*p0(fi)/a0(l,fi)-c43_1*Mr(fi,j_1)/kappa-c43_3*Mr3_pole(fi,j_3)/kappa);

d48=hx*(-c44_1*Mr(fi,j_1)-c44_3*Mr3_pole(fi,j_3)); %5

d51=2*(1-u2л2/a2(y,fi)+1-u0(fi)л2/a0(l,fi)); d52=2*(-u2*v2/a2(y,fi)-u0(fi)*v0(fi)/a0(l,fi));

d53=2*(-u2*w2/a2(y,fi)-u0(fi)*w0(fi)/a0(l,fi)); d54=0;

d55=hx*(-u2*v2/a2(y,fi)-rx2(y,fi)*(1-u2Л2/a2(y,fi))-u0(fi)*v0(fi)/a0(l,fi)-rx2(y,fi)*(1-u0(fi)л2/a0(l,fi))-c11_2*Mr(fi,j_2)-c11_3*Mr3_pole(fi,j_3));

d56=hx*(1-v2Л2/a2(y,fi)+rx2(y,fi)*u2*v2/a2(y,fi)+1-v0(fi)л2/a0(l,fi)+rx2(y,fi)*u0(fi)*v0(fi)/a0(l,fi)-c12_2*Mr(fi,j_2)-c12_3*Mr3_pole(fi,j_3));

d57=hx*(-v2*w2/a2(y,fi)+rx2(y,fi)*u2*w2/a2(y,fi)-v0(fi)*w0(fi)/a0(l,fi)+rx2(y,fi)*u0(fi)*w0(fi)/a0(l,fi)-c13_2*Mr(fi,j_2)-c13_3*Mr3_pole(fi,j_3)); d58=0; %6

d61=2*kappa*(u2*p2/a2(y,fi)+u0(fi)*p0(fi)/a0(l,fi));

d62=0;

d63=0;

d64=2*(1+1);

d65=hx*kappa*((v2-rx2(y,fi)*u2)*p2/a2(y,fi)+(v0(fi)-rx2(y,fi)*u0(fi))*p0(fi)/a0(l,fi)-c21_2*Mr(fi,j_2)-c21_3*Mr3_pole(fi,j_3)); d66=0; d67=0;

d68=hx*(-2*rx2(y,fi)); %7 d71=0;

d72=2*kappa*(u2*p2/a2(y,fi)+u0(fi)*p0(fi)/a0(l,fi));

d73=0;

d74=0;

d75=0;

d76=hx*kappa*((v2-rx2(y,fi)*u2)*p2/a2(y,fi)+(v0(fi)-rx2(y,fi)*u0(fi))*p0(fi)/a0(l,fi)-c32_2*Mr(fi,j_2)-c32_3 *Mr3_pole(fi,j _3)); d77=0;

d78=hx*(1+1); %8 d81=0; d82=0;

d83=2*kappa*(p2*u2/a2(y,fi)+p0(fi)*u0(fi)/a0(l,fi));

d84=0;

d85=0;

d86=0;

d87=hx*kappa*((v2-rx2(y,fi)*u2)*p2/a2(y,fi)+(v0(fi)-rx2(y,fi)*u0(fi))*p0(fi)/a0(l,fi)-c43_2*Mr(fi,j_2)-c43_3*Mr3_pole(fi,j_3));

d88=hx*(-c44_2*Mr(fi,j_2)-c44_3*Mr3_pole(fi,j_3));

D = [dll d12 d13 d14 d15 d16 d17 d18; d21 d22 d23 d24 d25 d26 d27 d28; d31 d32 d33 d34 d35 d36 d37 d38; d41 d42 d43 d44 d45 d46 d47 d48; d51 d52 d53 d54 d55 d56 d57 d58; d61 d62 d63 d64 d65 d66 d67 d68; d71 d72 d73 d74 d75 d76 d77 d78; d81 d82 d83 d84 d85 d86 d87 d88];

%Fij if st==0 if l==1

F11=2*hx*(-v1/r1(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_1+c11_3)*(M1 (fi,tt+1)+M1_3_0(fi,tt+1))+(c12_1+c12_3) *(M1(fi,tt+2)+M1_3_0(fi,tt+2))+( c13_1+c13_3)*(M2(fi,tt+3)+M2_3_0(fi,tt+3)))-d15*R1(1,1)-d16*R1(2,1)-d17*R1(3,1)+d11*u1+d12*v1+d13*w1;

F21=-hx*(c21_1+c21_3)*(M1(fi,tt+1)+M1_3_0(fi,tt+1))-d25*R1(1,1)-d28*R1(4,1)+d21*u1+d24*p1;

F31=2*hx*(wlл2*kappa*p1/(a1(y,fi)*r1(y,fi))+w0(fi)л2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_1+c32_3)*(M1(fi,tt+2)+M1_3_0(fi,tt+2))-d36*R1(2,1)-d38*R1(4,1)+d32*v1;

F41=2*hx*(-v1*w1*kappa*p1/(a1(y,fi)*r1(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_1+c43_3)*(M2(fi,tt+3)+M2_3_0(fi,tt+3))+(1/r1(y,fi)+1/r3(y,fi))*(M1(fi,tt+4)+M1_3_0(fi,tt +4)))-d47*R1(3,1 )+d43*w1;

F51=2*hx*(-v2/r2(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_2+c11_3)*(Ml(fi,tt+9)+M1_3_0(fi,tt+1))+(c12_2+c12_3)*(M1(fi,tt+10)+M1_3_0(fi,tt+2))+ (c13_2+c13_3)*(M2(fi,tt+11)+M2_3_0(fi,tt+3)))-d55*R2(1,1)-d56*R2(2,1)-d57*R2(3,1)+d51*u2+d52*v2+d53*w2;

F61=-hx*(c21_2+c21_3)*(M1(fi,tt+9)+M1_3_0(fi,tt+1))-d65*R2(1,1)-d68*R2(4,1 )+d61*u2+d64*p2;

F71=2*hx*(w2л2*kappa*p2/(a2(y,fi)*r2(y,fi))+w0(fi)л2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_2+c32_3)*(M1(fi,tt+10)+M1_3_0(fi,tt+2))-d78*R2(4,1)-d76*R2(2,1)+d72*v2;

F81=2*hx*(-v2*w2*kappa*p2/(a2(y,fi)*r2(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_2+c43_3)*(M2(fi,tt+11)+M2_3_0(fi,tt+3))+(1/r2(y,fi)+1/r3(y,fi))*(M1(fi,tt+12)+M1_3_0(fi, tt+4)))-d87*R2(3,1)+d83*w2; end;

if l>1

F11=2*hx*(-v1/r1(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_1+c11_3)*(M1 (fi,tt+1)+MM1_3(fi,1))+(c12_1+c12_3) *(M1(fi,tt+2)+MM1_3(fi,2))+(c13_1 +c13_3)*(M2(fi,tt+3)+MM2_3(fi,3)))-d15*R1(1,1)-d16*R1(2,1)-d17*R1(3,1)+d11*u1+d12*v1+d13*w1;

F21=-hx*(c21_1+c21_3)*(M1(fi,tt+1)+MM1_3(fi,1))-d25*R1(1,1)-d28*R1(4,1)+d21*u1+d24*p1;

F31=2*hx*(wlл2*kappa*p1/(a1(y,fi)*r1(y,fi))+w0(fi)л2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_1+c32_3)*(M1(fi,tt+2)+MM1_3(fi,2))-d36*R1(2,l)-d38*Rl(4,1)+d32*v1;

F41=2*hx*(-v1*w1*kappa*p1/(a1(y,fi)*r1(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_1+c43_3)*(M2(fi,tt+3)+MM2_3(fi,3))+(1/r1(y,fi)+1/r3(y,fi))*(M1(fi,tt+4)+MM1_3(fi,4)))-d47*R1(3,1)+d43*w1;

F51=2*hx*(-v2/r2(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_2+c11_3)*(Ml(fi,tt+9)+MM1_3(fi,1))+(c12_2+c12_3)*(M1(fi,tt+10)+MM1_3(fi,2))+(c13_ 2+c13_3)*(M2(fi,tt+11)+MM2_3(fi,3)))-d55*R2(1,1)-d56*R2(2,l)-d57*R2(3,1)+d51*u2+d52*v2+d53*w2;

F61=-hx*(c21_2+c21_3)*(M1(fi,tt+9)+MM1_3(fi,1))-d65*R2(1,1)-d68*R2(4,1 )+d61*u2+d64*p2;

F71=2*hx*(w2A2*kappa*p2/(a2(y,fi)*r2(y,fi))+w0(fi)A2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_2+c32_3)*(M1(fi,tt+10)+MM1_3(fi,2))-d78*R2(4,1)-d76*R2(2,1)+d72*v2;

F81=2*hx*(-v2*w2*kappa*p2/(a2(y,fi)*r2(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_2+c43_3)*(M2(fi,tt+11)+MM2_3(fi,3))+(1/r2(y,fi)+1/r3(y,fi))*(M1(fi,tt+12)+MM1_3(fi,4)) )-d87*R2(3,1)+d83*w2; end; end;

%расчет точек на обтекаемом теле if st==1; tt=36; if l==1

F11=2*hx*(-v1/r1(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_1+c11_3)*(M1(fi,tt+1)+M1_3_0(fi,tt+1))+(c12_1+c12_3)*(M1(fi,tt+2)+M1_3_0(fi,tt+2))+( c13_1+c13_3)*(M2(fi,tt+3)+M2_3_0(fi,tt+3)))-d15*R1(1,1)-d16*R1(2,1)-d17*R1(3,1)+d11*u1+d12*v1+d13*w1;

F21=-hx*(c21_1+c21_3)*(M1(fi,tt+1)+M1_3_0(fi,tt+1))-d25*R1(1,1)-d28*R1(4,1)+d21*u1+d24*p1;

F31=2*hx*(w1A2*kappa*p1/(a1(y,fi)*r1(y,fi))+w0(fi)A2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_1+c32_3)*(M1(fi,tt+2)+M1_3_0(fi,tt+2))-d36*R1(2,1)-d38*R1(4,1)+d32*v1;

F41=2*hx*(-v1*w1*kappa*p1/(a1(y,fi)*r1(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_1+c43_3)*(M2(fi,tt+3)+M2_3_0(fi,tt+3))+(1/r1(y,fi)+1/r3(y,fi))*(M1(fi,tt+4)+M1_3_0(fi,tt +4)))-d47*R1(3,1 )+d43*w1;

F51=2*hx*(-v2/r2(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_2+c11_3)*(M1(fi,tt+5)+M1_3_0(fi,tt+1))+(c12_2+c12_3)*(M1(fi,tt+6)+M1_3_0(fi,tt+2))+( c13_2+c13_3)*(M2(fi,tt+7)+M2_3_0(fi,tt+3)))-d55*R2(1,1)-d56*R2(2,1)-d57*R2(3,1)+d51*u2+d52*v2+d53*w2;

F61=-hx*(c21_2+c21_3)*(M1(fi,tt+5)+M1_3_0(fi,tt+1))-d65*R2(1,1)-d68*R2(4,1 )+d61*u2+d64*p2;

F71=2*hx*(w2A2*kappa*p2/(a2(y,fi)*r2(y,fi))+w0(fi)A2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_2+c32_3)*(M1(fi,tt+6)+M1_3_0(fi,tt+2))-d78*R2(4,1)-d76*R2(2,1)+d72*v2;

F81=2*hx*(-v2*w2*kappa*p2/(a2(y,fi)*r2(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_2+c43_3)*(M2(fi,tt+7)+M2_3_0(fi,tt+3))+(1/r2(y,fi)+1/r3(y,fi))*(M1(fi,tt+8)+M1_3_0(fi,tt +4)))-d87*R2(3,1)+d83*w2; end;

if l>1

F11=2*hx*(-v1/r1(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_1+c11_3)*(M1 (fi,tt+1)+MM1_3(fi,1))+(c12_1+c12_3) *(M1(fi,tt+2)+MM1_3(fi,2))+(c13_1 +c13_3)*(M2(fi,tt+3)+MM2_3(fi,3)))-d15*R1(1,1)-d16*R1(2,1)-d17*R1(3,1)+d11*u1+d12*v1+d13*w1;

F21=-hx*(c21_1+c21_3)*(M1(fi,tt+1)+MM1_3(fi,1))-d25*R1(1,1)-d28*R1(4,1)+d21*u1+d24*p1;

F31=2*hx*(w1A2*kappa*p1/(a1(y,fi)*r1(y,fi))+w0(fi)A2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_1+c32_3)*(M1(fi,tt+2)+MM1_3(fi,2))-d36*R1(2,1)-d38*R1(4,1)+d32*v1;

F41=2*hx*(-v1*w1*kappa*p1/(a1(y,fi)*r1(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_1+c43_3)*(M2(fi,tt+3)+MM2_3(fi,3))+(1/r1(y,fi)+1/r3(y,fi))*(M1(fi,tt+4)+MM1_3(fi,4)))-d47*R1(3,1)+d43*w1;

F51=2*hx*(-v2/r2(y,fi)-v0(fi)/r3(y,fi))-hx*((c11_2+c11_3)*(M1(fi,tt+5)+MM1_3(fi,1))+(c12_2+c12_3)*(M1(fi,tt+6)+MM1_3(fi,2))+(c13_2 +c13_3)*(M2(fi,tt+7)+MM2_3(fi,3)))-d55*R2(1,1)-d56*R2(2,1)-d57*R2(3,1)+d51*u2+d52*v2+d53*w2;

F61=-hx*(c21_2+c21_3)*(M1(fi,tt+5)+MM1_3(fi,1))-d65*R2(1,1)-d68*R2(4,1 )+d61*u2+d64*p2;

F71=2*hx*(w2A2*kappa*p2/(a2(y,fi)*r2(y,fi))+w0(fi)A2*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-hx*(c32_2+c32_3)*(M1(fi,tt+6)+MM1_3(fi,2))-d78*R2(4,1)-d76*R2(2,1)+d72*v2;

F81=2*hx*(-v2*w2*kappa*p2/(a2(y,fi)*r2(y,fi))-v0(fi)*w0(fi)*kappa*p0(fi)/(a0(l,fi)*r3(y,fi)))-

hx*((c43_2+c43_3)*(M2(fi,tt+7)+MM2_3(fi,3))+(1/r2(y,fi)+1/r3(y,fi))*(M1(fi,tt+8)+MM1_3(fi,4)))-d87*R2(3,1)+d83*w2; end;

d11_st(y,fi)=Gx 1(y, fi)/(2*(x_start+(kt+ 1)*hx)); d12_st(y,fi)=-1;

d13_st(y,fi)=Gx1(y,fi)A2*sin(2*(fi-1)*hf)/(2*(x_start+(kt+1)*hx)); %(1/r)*dG/dfi F11=0;

D =[d11_st(y,fi) d12_st(y,fi) d13_st(y,fi) 0 0 0 0 0; d21 d22 d23 d24 d25 d26 d27 d28; d31 d32 d33 d34 d35 d36 d37 d38; d41 d42 d43 d44 d45 d46 d47 d48; d51 d52 d53 d54 d55 d56 d57 d58; d61 d62 d63 d64 d65 d66 d67 d68; d71 d72 d73 d74 d75 d76 d77 d78; d81 d82 d83 d84 d85 d86 d87 d88];

end;

%конец расчета стенки

F =[F11;F21;F31;F41;F51;F61;F71;F81];

X = D\F;

u0(fi)=X(1,1);

v0(fi)=X(2,1);

w0(fi)=X(3,1);

p0(fi)=X(4,1);

fprintf(fid4,'%5.4f u0(fi), v0(fi), w0(fi), p0(fi)); end; %fi - расчет одной точки сразу по всем углам fclose(fid4);

%создание SP_3 - матрица для расчета производных по углам, 2 итерация

fid5=fopen('splain.txt', 'r'); SP=fscanf(fid5,'%g'); for s=1:N

for e=1:4

SP_3(s,e)=SP(e+(s-1)*4);

end; end; fclose(fid5);

%сплайны для четной функции,ММ1_3-матрица производных в точке 3 for e=1:4 for s=2:(N-1)

if s==2 Alpha_ch_3(s,e)=0; Beta_ch_3(s,e)=0; end;

Alpha_ch_3(s+1,e)=-1/(Alpha_ch_3(s,e)+4);

Beta_ch_3(s+1,e)=Alpha_ch_3(s+1,e)*(Beta_ch_3(s,e)-3*(SP_3(s+1,e)-SP_3(s-1,e))/hf); end; end;

for e=1:4

MM1_3(N,e)=0; for s=N:-1:2

MM1_3(s-1,e)=Alpha_ch_3(s,e)*MM1_3(s,e)+Beta_ch_3(s,e); end; end;

%конец сплайны для четной функции

%сплайны для нечетной функции, М2_3-матрица производных в точке 3 for e=3:3

for s=2:(N-1)

if s==2 Alpha_nech_3(s,e)=-0.5; Beta_nech_3(s,e)=3*SP_3(s,e)/(2*hf); end;

Alpha_nech_3(s+1,e)=-1/(Alpha_nech_3(s,e)+4);

Beta_nech_3(s+1,e)=Alpha_nech_3(s+1,e)*(Beta_nech_3(s,e)-3*(SP_3(s+1,e)-SP_3(s-

1,e))/hf);

end; end;

for e=3:3

MM2_3(N,e)=(-Beta_nech_3 (N,e)-3*SP_3(N- 1,e)/hf)/(2+Alpha_nech_3(N,e)); for s=N:-1:2

MM2_3(s-1,e)=Alpha_nech_3(s,e)*MM2_3(s,e)+Beta_nech_3(s,e); end; end;

%конец сплайны для нечетной функции

end; %l 3 итерации

%записываем результаты в файл fid6=fopen('splain.txt', 'r'); A_1=fscanf(fid6,'%g'); for s=1:N for e=1:4

A_2(s,e)=A_1(e+(s-1)*4); end;