К вопросу о снижении плотности в набегающем потоке перед летательным аппаратом

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

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

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

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

свою очередь снижает сопротивление тела в целом.

Постановка задачи. Рассмотрим математическую

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

мерном виде:

1^ + — ( р и ) + — ( р V ) = о,         (1)

д t дх      ду д и д и д u дp дтхх дт ху р —+ р и — + р v— + ^---—--- = 0, (2)

д t     дх     ду  дх   дх    ду д V     д v     д v  дp  дТхУ  дт р—+ р и — + р V— + _'  —y - ^y = 0, (3)

д t     дх    ду  дy   дх    ду д e    д e    д e р—+ р и  +р v + д t     дх     ду

+Р

д и д v

--+-- д х д у

= dQ д t

д q д q

- ^ - ^ + Ф , д х   д у

Р = ( Y- 1) Р ^е , Ц = c з ^{е Ю} , 0,5 1.

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

Ф = —

Re

2 ^f д и ) ²

2 f d v

+--

\

+

3 \ д х ,

3 д у

)

д v д и

⁺г д х д у

2 д и д v

+ —

.

3( д х д у ,

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

^Т ху

Т

xx

- f -у \

• 2    ц   д и дv

  --2---.

• 3    ^Re   д х д у

ц ^f д и д v

--+— ,

^Re д у д х

^Т уу

2 ц

3Re

д v д и

\

2---

д у д х

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

γ ∂e           γ    ∂e qx =--Ц , qy =--Ц ,

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

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

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

и ^ 0 при R ^ ^ , R = д/ x ² + у ² .

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

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

^Т xx n x ^{+ Т} xy n y = Pn x ^- P ext n x ^{на (0,} t fin ) ^хГ ’

Т уу П у +Т xy n x = pn y - P ext n на (0, t fin ) хГ ,

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

V e • п = 0 на Г .                  (7)

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

м ^ п >  0 на [0, t f_ш] xГ .                (8)

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

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

tor = {tk : tk = kТ, k = 0,..., m}, а в замыкании области Q - равномерную разностную сетку с шагом h = 1/ n :

Q h = {( x i , У j ): x i = ih , i = 0,..., n ; У j = jh , j = 0,..., n } , и обозначим

Q h = Q h nQ , Г = Q h nF.

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

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

- ^р k + 1, i , j ^-р k , i , j ⁺ Т

⁺ - , ⁽ u k , i + 1, j ^р k + 1, i + 1, j

2 h

-

u k , i - 1, j ^P k + 1, i - 1, j ⁺

^{+ V}k , i , j + 1 ^P k + 1, i , j + 1

^Vk , i , j - 1 ^р k + 1, i , j - 1 )     °.

-

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

“^P k + 1, i , j u k + 1, i , j

-

1 1/2     1/2         ,

Т ^P k + 1, i , j ^P k , i , j u k , i , j ⁺

⁺ 4 h ^{(( р} k + 1, i + 1, j u k , i + 1, j ^+P k + 1, i , j u k , i , j ) uk + 1, i + 1, j

-

-

'^{( p} k + 1, i - 1, j u k , i - 1, j ^{+ P} k + 1, i , j u k , i , j ) u k + 1, i - 1, j ) ⁺

⁺ 4 h ^{(( р} k + 1, i , j + 1 ^vk , i , j + 1 ^{+ р} k + 1, i , j^vk , i , j ) u k + 1, i , j + 1

-

^{( p} k + 1, i , j - 1 ^Vk , i , j - 1 ^{+ P} k + 1, i , j^Vk , i , j ) u k + 1, i , j - 1 ) ⁺

⁺ 3 h 2 ⁽ u k + 1’ i ■ j

⁺ 3 h 2 ⁽ u k + 1, i ’ j

uk+1, i-1, j )(ц k, i, j + Ц k, i-1, j ) + uk+1,i+1,j )(цk,i,j + Цk,i+1,j ) +

⁺ 12 ^{( Ц} k , i + 1, j ^{( V}k + 1, i + 1, j + 1

6 h

-

^Vk + 1, i + 1, j - 1 )

-

^Ц k , i - 1, j ^{( v}k + 1, i - 1, j + 1

-

^Vk + 1, i - 1, j - 1 )) ⁺

⁺ - ,2 ⁽ u k + 1, i , j u k + 1, i , j - 1^{)( ц} k , i , j ^+Ц k , i , j - 1 ) ⁺ 2 h

⁺ Д7Г⁽ u k + 1, i,j - u k + 1, i , j + 1^{)( p} k , i,j ^+Ц k , i ’ j + 1 ) - 2 h

- 4 h 2 ^{( Ц} k , i , j + 1 ^{( V}k + 1, i + 1, j + 1 - ^Vk + 1, i - 1, j + 1 ) -

"^Ц k , i , j - 1 ^{( V}k + 1, i + 1, j - 1

- ^Vk + 1, i - 1, j - 1^{)) +}

⁺ :1 ⁽ P k ,i + 1,j

2 h

- ^pk , i - 1, j ) = o.

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

^{— Р} k + 1, i , j v k + 1, i , j

T

^^^^^™

TP k + 1, i , j P k , i, j v k , i , j ⁺

^{+ (} v k + 1, i , J ^- v k + 1, i - 1, J ^{+ u} k + 1, i , J + 1 ^{- u} k + 1, i , J ^{)2 +}

⁺ 4 h ^{(( p} k + 1, i + 1, j u k , i + 1, j ^{+ P} k + 1, i , j u k , i , j ) v k + 1, i + 1, j

^{+ (} v k + 1, i + 1, j    v k + 1, i , j ^{+ U}k + 1, i , j   ^Uk + 1, i , j - 1 ) ⁺

^^^^^™

—I

"^{( p} k + 1, i - 1, j u k , i - 1, j ^{+ P} k + 1, i , j u k , i , j ) v k + 1, i - 1, j ) ⁺

⁺ 4 h ⁽⁽ P k + 1, i , j + 1 v k , i , j + 1 ^{+ P} k + 1, i , j v k , i , j ) v k + 1, i , j + 1

^^^^^™

—I

'^{( p} k + 1, i , j - 1 v k , i , j - 1 ^{+ P} k + 1, i , j v k , i , j ) v k + 1, i , j - 1 )

^^^^^^e

^{+ ( V}k + 1, i , j    v k + 1, i - 1, j ^{+ U}k + 1, i , j   ^Uk + 1, i , j - 1 ) ) ⁺

⁺--- 3"^P k , i , J ^{(( u} k + 1, i + 1, J ^{- u} k + 1, i , J ^- v k + 1, i , J + 1 ⁺ v k + 1, i , J ^{)2 +}

6 h

^{+ ( u} k + 1, i + 1, J - ^uk + 1, i , J - v k + 1, i , J ⁺ v k + 1, i , J - 1 ^{)2 +}

^{+ ( u}k + 1, i , j ^{- u}k + 1, i - 1, j ^- v k + 1, i , j ⁺ v k + 1, i , j - 1 )   ⁺

4 h^{1 ( P} k , i + 1, j ⁽ u k + 1, i + 1, j + 1

^^^^^^.

u k + 1, i + 1, j - 1 )

^^^^^^^e

^{+ ( u}k + 1, i ,j

-

^uk + 1, i - 1, j

-

v k + 1, i , j + 1 ⁺ v k + 1, i , j ) ^)-

^P k , i - 1, j ⁽ u k + 1, i - 1, j + 1

^^^^^^.

u k + 1, i -м-D ) ⁺

Линеаризация проведена таким образом, что выпол

⁺ 2 h 2 ⁽ v k + 1, i , j    vk + 1, i - 1, j ^{)( p} k , i , j ^{+ P} k , i - 1, j ) ⁺

+ X( v

2 h 2 ^{k + 1,} i , ^j

^- V k + 1, i + 1, j ⁾⁽ P k , i , j ⁺P k , i + J ⁺

+ Ал, 3h ²1 ^{k + 1,} i , J

⁺ TTT⁽ v k + 1, i ,j

^- v k + 1, i , j - 1^{)( p} k , i , j ^{+ p} k , i , j - 1 ) ^{+ -} V k + 1, i , J + 1⁾⁽ P k , i , J ⁺ P k , i , j J ⁺

⁺ 6 h 2 ^{( P} k , i , j + 1 ⁽ u k + 1, i + 1, j + 1

^^^^^^.

u k + 1, i - 1, j + 1 )

^^^^^^^e

^P k , i , j - 1 ⁽ u k + 1, i + 1, j - 1

u k + 1, i - 1, j - 1 )) ⁺

⁺ 2 h ⁽ p k , i , j + 1 Pk , i , j - 1 ) °’

Уравнение (4) для внутренней энергии в свою очередь преобразуется к виду

1          _1

— P k+1, i, jek+1, i, j — P k, i, jek, i, j + t              t

⁺ , ⁽ e k + 1, i + 1, j ^P k + 1, i + 1, j u k , i + 1, j 2 h

^^^^^™

e k + 1, i - 1, j ^P k + 1, i - 1, j u k , i - 1, j ) ⁺

⁺ 2 h ⁽ e k + 1, i , J + 1 ^P k + 1, i , J + 1 v k , i , J + 1    e k + 1, i , J - 1 ^P k + 1, i , J - 1 v k , i , J - 1 ) ⁺

⁺ 2 h p k , i , J ⁽ u k ^{+ 1}, i ^{+ 1}, J

^- u k + 1, i - 1, J ⁺ v k + 1, i , J + 1 ^- v k + 1, i , J - 1 ) ⁺

няется закон сохранения массы и энергии на сеточном уровне для всей системы в целом, т е. разностные уравнения и после линеаризации согласованы в смысле законов сохранения. Нетрудно заметить, что уравнения (9)^(12) распадаются согласно физическим законам сохранения. Так, система уравнений (2) на ( к + 1) шаге по времени разрешается относительно искомых P k + 1, i , j , поскольку коэффициенты u k , i , j и V k , i , j известны с предыдущего шага. Аналогично система алгебраических уравнений (10), (11) на ( к + 1) шаге по времени разрешается относительно искомых величин u k ₊ 1, i , j , V k ₊ 1, i , j , поскольку ^P k + 1, i , j , P k , i , j , u k , i , j , V k , i , j являются известными функциями. В системе алгебраических уравнений (12) известными являются ^P k + 1, i , j , ^u k + 1, i , j , ^V k + 1, i , j , p k , i , j , а также u k , i , j , V k , i , j , динамический коэффициент вязкости P k , i , j , который во всех уравнениях берется с предыдущего шага по времени. Таким образом, возможно раздельное рассмотрение системы уравнений (9)^(12) на каждом шаге по времени, что существенно сокращает порядок системы алгебраических уравнений в каждом отдельно взятом случае.

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

⁺₀ 7 2 ⁽ e k + 1, i , j    e k + 1, i - 1, j ^{)( p} k , i , j ^+P k , i - 1, j ) ⁺

2 h

⁺     2 ⁽ e k + 1, i , j ^- e k + 1, i + 1, j ⁾⁽P- k , i , j ^+P k , i + 1, j ) ⁺

^а P ^{, j} p i , j ^+а P^{+ 1, j} p i + 1,j ^+a P^{- 1, j} p i - 1,j ^{+ +a} P ^{, j + 1} p i , j + 1 ^+a P ^{, j - 1} p i , j - 1 = F P ^,j.

2 h ^{2 (} e ^{k + 1} , i,^j

e k + 1, i , j - 1 ^{)( P} k , i , j ⁺ (^ k , i , j - 1 ) ⁺

^( e.

2 h ² V ^{k + 1,} - . ^j

^- e k + 1, i , j + 1⁾⁽i ! k , i ,j ^+p k , i , j + 1 ) =

Q i , J ⁺ T 7 2 ^P k , i , j ⁽ u k + 1, i + 1, j    u k + 1, i , j ) ⁺

3 h

⁺ - ,2 ^P k , i , j ⁽ u k + 1, i , j

3 h

^- u k + 1, i - 1, j ^{)2 +}

⁺ 3 h 2 ^ k , i , J ⁽ v k + 1, i , J + 1

-

v k + 1, i , j ^{)2 +}

⁺ .,2 ^P k , i , j ⁽⁽ v k + 1, i + 1, j 4 h

⁺ -7 7 2 ^P k , i , j ^{( V}k + 1, i , j    v k + 1, i , j - 1 ) ⁺

3 h

^- v k + 1, i , j ^{+ u}k + 1, i , j + 1

^{- u}k + 1,i , j ^{)2 +}

V i = 1,..., n - 1, V j = 1,..., n - 1.

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

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

A h Р h = F p ,                    (14)

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

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

^а u ^J u - , j ^{+ а} u⁺¹’ ^J u - + 1, j ^+а u^-1’ ^J u - - 1, j ⁺

^+a u^{J + 1 u}i , j + 1 ^{+ а} u^{J 1 u}i , j - 1 ⁺                   (1₅)

+B i ^{+ 1} , ^{j +1}v        + B i ^{+ 1} , ^{j -1}v       +

^{+ e} u      ^Vi + 1, j + 1 ^{+ e} u     ^Vi + 1, j - 1 ⁺

+B i ^- *, J ⁺¹v        + B i ^- *, J ^-1v        = F i ,J

^{+ e} u      ^Vi - 1, j + 1 ^{+ e} u      ^Vi - 1, j - 1 F u ,

^а *, j v , , j ^{+ а} v⁺¹’ j V i + 1, j ^{+ а} t"^1, j v - - 1, j ⁺

^{+ a} v j ^{+ 1} V i , j + 1 ^+a v j ^{- 1} Г , j - 1 ⁺

^+e v^{1J+ 1} U i + 1,j + 1 ^{+ в} ^’ j - U i + 1,j - 1 ⁺             (16)

+в Г^Ь j ^{+ 1} u . - 1, j + 1 + в^ j ^{- 1} u . - 1, j - 1 = F t ’ j ,

V i = 1,..., n - 1, V j = 1,..., n - 1 , где u i , j и v i , j - искомые компоненты скорости в узлах сетки; а и в - коэффициенты при неизвестных.

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

A ^h u h + B ^ v h = F ^h , B_v h u h + A_v h v h = F h , (17) где A h , A ^ - пятидиагональные матрицы; B ^h , B ^, - четырехдиагональные матрицы. Здесь вдоль главной диагонали матрицы располагаются коэффициенты а _u и а _v , а вдоль главной диагонали верхнего правого и нижнего левого квадрантов - коэффициенты в _u и в _v .

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

^а e e , j ^+а e^+1, j , + 1, j ^+а i - ^1, e - 1, j +

^+а e^{j + 1} e , j + 1 ^+а e^{j - 1} e , j - 1 = F ^{, j} ,             ⁽¹⁸⁾

V i = 1,..., n - 1, V j = 1,..., n - 1.

Во всех случаях уравнения для граничных точек и узлов, т. е. для таких точек, у которых i = 0, i = n , j = 0 и j = n , также входят в систему, но коэффициенты при неизвестных, чьи индексы выходят за границу интервала [0, п ], обнуляются.

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

A h e h = F ,                 (19)

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

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

x m + 1 = x m ^- D ^{1 (} A h x m ^- F h ), где A h x = F h - исходная система уравнений; D ^{- 1} -итерационный множитель, совпадающий с соответствующим элементом главной диагонали матрицы.

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

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

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

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