Решение прямой задачи для плоского сопла Лаваля релаксационным численным методом по схеме мурмена-коула Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ И А Г И

Том XII

198 1

М з

УДК 532.525

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

Э. Г. Шифрин, М. А. Шулаков

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

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

В настоящее время имеется ряд методов численного решения прямой задачи сопла Лаваля, например, методы [3, 4], использующие схемы установления [5, 6]. К ним следует отнести также схемы работ [7, 8]. Кроме того, в приближенной постановке решение прямой задачи можно находить среди множества решений обратной задачи [9].

Однако методы установления являются медленно сходящимися. Существенно увеличить скорость сходимости при численном решении уравнений эллиптического и смешанного типа позволяет использование метода верхней релаксации [10], а также разработанных в последнее время попеременно-треугольного метода [11] и методов приближенной факторизации [12, 13].

В настоящей работе в развитие метода [1], апробированного ( на расчетах внешнего трансзвукового обтекания в рамках уравнений малых возмущений, предлагается метод решения прямой задачи, обладающий следующими характерными чертами:

1. Задача сведена к краевой задаче для одного функционально-дифференциального урайнения второго порядка. По сравнению с методами [3, 4, 7, 8], где решается система уравнений, это существенно экономит машинную память.

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

3. Используемая аппроксимация дифференциального уравнения имеет погрешность второго порядка на гладких решениях.

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

5. Применение метода блочной релаксации [10] в сочетании с многосеточным методом [14] для решения системы алгебраических уравнений на каждой итерации обеспечивает достаточно высокую скорость сходимости итерационного процесса.

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

(1)

д'о й2 йф

где М — число М, X — коэффициент скорости, угол наклона вектора скорости, ки й2 — коэффициенты Ламэ системы координат

Для коэффициентов /ц и /г2 справедливы соотношения:

(здесь х — показатель адиабаты, — произвольные функ-

ции от ср и <]>, определяющие нумерацию линий соответствующего семейства).

Выберем (ср) = Х(®, фс) = Л(<р), тогда на стенке сопла (<|> = фс) = 1 и, следовательно, <р — длина дуги стенки сопла. Положим далее /*;2 (ф) = 1.

Дифференцируя перекрестно по <р и ф уравнения системы (1) и исключая р, получим

X2 — 1 . ] _ Г Л (?) (X)

Л (<р) Хт*^-1) (X) 9 или в окончательном виде

і+1 ГХ4 - 1

ф

Х(Х2-1) Х9¥ - Л2 (?) Хт-1 (X) Хфф + [^1 + 2] Х2? +

н[

X2 — 1 т (X)

Л2 (?) (X) 4 - ЦХ* - 1) М. х¥ = 0. (2)

Л (9)

Уравнение (2) имеет смешанный эллиптико-гиперболический тип. Математическая формулировка прямой задачи сопла Лаваля

бесконечной длины для уравнения (2) будет следующей: в полу-полосе G (— оо <р С; 0 < ф •<:' фс) найти решение, удовлетворяющее условиям:

Хф = 0 при 6 = О,

1 k (?) ,

Ч = ^-1) — при ф = фс,

lim Х=:>>ос< 1,

ср ——оо

>->1 при о = С,

где k (?) — кривизна контура сопла.

Постоянная С определяется как значение ? на выходе из сопла.

Условие для Хф при ф = фс, полученное из второго уравнения системы (1), должно выполняться во всех точках контура, в которых кривизна у. (с?) непрерывна.

Для сопла конечной длины, которое и рассматривается в настоящей работе, вместо условия lim к -< 1 ставится условие

ср->—оо

JTW^5fX? = 1,(W ПРИ V = B<C’

где у] (Ф) кривизна ликии равного потенциала со = В. Это условие фактически задает распределение угла наклона вектора скорости р на входе в сопло.

При численном решении уравнения (2) используются идеи метода, предложенного^ в работе [1] для решения квазилинейного уравнения малых возмущений Ф^ФЛГЛ-=ФУГ Метод предполагает использование двух конечно-разностных схем аппроксимации производных для областей эллиптичности и гиперболичности решаемого дифференциального уравнения. Аналогично работе [1] в области эллиптичности уравнения (2) применяются центральные разностные отношения второго порядка точности для производных по f и ф, а для области гиперболичности — односторонние разностные отношения для производных по <р

1 __2\\ /— i, j + 2, j h-з. /

— (Д<Р)3

.1 34 i ~ Щ-l,- i + h-2, j

■ ~ 2Д?

и центральные разностные отношения для производных по ф (шаблоны для производных показаны на рис. 1), где i и / — соответственно номера линий <р = const и ф = const (1 -< iN\ 1 j^CjN).

■г ' ; 1

Рис. 1

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

2(1/ + £>£/)Х(-н> - 5/(№+_!) + №+!)) =

= + (1 _ ш) [2 (V + пи) т- 1/(Х(«)_1 + х(.»)+1)],

где

I *-f 1

к—X?

[D = 1; г'==.Г.(Дф)а + £/(хВД + ^)1./) при. Х{?/« 1,.

[D = - 1;.■.2 = W(ЦУ -U(5ХВД-Ч4№?} + ХЙ&) при Х$ > 1,

1F =

xf <*> -1 7|------------------f-2

'2 (л)

■Am

V ,2 (n) « ■ '

X?,<") — 1

:(XW.) 1 ~J''? Ь(Х<«).) +2.

}(n)

;(«) 2(n) ___ -.4 <?(l.jN) i(ra)

Л!. / I i» I 1) , (n) V *

A«,//V

Ш — параметр релаксации, принимающий значения

1<<о<2 при Х< 1,

to < 1 при X > 1,

п — номер итерации.

Полученная система решается методом прогонки вдоль линий <р = const. Решение начинается с линии <з> = В и заканчивается на линии ф = С, где ХС/.>1. Используемые конечно-разностные формулы обеспечивают устойчивость прогонки [15].

Число итераций, необходимое для получения решения с заданной точностью, определяется в основном скоростью сходимости в области эллиптичности. Последняя, в свою очередь, зависит от значения параметра релаксации со для этой области, которое при численном решении задачи выбирается равным оптимальному значению со0 при решении задачи Дирихле для уравнения Пуассона в прямоугольнике [11]

1 + "УЛШШ (2 — -A-mjn) где Amin — минимальное собственное значение задачи:

2 (Дф)3 sin2 — + 2 (Дер )2 Sin2 ~

Amin =-----------------------------;

2(A+)2sin2— + (Дф)2

здесь ik — номер линии f, выходящей из точки контура сопла, в которой у = ymin.

Значение параметра релаксации <о для точек расчетной области, где величина коэффициента скорости Х>1, выбирается из условия сходимости итерационного процесса, т. е. ш<1 (в рассматриваемых примерах со = ш0/2).

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

расчетной сетки Д<р и ДФ, т. е. испоЛьзуется многосеточный метод решения.

Начальное распределение коэффициента скорости X на линиях у = const задается в виде многочлена Х(Ф) = Ах ф4 + А,, коэффициенты которого определяются из условий:

Фс

/ (Л, Г+ А2) (Х0) d'b = Х0 ^<*-1) (Х0) Фс,

О

(А\ ф4 -М2)ф = ПРИ Ф“Фс.

где фс = Уо:^оо (Х00), у0— ордината входного сечения сопла, Х0 — значение коэффициента скорости, вычисленное по одномерной теории (Х00 —во входном сечении сопла).

Окончательно распределение коэффициента скорости X вдоль линий ср = const можно представить как

ч±) , *.Шс Г/ ^ V 1 I

В процессе решения задачи реализация граничного условия на стенке сопла (Ф = Фс) производится в виде

Х(,п)(ср ф ) = -___1М__________

Входящая в выражение для Хф (<р, Фс) кривизна контура сопла х(<р) в равноотстоящих точках определяется заранее при аналитическом или табличном задании контура.

Положение граничной линии тока (ty = Фс) должно соответствовать контуру сопла в физической плоскости ху. Для этого в процессе итераций необходимо поддерживать постоянство ординаты входного сечения сопла у0, которая определяется из условия *-/

При интегрировании .методом трапеций

Kj = ~U=\, j=jN); /С,= 1(1</</Л0.

Из полученного выражения для у0 определяется величина ДФ на следующей итерации, которая зависит от номера итерации и стремится к некоторому предельному значению.

На каждой итерации расчет ведется до линии <рг = const, для которой Х,-,;-> 1 и положение которой может меняться в процессе счета. Если в расчет включается' новая линия <р.+1, то начальное распределение коэффициента скорости на ней задается в виде

h+i,j = 2h,i—h-i,j'

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

1) X/f = XB = const на входном участке;

2) монотонный разгон от ХК = \В до ^=1 на участке с углом наклона вектора скорости р = const;

3) X^. = 1 на участке до критического сечения сопла.

/л?

it I»'

(Х)]|

Ц + О (Дф)2.

~6 ~й -4і -З ~1 ~1хО Рис. 2

При этом реализуется прямая звуковая линия в критическом сечении.

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

На рис. 1 и 2 приведены контуры двух сопл с углами {3 = 60° и 90° соответственно (кривая 3), начальные условия при их профилировании (кривая 1) и результаты численного решения прямой задачи для этих сопл в виде распределений коэффициента скорости X вдоль оси (кривая 2) и стенки сопла (нанесены точками). При профилировании и решении прямой задачи сопла Лаваля течение предполагалось изоэнтропическим с показателем адиабаты х=1,4.

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

Максимальное отличие распределений коэффициента скорости X вдоль контура сопла с р = 60° (см. рис. 1) при использовании расчетной сетки 77X29 (69X29 —в дозвуковой области) составляет 1,6%, а среднее отличие по контуру, исключая точки разрыва второй производной, составляет 0,58%. Среднее отличие и среднеквадратическое отклонение распределений коэффициента скорости в критическом сечении сопла равны 0,54 и 0,59% соответственно.

Установление решения с погрешностью г<10~7 (где г —максимальная из невязок в четырех контрольных точках), при использовании на начальном этапе расчетной сетки 39X15 с последующим уменьшением шагов А? и Д’1> вдвое, достигается примерно за 450 итераций (из них 300 итераций на сетке 39 X 15), при этом время счета на ЭВМ БЭСМ-6 составляет примерно 7,5 мин.

ЛИТЕРАТУРА

1. Мурмен Е., Коул Дж. Расчет плоских установившихся трансзвуковых течений. .Ракетная техника и космонавтика", 1971, № 1.

2. П о д с ы п а н и н а Н. А., Шифрин Э. Г. Об одном методе профилирования коротких плоских сопл. „Изв. АН СССР, МЖГ“, 1975, № 1.

3. Иванов М. Я., Крайко А. Н. Численное решение прямой задачи о смешанном течении в соплах. „Изв. АН СССР, МЖГ“, 1969, № 5.

4. Киреев В. И., Лифшиц Ю. Б., Михайлов Ю. Я. О решении прямой задачи сопла Лаваля. .Ученые записки ЦАГИ', т. 1, № 1, 1970.

5. Г о д у н о в С. К. О численном решении краевых задач для систем линейных обыкновенных дифференциальных уравнений. ,УМН‘, XVI, № 3, 1961.

6. Бабенко К. И., Воскресенский Г. П. Пространственное обтекание гладких тел идеальным газом. М., .Наука", 1964.

7. Л а в а л ь П. Нестационарный метод расчета трансзвуковых течений в соплах. В сб. „Численные методы в механике жидкостей*. М., „Мир\ 1973.

8. Васенин И. М., Рычков А. Д. Численное решение задачи о смешанном осесимметричном течении газа в некоторых криволинейных областях методом установления. „Изв. АН СССР, МЖГ“, 1971, № 1.

9. П и р у м о в У. Г. Расчет течения в сопле Лаваля. „Изв. АН СССР, МЖГ», 1967, № 5.

10. В а з о в В., Форсайт Дж. Разностные методы решения дифференциальных уравнений в частных производных. М., Изд. иностр. лит-ры, 1963.

11. Самарский А. А., Николаев Е. С. Методы решения

сеточных уравнений. М., „Наука", 1978. \

12. Бальхауз В., Джеймсон А.^Дж. Альберт. Неявный метод приближенной факторизации для рёшения стационарных трансзвуковых задач. „Ракетная техника и космонавтика", 1978, № 6.

13. Холст Т., Бальхауз В. Консервативные методы быстрого расчета полного уравнения для потенциала скорости в применении к трансзвуковым течениям. „Ракетная техника и космонавтика", 1979, № 2.

14. Ф е д о р е н к о Р. П. О скорости сходимости одного итерационного процесса. „Ж. вычисл. матем. и матем. физ.“, т. 4, № 3, 1964.

15. Годунов С. К., Рябенький B.C. Введение в теорию разностных схем. М., „Физматгиз", 1962.

Рукопись поступила 18jl 1979 г.

\