О решении прямой задачи сопла Лаваля

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

Томі 1970 М 1

УДК 533.6.011.55

О РЕШЕНИИ ПРЯМОЙ ЗАДАЧИ СОПЛА ЛАВАЛЯ

В. И. Киреев, Ю. Б. Лифшиц, Ю. Я. Михайлов

Прямая задача сопла Лаваля решается методом установления. Для его реализации предлагается неявная разностная схема. Приводятся примеры расчетов; один из них сравнивается с известным решением обратной задачи [5].

Исследование течений в соплах Лаваля имеет более чем шестидесятилетнюю историю. Первая работа по этому вопросу принадлежит Т. Мейеру [1] и относится к 1908 г. В ней потенциал скоростей был разложен в ряд по степеням декартовых координат в окрестности точки пересечения звуковой линии с осью симметрии. Полученное решение давало «расчетный» режим работы сопла Лаваля с дозвуковым полем скоростей до горла канала и сверхзвуковым — позади него.

Свойства трансзвуковых течений в соплах в окрестности звуковой линии хорошо известны.

Методы расчета течения в целом в сопле Лаваля разработаны в гораздо меньшей степени. Первые результаты в этом направлении получены Ф. И. Франклем [2]. Им построено течение в подводящей части плоского сопла, образованного двумя плоскостями, наклоненными под некоторым углом друг к другу. Границы области влияния, состоящие из прямых линий и двух характеристик, одна из которых вырождается в точку, легко строятся в плоскости годографа скоростей. Полученная задача Трикоми для уравнения С. А. Чаплыгина имеет единственное решение, которое находится в виде суммы частных решений.

Расчет поля потока через сопло произвольной заданной формы рассматривался в работе Я. И. Алихашкина, А. П. Фаворского и П. И. Чуш-кина [3] и в работе А. П. Фаворского [4]. В них точные уравнения газодинамики интегрировались методом интегральных соотношений, предложенным А. А. Дородницыным.

Обратная задача сопла Лаваля, состоящая в определении его стенок по заданному распределению параметров вдоль оси симметрии, решена У. Г. Пирумовым [5].

Ниже исследуется возможность расчета течения во входной части плоских и осесимметричных сопел при помощи метода установления на основе неявной разностной схемы. Для построения схемы используется конечно-разностный метод К. И. Бабенко, Г. П. Воскресенского, А. Н. Любимова и В. В. Русанова [6]. Разностная аппроксимация этого метода берется только для части уравнений, остальные уравнения аппроксимируются иным образом с учетом того, что в состав границ входят линии тока.

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

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

дZ , . дZ . „ дZ х с „ ,

— 4- Л-----[-В— + Р = 0, р* = />.

дЬ дх дг

Здесь А и В — матрицы

А —

и xpa2 0 V 0 xpa2

x~1 p~1 » 0 , 5 = 0 V 0

0 0 u 1p~1 0 V

P vr-1 xp a2v

Z = и , F = 0

V 0

a Z и F— векторы

Давление р, плотность р, скорость звука а и компоненты векторов скорости и и у отнесены к своим критическим значениям; х — показатель адиабаты Пуассона; v = 0 и 1 соответственно для плоских и осесимметричных течений.

Смешанная задача для гиперболической системы (1.1) формулируется следующим образом. На оси симметрии сопла и стенке r = R(x) зададим условие непротекания. В бесконечности вверх по потоку все параметры будем считать постоянными по сечению, но сами они вычисляются в процессе расчета. Невозможность определить их априори Следует из того факта, что в расчетном режиме через сопло Лаваля может быть пропущено в единицу времени только вполне определенное количество газа, которое невозможно заранее задать. Этот вопрос подробно обсуждался в работе [7]. В качестве начальных условий можно взять любое поле, полученное из достаточно разумных соображений, например, по гидравлической теории.

Численное решение сформулированной задачи будем проводить в переменных 1 = 1 (х, г), Т] = Т) (х, г).

Это преобразование осуществляет отображение внутренности заданного сопла на внутренность фиксированного прямоугольника. В переменных g и т] выбирается сетка с постоянным шагом hi по £ и hi по tj. Поэтому преобразование должно удовлетворять определенным требованиям, связанным с удобством проведения расчетов, а также с требуемой точностью решения.

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

l = r/R{x).

Преобразование г] (х, г) выбирается так, чтобы в горле сопла, где производные искомых функций велики, была достаточно густая сетка. Кроме того, оно должно учитывать асимптотическое проведение решения в дозвуковой части при х ->оо. Закон затухания возмущений в этой области определяется формой стенок [7].

2. Численный метод основан на замене исходной системы (1.1) эквивалентной, в которой каждое уравнение есть условие совместности для характеристической поверхности, проходящей через линию t = const,

г) — const. Для ее получения достаточно умножить (1.1) слева на матрицу

1 xpasincp xpacos«

Q= 1 —y.pasincp —xpacos<

0 cos f —sin<p

здесь через 1 -

sin ср = \х

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

В переменных ?, 1\ полученная таким образом система уравнений принимает вид

+,ЪВ)^ -QF,

dt дк <?*)

где Л—диагональная матрица, элементами которой являются временные составляющие характеристических нормалей.

При нормировке (2,1) имеем „

Лп = -)- 1)\г 4* а, (S* 4- Ег)1,

Л22 = и\х — а (S* +-?г)1/2>

• V , +

Если разностную аппроксимацию [6] применить для всей исходной системы (1.1) или эквивалентной ей системы (2.2), то в результате получилась бы разностная схема, по которой можно проводить расчеты только тогда, когда элементы матрицы Л не меняют знака. Такое ограничение вызвано тем, что в этом методе используется четырехточечная неявная схема для аппроксимации левой части (2.2). Более того, для успешного применения четырехточечной неявной схемы необходимо, чтобы элементы матрицы Л не были малыми по абсолютной величине. Для системы (2.2) это особенно ясно, если предположить, что матрица Q постоянна или слабо меняется.

В силу граничных условий непротекания

и%х + vlr—0 при Е = 0 и £=1

и непрерывности Азз имеем, что Лзз->- 0 при | -» 0 или £ -»1. Поэтому для аппроксимации левой части третьего уравнения системы (2.2) возьмем шеститочечную неявную схему. Вследствие (2.3) такая аппроксимация не требует дополнительных краевых условий при 1 = 0 и |=1. По существу это вызвано тем, что границы являются характеристическими поверхностями. Для двух же первых уравнений сохраним четырехточечную аппроксимацию. Обобщение на вторую пространственную переменную выполним аналогично работе [6].

Таким образом, для двух первых уравнений системы (2.2) имеем

д^ "t" ' 1_ .„я + 1 — П-\-\ ггП <7/1 \ ’

... — _ v^ot+1, ! ~г т, I ¿m+l, I ¿-т, I)

Ot /т+1/2 2х

- ^ (Zm+l, i+1 - 2Z”+1,, 4 Znm+1. /_! 4- Z" , ¡+1 - 2Z". , 4- Zm, ,_i),

X

(1г)Л+,/2 =7-[*(2»1и-0+Р(Й+,,,-2;;)],

\ oí / m+1/2, í

я+1/2 ь ,, , ,

2 Г ~ i 7/2+1 yn+1 <7/1 + 1 \ ,

— ~— Ia v¿m+l. í+1 4" ¿m, ¡+1 — ¿m+1, í-1 — ¿m, г-l) 4 m+1/2, i 4t

+ Р(^я+1,1+1 + 1+1 —2«+l,l-i - Z¡|, j-l)],

• Z"V/;l: = -j

(m«=0, 1....AÍ- 1; /=1,2,.../.).

Для последнего уравнения { dZ '\п+1/2 1 / 7п+1 -7Л \ °2 ^2 / 7» о 7п . \ 7п . \

( —г— )-------------{¿т, I ¿т, // \¿*m, Z+1 “Г ¿т, 1—\)ч

V dt )т, і х *

Г “Т+1/2 =“[а (z”+i-1 - z"->- +р (Z™+i-¡

V Я /Я, г 2г <?z y+i/2

Ч <?•»] /от. г

Z/I+1/2 1 / 7П+І , <7/1 ч

m, I — ~Z~

(

^ MZ^I+I - Z"t/-i) + p (Z". i+1 - z".,_!)],

2x .

(m = 0, 1,...,,M; Z=1,2,...,I),

здесь £г = х//гг, <x>p, a + p=l, >0;

M к L~ числа интервалов по 5 и tj; х —шаг по времени (t = т).

3. После подстановки формул (2.2) и (2.5) в систему уравнений (2.2) получим систему нелинейных разностных уравнений. Для ее решения воспользуемся итерационным процессом. Точно так же» как и в [6], все величины со слоя п -f 1, входящие в коэффициенты уравнений и их правые части, возьмем с і-й итерации. Тогда для расчета искомых величин на (/ + 1)-й итерации на каждом луче т) = const получается замкнутая система линейных разностных уравнений, которую запишем следующим образом:

amZm-i + bmZm + cmZm+1=dm

Здесь через ат, Ьт и ст обозначены матрицы;

—1/2 (xpa \ Sin cpm—1/2 (xpaXcos cp )m_

= 0 0 0

0 —

1 (xpa sin cp)m_ 1,2 (xpa COS <?)m-1/2

bm = 1 — (xpa sin cp)m+l/2 — (xpa cos cp)m+1/2

0 (2 cos cp)m - (2 sin cp)m

0 0 0

V-m+112 —( *pa¡A Sin ф)т+1/2

0

— (xpaji. cos<p)m+i/2 — («А, Л33 sin ср)от = (14- 2a.kx Л,2)/(1 — 2a.k1 Л2г).

(akx Л33 cos cp)m

X = (1 —2&k1 Ли)/(1 4- 2aftj Лп), jx >

В правую часть dm включены все члены, не зависящие от

величин на (i 4- 1)-й итерации.

Все уравнения на луче = const, которые не вошли в систему

(3.1), совместно с уравнениями (2.3) возьмем в качестве граничных

условий

goZo + SiZi = d0,

ём-1 Zm-i + ём Zm — dM-

Матрицы g¡ в этих формулах даются выражениями

1 — (хрд Sin <p)i/2 — (xpa cos <p)i/2

go = 0 (sin <p)0 (cos <p)0

0 (cos (p)0 — (sin ep)a

it

g 1 =

Єм-і'

V-m — (xpap. sin 9)1/2 — (xpajA cos <p)i/2

0 0 0

00 o

Хлі-1/2 (xpaX sin <Р)лі-1/2 0 0

0 0

(xpaA COS f)M-l/2 0 0

1 (xpa sin <р)ж-1/2 (xpa eos <?)M-V2 gM= 0 (sincp )m (cosfbi

0 (cos<f)M — (sin<pV Краевая задача (3.1) — (3.3) будет хорошо обусловленной (см. [8]), если |Х|<1 и |fJ.|<l. Для выполнения этих условий необходимо и достаточно, чтобы

Лп^>0, Л22<^0

для всех рассматриваемых значений i и rj. Указанные неравенства имеют место, если поток вдоль координаты g дозвуковой; этого всегда можно добиться соответствующим выбором вида функции | (х, г).

Решение задачи (3.1) — (3.3) осуществляется методом прогонки, который детально исследован в работе И. Д. Сафронова [9].

Последний луч I = L выбирается так, чтобы нормальная составляющая вектора скорости на нем всегда была больше скорости звука. Производные по г) на этом луче берутся по формулам односторонних разностей.

На луче 1 — 0 считается, что и = 0, а и и р связаны интегралом Бернулли и постоянны по Для определения их значений на каждом шаге по t будем требовать, чтобы расход газа в сечении 1 = 0 и в горле сопла был одинаковым. Начальный луч с / = 0 выбирается на некотором достаточно большом расстоянии вверх по потоку от горла сопла.

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

На фиг. 1 построены линии уровня приведенной скорости % — const поля течения в круглом сопле, образующая которого задавалась формулой

R W = i—£L + T*f

где y :

ch2x

= 0 при х<0 и у = 0,2 при х>0.

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

На фиг. 2 сплошными линиями нанесены кривые % = const в сопле, расчет которого проведен в [5]; пунктиром отмечены линии, полученные в [5], а кружочками показаны результаты расчета течения в выхлопной части методом характеристик. Некоторое расхождение объясняется, вероятно, неточным заданием контура сопла.

Специально проведенный численный эксперимент показал слабое влияние положения начального луча I = 0 на течение в горле сопла.

* *

*

ЛИТЕРАТУРА

1. Meyer Т. N. Uber zweidimensionale Bewegungvorgenge in einem Gaz das mit Überschallgeschwindiskeit strömt. Forschung heit, Hf. 72, 1908.

2. Франкль Ф. И. Истечение сверхзвуковой струи из сосуда с плоскими стенками. ДАН СССР, т. 58, № 3, 1947.

3.. А л и х а ш к и и Я. И., Фаворский А. П., ЧушкинП. И. О расчете течения в плоском сопле Лаваля. Журнал вычислительной математики и математич. физики, т. 3, № 6, 1963.

4. Фаворский А. П., Расчет сопел Лаваля. Журнал вычислительной математики и математич. физики, т. 5, № 5, 1965.

5. П и р у м о в У. Г. Расчет течения в сопле Лаваля. Изв. АН СССР, Механика жидкости и газа, 1967, № 5.

6. Бабенко К- И., Воскресенский Г. П., Любимов A. H., Русанов В. В. Пространственное обтекание гладких тел идеальным газом. М., «Наука», 1964.

7. Лифшиц Ю. Б., Рыжов О. С. О вариации расхода газа в расчетном режиме работы сопла Лаваля. Журнал вычислительной математики и математич. физики, т. 6, № 2, 1966.

8. Рябенький В. С. Необходимые и достаточные условия хорошей обусловленности краевых задач для системы обыкновенных разностных уравнений. Журнал вычислительной математики и математич. физики, т. 4, № 2, 1964.

9. Сафронов И. Д. О методе прогонки для решения разностных краевых задач. Журнал вычислительной математики и математич. физики, т. 4, № 2, 1964.

Рукопись поступила 21 /IV 1969 г.