CC BY
39
Вычислительные технологии
Том 12, Специальный выпуск 4, 2007
МОДЕЛИРОВАНИЕ СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ ГАЗА В КАНАЛЕ
В.М. Ковеня, А.Ю. Слюняев Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: kovenyaOict. nsc. ru, viimyOngs. ru
On the basis of the splitting method on physical processes and spatial directions, supersonic flows of gas in a channel with an injection of gas from a part of its surface, are numerically investigated. The influence of the Mach and Reynolds numbers on flow characteristics is investigated.
Введение
Интенсивное развитие вычислительных систем и повышение их быстродействия позволяют проводить исследования все более сложных физических процессов. Однако стремление к более полному исследованию изучаемых явлений и соответственно использование наиболее полных математических моделей приводят к возрастанию вычислительных затрат. Таким примером являются задачи аэродинамики, исследование которых на начальных этапах проводилось с помощью системы уравнений Эйлера и других приближений, а затем с развитием вычислительной техники — с использованием полной системы уравнений Навье — Стокса. Численное решение уравнении Навье Стокса и сегодня представляет большие трудности, что обусловлено нелинейностью исходных уравнений, наличием областей больших градиентов и других особенностей, возникающих при определенных параметрах и режимах гидродинамических течений. Как следствие, это вызывает необходимость разработки и создания специальных численных методов решения этих уравнений. Хотя к настоящему времени разработано много численных алгоритмов и специальных комплексов программ (см, например, [1-8]), проблема создания и применения эффективных численных методов остается открытой.
В настоящей работе используется неявная разностная схема, предложенная ранее авторами [9], основанная на специальном расщеплении исходных многомерных операторов по физическим процессам и пространственным направлениям. Расщепление исходного многомерного оператора по направлениям позволяет свести решение исходной задачи к последовательности одномерных аналогов. Дальнейшее расщепление каждого одномерного оператора по физическим процессам сводит реализацию схемы на дробных шагах к решению одномерных скалярных уравнений.
В настоящей работе дано краткое изложение алгоритма решения уравнений Навье — Стокса сжимаемого теплопроводного газа и представлены некоторые результаты расчетов плоских течений вязкого газа в канале. Расчеты проведены в предположении, что
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
с части канала вдувается газ с трансзвуковой скоростью. Исследовано влияние чисел Маха и Рейнольдса на картину течения, получены основные закономерности течения, зависящие от параметров набегающего потока.
1. Исходные уравнения и численный алгоритм
Для численного решения уравнений Навье — Стокса сжимаемого теплопроводного газа в канале в качестве базовой выбрана декартова система координат х1; х2, Однако при проведении расчетов ее использование не всегда удобно. Поэтому в расчетной односвяз-ной области О (рис, 1), в которой будем отыскивать численное решение задачи, введем невырожденное преобразование координат:
д. = qj (х*), где;, г = 1, 2, (1)
для которого существует его обратное преобразование:
х. = х.)•
Тогда система уравнений Навье — Стокса в новых переменных qj• может быть записана в консервативной форме (см. [9, 10]):
Л
д и
-= Л¥, Л¥ = - .
д1 ^ %
(2)
Здесь
и = -
(9 \
ру 1 РУ 2
Vе /
3 = det
( ддг дд2 \
дх1 дх1
дд\ дд2
\ дх2 дх2 )
£
= (А^2 - ql2q1
— соответственно вектор искомых функций и якобиан преобразования координат, а
( ри.и
Wj
1 1
\
ру 1 и. + q71 р - С1
ру2 и. + qj2р - С? V (Е + р)и. + Qj - Я. )
— векторы потоков в направлении qj•, где соответствующие коэффициенты равны
2
1=1
2 ,дТ 2
т= 1
2 / ду' ду' а) = 5}Л<Иуу + /X Е ^ + д] ^
1=1 т=1
2 2 ду
Л = ^ - 2/(3/х).
.=1 г=1 Щ1
Выбирая в качестве искомых функций вектор f = И, представим систему уравнений (2) в недивергентном виде:
2
Матричные операторы Б^- включают газодинамические члены уравнений и часть вязких членов, содержащих лишь повторные производные по каждому направлению qj•, а вектор Е содержит оставшиеся таены уравнений. Выберем f = (р, у1,у2,р)т. Тогда матричные операторы для данного вектора f примут вид
Бj
где
дЯо 3 0
0
0
о
д
0
0 0
0
\
д
^ дqj
■Е
т=1
с2 — j дqj
а1 —
j дqj
а2 — j дqj
з —1,1 о д
I = 1, 2,
2
га - 1 — э --
^з - 2^
т=1
4 = qj•/рп, е{ = •
В преобразованной расчетной области (0 < qj < 1) введем равномерную разностную сетку по каждому пространственному направлению qj• и аппроксимируем исходные операторы Б^^ ^^^^^^^^^^ операторами Б^ с порядком 0(Нк), Введем расщепление операторов Б^ на сумму операторов Б^ в виде
(4)
Бjh
где
Б
jh
0 0 0 0 \
0 0 0 а] А,
0 0 0 а2Л,-
0 с1 Л • cj Лj с2Л • cj Лj ипЛ Uj Лj _ /
Б
jh
( Лj ип 0
0 0
UjnЛj - С]н
0
0
0 0
0
— С 2
jh
0 0
0 0
Систему уравнений (2) или (3) аппроксимируем подобно [9] разностной схемой:
2 £П+1 _ £П
п
¿=1
(I + га В}„)(1 + таВ%)-= -(А"1)^.
(5)
Схеме (5) эквивалента схема в дробных шагах:
Сп = — (А-1)^, (I + гаБ1ь)Сп+1/4 = Сп,
(I + таБ2ь)Сп+1/2 = Сп+1/4,
(I + гаБ2ь)Сп+3/4 = Сп+1/2, (I + гаБ2ь)Сп+1 = Сп+3/4, £п+1 = £п + тС'п+1
1
р
0
1
2
Она аппроксимирует нестационарные уравнения Навье — Стокса с порядком 0(т + h2), а при установлении уравнения (2) в дивергентной форме — с порядком 0(h2) и в линейном приближении безусловно устойчива при а > 0.5, Реализация схемы (6) па дробных шагах сводится к трехточечным скалярным прогонкам [9]. При построении разностной схемы дифференциальные операторы d/dxj в оператоpax Bj аппроксимированы разностными операторами Aj или Aj с порядком 0(hk) по формулам
если Vj < 0, то Aj = Aj = Aj_,
если Vj > 0, то Aj = Aj_, Aj = Aj+
при k = 1, и по формулам
для k = 2, Здесь Лj± = ^(1 — T±1)/hj, T±fl = fl±1 — оператор сдвига. Компоненты вектора Wh в схемах (5) или (6) аппроксимированы центральными разностями. Как известно, симметричная аппроксимация первых производных приводит к осцилляциям в численном решении. Для их подавления вводится сглаживающий оператор по каждому направлению qj второго порядка малости подобно [9]
* гл Ql+1 — Ql-1 ■ ! \ 2 Ql+1 — 2Ql + Ql-1 jQl =-2h---а gn 1-2h~i-'
( \Qi+i -2Ql+Ql_l\ e\=\ IQi+i-Qil + lQi-Qi-il' (7)
[ о^лп |Ql+1 — Ql| + |Ql — Ql-11 = 0.
Отметим, что введение сглаживания по каждому направлению вида (7) приводит к некоторому понижению устойчивости схемы,
2. Численные расчеты
Апробирование и тестирование описанного выше алгоритма проведены в работе [9], Ниже приведены результаты численных расчетов течений вязкого газа в канале при наличии вдува газа с части нижней стенки при различных числах Маха и Рейнольдса, Расчетная область (рис, 1), в которой находилось решение задачи, представляет собой прямоугольник R G {—0.5 < x1 < 6,0 < x2 < 2}, а внутри него находится канал. Нижняя стенка канала отнесена от границы расчетной области на расстояние 0,5 и
x1 = 0.0
x1 = 0.5
Стационарное решение задачи находилось в приближении уравнений Навье — Стокса вязкого сжимаемого теплопроводного газа в преобразованных координатах (2) по разностной схеме (6), Преобразование координат задавалось в виде
x1 — a' в
qi = 1 I ' 92
V1(x1) — 1
У2(x1) — V1(x1)
Здесь у1(х1^ у2(х1) — уравнения нижней и верхней стенок канала; 9 ж в — коэффициенты сгущения узлов сетки к границам расчетной области; а и Ь — левая и правая границы расчетной области.
в
% Г3
Рис. 1. Расчетная область
На границе области при х1 = -0.5 задавался невозмущенный сверхзвуковой поток, параллельный оси х1; на линиях х2 = 0, -0.5 < XI < 0.0 и х2 = 2, -0.5 < XI < 0.5 задавались условия симметрии, па стенках капала — условия прилипания и тепловой изоляции Vj = 0 , дТ/дх2 = 0, а на выходе из канала — мягкие условия дт¥/дхт = 0, где т = 1 или 2. На нижней стенке канала задавался источник подачи газа, начинающийся в точке х = 1.0 шириной 0.05 единицы. Скорость вдува газа полагалась звуковой, т. е. МЕД = 1.0, скорость набегающего потока варьировалась в пределах М =1.5... 2.5, числа Ле = 103 ... 104, 7 = 1.4 Рг = 0.72, и = 0.76.
Для удобства численного интегрирования и анализа результатов исходные уравнения Навье — Стокса приведены к безразмерному виду, где в качестве безразмерных параметров выбраны следующие величины: высота половины канала Ь, скорость набегающего потока щ, плотность газа р0, температура Т0. Тогда уравнения (2) или (3) в криволинейной системе координат в безразмерном виде сохраняют ту же форму, а при вторых производных в уравнениях движения и энергии появляются безразмерные параметры — числа Рейнольдса и Прандтля [8|, Стационарное решение задачи находилось методом установления, критерием установления полагалось выполнение условия (рп+1 — рп)/(Грп) < 10-4 во всех узлах расчетной области. Расчеты проводились на сетке, содержащей 521 х 161 узлов. Установление к стационарному решению достигалось за ~ 5500 итераций при использовании в качестве начального распределения несогла-
х1 = -0.5 0 < х2 < 2 параметры набегающего потока, а в расчетной области задавались нулевая скорость и постоянные значения плотности и температуры). При задании в качестве начального приближения решения задачи при других параметрах (например, при другом числе Рейнольдса или числе Маха) число итераций значительно (в несколько раз) сокращалось.
2.1. Расчеты с различными числами Маха
В первой части работы проведена серия расчетов сверхзвукового течения газа в канале при различных числах Маха и фиксированном числе Яе = 104, На рис. 2 и 3 представлены распределения поля скорости и изолинии давления для чисел Маха 1.5 и 2.5. Хорошо видно, что изменение числа Маха в набегающем потоке существенно влияет
на картину течения в канале. Как следует из теории ударных волн, наклон головного скачка уплотнения для случая М = 1.5 (см. рис. 2, а и 3, а) существенно меньше, чем для М = 2.5 (рис. 2, б и 3, б). В первом случае головной скачок уплотнения, возникший около источника вдува, падает на верхнюю стенку и отражается от нее (точнее, от пограничного слоя), затем отражается от нижней стенки канала. Во втором случае скачок уплотнения вызывает отрыв пограничного слоя, возникают два скачка: один инициирован отрывом пограничного слоя, а второй — отражением от пограничного слоя. Далее они также взаимодействуют с потоком в канале и, поскольку они имеют значительно меньший угол наклона, для рассматриваемой расчетной области не успевают взаимодействовать с нижней стенкой.
От числа Маха набегающего потока существенно зависят и размеры областей отрыва пограничного слоя на нижней и верхней стенках канала. При малых числах Маха интенсивность головного скачка недостаточна для отрыва пограничного слоя у верхней стенки и отрыв возникает лишь при М> 1.5. На рис. 4 представлено распределение коэффициентов поверхностного трения для нижней и верхней стенок канала, это подтверждающие. Как и ожидалось, наименьшие размеры отрывной зоны в окрестности вдува наблюдаются для меньшего числа Маха набегающего потока М = 1.5, наибольшие — для М = 2.5 (рис. 4, а).
Рис. 2. Поле скорости для М = 1.5 (а) и М = 2.5 (б)
Рис. 3. Изолинии давления для М = 1.5 (а) и М = 2.5 (б)
Рис. 4. Коэффициент трения на нижней (а) и верхней (б) стенках
На рис. 5 представлены распределения полей скорости и области отрыва пограничного слоя около источника вдува для двух значений чисел Маха. Видно, что размеры отрывных зон значительно различаются, что связано с интенсивностью потока газа после прохождения головных скачков. При оценке размеров отрывной зоны для верхней стенки канала можно сказать, что для М = 1.5 интенсивности ударной волны не хватает для отрыва пограничного слоя (см. рис. 4, б, линия 1). Начиная с М = 1.75 и для всех последующих значений наблюдается отрыв пограничного слоя. Наибольшая ширина отрывной зоны, как и для нижней стенки, наблюдается для максимального числа Маха М = 2.5.
На рис. 6 представлено распределение чисел Маха в продольном и поперечном сечениях канала. Отметим, что поток газа после прохождения через головной ударный скачок на уровне У = 0.5 тормозится до примерно одинаковой скорости независимо от скорости набегающего потока (различие скоростей потоков составляет не более 6 %). В
х = 1.5
Распределение давления на нижней и верхней стенках канала приведено на рис. 7. Кривые на рис. 7, б подтверждает, что изменение скорости набегающего потока (числа Маха) приводит к изменению угла наклона головного скачка. При увеличении скорости
Рис. 5. Линии тока у источника вдува при М = 1.5 (а) иМ = 2.5 (б)
Рис. 6. Распределение чисел Маха в сечениях У = 0.5 (а) и X = 1.5 (б)
Рис. 7. Распределение давления на нижней (а) и верхней (б) стенках
набегающего потока за источником вдува происходит сильное разряжение газа (см. рис. 7, а). Это связано с ускорением течения газа над отрывной зоной и соответственно ускорением истечения газа из отрывной зоны.
2.2. Течения с различными числами Рейнольдса
В следующей серии расчетов проведено моделирование течения газа в канале воздухозаборника со вдувом газа с части поверхности для различных чисел Рейнольдса набегающего потока при фиксированном числе М = 2.0. Число Ие последовательно задавалось равным 103, 5 • 103 и 104.
На рис. 8 представлены поля скоростей и линии равного уровня давления для всех исследуемых чисел Рейнольдса. Видно, что с уменьшением числа Рейнольдса толщина пограничного слоя на стенках возрастает и от носков нижней и верхней пластин возникают все более сильные скачки уплотнения. В окрестности нижней стенки первый головной скачок и пограничный слой взаимодействуют с головным скачком, возникшим в области вдува. За отрывной зоной в области вдува возникает висячий скачок, который взаимодействует с отраженными скачками. Дальнейшее их взаимодействие приводит к сложной картине течения, которая, конечно, зависит от чисел Рейнольдса. Отрывные зоны у верней стенки канала, возникшие после падения головного скачка, приводят к перестройке течения и возникновению вторичных отраженных скачков.
Видно (см. рис. 8, в), что с уменьшением чисел Рейнольдса толщина пограничного слоя и область с дозвуковыми скоростями возрастают и скачки уплотнения отража-
а (Ие = 104)
в (Ие = 103)
Рис. 8. Поля скоростей (слева) и изолинии давления (справа)
Рис. 9. Распределение ж-компоненты скорости в сечении Y = 0.0125
ютея уже не от стенок канала, а от пограничного слоя. Повышение числа Рейпольдса также приводит к уменьшению зоны отрыва па верхней стойко капала и увеличению зоны отрыва пограничного слоя па нижней стопке капала. Об этом факте также можно судить по рис. 9. Скорость течения (линия 1) значительно падает и после прохождения через скачок уплотнения сохраняет только примерно 5 % от своей первоначальной величины. Это и говорит о том, что рассматриваемая область находится в пограничном с л оо. Совеем другую картину можно наблюдать для максимального числа Рейпольдса, Здесь после прохождения скачка уплотнения и присоединения пограничного слоя поток быстро набирает свою первоначальную скорость.
Выводы
Для расчета течений газа в канале типа воздухозаборник с вдувом газа с части ого поверхности применена экономичная схема расщепления (но физическим процессам и пространственным направлениям) второго порядка аппроксимации. Простота реализации схемы на дробных шагах, возможность варьирования шагов сетки и параметров потока в широком диапазоне (на основании большого запаса устойчивости схемы), экономичность алгоритма позволяют получать решения сложных задач при небольших затратах ресурсов ЭВМ,
Проведено моделирование вязких сверхзвуковых течений в канале со вдувом газа с части ого поверхности при различных числах Маха и Рейпольдса набегающего потока, Исследовано влияние параметров потока на картину течений в канале и получены основные закономерности их изменения.
Список литературы
[1] Флетчер К. Численные методы в динамике жидкостей. В 2-х т. М.: Мир, 1991.
[2] Численное решение многомерных задач газовой динамики / С.К. Годунов, A.B. Забродин, М.Я. Иванов и др. М.: Наука, 1989.
[3] Ковеня В.М., Яненко H.H. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, Сиб. отд-ние, 1981.
[4] Computer and Fluids, Special Issue // 6th Intern. Symp. on CFD. 1998. Vol. 27, N 5 6.
[5] Computational Fluid Dynamics J. // Special Issue dedicated to prof. H. Daiguji. 1999. Vol. 8, N 2.
[6] Computational Fluid Dynamics J. // 4th Asian Workshop on CFD. 2004. Vol. 13, N 2.
[7] Computer and Fluids // Special Issue dedicated to prof. R. Peuret. 2002. Vol. 31, N 4-7.
[8] Лойцянекий Л.Г. Механика жидкости и газа. М.: Наука, 1970.
[9] Ковеня В.М. Разностные методы решения многомерных задач: Курс лекций. Новосибирск: Изд-во НГУ, 2004.
[10] Ковеня В.М., Слюняев А.Ю. Модификации алгоритмов расщепления для решения уравнений газовой динамики и Навье — Стокса // Вычисл. технологии. 2007. Т. 3, 12.
Поступила в редакцию 1 декабря 2007 г.
