Численное исследование пространственных ламинарных течений вязкого химически реагирующего газа около затупленных тел

Забарко Д.А. (d9454@mail.ru) ФГУП НПО машиностроения

На основе метода установления решения по времени с использованием принципа разделения задачи обтекания по физическим процессам разработан эффективный алгоритм численного интегрирования уравнений Навье-Стокса, позволяющий

При гиперзвуковых скоростях полета неравновесные процессы в ударном слое оказывают влияние на структуру и характеристики поля течения около летательных

аппаратов. Моделирование реальных свойств воздуха в существующих гиперзвуковых экспериментальных установках зачастую оказывается практически невозможным. Поэтому актуальной задачей является теоретическое исследование неравновесных гиперзвуковых течений газа.

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

Классический подход к решению задач обтекания состоит в разделении течения на внешнее невязкое течение и вязкий пограничный слой, с которым связана проблема нахождения условий на его внешней границе. B рамках уравнений Навье-Стокса решение задачи проводится для всей области течения и указанные трудности не возникают.

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

рассчитывать

гиперзвуковые

неравновесные

течения

около

затупленных тел.

Уравнение неразрывности для смеси имеет вид:

др dt

+ div

( -л

pV V J

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

р

d dt

dp ~dt

+ div

V V

{ d v

grad-

2 1 -- - ---V x rot V

2 - - -/и V div V - q

В приведенных уравнениях Здесь h - энтальпия газа, p - давление, р - плотность,

V - вектор скорости, / - динамическая вязкость, q - тепловой поток, S -тензор скоростей деформации, Div - дивергенция тензора.

Набегающий сверхзвуковой поток рассматривается как смесь, состоящая из 23,3% кислорода (O2) и 76,7% азота (N2). Обтекающий тело газ состоит из семи компонентов: O,

N, NO, O2, N2, NO + , e-, а в качестве возможных химических реакций, протекающих в ударном слое при высоких температурах, принимаются следующие [7]:

Здесь M любой из 6 рассматриваемых компонентов, являющихся катализаторами, e-электронная компонента.

Введем следующую безразмерную концентрацию

где ci - массовая концентрация i -го компонента, mt - молекулярный вес i -го компонента,

тх - молекулярный вес газа набегающего потока.

Будем ассоциировать индекс i = 1, 2, 3, 4, 5, 6 с компонентами O, N, NO, O2, N2,

NO + соответственно. Условия сохранения атомарного состава и квазинейтральности смеси запишутся в виде:

г/4 = 0.21 - 0.5(п +Пз +П6), П = 0.79 - 0.5(/2 +Пз + П),

где г/7 - концентрация электронной компоненты.

Коэффициент вязкости неравновесной смеси рассчитывается по формуле Буденберга-Уилки [1]:

i=1 V k*i )

где xi

мольная концентрация i -го компонента, m - молекулярный вес смеси,

Г 1 1

Г mm

1 + -

V m, )

Коэффициенты вязкости в (i =1,--,6) чистых (однокомпонентных) газов вычисляются по формуле

В = 2.6693 10-6

mi T

а2 Ц(22)

где T - температура, а для вычисления интегралов столкновений Д22-* применяется следующая формула:

/ \ -0.1472

1.157Г kT

Д(2,2 )

Здесь k = 1.38 Ю-23 Дж/К - постоянная Больцмана, сг - диаметры столкновений, si -характеристический потенциал столкновений i -го компонента. Данные о параметрах сг и si приведены в таблице 1.

Таблица 1

Характеристические потенциалы и диаметры столкновений

ci m

Тепловой поток к поверхности тела обусловлен теплопроводностью и энергией диффундирующих компонентов, поэтому

6 и

q = -Я grad Т - пг grad уг

г=1 Sc

где у. =-, h - энтальпия г -го компонента, Sc. =- D - числа Шмидта (для нейтральных

х р

компонентов принимают равными Sc7 =0.5, а для заряженных - Sc7 =0.25, в соответствии с моделью амбиполярной диффузии), Л - коэффициент теплопроводности, D7 - эффективный

коэффициент диффузии г -го компонента.

Будем рассматривать газ в ударном слое как бинарную смесь атомов и молекул одного вещества, учитывая бинарную диффузию с коэффициентами, определяемыми числом Льюиса. Такое предположение хорошо выполняется для воздуха ввиду близости кинетических коэффициентов и атомных весов кислорода и азота. Кроме того, в работе [10] показано, что концентрации компонентов слабо зависят от степени упрощения диффузионной модели для рассматриваемых режимов обтекания.

С учётом сделанных предположений имеем

q=-grad h+-S (Lei- 1)h grad y - p- S y grad ,

Pr Pr i=1 Pr i=1

Л

где Pr = / = 0.72 - число Прандтля, cp - удельная теплоемкость смеси при постоянном

давлении, Le7 =--числа Льюиса, evi - колебательная энергия молекул г -го компонента

(для атомов ev1 = ev 2 = 0 ),

h = S hy = £ y(cpRAT + eu + h0 ) = -*-p 7=7 7=7 У* -1 р

y* -1 р

статическая энтальпия, h°- энтальпия образования г-го компонента (табл.2), у*

эффективный показатель адиабаты.

Таблица 2

Энтальпии образования компонентов

Распределение энергий по поступательным и вращательным степеням свободы принимается равномерными, т.е. молярные теплоемкости cpi при постоянном давлении

равны 5/2 для атомов и 7/2 для молекул. Колебания молекул считаются возбужденными равновесно, или по модели Лайтхилла. В работе [11] показано слабое различие между моделью Лайтхилла и строгой равновесной моделью учёта колебательной энергии молекул, вследствие того, что она составляет лишь небольшую часть от величины полной энтальпии. Колебательная энергия молекул для модели Лайтхилла принимается равной половине максимально возможного равновесного значения при данной температуре, т. е.

evr =

i =3,..Д

где Ra - универсальная газовая постоянная. При равновесной модели энергия колебаний двухатомной молекулы, просуммированная по всем колебательным уровням, выражается формулой [8]

где 6v - характеристическая температура (см. таблицу 3).

Таблица 3

Характеристические температуры колебаний молекул 9V, K

Уравнение сохранения массовой концентрации i -го компонента имеет вид

i = 1,2,3,6,

+ -div Ut = Wt dt p

где Ui = pDi gradf. - массовый поток i -го компонента вследствие диффузии, Wi массовая скорость образования i -го компонента в единице объема за единицу времени.

Выражения для источниковых членов реакций Щ ]

г с

имеют вид:

#1 = Р11 + Р13 + Р14-p15 -p>16> = Р22 + Р13 + Р15 -p14 -p16> #3 = -p13 + Р14 + Р15, #6 =р16,

где функции ptj определены следующим образом:

у4 -ру1

Р13 = k3р(к3у3 -ру1у2 ), р14 =-k4р(к4у1у3 -У2у4 ),

Р15 = k5p(k5y1y5 -у2у3 ) ,

Р22 Р16

:2k 2Р(К:

у5 -ру

здесь: у

р = рр-

размерная концентрация i -го компонента [ [yt]

6у1у2 -У6

m = 28.8-

размерная плотность

см3

константы скоростей обратных

реакций, Кг - константы равновесия.

Константы, используемые в функциях ptj, имеют следующие размерности:

[К, ]

[Kj ] - безразмерные,

где i = 1,2,3, j = 4,5,6.

Вводятся следующие безразмерные величины

[k]=

с

[kj]

моль сJ

Р р =Р Т = -

без без

V6e3p

x y z t

x = - y = z = - t = -

безр j безр j безр j безр j

- W =

Щ mгл1Р-Р-

где L - характерный линейный размер. Далее индекс безр опускается.

Дополняют систему уравнений сведения о константах равновесия и скоростях химических реакций [1,5] (табл. 4), и безразмерное уравнение состояния, которое принимает

вид

г

г

Таблица 4

Константы скоростей химических реакций

Константа скорости обратной реакции k

№ реакции

Катализатор

Константа равновесия К

N, N0

К, = 1.2 Ю3 Т

-0.5

595001

Т

3 1015Т-0.5

2k1 9k,

0 2, 0, N0

К 2 =18exp

113000

Т

k2 = 1.056 1016 T

2.47k 2 2.15 105 T lk

02 , N2

N0 , 0 , N

К 3 =

4 exp

755001

k3 =

9.75 1019T 2 20k3

К 4 = 0.24exp

16127

k 4 = 1.3 -1010 T exp

35731

К5 = 4.5exp

380001

k5 = 1.6 Ю13

A = 1.4 10-8 T + +1.2 10-12 T2 + 1.4 10-16 T3

325001

A-1011

k6 =1.8 1021T

Сформулируем основные положения, при которых исследуется рассматриваемая

задача:

1) набегающий на затупленное тело поток газа сверхзвуковой, однородный и невозмущённый (у = 1.4);

2) течение во всей возмущённой телом области - ламинарное и симметричное относительно плоскости, проходящей через ось симметрии тела;

3) уравнения Навье-Стокса справедливы для описания течения во всей возмущённой телом области;

4) из вязкостных эффектов учитывается динамическая вязкость и . Условия симметричности течения имеют следующий вид:

du = dv = др = dT = w = 0 dx dx dx dx

где u, v, w - компоненты вектора скорости V, переменная x соответствует координатной линии, перпендикулярной плоскости симметрии.

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

uw = vw = ww = 0 .

w w

На стенке задается температура Tw = Tw = const или изменение температуры по

закону, близкому к линейному. Индекс w соответствует параметрам течения газа, задаваемым на стенке.

Для некаталитической поверхности стенки граничные условия для концентрации компонентов (атомов O и N, молекул NO) имеют вид

а для заряженных частиц используется условие идеально каталитической поверхности Для обеспечения устойчивости расчёта налагалось ещё одно условие:

fd p

Переход набегающего потока через ударную волну происходит замороженным образом, т.е. концентрации компонентов здесь совпадают с соответствующими значениями концентрации невозмущенного течения. Для модели Лайтхилла это означает, что на ударной волне со стороны обтекаемого тела у*=1.33333.

Для задания тонкой головной ударной волны ( внешней границы области интегрирования) необходимы априорные данные о структуре поля течения. Для больших чисел Рейнольдса (Reш>103) можно пренебречь влиянием структуры тонкой головной ударной волны на течение вниз по потоку и принять, что ударная волна является поверхностью разрыва газодинамических параметров, на которой выполняются нестационарные условия Рэнкина-Гюгонио:

p(Vn -D)=px(Vn-D),

p + p(Vn - D )2 = px+px(Vn- D )2, h(p,p) +1 (Vn -D)2 = K(px,pj + 2D)

->

где D - скорость распространения волны по частицам газа, Vn,Vт - проекции вектора скорости на нормаль и касательную плоскость к поверхности ударной волны.

Математическая модель

Расчёт течений вблизи головных частей затупленных тел будем проводить в нормированной сферической системе координат x1=, x2=p, xera.).

R - RT (в, p)

Rb ((,в,р)-Rt (в,р)

R = R- (в<р) + (Rb (t,e,p)-R- (в,р)) .

x = R sinesinp, y = -R sinecosp, z = z0 - R cose,

0 <£< 1, 0 <p<, 0 <в<ж/2.

Рис. 1.Сферическая система координат.

Рис. 2 Нормированная система координат.

Здесь RT = RT (e,p), RB = RB (t,e,p) - уравнения, задающие поверхности тела и ударной

волны соответственно, причем RB определяется в процессе решения задачи.

Основой разработанного численного метода является расщепление задачи обтекания по физическим процессам.

I этап. Интегрируется невязкая часть системы, включающая в себя уравнения движения Эйлера, неразрывности и энергии только с невязкими членами:

U+F£ + Ев+G-+H = 0,

где

- - U = U

U =

f = £RF+e E+ф+ф, E = E,

H = H+r E+f-P +£-U,

P + puR

F = puRug

pURU-

(e+p k

- - G = G.

\pue

pURUg

p+pug

pugu- (E + P )u

puRu

sing

llgll -

P + pu- (E + P )u-

H = 1

[2puR + pugctgg

p((u - - u- )+ 3puRug p(ul - u- )+ 3puRug I pu-((2ugctgg + 3uR )

(E + P )(2uR + ugctgg)

uR, u

физические компоненты скорости в сферической системе координат,

E = p

e + -

полная внутренняя энергия единицы объема, e

. 1 £ у-1 p

внутренняя

энергия единицы массы газа, индексы при переменной £ обозначают дифференцирование по соответствующей координате.

Невязкая часть системы уравнений решается при помощи явной двухшаговой конечно-разностной схемы Мак-Кормака.

Введем в пространстве t,£,-,g сетку с шагами: At, А£, А- = л/K, Ag, обозначим

координаты узлов сетки tJ = jAt,gm = mAg, -к = kA-, а любую функцию f в точке как /г!тк, при этом 0 < n < N,0 < m < M,0 < k < K , где N, M, K - число интервалов в направлениях g, - соответственно.

Наличие узких областей с большими градиентами параметров потока, таких как пограничный слой вблизи тела, зона размазанной ударной волны при низких числах Рейнольдса, не позволяет получать надёжные количественные результаты при