К расчету стационарных сверхзвуковых течений с учетом реальных свойств газа Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ ЦАГ И Т о м IX 197 8

М 4

УДК 533.6.011

К РАСЧЕТУ СТАЦИОНАРНЫХ СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ С УЧЕТОМ РЕАЛЬНЫХ СВОЙСТВ ГАЗА

А. А. Голубинский, А. П. Косых

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

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

Методом С. К. Годунова получены решения большого числа задач внешней и внутренней аэродинамики в предположении о том, что рассматриваемый газ совершенный (х = cp/cv = const, см. библиографию в [4]). Учет реальных свойств газов, находящихся в термодинамическом равновесии, может существенно изменить картину течения. Поэтому имеет смысл приспособить широко применяемый метод к расчету равновесных течений газа. Для расчета таких течений, кроме численного интегрирования системы дифференциальных уравнений газодинамики, необходимо рассчитывать термодинамические функции газа. Этот расчет можно проводить, решая систему нелинейных алгебраических уравнений химического равновесия (такой подход требует больших затрат машинного времени) или используя табличное задание термодинамических функций. В данной работе реальные свойства газов учитывались с помощью таблиц и аппроксимаций, предложенных в [5], и созданной на их основе С. В. Пироговой стандартной программы.

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

*)

¡0

=— ударные волны;

— границы веера волн разрежения; - — тангенциальные разрывы

Фиг. 1

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

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

— ударная волна — тангенциальный разрыв — ударная волна (фиг. 1, а);

— волна разрежения — тангенциальный разрыв — волна разрежения (фиг. 1, б);

— ударная волна —тангенциальный разрыв — волна разрежения (фиг. I, в, г).

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

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

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

1. Расчет течения за двумя косыми скачками уплотнения. Запишем систему уравнений газодинамических законов сохранения для потоков / и 2 (см. фиг. 1, а) в общей форме, обозначив через р, р, V, Л давление, плотность, модуль вектора скорости потока, удельную энтальпию соответственно. Для каждого из потоков двумерного течения используем уравнение неразрывности, уравнение сохранения количества движения при прохождении через косой скачок уплотнения и уравнение сохранения энергии:

дексом „п“ отмечены проекции скоростей потоков на нормали к скачкам, параметры с индексом „¿=1“ относятся к первому потоку, а с „¿ = 2“ — ко второму.

Вводя в рассмотрение угол {3 между косым скачком и направлением набегающего потока, эту систему уравнений можно привести к следующему виду:

ных уравнений (1) и (2).

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

Параметров. Рск\у Рек 1» Ріт ^скії Рск2і Рск2> ^ск2-

Воспользуемся табличным заданием термодинамических функций: энтальпии Л = Л(/7, р) и эффективного показателя адиабаты -х==х(/7, р). Тогда систему (1), (2) можно дополнить соотношениями для значений энтальпии за скачками

Кроме того, из условий совпадения направлений потоков и знЭ' чений давления за скачками следуют соотношения:

ТОг = (рГСКп)г;

^СЕлЧ-^Ск) > ^ 2,

где индексом „ск“ отмечены значения параметров за скачками, ин-

(1)

(2)

Лек 2 = ¿2 + \ VI 5ІП2 ¡32 [і - (-ІІ-)2' .

В дальнейших рассуждениях будем исходить иэ преобразован-

Лск 1—к(Рск1> Рек і)> Лск 2 — Ь(Рск 2> Рсь г)1

(3)

ІБ(0І — 8і) = (62 + 5г)> Рек 1 =Рск 2 =Р,

(4)

где 0! и 02 '-г- углы между набегаюьцими потоками 1 и 2 и направлением неподвижной оси х (угол/считается положительным, если он отсчитываетсЯтдрввдвчасовой7^стрелки); §! и &2 — положительные приращения углов поворота^потоков относительно направлений набегания при переходе через скачок.

Из уравнений сохранения массы определим и й2:

Используя уравнения (1) — (5), составим замкнутую систему семи нелинейных алгебраических уравнений относительно семи неизвестных р, рСК1, р„ ЛСК1, рск2, ра, /гск2:

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

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

которое приближенно выполняется и для слабых скачков уплотнения (С —тангенс угла наклона вектора скорости к оси х, М —

Рек/

Р/

/= 1, 2.

(5)

^ск і—Л {р, рек і);

р = р, + ь УІ8ІП*Р2(і—-£-)

\ Рек 2 1

(6)

Лек 2 — Л (р> Рек г)>

(7)

л±Ц-£^м*-1 Лр=о,

р V2

(8)

число М). Переходя от соотношения (8) к конечным разностям, для 1 и 2 потоков получим соотношения:

2

«Ж-Л) = 0. М?-1, ¿=1, 2, (9)

4 ч ^ ^ Р 1*1

где Са = 0!, С2 = взГ

Разрешив эту линейную систему относительно р, получим некоторое „нулевое“ приближение для давления за скачками

р = р(°) = (а1р1 + а2р2 + Сх - С ¡¡Ж«! + а2). (10)

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

ря - рй„ рг, (11)

рщ - А йя, л2>>. <12>

В первоначальном варианте, реализованном на ЭВМ, для решения каждой из подсистем в (6) был применен метод простой итерации по параметру р. Все выкладки и рассуждения проведем для первого потока, для второго потока они будут аналогичны. Цепочка вычислений для первого потока была взята следующая:

р1+1±\р(0) ( Аск1 = Л(/>(0), Рс°к 0,

Р(0) - Рс°к\= *+ Л Р1 - ~(0) __ . , _1_ р<°> -Рг (! , _рЛ (13)

^Л+;( ) [Аск1 —«!+ 2 Р1 \ р£°\] '

Здесь значение плотности р(°>1 получено из адиабаты Гюгонио для совершенного газа. Сравнивая полученные значения энтальпии

Лек 1 и Лск1 и если они с заданной точностью в (е>0, наперед заданная малая величина) не совпадают, то итерации продолжаем по формуле

рО>■ = Л^<0>- >"■)-*■ р.— Л ’р.. (14)

^ т <'“-*> )

Зная р*1^, находим А(ск 1 = А(р{0), рйО, Ъ[к1 = Ь{ри Р1) + р( ]~р1 + и т- д-

Итерации продолжаются до момента совпадения ЛЙ1 и ЛсК\ с с заданной точностью (/=1, 2, 3, . . .). Аналогично проводятся итерации для второй подсистемы из системы (6). Как показали пробные расчеты задачи о столкновении двух симметричных потоков воздуха (высота Н — 5 км) с числами М1 = М2 = 5 и углами 61 = = — 62, которые изменялись от 0 до 30°, для значений | 61,215° вышеописанный итерационный процесс сходился. Начиная же с углов 101,2 [ ]> 15° алгоритм оказался быстро расходящимся. После тщательного анализа получаемых приближений и поведения функции И{р, р) и ее частной производной Ар(р, р) мы убедились, что правая

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

Наиболее подходящим и работоспособным итерационным способом решения подсистем оказался метод хорд

(Ш1) _а(» У „<*> _0<г+1Л РЙЧ»-РЙ, +] .. /=0.1,2................................... (15)

№,)-ВД’)-(їй,--*Й.)

где р^, получено, как и в (13), а

0(1) =о(°) -I — о(°> /У=102

"ск 1 Гск1 Т- ДГ рск1, ¿V їй .

Значение плотности считается достигнутым, как только

0(?+1> _ л(') I

Рек 1 Рек 1 I ^

- <8.

£>(*>

Рек 1

Определив рск1 и /гск1 по заданному /?<°>, находим сі^ и С по формулам

Г „ 1/2 /1 _ -£і-

Ч;______рс

р(0) — Рі

1 +(®Єі (1-“р^7)(СІ8р1 + с*вр» РекіЛ

/ Рі у\{1

4. а 1 / \ Рск 1 / 1

Рі -------------

Г<«

(0)

(16)

Аналогичным образом проводятся расчеты и для второго потока и определяются значения Рск2, Лск2, ctg р2, ¡¡(в0). Через С’ и С<0) обозначены нулевые приближения тангенсов углов наклона первого и второго потоков за скачками уплотнения.

Итак, после прохождения двух расчетных цепочек (11), (12) по значению /?(°> получим

/» + дс(0) = с10) —^0)-

Зададимся теперь некоторым значением давления /?(1) (первое приближение), близким к р(°\ Например, р(1> можно рассчитать по итерационным формулам для скачка уплотнения в совершенном газе, как и в работах [3, б]. Зная давление р^\ проводим, вычисления для первого потока по формулам (13), (15), (16) и по таким же формулам для второго потока. Тогда получим

р(1) + ас(1) = с^1) — с^)-

По двум приближениям для уравнения ДС(/?) = 0 нахо-

дим последующие приближения, применяя метод хорд. Приведем итерационное соотношение, по которому очередное приближение определяется по двум предыдущим приближениям (р(А), ДС(Й)) и (//*+1), д^(*+1>);

р№1=р,,, + дС«.^_г___ *=0, 1,2, ...

Итерации проводятся до тех пор, пока не будет выполнено условие

р(*+1> _;,(*)

р(к)

<е.

Таким образом, если „внутренние“ (по параметру р) и „внешт ние“ (по р) итерации сошлись, то поставленная в начале данного раздела задача решена.

2. Расчет течения с двумя центрированными волнами разрежения. Рассмотрим Случай (фиг. 1, б) возникновения двух центрированных волн разрежения при взаимодействии потоков 1 и 2. Соответствующая этому случаю автомодельная задача с заданными начальными данными рь рь Ух, С, — i%Ъu Pг■, Рг> ^2, = состоит

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

¿С 1 -К; <1р р V3

/М2

Йр

<1р

1

Р

*(/>, Р) Р ’

V 2

■к(р, ?)

VI

(17)

(знак „ — “ в первом уравнении и индекс „/=1“ относятся к потоку 1, а знак „ + “ и „¿ = 2“ — к потоку 2; функции ¡г(р, р) и х(/?, о) заданы, как и в задаче со скачками уплотнения, табличным способом).

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

Решение систем дифференциальных уравнений (17) с заданными начальными условиями и с дополнительным условием на контактном разрыве эквивалентно нахождению некоторой точки пересечения (р, С) интегральных кривых С = Сг(/*) соответствующих систем дифференциальных уравнений для потоков 1 и 2.

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

Систему уравнений для схемы Эйлера можно записать в виде

р=,р + Ар; I = с (р) + А/> - * + / М2 (р)— 1;

р = р(р) + Ьр рт^р] р); * = *(?, р); л = А(Яр);

У2 ¡Г .

Р.(А) = РГ, ^(А) = ^Г, С(/>,) = Сг; М2(/7г)

^?Р*

р/)

; ¿=1, 2,

(18)

где Л/? — шаг интегрирования, а величины параметров на новом шаге обозначены знаком и—м.

Для схемы второго порядка точности на каждом шаге интегрирования проводится пересчет найденных из (18) величин С, р, \ Л, V, М по формулам

р(р + *р) = р(р) +

&р

Р (Р)

+

2 I Р*(Р, Р) 1 (р+&р)

* (р + Д/7) = х \р + Д/7, р (р + Д/7)];

А (р + Д/>) = Л[/> + Д/7, р (/7 + Д/7)];

^2 (/7 + Д/7) = V? + 2 [А, - А (р + Д/7)]; М2(Р + До)= ^(Р±.^-Р^ + АР) .

'Р^аР> {р + Ар).ъ(р + Ар) >

Др/1+С8(р) ,__________

с (/, + Л/0=с (/7)+-£ кмV) -1+

г );

+

\т-УЧр)

1+ь УЩр+Щ-і)>

р(р + Ьр)-У1(р+Ьр) і= 1, 2.

(19)

В данной работе система (17) интегрировалась сначала с некоторым начальным шагом Д/7. Например, можно выбрать начальный шаг по формуле

Д/7

Р{0) — Рі N

N>1, і=1, 2,

(20)

где /7<0) вычислено по формуле (10).

Затем находилось приближенное значение давления /7^°>г на

контактной границе. Далее шаг интегрирования дробился -^-Дру

—Д/7, ... и интегрирование системы (17) повторялось. Итерационный процесс получения /7^>г, р%\, . . . заканчивался при

(т+1) _ (т)

г К. г______“к. Г

У к. г

<е, т = 0, 1, 2, . .

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

Иллюстрация описанного алгоритма дана на фиг. 2, где показана сходимость численных решений, полученных для симметрично разбегающихся потоков совершенного газа, когда интегрирование проводилось по двум различным разностным схемам. На этом графике и далее полагается, что газодинамические параметры обезразмерены: плотность отнесена к плотности набегающего потока р^, давление — К рга 1/шах, компоненты скорости — К 1Лпах (Ушах — максимальная скорость). Наиболее эффективной схемой следует признать схему с пересчетом, и с увеличением углов

поворота потоков целесообразно применять только схему с пересчетом. Измельчение шага Ар влияло на точность расчета в меньшей мере, чем пересчет параметров. Для сравнения результатов отрезок интегрирования разбивался на 5, 10, 20 шагов. Следует заметить, что по количеству итераций созданный алгоритм соответствует алгоритмам для совершенного газа [3], [6], но из-за численного интегрирования систем было заметно возрастание требуемого машинного времени для расчета течения разрежения.

3. Расчет течения с ударной волной и волной разрежения. Примерная картина взаимного влияния двух потоков для случая, когда возникают одновременно скачок уплотнения и волна разрежения, показана на фиг. 1, в и г. В этом случае возможно появление скачка уплотнения как в первом, так и во втором потоках.

0

/—точное решение: 2—расчет, схема с пересчетом; 3—расчет, схема без пересчета (Ы—число шагов интегрирования)

Фиг. 2

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

Для расчета рассматриваемой задачи необходимо решать совместно одну из двух подсистем (6) и систему (17). Как и для задачи п. 1, из линейного приближения задается /А°>. С помощью алгоритма (13), (15), (16) определяется Сн0) (фиг. 1, в) или С(в0) (фиг. 1, г), а по заданному и по значению давления в набегающем потоке находится начальный шаг интегрирования для течения разрежения по формуле (20). Затем проводится интегрирование по р системы (17) для соответствующего потока до конечного значения /7<°> и определяется С<°> (см. фиг. 1, в) или С(н0) (см. фиг. 1, г).

Определяется разность

дс(0) = ;(в0) - с£,0)

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

По двум известным точкам (/?<°>, ДС(0)) и (р(1), ДС(1)), используя метод секущих, определяются следующие приближения:

р(т) _ Р(т+1) Г

,(.«> = рт + дс<* , т - 0, 1, 2, .... (21)

Итерации продолжаются до совпадения двух последовательных приближений с заданной точностью

Замечание 1. Следует подчеркнуть важность согласования точности интегрирования системы дифференциальных уравнений в волне разрежения и точности, с которой находится р: точность итераций (21) не должна быть выше точности интегрирования системы. В противном случае процесс нахождения корня Д£(/?) = 0 может и не быть сходящимся.

х-1

Замечание 2. Использование переменной р 2х вместо р (см. [3]), позволяющей повысить точность расчетов в совершенном газе, в случае равновесного течения невозможно из-за зависимости * от р и р.

Замечание 3. Алгоритмы расчета течения за скачками уплотнения и в веере волн разрежения легко изменить для расчета течения у стенки.

4. Результаты расчетов. Рассмотренные выше алгоритмы были использованы при решении модельных задач о взаимодействии двух равномерных сверхзвуковых потоков совершенного газа (х=1,4) и воздуха, находящегося в состоянии термодинамического равновесия. Численные эксперименты показали надежность работы алгоритмов и созданных на их основе вычислительных программ. Многочисленные проверки и сравнения получаемых решений с результатами расчетов по ранее созданной программе для совершенного газа [6] дали совпадение с заданной высокой точностью.

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

с возникновением двух скачков уплотнения. Расчеты проводились для |С1| = |С2|, которые изменялись от tg 2,5° до, 30°, газ был совершенный (х=1,4), и в том же диапазоне изменения углов поворота потоков проведены расчеты для равновесного воздуха (высота Н=5 км). Влияние физико-химических свойств воздуха заметно проявилось лишь для чисел Мг = М2 = 20 и углов поворота больше десяти градусов (фиг. 3). На этой же фигуре для сопоставления нанесены данные расчетов параметров за скачками уплотнения в воздухе [7]. Незначительные отличия в поведении плотности, которые видны на графиках, можно объяснить различной методикой учета равновесных процессов в воздухе.

У—совершенный газ: 2 — воздух; 3—воздух, расчет [4]

Фиг. 3

На фиг. 4 показаны результаты расчета течения разрежения в совершенном газе и воздухе. Для воздуха в качестве начальных параметров в первом и во втором потоках было взято решение за косым скачком уплотнения с углом поворота потока, равного 25° на высоте 5 км.

Исследование алгоритмов было проведено также на примерах численного расчета элементов плоского сверхзвукового течения. Все такие течения создавались при столкновении двух плоских потоков. В направлении у в поле течения равномерно располагалось 400 ячеек расчетной сетки, граница двух потоков в начальном сечении х = 0 находилась между 200 и 201 ячейками (см. схему на фиг. 5). На фиг. 5 представлены данные численного расчета течения с двумя симметричными скачками уплотнения, где изображено изменение плотности р по автомодельной переменной у/х в двух

/—совершенный газ; 2—воздух

Фиг. 4

/—скачки; 2—линии Мк=1; 5—воздух; 4—совершенный газ Фиг. 6

сечениях х = 1 и 2 (показана только часть симметричной картины течения в потоке 2). Расчетные кривые для течения совершенного газа и воздуха проведены пунктиром, а соответствующие им точные решения, полученные итерационным путем, нанесены сплошными линиями. С увеличением числа ячеек, располагающихся в сжатом слое газа, численные решения все точнее описывают структуру течения.

В заключение рассмотрено коническое течение у треугольного в плане крыла, образованного четырьмя плоскими гранями (фиг. 6). На этой фигуре представлены обработанные результаты расчета пирамидального крыла с углами в плоскости симметрии ш1 = ш2 = = 7° 30', углом стреловидности по передней кромке Х = 70°, углом атаки а = 30° и числом М набегающего потока Моо = 20. Для совершенного газа полагалось х=1,4, а в воздухе параметры набегающего потока соответствовали высоте Н = 5 км. На фиг. 6 показаны форма поперечного сечения крыла, следы ударных волн и конически звуковых поверхностей Мк=К У1 + У%/а=1 (У в, У?— компоненты скорости газа в сферической системе координат с началом координат в вершине крыла, а —местная скорость звука). Очевидно, что наиболее существенно влияние физико-химических свойств воздуха сказалось на росте плотности газа, и, как следствие, заметно значительное поджатие головной ударной волны к поверхности тела.

1. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики. Математ. сб., т. 47, № 3, 1959.

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

3. Иванов М. Я., Крайко А. Н. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений. Ж. вычисл. матем. и матем. физ., т. 12, № 2, 1972.

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

5. Дьяконов Ю. Н. Пространственное обтекание затупленных тел с учетом равновесных физико-химических реакций. ДАН СССР, т. 157, № 4, 1964.

6. Косых А. П., М и н а й л о с А. Н. Исследование методов сквозного счета для задач сверхзвуковой аэродинамики. Ученые записки ЦАГИ, т. 7, № 1, 1976.

7. Нерсесов Г. Г., Бабинова JI. С. Диаграммы параметров потока за скачком уплотнения для воздуха, углекислого газа, кислорода, азота и аргона. Труды ЦАГИ, вып. 1043, 1966.

Рукопись поступила 2jIX 1977