автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Математическое моделирование пространственных течений газа в соплах

МОСКОВСКИЙ ОРДЕНА ЛЕНИНА, ОРДЕНА. ОКТЯБРЬСКОЙ РЕВОЛЮЦИИ И ОРДЕНА ТРУДОВОГО КРАСНОГО ЗНАМЕНИ ГОСУДАРСТВЕННЫЙ

УНИВЕРСИТЕТ ИМЕНИ М.В. ЛОМОНОСОВ А ФАКУЛЬТЕТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И КИБЕРНЕТИКИ

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

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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННЫХ ТЕЧЕНИЙ ГАЗА

В СОПЛАХ

Специальность 05.13.16 - "Применение вычислительной техники, математического моделирования и математических методов в научных исследованиях"

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Москва - 1998

Диссертация выполнена на факультете вычислительной математики и кибернетики Московского государственного университета им.М В .Ломоносова.

Научный руководитель:

доктор физико-математических наук, профессор Росляков Геннадий Степанович.

Официальные оппоненты:

доктор технических наук, чл.-корр. РАН, профессор Пирумов Ульян Гайкович,

доктор физико-математических наук, профессор Тишкин Владимир Федорович.

Ведущая организация:

НПО Энергомаш им.академика В.ПТлушко.

Защита состоится

оиш2

га К.053.05.

1999г. в

1Н

2.0

часов на

заседании диссертационного совета К.053.05.87 в Московском государственном университете им.М.В.Ломоносова, на факультете вычислительной математики и кибернетики (119899, Москва, Воробьевы горы, МГУ, 2 учебный корпус, факультет ВМиК, ауд. 685).

С диссертацией можно ознакомиться в библиотеке МГУ во 2 учебном корпусе.

Автореферат разослан'

1999г.

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

В.М.Говоров

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

Научная новизна. Разработан модифицированный полностью неявный метод (МПНМ) решения системы алгебраических уравнений, возникающих при симметричной аппроксимации на 19-точечном пространственном шаблоне краевой задачи для линейного дифференциального уравнения в частных производных эллиптического типа. Впервые МПНМ применен для расчета пространственных безвихревых течений газа в соплах. Такие течения описываются квазилинейным уравнением потенциала скорости смешанного типа, поэтому перенос численных методов, разработанных для линейных эллиптических уравнений, на уравнение потенциала не является тривиальным. К основным достоинствам нового алгоритма относятся высокая скорость сходимости и приемлемая точность. Для расчета смешанных течений, содержащих большие сверхзвуковые области, предложен алгоритм, представляющий собой комбинирование МПНМ и метода приближенной факторизации (метода ПФ). С помощью реализованных алгоритмов показана возможность построения сопел с прямоугольным выходом, по тяговым характеристикам близких к осесимметричным соплам.

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

Апробация работы. Основные результаты работы докладывались на VIII школе-семинаре "Современные проблемы аэрогидродинамики" (Севастополь, сентябрь 1996), на XXI научных чтениях по хосмонавтике (Москва, январь 1997), на Международной конференции студентов и аспирантов по фундаментальным наукам "Ломоносов" (Москва, апрель 1997).

Структура и объем работы. Диссертация состоит га введения, двух глав (первая глава содержит два параграфа, вторая - три), заключения и списка литературы из 78 названий. Общий объем, включая 32 рисунка, составляет 117 страниц.

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

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

Исследование пространственных течений газа в соплах представляет большой научный и практический интерес. В ряде случаев довольно точные предсказания можно получить, используя модель безвихревого течения. Такое течение описывается одним уравнением потенциала скорости, при этом все газодинамические функции выражаются через функцию потенциала. В последние годы рядом авторов были разработаны весьма эффективные численные методы решения уравнения потенциала, получившие название релаксационных методов. К ним относятся метод последовательной верхней релаксации (метод ПВРЛ) и различные варианты метода приближенной факторизации (метод ПФ). Сейчас известно несколько разновидностей метода ПФ, которые отличаются друг от друга видом одномерных операторов. Методы ПВРЛ и ПФ стали применяться американскими исследователями около 20 лет назад для численного решения внешних задач аэродинамики. Первые работы по использованию релаксационных методов для исследования задач внутренней газодинамики были выполнены в нашей стране: Шифриным Э.Г., Шулановым М.А.1; Ивановым М.Я., Корецким В.В.2; Веретенцевым В.А., Меченовой В.А., Росляковым Г.С.3. Кроме методов ПФ для решения уравнения потенциала скорости, хорошо зарекомендовали себя так называемые полностью неявные методы. Пионерская работа была опубликована Стоуном4, который применил, разработанный им метод и названный полностью неявным (ПНМ) для решения системы алгебраических уравнений, получающихся при дискретизации линейного дифференциального уравнения в частных производных эллиптического типа с двумя независимыми переменными на 5-точечном шаблоне. Шнейдером и Зеданом5 этот метод был распространен на 9-точечную схему и получил название модифицированного полностью неявного метода (МПНМ). За счет модификации удалось достичь ускорения сходимости в 2-3 раза. Меченова В.А. и Росляков Г.С.6 перенесли МПНМ на полное уравнение потенциала скорости смешанного типа для исследования двумерных течений газа в соплах. МПНМ был распространен и на пространственный случай7. В этой публикации авторы приводят свои расчетные формулы для 7-точечного шаблона и результаты проведенных по этой схеме расчетов трехмерного стационарного уравнения теплопроводности. Наиболее общим в этом классе методов является модифицированный полностью неявный метод решения системы алгебраических уравнений, получающихся при дискретизации линейного дифференциального уравнения в частных производных эллиптического типа с тремя независимыми переменными на 19-точечном шаблоне. На его основе можно легко получить расчетные формулы для 7-точечной схемы в трехмерном случае и 9-

1 Шнфрин Э.Г., Шуланов М.А. Решение прямой задача дпх плоского сопла Лаваля релаксационным численным методой по схеме Мурмена-Коула. //Учен. зап. ЦАГИ -1981. -Т. 12, №3.

2 Иванов МЯ., Корсцкнй ВВ. Расчет течений в двумерных и пространственных соплах методом праблнжеаной факторизации. //ЖВМ к МФ. -1985. -Т.25, №9.

1 Веретенцев В.А., Меченова В. А., Росляков Г.С. Численное исследование течений газа в каналах и соплах на основе уравнений для потенциала. //Нестационарные течения газов с ударными волнами. -Ленинград: Издательство ЛФТИ, 1990.

4 Stone H.L. Iterative Solution of Implicit Approximation of Multidimensional Partial Differential Equations. //Siam Journal of Numerical Analysis.-1968,-V.5, Sept

5 Schneider G., Zedan M. A Modified Strongly Implicit Procedure for the Numerical Solution of Field Problems. //Numerical Heat Transfer. -1981. - V.4, Jan.

6 Веретенцев В.А., Меченова В.А., Росляков Г.С. Укагработа,

7 Зедан М, ЦГнейдер Г.Э. Модифицированный полностью неявный метод расчета трехмерных температурных полей. //Аэрокосмическа* техника. -1983. -Т,1, №11.

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

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

Симметричная аппроксимация эллиптического уравнения

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

+ Д*М?\./-и+Д!м?'|+и-и + А/ЛЙ-иЛ +Д* мРц* +Д™.*й-1./+и + ^

М^и-и^ + 1 +

где

Д.у.» > > Д../> > ■ Д./.» > Д,у.» > Д../* • Д.;.* > Д./.» » Д'.» > Д./> > Д./.*»Д./.» > Д,/,* > А> Д¡Уд»Д' ДГл* > АЗ.»»Qt.it

- известные величины, <р - неизвестные. В матричной форме систему можно переписать в следующем виде: Аф = д. Уравнения в этой системе расположены в порядке возрастания значений i,j,k таким образом, что к меняется от 0 до К, у меняется от О до J, i меняется от 0 до Т, при этом цикл по г вложен в цикл по ], который, в свою очередь, вложен в цикл по к. Вид матрицы А показан на рис.1. Матрицу А можно представить в виде произведения двух треугольных матриц: нижней V и верхней 1Г, то есть А = 1ХГ. Каждая го них является квадратной (Г++ ЩК +1) х (/ ++- 1)(Х +1) -матрицей. Однако в общем случае V и £/' имеют (/+1)(./ + 2)+1 ненулевых диагоналей каждая. Ненулевые элементы V расположены между Аы -диагональю и Ар -диагональю, ненулевые элемента II' — между А' -диагональю и А1" -диагональю (рис.1). Прямое разложение матрицы А в произведение нижней и верхней треугольных матриц, не ^учитывающее разреженную структуру матрицы А, не является.эффективным по затратам машинного времени. Построим на основе матрицы А модифицированную матрицу А' в некотором смысле близкую матрице А. Затем прибавив к обеим частям системы (1) вектор-столбец А'ф„ получим уравнение А'ф = А'ф-(Аф-д), которое будет служить основой для организации итерационного процесса. Модифицированная матрица Л должна быть пред ставима в виде произведения сильно разреженных нижней Ь и верхней {/ треугольных матриц. Для этого потребуем, чтобы матрицы I и [/ имели по 10 ненулевых диагоналей каждая и все их ненулевые элементы располагались на местах ненулевых элементов исходной матрицы А. Структура матриц £ и и показана на рис.2, а и б соответственно. Кроме того, потребуем, чтобы ненулевые элементы

матрицы А и расположенные на тех же самых местах элементы матрицы А' были соответственно равны. Из этих условий получаем следующие уравнения для определения элементов матриц Ь и 17:

а!.и = 4>

^¡.М + аи.к°и1-1.к~1 = >

^и* +сшиил =

еи.1 + + с1,мР1.м-1 + ^./.Л+им = .

Лм = А/*»

'¡.М + А.у.Л-и,* ! + + Л;.* А ц-и ~ '

+суЛЛ +

ти.Л.ц ^АлЛШ-! +К1*Р,+и-и=А,'м>

тш°ш +е,.,Л/ии = -Са-

+ 1 ~ 4,/Л >

+ЛлЛии ^ + «..М^-У-и + = 4'.*, = 4^-

Из уравнений (2) можно получить явные рекуррентные формулы для определения

величин ам1г, /и>) 1я,м> о,м,

Г'.М> *>.м> и,./л> . поскольку к моменту вычисления элементов матриц

Ь,и с индексами г,],к элементы матрицы С/ с индексами /,у - -1, (- \,],к-1;

/ + + У — 1,_/' — 1,А; i + \,j-\,k■, / — 1,уже

известны. Модифицированная матрица А'= Ы], структура которой показана на рие.З, кроме 19 ненулевых диагоналей исходной матрицы содержит 24 дополнительные ненулевые диагонали, обозначенные через я) Выражения для

дополнительных элементов матрицы А' имеют следующий вид:

~ ^у/нм-ь ^(./Л =г1.М°|.)+и-1> и-Ь п1.],к~е1.1,кг1.1*и-1>

= +¡1.1^1-1.,*, (3)

Очевидно, что дополнительные элементы все сразу нулю не равны и, следовательно, матрица Л' не может быть идентичной матрице Л.

Построенная таким образом матрица А' легко разложима в произведение нижней и верхней сильно разреженных треугольных матриц, элементы которых вычисляются по рекуррентным формулам. Поэтому итерационный метод решения системы (1), имеющий вид: А'5"'1 = Щ", будет эффективным по затратам машинного времени на одну итерацию (здесь 3"*' -ф" - вектор разности, Я" Аф" -вектор невязки, Д - итерационный параметр). Однако по скорости сходимости этот итерационный процесс может оказаться неэффективным. Поэтому построим на основе матрицы А не одну модифицированную матрицу А', а семейство новых матриц {Л1}, конкретная представительница которого зависела бы от набора параметров. Меняя эти параметры (то есть выбирая разные матрицы из семейства {Л'}), можно было бы оказывать влияние на скорость сходимости итерационного процесса А'б"'1 = Щ".

Рассмотрим семейство матриц {А'}; построенных на основе матрицы А системы (I) следующим образом. К правой части каждого уравнения (2) прибавим линейную комбинацию величин Тогда с учетом (3) получим следующие

соотношения (для экономии места запишем лишь первое и последнее уравнения):

+ +

+ 1.М-1 +Лм°м./и) + % +К»гм.1-(41)

Г=1

+ (4„ )

Коэффициенты меняются при переходе от одного узла сетки к

другому, однако ограничимся случаем, когда они не зависят от номера узла (»,/,£)•

Каждое уравнение (4) линейно относительно неизвестных

Поэтому из уравнений (4,)-(4!1) можно получить явные аналитические выражения для этих величин. Полученные значения а,.1Л>^>>'>fi.ijfSi.jj.yh.ijoh.ij: зат€м используются в (410

19 ) Для

вычисления Таким образом, каждая

матрица из {А'} легко разложима в произведение нижней £ и верхней и сильно разреженных треугольных матриц, структура которых изображена на рис.2, а и б соответственно. Конкретная представительница семейства матриц {Л'} определяется заданием 456-компонентного вектора (л^,...,^,...,^,...,^) (здесь 456=19x24). На практике при выборе матрицы А' для конструирования итерационного метода решения системы линейных алгебраических уравнений (1) диапазон допустимых значений ограничен и зависит от конкретной задачи. Подбор комбинации 456 параметров, которая обеспечила бы наиболее быструю скорость сходимости итерационного процесса, является чрезвычайно сложной задачей. Кроме того, в общем случае этот подбор нужно производить на каждой итерации. Упростим задачу, построив подсемейство семейства матриц {А'}, элементы'которого зависят от одного параметра айв некотором смысле являются близкими матрице А.

На рис.4 изображен разностный шаблон, определяемый модифицированной матрицей А'. Кружочками обозначены узлы, соответствующие исходному сеточному уравнению, крестиками - узлы, вовлекаемые дополнительными ненулевыми элементами матрицы А'. Чтобы уменьшить влияние этих новых узлов, применим так называемый прием частичного исключения дополнительных членов. Для этого значение функции <р в узле, соответствующем дополнительному элементу матрицы А', выразим с помощью разложения в ряд Тейлора через значения <р в узлах, принадлежащих 7-точечному шаблону. Здесь используется 7-точечный шаблон для того, чтобы окончательные формулы для 19-точечной схемы включали формулы для 7-точечной схемы как частный случай. Пренебрегая членами второго и более высоких порядков, имеем:

<Рх = и ~2<Рш + 9м.ц, + 9,.н* = А»

= <РиХ!Л-\ й ~2<Ри.к + 2<Рм,1* + Фш т. = А> 9г =Р(-г/+ик-1 »-Зру* + 9и,лк +9,.м-1 =А»

9* ~ 9,- 1.,чи-1 * -2 Фи* + <Р,л.1м + 9,.,** + 9,.1.1,-1 = Д.. <Ръ ~ Рм.ц\.к-\ * -2<Ри* + 9м.ц. + <Р,,1*.к + <Рш-1 = А. 9с + и +<Ри.кл =А.

97 = 9,-¡./+2л и -Щ.,.* + 9,-+ + 9i.ij.-i = А. 9г = 9,.1>гк 1 +Й.Й-1 = А»

% =9м.у,2.к-1 ~~3<Ри.к +9м.1* +2<Рил* +9>,1,к-\ = А.

9т = 91.1-и, м -9ш + 29им.к = До- Рп = и + 2<Р..и> + = А1,

¥>12 = 9,-2.11, М ~9i.lt + 29,Л.,.к = А>. Из = и -9',.,л + = Аз.

Рн = .»V « + 2рм м + к = Д4, <ри = « -р,,, + 29,.,<1* = А5,

Ри, =Р,-и-1М1 а ~ъ9ил +9,-111. +29и-» *9>.1м\ = Аб.

Р)7 = 9,.]-Хк* ~ ~29ш + 29,.1-х.к + Рим, = Ат.

Pin = 9>I-I.hm¡ « +2 + + ÍY/.ttl = A»

?>20 = W-U-l.i+l * -2<P,.M +<Р,-Ц.к +<Pl./ lk I = D20,

P21 и + +P¡.y-U + = Al.

^22 = 5¡>i+2J-u+I ю + + <p,,,M 1 = Ав.

<Pn = Pi-2.ijc» a ~2<?,.1* 4 2pM /Jt + <pl/M, = Da,

P24 = <P,-wm\ ~2PíJ* + Pi-ш + + ftjjw = A»-

Рассмотрим следующую модификацию схемы (1):

u-i +4!"jtPfj,u-i + 4"t9Vu-i> +

+ Kii'Pi.i^.t +4"» PMJ-U +a,"Í.I.9>Í-W +ÁiT.,*<Pt.jt +Ai'.j.k9M./jt + 4"*<эи./<и +

+ +aIJ¿9>U-u+I +АЪл<Рим\ + KjíVmjmi +

2-1

где а - параметр. В этом разностном уравнении члены Dt (А = 1,...Д4) зависят от значений функции tp в узлах 7-точечного шаблона. Следовательно, прежде чем строить L U -разложение модифицированной матрицы Л, необходимо соответствующим образом перегруппировать члены. После такой перегруппировки уравнения принимают следующий вид:

tfjjtVi./->.*-i+ + ¿W<P,j.t-i + + i.t-1 + 4"t?Vi.í-u +

+ + + ¿ti.t'Pi.i-i.M + ¿fi.M-w* + ¿w<Pi.ji.»« + 4 j.kPw../.»-«+

24

¿Ul

где надчеркнугые коэффициенты вычисляются по формулам:

—ь '

¿¡.М = A?JJt -a^icfj,,,

V A=i6 Í-'ii )

Лш -ai +*,«)+2 +

\J=3.19.23 ¿-12,U if 7,1« '

ÜÍM = 4^ +«f Efe+2^; +3 +3<j)+

УЛ-1.4 ,1=7.16

Д=10.13 i-19,22 /

з'м=4v -4 +2^)4-2 +5Х.Д

Я-1.5.21 «

= 4?,.* -« —/ 24

И.И

Теперь можно записать модификацию уравнений (2): +а1.,Л°и 1.4 1 = Д./.* »

ил ш у и '

—Ь 9

Д=1

+ = Д.»>

/■ нал

-а + + ,

V -1=1« ¿=19 )

Км =

К1* + = -

-«С •

V ^3.19.23 ¿=П,И У

,1=10.13 .1=19.22 У

ти*"ш +си*уш-1 М-м".^./.*-! +««Ач» +КисРм.н.к: = Ллм = ■

Л=9.18

^А/Аии =

Щ./Л/.1 + Н.,= (5)

—/ м тихиих +8илн'и1.» +/илум.м = Ли' = Лбл

= л*

Уравнения (5) имеют форму (4), то есть каждое из них линейно относительно неизвестных а11МЬ>]М сш,(111Л,,,«¡.^Кц> 1ил - Поэтому из первых девяти уравнений (5) с учетом соотношений (3) можно получить явные формулы для определения этих величин:

а|.М = Д.Л» > = Д.» ~ а1.1.к°1.1-\*-1'

Ж^к,

■ft.lt — Д!1*

кл

+ ~ fi.ltЛ-Т .J-l.lt ~~

К

+ЛуЛ2<>/-и-и +2/м.ьи +и|-и-и)+

Когда вычислены, по формулам (3)

вычисляют дополнительные элементы. Затем полученные значения используются для

вычисления Явные формулы для

определения этих величин таковы: . ,

т1.!.к = Д'¡.к ~а1.]*М'и-1Л-1 ~Си.кии*-1 1./Л ~

~еШ*и+и-1 ~ ~Si.ixPi.lM-К,^.!-^ ~ lj.lt +

+а

.#«1.4 ^7.1«

/»-10.13 0-19.22

П1.].к ~ 81.,^.;-и \jxPM.i-1Л

(1=1.5.21 ¿=11.13 Д-9,18

= ~еиЛ1+1,/Ы -hj.kPvA.yk )!ти* •

Ри,к = {^л -еШии+и-1 ~Ь.1.кГ1-х.Цк ~

-а +

/1=3 ¿.7

= К"* ~¿¡.и^ы.!* 1 ~е1.иУ1.1ЛЛ-\\1п1,1Л >

'¡.Л* = ~ >

"ш = - Яи^и-и - - •

(7)

Применив формулы (б), (3), (7) ко всем узлам расчетной сетки, определим все элементы нижней £ и верхней Г/ треугольных матриц, произведение которых есть модифицированная матрица А'.

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

лЫ _ АЬ* _ дЫ _ дЪн _ ажш _ 1г* __

~ — Л./,* ~ "1.1.к ~ ~ К].к ~

-А** = А" - А* -А* - А*' - А* -О

откуда следует

«>.м = ьш = Л.,.к = гш = *и.к = ~ О

и

«А — - .я-6 — - _10 _ „12 _

71 I.¡.к - - - м - - 7СШ -

_ „13 _ 15 _ „16__19 _ „22 _ „24 _ п

Построенная таким образом матрица А' легко разложима в произведение нижней и верхней сильно разреженных треугольных матриц. Поэтому итерационный метод решения системы линейных алгебраических уравнений, имеющий вид: А16*" = , будет эффективным по затратам машинного времени на одну итерацию (здесь 5"*х -ф"*х -ф" - вектор разности, Я"=д-Аф" - вектор невязки, р -итерационный параметр). Важно подобрать такую комбинацию параметров а и /}, при которой этот итерационный процесс был бы эффективным и по скорости сходимости. На практике диапазон допустимых значений а и /7 ограничен и зависит от конкретной задачи.

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

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

д( гдЛ д( др\ д ( 1

дх) <> дг) Зе\ г де) Отнесем линейные размеры к некоторой характерной длине /., скорость и плотность — к их критическим значениям а. и р., давление - к р.а1, тогда из интеграла Бернулли и из условия сохранения энтропии имеем зависимость плотности от потенциала:

и зависимость давления от плотности (и, следовательно, от потенциала):

Г

Здесь у - показатель адиабаты.

На граничных поверхностях физической области течения определенным образом ставятся дополнительные условия. Полагаем, что на входе вектор скорости ортогонален начальному сечению х = а. Вследствие этого величина р на входе постоянна, для определенности: <р = 0. На поверхности сопла г=/(х,е) ставим

условие непротекания: <рг - На плоскостях симметрии е=0 и е- с'

ставится условие симметрии: <рс ~ 0.

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

£ = £(*), 7 =-,в—в(в) осуществим преобразование исходной области

физического пространства в прямоугольный параллелепипед (ограниченный плоскостями | = |(я), € = »7 = 0, 7 = 1, в = 0(0), в = в(е*), здесь х~Ъ -координата выходного сечения сопла) в пространстве, где будем проводить численный расчет на равномерной сетке. Использовалось такое преобразование координат, при котором расчетная сетка имела бы сгущение в физических координатах в областях сильного изменения газодинамических параметров. Уравнение потенциала скорости в новых переменных можно записать в следующей форме:

д4{ з J м з ) дек. з ;

Здесь

и = <р& - 9Л У = у+<Р,,

Выражение для плотности в новых переменных £ 7, в имеет вид

Г Г Г) V" (9)

Н^г-^гК^^т' (10)

где г - .

Граничные условия в переменных £,т],в будут выглядеть следующим образом:

а) на входе (4 = 4(a)): ?> = 0, (ll,)

б) на поверхности сопла (/7 = 1): F = О, (ll2)

в) на плоскостях в - 0(0) и 0 = в(е*): W = 0. (ll3) Требуется найти функцию p{4,jj,ff), определенную на множестве [4(a)4(b),0<ij<. 1,0(0) <&<9(е')} и удовлетворяющую уравнению (8), (9), (10) и граничным условиям (11).

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

Для уравнения (8) консервативную схему второго порядка аппроксимации можно записать в виде:

F , -F .

+— G , -ё , +— Atj\ и*-* I Ав

U.t+- i.M—

Я ,-Я. , =0, (12)

X рги ~ ртУ ~ (ЯГ

где г = --, Сг = -—,Я = -— - потоки. Величины потоков вычисляются через

3 3 3

производные от потенциала с помощью соотношений (9). Производные, входящие в (9), аппроксимируются со вторым порядком с помощью центральных разностей.

Уравнение потенциала имеет эллиптический тип при дозвуковом характере течения, гиперболический - при сверхзвуковых скоростях и параболический - на звуковой поверхности. Схема с центральными разностями (12) приводит к неустойчивым алгоритмам в областях сверхзвукового течения, поскольку передает возмущения вверх по потоку. Для обеспечения устойчивого счета в сверхзвуковой области вводилась искусственная сжимаемость в виде поправок к плотности. Модификация плотности вида

где - разностный оператор первой производной назад по направлению V -коэффициент искусственной сжимаемости и М — число Маха, вносит в схему ошибку первого порядка1. Введение искусственной сжимаемости эквивалентно использованию односторонних разностей против потока при аппроксимации уравнения в сверхзвуковой зоне течения. Можно вводить модифицированную плотность не по одному направлению 4, а по всем трем направлениям 4,4,0- В некоторых расчетах по продольной координате 4 использовалось другое выражение для р:

Р 1 =Р X +

1 Мэрией Э., Саут Дж., Хафез М, Применение методов качественной сжимаемости да» численного решения полного уравнения потенциала в трансзвуковом диапазоне скоростей. //Ракетна« техника я космонавтика. -1979. -Т.17, №8.

+ктах-ш,1-АГ

(Д£)2

I

\-М'\ )-<5Т

8\8:р ,

* * 14-l.il

Эта модификация плотности, которая вносит в схему ошибку второго порядка, была предложена Веретенцевым В. АЛ

Поверхность сопла задавалась уравнением (в декартовых координатах)

чЧ*) г

У Л ^

лМ

где уЛх) ~ образующая сопла в плоскости 2 - 0, г+(х) - образующая сопла в плоскости у = 0, а>(х) - показатель "суперэллиптичности" сопла. Если у(х) з г(х) и а>(х) э 2, то сопло осесимметричное.

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

0,625; а = -4,0; 6 = 1,0. На рис.6

(рис.5): 9а =~'А = = 3,125; Д, = 1,0; Л,

представлены распределения числа Маха на контуре (кривая 1) и оси (кривая 2) сопла (сетка 65x43x5), крестиками обозначены результаты названной работы. На рис.7 изображено распределение значений расхода по поперечным сечениям сопла. Кривая 1

соответствует сетке 33x22x22 ^0 5 е <, ^, кривая 2 - 42x22x17 < е 2 , кривая 3 -

50x30x10

< с <, , кривая 4 - 65x43x5 < £ 5 . В качестве второго примера было исследовано осесимметричное сопло3, имеющее следующие параметры (рис.5):

50 = —= —;Д,, = 2,2; Д, = 0,7;/?, = 0,625;а = -3,0;6 = 0,9. Расчеты проводились на 4 12

тех же самых равномерных сетках, что и в названной работе: 20x6x6 (сетка 1), 30x9x9 (сетка 2), 40x12x12 (сетка 3). Для уменьшения нормы невязки в 1000 раз (относительно нулевой итерации) на сетке 1 потребовалось 62 итерации (в названной работе -приблизительно 140), на сетке 2-108 итераций (в названной работе - приблизительно 200), на сетке 3 - 194 итерации (в названной работе - приблизительно 500). Проводились дополнительные расчеты на более мелких (по осевому и радиальному направлениям) сетках: 50x20x5, 65x43x5. На рис.8 показано поведение расхода. Кривая

1 соответствует сетке 40x12x12

кривая 2 - 50x20x5 0 ^ е <

=4

40 )

кривая 3

■ 65x43x5 0 < е < — .С увеличением числа узлов сетки разброс значений расхода по 40]

1 Веретенцев В.А. Численное исследование течений газа в соплах сложной геометрии. //Кандидатская диссертация. - М.: Моск. университет, 1991.

2 Меченова В.А., Росляков Г.С. Применение полностью неявных методов дня решения прямой задачи теории сопла. //Численные методы в математической физике. - М.: Издательство МГУ, 1996.

3 Корецкий В.В., Любимо» ДА Модифицированный метод приближенной факторизация для расчета потенциальных пространственных течений в каналах. //ЖВМ я М4>. -1990. - Т.30. №10.

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

После исследования осесиммегричиых задач проводился расчет течения газа в сопле, имеющем одну плоскость симметрии1. В полуплоскости у>0 плоскости симметрии 2 = 0 образующая сопла у+(х) имеет следующие параметры (рис.5):

¿>0 =—=— = 3,0;^ =1,0;^ = 1,0; а=-5,0; £ = 0,4, а в полуплоскости ><0 6 12

параметры образующей у (х) точно такие же, кроме одного: <90 = —. В плоскости хг

4

образующая (х)еу (г). Уравнение поверхности сопла имеет вид:

= 1. Расширяющаяся часть такого сопла геометрически является

и±(*); V *♦(*),!

осесимметричной, однако течение в этой области имеет пространственный характер. Неравномерность потока распространяется и на цилиндрическую часть сопла. Графики боковой силы (рис.9) и момента (рис.10) демонстрируют хорошее качественное совпадение с результатами названной работы (они обозначены крестиками), полученные методом установления.

Некоторое представление о структуре течения внутри сопла с прямоугольным выходом дает следующий иллюстративный пример. В плоскости >> = 0 образующая г+(х) имеет параметры (рис.5):

.?„= — — ;Л„ =2,0;^ =2,0;Л2=2,0;а = -4,0;£ = 2,0. В плоскости г = 0 6 12

образующая >,+(х) = 0,5г+(г). Функция "суперэллиптичности": <а(х)^2 при »¿0 и <»(г) = 2хг +2 при х > 0. На рис.11 изображены линии уровня давления в выходном сечении такого сопла. Имеет место выраженное сгущение этих линий в окрестности е = агс%2, и можно ожидать, что с ростом х газодинамические параметры либо их производные будут терпеть разрыв на поверхности е = агс<& 2.

Также было исследовано две группы новых, используемых на практике сопел. Для таких сопел в плоскости у = 0 на участке а<,х<. О образующая 7«. (ж) имеет

параметры (рис.5): 9а =— =ЗД25;Л =1ДЛ, = 0,625;а = -3,85, а при г>0

■•■>.- 4

образующая г+(х) задается набором значений полученным в результате

решения задачи профилирования, по которому для функции строится

интерполирующий кубический сплайн. В плоскости г = 0 образующая (х) з г+(х). Всего было исследовано четыре сопла, которые отличались друг от друга наборами ■ -{хх>гхУм 11)111 ВИД°М функции о(х): сопла разных групп отличались между собой наборами {^.г^}^, внутри одной группы сопла отличались видом функции а>(х) -для первого сопла каждой группы (сопла 1.1 и 2.1) функция а>(х) бралась в виде <»(*)= 2 при х¿0 и <э(х)з2х2 + 2 при х > О, сопла 1.2 и 2.2 выбирались осесимметричными (о>(х) = 2). Представляет ■ интерес сравнение распределений импульса в сверхзвуковой части для осесимметричного сопла и для сопла, выходное сечение которого близко к квадратному. Рис.12 иллюстрирует зависимость импульса от

' Дворецкий В.М. К жсследокшию пространственных смешанных течений в соплах с несимметричным входом. //Известия АН СССР. МЖГ. -1975. - №2.

площади поперечного сечения сопла (кривая 1 соответствует соплу 1.1, кривая 2 -соплу 1.2, кривая 3 - импульс, вычисленный по одномерной теории). Как уже отмечалось, сверхзвуковые контуры осесимметричных сопел 1.2 и 2.2 построены в результате решения задачи профилирования и, тем самым, они являются в некотором смысле "хорошими" соплами. Проведенные исследования показали, что можно сконструировать сопла с прямоугольным сечением (экономично заполняющие выход) по тяговым характеристикам близкие к осесиммегричным.

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

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

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

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

/ \»М / \»00 сопла имеет вид: ——— + —-— =1, где образующие

1,5-0,5ат|-+—\о<:с< 10 . ч „„ ТТ

10/ а ю(х) з 20. Число Маха в

2,10 < х < 20

Л О) = *,(*) =

1 Иванов М.Я., Корецюай В .В. Расчет течений в двумерных я пространственных соплах методом приближенной факторизации. //ЖВМ иМФ. -1985. - Т.25, №9.

2 Мечснова В.А., Росляков Г. С. Применение полностью неявных методов дая решения пр»мой задачи теории сопла. //Численные методы в математической физнхе. - М.: Издательство МГУ, 1996.

минимальном сечении полагалось равным единице. На рис.13 изображены линии уровня давления в плоскости у = О (в плоскости г = 0 картинка получилась аналогичной). Расчеты методом ПФ проводились также и для новых сопел с квадратным выходом, используемых на практике. Сравнивались результаты, полученные МПНМ и методом ПФ.

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

В заключении сформулированы основные результаты диссертации. К ним относятся:

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

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

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

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

Рис,4

Рис.6

Рис.7 Рис. В

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

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ ОПУБЛИКОВАНЫ В РАБОТАХ:

1. Росляков Г.С., Федоренко В.В. Исследование пространственного течения газа в сопле полностью неявным методом. //Методы математического моделирования. -М.: Издательство факультета ВМиК МГУ, 1998. - С.76 - 86.

2. Росляков Г.С., Федоренко В.В. Исследование течения газа в пространственном сопле полностью неявным методом. //Известия АН СССР. МЖГ. - 1997. - №4. -Материалы УШ школы-семинара "Современные проблемы аэрогидродинамики" (Севастополь, 4-13 сентября 1996). - С.179 -191.

3. Федоренко В.В. Исследование пространственных сопел. //Материалы XXI научных чтений по космонавтике. - М., 1997. - С.15 —16.

4. Федоренко В.В. Исследование пространственных сопел. //Материалы Международной конференции студентов и аспирантов по фундаментальным наукам "Ломоносов". -М.: Издательство МГУ, 1998. - С.204.

5. Федоренко В.В. Применение релаксационных методов к расчету пространственных течений газа в соплах. //Депонированная рукопись. Деп. в ВИНИТИ. -26.10.98, №3075-В98.

МОСКОВСКИЙ ОРДЕНА ЛЕНИНА, ОРДЕНА ОКТЯБРЬСКОЙ РЕВОЛЮЦИИ И ОРДЕНА ТРУДОВОГО КРАСНОГО ЗНАМЕНИ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМЕНИ М.В.ЛОМОНОСОВА ФАКУЛЬТЕТ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И

КИБЕРНЕТИКИ

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

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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННЫХ

ТЕЧЕНИЙ ГАЗА В СОПЛАХ

Специальность 05.13.16 - "Применение вычислительной техники, математического моделирования и математических методов в

научных исследованиях"

Диссертация на соискание ученой степени кандидата физико-

математических наук

Научный руководитель - доктор физико-математических наук, профессор Г.С.Росляков

Москва - 1998

ОГЛАВЛЕНИЕ

Введение..............................................................................3

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

1.1. Модифицированный полностью неявный метод, использующий 19-точечную аппроксимацию............................................12

1.2. Процедура расчета для 7-точечной схемы.............................34

Глава 2. Исследование смешанных пространственных течений газа в

соплах различных конфигураций.................................45

2.1. Математическая постановка задачи.....................................45

2.2. Исследование пространственных течений модифицированным полностью неявным методом.............................................54

2.2.1. Методические расчеты...........................................64

2.2.2. Течения в соплах с прямоугольным выходом...............83

2.3. Исследование пространственных течений методом приближенной факторизации............................................92

Заключение.......................................................................108

Библиография....................................................................109

ВВЕДЕНИЕ

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

Поиск аналитических решений полной системы дифференциальных уравнений в частных производных, описывающих течение газа внутри сопла на основе модели Эйлера, затруднен в силу нелинейности уравнений и их смешанного типа. Тем не менее, аналитические методы сыграли и продолжают играть большую роль при качественном исследовании течений газа в соплах. Были развиты асимптотические методы исследования течений в трансзвуковой области сопла: метод малых возмущений и разложение решения в ряд в окрестности прямолинейной звуковой линии. Изложение этих результатов и библиографические данные содержатся в [49, 50]. В монографии [56] приведено решение методом малых возмущений пространственных нестационарных уравнений для трансзвуковой области сопла. В работе [46] метод малых возмущений применен для

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

Область применимости аналитических методов ограничена. По этой причине, начиная с 50-х годов, для исследования газодинамических процессов в соплах стали применяться численные методы. Первые расчеты течений газа в соплах были связаны с задачей профилирования сверхзвуковой части осесимметричных и плоских сопел с угловой точкой и плоской поверхностью перехода через скорость звука. Для расчетов использовался метод характеристик [28, 30, 51, 53]. В последующие годы были разработаны вычислительные алгоритмы метода характеристик для более сложных типов течений: течений с химическими реакциями, двухфазных и слоистых течений [26, 27, 44, 50, 59]. Большой фактический материал по соплам с угловой точкой содержится в [38, 50].

Были разработаны методы профилирования сопел, реализующих на выходе заданный неравномерный [31, 34], в частности неизоэнтропический поток. Метод характеристик был применен для профилирования сверхзвуковой части пространственного сопла [6].

Сопло является основным элементом любого реактивного двигателя, эффективность которого во многом зависит от степени совершенства сопла. В реактивных летательных аппаратах цена единицы удельного импульса сопла достаточно высока, поэтому большой интерес представляет выбор оптимального контура, обеспечивающего максимум тяги при тех или иных ограничениях. Многими авторами разрабатывались вариационные принципы построения оптимальных сопел [33, 58, 68]. В работе [7] решена задача профилирования оптимальной сверхзвуковой части

пространственного сопла. Принципиально новый подход к профилированию, позволяющий получить контур сопла во всей до-, транс- и сверхзвуковой области, базируется на решении обратной задачи [31,45,47].

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

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

созданию новых схем, получивших название сеточно-характеристических методов [37]. Методика выделения разрывов для расчета пространственных течений применялась в работах [14, 15, 29, 66, 74].

Дополнительные трудности возникают, когда в потоке имеет место интерференция разрывов. По этой причине рядом авторов были разработаны методы, получившие название методов сквозного счета. Они имеют однородную структуру во всей расчетной области, независимо от того, имеются в потоке разрывы или нет. При этом в разностном решении поверхности разрывов получаются в виде узких областей больших градиентов газодинамических параметров. При расчете разрывных решений обычно используются консервативные формы записи газодинамических уравнений, то есть уравнения в виде законов сохранения, и консервативные (дивергентные) разностные схемы. Для нестационарных уравнений оригинальную консервативную разностную схему первого порядка точности разработал С.К.Годунов [12, 13]. В дальнейшем эта схема была перенесена и на стационарные уравнения [22, 23]. Метод сквозного счета для исследования пространственных течений применялся в [4, 52].

Долгое время для расчета трехмерных течений использовался так называемый метод установления [5, 16, 17, 18, 24, 25, 35, 60]. Идея метода основана на том, что нестационарные уравнения газовой динамики являются гиперболическими независимо от того до-, транс-или сверхзвуковым является поток. Для нестационарных уравнений решается краевая задача с граничными условиями стационарной задачи. Искомое стационарное решение получается как предел, к которому стремится нестационарный процесс с ростом г. Несмотря на то, что размерность задачи увеличивается на единицу, такой прием во многих случаях оправдан.

Во многих прикладных задачах не удается описать течение газа, используя лишь модель идеального газа. В реальном течении имеет место одновременное протекание физико-химических процессов, различных по своей природе. Вопросы методологии расчета и результаты исследований таких течений рассмотрены в [47 - 50, 57, 59].

Систематическое изложение имеющихся в настоящее время сведений о течениях газа в соплах дано в монографии [49]. Эта монография представляет собой существенно переработанный, в отдельных аспектах сокращенный вариант первой монографии авторов [50], дополненный изложением новых результатов. Обе книги в определенной степени дополняют друг друга. Кроме перечисленных работ обширный фактический материал содержится в справочнике [1] под редакцией В.П.Глушко.

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

зарекомендовавшие себя методы [2, 3, 42]. К релаксационным методам решения уравнения потенциала скорости относятся метод последовательной верхней релаксации (метод ПВРЛ), подробно описанный в [65], и различные варианты метода приближенной факторизации (метод ПФ). Сейчас известно несколько разновидностей методов факторизации: ПФ1, ПФ2, ПФЗ, которые отличаются друг от друга видом одномерных операторов. Методы факторизации стали применяться американскими исследователями около 20 лет назад для численного решения внешних задач аэродинамики [64, 65, 69 - 73, 75, 76].

Первые работы по использованию релаксационных методов для исследования задач внутренней газодинамики были выполнены в нашей стране [21, 67]. Авторы публикации [21] впервые применили метод ПФ для расчета потенциальных осесимметричных течений газа в соплах. При этом для обеспечения устойчивости в сверхзвуковой зоне использовались односторонние разности. Дальнейшее развитие методы ПФ получили в работах [8, 9, 10]. В этих публикациях представлены результаты расчетов двумерных потенциальных течений газа в соплах (для обеспечения устойчивости в сверхзвуковой зоне использовалась модифицированная плотность). Важные результаты были получены авторами работы [11], которые впервые распространили метод ПФ на случай двумерного непотенциального течения газа. В результате проведенных исследований оказалось, что методы ПФ на порядок более эффективны методов установления. Однако в случае пространственных течений методы ПФ имеют медленную скорость сходимости по сравнению с осесимметричным случаем (по числу итераций). Таким образом, проблема создания быстрых методов расчета пространственных течений газа в соплах, позволяющих конструкторам работать в режиме вычислительного эксперимента, остается актуальной.

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

Первая глава диссертации посвящена полностью неявным методам решения системы алгебраических уравнений, возникающих при аппроксимации краевых задач для эллиптических уравнений. Пионерская работа была опубликована Стоуном [78], который применил разработанный им метод и названный полностью неявным (ПНМ) для решения системы алгебраических уравнений, получающихся при дискретизации линейного дифференциального уравнения в частных производных эллиптического типа с двумя независимыми переменными. При этом использовался 5-точечный шаблон. В дальнейшем [77] метод Стоуна был распространен на 9-точечную схему и получил название модифицированного полностью неявного метода (МПНМ). При использовании 5-точечной аппроксимации МПНМ неэквивалентен ПНМ Стоуна. За счет модификации удалось достичь ускорения сходимости в 2-3 раза. В работах [9, 39, 40] МПНМ был перенесен на полное уравнение потенциала скорости смешанного типа для исследования двумерных течений газа в соплах. Затем модифицированный полностью неявный метод был распространен на пространственный случай - в публикации [20] авторы приводят свои расчетные формулы для 7-точечного шаблона и результаты проведенных по этой схеме расчетов трехмерного стационарного уравнения теплопроводности.

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

точечном шаблоне. На его основе можно легко получить расчетные формулы для 7-точечной схемы в трехмерном случае и 9-точечной схемы в двумерном случае. Однако до сих пор ни в отечественной, ни в зарубежной доступной литературе полное описание МПНМ, содержащее расчетные формулы для случая использования 19-точечной аппроксимации, не было опубликовано. Диссертация устраняет этот пробел.

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

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