CC BY
22
УДК 533.6.011.5+532.582.33
Д. А. Забарко, В. П. Котенев
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ЛАМИНАРНЫХ ТЕЧЕНИЙ ВЯЗКОГО ХИМИЧЕСКИ РЕАГИРУЮЩЕГО ГАЗА ОКОЛО ЗАТУПЛЕННЫХ ТЕЛ
На основе метода установления решения во времени с использованием принципа разделения задачи обтекания по физическим процессам разработан эффективный алгоритм численного интегрирования уравнений Навье-Стокса, позволяющий рассчитывать гиперзвуковые неравновесные течения около затупленных тел.
При гиперзвуковых скоростях полета неравновесные процессы в ударном слое оказывают влияние на структуру и характеристики поля течения газа около летательных аппаратов. Моделирование реальных свойств воздуха в существующих гиперзвуковых экспериментальных установках часто оказывается практически невозможным. Поэтому актуальной задачей является теоретическое исследование неравновесных гиперзвуковых течений газа.
Постановка задачи. Классический подход к решению задач обтекания состоит в разделении течения на внешнее невязкое течение и вязкий пограничный слой, с которым связана проблема нахождения условий на его внешней границе. B рамках уравнений Навье-Стокса решение задачи проводится для всей области течения и указанные трудности не возникают.
Для получения стационарного решения задачи обтекания тела равномерным сверхзвуковым потоком используется метод установления решения по времени для нестационарной системы уравнений Навье-Стокса
р— = — gradp — 2 grad(p divV) + 2 Div(pS). dt 3
Уравнение неразрывности для смеси имеет вид
I + div(pV) = 0.
Уравнение энергии записывается следующим образом:
d / V2 \ dp
pJt\h + tJ = т +
f f V2 1- Л 2 - - \ + div ( 2p i grad— — - V x rotV j — - pV divV — q j.
В приведенных уравнениях к — энтальпия газа, р — давление, р — плотность, V — вектор скорости, р — динамическая вязкость, д — тепловой поток, £ — тензор скоростей деформации, Б1у — дивергенция тензора.
Набегающий сверхзвуковой поток рассматривается как смесь, состоящая из 23,3% кислорода (02) и 76,7% азота (N2). Обтекающий тело газ состоит из семи компонентов: О, N N0, 02, N0+, е-, а в качестве возможных химических реакций, протекающих в ударном слое при высоких температурах, принимаются следующие [7]:
1) 02 + М ^ 20+ М,
2) N + М ^ 2^М,
3) N0 + M ^ N + 0 + М,
4) N0 + 0 ^ О2 + N
5) N2 + 0 ^ N0 + N
6) N + 0 ^ N0+ + е-.
Здесь М — любой из шести рассматриваемых компонентов, являющихся катализаторами, е- — электронный компонент.
Введем следующую безразмерную концентрацию:
Пг = -,
тг
где сг — массовая концентрация ¿-го компонента, тг — молекулярный вес ¿-го компонента, тж — молекулярный вес газа набегающего потока.
Будем ассоциировать индекс ¿ = 1, 2, 3, 4, 5, б с компонентами О, N N0, 02, N0+ соответственно. Условия сохранения атомарного состава и квазинейтральности смеси запишутся в виде
П4 = 0,21 - 0,5(т + Пз + Пб), П5 = 0,79 - 0,5(^2 + Пз + Пб), Пб = П7,
где — концентрация электронного компонента.
Коэффициент вязкости неравновесной смеси рассчитывается по формуле Буденберга-Уилки [1]:
б / х -1
р = ^ Рг хА Хг + йгк Хк
г=1 ^ к=г
где хг = сгт/тг — мольная концентрация ¿-го компонента, т — моле-
кулярныи вес смеси,
Gik —
Mk
! M Vi шЛ
тг)
2,/2(1 + m)
V mk
г — к.
Коэффициенты вязкости , i — 1,..., 6, чистых (однокомпонент-ных) газов вычисляются по формуле
6 v miT M — 2,6693 • W-6^-^:
где Т — температура, а для вычисления интегралов столкновений П»2'2)* применяется следующая формула:
Пг(2'2)* — 1,157^ kpj
-0,1472
Здесь к = 1,38 • 10-23 Дж/К — постоянная Больцмана, а» — диаметр столкновений, е» — характеристический потенциал столкновений ¿-го компонента. Данные о параметрах а» и е» приведены в табл. 1.
Таблица 1
Характеристические потенциалы и диаметры столкновений
i O N NO, NO+ O2 N2
т ® 106,7 71,4 116,7 106,7 71,4
<Л(А) 3,050 3,298 3,492 3,467 3,798
2
Тепловой поток к поверхности тела обусловлен теплопроводностью и энергией диффундирующих компонентов, поэтому
6
д = -А grad Т - V ~^кг grad 7», »=1 Бс»
где 7» = п/тж, к» — энтальпия ¿-го компонента, Бс» = рБ»/р — числа Шмидта (для нейтральных компонентов Бс» = 0,5, а для заряженных — Бс» = 0,25 в соответствии с моделью амбиполярной диффузии), А — коэффициент теплопроводности, Б» — эффективный коэффициент диффузии ¿-го компонента.
Будем рассматривать газ в ударном слое как бинарную смесь атомов и молекул одного вещества, учитывая бинарную диффузию с ко-
эффициентами, определяемыми числом Льюиса. Такое предположение хорошо выполняется для воздуха ввиду близости кинетических коэффициентов и атомных весов кислорода и азота. Кроме того, в работе [10] показано, что концентрации компонентов слабо зависят от степени упрощения диффузионной модели для рассматриваемых режимов обтекания.
С учетом сделанных предположений имеем
6 6 -q = т^т grad h + (Lei - grad Yi - p ^ Y grad em,
i=1 i=1
где Pr = pcp/X = 0,72 — число Прандтля, cp — удельная теплоемкость смеси при постоянном давлении, Lei = Pr/ Sc — числа Льюиса, evi — колебательная энергия молекул i-го компонента (для атомов
evi = eV2 = 0),
66
Оч Y* Р
h = ^ hiYi = ^2 Yi(cpiRAT + evi + h0)
сргКл1 +
*=1 *=1 ^ - 1 Р — статическая энтальпия, к0 — энтальпия образования ¿-го компонента (табл. 2), 7* — эффективный показатель адиабаты.
Таблица 2
Энтальпии образования компонентов
Компонент h0, кал/г h0, кал/моль
O 3680 5,8990-104
N 8030 1,1250-105
NO 715 2Д810-104
O2 0 0
N2 0 0
NO+ 8020 2,406-105
Распределения энергий по поступательным и вращательным степеням свободы принимаются равномерными, т. е. молярные теплоемкости Срг при постоянном давлении равны 5/2 для атомов и 7/2 для молекул. Колебания молекул считаются возбужденными равновесно или по модели Лайтхилла. В работе [14] показано слабое различие между моделью Лайтхилла и строгой равновесной моделью учета колебательной энергии молекул вследствие того, что эта энергия составляет лишь небольшую часть величины полной энтальпии. Колебательная энергия
молекул для модели Лайтхилла принимается равной половине максимально возможного равновесного значения при данной температуре, т. е.
ЯлТ
ev( —
2
г — 3,..., 6,
где Ял — универсальная газовая постоянная. Для равновесной модели энергия колебаний двухатомной молекулы, просуммированная по всем колебательным уровням, выражается формулой [8]
ЯлOv
ev(
Jv /T _ 1'
где — характеристическая температура (см. табл. 3).
Таблица 3
Характеристические температуры колебаний молекул в„, K
NO, NO+ O2 N2
2690 2230 3340
Уравнение сохранения массовой концентрации ¿-го компонента имеет вид
^ + 1 div и» = IV» ¿ = 1, 2,3, 6, at р
где и» = —рБ» grad 7» — массовый поток ¿-го компонента вследствие диффузии, IV» — массовая скорость образования ¿-го компонента в единице объема за единицу времени.
„ • (Гт • 1 Г моль"
Выражения для источниковых членов реакций IV» I \уУ»\ = -
имеют вид г с
Щ = ^11 + ^13 + £14 — £15 — £16,
^Из = £22 + £13 + £15 — £14 — £16, Щ = — £13 + <£14 + £15,
Щ = <16,
где функции <»3 определены следующим образом:
£11 = 2к1р (К174 — р072), <15 = к5р (К57175 — 7273),
£13 = к3р (К373 — Р7172), £22 = 2к2р (К75 — р7г),
£14 = —к4р (К47173 — 7274), £16 = к6р (^671-72 — 7е);
здесь 7» = п»/тж — размерная концентрация ¿-го компонента ([7»] = = [моль/г, тж = 28,8г/моль), р = ррж — размерная плотность
([Рго] = [г/см3]), кг — константы скоростей обратных реакций, Кг — константы равновесия.
Константы, используемые в функциях ^, имеют следующие размерности:
[Кг] =
моль" [k ] = [■ см6 ]
- см3 . , [ki] = .моль2 • с.
[Kj ] — безразмерные, [kj ] =
см3
l-моль • CJ
где i = 1, 2,3, j = 4, 5, 6.
Вводятся следующие безразмерные величины:
рбезр
р
Ро
рбезр
Р
_ T
T безр m ,
Ро
^безр
рбезр
Р
V
\/Рю/рО
хбезр
X
V
убезр
V
^безр
V
. = fpö ф = WfimiVpOp0
^безр г А / , ф i безр
4P»
ь
где Ь — характерный линейный размер.
Дополняют систему уравнений сведения о константах равновесия и скоростях химических реакций [1, 5] (см. табл.4) и безразмерное уравнение состояния, которое принимает вид
Р = PT~
m
Сформулируем основные положения, при которых исследуется рассматриваемая задача:
1) набегающий на затупленное тело поток газа — сверхзвуковой, однородный и невозмущенный (7 = 1,4);
2) течение во всей возмущенной телом области — ламинарное и симметричное относительно плоскости, проходящей через ось симметрии тела;
3) уравнения Навье-Стокса справедливы для описания течения во всей возмущенной телом области;
4) из вязкостных эффектов учитывается динамическая вязкость р.
Условия симметричности течения имеют следующий вид:
ди дь др др дТ дх дх дх дх дх
где и, ь, ы — компоненты вектора скорости V, переменная х соответствует координатной линии, перпендикулярной плоскости симметрии.
= w = 0,
Р
00
OG
z
Таблица 4
Константы скоростей химических реакций
№ реакции
Катализатор
Константа равновесия K
Константа скорости обратной реакции k
1
N, NO N2 O2 O
Kl = 1,2 • 103T-°>5x 59500 1
x exp < - -
k1 = 3 • 1015T-°'5 2k1 9k1 25k1
O2, O, NO N2 N
K2 = 18 exp <-
( 113000 \
16^-°,5
k2 = 1,056 • 1016T 2,47k2 2,15 • 105T-1k2
O2, N2 NO, O, N
K3
75500 4 exp --—\
-
ka = 9,75 • 1019T-3/2 20k3
16127
K4 = 0,24 exp j--— f
k4 = 1,3-1010 T expj - 3573}
K5 = 4,5 exp -
-
38000 1
k5 = 1,6 • 10
13
K6 = A1011 x
( 325001
X 6XP\ ^J kfi = 1,8 • 1021T-3/2
A =1,4 • 10-8T+ + 1,2 • 10-12T2 + +1,4 • 10-16 T3
В качестве граничных условий на поверхности обтекаемого тела для компонентов скорости задаются условия прилипания:
Uw = vw = Ww = 0.
На стенке задается температура Tw = TW = const или изменение температуры по закону, близкому к линейному. Индекс "w" соответствует параметрам течения газа, задаваемым на стенке.
Для некаталитической поверхности стенки граничные условия для концентрации компонентов (атомов O и N, молекул NO) имеют вид
^ w 0'
а для заряженных частиц используется условие идеально каталитической поверхности
(ce)w = 0.
Для обеспечения устойчивости решения налагалось еще одно усло-
вие:
тт) =0-
dnj w
2
3
4
5
Рис. 1. Сферическая система координат
Переход набегающего потока через ударную волну происходит "замороженным" образом, т. е. концентрации компонентов здесь совпадают с соответствующими значениями концентрации невозмущенного течения. Для модели Лайтхилла это означает, что на ударной волне со стороны обтекаемого тела % = 1,33333.
Для задания тонкой головной ударной волны ("внешней" границы области интегрирования) необходимы априорные данные о структуре поля течения. Для больших чисел Рейнольдса (Яе^ > 103) можно пренебречь влиянием структуры тонкой головной ударной волны на течение вниз по потоку и принять, что ударная волна является поверхностью разрыва газодинамических параметров, на которой выполняются нестационарные условия Рэнкина-Гюгонио:
р(Уп - Б) = рх(Уп>х - Б), Р + р(Уп - Б)2 = р^ + рх(Уп>х - Б)2,
%,р) + 2(Кг - Б)2 = ) + 2(Уп,ж - Б)2,
К = Ут><х,
где Б - скорость распространения волны по частицам газа, Уп ,УТ — проекции вектора скорости на нормаль и касательную плоскость к поверхности ударной волны.
Математическая модель. Расчет течений вблизи головных частей затупленных тел будем проводить в нормированной сферической системе координат х1 = £, х2 = у, х3 = в (рис. 1, 2):
Z =
R - Rt(6, у)
Rb (t,6,y) - Rt (6, у)'
то есть
R = Rt(6, у) + £(Rb(t, 6, у) - Rt(6, у)),
x = R sin 6 sin у, y = — R sin 6 cos у, z = z0 — R cos 6,
0 < £ < 1, 0 < у < n, 0 < 6 < n/2.
Рис. 2. Нормированная система координат
Здесь RT = RT(#,<), RB = Rb (t, — уравнения, задающие поверхности тела и ударной волны соответственно, причем RB определяется в процессе решения задачи.
Основой разработанного численного метода является расщепление задачи обтекания по физическим процессам.
I этап. Интегрируется "невязкая" часть системы, включающая в себя уравнения движения Эйлера, неразрывности и энергии только с "невязкими" членами:
Ut + F + E, + G ^ + H = 0,
где
U/ = J, F = ÇrF + Ç, E + Ç^G + ÇtJ, E = E, G = G,
H = H + ÇR F + £ + Ç^ G + &J,
pu, pnBue
2
P + Pn2 , pu, n^ (E + p) n, _
U =
P
PUR
pue E
F =
PUR
p + PuR
pURUe
PUR U^
(E + p) UR
E = R
G=
1
R sin в
pu^ 2p ur + PUe ctg в
P uRu^ ■ h=R P (2uR - u2 - U2) + 3pURUe
PUe U^ p (u2 - u2) + 3pURUe
P + P u2 pu^ (2ue ctg в + 3ur)
_ (E + p) u^ _ (E + p) (2ur + ue ctg в)
uR, uv, ua — физические компоненты скорости в сферической системе координат,
V2 ~
E = р
e+т
1
Р
вну-
— полная внутренняя энергия единицы объема, е =
1 - 1 р
тренняя энергия единицы массы газа, индексы при £ обозначают дифференцирование по соответствующей координате.
"Невязкая" часть системы уравнений решается при помощи явной двухшаговой конечно-разностной схемы Мак-Кормака.
Введем в пространстве Ь, £, у, в сетку с шагами ДЬ, Д£, Ду = п/К, Дв, обозначим координаты узлов сетки Ь = ]ДЬ, вт = тДв, ук = кДу, а любую функцию f в точке обозначим как , при этом 0 < п < N, 0 < т < М, 0 < к < К, где И,Ы,К — число интервалов в направлениях £, в, у соответственно.
Наличие узких областей с большими градиентами параметров потока, таких как пограничный слой вблизи тела, зона размазанной ударной волны при низких числах Рейнольдса, не позволяет получать надежные количественные результаты при использовании равномерных сеток. Для адекватного учета поведения параметров в областях больших градиентов вводится сгущение сетки в радиальном направлении вида
£(n) =
Г 1
N,
1 1 1 1
тг +
если n > N„
N ' AN + AN„ N sinf -П(2П — NY), если n < N„
2
2
2N
где А - произвольное целое число, задающее отношение наибольшего шага к наименьшему, N — номер сеточного узла, с которого вводится неравномерная сетка.
Предиктор разностной схемы имеет вид
Uj+1 = Uj — At
n,m, k n ,m, k
uj - иj
n+1 ,m, k n ,m, k
A£
Ej - Ej Gj - Gj
n,m+1,k n,m,k n,m,k+1 n,m,k Uj
+ A6 1 Ay n,m,k а корректор записывается следующим образом:
f
Uj+1 = 11 Uj + Uj+1 — At
n,m,k
n,m,k
Uj+1 _ Fj
n,m,k n— 1,m,k
?j+1
At
Ej+1, -Ej+1
Gj+1 _ gj+1
n,m,k n,m-1,k n,m,k n,m,k-1 Uj+1
+ » ч + » + 1 n,m,k
Ad
Ay
II этап. Решение уточняется с использованием "вязкой" части уравнений Навье-Стокса и энергии:
0
Ut = K, K =
- 2 grad (ß divK/) +2 Div (^S?
/ V2 1 - -Л 2 - -div ( 2ß I grad—--- V x rotK I - - ßV divK - q
"Вязкая" подсистема интегрируется при помощи метода прогонки:
и}]+2 _тт3+1 _ К п,т,к п,т,к _ К3+2
1 = ' ДЬ = '
Для решения "вязкой" подсистемы уравнений используется неявная разностная схема. Для ее описания представим "вязкую" часть уравнений количества движения и энергии в виде
ЗГ л д / ЗГ\ л д / ЗГ4
Ж = Алд£ + д£ {^ж 1 +
л д . ,
+ Aq^z ßte
dF
- за + вм ж + +СГ+Б
где А, В, С — коэффициенты при вторых производных, первых производных и самой функции соответственно, слагаемое Б содержит все оставшиеся члены уравнений со смешанными производными и свободными членами. Используя выражения конечно-разностных аналогов производных для сетки с переменным шагом, запишем разностный аналог уравнения, оставляя в качестве неизвестных целевые функции на новом временном слое для каждого уравнения в радиальном направлении. Все значения остальных функций считаются известными и берутся с текущего временного слоя.
Для получения решения на {г + 2)-м слое при заданных граничных условиях интегрирование одномерных уравнений сводится к последовательному решению отдельных разностных уравнений с трехдиаго-
нальной матрицей методом прогонки. Эти уравнения имеют стандартный вид:
An/n+i — Bnfn + Cnfn-i = Dn, 0 < n < N,
где в качестве целевой функции / рассматривается любая из неизвестных функций, /0 и fN заданы или определяются из граничных условий.
Таким образом, в "вязкой" части системы последовательно решаются разностные уравнения движения, записанные неявно относительно искомых параметров. Затем с использованием полученного поля скоростей решаются уравнения энергии и диффузии. Применение описанной процедуры уменьшает ограничение на шаг интегрирования по времени, связанное с "вязким" критерием устойчивости.
III эта п. Для интегрирования уравнений диффузии метод прогонки применяется в сочетании со специальной неявной разностной схемой, используемой для аппроксимации источниковых членов. При этом сначала при помощи неявного метода Ньютона совместно решается система разностных уравнений с источниковыми членами:
р^ = wWi, г = 1, 2, 3, 6. dt г' ' ' '
IV этап. Затем учитывается диффундирование компонентов смеси:
р^ + divUi = 0. dt
Уравнения последней системы решаются последовательно и независимо друг от друга методом прогонки.
Предложенный метод расщепления уравнений диффузии позволяет рассчитывать химически неравновесные течения в широком диапазоне значений параметров набегающего потока вплоть до высот полета, где реализуется околоравновесный режим протекания химических реакций.
Для исследования течений газа около удлиненных тел используется разделение всей области интегрирования на ряд взаимно перекрывающихся подобластей (рис. 3) и последовательный расчет в каждой из них. Сначала задача решается в окрестности затупления. Затем центр сферической системы координат переносится по оси тела. Выстраивается новая расчетная область, где на левой границе области задаются "жесткие" граничные условия из уже рассчитанной области, на выходных границах — "мягкие" граничные условия вида линейной экстраполяции искомых функций, на теле — условия прилипания и температура стенки, на ударной волне — нестационарные соотношения
Рис. 3. Разделение физической области интегрирования
Рэнкина-Гюгонио. Решение в полученной области устанавливается, и описанная процедура построения новой расчетной сетки повторяется. Описанное разделение области возможно в силу слабой передачи возмущений вверх по потоку при обтекании тел сверхзвуковым набегающим потоком вязкого газа [6] и позволяет проводить расчет длинных затупленных тел до 100 калибров и более.
Из-за большой сложности уравнений Навье-Стокса невозможно получить аналитическое выражение устойчивости для схемы Мак-Кормака. Однако практическое применение предложенного численного метода показало возможность использования, без потери устойчивости счета, эмпирической формулы
Д£ = аД^кФЛ,
где о — коэффициент запаса (о = 0,8); Д£КФЛ определяется по критерию Куранта-Фридрихса-Леви для линейных гиперболических уравнений в частных производных [9].
Перед очередным шагом по времени для всех точек сетки рассчитывается минимальный шаг интегрирования Д£, который используется для получения решения на следующем временном слое.
Анализ результатов. В настоящем пункте приводятся результаты расчетов осесимметричного обтекания сферически затупленных конусов, а также сравнения с экспериментальными и расчетными данными других авторов [2, 11-13]. Значения параметров для рассмотренных режимов приведены в табл. 5. В первых двух случаях рассматривается обтекание сферически затупленного конуса на высотах полета 45 км и 75 км. Температура на стенке Т№ задается по кусочно линейному закону: сначала на сфере и далее по длине конуса до 1,4 калибра Т№ = Т№1, затем происходит линейное уменьшение температуры до Тте2 для сечения, находящегося на расстоянии от носка, равном 101,5 калибра, и далее Т№ = Тте2. Третий режим соответствует течению на затупленном конусе с углом полураствора в = 9°, углом атаки а = 0° для высоты полета 70,104 км.
Таблица 5
Параметры расчетных случаев
№ H, км V», м/с тж,, к а, град ©, град Re», м 1 Rn, м Tw,K
1 45 4500 13,64 271,28 0 11,02 517218 0,04 Tw1 = 2300 TW2 = 1100
2 75 6807 23,89 202,06 0 9,5 22548 0,04 TW1 = 500 Tw2 = 300
3 70,104 8059,4 27,18 218,79 0 9 51454 0,1524 1000
На рис. 4-11 приведены результаты расчетов режимов 1, 2 и их сравнение с данными, полученными при использовании метода, изложенного в работе [2]. Число Стантона имеет вид
St =
P'X1V'X1(H,X К)
föT
где q = Л —— — тепловой поток, Л — коэффициент теплопровод-
\оп/
ности, =
Y P V
I 1 го . Vc
+---полная энтальпия газа в однородном
2
7 - 1 Р набегающем потоке.
Сравнивая расчеты, проведенные для неравновесного течения и совершенного газа, можно отметить, что отличия тепловых потоков для первого режима (см. рис.4) достигают 10% в области "ложки", увеличиваясь до 20% далее по длине конуса. Для второго режима
Рис. 4. Распределение числа Стантона:
режим течения 1; высота полета 45 км; число Маха М = 13,64; • — данные по методу работы [2]; ▲ — неравновесный расчет; □ — расчет по модели совершенного газа
2
2.4Е-02 1.8Е-02 Й1-2Е-02 6.0Е-03 О.ОЕ+ОО
О 20 40 60 80 100 120 140 Калибры
Рис. 5. Распределение числа Стантона:
режим течения 2; высота полета 75 км; число Маха М = 23,89; • — данные по методу работы [2]; ▲ — неравновесный расчет; □ — расчет по модели совершенного газа
(см. рис. 5) отличия тепловых потоков незначительны (до 5 %) вследствие близости режима к "замороженному" и, соответственно, малости концентраций химических компонентов в пристеночной области ударного слоя. На рис. 4,5 также приводятся сравнения числа Стантона с результатами, полученными по методу работы [2]. Наблюдается хорошее согласование неравновесных расчетов (отличия не превышают 5 %), выполненных по предложенному численному методу и методу автора работы [2].
Характерной особенностью неравновесных режимов течения является наличие явно выраженной "ложки", в отличие от модели совершенного газа, когда данная "ложка" отсутствует [4]. Тепловой поток для неравновесного течения оказывается на 20 % меньше теплового потока, рассчитанного без учета реальных свойств газа.
На рис. 6-11 приведены безразмерные концентрации компонентов смеси Пг, плотность и температура поперек ударного слоя для сечения, находящегося на расстоянии от носка, равном 10 калибрам. Качественное согласование позволяет судить о правильности расчета химических концентраций в структуре ударного слоя. Количественные отличия рассматриваемых параметров обусловлены, скорее всего, разными моделями химических реакций сравниваемых методик, включающих в себя константы скоростей химических реакций и граничные условия. При этом отмечается хорошее согласование концентраций (атомарный кислород), достигших сравнительно больших величин. Несколько большее различие наблюдается в концентрациях компонентов, величины которых малы. Для сравнения приведены концентрации, рассчитанные с применением грубой и подробной сеток. Первая расчетная
Рис. 6. Атомарный кислород О:
• — данные по методу работы [2]; ▲
Рис. 7. Атомарный азот N:
• — данные по методу работы [2]; ▲
сетка 1 (40 узлов); □ — сетка 2 (80 узлов) сетка 1 (40 узлов); □ — сетка 2 (80 узлов)
0.0 0.2 0.4 0.6 0.8 ъ,
Рис. 8. Окись азота NO:
• — данные по методу работы [2]; ▲ — сетка 1 (40 узлов); □ — сетка 2 (80 узлов)
0.0 0.2 0.4 0.6 0.8 Е,
Рис. 9. Ионы окиси азота NO+:
• — данные по методу работы [2]; ▲ — сетка 1 (40 узлов); □ — сетка 2 (80 узлов)
сетка имеет N = 40 узлов поперек ударного слоя, параметры сгущения А = 5, N = 40, параметры сетки по координате в: М = 46, Дв = 2°. Вторая расчетная сетка имеет N = 80 узлов, параметры А = 50, N = 80, М = 46, Дв = 2°. Шаги грубой и подробной сеток в пристеночной области различаются в 17 раз. На рис. 6-10 показано слабое влияние подробной сетки на распределение концентраций поперек ударного слоя. Ввиду того, что в программе автора работы [3] температура стенки принималась постоянной (Т№ = 2300 К), были также проведены соответствующие расчеты. На рис. 10 показано, что изменение температуры стенки на 100 градусов практически не влияет на распределение электронов в ударном слое. Сравнение температурного профиля для совершенного и реального газов (рис.11) показывает значительное отличие распределения температур поперек ударного слоя в пристеночной области, где наблюдаются значительные концентрации химических компонентов.
Рис. 10. Распределение электронов e-:
• — данные по методу работы [2], ♦ — данные по методу работы [3]; ▲ — сетка 1 (40 узлов); □ — сетка 2 (80 узлов); х — сетка 2 (80 узлов) при Т„ = 2300 К
Рис. 11. Распределение температуры:
• — данные по методу работы [2]; ▲ — сетка 1 (40 узлов); □ — сетка 2 (80 узлов); х - сетка 1 (80 узлов) для совершенного газа
Третий рассмотренный случай (см. табл. 5) соответствует обтеканию аппаратов серии ИАМ-С, на которых были получены натурные данные [11-13]. На аппаратах этой серии проводились, среди прочих, измерения профиля электронной концентрации в ударном слое с помощью гребенки Ленгмюра. Сравнение теоретической концентрации электронов с экспериментальными данными показывает их различие почти в два раза на расстоянии 8,2 калибров от носка (рис. 12). Тем не менее, порядок величины количества частиц в единице объема учитывается верно, что позволяет использовать результаты для прогнозирования радиосвязи с летательным аппаратом. Отличие объясняется приближением данного режима к "замороженному". Кроме того, как показано в работе [15], для подобных режимов значительное влияние на электронную концентрацию оказывает учет эффектов скольжения
1.Е+12 1.Е+11 ^ 1.Е+10 ^ 1.Е+09 'О) 1.Е+08 1.Е+07 1.Е+06
О 4 8 12 16
Расстояние по нормали, см
Рис. 12. Распределение электронов e-:
расчетный режим 3; Н = 70,104км; М = 27,18; калибр 8,2; • - БЛМ-С [11-13]; ■ — неравновесный расчет
газа в граничных условиях на поверхности тела. Эффекты скольжения значительно изменяют температурный профиль и, соответственно, концентрацию химических компонентов. Различие электронной концентрации, с учетом и без учета скольжения на поверхности тела, достигает двух порядков. Особенностью данного расчета является необходимость использования полных уравнений Навье-Стокса. Упрощенные уравнения (в приближении тонкого слоя) для одинаковых условий расчетов дают большую погрешность при определении параметров на теле, в частности тепловых потоков.
Выводы. Математические модели, основанные на уравнениях Навье-Стокса, позволяют сократить количество дорогостоящих экспериментов по определению аэродинамических, тепловых и прочностных характеристик изделий.
Анализ полученных результатов свидетельствует о возможности применения предложенного метода для расчета параметров химически неравновесного вязкого ударного слоя, включая характеристики прибортовой плазмы. Затраты машинного времени остаются приемлемыми для использования программного комплекса в инженерной практике.
СПИСОК ЛИТЕРАТУРЫ
1. Агафонов В. П., Вертушкин В. К., Гладков А. А., Полянский О. Ю. Неравновесные физико-химические процессы в аэродинамике. - М.: Машиностроение, 1972. - 344 с.
2. Власов В. И. Метод расчета вязкого ударного слоя с учетом неравновесных физико-химических процессов // Космонавтика и ракетостроение. - 1997. -№ 11.-С. 5-12.
3. В о р о н к и н В. Г. Неравновесный вязкий ударный слой на притупленных конусах // Изв. АН СССР. МЖГ. - 1979. - № 6.
4. Забарко Д. А., Котенев В. П. Численное исследование течений вязкого химически реагирующего газа около затупленных тел в рамках уравнений Навье-Стокса // Космическая наука и технология (Киев, НКАУ). - 2005. - Т. 11, № 1. - С. 36-42.
5.Кань Сань-Вук. Неравновесное ионизованное гиперзвуковое течение около затупленного тела при низких числах Рейнольдса // РТК. - 1970. - Т. 8, № 7. - C. 98-105.
6. Кокошинская Н. С., Павлов Б. М., Пасконов В. М. Численное исследование сверхзвукового обтекания тел вязким газом. - М.: Изд-во МГУ, 1980.-248 с.
7. К о т е н е в В. П., Сахаров В. И., Т и р с к и й Г. А. О расчете сверхзвукового обтекания пространственных затупленных тел химически неравновесным потоком газа // ЖВМиМФ. - 1987. - Т. 27, № 6. - C. 411-415.
8. Л у н е в В. В. Гиперзвуковая аэродинамика. - М.: Машиностроение, 1975. -327 с.
9. Андерсон Д., Таннехил Д ж., Плетчер Р. Вычислительная гидромеханика и теплообмен. В 2-х т. - М.: Мир, 1990.
10. Blottner F.G. Nonequilibrium Laminar Boundary Layer Flow of Ionized Air // AIAA J. - 1964. -V. 2, № 11. - P. 1921-1927.
11.Evans J. S., Schexnayder C. J., Huber P. W. Boundary Layer electron profiles for high-altitudes entry of a blunt slender body // AIAA J. - 1973. -V. 11, № 10. - P. 1371-1372.
12. H u b e r P. W., Evans J. S., Schexnayder C. J. Comparison of theoretical and flight-measured ionization in a blunt body re-entry flow field // AIAA J. - 1971. -V. 9, № 6. -Р. 1154-1162.
13. K a n g S. W., Jones W. Z., Dunn M. G. Theoretical and measured electron-density distributions at high altitudes // AIAA J. - 1973. - V. 11, № 2. -Р. 141-149.
14. R a k i c h J. V., Bailey H. E, Park C. Computation of Nonequilibrium, Supersonic Three-Dimensional Inviscid Flow over Blunt-Nosed Bodies // AIAA J. -1983. -V. 21, № 6. -Р. 834-841.
15. S w a m i n a t h a n S., K i m M. D., L e w i s C. H. Nonequilibrium Viscous Shock-Layer Flows over Blunt Sphere-Cones at Angle of Attack // J. of Spacecraft and Rockets. - 1983. - V. 20, № 4. - Р. 331-338.
Статья поступила в редакцию 13.05.2005
Котенев Владимир Пантелеевич окончил МГУ им. М.В. Ломоносова в 1978 г. Начальник отдела ФГУП "НПО машиностроения" (г. Реутов). Автор более 20 научных работ в области прикладной математики, численных и аналитических методов исследования течения газа при обтекании поверхностей летательных аппаратов.
V.P. Kotenev graduated from the Lomonosov Moscow State University in 1978. Chief of department of "NPO Mashinostroenie" (town Reutov Moscow region). Author of more than 20 publications in the field of applied mathematics, numerical and analytical methods of study of the gas flow around surfaces of flying vehicles.
Забарко Дмитрий Александрович окончил МГТУ им. Н.Э. Баумана в 1995 г. Ведущий инженер ФГУП "НПО машиностроения" (г. Реутов). Специализируется в области прикладной математики, численных методов и программного обеспечения аэрогазодинамических расчетов.
D.A. Zabarko graduated from the Bauman Moscow State Technical University. Leading engineer of "NPO Mashinostroenie" (town Reutov Moscow region). Specializes in the field of applied mathematics, numerical methods and software of aero- and gas-dynamic calculations.
