CC BY
94
С использованием предложенного алгоритма оценки качества системы наблюдения можно поставить обратную задачу обоснования структуры наблюдательной сети, а именно взаимного расположения источников выбросов и постов наблюдения. Целью такой задачи может быть максимизация вероятности обнаружения нарушения любым источником загрязнения.
Выводы, которые можно сделать из результатов имитационного моделирования, показывают возможность определения источников залповых выбросов. В ряде неблагоприятных случаев возможен пропуск события при существовании одной станции мониторинга, однако, наличие трех и более станций существенно уменьшает вероятность пропуска. Результаты имитационного моделирования создают основу для решения задачи обоснования структуры и параметров наблюдательной сети при заданных ограничениях (количество постов, вероятность обнаружения и т.д.)
Библиографический список
1. Овчаренко А.В., Белозеров А.Б., Березина С.В., Зобнин Б.Б., Уткин В.И., Черванев В.О. Отчет о научно-исследовательской работе РФФИ - Урал № 04-01-95096 «Динамические модели загрязнения воздушной среды в городской агломерации Екатеринбург-Первоуральск, Ревда, Дегтярск, Арамиль, Кольцово, Березовский, Среднеуральск». Екатеринбург: ИГФ УрО РАН, 2005. 51 с.
2. Овчаренко А.В., Белозеров А.Б., Березина С.В., Зобнин Б.Б., Уткин В.И., Черванев В.О. Методика создания динамической модели загрязнения воз-
душного бассейна Екатеринбурга // Геологическая интерпретация гравитационных, магнитных и электрических полей: Материалы 33 сессии Междунар. семинара им. Д.Г. Успенского. Екатеринбург: ИГФ УрО РАН, 2006. С. 264-267.
3. Зобнин Б.Б., Черванев В.О. Модель наблюдения совокупности природного и техногенного процессов // Глубинное строение, геодинамика, мониторинг, интерпретация геофизических полей: Материалы 3 науч. чтений Ю.П. Булашевича. Екатеринбург: ИГФ УрО РАН, 2005. С. 124-126.
УДК 533.6.011
ИССЛЕДОВАНИЕ И РАЗРАБОТКА ДВУМЕРНЫХ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ
ПУЛЬСИРУЮЩИХ ДЕТОНАЦИОННЫХ ДВИГАТЕЛЕЙ ДЛЯ ПЕРСПЕКТИВНЫХ ЛЕТАТЕЛЬНЫХ АППАРАТОВ
Г.П. Шиндяпин*, С.П. Шевырев, В.А. Поршнев
Саратовский государственный университет,
* кафедра вычислительного эксперимента в механике, кафедра прикладной информатики E-mail: [email protected]
В настоящей статье разрабатывается методика расчета параметров и характеристик пульсирующего детонационного двигателя с диффузором и эжектором на различных режимах и этапах его работы.
Ведется построение численных методов решения задачи по определению параметров течения продуктов детонации и воздуха в детонационном двигателе в различные моменты времени. Определяется суммарный импульс детонационного двигателя с диффузором и эжектором за один такт работы. Проведены пробные расчеты для математической модели детонационного двигателя с диффузором и эжектором на различных этапах его работы.
Research and Development of Bidimentional Mathematical Models of Pulsing Detonation Engines for Perspective Flying Devices
G.P. Shindyapin, S.P. Shevirev, V.A. Porshnev
In the present paper the design procedure of parameters and performances of the pulsing detonation engine with a diffuser and the ejector on different conditions and stages of its operation develops. Construction of numerical methods of a solution of a commitment by definition parameters of flow of yields of a knocking and air to the detonation engine in different instants is carried on. Determining of a general pulse of the detonation engine with a diffuser and the ejector for one clock tick of operation.
Trial accounts of a mathematical model of the detonation engine with a diffuser and the ejector at different stages of his operation are conducted.
Детонационные двигатели привлекают внимание разработчиков практически всех видов транспорта своей высокой производительностью. Создаются различные модели таких двигателей [1]-[3].
В [3] был дан обзор состояния исследований в разработке детонационных двигателей, построена математическая модель импульсного детонационного двигателя в квазидвумерной постановке, приведены результаты расчета параметров течения продуктов детонации и воздуха в случае наличия диффузора.
В настоящей статье исследования продолжены на случай плоской и осевой симметрии. Также добавлен эжектор, функция которого - засасывать окружающий воздух. Сложная геометрия иссле-
© Г.П. Шиндяпин, С.П. Шевырев, В.А. Поршнев, 2007
дуемого комплекса заставляет численно решать краевые задачи, зависящие от двух пространственных переменных и от времени.
Рассмотрены следующие задачи:
- разработка методики расчета параметров и характеристик пульсирующего детонационного двигателя с диффузором и эжектором на различных режимах и этапах его работы;
- разработка численных методов решения задачи по определению параметров течения продуктов детонации и воздуха в детонационном двигателе в различные моменты времени. Определение суммарного импульса детонационного двигателя с диффузором и эжектором за один такт работы;
- проведение пробных расчетов для математической модели детонационного двигателя с диффузором и эжектором на различных этапах его работы.
В первом разделе представлена методика расчета параметров и характеристик пульсирующего детонационного двигателя, опирающаяся на известные теоретические положения [4] и разработанные авторами математические модели [5], методы расчета [6] параметров детонационных двигателей. Методика расчёта содержит два этапа:
- этап распространения фронта детонации по камере (I этап). Этот этап исследован в [3];
- этап течения продуктов детонации и воздуха после завершения процесса детонации (II этап). Этот этап связан с распространением ударной волны и линии контактного разрыва в диффузорном и эжекторном участках и исследуется численным методом.
Методика позволяет рассчитать параметры и характеристики пульсирующего детонационного двигателя с диффузором и эжектором на различных режимах и этапах его работы, а также суммарный импульс давления на дно камеры и стенки диффузора и эжектора двигателя за один такт работы.
Во втором разделе разработан численный метод расчета параметров течения продуктов детонации и воздуха на II этапе работы двигателя после завершения процесса детонации. Сформулирована краевая задача для дифференциальных уравнений в частных производных, описывающих движение продуктов детонации и воздуха, решение которой строится численным методом Давыдова [7], [8], обобщенным на случай треугольной сетки [9].
Обсуждаются результаты расчетов при проведении вычислительного эксперимента по разработанной методике на этапах I и II. На всех этапах работы двигателя определяются: время протекания процесса; давление на дно камеры, стенки диффузора и эжектора; импульсы давления на дно камеры, стенки диффузора и эжектора в различные моменты времени, а также суммарный импульс для двигателя с диффузором и эжектором. Результаты вычислений иллюстрируются таблицами и графиками. Приведены сравнения с натурными экспериментами авторов. Обсуждаются физические закономерности и особенности работы детонационного двигателя.
1. АНАЛИЗ ОБЛАСТЕЙ ТЕЧЕНИЯ ГАЗОВ И ПРОДУКТОВ ДЕТОНАЦИИ В ДЕТОНАЦИОННОМ ДВИГАТЕЛЕ В РАЗЛИЧНЫЕ МОМЕНТЫ ВРЕМЕНИ. РЕЖИМЫ И ЭТАПЫ РАБОТЫ ДВИГАТЕЛЯ
1.1. Постановка задачи для импульсного детонационного двигателя с диффузором и эжектором
Детонационный двигатель представляет собой (рис. 1, а) цилиндрическую камеру постоянного сечения, заполненную при давлении, равном давлению окружающей среды, взрывчатым веществом, соединенную с диффузором и эжектором, заполненными при том же давлении нейтральным газом (воздухом).
Длина цилиндрической камеры 11, диаметр ¿1, длина диффузорного участка 10, угол полураствора диффузора а. Эжектор характеризуют длины 12, ¿3, диаметр ¿2, угол полураствора в. В камере находится взрывчатая смесь с давлением р1, плотностью р1. Скорость смеси = 0. Область камеры, где находится невозмущенная взрывчатая смесь, далее называется областью 1.
Давление воздуха в диффузорном участке р0, плотность р0, скорость ь0 = 0. Назовем эту область с невозмущенным воздухом областью 0.
Ось ОХ направлена вдоль оси камеры, а ее начало совпадает с дном камеры.
У дна камеры в момент времени £ = 0 инициируется детонационная волна (см. рис. 1), которая распространяется по взрывчатой смеси в режиме Чепмена-Жуге. После того, как в момент времени £4 фронт волны детонации закончит движение по взрывчатой смеси, по воздуху от сечения х = 11 пойдет ударная волна. Подвижная граница, разделяющая продукты детонации и воздух, будет представлять собой контактный разрыв, на котором плотность, температура и химический состав будут испытывать разрыв. Ударная волна, распространяясь по диффузорному участку, придет в сечения х = 11 + 10 и далее будет происходить вытекание воздуха из диффузора. Когда контактный разрыв достигнет сечения х = 11 + 10 + 12, начнется вытекание продуктов детонации из детонационного двигателя.
Области: 0 — невозмущенный воздух в диффузоре и эжекторе; 1 — невозмущенная взрывчатая смесь в камере; 2 — простая волна разрежения в продуктах детонации; 3 — область покоя в продуктах детонации; 4 — область разрежения в продуктах детонации; 5 — область воздуха, прошедшего через ударную волну
Фронты: Б — фронт волны детонации; — задний фронт волны разрежения; Б — фронт ударной волны; ^2 — фронт волны разрежения, движущийся в отрицательном направлении оси ОХ; С —линия контактного разрыва между воздухом и продуктами детонации
Время: £ = 0 — начало детонации; £ = £і — завершение детонации; £ = £2 — выхода ударной волны из диффузора; £ = — выхода контакт-
ного разрыва из двигателя; £ = £* — прихода волны разрежения в сечение х = 0; £ = —
при котором на выходе из дифузора скорость продуктов детонации обратиться в нуль
I этап
СІ1
І = 0
І = і і
ь 3 2 Т 1 і 1
0 < І < Іі
о
II этап
Задача состоит в том, чтобы определить изменение параметров течения продуктов детонации и газа в детонационном двигателе с течением времени и определить суммарный импульс тяги двигателя за счет действия давления продуктов детонации на дно его камеры, а затем воздуха и продуктов детонации на стенки диффузора и эжектора.
1.2. Области течения продуктов детонации и воздуха на этапе распространения фронта детонации по камере (I этап)
На рис. 1 изображены режимы Ь, c, d, e, f, g, И) и этапы (I, II) работы импульсного детонационного двигателя в различные моменты времени.
После того как от дна камеры начнет распространяться детонационная волна Б (рис. 1, Ь,
0 < £ < £1), в детонационном двигателе возникнут четыре области. Кроме областей покоя 0, 1 за фронтом детонационной волны Б возникнет область 2 — область простой волны разрежения в продуктах детонации. Задний фронт ^1 простой волны разрежения будет отделять область простой волны 2 от области покоя 3 и распространяться в положительном направлении оси ОХ. На рис. 1, Ь стрелками указаны направления распространения фронтов.
В момент времени £1 область 1 исчезнет (рис. 1, с), так как фронт детонации достигнет сечения х = ¿1, процесс детонации закончится. Завершится I этап - этап детонации.
І і < І < І2
І = І 2 f
І = ІС g
І = І*
1 р ^ А ■'""I |
3 2 ~с |4 _1 5 I 0 і I І I
—► ►- І I
І = Ік
Рис. 1. Режимы и этапы работы детонационного двигателя:
I - этап распространения фронта детонации по камере
II - этап течения продуктов детонации и воздуха после завершения детонации
a
d
є
h
1.3. Области течения продуктов детонации и воздуха после завершения процесса детонации (II этап)
В момент времени ¿1 в сечении х = ¿1 возникает произвольный разрыв, который сразу распадается (рис. 1, й). Возникает ударная волна 5, которая движется по воздуху в диффузорном участке двигателя, возникают контактный разрыв С, разделяющий продукты детонации и воздух, и волна разрежения ^2 в продуктах детонации, которая движется от контактного разрыва в сторону дна камеры.
При ¿1 < £ < £2 по воздуху будет распространяться ударная волна, за которой образуется область 5 воздуха, прошедшего через ударную волну. Далее движется контактный разрыв С, разделяющий воздух и продукты детонации. За контактным разрывом находится область 4 волны разрежения, которая возникает за счет распада произвольного разрыва. Фронт волны разрежения ^2 движется по области 2 в продуктах детонации в отрицательном направлении оси Ох.
Область 2 после момента времени £1 постепенно будет уменьшаться и затем в момент £ = £р исчезнет. В момент £ = £2 ударная волна Б достигнет выходного сечения диффузора х = 11 + , а
в момент £ = £т — выходного сечения двигателя х = 11 + ¿0 + ¿2, и воздух начнет вытекать через выходное сечение двигателя (рис. 1, е). Область 0 невозмущенного состояния воздуха исчезнет.
В момент времени £ = £3 (рис. 1, /) контактный разрыв С достигнет выходного сечения диффузора х = 11 + ¿0, а в момент £ = £с, сечения х = 11 + ¿0 + ¿2 и продукты детонации начнут вытекать через выходное сечение двигателя. Область 5 исчезнет.
На рис. 1, £ изображена картина течения при £ = £* > £с для момента времени, когда область 3 исчезнет, то есть фронт ^2 достигнет дна камеры. При £ > £* в камере, в диффузоре и эжекторе продукты детонации движутся в положительном направлении оси ОХ в сторону выходного сечения.
Расчет движения продуктов детонации при £ > £* (рис. 1, Н) будет производиться до момента времени £ = £к, при котором скорость продуктов детонации на выходе из сопла обратится в ноль.
1.4. Представление областей течения газов в детонационном двигателе
На рис. 2 изображены основные линии слабых и сильных разрывов и различные области течений в детонационном двигателе. Обозначения областей и линий разрывов соответствуют описаниям процессов п. 1.1—1.3, а также режимам и этапам работы детонационного двигателя (рис. 1, а-Н).
Этап I использует классическое решение для одномерного нестационарного случая [1] и полностью описан в статье [3].
2. ПОСТАНОВКА КРАЕВОЙ ЗАДАЧИ ДЛЯ РАСЧЕТА ТЕЧЕНИЯ НА ЭТАПЕ II
2.1. Системы дифференциальных уравнений в частных производных, описывающие течения на этапе II
Второй этап работы двигателя начинается с момента £1 , когда в сечении х = 11 возникает ударная волна Б, которая движется по воздуху и приводит его в движение, и контактный разрыв С, который разделяет продукты детонации, и воздух, прошедший через ударную волну. Контактный разрыв ведет себя как поршень. Он поджимает воздух и создает ударную волну и в то же время в продуктах детонации создает волну разрежения с фронтом ^2, движущимся в отрицательном направлении оси Ох (рис. 1, й).
Движение продуктов детонации в камере и диффузоре можно описать системой уравнений газовой динамики идеального газа в двумерной постановке (зависимость от двух пространственных переменных и от времени) [10].
Выпишем эти системы для случаев плоской и осевой симметрии.
Плоская симметрия:
др д (ри) д(рг’)
+
+
дх ду д (ри) д (рии) д (р + ри2)
д£
+
дх
+
ду
= 0,
= 0,
д (ри) д(р + ри2) д (рии)
д£ + дх + ду ’ д (ре) д ((е + р)и) д((е + р)г>)
д£
+
р = (к — 1)ре,
Е = е +
дх и2 + V2 2 '
+
ду
= 0,
(1)
Здесь £, х, у — независимые переменные, р — плотность, и, v — компоненты вектора скорости вдоль осей х и у соответственно, е — внутренняя энергия единицы массы газа, р — давление, Е — полная энергия единицы массы газа, к-отношение удельных теплоемкостей.
Осевая симметрия:
рv
др + д (ри) + д ^) д£ дх дг г
д ^) + д(р^) + д (р + рv2) рv2
д£ дх дг г
р = ( к - 1 ) ре,
д (ри) д (р + ри2) д(р^)
д£ + дх + дг д (ре) д ((е + р)и) д ((е + p)v) (е + р^
д£ + Е = е +
дх и2 + v2 2
+
дг
(2)
Здесь £, х, г — независимые переменные. Остальные параметры те же, что и в плоском случае.
2.2. Граничные условия
Решения систем (1), (2) должны удовлетворять граничным условиям: на дне камеры, на стенках камеры, на стенках диффузора, на эжекторе и в выходном сечении двигателя. На поверхности детонационного фронта, на контактном разрыве и на ударной волне граничные условия не ставятся, так как будет иметь место сквозной счет.
На дне, стенках камеры, а также на стенках эжектора выполняется граничное условие для скорости: нормальная к дну и стенкам компонента скорости равна нулю. То же самое выполняется на оси (плоскости) симметрии. На остальных внешних границах производные по нормали от всех интегрируемых функций равны нулю.
2.3. Начальные условия
Состояние продуктов детонации и воздуха в момент £ = ^ , т.е. в конце I этапа, служит начальным условием для II этапа. Начальные условия имеют следующий вид: v = Vl (х), р = р1 (х), р = р1(х) при £ = £1, 0 < х < 11; v = Vo, р = ро, р = ро при £ = £1, 11 < х < 11 +10.
Здесь v1 (х), р1 (х), р1(х) — функции распределения скоростей, плотностей и давлений в трубе в момент прихода фронта детонации Б в сечение х = 11. Эти функции определяются по известному аналитическому решению при £ = £1 из предыдущей статьи [3].
Отметим также, что хорошей и быстрой аппроксимацией полного начального условия является задание соответствующего энергоперепада между цилиндрической частью камеры и всем остальным пространством. Через небольшое число шагов по времени решение краевой задачи выходит на нужный режим.
2.4. Определение суммарного импульса давления на дно камеры, стенки диффузора и стенки эжектора
В результате численного решения краевой задачи, описанной ниже, в части 2 отчета, определяются: поля газодинамических параметров внутри и в окрестности детонационной камеры с диффузором и эжектором, давление р(0, у, £), действующее на дно камеры, распределение давления р на участке диффузора 11 < х < 11 + 10 и на участке эжектора, который составляет с детонационной камерой единое целое.
Теоретическое значение суммарного размерного импульса давления /С (за промежуток времени 0 < £ < Ьк) складывается из донного размерного импульса /1, дополнительного размерного импульса /2 за счет диффузорного участка и дополнительного размерного импульса /3 за счет эжекторного участка
/с = /1 + /2 + /3, (3)
г
где
tk di/2
/l =
FPi (y)(p(0,y,t) - po) dydt,
0 0
tk
/2 =
PD(x, y(x), t) — (x, y(x), t)FPD(x, y(x)) tan(a) dxdt,
(4)
(5)
0 li+lc tk
/з =
(х, у(х), £) — рэ(х, у(х), £)ЕРэ(х, у(х)) tan(в) йх^£. (6)
О /1+/0 —/3
Здесь ЕР^у) = 2пу, 0 < у < /2, ЕРи(х,у(х)) — периметр сечения канала на диффузорном
участке, ЕРэ(х, у(х)) — периметр сечения канала на эжекторном участке в точке х, ЕРд (х) = п(^1 + + 2(х — 11) tan(a)), 11 < х < 11 + 10, ЕРэ(х) = п(^2 + 2(11 + 10 — х) tan(в)), 11 + 10 — 13 < х < 11 + 10, рИ(х,у(х),£), р^(х,у(х) — давление на внутренней и внешней стенках диффузора, рэ(х, у(х),£), рВ(х,у(х),£) — давление на внутренней и внешней стенках эжектора.
Внутренние интегралы в (4), (5) и (6) вычисляются по методу трапеций, внешние интегралы вычисляются по методу левых прямоугольников.
3. ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ
3.1. Описание численного метода крупных частиц (метода Давыдова)
Численные методы в поставленных в работе задачах применяются по причине сложной нелинейной природы исходных систем уравнений, которые не поддаются аналитическому решению.
Статья [3] дает решение подобных задач при отсутствии эжектора как функцию от одной пространственной переменной х и от времени £. Этого достаточно, чтобы исследовать течение на первом этапе аналитически, а на втором этапе — численно с помощью метода Лакса-Вендроффа. Учет двух измерений (х, у) по пространству в настоящем исследовании заставляет применять более сложные системы уравнений и более мощные численные методы для их решения на II этапе. К числу таких методов относится метод крупных частиц (метод Давыдова) [7]-[10]. Расчетная область разбивается на простые геометрические фигуры (прямоугольники, треугольники и т.д.), прилегающие друг к другу. Газ, попавший в каждую такую фигуру, объявляется крупной частицей. Решение получается продвижением небольшими шагами по времени. Основной идеей этого метода является идея расщепления «по физическим факторам», когда сначала «работают» силы давления, затем происходит перетекание газа из одной крупной частицы в другую, и в конце происходит объединение результатов двух первых этапов в исходные уравнения, но уже на новом временном шаге.
Две системы (1), (2) будут численно интегрироваться, используя метод Давыдова (крупных частиц) [7], [8], [10], модифицированный на случай произвольной треугольной сетки [9].
Выпишем разностные схемы метода Давыдова для случая плоской симметрии. Для случая осевой симметрии разностные схемы аналогичны.
Расчетное поле с помощью алгоритма триангуляции разбивается на множество треугольников. Достоинство треугольной сетки — в ее быстрой адаптации к произвольным границам.
Все вычисляемые газодинамические параметры (плотность, скорость, полная энергия, давление) относятся к геометрическим центрам треугольников.
Расчет одного временного шага разбивается на три этапа. На первом, эйлеровом этапе, пренебрегают дивергентными членами. Укороченные дифференциальные уравнения эйлерова этапа
ди др
Ж + дх = 0 | + др = о.
д£ ду
дЕ дри др-и _ д£ + дх + ду
дискретизируются на треугольной сетке. В частности, члены с давлением можно заменить через контурные интегралы (применяется теорема Грина и теорема о среднем [11]):
dp _ fA pdy
dx §A x dy ’
dp
dy
§A pdx
fA x dy ’
dpu + dpv _ §A pu dy — §A pv dx dx dy §A x dy
Контурные интегралы при численных расчетах обходятся по сторонам треугольника против часовой стрелки:
3
E(pj+i + Pj )(у^+і - yj) dp _ j=i_________________________
dx з
E(xj + i + xj )(yj + i - yj) j=1
-E (Pj+i + Pj )(xj+i - x) dp j=i
dy
E(xj+i + xj )(yj+i - yj) j=i
E {(Pj+iuj+i + Pju)(yj+i- yj)- (Pj+ivj+i + Pjvj)(xj+i- xj)}
dpu dpv j=i
dx dy 3
y E (Xj+i + Xj)(yj+i - yj)
j=i
Здесь p4 = рь x4 = xi; y4 = yi; pi; p2, p3 — значения давления в трех соседних треугольниках;
(xi,yi), (x2,y2), (x3,y3) — координаты центров этих треугольников.
Внешние границы и границы на теле заменяется виртуальными ячейками, чтобы сохранить единую схему расчета. Все параметры в виртуальных ячейках, кроме скорости, заменяются аналогичными параметрами из соседних реальных ячеек, нормальная компонента скорости меняет знак на противоположный на стенках.
После аппроксимации производных с давлением величины с «волнами» определяются обычным образом, как это принято в классическом методе Давыдова [7], [8]:
иП = иП - р(dp/dx)n, pi
vn = vn - pf (dp/dy)n, (7)
pi
Ef = Ef------f [(dpu/dx)n + (dpv/dy)”] .
Pi
На лагранжевом этапе учитывают то, что было опущено на эйлеровом этапе, то есть члены, описывающие перенос массы, импульса, энергии.
Выберем произвольную сторону треугольника, например, отрезок (ж,уі) — (ж2, у2). Выпишем формулу для потока массы через выбранную границу:
IIП I 7/П
Л,™ n »(норм) + 7к1,2(норм) .. w ^
AMK1j2 = Рр---------- -----о--------- AtL((xi,yi), (х2,У2)),
2
(8)
где АМ_П12 — поток массы через сторону (^1 ,у1) — (х2,у2); рп — плотность в текущем треугольнике, если поток вытекает из него, и плотность в соседнем через указанную сторону треугольнике, если поток втекает; иП(норм) — нормальная составляющая вектора скорости для текущего треугольника; 2(норм) — нормальная составляющая вектора скорости для соседнего треугольника;
Е((х1 ,У1), (Х2,У2)) = л/(х1 — х2)2 — (У1 — У2)2 — длина стороны (Х1,У1) — (Х2,У2); А£ — шаг по вре-
мени.
На заключительном этапе
pl,+i = рП + !>m;/s,,
j=i
U n+i ___ і ~n ^n
(E” pn + £am“E”;/s,)/p
n+i
j=i
(9)
^n+i _ (et pn
(e” pn + £АМ"5”/$)/p
n+i
j=i
E,n+i _ (E”p” + £ amj>.e;/s,)/p:
j=i
n+i
где — площадь текущего (¿-того) треугольника, р — нижний индекс у переменных ип, ирП, -Ерп, равен номеру текущего треугольника, если поток вытекает из него через указанную сторону, в противном
случае, этот индекс равен номеру соседнего треугольника с общей указанной стороной. В процессе триангуляции номера соседних треугольников не обязательно образуют упорядоченную последовательность, поэтому, хотя задача двумерная по пространству, целесообразно применять векторную запись для хранения искомых газодинамических параметров.
На теле ставится условие непротекания:
гр.(норм) _ ивирт.(норм) •
(10)
На внешних границах ставятся обычные условия «мягкого» сопряжения (равенство нулю производных от интегрируемых функций по нормали к границе).
Таким образом, формулы (7)-(10) представляют метод Давыдова для случая произвольной треугольной сетки.
Создано три комплекса программ, один из которых выстраивает геометрию конкретного обтекаемого тела, второй — численно решает краевые задачи, третий — обрабатывает результаты (изолинии, параметры на обтекаемом теле и т.д.)
На рис. 3 и в таблице представлены результаты вычислений превышения донного давления над давлением окружающей среды во времени (0 < Ь < Ьь). Время завершения решения краевой задачи Ьь. Угол а = 15°, длина диффузорного участка 1о = 0.050 м, длина цилиндрического участка ¿1 = 0.450 м; длина эжекторного участка 13 = 0.400 м; угол в = 15°, Ь2 = 120.
1000 * t
Рис. 3. Превышение донного давления
Время Давление
0.0 12.512
0.1 11.901
0.2 11.788
0.3 8.213
0.4 4.110
0.5 2.002
0.6 1.015
0.7 0.478
0.8 0.311
0.9 0.199
1.0 0.256
1.1 0.214
1.2 0.161
1.3 0.144
1.4 0.123
1.5 0.100
и
На рис. 4. представлены зависимости суммарного интегрального импульса 1С от времени: 1 — без эжектора; 2 — с эжектором; 3 — отдельные компоненты суммарного интегрального импульса: 3.1 — донный импульс, 3.2 — диффузорный импульс, 3.3 — эжекторный импульс.
Рис. 4. Размерные интегральные импульсы
При сравнении результатов рис. 4 с результатами статьи [3] эжектор повышает это отношение на 19%, причем в начальные моменты времени — отношение меньше (кривые 1 и 2), чем без эжектора (через эжектор уходит газ за счет перепада давления), зато в конце расчета эжектор начинает выполнять свою основную функцию — засасывать воздух из окружающей среды.
На рис. 5 представлены поля давления в детонационном двигателе в фиксированные моменты времени £ = 0.10£&, 0.25£&, 1.00£&. Как видно из приведенных результатов, поток существенно перестраивается при изменении времени. Если в начальные моменты времени (£ = 0.1£&, 0.25£&) масса газа выбрасывается через эжектор, то в последующие моменты (£ > 0.5£к) воздух из окружающей среды засасывается через эжектор, что создает дополнительный положительный суммарный импульс /3 эжектора (для нашего варианта расчёта 19% при £ = £&). Ввиду наличия оси симметрии расчёты велись для половины камеры.
Рис. 5. Поля давления в детонационном двигателе в различные моменты времени: а — Ь = 0.10 Ьк, Ь — Ь = 0.2Ык, с — Ь = 1.00Ьк
3.2. Физические закономерности, отмеченные в результате расчетов
В результате аналитического (этап I) и численного (этап II) решения были получены распределения параметров продуктов детонации и воздуха после завершения детонации (£ = £1), определены при численном решении положение контактного разрыва и ударной волны в фиксированные моменты времени (для выделенных номеров итераций), а также значения давления продуктов детонации на дно камеры, донный и суммарный размерный импульсы. Определены зависимости этих импульсов по времени за один цикл работы двигателя, характеризующие работу отдельных элементов (цилиндрической, диффузорной и эжекторной частей) двигателя.
Вычисление каждого варианта заканчивалось, как только значение скорости в выходном сечении диффузора становилось отрицательным.
На рис. 3, иллюстрирующем превышение давления в донной части над давлением окружающей среды (в расчетах — атмосферным), начальный участок превышения представляет собой слабо изменяющуюся линию, так как соответствует области покоя продуктов детонации (с постоянным давлением). После прихода в эту область волны разрежения давление у дна камеры начинает падать и становится ниже давления окружающей среды.
Как видно из графика для полного интегрального импульса, а также составляющих полный импульс донного, диффузорного и эжекторного импульсов (рис. 4), он имеет качественно тот же вид, что и в одномерном нестационарном случае. Добавление третьего слагаемого неоднозначно влияет на вид графика, хотя в конечные моменты времени, превышающие 0.5£к расчета, эжектор начинает выполнять свою основную функцию — засасывать воздух из окружающей среды.
На рис. 5 изображены поля давления р/р0 в три момента времени (£ = 0. 10£к, 0.25£к, 1.00£&). Первые два момента иллюстрируют отрицательное действие эжектора, когда избыточное давление выдавливает газ наружу. На рис. 5, с, газ начинает засасываться.
На рис. 6 приведены экспериментальные результаты авторов (снимки, запечатлевшие моменты натурного эксперимента с макетом детонационного двигателя в режиме выхлопа), показывающие удовлетворительное совпадение численного и натурного экспериментов.
!
*•
Рис. 6. Результаты натурного эксперимента
ЗАКЛЮЧЕНИЕ
В статье приведены результаты теоретических и экспериментальных исследований по разработке математических моделей пульсирующих детонационных двигателей, представляющих значительный практический интерес для решения задач авиации, космонавтики, транспорта.
Для двигателя, представляющего цилиндрическую камеру с диффузором и эжектором, разработана методика расчета параметров и характеристик пульсирующего детонационного двигателя на различных режимах и этапах его работы.
На I этапе работы (этап детонации) построено аналитическое решение, на II этапе - численное решение задачи. Составлен комплекс программ для расчета течения продуктов детонации и воздуха, в результате работы которого определены параметры течения (поля давления и чисел Маха) в различные моменты времени, составляющие импульса давления, характеризующие давления в донной части, в эжекторе и диффузоре; суммарный импульс давления детонационного двигателя.
Таким образом, в работе разработаны методика расчета и математическая модель пульсирующего детонационного двигателя, позволяющие исследовать закономерности и особенности детонационных процессов. Результаты могут быть использованы при разработке и проектировании лабораторных образцов пульсирующих детонационных двигателей.
Библиографический список
1. Ремеев Н.Х. и др. Численное моделирование и экспериментальное исследование рабочего процесса в детонационной камере сгорания // Химическая физика. 2003. Т. 22, № 8. С. 45-56.
2. Зубов В.И. и др. Численное моделирование пространственных турбулентных течений газа в сложных сопловых устройствах // Журн. выч. мат. и мат. физ. 2005. Т. 45, № 10. С. 1871-1885.
3. Шиндяпин Г.П., Сапунков Я.Г., Шевырев С.П., Поршнев В.А., Федорец Н.В., Федорец В.Н. Проблемы математического моделирования детонационных двигателей для перспективных летательных аппаратов // Изв. Сарат. ун-та. Новая серия. Математика. Механика. Информатика. 2006. Т. 6, вып. 1/2. С. 74-85.
4. Станюкович К.П. Неустановившиеся движения сплошной среды. 2-е изд. М.: Наука, 1971.
5. Сапунков Я.Г., Шиндяпин Г.П., Поршнев В.А., Федорец B.Н. Математическая модель детонационного двигателя // Математика. Механика: Сб. науч. тр. Са-
ратов: Из-во Сарат. гос. ун-та, 2001. Вып. 3. С. 178-181.
6. Сапунков Я.Г., Шиндяпин Г.П., Поршнев В.А., Федорец В.Н. Расчет гидродинамических параметров в импульсном детонационном двигателе с диффузором // Тяжелое машиностроение. М., 2003. № 11. С. 7-10.
7. Математическая энциклопедия. Т. 3. М.: Сов. эн-цикл., 1982.
8. Математический энциклопедический словарь. М.: Сов. энцикл., 1988.
9. Шевырев С.П. Разностные схемы метода Давыдова на произвольной сетке // Математика. Механика: Сб. науч. тр. Саратов: Изд-во Сарат. ун-та, 2005. С. 205209.
10. Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Наука. Гл. ред. физ.-мат. лит., 1984.
11. Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления. М.: Наука. Гл. ред. физ.-мат. лит., 1970. Т. III.
