Моделирование сверхзвуковых течений газа в канале Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 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 г.