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

УДК 621.793+ 629.78

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

© 2013 А.И. Рязанов

Самарский государственный аэрокосмический университет имени академика С.П. Королёва (национальный исследовательский университет)

Поступила в редакцию 02.12.2013

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

Ключевые слова: ракетная камера, металлизатор, задача Стефана, численное решение, алгоритм расчета, интерфейс программы.

Ракетный двигатель является эффективным генератором энергии. Горение топлива происходит в ограниченном объеме при давлении многократно превышающем атмосферное. Высокотемпературные продукты сгорания, приближаясь к скорости звука, истекают через небольшое критическое сечение. В этом сечении концентрация энергии, запасенная в рабочем теле, максимальна и уступает лишь показателям, достигаемым в ядерной промышленности. Область сосредоточения энергии может быть использована для плавки металлов, а поток продуктов сгорания для транспортировки расплава. Подобная задача возникает при создании металлизаторов - устройств, наносящих покрытия из материалов в форме проволоки. На кафедре механической обработки материалов СГАУ был спроектирован [1], изготовлен и испытан [2] металлизатор на газообразной топливной паре: пропан, воздух. Устройство предназначено для использования в наземных условиях. Материал наносимого покрытия - алюминий. Разработка базируется на опыте и конструктивных решениях, применяемых при создании ракетных двигателей малой тяги [1]. Проволока подается в область критического сечения, где теплоотдача рабочего тела максимальна (рис. 1). Она может быть представлена как цилиндр конечных размеров, поступательно движущийся в среде, с которой осуществляется конвективный и радиационных теплообмен. В случае постоянной скорости подачи устанавливается равенство между объемом внесенной в камеру проволоки и объемом получае-

Рязанов Александр Ильич, ассистент кафедры механической обработки материалов. E-mail: [email protected]

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

Аналитическое решение задачи Стефана в классической постановке на данный момент никем не получено. Существует вариант численного решения задачи нагрева и плавления подвижной проволоки внутри камеры сгорания металлиза-тора методом конечных объемов [3]. Расчетная область разбивается на сетку конечных объемов. Каждому узлу сетки соответствует свое уравнение баланса удельных энергий:

Е = Е + Е + Е + Е

(1)

На рис. 2 показана схема распространения энергий для одного конечного объема (КО).

Рис. 1. Принципиальная схема процесса плавления проволоки: 1 - металлизатор; 2 - проволока; 3 - поток продуктов сгорания уже с частичками расплава

1Г

X

Рис. 2. Схема распространения энергий для конечного объема:

--связи между узлами КО,

------границы КО.

sv =■

V дТ _dx_

V

SvL V

количество теп-

s, =

1 dqr

f

Ar dr

sx =

1 dq

Л

A dx

- количество

V у

тепла, проходящее расстояние между соседними в радиальном (осевом) направлении узлами сетки за время т , отнесенное к площади перпендикулярной направлению его распространения Аг.

*

s=

1

q.

r

- дополнительное количество

s = s + s д = — + qр

вн конв рад v V

количество

а

s„„„e = ■

' Аа ■ (ТОС ТПР )

V

(2)

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

Л

р ■ с ■ А

Ь прив 0 а

Т

1 г1

Т

ИЛ

100

100

s рад =

V

•,(3)

где ТОТ - температура облекающего проволоку тела, ТПР - температура конечного объема проволоки, С0 - коэффициент лучеиспускания абсолютно черного тела, £прив - приведенный коэффициент излучения для системы двух тел равный:

рприв

1

(4)

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

ПР

ЮТ

-1

v' кс

где - <!;ПР и £КС степени черноты проволоки и камеры сгорания, АОТ - площадь облекающего тела (камеры сгорания).

Тепловые потоки , дг в металле распространяются по закону Фурье:

1 л dT , А дТ

qx = Ax, qr = Ar.

дх дг

(5)

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

Источник тепла, а значит и евн присутствует только в объемах расположенных на поверхности цилиндра, т.е. контактирующих с внешней средой. Тепловой поток приходящий извне к поверхностному объему заменяется равным ему тепловым потоком, вырабатываемым внутренним источником тепла находящемся в этом же объеме [4].

Задача сводиться к нахождению температуры в каждом узле сетки. Для нахождения численного решения заменим частные производные температуры в точке (I .Дх, ] -Аг ) через разностные соотношения по формулам конечных разностей (для г аналогично):

дТ . _11

дх

Т - Т1 .

1.1 1 -1.1 , гг -;

Дх

д 2Т

Т 1 . - т

1 -1.1 1,3

т.,, . - т 1+1,1 1J-+Р;

(6)

(7)

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

дх2 Дх2 Дх2

где р ...р4 - остаточные члены, стремящиеся к нулю при стремлении к нулю Дх и Дг, которые можно отбросить.

Для радиационной составляющей сделаем замену, переменную Tt. возведенную в четвертую степень аппроксимируем кусочно-заданной

(\ 4

Т ) = K1 ■ Т + K2.

Для температуры об4лекающего тела замены не требуется, т.к. (ТОТ ) = const.

K1 = (Т1- ТТ2) , K2 = (Т )4 - K1 ■ Т. (8) Т1 Т 2

При малой разности температур Тх — Т2 на границах участка кусочно-заданной функции погрешность аппроксимации будет также мала. Сделав замены и подставив все приведенные зависимости в уравнение баланса удельных энергий (1) получим линейное уравнение (9) для одного узла сетки конечных объемов.

Распределение температур по узлам можно получить, решив систему из таких уравнений. Используем для этого один из разновидностей метода Гаусса - метод так называемых ЬИ - разложений. Система приводится к виду А ■ X = В, где А - квадратная матрица коэффициентов размера п, X - вектор неизвестных размера п, В - вектор правых частей уравнения размера п. Неизвестными являются температуры в узлах. Вектор В содержит известные граничные условия и формируется из членов уравнений соответствующих пограничным узлам, которые контактируют с продуктами сгорания (рис. 3). Матрица А представляется в виде произведения двух треугольных матриц специального вида А = Ь■и. Матричные преобразования требуют существенных затрат времени. Однократный расчет поля температур подразумевает десятки итераций поиска коэффициентов кусочно-заданной функции и итерации подбора скорости подачи проволоки исходя из фиксированного положения вершины конуса границы плавления (рис. 1).

сР-p-Vn

т - T Т - Т

• ,j _ х • • ,j

^Р Н • ПР . " .2

Ах Дг

Т. .+, - Т. . Х т .+, - Т. . Т. , . - Т. .

Дг г. Дг

Дхz

+х

Т1+1,j - Т1. , аа ' Aa,i,j ' (ТОС Т1. )

Дх

Vi

(9)

Е . с . A

Ьприв ,i,j 0 a (гр \4

—^^—:—г' (1оТ)

100* V..' Ср 'р Е ' С ' A

Ьприв ,i, j 0 a,i, j

1004 'V,j' Ср 'P

(K1T, j. + K 2)

временных средств программирования. В среде разработки Delphi 7 создана авторская программа PLAVCA для проведения теоретических исследований образования расплава в ракетной камере под воздействием потока продуктов сгорания. Интерфейс программы прост и интуитивно понятен. Пользователь имеет возможность задать все условия задачи: геометрические параметры расчетной области, размер и количество конечных объемов, физические параметры потока продуктов сгорания и расплавляемого металла. Введенные однажды начальные условия можно сохранять в файл и загружать при новом запуске программы. Расчет проводится как для одного квазистационарного состояния системы стержень - рабочее тело, так и для серии состояний получаемых при изменении одного из параметров задачи с некоторым шагом. Такими параметрами могут быть: коэффициент избытка окислителя, давление в ракетной камере, массовый расход топлива, диаметр критического сечения, диаметр стержня. Результаты расчета представляются в виде распределения температур на каждой итерации и сохраняются в файле формата Excel. Также сохраняется файл с основными результирующими физическими величинами, полученными при достижении равновесия системы.

Программа PLAVCA располагается в единственном файле plavca.exe, размер которого не превышает 1 Mb. Обязательным условием работы программы является наличие на персональном компьютере операционной системы семейства WINDOWS и компонента Excel программного пакета Office. Данное программное обеспечение, согласно статистическим данным, установлено на 90-95% вычислительных машин. Основная нагрузка во время выполнения расчетов приходится на процессор. Чем выше тактовая частота процессора, тем выше скорость вычислений и тем более подробное разбиение расчетной области можно задать. Для выполнения программы необходимо запустить командный файл plavca.exe.

Моделирование плавления осуществляется путем удаления контрольных объемов с поверхности раздела фаз. Критерием для удаления является температура конечного объема, превышающая температуру плавления материала. Сетка перестраивается, снова проводятся матричные преобразования, и расчет повторяется заново. Таким образом, итеративно подбирается линия раздела фаз.

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

Рис. 3. Сетка, включающая все девять отличных конечных объемов, включая пограничные

Т,К -1060 1040 1020 -1000 980 -960 -940 920 -900

С

Рис.

] Б00 1000 1Б00 2000 2Б00 3000 3500 кол-но КО

4. Зависимость температур ТКОН (внизу) и ТКОН+РАД (вверху) от количества конечных объемов

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

Задача называется поставленной корректно, если для любых значений исходных данных из некоторого класса ее решение существует, единственно и устойчиво. Задача получения расплава в ракетной камере поставлена корректно. Рассматривается квазистационарная стадия плавления проволоки, когда искомое распределение температур в материале постоянно во времени, т.е. дТ/ дт ^ 0. Заданы механизмы протекания тепловых потоков внутри материала и известна мощность источника энергии (уравнения 2, 3). Система линейных алгебраических уравнений, например, составляется всегда таким образом, чтобы количество уравнений совпадало с количеством неизвестных. Значит, решение всегда существует и оно единственно.

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

вания в рассматриваемой модели достигается увеличением количества конечных объемов в расчетной области. Точность определения скорости подачи проволоки повышается уменьшением шага итерации. Детальность сетки конечных объемов и минимальные значения шагов итерации ограничиваются только производительностью ЭВМ. Необходимо находить компромисс между точностью расчета и затратами времени на него. Проведем серию расчетов и посмотрим, как изменяется температура в наиболее тепло-напряженном узле сетки при увеличении количества конечных объемов (рис. 4).

Температура узла стабилизируется и перестает существенно зависеть от подробности сетки. Допустимая точность может быть получена уже при 2500-3000 конечных объемах в расчетной области. Стремление температуры к предельному значению подтверждает сходимость численного метода.

Устойчивость численного решения подтверждается монотонностью возрастания графиков на рис. 4. Решение непрерывно зависит от исходных данных и изменяется пропорционально их изменению. К неустойчивости решения могло бы привести накопление погрешностей округлений в процессе вычислений. Однако при программной реализации модели получения расплава в ракетной камере, для дробных чисел использовался действительный тип данных Extended. Число значащих разрядов для данного типа достигает 20-ти. Это позволило исключить влияние погрешностей округления.

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

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

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

1. Рязанов, А.И., Первышин А.Н. Стенд для исследований внутрикамерных рабочих процессов ракетных

двигателей малой тяги // Вестник Самарского государственного аэрокосмического университета им. С.П. Королева. 2011. №3 (27). Ч. 4. С 97-102.

2. Рязанов, А.И., Первышин А.Н. Рабочий процесс метал-лизатора на газообразном топливе / / Вестник Самарского государственного аэрокосмического университета им. С.П. Королева. 2011. №3 (27). Ч. 4. С. 103-109.

3. РязановА.И. Численное решение задачи Стефана для подвижного цилиндра в струе продуктов сгорания с учетом радиационного теплообмена. // Четвертая Российская Национальная Конференция по Теплообмену. М..: 2006. Т. 7. С. 124-127.

4. Пантакар С.В. Численное решение задач теплопроводности и конвективного теплообмена при течении в каналах [пер. с англ. Е.В. Калабина; под ред Г.Г. Янькова]. М.: Издательство МЭИ, 2003. 312 с.

PROGRAM IMPLEMENTATION OF THE NUMERICAL SOLUTION OF THE STEFAN PROBLEM FOR ROD INTRODUCED INTO THE ROCKET CHAMBER'S COMBUSTION FLOW

© 2013 A.I. Ryazanov

Samara State Aerospace University named after Academician S.P. Korolyov (national research university)

In this work, a theory grounded numerical solution of the Stefan problem for rod introduced into the flow of combustion products.Describes the approaches used in the software implementation of this method. With help of the program created by the author proved the correctness of the problem and the stability of solutions in a wide range of initial conditions. Is confirmed the accuracy and convergence of the numerical method. Keywords: rocket chamber, metallizer, Stefan problem, numerical solution, calculation algorithm, program interface.

Alexandr Ryazanov, Assistant Lecturer at the Machining Materials Department. E-mail: [email protected]