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



УДК 519.688

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

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

Используемый в настоящее время метод расчета температуры поверхности для гиперзвуковых летательных аппаратов (ГЛА) может давать значительные погрешности, поскольку температура аэродинамического нагрева конструкций может быть настолько высока, что потребуется применение теплозащитных материалов, вбрасывающих в набегающий поток смесь многокомпонентных газов. В статье предложен новый численный метод решения полной сопряженной задачи аэрогазодинамики и внутреннего теплообмена в конструкциях перспективных ГЛА, основанный на итерационном решении трех типов самостоятельных задач: задачи газодинамики для идеального газа, задачи динамики вязкого газа в рамках полных динамических уравнений Навье—Стокса для трехмерного пограничного слоя, уравнения теплопроводности для оболочки ГЛА. Предложены алгоритмы численного решения этих задач в криволинейных неортогональных координатах. Выполнена модификация конечно-разностных схем типа TVD для их применения на геометрически-адаптивных сетках. Представлены результаты моделирования обтекания фрагмента носовой части модельного ГЛА эллипсоидальной формы. Проведено сравнение результатов расчета температуры для случаев адиабатической стенки и с учетом теплообмена между газом и стенкой, показавшее важность учета теплообмена при расчете температуры поверхности ГЛА.

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

Modeling of coupled aerogasdynamics and heat transfer processes on the thermal protection surface of a future hypersonic aircraft

Yu.I. Dimitrienko, A.A. Zakharov, M.N. Koryakov, E.K. Syzdykov

To create future high-speed hypersonic aircrafts, it is necessary to predict with sufficient accuracy heat flows and temperature fields on the surface of the

Исследования выполнены при поддержки гранта Президента РФ МК-3218.2013.8.

ДИМИТРИЕНКО Юрий Иванович

(МГТУ им. Н.Э. Баумана)

DIMITRIENKO Yuriy Ivanovich

(Moscow, Russian Federation, Bauman Moscow State Technical University)

ЗАХАРОВ Андрей Алексеевич (МГТУ им. Н.Э. Баумана)

ZAKHAROV Andrey Alekseevich

(Moscow, Russian Federation, Bauman Moscow State Technical University)

КОРЯКОВ Михаил Николаевич (МГТУ им. Н.Э. Баумана)

KORYAKOV Mikhail Nikolaevich

(Moscow, Russian Federation, Bauman Moscow State Technical University)

СЫЗДЫКОВ Елтуган Кимашевич

(ОАО «ГосМКБ «Радуга» им. А.Я. Березняка»)

SYZDYKOV Eltugan Kimashevich

(Dubna, Moscow Region, Russian Federation, JSC «Raduga» n. a. A.Ya. Bereznyak)

*

aircraft caused by aerodynamic heating. The surface temperature is currently calculated in two steps: the calculation of the aerogasdynamic flow parameters of a "cold" surface and the calculation of the surface temperature of the real structure. The method that is currently used for calculating the surface temperature of hypersonic vehicles can give significant errors because of the intensive aerodynamic heating of ablative beat shield materials throwing the multicomponent gas mixture into the incident flow. In this paper, a new numerical method for solving the coupled problem of aerogasdynamics and internal heat transfer in future hypersonic aircraft structures is proposed. The method implies an iterative solution of the following three independent problems: ideal gas dynamics, viscous gas dynamics on the basis of Navier-Stokes equations for a three-dimensional boundary layer, and heat conduction in a hypersonic aircraft shell. Algorithms for the numerical analysis of these problems in non-orthogonal curvilinear coordinates are developed. The TVD-type finite difference schemes were modified to be applied to geometrically adaptive grids. The results of simulation of the flow about an ellipsoidal nose segment of a model hypersonic aircraft are presented. The results of calculation of the adiabatic wall temperature in the case of heat transfer between the gas and the wall showed that the heat transfer must be taken into account when calculating the temperature distribution on the surface of a hypersonic aircraft. The presented numerical method can be applied to the analysis of aerodynamics of hypersonic aircrafts.

Keywords: aerogasdynamics, hypersonic aircraft, coupled problems, heat transfer, numerical modeling, TVD scheme, adaptive mesh, thermal protection, ablative beat shield materials.

ОЛдной из важнейших технических проблем, возникающих при создании гиперзвуковых летательных аппаратов (ГЛА), является преодоление теплового барьера: создание материалов и технологий, обеспечивающих защиту силовых конструкций ГЛА от воздействия высоких температур набегающего газового потока [1—7]. Современные теплозащитные системы пассивного типа ограничены по времени их применения при высоких температурах, активные же системы теплозащиты обладают меньшей надежностью и стоят на много дороже. В этой связи задача повышения точности расчета тепловых режимов на обтекаемых поверхностях ГЛА является важной составной частью проблемы преодоления теплового барьера.

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

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

значительных вычислительных ресурсов и позволяет более точно определять поля газодинамических параметров вблизи поверхности летательных аппаратов. В предложенном подходе используются конечно-разностные схемы высокого порядка точности с малой схемной диффузией [11 — 13].

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

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

Ы Эре

ы

+ V■рv = 0; + 0 v + pE) = 0; + 7-[(ре + р^] = 0,

(1)

где р — плотность газа; г — время; v — вектор скорости; р — давление, р = рЯ9 (Я — газовая постоянная, 9 — температура газа); E — метрический тензор; е — плотность полной энергии газа, е = е + Н / 2 (е — плотность внутренней энергии, е = сг9; сг — теплоемкость при постоянном объеме; = v - v — квадрат модуля скорости); V — набла-оператор в криволинейных

(адаптивных) координатах X' [14], У = г '

,■ Э

ЭХ'

Граничные условия для этой системы имеют вид:

♦ на жесткой стенке

v - п = 0, (2)

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

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

р = ри; v = vм; р = ри. (3)

Здесь рм, рм — заданные значения плотности и давления; vм — заданный вектор скорости набегающего потока;

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

^р = 0; у - „ = 0; ^ = 0; ^ = 0, (4)

Эп Эп Эп

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

На сверхзвуковой выходной границе условия не задаются.

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

г =0: р(0,x) = р0(x); v(0,x) = v 0 (x); р(0, x) = р0 ^). (5)

Здесь р0(x), р^) — заданные распределения плотности и давления; v ) — заданное распределение вектора скорости; x — радиус-вектор точки области течения газа.

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

Эр

+ = 0;

Эг

ЭрV + V - ^ 0 v + pE - Tv )= 0; (6) + У-[(ре + р^ — Tv - v + q] = 0.

Эг Эре Эг

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

^ = ц 1(У-v ^ + ц 2(У0 v + 70 vт); (7)

q = —АУ9, (8)

где ц 1 и ц 2 — коэффициенты вязкости; А — коэффициент теплопроводности газа. Для расчета коэффициента вязкости ц2 применялась формула из работы [15]:

ц 2 = ц

/ \0.5/

I 9 \ /

\90 /

1+9.)

9 /

4 9 \

1+ ^

9 0 /

(9)

где 90 = 288,15 К; ^ = 111 К; ц 0 = 1,7894 -10-5 Па -с.

Газ предполагался вязко-несжимаемым

2

(объемно невязким), поэтому ц 1 = — 3ц2. Граничные и начальные условия к системе (6):

♦ на жесткой стенке

v = 0, - Х70-п = , где — тепловой поток к стенке;

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

Р = Р., v = v ., р = р.;

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

^ = 0, v - п = 0,

дП

дУ,

~ др „ = 0, — = 0; дП дП

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

д'У „ др — = 0, — = 0. дП дП

(10)

(11)

(12)

(13)

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

[р] = 0; М = 0; п - Ty-т7 = 0; [0] = 0. (14)

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

Начальные условия к системе (6) имеют вид (5). В области ¥3, соответствующей конструкции ГЛА, рассматривается уравнение теплопроводности

д0

Р,С, — = X 5 А0,

' ' д1 '

(15)

где XЛ — коэффициент теплопроводности материала конструкции, который предполагается изотропным; р^ — плотность материала; ев — удельная теплоемкость; А — оператор Лапласа. Граничные условия на твердой стенке — поверхности контакта высокоскоростного потока и теплозащитной конструкции имеют следующий вид:

Е „: X, 70^ - п = Х70^ - П + в^ а04 - в, а04„,

[0] = 0, (16)

где 0 ж — температура твердой стенки (совпадает с температурой газа на этой стенке); 0е — температура внешней поверхности пограничного слоя; 70Л — градиент температуры на твердой стенке со стороны конструкции; 70 — градиент температуры со стороны пограничного слоя газа; в и в Л — интегральные коэффициенты излучения нагретого газа и твердой поверхности; а — коэффициент Стефана—Больцмана. Физико-химические превращения теплозащит-

ного материала (унос, плавление, термодеструкция) в данной работе не учитывались.

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

Е • 70-п = 0.

(17)

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

1=0: 0 = 00.

(18)

Разработка метода решения сопряженной задачи. Для решения сформулированной выше сопряженной задачи предложен следующий метод: вводится цикл по «медленному» времени 1 = 1 /10, соответствующему процессу распространения тепла в стенке конструкции (10 — характерное время нагрева конструкции). Внутри этого цикла вводится «быстрое» время т = 1 / 1у (1у — характерное время установления газового потока). Для каждого фиксированного момента «медленного» времени 1{ тепловой поток на твердой стенке ^^ = X Л 70 Л - п, вообще говоря, неизвестный, полагается фиксированным, тогда для уравнений газовой динамики на твердой стенке из уравнений (16) рассматривается только первое:

Е „: Х70е - п = д, + в, а04„ - в^ а04е (19)

и система уравнений Навье—Стокса (6) отделяется от уравнения теплопроводности (15) на одном шаге «медленного» времени.

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

Е е : Р= Ре , v - П = 0, v-т7 = У е1 , 0 = 0 е . (20)

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

Тепловой поток п+1 на твердой стенке на очередном (п + 1)-м временном шаге рассчитывается по следующему алгоритму. Согласно этому методу, сначала ищется численно-аналитическое решение уравнения теплопроводности (15) только для главных членов теплового потока У0Л • п в направлении по нормали к нагреваемой поверхности (тепловыми потоками в касательной плоскости У0 • т 1 пренебрегает-ся) и только со вторым граничным условием (16) (с заданной температурой поверхности). Тогда безразмерное уравнение теплопроводности (15) с граничными и начальным условиями (16)—(18) принимает вид

д 0 ^ д2 0 „ _ , — = Бо—0 < х <1; дг дх2

д 0

х = 0 : 0 = 0^; х =1 — = 0; I = 0: 0 = 1,

дХ

Т7 ^ s* 0

— параметр Фурье; H — толщи-

PsCsH_

на оболочки; 0 = 0/00 — безразмерная температура, 0w = 0w / 00; Х = x / H — безразмерная

координата по толщине оболочки. В силу линейности задачи (21) ее решение — безразмерная температура 0(Х,t)и тепловой поток

- д 0 - „ q = (x ,t) являются линейными функциями

дХ

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

лового потока qs = — (0,1) на нагреваемой по-дх

верхности в момент времени t =1 можно представить в виде qs = g(Fo)(0w -1), где g(Fo) — некоторая функция от параметра Фурье,

(21) (р pv1 / -к \ pv | pvkv1+ pS k 1

1, U = pv2 VI и pvkv2 + pSk2

pv3 pvkv3 + pSk3

толщи- vk (pe + p) ,

) д 0(0,1) -g(Fo)= — /(0w -1).

дХ

(22)

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

qs = а(0„ - 00); а = X^(Бо)/ Н, (23)

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

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

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

—+pí =0,

д*

дХ'

(24)

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

(25)

Pi =

дХ''

дХ1

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

матрицы преобразования декартовых координат х' в адаптивные X'; V' — компоненты вектора скорости в декартовом базисе ёг.

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

тт — 1 U2 —2 U 3 — 3 U4

p = U1; v1 = v2 v3 =

H 1 Ux Ux Ux

p = (Y-1)

U -U2 +U2 +U42

2U1

(26)

Для решения систем дифференциальных уравнений (24) использовался численный метод ТУО 2-го порядка аппроксимации. Для генерации конечно-разностной почти регулярной сетки для области Ух+ У2+ ¥3 вводилась специальная криволинейная (адаптивная) система координат X', согласованная с границей анализируемой области [11].

Рассмотрим способ получения ТУО схемы 2-го порядка аппроксимации, основанный на методе Хартена (метод модифицированного

потока), который заключается в применении ТУО схемы 1-го порядка аппроксимации к соответствующим образом модифицированным по-

~ ( 1 \ А1

токам функций: / = / + - ^ I, где X = —, а g — \ X ) АХ

модифицирующая добавка (будет введена ниже). Тогда конечно-разностные соотношения, аппроксимирующие уравнения Эйлера (24), можно представить в виде явной схемы [12]:

ип+1 = ип - А1

П п

V1,П — V1,П

т п+1/2 т П—1/2 2,п 2,П

- Хп

V1,П

Т П+1/2

. V1,П т п-1/2

X 1п _х 1 п

V1,П - V1,П ^

' п+1/2 ' п-1/2

- А1

\г2,П _Л/"2,п

' п+1/2 ' п-1/2

X ЗП - X:

\г 2 ,п _Д7 2 ,П

' п+1/2 ' п-1/2

X 1п - Х1 П

\

V 2 П

V 2 П

'Лп

- (27)

V 2 , П _V 2 , П

У п+1/2 У п-1/2

X 3

V 3 ,П

'Лп

- А1

V 3'П - V3,П

' п+1/2 ' п-1/2

X 1п - X1 П

V3 п - V3 ,П

' п+1/2 ' п-1/2

\

V3 п - V3 ,П ^

' п+1/2 ' п-1/2

2 П Л*п

2 П

'Лп

X

3 П

'Лп

где — столбцы численного потока в /-м

направлении

V' П Т п+1/2

+

= ^/ П + Vп/+П )Р1 п +

2

1

2^

И /

/' ■"'п+1/2 ^ п+1/2

Ф

/ П

(28)

где V/ П = V/ X , 1П); X/ = А1 / (Xп + - Xп);

+

^, /=1; Я, / = 2;

, / = 3;

И

п+1/2

= (И п++ + И п,П)/ 2,

Ип,П = И' (Xn ,1П), И/ — матрица правых собственных векторов матрицы О/ =дV/ / ди [12]. Вектор численной диссипации Фп+1/2 состоит из следующих компонентов [12]:

1

Ф5 ,п +1/2 2 ,п +1/2)(&5 ,п + gs ,п + )

);

^ ,п +1/2 •

2 Т V5 ,п +1/ 2 / \0 5 ,п ,п +1/2 + У 5 ,п +1/2 )а5 ,п +1/2 ;

gs ,п = тттоё(а 5 ^-1/2,а'

т1птоё(х, у ) =

= 0,5^п( х)+ 81яп( у )]тт(| х | ,| у|);

,п +1/2 ,п + - gs ,п )

(29)

У

5 ,п +1/2

2а

а

5 ,п +1/2

5 ,п +1/2

* 0;

0,

а5 ,п+1/2

,п +1/2 аз ,п +1/2Р5 ,п +1/2, а

= а'(X ,1П) —

5 ,п 5 У п ' '

Здесь а/

собственные значения матрицы С, которые в декартовой системе координат имеют вид

а ' = у/

а; а2 = у'

I'

тт'

а 3 = у';

а4' = у/; а5' = у ! + а,

где а — локальная скорость звука.

(30)

Величины

а.

^п+1/2, представляющие перепады вектора газодинамических переменных иП

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

а

п+1/2

=к:

/,П

1,п +1/2 5

* /,П \-1 +1/2

а 'П )т = • • а5,п +1/2 )

= (И п++1/2 )-1 (и

П

+

и5,п +1/2 -

■и П).

(31)

Здес1> (Ип+^Г1 = 0,5((Ип+ )-1 +(Ип,П)-1), (Ип,П)-1 = = (И' )-1 (Xn ,1П) — матрицы левых собственных векторов матрицы О', выражения для них приведены в работе [12]. Функция численной вязкости у определяется по формуле

) =

г > у;

(г2 + У2) / 2у, г < у,

(32)

где у — параметр диссипации.

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

Численный метод решения задачи газодинамики для вязкого газа. После решения во всей области У1+ У задачи газодинамики для идеаль-

п

П

П

/ П

П

П

/п

ного газа в области У2 погранслоя решается задача для вязкого газа. Для этого система (6) записывается в полных дифференциалах в адаптивных координатах [11]:

дг

р|+у'рк

д

дХк

(р))

=Рк

д-'

дХ

к '

д-к

д-к

+ р-'Р" дг р ' дХп

+ РП

др

~дХп

_ рп

дХп

д-! ) ! ' дХ

(33)

'де

р — + р-'Р! „ дг ' дХ

де пп д-' — + рРП —

дХп

= Рп

д

дХп

ХРк

де

дХк

+ м>

где Ифр _ цр8т8зр + ц2(8тз8'р + 5тр8'р) — компоненты тензора вязкости газа; 8' — символы

Кронекера; ^* = ИрфРк:Р1

д-!

д-„

к р дХ3 дХ'

— функ-

ция диссипации.

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

дU -р — = А U + F.

дг

(34)

Здесь U — координатный столбец неизвестных функций; F — координатный столбец правой части уравнений; А — матрица (4x4) дифференциальных операторов:

(0 I

0 (А1 0

-1)

U =

-2 -3

F =

0

0

/

А=

)

0 А

(35)

4 у

где дифференциальные операторы

Л _ РП

Л = р

д

дХп д

! ' дХ3

дХп

Х8 'Рк-^г ' дХк

(36)

Факторизуем диагональные операторы по трем координатным направлениям следующим образом: Лрр = Лр1 + Лр2 + Лр3 + Лр, р = 1...4, где

д

Л = Р а-

Лра Р дХ'

И® Ра — р ' дХ'

Л = Ра-

Л4 а Р' дХр

д

Х8 ''Ра-

' дХ'

], в = 1,2,3, и ; по индексам, обо-

значенным греческими буквами, суммирование не производится. Операторы Лр содержат оставшуюся часть операторов Лр. Обозначим

ЛР =

Л

р1

0 Л

0 )

— диагональные матрицы (4x4)

и матрицы А = А — ( Ар + А 2 + А 3), содержащие оставшиеся части дифференциальных операторов.

Для интегрирования системы (34) воспользуемся методом дробных шагов, который состоит из четырех шагов: шаг 1:

р

п+2/6 (U п+3/6

п+2/6)_ М Ар • Uп+3/6, (37)

где Uп+2/6 — значения, полученные на этапе решения системы Эйлера; шаг 2:

р

шаг 3: р

п+3/6 (u п+4/6 — u п+3/6 )_ М а • U п+4/6.

и+4/6(ии+5/6— и"^1/6)_ М А U •

(38)

(39)

шаг 4:

п+2/6 (Пп+1 — П п+5/6)_ М ((Ар + А 2 + А з)х

р

х и

п+5/6

+ А- и

п+2/6

2

+ Fn+2/6 ^

(40)

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

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

0п +1 (лп Г\п+1 ллп+р I г\ п+1

к — 0к _ро 0к+1 — 20к + 0к-1

М к2

(41)

р

т

00П = 0«, 0

Г\П

0М-V-

где количество расчетных узлов по толщине оболочки М +1 Полученная из (41) система алгебраических уравнений решается методом прогонки. Формула (22) для нахождения безразмерного коэффициента теплообмена g(Fo) в конечно-разностной аппроксимации имеет следующий вид:

0Г1 - 0«

g(Fo)п+l = 0Г «

(42)

(0« -т

Расчет параметров теплообмена. После решения сопряженной задачи (1)—(18) вычисляется тепловой поток к поверхности конструкции де, коэффициент теплообмена ат и число Нуссельта Ми:

д, = а(0« -00); ат = д8/(0е -0«); Ми = ат Ь / X,

где Ь — характерный размер обтекаемого тела.

Результаты численного решения. На рисунках 1—7 представлены результаты численного моделирования обтекания фрагмента корпуса модельного летательного аппарата, имеющего эллипсоидальный профиль (поверхность которого образована двумя половинами эллипсоидов вращения с разными значениями полуосей, см. рис. 1), гиперзвуковым потоком газа. В настоящее время такая форма считается одной из наиболее аэродинамически целесообразных для гиперзвуковых полетов [5, 16]. Начальные условия для задачи и скорость набегающего потока (граничное условие на входной границе) задавались следующими: р = 0,195 кг/м3, ух = уу = 0 м/с, уг =1 800 м/с, р = 12 346 Па, М = 6, #=15 км. Значения размерных параметров задачи принимались следующими: 00 = 293 К, б1 = 0,8, е2 = 0,3, X5 = 0,3 Вт/(м-К), р, =1 800 кг/м3, с, = 0,8 кДж/(кг-К).

Для выявления особенностей течения в наиболее ответственных зонах (носовой части и на кромках ГЛА) необходимо использовать достаточно мелкие сетки, что обусловливает существенное увеличение времени счета и объема машинной памяти, как для проведения непосредственно вычислений, так и визуализации результатов расчетов. Для расчетов использовался вычислительный комплекс Научно-образовательного центра «Суперкомпьютерное ин-

Рис. 1. Носовая часть модельного варианта ГЛА с эллипсоидальным профилем. Изображены линии, вдоль которых проводится анализ расчетов:

1 — верхняя образующая; 2 — нижняя образующая;

3 — кромка

(Полноцветную версию см. http://www. 1/уи/та$Ь.Ьт$Ш.ги)

женерное моделирование и разработка программных комплексов» (НОЦ «СИМПЛЕКС») МГТУ им. Н.Э. Баумана, состоящий из 80 процессоров с общей памятью. Для построения адаптивных высококачественных конечно-разностных почти регулярных сеток использовался специализированный генератор сеток 3БЯМО, разработанный в НОЦ «СИМПЛЕКС». На рисунке 1 показаны три характерные линии на поверхности теплозащитной конструкции ГЛА, вдоль которых анализировались значения решения параметров потока и теплообмена на поверхности ГЛА: вдоль верхней образующей корпуса, вдоль нижней образующей и по кромке ГЛА. Некоторые из результатов расчетов параметров аэрогазодинамического потока, обтекающего поверхность модельного ГЛА в областях У1+У2, представлены на рис. 2, 3. На этих рисунках изображены значения параметров потока: скорости, давления, плотности, температуры в 3-мерном сечении расчетной области потока У1+ У вертикальной плоскостью.

Численное решение с хорошим разрешением воспроизводит головной скачок уплотнения, который образуется в критической точке ГЛА и формирует ударную волну, окружающую поверхность корпуса ГЛА. Метод ТУБ позволяет получить достаточно четкий фронт ударной волны практически на всем протяжении корпуса ГЛА, без размазывания ее на большом удалении от критической точки, что наблюдается обычно при использовании немонотонных разностных схем. Максимумы плотности, давления и температуры находятся около критической точки носка аппарата, максимум температуры

0,0376 368 737 1,11е + 03 1,47е+03

184 553 921 1,29е+03 1,66еН-03

а

7,Обе+03 1,Зе + 05 2,53е+05 3,76е + 05 4,98е + 05 6,85е+04 1,91е+05 3,14е+05 4,37е+05 5,6е+05

б

0,0305 0,242 0,453 0,664 0,875

0,136 0,347 0,559 0,77 0,981

в

Рис. 2. Распределение продольной компоненты скорости [м/с] (а), давления [Па] (б) и плотности

[кг/м3] (в) набегающего газового потока в окрестности поверхности ГЛА в области Ур+У2 (Полноцветную версию см. izvuzmash.bmstu.ru)

170 544 918 1,29е+03 1,б7е + 03

357 731 1,11е+03 1,48е+03 1,85е+03

а

170 544 918 1,29е+03 1,67е+03

357 731 1,11е+03 1,48е+03 1,85е+03

б

Рис. 3. Распределение температуры [К] набегающего газового потока в окрестности поверхности ГЛА в области Ур+ У2:

а — для адиабатической стенки; б — для стенки с учетом теплообмена

(Полноцветную версию см. ЫХр://^"^. izvuzmash.bmstu.ru)

при выбранных условиях численного эксперимента достигает значения 2 040 К.

На рисунке 3 представлены сравнительные картины распределения температуры набегающего потока для двух различных вариантов расчета: для модели адиабатической стенки, когда внешняя поверхность корпуса ГЛА считалась теплоизолированной, и с учетом теплообмена между газовой средой и оболочкой в момент времени 50 с после начала расчета, т.е. на основе полного решения сопряженной задачи. Во втором случае температура в среднем по поверхности тела на 25% ниже, чем в первом, что

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

Распределение числа Рейнольдса по безразмерной продольной координате X, ориентированной по продольной оси ГЛА, приведено на рис. 4, а. Число Рейнольдса при выбранных значениях параметра набегающего потока достигает максимального значения на кромке ГЛА на некотором расстоянии (X = 0,06) от критической точки, и составляет 1,7-108, что свидетельствует о возможности реализации турбулентного режима течения потока в окрестности поверхности ГЛА. Максимальные значения Яе по верхней и нижней образующей линии корпуса располагаются вблизи критической точки ГЛА. Вследствие отличия геометрических параметров верхней и нижней части корпуса ГЛА, кривые значений Яе(Х) отличаются для этих частей и их значения существенно меньше, чем на кромке ГЛА.

Распределение числа Нуссельта по верхней и нижней образующей корпуса и вдоль кромки

Рис. 4. Распределение числа Рейнольдса (а) и числа Нуссельта (б) по продольной координате X для верхней, нижней частей поверхности и по соединяющей их кромке ГЛА:

1 — верхняя образующая; 2 — нижняя образующая; 3 — кромка

показано на рис. 4, б. Максимальные значения № достигаются на кромке для всех значений X, кроме того, максимальное абсолютное значение № достигается на некотором удалении от критической точки, примерно на том же расстоянии (Х = 0,06), что и число Яе, что обусловлено относительно малым радиусом затупления рассмотренной поверхности ГЛА. Соответствующие кривые распределения коэффициента теплообмена представлены на рис. 5.

0,24 0,35 0,46 0,58 0,69 Продольная координата, Т

Рис. 5. Распределения коэффициента теплообмена по продольной координате X для верхней, нижней частей поверхности и по соединяющей их кромке ГЛА:

1 — верхняя образующая; 2 — нижняя образующая; 3 — кромка

Зависимость коэффициента теплообмена от угловой координаты в разных сечениях по X приведена на рис. 6. В верхней части ГЛАтеплооб-мен более интенсивный, чем в нижней его части, что обусловлено различием геометрической формы верхней и нижней частей поверхности ГЛА.

Распределение температуры по телу, рассчитанное на основе решений полной сопряжен-

90 105 120 135 150 165 180 195 210 225 240 255 Угловая координата, град

Рис. 6. Распределение коэффициента теплообмена по угловой координате ГЛА для разных значений продольной координаты X:

1 — X =0,118; 2 — X = 0,240; 3 — X = 0,360; 4 — X = 0,600; 5 — X = 1,000

МАШИНШТРШНИ

ной задачи с учетом теплообмена между теплозащитным материалом корпуса ГЛА и пограничным слоем набегающего потока ГЛА, изображено на рис. 7. В критической точке значение температуры достигает 1 600 К, по мере удаления от критической точки температура монотонно убывает, однако остается достаточно высокой — ее значение на максимальном удалении составляет примерно 800 К для кромок и 1 000 К для верхней части ГЛА, обладающей большим значением угла конусности.

Рис. 7. Распределение температуры

по продольной координате Zдля верхней, нижней

частей поверхности и по соединяющей их кромке ГЛА:

1 — верхняя образующая;

2 — нижняя образующая; 3 — кромка

Выводы

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

2. Проведено численное моделирование обтекания фрагмента носовой части модельного ГЛА эллипсоидальной формы.

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

Литература

[1] Anderson J.D. Hypersonic and High-Temperature Gas Dynamics. American Institute of Aeronautics and Astronautics, Virginia, Reston, 2006. 232 p.

[2] Тирский Г.А., ред. Гиперзвуковая аэродинамика и тепломассообмен спускаемых космических аппаратов и планетных зондов. Москва, ФИЗМАТЛИТ, 2011. 548 с.

[3] McNamara J.J., Friedmann P.P. Aeroelastic and Aerothermoelastic Analysis of Hypersonic Vehicles: Current Status and Future Trends. 48th AIAA/ASME/ASCE/ÄHS/ASC Structures, Structural Dynamics, and Materials Conference, 23—26 April 2007,

Honolulu, Hawaii, available at: http:// www. mecheng. osu.edu/lab/cael/sites/default/files/AIAA—2007—2013 (accessed 20 November 2013 ). 55 p.

[4] Crowell A.R., McNamara J.J., Miller B.A. Hypersonic Aerothermoelastic Response Prediction of Skin Panels Using Computational Fluid Dynamic Surrogates. ASDJournal, 2011, vol. 2, no. 2, pp. 3—30.

[5] Dimitrienko Yu.I. Termomechanics of Composites under High Temperatures. Dordrecht, Boston, London, Kluwer Academic Publishers, 1999. 347 p.

[6] Dimitrienko Yu. Heat-Mass-Transport and Thermal Stresses in Porous Charring Materials. Transport in Porous Media, 1997, no. 27(2), pp. 143-170.

[7] Dimitrienko Yu.I. Thermomechanical behaviour of composite materials and structures under high temperatures: 1. Materials. Composites Part A: Applied Science and Manufacturing, 1997, no. 28(5), pp. 453-461.

[8] Димитриенко Ю.И., Захаров АА., Коряков М.Н., Сызды-ков Е.К. Численное решение сопряженной задачи аэрогазодинамики и внутреннего теплопереноса в конструкциях гиперзвуковых летательных аппаратов. Вестник МГТУим. Н.Э. Баумана. Сер. Естественные науки. Спец. вып. № 6. Моделирование и исследование физических и технических систем, 2012, с. 84—101.

[9] Димитриенко Ю.И., Минин В.В., Сыздыков Е.К. Численное моделирование процессов тепломассопереноса и кинетики напряжений в термодеструктирующих композитных оболочках. Вычислительные технологии, 2012, т. 17, № 2, с. 43—59.

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

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

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

[13] Гильманов А.Н. Методы адаптивных сеток в задачах газовой динамики. Москва, ФИЗМАТЛИТ, 2000. 248 с.

[14] Димитриенко Ю.И. Тензорное исчисление. Москва, Высшая школа, 2001. 575 с.

[15] Калугин В.Т., ред. Аэродинамика. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010. 687 с.

[16] Братчев А.В., Забарко ДА., Ватолина Е.Г, Коробков АА, Сахаров В.И. Вопросы теплотехнического проектирования перспективных гиперзвуковых летательных аппаратов аэробаллистического типа. Известия института инженерной физики, 2009, т. 2, № 12, с. 42—49.

References

[1] Anderson J.D. Hypersonic and High-Temperature Gas Dynamics. American Institute of Aeronautics and Astronautics, Reston, Virginia, 2006. 232 p.

[2] Giperzvukovaia aerodinamika i teplomassoobmen spuskaemykh kosmicheskikh apparatov i planetnykh zondov [Hypersonic aerodynamics and heat transfer reentry spacecraft and planetary probes]. Ed. Tirskii G.A. Moscow, FIZMATLIT publ., 2011. 548 p.

[3] McNamara J. J., Friedmann P.P. Aeroelastic and Aerothermoelastic Analysis of Hypersonic Vehicles: Current Status and Future Trends. 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 23-26 April 2007. Honolulu, Hawaii, available at: http:// www mecheng.

osu.edu/lab/cael/sites/default/files/AIAA—2007—2013 (accessed 20 November 2013). 55 p.

[4] Crowell A.R., McNamara J.J., Miller B.A. Hypersonic Aerothermoelastic Response Prediction of Skin Panels Using Computational Fluid Dynamic Surrogates. ASDJournal, 2011, vol. 2, no. 2, pp. 3—30.

[5] Dimitrienko Yu.I. Termomechanics of Composites under High Temperatures, Dordrecht, Boston, London, Kluwer Academic Publishers, 1999. 347 p.

[6] Dimitrienko Yu. Heat-Mass-Transport and Thermal Stresses in Porous Charring Materials. Transport in Porous Media, 1997, no. 27(2), pp. 143-170.

[7] Dimitrienko Yu.I. Thermomechanical behaviour of composite materials and structures under high temperatures: 1. Materials. Composites Part A, Applied Science and Manufacturing, 1997, no. 28(5), pp. 453-461.

[8] Dimitrienko Iu.I., Zakharov A.A., Koriakov M.N., Syzdykov E.K. Chislennoe reshenie sopriazhennoi zadachi aerogazodinamiki i vnutrennego teploperenosa v konctruktsiiakh giperzvukovykh letatel'nykh apparatov [Numerical solution of the dual problem aerogasdynamics and internal heat transfer in hypersonic aircraft designs]. Vestnik MGTUim. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman Moscow State Technical University. Ser. Natural Sciences]. Spetsial'nyi vypusk no. 6 Modelirovanie i issledovanie fizicheskikh i tekhnicheskikh system [Special issue no. 6 Simulation and study of the physical and technical systems]. 2012, pp. 84-101.

[9] Dimitrienko Iu.I., Minin V.V., Syzdykov E.K. Chislennoe modelirovanie protsessov teplomassoperenosa i kinetiki napriazhenii v termodestruktiruiushchikh kompozitnykh obolochkakh [Numerical modeling of heat-mass-transfer and stress kinetics in thermodestructing composite shells]. Vychislitel'nye tekhnologii [Computational Technologies]. 2012, vol. 17, no. 2, pp. 43-59.

[10] Dimitrienko Iu.I., Zakharov A.A., Koriakov M.N. Model' trekhmernogo pogranichnogo sloia i ee chislennyi analiz [Three-dimensional model of boundary layer and its numerical analysis]. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman Moscow State Technical University. Natural Sciences]. 2011, special issue, pp. 136—150.

[11] Dimitrienko Iu.I., Kotenev V.P., Zakharov A.A. Metod lentochnykh adaptivnykh setok dlia chislennogo modelirovaniia vgazovoi dinamike [Adaptive mesh tape method for numerical simulation of gas dynamics]. Moscow, FIZMATLIT publ., 2011. 286 p.

[12] Dimitrienko Iu.I., Koriakov M.N., Zakharov A.A., Syzdykov E.K. Razvitie metoda lentochno-adaptivnykh setok na osnove skhem TVD dlia resheniia zadach gazovoi dinamiki [Development of the band-adaptive grids based on TVD schemes for solving problems of gas dynamics]. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman Moscow State Technical University Natural Sciences]. 2011, no. 2, pp. 87—97.

[13] Gil'manov A.N. Metody adaptivnykh setok v zadachakh gazovoi dinamiki [Methods of adaptive grids in problems of gas dynamics]. Moscow, FIZMATLIT publ., 2000. 248 p.

[14] Dimitrienko Iu.I. Tenzornoe ischislenie [Tensor calculus]. Moscow, Vysshaia shkola publ., 2001. 575 p.

[15] Golubev A.V., KaluginV.T, Lutsenko A.Iu., GolubevA.G., Kalugna V.T. Aerodinamika [Aerodynamics]. Moscow, Baumana Press, 2010. 687 p.

[16] Bratchev AV, Zabarko D.A., Vatolina E.G., Korobkov AA, Sakharov V.I. Voprosy teplotekhnicheskogo proektirovaniia perspektivnykh giperzvukovykh letatel'nykh apparatov aeroballisticheskogo tipa [Questions teplotehnicheskogo design perspective aeroballistic hypersonic aircraft type]. Izvestiia institu-ta inzhenernoi fiziki [Proceedings of the Institute of Engineering Physics]. 2009, vol. 2, no. 12, pp. 42-49.

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

Информация об авторах

ДИМИТРИЕНКО Юрий Иванович (Москва) — директор Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов», доктор физико-математических наук, профессор, зав. кафедрой «Вычислительная математика и математическая физика». МГТУ им. Н.Э. Баумана (105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1, e-mail: dimit.bmstu@gmail.com).

ЗАХАРОВ Андрей Алексеевич (Москва) — старший научный сотрудник Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов», доцент кафедры «Вычислительная математика и математическая физика». МГТУ им. Н.Э. Баумана (105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1, e-mail: azaharov@bmstu.ru).

КОРЯКОВ Михаил Николаевич (Москва) — младший научный сотрудник Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов», ассистент кафедры «Вычислительная математика и математическая физика». МГТУ им. Н.Э. Баумана (105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1).

СЫЗДЫКОВ Елтуган Кимашевич (Московская область) — кандидат технических наук, первый заместитель Генерального директора ОАО «ГосМКБ «Радуга» им. А.Я. Березняка» (141980, Московская область, Дубна, Российская Федерация, ул. Жуковского, д. 2а, e-mail: raduga@dubna.ru).

Information about the authors

DIMITRIENKO Yuriy Ivanovich (Moscow) — Director of Scientific Educational Center «Supercomputer Engineering Modeling and Development of Software Systems», Dr. Sc. (Phys. Math.), Professor, Head of «Computational Mathematics and Mathematical Physics» Department. Bauman Moscow State Technical University (BMSTU, building 1, 2-nd Baumanskaya str., 5, 105005, Moscow, Russian Federation, e-mail: dimit.bmstu@gmail.com).

ZAKHAROV Andrey Alekseevich (Moscow) — Senior Researcher of Scientific Educational Center «Supercomputer Engineering Modeling and Development of Software Systems», Associate Professor of «Computational Mathematics and Mathematical Physics» Department. Bauman Moscow State Technical University (BMSTU, building 1, 2-nd Baumanskaya str., 5, 105005, Moscow, Russian Federation, e-mail: azaharov@bmstu.ru).

KORYAKOV Mikhail Nikolaevich (Москва) — Junior Researcher of Scientific Educational Center «Supercomputer Engineering Modeling and Development of Software Systems», Assistant of «Computational Mathematics and Mathematical Physics» Department. Bauman Moscow State Technical University (BMSTU, building 1, 2-nd Baumanskaya str., 5, 105005, Moscow, Russian Federation).

SYZDYKOV Eltugan Kimashevich (Московская область) — Cand. Sc. (Eng.), First Deputy Director General of JSC «Raduga» n. a. A.Ya. Bereznyak (Zhukovskogo str., 2a, 141980, Dubna, Moscow region, Russian Federation, e-mail: raduga@dubna.ru).