CC BY
31
Том XXIII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 1992
№ 2
УДК 533.6.011.55 : 532.582.3
К ВОПРОСУ О НАЛИЧИИ ОБЛАСТЕЙ НЕУСТОЙЧИВОСТИ СТАЦИОНАРНЫХ ТЕЧЕНИЙ С УДАРНЫМИ ВОЛНАМИ В СОВЕРШЕННОМ ГАЗЕ (ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ)
И. В. Савин, И. И. Шитиков
Приведены результаты вычислительного эксперимента по проверке {с позиции научного принципа повторяемости) существования якобы обнаруженных областей неустойчивости стационарных течений совершенного газа с сильной ударной волной перед тупым телом в определенном диапазоне значений чисел Маха набегающего потока и показателя адиабаты.
Проведены расчеты обтекания конуса со сферическим затуплением с помощью метода установления по времени, путем численного интегрирования нестационарных уравнений Эйлера, записанных в дифференциальной форме в сферической системе координат с подвижным центром. Для конечно-разностной аппроксимации использовалась неконсервативная форма явной двухшаговой схемы Маккормака. Основное внимание в работе уделяется процессу установления по времени к стационарному решению. '
Проблема теоретического изучения устойчивости ударной волны уже на протяжении многих лет остается в сфере внимания газодинамиков и механиков. В последние годы отмечается некоторая активизация работ в этом направлении, в печати стали появляться все новые и новые публикации.
Среди них, в частности, встречаются такие, которые по сути своей, с одной стороны, являются весьма интересными, но, с другой стороны, изложенные в них результаты настораживают и требуют проверки. К ним можно отнести работы [1] и [2].
В работе [1] В. П. Федосов проводит теоретический анализ уравнений Эйлера с целью выяснения причин появления неустойчивости головных скачков высокой интенсивности и делается вывод о том, что при определенном фиксированном значении числа Маха набегающего потока существует критическое значение показателя адиабаты хкр, при котором головной скачок еще остается устойчивым, а для состояний газа с меньшими значениями х скачок становится уже неустойчивым.
В публикации [2] Г. А. Тарнавский путем прямого численного моделирования, исследуя процесс установления головной ударной волны около затупленного конуса, определяет границы области устойчивости скачка по значению х при заданном числе М. Так, отмечается, что для М„=10 неустойчивость головного скачка с уменьшением х наступает в интервале значений 1,15Сх< 1,2, а для М», = 5 критическое значение хкр, ниже которого происходит деструкция головной ударной волны, лежит в диапазоне 1,05 < <х < 1,10.
Несомненный интерес представляет работа [3], в которой получено интегро-дифференциальное уравнение, описывающее поведение поверхности скачка уплотнения и найдено в явном виде решение этого уравнения. Одним из следствий анализа полученного в явном виде решения является вывод о том, что в совершенном газе ударная волна абсолютно устойчива относительно малых возмущений. Утверждение относительно устойчивости ударной волны в совершенном газе полностью соответствует выводам, которые следуют из ранее проведенных исследований [4—6].
Лишний раз в справедливости этого вывода авторы убедились, выполняя с помощью метода установления по времени численные расчеты по изучению обтекания совершенным газом с различными показателями адиабаты х= 1,4; 1,3; 1,2; 1,1; 1,05 при числах Маха набегающего потока М», = 5 и 10 конуса с полууглом раствора 0 = 25° и сферическим затуплением. Одной из целей проведенного исследования была проверка с позиций научного принципа повторяемости результатов работы [2].
Отметим, что подобные расчеты с вариацией значения показателя адиабаты х выполнялись и ранее, и не одним автором (см., например, [7]). Как по результатам вновь проведенных исследований, так и по результатам, полученным ранее другими авторами, можно сделать вывод, что результаты численного эксперимента, представленные в работе [2], не могут являться подтверждением выводов, сделанных в работе [1].
1. Были проведены расчеты обтекания совершенным газом под нулевым углом атаки осесимметричного тела вращения, представляющего собой конус с полууглом раствора 0 = 25° и сферическим затуплением, при числах Маха набегающего потока = 5 и 10 и варьировании показателя адиабаты в пределах ] ,05 ^ х ^ 1,4.
Расчеты проводились с помощью комплекса АРГОЛА [8], который позволяет получить решение, используя метод установления по времени, путем численного интегрирования нестационарных уравнений Эйлера, записанных в дифференциальной форме в сферической системе координат с подвижным центром. Для конечно-разностной аппроксимации используется неконсервативная форма явной двухшаговой схемы Маккормака.
Счетная область, для которой строится разностная сетка, ограничивается поверхностью тела, для которой ставится условие непротекания; внешней границей, являющейся головной ударной волной, положение которой определяется в процессе установления граничными условиями Ренкина — Гюго-нио; и выходной границей, выбираемой заведомо лежащей в зоне сверх; звукового течения и не требующей постановки граничных условий. С целью наблюдения процесса установления к стационарному решению при фиксированном значении числа Маха и варьировании показателя адиабаты расчеты выполнялись на разностной сетке, имеющей 10 интервалов (ЫЫ) вдоль лучей от тела к волне и 17 или 34 интервала вдоль поверхности тела от оси симметрии к выходной границе (Л^/7) (рис. 1). Для уточнения газодинамических параметров внутри поля и оценки влияния вариации числа узлов в расчетной области на процесс установления проводились расчеты на сетке с параметрами NN X NF — 20 X 34.
2. Основное внимание в численном эксперименте уделялось процессу установления по времени. Сходимость результатов во время расчетов контролировалось проверкой выполнения законов сохранения энергии (в виде константы Бернулли) в каждой расчетной точке, массы и импульса — интегрально, стационарности положения головного скачка, гладкости всех газодинамических функций.
Процесс установления изучался на примере эволюции значения плотности р = р/р,„ к стационарному решению в следующих контрольных точках в поле течения: 1 и 2 — лежащих соответственно на теле и скачке на оси симметрии; 3 и 4 — лежащих на поверхности тела и головном скачке уплотнения в месте их пересечения с лучом, выходящим из центра сферического затупления и перпендикулярным оси симметрии (см. рис. 1). Плот-
—=^о д)
1)
Рис. I
ность выбиралась в качестве контрольной функции потому, чтО она, как известно, является одним из наиболее медленно устанавливающихся газодинамических параметров.
Для ускорения процесса установления комплекс АРГОЛА позволяет осуществить следующие вычислительные технологические приемы.
Первый: при переходе от варианта к варианту по х перед началом расчета в стартовом приближении для узлов, находящихся на поверхности тела, проводился пересчет значения плотности, исходя из условия сохранения энтропии при вновь заданных параметрах набегающего потока вдоль линии тока, приходящей в критическую точку и «омывающей» всю поверхность тела. Значение энтропии в случае совершенного газа для такой линии тока нетрудно вычислить, если учесть, что головной скачок уплотнения при обтекании затупленных тел локально, около того места, где его пересекает линия тока, затем попадающая в критическую точку, можно рассматривать прямым.
После пересчета в стартовом поле значений плотности на теле проводилось сглаживание поля плотности вдоль лучей в направлении от тела к ударной волне по типу схемы Лакса: <р, = (ф,_| + 2<р, + <р,+ 1)/4 при условии сохранения плотности на волне.
С подготовленного изложенным выше способом стартового поля затем начиналось установление по времени.
Второй: в начальный период установления в дозвуковой части течения задается принудительное втекание газа через поверхность тела, которое автоматически отключается при выходе на стационарное решение поверхности ударной волны.
Ниже в тексте при описании вариантов расчета обтекания затупленного конуса без особых оговорок, везде подразумевается, что при переходе от варианта к варианту использовался первый прием.
Все расчеты, в которых варьировались значения показателя адиабаты как для Моо=10, так и для М,* = 5, выполнялись при однотипном переходе от варианта к варианту. В качестве начального поля для варианта с каким-либо фиксированным значением х бралось поле, полученное для ближайшего и большего значения х, например для варианта с х= 1,3 стартовым полем служило поле течения, рассчитанное для варианта со значением х = 1,4, и т.д.
Для получения решения при х= 1,4 начальным приближением выбиралось поле течения около сферы с аналогичными параметрами набегающего потока.
Отметим, что ниже при анализе результатов расчетов газодинамические параметры представляются в безразмерном виде и обезразмеривание произведено следующим образом:
р = —--------давление; р = —-------плотность,
Р« К» р~
где Ктах определяется из соотношения Бернулли
V2 о V2
9 тах К "» ао
~2~~ х-1 Р„ Г'
а индекс «оо» относится к параметрам набегающего потока.
На рис. 2 отражен процесс установления плотности для вариантов расчета с числом Маха набегающего потока М», = 10 и вариацией показателя адиабаты х= 1,2; 1,1 и 1,05.
Комплекс АРГОЛА позволяет осуществлять продвижение по времени с минимальным во всем расчетном поле шагом, определяемым из условия Куранта — Фридрихса — Леви. Однако на практике расчет проводится с несколько меньшим шагом, который определяется заданием коэффициента запаса СОЕР, меньшего единицы. В серии расчетов для М00= 10 коэффициент СОЕР равнялся 0,7, расчетная сетка бралась размером 10 X 17 (ЫР/Х Л^). Для этой серии расчетов установление решения с точностью не менее 0,5% для всех газодинамических параметров происходило в основном за первые 500 итераций по времени, исключая вариант х= 1,05, причем, как и следовало ожидать, за скачком быстрее, чем на теле. Продвигаясь в расчетах от варианта к варианту в сторону уменьшения х, интегральная погрешность в законах сохранения массы и импульса на счетных лучах, близких к выходной границе, увеличивается, напррмер для М.» = 10 и х= 1,05 она достигает 15%. Это говорит о том, что, хотя по мере уменьшения х головной скачок и приближается к телу (см. рис. 1), градиенты параметров течения в направлении, перпендикулярном поверхности головного скачка, возрастают, а область их изменения существенно уменьшается. Следовательно, расчетная сетка 10X17 становится достаточно грубой и для достижения более высокой точности результатов требуется переход на более мелкую или неравномерную расчетную сетку. Поэтому для варианта М„ = 10, х= 1,1 специально был осуществлен переход на сетку 20 X 34 и пройдено 3000 итераций по времени — погрешность результатов по всем контрольным параметрам не превышала 1%.
х-/,2
Численное решение с параметрами набегающего потока М,* =10, х = = 1,05 (см. рис. 2) получено с использованием процедуры, ускоряющей процесс установления плотности на теле. Выход решения на установившееся значение носит более выраженный колебательный и длительный характер. Если плотность за скачком (кривая 1) достигла своего стационарного значения через 500 шагов по времени, то на теле (кривая 2) установление достигается примерно через 2200 итераций. Затянувшийся процесс установления в данном случае, по-видимому, связан с неудачным заданием параметра, отвечающего за величину «втекания» газа через поверхность. Когда после первых 1000 итераций этот параметр был подкорректирован, выход решения на установление существенно ускорился.
На рис. 3 представлен процесс установления значений безразмерной плотности р в контрольных точках 1, 2 и 3, 4 для 1АХ = 5 и сетках с параметрами X Л7/г= 10 X 34 и 20X 34 с постоянным шагом по времени.
При переходе от варианта к варианту по х для ускорения процесса установления использовался первый из описанных выйе вычислительных тех-нологическйх приемов. Расчет каждого варианта с различными х выполнялся в несколько этапов по 500 шагов. После каждой серии в 500 шагов проводилось обновление шага по времени, а именно: при СОЕР = 0,8
выбирался минимальный во всем поле шаг по времени, вычисленный из условия Куранта — Фридрихса — Леви. Таким образом, каждый из этапов рассчитывался со своим (в общем случае) шагом по времени, однако, как показал анализ значений шага на каждом из этапов (особенно после первых 500 шагов, когда наблюдается относительно бурный процесс установления), они отличались друг от друга незначительно.
В случае варианта с х = 1,3 установление с простой машинной точностью (т. е. когда переставали изменяться пять точных цифр после запятой) для контрольных точек / и 2 наступило фактически после 1500 шагов по времени, а установление с точностью 0,1% — к 500-й итерации. Что касается контрольных точек 3 и 4, то установление четырёх значащих цифр после запятой в точке 3 произошло к 1600-й итерации, а в точке 4 — к 1300-й итерации.
С уменьшением значения х наблюдалась общая тенденция, когда граница установления сдвигается в сторону увеличения количества шагов по времени. Так, в случае варианта с х= 1,1 установление четырех значащих цифр после запятой в значении плотности в точках / и 2 произошло уже
к 3500-й итерации, в точке 3— к 4350-й итерации, а в точке 4— к 3500-й. Установление по времени решений для х = 1,4, 1,3, 1,2, 1,1 осуществлялось на расчетной сетке с параметрами NN X NF= 10 X 34.
Следует отметить, что попытка в рамках сложившейся методики на этой сетке перейти от варианта с х = 1,1 к варианту с х = 1,05 оказалась неудачной. При этом наблюдалось следующее явление: при относительно установившихся значениях плотности в контрольных точках на волне {2 к 4) и в критической точке /, около поверхности конуса в районе 28-го и 29-го лучей возникали осцилляции плотности, амплитуда которых в процессе установления возрастала, и, в конце концов, когда плотность принимала отрицательное значение, происходила аварийная остановка программы. Наступал так называемый «развал» решения.
Поэтому было решено осуществить переход от х= 1,1 к х = 1,05 более «плавно», т. е. через промежуточное решение, полученное для х= 1,07 на расчетной сетке NN X N1? = 20 X 34.
Процесс установления значения плотности в контрольных точках для варианта расчета с х= 1,05 представлен на рис. 3. Как показывает анализ решения при х= 1,05, при ‘ отсутствии явных признаков катастрофических явлений, которые могли бы повлечь за собой «развал» решения, когда в целом поля газодинамических параметров гладки (рис. 4), наблюдается следующее. Установление значения плотности р в контрольной точке 2 с простой точностью пяти значащих цифр после запятой произошло к 2000-му шагу, установление с точностью четырех цифр в контрольной точке 4 (на волне) произошло к 4500-му шагу. Однако установление в контрольных точках
■Р
№-22
V
0,14
0,12
0,10
0,08
0,06
0
5
10
15
NN
6)
Рис. 4
на поверхности тела значительно затягивается, так, к 6000-му шагу произошло установление лишь трех значащих цифр после запятой, далее при установлении наблюдается медленное изменение пятого знака после запятой, влекущее за собой еще более медленное изменение четвертого знака. Настораживающим является то, что проверка выполнения закона сохранения интеграла Бернулли после очереднего этапа дает увеличение, хотя и незначительное, ошибки в узлах на поверхности тела и в узлах, ближайших к ним.
По-видимому, такую ситуацию следует объяснить тем, что на процессе установления уже начинают сказываться не физические явления, которые моделируются в численном эксперименте и установление которых в основном закончилось, а чисто программные эффекты. Ведь каждую программу можно рассматривать как инструмент, предназначенный для исследования какого-либо физического явления, и, в принципе, возможна такая экстремальная ситуация, когда точности инструмента начинает не хватать и начинает ощущаться влияние самого инструмента на исследуемое явление. Как показал опыт с вариацией параметров сетки, в данном численном эксперименте с х=Л,05 имеет место сходное явление, когда Для такого низкого значения к реализованный в программе алгоритм учета граничных условий непротекания на теле для используемой расчетной сетки не обеспечивает, по-видимому, достаточной точности.
В заключение следует остановиться на одной интересной и, на первый взгляд, неочевидной газодинамической особенности, попутно наблюдавшейся при вариации х. На рис. 1 представлено положение головной ударной волны в зависимости от значения показателя адиабаты. Характерным является то, что с уменьшением показателя адиабаты отход ударной волны при фиксированном числе Маха уменьшается. Кроме того, звуковые линии (М = 1) при различных местах своего начала на волне «дружно» заканчиваются почти в одном и том же месте на сферической части затупленного конуса. Из' известного газодинамического соотношения, связывающего критическое давление, соответствующее звуковой линии (М = 1), с давлением торможения, имеем
X
£*-= (-^) *+1 » 0,525 -Ь 0,605, если х = 1,4 1
(\\т—= eTY
\х— 1 Ра /
С другой стороны, распределение давления вдоль поверхности тел гладкой формы можно приближенно определять по формуле Ньютона:
Р*/Ро = COS2 со*.
Из этих двух соотношений вытекает, что положение звуковых точек, определяемое углом ш*, для тел с гладкой формой поверхности слабо зависит от условий обтекания и от показателя адиабаты х. Например, для сферы и цилиндра центральный звуковой угол <о* = 36 -г- 41 ° в широком диапазоне условий обтекания при М^ ^6 [9]. Поэтому при численной дискретной аппроксимации поверхности звуковые точки на теле для вариантов с разными х с точностью до ячейки будут совпадать, что и наблюдается на рис. 1.
Подводя итоги, следует сказать, что авторы данной статьи лишний раз убедились в справедливости выводов работ [3—6], из которых следует, что ударная волна в совершенном газе, характеризующемся введенным в работе
rci • dV I Р — Р» 1 /,, 1 \
[5J параметром I *= —\я у _у = — — ^ V= —j , всегда устойчива.
Что касается результатов работы [2], то они не убеждают в том, что причиной развала решения является неустойчивость головного скачка уплотнения. Скорее всего, как показывает внимательное изучение графического материала этой работы, более вероятной причиной развала решения является неудачная реализация краевых условий на поверхности тела.
ЛИТЕРАТУРА
1. Федосов В. П. О наличии областей неустойчивости течений с ударными волнами в совершенном газе // Препринт ИТПМ СО АН СССР.—
1985, № 21. Ч. 1.
2. Тарнавский Г. А. О наличии областей неустойчивости течений с ударными волнами в совершенном газе // Препринт ИТПМ СО АН СССР.—
1985, № 22. Ч. 2.
3. М и х а й л о в Ю. Я., Ч и н и л о в А. Ю. Построение решения линейной задачи об устойчивости ударной волны // Изв. АН СССР, МЖГ.—
1988, № 4.
4. Д ь я к о н о в С. П. Об устойчивости ударных волн // ЖЭТФ.—1954.
Т. 27, вып. 3 (9).
5. Иорданский С. В. Об устойчивости плоско# стационарной ударной волны // ПММ.— 1957. Т. 21, вып. 4.
6. Конторович В. Н. К вопросу об устойчивости ударных волн //
ЖЭТФ.— 1957. Т. 33, вып. 6(12).
7. Благосклйнов В. И., Минайлос А. Н. Обтекание кругового цилиндра сверхзвуковым потоком совершенного газа // Ученые записки ЦАГИ.— 1972. Т. 3, № 2.
8. Михайлов Ю. Я., Юмашев В. М. Комплекс АРГОЛА: автоматизированный расчет гиперзвукового обтекания ЛА // Материалы VII Всесоюзного семинара по комплексным программам математической физики.— Новосибирск: 1982.
9. Лунев В. В. Гиперзвуковая аэродинамика.— М.: Машиностроение,
1975.
Рукопись поступила 2/Х 1990 г.
