МЕТОДИКА, АЛГОРИТМ И ПРОГРАММНЫЙ КОМПЛЕКС ПО ОПРЕДЕЛЕНИЮ ТЕПЛОВОГО СОСТОЯНИЯ ОХЛАЖДАЕМЫХ МИКРОРАКЕТНЫХ ДВИГАТЕЛЕЙ Текст научной статьи по специальности «Физика»

Электронный журнал «Труды МАИ». Выпуск № 78

www.mai.ru/science/trudy/

УДК 536.21

Методика, алгоритм и программный комплекс по определению теплового состояния охлаждаемых микроракетных двигателей

Формалёв В.Ф.,* Колесник С.А.**

Московский авиационный институт (национальный исследовательский университет), МАИ, Волоколамское шоссе, 4, Москва, A-80, ГСП-3, 125993, Россия

*e-mail: formalev38@yandex. ru **e-mail: sergey@oviont.com

Аннотация

В работе приведена комплексная физико-математическая модель, алгоритм и описание программного комплекса по определению теплового состояния блока охлаждаемых ракетных микродвигателей, используемых для управления космическими аппаратами. Изложен метод численно-аналитического решения сопряженных задач теплообмена в гидрогазодинамических каналах и трехмерной нестационарной теплопроводности в многосвязных областях. Предложен метод погружения с использованием метода дробных шагов численного решения сопряженных задач теплообмена и теплопроводности в многосвязных областях с негладкими границами. Разработан и описан программный комплекс по модульному принципу на языке ФОРТРАН. Получены и проанализированы результаты по тепловым потокам в охлаждаемых каналах и газодинамическом сопле, а также трехмерные нестационарные температурные поля.

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

Введение

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

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

- об идеальном течении газа в сопле микродвигателя;

- о течении и теплообмене в тепло-газодинамическом пограничном слое;

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

- о нестационарном трехмерном теплопереносе в блоке газодинамических микродвигателей.

В силу малых размеров блока микродвигателей (длина камеры сгорания с соплом - не более 15 мм) первые два пункта, по-видимому, можно объединить уравнениями тепло-газодинамики на основе уравнений Навье-Стокса.

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

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

В такой комплексной постановке задача решается впервые, поскольку по отдельным компонентам имеется значительное число работ, среди которых следует отметить работы Шлихтинга Г. [1], Формалева В.Ф. [2, 3], Михеева М.А. [4], Галицейского Б.М., Совершенного В.Д., Формалева В.Ф. [5], Формалева В.Ф., Кузнецовой Е.Л. [6], а также Авдуевского В.С. и др. [7].

Поэтому данная работа по решению комплексной проблемы сопряженного теплообмена между гидро- и газодинамическими течениями и телами сложной

геометрической формы является весьма актуальной и востребованной проектировщиками блоков микроракетных двигателей.

1. Физико-математическая постановка

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

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

потоков

а.

, (Т, - Т, )-£ аТ 4 +£ аТ4 -Л

1 V е1 м1 / w1 w1 е е

дТ_ дп

= 0, (х, y, г) е 5^1,

(1)

м1

где а ,, Т, - коэффициент теплоотдачи от газодинамического потока к стенке и

эффективная температура пограничного слоя; , £е

степени черноты

поверхности м1

и газа, а

- постоянная Стефана-Больцмана

а = 5,67 • 10"

Ад

I 2 • К4

Л - теплопроводность материала корпуса, Т -

температура, п - направление внешней нормали (от границы в сторону

газодинамического потока.

Рис. 1. Отдельный ракетный микродвигатель в блоке с охлаждающими

каналами

2. Уравнение теплопереноса в трехмерной многосвязной области с двумя цилиндрическими каналами для охладителя

\Т)р— = — д? дх

л—

V дх у

+

ду

V дУ у

+

А д2

лдТ

V д 2 у

, (х, у, 2)е V

(2)

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

направлении оси 2

дТ

а 9 (Т9 - Т 9) + Л —

w2 V е2 w2 / ~

дп

0, (x,У,2)е ^2.

(3)

w 2

4. Аналогичное соотношение на границе второго канала для охладителя с

осью, проходящей через точку О2 в направлении оси 2

дТ

^3 (Те3 - ^3 ) + Л —

дп

0, ( х 2 )е Swз.

(4)

В соотношениях (3) и (4) п - внешняя нормаль в направлении от границы тела к охладителю. Нормальные производные температуры определяются

соотношением

дТ дТ , ч дТ , ч дТ , ч

— = — cos (п, х) л--cos (п, у) л--cos (п, г).

дп дх У ' ду К ' дг К '

(5)

5. Параметры теплообмена а№1, Те1 от газа к границе м! определяются соотношениями [7].

0,8 0,8 0,2

а^= 0,0296 Pr-0,57 • Кт р 'Щ К

Т 0,2

х,\

уо

(6)

1 + а к -1 , 2

е1 0

1 + а

, а =- 1 2, г = 0,88;

2

(7)

где

Кт =

х ч 0,4+0,2ехр(-га)

Т

V Те1 У

/1 \°,11 (1 + аг ) ;

\р^1и1К 54 ^

х/— хп л

уо уо 0

р м.Я

5/4

в плоской камере сгорания равна ее толщине; х - текущее расстояние; щ -скорость на оси газодинамического канала определяется ниже; М1 = щ / а; а -скорость звука; Рг - число Прандтля (для газа его можно принять 0,7-0,85);

к = ср / су - постоянная адиабаты; ер, еу - теплоемкости при постоянном давлении и

р ' v

постоянном объеме; Т0 - температура торможения в камере сгорания.

6

х

6. Определение параметров теплообмена ак2, Те2 на границах м2 и м3. Стационарное ламинарное течение газообразного охладителя (число Прандтля Рг = 1) описывается следующей краевой задачей (без инерционных сил):

дг

г ди г ц—

V дг.

= г

др ди дг' дг

г=0

0 и|г= 0

(8)

решением которой является функция

и

( г ) = -

др_ 8ц дг:

(9)

в которой падение давления др / дг определяется формулой

-др = %•£!*— о = 2 я дг 2В

1'

(10)

а коэффициент сопротивления % определяется выражением

Яе V

и =

ср

Я[_ др 8ц дг

(11)

Соотношения для теплового пограничного слоя имеют вид

дТ

дТ Л д

сри--V сру— =--

дг дг г дг

V дг У

+ Л

д 2Т

дТ дг

= 0, Т| = Т ,.

' \г=К1 м 2

(12)

Используя продольный компонент скорости и из выражения (9) и введя новые переменные

а = ^^ % = ч = ±

Т9 - Т.. Я Я

е2 ао 1 1

получим

+ 1д^ = Ре(1 (14)

д%2 % д% Х } д% дг!2

и ср • О

где Ре = —--число Пекле, а = Л / ср. Для числа Ре > 10 решение задачи (12)

а

представляется в виде ряда

3 = 0,82ехр (-14,6 ^ ) + 0,098ехр (-88,2А ) + 0,019ехр (-212 А ) +...

А = . (15)

и О

ср

Выражение (13) определяет среднюю температуру в сечении г.

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

аппроксимируются зависимостями

Г , 2Ре • Я ^ 3,66 аёу -1 < 12,

Ыи = \ Р (16) , ч 2 Ре • Я 1,61(2Ре • Я / Р) аёу -1 > 12,

Р

Лт тт а 7 • О

где Ыи - число Нуссельта, Ыи = —^-.

Л

В соотношениях (8)-(16) введены следующие обозначения: г - радиальная переменная, отсчитываемая от оси канала; ц,у - динамическая и кинематическая вязкости; г - продольная переменная г е(0; Р); О = 2 • Рх; и, V - продольный и поперечный компоненты вектора скорости; Яе - число Рейнольдса; плотность р, динамическую вязкость ц и теплопроводность Л можно определить по формуле Сезерленда

р Т А 1 + 216/ Т

г „ л1/2

Т

Т V Т0 У

± = (17)

А> А

Р0 Т А 1 + 216 / Т где р0, а , Д0 - плотность, динамическая вязкость, теплопроводность при некоторой фиксированной температуре, например, при Т0 = 300Е.

Таким образом, соотношения (16) определяют коэффициент теплоотдачи аш на стенках охлаждающего канала с радиусом поперечного сечения Кх, а соотношения (13), (15) - температуру Те2 охладителя на границе »2 этого канала. Аналогичные соотношения можно выписать для охлаждающего канала с радиусом поперечного сечения Я2.

Для определения давления р (г) вдоль канала в зависимости от давления во входном сечении канала р (г = 0 ) = р0 можно воспользоваться соотношениями (10), (11)

64

р(иср )2

р (г ) = р, -— г, (18)

у ; 0 Яе 2й

где (ис^) - скорость на входе в канал.

2. Метод численного решения сопряженного теплообмена и

теплопроводности

Поставленная задача является совместной задачей теплообмена между

вязкими газодинамическими течениями и теплопроводности в трехмерных

9

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

Для решения таких сложных задач сопряженного теплообмена и теплопроводности в нестационарном случае предлагается метод погружения многосвязной области с негладкой границей в область простейшей геометрической формы, в данном случае в параллелепипед с размерами x = AL, y = O2F z = L (рис. 1), и дальнейшим использованием итерационного процесса на каждом временном слое. В этом случае параметры теплообмена awl, Те1, aw2, Те2, aw3, ТеЪ на первой

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

На область классической формы, полученную в результате погружения, накладывается конечно-разностная сетка с равномерным шагом h3 в направлении

оси Oz и неравномерными шагами h1;, i = 1,4 и h2 j = 1,4, учитывающими

различные подобласти, на которые разбивается вся расчетная область.

Для решения всей существенно нелинейной задачи используется конечно-

разностный метод дробных шагов Яненко Н.Н. с модификацией сквозного счета

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

Конечно-разностная схема метода дробных шагов с модификацией для сквозного счета с неравномерными шагами имеет вид

-¿+1/3

' Р

А А

т

к + к

1г 1г

п ¿+1/3 . п ¿+1/3

4-1А + ЛА

2К

в

¿+1/3 1]1

Т

к+1/3

А ,

п ¿+1/3 , п ¿+1/3

4'/ + 4+1 л

2ки

в

¿+1/3 _ грк+иъ

А Т'+1А ,

, = 1,1, а = 1,3, / = 1, IV;

(19)

грк+2/3 грк+1/3 -¿+2/3 1ч1 — у1

у I ' Р

Т

2

к ■ + к

2 А 2 А

пк+2/3 . пк+2/3

А-- и + А-;

2к

в

¿+2/3 _ /^¿+2/3

4-1/ А

)-

2А

^¿+2/3 2

+ ^+1/ / грк+213 _ грк+213 ' '(1,}1 Т4+1/

2к

2 А

,, = 1,1, 4 = 1,3, / = 1, V;

(20)

грк+1 грк+2/3 -¿+1 1ч1 ~ 1г)1 1

'р-= г

т к3

пк+1 л ¿+1

4,7-1 + 4

Чу7

2 к

(грк+\ грк+\ \

Т4/-1 - Т4/ )-

п ¿+1 л ¿+1

47 +4/+1

2 к

(грк+\ грк+\ \

ТА - Т

А+1 ,/

, = 1,1, 4 = 1,3, / = 1, V;

(21)

где теплофизические характеристики /к+п'/3 {ск+"'13, Лк+п/3|, и = 1,2,3 на верхних

временных дробных слоях вычисляются с помощью экстраполяционных формул по температурам на двух предыдущих временных слоях

к+п/3 _ 0 гк+(п-1)!Ъ гк+(п-2)/3

ук+п/з _

г

(22)

Пространственные шаги \г, к2г берутся слева от узла I, ], I, а шаги , Л2г -

справа от этого узла.

Аппроксимация краевых условий (1), (3), (4) осуществляется неявным образом, включая нелинейные лучистые тепловые потоки. При этом в конечно-разностной аппроксимации краевых условий возникают, так называемые, консервативные слагаемые, учитывающие аккумуляцию тепловой энергии в граничных узлах [6]. Для граничного условия (1) такая аппроксимация имеет вид

а

(грк+ИЪ грк+ИЪ \ 1 (/е! )~Л

,к+2/3 /г^к+2/3 w1

¿+2/3

АТ

¿+2/3

Ап

(гг^к+2/3 \4 /грк+2/3 \4

(^1 ) + ^(Т,1 ) =

w1

Су1Р ' Ап^1 ( грк+2/3 _ грк+1/3 \

Tw1 I,

т 4 7

(23)

где Апм - половина шага вдоль нормали п, примыкающего к границе w1.

Аппроксимация производной температуры в нормальном направлении следует из соотношения (5)

АТ

к+2/3

Ап

w1

грк+2/3 грк+2/3 Т+1 ]1 ~ Т,-1 ]1

2 (Л )„

грк+2/3 грк+2/3

I \ , ТЧ+11 — Тч1 сов ( п, х ) + ——;—^--сов (п

( Л2 ),,

Т — Т

(п, у)+ ^— 1]1—1 сов (п, 7). (24)

Алгоритм численного решения содержит следующие пункты:

-определяются параметры теплообмена аы, Те1, ак2, Те2, а№3, Те3, используя распределение температур на к -м временном слое;

- на первом дробном шаге по времени (слой к +1/3) выполняются скалярные прогонки в направлении оси Ох с использованием подсхемы (19) во всей расчетной области; в результате получается температурное поле Т^+1/3;

- на втором дробном шаге по времени (слой к + 2/3) выполняются скалярные прогонки в направлении оси Оу с использованием подсхемы (20) и распределения

Т^+1/3; при этом используется метод, изложенный в [6] с неявной аппроксимацией

лучистого потока и решением в квадратурах на границе w1 нелинейного уравнения

~ гтк+2/3

четвертой степени; в результате получается температурное поле ;

- на третьем дробном шаге по времени (слой к +1) выполняются скалярные прогонки в направлении оси Ог по подсхеме (21) во всей расчетной области, в результате чего получается решение Т,к+1 на (к +1) -м временном слое;

- итерационный процесс повторяется, начиная с первого пункта уточнением параметров теплообмена с новым распределением температур; этот процесс повторяется не более двух раз.

3. Программный комплекс

По изложенному методу и алгоритму по модульному принципу разработан программный комплекс на языке ФОРТРАН по расчету тепловых потоков от высокотемпературного газодинамического пограничного слоя, от пограничного слоя

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

На рис. 2 приведена иерархическая структура программного комплекса, где в вершинах отмечены имена программных модулей, а стрелками обозначены направления вызова программных модулей.

Наиболее трудоемкой программой является управляющая программа ЦР,

находящаяся в вершине иерархической структуры.

Она описывает переменные и массивы, формирует многочисленные входные

данные, обращается и выполняет программные модули второго уровня УРЯ, ХРЯ,

14

ZPR, AW1, AW23, MNK, формирует многомерные массивы выходных данных и выводит их на печать.

Программные модули второго уровня XPR, YPR, ZPR вычисляют коэффициенты матриц систем алгебраических уравнений (СЛАУ) и их правые части, по ним определяются прогоночные коэффициенты и распределения температур в скалярных прогонках в направлении осей Ox, Oy, Oz соответственно.

Программный модуль AW23 вычисляет параметры теплообмена а2, Те2, а3, Те3 на стенках двух охлаждаемых каналов в зависимости от переменной z. Модуль AW1 вычисляет параметры теплообмена на границе w1 с газодинамическим пограничным слоем в зависимости от переменных x и z и температуры Twl (x, z).

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

Модули третьего уровня обеспечивают модули второго уровня XPR, YPR, ZPR необходимой информацией.

Модуль YW вычисляет прогоночные коэффициенты и температуры в узлах на границе w1, модуль DXY формирует дробные шаги по пространственным переменным, в зависимости от местонахождения узла, модуль INTERP интерполирует входные табличные функции по заданным значениям аргументов. Модули XW и ATB формируют коэффициенты СЛАУ на границе w1 вдоль оси Ox.

Модули AZET, ALPHA, ITER, AGI, PSI формируют газодинамические характеристики для модуля AW1.

Входные данные в программный комплекс:

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

- температура и давление торможения в камере сгорания;

- давление и скорость охладителя на входе в каналы охлаждения, а также теплофизические характеристики охладителя с, р, ¡и, Л;

- зависимости от температуры теплофизических характеристик материала корпуса Л(Т), с (Т), р.

Выходные данные: параметры теплообмена аы(х,г,г), Те1 (х,г,г), ак2(г,г), Те2 (г,г), а№3 (г,г), Те3 (г,г), распределение температур Т(х,у,г,г).

4. Результаты численного решения

На рисунках 3-5 приведены некоторые результаты о тепловом состоянии микроракетного двигателя, полученные с помощью программного комплекса. Пакет плоских микродвигателей длиной 14 мм с расстоянием в критическом сечении 0,4 мм и расстоянием на среде сопла 5,0 мм функционирует не более 10 секунд.

Рис. 3. Параметры теплообмена на границе и! в момент времени ^ = 2,6 с.(

----;

)

На рисунке 3 приведены параметры теплообмена ^, Тл на границе и!, омываемой высокотемпературным газодинамическим потоком, из которого видно резкое падение теплового потока - Ты) в выхлопном сопле за счет

значительного падения давления и плотности газа.

Рисунок 4 демонстрирует величину нагрева газообразного охладителя при течении его в левом канале в прямом направлении, а затем при течении того же охладителя в обратном направлении в правом канале. За время, равное 5,6 с, охладитель в обоих каналах нагревается до температуры 800 К, то есть охладитель снимает значительное количество тепловой энергии, поступающей в корпус от газодинамического потока.

Те2, К

ТеЪ, К

1 = 2,6с

400

800

600

200

0

О

г П.

I

Рис. 4. Температура охладителя в левом канале при течении его в прямом направлении и в правом канале - при течении в обратном направлении.

На рисунке 5 приведены результаты расчетов распределения температурного распределения в критическом сечении микроракетного двигателя. Максимальная температура на границе мЛ в критическом сечении, достигаемая за 5,6 с. равна 8500С, что еще может выдержать стальной корпус, правда при резком снижении прочности.

Рис. 5. Распределение температуры в критическом сечении корпуса

микродвигателя.

Заключение

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

2. Разработаны методы и алгоритмы определения параметров теплообмена на свободных границах многосвязной области, каковой является корпус пакета микроракетных двигателей.

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

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

Работа выполнена при финансовой поддержке Министерства образования и науки Российской Федерации, проект № 118, руководитель Формалев В.Ф.

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

1. Шлихтинг Г. Теория пограничного слоя. - М.: Наука, 1969. 742 с.

2. Формалев В.Ф. Методика, алгоритм и программный комплекс решения задач о тепловом состоянии теплозащитных композиционных материалов при аэродинамическом нагреве // Электронный журнал «Труды МАИ», 2014, № 72: http://www.mai.ru/science/trudy/published.php?ID=47581 (дата публикации 27.01.2014).

3. Формалев В.Ф., Воробьева О.Ф. Метод переменных направлений с экстраполяцией численного решения задач теплопроводности с тензором теплопроводности и конвективными членами // Вестник Московского авиационного института. 1998. Т. 5. № 1. С. 41-48.

4. Михеев М.А., Михеева И.М. Основы теплопередачи. - М.: Энергия, 1973. -

320 с.

5. Галицейский Б.М., Совершенный В.Д., Формалев В.Ф. Тепловая защита лопаток турбин. - М.: Изд-во МАИ, 1996. - 356 с.

6. Формалев В.Ф., Кузнецова Е.Л. Тепломассоперенос в анизотропных телах при аэрогазодинамическом нагреве. - М.: Изд-во МАИ, 2010. - 308 с.

7. Авдуевский В.С., Галицейский Б.М., Глебов Г.А. и др. Основы теплопередачи в авиационной и ракетно-космической технике. - М.: Машиностроение, 1992. - 552 с.