К вопросу о снижении плотности в набегающем потоке перед летательным аппаратом Текст научной статьи по специальности «Физика»

УДК 533.6.011

В. В. Шайдуров, Г. И. Щепановская

К ВОПРОСУ О СНИЖЕНИИ ПЛОТНОСТИ В НАБЕГАЮЩЕМ ПОТОКЕ ПЕРЕД ЛЕТАТЕЛЬНЫМ АППАРАТОМ

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

Расчет сопротивления движущихся тел в среде со сверхзвуковой скоростью является одной из основных проблем аэродинамики больших скоростей. В работах [1; 2] приведены результаты по полному сопротивлению звездообразных конфигураций в равномерном сверхзвуковом потоке с учетом всех его составляющих: волнового, вязкого и донного сопротивлений. Изучение тел с лучевой структурой показало большие возможности существенно пространственных форм в решении вопроса снижения волнового сопротивления. И наоборот, вязкая составляющая сопротивления таких тел по сравнению с эквивалентной по объему осесимметричной формой значительно больше. Вопрос об оптимальности звездообразных конфигураций по сравнению с осесимметричным эквивалентом решается с помощью многочисленных вычислительных экспериментов по параметрическому анализу [1].

В последние годы исследователями проявляется значительный интерес к изучению влияния неравномерностей распределения газодинамических параметров в потоке на сопротивление тел при их движении со сверх- и гиперзвуковыми скоростями. Этот интерес связан с тем, что иногда с помощью относительно небольших затрат энергии или вещества можно изменить структуру течения. Особое внимание уделяется решению задачи активного управления обтеканием тела посредством энергетического и силового воздействия на набегающий поток, в частности посредством локального воздействия на сверхзвуковые течения перед телом. Для технической реализации предлагается использование лазерного и СВЧ-излучения, электрического разряда. Уменьшение аэродинамического сопротивления связывается с уменьшением плотности газа в набегающем потоке, что подтверждено расчетами и непосредственными измерениями [3; 4]. Экспериментально найдено, что при использовании мощного импульсного оптического разряда перед телом (конусом, сферой), обтекаемым сверхзвуковым потоком, его аэродинамическое сопротивление уменьшается при увеличении частоты следования импульсов лазерного излучения. В работах [3; 4] нестационарные явления, возникающие в объеме газа при мгновенном локальном подводе энергии, рассматриваются как правило, на основе решения уравнений Эйлера. Для численного расчета двумерных нестационарных течений развивается метод Годунова на основе TVD- и ENO-схем. Показано значительное влияние воздействий, связанных с нелинейной природой явлений. Отмечается, что невязкий механизм перестройки сверхзвукового течения отличает бесконечная длина разреженного канала [4].

В данной статье предлагается математическое и численное моделирование нестационарного распространения

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

Постановка задачи. Рассмотрим математическую модель нестационарного распространения локального импульса энергии большой мощности в вязком теплопроводном газе [5; 6]. Для этого выпишем полную систему уравнений вязкого теплопроводного газа в декартовой системе координат для двумерного случая в безразмерном виде:

^ + — (ры) + — (ру) = 0,

ді дх ду

ды ды ды др

р +р„_ + ру_ + д!.

ді дх ду дх

ду ду ду др

р — + ры — + ру— + г

дт дт

дх

дт

ді

дх ду ду дх

де де де

р—+ры +ру +

дх ду

- = 0, ду дт

■—УУ- = 0,

ду

(1)

(2)

(3)

ді

+Р

ды ду

----+-----

дх ду

д£

ді

дд

____х_

дх

ю

дд

(4)

■У+Ф,

дУ

Р = (У- 1)ре, Д = с3е , 0,5 <ю< 1.

(5)

Здесь искомыми величинами являются плотность р(/, х, у) , компоненты вектора скорости и(/, х, у), у(/, х, у), внутренняя энергия единицы массы е(/, х, у), давление р(11, х, у) и динамический коэффициент вязкости Ц. Диссипативная функция ф определяется как положительная функция:

ч2 / \2 с \2 с \2

2 ( ду

+

Д

Ф= —

Яе

ды

дх

ду

дУ

ду ды дх ду

ды

дх

ду

ду

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

Т =-

2 Д

3

Яе

т =-Д-ху Яе

ды ду

----+-----

ду дх

ды ду

2------------

дх ду

2 Д Туу = 7 Ке

ду ды

2------------

ду дх

хх

ху

а компоненты вектора теплового потока как

у де у де

Чх =--------Ц , Чу =-------------Ц ,

Яе • Рг дх Яе • Рг ду

где Q(t, х, у) - заданное распределение импульса энергии (мощность источника); Яе,Рг - числа Рейнольдса и Прандтля соответственно; У - отношение удельных теплоемкостей.

Начальные и краевые условия. Начальные условия определяет покоящийся газ с соответствующей плотностью и температурой (и = у = 0, р = 1, е = ед при t = 0), давление и динамический коэффициент вязкости находятся согласно уравнениям (5). Мощность и положение источника в начальный момент времени соответствуют заданному режиму.

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

и ^ 0 при Я ^ гс, Я = /х 2 + у 2 .

Аналогично все искомые газодинамические величины должны стремиться к состоянию в покоящемся газе в бесконечном удалении от источника.

Будем искать приближенное решение для и,у, р, р, е, Ц на прямоугольнике [0, 7 ]х[0,1] х [0,1]. На границе Г прямоугольника О = (0,1) х (0,1) для всех t е [0, tgn ] задаются следующие граничные условия:

тххпх + тхупу = Рпх - РехГпх на (0, tfln ) х Г

ТууПу + ТхуПх = РПу - Рехна (0, tf1n ) X Г,

(6)

Ок = Ок пА, Г = Ок пГ.

Далее проведем дискретизацию двумерной задачи (1)...(5) по пространству методом конечных элементов. При этом применен метод Бубнова-Галеркина для вариационной постановки в комбинации с методом линий. В качестве пробных функций рассмотрены кусочно-билинейные функции на прямоугольниках. Для вычисления интегралов использована формула метода трапеций и ее многомерные аналоги. Этот подход дает схемы второго порядка точности по пространству.

Линеаризация вариационно-разностных уравнений. Полученные разностные уравнения в результате аппроксимации дифференциальной задачи являются нелинейными относительно искомых функций. Проведем следующую линеаризацию по времени. На (к + 1) шаге по времени зададим неизвестные функции-коэффициенты нелинейных членов их значениями с предыдущего шага. В результате уравнение (1) аппроксимируется следующим образом:

1

~рк+1,2, / - рк,2, / +

Т

1

+ 2Ь (ик+1,1Рк+1,2+1,j ик,2-1,j рк+1,2-1,j + (9)

+Ук,2,1+1рк+1,2,1+1 - Ук,2,1 -1рк+1,2,1 -1) = 0-Уравнение (2) для составляющей и вектора скорости преобразуется к виду

1 1 1/2 1/2

~Рк+1,2,1ик+1,2,1 -"ТРк+1,2,1 рк,2,1ик,2,1' +

+ ~^^((р к+1,2+1,1ик ,2+1,1 + рк+1,2,1ик,2,1 )ик+1,2+1,1 -

4к

где п(х, у) = (пх (х, у), Пу (х, у)) - вектор единичной внешней нормали к Г в точке (х, у), рех^, х, у) - заданное внешнее давление на границе расчетной области. Эти краевые условия являются неотражающими с вычислительной точки зрения, поскольку они пропускают импульс через границу Г без влияния на течение внутри области. Выполнение условий (6) определяет задание условия Неймана на внутреннюю энергию на границе расчетной области:

Уе • п = 0 на Г. (7)

Кроме того, для корректной постановки задачи предполагается, что поток на Г направлен вне области О и выполняется следующее условие:

и • п > 0 на [0, tfln]хГ. (8)

Представленные условия естественным образом вытекают из интегральных законов сохранения массы, импульса и энергии [7].

Аппроксимация дифференциальной задачи. Введем равномерную разностную сетку по времени с шагом Т = Тйп / т :

шг = : tk = кт, к = 0,..., т},

а в замыкании области О - равномерную разностную сетку с шагом к = 1/ п :

Оh = {(хг,у1): х, = 2к, 2 = 0,...,п; у1 = ;к, 1 = 0,...,п} и обозначим

(р к+1,2-1,1ик ,2-1,1 + р к+1,2,1ик ,2,1 )ик+1,2-1,1) + + ((рк+1,2,1+1Ук,2,1+1 + рк+1,2,1Ук,2,1 )ик+1,2,1+1 '

-(рк+1,2,1 -1Ук ,2,1 -1 + р к+1,2,1Ук ,2,1 )ик+1,2,1 -1 ) + _2_ зк

+ -з ;_2 (ик+1,2,1 ик+1,2-1,1 )(цк,2,1 + Цк,2-1,1 ) +

+ 3к2 (ик+1,2,1 ик+1,2+1,1 )(Цк,2,1 +Цк,2+1,1 ) +

зк (10)

+ ,, 2 (цк,2+1,1 (ук+1,2+1,1+1 - Ук+1,2+1,1-1 ) -6к

-Цк,2-1,1 (ук+1,2-1,1+1 - Ук+1,2-1,1 -1 )) +

+ ^(

2к

1)(Ц к ,2,1 +Цк ,2,1 -1) +

+ 2к2 (ик+1,2,1 ик+1,2,1+1)(цк,2,1 +Цк,2,1+1)

л 1.2 (цк,2,1+1 (

4к

1 +1 (Ук+1,2+1,1 +1 Ук+1,2-1,1+1,

-1,1+1 ) "

Цк ,2,1-1 (ук+1,2+1,1 -1 Ук+1,2-1,1 -1 )) +

+ 2к (Рк,2+1,1 - Рк,2-1,1 ) = 0.

Уравнение (3) количества движения для составляющей у вектора скорости будет выглядеть следующим образом:

Тр*+1,л+1>,. 1-Трк+и 1 р^л] + +(у*+и,]- Ук+и-и + Мк+1и+1- ик+1и)2 +

+(у+1,1+1, ]— +1,/, ] +и*+1,/, ]- и*+1,/, ]-1) +

+(ук+1,/, ] — Ук+1,/-1, ] + ик+1,/, ] - ик+1,/, ]-0 ) +

+ 4^ ((р*+1,/+1,]ик,/+1,] +рк+1,/,]ик,/,] )ук+1,/+1,] —(рк+1,/-1,]ик,/-1,] + рк+1,/,]ик,/,] )ук+1,/-1,] ) + + А 1„ ((рк+1,/,]+1Ук,/,]+1 + рк+1,/,]Ук,/,] )ук+1,/,]+1

+------2 д к,/,] ((ик+1,/+1,] - ик+1,/,] - Ук+1,/,]+1 + Ук+1,г,] )

6к

2+

,1 \\| к+1,/,]+1 к,/,]+1 1 к+1,/,] к,/,] / к+1,/,]+1 * \ 2

+ (ик+1,/+1,] - ик+1,/,] - Ук+1,/,] + Ук+1,/,]-1) +

-(рк+1,/, ]-1ук ,,■, ]-1+рк+1,/, ]ук ,,■, ] )у+1,/, ]-1) - +(ик+1,,- ]- ик+1,/-1, ]- у+1,/, ]+ук+1,,-, ]-1 )2+

+(ик+1,/,] - ик+1,/-1,] - Ук+1,/,]+1 + Ук+1,/,] ) У Линеаризация проведена таким образом, что выполняется закон сохранения массы и энергии на сеточном уровне для всей системы в целом, т. е. разностные уравнения и после линеаризации согласованы в смысле законов сохранения. Нетрудно заметить, что уравнения (9).. .(12) распадаются согласно физическим законам сохранения. Так, система уравнений (2) на (к + 1) шаге по времени разрешается относительно искомых рк+1/,] , поскольку коэффициенты ик] и у*] известны с предыдущего шага. Аналогично система алгебраических уравнений (10), (11) на (к + 1) шаге по времени разрешается относительно искомых величин ик+!/] , У^+у,] , поскольку Рк+1/] , Рк] , ик] , Ук] являются известными функциями. В системе алгебраических уравнений (12) 1 известными являются Рк+1,/,], ик+1,/,] > Ук+1,/,] > Рк,/,] , а

+^( Рк,/, ]+1 - рк I ]-1) = 0. также ик/], у*/] , динамический коэффициент вязкости

2к дк,1,] , который во всех уравнениях берется с предыду-

Уравнение (4) для внутренней энергии в свою оче- щего шага по времени. Таким образом, возможно раз-редь преобразуется к виду дельное рассмотрение системы уравнений (9).. .(12) на

1 р е -1 р е + каждом шаге по времени, что существенно сокращает

Кк+1,/,]Кк+1,/,] г к,/,]Кк,/,] т „

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

1 ^ - •(д к,. +1,] (ик+1,/+1,]+1 ик+1,/+1,]-1 )

-д к,/-1, ] (ик+1,/-1,]+1 ик+1,/-1 1, ]-1)) +

2 + Ук+и ,] - Ук+1,/-1,] )(дк,1,] + дк,/-1, ], + (11)

+2к^( Ук+1,/ ,] - Ук+1,/+1,] )(дк,/,] + дк,/+1,] , +

+3кИ Ук+1,1 , ] Ук+1,/, ]-1 )(д к,/, ] + дк ,/,]-1, +

2 + Ук+1/ ,] - Ук+1,/,]+1 )(дк,/,] + д к ,/,]+1 у +

6к -(Дк, /,]+1 (ик+1,/+1,]+1 ик+1,/-1,]+1 ) "

-д к,/, ]-1 (ик+1,/+1, ]-1 ик+1,/- 1,]-1)) +

1 / \ отдельно взятом случае.

+ 7ТГ (ек+1,/+1,]рк+1,/+1,]ик,1+1,] - ек+1,/-1,]рк+1,/-1,]ик,/-1,] ) +

2^ к+1,/+1, + ,'+1,] к,/+1,] к + ,'-1, +1,/-1,] к /-1, Сеточный аналог закона сохранения массы. Система

1 линейных алгебраических уравнений (9) на (к + 1) шаге

+ ^ (ек+1,', ]+1 р к+1,', ,/, - ек+1,', ]-1р к+1,', ]-1Ук,/, ]-1) + по времени может быть записана в следующем виде (опу-

стим индекс к):

1 +

2^к,1,] ч-^1^1,] 1,^1,] ■ ■ ^1,.,^1 ■ ^1,,ар]р/,] + ар+1,]р/+1,] + ар-1,]ру_1,] +

V1', ]+1р +а', ]-1р = Р/,]

Рк,/, (ик+1,/+1, - ик+1,/-1,] + Ук+1,/, ]+1 - Ук+1,

+ ^ уек+1,/, ] - ек+1,/-1,] )(д%к, I, ] +д к,/-1, ] )+

+ ^ уек+1,/, ] - ек+1,/+1,] )(дк, /, ] + к,/+1, ] )+

+ ек+1,/, ] - ек+1,/,]-1 )(д%к, /, ] + 1% к,/, ]-1 )+

ек+1,/, ] - ек+1,/,]+1 )(дк, /, ] + к,/, ]+1 )=

= б'- 1 + Зк2 дк,/,] (ик+1,/+1,] ' - ик+1,/,] )2 +

(13)

+ар,]+1р/,]+1 +ар,1р/,]-1 = Fр', ,

V / = 1,..., п -1, V ] = 1,..., п-1.

Здесь р/,] - искомые значения плотности в узлах сетки, а - коэффициенты при неизвестных.

Запишем уравнение (13) в матричном виде:

Лкр к = ^, (14)

лк

где Л - сильно разреженная матрица, имеющая отличные от нуля элементы на пяти диагоналях, включая главную диагональ.

Сеточный аналог закона сохранения импульса. Сис-+ 1 д (и .. )2 + тема линейных алгебраических уравнений (10), (11) мо-

о/2 Дк,/,]\ик+1,/,] ик+1,/-1,]) г-

3к ^ ^ жет быть записана в следующем виде:

+д (У - У )2 + (12) < 1ии] + а1+1']и/+1,]+<х’4-1,] +

3к^ к,/,] +1,/,]+1 к+1,/,] / ~ У У ^

+ аи]+1и/,]+1 + а«] 1и/,]-1 + (1 5)

( )2 + (15)

о .2 Дк,/, ] (Ук+1,/, ] Ук+1,/, ]-1 ) +0!+1, ]+1 + в /+1, ]-1 +

3к +1,]+^ У*и У1+1,]-1

+ 4^ Дк’1’] ((Ук+1>/+1>] - Ук+Ш + ик+1,/,у+1 - ик+1,/,] )2 + +Ри_1']+1 У/-1,]+1 + РГ1,]-Ч-и-1 = Р'и,]

33

аУ Ч ]+аУ+1' 'у+1, ] + аГ‘' -1, ] +

+а!;]+Ч]+1+^-%-1 +

+РУ+1,1+4+1,]+1+Г’14+1,]-1 + (16)

+РГи+4-1,1+1+РГи-4-1,1-1 = К -

V / = 1,..., п-1, V ] = 1,..., п -1,

где и/, и У/, - искомые компоненты скорости в узлах

сетки; а и в - коэффициенты при неизвестных.

Особенностью уравнений (15), (16) является то, что их решение необходимо искать одновременно для обеих компонент скорости и и У . В данном случае матричный вид системы (15), (16) будет следующий:

Лкик + бУ = Кк, вкик + ЛУ = Кук, (17)

где Лик, Лук - пятидиагональные матрицы; Бгк, Б^ - четырехдиагональные матрицы. Здесь вдоль главной диагонали матрицы располагаются коэффициенты аи и аУ , а вдоль главной диагонали верхнего правого и нижнего левого квадрантов - коэффициенты ви и .

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

а',] е + а'+1,]е + а'-1,]е +

^е ,] Т ^е +1,] Т ^е -1,] Т

+ае-]+1е<,]+1 + аГЧ,1.1 = К", (18)

V / = 1,..., п -1, V ] = 1,..., п -1.

Во всех случаях уравнения для граничных точек и узлов, т. е. для таких точек, у которых / = 0, / = п, = 0 и

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

неизвестных, чьи индексы выходят за границу интервала [0, и], обнуляются.

Система линейных алгебраических уравнений (18) может быть записана также в виде

Лекек = Кек, (19)

к

где Ле - также матрица, имеющая отличные от нуля элементы на пяти диагоналях, включая главную диагональ.

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

хт+1 = хт - Б 1 (Л^хт - X

где Лкх = Кк - исходная система уравнений; Б_1 - итерационный множитель, совпадающий с соответствующим элементом главной диагонали матрицы.

Рис. 1. Распределение плотности как функции пространственной переменной для моментов времени t = 0; 0,4; 0,8; 1,2; 1,6; 2,0

Среднее количество итераций, необходимое для получения решения при n=300, составляет в зависимости от начальных данных от 15 до 100.

Вычислительный эксперимент. Приведем расчеты газодинамических величин для характерных чисел Маха M = 4 и Рейнольдса Re = 105 [8] (рис. 1, 2).

Таким образом, отметим, что в данной работе для уравнений вязкого теплопроводного газа впервые получены граничные условия (6) и многосеточный метод, позволяющие существенно уменьшить время счета при достижении четвертого порядка точности для приближенного решения. Кроме того, полученный программный продукт используется для решения важной задачи построения воздушного разреженного коридора в набегающем потоке перед летательным аппаратом путем использования периодических источников энергии большой мощности.

Работа выполнена при финансовой поддержке программы «Университеты России» (проект № УР.03.01.101).

Библиографический список

1. Щепановский, В. А. Газодинамическое конструирование в проблеме сопротивления / В. А. Щепановский, Г. И. Щепановская // Газодинамическое конструирование в проблеме сопротивления. Новосибирск : Наука. Сиб. предприятие РАН, 1998. 216 с.

2. Щепановский, В. А. Вычислительное моделирование воздушно-космических систем / В. А. Щепановский, Г. И. Щепановская // Модели, методы, технологии. Т. 1. Новосибирск : Сиб. изд. фирма «Наука» РАН, 2000. 232 с.

3. Коротаева, Т. А. Пространственное сверхзвуковое обтекание заостренного тела при подводе энергии перед ним / Т. А. Коротаева, В. М. Фомин, А. П. Шашкин // ПМТФ. 1998. Т. 39, №» 5. С. 116-121.

4. Зудов, В. Н. Сверхзвуковое обтекание импульсного источника энергии большой мощности / В. Н. Зудов ; ИТПМ. Новосибирск, 2000. 5 с.

5. Шайдуров, В. В. Математическое моделирование нестационарного распространения импульса энергии большой мощности в вязком теплопроводном газе / В. В. Шайдуров, Г. И. Щепановская // Вычислительные технологии. 2001. Т. 6. Ч. 2. Спец. вып. С. 693-698.

6. Shaidurov, V. V. Solution to viscous heat-conductive gas equations based on multiprocessor computer system /

Рис. 2. Распределение внутренней энергии как функции пространственной переменной для моментов времени t = 0; 0,4; 0,8; 1,2; 1,6; 2,0

V. V. Shaidurov, G. I. Shchepanovskaya // Proceedings of the International Conference on Computational Mathematics. Novosibirsk : ICM & MG Publisher, 2002. Pt. 1. P. 83-87.

7. Шайдуров, В. В. Математическое и численное моделирование нестационарного распространения импульса энергии большой мощности в вязком теплопроводном газе. Ч. 1. Математическая постановка задачи / В. В. Шайдуров, Г. И. Щепановская ; Ин-т вычисл. моде-

лирования Сиб. отд-ния Рос. акад. наук. Красноярск, 2003. S0 с. Деп. в ВИНИТИ 24.10.03, № 1860-В2003.

8. Малышев, А. В. Повышение точности приближенного решения при расчете плотности в тепловом пятне / А. В. Малышев, В. В. Шайдуров, Г. И. Щепановская // Вычисл. технологии. 2003. Т. 8. Регион. вестн. Востока. 2003. J№ 3. (Совмест. вып. Ч. 2). С. 178-190.

V. V. Shaidurov, G. I. Shchepanovskaya

THE PROBLEM OF DENSITY DECREASE FOR FILLING FLOW IN FRONT OF AIRCRAFT

Mathematical and numerical modeling of impedance decrease of aircrafts in supersonic flow is realized.