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

УДК 517.958:533.7

ЧИСЛЕННЫЙ АЛГОРИТМ РАСЧЕТА ДОЗВУКОВЫХ ТЕЧЕНИЙ ВЯЗКОГО СЖИМАЕМОГО ГАЗА

Т. Г. Елизарова*), М. Е. Соколова

(.кафедра математики) E-mail: [email protected]

Предложен и протестирован новый численный алгоритм расчета дозвуковых течений, основанный на системе квазигазодинамических уравнений. Особенностями алгоритма являются способ введения регуляризации и простота постановки граничных условий. Алгоритм протестирован на задачах течения в каналах с внезапным расширением и сужением.

Введение и система уравнений

В работах [1, 2] был предложен и протестирован новый численный алгоритм расчета газодинамических течений со сверхзвуковыми скоростями. Алгоритм основан на системе квазигазодинамических (КГД) уравнений, которые отличаются от уравнений Навье-Стокса (НС) дополнительными диееипатив-ными слагаемыми. Наличие этих слагаемых позволяет строить эффективные численные алгоритмы расчета течений вязкого сжимаемого газа.

В настоящей работе алгоритм [1, 2] адаптирован для расчетов дозвуковых течений газа. Предложен способ построения неотражающих условий на свободных границах.

Запишем систему КГД уравнений [1] для плоского случая:

др ^ д]тх ^ д]ту т

= о,

дх ду д(рих) d(jmxux) d(jmyux) dp

(1)

at

дх (ЯГ,

ш.

ду

ух

дх

дх ду '

д (риу) + d(jmxuy) + d(jmyUy) dp

ду ду

dt

дх Ш

ху

ш.

УУ

дх ду ' д Е + d(jmxH) d(jmyH) + dqx + dq]L ду дх ду

дх

dt

д д

(2)

(3)

(4)

Здесь р — плотность газа, их и иу — проекции вектора скорости и на оси Ох и Оу соответственно, р — давление, Т — температура, -у — показатель адиабаты, Я — газовая постоянная, Е и Н — полная энергия единицы объема и полная удельная энтальпия, которые определяются по формулам

Е = р-

р

7-1'

Н =

(Е + р) Р

О

Институт математического моделирования РАН.

Компоненты вектора плотности потока массы j вычисляются по формулам

Зтпх = Р{^х Зтпу = Р{^у

~д(ри2х) д(рихиу) др

т

Wx = -р

дх ду дх

д(рихиу) 9(pUy) др

дх

ду ду

Компоненты тензора вязких напряжений П определяются как

П-жж = Пхх + ихгох + Н, ,

^хх — 2?]

дих дх

-rjdivu,

Пху = П%? + ихт*

ПNS \\хч п

ху ух I

ди„ дио

дх ду

(5)

п ух — Пух + uywx,i

^уу — ^-уу* "I" uywy "I" R ' =

Здесь П^, ПхУ , П^,5, П^ — компоненты тензора вязких напряжений НС, а выражения для юх, ю*, К* и дивергенция вектора скорости сНуи задаются формулами

wy = t

PUx

PUx/

дих дх

диу дх

риу

риу

дих ду

диу ду

R* = т

др

' дх

др

их— + 7i>div и

Уду

др дх

др ду_

divu =

диг диь

дх ду

Компоненты вектора теплового потока q вычисляются как

Rq = тр

Lt •

;__f Р

• 1 дх \р д

PUxTT дх

a ~aNS Чу — Чу

uyRq,

д

7

-1 ду д

PUy

ду

(6)

где НС-слагаемые q, формулам

NS х

<£s =

-X

дт

дх'

и

qNS =

Чу

определяются по

дт

-УС-

ду '

Коэффициенты динамической вязкости г], теплопроводности ус и релаксационный параметр г связаны соотношениями

и>

V = Voo

т т~

-L ОТ:

ус =

7 R

1

(7- 1)Рг

V,

т =

р Sc

Г].

Здесь jjoo — известное значение г) при температуре Tqo , Рг и Sc - числа Прандтля и Шмидта соответственно, ш — заданный показатель степенной зависимости из промежутка [0.5; 1]. Система уравнений (1)-(4) дополняется начальными и краевыми условиями.

Уравнения приводятся к безразмерному виду с помощью замены:

x^x/h, y^y/h, t-ttuoc/h, p^p/poo, Ux ^ Ux/Uoo, Uy^Uy/Uoo,

p^p/(pocu2oo), E Е/(Poo^o), II -H- ll/ll'i-где h — линейный размер (высота уступа канала), Poo > и0о, и Tqo плотность, скорость и температура в невозмущенном потоке.

Для численного решения начально-краевой задачи используется явная по времени разностная схема. Разностный алгоритм и способ задания граничных условий (введение системы фиктивных ячеек) аналогичны описанным в [1, 2].

Введение регуляризации

Для стабилизации численного алгоритма при расчете дозвуковых течений в релаксационный параметр г вводится искусственная диссипация (добавка), пропорциональная шагу пространственной сетки:

где Ма = «оо/соо и Re= («ооРоо^/^оо — числа Маха и Рейнольдса соответственно, hxy = + Щ,

hx> hy шаги разностной сетки по пространству, а — численный коэффициент порядка единицы, который мы определяем подбором.

Оценочную формулу для релаксационного параметра г можно получить, предварительно оценив безразмерные значения давления и температуры как

РооТоо Роо 2

Роо = - = -Coo ~

7 7 ^ Роо ^ РоО ^ 1

7 \Ма/ ~ 7Ma2 ~ 7Ма2'

7Роо 7 1

Т =

,2'

Роо 7 Ма Ма где безразмерные плотность и скорость равны единице (роо = 1 > «оо = !)•

В результате оценочная формула для г имеет вид

1 7Ма2

г =

ahxy Ма.

Ые Эс

Формально искусственную вязкость можно считать малой, если

Ма

у Ке '

В данном алгоритме в отличие от алгоритма расчета сверхзвуковых течений [2] искусственная диссипация вводится только в КГД слагаемые. Благодаря этому искусственная диссипация отлична от нуля только в поле течения, а на границах расчетной области,включая твердые стенки, эти дополнительные слагаемые исчезают. Это можно показать, записав выражение для компоненты тензора вязких напряжений Пху (5) и тепловой поток qx (6) на вертикальной стенке. Действительно, учитывая «условие непротекания» для скорости их = 0, выражения для теплового потока и коэффициента трения на стенке (5)-(6) принимают вид

9Т П -П^-гЛ"

11ху — *-*-Ху — Ч

NS

Qx = qx =

дх'

дх

Постановка неотражающих граничных

условий и расчет течения в канале

В расчетах дозвуковых течений возникает проблема построения неотражающих граничных условий на свободных границах. Традиционно применяемые условия, основанные на римановых инвариантах для уравнений Эйлера, представляются нефизичными для расчета вязких течений и трудоемкими в численной реализации. Обзор условий на свободных границах приведен, например, в [3] и [4].

В рамках КГД алгоритма для постановки условий на свободных границах удается использовать простые и естественные граничные условия, аналогичные условиям для течений вязкой несжимаемой жидкости [5, 6]. Опишем постановку таких граничных условий на примере задачи о течении дозвукового потока вязкого сжимаемого газа в канале с внезапным расширением. Пусть газ втекает со стороны уступа (рис. 1).

Предположим, что на входе в канал течение представляет собой параболу Пуазейля [7]. Для течения

Рис. 1

вязкой несжимаемой жидкости (р = const) решение Пуазейля имеет вид

их(у) = ~ ~ у^

р

-o-i)

(7)

~Т)Р 1 + ¿Р2,

где р1 — давление в точке с координатами (Л, 0), а р2 - давление в точке (Ц), Л и / — высота и длина уступа, Н и Ь — высота и длина канала соответственно.

В случае Н = 2 и К — 1 (рис. 1), полагая, что средняя скорость во входном сечении гг^ = 1, получаем:

«Л) = -«2-„)(!-,). ! = (8)

Граничные условия на входе (8) дополняем выражениями для двух других газодинамических параметров: р= 1 и г^ = 0.

На выходной границе также по аналогии со случаем несжимаемой жидкости [5, 6] используем так называемые «мягкие» граничные условия (или условия «сноса») — равенство нулю производных для всех переменных за исключением давления, которое поддерживаем постоянным р= 1/(7 Ма2).

Эти условия будем использовать в качестве неотражающих условий на входной и выходной границах расчетной области. На твердых стенках ставятся условия «прилипания» и «непротекания» совместно с условием адиабатичности для температуры [2].

В качестве начальных условий используются параметры невозмущенного потока. Рассматривается

Таблица 1

Расчет течения в канале с внезапным расширением

Ее Ма а Сетка Нх — Ьу ей Невязка Число шагов Время ь3

100 0.1 0 100 х 20 0.1 10-4 2 • М-2 201500 20 4.7

100 0.1 0.01 100 х 20 0.1 10-4 4. Ю-4 1 201500 111 5.0

100 0.1 0.1 100 х 20 0.1 М-3 10-4 86145 86 5.1

100 0.1 0.2 100 х 20 0.1 М-3 10-4 80 333 80 5.0

100 0.1 0.5 100 х 20 0.1 М-3 10-4 80 000 80 5.1

100 0.1 1.0 100 х 20 0.1 Ю-3 10-4 55133 55 5.1

100 0.1 1.5 100 х 20 0.1 м-3 10-4 45 413 45 5.0

100 0.1 2.0 100 х 20 0.1 5-Ю"4 10-4 75 405 38 5.1

100 0.1 5.0 100 х 20 0.1 10-4 10-4 219 849 22 4.8

100 0.1 0.5 100 х 20 0.1 М-3 10-4 80 000 80 5.1

100 0.1 0.5 200 х 40 0.05 10"3 10-4 80 832 81 5.0

200 0.1 0.5 140 х 20 0.1 М-3 10"» 170 623 171 8.50

200 0.1 0.5 280 х 40 0.05 М-3 10"» 203 450 204 8.35

300 0.1 0.5 160 х 20 0.1 М-3 10-4 238 033 238 10.10

300 0.1 0.5 320 х 40 0.05 10"3 10-4 140 000 140 9.95

400 0.1 0.5 200 х 20 0.1 М-3 10-4 535 000 535 11.40

400 0.1 0.5 400 х 40 0.05 М-3 10-4 626 000 626 12.70

100 0.5 0.5 100 х 20 0.1 м-3 10-4 95142 95 4.9

100 0.2 0.5 100 х 20 0.1 м-3 10-4 103 775 104 5.0

100 0.1 0.5 100 х 20 0.1 10"3 10-4 80 000 80 5.1

100 0.05 0.5 100 х 20 0.1 м-3 10-4 39 774 40 5.0

100 0.01 0.5 100 х 20 0.1 10-4 5-Ю"6 271 066 27 4.9

течение воздушного потока при нормальном давлении. Значения молекулярных параметров для воздуха: 7= 1.4, Рг = 0.737, 8с = 0.737, и; = 0.74.

Расчеты проведены на равномерных пространственных сетках для чисел Рейнольдса Ые = 100, 200, 300, 400 и чисел Маха Ма = 0.5, 0.2, 0.1, 0.05, 0.01. Стационарное решение находилось методом установления. Результаты расчетов течения в канале с внезапным расширением приведены в табл. 1. На рис. 1 изображен процесс установления течения для варианта Ые = 300, а = 0.5, Ма = 0.1. Построены распределения плотности и приведены линии тока для времен ¿= 1, 8, 15, 22, 238. Видно, что возмущения свободно пересекают выходную границу области.

Для рассмотренных течений градиенты плотности пропорциональны 1/Ма2 [7], что позволяет проводить сопоставление решения с расчетами, выполненными в приближении вязкой несжимаемой жидкости [5, 6]. Полученная в расчетах длина отрывной зоны Ья для всех вариантов хорошо соответствует эталонным результатам. При измельчении шага пространственной сетки длина отрывной зоны становится ближе к эталонному решению (табл. 2).

Длина отрывной зоны практически не зависит от числа Маха, когда 0.01 < Ма < 0.5. При дальнейшем уменьшении (Ма < 0.01) значительно замедляется

Таблица 2

Длина отрывной зоны Ь8 при разных Не

Жидкость Газ

Ее Нх = 0.05 Нх =0.1 Нх = 0.05

100 5.0 5.1 5.00

200 8.2 8.5 8.35

300 10.1 10.1 9.95

400 14.8 11.3 12.70

скорость сходимости, поскольку время счета обратно пропорционально числу Маха: £ ~ 1/ Ма.

Оптимальный диапазон значений для регуляри-зирующего параметра составляет 0.5 < а < 2. Этот диапазон соответствует наилучшей точности решения и минимальному числу шагов до сходимости (табл. 1).

В качестве второго примера рассматривалось течение воздушного потока в канале с внезапным сужением. По сравнению с предыдущим случаем изменяется направление течения — газ втекает справа налево, размер области и уступа остаются прежними. Проведены расчеты Ые = 100 и йе = 400 при Ма = 0.1 на двух пространственных сетках.

При Ые = 100, а = 0 устойчивое решение получить не удается. При а < 0.2 и а > 0.5 сходимость и устойчивость алгоритма ухудшается. В диапазоне

Ее = 400, а = 0.5

0 2 4

х

Рис. 2

0.1 < а < 1.0 перед передней стенкой цилиндра образуется небольшой вихрь. При числе Рейнольдса Ые = 400 образуются два вихря: на торце и на боковой стенке уступа. На рис. 2 приведено распределение плотности и линий тока для варианта Ые = 400, а = 0.5.

При а = 0.5 полученные результаты расчетов для Ые = 100 и Ые = 400 хорошо согласуются с данными [8], где приведен двумерный расчет ламинарного течения вязкой несжимаемой жидкости в каналах различной формы в ж, у-геометрии. Расчеты выполнены в переменных «функция тока» — «вихрь скорости». На входе задан профиль Пуазейля.

Для Ые = 100 оптимальный диапазон составляет 0.1 < а < 1, а для Ые = 400 область значений параметра а сужается до 0.5 < а < 1. Значения а > 1 не исследовались.

Выводы

В работе построен и протестирован эффективный численный метод расчета дозвуковых течений, основанный на КГД уравнениях. Основными особенностями метода являются способ введения искусственной диссипации, который позволяет рассчитывать тепловые потоки и коэффициент трения без искажений, а также естественный способ задания неотражающих граничных условий на свободных границах, который позволяет избежать применения условий, основанных на характеристиках для уравнений Эйлера.

Для рассмотренных задач параметр регуляризации целесообразно выбирать в диапазоне 0.5 <а< 1.5.

С уменьшением числа Маха скорость сходимости построенного метода замедляется, что является естественным при использовании полных уравнений газовой динамики для расчета дозвуковых течений. Для преодоления этого недостатка существуют специальные методы, например, LM-приближение [9].

Работа выполнена при финансовой поддержке РФФИ (грант 01-01-00061). Авторы благодарны Ю. В. Шеретову за сотрудничество.

Литература

1. Шеретов Ю.В. // Применение функционального анализа в теории приближений. Тверь, 2001. С. 191.

2. Елизарова Т.Г., Соколова М.Е. // Вестн. Моск. ун-та. Физ. Астрон. 2004. №1. С. 10 (Moscow University Phys. Bull. 2004. N 1. P. 12).

3. Ильгамов M.A., Гильманов A.H. Неотражающие условия на границах расчетной области. М., 2003.

4. Дородницын Л.В. // ЖВМ и МФ. 2002. 42, №4. С. 522.

5. Елизарова Т.Г., Калачинская И.С., Шеретов Ю.В., Шиль-ников Е.В. // Сб. Прикладная тематика и информатика. Тр. факультета ВМиК. М„ 2003. №14. С. 85.

6. Armaly B.F., Li A., Nie J.H. // J. of Thermophysics and Heat Transfer. 2002. 16, N 2. P. 222.

7. Лойцянский Л.Г. Механика жидкости и газа. М., 1987.

8. R.W. Меи, Л Plotkin. // AIAA Paper. 1986. N 86. P. 110.

9. Лапин Ю.В., Стрелец М.Х. Внутреннее течение газовых смесей. М., 1989.

Поступила в редакцию 13.04.04