Численное решение сопряженной задачи аэрогазодинамики и внутреннего теплопереноса в конструкциях гиперзвуковых летательных аппаратов Текст научной статьи по специальности «Математика»

УДК 519.688

Ю. И. Д и м и т р и е н к о, А. А. З а х а р о в, М. Н. К о-р я к о в, Е. К. С ы з д ы к о в

ЧИСЛЕННОЕ РЕШЕНИЕ СОПРЯЖЕННОЙ ЗАДАЧИ АЭРОГАЗОДИНАМИКИ И ВНУТРЕННЕГО ТЕПЛОПЕРЕНОСА В КОНСТРУКЦИЯХ ГИПЕРЗВУКОВЫХ ЛЕТАТЕЛЬНЫХ АППАРАТОВ

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

E-mail: dimit@serv.bmstu.ru

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

Наличие автоматизированных программных комплексов для расчета аэродинамики конструкций гиперзвуковых летательных аппаратов (ЛА) с учетом теплообмена в широком диапазоне изменения геометрии аппарата и режимов полета дает возможность проводить широкомасштабное математическое моделирование реальных прототипов изделий. В настоящее время существует значительное число коммерческих программных комплексов (ANSYS CFX, FLUENT, STAR CD, FlowVision и др.), предназначенных для моделирования аэродинамики и теплообмена ЛА, однако каждый из них реализует один или несколько численных методов решения уравнений газодинамики и тот или иной алгоритм расчета теплообмена в окрестности поверхности ЛА. В связи с этим в настоящее время продолжают интенсивно разрабатываться новые вычислительные технологии в области сверх- и гиперзвуковой аэродинамики.

Чаще всего [1-4] при численных расчетах аэродинамики ЛА температуру на поверхности тела находят либо из условия «холодной

стенки», когда температура или тепловой поток на поверхности являются заданными величинами; либо из условия теплоизолированной (адиабатической) стенки, когда предполагают отсутствие обмена теплоты между газом и стенкой. Прямые методы совместного решения задачи газодинамики и теплообмена применительно к конструкции также существуют, однако, как правило, они неэффективны с вычислительной точки зрения. Дело в том, что в такой сопряженной задаче аэродинамики и внутреннего теплообмена существуют два различных временных масштаба: установления газового потока и прогрева стенки конструкции ЛА, которые различаются на несколько порядков. Поэтому совместное решение этих задач является весьма затратным с точки зрения машинного времени. Использование же постановок с априорно установившимся режимом течения не всегда обеспечивает заданную точность решения, особенно для ЛА сложной формы с наличием нескольких точек торможения потока. Кроме того, нелинейное граничное условие теплового баланса между газовым потоком и корпусом ЛА препятствует применению унифицированных вычислительных алгоритмов для совместного решения задач в газовой области и в твердом теле.

В настоящей работе предложен новый алгоритм сопряженного решения задачи аэротермодинамики и внутреннего теплообмена, разработано программное обеспечение для численной реализации этого алгоритма и проведена проверка этого алгоритма на модельной конструкции ЛА. Алгоритм основан на прямом численном моделировании аэротермодинамики с использованием модели 3-мерного пограничного слоя [5] и специального численного алгоритма решения уравнения теплопроводности в конструкции летательного аппарата. Данный подход не требует значительных вычислительных ресурсов и постоянного повышения эффективности вычислительных технологий генерации адаптивных сеток [6]. В предложенном подходе используются конечно-разностные схемы высокого порядка точности с малой схемной диффузией [1, 7-12].

Математическая постановка сопряженной задачи аэротермодинамики и внутреннего теплообмена. Рассмотрим носовую часть конструкции ЛА, обтекаемую гиперзвуковым газовым потоком. Будем рассматривать три характерные области: У1 - область высокоскоростного течения идеального нетеплопроводного газового потока; У2 - область пограничного слоя, в котором решаются полные динамические уравнения Навье - Стокса для теплопроводного газа; У3 - область, соответствующая стенке конструкции ЛА.

В области У1 имеет место система уравнений Эйлера, описывающая поведение идеального нетеплопроводного газа:

P + Vp = 0, dt

Р + У(ру ®v + рЕ) = 0, (1)

^ + 4^ + р) ) = 0, где р - плотность газа; I - время; V - вектор скорости; р - давление

( р Л

р = р —0 в ; Я0 - универсальная газовая постоянная; ц - молекуляр-V )

ная масса газа; Е - метрический тензор; е - плотность полной энергии

I I2

V

газа (а = сув + —*—, здесь су - теплоемкость при постоянном объеме,

2 . д в - температура газа, |VI2 = V • V - квадрат модуля скорости); У = гг-- -

дХ'

набла-оператор в криволинейных (адаптивных) координатах X' [7]. Граничные условия к этой системе имеют вид

V • п = 0 (2)

на жесткой стенке, где п - вектор внешней нормали к поверхности;

Р = Р^ = ^ Р = Р^ (3)

на сверхзвуковой входной границе, где рж, рж - заданные значения плотности и давления, vtx - заданный вектор скорости набегающего потока;

|Р = о, V • п = 0, ^ = 0, дг = 0 (4)

дп дп дп

на границе симметрии, где V = V • т 1 - касательные к границе симметрии составляющие скоростей газа, т1 - единичные векторы в касательной плоскости, = 1, 2; на сверхзвуковой выходной границе условия не задаются.

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

I = 0: р(0, х) = Р0(х), v(0, х) = Vо(x), р(0, х) = р0(х), (5)

где р0(х), р0(х) - заданные распределения плотности и давления, v0(x) - заданное распределение вектора скорости, х - радиус-вектор точки области течения газа.

В области пограничного слоя У2 рассматриваются уравнения Навье - Стокса, описывающие движение вязкого теплопроводного газа:

dp + Vpv = 0, dt

dpv

~еГ

dps

~dt

+ V(pv ®v + pE-Tv ) = 0, + V((ps + p) v -Tv -v + q) = 0.

К этой системе добавляются соотношения для тензора вязких напряжений и вектора теплового потока:

Т = ^(Уу)Е + /л2(у®\ + У®ут),

Я = -XV в,

где цх и /л2 - коэффициенты вязкости, X - коэффициент теплопроводности газа. Для расчета коэффициента вязкости ¡л2 применяли формулу из работы [10]:

/и2 =ßo

/ \ 0.5 Г 0Л

V0 J

1Л

0

1+00 0

0 J

в0 = 288,15 K, вх = 111 K, = 1,7894-10-5 Па • с.

Газ предполагался вязко-несжимаемым (объемно невязким), по-2

этому // =

Граничные и начальные условия к системе (6) имеют вид

v = 0, -XV в • n = qe на жесткой стенке, где qe - тепловой поток к стенке;

Р = Р^ v = ТО5 ^ X ад

на сверхзвуковой входной границе;

др дп

= 0, v ■ n = 0,

dvT

_Ti

дп

= 0,

dp дп

= 0

на границе симметрии;

dv = 0, dp = 0

дп

дп

на сверхзвуковой выходной границе.

На внешней поверхности пограничного слоя Ъ являющейся границей раздела движения идеального и вязкого потоков, формулируются следующие условия непрерывности:

[р] = 0, [V] = 0, п гл = 0, [в] = 0. (7)

Здесь [в] - скачок функций на внешней поверхности пограничного слоя.

Начальные условия к системе (6) те же, что и к системе (1).

В области У3 стенки ЛА рассматривается уравнение теплопроводности

дв

РС дв = К дв, (8)

дг

где Хя - коэффициент теплопроводности в твердом теле, р - плотность материала, св - удельная теплоемкость, Д - оператор Лапласа.

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

£: X Ув .п = ХУв .п + еав4 - е ав4, [в] = 0, (9)

ю я я g g е я w, 1--> > V ^

где в№ - температура твердой стенки (совпадает с температурой газа на этой стенке), ве - температура внешней поверхности пограничного слоя, Увя - градиент температуры на твердой стенке со стороны конструкции, Ув^ - градиент температуры со стороны пограничного слоя газа, е и - интегральные коэффициенты излучения нагретого газа и твердой поверхности, а - коэффициент Стефана - Больцмана. Физико-химические превращения материала стенки (унос, плавление, термодеструкция) не учитывались.

На внутренней поверхности конструкции ЛА выполняется условие теплоизоляции:

У в. п = 0. (10)

Начальное условие для уравнения (8) имеет вид

г = 0: в = в0. (11)

Разработка метода решения сопряженной задачи. Для решения сформулированной ранее сопряженной задачи предложен следующий метод: вводится цикл по «медленному» времени г = Х/Х0, соответствующему процессу распространения тепла в стенке конструкции, где г0 - характерное время нагрева конструкции. Внутри этого цикла вводится «быстрое» время т = Х/Ху, где гу - характерное время установления газового потока. В связи с тем, что для каждого фиксированного момента медленного времени ti тепловой поток на жесткой стенке ^ = ХУвя. п неизвестен, его полагают фиксированным; тогда для уравнений газовой динамики на жесткой стенке из уравнений (9) рассматривается только первое:

Ъ : ХУв .п = а + еав4 - е ив4,

ю g * s s V g е '

система уравнений газодинамики и (6) отделяется от уравнения теплопроводности (8) на одном шаге медленного времени.

Согласно модели 3-мерного пограничного слоя [5], уравнения идеального газа (1) и вязкого газа (6) также разделяются: решение уравнений (1) ищется во всей области У1 + У2 течения газового потока с граничными условиями (2) на жесткой стенке, затем полученное решение идеального потока на жесткой стенке для плотности, касательных компонент скорости и температуры (р Уе1 = \тр ве), переносится на внешнюю поверхность пограничного слоя и вместо условий (7) формулируются следующие условия для системы уравнений (6) 3-мерного пограничного слоя:

Ъ: р = Pe, V.п = а V.т} = в = ве. (12)

Далее осуществляется решение системы уравнений (6) в области У2 по «быстрому» времени до установления с условием для теплового потока на жесткой стенке. После этого осуществляется переход к следующему моменту -+1 медленного времени.

Тепловой поток п+1 на жесткой стенке на очередном (п + 1) временном шаге рассчитывается с помощью специального метода, согласно которому сначала ищется численно-аналитическое решение уравнения теплопроводности (8) только для главных членов теплового потока . п в направлении по нормали к нагреваемой поверхности (тепловыми потоками в касательной плоскости Увт1 пренебрегают) и только со вторым граничным условием (9) (с заданной температурой поверхности). Тогда безразмерное уравнение теплопроводности (8) с граничными и начальным условиями (9) - (11) принимает вид

дв _ дв _ _ 1 = Бо—-, 0 < х < 1;

д дх (13)

— — дв — х = 0: в = в*; х = 1: — = 0; Т = 0: в = 1,

дх

- в

где Бо = —^^ - параметр Фурье, Н - толщина оболочки, в = - -

РСН - в Х в0

безразмерная температура, в = --, х = - - безразмерная координата

* в0 Н

по толщине оболочки. В силу линейности задачи (13) ее решение -

- дд безразмерная температура в(-, 1 ) и тепловой поток а = — (х, 1),

дх

являются линейными функциями от входных данных задачи - от

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

дв

тока ^ =— (0,1) на нагреваемой поверхности в момент времени

дх

г = 1 можно представить в виде ^ = g(Fo)( в№ - 1), где g(Fo) - некоторая функция от параметра Фурье, которую находят из формулы

g (Fo) =

дв_ дх

(0,1)

(14)

в -1)

Возвращаясь к размерным величинам, для теплового потока получают формулу Ньютона:

^ (Бо)

q = a(ßw а =

H

(15)

где а - вычисляемый коэффициент теплообмена для жесткой стенки. Подставляя далее выражение (15) в граничное условие (12), замыкаем итерационный цикл по «быстрому» времени.

Для цикла по медленному времени также используют решения (14), (15), полученные при различных значениях параметра Фурье (изменение значений Бо определяется изменением значений характерного медленного времени г0).

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

дУ к

= 0, (16)

dU D¡ — + P dt

дХ1

где обозначены координатные столбцы неизвестных функций

f pVk ^

^k—1. „ ski

f P ^

U =

pv1 pv2 pv3

ps

Vk =

pv v1 + pök pvkv2 + pSk 2 pvkv3 + pSk 3 V (ps + p)

dX1

- компоненты обратной якобиевой матрицы преоб-

Здесь Р'3 =

1 дх1

разования декартовых координат х3 в адаптивные X', -' - компоненты вектора скорости в декартовом базисе -'.

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

TT —1 U2 —2 U3 —3 U4

p = Ui, V =-=2-, V =-=3 V = -=f-,

Ui

U

U

p = (Y-1)

" _ u 2 + U2 + U2

U 2

5 2U

Для решения систем дифференциальных уравнений (16) использовался численный метод ТУО 2-го порядка аппроксимации [1, 7-9]. Для его применения необходимо ввести разностную регулярную сетку, в узлах которой будут определяться конечно-разностные аппроксимации производных. Для кубической формы области задание сетки на поверхности представляет собой тривиальную задачу. Однако в более сложных случаях, когда форма области не может быть рассмотрена как кубическая, вводится специальная криволинейная (адаптивная) система координат X1, согласованная с границей рассматриваемой области [7, 11].

Рассмотрим способ получения ТУО схемы 2-го порядка аппроксимации, основанный на методе Хартена [9] («метод модифицированного потока»), который заключается в применении ТУО схемы 1-го порядка аппроксимации к соответствующим образом модифициро-

~ 1 — ванным потокам функций: / = (/ + ), где л = , g - модифицирующая добавка (будет введена ниже). Тогда конечно-разностные соотношения, аппроксимирующие уравнения Эйлера (16), можно представить в виде явной схемы [7]:

U+1 = un -At

f — 1,n

Vn+1/2

n-1/2

П+1/2

П-1/2

' n+1/2

— 1,n ^

' Vn-1/2

-At

-At

f — 2,n

Vn+1/2

V

^ 2,n

' Vn-1/2

Xi'n - Xi'n

Fn n

X2'n - X2'n

Rn n

\r3,n \r3,n XUn - Xn

— 2,n

Vn+1/2

^ 2,n

' Vn-1/2

^ 2,n

Vn+1/2

— 2,n ^

' Vn-1/2

XFn - Xn

X l,n

V

( — 3,n

Vn+1/2

1 ,n

XRn - Xn

Rn n

— 3,n

' Vn-1/2

^ 3,n

Vn+1/2

^ 3,n

' Vn-1/2

Ч - Xn

— 3,n

Vn+1/2

3,n

J

— 3,n \

Vn-1/2

XF" - Xn

1 ,n

XR,n - Xn

Rn n

2,n

4 - Xn

3,n

J

X = At/(Xn+-Xn), n+ = \

Rn

U„

i = 1;

i = 2; n- =

i = 3;

B

L >

D,

i = 1; i = 2;

i = 3;

где V п+1/2 - столбцы численного потока в 1-м направлении, которые определяются по следующей формуле:

— i ,n

V

f,n

7+1/2^П+1/2.

= !( Vкп + Vп ) р' + 1» ф

7+1/2 2 ^ ^^ 7+ ^ к,П 2 п+т к

Здесь V- = V' (Хп, г"), уп = У (Хп, г"), = (»£ + )/2, К-П" = Я' (X , гп), а Я' - матрица правых собственных векторов матрицы О' = дУ' / ди:

R1

1

R2

v1 - с

—2 v

—3

v

H - cv1 1

—1 v

v2 - с

—3

v

1

—1 v

—2 v

—3

v

2 1

—1 v

—2 v

—3

v

0 0 -1 0

0 0 0 1

—2 —3

-v v

—2

R3

H - cv

—1 v

—2 v

v3 - с

H - cv2

0 1 0 0

—1 v

0 0 0 -1

—3

-v

1

v1 + с

—2 v

—3

v

H + cv1 1

—1 v

v2 + с

—3

v

—2

—1 v

—2 v

—3

v

0 0

-1 0

0 1

0 0

—1 —2 v v

H + ov

—1 v

—2 v

v3 + с

H + су3

где с - локальная скорость звука, Н - полная энтальпия.

Вектор численной диссипации ФП+"/2 имеет следующие компоненты [1, 7]:

p

s,n+H2

= 1 Ж,п+1/2 ) • (j + gsп+ ) - Ж,n+1/2 + Г1 ,п+1/2 №

2 1 4 s,П + 1/ 2S N«-7 s, j ' о s,n +

gSs ,n = minmodK ,„+1/2);

min mod(x, y) = 0,5[sign(x) + sign(y)] min(| x |,| y |);

Ж )( g s ,n+

s ,7+1/2/ s,n+1/2'

Y

s,n+1/2

2a

s,n+1/2

0,

^,П+1/2 * 0,

^,П+1/2 * 0,

где +1/2 = а1п+У2Р\п+1/2, а <П = «1 (ХП> ) - с0бственные значения матрицы О, которые в декартовой системе координат имеют вид

«1 = у1 — с; «2 = у1; « = у1; @'4 = у1; « = у1 + с. Величины агкпц+1/2, представляющие перепады вектора газодинамических переменных Ц в характеристической форме, определяются из следующих соотношений:

а- = (а;;+1/2,...,а5;:+1/2)г = 2)"Щ+ — Ц).

Здесь (К£ш)—1 = 0,5(( )—1 + (Я-)—1), а (К^и)—1 = (К1')—^, ^) - матрицы левых собственных векторов матрицы О1, которые имеют следующий

(R1)-1 =

(R2)-1 =

(R3)-1 =

вид:

(Y- 1)|v|2 +2 (v1 (1 -Y) v1 - с (1 -у) v2 (1 -Y) v3 Y -1

4с2 2с2 2с2 2с2 2с2

2 с2 - (Y- 1)|v|2 (Y -1)v1 (Y -1)v2 (Y -1)v3 Y -1

4с2 с2 с2 с2 с2

—2 v 0 -1 0 0

—3 v 0 0 1 0

(y - 1)|v|2 -2 с¥1 (1 -Y) v1 + с (1 -y) v2 (1 -Y) v3 Y -1

4с2 2с2 2с2 2с2 2с2 .

(y- 1) 1 v |2 +2 ov2 (1 -y) v1 (1 -y)v2 - с (1 -y) v3 y -1

4с2 2с2 2с2 2с2 2с2

2 с2 - (y- 1) 1 v |2 (y-1) v1 (y-1) v2 (y-1) v3 y-1

4с2 с2 с2 с2 с2

—1 -v 1 0 0 0

—3 v 0 0 -1 0

(y- 1) 1 v |2 -2 су2 (1 -y) v1 (1 -у) v2 + с (1 -y) v3 y-1

4с2 2с2 2с2 2с2 2с2

(y- 1) 1 v |2 +2 dv3 (1 -y) v1 (1 -y) v2 (1 -y) v3 - с y-1

4с2 2с2 2с2 2с2 2с2

2 с2 - (y- 1) 1 v |2 (y-1) v1 (y-1) v2 " (y-1) v3 y-1

4с2 с2 с2 с2 с2

v1 -1 0 0 0

—3 -v 0 1 0 0

(y- 1) 1 v |2 -2 ^3 (1 -y) v1 (1 -y) v2 (1 ) v3 + с y-1

4с2

2с2

2с2

2с2

2с2

Функция численной вязкости у определяется по формуле

П 2 |, 2 >8,

^(2) = Ь 2 2ч /о

[(2 + 8 ) / 28, Z <8,

где £ - параметр диссипации. При £ = 0 диссипация будет наименьшей, чем больше £, тем более диссипативной является разностная схема.

Численная аппроксимация граничных условий (2)-(4) осуществлялась на основе метода фиктивных ячеек, который, как и используемые разностные схемы, обеспечивает второй порядок точности. Для расчетной области V вводились дополнительные внешние слои фиктивных ячеек и разностные операторы для граничных узлов (рис. 1). Таблица поясняет аппроксимацию граничных условий в фиктивных узлах.

Разностные операторы для граничных условий

Рис. 1. Добавление фиктивных ячеек

Жесткая стенка PZ(n)' n vtiZ(n^ 1 1 2' vn n vnZ (n)' Pn PZ(n)

(vn < 0, jvj > a): p = 2p — p * vj = 2vj — vj * ' n 'н ' Z(nY n n Z(n)' Pn = 2Рн — pm

(vn > 0, jvj > a): pn = °Z(n)' vn = vJZ(n); p n ~ pZ(n)

V > 0, jvj < a): Pn = Pzw vn = vW Pn = 2Pg — PZ(n)

Граница симметрии 3 h - h Y< hA(11)+ = 0; Л y 1 Yl i=1 XA(n)+ A(n)- vnn vnZ(n)' h vtj, p};i I = 1, 2

Численный метод решения задачи газодинамики для вязкого газа. После решения во всей области V1 + V2 задачи газодинамики для идеального газа, в области V2 пограничного слоя решается задача для вязкого газа. Для этого систему (6) записывают в полных дифференциалах:

р

d (1 ^

dt

= Vv,

\Pj

P^v = -Vp + VTv, dt

de т

р- = -pVv -Vq + TvV® vr, dt

где е - внутренняя энергия газа.

Затем записывают систему (17) в адаптивных координатах в недивергентной форме:

р

1

+ v1Pk

Vdt P

Птгк

д

dXk

Птгк

(- и

= P

k dv1

dXk ' „ д (

dvk dvk dp

p^ + p^P" — + P"^- = P -

dt 1 dX" 1 dX" 1 dX"

pde+PviP"^e_+pP" =p"A_ dt 1 dX" 1 dX" 1 dX"

dv^ dv

+Mpqk,P,sPm —--q-,

1 k p dXs dXm

M1k],Ps-

v 1 dX J

Äpk d6 ^ ■

(18)

dXk

гдеМтрр = цх6т61 + /л2(ёт16№ + 6тр6]р), 61 - символы Кронеккера.

Для решения системы (18) используют метод расщепления по физическим процессам. На первом шаге интегрируют систему без учета вязких членов (т. е. систему уравнений Эйлера) методом ТУО. На втором шаге учитывают вязкие члены, не рассматривая конвективные и уравнение неразрывности. Для второго шага система будет выглядеть следующим образом:

dv

р— = P" dt

d

de D P— = P 1 dt 1 dX

dX" d (

M1kkps

dv1 ^

1 dXs

ÄP

k de

- t dv1 dvq

—-i + Mpqk,P\Pm —--^.

dXk ) 1 k p dXs dXm

(19)

Для интегрирования системы (19) используют метод дробных шагов, который состоит из четырех шагов [7, 12]. Шаг 1:

рЯ+2/6(Цп+3/6 Цп+2/6 ) = At д Цп+3/6

где Цп+ 6 - значения полученные на этапе решения системы Эйлера.

Шаг 2: Шаг 3: Шаг 4:

Здесь

pn+3/6(Un+4/6 _ иn+3/6) _ AAt л UJn+4/6

p«+4/6(jjn+5/6 üjn+4/6) _ At Д Ujn+5/6

2 3

pn+2/6(U "+1 - U n+5/6) _ At[(A1 + л2 + A3)U n+5/6 + + (A12 + A13 + A23)U n+2/6 + Fn+2/6].

U =

С —1 ^ —2

—3

vn)

F =

0 0 0

д—i д—

M pqk ps pm __

дХ* dXm )

Разностные операторы имеют следующий вид:

А.=ЛГл-д

дх v дх

л,=jl с л л

11 дх1 I дх1

i * j,

где А - некоторая функция.

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

Численное решение уравнения теплопроводности. Для решения одномерной задачи теплопроводности (13) используют пошаговый метод с неявной разностной схемой:

О n+1 Оп О n+1 о О n+1 I О n+1

nk -nk = р0 nk+1 - 2nk + nk-1

А ~t

П = П,

h2

nn = nn

M M-1 •

(20)

Полученная из (20) система алгебраических уравнений решается методом прогонки. Формула (14) для нахождения безразмерного коэффициента теплообмена ^Го) в конечно-разностной аппроксимации имеет следующий вид:

£ (Ро)«+1 =

П - 1)h

Результаты численного решения. На рис. 2-6 представлены результаты численного моделирования обтекания фрагмента корпуса модельного летательного аппарата гиперзвуковым потоком газа на высоте 15 км. Начальные условия и скорость набе-

0.0305 0.136 0.242 0.347 0.453 0.559 0.664 0.77 0.875 0.981 1.09

Рис. 2. Распределение плотности газового потока в окрестности поверхности ГЛА (области У1 + К2), кг/м3

■0.0376 184 368 553 737 921 1.11е+03 1.29е+03 1.47е+03 1.669+03 1.849+03

Рис. 3. Распределение продольной компоненты скорости газового потока в окрестности поверхности ГЛА, м/с

Рис. 4. Распределение давления газового потока в окрестности поверхности ГЛА, Па

temperature

170 357 544 731 918 1.11е+03 1.29е+03 1.48е+03 1.67е+03 1.85е+03 2.04е+03

Рис. 5. Распределение температуры газового потока в окрестности поверхности ГЛА для адиабатической стенки, K

170 357 544 731 918 1.11 е+03 1.29е+03 1.48е+03 1.67е+03 1.85е+03 2.04е+03

Рис. 6. Распределение температуры газового потока в окрестности поверхности ГЛА для стенки с учетом теплообмена, К

гающего потока (граничное условие на входной границе) задавались равными друг другу: р = 0,195 кг/м3; ух = = 0 м/с; у2 = 1800 м/с; р = 12346 Па. Значения рамерных параметров задачи принимались следующими: в0 = 293 К; е1 = 0,8 и е2 = 0,3; ^ = 0,3 Вт/мК; рз = 1800 кг/м3;

= 0,8 кДж/кгК.

Для выявления особенностей течения в наиболее ответственных зонах необходимо использовать мелкие сетки, что приводит к существенному увеличению времени счета. Проведем оценку вычислительной сложности описанного алгоритма для решения задачи газовой динамики. Так, если размер сетки - 1,5106, то число операций для решения уравнений Навье - Стокса - 65 103, а число временных слоев до установления - 1104. В итоге получаем, что необходимо выполнить около 1015 арифметических операций.

Численное решение с хорошим разрешением передает головной скачок уплотнения в критической точке ЛА, максимум плотности, давления и температуры приходится на критическую точку, в которой температура достигает 2000 К. На рис. 5 и 6 показано сравнение температуры при адиабатической стенке и с учетом теплообмена между газовой средой и оболочкой в момент времени 50 с после начала движения. Во втором случае температура в среднем по поверхности тела получается на 25 % ниже, чем в первом, что свидетельствует о важ-

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

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

Исследования выполнены при частичной поддержке грантов Президента РФ МК-2498.2011.8, МК-3150.2012.8. Результаты расчета были получены с использованием суперкомпьютера СКИФ МГУ «Чебышев».

СПИСОК ЛИТЕРАТУРЫ

1. Г и л ь м а н о в А. Н. Методы адаптивных сеток в задачах газовой динамики. М.: Физматлит, 2000. - 248 с.

2. Ч и с л е н н о е решение многомерных задач газовой динамики / С.К. Годунов и др. М.: Наука, 1976. - 400 с.

3. С а м а р с к и й А. А., П о п о в Ю. П. Разностные методы решения задач газовой динамики. М.: Едиториал УРСС, 2004. - 424 с.

4. Р о у ч П. Вычислительная гидромеханика. М.: Мир, 1980. - 612 с.

5. Д и м и т р и е н к о Ю. И., З а х а р о в А. А., К о р я к о в М. Н. Модель трехмерного пограничного слоя и ее численный анализ // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. Спец. вып. 2011. - С. 136-149.

6. Т е х н о л о г и я построения разностных сеток / В.Д. Лисейкин и др. Новосибирск: Наука, 2009. - 414 с.

7. Д и м и т р и е н к о Ю. И., К о т е н е в В. П., З а х а р о в А. А. Метод ленточных адаптивных сеток для численного моделирования в газовой динамике. М.: Физматлит, 2011. - 286 с.

8. Р а з в и т и е метода ленточно-адаптивных сеток на основе схем TVD для решения задач газовой динамики / Ю.И Димитриенко и др. // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2011. - № 2. - С. 87-97.

9. H a r t e n A. High Resolution Schemes for Hyperbolic Conservation Laws // J. Comp. Phys. - 1983. - Vol. 49. - P. 357-393.

10. К р а с н о в Н. Ф. Аэродинамика. В 2 т. М.: Высш. школа, 1980. - Т.1. - 495 с.; Т. 2. - 416 с.

11. Р а з р а б о т к а программного обеспечения для математического моделирования в задачах сверхзвуковой аэрогазодинамики перспективных летательных аппаратов / Димитриенко Ю.И. и др. // Супервычисления и математическое моделирование: Тр. Междунар. сем. Саров, 2010. - С. 148-155.

12. Я н е н к о Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, сиб. отд., 1967. - 197 с.

Статья поступила в редакцию 30.05.2012