Применение ENO-схемы (essentially nonoscillatory) для моделирования течения многокомпонентной газовой смеси Текст научной статьи по специальности «Математика»

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

Том 12, Специальный выпуск 4, 2007

ПРИМЕНЕНИЕ ЕХО-СХЕМЫ (ESSENTIALLY NONOSCILLATORY) ДЛЯ МОДЕЛИРОВАНИЯ ТЕЧЕНИЯ МНОГОКОМПОНЕНТНОЙ ГАЗОВОЙ СМЕСИ

А. О. Бекетаева, А. Ж. Найманова Институт математики Министерства образования и науки,

Алматы, Казахстан e-mail: ked@math.kz

This paper addresses a numerical investigation of the plane supersonic turbulent multi-component flow with the perpendicular injection of a hydrogen from the slots located symmetrically on the lower and the upper walls of a channel. For the mathematical modeling of such flows the full Navier — Stokes equation for multi-component flow fields were used. On the basis of the ENO-scheme the computational code has been developed. We study the influence of the Mach number on the length of separation region and the penetration depth of a hydrogen jet into a supersonic air crossflow.

Введение

Обтекание струй сверхзвуковым потоком исследовано в основном экспериментально, при этом многие закономерности, такие как образование системы скачков уплотнения, отрывных зон в области вдува струи, изучены качественно [1]. В теоретических работах по моделированию такого рода течений в основном исследуется одноатомный газ и практически отсутствуют исследования влияния на схему течения таких важных с практической точки зрения параметров, как число Маха, степень нерасчетно-сти, ширина щели. Так, например, в работе [2] рассмотрен поперечный вдув струи в сверхзвуковой поток, изучено влияние различных моделей турбулентности на характеристики течения. В [3] для аналогичной задачи используется адаптированная сетка. В работах по исследованию многокомпонентного газа (например, в [4, 5|) для решения уравнений Навье — Стокса применены традиционные разностные методы, где для обеспечения устойчивости используется искусственная вязкость. На сегодняшний день для решения уравнений Навье — Стокса применяются квазимонотонные консервативные схемы повышенного порядка аппроксимации без введения искусственных диссипатив-ных членов, такие как TVD-схемы (Total Variation Diminution Schemes) [2]. При этом их основной недостаток состоит в том, что в окрестности разрыва решения порядок точности понижается до первого, в результате скачки уплотнения, возникающие в течении,

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

Lh и»

Рис. 1. Схема точения

могут быть сильно размазаны. Один из способов решения этой проблемы — использование копечпо-разпоетпых схем высокого порядка точности, таких как ЕХО-ехемы (Essentially Xonoseillatory Schemes) |6-8|.

В настоящей работе на основе разработанной методики решения уравнений Навье — Стокса для течения многокомпонентной газовой смеси численно моделируется плоское сверхзвуковое турбулентное течение воздуха при наличии поперечного вдува водорода со стенок канала. При исследовании для удобства вычисления рассматривается вдув струи с нижней стенки. Схема течения показана на рис. 1. Изучается влияние числа Маха потока на длину отрывной зоны и дальнобойность вдуваемой водородной струи.

1. Постановка задачи

Исходной является система двумерных оередненных но Рейнольдсу уравнений Навье — Стокса для многокомпонентного газа, записанная в безразмерной форме в общепринятых обозначениях:

9и + а(Е-Е.) + а(Е-Р.)=01 (1)

■'v )

dt ' dx ' dz

\T

U = (p, pu, pw, Et, pYk) E = (pu, pu2 + p, puw, (Et + p) u, puYk)T F = (pw, puw, pw2 + p, (Et + p) w, pwYk)

T

Ev (0, Txx, Txz, uTxx + uTxz OX? Jkx) ,

Fv (0, Txz, Tzz, wTxz + wTzz qz, Jkz)

T

N i T

Et =-~j\f2~ X! Ykhk - P+oP ("2 + w2) 1 hk = h°k + / c.pkdT, cpk = Cpk/W. (3)

7°M<o k=1 2 J

To

Здесь молярные удельные теплоемкости Срк определяются по экспериментальным данным при помощи полиномиальной интерполяции четвертого порядка по температуре; численные значения эмпирических констант а, к приведены в таблице Л А КАР [9]; Ук — массовая концентрация к-й компоненты, к = 1,..., N оде N — число компонентов смеси газов; тхх, тгг, тхх = тхх — тензоры вязких напряжений (для определения турбулентной вязкости используется модель Болдуина — Ломакса); дх, и ]хк, .]гк — тепловые и диффузионные потоки (диффузионные потоки вычисляются по закону Фика [10]).

В качестве определяющих параметров при обезразмеривании приняты параметры па входе (ите, рте, Тте), давление и полная энергия отнесены к значению ртеи2те, удельная энтальпия отнесена к Я0Тте/Ште, молярные удельные теплоемкости — к Л0, характерным параметром длины является ширина щели.

На входе задаются следующие параметры потока:

= Р = Рте, Т = и = Мо,

7те>Л0 Тте

Ж* '

т = 0, Ук = Укте, х = 0, 0 < г < Н;

во входном сечении вблизи стенки задается пограничный слой, продольная составляю-и

На щели предполагается:

№к = Шк о, р = проа, Т = Т0, ъи = Мо/У°,

и = 0, Ук = Ук0, г = 0, г = 0, Ьъ < х < Ьъ + К.

Здесь п = р0/рте — степень нерасчетпости, р0 — давление в струе, рте — давление в потоке. На нижней стенке задается условие прилипания и теплоизоляции, на верхней гра-нпце - условие симметрии; па выходной границе задаются условия неотражения [11].

2. Метод решения

Система уравнений (1) решается численно, при этом аппроксимация конвективных слагаемых производится на основе ЕМО-схемы третьего порядка точности. Способ построения этой схемы будет показан на примере одномерного невязкого случая, который обобщается на двумерный случай. Для этого рассматривается следующая задача Коши:

д и д/ (и) ди . ди м . .

Ж + = ж + »(«.») = "<>. м

д / (и)

где А = —---матрица Якоби.

ди

К

начальные данные функции и(х, заменяются кусочно-постоянными, а именно его

интервалов возникает задача о распаде произвольного разрыва для малого отрезка по времени Д£ = (£п+1 — ¿п). Шаг по времени должен выбираться настолько малым, чтобы решения, полученные для каждых двух соседних стыков, не взаимодействовали. Чтобы

совершить следующий шаг по времени, необходимо снова иметь кусочно-постоянные функции. Таким образом, при найденных V0 = и0(ж) рассматриваемая задача Коши в полосе < £ < Ьп+\ для Л > 0 примет вид

дш дш

+ = 0, ю(х,гп) = Щх-,ь)™)

(5)

тогда

ш(ж, £) = — Л£; ш"),

х1 + 1/2

1

(6)

^'-1/2

где w = R-1vh — инвариант Римана.

Из (6) следует, что необходимо знание функции R(•; ш"). Согласно методу ЕР [6] примем ее в виде

й

Д(ж;гу) = —Нт(х; Ш), аж

где Нт(ж; Ш) есть кусочно-полиномиальная функция степени т; Ш(ж) — первообразная для функции ш(ж).

Тогда (6) запишется как

-га+1 _

1

гЗ + 1/2

1

^ =1

Я (х- Ь))<1х = - (Нт (ху+1/2 - А£; Ш) - Нт (х^1/2 - Лí; Ш)) . (7)

Х1-1/2

В качестве полинома Нт(ж; Ш) принимается формула Ньютона третьей степени, алгоритм построения которого приведен в работе [8].

Л > 0 Л < 0

+

"+1 = ш" _ а+Д_ш" - а+Д_

г«/ = щ

1- а+

(minmod (Д_ш", Д+ш")) +

-—За+ + (т (А_А_ц)?; А-А+гВ?)) если |Д_<| < |Л+го

~ 1

6

(т (Д_ Д+ш"; Д+Д+ш")) если |Д_ш"| > | Д+г«

Л

—а_Д+ш" + а_Д+ (а_)2 — 1 „

1 +

(minmod (Д_ш", Д+г«")) —

2 + 3а_ + (а_)2 ( .

6

(т. (Д_Д_ш"; Д_Д+ш")) если | Д_г«n| < |Д+г«

о I

6

(т(Д_Д+г«"; Д+Д+ш")) если |Д_ш"| > | Д+г«

где а± =

А г к 5

Л

±

Л ±|Л|

Д + Wj = — Д — = —

. , ,л , а, если а < О, т(а, о) = < , ,

I о, если |а| > |о| ;

(8)

2

штшо^а, Ь)

в шт(|а| , |Ь|), если sgn (а) = sgn (Ь) = 5, 0, если иначе.

После осуществления перехода от переменных -шп — инвариантов Римапа — к переменным уп схема (8) будет иметь вид

п+1 _ 0

^ ' V/ -.К

1-1,2*-*? - %А7+1/2= о,

к

(9)

где — модифицированный поток на узловых точках состоящий из исходного конвективного вектора fj и добавочных членов высокого порядка точности, а именно:

е т 10

1 0 + Е, + Б,

при этом векторы Е, и Б, определяются следующим образом:

Е, = шinшod (-Ё,-+1/2, ,

Б ( т(Д_О_1/2, Д+О0-1/2), 0 I М.Д-оО,+1/2, Д+оО,+1/2),

Б,

т(Д_О0+1/2, Д+О0+1/2), ш(Д_Д-1/2, Д+ОО0-1/2),

если |Д_V,1 < |Д+V,1

если |Д_vj I > |Д+vj |

если |A_Vj| < |Д+V, |

если 1 > |A+Vj |

для А > 0, для А < 0,

где

Щ+1/2 = ^П (Л-+1/2) |Л'+1/2|) А + ^'/2 А-+1/2 = 8ёП (Л+1/2) (^|Л+1/2|) "

Д + 10 /6,

О0+1/2 = ( А,+1/2)

т | | /Т I

Д + ^/6,

+ Л = /, А1*1 = //А : // 1 // (-^—- ) // 1: / — единичная матрица, матрицы А1*1

рассматриваются как А± = /А.

Из (9) следует, что величина модифицированного потока 1т совпадает с выражением потока, полученным в работе [12] на основе построения решений вдоль характеристических направлений. Аналогично этой работе уравнения (9) можно формально рассматривать как одностороннюю схему с разностями против потока первого порядка аппроксимации для следующего уравнения:

^ - д 1т

0.

(10)

Для более точного учета течения в областях больших градиентов, т. е. в пограничном слое, вблизи стенки и на уровне щели, вводится сгущение сетки с помощью простых преобразований вида

£ (х) = К + — агвЬ т

X

1 ^ вЬ(тК)

П(г) = Н

(в +1) - (в - 1)

/3 + 1 /3-1

а

/3+1 /3-1

1-2-

а

+1

2

К = — 1п 2т

где в т _ коэффициенты сгущения, в, т > 1; а — высота расчетной области в новой системе координат; жс — точка, относительно которой производится сгущение. Тогда система уравнений (1) в координатах (С, п) примет вид

дУ | дЕ | дР = | | |

д£ дС дп дС дС дп дп

(П)

где и = и//, Е = СхЕ// / = п^Г// = СхЕ„2// = СхЕ„т// = п^Г„2// = п*// / = д (С, п) /д (ж, г) — якобиан преобразования; Е„ = + Е„т и Р = + Р„т — диффузионные члены, содержащие вторые (£„2, и смешанные

производные (Е„т,

)

дЦ_

+ А+ +

+ (в + + В -

дп

д(Е„2 + Е„

дС

г

дп

0. (12)

Здесь модифицированные векторы потоков имеют вид Ет = Е + Е^ + и ¥т =

Б + Еп +

Одношаговая конечно-разностная схема для интегрирования по времени системы (12) будет иметь следующий вид:

Ди"+1 + Д£

V / \ ) дг]

д (Е"2+1 + Е"т) д (/Р^1 + Щ

дС

дп

0.

(13)

Численное решение осуществляется согласно принципу расщепления методом матричной прогонки [13]. На каждом шаге по времени в системе уравнений (13) конвективные слагаемые линеаризуются с помощью замены их значений с (п + 1)-го слоя разложениями в ряд Тейлора с известными значениями на предыдущем слое с номером п: Е"+1 = Сх = пх оде В = дГ/ди _ матрица Якоби. После

этого они записываются как односторонние схемы с разностями против потока первого порядка аппроксимации аналогично (9). Диффузионные члены аппроксимируются с точностью второго порядка, а смешанные производные — по явной схеме со вторым порядком точности.

3. Анализ результатов

Расчет производился на разнесенной по пространственным координатам сетке размером 201 х 101 с диапазоном параметров 2 < < 6 2 < п < 11. Расчет температуры производился итерационно методом Ньютона — Рафсона.

Для проверки метода численно рассчитывалось течение воздуха = 3.75, = 478 К) с перпендикулярным вдувом звуковой струи воздуха (М0 = 1 п = 10.29) че-0. 1

5 -

Р Р.,

4 -

5 -4 -3 -2 -1 0 1 2 3 4 5

|14|; кривые 13

соответственно для схем Бима — Уорминга, Лакса — Фридрихса, ЕХО-схемьт

в

Рис. 3. Изолинии: а — изобары, б — распределение местного числа Маха, в — массовая концентрация водорода

и

2.0 2.5 3.0 3.5 "Но 4.0

Рис. 4. Зависимость длины отрывной зоны от числа Маха потока

Видно, что результаты, полученные с использованием ЕЬЮ-схемы, лучше согласуются с экспериментом, тогда как данные расчета с искусственным сглаживанием очень сильно занижены.

Анализ влияния числа Маха потока рассмотрен на примере истечения звуковой струи водорода с нерасчетностью п = 4 через щель шириной К = 0.3 см в поток воздуха с параметрами = 2... 4, Рг = 0.9, Ие = 104. На рис. 3, а представлено распределение давлений для двух значений числа Маха потока. Как следует из рисунка, общая картина ударно-волновой структуры для этих значений схожа. Перед струей наблюдаются головной, косой и замыкающий скачки уплотнения, которые, пересекаясь в одной точке, образуют сложную А-образную систему скачков уплотнения. Видно, что с увеличением числа Маха потока угол наклона головного скачка уплотнения существенно уменьшается. Это объясняется тем, что с ростом числа Маха потока происходит ускорение основного потока, соответственно для расширяющейся струи появляется добавочное сопротивление. При этом размеры бочки, которую можно наблюдать из рис. 3, б, где представлено распределение местного числа Маха (М = ^/и2 + т2/с, с — местная скорость звука), заметно уменьшаются.

Влияние числа Маха на дальнобойность вдуваемого водорода можно увидеть из распределения массовой концентрации водорода (рис. 3,в). Численные эксперименты показывают, что при числе Маха потока = 2 линия однопроцентной концентрации водорода максимально поднялась на высоту 1.8 см от стенки канала, в то время как для = 4 максимальная высота однопроцентной концентрации водорода располагается на высоте 1.125 см. При этом четко наблюдается уменьшение линии перетекания водорода перед струей. В соответствии с этим увеличивается не только дальнобойность вдуваемой струи (глубина проникновения водорода на рис. 3, в), но и длина отрывной зоны перед струей. Так, на рис. 4 приведена зависимость длины отрывной зоны перед струей от числа Маха потока. Из графика следует, что с увеличением числа Маха потока длина отрывной зоны линейно уменьшается.

Описанная в работе численная модель на основе ЕЬЮ-схемы позволяет моделировать течение многокомпонентных газовых смесей и учитывать такие сложные физические явления, как возникновение скачков уплотнения и вихревых зон. Изучено влияние числа Маха потока на длину отрывной зоны и дальнобойность вдуваемой водородной струи. Разработанная модель может быть распространена на пространственный случай и горение газовоздушных смесей.

Список литературы

[1] Авдуевский B.C., Медведев К.И. Взаимодействие сверхзвукового потока с поперечной струей, вдуваемой через круглое отверстие в пластине // Изв. АН СССР. Механика жидкости и газа. 1970. Т. 37, № 5. С. 193-197.

[2] Grasso F., Magi V. Simulation of transverse gas injection in turbulent supersonic air flows // AIAA J. 1995. Vol. 33, N 1. P. 56-62.

[3] Ramakrishnan R., Singh D.J. Modeling scramjet combustor flowfields with a grid adaptation scheme // AIAA J. 1994. Vol. 32, N 5. P. 930-935.

[4] Вайднер Э.Х., Драммонд Д.Ф. Расчетное исследование схемы вдува топлива в канал ПВРДсг с последовательным расположением щелей // Аэрокосм, техника. 1983. Т. 1, № 5. С. 103-111.

[5] Шунь Дж. Ш., Юнь С. Численное исследование течений с химическими реакциями на основе LU-факторизованной схемы, построенной методом симметричной последовательной верхней релаксации // Аэрокосм, техника. 1990. № 510. С. 102-113.

[6] Harten A., Engquist В., Osher S., Chakravarthy S.R. Uniformly high order accurate non-oscillatory schemes. Ill //J. of Сотр. Phys. 1987. Vol. 71. P. 231-303.

[7] Harten A., Osher S. Uniformly high-order accurate non-oscillatory schemes. I // SIAM. J. Numer. Analysis. 1987. Vol. 24, N 2. P. 279-309.

[8] Harten A., Osher S., Engquist В., Chakravarthy S.R. Some results on uniformly high-order accurate essentially non-oscillatory schemes // Appl. Num. Math. 1986. Vol. 2. P. 347-377.

[9] Kee R.J., Rupley F.M., Miller J.A. CHEMKIN-II: a Fortran Chemical Kinetic Package for the Analysis Of Gas-phase Chemical Kinetics. SANDIA Report SAND89-8009. 1989.

[10] Лапин Ю.В., Стрелец M.K. Внутренние течения газовых смесей. М.: Наука, 1989. 366 с.

[11] poinsot T.J., Lele S.K. Boundary conditions for direct simulation of compressible viscous flow // J. of Comput. Phys. 1992. Vol. 101. P. 104-129.

[12] Yang J.Y. Third-order nonoscillatory schemes for the Euler equations // AIAA J. 1991. Vol. 29, N 10. P. 1611-1618.

[13] Бекетаева A.O., Найманова А.Ж. Численное моделирование сверхзвукового течения с поперечным вдувом струй // Прикл. механика и техн. физика. 2004. Т. 45, № 3. С. 367-374.

[14] Chenault C.F., Beran P.S. К-е and Reynolds stress turbulence model comparisons for two-dimensional injection flows // AIAA J. 1998. Vol. 36, N 8. P. 1401-1412.

Поступила в редакцию 18 мая 2007 г.