CC BY
59
________УЧЕНЫЕ ЗАПИСКИ и А Г И
Т о м XI 19 80
Л* 4
УДК £32.525.2:532.517
МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ О РАСПРОСТРАНЕНИИ СВЕРХЗВУКОВОЙ НЕДОРАСШИРЕННОЙ ТУРБУЛЕНТНОЙ СТРУИ В СПУТНОМ СВЕРХЗВУКОВОМ ПОТОКЕ
В. И. Копченое
В рамках упрощенной системы уравнений Навье — Стокса с привлечением дифференциального уравнения для турбулентной вязкости рассматривается задача о распространении сверхзвуковой недорасширенной турбулентной струи в спутном сверхзвуковом потоке. Для численного интегрирования системы уравнений построена полунеявная двухшаговая конечно-разностная схема.
1. Экспериментальному исследованию сверхзвуковых вязких недорасширенных струй, истекающих в спутныи сверхзвуковой поток, посвящены работы [1, 2|. В силу эллиптичности полной системы уравнений Навье — Стокса, описывающей стационарное течение, численное решение задачи с использованием современных ЭВМ получено в [3] лишь в сравнительно небольшой области, примыкающей к срезу сопла. В то же время протяженность первой „бочки“ в ряде случаев составляет десятки и даже сотни радиусов среза сопла. Поэтому заслуживает внимания подход [4], основанный на численном интегрировании упрощенной системы уравнений Навье — Стокса, являющейся параболической по продольной координате при сверхзвуковом компоненте скорости в указанном направлении. Методы расчета ламинарных струй в рамках упрощенной системы уравнений развивались в [4—7|. С применением численного метода [4—6| решена задача о турбулентной струе [8]. При этом система уравнений для осредненных параметров замыкалась с помощью приближенной интерполяционной формулы для турбулентной вязкости. В настоящей работе используется дифференциальная модель для турбулентной вязкости, предложенная в [9].
Основные требования к численному методу диктуются особенностями газодинамической структуры течения, схема которого
изображена на рис. !. Наличие достаточно сложной системы скачков уплотнения делает весьма привлекательными методы сквозного счета. Последние нашли широкое применение наряду с методами, включающими аккуратное построение разрывов [10]. Однако разрешающая способность конечно-разностных схем сквозного счета в значительной мере зависит от свойства монотонности, о чем свидетельствуют результаты, представленные, например, в [11]. Вместе с тем необходимо, чтобы погрешности аппроксимации были малы в областях, где существенна вязкость. Опыт использования модели турбулентности (см. [12]) показывает, что более предпочтительными являются монотонные конечно-разностные схемы, а точность модели в ряде случаев допускает применение методов первого порядка.
В настоящей работе за основу принята конечно-разностная схема сквозного счета для стационарных сверхзвуковых течении [13], которая ранее была распространена на случай вязких ламинарных струй в [7]. Следует отметить, что в методах, предложенных в 14—6, 8] и в [14], скачки уплотнения либо выделяются, либо возникает необходимость в искусственном размазывании ступенчатого распределения давления в сечении среза сопла даже при умеренных степенях недорасширения.
2. Рассмотрим плоскую или осесимметричную сверхзвуковую струю, истекающую в спутный сверхзвуковой поток. Пусть X н у — оси прямоугольной системы координат (ось х направлена по оси симметрии в осесимметричном случае или лежит в плоскости симметрии в плоском случае). Газ считается вязким и теплопроводным, а режим течения — турбулентным. Для описания осредненного потока воспользуемся упрощенной системой уравнений Навье — Стокса (см., например, [7|), которая может быть представлена в следующем виде:
др</ , др1>
дх
д\
У ’
5, - р«0
д (рца + р) ^ дриу дзх ,
дх ду ду У
д(ли ' д ((/Ч- + р)
дх
дьи (2/ 4 а>3) дх
ду
Ору (2/ 4 а,а)
ду
даУ , ,
ду' у
= о
ду
Ь/
2д — ?у (2; 4- а’г)
(2.1)
где
рг 4-1* ди Re ду ’
2 ?* - н-
3 Ре
/о ду . у \ , _ ,, рг - [* I ду_______________1<_\
I ду у)’ у " Ие I ду у /’
1 / ;* \ д1 , рг 4 !* Г ди 2 /п ду .у \1
~ Ре ( ~Рг Р'т ) ду _Г Ие [“ Ну ~~ Т С' ( 17 "•/7)] '
(2.2)
Система уравнений (2.1) и (2.2) намыкается с помощью дифференциального уравнения для турбулентной вязкости в приближении пограничного слоя (9]
О'АП
дх
л. д'*т —
ду
_ д / рм + р 0« \ ] / РХ«+Ц дг \
ду \ Яе ду ) у \ Ре ду )
ди
ду
дх
(2.3)
Здесь ю—модуль вектора скорости го, проекции которого на оси л: и у соответственно и и V; р — давление, ¡> — плотность, удельная энтальпия, ц — динамический коэффициент вязкости, / = О или 1 в плоском и осесимметричном случае соответственно,
/—головной скачок уплотнения; 2—разделительная линия гока; 3—висячий скачок* углотнения; «/—отраженный скачок уплотнения: о—границы слоя
смешения: 6—веер волн разрежения
Рис. 1
е —турбулентная вязкость. Отметим, что а, х, ; — полуэмпиричес-кие постоянные, присущие модели турбулентности, используемой в работе (вообще говоря, величина г — функция отношения г к кинематическому коэффициенту ВЯЗКОСТИ V).
Все величины считаются безразмерными. Последнее достигается следующим образом. Пусть /.*, ¡>*, ги# — характерные размерные значения длины, плотности, скорости. Тогда все координаты относятся к /„*, плотность — к р*, компоненты скорости — к го», давление — к и& удельная энтальпия / — к ■£'?, динамический коэффициент вязкости 1* — к своему значению на срезе сопла, V и г — к Числа Ие и Рг определяются по характерным размерным величинам. В дальнейшем за принимается радиус выходного сечения сопла, а в качестве р*, выбраны критические
значения плотности и скорости в сопле.
Удельная внутренняя энергия е, удельная энтальпия I, скорость звука а — известные функции давления и плотности. Для совершенного газа с показателем адиабаты ?
а =/^7р, е = р!(1— 1) р, ¿ = е + />/р*= ?/>/(?— 1>Р-
Предполагается, что проекция и вектора скорости на ось л: всюду в рассматриваемой области превосходит скорость звука а.
Тогда исходная система уравнений параболична в направлении х. В начальном сечении должны быть заданы распределения компонентов скорости и и V, давления, плотности, турбулентной вязкости или эквивалентных им параметров. При у -* оо все параметры стремятся к своим невозмущенным значениям в спутном потоке. В расчетах указанные граничные условия сносятся на некоторую, достаточно удаленную от оси линию у(х) = У, — такую, чтобы не вызывать существенных возмущений в поле течения. На оси выставляются условия симметрии.
3. В работе [7| для численного интегрирования системы уравнений (2.1) —(2.2) была предложена явная конечно-разностная схема первого порядка точности. В процессе выполнения настоящего исследования эта схема была опробована применительно к турбулентным снутным струям. Обладая рядом достоинств, связанных с расчетом разрывов, метод имеет существенный недостаток. Вследствие явной формы конечно-разностного представления „диссипативных“ членов, содержащих вторые производные от параметров в поперечном направлении, при малых эффективных числах Рейнольдса !?£',= Ре/(г + ч) проявляется зависимость шага интегрирования кх от Рег
Соображения, приведенные в |7], и результаты методических расчетов показывают, что для Ие,— 100 шаг интегрирования в направлении х приходится уменьшать в некоторых случаях в 20 раз по сравнению с шагом, выбранным из условия устойчивости для
невязкой задачи. Последнее имеет вид Нх—V М* — 1 Ау, где Мх — число М, определенное по продольной составляющей скорости, а Ау —размер разностной сетки в направлении у. В турбулентном слое смешения величина £ может достаточно сильно нарастать по длине струи. В результате машинное время, потребное для расчета турбулентных струй, может значительно превосходить время решения аналогичной задачи в рамках идеального газа при одинаковых прочих определяющих параметрах. Отметим, что ограничение на шаг 1гх для невязкой задачи является естественным, так как отражает связь масштабов неоднородности идеального сверхзвукового потока в продольном и поперечном направлениях. Модификация метода [7], предпринятая в настоящей работе, позволяет избежать зависимости шага интегрирования Нх от числа Ие,. При ‘ построении нолунеявной двухшаговой конечно-разностной схемы использовался опыт применения метода’„предиктор—корректор“ в задачах одномерной нестационарной газовой динамики (15].
Модифицированный метод включает следующие элементы. На первом этапе с помощью неявной конечно-разностной схемы определяются вспомогательные величины на новом слое. На втором полушаге осуществляется пересчет по конечно-разностной схеме [7], однако для аппроксимации вязких членов используются вспомогательные распределения параметров на новом слое. Будем применять обычный способ (см., например, [7]) обозначения и нумерации величин, осредненных по сторонам элементарной четырехугольной ячейки, изображенной на рис. 1. Параметрам на известном слое с абсциссой х0 присвоим нижние полуцелые индексы, а на слое х0 + Их — верхние. В качестве примера представим конечно-разностные соотношения, следующие из уравнения количества движения в проекции на ось х:
первый полушаг
(/?Ю* + (ЯЮ*_, ик - ик_х
2
~~Рк—12
второй полушаг
X
(3.2)
2
2
Здесь нижний индекс р приписывается вспомогательным величинам, определяемым на первом полушаге; большими буквами обозначены осредненные значения параметров на верхней и нижней сторонах ячейки; Ну.к-т и Лу-12 — высоты боковых граней ячейки, а йу,т = (Л5~12-{-Л5+,г)/2. „Большие“ величины Н,и, V, Р, как и в [13|! находятся из решения плоской задачи о взаимодействии двух равномерных сверхзвуковых потоков. Параметры на новом слое с целыми индексами вычисляются, как полусуммы соответствующих величин с полуцелыми индексами. Например, принимается, что =* = (г*-12 + г*+Ч2)/2.
Получающаяся на первом полушаге алгебраическая система уравнений для ир решается методом скалярной прогонки. Заметим. что в уравнении движения в проекции на ось х содержится заранее неизвестная величина др/дх. Поэтому приходится использовать итерационный процесс (нижний индекс / указывает, что значение функции определяется на предыдущей итерации). Аналогичным образом рассматриваются остальные уравнения системы
(2.1) — (2.3) за исключением уравнения неразрывности. Алгебраическая система, отвечающая уравнению для турбулентной вязкости, получается нелинейной относительно неизвестных величин £* из-за диффузионного члена. От нелинейности можно избавиться, если аппроксимировать коэффициенты при соответствующих производных в уравнении (2.3) с помощью величин г, и р(., найденных на предыдущей итерации. Эти же значения г( и р/ используются для представления в конечных разностях „диссипативных“ членов
(2.2). В результате выполнения первого полушага рассчитываются распределения параметров гр, ир, Vр, ¡р.
Второй этап основан на алгоритме, предложенном в [7), с изменениями, касающимися конечно-разностной записи ад., ау, 3^, д.
Так, например, в уравнении (3.2) величина I* определяется следующим образом:
1 — Ц (Рг)»+1/2 ~ 'Хк + Ц2 ~ (?е)*-1 2 + :хЬ-1-2
“*4-1/2 ~ “*-12 Лу. т, к
2
где 12<6<1.
Численное интегрирование уравнения неразрывности осуществляется но методу [13|. Найденные по результатам второго этапа значения параметров р, г, р используются на следующей итерации в качестве р,, а также на первом неявном полушаге при вычислении в конечных разностях производной др/дх. Сходимость итераций контролируется по двум относительным величинам 5, = | (р — р,)1р| и ?_>'=’ (Е—е|)/*1- В качестве начального приближения для р, и г, выбираются значения параметров р и е с известного слоя, а производная дрдх заменяется конечно-разностным выражением
Построенная разностная схема аппроксимирует исходную систему дифференциальных уравнений с первым порядком, она устойчива при выполнении условия того же типа, что и для невязкой задачи. Как показывают примеры расчетов, на разрывах решение ведет себя монотонным образом. Вместе с тем модифицированная схема имеет ограничения по сравнению с явной схемой \1\. Во-первых, приходится запоминать больше информации, например, значения параметров, найденные на промежуточном этапе. Во-вторых, необходимо совершать итерации, что в сочетании с возрастанием числа арифметических операций приводит к увеличению времени решения задачи. Кроме того, при достаточно больших степенях нерасчетности результаты не могут быть получены с помощью модифицированной схемы из-за расходимости итераций. Поэтому в разработанном алгоритме предусмотрена возможность переключения с явной схемы на двухшаговую и наоборот, что позволяет решать задачу в достаточно широком диапазоне определяющих параметров.
4. Перейдем к описанию некоторых результатов расчета. Пусть осесимметричная сверхзвуковая струя вытекает в спутный сверхзвуковой поток. В начальном сечении каждый из потоков равномерен и параллелен оси х. В качестве основных определяющих параметров приняты числа М на выходе из сопла М0 и в набегающем потоке Мг, степень нерасчетности п — ра1ре, отношение температур Та/Те. Индексы а не служат для обозначения величин на срезе сопла и в спутном потоке. Расчет проводится с выделением линии тока у = ус(дс), выходящей с кромки сопла. В области между осью струи и линией тока ус(*) используется равномерное разбиение на ячейки. В полосе _ус У размер ячейки линейно
увеличивается но направлению к верхней границе. Постоянные, входящие в уравнение (2.3), выбраны в соответствии с результатами работы [9|: а = 0,2, х = 2, ^ = 0,5. Принимается, что Ргг=0,7.
Рис. 2
На рис. 2 представлены распределения давления по поперечной координате для случая п — 5, М„ = 3, М, = 4, Та/Те = 1 на расстоянии л-=15 и 40 от среза сопла (все координаты относятся к радиусу выходного сечения сопла». Число Re0, определенное по параметрам на срезе сопла, равно 104. В расчетах варьировалось количество ячеек. Сплошной кривой соответствует конечно-разностная сетка, для которой на область между осью струи и линией тока _у = ус(х) приходится 40 ячеек, штриховой — 80, кривой, обозначенной точками, —20. Указанное число ячеек составляет приблизительно 14 часть от их общего количества. Наибольшее различие результатов наблюдается в окрестности головного скачка уплотнения. При уменьшении размера ячейки пространственное размазывание скачка уменьшается. Аналогичные выводы следуют также из рассмотрения распределений продольного компонента скорости и по координате у.
Результаты, изображенные на рис. 3, получены для сверхзвуковой струи с Ма = 2,56. распространяющейся в снутном потоке с М, = 3,1. Отношения статических давлений и полных температур на срезе сопла и в набегающем потоке равны соответственно 4,7 и 1,6. Число Re0 принимается равным 10\ На внешней и внутренней стенках сопла имеются пограничные слон толщиной 1,3 и 0,4. Пограничный слой моделируется следующим образом. Профиль скорости задается по степенному закону с показателем 1/7. Счи-
тается, что поля перепадов температуры торможения и относительной скорости подобны, а величина турбулентной вязкости выбирается в соответствии с профилем скорости по данным, представленным в |16|. На рис. 3 показано распределение относительной избыточной температуры торможения Г0-= (Г0 — Т0е),(Т0а—Т0е) по координате у в сечениях х = 14 и 24. Там же приведены соответствующие распределения величины р[>!р0е, где р0— давление торможения за прямым скачком уплотнения в случае, когда параметры перед скачком совпадают с параметрами в рассматриваемой точке. Кружками изображены данные из экспериментальной работы [2].
В методе расчета предусмотрена возможность, когда молекулярные массы и теплоемкости газа в струе и в спутном потоке различаются. При этом к системе уравнений (2.1) —(2.3) добавляется уравнение для массовой концентрации газа, вытекающего
нз сопла. Кроме того, ряд дополнительных членов присутствует в уравнении энергии (см., например, (17|). Необходимо также учитывать переменность ?, в частности, при решении автомодельной задачи о взаимодействии двух сверхзвуковых потоков. В качестве примера на рис. 4 проиллюстрировано влияние показателя адиабаты 7 на положение линии тока, выходящей с кромки сопла. Расчеты проведены для воздуха при определяющих параметрах задачи Ма = 3, М, = 4, д=10, TjTe=\, Re0 = 10\ Числа Прандтля Рг и Шмидта Se, отвечающие как молекулярному, так и турбулентному переносу, выбраны равными 0,7.
Следует отметить, что подход, основанный на использовании упрощенных уравнений Навье — Стокса (2.1) — (2.2), записанных в декартовой системе координат, связанной с осью струи, справедлив при сравнительно небольших углах наклона разделительной линии тока к оси. Один из возможных путей преодоления указанного ограничения рассмотрен в (18], где строится криволинейная система координат и для нее выводятся упрощенные уравнения Навье — Стокса.
Остановимся на подробностях вычислительного характера. Программы для ЭВМ „БЭСМ-6“ составлены на языке „АЛГОЛ-бО“. Время счета типичных вариантов на расстояния до 40 радиусов среза сопла не превышало 30 минут. При этом получается не менее чем двух-трехкратный выигрыш по времени в сравнении с полностью явной схемой.
В заключение автор благодарит А. Н. Крайко, высказавшего ряд основных соображений по методу решения задачи, А. Н. Се-кундова, давшего исчерпывающие консультации по особенностям
модели турбулентности, М. Я. Иванова, любезно предоставившего вариант программы расчета вязких ламинарных струй но явной конечно-разностной схеме. Автор выражает признательность И. П. Смирновой, оказавшей большую помощь при проведении сопоставления с экспериментом.
ЛИТЕРАТУРА
1. А в д у е в с к и й В. С., И в а н о в А. В., Карпман И. М.,
Т р а с к о в с к и й В. Д., Юделович jM. Я. Теченне в сверхзвуковой вязкой иедорасширенной струе. .Изв. АН СССР, МЖГ*.
1970, Лй 3.
2. Авдуевский В. С., Иванов А. В., К а р п м а н И. М., Трасковский В. Д., Юделович М. Я. Структура турбулентных недорасширенных струй, вытекающих в затопленное пространство и спутный ноток. „Из'в. АН СССР, МЖГ‘, 1972, № 3.
3. Ковалев Б. Д., Мышеи ков В. И. Расчет сверхзвуковой струи, истекающей в спутный поток. .Ученые записки ЦАРИ", т. 9, .V 3, 1978.
4. Б о и д а р е в Е. Н., Г о р и н а А. Н. Решение задачи о сверхзвуковой ламинарной нерасчетной струе в спутном потоке разностным методом. „Изв. АН СССР, МЖГ*, 1968, Аг» 4.
5. Бондарев Е. Н., Гущин Г. А. Пространственное взаимодействие струй, распространяющихся в спутном сверхзвуковом потоке. „Иьв. АН СССР, МЖГ", 1972, № 6.
6. Бондарев Е. Н., Л и с и ч к о И. Д. О влиянии вязкости на течение иедорасширенной струи, распространяющейся в спутном сверхзвуковом потоке. .Изв. АН СССР, МЖГ‘, 1973, Jsfe 2.
7. Иванов М. Я., К рай ко А. Н. К численному решению задачи о нерасчетном истечении сверхзвуковой струи вязкого газа в спутный сверхзвуковой поток. Сб. .Численные методы механики сплошной среды“, т. 6, № 2, 1975.
8. Бондарев Е. Н., Л и с и ч к о И. Д. Распространение иедорасширенной турбулентной струи в спутном сверхзвуковом потоке. .Изв. АН СССР, МЖГ*, 1974, № 4.
9. Абрамович Г. Н., Крашенинникове. Ю, С е к у н-д о в А. Н. Турбулентные течения при воздействии объемных сил и неавтомодельности. М., Машиностроение, 1975.
10. А ш р а т о в Э. А., Волконская Т. Г., Росляков Г. С.,
У с к о в В. И. Исследование сверхзвуковых течений газа в струях.
Сб. .Некоторые применения метода сеток в газовой динамике", М., Изд-во МГУ, 1974, вып. 6.
11. Виноградов В. А., Иванов М. Я., Край ко А. Н.
О применении разностных схем .сквозного счета* для расчета течений с ударными волнами. Сб. „Численные методы механики сплошной среды", 1977, т. 8, № 2.
12. Г и невский А. С., И осе л е вич В. А., Колесников А. В., Л а п и н Ю. В., П и л и и е н к о В. Н., С е к у н д о в А. Н. Методы расчета турбулентного пограннчного слоя. "Итоги науки и техники. Серия механика жидкости и газа, М., ВИНИТИ, 1978, т. 11.
13. Годунов С. К., Забродин А. В., Иванов М. Я.,
Край ко А. Н„ Прокопов Г. П. Численное решение многомерных задач газовой динамики. М., Наука, 1976.
14. Jenkins R. V. Mixing and combustion of an underexpanded H2 jet in supersonic flow. A1AA Paper, N 76-610.
15. Алалыкин Г. Б., Годунов С. К., Киреева И. Л., Плинер Л. А. Решение одномерных задач газовой динамики в подвижных сетках. М., „Наука", 1970.
16. Хинце И. О. Турбулентность. М., Физматгиз, 1963.
17. Лойця некий Л. Г. Механика жидкости и газа. М., „Наука", 1973.
18. Борисов Н. Ф., Сыровой В. А. Об уравнениях вязких сверхзвуковых струй с большой степенью нерасчетности. .Изв.
АН СССР, МЖГ", 1977, Лё 2.
Рукопись поступила 14/V 1979 г.
