УЧЕНЫЕ ЗАПИСКИ ЦАГИ Том XVIII 1987
№ 3
УДК 532.525.011.5
О ПРОФИЛИРОВАНИИ плоского СВЕРХЗВУКОВОГО СОПЛА, ОБЕСПЕЧИВАЮЩЕГО РАВНОМЕРНЫЙ ПОТОК В ВЫХОДНОМ СЕЧЕНИИ
В. Н. Котеров, И. Л. Осипов, В. П. Пащенко
Предлагается метод профилирования короткого плоского сверхзвукового сопла по распределению скорости, заданному на неизвестной до решения образующей сопла. При этом требуется равномерность потока в выходном сечении. Уравнения газовой динамики решаются в переменных «функция тока—потенциал скорости» конечно-разностным способом, заключающимся в комбинации итерационного метода и метода характеристик.
Несмотря на интенсивное развитие численных методов решения прямой задачи сверхзвукового сопла (см., например, [1—3]), задача о профилировании сопла, обеспечивающего равномерный поток в выходном сечении, по-прежнему привлекает значительное внимание, в частности, в связи с развитием новой техники [4, 5].
Задачу профилирования сверхзвукового сопла обычно разделяют на две независимые проблемы: построение дозвуковой [6, 7] и сверхзвуковой (например, [8]) частей сопла, ограничиваясь тем самым классом сопл с плоской звуковой линией и угловой точкой. К существующим методам сквозного профилирования следует, по-видимому, отнести лишь подходы, основанные на задании параметров потока на известных до решения задачи линиях тока, например, на оси сопла [9]. В подобных постановках задача профилирования сопла формулируется как задача Коши, которая в дозвуковой области течения обладает условной корректностью. Кроме того, при этом подходе трудно обеспечить безударное и безотрывное течение у стенЪк короткого сопла, а также добиться равномерного течения на выходе без чрезмерного удлинения сверхзвукового участка сопла.
В данной статье излагается методика построения короткого плоского сверхзвукового сопла, основанная на решении уравнений, описывающих потенциальное течение невязкого газа в переменных «функция тока — потенциал скорости» [7].
1. Рассмотрим безвихревое стационарное течение идеального совершенного газа в схематически изображенном на рис. 1, а плоском симметричном сопле. Пусть при —сю поток газа приближается к равно-
Рис. 1
мерному с числом М0<1. Требуется определить профиль сопла (АВС) с известной полушириной выходного сечения, обеспечивающий равномерный поток с заданным числом М1>1 в сечении х = х^ Ниже везде индекс «О» относится к параметрам, определенным при х-*—°° (при практической реализации методики — параметрам на некоторой расположенной достаточно далеко вверх по потоку линии х=хг, а индекс «1» — к параметрам на срезе сопла при х=х1).
Если в качестве обезразмеривающих параметров выбрать максимальную скорость звука ат и давление торможения р*, то зависимости модуля скорости газа Ш и удельного расхода б от местного числа М примут вид
янм)= ■ О)
[1+0,5 (т-1) М2]
[1 +0,5(7— 1)М2]^
где у — отношение удельных теплоемкостей газа.
Погонный расход газа т через плоское сопло дается выражением
от =■ 2Т /?* Ь/ат,
где величина Ь, имеющая размерность длины, может быть рассчитана по данным в сечениях х=хг или х — хи где реализуются равномерные потоки
1 = (М0) = ЬхО (Мг) . (2)
Введем безразмерную функцию тока ч]? и потенциал скорости ср:
й$ = Ь~1 в (М) (сое & йу -- йх) ; йу = С (?) V/ (М) (сов Ъ йх вш & йу)
,!
(3)
где '& — угол наклона вектора скорости к оси х. Пусть на оси сопла ■ф = 0, тогда, согласно (2), стенке сопла соответствует значение г|), равное 1. Гладкую функцию С(ф) можно выбирать произвольным образом, например, ее удобно ввести в форме
С (*) = £/(?)/ИР (М) \^
(4).
где (М) |ф=1—заданное распределение скорости на образующей сопла, а и (ф) — произвольная гладкая функция. Если положить О (<р) = = 1/Ь*, где Ь*— предполагаемая длина сопла, то вдоль стенки, вследствие (3) и (4), будет выполняться соотношение й(р = с11/Ь* (й.1—дифференциал длины стенки) и, таким образом, устанавливается линейная связь между величинами ф и /. Если же выбрать [/(ф)=^соп$1:, то эта связь становится нелинейной. При численном решении рассматриваемой задачи это позволяет получать неравномерную разностную сетку в физических переменных х, у, сохраняя равномерную дискретизацию расчетной области в плоскости ф, г|э. В настоящей статье выбираем ¿/(Ф) = 1ДЛ
Отметим, что отображение плоскости (х, у) на плоскость (ф, г|з), во-первых, всегда однозначно и, во-вторых, трансформирует исходную криволинейную область в прямоугольник, что несколько упрощает численное решение задачи.
В независимых ортогональных переменных ф, ч|з уравнения неразрывности и безвихренности записываются в следующей форме:
— = - Л(М)— ;
^■=.В (М) ^ ;
(5)
(6)
А( М) =
(I — М3)/.
I* № О (М)^-2 IV (М) ’
В( М) =
\У (М)
М2Д(М)Т I. ’
О (М) = 2 [1 + 0,5 (т - 1) М2]!“1 .
Из уравнений (5) и (6) можно получить уравнение второго порядка относительно числа М:
Система (5) — (6) и соответственно уравнение (7) являются эллиптическими при М<1 и гиперболическими при М>1. В последнем случае тангенсы углов наклона С±-характеристик к оси ф и соот-
ношения вдоль этих характеристик представляются следующим образом:
У - В (М)/Л (М) = Ф (М) , )
(8)
= + V — А (М)-В (М) ¿Ш. ]
Рассматриваемая область течения в плоскости ф, г); схематически изображена на рис. 1,6. Точками помечена звуковая линия ВО, слева от которой поток газа имеет дозвуковую скорость. Зону сверхзвукового течения разобьем на две части криволинейной характеристикой С_, выпущенной из точки £), которая, находится на оси сопла и в которой М(£))=М1. В области АРОЕ формулируется задача профилирования
дозвукового АВ и «разгонного» ВР участков сопла по заданному распределению числа Маха Мс (/) на стенке ЛF [7]. Это достигается посредством решения уравнения (7) с последующим определением #(ф, я|)) из уравнения (5) — в том числе и на неизвестном профиле АР. В «выравнивающем» участке сопла, ограниченном справа прямолинейной характеристикой С+, на которой § = 0 и М=Мц, поставим задачу Гурса [8] для характеристической системы (8). Следовательно, форма сопла полностью определяется распределением Мс(ф) =Мс[/(ф)] на «ускоряющей» части стенки Л/\
Переход в физическую плоскость (х, у) при рассчитанных полях М(ф, -ф) и ^(ф, гр) осуществляется по формулам, обратным выражениям (3):
здесь фг — эквипотенциаль, принимаемая при вычислениях за входной срез сопла, а хт—координата точки на оси сопла, соответствующая значению ф = фг (см. рис. 1, а).
2. Решение поставленной задачи осуществляется методом конечных разностей. Дискретизация расчетной области ACHE (см. рис. 1, б) проводится равноотстоящими линиями ф = const и г(з = const; шаги разностной сетки обозначаются ниже через h(p и hty соответственно. В дозвуковой области уравнение (7) аппроксимируется на пятиточечном шаблоне «крест» (рис. 2, а) центральными разностями второго порядка точности. В «разгонном» участке BFDG (рис. 1,6) используется семиточечный несимметричный шаблон (рис. 2,6). Отказ от аппроксимации
Г sin »
1 J G (М) ^ ;
о
(9)
о
L,J i.j+t
i't.J
в)
/' \
Рис. 2
уравнения (7) со вторым порядком точности по переменной ф связан с условной устойчивостью и немонотонностью известных нам схем подобного типа в сверхзвуковой области. Наиболее употребительная схема Мермена [10], имеющая второй порядок аппроксимации вдоль потока, была разработана для расчета обтекания тонких профилей, которое характеризуется малыми градиентами давления в трансзвуковой области. Применение схемы [10] к расчету внутренних течений [3] сочеталось с тщательным выбором начального приближения и небольшими значениями величины М — 1 (расчет заканчивался при выполнении на оси сопла условия М=1). Использование «косого» шаблона первого порядка аппроксимации по ф не обнаруживает неустойчивости даже при больших значениях М1. Наличие же семи точек в шаблоне обусловлено тем, что аппроксимация второго члена в уравнении (7) при М>1 производится центральными разностями на слоях ф^ и ф,, а затем результаты суммируются с весами, равными р и 1—р соответственно, где
г==_}_ 0<г<1 ,
Ф АТ } (10)
1, г > 1 .
Параметр г представляет собой долю шага йф, отсекаемую С±-харак-теристиками на линиях грг-1 и гр+1 (рис. 2,6), а величина Ф дается формулой из (8). Из (10) следует, что при г> 1 (т. е. для больших значений М) семиточечный шаблон вырождается в пятиточечный и, следовательно, схема расчета становится явной. Это позволяет избежать асимптотической неустойчивости, присущей неявным схемам [11], а введение весовой аппроксимации (при г<1) обеспечивает плавный переход от неявной к явной схеме при увеличении М.
В выравнивающем участке сопла РСВ применяется обратный метод характеристик второго порядка точности [12], основанный на соотношениях (8) (см. также рис., 2, в).
3. Решение разностных уравнений проводится итерационным способом. Начальное приближение М во всей расчетной области полагается равным распределению Мс(ф) вдоль всей стенки АС. На каждой текущей итерации рассматривается система линейных алгебраических уравнений, полученная заменой функций А (М.) и В(М) [см. (7)] их значениями с предыдущей итерации, а само распределение М с предшествующей итерации привлекается для определения способа аппроксимации уравнения (7) в каждом узле сетки. Для разрешения этой системы применяется метод релаксации с прогонкой по линиям Ф = сопэ1 [13]; прогонки проводятся последовательно в направлении возрастания координаты ф и дают в результате некоторое поле М' (ф, о|:). При этом на оси сопла и на входе используется условие дМ/дп = 0 (п — нормаль), которое является следствием уравнений (5), (6) и ра-> венства нулю угла 'О на этих границах. Новое распределение числа М(А формируется следующим образом:
М<*) = м<*-1> + ш (М«6-1), я'*“1)) (М' — М(*_1)) . (11)
Здесь верхний индекс указывает номер итерации, а через со обозначается параметр релаксации [13], зависящий от величины М(А_1)и итерационной ошибки я:
5(*) = | 1 — !> | .
Функция со имеет вид:
1, М(*_1)>1
Ь9_ м(*-1) < 1.
И +2
При этом представлении о» итерационная ошибка уменьшается
монотонно, т. е. без скачков, характерных для расчетов методом верхней релаксации с постоянным или кусочно-постоянным поведением со.
Процедура последовательных прогонок заканчивается определением точки й на оси сопла, в которой М(й) (О) = Затем в области AK.DE (ем. рис. 1,6) проводится интегрирование линейного относительно # уравнения (5) в направлении поперек потока, начиная от оси, где 8- ^0. По рассчитанным на данной итерации полям М и Ф выстраивается С--характеристика .О/7 и вычисляются распределения М и й на £>/\ Наконец, с помощью условий однородности потока на прямолинейной С+-характеристике СЕ> рассчитываются значения М и Ф в расчетных узлах, содержащихся в криволинейном треугольнике /-’СД путем решения задачи Гурса. В частности, определяется поведение числа М вдоль участка ЕС и тем самым в выравнивающей части профиля сопла корректируются заданное начальное распределение Мс(ф) и длина сопла ¿*.
Итерационная процедура заканчивается при выполнении условия
max (sft) < в,
где е — заданная точность сходимости, а максимум берется по всей расчетной области.
4.. Предложенный выше алгоритм был реализован в виде программы на ЭВМ БЭСМ-6, написанной на языке АЛГОЛ-бО. На рис. 3 и 4 представлены некоторые результаты проектирования плоских сопл для ■у =1,48. Распределение Мс(<р) вдоль стенки во всех вариантах задавалось по формуле
Мс(?) =
[(1—М0)е*-А + М0 —ср<0, 1+0?, ?>0; A = 0,5-L*u 0 = 1,8-¿ь
где /-1 —длина дозвукового участка сопла, а величиной д определяется скорость возрастания числа М в «разгонной» части.
На рис. 3, а приводится рассчитанный контур сопла с М0=0,11, М1 = 2, ¿1=6, ¿* = 26. При этом размер сетки составлял =
= 20X200 ячеек, а машинное время, необходимое для достижения итерационной точности е= 10~4, — 15 мин.
Поскольку, как уже отмечалось, исходное распределение [в данном случае (12)] в сверхзвуковой части сопла автоматически корректируется и тем самым постановка задачи изменяется в ходе итерационного процесса, то зависимость итогового решения от количества узлов сетки довольно сложна и, в некотором роде, нетрадиционна. В связи с этим был проведен ряд методических расчетов по профилированию сопл с очень короткой дозвуковой частью, чтобы возможно большая часть узлов попала на разгонный участок.
Рис. 4
В таблице сведены данные указанных расчетов: длины сверхзвуковых частей сопл с ¿* = 6, Мо = 0,8, ¿* = 1, М.1 = 2, ¿1 = 0,75:
Сетка ¿2 ¿3
5X50 3,160 3,160
10X100 2,853 3,103
15X150 2.685 3,071
В первой колонке таблицы дается размер сетки (первая цифра означает число ячеек поперек сопла), во второй — длина Ь2 сверхзвукового участка образующей, полученная при тривиальном начальном приближении М(<р, ф)=Мс(<р), а в третьей — та же длина Ь3 для расчетов, в которых в качестве начального приближения использовалось численное решение задачи на более крупной сетке (указанной в предыдущей строке таблицы).
/
Рис. 5 о
Резкое отличие длин Ь2 от Ь3 объясняется тем, что, несмотря на идентичное распределение Мс(ф) для всех сеток, итоговое распределение Мс(ф) в сверхзвуковой части сопла отличается от исходного (см. рис. 3,6). На рис. 3,6 сплошной линией нанесен исходный профиль Мс(ф), а штриховой — окончательный для сетки 5x50. Вертикальная линия ф = ф* = 0,134 отделяет разгонный участок от выравнивающего. Заметно расхождение графиков уже в разгонном участке. Это является следствием итерационного перемещения С_-характеристики Т7/) (см. рис. 1,6) с соответствующим «исправлением» разгонного участка. Если же начальное приближение выбирается посредством предварительного расчета на более крупной сетке, то результаты оказываются близкими (см. третью колонку таблицы) и различие между соседними вариантами уменьшается при измельчении сетки. Это происходит из-за постепенного изменения постановки задачи: в качестве исходного распределения Мс(ф) в каждый расчет закладывается распределение Мс(ф), полученное в результате расчета на предыдущей, более крупной сетке.
На рис. 4 представлено сравнение сверхзвуковых частей сопл с М! = 3 (у=1,4), реализующих равномерный выходной поток. Сплошная линия соответствует контуру сопла с минимально возможной длиной — с угловой точкой [14], штриховая — соплу со скругленным по дуге окружности изломом [15], а штрихпунктирная — рассчитанному по предложенной методике. Отличие длин последних двух сопл • составило — 17%, при этом итоговое распределение Мс(ф) (верхняя кривая на рис. 4) достаточно благоприятно.
В заключение отметим, что данный метод легко обобщается для проектирования осесимметричных сопл на основе метода [7]. На рис. 5 сплошной линией нанесен профиль осесимметричного сопла, испытанного в работе [16], а крестиками — контур стенки, восстановленный расчетом по измеренному в этой работе распределению давления вдоль стенки. Штриховой линией обозначена полученная в результате расчета звуковая поверхность. Отличие двух профилей — реального и восстановленного— лежит в пределах гарантированной в работе [16] точности измерений.
ЛИТЕРАТУРА
1. Киреев В. И., Лифшиц Ю. Б., Михайлов Ю. Я- О решении прямой задачи сопла Лаваля. — Ученые записки ЦАГИ, 1970, т. 1, № 1.
2. Годунов С. К., Забродин А. В., Иванов М. Я-, К р а й-к о А. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики. — М.: Наука, 1976.
\
\
3. Ш и ф р и н Э. Г., Ш у л а к о в М. А. Решение прямой задачи для плоского сопла Лаваля релаксационным численным методом по схеме Мур-мена—Коула. — Ученые записки ЦАГИ, 1981, т. 12, № 3:
4. Лосев С. А. Газодинамические лазеры. ■—М.: Наука, 1977.
5. Бабаев И. К., Басов Н. П., Г л о т ов Е. П. Расчет оптимальных режимов генерации непрерывных сверхзвуковых электроионизационных СО-лазеров. — ДАН СССР, 1983, т. 270, № 3.
6. П о д с ы и а н и н а Н. А., Ш и ф р и н Э. Г. Об одном методе профилирования коротких плоских сопл. — Изв. АН СССР, МЖГ, 1975, № 1.
7. Осипов И. Л. Численный метод построения двумерных сопл. — Изв. АН СССР, МЖГ, 1979, № 2.
8. Ш мыглевский Ю. Д. Некоторые вариационные задачи газовой динамики. — М.: Труды ВЦ АН СССР, 1963.
9. Пир у м о в У. Г. Расчет течения в сопле Лаваля.—Изв. АН СССР, МЖГ, 1967, № 5.
10. Murman Е. М., Cole J. D. Calculation of plane steady transonic flows. —AIAA J, 1971, vol. 9, N 1.
11. Самарский А. А. Теория разностных схем. — М.: Наука, 1977.
12. Кацкова О. Н., Чушкин П. И. О единой схеме численного метода характеристик.-—ДАН СССР, 1962, т. 154.
13. Вазов В., Форсайт Д. Разностные методы решения дифференциальных уравнений в частных производных.—М.: Изд. иностр. лит., 1963.
14. Кацкова О. Н., Наумова И. Н., Шмыглевский Ю. Д., Шулишнина Н. П. Опыт расчета плоских и осесимметричных сверхзвуковых течений методом характеристик. — М.: ВЦ АН СССР, 1961.
15. Крайко А. Н., Шелом овский В. В. Сравнение двух методов профилирования контуров сверхзвуковых частей сопл, реализующих равномерный поток. — Ученые записки ЦАГИ, 1981, т. 12, № 4.
16. В а с k L. Н., М a s s i е г P. F., G i е г Н. L. Comparison of measured and predicted flows through conical supersonic nozzles with emphasis on the transonic region. — AIAA J., 1965, vol. 3, N 9.
Рукопись поступила 6/Xl 1985